A FFLUX Water Model: Flexible, Polarizable and with a Multipolar Description of Electrostatics

FFLUX molecular simulations driven by atomic Gaussian process regression models

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.

Density functional theory study of crown ether-magnesium complexes: from a solvated ion to an ion trap

Metal ion detection rests on host-guest recognition. We propose a theoretical protocol for designing an optimal trap for a desired metal cation. A host for magnesium ions was sought for among derivatives of crown ethers 12-crown-4, 15-crown-5, and 18-crown-6. Mg-crown complexes and their hydrated counterparts with water molecules bound to the cation were optimized using density functional theory. Based on specific geometric criteria, Interacting quantum atoms analysis and density functional theory-based molecular dynamics of Mg-crown complexes immersed in water, crown ethers for optimal accommodation of Mg2+ in aqueous solution were identified. Selectivity of the chosen crowns towards Na+, K+, and Ca2+ ions is addressed.

Interacting Quantum Atoms and Multipolar Electrostatic Study of XH···π Interactions

The interaction energies of nine XH···π (X = C, N, and O) benzene-containing van der Waals complexes were analyzed, at the atomic and fragment levels, using QTAIM multipolar electrostatics and the energy partitioning method interacting quantum atoms/fragment (IQA/IQF). These descriptors were paired with the relative energy gradient method, which solidifies the connection between quantum mechanical properties and chemical interpretation. This combination provides a precise understanding, both qualitative and quantitative, of the nature of these interactions, which are ubiquitous in biochemical systems. The formation of the OH···π and NH···π systems is electrostatically driven, with the Qzz component of the quadrupole moment of the benzene carbons interacting with the charges of X and H in XH. There is the unexpectedly intramonomeric role of X-H (X = O, N) where its electrostatic energy helps the formation of the complex and its covalent energy thwarts it. However, the CH···π interaction is governed by exchange-correlation energies, thereby establishing a covalent character, as opposed to the literature's designation as a noncovalent interaction. Moreover, dispersion energy is relevant, statically and in absolute terms, but less relevant compared to other energy components in terms of the formation of the complex. Multipolar electrostatics are similar across all systems.

Construction of a Gaussian Process Regression Model of Formamide for Use in Molecular Simulations

FFLUX, a novel force field based on quantum chemical topology, can perform molecular dynamics simulations with flexible multipole moments that change with geometry. This is enabled by Gaussian process regression machine learning models, which accurately predict atomic energies and multipole moments up to the hexadecapole. We have constructed a model of the formamide monomer at the B3LYP/aug-cc-pVTZ level of theory capable of sub-kJ mol^-1 accuracy, with the maximum prediction error for the molecule being 0.8 kJ mol^-1. This model was used in FFLUX simulations along with Lennard-Jones parameters to successfully optimize the geometry of formamide dimers with errors smaller than 0.1 Å compared to those obtained with D3-corrected B3LYP/aug-cc-pVTZ. Comparisons were also made to a force field constructed with static multipole moments and Lennard-Jones parameters. FFLUX recovers the expected energy ranking of dimers compared to the literature, and changes in C═O and C-N bond lengths associated with hydrogen bonding were found to be consistent with density functional theory.

An Imbalance in the Force: The Need for Standardized Benchmarks for Molecular Simulation

Force fields (FFs) for molecular simulation have been under development for more than half a century. As with any predictive model, rigorous testing and comparisons of models critically depends on the availability of standardized data sets and benchmarks. While such benchmarks are rather common in the fields of quantum chemistry, this is not the case for empirical FFs. That is, few benchmarks are reused to evaluate FFs, and development teams rather use their own training and test sets. Here we present an overview of currently available tests and benchmarks for computational chemistry, focusing on organic compounds, including halogens and common ions, as FFs for these are the most common ones. We argue that many of the benchmark data sets from quantum chemistry can in fact be reused for evaluating FFs, but new gas phase data is still needed for compounds containing phosphorus and sulfur in different valence states. In addition, more nonequilibrium interaction energies and forces, as well as molecular properties such as electrostatic potentials around compounds, would be beneficial. For the condensed phases there is a large body of experimental data available, and tools to utilize these data in an automated fashion are under development. If FF developers, as well as researchers in artificial intelligence, would adopt a number of these data sets, it would become easier to compare the relative strengths and weaknesses of different models and to, eventually, restore the balance in the force.

Binding Energy Partition of Promising IRAK-4 Inhibitor (Zimlovisertib) for the Treatment of COVID-19 Pneumonia

