Ordering and Dynamics for the Formation of Two-Dimensional

Article pubs.acs.org/JPCC

Ordering and Dynamics for the Formation of Two-Dimensional Molecular Crystals on Black Phosphorene Titas Kumar Mukhopadhyay and Ayan Datta* Department of Spectroscopy, Indian Association for the Cultivation of Science, Jadavpur Kolkata 700032, West Bengal, India S Supporting Information *

ABSTRACT: Phosphorene, a monolayer of elemental phosphorus, has emerged as the most important 2D atomic crystal post graphene. As a consequence of its suitable direct band gap, it has high carrier mobility, thereby making it suitable for potential applications in nanoelectronic devices. In this article, we have performed molecular dynamics simulations for the first time to unveil the potential of black phosphorene (BP) as a platform to fabricate two-dimensional crystals through van der Waals epitaxy of organic molecules. Eight molecules (benzene, 1,3,5-trimethylbenzene, 1,3,5-trifluorobenzene, 1,3,5-trihydroxybenzene, 7,7′,8,8′-tetracyanoquinodimethane (TCNQ), pentacene, coronene, and fullerene) of various size, shape, and polarity are shown to form highly periodic 2D monolayers on BP. We show that self-assembly proceeds through a hierarchy of steps, beginning with adsorption on the surface followed by interplay between the kinetic energy and intermolecular interactions. Crystalline 2D layers are obtained through the minimization of van der Waals interactions for nonpolar and electrostatic interactions for polar molecules. Long-range attractive electrostatic interactions contribute significantly toward the stabilities of the superlattices, thereby making the melting points of polar crystals much higher compared to their nonpolar analogs (Tm (polar) > Tm (nonpolar)). Interestingly, although nonpolar fullerene is shown to form highly stable superlattices due to extensive intermolecular van der Waals interactions, the underlying BP monolayer helps to orient the molecules and, hence, indirectly assists molecular ordering. It also contributes to the overall stability of the molecular crystals through interfacial adsorption. A comparison of self-assembly patterns over BP vis-a-vis graphene and hexagonal boron nitride shows that the nonaromatic nature and undulation of black phosphorene do not affect the intrinsic stability of the assemblies. These findings offer new avenues for generating new molecular doped crystalline superlattices for semiconductor industry.

■

gained remarkable research interest.9,10 Additionally, they provide atomically thin surfaces which have made them one of the most promising playgrounds to study 2D thin-film formation owing to the potential applications in optoelectronic devices such as organic light-emitting diodes (OLEDs), organic photovoltaics (OPVs), and organic field effect transistors (OFETs).11 Prototypical among them, epitaxial graphene supported on metal surfaces such as Ir(111) and Ru(0001) has been extensively studied for the formation of organic self− assembled monolayers because of its chemical inertness which strongly reduces the influence of the metal underneath.12,13 In the past few years, a large variety of molecules such as 3,4,9,10perylene tetracarboxylic dianhydride (PTCDA), perylene tetracarboxylic diimide (PTCDI), tetracyanoquinodimethane (TCNQ), tetrafluoro-tetracyanoquinodimethane (F4TCNQ), pentacene, copper phthalocyanine, fullerene, dioctylbenzothienobenzothiophene (C8-BTBT), long-chain alkanes such as CnH2n+2 (n = 8, 16, 24, 34), haloalkanes, cholesterol and cholesteryl esters have been assembled into layered twodimensional crystalline state as monolayers on graphene.14−24

INTRODUCTION Molecular self-assembly is a process in which randomly placed molecules initially in a disordered state organize into distinctly patterned aggregates. Apart from the ceaseless quest to decipher the mechanism behind such order out of chaos, selfassembly of nano- and mesoscale objects such as organic πconjugated systems, quantum dots, nanotubes, nanowires, and nanoscrolls are of the utmost importance as they have potentials to perform specific functions. In particular, fabrication of molecular ordered and superordered architectures on solid surfaces has become a major challenge in the realm of nanotechnology and semiconductor industry.1−5 Traditionally, self-assembled 2D networks of organic molecules are prepared over metal surfaces such as meso-tetramesitylporphyrin over Cu(100), phthalocyanin, and its derivatives (e.g., subphthalocyanin) on Ag(111) and Cu(111) and phenyl-C61-butyric acid on Au(111) among many others.6−8 However, these metallic surfaces strongly interact with the adsorbed molecules, and molecular ordering is largely dominated by surface−molecule interactions. Therefore, there has been an urge to explore materials which can efficiently decouple the interactions between metallic substrate and molecules. In recent years, the family of two-dimensional materials with their intriguing thermal, structural, and electronic properties has © 2017 American Chemical Society

Received: March 16, 2017 Revised: April 24, 2017 Published: April 25, 2017 10210

DOI: 10.1021/acs.jpcc.7b02480 J. Phys. Chem. C 2017, 121, 10210−10223

Article

The Journal of Physical Chemistry C

form stable 2D periodic networks over BP. The mechanism of self-assembly is investigated and analyzed for both polar and nonpolar molecules. Extensive free energy calculations are performed to determine the thermodynamic preference of 2D lattice formation. We also investigate the temperature-dependent stability of the lattices and estimated their formation as well as melting temperatures. The effect of the presence of undulation and lack of conjugation in phosphorene toward molecular ordering is compared with graphene as well as hexagonal boron nitride. Multilayer surface effects are demonstrated to have a minor influence on the assembly process.

