Towards exact molecular dynamics simulations with machine-learned force fields

We develop a combined machine learning (ML) and quantum mechanics approach that enables data-efficient reconstruction of flexible molecular force fields from high-level ab initio calculations, through the consideration of fundamental physical constraints. We discuss how such constraints are recovered and incorporated into ML models. Specifically, we use conservation of energy - a fundamental property of closed classical and quantum mechanical systems -- to derive an efficient gradient-domain machine learning (GDML) model. The challenge of constructing conservative force fields is accomplished by learning in a Hilbert space of vector-valued functions that obey the law of energy conservation. We proceed with the development of a multi-partite matching algorithm that enables a fully automated recovery of physically relevant point-group and fluxional symmetries from the training dataset into a symmetric variant of our model. The developed symmetric GDML (sGDML) approach is able to faithfully reproduce global force fields at the accuracy high-level ab initio methods, thus enabling sample intensive tasks like molecular dynamics simulations at that level of accuracy.

Machine learning (ML) force fields are revolutionizing molecular dynamics (MD) simulations as they bypass the computational cost associated with ab initio methods but do not sacrifice accuracy in the process. In this work, the GPyTorch library is used to create Gaussian process regression (GPR) models that are interfaced with the next-generation ML force field FFLUX. These models predict atomic properties of different molecular configurations that appear in a progressing MD simulation. An improved kernel function is utilized to correctly capture the periodicity of the input descriptors. The first FFLUX molecular simulations of ammonia, methanol, and malondialdehyde with the updated kernel are performed. Geometry optimizations with the GPR models result in highly accurate final structures with a maximum root-mean-squared deviation of 0.064 Å and sub-kJ mol-1 total energy predictions. Additionally, the models are tested in 298 K MD simulations with FFLUX to benchmark for robustness. The resulting energy and force predictions throughout the simulation are in excellent agreement with ab initio data for ammonia and methanol but decrease in quality for malondialdehyde due to the increased system complexity. GPR model improvements are discussed, which will ensure the future scalability to larger systems.

Stable molecular dynamics simulations of halide perovskites from a temperature-ensemble gradient-domain machine learning approach

MDWiZ: A platform for the automated translation of molecular dynamics simulations

Machine-learning (ML) techniques have drawn an ever-increasing focus as they enable high-throughput screening and multiscale prediction of material properties. Especially, ML force fields (FFs) of quantum mechanical accuracy are expected to play a central role for the purpose. The construction of ML-FFs for polymers is, however, still in its infancy due to the formidable configurational space of its composing atoms. Here, we demonstrate the effective development of ML-FFs using kernel functions and a Gaussian process for an organic polymer, polytetrafluoroethylene (PTFE), with a data set acquired by first-principles calculations and ab initio molecular dynamics (AIMD) simulations. Even though the training data set is sampled only with short PTFE chains, structures of longer chains optimized by our ML-FF show an excellent consistency with density functional theory calculations. Furthermore, when integrated with molecular dynamics simulations, the ML-FF successfully describes various physical properties of a PTFE bundle, such as a density, melting temperature, coefficient of thermal expansion, and Young's modulus.

Source rocks contain significant volumes of hydrocarbon fluids trapped in kerogen, however their effective recovery is challenged due to amplified fluid-wall interactions and the nanopore confinement impact on fluid composition. Enhanced oil production can be achieved by modifying the existing molecular forces in kerogen pore-network by using custom-designed targeted chemistry technologies. Our objective is to show how the transport of hydrocarbons in kerogen and its recovery can be altered with the delivery of microemulsion nanodroplets into the pore network. This is done by using computational chemistry and molecular dynamics (MD) simulations. Molecular dynamics simulation is used to generate a 3D model replica of a Type II kerogen representative of source rocks located in Delaware and Midland basins in the United States. Oil phase saturated kerogen is modeled as consisting of nine different types of molecules: dimethyl naphthalene, toluene, tetradecane, decane, octane, butane, propane, ethane and methane. The delivered microemulsion is an aqueous dispersion of solvent-swollen surfactant micelles. The solvent and nonionic surfactant present in the microemulsion are modeled as d-limonene and dodecanol heptaethyl ether (C12E7), respectively. Molecular dynamics simulation experiments include two steps: (i) the injection of microemulsion treatment fluid into the oil-saturated kerogen pore-network, and (ii) transient flow-back of the oil-chemical mixture in the pore network. The utilized 3D kerogen models were developed based on a representative oil sample composition (H, C, O, S, N) from the region. Simulation results show that microemulsions can affect the reservoir via two different mechanisms. During the injection, microemulsion nanodroplets that enter the nano-capillaries of pore network disperse in the liquid present in the pore space under the influence of pore walls. The solvent dissolves in the oil phase and alters the physical and transport properties of the phase, while the surfactant molecules modify the wettability of the solid kerogen surfaces. The recovery effectiveness of heavier oil fractions is improved compared to the recovery effectiveness achieved with surfactant micelles without the solubilized solvent. New 3D kerogen models are presented using atomistic modeling and molecular simulations. These models possess important chemical and physical characteristics of the organic matter of the source rock. Molecular dynamic experiments indicate that solubilized solvent and surfactant are delivered as part of a microemulsion droplet and are expected to aid the mobilization of oil present within kerogen.

