Methodology of Parameterization of Molecular Mechanics Force Field From Quantum Chemistry Calculations using Genetic Algorithm: A case study of methanol

In molecular dynamics (MD) simulation, force field determines the capability of an individual model in capturing physical and chemistry properties. The method for generating proper parameters of the force field form is the key component for computational research in chemistry, biochemistry, and condensed-phase physics. Our study showed that the feasibility to predict experimental condensed phase properties (i.e., density and heat of vaporization) of methanol through problem specific force field from only quantum chemistry information. To acquire the satisfying parameter sets of the force field, the genetic algorithm (GA) is the main optimization method. For electrostatic potential energy, we optimized both the electrostatic parameters of methanol using the GA method, which leads to low deviations of between the quantum mechanics (QM) calculations and the GA optimized parameters. We optimized the van der Waals (vdW) parameters both using GA and guided GA methods by calibrating interaction energy of various methanol homo-clusters, such as nonamers, undecamers, or tridecamers. Excellent agreement between the training dataset from QM calculations (i.e., MP2) and GA optimized parameters can be achieved. However, only the guided GA method, which eliminates the overestimation of interaction energy from MP2 calculations in the optimization process, provides proper vdW parameters for MD simulation to get the condensed phase properties (i.e., density and heat of vaporization) of methanol. Throughout the whole optimization process, the experimental value were not involved in the objective functions, but were only used for the purpose of justifying models (i.e., nonamers, undecamers, or tridecamers) and validating methods (i.e., GA or guided GA). Our method is general and can be extended to any type of small organic molecule systems, as well as other descriptive polarizable force field.


INTRODUCTION
Molecular mechanics (MM) modeling have been extensively employed in the research of chemistry, biochemistry, and condensed-phase physics, due to its capability of providing atomistic resolution data within relatively small to moderate computational cost. The accuracies of the MM calculation results relies on several important factors, such as sampling techniques, the underlying function forms and the parameters employed in the potential, with the latter two being the indispensable components of MM force field. Force field determines the capability of a certain model in capturing physical and chemistry properties in molecular dynamics (MD) simulations. The construction of a force field and its parameterization are non-trivial processes that require thoroughly understanding of the underlined physics and chemistry, carefully designing of parameterization protocols and appropriately applying of the optimization techniques. While several force field parameterization strategies have been widely employed and well-documented elsewhere, 1 these approaches rely more or less on both experimental and quantum mechanics (QM) calculated properties. The present study intends to extend the existing scopes of force field developing strategies, to explore the possibility to use only QM calculation approaches for force field parameterization.
In this study, we probed the procedure for parameterization of the polarizable force field in the form of Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) based on QM calculations. Comparing with non-polarizable force fields, which describes the system dipole moments in an averaged and thus fixed fashion, polarizable models were designed in a way, so that they can effectively describe the capability of chemistry species to change electronic distribution under the influence of its electrostatics environment. 2 The polarizable models represented a significant improvement in the field of theoretical chemistry since last decade and had led to important scientific insights such as in the study of ionic liquids and ligand binding to biological macromolecules. 3 Knowledge on efficient and accurate parameters development of polarizable models will lead to more discoveries in biophysics, biochemistry and material science. There are several polarizable force fields available for simulating various chemical and biomolecular systems, such as the Drude polarizable force field that utilizes charges on springs, 2a, 2b, 4 the CHARMM charge equilibration (CHEQ) model that employs fluctuating charges, 5 the AMBER ff02 using partial charges with inducible dipoles, 6 and the AMOEBA force field that incorporates contributions from monopoles (charges), dipoles and quadrupoles, that are of higher order multipoles. 2c, 7 The AMOEBA polarizable force field has illustrated its capability of providing reasonable predictions of interaction energy and structural properties of small molecules and clusters in gas phase, as well as the physicochemical properties of bulk phase, such as density, self-diffusivity and static dielectric constants. 3e In order to improve the ability to describe the physical models, AMOEBA polarizable force field invokes the utilization of a vast number of variables at the atomic site. The increments of the dimension of variables lead to a series of challenges in parameter searching and thus have been a focus of the present study. To address this issue, we implement the genetic algorithm (GA) to approach the global optimization regarding parameterizations. GA is an evolutionary algorithm that mimics the process of natural selection. 8 For traditional parameters optimization, gradient-based algorithms, such as the steepest descent method, conjugate gradient method, etc. are the standard approach. 9 However, for the gradient-based method, the optimal result depends strongly on the existing knowledge of the objective function. For example, the convergence of the optimal parameters is not guaranteed when the surface of the objective function is considerably rugged or non-differentiable, or when the initial value deviates significantly from the global minimum. In parameters searching of force field, the value range of parameters exponentially expands as the complexity of the force field functional form increases. To tackle the optimization algorithm issue, statistical search heuristic (such as GA, simulated annealing, etc.) is much needed to drive the parameter space towards the optimal region for the objective function. One successful example is using the GA to determine the parameters for the reactive force field (ReaxFF) for materials system. 10 Another issue to address in this study is the selection of reference data. Force fields that have been extensively employed and validated are often based on reference data from both experiments and theoretical calculations. For example, the development of almost all the widely employed water models 2c, 11 are based on achieving a balance between the high level ab initio QM calculated dataset and condensed phase properties of bulk water measured by experimentalists. Calibration of force field parameters can be achieved by minimizing the difference of the desired properties between the MM calculations and QM calculations/experimental data. 12 For small molecules and clusters in gas phase, QM calculations are often straightforward and within reasonable computational cost. While the evaluation of the bulk phase properties from QM calculations is often computationally unfeasible, thus referring to experimental data (such as densities, dielectric constants, etc.) is necessary. However, for many molecules and especially mixtures system, experimental data are either unavailable or have significant discrepancies between different studies, due to difficulties in achieving precise measurements. Thus it is necessary to explore the possibility to select reference data of small molecules or clusters solely from theoretical calculations, i.e. ab initio QM calculations, as the training dataset for force field development to make the predictive calculation in computational modeling.

