Molecular Dynamics Simulations Using a CapacitanceâPolarizability

Page 1 of 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

The Journal of Physical Chemistry

Molecular Dynamics Simulations using a Capacitance−Polarizability Force Field Xin Li,*,† Hans Ågren† †

Division of Theoretical Chemistry and Biology, School of Biotechnology, KTH Royal Institute

of Technology, SE-10691 Stockholm, Sweden KEYWORDS. Molecular dynamics, capacitance−polarizability force field, gold nanoparticles.

ABSTRACT. We present molecular dynamics (MD) simulations using a capacitance– polarizability force field. This force field allows an atomistic description of charge migration within a particle and hence the image charge effects at the interface of such a particle. By employing atomic capacitance and polarizability as the key parameters that describe fluctuating charges and dipoles we can thus explore the effect of charge migration on the structural dynamics. We illustrate the method by exploring gold nanoparticles in aqueous solutions and compare with previous simulation work. We reach the conclusion that the capacitance– polarizability force field MD method serves as a promising tool for simulating gold–water systems, indicating probable extensions to other metal solutions and for studies of more complicated systems provided that a proper parameterization of the capacitance force field can be made. For the particular system studied, it is found that the water molecules interact with the

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 33

surface through oxygen atoms, leading to more hydrogen bond donors than acceptors at the gold−water interface. A prominent shoulder peak is found in the radial distribution of oxygen atoms with respect to the gold surface, owing to the fact that the oxygen atoms adsorbed at the on-top sites of the gold nanoparticle. The surface of the aqueous gold nanoparticle carries negative charge, which is balanced by the positive charge in the second outermost layer.

Introduction Noble metal nanoparticles are valuable not only because of their rareness, but also owing to the salient and intriguing physicochemical properties of their interfaces in contact with other chemical species. Nowadays, noble metal surfaces and nanoparticles have been widely employed as an important part of functional materials with many applications in catalysis, bioimaging and biosensing.1-6 In particular, these nanoparticles can provide ordered orientation of the adsorbed chromophores and a local field enhancement for optical response7-9 − two characteristics that offer combined functionality and convenient detection of events taking place in the adsorbent. Rational design and development of such novel functional materials with desirable properties calls for theoretical modelling of the underlying interaction between the metal nanoparticles and the surrounding environment. Unfortunately, the interface of noble metal nanoparticles is in general large, complex and coupled with the almost-free motions of the electrons, which give rise to the image charge effects at the metal surfaces. This poses challenges in theoretical modelling of the interface formed by the noble metal nanoparticles and the non-metal surroundings.


2

Page 3 of 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


We have recently introduced linear and nonlinear response theory with heterogeneous environments in a multiscale quantum mechanical/classical mechanical (QM/MM) framework,1013

which allows simulations of optical properties of organic chromophores physisorbed on metal

surfaces. It is our aim to further develop a molecular dynamics (MD) force field to allow efficient sampling of configurations in such calculations. Owing to the complex interaction at the metal–nonmetal interfaces, quantum chemical calculations can provide useful information for the development of a practical force field for the purpose. For instance, Cicero and co-workers14 reported ab initio MD simulations of an aqueous Au(111) surface. It was found that the oxygen atoms of the water molecules adsorbed at the on-top site of the gold atoms, leading to a prominent shoulder peak in the distribution function of the oxygen atom. The deficiencies in some available classical force fields were also pointed out for simulating gold−water interfaces, that is, i) the orientation of hydrogen atoms at the surface, and ii) the adsorption site of oxygen atoms on the surface. It was reported that the GolP force field15,16 gives an erroneous orientation of adsorbed water where a hydrogen atom is pointing to the gold surface, due to the absence of van der Waals interaction between the gold and the hydrogen atoms. Owing to the same reason the METAL force field proposed by Heinz and co-workers17 cannot describe the hydrogen bond structures at the gold surface. In addition, the METAL force field was found to predict hollow adsorption sites for oxygen atoms on gold surface, which is inconsistent with quantum chemical calculations. Recently, Walsh and co-workers developed the GolP-CHARMM force field18 for gold surfaces with parameterized image charge effects. This force filed overcomes the aforementioned problems in the orientation of hydrogen atoms and adsorption of oxygen atoms, and explicitly takes into account the polarization effects on gold surfaces by the surrounding molecules. However, the applicability of this force field is limited to specific surfaces, namely


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 33

