Force Field Simulations Research Articles

Accurate Force Fields for Atomistic Simulations of Oxides, Hydroxides, and Organic Hybrid Materials up to the Micrometer Scale.

The simulation of metals, oxides, and hydroxides can accelerate the design of therapeutics, alloys, catalysts, cement-based materials, ceramics, bioinspired composites, and glasses. Here we introduce the INTERFACE force field (IFF) and surface models for α-Al2O3, α-Cr2O3, α-Fe2O3, NiO, CaO, MgO, β-Ca(OH)2, β-Mg(OH)2, and β-Ni(OH)2. The force field parameters are nonbonded, including atomic charges for Coulomb interactions, Lennard-Jones (LJ) potentials for van der Waals interactions with 12-6 and 9-6 options, and harmonic bond stretching for hydroxide ions. The models outperform DFT calculations and earlier atomistic models (Pedone, ReaxFF, UFF, CLAYFF) up to 2 orders of magnitude in reliability, compatibility, and interpretability due to a quantitative representation of chemical bonding consistent with other compounds across the periodic table and curated experimental data for validation. The IFF models exhibit average deviations of 0.2% in lattice parameters, <10% in surface energies (to the extent known), and 6% in bulk moduli relative to experiments. The parameters and models can be used with existing parameters for solvents, inorganic compounds, organic compounds, biomolecules, and polymers in IFF, CHARMM, CVFF, AMBER, OPLS-AA, PCFF, and COMPASS, to simulate bulk oxides, hydroxides, electrolyte interfaces, and multiphase, biological, and organic hybrid materials at length scales from atoms to micrometers. The nonbonded character of the models also enables the analysis of mixed oxides, glasses, and certain chemical reactions, and well-performing nonbonded models for silica phases, SiO2, are introduced. Automated model building is available in the CHARMM-GUI Nanomaterial Modeler. We illustrate applications of the models to predict the structure of mixed oxides, and energy barriers of ion migration, as well as binding energies of water and organic molecules in outstanding agreement with experimental data and calculations at the CCSD(T) level. Examples of model building for hydrated, pH-sensitive oxide surfaces to simulate solid-electrolyte interfaces are discussed.

Cathode Atomic Structures and Their Electrolyte Interfaces in Lithium‐sulfur Batteries

AbstractLithium‐sulfur batteries (LSBs) are a promising alternative to conventional battery technologies due to their high theoretical energy density. However, the insulating nature of sulfur and the polysulfide shuttle effect present practical challenges to achieving this. This study employs molecular dynamics simulations with a ReaxFF reactive force field to investigate the structure of LSB cathode‐electrolyte interfaces, with cathode materials made from different sized graphene flakes and elemental sulfur (S8). Smaller graphene flakes were seen to bond more with the sulfur, while structures with larger graphene flakes featured fewer C−S bonds and more long‐chain sulfur species. Observing the interface between the cathode and an organic electrolyte revealed that cathode structures with larger graphene flakes could orientate to shield the long‐chain sulfur and lithium polysulfide species from interacting with and dissolving into the electrolyte. Additional simulations were performed to investigate cathode materials with regions of crystalline α‐S8. However, it was found that the ReaxFF force field used for LSB systems cannot effectively model crystalline sulfur under constant pressure conditions. Further development of reactive force fields for simulation of graphene flake and sulfur electrodes in the presence of electrolyte are therefore required to explore the full range of likely structures using these methods.

Novel Computational Chemistry Infrastructure for Simulating Astatide in Water: From Basis Sets to Force Fields Using Particle Swarm Optimization.

