為非線性自適應力批配演算法打下基礎：具效益、系統化且從第一原理延伸的簡單水模型重參數化

全文

(1)國立台灣師範大學化學系碩士論文. 指導教授：蔡明剛博士. 為非線性自適應力批配演算法打下基礎：具效益、系統化且從第一原理延伸的簡單水模型重參數化. Laying the Foundations of a Non-Linear Adaptive Force Matching Algorithm: An Economic, Systematic, & Ab Initio Derived ReParameterization of a Simple Water. 研究生：黃怡碩. 中華明國一百零五年七月.

(2) Table of Contents Table of Figures ....................................................................................... 3 中文摘要 ................................................................................................. 4 Abstract ................................................................................................... 5 Chapter 1 Introduction ............................................................................. 6 Chapter 2 Methodology ........................................................................... 9 2-1. The Adaptive Force-Matching Method ....................................... 9 2-1-1. Step 1: Force Field Validation .......................................... 9 2-1-2. Step 2: Molecular Dynamic Simulation .......................... 10 2-1-3. Step 3: Snapshot Sampling. ............................................ 11 2-1-4. Step 4: QM/MM Force Calculations. ............................. 12 2-1-5. Step 5: Data Weighting ................................................... 12 2-1-6. Step 6: Force Matching .................................................. 13 2-2. Radial Distribution Function .................................................... 14 2-3. Potential of Mean Force ........................................................... 15 Chapter 3 Results & Discussion ............................................................. 17 3-1. The Generated Force Fields...................................................... 17 3-2. Radial Distribution Function .................................................... 23 3-2-1. Oxygen-Oxygen RDF .................................................... 23 3-2-2. Hydrogen- Hydrogen RDF ............................................. 25 3-2-3. Oxygen-Hydrogen RDF ................................................. 26 3-3. Cluster Energy Calculation ...................................................... 27 3-3-1. Energy Component Analysis .......................................... 28 3-4. Potential of Mean Force ........................................................... 29 1.

(3) 3-5. Slab Integrity ........................................................................... 33 Chapter 4 Conclusion ............................................................................. 35 Chapter 5 References ............................................................................. 36. 2.

(4) Table of Figures Figure 3-1-1. AFMM trajectory of L-J parameters .................... 19 Figure 3-1-2. Second AFMM trajectory of L-J parameters ....... 21 Figure 3-2-1. O-H RDF. ............................................................ 24 Figure 3-2-2. H-H RDF .......................................................... 25 Figure 3-2-3. O-H RDF. .......................................................... 26 Figure 3-3. The 17-water-clusters chosen ............................... 27 Figure 3-4. Potential of mean force free energy graph. ........... 30 Figure 3-4. Molecular simulations for slab. .............................. 31 Figure 3-5. Hydrogen atom linear density chart ...................... 33. 3.

(5) 中文摘要. 非線性的力參數在參數化是目前計算化學的重大難題，由於非線性的最佳化方法可以產生的錯誤極多。雖然我們可以選擇用線性的最佳化方法，但是如果用在非線性的力參數上，會有系統上的誤差。本文希望提供一個升級過的自適應力批配演算法 (Adaptive Force Matching Method, AFMM)來處理此問題，並以液態水的分子模擬計算作為驗證。初步測驗，包含徑向分布函數以及平均力勢計算實驗，顯示結果有部分地被改良。. 本文產生的簡單點電荷模組力參數，提供更正確的第一溶劑殼，以及更準確地描述成對力。然而，這些優勢也帶來了一些缺點，本文產生的力參數無法提供正確的第二溶劑殼，也同時完全放棄多體效應。目前已有文獻提出簡單點電荷模組無法完整的描述液態水的性質，例如第一溶劑殼及第二溶劑殼，本研究試圖在此方面提供更徹底的探討。. 關鍵字：力批配，第一原理，分子模擬，液態水. 4.

(6) Abstract On the subject of force field re-parameterization, the Annual Reports in Computational Chemistry, Volume 10 stated that “Non-linear optimization is a grand challenge.” Indeed, the complication and errors of such a non-linear method is not easily fixed. Although linear methods can be used instead, this is at the cost of less precise results. To tackle this issue, an upgraded Adaptive Force Matching Method (AFMM) is proposed based on testing with liquid water model. Preliminary results from radial distribution functions and potential of mean force calculations suggest a partially improved result from previous experiments.. It would seem that the generated simple point charge water force field from AFMM, has a better first solvation shell, and a better pair-wise force description. This, however, came at the cost of a poorer second solvation shell fit, and a forfeit of many-body effects. It has been shown that a simple point charge model of water cannot effectively cover both solvation shell, and this work may provide additional insight on the matter.. Keyword: Force Matching, Ab Initio, Molecular Dynamics Simulation, Liquid Water. 5.