Au(111) and Au(100) surfaces, owing to the employment of virtual interacting sites. Simulations of nanoparticles are therefore not straightforward with this force field. On the other hand, the reactive force field (ReaxFF) for Au–S–C–H systems19 provides a reasonable description of thiol-protected gold surface, which, unfortunately, lack universal applicability due to the systemdependent interacting parameters and the tedious parameterization procedure for introducing parameters for particular atoms. A practical direction for the development of a proper force field for noble metal surfaces and nanoparticles is to employ different interacting potential functions for the metallic and non-metallic components and treat the interaction between them carefully. Towards proper treatment of the almost-free motions of the electrons in metallic materials, Gaussian-distributed charges have been employed in the simulation of electrostatic interaction in water/electrode systems.20 During the simulations the charges are determined at each step by minimizing the potential energy of the system and maintaining the equipotential character of the metallic electrode. In practice the fluctuation of the induced charges can also be coupled with the motions of the atoms by means of an extended Lagrangian technique. A similar approach21 has also been applied in the simulations of ionic liquid in contact with metallic and graphite electrodes.22 Other methods describing the charge migration within molecules and/or nanoparticles have also been reported, for instance the charge equilibration (QEq) approach23 and the electronegativity equalization method (EEM).24 In particular, the EEM approach, which has been implemented in ReaxFF,25 employs electronegativity and hardness as parameters that determine the atomic partial charges at a given geometry. Based on the electronegativity equalization principle, a system of linear equations can be solved under the constraint of the net charge of the molecule. This charge-fluctuating approach has been successfully applied in the simulations of gold−thiol systems.19 Recently, Mayer further incorporated induced dipoles into


4

Page 5 of 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


the system and proposed the charge–dipole model26 for carbon nanotubes. This model was further developed by Jensen and co-workers as the capacitance–polarizability interaction model27,28 for gold and silver nanoclusters. In this paper we present development of force field for simulations of metal nanoparticles in the presence of surrounding non-metallic molecules by adopting a similar treatment of electrostatic interaction to that reported by Jensen and co-workers, where atomic capacitance and polarizability are employed as key parameters that describe the energy cost of creating induced charges and dipoles at the metal atoms. The atomic capacitance together with polarizability describes the response of the metal nanoparticle to an external electric field, the physical basis of which differs from the partial charges used in common parameterization of MD force fields. In the absence of external electric field, each metal atom in an electrically neutral nanoparticle carries zero charge, while the partial charges in an organic molecule persist due to its own electronic structure. We therefore aim to treat the electrostatic response of the metallic part using the capacitance–polarizability interaction model, while the non-metallic part is described by conventional force fields.

Theory The capacitance–polarizability interaction model27 describes the induced dipoles and charges of an arbitrary-shaped metal surface/nanoparticle in the presence of external electric field and potential. The induced dipoles and charges are obtained by solving the matrix problem  A  T − M  0 

− M 0  p   E      − C 1  q  =  V  1 0  λ   Q 

(1)


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 33

where p is the induced dipoles, q is the induced charges, λ is the Lagrangian multiplier, E is the electric field from the surrounding non-metallic molecules plus external electric fields, V is the corresponding electrostatic potential, and Q is the net charge of the nanoparticle. The sub matrices A, M and C describe the dipole−dipole, dipole−charge and charge−charge interaction among metal atoms, which are expressed as

(i ≠ j ) M ij ,α = Tij(,1α) (i ≠ j ) Cij = Tij( 0) (i ≠ j )

1 Aii ,αβ = α i−,αβ ; Aij ,αβ = −Tij(,2αβ)

M ii ,α = 0; Cii = ci−1 ;

(2)

The charge−charge, dipole−charge and dipole−dipole interaction between different metal atoms are 1 r Tij( 0) = erf  , r R Tij(,1α) = ( 2) ij ,αβ

T

r = rij , R = Rijqq

 r 2  rα   r  2r   − 2  , erf − exp    r 3   R  πR  R 

r = rij , R = Rijpq

3rα rβ − δ αβ r 2   r   r2  2r   − 2   erf exp = −   R   r5 π R  R   4rα rβ  r2   − 2 , r = rij , R = Rijpp − exp 3 2 πR r  R 

(3)

where the R parameters are obtained from the individual atomic capacitance and atomic polarizability parameters, c and α


6

Page 7 of 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


Rijqq = ( Riq ) 2 + ( R qj ) 2 Rijpq = ( Rip ) 2 + ( R qj ) 2 Rijpp = ( Rip ) 2 + ( R jp ) 2 Riq =

2

π

(4)

ci 1/ 3

 2 αi   Ri =   π 3   p

The above formulae are derived based on Gaussian distributions of induced dipoles and charges, as proposed by Mayer.26 During MD simulations Eq. (1) is solved at each step, using the biconjugate gradient stabilized (BiCGSTAB) method.29 The BiCGSTAB method is an iterative approach, such that the solved charges and dipoles at the previous MD step can be used as the initial guess for the following MD step. The computational complexity of solving Eq. (1) is O(N2), and future works will be focused on boosting the computational efficiency by taking advantage of the sparsity of the matrix, spatial distributions of the charges and extended Lagrangian dynamics. The BiCGSTAB approach has been parallelized and the computational time can be reasonably reduced through parallelization (see Figure S1 in the Supporting Information). The forces from the capacitance–polarizability interaction model can be calculated by differentiating the potential energy