Other 2D materials such as hexagonal boron nitride and molybdenum disulfide have also been explored to fabricate 2D molecular crystals.23,25,26 One of the latest additions to the family of two-dimensional materials is atomically thin layers of phosphorus, known as phosphorene.27 Elemental phosphorus is known to exist in several polymorphic forms, among which the white, red, and black forms are widely known. The structure of black phosphorene (BP) is unique, consisting of hexagonal units of phosphorus atoms situated in the same plane, but each of these atoms is covalently bound to three neighboring other P atoms. This results in a puckered shape of the P6 rings with out-ofplane ridges, unlike graphene.28 As with other van der Waals solids, intralayer binding in black phosphorus is dominated by weak dispersion forces which offer exciting strategies to produce mono- to few layers of BP through mechanical or liquid phase exfoliation (LPE) techniques.29−31 In 2014, Zhang and co-workers fabricated a BP field effect transistor which displayed charge carrier mobility as high as 103 cm2 V−1 s−1.32 Koenig et al. demonstrated a few-layer BP FET on Si/SiO2 with a drain current modulation of 103 and low-temperature ON− OFF ratio exceeding 105.33 Subsequently, theoretical investigations have unraveled some remarkable mechanical and electronic properties of this material. For instance, in contrast to graphene, a monolayer of BP behaves as a p-type semiconductor possessing a direct bandgap of 2.0 eV, which reduces with increasing number of layers, converging to 0.31− 0.36 eV for the bulk phase.34,35 Li et al. demonstrated that edge-hydrogenated black phosphorene nanoribbons can also behave as direct band gap semiconductors.36 In addition, BP undergoes a series of phase transitions and even behaves as a superconductor at high pressures.37,38 It has been predicted that a deforming strain normal to the plane of BP can reduce the bandgap, thereby inducing a semiconductor−metal transition.39 Moving one step ahead, researchers have suggested molecular doping of BP can be an effective way to tune the electronic properties as has been established for other 2D materials.40−42 Recently, the interactions between BP and donor (tetrathiofulvalene (TTF) and its derivatives) as well as acceptor molecules (tetracyanoquinodimethane (TCNQ) and tetracyanoethylene (TCNE)) were investigated through density functional theory (DFT) which reveal that the band gap of BP is strongly reduced because of the charge transfer between BP and the dopant molecules. BP doped with strong donors such as tetraamino-tetrathiofulvalene is found to behave as an n-type semiconductor, while acceptor dopants are capable of converting BP into a p-type semiconductor.43−45 Clearly, doping BP through an epitaxial growth of noncovalently bound organic molecules is of seminal importance for designing BP-based materials. However, in spite of several experimental studies of molecular ordering over graphene, the same on the BP surface is yet to be achieved. A molecular level understanding is essential to scrutinize the potential of BP to act as a platform for the growth of molecular patterns. Molecular dynamics (MD) simulations can provide in-depth qualitative and quantitative information about these molecularscale events. Recently, it has been employed to study the selfassembly of small organic molecules and fullerenes over graphene as well as graphyne.46−48 However, to the best of our knowledge, there are no previous reports for the study of molecular self-assembly over BP and the factors involved therein. In this paper, eight different organic molecules with varying size, shape, polarity, and functional groups are shown to

■

COMPUTATIONAL DETAILS To study the self-assembly of small organic molecules on the BP surface, we have taken a BP monolayer having the dimension of 10 nm × 10 nm in the XY plane, which was placed at the center of a simulation box of dimension 10 × 10 × 20 nm3 with periodic boundary conditions applied in all three directions. We performed simulations with different surface coverage for each molecule. However, the results presented here correspond to intermediate surface coverage (50−60%) where sufficient space is available for the molecules to undergo thermal motion while avoiding overlap with periodic images. Each system was subject to energy minimization for 104 steps using the conjugate gradient method to eliminate unfavorable contacts between atoms followed by equilibration for 5 ns using the canonical (NVT) ensemble at 50 K. Following this, each system was heated to the desired temperature at a heating rate of 50 K/ns, and production simulations were performed for each system using the NVT ensemble at least for 200 ns. Nevertheless, most of the results presented in this article are up to 100 ns, indicating that no significant changes were observed after that instant of time. The 2D lattice parameters were calculated for the last 10 ns of our 200 ns production simulations to get a high quality ordered structure. During all the simulations each atom of the BP sheet was harmonically constrained with a force constant of 100 kcal/mol/Å2. For each system, we have taken three different initial configurations consisting of the same number of substrate molecules. The simulation protocol mentioned above was then applied to all the initial configurations, and two parallel simulations were run (total of six simulations per system) starting from the same configuration, for each system. Each structure produced very similar 2D ordering providing us the necessary confidence on the quality of our simulations. We have used Langevin dynamics with a damping coefficient of 5 ps−1 to maintain isothermal conditions in all of our simulations.49 Full periodic electrostatic interactions of the system were calculated engaging the particle mesh Ewald method with 1 Å grid.50 Classical equations of motions were integrated using the Velocity Verlet algorithm, and a 2 fs time step was used for the propagation of the system. SHAKE algorithm was used to hold rigid covalent bonds involving hydrogen atoms.51 A switching function was applied for the calculation of all Lennard-Jones interactions within the distance range 10−12 Å, and the cutoff for nonbonded interactions was set to 12 Å. Atomic coordinates were stored after every 2 ps for the trajectory analyses. We have used Packmol utility for the preparation of initial configurations.52 NAMD 2.10 was used for all MD simulations and VMD for visualization, and our in house Tcl scripts as well as VMD plug-ins were used for the analysis of data.53,54 10211


Article


Table 1. Unit Cell Parameter, Formation Temperature (Tf), Melting Temperature (Tm), ΔGfree, ΔGassembly, and ΔΔGrel assembly of the 2D Lattices for the Molecules 1−8 over Black Phosphorene (BP)a entry

molecule

surface

1

benzene

BP

2

1,3,5-trimethylbenzene

BP

3

1,3,5-trifluorobenzene

BP

4

1,3,5-trihydroxylbenzene

BP

5

TCNQ

BP

BP (trilayer) graphene

h-BN

6

pentacene

BP

BP (trilayer) graphene

h-BN

a

7

coronene

BP

8

fullerene

BP

unit cell parameters a = b = 6.51 ± 0.26 Å α = 121.8° ± 4.6° β = 60.2° ± 3.0° a = 8.69 ± 0.25 Å b = 7.05 ± 0.15 Å α = 108.0° ± 2.4° β = 71.9° ± 2.1° a = b = 7.11 ± 0.11 Å α = 116.2° ± 3.4° β = 65.1° ± 2.1° a = b = 7.34 ± 0.18 Å α = 119.6° ± 3.3° β = 60.2° ± 3.0° a = 12.44 ± 0.28 Å b = 9.18 ± 0.20 Å α = 133.0° ± 2.5° β = 47.8° ± 2.0° similar to above a = 12.40 ± 0.37 Å b = 9.20 ± 0.24 Å α = 132.9° ± 2.1° β = 47.2° ± 1.7° a = 12.40 ± 0.18 Å b = 9.21 ± 0.21 Å α = 133.1° ± 0.9° β = 47.7° ± 2.5° a = 16.50 ± 0.20 Å b = 7.11 ± 0.11 Å α = 104.3° ± 1.9° β = 76.1° ± 1.8° similar to above a = 16.54 ± 0.05 Å b = 8.10 ± 0.15 Å α = 116.8° ± 1.0° β = 64.0° ± 2.3° a = 16.51 ± 0.15 Å b = 8.20 ± 0.09 Å α = 118.6° ± 1.2° β = 61.2° ± 1.1° a = b = 11.55 ± 0.14 Å α = 120.6° ± 1.2° β = 59.4° ± 1.0° a = b = 10.14 ± 0.21 Å α = 121.7° ± 2.7° β = 58.4° ± 1.1°