(7) Chapter 1 Introduction The Adaptive Force-Matching Method (AFMM) can reliably generate molecular mechanical (MM) force field (FF) parameters, without empirical dynamics observations, for most given systems, and at a remarkably low cost.. [1] [2] [3] [4] [5] [6] [7]. This is achieved by fitting the parameters of a. molecular mechanical force field, en masse, to a data sample of ab initio calculated forces of a small but representative system, the structure of which are in turn generated by prior force field parameters via molecular dynamic (MD) simulation. Once the fitting is done, a new set of parameter is conceived, and thus can be used again in another MD simulation, whose structure can be then sampled for a new group of ab initio force calculation, facilitating the fitting of a fresh set of, ideally improved, parameters. Through iteration, a viable force field parameter set will eventually be generated, even if the individual improvement over each cycle is limited.. Such a force field parameter set is crucial to the efficacy of its corresponding molecular mechanical model for molecular dynamic simulations, to the point where many refer to such simulations as of Force Field Molecular Dynamics (FFMD). Due to the classically physical nature of FFMD simulations, they are arguably the only viable way of computational studies where a large system is unavoidable, such as the case of heterogeneous catalysis systems where a huge simulation box has to be used to depict the nature of surfaces.. 6.

(8) While there are many details in the mechanisms and history of AFMM, the most relevant part for our purpose right now is the curve fitting process found in the core of AFMM, said process, to our knowledge, has always been chosen to be suited for statistically linear equation systems, e.g. Singular Value Decomposition (SVD) or QR Decomposition. [4] [5] [6] This research direction decision primarily stems from the desire to avoid complications that may arise from using non-linear fitting methods, even though the employment of fitting processes assuming otherwise will likely damage the epistemological integrity of AFMM when the force field model in question is not really linear.. Therein lies the central problem: most force fields models do not have statistical linearity, and using linear fitting regardless would yield an invalid result in worst case scenarios. Admittedly, the results generated by AFMM so far (i.e. with linear fitting) has been acceptable — meaning that the assumption of linearity is a viable approximation in certain cases, bearing in mind that only simpler force field models have been parameterized this way; nonetheless, if a nonlinear fitting process can be utilized in AFMM, such a ‘pseudo-linearity’ approach ultimately may not be needed.. Our main contribution thus is upgrading AFMM with a nonlinear fitting method, to hopefully create a more legitimate methodology and to broaden its applicability. In doing so we employed various solutions to overcome aforementioned complications related to nonlinear fitting. 7.

(9) To that end, the primary subject of our study is liquid water. Other then it’s physical ubiquitousness and importance, it is an oft studied system, with ample literary references. The apparent simplicity of the water molecule also lend to less complicated modelling, making it a suitable option as a starting point for building a re-parameterization method.. 8.

(10) Chapter 2 Methodology 2-1. The Adaptive Force-Matching Method The Adaptive Force-Matching Method (AFMM) is an automated iterative process of six steps aiming to render a force field parameter set derived from ab initio level calculations. Every six steps is grouped as a “generation”, and every generation observes its own associated force field while providing the next generation a new, ideally improved, force field. With this work, our goal is to model water with a flexible single point charge model via AFMM.. 2-1-1. Step 1: Force Field Validation A force field model and parameter set is provided and checked, resulting in a new generation and possibly halting the circular process with completion declared.. At the very beginning, an initial guess force field is manually provided and no automated validation is needed. Said initial guesses can likely be acquired from previous literature, though an educated guess should also suffice.. If the parameters provided came from the last generation, then its convergence with the parameters of previous generations will be checked. 9.

(11) Upon convergence, AFMM is finished and the latest force field, the final generation, is declared as the result; otherwise the provided force field will be used in the now-current generation.. With this work, we have employed Wang Feng’s initial-guess parameters as ours. [4] 1. 1. 2. 2. 2. 𝐸 = ∑ [ 𝑘𝑟 (𝑟𝑖 − 𝑟𝑒 )2 + 𝑘𝜃 (𝜃𝑗 − 𝜃𝑒 ) +. 𝐴 𝑟𝑘. 12 −. 𝐶 𝑟𝑘. 6 +. 𝑞𝑚 𝑞𝑛 𝑟𝑚𝑛. ]. Where 𝑘𝑟 is the linear force constant, 𝑟𝑒 is the equilibrium bond length, 𝑘𝜃 is the angle force constant, 𝜃𝑒 is the equilibrium bond angle, 𝐴 and 𝐶 are respectively the Lennard-Jones repulsion and attraction constant for oxygen,. 𝑞𝑚 𝑞𝑛 𝑟𝑚𝑛. is the columbic energy term, and 𝑟𝑖 , 𝜃𝑗 , 𝑟𝑘 , and 𝑟𝑚𝑛 are. respectively the bond length, bond angle, oxygen-oxygen intermolecular distance, and partial charge distance arguments.. In this work, the halting is manually initiated and executed, and is based on the preliminary observation of convergence in the radial distribution function across each generation.. 2-1-2. Step 2: Molecular Dynamic Simulation A molecular dynamics simulation of the system in question is run with the current generation force field. In doing so a thermodynamic population of possible configuration states of the system is generated.. This work uses the TINKER molecular modeling software for its 10.

(12) molecular dynamic calculations, [8] The MD model is a 18.6445 Å periodic water cube with 216 water molecules, all inter molecular forces is cut-off at 9.0 Å . Lennard-Jones potential and Ewald Summation [9] are both utilized. The epsilon combination rule for Lennard Jones potential is geometric, meaning Oxygen-Hydrogen interactions will not include van der Waals forces. All simulation is done under canonical ensemble with Nosé-Hoover thermostat. [10]. (tau=0.2) set at 300K. Each simulation lasts 40 simulated. picoseconds with a time-step of 0.5 femtosecond, structures are saved every femtosecond.. 2-1-3. Step 3: Snapshot Sampling. With the molecular dynamic trajectory generated from the previous step, one can provide statistically significant samples of snapshots.. Additionally, each time frame can provide multiple samples by focusing on different molecules. As step 4 involves QM/MM calculation, molecules “in focus” will be calculated ab initio, while those “out of focus” but not “out of frame” will be handled by molecular dynamic calculation using the force field of the current generation; which is why the former is dubbed the QM region, and the latter MM region. Further details on the QM/MM model will be covered in the next section. In this work, we picked the very last structural frame of the MD trajectory created from the previous step as the basis for all of the samples 11.

