Interplay between Polarizability and Hydrogen Bond Network of Water

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

A: Molecular Structure, Quantum Chemistry, and General Theory

The Interplay Between Polarizability and Hydrogen Bond Network of Water: Reparametrizing the Flexible Single Point-Charge Water Model by the Nonlinear Adaptive Force Matching Approach I-Shou Huang, and Ming-Kang Tsai J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b12726 • Publication Date (Web): 25 Apr 2018 Downloaded from http://pubs.acs.org on April 26, 2018

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 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 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.

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

The Journal of Physical Chemistry

The Interplay between Polarizability and Hydrogen Bond Network of Water: Reparametrizing the Flexible Single Point-Charge Water Model by the Nonlinear Adaptive Force Matching Approach I-Shou Huang and Ming-Kang Tsai* Department of Chemistry National Taiwan Normal University, Taipei 11677, Taiwan Email: [email protected] TEL: 886-2-77346217

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 31

Abstract

An adaptive force matching (AFM) scheme using the nonlinear optimization to re-parameterize the three-site, flexible and polarizable single-point-charge (SPC) water model is reported. We compare the radial distribution functions of the intermolecular oxygen-oxygen, oxygen-hydrogen and hydrogenhydrogen distances with the recent scattering experiments, the previous AFM-fitting water model – MP2f and the atomic multipole expanded AMOEBA model. Our non-polarizable SPC-3f(0) model captures the feature of the first solvation shell of bulk water. With the ad hoc inclusion of the isotropic polarizability, the polarizable SPC-3f(0.6) water model recovers the many-body effect of the second solvation shell. In the n-body decomposition analysis, SPC-3f(0) model predicts the best agreement with MP2/aug-cc-pVTZ calculations with using the low-dimensional (H2O)4-ring and (H2O)6-ring clusters. For the comparison using three-dimensional (H2O)6-prism and (H2O)16-4444a clusters, SPC-3f(0.6) predicts the consistent results as those of AMOEBA and MP2 levels. For simulating the water-clusterdominant system like supercritical water, SPC-3f(0) well characterizes the combination mode of bending and stretching at 5300 cm-1.


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


Introduction Water, the universal and life-vital solvent, plays a critical role in numerous chemical and biological reactions in nature. Representing such polar medium through the classical viewpoint of intermolecular interactions has been one of the core developments in the history of molecular mechanics. The representations of multiple interaction sites, ranging from one to six sites, have all been introduced in the parameterization of classical water models.1-2 The one-site soft sticky dipole (SSD) water model reported by Tan et al. contains only one interaction site at the centre of mass with a tetrahedral-coordinated sticky potential to restrict the tetrahedral coordination of neighbouring molecules.3 A two-site model was introduced by Dyer et al. using site-renormalized molecular fluid theory in the simulation of dielectric and solvation properties of bulk water.4 The models using three interaction sites centred at the atoms of water molecules are commonly seen in the molecular dynamics (MD) simulations, and these 3-site models can be further classified as the rigid5-7 or flexible8-13 cases. Using more than three sites representing water molecules has also been the focus of water model development, namely TIP4P,6 TIP5P,14 AMOEBA,15 TTM16-18 and others.19-23 The recent development in simulating the inter-water interactions in terms of the explicit many-body potential energy functions, being summarized by the latest review on water models,2 has ultimately converged the accuracy of water models as being demonstrated by MB-pol potential.2, 24-27 However, the computational cost increases rapidly as the number of interaction sites used to represent the water monomer, and the number of water molecules in the simulations, particularly for the interfacial-heterogeneous simulations. Maintaining an ideal balance between the computational cost and scientific accuracy is a subject of great research interest.



Page 4 of 31