METHODS
The principal goal of this work is to establish guidance for constructing a framework of force field parameters development using primarily QM calculation data. In this initial attempt, we selected liquid methanol as the model system. Methanol represents a molecule of significant importance to organic liquids and polymers that contain bulky alkyl moieties. The hydroxyl group leads the molecule to the polar nature. The formation of the methanol clusters primarily comes from the contribution of hydrogen bond networks. 13 These ubiquitous aspects make methanol as an excellent candidate for designing and examination of force field development.
Ren et. al. suggested a general protocol of parameter development for the AMOEBA polarizable force field. 7,14 In the AMOEBA model, the electrostatic potential explicitly considers the effect of atomic multipoles (charge, dipoles, and quadrupoles) and polarization effect on induced dipole. The van der Waals (vdW) term is modeled as a buffered 14-7 potential proposed by Halgren. 15 Compared to the 12-6 potential, the buffered 14-7 potential has the advantage of optimizing structures with initial crude geometries. For hydrogen atoms, an additional parameter (the so-called reduction factor) is used to scale the position of the hydrogen atom interacting site along the corresponding covalent bond. It is meant to reflect the degree of which hydrogenic electron density displaced toward the heavy atom when covalent bonding takes place. 3e The full potential energy function of the AMOEBA force field has been described in full detail elsewhere. 2c, 14 The general energy expression of the AMOEBA force field is written in equation (1), where parameterizations of electrostatic potential energy E ESP and vdW interaction energy E vdW terms are studied in this work. Bonded interaction, including bond, angle, dihedral and torsional parameters, were kept consistent with that in the original AMOEBA force field (amoeba09.prm). 7 In Ren et. al.'s protocol, the parameters of multipoles and torsional bonded terms are optimized using QM results. However, besides of taking QM calculations as the reference, the protocol adapts data from existing database of experimental values, which are not entirely available for every desired system. We followed their approaches with the context of parameterizing standard AMOEBA force field for atomic multipoles, polarizabilities and vdW interaction by GA methods using solely results from QM calculations, which are available for the methanol system within given computational resources.
All the QM calculations were performed with Gaussian 09 16 package of electronic structure programs. The second order Møller-Plesset perturbation theory 17 (MP2) was employed, with the basis set superposition error (BSSE) correction. 18 This level of theory and basis set to predict energetics and structural properties has been widely verified by previous studies. 19 The electronic density calculation of methanol monomer is initially used for calibrating the magnitude of the multipole moments at the atomic positions.
Combining with the Distributed Multipole Analysis (GDMA 2.2) tool developed by Anthony Stone, 20 the Tinker 2c, 21 package is used to optimize the permanent atomic multipoles by fitting the electrostatic potential measured on the Connolly surface of the methanol monomer. Table S1 in the Supporting Information is showing the electrostatic potential results for methanol monomer using different basis sets.
In the parameterization process, the training dataset for the optimization should contain a good sampling of the possible structural configurations and their respective energies. For instance, to predict the polarizable effect between molecules, methanol dimers with a continuous range of electrostatic potential energies ( E ESP ) should be included. In AMOEBA, to compute clusters electrostatic potential energy ( E ESP ) including the polarizable effect between methanol molecules, 44 independent parameters are needed. Those parameters are monopole (q), dipole (µ x , µ y , µ z ), quadrupole-a traceless and systemic matrix-(Q xx , Q yx , Q yy , Q zx , Q zy ), atomic polarizability (α) and Thole's description of damping factor (a) for four types of atoms (O, H(-O), C and H(-C)) in methanol (The detailed description of the parameter form can be found in somewhere else 2c, 14 ). We optimized all 44 independent parameters listed above, using an extensive training dataset taking electronic density calculations for 4943 methanol dimers. The methanol dimer configurations were optimized via Guassian09 program at the MP2 level using 6-311G(d, p) basis set. The electrostatic potential energy ( E ESP ) of 4943 dimers was computed using the same level of theory and basis set. The detailed description of optimization for electrostatic parameters is in the Supporting Information Section S2.
To ensure an adequate representation of the various possible coordination environments and cluster sizes as well as the energy landscape, we sampled three training datasets containing different methanol homoclusters, which consisting of 9, 11, and 13 molecules referred as nonamer, undecamer and tridecamer in the descriptions below, respectively, to calibrate the vdW interaction energy ( E vdW ). For methanol, there are vdW potential depth (ε) and minimum energy distance (R * ) for four types of atoms, and additional reduction factor (λ) for two types of hydrogen atoms, which constitutes 10 independent parameters to optimize. We optimized these 10 independent parameters using interaction energy ( ΔE ) results of the sampled 502 nonamers, 157 undecamers and 94 tridecamers at the MP2/6-31G (d, p), respectively. For instance, the nonamer was sampled through randomly placing 9 C atoms in a 10 Å × 10 Å × 10 Å computational supercell, where one C atom placed in the center and surrounded by other 8 C atoms with relative distance at least of 3.25 Å, the minimum distance of carbon bond networks in methanol, 22 to each other. Then, methanol with random orientations will be assigned to the position where the C atoms sit. The methanol nonamer configuration is optimized in Gaussian09 using the eigenvalue-following algorithm 23 with constraints of fixed C atoms positions. The relaxed nonamer configurations and their interaction energies ( ΔE ) between the center methanol and surrounding methanols are then employed in the training dataset for optimizing the vdW parameters. This procedure was applied to undecamers and tirdecamers as well for getting the relaxed configuration and the interaction energy, respectively. Figure 1 (a), (b) and (c) are showing example configurations of nonamer, undecamer and tridecamer. In the AMOEBA force field, except for the electrostatic potential energy ( E ESP ), vdW interaction energy ( E vdw ) is the only term consisting of interaction energy between molecules. The vdW interaction energy ( E vdw ) can be calibrated from the interaction energy ( ΔE ) by separating apart molecules, with the assumption that the deformation energy of each molecule upon binding is reasonably small. 7 Through sampling configurations widely and calculating the coordinated interaction energies, we consider the vdW interaction by taking the effect of electrostatic interaction inclusively. The interaction energy ( ΔE ) is calculated as the total energy of the cluster subtracts the total energy of the center methanol molecule and that of the surrounding methanol clusters, as shown in equation (2). The BSSE correction is applied for getting the ΔE of every configuration. The MP2/6-31G (d, p) computed interaction energies of the 502 nonamer, 157 undecamer and 94 tridecamer configurations exhibit a Gaussian distribution as shown in figure 1 (d), (e) and (f), which signify the sampling of those configurations are adequate and in the equilibrated states.
Using the training dataset described above, we optimized the non-bonded parameters in AMOEBA force field by employing genetic algorithm; the procedure is outlined in Figure 2. Taking the optimization of vdW parameters as an illustration, we begin the optimization process by generating a population of N p = 120 parameter sets randomly, such that their values lie within physically allowable limits (parameter search ranges are listed in Supporting Information Table S2). Each set of parameters in this population is called a member. For each member i, we compute the interaction energies for all structures in the training dataset using MD simulation package Tinker 2c, 21 and evaluate the objective function Δ i given by equation (3), where ΔE j MP 2 and ΔE j AMOEBA i are the MP2 calculated and the AMOEBA calcuated interaction energies for the structure j in the training dataset, respectively. Since the electrostatic potential energy ( E j ESP ) is computed for the structure j to represent the atomic electrostatic potential precisely, the only uncertain term in the interaction energy is the vdW interaction energy ( E j vdW ), represented by a member i, a set of 10 independent parameters. Here, members are then ranked in an ascending order of Δ i . After the ranking, non-linear roulette wheel selection 24 will be performed to select the top 60% members, i.e., the ones with lowest values of Δ i , which are then subjected to genetic operations: crossover with crossover-rate 3% and mutation. These mutations introduce sufficient diversity into the population, and the non-linear selection scheme helps to avoid premature convergence of the GA run. After the genetic operations, both the old and the new members are ranked by their Δs. The best N p parameter sets (members) were then chosen to constitute the next generation. Such an optimization routine ensures that only satisfactory parameter sets survive after each generation; upon repeating this workflow for sufficient generations and sampling viable regions in the parameters space, we perform three separate GA runs starting with different random populations. From each of the converged GA run, we choose the final parameter set corresponding the lowest Δ.
To overcome one of the serious shortcomings of MP2 theory, a noticeable overestimation of the interaction energy, 25 which plays a major role in stability and formation of condensed phase molecular structure, we improved the form of objective function Δ in the GA program. We introduced a penalty parameter δ in the GA optimization procedure as the unexplained discrepancy between the AMOEBA results and MP2 calculations, as shown in equation (4), which obeys the design concept of guided genetic algorithm (guided GA). 26 This penalty parameter δ never participates as the form of AMOEBA force field in MD simulations, adds up to 11 independent parameters as a subset of i th member in the GA optimization process for vdW parameters.