(13) of that generation, with 96 different QM region chosen from it. Our QM region has 7 water molecules, which was deemed sufficient by Wang Feng after special weighing (covered in step 5), and contains 6 closest water molecule to the central water molecule, or the focal point, of that QE region. In every generation the 96 water molecule closest to an arbitrary origin is chosen to be the focal points of their corresponding QM Region. All water molecules not included in the QM region is included in the MM region.. 2-1-4. Step 4: QM/MM Force Calculations. The provided samples of snapshots will be calculated and the forces applied to the QM region atoms is recorded. Because the calculation of each sample is not dependent on each other, this step can be perfectly parallelized.. In this work, Gaussian 09 is used for QM/MM force calculations, with the QM part calculated with MP2 [11]/AUG-cc-pVDZ [12] [13], and the MM part with the force field of current generation along with electronic embedding.. 2-1-5. Step 5: Data Weighting The parallelized computational results of the previous step is assembled into a collection of ab initio forces linked to relevant configurational snapshot. Atoms deemed to be more or less representative 12.

(14) of the system will be weighted accordingly.. In this work we utilize Wang Feng’s weighing rule, [4] where a water molecule’s weight is dependent on how many of its first solvation shell neighbors are included in the QM region. If there is three or four such neighbor included in the QM region, then the relevant water molecule will be weighted as standard. With only two neighbors, the data is halfdiscounted in weighting. With one or less the molecule is completely ignored.. 2-1-6. Step 6: Force Matching Finally, the force data from step 4 will be used as the target for the new force fields’ re-parameterization via fitting algorithms, using the relevant configurational snapshot (as chosen in step 3 and employed in step 4) as context, and with considerations from the weighting of step 5. Once the fitting is done, a new set of force field is created, step 1 is revisited, and in the process starting a new generation.. In this work we employ a non-linear fitting process, L-BFGS-B, provided by the Python module SciPy and NumPy.. 13.

(15) The objective function has three terms, each covers atomic forces, inter-molecular forces, and inter-molecular torque.. = √∑ 𝑤𝛼 |. 𝑓𝑖𝑡. 𝑓𝛼𝑎𝑖 −𝑓𝛼. √∑ 𝑢𝑖 |. 𝑓𝛼𝑎𝑖. |2 / ∑ 𝑤𝛼 𝑛𝛼. 𝑓𝑖𝑡. 𝜏𝑖𝑎𝑖 −𝜏𝑖 𝜏𝑖𝑎𝑖. + √∑ 𝑤𝑖 |. 𝑓𝑖𝑡. 𝑓𝑖𝑎𝑖 −𝑓𝑖 𝑓𝑖𝑎𝑖. |2 / ∑ 𝑤𝑖 𝑛𝑖 +. |2 / ∑ 𝑢𝑖 𝑛𝑖. 𝑤𝛼 , 𝑤𝑖 , 𝑢𝑖 is the weight assigned from Step 5, 𝑛𝛼 and 𝑛𝑖 are the amount of atoms and molecules considered.. The input data is collect from the previous 4 generations of configurations and ab initio calculations, this is done to give the algorithm some inertia, and to prevent extreme fluctuation in parameter values across generations.. 2-2. Radial Distribution Function A Radial Distribution Function (RDF) of a configuration describes the probability of a specific type of atom being found around another type of atom. RDF is frequently used to represent the structural property of fluids, and thus we employ it to validate the force field of this work. An RDF can be easily extracted from molecular dynamic simulations by integrating data over the duration of a simulation of a fluid after equilibrium is reached. RDF can also be interpreted from physical experiments such as X-ray or 14.

(16) neutron scattering, though the operating principles of these methods is beyond this work.. On the subject of RDFs of water, it is a topic of study for over a couple of decades, [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] and during this period of time, the pinnacle of physical experiments was conducted and released back in 2000 by A. K. Soper.. [33]. A more recent. revision with updated methodology, however, has prompt A.K. Soper to publish a rather different RDF results in 2013, and has accordingly and aptly asked the following question in the subtitle of the paper: “Is There Anything We Can Say for Sure?” [34]. It is one of this work’s aim to contribute on this topic by providing additional certainty in the structure of water, via an updated ab initio method, and if not immediately and directly, eventually and indirectly as a stepping stone through some of its innovation.. 2-3. Potential of Mean Force This work employs one final in silico experiment in the form of potential of mean force (PMF) experiments. This type of experiment is designed to find the free energy required for a selected group of atoms to go through a specific path, spatial or chemical. This is done through the summation or integration of the mean forces the subject is influenced by, throughout different points of its trajectory. As integration of force over 15.

