Recent applications and developments of charge equilibration force fields for modeling dynamical charges in classical molecular dynamics simulations

Unifying Charge-Flow Polarization Models

We show that several models where electric polarization in molecular systems is modeled by charge-flow between atoms can all be considered as different manifestations of a general underlying mathematical structure. The models can be classified according to whether they employ atomic or bond parameters and whether they employ atom/bond hardness or softness. We show that an ab initio calculated charge response kernel can be considered as the inverse screened Coulombic matrix projected onto the zero-charge subspace, and this may provide a method for deriving charge screening functions to be used in force fields. The analysis suggests that some models contain redundancies, and we argue that a parameterization of charge-flow models in terms of bond softness is preferable as it depends on local quantities and decay to zero upon bond dissociation, while bond hardness depends on global quantities and increases toward infinity upon bond dissociation.

Harnessing Deep Learning for Optimization of Lennard-Jones Parameters for the Polarizable Classical Drude Oscillator Force Field

The outcomes of computational chemistry and biology research, including drug design, are significantly influenced by the underlying force field (FF) used in molecular simulations. While improved FF accuracy may be achieved via inclusion of explicit treatment of electronic polarization, such an extension must be accompanied by optimization of van der Waals (vdW) interactions, in the context of the Lennard-Jones (LJ) formalism in the present study. This is particularly challenging due to the extensive nature of chemical space combined with the correlated nature of LJ parameters. To address this challenge, a deep learning (DL)-based parametrization framework is developed, allowing for sampling of wide ranges of LJ parameters targeting experimental condensed phase thermodynamic properties. The present work utilizes this framework to develop the LJ parameters for atoms associated with four distinct groups covering 10 different atom types. Final parameter selection was facilitated by quantum mechanical data on rare-gas interactions with the training set molecules. The chosen parameters were then validated through experimental hydration free energies and condensed phase thermodynamic properties of validation set molecules to confirm transferability. The ultimate outcome of utilizing this framework is a set of LJ parameters in the context of the polarizable Drude FF, which demonstrated improvement in the reproduction of both experimental pure solvent and crystal properties and hydration free energies of the molecules compared to the additive CHARMM General FF (CGenFF) including the ability of the Drude FF to accurately reproduce both experimental pure solvent properties and hydration free energies. The study also shows how correlations between difference in the reproduction of condensed phase data between model compounds may be used to direct the selection of new atom types and training set molecules during FF development.

Singular value decomposition analysis of the electron density changes occurring upon electrostatic polarization of water

In-depth elucidation of how molecules are electrically polarized would be one key factor for understanding the properties of those molecules under various thermodynamic and/or spatial conditions. Here this problem is tackled for the case of hydrogen-bonded water by conducting singular value decomposition of the electron density changes that occur upon electrostatic polarization. It is shown that all those electron density changes are approximately described as linear combinations of ten orthonormal basis “vectors”. One main component is the interatomic charge transfer through each OH bond, while some others are characterized as the atomic dipolar polarizations, meaning that both of these components are important for the electrostatic polarization of water. The interaction parameters that reasonably well reproduce the induced dipole moments are derived, which indicate the extent of mixing of the two components in electrostatic polarization.

The main features of the electron density changes that occur upon electrostatic polarization of water are elucidated by conducting singular value decomposition analysis of those changes.

Dividing the Periodic Box into Subdivisions with Their Surroundings to Accelerate Molecular Dynamics Simulation with High Accuracy

A major field of current research in chemistry and biology is the development of the tools that enable in situ analysis of complex systems. However, the long-time dynamics simulation for an extremely large system in solution is almost impossible by an all-atom force field combined with an explicit solvent model. The results show that the larger the periodic box is, the closer the properties of the system are to the experimental values. Therefore, how can we carry out simulations for systems that are fast, accurate, and large enough? A method of dividing the periodic box into subdivisions with their surroundings (DBSS) is presented here, and it clearly increases the computation speed without losing accuracy and enables the simulation of extremely large systems by strongly decreasing the dimension of the charge matrix. The DBSS method divides a single periodic box or unit in an extremely large system into several subdivisions with a suitable choice according to atomic coordinates. This method ensures that these subdivisions are always changing and allows the atoms to communicate with each other. Intermolecular communication is important for molecular properties and functions but is not possible with other fragment methods. The partial charges are calculated in each subdivision with an overlapping surrounding used to take hydrogen bond interaction between the subdivisions into account. This is an iterative process because the charge population will be recalculated at intervals during the dynamics simulations. Taking a water system as an example, each subdivision is extended by 4 Å as the surrounding. The computation time scales almost linearly with the size of the system, and the slope is small. MD simulations for several properties have been performed by the ABEEM-DBSS method. The results indicate that the ABEEM-DBSS method can accurately simulate the properties of water system, and the accuracy can reach or approach that of the experimental data or of other water potentials. Interestingly, the properties become closer to the experimental data as the sizes of the periodic box increase, further validating the need for the simulation of a large system and demonstrating the value of the DBSS method.