The technique of Fragment-Based Drug Design (FBDD) considers the interactions of different moieties of molecules with biological targets for the rational construction of potential drugs. One basic assumption of FBDD is that the different functional groups of a ligand interact with a biological target in an approximately additive, that is, independent manner. We investigated the interactions of different fragments of ligands and Interleukin-1 Receptor-Associated Kinase 4 (IRAK-4) throughout the FBDD design of Zimlovisertib, a promising anti-inflammatory, currently in trials to be used for the treatment of COVID-19 pneumonia. We utilised state-of-the-art methods of wave function analyses mainly the Interacting Quantum Atoms (IQA) energy partition for this purpose. By means of IQA, we assessed the suitability of every change to the ligand in the five stages of FBDD which led to Zimlovisertib on a quantitative basis. We determined the energetics of the interaction of different functional groups in the ligands with the IRAK-4 protein target and thereby demonstrated the adequacy (or lack thereof) of the changes made across the design of this drug. This analysis permits to verify whether a given alteration of a prospective drug leads to the intended tuning of non-covalent interactions with its protein objective. Overall, we expect that the methods exploited in this paper will prove valuable in the understanding and control of chemical modifications across FBDD processes.

Non-covalent interactions from a Quantum Chemical Topology perspective

About half a century after its little-known beginnings, the quantum topological approach called QTAIM has grown into a widespread, but still not mainstream, methodology of interpretational quantum chemistry. Although often confused in textbooks with yet another population analysis, be it perhaps an elegant but somewhat esoteric one, QTAIM has been enriched with about a dozen other research areas sharing its main mathematical language, such as Interacting Quantum Atoms (IQA) or Electron Localisation Function (ELF), to form an overarching approach called Quantum Chemical Topology (QCT). Instead of reviewing the latter's role in understanding non-covalent interactions, we propose a number of ideas emerging from the full consequences of the space-filling nature of topological atoms, and discuss how they (will) impact on interatomic interactions, including non-covalent ones. The architecture of a force field called FFLUX, which is based on these ideas, is outlined. A new method called Relative Energy Gradient (REG) is put forward, which is able, by computation, to detect which fragments of a given molecular assembly govern the energetic behaviour of this whole assembly. This method can offer insight into the typical balance of competing atomic energies both in covalent and non-covalent case studies. A brief discussion on so-called bond critical points is given, highlighting concerns about their meaning, mainly in the arena of non-covalent interactions.

Application of Quantum Chemical Topology Force Field FFLUX to Condensed Matter Simulations: Liquid Water

We present here the first application of the quantum chemical topology force field FFLUX to condensed matter simulations. FFLUX offers many-body potential energy surfaces learnt exclusively from ab initio data using Gaussian process regression. FFLUX also includes high-rank, polarizable multipole moments (up to quadrupole moments in this work) that are learnt from the same ab initio calculations as the potential energy surfaces. Many-body effects (where a body is an atom) and polarization are captured by the machine learning models. The choice to use machine learning in this way allows the force field’s representation of reality to be improved (e.g., by including higher order many-body effects) with little detriment to the computational scaling of the code. In this manner, FFLUX is inherently future-proof. The “plug and play” nature of the machine learning models also ensures that FFLUX can be applied to any system of interest, not just liquid water. In this work we study liquid water across a range of temperatures and compare the predicted bulk properties to experiment as well as other state-of-the-art force fields AMOEBA(+CF), HIPPO, MB-Pol and SIBFA21. We find that FFLUX finds a place amongst these.

Using atomic charges to model molecular polarization

We review different models for introducing electric polarization in force fields, with special focus on methods where polarization is modelled at the atomic charge level. While electric polarization has been included in several force fields, the common approach has been to focus on atomic dipole polarizability. Several approaches allow modelling electric polarization by using charge-flow between charge sites instead, but this has been less exploited, despite that atomic charges and charge-flow is expected to be more important than atomic dipoles and dipole polarizability. A number of challenges are required to be solved for charge-flow models to be incorporated into polarizable force fields, for example how to parameterize the models and how to make them computational efficient.

DL_FFLUX: A Parallel, Quantum Chemical Topology Force Field