(17) distance is work, assuming the absence of pressure-volume and other forms of work, the work we’ve calculated is the Gibbs/Helmholtz free energy. (Note that since pressure-volume work is the sole difference among Gibbs and Helmholtz free energy, they are indistinguishable here.). In this work, the force field acquired through AFMM is also tested by PMF experiments. Among the bulk of water of 729 molecules is a sole water molecule affixed to a specific point among a leaving trajectory. This 730th water molecule is referred to as the probing molecule, and for every single ångström that it moves on its trajectory, each of the configuration underwent its own MD simulation, which provides the mean force data needed to conduct PMF. All of the simulations was down with 2fstimesteps over 400ns under canonical ensemble, with 730 water molecules in periodic boundary condition of (27.96675 Å x83.90025 Å x27.96675 Å ), at a temperature of 300K.. To physically constrain a water molecule through a pre-planned path and make measurements is beyond the current capabilities of contemporary science, so this experiment can only be conducted in silico. Nonetheless high level ab initio calculations of PMF concerning water do already exist, and certain physical experiments, such as solvation experiments, can be used to cross-reference the calculated results.. 16.

(18) Chapter 3 Results & Discussion 3-1. The Generated Force Fields The following tables show the 40 generations of force fields generated by this work’s AFMM, each labeled GXX where XX is the number of the generation. The initial guess, G00, is the same as the ones used by Wang Feng (2008). We have chosen the tenth generation (G10) of this batch of force matching as our definitive result, as additional cycle may cause over fitting, and both the radial distribution function and the parameters themselves has been deemed sufficiently converged. For completeness sake however, we’ve ran the algorithm until the fortieth generation. Table 1. Generational force fields 𝑘𝑟. 𝑟𝑒. 𝑘𝜃. 𝜃𝑒. 𝐴. 𝐶. 𝑞𝑂. 𝑞𝐻. G00 1232.580. 0.957. 99.480. 104.390. 629.327. 625.502. -0.662. 0.331. G01 1241.182. 0.964. 84.642. 107.320. 689.924. 567.181. -0.797. 0.399. G02 1183.727. 0.964. 85.922. 107.386. 676.177. 486.275. -0.802. 0.401. G03 1161.741. 0.964. 87.783. 107.409. 671.411. 419.640. -0.814. 0.407. G04 1145.404. 0.963. 86.665. 107.501. 661.012. 360.593. -0.813. 0.407. G05 1120.199. 0.963. 87.679. 107.657. 661.447. 312.386. -0.825. 0.413. G06 1118.590. 0.963. 87.909. 107.680. 658.483. 269.927. -0.832. 0.416. G07 1121.824. 0.963. 87.140. 107.885. 661.043. 234.218. -0.839. 0.420. G08 1124.652. 0.962. 87.029. 108.043. 661.373. 202.889. -0.852. 0.426. G09 1114.342. 0.962. 85.748. 108.099. 658.480. 175.322. -0.853. 0.426. G10 1118.426. 0.962. 85.694. 108.226. 656.149. 151.565. -0.857. 0.429. 17.

(19) Table 2. Generational force fields (cont.) 𝑘𝑟 G11 G12 G13 G14 G15 G16 G17 G18 G19 G20 G21 G22 G23 G24 G25 G26 G27 G28 G29 G30 G31 G32 G33 G34 G35 G36 G37 G38 G39 G40. 1110.151 1118.589 1113.323 1112.273 1115.656 1108.327 1106.132 1105.337 1098.446 1088.405 1094.833 1087.198 1089.295 1095.699 1099.452 1092.206 1087.246 1076.581 1084.395 1115.313 1125.147 1127.240 1120.246 1111.794 1107.162 1110.032 1112.152 1099.175 1098.673 1112.948. 𝑟𝑒 0.962 0.963 0.963 0.963 0.963 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.963 0.963 0.963 0.963 0.962 0.962 0.962 0.962 0.962 0.962 0.962. 𝑘𝜃 84.913 87.011 86.905 86.470 85.828 84.739 85.779 86.112 86.566 86.656 85.774 86.140 85.647 85.712 86.103 86.522 87.448 84.882 86.608 86.396 87.074 88.814 88.905 88.251 87.732 87.043 85.353 85.156 84.791 85.282. 𝜃𝑒 108.123 107.958 107.890 107.913 108.024 108.232 108.274 108.264 108.230 108.124 108.214 108.171 108.256 108.239 108.224 108.170 108.014 108.190 108.105 108.140 108.179 108.058 108.037 108.062 108.064 108.100 108.168 108.266 108.277 108.222. 𝐴 648.784 649.353 643.899 633.852 632.241 634.749 636.390 636.530 636.100 637.413 642.803 651.674 652.519 657.083 659.703 653.794 648.675 643.065 642.047 640.624 644.292 641.770 647.447 651.196 651.924 649.798 632.388 626.335 613.266 605.382. 𝐶 174.027 201.037 231.161 215.093 186.039 161.434 139.986 121.245 115.540 117.071 135.752 153.974 157.526 182.531 189.760 163.599 150.194 129.508 136.858 157.855 182.797 210.662 219.219 231.879 200.925 173.722 148.419 127.918 109.618 94.320. 𝑞𝑂. 𝑞𝐻. -0.847 -0.843 -0.836 -0.833 -0.837 -0.842 -0.847 -0.847 -0.848 -0.848 -0.852 -0.851 -0.854 -0.853 -0.852 -0.853 -0.847 -0.846 -0.844 -0.844 -0.848 -0.844 -0.848 -0.849 -0.849 -0.854 -0.849 -0.847 -0.846 -0.844. 0.424 0.421 0.418 0.417 0.419 0.421 0.423 0.424 0.424 0.424 0.426 0.426 0.427 0.427 0.426 0.426 0.423 0.423 0.422 0.422 0.424 0.422 0.424 0.424 0.425 0.427 0.425 0.424 0.423 0.422. 18.