Tf (K)

Tm (K)

ΔGfree (kcal/mol)

ΔGassembly (kcal/mol)

ΔΔGrel assembly (kcal/mol)

100

140

−11.0

−13.3

−2.3

120

180

−16.3

−18.5

−2.2

120

180

−13.8

−17.4

−3.6

400

500

−12.8

−23.2

−10.4

350

500

−24.2

−35.0

−10.8

350 350

500 -

−21.7

−36.0 −34.3

−12.6

350

-

−25.5

−36.1

−10.6

120

180

−33.6

−40.4

−6.8

120 120

-

−31.3

−40.7 −37.5

−6.2

120

-

−36.2

−43.8

−7.6

120

160

−36.9

−42.6

−5.7

450

1300

−23.3

−48.7

−25.4

For pentacene and TCNQ, these values are also compared with respect to the 2D lattices formed over graphene and h-BN.

Potential of Mean Force (PMF) Calculations. The adaptive biasing force (ABF) method, as implemented in NAMD, was employed to determine the free energies in terms of potential of mean force (PMF).55 In this method, a locally estimated external biasing force is applied to a set of atoms along a proper reaction coordinate, so that the system can overcome sufficiently large potential energy barriers. Herein ΔGassembly and ΔGfree represent the potential of mean force (PMF) to detach a molecule vertically (along z axis) from the phosphorene surface, from the assembly, or free, respectively (Figure 6). For the calculation of ΔGfree and ΔGassembly, the

vertical distance (dz) between the center of mass of the BP surface and a single molecule was chosen as the reaction coordinates which ranged from 0 to 1.5 nm. During the calculation of ΔGfree, a single molecule was placed over the center of the surface, and for ΔGassembly, a molecule is chosen from the center of the already formed molecular assembly over BP to allow maximal molecule−molecule interactions. The reaction coordinates were divided into several equally spaced windows, each having width of 0.1 nm. Further, each of the windows was divided into small bins of 0.02 nm width, and a biasing force was applied once 5000 force samples were 10212


Article


Figure 1. Snapshots from the simulation of self-assembly of 136 benzene molecules at 100 K over a black phosphorene surface representing different stages of 2D lattice formation. (a) Initial structure after adsorption, (b) formation of self-assembled islands (shown by red circles), (c) growth of islands, (d) joining of the islands, (e) complete self-assembly of benzene, and (f) unit cell parameters of the 2D lattice of benzene, as obtained from the simulations.

respectively, according to Lorentz−Berthelot combining rules. We have used CHARMM General Force Field (CGenFF) parameters (version 3.0.1) to model all the aromatic molecules with zero penalty except for TCNQ that had a high (>50) CGenFF charge penalty.56 Therefore, we evaluated the CHelpG charges of TCNQ at the MP2/6-31G+(d,p) level of theory and simplified them by analogy in accordance with the parametrization protocol used for the generation of CGenFF.56 Phosphorus atoms in BP were assumed as uncharged LennardJones spheres, and the parameters as used by Blankschtein et al. were considered.57 Further details of heating simulations, parameters, and parameter modification protocols can be found in the Supporting Information.

collected for each bin. Finally, for each window, a 3 ns production ABF simulation was performed using an NVT ensemble at the same temperature used for the formation of self-assembly. Force-Field Parameters. The potential energy function used in the simulations has the mathematical form U (b , θ , χ , ϕ , rij) =

∑

Kb(b − b0)2 +

bonds

+

∑

K χ (1 + cos(nχ − δ)) +

dihedrals

+

∑

K θ(θ − θ0)2

angles

∑

Kϕ(ϕ − ϕ0)2

impropers

⎡⎛ min ⎞12 ⎛ R min ⎞6 ⎤ qq R ij ⎟ − 2⎜ ij ⎟ ⎥ + i j εij⎢⎢⎜⎜ ⎟ ⎜ ⎟ ⎥ r ϵrij ⎝ rij ⎠ ⎦ nonbonded ⎣⎝ ij ⎠

∑

■

RESULTS AND DISCUSSION To understand the molecular scale of events occurring during the epitaxial growth of organic thin films over BP, we have investigated the 2D ordering of eight different molecules with varying chemical and physical properties. Benzene (1), because of its simple structure and nonpolar nature, is considered as a benchmark to explore the mechanism of self-assembly. To compare and contrast with respect to benzene, we studied the self-assembly of three trisubstituted aromatic molecules, namely, 1,3,5-trimethylbenzene (mesitylene) (2), 1,3,5-trifluorobenzene (3), and 1,3,5-trihydroxybenzene (4), each possessing a 3-fold rotational symmetry. In addition, the selfassemblies of 7,7′,8,8′-tetracyanoquinodimethane (TCNQ)

Kb, Kθ, Kχ, and Kφ are force constants associated with bond stretching, angle bending, torsional angle, and improper angles, respectively. The subscript 0 denotes the equilibrium values of bond (b), angle (θ), and improper (φ), while n determines the periodicity of the dihedral potential in the interval [0,2π]. The nonbonded interactions were calculated as a sum of Lennard− Jones (LJ) “6-12” and Columbic terms wherein ϵ represents the effective dielectric constant; qi and qj are the partial atomic charges; rij is the distance between atoms i and j; Rijmin is the distance at the LJ minimum; and εij is the LJ well depth. The g eo m e t r ic m e an εij = εiεj a n d a r i t h m et i c m e a n R ijmin =

R imin + R jmin 2

were used for calculating εij and Rijmin, 10213


Article


Figure 2. Time evolution of the (a) radius of gyration of whole assembly and the orientational order parameter (S) of benzene (inset represents change in these quantities within the first 7 ns), (b) electrostatic, and (c) van der Waals interaction energies (per molecule) between benzene molecules, (d) van der Waals interaction energy (per molecule) between benzene and the black phosphorene surface, at 100 K. Note that in (b−d) the bold line corresponds to running averages of energies with an average length of 200 data points.

or the number density profile along the x-axis (Figure S2). The