Using the example of astatine, the heaviest naturally occurring halogen whose isotope At-211 has promising medical applications, we propose a new infrastructure for large-scale computational models of heavy elements with strong relativistic effects. In particular, we focus on developing an accurate force field for At- in water based on reliable relativistic density functional theory (DFT) calculations. To ensure the reliability of such calculations, we design novel basis sets for relativistic DFT, via the particle swarm optimization algorithm to optimize the coefficients of the new basis sets and the polarization-consistent basis set idea's extension to heavy elements to eliminate the basis set error from DFT calculations. The resulting basis sets enable the well-grounded evaluation of relativistic DFT against "gold-standard" CCSD(T) results. Accounting for strong relativistic effects, including spin-orbit interaction, via our redesigned infrastructure, we elucidate a noticeable dissimilarity between At- and I- in halide-water force field parameters, radial distribution functions, diffusion coefficients, and hydration energies. This work establishes the framework for the systematic development of polarization-consistent basis sets for relativistic DFT and accurate force fields for molecular dynamics simulations to be used in large-scale models of complex molecular systems with elements from the bottom of the periodic table, including actinides and even superheavy elements.

Derivation of a Force Field for Computer Simulations of Multi-Walled Nanotubes Using Genetic Algorithm. I. Tungsten Disulfide

A technique for constructing force fields based on the use of genetic algorithms is proposed, which is aimed at parameterization of potentials intended for computer simulation of polyatomic nanosystems. To illustrate the proposed approach, a force field has been developed for modeling layered modifications of WS2, including multi-walled nanotubes, the dimensions of which are beyond the capabilities of ab initio methods. When determining the potential parameters, layered polytypes of bulk crystals, monolayers, bilayers, and nanotubes of small diameters were used as calibration systems. The parameterization found was successfully tested on double-walled nanotubes, the structure of which was determined using density functional calculations. The obtained force field was used for the first time to model the structure and stability of achiral multi-walled nanotubes based on WS2. The interwall distances obtained from the simulation are in good agreement with the results of recent measurements of these parameters for existing nanotubes.

Derivation of a Force Field for Computer Simulations of Multi-Walled Nanotubes Using Genetic Algorithm. I. Tungsten Disulfide

Derivation of a Force Field for Computer Simulations of Multi-Walled Nanotubes Using Genetic Algorithm. I. Tungsten Disulfide

Theoretical models of staurosporine and analogs uncover detailed structural information in biological solution

Staurosporine and its analogs (STA-analogs) are indolocarbazoles (ICZs) compounds able to inhibit kinase proteins in a non-specific way, while present antimicrobial and cytostatic properties. The knowledge of molecular features associated to the complexation, including the ligand shape in solution and thermodynamics of complexation, is substantial to the development of new bioactive ICZs with improved therapeutic properties. In this context, the empirical approach of GROMOS force field is able to accurately reproduce condensed phase physicochemical properties of molecular systems after parameterization. Hence, through parameterization under GROMOS force field and molecular simulations, we assessed STA-analogs dynamics in aqueous solution, as well as its interaction with water to probe conformational and structural features involved in complexation to therapeutic targets. The coexistence of multiple conformers observed in simulations, and confirmed by metadynamics calculations, expanding the conformational space knowledge of these ligands with potential implications in understanding the ligand conformational selection during complexation. Also, changes in availability to H-bonding concerning the different substituents and water can reflect on effects at complexation free energy due to variation at the desolvation energetic costs. Based on these results, we expect the obtained structural data provide systemic framework for rational chemical modification of STA-analogs.

Phase transition behaviors of Al nanoparticles with low oxidation degree: A molecular dynamics study

In this study, the ReaxFF reactive force field is applied to simulate the melting and annealing behaviors of aluminum (Al) nanoparticles (ANP) by molecular dynamics (MD) simulations. Potential energy, specific heat and FCC lattice number are used to study the thermal behavior of the entire phase transition process. First, the melting and solidification process of 3[Formula: see text]nm cubic Al and Al2O3 models under constant-pressure/temperature ensemble (NPT) system was simulated without boundary, and the steady-state simulation was discussed at specific temperature points. The non-boundary simulation proved the effectiveness of the Reaxff reactive force field for ANP phase transition simulations. It is confirmed that the physical properties of Al and Al2O3 could be effectively reflected in the MD simulations. ANP with low oxidation degree cannot be completely ignored in engineering applications. The existence of oxide layer hinders the annealing behavior. Therefore, this study also carried out an additional simulation process of 6 nm ANP and its oxide particles in canonical ensemble (NVT). Those obtained results show that a small amount of oxide layer can greatly change the thermal properties of ANP. If the surface of ANP is isolated, then its melting point could be increased by at least 30[Formula: see text]K.