(20) We can see that, after the tenth generation, the parameters 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 LennardJones equation is already trying to depict: the Van der Waals force. Thus, these values wax and wane slightly across the generations, finding their “jurisdictions” somewhat fuzzy. Because the two parameters pertaining to the Lennard-Jones equation covers the repulsion and attraction of VdW force respectively, by charting a trajectory of all A and C values through the generations, we can more intrinsically understand the generational changes of the force fields.. 700.000 600.000 G00 500.000 400.000. C 300.000 200.000 G10. 100.000 0.000 600.000. 620.000. 640.000 660.000 A. 680.000. 700.000. Figure 3-1-1. AFMM trajectory of L-J parameters 19.

(21) We can see that the initially guessed parameters are noticeably away from the rest of the trajectory. After about ten generations, during which the parameter set serpentines as it gets over-corrected slightly every time, 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) for the sake of expedience, and 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. This is likely to occur if the initial guess for the parameter is too far off the final result, as it indeed is in this case.. To further validate the decision to use the tenth generation parameters, and this variant of AFMM in general, we ran a second batch of 40 generations using the same method, with 2 changes. First being the initial guess is replaced by the final result of Wang Feng’s force matching, the second being the introduction of a “final fit”, where all 40 generations, including the initial guess, are used to fit for the final 40th generation.. 20.

(22) Table 3. Second batch generational force fields. G00* G01* G02* G03* G04* G05* G06* G07* G08* G09* G10* G11* G12* G13* G14* G15* G16* G17* G18* G19* G20*. 𝑘𝑟 1095.620 1106.603 1119.880 1115.979 1112.878 1111.776 1106.486 1101.775 1104.689 1101.670 1104.206 1112.462 1114.445 1115.901 1114.314 1116.385 1115.440 1121.378 1115.421 1112.901 1109.757 300. 𝑟𝑒 0.960 0.961 0.961 0.961 0.961 0.961 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962 0.962. 𝑘𝜃 85.000 84.121 83.749 84.170 84.980 85.166 87.350 88.041 87.555 87.424 86.483 86.018 85.933 86.547 86.200 85.456 85.571 84.216 83.952 84.309 84.780. 𝜃𝑒 109.350 109.225 108.884 108.808 108.739 108.511 108.389 108.307 108.205 108.197 108.206 108.195 108.212 108.179 108.155 108.221 108.197 108.203 108.209 108.171 108.108. 𝐴 510.720 572.802 591.227 595.063 597.766 625.763 628.336 642.878 646.700 656.358 665.570 669.440 668.048 658.741 644.942 639.200 625.410 619.361 618.817 622.603 639.648. 𝐶 257.475 236.144 207.770 180.517 156.687 138.836 120.482 140.472 162.685 189.249 213.634 233.252 201.792 173.535 198.272 197.411 169.109 145.742 126.161 109.593 128.267. 𝑞𝑂 -0.854 -0.843 -0.850 -0.850 -0.849 -0.855 -0.854 -0.858 -0.856 -0.853 -0.850 -0.845 -0.849 -0.847 -0.847 -0.844 -0.843 -0.842 -0.840 -0.840 -0.840. 𝑞𝐻 0.427 0.421 0.425 0.425 0.424 0.428 0.427 0.429 0.428 0.427 0.425 0.423 0.424 0.424 0.423 0.422 0.421 0.421 0.420 0.420 0.420. 250 G00*. 200. G40*. C 150. G10 100 500. 550. 600 A. 650. 700. Figure 3-1-2. Second AFMM trajectory of L-J parameters 21.

(23) Table 4. Second batch generational force fields (cont.) 𝑘𝑟. 𝑟𝑒. 𝑘𝜃. 𝜃𝑒. 𝐴. 𝐶. 𝑞𝑂. 𝑞𝐻. G20* 1109.757 0.962. 84.780. 108.108 639.648 128.267 -0.840. 0.420. G21* 1112.563 0.963. 84.956. 108.152 640.010 148.152 -0.839. 0.419. G22* 1116.809 0.963. 84.375. 108.170 652.034 172.671 -0.840. 0.420. G23* 1116.323 0.962. 86.533. 107.978 657.247 200.178 -0.840. 0.420. G24* 1117.704 0.962. 88.033. 107.904 659.164 200.225 -0.843. 0.422. G25* 1100.133 0.962. 89.764. 107.815 665.410 232.293 -0.840. 0.420. G26* 1100.797 0.962. 90.710. 107.795 674.757 270.106 -0.841. 0.421. G27* 1093.768 0.962. 87.742. 107.896 669.671 237.768 -0.840. 0.420. G28* 1100.139 0.962. 84.298. 108.221 672.993 263.873 -0.841. 0.420. G29* 1109.671 0.962. 84.268. 108.174 678.789 251.480 -0.849. 0.424. G30* 1110.753 0.962. 83.361. 108.196 677.807 247.601 -0.853. 0.426. G31* 1114.685 0.962. 85.268. 108.176 670.205 213.223 -0.855. 0.428. G32* 1102.564 0.962. 85.949. 107.949 659.484 183.174 -0.845. 0.423. G33* 1104.828 0.963. 87.295. 107.964 656.577 167.388 -0.842. 0.421. G34* 1095.656 0.963. 88.186. 107.954 651.217 192.493 -0.840. 0.420. G35* 1092.737 0.963. 86.716. 108.081 656.094 223.102 -0.840. 0.420. G36* 1107.726 0.962. 86.322. 108.209 663.754 259.116 -0.845. 0.423. G37* 1110.414 0.962. 85.309. 108.209 664.400 224.510 -0.846. 0.423. G38* 1109.989 0.963. 83.163. 108.229 669.365 195.157 -0.846. 0.423. G39* 1119.049 0.963. 84.285. 107.996 658.459 167.628 -0.842. 0.421. G40* 1107.307 0.962. 84.196. 107.993 655.570 182.910 -0.841. 0.421. 22.

