Force Field Parametrization of Colloidal CdSe Nanocrystals Using an

Subscriber access provided by TRAKYA UNIVERSITESI

Article

Force Field Parametrization of Colloidal CdSe Nanocrystals Salvatore M. Cosseddu, and Ivan Infante J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01089 • Publication Date (Web): 18 Nov 2016 Downloaded from http://pubs.acs.org on November 22, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Force Field Parametrization of Colloidal CdSe Nanocrystals Salvatore Cosseddu# and Ivan Infante# Department of Theoretical Chemistry and Amsterdam Center for Multiscale Modeling (ACMM), VU University Amsterdam, De Boelelaan 1083, 1081 HV Amsterdam, The Netherlands

KEYWORDS. Force Field Parametrization, canonically-averaged force fitting, Monte Carlo, Simulated Annealing, Adaptive Rate Monte Carlo, Metropolis Monte Carlo, Born Oppenheimer Molecular Dynamics, Density Functional Theory, CdSe Nanocrystals

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 40

ABSTRACT

In a typical colloidal CdSe nanocrystal more than 50% of the atoms are located at the surface. These atoms can give rise to electronic traps that can deteriorate the performance of optoelectronic devices made of these nanomaterials. A key challenge in this field is thus to understand with atomistic detail the chemical processes occurring at the nanocrystal surface. Molecular dynamics simulations represent an important tool to unveil these processes, but its implementation is strongly limited by the difficulties of finely tuning classical force fields parameters, primarily caused by the unavailability of experimental data of these materials that are suitable in the parametrization procedures. In this work, we present a general scheme to produce force field parameters from first principles calculations. This approach is based on a newly-developed stochastic optimization algorithm called Adaptive Rate Monte Carlo, which is designed to be robust, accurate, easy-to-use and flexible enough to be straightforwardly extended to other nanomaterials. We demonstrate that our algorithm provides a set of parameters capable of satisfactorily describing non-stoichiometric CdSe nanocrystals passivated with oleate ligands akin to experimental conditions. We also demonstrate that our new parameters are robust enough to be transferable among crystal structures and nanocrystals of increasing sizes up to the bulk.


2

Page 3 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction. Colloidal semiconductor nanocrystals (NC) show unique optoelectronic properties deriving from large absorption cross-sections and bandgap tunability with size1,2 that allow them to find applications in a large variety of fields, such as solar cells,3 light emitting devices (LEDs),4 photodetectors, photocatalysis5 and bio-sensors.6 Despite the large amount of available experimental data and applications, research that seeks to improve NC engineering by understanding their underlying atomistic processes remains a big challenge. In this context, theoretical modeling is assuming a fundamental role in the interpretation of experimental data and in guiding the design of new materials by capturing nontrivial aspects of the NC unique chemistry, like surface processes,7 ligand passivation8–10 and charge carrier trapping.11–13 Binary II-VI and IV-VI NCs such as CdE and PbE (E= S, Se, Te), usually referred to as quantum dots (QD), are prototypical for the investigation of the properties of these nanomaterials.14 They exhibit non-stoichiometric crystalline structures with an excess of metal cations at the surface, charge compensated by a shell of surfactant organic anions employed in the colloidal NC synthesis.15–19 Theoretical investigations of these hybrid inorganic/organic systems require a dynamical modeling of the surface under realistic conditions (i.e. in the presence of solvent and surfactant molecules) to capture critical processes like surface reconstructions and ligand mobility. Fully classical molecular dynamics (MD) are, at present, the best option to explore the large time scales and phase space regions over which these events occur. In fact, the extremely high computational costs of ab-initio MD simulations strongly limit their applicability and, in the context of colloidal seminconductor nanocrystals, they have been used in the past only to simulate small NCs in vacuum, passivated by very small ligands, for very short timescales


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 40

(picoseconds).11,20–22 Fully classical MD, on the other hand, requires a targeted optimization of the potential energy functions (force fields, FF) to accurately describe the significant deviation of NC properties from the bulk.

Figure 1. Several CdSe NCs snapshots taken at the last point of MD trajectories. (top) CdSe NCs of 1.8 nm in diameter, sampled during 15 ps; (bottom) CdSe NCs of 2.3 nm in diameter, sampled during 20 ps. These simulations were carried out using: (a) ab-initio MD, (b) classical MD using the force field developed in this work, and (c) classical MD using the force field fitted by Rabani on bulk properties. In both classical MD simulations, the FFs were mixed with parameters from CHARMM FFs to describe the formate ligands. In case of (c), Oxygen charges were fixed to -0.606e to preserve charge neutrality of the nanocrystal and to mimic the charge transfer from the ligands, as seen in the DFT calculations.


4

Page 5 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Parameters for CdSe materials developed by Rabani for a similar Lennard-Jones FF have been used in the past to simulate CdSe NCs,23,24 despite their optimization did not focus on reproducing NCs properties but on fitting some experimental properties of their bulk counterparts (lattice and elastic constants, phonon dispersion relation, etc.). A similar approach was followed by Yalcin et el.25 that have produced better performing and transferable parameters for CdSe using a Buckingham force field. In these cases, the target NCs were simulated either in the absence of surfactant ligands or by covalently binding ligands on the surface with an arbitrary large force constant. In a semiconductor colloidal NCs, however, a great percentage of the atoms is passivated with organic ligands, strongly bound to the surface, that serve to stabilize and functionalize the NC. Consequently, explicit NC-ligand interactions need to be explicitly parameterized. In this contest, simple FFs based on combinations of Lennard-Jones (LJ) potentials and point charges are extremely attractive because of the large spatial and time scales they can reach, and because of the availability in the literature of a number of robust parameters to model the most common passivating organic ligands.26 To prove why a tailored FF parametrization of colloidal CdSe NCs is important, we compared the outcome of an ab-initio (DFT) MD simulations with a fully classical MD simulation constructed from Rabani’s23,24 FF parameters for the inorganic CdSe core and CHARMM FF parameters for the ligands26 (See also Section S1 in SI for further methodological details). We have chosen for this test a 1.8 nm CdSe QD in vacuum and passivated by formates (see below for details). This calculation reveals the inability of Rabani’s FF parameters to maintain the NC crystalline structure (Figure 1, top-right), which, on the other hand, is finely described by the abinitio MD model (Figure 1, top-left). Although for larger NCs, the inclusion of the solvent and longer ligands such as oleates somehow alleviate this issue, these calculations emphasize how


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 40