DL_FFLUX is a force field based on quantum chemical topology that can perform molecular dynamics for flexible molecules endowed with polarizable atomic multipole moments (up to hexadecapole). Using the machine learning method kriging (aka Gaussian process regression), DL_FFLUX has access to atomic properties (energy, charge, dipole moment, etc.) with quantum mechanical accuracy. Newly optimized and parallelized using domain decomposition Message Passing Interface (MPI), DL_FFLUX is now able to deliver this rigorous methodology at scale while still in reasonable time frames. DL_FFLUX is delivered as an add-on to the widely distributed molecular dynamics code DL_POLY 4.08. For the systems studied here (10³–10⁵ atoms), DL_FFLUX is shown to add minimal computational cost to the standard DL_POLY package. In fact, the optimization of the electrostatics in DL_FFLUX means that, when high-rank multipole moments are enabled, DL_FFLUX is up to 1.25× faster than standard DL_POLY. The parallel DL_FFLUX preserves the quality of the scaling of MPI implementation in standard DL_POLY. For the first time, it is feasible to use the full capability of DL_FFLUX to study systems that are large enough to be of real-world interest. For example, a fully flexible, high-rank polarized (up to and including quadrupole moments) 1 ns simulation of a system of 10 125 atoms (3375 water molecules) takes 30 h (wall time) on 18 cores.

Gaussian process models of potential energy surfaces with boundary optimization

A strategy is outlined to reduce the number of training points required to model intermolecular potentials using Gaussian processes, without reducing accuracy. An asymptotic function is used at a long range, and the crossover distance between this model and the Gaussian process is learnt from the training data. The results are presented for different implementations of this procedure, known as boundary optimization, across the following dimer systems: CO-Ne, HF-Ne, HF-Na⁺, CO₂-Ne, and (CO₂)₂. The technique reduces the number of training points, at fixed accuracy, by up to ∼49%, compared to our previous work based on a sequential learning technique. The approach is readily transferable to other statistical methods of prediction or modeling problems.

Systematic Comparison of the Structural and Dynamic Properties of Commonly Used Water Models for Molecular Dynamics Simulations

Water is a unique solvent that is ubiquitous in biology and present in a variety of solutions, mixtures, and materials settings. It therefore forms the basis for all molecular dynamics simulations of biological phenomena, as well as for many chemical, industrial, and materials investigations. Over the years, many water models have been developed, and it remains a challenge to find a single water model that accurately reproduces all experimental properties of water simultaneously. Here, we report a comprehensive comparison of structural and dynamic properties of 30 commonly used 3-point, 4-point, 5-point, and polarizable water models simulated using consistent settings and analysis methods. For the properties of density, coordination number, surface tension, dielectric constant, self-diffusion coefficient, and solvation free energy of methane, models published within the past two decades consistently show better agreement with experimental values compared to models published earlier, albeit with some notable exceptions. However, no single model reproduced all experimental values exactly, highlighting the need to carefully choose a water model for a particular study, depending on the phenomena of interest. Finally, machine learning algorithms quantified the relationship between the water model force field parameters and the resulting bulk properties, providing insight into the parameter-property relationship and illustrating the challenges of developing a water model that can accurately reproduce all properties of water simultaneously.

An electrostatic energy-based charge model for molecular dynamics simulation

The interactions of the polar chemical bonds such as C=O and N-H with an external electric field were investigated, and a linear relationship between the QM/MM interaction energies and the electric field along the chemical bond is established in the range of weak to intermediate electrical fields. The linear relationship indicates that the electrostatic interactions of a polar group with its surroundings can be described by a simple model of a dipole with constant moment under the action of an electric field. This relationship is employed to develop a general approach to generating an electrostatic energy-based charge (EEC) model for molecules containing single or multiple polar chemical bonds. Benchmark test studies of this model were carried out for (CH₃)₂-CO and N-methyl acetamide in explicit water, and the result shows that the EEC model gives more accurate electrostatic energies than those given by the widely used charge model based on fitting to the electrostatic potential (ESP) in direct comparison to the energies computed by the QM/MM method. The MD simulations of the electric field at the active site of ketosteroid isomerase based on EEC demonstrated that EEC gave a better representation of the electrostatic interaction in the hydrogen-bonding environment than the Amber14SB force field by comparison with experiment. The current study suggests that EEC should be better suited for molecular dynamics study of molecular systems with polar chemical bonds such as biomolecules than the widely used ESP or RESP (restrained ESP) charge models.

On the many-body nature of intramolecular forces in FFLUX and its implications