Nonaqueous Ion Pairing Exemplifies the Case for Including Electronic Polarization in Molecular Dynamics Simulations.

The inclusion of electronic polarization is of crucial importance in molecular simulations of systems containing charged moieties. When neglected, as often done in force field simulations, charge-charge interactions in solution may become severely overestimated, leading to unrealistically strong bindings of ions to biomolecules. The electronic continuum correction introduces electronic polarization in a mean-field way via scaling of charges by the reciprocal of the square root of the high-frequency dielectric constant of the solvent environment. Here, we use ab initio molecular dynamics simulations to quantify the effect of electronic polarization on pairs of like-charged ions in a model nonaqueous environment where electronic polarization is the only dielectric response. Our findings confirm the conceptual validity of this approach, underlining its applicability to complex aqueous biomolecular systems. Simultaneously, the results presented here justify the potential employment of weaker charge scaling factors in force field development.

On-surface synthesis of ballbot-type N-heterocyclic carbene polymers.

N-Heterocyclic carbenes (NHCs) are established ligands for metal complexes and surfaces. Here we go beyond monomeric NHCs and report on the synthesis of NHC polymers on gold surfaces, consisting of ballbot-type repeating units bound to single Au adatoms. We designed, synthesized and deposited precursors containing different halogens on gold surfaces under ultrahigh vacuum. Conformational, electronic and charge transport properties were assessed by combining low-temperature scanning tunneling microscopy, non-contact atomic force microscopy, X-ray photoelectron spectroscopy, first-principles calculations and reactive force field simulations. The confirmed ballbot-type nature of the NHCs explains the high surface mobility of the incommensurate NHC polymers, which is prerequisite for their desired spatial alignment. The delicate balance between mobility and polymerization rate allows essential parameters for controlling polymer directionality to be derived. These polymers open up new opportunities in the fields of nanoelectronics, surface functionalization and catalysis.

Mechanism investigations on co-pyrolysis of polyethylene and biomass using ReaxFF simulation and DFT computation

Co-pyrolysis of plastics and biomass is a promising solution for waste management and renewable energy utilization. There is still no consensus on the distributions of co-pyrolysis products and the reaction mechanisms between the polyethylene (PE) plastics and biomass components at the molecular level. In this work, the product distribution, co-pyrolysis mechanism, and electronic characteristic of PE and major biomass components such as cellulose (CE), hemicellulose (HC), and lignin (LG) were investigated at the atomistic level using reactive force field simulation and density functional theory computation. The co-pyrolysis of PE and biomass facilitates the release of non-carbon radicals, which enhances the reactivity of macromolecular fractions to high-quality oils. HC promotes the char degradation and oxygenated oil retention, whereas CE accelerates the decomposition of oxygen-free char to light hydrocarbons and inorganic gases. LG inhibits the conversion of oxygenated oil to CO/CO2 and retards the secondary cracking of oxygen-free oil. Alkyne is mostly generated from α-O-4 LG dimer, whereas olefin is mainly derived from γ-O-4 LG dimer. β-O-4 LG dimer accelerates the conversion of char to oil, while it decreases gas emissions. The electronic characteristics such as the Fukui function, frontier orbitals, and total electron density demonstrate the reaction sites of biomass monomers. For CE and HC, the hydroxyl groups on methylene are the nucleophilic sites, and the O atoms of hydroxyl on the pyran ring are prone to the electrophilic attack by PE. Except for O atoms in LG dimers, the unsubstituted C atoms on the benzene ring are inclined to the nucleophilic sites, whereas the C and para-C atoms located at hydroxyl and methoxy are susceptible to the electrophilic attack by PE. This work elucidates the differences in the co-pyrolysis reaction characteristics of various biomass structures, which provides scientific and theoretical guidance for the efficient resource exploitation of biomass and waste plastics.

Polarizable atomic multipole-based force field for cholesterol