(24) Because the initially guessed parameter set is much closer to the final result for the second batch, we did not observe any over-fitting-induced outliers. This second batch of AFMM results are otherwise similar to the first one: the first ten generations are corrected more noticeably, after which the trajectory wanders near a “valid region”. The final fit is shown to be almost identical with the tenth generation parameter set from the first batch of AFMM, which is rather reassuring.. 3-2. Radial Distribution Function To further investigate the physical implications of these generated parameters sets, molecular dynamic simulations identical to the ones used in the AFMM was used to generate water box trajectories. 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 (labeled G00, G01, G05, G10 respectively). For comparison purposes, the radial distribution function from the final result of Wang Feng’s adaptive force matching, and Soper’s physical experimental data are also included in the graphs.. 3-2-1. Oxygen-Oxygen RDF We start with the O-O RDF, observing that the first solvation shell is noticeably better described after a single cycle of AFMM. After which the difference among generations started to dwindle, by the tenth generation, 23.

(25) the RDF is mostly converged, much like the force field parameter sets themselves as previously noted. Compared to Wang Feng’s results, our water box is not over-structured at the first solvation shell. However, Wang Feng’s parameters fit the rest of curve much better. We believe this reaffirms the idea that a simple point charge model for water, simple isn’t capable of describing both solvation shell perfectly, as its parameters doesn’t contain terms regarding many-body effects. [35] Additionally, we believe this is the first simple point charge parameter set for water that specifically matches the physically observed first solvation shell like this. We thus speculate that this parameter set could be a better choice in molecular simulations involving water where the more immediate geometry is more relevant. In other words, perhaps this new parameter set can be used when over-structuring is highly undesirable. 3 2.5 2 1.5 g(r) 1 0.5 0 2. 3. 4 r/Å. 5. 6. Figure 3-2-1. O-H RDF. Green dotted line is the final result of Wang Feng’s force matching. Orange dashed line is the experimental results from Soper. While the blue lines are G00, G01, G05, and G10 respectively from darkest to lightest in brightness.. 24.

(26) 3-2-2. Hydrogen- Hydrogen RDF Here, we see a much poorer fit from our parameter set, and Wang Feng’s result seems to fit the first solvation shell more accurately. And, surprisingly, both parameter set seems to describe the second solvation shell equally well. We would like to argue that the difference from experimental values pertaining to the first solvation shell is relatively minor, assuming this is the trade-off for an overall non over-structured configuration. Additionally, we believe that the H-H RDF mostly reflects how the water in this model turn its hydrogen away in the presence of a partially positively charged atom or ion, which should be of less. g(r). importance when compared to the other aspect of the force field.. 1.4 1.2 1 0.8 0.6 0.4 0.2 0 1. 2. 3. 4. 5. 6. r/Å Figure 3-2-2. H-H RDF 4 Green dotted line is the final result of Wang Feng’s force matching. Orange dashed line is the experimental results from Soper. While the blue lines are G00, G01, G05, and G10 respectively from darkest to lightest in brightness.. 25.

(27) 3-2-3. Oxygen-Hydrogen RDF Finally, we have the O-H radial distribution function. We believe our parameter set excels in this regard, being able to match with experimental results for both solvation shells, including the trough between the peaks, with only slight over-structuring around the first solvation shell. We believe the O-H RDF is of a much greater importance when compared with the O-O and H-H RDF results, owning to the fact that many intricate properties regarding water is predicated upon the positional relationship between its hydrogen and oxygen atoms, examples include phenomenon such as hydrogen bonding, auto-ionization, protonation, cluster structures, and many more.. 2. g(r). 1.5 1 0.5 0 1. 2. 3. 4. 5. 6. r/Å Figure 3-2-3. O-H RDF. 5 Green dotted line is the final result of Wang Feng’s force matching. Orange dashed line is the experimental results from Soper. While the blue lines are G00, G01, G05, and G10 respectively from darkest to lightest in brightness.. 26.