FFLUX is a biomolecular force field under construction, based on Quantum Chemical Topology (QCT) and machine learning (kriging), with a minimalistic and physically motivated design. A detailed analysis of the forces within the kriging models as treated in FFLUX is presented, taking as a test example a liquid water model. The energies of topological atoms are modeled as 3N_atoms -6 dimensional potential energy surfaces, using atomic local frames to represent the internal degrees of freedom. As a result, the forces within the kriging models in FFLUX are inherently N-body in nature where N refers to N_atoms . This provides a fuller picture that is closer to a true quantum mechanical representation of interactions between atoms. The presented computational example quantitatively showcases the non-negligible (as much as 9%) three-body nature of bonded forces and angular forces in a water molecule. We discuss the practical impact on the pressure calculation with N-body forces and periodic boundary conditions (PBC) in molecular dynamics, as opposed to classical force fields with two-body forces. The equivalence between the PBC-related correction terms in the general virial equation is shown mathematically.

Water in Nanopores and Biological Channels: A Molecular Simulation Perspective

This Review explores the dynamic behavior of water within nanopores and biological channels in lipid bilayer membranes. We focus on molecular simulation studies, alongside selected structural and other experimental investigations. Structures of biological nanopores and channels are reviewed, emphasizing those high-resolution crystal structures, which reveal water molecules within the transmembrane pores, which can be used to aid the interpretation of simulation studies. Different levels of molecular simulations of water within nanopores are described, with a focus on molecular dynamics (MD). In particular, models of water for MD simulations are discussed in detail to provide an evaluation of their use in simulations of water in nanopores. Simulation studies of the behavior of water in idealized models of nanopores have revealed aspects of the organization and dynamics of nanoconfined water, including wetting/dewetting in narrow hydrophobic nanopores. A survey of simulation studies in a range of nonbiological nanopores is presented, including carbon nanotubes, synthetic nanopores, model peptide nanopores, track-etched nanopores in polymer membranes, and hydroxylated and functionalized nanoporous silica. These reveal a complex relationship between pore size/geometry, the nature of the pore lining, and rates of water transport. Wider nanopores with hydrophobic linings favor water flow whereas narrower hydrophobic pores may show dewetting. Simulation studies over the past decade of the behavior of water in a range of biological nanopores are described, including porins and β-barrel protein nanopores, aquaporins and related polar solute pores, and a number of different classes of ion channels. Water is shown to play a key role in proton transport in biological channels and in hydrophobic gating of ion channels. An overall picture emerges, whereby the behavior of water in a nanopore may be predicted as a function of its hydrophobicity and radius. This informs our understanding of the functions of diverse channel structures and will aid the design of novel nanopores. Thus, our current level of understanding allows for the design of a nanopore which promotes wetting over dewetting or vice versa. However, to design a novel nanopore, which enables fast, selective, and gated flow of water de novo would remain challenging, suggesting a need for further detailed simulations alongside experimental evaluation of more complex nanopore systems.

A FFLUX Water Model: Flexible, Polarizable and with a Multipolar Description of Electrostatics

Key to progress in molecular simulation is the development of advanced models that go beyond the limitations of traditional force fields that employ a fixed, point charge-based description of electrostatics. Taking water as an example system, the FFLUX framework is shown capable of producing models that are flexible, polarizable and have a multipolar description of the electrostatics. The kriging machine-learning methods used in FFLUX are able to reproduce the intramolecular potential energy surface and multipole moments of a single water molecule with chemical accuracy using as few as 50 training configurations. Molecular dynamics simulations of water clusters (25-216 molecules) using the new FFLUX model reveal that incorporating charge-quadrupole, dipole-dipole, and quadrupole-charge interactions into the description of the electrostatics results in significant changes to the intermolecular structuring of the water molecules. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.

Implementation of Geometry-Dependent Charge Flux into the Polarizable AMOEBA+ Potential

Molecular dynamics (MD) simulations employing classical force fields (FFs) have been widely used to model molecular systems. The important ingredient of the current FFs, atomic charge, remains fixed during MD simulations despite the atomic environment or local geometry changes. This approximation hinders the transferability of the potential being used in multiple phases. Here we implement a geometry-dependent charge flux (GDCF) model into the multipole-based AMOEBA+ polarizable potential. The CF in the current work explicitly depends on the local geometry (bond and angle) of the molecule. To our knowledge, this is the first study that derives energy and force expressions due to GDCF in a multipole-based polarizable FF framework. Due to the inclusion of GDCF, the AMOEBA+ water model is noticeably improved in terms of describing the monomer properties, cluster binding/interaction energy, and a variety of liquid properties, including the infrared spectra that previous flexible water models were not able to capture.