(5), pentacene (6), coronene (7), and fullerene (8) were also simulated, all of which have been previously shown to form thin films over graphene or other substrates.16,17,21 Molecules are classified into two categories, polar (TCNQ and 1,3,5trihydroxybenzene) and predominantly nonpolar (all others). However, here we will discuss the results taking one molecule from each category, and the rest are discussed in the Supporting Information. Figure S1 displays the unit cell structures, and Table 1 shows the unit cell parameters, formation, and melting temperatures of the eight 2D lattices over BP as obtained from our simulations after 200 ns, which are in excellent agreement with previous results over graphene as well as other substrates.2,16,17,48 General Mechanism of Self-Assembly and 2D Lattice Formation. Predominantly Nonpolar Molecules. Initially, 136 benzene molecules were placed randomly over a monolayer (10 nm × 10 nm) of BP at 100 K, and this disordered state was converted completely into a twodimensional periodic network. Figures 1(a)−1(e) show the snapshots of the simulation as viewed from the top of the phosphorene surface at various time intervals. The entire process of self-assembly can be divided into several distinct events which can be monitored through the time evolution of the radius of gyration (Rg) of the entire assembly (Figure 2(a))

radius of gyration is calculated as R g =

1 n

n

2 ∑i = 0 |ri2 − rcom |,

where n is the total number of atoms, and ri and rcom are the positions of the center of mass of the i-th atom and the whole assembly, respectively, at any instant of time. On the other hand, phase transition into an ordered arrangement is realized through the orientational order parameter S = 3/2⟨cos2 θ⟩ − 1/ 2 (Figure 2(a)), where S depicts the extent of order in the system and θ is the angle between the axis of a single benzene molecule and the preferred orientation of the molecules in the crystal. During energy minimization and equilibration periods, a significant fraction of the benzene molecules get adsorbed parallel to the BP surface. Between 0 and 20 ps of the production simulation, adsorption continues followed by random thermal movement resulting in collision with other molecules in the neighborhood. Following such collisions, molecules get attracted by weak intermolecular forces giving rise to several self-assembled nucleation islands (Figure 1b, red circles). This is characterized by a sharp decrease in Rg from 4 to 3.6 nm within 1.5 ns in Figure 2a (inset). Subsequently, these islands grow very fast as more and more molecules are added up through collision which is illustrated by another fall in Rg from 3.6 to 3.4 nm between 1.5 and 7 ns. Interestingly, the 10214


Article


Figure 3. (a−d) Snapshots representing structural changes observed through the modification of force-field parameters of the benzene− phosphorene system. Note that, Em−m and Em−s represent molecule−molecule and molecule−surface interactions, respectively.

Figure 4. (a−e) Snapshots from the simulation of self-assembly and 2D lattice formation of 75 TCNQ molecules over the black phosphorene surface at 350 K at different time instants. (f) 2D lattice parameters of TCNQ as obtained from our simulation at the end of 200 ns.

10215


Article


Figure 5. Time evolution of the (a) radius of gyration (Rg) and the order parameter (S) of the TCNQ assembly, (b) electrostatic and (c) van der Waals interaction energies (per molecule) between TCNQ molecules, and (d) van der Waals interaction energy (per molecule) between TCNQ and the black phosphorene surface, at 350 K. (e) Average charge density profiles along the x axis at different time instants showing increasing periodicity of charges.

Figure 2(a) inset also shows that, within the first 7 ns, the order parameter S rapidly increases up to 0.9, and interestingly, in this period the time evolutions of S and Rg are almost mirror images. Thus, the system rapidly attains periodicity through the formation of nucleation islands which behave as separate crystalline networks. Eventually, these islands join altogether (Figure 1d) yielding a complete assembled phase at about 60 ns (Figure 1e). This is reflected in the gradual decrease in Rg from 3.4 nm (∼7 ns) to 2.9 nm (∼60 ns). After 60 ns Rg attains a constant value, although the molecules within the assembly get more organized as observed from the fluctuation of the order parameter S up to 80 ns. A complete 2D crystal is obtained at 90 ns with lattice parameters 6.51 ± 0.26 Å and 121.8° ± 4.6° as shown in Figure 1(f).

In order to identify the driving force behind the above processes, the benzene−benzene and benzene−BP interaction energies are examined and shown in Figure 2(b−d). It is clear that, during the stages of adsorption, collision (0−20 ps), and island formation (up to 1.5 ns), the interaction energies reduce, suggesting that during this time period the system rapidly goes into a more stable state. After that, the intermolecular electrostatic interaction energy of benzene molecules becomes nearly constant with some minor fluctuations, but the benzene−BP interaction energy goes up until it becomes constant near 20 ns, although the change in the latter one is sufficiently low. However, a steady and larger decrease in the time evolution of the intermolecular van der Waals interaction energy is observed, which follows the same trend as the radius 10216


Article


After that, the structure remains practically unperturbed. Hence, we estimate ∼60 ns as the time required to form crystalline TCNQ films on BP. For polar molecules, formation of the initial assembled islands through adsorption and collision is guided by both van der Waals and Coulomb interactions simultaneously. Figure 5(e) represents the average charge density profile along a particular axis (x axis) during the selfassembly of TCNQ at different time instants which further support the above conclusions. It is seen that for the first 10 ns the charge density profile has no particular pattern due to random arrangements of charges. The beginning of the phase transition into an ordered arrangement is realized when periodicity develops in the charge density at ∼30 ns. Further simulation leads to an increase in the periodicity of the charges leading to sharp and alternate peaks in the density profile between 50 and 100 ns (Figure 5e) which indicates increasing crystallinity of the TCNQ molecular assembly. Hence, conversion of the islands into a periodic framework is achieved by minimization of long-range Coulomb forces, which being much stronger compared to van der Waals interactions can easily propagate through the molecular framework. Our conclusions are in excellent agreement with the lateral charge screening mechanism proposed by Tsai et al. for the formation of nucleation islands for polar molecules (such as F4TCNQ) over graphene.20 Interestingly, unlike benzene, the TCNQ−BP vdW interaction decreases with some fluctuations up to 50 ns, before approaching constancy. This is presumably due to larger intermolecular interaction between TCNQ molecules, which make the initial adsorption over BP somewhat difficult. However, as seen from Figure 5(b) and 5(d), the magnitude of the TCNQ−BP interaction is significantly smaller than the intermolecular electrostatic attractions between TCNQ molecules; hence, we can safely rule out the importance of molecule−surface interactions for the 2D ordering of polar molecules. In the force-field control strategy, when the molecule− molecule van der Waals interaction was halved with charges being the same, the lattice remained perfectly intact. However, when the TCNQ−phosphorene vdW interaction was reduced by 50%, the nature of the assembly changed with simultaneous loss of the structural periodicity due to the same stacking effect as observed for nonpolar benzene. Both the assembly and 2D lattice of TCNQ were completely disrupted within 100 ps after withdrawing the partial charges from our simulation, confirming the prime importance of electrostatic interaction behind the 2D lattice formation (Figure S17). This conclusion is valid for other polar molecules over BP as well. Role of the Black Phosphorene Surface. Interestingly, for all molecules, the intermolecular van der Waals interaction energies were positive throughout our simulations, suggesting a sterically repulsive environment in the 2D lattice, irrespective of their polarity. For polar molecules, such an unfavorable arrangement can be stabilized by highly negative Coulomb interactions, but for nonpolar systems, faintly weak electrostatic forces cannot overwhelm steric repulsion. Recently, Tsai et al. showed through DFT calculations that the interaction between polar F4TCNQ molecules is repulsive, and 2D island formation is thermodynamically unfavorable in the absence of a graphene surface.18 Therefore, a very relevant question would be to enquire about the underlying force which lead to a 2D order in spite of such difficulties. Comparing the energy scales in Figures 2(d) and 5(d), it is evident that the molecule−phosphorene interactions are highly negative for both polar (TCNQ) and

