Polarizable empirical force field for sulfur‐containing compounds based on the classical Drude oscillator model

Condensed‐phase computational studies of molecules using molecular mechanics approaches require the use of force fields to describe the energetics of the systems as a function of structure. The advantage of polarizable force fields over nonpolarizable (or additive) models lies in their ability to vary their electronic distribution as a function of the environment. Toward development of a polarizable force field for biological molecules, parameters for a series of sulfur‐containing molecules are presented. Parameter optimization was performed to reproduce quantum mechanical and experimental data for gas phase properties including geometries, conformational energies, vibrational spectra, and dipole moments as well as for condensed phase properties such as heats of vaporization, molecular volumes, and free energies of hydration. Compounds in the training set include methanethiol, ethanethiol, propanethiol, ethyl methyl sulfide, and dimethyl disulfide. The molecular volumes and heats of vaporization are in good accordance with experimental values, with the polarizable model performing better than the CHARMM22 nonpolarizable force field. Improvements with the polarizable model were also obtained for molecular dipole moments and in the treatment of intermolecular interactions as a function of orientation, in part due to the presence of lone pairs and anisotropic atomic polarizability on the sulfur atoms. Significant advantage of the polarizable model was reflected in calculation of the dielectric constants, a property that CHARMM22 systematically underestimates. The ability of this polarizable model to accurately describe a range of gas and condensed phase properties paves the way for more accurate simulation studies of sulfur‐containing molecules including cysteine and methionine residues in proteins. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Systematic benchmarking of large molecular dynamics simulations employing GROMACS on massive multiprocessing facilities

The influence of the total number of cores, the number of cores dedicated to Particle mesh Ewald (PME) calculation and the choice of single vs. double precision on the performance of molecular dynamic (MD) simulations in the size of 70,000 to 1.7 million of atoms was analyzed on three different high‐performance computing facilities employing GROMACS 4 by running about 6000 benchmark simulations. Small and medium sized systems scaled linear up to 64 and 128 cores, respectively. Systems with half a million to 1.2 million atoms scaled linear up to 256 cores. The best performance was achieved by dedicating 25% of the total number of cores to PME calculation. Double precision calculations lowered the performance by 30–50%. A database for collecting information about MD simulations and the achieved performance was created and is freely available online and allows the fast estimation of the performance that can be expected in similar environments. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

A polarizable,charge transfer model of water using the drude oscillator

The transfer of small amounts of charge between neighboring particles can be a significant part of interactions among particles. A model is developed for treating charge transfer (CT) combined with the Drude model for polarizability to create an efficient model for liquid water which includes both CT and polarizability. The model is shown to be accurate for a variety of liquid properties, including the density as a function of temperature and the dielectric constant. A new model for water with CT and polarization is developed and applied to the liquid. The inclusion of CT increases the accuracy of many properties, like the density as a function of temperature, indicating the importance of charge redistribution as induced by other particles. © 2016 Wiley Periodicals, Inc.     

High-performance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD

Incorporating the influence of induced polarization in large-scale atomistic molecular dynamics (MD) simulations is a critical challenge in the progress toward computations of increased accuracy. One computationally efficient treatment is based on the classical Drude oscillator in which an auxiliary charged particle is attached by a spring to each nucleus. Here, we report the first implementation of this model in the program NAMD. An extended Lagrangian dynamics with a dual-Langevin thermostat scheme applied to the Drude-nucleus pairs is employed to efficiently generate classical dynamic propagation near the self-consistent field limit. Large-scale MD simulations based on the Drude polarizable force field scale very well on massively distributed supercomputing platforms, the computational demand being only about 50-100% higher than for nonpolarizable models. As an illustration, a large-scale 150 mM NaCl aqueous salt solution is simulated, and the calculated ionic conductivity is shown to be in excellent agreement with experiment.     

Best bang for your buck: GPU nodes for GROMACS biomolecular simulations