1 1 VCPIM = q TCq + p T Ap − p T Mq + q T V − p T E 2 2

(5)

The derivative of the potential energy with respect to Cartesian coordinates is expressed as


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

∂VCPIM 1 T ∂C 1 ∂A ∂M ∂V ∂E = q q + pT p − pT q + qT − pT ∂r 2 ∂r 2 ∂r ∂r ∂r ∂r

Page 8 of 33

(6)

where only the explicit dependences of A, M, C, E and V on r exist, and the implicit dependences of p and q on r disappear owing to the stationary condition required by equation (1). The forces are obtained as the negative derivative of the potential energy

FCPIM = −

∂VCPIM ∂r

(7)

and the derivatives of the charge−charge, dipole−charge and dipole−dipole interaction are

∂Tij( 0) ∂rij ,ε ∂Tij(,1α) ∂rij ,ε

 r2  rε   r  2r  = − 3  erf   − exp − 2  , r = rij , R = Rijqq r  R πR  R   r 2  3rα rε − r 2δ αε   r  2r    =−  erf  R  − π R exp − R 2   r5     r2  + exp − 2 , π R 3r 2  R  4rα rε

∂Tij(,2αβ) ∂rij ,ε

r = rij , R = Rijpq

(8)

5rα rβ rε − r 2 (rα δ βε + rβ δ αε + rε δ αβ )   r2  4r 3   r   6r      =−  3erf  R  −  π R + π R 3  exp − R 2   r7      2 8r r r  r  + α β5 ε 2 exp − 2 , r = rij , R = Rijpp πR r  R 

where α, β and ε denote Cartesian components. For long-range electrostatic interaction the damped shifted force approach30 was employed, which is a variant of the Wolf summation method.31,32 In this DSF approach the electrostatic potential and electric field due to a point charge k is expressed as


8

Page 9 of 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


 erfc( wr )

ϕ ext ,k = qk  Eαext ,k

−

 erfc( wRc )  erfc( wRc ) 2 w exp(− w 2 Rc2 )  (r − Rc ) +  + 2  Rc Rc Rc π   

r   erfc( wr ) 2w exp(− w 2 r 2 )  erfc( wRc ) 2w exp(− w 2 Rc2 )   rα   = qk  + −  + 2 2  r r R R π π c c   r 

(9)

where Rc is the cutoff radius and w is the damping parameter. The derivatives of the electrostatic potential and the electric field with respect to Cartesian coordinates are

∂ϕ ext , k ∂rε ∂Eα ∂rε

ext , k

 erfc( wr ) 2 w exp(− w 2 r 2 )  erfc( wRc ) 2 w exp(− w 2 Rc2 )   rε   = qk  − − +  + 2 2  r r R R π π c c   r   3rα rε − r 2δ αε  4 w3 rα rε 2 wr 2 2  2 2  − erfc ( wr ) + exp ( − w r ) − exp ( − w r )   r5 π π r2     = qk   2 2 2  +  erfc( wRc ) + 2 w exp(− w Rc )  rα rε − r δ αε  2 3     R R r π c c    

(10)

While in this paper we focus on single metallic nanoparticle in aqueous environment, we would like to mention the possibility of treating multiple nanoparticles. The capacitance−polarizability interaction model described above allows charges to migrate within one nanoparticle, where the induced charges and dipoles are determined by minimizing the potential energy of the system (e.g. solving matrix equation (1)). In case of well-separated systems or multiple nanoparticles, equation (1) needs to be modified to take into account the additional metallic clusters. For instance, if there are two nanoparticles, two Lagrangian multipliers will be needed to constrain the net charges, forming a matrix equation of 4n+2 dimension (instead of 4n+1 as in equation (1) where n is the number of metallic atoms)


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

 A  T − M  0   0 

−M

0

−C

X1

T 1 T 2

X

0

X

0

0  p   E      X 2  q   V  = , 0  λ1   Q1      0  λ2   Q2 

1   M 1 X1 =  , 0 M   0  

Page 10 of 33

0   M 0 X2 =   1 M   1  

(11)

Here λ1 and λ2 are the Lagrangian multipliers, Q1 and Q2 are the net charges of the nanoparticles, and the non-zero elements in vectors X1 and X2 correspond to the metal atoms in the first and the second nanoparticle, respectively. This matrix equation can be easily expanded to describe multiple nanoparticles, where the net charge of each nanoparticle is constrained. This provides a correct description for multiple nanoparticles and eliminates spurious long-range motion of charges that may take place when only one Lagrangian multiplier is used. Our simulation code CapacMD, developed for molecular dynamics simulations using capacitance−polarizability force field, is available on the GitHub platform https://github.com/recoli/CapacMD/ and on our webpage http://www.theochem.kth.se/~lixin/capacmd/.

