POLARIZABLE FORCE FIELD FOR MOLECULAR DYNAMICS SIMULATIONS OF SILVER NANOPARTICLES

Contact of silver metal surfaces with water, ions and organic ligands experiences induced charges, leading to attractive polarization. These forces play an important role at inorganic/organic interfaces and complement other non-bonded surface interactions. Despite the importance of these interactions, it, however, remains difficult to implement polarization effects to classical molecular dynamics (MD) simulations. In this contribution, we first present an overview of two popular polarizable models, such as Drude oscillator and the rigid rod model, which are utilized to mimic the polarizability of bulk metals. Second, we implemented the rigid rod model to the polarizable force field (FF) for a silver atom, which was further adapted for atomistic MD simulations of silver nanoparticles (AgNPs) composed of 1397 atoms. In our model, induced charge polarization is represented by the displacement of a charge-carrying virtual site attached rigidly to an original Ag atom. To explore the role of polarization, we compared the performance of the classical nonpolarizable FF and the new polarizable model in the MD simulations of adsorption of water and ions onto quasispherical AgNP and the flat crystalline silver surface. The analysis of the radial distribution function of Ag-Ag atoms demonstrated that the introduction of the polarization effect had minor effects on face-centered cubic (fcc) packing of silver atoms of bare and water-solvated AgNPs. We found that the polarizable FF causes some increase in attractive interactions between the silver surface and water molecules and Na ions. As a crucial test of the developed polarizable model, the structure of adsorbed interfacial water molecules was analyzed. Our data suggest that the environment-induced polarization of the silver surface contributes significantly to the structure of adsorbed interfacial water layers and it also plays an important role in the adsorption of positive ions. However, it was also found out that the polarization effect has a rather short-range effect, so that a minor contribution of silver polarization was seen for adsorption of water molecules and ions from distant solvation shells.

Silver metal surfaces in contact with the aqueous environment experience attractive polarization owing to metal-induced charges [31]. This fundamental physical interaction complements Van der Waals surface interactions; however, despite recent progress in this field, it remains difficult to implement polarization effects to classical MD simulations [32]. A much less considered contribution to the binding strength is a polarization of the metal surface through the typically polar and often ionic environment created by solvents, surfactants, peptides or DNA. Therefore, the development of advanced computational models for complex hybrid organic/inorganic nanomaterials is required [33][34][35][36]. https://doi.org/10.26565/2220-637X-2019-32-03 The main goal of our study is the implementation of the polarizable Ag model to the existing nonpolarizable FF, which were recently used for MD simulations of silver nanoparticles and interfaces [8,19,24,[37][38]. Besides, the influence of polarization effects on the adsorption of water molecules and ions at the silver/water interface is considered.