Cholesterol is one of the essential component of lipid in membrane. We present a polarizable atomic multipole force field (FF) for the molecular dynamic simulation of cholesterol. The FF building process follows the computational framework as the atomic multipole optimized energetics for biomolecular applications (AMOEBA) model. In this framework, the electronics parameters, including atomic monopole moments, dipole moments, and quadrupole moments calculated from ab initio calculations in the gas phase, are applied to represent the charge distribution. Furthermore, the many-body polarization is modeled by following the same pattern of distributed atomic polarizabilities. Then, the bilayers composed of two typical phospholipid molecules, 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC), and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), in a range of different cholesterol concentrations are built and implemented by molecular dynamics (MD) simulations based on the proposed polarizable FF. The simulation results are statistically analyzed to validate the feasibility of the proposed FF. The structural properties of the bilayers are calculated to compare with the related experimental values. The MD values show the same trend of experimental values changes. Communicated by Ramaswamy H. Sarma

Data scheme and data format for transferable force fields for molecular simulation

A generalized data scheme for transferable classical force fields used in molecular simulations, i.e. molecular dynamics and Monte Carlo simulation, is presented. The data scheme is implemented in an SQL-based data format. The data scheme and data format is machine readable, re-usable, and interoperable. A transferable force field is a chemical construction plan specifying intermolecular and intramolecular interactions between different types of atoms or different chemical groups and can be used for building a model for a given component. The data scheme proposed in this work (named TUK-FFDat) formalizes digitally these chemical construction plans, i.e. transferable force fields. It can be applied to all-atom as well as united-atom transferable force fields. The general applicability of the data scheme is demonstrated for different types of force fields (TraPPE, OPLS-AA, and Potoff). Furthermore, conversion tools for translating the data scheme between .xls spread sheet format and the SQL-based data format are provided. The data format can readily be integrated in existing workflows, simulation engines, and force field databases as well as for linking such.

Improving machine learning force fields for molecular dynamics simulations with fine-grained force metrics.

Machine learning force fields (MLFFs) have gained popularity in recent years as they provide a cost-effective alternative to abinitio molecular dynamics (MD) simulations. Despite a small error on the test set, MLFFs inherently suffer from generalization and robustness issues during MD simulations. To alleviate these issues, we propose global force metrics and fine-grained metrics from element and conformation aspects to systematically measure MLFFs for every atom and every conformation of molecules. We selected three state-of-the-art MLFFs (ET, NequIP, and ViSNet) and comprehensively evaluated on aspirin, Ac-Ala3-NHMe, and Chignolin MD datasets with the number of atoms ranging from 21 to 166. Driven by the trained MLFFs on these molecules, we performed MD simulations from different initial conformations, analyzed the relationship between the force metrics and the stability of simulation trajectories, and investigated the reason for collapsed simulations. Finally, the performance of MLFFs and the stability of MD simulations can be further improved guided by the proposed force metrics for model training, specifically training MLFF models with these force metrics as loss functions, fine-tuning by reweighting samples in the original dataset, and continued training by recruiting additional unexplored data.

Toward a general neural network force field for protein simulations: Refining the intramolecular interaction in protein.

Molecular dynamics (MD) is an extremely powerful, highly effective, and widely used approach to understanding the nature of chemical processes in atomic details for proteins. The accuracy of results from MD simulations is highly dependent on force fields. Currently, molecular mechanical (MM) force fields are mainly utilized in MD simulations because of their low computational cost. Quantum mechanical (QM) calculation has high accuracy, but it is exceedingly time consuming for protein simulations. Machine learning (ML) provides the capability for generating accurate potential at the QM level without increasing much computational effort for specific systems that can be studied at the QM level. However, the construction of general machine learned force fields, needed for broad applications and large and complex systems, is still challenging. Here, general and transferable neural network (NN) force fields based on CHARMM force fields, named CHARMM-NN, are constructed for proteins by training NN models on 27 fragments partitioned from the residue-based systematic molecular fragmentation (rSMF) method. The NN for each fragment is based on atom types and uses new input features that are similar to MM inputs, including bonds, angles, dihedrals, and non-bonded terms, which enhance the compatibility of CHARMM-NN to MM MD and enable the implementation of CHARMM-NN force fields in different MD programs. While the main part of the energy of the protein is based on rSMF and NN, the nonbonded interactions between the fragments and with water are taken from the CHARMM force field through mechanical embedding. The validations of the method for dipeptides on geometric data, relative potential energies, and structural reorganization energies demonstrate that the CHARMM-NN local minima on the potential energy surface are very accurate approximations to QM, showing the success of CHARMM-NN for bonded interactions. However, the MD simulations on peptides and proteins indicate that more accurate methods to represent protein-water interactions in fragments and non-bonded interactions between fragments should be considered in the future improvement of CHARMM-NN, which can increase the accuracy of approximation beyond the current mechanical embedding QM/MM level.