The structural flexibility of DNA plays a key role in many biological processes of DNA, such as protein-DNA interactions, DNA packaging in viruses and nucleosome positioning on genomic DNA. Some experimental techniques have been employed to investigate the structural flexibility of DNA with the combination of elastic models, but these experiments could only provide the macroscopic properties of DNA, and thus, it is still difficult to understand the corresponding microscopic mechanisms. Recently, all-atom molecular dynamics (MD) simulation has emerged as a useful tool to investigate not only the macroscopic properties of DNA, but also the microscopic description of the flexibility of DNA at an atomic level. The most important issue in all-atom MD simulations of DNA is to choose an appropriate force field for simulating DNA. Very recently, a new force field for DNA has been developed based on the last generation force field of Amber bsc0, which was named Amber bsc1. In this work, all-atom MD simulations are employed to study the flexibility of a 30-bp DNA with the force fields of Amber bsc1 and Amber bsc0 in a comparative way. Our aim of the research is to examine the improvement of the new development of force field (Amber bsc1) in the macroscopic and microscopic properties of DNA, in comparison with the corresponding experimental measurements. All the MD simulations are performed with Gromacs 4.6 and lasted with a simulation time of 600 ns. The MD trajectories are analyzed with Curves+ for the last 500 ns, since the system reaches equilibrium approximately after ~100 ns. Our results show that the new force field (Amber bsc1) can lead to the improvements in the macroscopic parameters of DNA flexibility, i.e., stretch modulus S and twist-stretch coupling D become closer to experimental measurements, while bending persistence lengths lp and torsional persistence lengths C from the two force fields (bsc1 and bsc0) are both in good agreement with experimental data. Our microscopic analyses show that the microscopic structure parameters of DNA from the MD simulation with the Amber bsc1 force field are closer to the experimental values than those with the Amber bsc0 force field, except for slide, and the obvious improvements are observed in some microscopic parameters such as twist and inclination. Our further analyses show that the improvements in macroscopic flexibility from the Amber bsc1 force field are tightly related to the microscopic parameters and their fluctuations. This study would be helpful in understanding the performances of Amber bsc1 and bsc0 force fields in the description of DNA flexibility at both macroscopic and microscopic level.

Four commonly used molecular mechanics force fields, CHARMM22, OPLS, CVFF, and GROMOS87, are compared for their ability to reproduce experimental free energies of hydration (ΔGhydr) from molecular dynamics (MD) simulations for a set of small nonpolar and polar organic molecules: propane, cyclopropane, dimethylether, and acetone. ΔGhydr values were calculated by multiconfiguration thermodynamic integration for each of the different force fields with three different sets of partial atomic charges: full charges from an electrostatic potential fit (ESP), and ESP charges scaled by 0.8 and 0.6. All force fields, except for GROMOS87, give reasonable results for ΔGhydr · if partial atomic charges of appropriate magnitude are assigned. For GROMOS87, the agreement with experiment for hydrocarbons (propane and ethane) was improved considerably by modifying the repulsive part of the carbon-water oxygen Lennard-Jones potential. The small molecules studied are related to the chemical moieties constituting camphor (C10Hl6O). By invoking force-field transferability, we calculated the ΔGhydr for camphor. With the modified GROMOS force field, a ΔGhydr within 4 kJ/mol of the experimental value of −14.8 kJ/mol was obtained. Camphor is one of the largest molecules for which an absolute hydration free energy has been calculated by molecular simulation. The accuracy and reliability of the thermodynamic integration calculations were analyzed in detail and we found that, for ΔGhydr calculations for the set of small molecules in aqueous solution, molecular dynamics simulations of 0.8–1.0 ns in length give an upper statistical error bound of 1.5 kJ/mol, whereas shorter simulations of 0.25 nm in length given an upper statistical error bound of 3.5 kJ/mol. © 1997 by John Wiley & Sons, Inc.