Computational Details We take the gold−water interface as a prototypical system for illustration of our approach. The interaction in the nanoparticle−water system is divided into three parts: i) interaction within the metallic nanoparticle, ii) interaction between the metallic nanoparticle and the non-metallic region, and iii) interaction within the non-metallic region. For the interaction within the metallic nanoparticle, the quantum Sutton-Chen potential33 is known to offer an accurate description of


10

Page 11 of 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


many properties of metals. The interaction within the non-metallic region, which is water in this study, can be well described by the conventional TIP3P model.34 The most complex task is to properly describe the interaction between the metallic and non-metallic regions, which consists of electrostatic and van der Waals contributions. In this work, the electrostatic interaction is described by the capacitance−polarizability interaction model,27 which is introduced by Jensen and co-workers, and which is able to describe the almost-free motions of electrons within the nanoparticle. The model was recently generalized to heterogeneous environments and included in multiscale QM/MM property calculations by the present authors.10-13 The model consists of two parameters for gold nanoparticles, the atomic capacitance and the atomic polarizability, which are parameterized as cAu = 1.2159 a.u. and αAu = 39.5297 a.u., respectively.28 The implementation details of the electrostatic interaction and the corresponding forces are presented in the Theory section. For van der Waals interaction the Lennard-Jones (LJ) potential function has been widely used in classical MD simulations. The LJ potential consists of an r−6 term describing the dispersion interaction and an r−12 term describing the repulsion. A special character of the repulsive r−12 term is the steepness, which makes the particles behave similarly to hard spheres at repulsive interacting distances, and it has been reported that simply applying the LJ potential for gold−water interaction leads to hollow adsorption sites of oxygen atoms on gold surfaces.18 Efforts have been made to reproduce the correct adsorption site of water molecules on metal surfaces, by introducing virtual interacting particles on the surface18 or employing many-body interaction potential functions.20 These approaches, however, introduce extra complexity and hamper the flexibility and transferability of the force field for future applications in different systems. We thus aim to use a new type of van der Waals potential function with softer


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 33

repulsion, which properly describes water adsorption on gold and retains the simple pair-additive character of the potential function. We optimized the Au−O and Au−H van der Waals interaction potential functions, starting from the LJ parameters of the GolP-CHARMM force field.18 The potential energy curves from these parameters are shown in Figure 1. To enhance the softness of the repulsive region, we employed the Buckingham potential function and derived the interacting parameters from the available LJ parameters. The Buckingham parameters were refined such that the MD simulations provide comparable radial distribution functions and hydrogen-bond structures to those from ab initio simulations.14 From Figure 1 it can be seen that the repulsive part of the Buckingham potential is less steep than the LJ potential. Unfortunately, the Buckingham potential predicts hollow adsorption site of oxygen atoms on gold surface, which is inconsistent with previously reported quantum chemical calculations. We therefore further increased the softness of the repulsive region by scaling the LJ potential function by an error function with a modulation parameter RMod  C12 C6   r E ModLJ =  12 − 6 erf  r   RMod r

  

(11)

This modified LJ potential function has one extra parameter compared with conventional LJ potential function, and retains the r−6 behavior at long distances. Moreover, the softness of the repulsive region can be tuned by the RMod parameter. To our delight, the modified LJ potential function is able to reproduce the radial distribution functions, hydrogen bond structures, and adsorption site of oxygen atoms on gold surface. As shown in Figure 1, the finally optimized Au−O and Au−H van der Waals potential functions mainly show repulsive character, and are


12

Page 13 of 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


much softer than the conventional LJ and Buckingham potential functions. The van der Waals parameters used in our simulations are listed in Table 1.

Figure 1. Comparison among the Lennard-Jones, Buckingham and modified Lennard-Jones potential functions for (a) Au−O and (b) Au−H interaction.

Table 1. The van der Waals parameters for water−gold interaction. Interaction

C12 / kJ mol−1 nm12

C6 / kJ mol−1 nm6

RMod / nm

Au−O

3.625 × 10−6

2.522 × 10−3

2.25

Au−H

1.455 × 10−6

1.029 × 10−3

2.25

We illustrate our approach by carrying out MD simulation of an icosahedral gold nanoparticle consisting of 309 atoms, solvated by 3866 water molecules in a cubic box of 4.916 × 4.916 × 4.916 nm3. Periodic boundary conditions were applied during the simulation, and the short-range


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 33