Stochastic Constrained Extended System Dynamics for Solving Charge Equilibration Models

We present a new stochastic extended Lagrangian molecular dynamics solution to charge equilibration that eliminates self-consistent field (SCF) calculations, thus eliminating the computational bottleneck in solving the charge distribution with standard SCF solvers. By formulating both charges and chemical potential as latent variables and introducing a holonomic constraint that satisfies charge conservation, the SC-XLMD method accurately reproduces thermodynamic, dynamic, and structural properties within the framework of ReaxFF for a bulk water system and highly reactive RDX molecules simulated at high temperature. The SC-XLMD method shows excellent computational performance and is available in the publicly available LAMMPS package.

A fast and high-quality charge model for the next generation general AMBER force field

The General AMBER Force Field (GAFF) has been broadly used by researchers all over the world to perform in silico simulations and modelings on diverse scientific topics, especially in the field of computer-aided drug design whose primary task is to accurately predict the affinity and selectivity of receptor-ligand binding. The atomic partial charges in GAFF and the second generation of GAFF (GAFF2) were originally developed with the quantum mechanics derived restrained electrostatic potential charge, but in practice, users usually adopt an efficient charge method, Austin Model 1-bond charge corrections (AM1-BCC), based on which, without expensive ab initio calculations, the atomic charges could be efficiently and conveniently obtained with the ANTECHAMBER module implemented in the AMBER software package. In this work, we developed a new set of BCC parameters specifically for GAFF2 using 442 neutral organic solutes covering diverse functional groups in aqueous solution. Compared to the original BCC parameter set, the new parameter set significantly reduced the mean unsigned error (MUE) of hydration free energies from 1.03 kcal/mol to 0.37 kcal/mol. More excitingly, this new AM1-BCC model also showed excellent performance in the solvation free energy (SFE) calculation on diverse solutes in various organic solvents across a range of different dielectric constants. In this large-scale test with totally 895 neutral organic solvent-solute systems, the new parameter set led to accurate SFE predictions with the MUE and the root-mean-square-error of 0.51 kcal/mol and 0.65 kcal/mol, respectively. This newly developed charge model, ABCG2, paved a promising path for the next generation GAFF development.

Impact of electronic polarizability on protein-functional group interactions

Interactions of proteins with functional groups are key to their biological functions, making it essential that they be accurately modeled. To investigate the impact of the inclusion of explicit treatment of electronic polarizability in force fields on protein-functional group interactions, the additive CHARMM and Drude polarizable force field are compared in the context of the Site-Identification by Ligand Competitive Saturation (SILCS) simulation methodology from which functional group interaction patterns with five proteins for which experimental binding affinities of multiple ligands are available, were obtained. The explicit treatment of polarizability produces significant differences in the functional group interactions in the ligand binding sites including overall enhanced binding of functional groups to the proteins. This is associated with variations of the dipole moments of solutes representative of functional groups in the binding sites relative to aqueous solution with higher dipole moments systematically occurring in the latter, though exceptions occur with positively charged methylammonium. Such variation indicates the complex, heterogeneous nature of the electronic environments of ligand binding sites and emphasizes the inherent limitation of fixed charged, additive force fields for modeling ligand-protein interactions. These effects yield more defined orientation of the functional groups in the binding pockets and a small, but systematic improvement in the ability of the SILCS method to predict the binding orientation and relative affinities of ligands to their target proteins. Overall, these results indicate that the physical model associated with the explicit treatment of polarizability along with the presence of lone pairs in a force field leads to changes in the nature of the interactions of functional groups with proteins versus that occurring with additive force fields, suggesting the utility of polarizable force fields in obtaining a more realistic understanding of protein-ligand interactions.

Polarizable molecular dynamics simulations of ionic liquids: Influence of temperature control

Ionic liquids are an interesting class of soft matter with viscosities of one or two orders of magnitude higher than that of water. Unfortunately, classical, non-polarizable molecular dynamics (MD) simulations of ionic liquids result in too slow dynamics and demonstrate the need for explicit inclusion of polarizability. The inclusion of polarizability, here via the Drude oscillator model, requires amendments to the employed thermostat, where we consider a dual Nosé-Hoover thermostat, as well as a dual Langevin thermostat. We investigate the effects of the choice of a thermostat and the underlying parameters such as the masses and force constants of the Drude particles on static and dynamic properties of ionic liquids. Here, we show that Langevin thermostats are not suitable for investigating the dynamics of ionic liquids. Since polarizable MD simulations are associated with high computational costs, we employed a self-developed graphics processing unit enhanced code within the MD program CHARMM to keep the overall computational effort reasonable.