Review of Polarizable Force Fields
Metal surfaces acquire a charge in aqueous systems and these attracted surface ions form a layer balancing the oppositely charged counter ions. This structure is known as the electrochemical double layer (EDL). EDLs profoundly affect the physical and chemical behavior of such heterogeneous systems [39]. Unless there are a lot of common techniques, which involve the measurement of electrokinetic phenomena caused by the tangential movement of an electrolyte fluid relative to the charged surface, the molecular dynamics methods of charged NPs surfaces simulations are only arising. The investigations, which have already found their use in studying these systems, are summarized below.
Iori and Corni have recently introduced novel polarizable MD models for a bulk gold and gold surface (001) and (111) [40]. Their approach is based on the Drude oscillator model [32] in which a virtual interaction site is added to the real metal atom being connected by a harmonic bond with a force constant k and the null rest length l 0 =0 ( Fig. 1 left). To keep the system neutral, the authors offered to assign negative charge -q on the real metal atom and positive charge +q placed on the virtual site. The electrostatic interaction between these two charges on the same polarizable metal atom is not included in the energy and force calculations during the simulations. Additionally, a virtual site has no Van der Waals interactions with other atoms. Thus, the Drude oscillator model for metal requires three parameters: a force constant k, point charge q and a mass m, which help to mimic the behavior of the polarizable metal atom. Their combinations define the quantities that determine the behavior of the polarizable metal atom: the polarizability α (1) and the harmonic oscillator frequency ω (2) [40], respectively.
However, some disadvantage of this model is that the fictitious moving charge can collapse on other charges, which can lead to the system's crash. Attractive polarization occurs in response to exterior charges in the vicinity of a metal surface and it scales with the magnitude of atomic charges and with the length of multipoles, such as the distance between a real metal atom and a virtual site. Further improvement of the Drude oscillator model was given by the rigid rod model (Fig. 1 right) [32,40]. In the rod model, the rod always has a finite dipole moment and there is a fixed distance l 0 between the real gold atom and a virtual interaction site. Thus, the virtual site is forced to be at a given distance l 0 from the original metal atom, however, its orientation is still free. Similarly to the Drude oscillator model, the virtual site has no Van der Waals interactions, while it bears a positive charge q (−q is assigned to the original metal atom) and it has an atomic mass m. The rod model also requires three parameters: distance l 0 , partial charge q, a mass m. Unlike the Drude oscillator model, the rod model always has a finite dipole moment μ (3). To reproduce the physically meaningful polarizability of bulk metals, this model requires an orientational averaging (i.e., T > 0), so that the relation among the parameters and the characterizing physical quantities involve the temperature: the polarizability is also given by α as defined by (1) [32,40]. The typical values of these parameters for gold are T=300 K, q=0.3, l 0 =0.7 Å and m=0.5 a.m.u. [40]. This choice was confirmed by the necessity of making stable MD simulations and the fact that the dipole moment of every single rod should average to zero. It was found that q=0.3 gave the sufficient polarizability while not causing the system to freeze [40]. The MD simulations based on the Lennard-Jones 6-12 potential for gold and performed in an NVT ensemble, [40] demonstrated that the rod model is capable to reproduce the interaction energy between a charge and neutral metal slab.
( 4 ) The polarizable model for gold, suitable for biomolecular MD simulations, has recently been introduced in the GolP [41] and GolP-CHARMM [42][43] force fields, respectively. Recently, this approach has further been extended and implemented in new FF, AgP-CHARMM [44][45][46], suitable for simulations of the aqueous Ag(111) and Ag(100) interfaces, where abbreviations (111) and (100) refer to the Miller indices form a notation system in crystallography for planes in crystal lattices. AgP-CHARMM utilizes the rigid rod model to mimic the polarizability of the Ag atoms, where the total average dipole moment μ, is determined by (4), where T=300 K, q=0.308, l 0 =0.7 Å. Moreover, to ensure atop-site adsorption, both GolP-CHARMM and AgP-CHARMM introduce additional virtual interaction sites, rigidly placed between real metal atoms in crystal slabs, which allow overcoming some problems with the artificial orientation of hydrogen atoms of water molecules. Most recent extension of polarizable FF has been for the aqueous Pd(111) interface, referred to as PdP-CHARMM [47]. Unfortunately, the applicability of all these force fields is limited to specific (111) and (100) surfaces. Therefore, the tedious re-parametrization procedure is required before introducing new parameters for other crystal facets [48]. Therefore, MD simulations of spherical nanoparticles are not straightforward with this force field.

Molecular dynamics simulation setup
The preformed silver nanoparticles (AgNPs) were approximated by the perfect face-centered cubic (fcc) crystalline structure (Fig. 2). Two force-field (FF) models for AgNP were considered: (i) a nonpolarizable model for Ag, consisting of neutral atoms with the zero charge, (ii) a polarizable model for Ag, in which metal polarization was implemented by using the rigid rod approach. In both models, the repulsion and dispersion terms of nonbonded interactions between silver atoms were computed by using the Lennard-Jones 12−6 potential energy function (5), which describes the dependence of the potential interaction energy V LJ (r ij ) of two silver atoms as a function of the interatomic distance. The nonbonded interaction parameters Ag-Ag were based on non-polarizable INTERFACE-FF [31,36,49] and were validated in our recent works [8,19,24,37]. For MD simulations in vacuum, no of any rigid bonds and restraints were applied between silver atoms, so that the silver core crystalline structure was maintained by the Ag-Ag nonbonded LJ interactions.
Electrostatic interactions between two charged particles were calculated by the Coulomb potential (6). Therefore, during MD simulations the electrostatic interactions among the polarizable Ag atom dipoles were treated by Coulomb term allowing collective and cooperative motions of dipoles in response to external charges.
The description of the polarizable model for AgNP is represented as follows: Ag atom has q=-0.308, whereas the virtual site has q=+0.308 to make the model neutral. An alternative model with positive and negative charges altered between each other was also considered. In the original rigid rod model, the mass of the virtual site has been assigned to zero [40]. In our polarizable model, we assigned the minimal nonzero mass of 1.0, as required in the GROMACS package. The rigid rod was approximated with the rigid chemical bond of l 0 =0.07 nm and k b =1×10 12 kJ/molꞏnm 2 , respectively. The initial values of q and l 0 were taken from the original study [44]. The nonbonded interaction parameters Ag-Ag σ=0.2964 nm and ε=19.05865 kJ/mol were taken from our recent works [9,24]. However, both parameters σ and ε were assigned to zero for the virtual site. For MD simulations of polarizable AgNPs in water, the model was adjusted by adding the LJ interactions between Ag atoms and hydrogen atoms (HW) of SPC waters. The correction LJ interaction parameters Ag-HW were taken [50] as following: σ=0.30584 nm and ε=0.74323 kJ/mol, respectively.
AgNP was simulated according to the following procedure: (i) Steepest descent energy minimization was performed for 500 steps. (ii) Initial orientation of dipoles and atomic velocities were generated with the Maxwellian distribution at T=250 K. Productive MD sampling was carried out at the reference temperature of T=303 K. The water molecules are described by the simple point charge (SPC) model [51]. The parameters for Na + and Clions were used as implemented in the GROMOS 45a3 FF [52]. The reference temperature of T=303 K, which was kept constant using the Berendsen weak coupling scheme with the temperature coupling constant of  T =0.1 ps [53]. The cutoff distance of 0.8 nm was used for Lennard-Jones interactions. The MD simulation time step was 2 fs with the neighbor list updates every 10 fs. The MD simulations were carried out using the GROMACS set of programs, version 4.6.5 [54]. Molecular graphics and visualization were performed using VMD 1.9.2 [55].