targeting the NC surface and the ligand-NC interactions is crucial for the reliability of the FF. This is a great challenge because experimental data on the NC surface are scarce and not sufficiently detailed to be employed in the parameterization. For this work we thus developed an automated FF parameterization tool to retrieve FF parameters directly from ab-inito MD calculations (Born Oppenheimer Molecular Dynamics). This was used to generate a force field for CdSe specifically targeted on colloidal NCs. Central to this approach is a newly-developed stochastic optimization algorithm that we called Adaptive Rate Monte Carlo (ARMC). This method is based on an on-the-fly modification of the explored hypersurface borrowed from the Wang-Landau algorithm27 and is designed to be robust, accurate, easy-to-use, and governed by a small number of problem-independent algorithm parameters. In this manuscript, we firstly describe the overall parametrization scheme, then we present the ARMC algorithm, and finally we discuss the application on CdSe NCs by simulating a realistic sized NCs passivated with long alkylic chains in the presence explicit solvent molecules. We then compare the performance of the optimized force field with experimental data and available force fields from the literature.

2. Methodology and Results. Our FF parametrization scheme is based on the idea that a set of physical properties of interest computed using the fully classical model (molecular mechanics, MM) must converge to reference values calculated using higher-level quantum chemical methodologies (quantum mechanics, QM). In this work, the properties of interest are the radial distribution functions (RDFs) of model CdSe NCs computed using ab-initio MD simulations based on Density Functional Theory (DFT). Our parameterization scheme follows a three-step procedure (Scheme


6

Page 7 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1): (i) the first step consists of preparing a model for a reference (small) CdSe NC and computing, with ab-initio MD, the relevant radial distribution functions (RDFs); (ii) in the second step, a FF potential is chosen to run MM simulations and the ARMC algorithm (vide infra) is employed to optimize the FF parameters so that MM RDFs closely match the QM reference; (iii) in the final step, the transferability and scalability of the newly generated FF are tested on CdSe NCs of different sizes, on longer simulation time scales, and against experimental bulk properties of different crystal structures.

Scheme 1. Flow chart for the FF parameterization procedure using the ARMC scheme presented in this work. A more detailed flow chart is shown in the SI (Scheme S1).

2.1 Step 1: Preparation of the NC models. A non-stoichiometric Cd68Se55 NC with a diameter of 1.8 nm and a zinc-blende structure (Figure 1, top) was chosen as reference model . It exhibits a Cd:Se ratio of 1.1, in accordance to experimental evidence,17 and the excess of metal charge was compensated by adding 26 formate molecules. Formate was chosen to emulate the native oleate groups, commonly employed during the synthesis of these NCs.28


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 40

Figure 2. (a) Comparison between RDFs computed at different level of theory for the 1.8-nm CdSe NCs. (b) Comparison between the free energies computed with Equation (2) from the RDFs. (Blue lines) RDFs (or free energies) computed using classical FF potential using our optimized FF parameters;


8

Page 9 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(magenta line) RDFs (or free energies) computed with ab-initio MD. (Yellow line) RDFs (or free energies) computed with FF parameters taken from literature.23,24

The ab-initio MD simulations (QM), were carried out with Density Functional Theory (DFT)29 with a PBE exchange-correlation functional30 and a double-zeta basis set with polarization functions as implemented in the CP2k 2.5.1 quantum chemical package.31,32 This combination has been shown to be reliable in describing the structural properties of CdSe quantum dots.33 The integration step was chosen to be 1 fs and the velocity-rescaling thermostat was used to maintain a constant temperature of 300 K.34 Canonical simulations were run for 15 ps, with the initial 5 ps discarded as equilibration step. RDFs were computed for each pair of atoms of the inorganic CdSe core and between the inorganic core atoms and the anchoring oxygen atoms of the formate ligands bound to the NC surface. The relevant RDFs are shown in Figure 2 with magenta lines. Doubling the length of the simulation did not modify the RDFs.

2.3 Step 2: FF parametrization using the Adaptive Rate Monte Carlo (ARMC) scheme 2.3.1. Energy functions and parametrization approach In our models, NC and ligand-NC interatomic interactions are described by the effective pairwise-addictive potential: 𝑉 𝑟#$ =

4𝜀#$ #,$

𝜎#$ 𝑟#$

)*

𝜎#$ − 𝑟#$

,

+ #,$

𝑞# 𝑞$ 1 𝜖0 𝑟#$

where the first term is the Lennard-Jones 12-6 potential and includes both the Van der Waals interactions and the Pauli repulsions defined by the well-depth eij and distance sij; the second term is the Coulomb electrostatic interaction between the qi and qj point charges. This potential is one of the most popular choice to represent non-bonded interactions in MD simulations


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 40

because it provides a good compromise between performance and accuracy,35 allowing the simulation of realistic size systems with highly-parallelized MD software (i.g. GROMACS36 and NAMD37). Using only non-bonded interactions to describe our system means that all pairwise interactions inside the NC, and between the ligands and the NC, are considered as being primarily electrostatic (ionic), with the covalent components having a minor contribution adequately approximated within the non-bonded terms. Although this approximation may represent a limitation, it otherwise allows diffusion of ligands on the surface and better reproduces the dynamical nature of the NC surface, in agreement with experimental data. In this framework, our parameterization scheme focused only on the parametrization of the NC inorganic core and the ligand-NC interactions, while the ligand and solvent parameters were taken from the literature, including bonded terms described as in the CHARMM FF.26 The total number of parameters optimized in this work boiled down to the following set: (i) the LJ parameters (ε and σ) for the Cd-Cd, Se-Se, Cd-Se, Cd-O, and Se-O pairs; (ii) the point charges of Cd, Se in CdSe and of O and C of the carboxylate groups in the ligands. Overall, the total number of parameters to be optimized was 12. In addition, we also considered that CdSe NCs are non-stoichiometric and present an excess of positive charge stemming from the presence of a larger number of Cd atoms than Se during their synthesis. Usually this excess of charge is compensated by the carboxylate anion ligands anchoring the surface. This feature was taken into account in our optimization scheme by imposing the following constraints: q(Se) = −q(Cd) and q(C) = 2q(O) + 0.5q(Cd). The development of a robust and automatized approach for the construction of FF parameters primarily depends on the choice of an error function ξ that is a measure of the separation between the QM reference and the tested MM model. The FF development becomes therefore an


10

Page 11 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


optimization problem of ξ in the FF-parameter space. In this context, a major issue is notoriously posed by the point charges, which have a long-range nature difficult to capture.38 As anticipated, we chose to construct our ξ on the radial distribution functions (RDFs), i.e. the distribution of probabilities to find a pair of atoms i and j at a distance rij. This provides a description of nonbonded interactions over long ranges and allows an accurate collective optimization of LJ parameters and charges. In addition, our approach allows a simultaneous fitting of the free energy term providing, unlike alternative force matching procedures, an immediate and unambiguous evaluation of the accuracy of the FF in delivering accurate structural features. The choice of RDFs can be easily understood from standard statistical mechanics. At equilibrium, a molecular system is governed by its free energy, F, which can be written in terms of interatomic distances rij as: 𝐹 𝑟#$ = −𝑘5 𝑇𝑙𝑛 𝑔 𝑟#$