RESULTS
The vdW parameters are critical in the force field for getting the correct condensed phase structures, such as right intermolecular distance, through molecular dynamic simulations.
In this section, we present the optimization results of the vdW parameters from both GA and guided GA program by calibrating interaction energy between molecules. We follow with showing binding curves for methanol dimers using both GA and guided GA optimized parameters. Finally, we show the MD simulation results (i.e., density and heat of vaporization) of a condensed phase methanol system from both GA and guided GA optimized parameters using AMOEBA force field.

VDW parameters of nonamers, undecamers, and tridecamers from GA and Guided GA
The GA and guided GA programs were applied to optimize the vdW parameters using the interaction energies ( ΔE ) calculated at MP2/6-31G (d, p) with BSSE correction for 502, 157 and 94 methanol nonamers, undecamers and tridecamers, respectively. The excellent correlations (R = 0.970, 0.976, and 0.968) of the ΔE between MP2 calculations and AMOEBA with GA optimized parameters for nonamers, undecamers, and tridecamers, respectively, are shown in figure 3 (a), (b) and (c). Meanwhile with apparent overestimation of the MP2 calculation, figure 3 (d), (e) and (f) are also showing excellent correlations (R = 0.956, 0.965 and 0.946) of the ΔE between MP2 calculations and AMOEBA with guided GA optimized parameters for nonamers, undecamers, and tridecamers, respectively. There are in total six sets of vdW parameters from nonamers, undecamers, and tridecamers all using GA and guided GA optimization method, respectively.