cutoff radius was set to 1.2 nm. Long-range electrostatic interaction was taken into account by the damped shifted force approach30 with a damping parameter of 0.2 nm–1, which have been shown to give very good agreement with the Ewald summation in a recent benchmark study.35 A time-step of 2 fs was employed for the simulation, and the geometry of the water molecule was constrained by the RATTLE algorithm.36 The simulation was conducted for 1.5 ns under an NVT ensemble with T = 300 K, and the last 1-ns trajectory was employed in statistical analysis.

Results and Discussions A snapshot of the gold nanoparticle and surrounding water molecules is shown in Figure 2a, where the water molecules form intermolecular hydrogen bond networks. We calculated the average number of hydrogen bonds per water molecule with respect to the distance from the surface of the gold nanoparticle, as shown in Figure 2b. Noteworthy, the hydrogen bond donor shows a peak at around 2.1 Å from the nanoparticle surface, while very few acceptors appear at this distance. This indicates that the adsorbed water molecules are mostly acting as hydrogen bond donors, donating their hydrogen atoms to other water molecules, while directly interacting with the gold surface through the oxygen atoms. The average number of hydrogen bond donor shows a trough at around 3.2 Å, where the average number of hydrogen bond acceptor reaches its plateau. At this distance the interaction between water molecules and the gold surface becomes weaker, such that a water molecule may have one of its hydrogen atoms pointing at the gold atoms. The effect of the gold surface on the hydrogen bonds among water molecules becomes negligible at distances larger than 5 Å. The characteristics in the hydrogen bond


14

Page 15 of 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


network of adsorbed water molecules in general show good agreement with the previously reported ab initio MD simulations of an aqueous Au(111) surface.14

Figure 2. (a) A snapshot of a gold nanoparticle with surrounding water molecules. (b) Average number of hydrogen bonds per water molecule with respect to the distance from the surface of the gold nanoparticle.

The radial distributions of oxygen and hydrogen atoms in water molecules were calculated with respect to the surface of the gold nanoparticle, as shown in Figure 3a. Noteworthy, the radial density of oxygen shows a shoulder peak at ~2.5 Å, corresponding to the direct interaction between oxygen atoms and the surface gold atoms. This is in agreement with the previously reported ab initio MD study of aqueous Au(111) surface, and such a prominent shoulder peak is absent in other classical force fields for metal surfaces.14 The main peak of the radial distribution of oxygen atoms appears at ~3.3 Å, and the first trough is located at ~5.0 Å. The radial distribution of hydrogen atoms with respect to the gold surface shows a broader peak at around 3.1 Å, and no shoulder peak is present for hydrogen. It is the combined force field, consisting of


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 33

the capacitance−polarizability interaction model and the modified Lennard-Jones potential, that provides a proper description between the water molecules and the gold nanoparticle. The modified Lennard-Jones potential is mainly of repulsive character (Figure 1), while the attractive interaction largely arises from electrostatic interaction between the permanent charges of water and the induced charges and dipoles of the gold nanoparticle. Without the modified LennardJones potential, water molecules will crash into the nanoparticle; while without the capacitance−polarizability interaction model, the water−gold interaction will suffer from repulsive forces, leading to incorrect radial distributions of water molecules as shown in Figure S2. We also calculated the orientation of the water molecules in terms of two angles, the Au−O−M and the Au−O−H angles, where M denotes the midpoint of the two hydrogen atoms in one water molecule. As shown in Figure 3b, at the oxygen−gold distance of 2.5 Å the distribution of the Au−O−M angle shows a prominent peak at around 150°, reflecting the direct interaction between the oxygen atoms and the gold atoms at the surface of the nanoparticle. Moreover, the Au−O−H angles show preference at around 120°, where the O−H bonds are pointing away from the gold surface. At a larger oxygen−gold distance of 3.3 Å, the distribution of the angles show completely different patterns from those at 2.5 Å. Both the distributions of Au−O−M and Au−O−H become very broad, indicating that the orientation of the water molecules are subject to significant tumbling, since the interaction between the gold surface and the oxygen atoms is relatively weak at this distance.


16

Page 17 of 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


Figure 3. (a) Radial distributions of water molecules and (b) distribution of angles with respect to the surface of the gold nanoparticle. M denotes the midpoint of the two hydrogen atoms, and r denotes the distance (in Å) between the oxygen atom and the gold surface.

Since the induced charges of the gold atoms are available from the capacitance−polarizability interaction model, we examined the average charges in different layers of atoms. A snapshot showing the induced charges of the gold nanoparticle is shown in Figure 4a. The icosahedral gold nanoparticle consists of a central atom and four successive layers, and the averaged charge per atom in each layer is shown in Figure 4b. We can see that the average charge per atom in the surface layer (layer 4) has a small negative value and a very large standard deviation, indicating that the induced charges of the surface atoms are significantly fluctuating due to the presence of the surrounding water molecules. The second outermost layer (layer 3 in Figure 4b) has a smaller positive average charge per atom, which balances the negative charges in the surface layer. The standard deviation in the second outermost layer is much smaller than that of the surface layer. The inner layers are almost neutral with even smaller standard deviations, reflecting the equipotential character of the bulk of the gold nanoparticle.


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