The molecular dynamics simulation package GROMACS runs efficiently on a wide variety of hardware from commodity workstations to high performance computing clusters. Hardware features are well‐exploited with a combination of single instruction multiple data, multithreading, and message passing interface (MPI)‐based single program multiple data/multiple program multiple data parallelism while graphics processing units (GPUs) can be used as accelerators to compute interactions off‐loaded from the CPU. Here, we evaluate which hardware produces trajectories with GROMACS 4.6 or 5.0 in the most economical way. We have assembled and benchmarked compute nodes with various CPU/GPU combinations to identify optimal compositions in terms of raw trajectory production rate, performance‐to‐price ratio, energy efficiency, and several other criteria. Although hardware prices are naturally subject to trends and fluctuations, general tendencies are clearly visible. Adding any type of GPU significantly boosts a node's simulation performance. For inexpensive consumer‐class GPUs this improvement equally reflects in the performance‐to‐price ratio. Although memory issues in consumer‐class GPUs could pass unnoticed as these cards do not support error checking and correction memory, unreliable GPUs can be sorted out with memory checking tools. Apart from the obvious determinants for cost‐efficiency like hardware expenses and raw performance, the energy consumption of a node is a major cost factor. Over the typical hardware lifetime until replacement of a few years, the costs for electrical power and cooling can become larger than the costs of the hardware itself. Taking that into account, nodes with a well‐balanced ratio of CPU and consumer‐class GPU resources produce the maximum amount of GROMACS trajectory over their lifetime. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation,validation, and benchmarks

Presented is the implementation of the Drude force field in the open‐source OpenMM simulation package allowing for access to graphical processing unit (GPU) hardware. In the Drude model, electronic degrees of freedom are represented by negatively charged particles attached to their parent atoms via harmonic springs, such that extra computational overhead comes from these additional particles and virtual sites representing lone pairs on electronegative atoms, as well as the associated thermostat and integration algorithms. This leads to an approximately fourfold increase in computational demand over additive force fields. However, by making the Drude model accessible to consumer‐grade desktop GPU hardware it will be possible to perform simulations of one microsecond or more in less than a month, indicating that the barrier to employ polarizable models has largely been removed such that polarizable simulations with the classical Drude model are readily accessible and practical.     

MkVsites: A tool for creating GROMACS virtual sites parameters to increase performance in all-atom molecular dynamics simulations

The absolute performance of any all-atom molecular dynamics simulation is typically limited by the length of the individual timesteps taken when integrating the equations of motion. In the GROMACS simulation software, it has for a long time been possible to use so-called virtual sites to increase the length of the timestep, resulting in a large gain of simulation efficiency. Up until now, support for this approach has in practice been limited to the standard 20 amino acids however, shrinking the applicability domain of virtual sites. MkVsites is a set of python tools which provides a convenient way to obtain all parameters necessary to use virtual sites for virtually any molecules in a simulation. Required as input to MkVsites is the molecular topology of the molecule(s) in question, along with a specification of where to find the parent force field. As such, MkVsites can be a very valuable tool suite for anyone who is routinely using GROMACS for the simulation of molecular systems.     

A polarizable empirical force field for molecular dynamics simulation of liquid hydrocarbons

Electronic polarizability is usually treated implicitly in molecular simulations, which may lead to imprecise or even erroneous molecular behavior in spatially electronically inhomogeneous regions of systems such as proteins, membranes, interfaces between compounds, or mixtures of solvents. The majority of available molecular force fields and molecular dynamics simulation software packages does not account explicitly for electronic polarization. Even the simplest charge‐on‐spring (COS) models have only been developed for few types of molecules. In this work, we report a polarizable COS model for cyclohexane, as this molecule is a widely used solvent, and for linear alkanes, which are also used as solvents, and are the precursors of lipids, amino acid side chains, carbohydrates, or nucleic acid backbones. The model is an extension of a nonpolarizable united‐atom model for alkanes that had been calibrated against experimental values of the density, the heat of vaporization and the Gibbs free energy of hydration for each alkane. The latter quantity was used to calibrate the parameters governing the interaction of the polarizable alkanes with water. Subsequently, the model was tested for other structural, thermodynamic, dielectric, and dynamic properties such as trans/gauche ratios, excess free energy, static dielectric permittivity, and self‐diffusion. A good agreement with the experimental data for a large set of properties for each considered system was obtained, resulting in a transferable set of polarizable force‐field parameters for CH₂, CH₃, and CH₄ moieties. © 2014 Wiley Periodicals, Inc.     

Time-reversible always stable predictor-corrector method for molecular dynamics of polarizable molecules

An improved method for classic molecular dynamics of polarizable molecules is proposed. The method uses a predictor, one evaluation of the electrostatic field per integration step, and relaxation (damping). The self-consistent solution is approximated with error of the second order (with respect to the timestep). The time reversibility (long-time energy conservation) error is of the (2n - 1)th order, where n is the predictor length. The method is easy to implement, efficient, accurate, and suitable for any model of polarizability.     

Explicit proton transfer in classical molecular dynamics simulations