(28) 3-3. Cluster Energy Calculation To further investigate our force field parameters, we’ve put it through a water cluster single point calculation. Results from Wang Feng’s parameters and the MP2 ab initio method is provided for comparison. We’ve chosen 17-water clusters specifically as it’s about the size of second solvation shell.. Table 5. Cluster Energy Chart, all values in kcal/mol, ΔE is the energy difference from the most stable configuration. (kcal/mol) E_MP2. ΔE. This Work. ΔE. W., Feng. ΔE. Sphere. -1296.73226 0.000 -146.90070 0.163 -168.9233 0.641. 552’5. -1296.73064 1.015 -147.06370 0.000 -169.5644 0.000. 441’44. -1296.73016 1.319 -143.46910 3.595 -165.5717 3.993. L-shape. -1296.72951 1.729 -143.47160 3.592 -165.6175 3.947. Figure 3-3. The 17-water-clusters chosen 6. Both single point charge models, from this work and from Wang Feng’s, misidentified the most stable cluster configuration as 552’5, even 27.

(29) though the actually most stable configuration, the sphere, is a close second in both model. Similarly, both force field models misidentified the 441’44 configurations as the least stable of the group, while the truly least stable configuration, the L-shaped cluster, is a close second in both cases. However, we believe our parameters is a better result when compared to our Wang Feng’s. And most notably, under our model the energy difference between the sphere and 552’5 cluster is rather negligible. Overall, these results seem to reinforce our previous assertion that our parameter set better describes the intermolecular alignment of Oxygen and Hydrogen atoms in water. 3-3-1. Energy Component Analysis. We also noticed that our parameters seem to render much less absolute energy when compared to Wang Feng’s results. Further investigation reveals a possible explanation for the drop in our absolute energy result. If we were to break these components down,. [36]. the many-body effect. contributes about 20 kcal/mol in total absolute energy, which fits the description of our absolute energy dap. We assert that out water model omitted most if not all energy contributed from the many-body effect of water. We further hypothesize that this is a direct result of the updates we’ve added to our AFMM. We believe that, while this seems like a flaw, in a certain sense it’s an accurate and legitimate answer, because the functional form chosen, the simple point charge model, doesn’t have manybody effect itself. Whereas Wang Feng’s linear adaptive force-matching 28.

(30) seems to have accounted for all of the many-body effects, and reassigned them under the label of other forces, such as the Lennard-Jones energy term; Our AFMM seems to have filtered said many-body effect away during the process. We believe this is caused by the gradual and iterative nature of our non-linear fitting method, making it more sensitive to the discrepancy in the energy contribution of each functional term. We hypothesize that this quality of our AFMM could be instrumental in picking the best functional form for force field employed in molecular dynamic simulations, through its use in benchmarking. In the case of our simple point charge water model, the functional form chosen seems to be suboptimal, as implicated by the fact that the value of our objective function actually increased across the generations. Regardless, we believe our parameter set is making the most out of the simple functional form given, while following the physical meaning of the energy terms involved.. 3-4. Potential of Mean Force To further investigate the properties of our water model, we’ve conducted a self-solvation potential of mean force in silico experiment.. 29.

(31) 2.5. Free Energy/(Kcal/mol). 2 1.5 1 0.5 0 0. 5. 10. 15. 20. 25. 30. -0.5 Distance from center/Å Figure 3-4. Potential of mean force free energy graph. 7. The experimental self-solvation free energy of water is reported to be -6.32 kcal/mol, which is quite different from our water model’s -2 kcal/mol. [7]. However, according to in silico experiments, as the probing cluster. increase in size up to 18 molecule cluster, about the size of second solvation shell, the average self-solvation free energy per fringe molecule starts to decline, and it just so happens to approach to -2 kcal/mol.. 30.

(32) Figure 3-4. Molecular simulations for slab. 8 Four of the molecular simulations used in the potential of mean force experiment. The configurations are stack over each other, and are tagged by color.. We believe this is further evidence suggesting that our AFMM have successfully rooted out many-body effects in the parameter set, because as the probing cluster grow in size, the many-body effect among the probing and non-probing molecules starts to rapidly decrease. While molecules at the fringe of the cluster will contribute the bulk of the intermolecular energy between cluster and bulk. To elaborate, the Jaguar model gave the 31.

(33) result of -36.48 kcal/mol self-solvation free energy for an 18 molecule water cluster, and since all of these water molecules are on the surface, we have about -2 kcal/mol self-solvation free energy per fringe molecule. We believe this suggests that our water model is one that exhibit no bulk water properties, and as such could be a better description for water found on the interface between liquid water and a protein, for example.. 32.

(34) 3-5. Slab Integrity To validate the legitimacy of our potential of mean force experiment, we’ve checked for the structural integrity of our water slab.. 0.07. Hydrogen Density/(1/Å 3). 0.06 0.05. 0.04 0.03 0.02 0.01 0 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 x/Å. Figure 3-5. Hydrogen atom linear density chart 9 Binned every .125Å , sampled 2000 frames, once every 100 2fs-timestep over 400ns. The blue line is the observed results, while the orange one is the experimental result for the bulk region.. We can see that our water model indeed displays condensed phase properties, with a rapid drop in density near the surface, and a flat plateau in the bulk region. Note however, without the constraint of an accurate 33.

(35) periodic boundary condition, the density is only two thirds of that of experimental values. This is coherent with the weakened Lennard-jones attraction presented in our parameter set, and the fact that the boundary condition is a part of the AFMM conditions. Our model is not fully applicable in this case because it wasn’t optimized for this purpose. One of the greatest advantages of AFMM is its ability to produce force fields tailored to specific configurations and situations, but this at the cost of general transferability.. 34.