Introduction: Milk kefir, a fermenting milk made with kefir grains, has shown potential in promoting apoptosis, regulating the cell cycle, and reducing tumor growth in breast cancer cells. This study aimed to investigate the stability and potential of milk kefir metabolites as inhibitors of breast cancer growth by interacting with the estrogen receptor alpha (ER-α), a key protein involved in breast cancer cell proliferation. We used computational methods, specifically molecular docking simulations with AutoDock and molecular dynamics (MD) simulations with Gromacs, to analyze how these metabolites bind to ER-α. Methods: A combination of molecular docking and MD simulations was used to explore how metabolites derived from milk kefir interact with ER-α, a crucial target in breast cancer therapy. The methodology included multiple stages: preparation of target proteins, preparation and screening of the metabolites, geometry optimization, molecular docking, and MD simulations. Results: The molecular docking simulations of 43 metabolites revealed three promising candidates: 2-Methyl (S35), benzeneethanol (S42), and 2,6-dimethyl-4-heptanone (S54), with binding affinities (ΔG) of -5.08, -5.06, and -4.90 kcal/mol, respectively. MD simulations further showed that the selected metabolites stabilized the ER-α-metabolite complex, with the 2,6-dimethyl-4-heptanone (S54) metabolite demonstrating the most negative total MM-GBSA energy value (ΔG = -22.98 kcal/mol), indicating a strong and stable binding interaction. Conclusion: 2,6-Dimethyl-4-heptanone, a metabolite from milk kefir, showed promising potential as a candidate for further development as a breast cancer treatment, offering a novel alternative to conventional therapies.

Molecular dynamics (MD) simulations based on atomic models play an important role in the drug-discovery process to screen molecules, estimate binding free energies, and optimize lead compounds in chemical space. Accurate computations of thermodynamic and kinetic properties using MD simulations are highly dependent on the accuracy of the underlying atomic force field. In this context, going beyond the nonpolarizable fixed-charge model by accounting explicitly for induced polarization is highly desirable. The CHARMM polarizable force field based on classical Drude oscillators, in which an auxiliary charged particle is attached via a harmonic spring to its parent nucleus, offers both a computationally convenient and rigorous framework to model explicitly induced electronic polarization in MD simulations. For any molecule of interest, electrostatic partial charges, atomic polarizabilities, and Thole shielding factors, as well as bonded parameters can either be determined from ab initio calculations or ascribed from the knowledge-based library of the CHARMM Generalized force field. While this approach is fairly reliable in general, it is well understood that the overall accuracy of the models with respect to thermodynamic properties such as bulk density, enthalpies, and solvation free energies is particularly sensitive to the nonbonded Lennard-Jones (LJ) parameters. In the present study, we systematically refined the set of LJ parameters for the atom types available in the Drude force field to best match the experimental thermodynamic properties for 416 small druglike organic molecules. To further test the transferability of the optimized parameters, the hydration free energy of 372 molecules was computed. The calculations resulted in a small average error of 0.46 kcal/mol and a Pearson R of 0.9, representing a significant improvement over the additive GAFF force field in our previous study, where an average error of ∼2 kcal/mol was obtained. Such an improvement is consistent with the ability of the polarizable Drude model to more accurately model interactions in different environments. The effort provides a roadmap for the global optimization of force field parameters using experimental data. It is hoped that the present effort will further the application of the Drude polarizable force field in molecular simulations including drug design and discovery.