Reproducing the information produced by ab initio calculations by classical force fields has been one of the popular means in parameterizing such classical representations. Force matching (FM) method reported by Ercolessi and Adames is one of the pioneer studies for demonstrating such QM-to-MM mapping process for aluminium potential.28 Izvekov et al. used the spline representation for the short-range interatomic forces of the water model in order to maintain the linearity of the objective function, followed by large configuration sampling to reach an overdetermined system, and solved such a system by the singular value decomposition algorithm.29 The required ingredients such as adapting spline expression and collecting large configuration sampling could present mathematical challenges for further generalizing the application of FM. Akin-Ojo et al. introduced the adaptive force match (AFM) approach with using an efficient QM/MM sampling scheme and the combinatorial parameter set to well generate the parameters of a simple point-charge (SPC) water model at a remarkably low computational cost.30-32 While there are many details in the mechanisms and history of FM, the most concerned part for the purpose of this study is the fitting process found in the core of AFM, said process, to our knowledge, has always been chosen to be suited for statistically linear equation systems, like Singular Value Decomposition (SVD) or QR Decomposition. The central problem here is that most force fields models do not have statistical linearity. Therefore, we employ a nonlinear fitting algorithm in the current AFM study to re-parameterize the SPC water model, one of the most popular and computational cost-effective water model, and attempt to create a more legitimate methodology and to broaden its applicability.

Methodology (a). The Adaptive Force-Matching method


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


AFM is an automated iterative process of five steps aiming to render a force field parameter set derived from ab initio level calculations as shown in Fig. 1. Every six steps are grouped as a “generation”, and every generation observes its own associated force field while providing the next generation a new, ideally improved, force field. Upon minimizing the difference between force field and ab initio levels, the circular process will be halted if the convergence criteria is reached. In this study, we aim to simulate bulk water with a flexible single point charge model via AFM.

Figure 1. Schematic representation of adaptive force-matching method (AFM). With this work, the intermolecular interactions between water molecules are described as the following equation:

∑ ∑


Eq. (1)


Page 6 of 31

where kr is the stretching force constant, re is the equilibrium bond length, kθ is the bending force constant, θe is the equilibrium bond angle, A and C are respectively the repulsive and attractive coefficients for the Lennard-Jones potentials (LJ) on oxygen atoms. qm and qn denote the atomcentred charges of hydrogen or oxygen atoms (2*qH + qO = 0 is constrained in this study). ri, θj, rk and rmn are the arbitrary bond length, bond angle, oxygen-oxygen distance and intermolecular charge-charge distance, respectively.

(b). Force field molecular dynamics simulations Following each iterative prediction of force field parameter set, we conduct a molecular dynamics (MD) simulation of (H2O)216 under canonical ensemble using the LAMMPS Molecular Dynamics Simulator, and the cubic supercell is fixed at 18.6445 Å3.33-34 All intermolecular Lennard Jones and columbic interactions are calculated with the cut-off at 9.0 Å, and Ewald Summation is utilized to describe the long range Columbic interactions. Nosé-Hoover thermostat (τ = 0.2) is adopted to control the temperature at 300K.35 Each simulation lasts 40 simulated picoseconds with a time-step of 0.5 femtosecond, structures are saved every femtosecond. The standard Velocity-Verlet integrator is chosen for the MD simulations.

(c). Snapshot Sampling We choose the very last structural frame of the 40-ps MD trajectory as the basis for all of the samples of that generation, with 96 different QM regions chosen from it. Our QM region contains 7 water molecules: one water molecule is randomly selected as the center of QM region along with its 6 closest neighbors. Each QM region is surrounded by 209 water molecules described at MM level.


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


(d). QM/MM Force Calculations The collected geometries of (H2O)216 will be calculated by ONIOM method and the forces applied to the QM region atoms are recorded. Because the calculation of each snapshot is not dependent on each other, this step can be perfectly parallelized. In this work, Gaussian 09 is used for ONIOM force calculations,36-37 with the QM part calculated at MP2/aug-cc-pVDZ level under the embedded MM charge, and the MM part using the force field of current generation.

(e). Data Weighting The QM atoms deemed to be more or less representative of the system will be weighted accordingly. We utilize the weighing rule reported by Akin-Ojo et al.30 where a water molecule’s weight is subject to the number of neighbors in its first solvation shell described by the QM region. If there is three or four such neighbour included in the QM region, then the relevant water molecule will be weighted as 100%. With only two neighbours, this molecule is halfdiscounted in weighting. With one or less neighbours, this molecule is completely ignored. The weighting for the atoms in a molecule automatically inherits its molecular weighting.