Results and discussion
Isolated Silver Nanoparticles MD simulations of silver nanoparticles of 1397 atoms were carried out by the using of polarizable and nonpolarizable FFs in vacuum. To investigate the accuracy of the nanoparticles' structure ( Fig. 2a) and to compare their behavior, we used the dependences shown below.
Radial Distribution Function. To characterize the structure, atom packing and the long-range order of Ag atoms in AgNPs, the pair radial distribution function (RDF) g(r) was calculated between all silver atoms. As can be seen in Fig. 2b, the RDF plots of Ag-Ag interactions of both models are characterized by narrow and sharp peaks and have the same peak positions. The sharp RDF peaks are indicative of small-amplitude vibrations of Ag atoms about the lattice site position and very low diffusive movements over the local region of the crystal. These results demonstrate that the polarizable model represents well quasi-spherical AgNPs with face-centred cubic fcc crystalline structure [56][57].
Root Mean Square Displacements. For further characterization of the dynamic stability of AgNPs, their structures were analyzed by calculating the root mean square displacement (RMSD) of Ag atoms as a function of time. The RMSD of all Ag atoms was calculated with respect to their position in a perfect fcc structure by using the GROMACS utility g_rms, which carries out least-square fitting the AgNP structure (τ 2 ) to its initial perfect structure (τ 1 = 0) (7).

   
where N and r i (τ) are the number of Ag atoms, and the position Ag atom i and its reference position at time τ [54,58].

Silver Nanoparticles in Aqueous Solution
On the next step, the performance of the polarizable FF was tested on AgNPs in the aqueous environment. For this purpose, AgNP 1397 was placed in the cubic box of 6.5 nm size and solvated by 7880 SPC water molecules. Moreover, to study the effect of the polarization on adsorption of positive and negative ions on AgNP the additional MD simulations were carried out in the presence of 2 M of NaCl. Fig. 3a shows the MD snapshot of the solvated AgNP 1397 in the presence of NaCl.
Adsorption of water molecules onto AgNP The adsorption of water molecules on AgNP was analyzed by calculating of the RDF between Ag and water oxygen atom OW, respectively. The comparison of the RDF plots of g(r) Ag-OW calculated for AgNP 1397 modelled with the nonpolarizable and polarizable MD models shows very similar results (Fig. 3b). These findings suggest that the metal polarization does not perturb significantly the strength of the water adsorption.
Adsorption of ions onto AgNP. The further comparison of the nonpolarizable and polarizable FFs was performed by the analyzing of adsorption of ions onto the silver surface. Fig. 3c shows the time evolution of the number of contacts between AgNP and Na + and Clions, respectively. As can be seen, the striking difference was observed for the adsorption of positive Na + ions as compared to that of the negative Clones. The use of the polarizable FF significantly favors the adsorptions of the Na + ions onto AgNP. Such influence of the salt is consistent with by the formation of the electrochemical double layer (EDL) between polarizable surface of silver atoms and counter ions.