VDW parameters from Dimers
To validate the transferability of the GA and guided GA optimized AMOEBA parameters from nonamers, undecamers, and tridecamers for gas-phased methanol molecules, we use the six sets of parameters to calibrate the binding energy curves of methanol dimers. We fully optimized 141 methanol dimer configurations at MP2/6-31G (d, p) with BSSE correction with the constraint of pre-defined distance between carbon atoms in the two methanol molecules. Figure 4 is showing the binding energy of methanol dimers as a function of the distance between two carbon atoms.
We obtained the binding energy curves for the 141 configurations of methanol dimer using the three sets of vdW parameters from nonamers, undecamers, and tridecamers using GA optimization method, respectively. Figure 5 (a), (b) and (c) are showing the good agreement of the binding energy between the MP2 calculation and the GA optimized parameters from nonamers, undecamers, and tridecamers, respectively. We also obtained the binding energy curves for those dimers using the other three sets of vdW parameters using guided GA optimization method from nonamers, undecamers, and tridecamers. Figure 5 (d), (e) and (f) are showing the binding energy curves of methanol dimers from the MP2 calculations and from the guided GA optimized vdW parameters by nonamers, undecamers, and tridecamers, respectively. From the lower row of figure 5, we can see discrepancies of binding energy between the gas-phased MP2 calculations and the guided GA optimized parameters from clusters, where the overestimation of MP2 calculations is conspicuous near the lowest energy curve. The shifting down of the binding energy curves from the guided GA optimized parameters implies the interaction energy of the same level of QM theory (i.e., MP2/6-31G(d, p) + BSSE) calculations from gas-phased molecules leads to the overestimation of real potential energy surface (PES). This discrepancy of interaction energy from MP2 calculations is commonly seen as the incompleteness of the basis sets for the particular level of MP2 theory and the insufficiency sampling of configurations in the parameterization of force field. 19,25