of gyration (see Figure 2(a) and 2(c)). In fact, in our simulations with other nonpolar molecules (Figures S3 to S16) we observed a similar decrease in the molecule−molecule vdW interaction energies between adsorbed molecules. Pentacene and mesitylene show a small increase in the intermolecular electrostatic interaction energies as the assembly became more ordered. Therefore, for nonpolar molecules, the phase transition into a two-dimensional ordered state is predominantly governed by intermolecular vdW interactions rather than Coulomb forces. To get further insights, we adopted a force-field control protocol. A set of simulations were ran with the lattice structure of benzene molecules obtained at the end of the 200 ns production simulation with modified force field parameters. If all charges are eliminated, the original assembled structure (Figure 3(a)) remains almost unaltered (Figure 3(b)). On the other hand, when benzene−BP interaction was halved or benzene−benzene interaction was doubled, in either of the cases the benzene molecules became perpendicular to the surface within few nanoseconds, and the packing of the 2D lattice changes completely (see Figure 3(c)). Nevertheless, complete loss of the assembly took place within only 1 ns when benzene−benzene vdW interactions are halved (Figure 3(d)). Therefore, we conclude that, upon adsorption on the phosphorene surface, the herringbone packing of benzene is lost as the benzene−surface vdW interaction overweighs the stabilization gained through intermolecular interactions. However, if the benzene−BP interaction becomes less than the benzene−benzene interaction, the π stacking of benzene over phosphorene becomes more favored, and the planarity of the 2D lattice is lost. This explains why the plane of aromatic molecules such as pentacene and C8-BTBT becomes perpendicular in the presence of feebly interacting surfaces such as MoS2 and SiO2, leading to quasifreestanding vdW heterostructures.11,23 We obtained similar results for other nonpolar molecules (Figure S17) indicating that they follow the same general mechanism of self-assembly over BP. The mechanism of the planar molecules being “upright” is discussed in greater detail in the Supporting Information file (Figures S42 and S43), taking pentacene as an example. Predominantly Polar Molecules. We have studied the selfassembly followed by 2D lattice formation of 75 TCNQ molecules over BP at 350 K, and Figure 4(a−e) shows the snapshots of the simulation. All the events observed in the selfassembly of benzene (e.g., adsorption, collision, island formation, and their eventual coalescence) can also be traced for TCNQ. However, in this case, assembly formation is extremely fast as within only 5 ns Rg decreases from 3.8 to 3.35 nm (Figure 5a), and the order parameter (S) undergoes a sharp increase up to 0.9. After that, the observed changes are minor, and the molecules in the assembly organize themselves to achieve a thermodynamically more stable planar two-dimensional lattice, as observed by a small decrease and fluctuations in Rg between the time intervals of 5 and 60 ns. Figure 4(f) represents the unit cell of the TCNQ 2D lattice obtained from our simulations and the corresponding unit cell parameters. Figure 5(b−d) shows that for the first ∼7−8 ns there is a decrease in all the interaction energies and, in particular, the electrostatic component. It suggests the completion of the initial stages of adsorption and island formation. Subsequently, the van der Waals interaction energy among the TCNQ molecules remains nearly constant, although the electrostatic interaction energy undergoes a gradual decrease up to 60 ns. 10217


Article


Figure 6. (a) Snapshots representing initial and final structures used for the calculation of ΔGfree and ΔGassembly, respectively. The particular molecule, for which PMF calculations are performed, is shown in red. (b, c) PMF curves for benzene and TCNQ, respectively. Note that, unlike that for a free single molecule (red lines), the PMF curves representing detachment of a molecule from the assembly (black lines) have secondary minima at dz ≈ 8 Å, indicating stabilization due to π-stacking with adjacent molecules.

TCNQ, the overall direction of the molecular arrangement was random. For instance, among our six simulations for pentacene, we observed the long axis of the molecules to be oriented along the armchair direction in three cases and along the zigzag orientation in two of them. In the other simulation the long axis was oriented in between the above two. Thus, we can safely conclude that during self-assembly weakly adsorbed molecules easily surpass the small energy barriers between different adsorption geometries due to sufficient kinetic energies. In addition, free energies are calculated as ensemble averages through extensive sampling between all possible adsorption sites in phase space. Therefore, we believe that ΔGfree and ΔGassembly are not very sensitive to the site preferences as they are calibrated incorporating the probabilities of occurrences of different adsorption sites. To further test this assumption, we choose two different benzene molecules from the assembly and show their PMFs to be very similar (Figure S41). From the PMF profiles it is clear that it is more difficult to detach a molecule from an assembled state over phosphorene rather than a free molecule. Consequently, ΔΔGrel assembly is negative for rel all the molecules. A more negative value of ΔΔGassembly indicates higher preference for the molecule in the twodimensional ordered state. Table 1 shows that among the monocyclic aromatic molecules (1−4) ΔGassembly is relatively larger for the predominantly polar molecule trihydroxybenzene (4) than their nonpolar counterparts (1−3), although ΔGfree is similar for them. As a result, for polar molecules the self-assembled