Page 18 of 33

Figure 4. (a) Surface charges of the gold nanoparticle represented by red (positive charge) and blue (negative charge) colors. (b) Average charge per atom and standard deviation for different layers of the gold nanoparticle; layer 0 denotes the central atom and layer 4 corresponds to the surface. (c) Time evolution of total surface charge of the gold nanoparticle (gray line) and the smoothed curve by adjacent averaging (blue line). (d) Distribution of positive and negative induced charges in the surface layer of the gold nanoparticle.

Interestingly, the gold nanoparticle accumulates negative charges at the surface and positive charges in the second outermost layer. We have calculated the time evolution of the total charges


18

Page 19 of 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


of the surface, as shown in Figure 4c. It can be seen that the surface charge is fluctuating with a statistical average of −0.228 ± 0.081 e, as calculated from the last 1-ns trajectory. We further calculated the distributions of positive and negative induced charges of the atoms in the surface layer of the nanoparticle, as shown in Figure 4d. Clearly, the population of negative charges is larger for small charges with |q| < 0.056 e, while the population of positive charges becomes larger for more significant charges. This is in accordance with the snapshot shown in Figure 4a, where several gold atoms carrying relatively large positive charges can be identified, and where more atoms were found to carry moderate negative charges. This is attributed to the larger negative charges of the oxygen atoms in water molecules, which are outnumbered by the hydrogen atoms which carry smaller positive charges. Finally, we inspected the adsorption site of water molecules on the surface of the gold nanoparticle. The projection of oxygen and hydrogen atoms on the facet of the icosaheral nanoparticle is plotted in Figure 5. We can see that the projection of oxygen atoms appears at the position of the surface gold atoms, confirming that the oxygen atoms prefer on-top adsorption sites. The projection of hydrogen atoms, on the other hand, exhibits distributions around the gold atoms. This is in accordance with the spatial distribution function of water molecules on Au(111) surface reported in an ab initio MD simulation study.14 Therefore, correct adsorption sites of water molecules are obtained by our approach, which employs a capacitance−polarizability force field and a modified Lennard-Jones potential function, without using any virtual atoms. This offers transferability of the force field in future parameterization and simulation of more complicated organic−inorganic systems where irregular metallic surfaces are interacting with complex environment consisting of solvents and biochemical species.


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 33

Figure 5. Projection of (a) oxygen atoms and (b) hydrogen atoms within 2.6 Å from the facet of the icosahedral gold nanoparticle.

Conclusions In summary, we have presented molecular dynamics simulations by means of a capacitance−polarizability force field and used a gold nanoparticle in aqueous solution as an illustrating example of the method. The force field divides the water−gold interaction into three categories, which are tackled by different potential functions. The TIP3P water model, the quantum Sutton-Chen potential, the capacitance−polarizability interaction model, and a modified Lennard-Jones potential compose the force field, providing promising results for the interfacial structure of aqueous gold nanoparticles. In particular, the oxygen atoms were found to adsorb at the on-top sites of the gold atoms, leading to more hydrogen bond donors than acceptors in the vicinity of the gold surface. This is also reflected by the prominent shoulder peak, located at around 2.5 Å, in the radial distribution of oxygen atoms with respect to the gold surface. The main peak of the radial density of oxygen appears at around 3.3 Å, where the interaction between


20

Page 21 of 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


the oxygen atoms and the gold surface becomes weaker, as reflected by the broad distribution of the Au−O−M and Au−O−H angles, where M is the midpoint of the two hydrogen atoms. Moreover, the surface of the aqueous gold nanoparticle was found to carry negative charges, which are balanced by the positive charges from the second outermost layer of atoms. Although the oxygen atoms induce relatively large positive charges at the surface gold atoms, the hydrogen atoms, which are greater in number, lead to a larger population of moderate negative charges at the surface of the gold nanoparticle. Through comparison with a previous ab initio molecular dynamics study, we conclude that the capacitance−polarizability force field presented in this work offers an accurate simulation method for gold−water interfaces, and that with proper parameterization it probably will enable future simulations of more complex systems involving biochemical species and irregular metallic surfaces.

ASSOCIATED CONTENT AUTHOR INFORMATION

Corresponding Author *E-mail: [email protected]

Notes The authors declare no competing financial interest. ACKNOWLEDGMENT The authors thank Dr. Sai Duan (KTH, Sweden) for many helpful and stimulating discussions. The Swedish National Infrastructure for Computing (SNIC) is acknowledged for providing


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 33