Describing Molecular Polarizability by a Bond Capacity Model

We propose a bond capacity model for describing molecular polarization in force field energy functions at the charge-only level. Atomic charges are calculated by allowing charge to flow between atom pairs according to a bond capacity and a difference in electrostatic potential. The bond capacity is closely related to the bond order and decays to zero as the bond distance is increased. The electrostatic potential is composed of an intrinsic potential, identified as the electronegativity, and a screened Coulomb potential from all other charges. The bond capacity model leads to integer fragment charges upon bond dissociation and displays linear scaling of the polarizability with system size. Bond capacity parameters can be derived from reference molecular polarizabilities, while electronegativity parameters can be derived from reference atomic charges or a reference molecular electrostatic potential. Out-of-plane polarization for planar systems is modeled by off-nuclei charge sites. The model is shown to be able to reproduce anisotropic reference molecular polarizabilities with an accuracy of ∼10% using a limited set of bond capacity parameters and can describe both inter- and intramolecular polarization.

Toward Prediction of Electrostatic Parameters for Force Fields That Explicitly Treat Electronic Polarization

The derivation of atomic polarizabilities for polarizable force field development has been a long-standing problem. Atomic polarizabilities were often refined manually starting from tabulated values, rendering an automated assignment of parameters difficult and hampering reproducibility and transferability of the obtained values. To overcome this, we trained both a linear increment scheme and a multilayer perceptron neural network on a large number of high-quality quantum mechanical atomic polarizabilities and partial atomic charges, where only the type of each atom and its connectivity were used as input. The predicted atomic polarizabilities and charges had average errors of 0.023 ų and 0.019 e using the neural net and 0.063 ų and 0.069 e using the simple increment scheme. As the algorithm relies only on the connectivities of the atoms within a molecule, thus omitting dependencies on the three-dimensional conformation, the approach naturally assigns like charges and polarizabilities to symmetrical groups. Accordingly, a convenient utility is presented for generating the partial atomic charges and atomic polarizabilities for organic molecules as needed in polarizable force field development.

Inertial extended-Lagrangian scheme for solving charge equilibration models

The inertial EL/SCF method is developed to solve charge equilibration models for molecular dynamics, reducing the number of SCFs by 50–80% at each time step when compared to a conjugate gradient SCF solver and tested on diverse reactive systems.

On flexible force fields for metal-organic frameworks: Recent developments and future prospects

Classical force field simulations can be used to study structural, diffusion, and adsorption properties of metal-organic frameworks (MOFs). To account for the dynamic behavior of the material, parameterization schemes have been developed to derive force constants and the associated reference values by fitting on ab initio energies, vibrational frequencies, and elastic constants. Here, we review recent developments in flexible force field models for MOFs. Existing flexible force field models are generally able to reproduce the majority of experimentally observed structural and dynamic properties of MOFs. The lack of efficient sampling schemes for capturing stimuli-driven phase transitions, however, currently limits the full predictive potential of existing flexible force fields from being realized. This article is categorized under: Structure and Mechanism > Computational Materials ScienceMolecular and Statistical Mechanics > Molecular Mechanics.

Machine Learning Force Field Parameters from Ab Initio Data

Machine learning (ML) techniques with the genetic algorithm (GA) have been applied to determine a polarizable force field parameters using only ab initio data from quantum mechanics (QM) calculations of molecular clusters at the MP2/6-31G(d,p), DFMP2(fc)/jul-cc-pVDZ, and DFMP2(fc)/jul-cc-pVTZ levels to predict experimental condensed phase properties (i.e., density and heat of vaporization). The performance of this ML/GA approach is demonstrated on 4943 dimer electrostatic potentials and 1250 cluster interaction energies for methanol. Excellent agreement between the training data set from QM calculations and the optimized force field model was achieved. The results were further improved by introducing an offset factor during the machine learning process to compensate for the discrepancy between the QM calculated energy and the energy reproduced by optimized force field, while maintaining the local "shape" of the QM energy surface. Throughout the machine learning process, experimental observables were not involved in the objective function, but were only used for model validation. The best model, optimized from the QM data at the DFMP2(fc)/jul-cc-pVTZ level, appears to perform even better than the original AMOEBA force field (amoeba09.prm), which was optimized empirically to match liquid properties. The present effort shows the possibility of using machine learning techniques to develop descriptive polarizable force field using only QM data. The ML/GA strategy to optimize force fields parameters described here could easily be extended to other molecular systems.