Hydrophobic charge-induction chromatography (HCIC) with 4-mercaptoethyl-pyridine (MEP) as the ligand is a novel technology for antibody purification. In the present work, the molecular simulation methods were used to investigate the interactions between MEP ligand and Fc fragment of IgG (Fc-A). Six ligands with different structures of spacer arm were studied with molecular docking and dynamics simulation at neutral and acidic pH. The binding modes and the interaction energies were analyzed. The results indicated that all ligands tested could bind into the selected pocket on the C(H2) domain of Fc-A at neutral pH. The pyridine ring on the top of MEP ligands acts as a major role to provide the hydrophobic association and hydrogen bond for the ligand-IgG binding; meanwhile, the sulfone group on the spacer arm might form the additional hydrogen bond and enhance the binding of ligand onto the surface of IgG. The replacements of thioether sulfur atom on the spacer arm with either nitrogen or oxygen atom seem to have little influence on the binding. The influences of pH on the ligand-IgG interactions were also studied with the molecular dynamics simulation. It was found that MEP ligands would departed from the surface of Fc-A at low pH due to the electrostatic repulsion. The ligands with a sulfone group on the spacer arm would weaken the electrostatic repulsion and need more acidic conditions for the departing of ligand. The molecular simulation results were in agreement with some experimental observations, which would be useful to elucidate the molecular mechanism of HCIC and design a novel ligand to improve the efficiency of antibody separation.

This work provides an immersive visualization of oil & gas relevant systems and enhanced oil recovery (EOR) processes at the nanoscale by coupling molecular dynamics (MD) simulations with gamming virtual reality (VR) technologies. The main objective is to understand oil/brine/rock interfaces at molecular level and identify the underlying EOR mechanisms at the atomic scale. Within this immersive experience, the user can directly interact and enhance its perception of atomic environment for EOR applications. The experiences cover nano-EOR, nano-IOR and low-salt processes at nanoscale based on MD calculations of nanoparticles at oil-brine interfaces, oil-brine at silica nanopores and calcite-brine-oil interfaces, respectively. The MD simulations are performed with the Lammps package. The visualizations were done with an HTC Vive and Oculus Rift virtual reality headsets with the Nomad VR and Unitymol software. For the Nomad VR, the trajectories are previously saved from a Lammps molecular dynamics simulation, whereas, for the Unitymol, the simulation with Lammps is performed on-the-fly through the iMD (interactive MD) plugin. The user can visualize and navigate through the trajectories using the Nomad VR. Furthermore, the Unitymol also allows the user apply forces on selected molecules on real time during the VR experience. As means of comparison, the visualization was also performed with cell-phone based VR headsets with the Nomad VR application. Our demonstrations show that VR combined with molecular simulations can be an interesting and attractive way to improve the perception of the nanoscale for the general public. Additionally, it is an emergent tool to characterize, improve the understanding and provide molecular insights about nanosystems and the EOR methods, and also to be integrated with on-going digitalization processes within the oil & gas industry.

Molecular dynamics (MD) simulations are extensively used in the study of the structures and functions of proteins. Ab initio protein structure prediction is one of the most important subjects in computational biology, and many trials have been performed using MD simulation so far. Since the results of MD simulations largely depend on the force field, reliable force field parameters are indispensable for the success of MD simulation. In this work, we have modified atom charges in a standard force field on the basis of water-phase quantum chemical calculations. The modified force field turned out appropriate for ab initio protein structure prediction by the MD simulation with the generalized Born method. Detailed analysis was performed in terms of the conformational stability of amino acid residues, the stability of secondary structure of proteins, and the accuracy for prediction of protein tertiary structure, comparing the modified force field with a standard one. The energy balance between alpha-helix and beta-sheet structures was significantly improved by the modification of charge parameters.

In recent years, short peptide self-assembled materials, prepared under the control of the thermolysin catalyst, have been investigated extensively and shown to acquire various morphologies and fun...

New force fields for molecular dynamics (MD) simulation of aqueous zwitterionic amino acid simulations were developed. These were especially designed to calculate activity coefficient of water in amino acid solutions with high accuracy. For example, aqueous solutions of the following amino acids were considered: glycine, alanine, α-aminobutyric acid, α-aminovalerianic acid, valine and leucine. The force fields were obtained by quantum chemical calculations using B3LYP/6-31G and MP2/6-311(d,p) model theories in combination with the Merz–Kollmann–Singh scheme. To further increase the accuracy of the force field, a polarised continuum was considered in all quantum chemical calculations. Water activity coefficients obtained from MD using different all-purpose literature force fields, namely, OPLS, AMBER ff03 and GROMOS 53A6 as well as experimental data are compared with the results utilising the new force field. The new force field is shown to give better results compared with experimental data than existing force fields.