nonpolar (benzene) molecules due to adsorption. In fact, we observed a steep decrease in molecule−phosphorene interaction energies during energy minimization, equilibration, and the first few nanoseconds of most of the production simulations. Evidently, adsorptions of the molecules act as the first trigger toward self-assembly. Nevertheless, it is not clear how surface affects the nature and pattern of 2D ordering once molecules are adsorbed. For a better understanding, the gain in stabilization energy for the incorporation of a free molecule within the assembly over the BP surface is quantified (see Figure 6(a) for schematic representation of free energy changes). We define the relative free energy preference for self-assembly for a molecule as ΔΔGrel assembly = ΔGassembly − ΔGfree. Table 1 presents all the free energy values for the eight molecules under study, and two representative PMF curves for benzene and TCNQ are shown in Figure 6(b) and 6(c) (see Figures S17 and S18 for other PMF curves). It is worthwhile to mention that 2D surfaces have multiple adsorption sites and hence have different local minima and exhibit different binding energies.45 Indeed, in our simulations performed at low temperatures ( phosphorene > graphene. This may be attributed to the difference in van der Waals interactions between the surface and the molecules, which follow the same trend (Figures S23 (a) and S27 (a)). However, the stabilization gained for a rel molecule through self-assembly (ΔΔGassembly ) displays the order h-BN (−7.6 kcal/mol) > phosphorene (−6.8 kcal/mol) > graphene (−6.2 kcal/mol) for nonpolar pentacene, while for polar TCNQ, exactly the opposite trend is found. Therefore, for nonpolar molecules, strongly adsorbing surfaces favor selfassembly more as intermolecular interactions are too weak to stabilize the lattice alone. Contrarily, for polar networks, molecular alignments through strong long-range electrostatics require larger penalty when the molecular motions are restricted due to strong adsorption on the surface. As a result weakly interacting surfaces are preferred for self-assembly. For instance, the 2D lattice of TCNQ is most stable over graphene even though ΔGfree is the least. Nevertheless, ΔΔGrel assembly for all the molecules are similar in graphene, h-BN and BP, suggesting that the 2D surfaces only contribute toward the overall stability of the assembly but cannot influence the fundamental nature of the molecular ordering. Though we have considered the formation of molecular assemblies over monolayers of BP, in experiments the situation is more complicated as the surfaces over which epitaxial growth is carried out are usually few layers in thickness. To better understand the multilayer effects, we have again chosen 10220


Article

The Journal of Physical Chemistry C pentacene as well as TCNQ and simulated their 2D ordering over a three-layer BP surface (Figures S29 and S31). It is seen that the formation (Tf) and melting (Tm) temperatures for the molecular assemblies remain the same, and 2D lattices are formed with similar lattice parameters as observed in the case of a monolayer. The time evolutions of the interaction energies between molecule and phosphorene as well as the intermolecular interaction energies for a trilayer and monolayer phosphorene surface are also identical (Figure S30 and S32). Additionally, Figure S33 (a) shows that ΔGassembly from a trilayer of phosphorene is only −1.0 and −0.3 kcal/mol higher compared to that from a monolayer for TCNQ and pentacene, respectively. Hence, multilayer surface effects are practically negligible, and they can neither influence molecular ordering nor the stability of the lattice.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +91-33-24734971. ORCID

Ayan Datta: 0000-0001-6723-087X Notes

The authors declare no competing financial interest.

■

■

ACKNOWLEDGMENTS AD thanks DST and BRNS for partial funding. We thank CRAY supercomputer and IBM P7 cluster for computational facilities.

CONCLUSION In summary, we have ventured into black phosphorene to evaluate its potential as a planar scaffold for the fabrication of two-dimensional molecular crystals with high periodicity. We selected a set of polar and nonpolar molecules with different shapes and sizes and successfully simulated their self-assembly and 2D lattice formation over black phosphorene substrate. It is shown that molecular ordering is always preceded by rapid adsorption on the surface in a favorable geometry. Once adsorption is complete, nonpolar molecules organize themselves through the minimization of intermolecular van der Waals interaction energies, while polar molecules assemble into a geometry where intermolecular electrostatic contacts are highly stabilizing. Our free energy estimations indicate that, for planar molecules adsorbed parallel to the surface, long-range Coulomb forces can stabilize the lattice more than weak van der Waals forces. Indeed, estimation of melting temperatures show that nonpolar crystals are thermally less stable compared to polar networks. However, the lattice of nonplanar and bulky fullerene was highly stable as it provides larger contact area to interact with neighbors. It seems that these stable 2D assemblies may be used for the passivation of BP. Also, it is confirmed that the tendency for adsorption has no correlation with the stability of the lattices which clearly demarks the active influence of the BP surface on self-assembly. Nevertheless, the surface governs the orientation of molecular ordering through interfacial interactions. Further comparison of self-assembly over graphene and h- BN establishes that the nonaromatic nature and the presence of buckling in BP do not affect the fundamental nature of molecular ordering, although overall stabilities of the lattices were different for the three substrates. In addition, multilayer effects were shown to be negligible. Therefore, the present paper clearly demonstrates that BP might be used for the epitaxial growth of organic thin films, and it simultaneously disentangles the various factors influencing the nature and stability of the 2D lattices grown over black phosphorene. We believe this work would serve as an a priori guideline for experimentalists to design specific self−assembled networks over black phosphorene for a variety of interesting applications involving heteroepitaxial thin films.

■

Unit cell structures, lattice parameters, time evolution of number density profile, analysis of self−assembly on BP, additional simulations at different coverages, force field modification simulations, PMF curves, multilayer effects and details of force−field parameters (PDF)

■

REFERENCES

(1) Rosei, F.; Schunack, M.; Naitoh, Y.; Jiang, P.; Gourdon, A.; Laegsgaard, E.; Stensgaard, I.; Joachim, C.; Besenbacher, F. Properties of Large Organic Molecules on Metal Surfaces. Prog. Surf. Sci. 2003, 71, 95−146. (2) Otero, R.; Gallego, J. M.; Vázquez de Parga, A. L.; Martín, N.; Miranda, R. Molecular Self − Assembly at Solid Surfaces. Adv. Mater. 2011, 23, 5148−5176. (3) Forrest, S. R. Ultrathin Organic Films Grown by Organic Molecular Beam Deposition and Related Techniques. Chem. Rev. 1997, 97, 1793−1896. (4) Yang, J.; Yan, D.; Jones, T. S. Molecular Template Growth and Its Applications in Organic Electronics and Optoelectronics. Chem. Rev. 2015, 115, 5570−5603. (5) Poelking, C.; Tietze, M.; Elschner, C.; Olthof, S.; Hertel, D.; Baumeier, B.; Wurthner, F.; Meerholz, K.; Leo, K.; Andrienko, D. Impact of Mesoscale Order on Open-Circuit Voltage in Organic Solar Cells. Nat. Mater. 2014, 14, 434−439. (6) Ecija, D.; Trelka, M.; Urban, C.; Mendoza, P. D.; Mateo-Martí, E.; Rogero, C.; Martín-Gago, J. A.; Echavarren, A. M.; Otero, R.; Gallego, J. M.; et al. Molecular Conformation, Organizational Chirality, and Iron Metalation of Meso-tetramesitylporphyrins on Copper (100). J. Phys. Chem. C 2008, 112, 8988−8994. (7) Cheng, Z. H.; Gao, L.; Deng, Z. T.; Jiang, N.; Liu, Q.; Shi, D. X.; Du, S. X.; Guo, H. M.; Gao, H. J. Adsorption Behavior of Iron Phthalocyanine on Au (111) Surface at Submonolayer Coverage. J. Phys. Chem. C 2007, 111, 9240−9244. (8) Écija, D.; Otero, R.; Sánchez, L.; Gallego, J.; Wang, Y.; Alcamí, M.; Martín, F.; Martín, N.; Miranda, R. Crossover Site-Selectivity in the Adsorption of the Fullerene Derivative PCBM on Au (111). Angew. Chem., Int. Ed. 2007, 46, 7874−7877. (9) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. Science 2004, 306, 666−669. (10) Novoselov, K. S.; Jiang, D.; Schedin, F.; Booth, T. J.; Khotkevich, V. V.; Morozov, S. V.; Geim, A. K. Two-dimensional Atomic Crystals. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 10451. (11) Huang, H.; Huang, Y.; Wang, S.; Zhu, M.; Xie, H.; Zhang, L.; Zheng, X.; Xie, Q.; Niu, D.; Gao, Y. van Der Waals Heterostructures Between Small Organic Molecules and Layered Substrates. Crystals 2016, 6, 113. (12) De Parga, A. V.; Calleja, F.; Borca, B.; Passeggi, M. C. G., Jr; Hinarejos, J. J.; Guinea, F.; Miranda, R. Periodically Rippled Graphene: Growth and Spatially Resolved Electronic Structure. Phys. Rev. Lett. 2008, 100, 056807. (13) Borca, B.; Barja, S.; Garnica, M.; Sánchez-Portal, D.; Silkin, V. M.; Chulkov, E. V.; Hermanns, C. F.; Hinarejos, J. J.; de Parga, A. V.;

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b02480. 10221


Article

The Journal of Physical Chemistry C Arnau, A.; Echenique, P. M.; Miranda, R. Miranda, R. Potential Energy Landscape for Hot Electrons in Periodically Nanostructured Graphene. Phys. Rev. Lett. 2010, 105, 036804. (14) Wang, Q. H.; Hersam, M. C. Room-Temperature MolecularResolution Characterization of Self-Assembled Organic Monolayers on Epitaxial Graphene. Nat. Chem. 2009, 1, 206−211. (15) Huang, H.; Chen, S.; Gao, X.; Chen, W.; Wee, A. T. S. Structural and Electronic Properties of PTCDA Thin Films on Epitaxial Graphene. ACS Nano 2009, 3, 3431−3436. (16) Pollard, A. J.; Perkins, E. W.; Smith, N. A.; Saywell, A.; Goretzki, G.; Phillips, A. G.; Argent, S. P.; Sachdev, H.; Muller, F.; Hufner, S.; et al. Supramolecular Assemblies Formed on an Epitaxial Graphene Superstructure. Angew. Chem., Int. Ed. 2010, 49, 1794−1799. (17) Chen, W.; Chen, S.; Qi, D. C.; Gao, X. Y.; Wee, A. T. S. Surface Transfer P-Type Doping of Epitaxial Graphene. J. Am. Chem. Soc. 2007, 129, 10418−10422. (18) Tsai, H.; Omrani, A. A.; Coh, S.; Oh, H.; Wickenberg, S.; Son, Y.; Wong, D.; Riss, A.; Jung, H. S.; Nguyen, G. D.; et al. Molecular Self-Assembly in a Poorly Screened Environment: F4TCNQ on Graphene/BN. ACS Nano 2015, 9, 12168−12173. (19) Lee, W. H.; Park, J.; Sim, S. H.; Lim, S.; Kim, K. S.; Hong, B. H.; Cho, K. Surface-Directed Molecular Assembly of Pentacene on Monolayer Graphene for High-Performance Organic Transistors. J. Am. Chem. Soc. 2011, 133, 4447−4454. (20) Xiao, K.; Deng, W.; Keum, J. K.; Yoon, M.; Vlassiouk, I. V.; Clark, K. W.; Li, A. P.; Kravchenko, I. I.; Gu, G.; Payzant, E. A.; et al. Surface-Induced Orientation Control of CuPc Molecules for the Epitaxial Growth of Highly Ordered Organic Crystals on Graphene. J. Am. Chem. Soc. 2013, 135, 3680−3687. (21) Cho, J.; Smerdon, J.; Gao, L.; Suzer, O.; Guest, J. R.; Guisinger, N. P. Structural and Electronic Decoupling of C60 from Epitaxial Graphene on SiC. Nano Lett. 2012, 12, 3018−3024. (22) MacLeod, J. M.; Rosei, F. Molecular Self-Assembly on Graphene. Small 2014, 10, 1038−1049. (23) He, D. W.; Zhang, Y. A.; Wu, Q. S.; Xu, R.; Nan, H. Y.; Liu, J. F.; Yao, J. J.; Wang, Z. L.; Yuan, S. J.; Li, Y.; et al. Two-Dimensional Quasi-Freestanding Molecular Crystals for High-Performance Organic Field-Effect Transistors. Nat. Commun. 2014, 5, 5162. (24) Yang, T.; Berber, S.; Liu, J. F.; Miller, G. P.; Tomanek, D. SelfAssembly of Long Chain Alkanes and Their Derivatives on Graphite. J. Chem. Phys. 2008, 128, 124709. (25) Zhang, Y.; Qiao, J.; Gao, S.; Hu, F.; He, D.; Wu, B.; Ji, W.; et al. Probing Carrier Transport and Structure-Property Relationship of Highly Ordered Organic Semiconductors at the Two-Dimensional Limit. Phys. Rev. Lett. 2016, 116, 016602. (26) Jariwala, D.; Howell, S. L.; Chen, K. S.; Kang, J.; Sangwan, V. K.; Filippone, S. A.; Hersam, M. C.; et al. Hybrid, Gate-Tunable, van der Waals p−n Heterojunctions from Pentacene and MoS2. Nano Lett. 2015, 16, 497−503. (27) Reich, E. S. Phosphorene Excites Materials Scientists. Nature 2014, 506, 19. (28) Crichton, W. A.; Mezouar, M.; Monaco, G.; Falconi, S. Phosphorus: New in situ Powder Data from Large Volume Apparatus. Powder Diffr. 2003, 18, 155−158. (29) Du, Y.; Ouyang, C.; Shi, S.; Lei, M. Ab initio Studies on Atomic and Electronic Structure of Black Phosphorus. J. Appl. Phys. 2010, 107, 093718. (30) Castellanos-Gomez, A.; Vicarelli, L.; Prada, E.; Island, J. O.; Narasimha-Acharya, K. L.; Blanter, S. I.; Groenendijk, D. J.; Buscema, M.; Steele, G. A.; Alvarez, J. V.; et al. Isolation and Characterization of Few-Layer Black Phosphorus. 2D Mater. 2014, 1, 025001−025020. (31) Yasaei, P.; Kumar, B.; Foroozan, T.; Wang, C.; Asadi, M.; Tuschel, D.; Indacochea, J. E.; Klie, R. F.; Salehi-khojin, A. HighQuality Black Phosphorus Atomic Layers by Liquid-Phase Exfoliation. Adv. Mater. 2015, 27, 1887−1892. (32) Li, L.; Yu, Y.; Ye, G. Y.; Ge, Q.; Ou, X.; Wu, H.; Feng, D.; Chen, X. H.; Zhang, Y. Black Phosphorus Field-effect Transistors. Nat. Nanotechnol. 2014, 9, 372−377.

(33) Koenig, S. P.; Doganov, R. A.; Schmidt, H.; Neto, A. H. C.; Ozyilmaz, B. Electric Field Effect in Ultrathin Black Phosphorus. Appl. Phys. Lett. 2014, 104, 103106. (34) Churchill, H. O. H.; Jarillo-Herrero, P. J. Two Dimensional Crystals: Phosphorus Joins The Family. Nat. Nanotechnol. 2014, 9, 330−331. (35) Liu, H.; Neal, A. T.; Zhu, Z.; Luo, Z.; Xu, X.; Tomanek, D.; Ye, P. D. Phosphorene: An Unexplored 2D Semiconductor with a High Hole Mobility. ACS Nano 2014, 8, 4033−4041. (36) Li, W.; Zhang, G.; Zhang, Y.-W. Electronic Properties of EdgeHydrogenated Phosphorene Nanoribbons: A First-Principles Study. J. Phys. Chem. C 2014, 118, 22368−22372. (37) Wittig, J.; Matthias, B. T. Superconducting phosphorus. Science 1968, 160, 994. (38) Okajima, M.; Endo, S.; Akahama, Y.; Narita, S.-I. Electrical Investigation of Phase Transition in Black Phosphorus Under High Pressure. Solid State Commun. 1984, 49, 879. (39) Rodin, A. S.; Carvalho, A.; Neto, A. H. C. Strain-Induced Gap Modification in Black Phosphorus. Phys. Rev. Lett. 2014, 112, 176801. (40) Willke, P.; Amani, J. A.; Sinterhauf, A.; Thakur, S.; Kotzott, T.; Druga, T.; Weikert, S.; Maiti, K.; Hofsäss, H.; Wenderoth, M. Doping of Graphene by Low-Energy Ion Beam Implantation: Structural, Electronic, and Transport Properties. Nano Lett. 2015, 15, 5110− 5115. (41) Liu, H.; Liu, Y.; Zhu, D. Chemical Doping of Graphene. J. Mater. Chem. 2011, 21, 3335−3345. (42) Tang, Q.; Zhou, Z.; Chen, Z. Molecular Charge Transfer: A Simple and Effective Route to Engineer the Band Structures of BN Nanosheets and Nanoribbons. J. Phys. Chem. C 2011, 115, 18531− 18537. (43) Zhang, R.; Li, B.; Yang, J. A First-Principles Study on Electron Donor and Acceptor Molecules Adsorbed on Phosphorene. J. Phys. Chem. C 2015, 119, 2871. (44) Jing, Y.; Tang, Q.; He, P.; Zhou, Z.; Shen, P. Small Molecules Make Big Differences: Molecular Doping Effects on Electric and Optical Properties of Phosphorene. Nanotechnology 2015, 26, 095201−095209. (45) Chowdhury, C.; Jahiruddin, S.; Datta, A. Pseudo-Jahn−Teller Distortion in Two-Dimensional Phosphorus: Origin of Black and Blue Phases of Phosphorene and Band Gap Modulation by Molecular Charge Transfer. J. Phys. Chem. Lett. 2016, 7, 1288−1297. (46) Feng, J.; Ding, H.; Ma, Y. Self − Assembly of Fullerenes and Graphene Flake: A Molecular Dynamics Study. Carbon 2015, 90, 34− 43. (47) Ozmaian, M.; Fathizadeh, A.; Jalalvand, M.; Ejtehadi, M. R.; Allaei, S. M. V. Diffusion and Self-Assembly of C60 Molecules on Monolayer Graphyne Sheets. Sci. Rep. 2016, 6, 21910. (48) Zhao, Y.; Wu, Q.; Chen, Q.; Wang, J. Molecular Self-Assembly on Two-Dimensional Atomic Crystals: Insights from Molecular Dynamics Simulations. J. Phys. Chem. Lett. 2015, 6, 4518−4524. (49) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant-Pressure Molecular-Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177− 4189. (50) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - an N· Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (51) Andersen, H. C. Rattle - A Velocity Version of the Shake Algorithm for Molecular-Dynamics Calculations. J. Comput. Phys. 1983, 52, 24−34. (52) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157− 2164. (53) Kale, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. NAMD2: Greater Scalability for Parallel Molecular Dynamics. J. Comput. Phys. 1999, 151, 283−312. (54) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. 10222


Article

The Journal of Physical Chemistry C (55) Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive Biasing Force Method for Scalar and Vector Free Energy Calculations. J. Chem. Phys. 2008, 128, 144120. (56) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Field. J. Comput. Chem. 2010, 31, 671−690. (57) Sresht, V.; Padua, A. A. H.; Blankschtein, D. Liquid Phase Exfoliation of Phosphorene: Design Rules From Molecular Dynamics Simulations. ACS Nano 2015, 9, 8255−8268. (58) Reddy, C. D.; Yu, Z. G.; Zhang, Y. Two-Dimensional van der Waals C60 Molecular Crystal. Sci. Rep. 2015, 5, 12221. (59) Mukhopadhyay, T. K.; Datta, A. Deciphering the Role of Solvents in the Liquid Phase Exfoliation of Hexagonal Boron Nitride: A Molecular Dynamics Simulation Study. J. Phys. Chem. C 2017, 121, 811−822.

10223


Ordering and Dynamics for the Formation of Two-Dimensional

Recommend Documents