Force Field Development for Lipid Membrane Simulations

With the rapid development of computer power and wide availability of modelling software computer simulations of realistic models of lipid membranes, including their interactions with various molecular species, polypeptides and membrane proteins have become feasible for many research groups. The crucial issue of the reliability of such simulations is the quality of the force field, and many efforts, especially in the latest several years, have been devoted to parametrization and optimization of the force fields for biomembrane modelling. In this review, we give account of the recent development in this area, covering different classes of force fields, principles of the force field parametrization, comparison of the force fields, and their experimental validation. This article is part of a Special Issue entitled: Biosimulations edited by Ilpo Vattulainen and Tomasz Róg.

Refined OPLS all-atom force field for saturated phosphatidylcholine bilayers at full hydration

We report parametrization of dipalmitoyl-phosphatidylcholine (DPPC) in the framework of the Optimized Parameters for Liquid Simulations all-atom (OPLS-AA) force field. We chose DPPC as it is one of the most studied phospholipid species and thus has plenty of experimental data necessary for model validation, and it is also one of the highly important and abundant lipid types, e.g., in lung surfactant. Overall, PCs have not been previously parametrized in the OPLS-AA force field; thus, there is a need to derive its bonding and nonbonding parameters for both the polar and nonpolar parts of the molecule. In the present study, we determined the parameters for torsion angles in the phosphatidylcholine and glycerol moieties and in the acyl chains, as well the partial atomic charges. In these calculations, we used three methods: (1) Hartree-Fock (HF), (2) second order Møller-Plesset perturbation theory (MP2), and (3) density functional theory (DFT). We also tested the effect of the polar environment by using the polarizable continuum model (PCM), and for acyl chains the van der Waals parameters were also adjusted. In effect, six parameter sets were generated and tested on a DPPC bilayer. Out of these six sets, only one was found to be able to satisfactorily reproduce experimental data for the lipid bilayer. The successful DPPC model was obtained from MP2 calculations in an implicit polar environment (PCM).

Polarizable empirical force field for hexopyranose monosaccharides based on the classical Drude oscillator

A polarizable empirical force field based on the classical Drude oscillator is presented for the hexopyranose form of selected monosaccharides. Parameter optimization targeted quantum mechanical (QM) dipole moments, solute-water interaction energies, vibrational frequencies, and conformational energies. Validation of the model was based on experimental data on crystals, densities of aqueous-sugar solutions, diffusion constants of glucose, and rotational preferences of the exocylic hydroxymethyl of d-glucose and d-galactose in aqueous solution as well as additional QM data. Notably, the final model involves a single electrostatic model for all sixteen diastereomers of the monosaccharides, indicating the transferability of the polarizable model. The presented parameters are anticipated to lay the foundation for a comprehensive polarizable force field for saccharides that will be compatible with the polarizable Drude parameters for lipids and proteins, allowing for simulations of glycolipids and glycoproteins.

The Polarizable Atomic Multipole-based AMOEBA Force Field for Proteins

Development of the AMOEBA (Atomic Multipole Optimized Energetics for Biomolecular Simulation) force field for proteins is presented. The current version (AMOEBA-2013) utilizes permanent electrostatic multipole moments through the quadrupole at each atom, and explicitly treats polarization effects in various chemical and physical environments. The atomic multipole electrostatic parameters for each amino acid residue type are derived from high-level gas phase quantum mechanical calculations via a consistent and extensible protocol. Molecular polarizability is modeled via a Thole-style damped interactive induction model based upon distributed atomic polarizabilities. Inter- and intramolecular polarization is treated in a consistent fashion via the Thole model. The intramolecular polarization model ensures transferability of electrostatic parameters among different conformations, as demonstrated by the agreement between QM and AMOEBA electrostatic potentials, and dipole moments of dipeptides. The backbone and side chain torsional parameters were determined by comparing to gas-phase QM (RI-TRIM MP2/CBS) conformational energies of dipeptides and to statistical distributions from the Protein Data Bank. Molecular dynamics simulations are reported for short peptides in explicit water to examine their conformational properties in solution. Overall the calculated conformational free energies and J-coupling constants are consistent with PDB statistics and experimental NMR results, respectively. In addition, the experimental crystal structures of a number of proteins are well maintained during molecular dynamics (MD) simulation. While further calculations are necessary to fully validate the force field, initial results suggest the AMOEBA polarizable multipole force field is able to describe the structure and energetics of peptides and proteins, in both gas-phase and solution environments.