VDW parameters for Condensed Phase Methanol
To verify the ability of our parameters of predicting the condensed phase properties, sets of MD simulations using the developed parameters of AMOEBA force field were performed. Here, we applied the six sets of optimized parameters to perform MD simulations to the condensed phase system consisting of 344 methanol molecules. For getting the density (ρ) and heat of vaporization (H vap ), we used Tinker package (omm_dynamics program --GPU version of Tinker dynamics program 27 ) to run MD simulation. NPT ensemble simulations were performed at T = 298.15 K and P = 40 atm with an integration time step of 1 fs. With three different random seeds for initiating velocity on the system, the averaged value of density and heat of vaporization were obtained. Table 1 shows the result of density and heat of vaporization from Tinker MD simulations with the GA optimized AMOEBA force field parameters. Among all the results, the simulated density from tridecamers is the best, while the winner of the heat of vaporization is from undecamers. However, from Table 1, it is clear that the value of density and heat of vaporization are nowhere close to the experimental value.  Table 2 shows the results of density and heat of vaporization from Tinker MD simulations with the guided GA optimized AMOEBA force field parameters. The simulated density and heat of vaporization from tridecamers are much closer to the experimental values. We can see that using the guided GA, which introduces a penalty parameter δ, eliminates the overestimation of interaction energy from MP2 calculation. The guided GA optimized parameters describe better condensed phase properties (i.e., density and heat of vaporization) than the GA optimized parameters in MD simulation. Among all the results, the simulated density and heat of vaporizations from more number of methanol molecule clusters are better. MP2 calculations of interaction energy of tridecamers and undecamers are more descriptive than nonamers for optimizing vdW parameters, due to more detailed information can be calibrated as more molecules in the solvation shell of methanol involved.

CONCLUSIONS
Our study showed that the feasibility to predict experimental condensed phase properties (i.e., density and heat of vaporization) of methanol through problem specific (here, we use AMOEBA) force field from only quantum chemistry information. To acquire the satisfying parameter sets, the genetic algorithm (GA) is the main optimization method. For electrostatic potential energy ( E ESP ), we optimized the both the multipoles and polarizability, Thole damping factors of methanol using the GA method, which leads to lower deviations of E ESP between QM calculations and the GA optimized parameters than that from amoeba09.prm. We optimized the vdW parameters both using GA and guided GA methods by calibrating interaction energy ( ΔE ) of various methanol homo-clusters, such as nonamers, undecamers, or tridecamers. Excellent agreement between the training dataset from MP2 calculations and GA optimized parameters can be achieved. However, only the guided GA method, which eliminates the overestimation of interaction energy from MP2 calculations in the optimization process, provides proper vdW parameters for MD simulation to get condensed phase properties (i.e., density and heat of vaporization) of methanol. Throughout the whole optimization process, the experimental value were not involved in the objective functions, but were only used for the purpose of justifying models (i.e., nonamers, undecamers, or tridecamers) and validating methods (GA or guided GA). We conclude that the main difficulty of parameterizing the force field for liquid phase methanol from solely QM calculations could be coming from the reliability of the training dataset (i.e., MP2 calculations), as it is well known that different basis sets of MP2 could lead to different accuracy of QM calculations, the conspicuous overestimation of interaction energy from MP2 calculations. Some of the other higher level of QM calculations (e.g., MP3/4, CCSD(T), etc.) can only deal with up to certain amount of atoms/electrons. 29 It is ambiguous to study the property of condensed phase properties from a few molecules. In fact, we tried to use the interaction energy ( ΔE ) calculated for methanol tridecamers from another type of QM calculation, Symmetry-Adapted Perturbation Theory (SAPT), as the training dataset using GA optimization, which yields a lousy performance of the condensed phase properties. However, our method was able to utilize the QM calculations (i.e., MP2 calculations) from the limited calculation resource to tune the force field parameters for reasonable agreement between the MD simulated values and experimental data. The method is general and can be extended to any type of small organic molecule systems, as well as other descriptive polarizable force field.