We present Hydrogen Dynamics (HYDYN), a method that allows explicit proton transfer in classical force field molecular dynamics simulations at thermodynamic equilibrium. HYDYN reproduces the characteristic properties of the excess proton in water, from the special pair dance, to the continuous fluctuation between the limiting Eigen and Zundel complexes, and the water reorientation beyond the first solvation layer. Advantages of HYDYN with respect to existing methods are computational efficiency, microscopic reversibility, and easy parameterization for any force field. © 2014 Wiley Periodicals, Inc.     

Optimization of parameters for molecular dynamics simulation using smooth particle-mesh Ewald in GROMACS 4.5

Based on our critique of requirements for performing an efficient molecular dynamics simulation with the particle-mesh Ewald (PME) implementation in GROMACS 4.5, we present a computational tool to enable the discovery of parameters that produce a given accuracy in the PME approximation of the full electrostatics. Calculations on two parallel computers with different processor and communication structures showed that a given accuracy can be attained over a range of parameter space, and that the attributes of the hardware and simulation system control which parameter sets are optimal. This information can be used to find the fastest available PME parameter sets that achieve a given accuracy. We hope that this tool will stimulate future work to assess the impact of the quality of the PME approximation on simulation outcomes, particularly with regard to the trade-off between cost and scientific reliability in biomolecular applications.     

FFParam: Standalone package for CHARMM additive and Drude polarizable force field parametrization of small molecules