+𝐶 2

where g(rij) is the RDF, kB the Boltzmann constant, T the temperature and C a constant.39 By definition, equation (2) corresponds to the potential of mean forces (PMF), i.e. the interactions between the two atoms i and j held at distance rij, with all the remaining interactions canonically averaged. From (2) it is therefore possible to construct an error function between the reference ab-initio MD and the fully classical MD that allows capturing the correct free energy of the system. Both QM and MM free energies have an error ξ from the true free energy value: Ftrue = FMM + ξMM = FQM + ξQM. A relative error ξQM-MM of the MM from the QM reference can then be conveniently defined as the free-energy difference: ∆𝐹 𝑟#$ = 𝐹== 𝑟#$ − 𝐹>= 𝑟#$ = 𝜉>= − 𝜉== = 𝜉>=@== 3


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 40

Equation (3) is the mathematical form for the obvious necessity of accurate QM references, i.e. the quality of the optimized FF is primarily determined by the accuracy of the ab-initio MD simulations. A similarly defined free-energy difference has been used in the literature to approximate an energy error, iteratively applied to correct the interatomic potential in the Iterative Boltzmann Inversion method proposed by Reith et al.40,41 In our approach, ξQM-MM is used to construct a collective error function that is minimized. This is done by rewriting (2) and (3) as a function of ∆g(rij) = gMM(rij) − gQM(rij), the distances between MM and QM RDFs, as follow: 𝜉>=@== = −𝑘5 𝑇𝑙𝑛 1 +

∆𝑔 𝑟#$

4

𝑔>= 𝑟#$

As expected, the error ξQM-MM goes to zero as ∆g(rij) → 0. This allows us to write an auxiliary error function ℇ for collectively minimizing the ∆g(rij) for n atom pairs ij of interest, over arbitrary ranges (0, rmax) of rij: DEFG

ℇ>=@== =

∆𝑔C 𝑟#$ C

*

5

DHI JK

where ℇ>=@== is function of the FF parameters. Its form as a sum of Euclidean distances of QM-MM RDFs is particularly suitable for the optimization procedure. Figure 3 shows how the presented approach is stable in optimizing charges. The error function ℇ>=@== was tested on different sets of pairs of atoms to compute the RDFs (Figure 3a): P1 includes only pairs of atoms stemming from the NC, i.e. Cd–Cd, Se–Se, and Cd–Se; ii) P2 additionally includes atom pairs between NC atoms and the anchored oxygen atoms, i.e. Cd–O and Se–O; iii) P3 considers also the O–O RDF. Figure 3 depicts the complete surfaces of ℇ>=@== computed firstly by varying q(Cd) in the range (0,2) (Figure 3b) and, secondly, q(O) in (−1,0) (Figure 3c), while the remaining parameters were maintained to almost optimized values.


12

Our computed surfaces also demonstrate that a good estimation of the charges is obtained irrespective of the P size. Because RDFs used to build ℇ>=

== are

correlated, a larger number

does not necessarily deliver additional information. On the other hand, larger P sets provide better separation between the noise and the true underlying surface that can be greatly beneficial in many practical cases. Further details on how to noise was filtered in the productive runs are provided in the SI (Section S4 and Figure S1).

Cd Se O Cd Cd,Cd Cd,Se Cd,O Se Se,Se Se,O O O,O

a)

P1

b) εQM-MM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


P2

P3

600

400

P1 P2 P3

200

0 0.0

0.5

1.0

1.5

2.0

q(Cd) (e)

c) εQM-MM

Page 13 of 40

250

200

150

P2 P3

100 −1.00

−0.75

−0.50

−0.25

0.00

q(O) (e) Figure 3. Panel a) reports different P sets of pairs of atoms type used to build the error function


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 40

ℇ>=@== . P1 includes all pairs of NC atoms; P2 additionally comprises pairs of NC atoms and the attached oxygen atoms (i.e. all atom pairs involved in the parameterized LJ interactions); and P3 further includes the RDFs between oxygen atoms. Panel b) and c) reports ℇ>=

== surfaces

computed by

varying q(Cd) in the range (0,2) and q(O) in (−1,0).

2.3.2 Adaptive Rate Monte Carlo To perform global optimizations of ℇ>=

== ,

we developed an algorithm driven by the will to

create a robust algorithm that is able to provide high-quality solutions, that requires a minimal interference from the final user and that is easy to implement. Monte Carlo methods are iterative algorithms in which, at each iteration I, a new set of parameters SI (state) is randomly generated (move), tested and accepted according to some probability. Differently from many conventional approaches, our Adaptive Rate Monte Carlo (ARMC) is not built upon the Metropolis-Hastings rule,42,43 but on a simpler binary probability based on the optimization problem. In our case: 𝑝 𝑖←𝑗 =

0𝑖𝑓∆ℇ>=@== > 0 6 1𝑖𝑓∆ℇ>=@== < 0

where ∆ℇ>=@== = ℇ# −ℇ$ , i.e. the difference of ℇ>=@== at subsequent iterations j and i of the optimization procedure. Equation (6) corresponds to Metropolis-Hastings probability when T = 0. A simple application of (6) rapidly descends to the closest local minimum. Global exploration is obtained by incrementing ℇ$ in a systematic way by a factor 𝜙 > 0 at each iteration in which the state Sj is visited (ℇ$,WX) = ℇ$,W + 𝜙), in a fashion ideally analogous to successful biased methods for sampling Boltzmann distributions, most notably the Wang-Landau algorithm.27 In this way, the underlying hypersurface is continuously modified, retaining memory of the visited


14

Page 15 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


