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
18 Table 2. Generational force fields (cont.)
𝑘𝑟 𝑟𝑒 𝑘𝜃 𝜃𝑒 𝐴 𝐶 𝑞𝑂 𝑞𝐻
19 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 Lennard-Jones 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.
Figure 3-1-1. AFMM trajectory of L-J parameters G00
20 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.
21 Table 3. Second batch generational force fields
Figure 3-1-2. Second AFMM trajectory of L-J parameters G00*
22 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
23 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,
24 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.
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.
0
25
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 importance when compared to the other aspect of the force field.
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.
0
26
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.
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.
0 0.5 1 1.5 2
1 2 3 4 5 6
g(r)
r/Å
27
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
28 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 many-body effect itself. Whereas Wang Feng’s linear adaptive force-matching
29 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.
30 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.
-0.5 0 0.5 1 1.5 2 2.5
0 5 10 15 20 25 30
Free Energy/(Kcal/mol)
Distance from center/Å
31 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
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.
32 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.
33
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.
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
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 Hydrogen Density/(1/Å3 )
x/Å
34 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.
35
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 molecule-molecule 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 many-body 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.
36
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.
10. Nosé, S. , Journal of Chemical Physics 1984, 81, 511–519.
11. Møller, C.; Plesset, M. S. , Phys. Rev. 1934, 46, 618–622.
37
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.