Optimized Quantum Drude Oscillators for Atomic and Molecular Response Properties.

The quantum Drude oscillator (QDO) is an efficient yet accurate coarse-grained approach that has been widely used to model electronic and optical response properties of atoms and molecules as well as polarization and dispersion interactions between them. Three effective parameters (frequency, mass, and charge) fully characterize the QDO Hamiltonian and are adjusted to reproduce response properties. However, the soaring success of coupled QDOs for many-atom systems remains fundamentally unexplained, and the optimal mapping between atoms/molecules and oscillators has not been established. Here we present an optimized parametrization (OQDO) where the parameters are fixed by using only dipolar properties. For the periodic table of elements as well as small molecules, our model accurately reproduces atomic (spatial) polarization potentials and multipolar dispersion coefficients, elucidating the high promise of the presented model in the development of next-generation quantum-mechanical force fields for (bio)molecular simulations.

Insights into the Differential Composition of Stem-Loop Structures of Nanoviruses and Their Impacts.

Multipartite viruses package their genomic segments independently and mainly infect plants; few of them target animals. Nanoviridae is a family of multipartite single-stranded DNA (ssDNA) plant viruses that individually encapsidate ssDNAs of ~1 kb and transmit them through aphids without replication in aphid vectors, thereby causing important diseases in host plants, mainly leguminous crops. All of these components constitute an open reading frame to perform a specific role in nanovirus infection. All segments contain conserved inverted repeat sequences, potentially forming a stem-loop structure and a conserved nonanucleotide, TAGTATTAC, within a common region. This study investigated the variations in the stem-loop structure of nanovirus segments and their impact using molecular dynamics (MD) simulations and wet lab approaches. Although the accuracy of MD simulations is limited by force field approximations and simulation time scale, explicit solvent MD simulations were successfully used to analyze the important aspects of the stem-loop structure. This study involves the mutants' design, based on the variations in the stem-loop region and construction of infectious clones, followed by their inoculation and expression analysis, based on nanosecond dynamics of the stem-loop structure. The original stem-loop structures showed more conformational stability than mutant stem-loop structures. The mutant structures were expected to alter the neck region of the stem-loop by adding and switching nucleotides. Changes in conformational stability are suggested expression variations of the stem-loop structures found in host plants with nanovirus infection. However, our results can be a starting point for further structural and functional analysis of nanovirus infection. IMPORTANCE Nanoviruses comprise multiple segments, each with a single open reading frame to perform a specific function and an intergenic region with a conserved stem-loop region. The genome expression of a nanovirus has been an intriguing area but is still poorly understood. We attempted to investigate the variations in the stem-loop structure of nanovirus segments and their impact on viral expression. Our results show that the stem-loop composition is essential in controlling the virus segments' expression level.

Ab Initio Derived Classical Force Field for Molecular Dynamics Simulations of ZnO Surfaces in Biological Environment.

Zinc oxide nanostructures are used in an ever increasing line of applications in technology and biomedical fields. This requires a detailed understanding of the phenomena that occur at the surface particularly in aqueous environments and in contact with biomolecules. In this work, we used ab initio molecular dynamics (AIMD) simulations to determine structural details of ZnO surfaces in water and to develop a general and transferable classical force field for hydrated ZnO surfaces. AIMD simulations show that water molecules dissociate near unmodified ZnO surfaces, forming hydroxyl groups at about 65% of the surface Zn atoms and protonating 3-coordinated surface oxygen atoms, while the rest of the surface Zn atoms bind molecularly adsorbed waters. Several force field atom types for ZnO surface atoms were identified by analysis of the specific connectivities of atoms. The analysis of the electron density was then used to determine partial charges and Lennard-Jones parameters for the identified force field atom types. The obtained force field was validated by comparison with AIMD results and with available experimental data on adsorption and immersion enthalpies, as well as adsorption free energies of several amino acids in methanol. The developed force field can be used for modeling of ZnO in aqueous and other fluid environments and in interaction with biomolecules.