states. A key feature of the method is that 𝜙 is adapted each ω iterations in order to maintain the average acceptance rate 𝑎 near a target rate at: 𝜙[\ = 𝜙 [@) \ ∙ 𝛾 _`C

ab @a(def)

𝑘 = 1,2,3, … , 𝑁 7

where Nω = M is the total number of iterations, γ > 1 regulates 𝜙 correction, and 𝑎([@)) is computed over all the ω iterations of the (k−1)th block, and can assume only values between 0 and 1. ARMC generates a non-markovian random walk driven by the on-the-fly modification of the underlying hypersurface, that is capable of promoting the exploration of new space regions. This means that ARMC is theoretically able to sample every state S as 𝑀 → ∞, given any 𝑎n > 0. This is clearly not necessary, and the actual value of at becomes important for practical purposes: for example, at=1 and at=0 correspond respectively to diffusive behavior and freezing in local minimum, similar to high and low T in Metropolis Monte Carlo (MMC) methods. All the values of at in between 0 and 1 might be seen as if, in MMC, T was dynamically adjusted to optimize the local exploration. The actual value of at determines the level of this optimization, providing an intuitive parameter for tuning the balance between accuracy and extent of the search. Additional important advantages of defining at in ARMC algorithm is that the search never stops, being designed to overcome the “freezing” issues of similar stochastic global optimization algorithms, such as MMC combined with Simulated Annealing. It is sometime convenient to automatically control at by applying an annealing-like procedure instead of setting a fixed value for the entire calculation. We label this as “simulated braking” (ARMC-SB). ARMC-SB should theoretically have advantages on plain ARMC in many cases, for instance when surface roughness and the distance between starting point and solution are


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 40

completely unknown (for instance when a random initial state is used). More detailed information on this approach is provided in the SI (Section S2).

2.3.3. Set-up of the ARMC Procedure At each step of the optimization procedure using the ARMC method, RDFs were computed by generating classical MD trajectories for the Cd68Se55 in a canonical ensemble using the same setup employed for the ab-initio simulations, i.e. 15 ps of dynamics, with the initial 5ps discarded as equilibration. Also in this case the integration step was chosen to be 1 fs and the velocityrescaling thermostat was used to maintain a constant temperature of 300 K. More details about these MM simulations, performed with GROMACS 5.0.4,44 are reported in section S1 of the SI. The MM RDFs were then compared to the reference ab-initio RDFs to compute ℇ>=

==

and

the ARMC algorithm was used to find the most optimal FF parameters that minimized ℇ>=

== .

In computing the RDFs, a rmax = 12 Å was used, that is the length of the LJ and real-space Coulomb cut-offs, together with a bin width of 0.05 Å. A careful choice of the latter (bin width) is important because it influences the numerical noise. After some testing, in the ARMC procedure, we set the ω intervals to 100, ϕ = 1 and γ = 2 for the productive runs. See section S3 in SI for further details on how we chose these parameters.

2.3.4. Performance of ARMC Exhaustive testing of ARMC and comparisons with similar methods are out of the scope of this manuscript, and will be topic of future publications. In this section, however, we show the performance of ARMC algorithm in finding the FF parameters using different settings: i) ARMC with different acceptance rates, and ii) 2 ARMC-SB. For sake of clarity, among the many


16

Page 17 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculations performed, we report only the results obtained with 60000 iterations with the P2 set and without treatment of the noise. These can be considered non-optimal starting conditions and were chosen to validate the algorithm for more difficult and general cases. Final solutions were assessed with 10 independent MD runs to filter out false solutions (outliers) due to the lack of noise filtering. The best minimum is simply labeled as ℇo#C , while the mean and the associated standard deviation over the 10 MD runs are labeled as ℇo#Cfp and s, respectively. A comparison is provided against Metropolis Monte Carlo Simulated Annealing (MMC-SA), one of the most popular stochastic optimization algorithm used in similar FF optimizations.42,43 Figure 4 depicts the different behaviors of ARMC, ARMC-SB and MMC-SA for the runs that resulted in the respective best minima. A summary of the final results is further presented in Table 1. All the best ℇo#C minima are provided by ARMC and ARMC-SB that lied between 96.1 and 92.6. Plain ARMC was found as being the fastest method to descend towards the low region in the FF-parameter space, requiring a small number of steps to find good solutions. In cases different from the one under investigation this might become important, for instance, when dealing with larger systems or when computational resources are limited. ARMC-SB usually provided the best solutions and delivered a better final local refinement. MMC-SA performed worse than ARMC on our problem and provided the largest solution variability among different runs. The best minimum found in the various test runs was 115.57, ~20% higher than ℇo#C of ARMC-SB, and was found scaling kBT from 10000 every 100 steps by 0.948. Results variability arose from “freezing” (being trapped) in local minima, which is common to many stochastic algorithms. Notably, setting of MMC-SA required more effort than ARMC-SB to optimize the cooling scheme for the given problem. FF parameters optimized with MMC-SA were unable to maintain the nanocrystalline structure of the quantum dot.


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. ARMC and ARMC-SB are compared by optimizing the error function ℇ>=

Page 18 of 40

== for

the P2

pair set without treatment of the noise. A maximum of 60000 iterations was carried out. Only the runs that provided to the best solution (minimum) are plotted for plain ARMC (at= 0.25), ARMC-SB (at : 1 to 0.002) and MMC-SA (kT from 10000 scaled by 0.948 every 100 steps).

Figure 5 shows the effect of the average acceptance rate, a key parameter in ARMC. The at rates are correctly maintained by both ARMC and ARMC-SB algorithm. In principle, higher at rates should guarantee broader sampling while lower values should theoretically offer a greater accuracy of the local search. In practice, for high at rates we observed an increase in the sampled space but not a dramatic change in accuracy. This is more likely due to the noisy nature of the optimized function. In Table 1, we also present results for different at values, which shows that differences are very small and within the noise uncertainty. However, when using plain ARMC, we usually found the best minima for at between 0.2 and 0.3, which appears a good compromise between sampling and accuracy.


18

Page 19 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Acceptance rate of different ARMC runs: i) at = 0.25, ii) at = 0.55 and iii) simulated braking (ARMC-SB), at from 1 to 0.002.

ARMC has additional important advantages. First of all, the algorithm is never trapped in a local minimum and continues the search until the very last iteration. As the number of iteration increases, it clearly becomes more difficult to improve the solution, however the “adapted” modifications 𝜙 prevent the freezing and widen the explored space. A second major advantage over similar methods is that ARMC and ARMC-SB can be easily suited to refinement of preoptimized FF-parameters, i.e. starting from a region already close to the best solution, by lowering at. Method

Run

ℇo#C

ℇo#Cfp

s

ARMC

at = 0.55

96.1

102.2

6.0

ARMC

at = 0.25

95.7

101.6

4.2

ARMC

at = 0.10

96.5

102.1

2.9

ARMC-SB

at = from 1 to 0.002

92.6

97.4

3.7

MMC-SA

kBT = from 10000

107.0

115.57

5.36


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 40

scaled by 0.948 every 100 steps Table 1. Comparison of the best minimum values of ℇo#C found with different methods in the same set of optimization runs (60000 iterations, P2 atom-pairs set). The set of FF parameters correspondent to these minima were tested with 10 trials for which ℇo#Cfp and standard deviation (σ) are also reported.

2.4 Step 3: Test of the FF Parameters on CdSe Nanocrystals The final optimized FF parameters for the fully passivated CdSe NC are reported in Table 2. These were obtained using ARMC-SB in two steps: (i) a global pre-optimization and (ii) a final refinement in which noise filtering was applied (see Section S4 and Scheme S1 in SI for more details). Atom-pairs set P3 was used to compute ℇ>=

== .

Convergence was overall good in

general and the largest uncertainties were observed for the e parameters because their variations affect more weakly the overall potential (Eq. (1)) than the oscillations of other FF parameters. Lennard-Jones Charges (e)

s(nm)

e(kJ/mol)

Cd

0.9768

Cd-Cd

0.1234

0.3101

Se

-0.9768

Se-Se

0.4852

0.4266

O

-0.4704

Cd-Se

0.2940

1.5225

C

0.4524

Cd-O

0.2471

1.8340

Se-O

0.3526

1.6135

Table 2. Force Field parameters optimized to describe crystalline and formate-NC interactions in fullypassivated CdSe NCs. ARMC-SB was employed during the optimization procedure using the P3 atompairs set to compute ℇ>=

== .

Full details about the procedure are reported in Section S4 of SI.


20

Page 21 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Unlike Rabani’s FF parameters that were fitted to bulk properties, our FF maintained the crystalline structure of the NC as in the QM calculations (Figure 1 and Video V1 in SI). Figure 2a shows the agreement between the RDFs computed with DFT and MM simulations. This agreement turns out to be exceptionally good considering the strong differences between the two computational models, while the computational speed of the fully classical MD was about 7 orders of magnitude faster than the ab-initio MD. The 15 ps trajectory was simulated with the QM-MD in 6 days on 48 cores, while only in 6 seconds with the MM-MD on 12 cores (Intel Xeon E5-2650 v2, 2.60 GHz). We also notice that the majority of the RDFs peaks, and corresponding minima when passing to free energies (Figure 2b), are well reproduced. The most noticeable deviation is in the Cd–Cd RDF, which is shifted to higher values (to the right) in the MM RDF. We explain this effect by observing that the Cd atoms close to the surface at Se-rich (111) facets exhibited outward fluctuations more marked in the MM model. This was associated to a slightly enhanced mobility of Se atoms, which is noticeable in the free energies plot (Figure 2b). These deviations are however small and it is therefore difficult to clearly estimate their real physical meaning. We also found a very good scalability of the FF with the nanocrystal size. Larger nonstoichiometric CdSe NC were prepared with similar Cd:Se ratios of ~1.1. In the first test, the scalability was tested on a 2.3-nm Cd152Se135 NC (Figure 1, bottom) simulated both at DFT and MM levels for 20 ps, with the initial 10 ps discarded as relaxation, in vacuum and in the presence of 34 passivating formates. In the second test, a Cd1012Se911 NC of realistic size, 4-nm, was simulated in the presence of 202 oleates, experimental ligands, both in vacuum and in dichloromethane at constant temperature (300 K) and pressure (1 atm).45 The preparation of the


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 40

latter model is rather involved and a step-by-step explanation is provided in the Appendix for researchers that are approaching this field.


22

Page 23 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. (a) Comparison between RDFs computed at different level of theory for the 2.3-nm CdSe NCs. (b) Comparsion between the free energies computed with Equation (2) from the RDFs. (blue lines) RDFs (or free energies) computed using classical FF potential using our optimized FF parameters; (magenta line) RDFs (or free energies) computed with ab-initio MD. (yellow line) RDFs (or free energies) computed with FF parameters taken from literature.23,24

In the case of the 2.3-nm NC, both the RDFs and free energies plots revealed a very good agreement between QM and MM (Figure 6). Surprisingly, the agreement is better than what seen for the smaller NC (Figure 2), in particular for the Cd–Cd RDF. The strong benefit of using a fully-classical model is easily understood by comparing the computational times for the simulation of this larger NC: only 10 seconds on 16 cores for the MM-MD against 11 days on 72 cores (E5-2690 v3, 2.6 GHz) for the ab-initio MD simulations. In the case of the large (4-nm) CdSe NC (Figure 7), our optimized FF parameters were able to maintain a stable NC core in the zinc-blende crystalline structure over time scales of hundreds of nanoseconds. The whole system was simulated for ~700 ns and no phase transformations to other crystalline structures were observed. Figure 7b shows the strong similarities among RDFs for the three different NC sizes considered in this work, despite the very different timescales of the MD simulations and the different environments (the 4-nm NC was passivated by oleates and solvated by dichloromethane). A trend towards more ordered structure is noticeable with increasing the NC size. Overall, the motions of the ligands at the surface are consistent with the outcomes (both QM and MM) from the smaller NCs. In addition, on this larger time scales, we noticed a significant diffusion of ligands between the facets. The behavior of CdSe NC on the μs time scales will be described in more details in future publications.


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 40

Figure 7. a) Starting simulation box for the MD simulation (~700 ns) of the Cd1012Se911 NC (~4.0

nm) at room temperature and pressure; the non-stoichiometric NC is passivated by 202 oleates, and solvated in dichloromethane. b) Comparison of RDFs computed for NCs of three different sizes (1.8, 2.3 and 4.0 nm); the 2 smaller NCs were passivated by formate and simulated in vacuum, the larger NC (4.0 nm) was passivated by oleates and simulated in dichloromethane at 1 atm.


24

Page 25 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We performed a further validation of the optimized FF parameters by computing, for the 1.8nm QD, the vibrational densities of states (vDOS), i.e. the power spectra of the velocity autocorrelation functions (VACF).46 Again, we compared our FF parameters with those available in the literature (Rabani) and with ab-initio MD results. Similarity among the different calculations is very high. Notably, Rabani’s FF parameters delivered good vDOSs despite being unable to properly provide a structural description of the NC. The largest deviations with ab-initio MD were observed in the oxygen vibrations associated with C-O stretching modes at 1230 and 1700 cm-1 (Figure 8), which primarily depends on the C-O bonded interactions as parameterized in CHARMM.


25


a)

b) 5

4 Spectrum (Arbitrary units)

0.9 Spectrum (Arbitrary units)

0.6

0.3

3

2

1

0.0

0 0

100

200

300

400

500

0

500

Frequency (cm−1)

1000

1500

2000

Frequency (cm−1)

c) FF: Rabani BOMD: This work FF: This work

0.9 Spectrum (Arbitrary units)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 40

0.6

0.3

0.0 0

100

200

300

400

500

Frequency (cm−1)

Figure 8. Comparison of vibrational density of states (vDOS) for (a) Se, (b) O and (c) Cd, obtained using ab-initio MD (green line), our optimized FF parameters (blue line), and Rabani’s FF paramters (orange line).

2.4 Bulk properties For practical purposes, a good description of the surface atoms is probably the most important feature of a classical FF for colloidal NCs where most of the relevant events are known to occur at the surface (reconstruction, creation of defects, ligand exchange, etc.). However, the success of such a FF parametrization also depends on its ability to provide a good representation of the


26

Page 27 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nanocrystal bulk counterpart and its transferability to different crystalline phases, e.g. zincblende and wurtzite for CdSe. In Table 3 we compared our parameters against experimental bulk properties (lattice parameters, elastic constants, and bulk moduli) and available FFs parameters by Rabani, and Yalcin et al.25 The latter uses a Buckingham potential for the Van der Waals interactions, unlike ours and Rabani’s FFs. Deviation from the experimental values are reported in Table 3 as the median absolute percentage error (MAPE).47 Our as-optimized ARMC parameters show a MAPE of 22.3%, which is better than Rabani’s (25.8%) but quite large compared to the 5.7% of Yalcin’s. This can be explained by noticing that NCs commonly employed in experiments (3-5 nm) expose a large fraction of atoms at the surface, directly bonded to surfactant ligands. In this framework, a significant charge transfer between the ligands and the NC takes place, which is implicitly represented in the QM reference calculations used in our fitting procedure. This effect is noticeable in the fitted ionic charges (q(Cd)=0.9768e) which are lower than the bulk charges fitted by Rabani (q(Cd)=1.18e) using the same FF. Remarkably, when we plug Rabani’s charges in our FF, we obtain bulk properties that exhibit the best MAPE (2.3%), indicating the high-quality of our fitted Lennard-Jones parameters. In other words, we could envision our parameters to be suitable for experimental sized NCs, which usually spans sizes up to 10 nm, while for larger NCs, when the number of atoms at the surface is much smaller than the number of core atoms charges should be adapted to match bulk values as fitted by Rabani. These results clearly demonstrate the quality of our approach that is capable of providing DFT quality result for small NCs and bulk macroscopic properties using a very simple FF potential.

Property

Exp.

48

ARMC ARMC (qCd=0.9768) (qCd=1.18)

Rabani23,24

Yalcin et al.25


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 40

Zinc Blende a

6.08

6.24

6.07

6.11

6.08

c11

66.7

41.6

65.3

86.2

61.6

c12

46.3

30.1

46.9

30.3

48.8

c44

22.3

19

30.5

36.8

20.3

B

53.1

34

53.0

48.9

53.1

Wurtzite a

4.30

4.44

4.32

4.34

4.32

c

7.01

7.01

6.91

7.02

6.94

c11

74.1

50.7

79.8

90.7

72.4

c12

45.2

28.4

43.8

30.4

47.9

c13

39.0

24.2

37.3

13.3

41.4

c33

84.3

52.1

82.1

122.0

72.5

c44

13.4

11.9

19.1

15.8

12.3

c66

14.5

13.6

18.0

30.1

12.2

B

53.1

50.2

53.1

45.6

53.1

23.2

2.3

25.8

5.7

Deviation (MAPE, %)

Table 3. Experimental and calculated bulk properties of CdSe: lattice constants a, b and c (Å), elastic constants cij and bulk moduli B (Gpa). GULP 4.4.2 was used for the calculations.49 Differently to our approach, FF of Rabani23,24 and Yalcin et al.25 were fitted to reproduce bulk properties, and the reported data are from the original manuscripts. Estimates of the deviation from the experimental data48 are provided as the median absolute percentage error (MAPE).


28

Page 29 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3. Conclusions. In this work, we have presented a FF parameterization of CdSe NCs with NC models that takes into account some specific aspects of the NC surface that have been experimentally demonstrated to exist only recently: (i) the presence of excess Cd on the NC surface that is charge compensated by ligands and (ii) the charge neutrality of the overall NC/ligands system. In this respect the FF parameterization has to be finely tuned to these experimental conditions. We thus developed an automated approach for parameterizing classical force fields from first principle calculations based on a stochastic global optimization algorithm that we called Adaptive Rate Monte Carlo (ARMC). Our optimization scheme is able to outperform similar algorithms present in the literature. Our approach is a Monte Carlo method based on an on-thefly modification of the explored surface and on the simple problem-independent parameter at, i.e. the “target average acceptance rate”. This generates an algorithm that, differently from similar approaches, always continues searching for better solutions (i.e. it is not affected by freezing and the solution is truly dependent on the number of iterations). Its application on real problems is easy and it is robust with respect to the choice of algorithm parameters. Depending on the purpose, ARMC performance can be increased by varying at according to some useful scheme. In our case, we found convenient to lower at in an annealing-like procedure (“simulated braking”, ARMC-SB). In this framework, ARMC can be considered as a serious alternative to Metropolis Monte Carlo for optimization, and can be easily associated with any layer designed to enhance MMC performances. It can be easily implemented in any Monte Carlo code and eventually straightforwardly parallelized by using multiple replicas that share the on-the-fly modification.


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 40

The quality of the refined FF parameters developed using our automatic procedure primarily depends on the accuracy of the ab-initio reference and on the significance of the sampled phase space. The application to the CdSe nanocrystal demonstrated how our method provides FF parameters able to both finely reproduce NC structural features as the ab-initio simulations and predict bulk properties by using a very simple interatomic potential (Lennard Jones and Coulomb). This FF will be used in future works to study the chemistry of CdSe QDs at the microscopic level, while applicability of the presented FF parameterization approach will be further explored on other systems, including different colloidal nanocrystals

APPENDIX: Our newly developed FF parameters were used to simulate large sized CdSe NC (4.0 nm) under realistic experimental conditions, i.e. oleate-capped and solvated in dichloromethane. The procedure to generate the simulation box is rather elaborated, thus we decided to provide a stepby-step procedure to create these models. A video that summarizes the following steps is provided in the SI (Video V4) Step 1. A Cd1012Se911 NC is cleaved out of the zinc-blend bulk to expose the (100), the Cd-rich (111), the Se-rich (111) and (101) staircase facets. The remaining parameters for the ligands were taken from CHARMM FF, while for the solvent parameters from Fennel et al.50 were used. Step 2. The distribution of the ligands at the NC surface is unknown from the experiments, therefore we placed 202 formates manually to the surface. The (100) facets were passivated first, because are composed by Cd atoms with the lowest coordination number; the remaining ligands were placed at Cd-rich (111) facets.


30

Page 31 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Step 3. The positions of the formate groups were optimized freezing the positions of NC atoms and using the conjugate gradient algorithm as implemented in GROMACS. Step 4. The hydrogen atoms of the formates were substituted with tails of oleic acid and a minimization was performed to remove the overlaps between nearby oleate tails, fixing the positions of the atoms of NC and carboxylates. The more robust minimization algorithm implemented in NAMD 2.9 (conjugate gradient and line search) was used for this step. Step 5. The system was solvated in a dichloromethane box of size 12.5×12.5×12.5 nm (16301 molecules). Relaxation of the solvated NC was performed in two steps: a) relaxation of solvent and tails with an energy minimization followed by 200 ps of MD simulation, in which atoms of NC and carboxylates were constrained; b) relaxation of oleates on the surface with 2 ns of MD simulation, in which only positions of the NC atoms were restrained with harmonic potentials. Step 6. A MD simulation was finally performed without any restraint. In this way, we were able to simulate experimentally relevant time scales (~700 ns) in a relatively short time (~85 ns/day on 96 Intel Xeon CPU E5-2690 v3, 2.60GHz). Part of the trajectory (50 ns) is shown in Video S3.

AUTHOR INFORMATION Corresponding Author * [email protected] Author Contributions


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 40

II and SMC sought the idea, while SMC developed the ARMC algorithm. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Funding Sources Ivan Infante would like to thank the Netherlands Organization of Scientific Research (NWO) for providing financial support within the Innovational Research Incentive (Vidi) Scheme. Notes Any additional relevant notes should be placed here. ACKNOWLEDGMENT The authors would like to thank Dr. Andrea Gabrieli for the useful discussions. I. I. and S.C. would like to thank the Netherlands Organization of Scientific Research (NWO) for providing financial support within the Innovational Research Incentive (Vidi) Scheme (Grant No. 723.013.002). This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.

Supporting Information: In the Supporting Information additional methodology and results are described. Videos of some MD simulations are also provided. This information is available free of charge via the Internet at http://pubs.acs.org


32

Page 33 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ABBREVIATIONS MD = Molecular Dynamics; MMC-SA = Metropolis Monte Carlo Simulated Annealing; ARMC = Adaptive Rate Monte Carlo; SB = Simulated Breaking; FF = Force Field; QM = Quantum Mechanics; MM = Molecular Mechanics; NC = Nanocrystals; DFT = Density Functional Theory REFERENCES (1)

Murray, C. B.; Norris, D. J.; Bawendi, M. G. Synthesis and Characterization of Nearly Monodisperse CdE (E = Sulfur, Selenium, Tellurium) Semiconductor Nanocrystallites. J. Am. Chem. Soc. 1993, 115 (19), 8706–8715.

(2)

Dabbousi, B. O.; Rodriguez-Viejo, J.; Mikulec, F. V; Heine, J. R.; Mattoussi, H.; Ober, R.; Jensen, K. F.; Bawendi, M. G. (CdSe)ZnS Core−Shell Quantum Dots: Synthesis and Characterization of a Size Series of Highly Luminescent Nanocrystallites. J. Phys. Chem. B 1997, 101 (46), 9463–9475.

(3)

Sargent, E. H. Colloidal Quantum Dot Solar Cells. Nat. Photonics 2012, 6 (3), 133–135.

(4)

Chen, O.; Zhao, J.; Chauhan, V. P.; Cui, J.; Wong, C.; Harris, D. K.; Wei, H.; Han, H.-S.; Fukumura, D.; Jain, R. K.; Bawendi, M. G. Compact High-Quality CdSe–CdS Core–shell Nanocrystals with Narrow Emission Linewidths and Suppressed Blinking. Nat. Mater. 2013, 12 (5), 445–451.

(5)

Roo, J. De; Driessche, I. Van; Martins, J. C.; Hens, Z. Colloidal Metal Oxide Nanocrystal Catalysis by Sustained Chemically Driven Ligand Displacement. Nat. Mater. 2016, 15 (January), 1–6.


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(6)

Page 34 of 40

Sperling, R. a; Parak, W. J. Surface Modification, Functionalization and Bioconjugation of Colloidal Inorganic Nanoparticles. Philos. Trans. A. Math. Phys. Eng. Sci. 2010, 368 (1915), 1333–1383.

(7)

Knowles, K. E.; Frederick, M. T.; Tice, D. B.; Morris-Cohen, A. J.; Weiss, E. A. Colloidal Quantum Dots: Think Outside the (Particle-in-a-)Box. J. Phys. Chem. Lett. 2012, 3 (1), 18–26.

(8)

Giansante, C.; Infante, I.; Fabiano, E.; Grisorio, R.; Suranna, G. P.; Gigli, G. “Darkerthan-Black” PbS Quantum Dots: Enhancing Optical Absorption of Colloidal Semiconductor Nanocrystals via Short Conjugated Ligands. J. Am. Chem. Soc. 2015, 137 (5), 1875–1886.

(9)

Azpiroz, J. M.; De Angelis, F. Ligand Induced Spectral Changes in CdSe Quantum Dots. ACS Appl. Mater. Interfaces 2015, 7 (35), 19736–19745.

(10)

Azpiroz, J. M.; Matxain, J. M.; Infante, I.; Lopez, X.; Ugalde, J. M. A DFT/TDDFT Study on the Optoelectronic Properties of the Amine-Capped Magic (CdSe)13 Nanocluster. Phys. Chem. Chem. Phys. 2013, 15 (26), 10996.

(11)

Voznyy, O.; Sargent, E. H. Atomistic Model of Fluorescence Intermittency of Colloidal Quantum Dots. Phys. Rev. Lett. 2014, 112 (15), 157401.

(12)

Voznyy, O.; Thon, S. M.; Ip, A. H.; Sargent, E. H. Dynamic Trap Formation and Elimination in Colloidal Quantum Dots. J. Phys. Chem. Lett. 2013, 4 (6), 987–992.

(13)

Hwang, G. W.; Kim, D.; Cordero, J. M.; Wilson, M. W. B.; Chuang, C.-H. M.; Grossman, J. C.; Bawendi, M. G. Identifying and Eliminating Emissive Sub-Bandgap States in Thin


34

Page 35 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Films of PbS Nanocrystals. Adv. Mater. 2015, 27 (30), 4481–4486. (14)

Smith, A. M.; Nie, S. Semiconductor Nanocrystals: Structure, Properties, and Band Gap Engineering. Acc. Chem. Res. 2010, 43 (2), 190–200.

(15)

Luther, J. M.; Pietryga, J. M. Stoichiometry Control in Quantum Dots: A Viable Analog to Impurity Doping of Bulk Materials. ACS Nano 2013, 7 (3), 1845–1849.

(16)

Morris-Cohen, A. J.; Frederick, M. T.; Lilly, G. D.; McArthur, E. A.; Weiss, E. A. Organic Surfactant-Controlled Composition of the Surfaces of CdSe Quantum Dots. J. Phys. Chem. Lett. 2010, 1 (7), 1078–1081.

(17)

Fritzinger, B.; Capek, R. K.; Lambert, K.; Martins, J. C.; Hens, Z. Utilizing Self-Exchange To Address the Binding of Carboxylic Acid Ligands to CdSe Quantum Dots. J. Am. Chem. Soc. 2010, 132 (29), 10195–10201.

(18)

Moreels, I.; Lambert, K.; De Muynck, D.; Vanhaecke, F.; Poelman, D.; Martins, J. C.; Allan, G.; Hens, Z. Composition and Size-Dependent Extinction Coefficient of Colloidal PbSe Quantum Dots. Chem. Mater. 2007, 19 (25), 6101–6106.

(19)

Moreels, I.; Fritzinger, B.; Martins, J. C.; Hens, Z. Surface Chemistry of Colloidal PbSe Nanocrystals. J. Am. Chem. Soc. 2008, 130 (45), 15081–15086.

(20)

Kilina, S.; Velizhanin, K. A.; Ivanov, S.; Prezhdo, O. V.; Tretiak, S. Surface Ligands Increase Photoexcitation Relaxation Rates in CdSe Quantum Dots. ACS Nano 2012, 6 (7), 6515–6524.

(21)

Kilina, S. V.; Neukirch, A. J.; Habenicht, B. F.; Kilin, D. S.; Prezhdo, O. V. Quantum


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 40

Zeno Effect Rationalizes the Phonon Bottleneck in Semiconductor Quantum Dots. Phys. Rev. Lett. 2013, 110 (18), 180404. (22)

Kilina, S. V.; Kilin, D. S.; Prezhdo, V. V.; Prezhdo, O. V. Theoretical Study of Electron– Phonon Relaxation in PbSe and CdSe Quantum Dots: Evidence for Phonon Memory. J. Phys. Chem. C 2011, 115 (44), 21641–21651.

(23)

Rabani, E. An Interatomic Pair Potential for Cadmium Selenide. J. Chem. Phys. 2002, 116 (1), 258.

(24)

Rabani, E. Structure and Electrostatic Properties of Passivated CdSe Nanocrystals. J. Chem. Phys. 2001, 115 (3), 1493.

(25)

Yalcin, A. O.; Fan, Z.; Goris, B.; Li, W.-F.; Koster, R. S.; Fang, C.-M.; van Blaaderen, A.; Casavola, M.; Tichelaar, F. D.; Bals, S.; Van Tendeloo, G.; Vlugt, T. J. H.; Vanmaekelbergh, D.; Zandbergen, H. W.; van Huis, M. A. Atomic Resolution Monitoring of Cation Exchange in CdSe-PbSe Heteronanocrystals during Epitaxial Solid–Solid– Vapor Growth. Nano Lett. 2014, 14 (6), 3661–3667.

(26)

Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell, A. D. CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2009, 31 (4), 671-690.

(27)

Wang, F.; Landau, D. P. Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States. Phys. Rev. Lett. 2001, 86 (10), 2050–2053.

(28)

Boehme, S. C.; Azpiroz, J. M.; Aulin, Y. V.; Grozema, F. C.; Vanmaekelbergh, D.;


36

Page 37 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Siebbeles, L. D. A.; Infante, I.; Houtepen, A. J. Density of Trap States and AugerMediated Electron Trapping in CdTe Quantum-Dot Solids. Nano Lett. 2015, 15 (5), 3056– 3066. (29)

Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press, New York, United States, 1989.

(30)

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865–3868.

(31)

Hutter, J.; Iannuzzi, M.; Schiffmann, F.; Vandevondele, J. Cp2k: Atomistic Simulations of Condensed Matter Systems. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (1), 15–25.

(32)

Lippert, B. G.; Parrinello, J. H. and M. A Hybrid Gaussian and Plane Wave Density Functional Scheme. Mol. Phys. 1997, 92 (3), 477–488.

(33)

Azpiroz, J. M.; Ugalde, J. M.; Infante, I. Benchmark Assessment of Density Functional Methods on Group II–VI MX (M = Zn, Cd; X = S, Se, Te) Quantum Dots. J. Chem. Theory Comput. 2014, 10 (1), 76–89.

(34)

Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126 (1), 14101.

(35)

Mackerell, A. D. Empirical Force Fields for Biological Macromolecules: Overview and Issues. J. Comput. Chem. 2004, 25 (13), 1584–1604.

(36)

Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 40

Comput. 2008, 4 (3), 435–447. (37)

Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781–1802.

(38)

Ivanov, M. V.; Talipov, M. R.; Timerghazin, Q. K. Genetic Algorithm Optimization of Point Charges in Force Field Development: Challenges and Insights. J. Phys. Chem. A 2015, 119 (8), 1422–1434.

(39)

McQuarrie, D. Statistical Mechanics, 1st Editio.; University Science Book: Sausalito, CA, 2000.

(40)

Doemer, M.; Liberatore, E.; Knaup, J. M.; Tavernelli, I.; Rothlisberger, U. In Situ Parameterisation of SCC-DFTB Repulsive Potentials by Iterative Boltzmann Inversion. Mol. Phys. 2013, 111 (22–23), 3595–3607.

(41)

Reith, D.; Pütz, M.; Müller-Plathe, F. Deriving Effective Mesoscale Potentials from Atomistic Simulations. J. Comput. Chem. 2003, 24 (13), 1624–1636.

(42)

Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220 (4598), 671–680.

(43)

Kirkpatrick, S. Optimization by Simulated Annealing: Quantitative Studies. J. Stat. Phys. 1984, 34 (5–6), 975–986.

(44)

Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism


38

Page 39 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from Laptops to Supercomputers. SoftwareX 2015, 1–2, 19–25. (45)

Parrinello, M.; Rahman, A. Crystal Structure and Pair Potentials: A Molecular-Dynamics Study. Phys. Rev. Lett. 1980, 45 (14), 1196–1199.

(46)

Cosseddu, S. M.; Khovanov, I. A.; Allen, M. P.; Rodger, P. M.; Luchinsky, D. G.; McClintock, P. V. E. Dynamics of Ions in the Selectivity Filter of the KcsA Channel. Eur. Phys. J. Spec. Top. 2013, 222 (10), 2595–2605.

(47)

Armstrong, J. S.; Collopy, F. Error Measures for Generalizing about Forecasting Methods: Empirical Comparisons. Int. J. Forecast. 1992, 8 (1), 69–80.

(48)

Adachi,

S.

Properties

of

Semiconductor

Alloys:

Group-IV,

III-V

and

II-VI

Semiconductors; John Wiley and Sons, Ltd., Chichester, United Kingdom, 2009. (49)

Gale, J. D.; Rohl, A. L. The General Utility Lattice Program ( GULP ). Mol. Simul. 2003, 29 (5), 291–341.

(50)

Fennell, C. J.; Li, L.; Dill, K. A. Simple Liquid Models with Corrected Dielectric Constants. J. Phys. Chem. B 2012, 116 (23), 6936–6944.


39


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 40

For Table of Contents use only

A Force Field Parametrization of Colloidal CdSe Nanocrystals Salvatore Cosseddu# and Ivan Infante#


40

Force Field Parametrization of Colloidal CdSe Nanocrystals Using an

Recommend Documents