Structure of AgNP in solution.
Finally, the performance of the new polarizable FF model for reproducing of the fcc structure of AgNP 1397 was summarized by a series of MD simulations in vacuum and in the aqueous solution in the absence and in the presence of NaCl as compared in Fig. 3d by the corresponding RMSD time traces. It is seen that, upon MD simulations of AgNP 1397 by the using the polarizable FF model, the RMSD values of silver atom positions were increased from 0.0189±0.0007 nm in vacuum up to 0.0215±0.0009 (curve 3) nm in bulk water. However, only a small additional increase of up to 0.0218±0.0009 (curve 4) was observed in the presence of 2 M NaCl.

Adsorption Configurations of Water Molecules on Silver Surface
The structure and properties of water molecules, adsorbed on well-defined metal surfaces, have been the subject of numerous experimental and theoretical investigations [59][60][61][62][63]. It these studies, the most favorable adsorption configurations of the water molecule on the silver surface were identified as shown in Fig. 4. The most stable structures were found to be atop configurations A and B, in which water is coordinated through the oxygen atom to the silver atom. The coordination through the hydrogen atom in configuration C was found to be unfavorable, so that this structure tends to rearrange to the above mentioned configurations A and B, respectively.
The adsorption configuration of water molecules on the silver surface can therefore be considered as a crucial test and validation of the polarizable Ag model for reproducing of the well-known interaction characteristics of water molecules onto the silver surface. For simplicity, the silver surface was approximated with the three fixed (111) layers, for which the two sets, nonpolarizable and polarizable Ag models were applied. The adsorption dynamics and preferred structures of water molecules onto the silver slab were simulated at the NVT ensemble as shown in Fig. 5a.  Figure 4. Schematic side view of the stable structures found for the water molecule on the silver surface adapted from [61].
Surprisingly, the analysis of the preferred adsorption configurations of water molecules on the silver surface modelled with the polarizable silver model revealed the large fraction of water molecules bound to the surface through atop configuration C (Insert in Fig. 5a). In contrast, such binding mode of water molecules was rarely observed for MD simulation with the nonpolarizable silver model. These findings were further confirmed by the RDF analysis of g(r) Ag-HW calculated between silver atoms and hydrogen atoms of water molecules (Fig. 5b), which revealed a small peak at 0.18 nm on the RDF curve (2) derived from the polarizable FF simulations. Figure 5. (a) Side-view snapshots of MD simulations of water molecules adsorbed on the silver surface as estimated with different FFs. Inset shows that the polarizable FF revealed that some fraction of adsorbed water molecules were found to be misoriented with one of its HW atoms pointing toward the silver surface. (b-c) Comparison of the RDF plots of g(r) Ag-HW (b) and g(r) Ag-OW (c) calculated between the silver atom and either hydrogen atoms (HW) or oxygen atoms (OW) of water molecules calculated with the nonpolarizable FF (1), polarizable FF (2), and polarizable FF with LJ-corrections on HW atoms.
Artificial orientation of hydrogen atoms of water molecules adsorbed onto the silver could be emerged as due to some overestimation of Coulomb attractions. This question was already discussed in literature and the following solution was suggested [42,44]. To shift water hydrogens from the surface and make the computational model more physical, there is a necessity to compensate the Coulomb attraction by the Lennard-Jones nonbonded repulsion interaction between HW and silver atoms [50]. Therefore, the polarizable silver model was further modified by adding the LJ interactions between Ag atoms and hydrogen atoms (HW) of SPC waters. This modified polarizable model is referred to as Polarizable Ag+LJ-for-H. The LJ interaction parameters Ag-HW were as following: σ=0.30584 nm and ε=0.74323 kJ/mol, respectively. As shown on the Fig. 5a-c the water adsorption modelled with the modified polarizable silver model becomes very similar to physical expectations. The RDF curve (3) of g(r) Ag-HW shows no artificial peaks at 0.18 nm observed before.
The comparison of the performance of the three different Ag models (Fig. 5) revealed that the implementation of the Ag polarization resulted in the limited influence on the adsorption behavior of water molecules near the Ag interface. The small changes in RDF peaks of g(r) Ag-OW in Fig. 5c demonstrated that the polarization effect of Ag has rather short-range effects, so that a minor contribution of the silver atom polarization was seen for adsorption of water molecules and ions form distant solvation shells. Recently, Menziani and coauthors [64] used MD simulations to study the adsorption of cytochrome c on various solvated silver surfaces by using the nonpolarizable and polarizable AgP-CHARMM model [44] and they have also reported small effects of the implemented Ag polarization on binding modes of cytochrome c. Taking into account the importance of the proper characterization of water molecule orientation and dynamics at the aqueous/metal interfaces, more sophisticated approaches have also been suggested to mimic polarizability of bulk metal and surfaces [50,65]. Towards the proper treatment of the polarizability of metallic gold, Li and Ågren have developed the capacitance-polarizability force field for noble metal surfaces and nanoparticles by employing of different interacting potential functions for the metallic and nonmetallic components [66]. In their force field, the water−gold interaction was divided into three categories, each of which is described by different potential functions: the TIP3P model for water, the quantum Sutton−Chen potential for inner gold atoms, and the capacitance−polarizability interaction model for outermost gold atoms, represented by a modified Lennard-Jones potential [66]. While this FF introduces extra complexity and hampers the practical flexibility and transferability for applications in different biomolecular systems, it provides some evidence that the surface of the aqueous gold nanoparticle may carry a negative charge, which is balanced by the positive charge in the second outermost gold layer.
Following the above arguments, we also verified the performance of the rigid rod polarizable model for AgNP by altering its positive and negative charges. The original rigid rod model was parameterized with a negative charge -q assigned at metal atom Me, whereas a positive charge +q fixed at a virtual site (Fig. 1). However, physical meaning of atomic polarizability is due to the displacement of electronic density under the influence of the external field, so that an alternative charging scheme of the rod model can also be valid, in which negative charge -q is placed at a virtual site and positive charge +q is fixed at metal atom Me, as shown in Fig. 6a. We carried out control MD simulations by using such alternative polarizable scheme. We found, however, that the charge reversing in the rigid rod model has a little effect on water adsorption at the silver surface. In particular, the structure of the first water layer and artificial water ordering were observed with MD simulations by using the original and reversed charging scheme of the rigid rod model, as seen in a snapshot and RDF plots of g(r) Ag-HW in Fig. 6b-c. Moreover, we also found that the charge reversing in the rigid rod model did not affect the adsorption of Na + and Clions (not shown), which suggests that due to the short-range polarizability effects of the silver surface these ions prefer residing in bulk solution.