computational resources for the project “Multiphysics Modeling of Molecular Materials”, SNIC 2014/11-31. REFERENCES (1)

Qian, X.; Peng, X.-H.; Ansari, D. O.; Yin-Goen, Q.; Chen, G. Z.; Shin, D. M.; Yang, L.;

Young, A. N.; Wang, M. D.; Nie, S. In Vivo Tumor Targeting and Spectroscopic Detection with Surface-Enhanced Raman Nanoparticle Tags. Nature Biotech. 2008, 26, 83–90. (2)

Cao, Y. C.; Jin, R.; Mirkin, C. A. Nanoparticles with Raman Spectroscopic Fingerprints

for DNA and RNA Detection. Science 2002, 297, 1536–1540. (3)

Murphy, C. J.; Gole, A. M.; Stone, J. W.; Sisco, P. N.; Alkilany, A. M.; Goldsmith, E. C.;

Baxter, S. C. Gold Nanoparticles in Biology: Beyond Toxicity to Cellular Imaging. Acc. Chem. Res. 2008, 41, 1721–1730. (4)

Palmal, S.; Jana, N. R. Gold Nanoclusters with Enhanced Tunable Fluorescence as

Bioimaging Probes. WIREs Nanomed. Nanobiotechnol. 2014, 6, 102–110. (5)

Eustis, S.; El-Sayed, M. A. Why Gold Nanoparticles Are More Precious than Pretty

Gold: Noble Metal Surface Plasmon Resonance and Its Enhancement of the Radiative and Nonradiative Properties of Nanocrystals of Different Shapes. Chem. Soc. Rev. 2006, 35, 209– 217. (6)

Zhang, Y.; Cui, X.; Shi, F.; Deng, Y. Nano-Gold Catalysis in Fine Chemical Synthesis.

Chem. Rev. 2012, 112, 2467–2505.


22

Page 23 of 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


(7)

Cesca, T.; Calvelli, P.; Battaglin, G.; Mazzoldi, P.; Mattei, G. Local-Field Enhancement

Effect on the Nonlinear Optical Response of Gold-Silver Nanoplanets. Opt. Express 2012, 20, 4537–4547. (8)

Kim, S.; Jin, J.; Kim, Y.-J.; Park, I.-Y.; Kim, Y.; Kim, S.-W. High-Harmonic Generation

by Resonant Plasmon Field Enhancement. Nature 2008, 453, 757–760. (9)

Prodan, E.; Radloff, C.; Halas, N. J.; Nordlander, P. A Hybridization Model for the

Plasmon Response of Complex Nanostructures. Science 2003, 302, 419–422. (10)

Li, X.; Rinkevicius, Z.; Ågren, H. Two-Photon Absorption of Metal-Assisted

Chromophores. J. Chem. Theory Comput. 2014, 10, 5630–5639. (11)

Rinkevicius, Z.; Li, X.; Sandberg, J. A. R.; Agren, H. Non-Linear Optical Properties of

Molecules in Heterogeneous Environments: A Quadratic Density Functional/Molecular Mechanics Response Theory. Phys. Chem. Chem. Phys. 2014, 16, 8981–8989. (12)

Li, X.; Rinkevicius, Z.; Ågren, H. Electronic Circular Dichroism of Surface-Adsorbed

Molecules by Means of Quantum Mechanics Capacitance Molecular Mechanics. J. Phys. Chem. C 2014, 118, 5833–5840. (13)

Rinkevicius, Z.; Li, X.; Sandberg, J. A. R.; Mikkelsen, K. V.; Ågren, H. A Hybrid

Density Functional Theory/Molecular Mechanics Approach for Linear Response Properties in Heterogeneous Environments. J. Chem. Theory Comput. 2014, 10, 989–1003. (14)

Cicero, G.; Calzolari, A.; Corni, S.; Catellani, A. Anomalous Wetting Layer at the

Au(111) Surface. J. Phys. Chem. Lett. 2011, 2, 2582–2586.


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

(15)

Page 24 of 33

Iori, F.; Di Felice, R.; Molinari, E.; Corni, S. GolP: An Atomistic Force-Field to Describe

the Interaction of Proteins with Au(111) Surfaces in Water. J. Comput. Chem. 2009, 30, 1465– 1476. (16)

Hoefling, M.; Monti, S.; Corni, S.; Gottschalk, K. E. Interaction of β-Sheet Folds with a

Gold Surface. PLoS ONE 2011, 6, e20925. (17)

Heinz, H.; Vaia, R. A.; Farmer, B. L.; Naik, R. R. Accurate Simulation of Surfaces and

Interfaces of Face-Centered Cubic Metals Using 12–6 and 9–6 Lennard-Jones Potentials. J. Phys. Chem. C 2008, 112, 17281–17290. (18)