Free Energy Calculations and Conformational Analysis of Dibenzo-30-crown-10 with Sm2+, Eu2+, and Three Halide Salts in THF Using the AMOEBA Force Field.

Crown ether complexes have been tailored for use in industrial separations of lanthanides (Ln) as a part of rare earth mining and refining. Dibenzo-30-crown-10 (DB30C10) is one of the most efficient complexants for the separation of rare earth mixtures based on the cation size. To understand the origin of this complexation, molecular dynamics (MD) simulations of DB30C10 have been performed using different combinations of divalent Sm and Eu and three halide salts Cl-, Br-, and I- in tetrahydrofuran (THF) solvent. DB30C10 was parameterized here for the polarizable atomic multipole optimized energetics for biomolecular simulation (AMOEBA) force field, and the existing parameters of THF, Sm2+, and Eu2+ were employed from our previous efforts. The large conformational fluctuations present in the DB30C10 systems were found to be dependent both on the identity of the lanthanide and halide complexes. For Cl- and Br- systems, there were no observed conformational changes at 200 ns, while in I- systems, there were two conformational changes with Sm2+ and one with Eu2+ within that same timeframe. In SmI2-DB30C10, there were three stages of conformational changes. In the first stage, the molecule is unfolded, in the second stage, the molecule is partly folded, and finally, in the third stage, the molecule is completely folded. Lastly, the Gibbs binding free energies of DB30C10 with SmBr2 and EuBr2 have been computed, which resulted in nearly identical ΔGcomp values for each lanthanide with Sm2+ being slightly more favorable. Considering the folding mechanism of the SmI2 system with DB30C10, the Gibbs binding free energies of DB30C10 and dicyclohexano-18-crown-6 (DCH18C6) with SmI2 were calculated separately and compared to probe their complexation affinities, in which the former was found to be more favorable.

A Reactive Force Field for Molecular Dynamics Simulations of Glucose in Aqueous Solution.

To expand the capabilities of reactive force field (ReaxFF) in simulations of biological processes involving glucose, in this work, using Metropolis Monte Carlo algorithm, new ReaxFF parameters for glucose have been developed to better describe the properties of glucose in water during molecular dynamics (MD) simulations. With the newly trained ReaxFF, the mutarotation of glucose in water can be better described, as suggested by our metadynamics simulations. In addition, the newly trained ReaxFF can better describe the distributions of the three stable conformers along the key dihedral angle of α-anomer and β-anomer. With better descriptions of hydration around glucose, the Raman and Raman optical activity spectra can be more accurately calculated. In addition, the infrared spectra obtained from simulations with the new glucose ReaxFF are more accurate than those obtained with the original ReaxFF. We note that although our trained ReaxFF performs better than the original ReaxFF, it is not generally applicable to all carbohydrates, which require further parametrization. We also find that the absence of explicit water molecules in the training sets may lead to inaccurate descriptions of water-water interactions around the glucose, implicating that it is necessary to optimize the water ReaxFF parameters together with the target molecule. The improved ReaxFF makes it possible to explore interesting biological processes involving glucose more accurately and efficiently.

A brief history of visualizing membrane systems in molecular dynamics simulations.

Understanding lipid dynamics and function, from the level of single, isolated molecules to large assemblies, is more than ever an intensive area of research. The interactions of lipids with other molecules, particularly membrane proteins, are now extensively studied. With advances in the development of force fields for molecular dynamics simulations (MD) and increases in computational resources, the creation of realistic and complex membrane systems is now common. In this perspective, we will review four decades of the history of molecular dynamics simulations applied to membranes and lipids through the prism of molecular graphics.