Accurate force-field (FF) parameters are key to reliable prediction of properties obtained from molecular modeling (MM) and molecular dynamics (MD) simulations. With ever-widening applicability of MD simulations, robust parameters need to be generated for a wider range of chemical species. The CHARMM General Force Field program (CGenFF, https://cgenff.umaryland.edu/ ) is a tool for obtaining initial parameters for a given small molecule based on analogy with the available CGenFF parameters. However, improvement of these parameters is often required and performing their optimization remains tedious and time consuming. In addition, tools for optimization of small molecule parameters in the context of the Drude polarizable FF are not yet available. To overcome these issues, the FFParam package has been designed to facilitate the parametrization process. The package includes a graphical user interface (GUI) created using Qt libraries. FFParam supports Gaussian and Psi4 for performing quantum mechanical calculations and CHARMM and OpenMM for MM calculations. A Monte Carlo simulated annealing (MCSA) algorithm has been implemented for automated fitting of partial atomic charge, atomic polarizabilities and Thole scale parameters. The LSFITPAR program is called for automated fitting of bonded parameters. Accordingly, FFParam provides all the features required for generation and analysis of CHARMM and Drude FF parameters for small molecules. FFParam-GUI includes a text editor, graph plotter, molecular visualization, and text to table converter to meet various requirements of the parametrization process. It is anticipated that FFParam will facilitate wider use of CGenFF as well as promote future use of the Drude polarizable FF.     

All‐atom polarizable force field for DNA based on the classical drude oscillator model

Presented is a first generation atomistic force field (FF) for DNA in which electronic polarization is modeled based on the classical Drude oscillator formalism. The DNA model is based on parameters for small molecules representative of nucleic acids, including alkanes, ethers, dimethylphosphate, and the nucleic acid bases and empirical adjustment of key dihedral parameters associated with the phosphodiester backbone, glycosidic linkages, and sugar moiety of DNA. Our optimization strategy is based on achieving a compromise between satisfying the properties of the underlying model compounds in the gas phase targeting quantum mechanical (QM) data and reproducing a number of experimental properties of DNA duplexes in the condensed phase. The resulting Drude FF yields stable DNA duplexes on the 100‐ns time scale and satisfactorily reproduce (1) the equilibrium between A and B forms of DNA and (2) transitions between the BI and BII substates of B form DNA. Consistency with the gas phase QM data for the model compounds is significantly better for the Drude model as compared to the CHARMM36 additive FF, which is suggested to be due to the improved response of the model to changes in the environment associated with the explicit inclusion of polarizability. Analysis of dipole moments associated with the nucleic acid bases shows the Drude model to have significantly larger values than those present in CHARMM36, with the dipoles of individual bases undergoing significant variations during the MD simulations. Additionally, the dipole moment of water was observed to be perturbed in the grooves of DNA. © 2014 Wiley Periodicals, Inc.     

Polarizable empirical force field for nitrogen‐containing heteroaromatic compounds based on the classical Drude oscillator

The polarizable empirical CHARMM force field based on the classical Drude oscillator has been extended to the nitrogen‐containing heteroaromatic compounds pyridine, pyrimidine, pyrrole, imidazole, indole, and purine. Initial parameters for the six‐membered rings were based on benzene with nonbond parameter optimization focused on the nitrogen atoms and adjacent carbons and attached hydrogens. In the case of five‐member rings, parameters were first developed for imidazole and transferred to pyrrole. Optimization of all parameters was performed against an extensive set of quantum mechanical and experimental data. Ab initio data were used for the determination of initial electrostatic parameters, the vibrational analysis, and in the optimization of the relative magnitudes of the Lennard‐Jones (LJ) parameters, through computations of the interactions of dimers of model compounds, model compound‐water interactions, and interactions of rare gases with model compounds. The absolute values of the LJ parameters were determined targeting experimental heats of vaporization, molecular volumes, heats of sublimation, crystal lattice parameters, and free energies of hydration. Final scaling of the polarizabilities from the gas‐phase values by 0.85 was determined by reproduction of the dielectric constants of pyridine and pyrrole. The developed parameter set was extensively validated against additional experimental data such as diffusion constants, heat capacities, and isothermal compressibilities, including data as a function of temperature. © 2008 Wiley Periodicals, Inc. J Comput Chem 2009     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.     

Charge distributions for molecular dynamics simulations from self-consistent polarization method

Partial atomic charges are important force field parameters. They are usually computed by applying quantum-chemical calculations and the assumed population scheme. In this study polarization consistent scheme of deriving a charge distribution inside solute molecule is proposed. The environment effect is explicitly taken into account by distributing solvent molecules around the solute target. The performed analysis includes a few computational schemes (HF, MP2, B3LYP, and M026X), basis sets (cc-pvnz, n = 2, 3, …, 6), and electrostatically derived charge distributions (KS, CHELP, CHELPG, and HLY). It is demonstrated that the environment effect is very important and cannot be disregarded. The second solvation shell should be included to achieve the charge convergence. Huge corrections to charge distribution are due to induction and dispersion. The B3LYP/cc-pvqz level of theory is recommended for deriving the charges within self-consistent polarization scheme.     

Fluctuating charge normal modes: An algorithm for implementing molecular dynamics simulations with polarizable potentials

A new method for performing molecular dynamics simulations with fluctuating charge polarizable potentials is introduced. In fluctuating charge models, polarizability is treated by allowing the partial charges to be variables, with values that are coupled to charges on the same molecule as well as those on other molecules. The charges can be efficiently propagated in a molecular dynamics simulation using extended Lagrangian dynamics. By making a coordinate change from the charge variables to a set of normal mode charge coordinates for each molecule, a new method is constructed in which the normal mode charge variables uncouple from those on the same molecule. The method is applied to the TIP4P-FQ model of water and compared to other methods for implementing the dynamics. The methods are compared using different molecular dynamics time steps.     

Molecular dynamics simulations of liquid methanol and methanol-water mixtures with polarizable models

A polarizable model for simulation of liquid methanol, compatible with the COS/G2 water model, has been developed using the Charge-on-Spring (COS) technique. The model consists of three point charges, with one polarizable center on the oxygen atom. The Lennard-Jones parameters on the oxygen atom together with the molecular polarizability were varied to reproduce the experimental heat of vaporization and density of liquid methanol at ambient conditions. We examined the energies of various methanol dimers in the gas phase and compared them with values obtained from ab initio calculations. The model was then used to study the thermodynamic, dynamic, structural, and dielectric properties of liquid methanol as well as of a methanol-water mixture. A microscopic picture of the structure of pure liquid methanol and of the methanol-water mixture is provided. Good agreement was found between the results from our model simulations and available experimental and ab initio calculation data. In particular, the experimental dielectric permittivity of 32 could be reproduced, which had been shown to be difficult when using nonpolarizable models.     

The distributed diagonal force decomposition method for parallelizing molecular dynamics simulations

Parallelization is an effective way to reduce the computational time needed for molecular dynamics simulations. We describe a new parallelization method, the distributed-diagonal force decomposition method, with which we extend and improve the existing force decomposition methods. Our new method requires less data communication during molecular dynamics simulations than replicated data and current force decomposition methods, increasing the parallel efficiency. It also dynamically load-balances the processors' computational load throughout the simulation. The method is readily implemented in existing molecular dynamics codes and it has been incorporated into the CHARMM program, allowing its immediate use in conjunction with the many molecular dynamics simulation techniques that are already present in the program. We also present the design of the Force Decomposition Machine, a cluster of personal computers and networks that is tailored to running molecular dynamics simulations using the distributed diagonal force decomposition method. The design is expandable and provides various degrees of fault resilience. This approach is easily adaptable to computers with Graphics Processing Units because it is independent of the processor type being used.