Wright, L. B.; Rodger, P. M.; Corni, S.; Walsh, T. R. GolP-CHARMM: First-Principles

Based Force Fields for the Interaction of Proteins with Au(111) and Au(100). J. Chem. Theory Comput. 2013, 9, 1616–1630. (19)

Järvi, T. T.; van Duin, A. C. T.; Nordlund, K.; Goddard, W. A. Development of

Interatomic ReaxFF Potentials for Au–S–C–H Systems. J. Phys. Chem. A 2011, 115, 10315– 10322. (20)

Siepmann, J. I.; Sprik, M. Influence of Surface Topology and Electrostatic Potential on

Water/Electrode Systems. J. Chem. Phys. 1995, 102, 511–524. (21)

Reed, S. K.; Lanning, O. J.; Madden, P. A. Electrochemical Interface between an Ionic

Liquid and a Model Metallic Electrode. J. Chem. Phys. 2007, 126, 084704.


24

Page 25 of 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


(22)

Vatamanu, J.; Borodin, O.; Smith, G. D. Molecular Insights into the Potential and

Temperature Dependences of the Differential Capacitance of a Room-Temperature Ionic Liquid at Graphite Electrodes. J. Am. Chem. Soc. 2010, 132, 14825–14833. (23)

Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular Dynamics

Simulations. J. Phys. Chem. 1991, 95, 3358–3363. (24)

Mortier, W. J.; Ghosh, S. K.; Shankar, S. Electronegativity-Equalization Method for the

Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315–4320. (25)

van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force

Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396–9409. (26)

Mayer, A. Formulation in Terms of Normalized Propagators of a Charge-Dipole Model

Enabling the Calculation of the Polarization Properties of Fullerenes and Carbon Nanotubes. Phys. Rev. B 2007, 75, 045407. (27)

Jensen, L. L.; Jensen, L. Electrostatic Interaction Model for the Calculation of the

Polarizability of Large Noble Metal Nanoclusters. J. Phys. Chem. C 2008, 112, 15697–15703. (28)

Morton, S. M.; Jensen, L. A Discrete Interaction Model/Quantum Mechanical Method for

Describing Response Properties of Molecules Adsorbed on Metal Nanoparticles. J. Chem. Phys.

2010, 133, 074103. (29)

van der Vorst, H. Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for

the Solution of Nonsymmetric Linear Systems. SIAM J. Sci. and Stat. Comput. 1992, 13, 631– 644.


25


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

(30)

Page 26 of 33

Fennell, C. J.; Gezelter, J. D. Is the Ewald Summation Still Necessary? Pairwise

Alternatives to the Accepted Standard for Long-Range Electrostatics. J. Chem. Phys. 2006, 124, 234104. (31)

Zahn, D.; Schilling, B.; Kast, S. M. Enhancement of the Wolf Damped Coulomb

Potential: Static, Dynamic, and Dielectric Properties of Liquid Water from Molecular Simulation. J. Phys. Chem. B 2002, 106, 10725–10732. (32)

Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. Exact Method for the Simulation

of Coulombic Systems by Spherically Truncated, Pairwise r–1 Summation. J. Chem. Phys. 1999, 110, 8254–8282. (33)

Cagin, T.; Kimura, Y.; Qi, Y.; Li, H.; Ikeda, H.; Johnsonb, W. L.; Goddard, W. A.

Calculation of Mechanical, Thermodynamic and Transport Properties of Metallic Glass Formers. Mater. Res. Soc. Symp. Proc. 1998, 554, 43–48. (34)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.

Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. (35)

McCann, B. W.; Acevedo, O. Pairwise Alternatives to Ewald Summation for Calculating

Long-Range Electrostatics in Ionic Liquids. J. Chem. Theory Comput. 2013, 9, 944–950. (36)

Andersen, H. C. Rattle: A "Velocity" Version of the Shake Algorithm for Molecular

Dynamics Calculations. J. Comput. Phys. 1983, 52, 24–34.


26

Page 27 of 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


Table of Contents Graphic and Synopsis


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

Fig. 1 178x76mm (300 x 300 DPI)


Page 28 of 33

Page 29 of 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


Fig. 2 168x75mm (300 x 300 DPI)



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

Fig. 3 183x76mm (300 x 300 DPI)


Page 30 of 33

Page 31 of 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


Fig. 4 177x142mm (300 x 300 DPI)



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

Fig. 5 242x106mm (300 x 300 DPI)


Page 32 of 33

Page 33 of 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


TOC 34x14mm (300 x 300 DPI)


Molecular Dynamics Simulations Using a CapacitanceâPolarizability

Recommend Documents

Molecular Dynamics Simulations Using a CapacitanceâPolarizability