(36) Chapter 4 Conclusion We’ve built an upgraded version of the adaptive force matching algorithm method, and have generated a simple point charge water model, based on ab initio calculations. Said model seems to have ignored the many-body effect presented in the ab initio calculations, but this is in accordance with the selected functional force of our force field. Radial distribution function suggests a better description of first solvation shell in terms of oxygen-oxygen distribution, which is essentially moleculemolecule distribution in the case of water. It also has a decent description of oxygen-hydrogen distribution, which we believed to be of the most importance.. In single point cluster energy calculation, we’ve discovered the primary featured of this new water model, its completely lack of manybody effect, instead of reassigning it to terms not designed for it.. Finally, in the potential of mean force experiment, we’ve realized the limits of the AFMM, the lack of transferability in its generated force fields. We believe our updated method is a success in general, and is excited for future work that make use of its ability to generate force fields that’s tailored for specific cases.. 35.

(37) Chapter 5 References 1. Pinnick, E. R.; Calderon, C. E.; Rusnak, A. J.; Wang, F. , Theoretical Chemistry Accounts 2012, 131, 1-11. 2. Ercolessi, F.; Adams, J. B. , Europhys. Lett. 1994, 26, 583-588. 3. Izvekov, S.; Parrinello, M.; Bumham, C. J.; Voth, G. A. , Journal of Chemical Physics. 2004, 120, 10896-10913. 4. Akin-Ojo, O.; Song, Y.; Wang, F. , The Journal of Chemical Physics 2008, 064108, 129. 5. Akin-Ojo, O.; Wang, F. , Journal of Physical Chemistry B 2009, 113, 1237-1240. 6. Akin-Ojo, O.; Wang, F. , Journal of Computational Chemistry 2011, 32, 453-462. 7. Asthagiri, D.; Pratt, L. R.; Ashbaugh, H. S. J. , Chem. Phys. 2003, 119, 2702. 8. Ponder, J. , Tinker [Online], http://dasher.wustl.edu/tinker/ (accessed 2016).. version. 7.1;. 9. Ewald, P. , Ann. Phys. 1921, 369, 253–287. 10. Nosé, S. , Journal of Chemical Physics 1984, 81, 511–519. 11. Møller, C.; Plesset, M. S. , Phys. Rev. 1934, 46, 618–622. 12. Kendall, R. A.; Dunning Jr, T. H.; Harrison, R. J. , J. Chem. Phys. 1992, 96, 6796-6806. 13. Woon, D. E.; Dunning Jr., T. H. , J. Chem. Phys. 1993, 98, 1358-1371. 14. Soper, A. K. , Physica 1986, 136B, 322. 15. Andreani, C.; Menzinger, F.; Ricci, M. A.; Soper, A. K.; Dreyer, J. , Phys. Rev. B 1994, 49, 3811. 16. Andreani, C.; Nardone, M.; Ricci, F. P.; Soper, A. K. , Phys. Rev. A 1992, 46, 4709. 17. Andreani, C.; Ricci, M. A.; Soper, A. K. , J. Chem. Phys. 1997, 107 , 214. 18. Santoli, G.; Bruni, F.; Ricci, F. P.; Ricci, M. A.; Soper, A. K. , Mol. Phys. 1999, 97, 777. 19. Ricci, M. A.; Nardone, M.; Ricci, F. P.; Andreani, C.; Soper, A. K. , J. 36.

(38) Chem. Phys. 1995, 102 , 7650. 20. Soper, A. K.; Silver, R. N. , Phys. Rev. Lett. 1982, 49, 471. 21. Soper, A. K.; Phillips, M. G. , Chem. Phys. 1986, 107, 47. 22. Chialvo, A.; Cummings, P. T. , J. Chem. Phys. 1994, 101, 4466. 23. Kuharski, R. A.; Rossky, P. J. , J. Chem. Phys. 1985, 82, 5164. 24. Root, J. H.; Egelsta, P. A.; Hime, A. , Chem. Phys. 1986, 109, 437. 25. Okhulkov, A. V.; Damianets, Y. N.; Gorbaty, Y. E. , J. Chem. Phys. 1994, 100, 1578. 26. Postorino, P.; Tromp, R. H.; Ricci, M. A.; Soper, A. K.; Neilson, G. W. , Nature 1993, 366, 668. 27. Löffler, G.; Schreiber, H.; Steinhauser, O.; Bunsen-Ges, B. , Phys. Chem. 1994, 98, 1575. 28. Pusztai, L. , Phys. Rev. B 1999, 60 , 11851. 29. Bellisent-Funel, M.-C.; Sridi-Dorrez, R.; Bosio, L. , J. Chem. Phys. 1996, 104 , 1. 30. Soper, A. K.; Bruni, F.; Ricci, M. A. , J. Chem. Phys. 1997, 106, 247. 31. Jedlovsky, P.; Bako, I.; Palinkas, G.; Radnai, T.; Soper, A. K. , J. Chem. Phys. 1996, 105 , 245. 32. Soper, A. K. , Chem. Phys. 1996, 202, 295. 33. Soper, A. K. , Chem. Phys. 2000, 258, 121. 34. Soper, A. K. , ISRN Physical Chemistry 2013, 2013, 67. 35. Nicolini, P.; Guàrdia, E.; Masia, M. , J Chem Phys. 2013, 139, 184111. 36. Wang, F.-F.; Deible, M. J.; Jordan, K. D. , J. Phys. Chem. A 2013, 117 , 7606–7611.. 37.

(39)