(f). Force Matching Finally, the collected and weighted QM forces, including the forces of the previous three generations, will be used as the QM reference in the objective function in order to reparameterize the new generation of force field via an optimization algorithms. The objective function (Obj) as listed below is used to describe the consistency between QM and MM levels, containing three terms, i.e. atomic forces (f), molecular forces (F) and molecular torques (τ).



∑!

$% ! "#!

&

#!%% " )∑! ∑'

$% & ! "#! " $% ' "*'

∑'

Page 8 of 31

$% ' "('

&

*'%% " )∑'

&

('%% " )∑'

$% & ' "*' "

$% & " ' "('

Eq. (2)

The superscripts QM and MM denote the theory level used to calculate the forces or torques. and

'

+

are the weights assigned on the atoms and molecules, respectively. ,+ and ,' are the

number of atoms and molecules considered. It should be noted that our objective function aims to maintain the equal balance between the atomic and molecular descriptions, however, the energetics contributed by each interaction term of Eq. (1) to the summation of atomic and molecular forces, and torques could be swappable or substitutable by each other. In this work we employ a non-linear optimization algorithm, L-BFGS-B, provided by the Python modules SciPy and NumPy to minimize Obj.

Results and Discussions (a). The Generated Force Fields Six batches of the AFM generated force fields are tabulated in the Tables S1-S6, each generation labelled as Bn(X) where n is the batch number and X denotes the number of the generation. The initial guess – Bn(0), n = 1-5, is the same initial guess used by Akin-Ojo et al.31 In the 10-generation history of B1-B4, all of the parameters seem to be monotonically converge while C of the LJ interaction reduces significantly (close to zero) starting from 625.502 kcal· mol-1· Å6. The reduction of C to negative value for LJ interaction was also observed by Nicolini et al. using rigid water models.38 The history for the parameters (A and C) of the Lennard-Jones interaction for all batches is summarized in Fig. 2. 8 ACS Paragon Plus Environment

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


For the B5 and B6, we restrict the change of each parameter to be less than 20% from its previous values at each generation during the courses of AFM. The 40-generation history of B5 and B6 are carried out where B5 starts from the same initial guess as B1-4 and B6 starts from the final results of Akin-Ojo et al.31 After the 10th generation, the parameters of B5 and B6 are mostly settled, though there are minor fluctuations among the later generations of parameters. Specifically, the bond length and angle of water defined by the force fields appears to be stable, but the bond strength, angular strength, Lennard-Jones parameters, and partial charges all mildly oscillate. We interpret this as the molecular stiffness along with the partial charges of the simple point charge model, attempting to account for what the Lennard-Jones equation is already trying to depict: the van der Waals force. Thus, these values wax and wane slightly across the generations. We can see that the initially guessed parameters are noticeably away from the rest of the trajectory. After about 10 generations in B5 and B6, it finally slows down and started to wander around (A,C) = (650,150). We believe this is evidence suggesting that convergence is achieved, and values of parameters around that region are all valid, positive results. This work picked the tenth generation parameter set (which first reached this region) to prevent over-fitting. Observe further, that nearing the end of this 40 generation of force matching, the final 3 points started to stray away from the valid region, we posit that this is the result of over-fitting finally manifesting. Such over-fitting could be caused the non-equipartitioning kinetic energy in the collected samples.



Page 10 of 31

Figure 2. The AFM history for the parameters of Lennard-Jones potential. (b). Radial Distribution Function To further investigate the physical implications of the generated parameters sets, molecular dynamic simulations of (H2O)1000 were used to generate water box trajectories with 1ns pre-equilibrium run before 1 ns production run. These trajectories were then analyzed for their radial distribution function, which contains valuable information regarding the structure of molecular configurations. Four generations of parameter sets have been chosen for this investigation, namely the initial guess, the first, fifth, and tenth generation parameter sets, each labelled B5(0), B5(1), B5(5), B5(10), respectively. For comparison purposes, the radial distribution function from the optimized result of Akin-Ojo et al.,30 and the radiation total scattering experiments39 are also included in the comparison. We notice that all of the RDFs are significantly improved after the first iteration of AFM. After which the difference among generations started to dwindle, by the tenth generation, the RDF is mostly converged. Compared to the previous results labelled in dotted-pink curve, our water box is generally less over-structured. We observe good matches for the first peak maximum of gOO and gOH as well as the second peak maximum of gHH and gOH. Despite the 10 ACS Paragon Plus Environment

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


equilibrium distance of OO and OH are well captured by the force field, the broadness of the first peaks of gOO and gOH are under-represented, which could be attributed to the harmonic stretching and bending described in the functional of force field. The first peak of gHH and the second peak of gOO are substantially under-estimated where the nature of many-body effect plays an important role to these structural features. We re-confirm the current non-linear approach could not redistribute the many-body effect and anharmonicity, through the iterative matching between ab initio and fore field calculations, to the pair-wise electrostatic and Lennard Jones interactions.



Page 12 of 31

Figure 3. The radial distribution functions for the intermolecular OO, HH and OH using B5(0), B5(1), B5(5) and B5(10) force fields are labelled from bright blue to dark blue, respectively. The pink-dotted curve denotes the optimized parameter of Akin-Ojo et al.31 The black and grey curves represent Soper’s and Skinner’s measurement, respectively.39-41 (c). Adding Polarization In order to take into account the many-body effect, we ad hoc include the polarization effect through the thermalized Drude oscillator model in LAMMPS. A similar water model was previously introduced by Lamoureux et al. as known as the SWM4-NDP model42-43 where a water molecule is described by 4-site interactions and oxygen-based Drude oscillator. In this study, we include polarization effect on top of B5(10) parameter set containing only 3 sites, in which the polarization represented by point charge of Thole’s dipole screening scheme44 with the Thole damping parameter at 2.6 (unitless). The polarizability (α, unit in Å3) determines the charge separation of the point charges. We increase α of water from 0.0 to 0.6 at 0.1 incremental, and simulate again the RDFs of (H2O)1000 at room temperature and density at 0.997 g/cm3.


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

Figure 4. The radial distribution functions for the intermolecular OO, HH and OH using SPC3f(α) with the pink-dotted curve denotes the optimized parameter of Akin-Ojo et al.31 The black and grey curves represent Soper’s and Skinner’s measurements, respectively.39-41 As the many-body effect is gradually included, the second peak of gOO is also gradually reproduced as α increases to 0.6 as shown in the gOO of Fig. 4. However, the recovery of the second peak of gOO also causes the over-structuring of the first peak. This is the same behaviour as being described by Akin-Ojo’s early simulations. The broadness of the first peaks of gOO and gOH are partially recovered, but the noticeable underestimation of width by polarizable B5(10) force field is in present. Further improvement could use anharmonic functions like Morse potential to recapture such broadness. The average pressures of the B5(10) constant-density simulations using α = 0 and 0.6 are estimated to be 7333(970) and 3802(832) atm, respectively, with the corresponding standard deviation listed in the parentheses. For the rest of discussions, we would denote the polarizable B5(10) force field as SPC-3f(α) – a flexible, polarizable and 3site simple point-charge model where α is the choice of polarizability.


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


(d). Analysing the Many-Body Effect The n-body decomposition procedure45 is adopted to elaborate the many-body effect using the finite size cluster models. The total binding energy of a given cluster can be expressed as the summation of the nth-body energy as 1

BE / 0! !23

where 0! is the nth-body energy and N is the size of the calculated fragment of the cluster. In this study, we calculate the first six terms using various water cluster models. The first three terms can be expressed as the following, and the last three terms are elaborated in the supporting materials: 1

03 /40 5 0678 59 '23

1=3

1

0& / / 40 5, 0 5 0 9 '23 ;2' / /

1

/ 40 5, , ? 0 5, 0 5, ? 0 , ? 0 5 0 0 ?9

'23 ;2'

Interplay between Polarizability and Hydrogen Bond Network of Water

Recommend Documents