Conclusions and Future Perspectives
The polarization forces at inorganic/organic interface play an important role in the adsorption of water molecules, ions and organic ligands [31]. While the most of structural and thermodynamic properties of adsorbed species at silver surfaces are well accounted with the existing empirical force field, however, the accuracy of the simulations still depends on the suitability of the force field utilizing the pair-wise Lennard-Jones interaction models. Therefore, it has been shown that some certain aspects of the structure and dynamics of the first adsorbed water layer require the use of the softer Morse and Buckingham potentials [50]. Alternatively, it has been suggested that the existing force field may be upgraded by implementing various approaches, mimicking metal polarization effects [32].
In this work, we implemented and validated the polarizable FF parameters for Ag suitable for atomistic MD simulations of silver nanoparticles (AgNPs) and bulk-Ag/water interfaces. Our model is derived from the rigid rod model [40], in which induced charge polarization is represented by the displacement of a charge-carrying virtual site, attached rigidly to a polarizable Ag atom. We compared the performance of the classical nonpolarizable FF and the new polarizable model of Ag by the modeling of adsorption of water molecules and NaCl ions onto quasi-spherical AgNP and the flat crystalline (111) silver surface. The analysis of the radial distribution function of Ag-Ag atoms demonstrated that the introduction of the polarization effect had minor effects on the fcc packing of silver atoms in both isolated and water solvated AgNPs, respectively. The implementation of the polarization effect led to some increase in attractive interactions between the silver surface and water molecules and Na + ions. As a crucial test of the developed polarizable model, the structure of adsorbed interfacial water molecules was analyzed. We found that the environment-induced polarization of the silver surface contributes significantly to the structure of the adsorbed interfacial water layer and plays an important role in the adsorption of positive ions. However, it was also demonstrated that the polarization effect of Ag is rather short-range, so that a minor contribution of the silver atom polarization was seen for adsorption of water molecules and ions from distant solvation shells [64]. Despite the short-range polarization effects, the proposed polarizable FF model for Ag can, however, have an important contribution for MD simulations of adsorption of charged organic ligands and polyelectrolytes due to Coulomb attractions. Moreover, the use of the polarizable Ag model might be promising for accurate estimation of the electrostatic potential at the silver surface. These studies are currently in progress in our team.