Polarizable force fields

Evaluating the Functional Importance of Conformer-Dependent Atomic Partial Charge Assignment

Physics-based methods such as protein-ligand binding free energy calculations have been increasingly adopted in early-stage drug discovery to prioritize promising compounds for synthesis. However, the accuracy of these methods is highly dependent on the details of the calculation and choices made while preparing the ligands and protein ahead of running calculations. During ligand preparation, researchers typically assign partial atomic charges to each ligand atom using a specific ligand conformation for charge assignment, often the input conformer. While it is a well-known problem that partial charge assignment is dependent on conformation, little investigation has explored the downstream effects of varied partial charge assignment on free energy estimates. Preliminary benchmarks from the Open Free Energy Project show that generating partial charges from different input conformers leads to variation of up to± 5 . 3 $$ \pm 5.3 $$  kcal/mol in calculated relative binding free energies due to variation in partial charges alone. In this study, we more systematically explore this issue, investigating it in smaller systems using absolute hydration free energy calculations to reduce the degrees of freedom and sources of statistical error as compared to larger protein-ligand systems. We investigate how differences in partial charge generation (such as those caused by input conformer choice, partial charge engine, and hardware) may lead to differences in calculated absolute hydration free energy (AHFE) values. We demonstrate that supplying different input conformers to a partial charge engine can result in atomic partial charge discrepancies of up to 0.681 e, resulting in differences in calculated AHFE of6 . 9 ± 0 . 1 $$ 6.9\pm 0.1 $$  kcal/mol. We find that even relatively small variations in partial charge assignment can result in notable differences in calculated AHFE. Thus, care should be taken when assigning partial charges to ensure reproducibility and accuracy of any resulting free energy calculations. We expect that these effects will be magnified in pharmaceutically relevant binding free energy calculations with additional degrees of freedom, more highly directional interactions than in water, and potentially more statistical error.

Shadow Molecular Dynamics for a Charge-Potential Equilibration Model

We introduce a shadow molecular dynamics (MD) approach based on the Atom-Condensed Kohn-Sham second-order (ACKS2) charge-potential equilibration model. In contrast to regular flexible charge models, the ACKS2 model includes both flexible atomic charges and potential fluctuation parameters that allow for physically correct charge fragmentation and improved scaling of the polarizability. Our shadow MD scheme is based on an approximation of the ACKS2's flexible charge-potential energy function, in combination with extended Lagrangian Born-Oppenheimer MD. Utilizing this shadow charge-potential equilibration approach mitigates the costly overhead and stability problems associated with finding well-converged iterative solutions to the charges and potential fluctuations of the ACKS2 model in an MD simulation. Our work provides a robust and versatile framework for efficient, high-fidelity MD simulations of diverse physical phenomena and applications.

Role of Electron Spin, Chirality, and Charge Dynamics in Promoting the Persistence of Nascent Nucleic Acid-Peptide Complexes

Primitive nucleic acids and peptides likely collaborated in early biochemistry. What forces drove their interactions and how did these forces shape the properties of primitive complexes? We investigated how two model primordial polypeptides associate with DNA. When peptides were coupled to a ferromagnetic substrate, DNA binding depended on the substrate's magnetic moment orientation. Reversing the magnetic field nearly abolished binding despite complementary charges. Inverting the peptide chirality or just the cysteine residue reversed this effect. These results are attributed to the chiral-induced spin selectivity (CISS) effect, where molecular chirality and electron spin alter a protein's electric polarizability. The presence of CISS in simple protein-DNA complexes suggests that it played a significant role in ancient biomolecular interactions. A major consequence of CISS is enhancement of the kinetic stability of protein-nucleic acid complexes. These findings reveal how chirality and spin influence bioassociation, offering insights into primitive biochemical evolution and shaping contemporary protein functions.

Isotropic Periodic Sum for Polarizable Gaussian Multipole Model

The isotropic periodic sum (IPS) method provides an efficient approach for computing long-range interactions by approximating distant molecular structures through isotropic periodic images of a local region. Here, we present a novel integration of IPS with the polarizable Gaussian multipole (pGM) model, extending its applicability to systems with Gaussian-distributed charges and dipoles. By developing and implementing the IPS multipole tensor theorem within the Gaussian multipole framework, we derive analytical expressions for IPS potentials that efficiently handle both permanent and induced multipole interactions. Our comprehensive validation includes energy conservation tests in the NVE ensemble, potential energy distributions in the NVT ensemble, structural analysis through radial distribution functions, diffusion coefficients, induced dipole calculations across various molecular systems, and ionic charging free energies. The results demonstrate that the pGM-IPS approach successfully reproduces energetic, structural, and dynamic properties of molecular systems with accuracy comparable to the traditional particle mesh Ewald method. Our work establishes pGM-IPS as a promising method for simulations of polarizable molecular systems, achieving a balance between computational efficiency and accuracy.

Shadow Molecular Dynamics with a Machine Learned Flexible Charge Potential

We present an extended Lagrangian shadow molecular dynamics scheme with an interatomic Born-Oppenheimer potential determined by the relaxed atomic charges of a second-order charge equilibration model. To parametrize the charge equilibration model, we use machine learning with neural networks to determine the environment-dependent electronegativities and chemical hardness parameters for each atom, in addition to the charge-independent energy and force terms. The approximate shadow molecular dynamics potential in combination with the extended Lagrangian formulation improves the numerical stability and reduces the number of Coulomb potential calculations required to evaluate accurate conservative forces. We demonstrate efficient and accurate simulations with excellent long-term stability of the molecular dynamics trajectories. The significance of choosing fixed or environment-dependent electronegativities and chemical hardness parameters is evaluated. Finally, we compute the infrared spectrum of molecules via the dipole autocorrelation function and compare to experiments to highlight the accuracy of the shadow molecular dynamics scheme with a machine learned flexible charge potential.

Simulating Ionized States in Realistic Chemical Environments with Algebraic Diagrammatic Construction Theory and Polarizable Embedding

Theoretical simulations of electron detachment processes are vital for understanding chemical redox reactions, semiconductor and electrochemical properties, and high-energy radiation damage. However, accurate calculations of ionized electronic states are very challenging due to their open-shell nature, importance of electron correlation effects, and strong interactions with chemical environment. In this work, we present an efficient approach based on algebraic diagrammatic construction theory with polarizable embedding that allows to accurately simulate ionized electronic states in condensed-phase or biochemical environments (PE-IP-ADC). We showcase the capabilities of PE-IP-ADC by computing the vertical ionization energy (VIE) of thymine molecule solvated in bulk water. Our results show that the second- and third-order PE-IP-ADC methods combined with the basis of set of triple-ζ quality yield a solvent-induced shift in VIE of -0.92 and -0.93 eV, respectively, in an excellent agreement with experimental estimate of -0.9 eV. This work demonstrates the power of PE-IP-ADC approach for simulating charged electronic states in realistic chemical environments and motivates its further development.

Performance Tuning of Polarizable Gaussian Multipole Model in Molecular Dynamics Simulations

Molecular dynamics (MD) simulations are essential for understanding molecular phenomena at the atomic level, with their accuracy largely dependent on both the employed force field and sampling. Polarizable force fields, which incorporate atomic polarization effects, represent a significant advancement in simulation technology. The polarizable Gaussian multipole (pGM) model has been noted for its accurate reproduction of ab initio electrostatic interactions. In this study, we document our effort to enhance the computational efficiency and scalability of the pGM simulations within the AMBER framework using MPI (message passing interface). Performance evaluations reveal that our MPI-based pGM model significantly reduces runtime and scales effectively while maintaining computational accuracy. Additionally, we investigated the stability and reliability of the MPI implementation under the NVE simulation ensemble. Optimal Ewald and induction parameters for the pGM model are also explored, and its statistical properties are assessed under various simulation ensembles. Our findings demonstrate that the MPI-implementation maintains enhanced stability and robustness during extended simulation times. We further evaluated the model performance under both NVT (constant number, volume, and temperature) and NPT (constant number, pressure, and temperature) ensembles and assessed the effects of varying timesteps and convergence tolerance on induced dipole calculations. The lessons learned from these exercises are expected to help the users to make informed decisions on simulation setup. The improved performance under these ensembles enables the study of larger molecular systems, thereby expanding the applicability of the pGM model in detailed MD simulations.

Revised 4-Point Water Model for the Classical Drude Oscillator Polarizable Force Field: SWM4-HLJ

In this work the 4-point polarizable SWM4 Drude water model is reparametrized. Multiple models were developed using different strategies toward reproduction of specific target data. Results indicate that no individual model can reproduce all the selected target data in the context of the present form of the potential energy function. The changes considered in the new models include, 1) variations in the gas phase dipole moment, 2) variations in the molecular polarizability, 3) variations of the distance between the oxygen and the M site, 4) variation of the oxygen Lennard-Jones (LJ) parameters, 5) introduction of a LJ potential to the hydrogen atoms, and 6) variations of the H-O-H angle. Detailed analysis is presented for 3 new water models from which a final model, SWM4-HLJ, is selected as the future default model for the Drude polarizable force field. The model maintains the gas phase dipole moment as the experimental value while the remaining listed terms were adjusted including a larger H-O-H angle (108.12). Compared to its predecessor, SWM4-NDP, the self-diffusion coefficient, water dimer properties, and water cluster energies are greatly improved. The temperature dependence of the density of the new model also performs better. Overall, the new SWM4-HLJ water model is a general improvement and a good balance between microscopic and bulk properties is achieved.

Theoretically grounded approaches to account for polarization effects in fixed-charge force fields

Non-polarizable, or fixed-charge, force fields are the workhorses of most molecular simulation studies. They attempt to describe the potential energy surface (PES) of the system by including polarization effects in an implicit way. This has historically been done in a rather empirical and ad hoc manner. Recent theoretical treatments of polarization, however, offer promise for getting the most out of fixed-charge force fields by judicious choice of parameters (most significantly the net charge or dipole moment of the model) and application of post facto polarization corrections. This Perspective describes these polarization theories, namely the "halfway-charge" theory and the molecular dynamics in electronic continuum theory, and shows that they lead to qualitatively (and often, quantitatively) similar predictions. Moreover, they can be reconciled into a unified approach to construct a force field development workflow that can yield non-polarizable models with charge/dipole values that provide an optimal description of the PES. Several applications of this approach are reviewed, and avenues for future research are proposed.

Anomalous Water Penetration in Al3+ Dissolution

The physicochemical characterization of trivalent ions is limited due to a lack of accurate force fields. By leveraging the latest machine learning force field to model aqueous AlCl3, we discover that upon dissolution of Al3+, water molecules beyond the second hydration shell are involved in the hydration process. A combination of scissoring of coordinating water is followed by synchronized secondary motion of water in the second solvation shell due to hydrogen bonding. Consequently, the water beyond the second solvation penetrates through the second solvation shell and coordinates to the Al3+. Our study reveals a novel microscopic understanding of solvation dynamics for the trivalent ion.

Diffusivity of CO2 in H2O: A Review of Experimental Studies and Molecular Simulations in the Bulk and in Confinement

An in-depth review of the available experimental and molecular simulation studies of CO2 diffusion in H2O, which is a central property in important industrial and environmental processes, such as carbon capture and storage, enhanced oil recovery, and in the food industry is presented. The cases of both bulk and confined systems are covered. The experimental and molecular simulation data gathered are analyzed, and simple and computationally efficient correlations are devised. These correlations are applicable to conditions from 273 K and 0.1 MPa up to 473 K and 45 MPa. The available experimental data for diffusion coefficients of CO2 in brines are also collected, and their dependency on temperature, pressure, and salinity is examined in detail. Other engineering models and correlations reported in literature are also presented. The review of the simulation studies focuses on the force field combinations, the data for diffusivities at low and high pressures, finite-size effects, and the correlations developed based on the Molecular Dynamics data. Regarding the confined systems, we review the main methods to measure and compute the diffusivity of confined CO2 and discuss the main natural and artificial confining media (i.e., smectites, calcites, silica, MOFs, and carbon materials). Detailed discussion is provided regarding the driving force for diffusion of CO2 and H2O under confinement, and on the role of effects such as H2O adsorption on hydrophilic confining media on the diffusivity of CO2. Finally, an outlook of future research paths for advancing the field of CO2 diffusivity in H2O at the bulk phase and in confinement is laid out.

Force field for halide and alkali ions in water based on single-ion and ion-pair thermodynamic properties for a wide range of concentrations

A classical non-polarizable force field for the common halide (F-, Cl-, Br-, and I-) and alkali (Li+, Na+, K+, and Cs+) ions in SPC/E water is presented. This is an extension of the force field developed by Loche et al. for Na+, K+, Cl-, and Br- (JPCB 125, 8581-8587, 2021): in the present work, we additionally optimize Lennard-Jones parameters for Li+, I-, Cs+, and F- ions. Li+ and F- are particularly challenging ions to model due to their small size. The force field is optimized with respect to experimental solvation free energies and activity coefficients, which are the necessary and sufficient quantities to accurately reproduce the electrolyte thermodynamics. Good agreement with experimental reference data is achieved for a wide range of concentrations (up to 4 mol/l). We find that standard Lorentz-Berthelot combination rules are sufficient for all ions except F-, for which modified combination rules are necessary. With the optimized parameters, we show that, although the force field is only optimized based on thermodynamic properties, structural properties are reproduced quantitatively, while ion diffusion coefficients are in qualitative agreement with experimental values.

Importance of Polarizable Embedding for Computing Optical Rotation: The Case of Camphor in Ethanol

Solvation effects on optical rotation are notoriously challenging to model for computational chemistry, as the specific rotatory power of a molecule can vary wildly going from apolar to polar or even protic solvents. To address such a problem, we present a polarizable embedding implementation of an electric and magnetic response property based on density functional theory and the AMOEBA polarizable force field, and apply such an implementation to the study of the optical rotation of camphor in ethanol. By comparing a continuum model, and electrostatic and polarizable embedding QM/MM models, we observe that accounting for the environment's polarization gives rise to not only a different quantitative prediction, in very good agreement with experiments for the QM/AMOEBA model, but also to a very different qualitative picture, with the values of the optical rotation computed along a classical molecular dynamics trajectory with electrostatic embedding being statistically uncorrelated to the ones obtained with the polarizable description.

Unraveling the mechanisms of triplet state formation in a heavy-atom free photosensitizer

Triplet excited state generation plays a pivotal role in photosensitizers, however the reliance on transition metals and heavy atoms can limit the utility of these systems. In this study, we demonstrate that an interplay of competing quantum effects controls the high triplet quantum yield in a prototypical boron dipyrromethene-anthracene (BD-An) donor-acceptor dyad photosensitizer, which is only captured by an accurate treatment of both inner and outer sphere reorganization energies. Our ab initio-derived model provides excellent agreement with experimentally measured spectra, triplet yields and excited state kinetic data, including the triplet lifetime. We find that rapid triplet state formation occurs primarily via high-energy triplet states through both spin-orbit coupled charge transfer and El-Sayed's rule breaking intersystem crossing, rather than direct spin-orbit coupled charge transfer to the lowest lying triplet state. Our calculations also reveal that competing effects of nuclear tunneling, electronic state recrossing, and electronic polarizability dictate the rate of non-productive ground state recombination. This study sheds light on the quantum effects driving efficient triplet formation in the BD-An system, and offers a promising simulation methodology for diverse photochemical systems.

Biomolecular dynamics with machine-learned quantum-mechanical force fields trained on diverse chemical fragments

The GEMS method enables molecular dynamics simulations of large heterogeneous systems at ab initio quality.

QMrebind: incorporating quantum mechanical force field reparameterization at the ligand binding site for improved drug-target kinetics through milestoning simulations

Understanding the interaction of ligands with biomolecules is an integral component of drug discovery and development. Challenges for computing thermodynamic and kinetic quantities for pharmaceutically relevant receptor-ligand complexes include the size and flexibility of the ligands, large-scale conformational rearrangements of the receptor, accurate force field parameters, simulation efficiency, and sufficient sampling associated with rare events. Our recently developed multiscale milestoning simulation approach, SEEKR2 (Simulation Enabled Estimation of Kinetic Rates v.2), has demonstrated success in predicting unbinding (koff) kinetics by employing molecular dynamics (MD) simulations in regions closer to the binding site. The MD region is further subdivided into smaller Voronoi tessellations to improve the simulation efficiency and parallelization. To date, all MD simulations are run using general molecular mechanics (MM) force fields. The accuracy of calculations can be further improved by incorporating quantum mechanical (QM) methods into generating system-specific force fields through reparameterizing ligand partial charges in the bound state. The force field reparameterization process modifies the potential energy landscape of the bimolecular complex, enabling a more accurate representation of the intermolecular interactions and polarization effects at the bound state. We present QMrebind (Quantum Mechanical force field reparameterization at the receptor-ligand binding site), an ORCA-based software that facilitates reparameterizing the potential energy function within the phase space representing the bound state in a receptor-ligand complex. With SEEKR2 koff estimates and experimentally determined kinetic rates, we compare and interpret the receptor-ligand unbinding kinetics obtained using the newly reparameterized force fields for model host-guest systems and HSP90-inhibitor complexes. This method provides an opportunity to achieve higher accuracy in predicting receptor-ligand koff rate constants.

MBX: A many-body energy and force calculator for data-driven many-body simulations

Many-Body eXpansion (MBX) is a C++ library that implements many-body potential energy functions (PEFs) within the "many-body energy" (MB-nrg) formalism. MB-nrg PEFs integrate an underlying polarizable model with explicit machine-learned representations of many-body interactions to achieve chemical accuracy from the gas to the condensed phases. MBX can be employed either as a stand-alone package or as an energy/force engine that can be integrated with generic software for molecular dynamics and Monte Carlo simulations. MBX is parallelized internally using Open Multi-Processing and can utilize Message Passing Interface when available in interfaced molecular simulation software. MBX enables classical and quantum molecular simulations with MB-nrg PEFs, as well as hybrid simulations that combine conventional force fields and MB-nrg PEFs, for diverse systems ranging from small gas-phase clusters to aqueous solutions and molecular fluids to biomolecular systems and metal-organic frameworks.

Drude Polarizable Lipid Force Field with Explicit Treatment of Long-Range Dispersion: Parametrization and Validation for Saturated and Monounsaturated Zwitterionic Lipids

Accurate empirical force fields of lipid molecules are a critical component of molecular dynamics simulation studies aimed at investigating properties of monolayers, bilayers, micelles, vesicles, and liposomes, as well as heterogeneous systems, such as protein-membrane complexes, bacterial cell walls, and more. While the majority of lipid force field-based simulations have been performed using pairwise-additive nonpolarizable models, advances have been made in the development of the polarizable force field based on the classical Drude oscillator model. In the present study, we undertake further optimization of the Drude lipid force field, termed Drude2023, including improved treatment of the phosphate and glycerol linker region of PC and PE headgroups, additional optimization of the alkene group in monounsaturated lipids, and inclusion of long-range Lennard-Jones interactions using the particle-mesh Ewald method. Initial optimization targeted quantum mechanical (QM) data on small model compounds representative of the linker region. Subsequent optimization targeted QM data on larger model compounds, experimental data, and dihedral potentials of mean force from the CHARMM36 additive lipid force field using a parameter reweighting protocol. The use of both experimental and QM target data during the reweighting protocol is shown to produce physically reasonable parameters that reproduce a collection of experimental observables. Target data for optimization included surface area/lipid for DPPC, DSPC, DMPC, and DLPC bilayers and nuclear magnetic resonance (NMR) order parameters for DPPC bilayers. Validation data include prediction of membrane thickness, scattering form factors, electrostatic potential profiles, compressibility moduli, surface area per lipid, water permeability, NMR T1 relaxation times, diffusion constants, and monolayer surface tensions for a variety of saturated and unsaturated lipid mono- and bilayers. Overall, the agreement with experimental data is quite good, though the results are less satisfactory for the NMR T1 relaxation times for carbons near the ester groups. Notable improvements compared to the additive C36 force field were obtained for membrane dipole potentials, lipid diffusion coefficients, and water permeability with the exception of monounsaturated lipid bilayers. It is anticipated that the optimized polarizable Drude2023 force field will help generate more accurate molecular simulations of pure bilayers and heterogeneous systems containing membranes, advancing our understanding of the role of electronic polarization in these systems.

Symmetrized Drude Oscillator Force Fields Improve Numerical Performance of Polarizable Molecular Dynamics

Drude oscillator potentials are a popular and computationally efficient class of polarizable models that represent each polarizable atom as a positively charged Drude core harmonically bound to a negatively charged Drude shell. We show that existing force fields that place all non-Coulomb forces on the Drude core and none on the shell inadvertently couple the dipole to non-Coulombic forces. This introduces errors where interactions with neutral particles can erroneously induce atomic polarization, leading to spurious polarizations in the absence of an electric field, exacerbating violations of equipartition in the employed Carr-Parinello scheme. A suitable symmetrization of the interaction potential that correctly splits the force between the Drude core and shell can correct this shortcoming, improving the stability and numerical performance of Drude oscillator-based simulations. The symmetrization procedure is straightforward and only requires the rescaling of a few force field parameters.

Gaussian Process Regression Models for Predicting Atomic Energies and Multipole Moments

Developing a force field is a difficult task because its design is typically pulled in opposite directions by speed and accuracy. FFLUX breaks this trend by utilizing Gaussian process regression (GPR) to predict, at ab initio accuracy, atomic energies and multipole moments as obtained from the quantum theory of atoms in molecules (QTAIM). This work demonstrates that the in-house FFLUX training pipeline can generate successful GPR models for six representative molecules: peptide-capped glycine and alanine, glucose, paracetamol, aspirin, and ibuprofen. The molecules were sufficiently distorted to represent configurations from an AMBER-GAFF2 molecular dynamics run. All internal degrees of freedom were covered corresponding to 93 dimensions in the case of the largest molecule ibuprofen (33 atoms). Benefiting from active learning, the GPR models contain only about 2000 training points and return largely sub-kcal mol^-1 prediction errors for the validation sets. A proof of concept has been reached for transferring the model produced through active learning on one atomic property to that of the remaining atomic properties. The prediction of electrostatic interaction can be assessed at the intermolecular level, and the vast majority of interactions have a root-mean-square error of less than 0.1 kJ mol^-1 with a maximum value of ∼1 kJ mol^-1 for a glycine and paracetamol dimer.

Regularized by Physics: Graph Neural Network Parametrized Potentials for the Description of Intermolecular Interactions

Simulations of molecular systems using electronic structure methods are still not feasible for many systems of biological importance. As a result, empirical methods such as force fields (FF) have become an established tool for the simulation of large and complex molecular systems. The parametrization of FF is, however, time-consuming and has traditionally been based on experimental data. Recent years have therefore seen increasing efforts to automatize FF parametrization or to replace FF with machine-learning (ML) based potentials. Here, we propose an alternative strategy to parametrize FF, which makes use of ML and gradient-descent based optimization while retaining a functional form founded in physics. Using a predefined functional form is shown to enable interpretability, robustness, and efficient simulations of large systems over long time scales. To demonstrate the strength of the proposed method, a fixed-charge and a polarizable model are trained on ab initio potential-energy surfaces. Given only information about the constituting elements, the molecular topology, and reference potential energies, the models successfully learn to assign atom types and corresponding FF parameters from scratch. The resulting models and parameters are validated on a wide range of experimentally and computationally derived properties of systems including dimers, pure liquids, and molecular crystals.

State averaged CASSCF in AMOEBA polarizable water model for simulating nonadiabatic molecular dynamics with nonequilibrium solvation effects

This paper presents a state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular application (AMOEBA) polarizable water model, which enables rigorous simulation of non-adiabatic molecular dynamics with nonequilibrium solvation effects. The molecular orbital and configuration interaction coefficients of the solute wavefunction, and the induced dipoles on solvent atoms, are solved by minimizing the state averaged energy variationally. In particular, by formulating AMOEBA water models and the polarizable continuum model (PCM) in a unified way, the algorithms developed for computing SA-CASSCF/PCM energies, analytical gradients, and non-adiabatic couplings in our previous work can be generalized to SA-CASSCF/AMOEBA by properly substituting a specific list of variables. Implementation of this method will be discussed with the emphasis on how the calculations of different terms are partitioned between the quantum chemistry and molecular mechanics codes. We will present and discuss results that demonstrate the accuracy and performance of the implementation. Next, we will discuss results that compare three solvent models that work with SA-CASSCF, i.e., PCM, fixed-charge force fields, and the newly implemented AMOEBA. Finally, the new SA-CASSCF/AMOEBA method has been interfaced with the ab initio multiple spawning method to carry out non-adiabatic molecular dynamics simulations. This method is demonstrated by simulating the photodynamics of the model retinal protonated Schiff base molecule in water.

Treating Polarization Effects in Charged and Polar Bio-Molecules Through Variable Electrostatic Parameters

Polarization plays important roles in charged and hydrogen bonding containing systems. Much effort ranging from the construction of physics-based models to quantum mechanism (QM)-based and machine learning (ML)-assisted models have been devoted to incorporating the polarization effect into the conventional force fields at different levels, such as atomic and coarse grained (CG). The application of polarizable force fields or polarization models was limited by two aspects, namely, computational cost and transferability. Different from physics-based models, no predetermining parameters were required in the QM-based approaches. Taking advantage of both the accuracy of QM calculations and efficiency of molecular mechanism (MM) and ML, polarization effects could be treated more efficiently while maintaining the QM accuracy. The computational cost could be reduced with variable electrostatic parameters, such as the charge, dipole, and electronic dielectric constant with the help of linear scaling fragmentation-based QM calculations and ML models. Polarization and entropy effects on the prediction of partition coefficient of druglike molecules are demonstrated by using both explicit or implicit all-atom molecular dynamics simulations and machine learning-assisted models. Directions and challenges for future development are also envisioned.

Computer Aided Development of Nucleic Acid Applications in Nanotechnologies

Utilization of nucleic acids (NAs) in nanotechnologies and nanotechnology-related applications is a growing field with broad application potential, ranging from biosensing up to targeted cell delivery. Computer simulations are useful techniques that can aid design and speed up development in this field. This review focuses on computer simulations of hybrid nanomaterials composed of NAs and other components. Current state-of-the-art molecular dynamics simulations, empirical force fields (FFs), and coarse-grained approaches for the description of deoxyribonucleic acid and ribonucleic acid are critically discussed. Challenges in combining biomacromolecular and nanomaterial FFs are emphasized. Recent applications of simulations for modeling NAs and their interactions with nano- and biomaterials are overviewed in the fields of sensing applications, targeted delivery, and NA templated materials. Future perspectives of development are also highlighted.

Application of Computational Biology and Artificial Intelligence in Drug Design

Traditional drug design requires a great amount of research time and developmental expense. Booming computational approaches, including computational biology, computer-aided drug design, and artificial intelligence, have the potential to expedite the efficiency of drug discovery by minimizing the time and financial cost. In recent years, computational approaches are being widely used to improve the efficacy and effectiveness of drug discovery and pipeline, leading to the approval of plenty of new drugs for marketing. The present review emphasizes on the applications of these indispensable computational approaches in aiding target identification, lead discovery, and lead optimization. Some challenges of using these approaches for drug design are also discussed. Moreover, we propose a methodology for integrating various computational techniques into new drug discovery and design.

Inclusion of High-Field Target Data in AMOEBA's Calibration Improves Predictions of Protein-Ion Interactions

The reliability of molecular mechanics simulations to predict effects of ion binding to proteins depends on their ability to simultaneously describe ion-protein, ion-water, and protein-water interactions. Force fields (FFs) to describe protein-water and ion-water interactions have been constructed carefully and have also been refined routinely to improve accuracy. Descriptions for ion-protein interactions have also been refined, although in an a posteriori manner through the use of "nonbonded-fix (NB-fix)" approaches in which parameters from default Lennard-Jones mixing rules are replaced with those optimized against some reference data. However, even after NB-fix corrections, there remains a significant need for improvement. This is also true for polarizable FFs that include self-consistent inducible moments. Our recent studies on the polarizable AMOEBA FF suggested that the problem associated with modeling ion-protein interactions could be alleviated by recalibrating polarization models of cation-coordinating functional groups so that they respond better to the high electric fields present near ions. Here, we present such a recalibration of carbonyls, carboxylates, and hydroxyls in the AMOEBA protein FF and report that it does improve predictions substantially─mean absolute errors in Na⁺-protein and K⁺-protein interaction energies decrease from 8.7 to 5.3 and 9.6 to 6.3 kcal/mol, respectively. Errors are computed with respect to estimates from van der Waals-inclusive density functional theory benchmarked against high-level quantum mechanical calculations and experiments. While recalibration does improve ion-protein interaction energies, they still remain underestimated, suggesting that further improvements can be made in a systematic manner through modifications in classical formalism. Nevertheless, we show that by applying our many-body NB-fix correction to Lennard-Jones components, these errors are further reduced to 2.7 and 2.6 kcal/mol, respectively, for Na⁺ and K⁺ ions. Finally, we show that the recalibrated AMOEBA protein FF retains its intrinsic reliability in predicting protein structure and dynamics in the condensed phase.

Classical Exchange Polarization: An Anisotropic Variable Polarizability Model

Polarizability, or the tendency of the electron distribution to distort under an electric field, often depends on the local chemical environment. For example, the polarizability of a chloride ion is larger in gas phase compared to a chloride ion solvated in water. This effect is due to the restriction the Pauli exclusion principle places on the allowed electron states. Because no two electrons can occupy the same state, when a highly polarizable atom comes in close contact with other atoms or molecules, the space of allowed states can dramatically decrease. This constraint suggests that an accurate molecular mechanics polarizability model should depend on the radial distance between neighboring atoms. This paper introduces a variable polarizability model within the framework of the HIPPO (Hydrogen-like Intermolecular Polarizable Potential) force field, by damping the polarizability as a function of the orbital overlap of two atoms. This effectively captures the quantum mechanical exchange polarization effects, without explicit utilization of antisymmetrized wave functions. We show that the variable polarizability model remarkably improves the two-body polarization energies and three-body energies of ion-ion and ion-water systems. Under this model, no manual tuning of atomic polarizabilities for monatomic ions is required; the gas-phase polarizability can be used because an appropriate damping function is able to correct the polarizability at short range.

Extension of the CHARMM Classical Drude Polarizable Force Field to N- and O-Linked Glycopeptides and Glycoproteins

Molecular dynamic simulations are an effective tool to study complex molecular systems and are contingent upon the availability of an accurate and reliable molecular mechanics force field. The Drude polarizable force field, which allows for the explicit treatment of electronic polarization in a computationally efficient fashion, has been shown to reproduce experimental properties that were difficult or impossible to reproduce with the CHARMM additive force field, including peptide folding cooperativity, RNA hairpin structures, and DNA base flipping. Glycoproteins are essential components of glycoconjugate vaccines, antibodies, and many pharmaceutically important molecules, and an accurate polarizable force field that includes compatibility between the protein and carbohydrate aspect of the force field is essential to study these types of systems. In this work, we present an extension of the Drude polarizable force field to glycoproteins, including both N- and O-linked species. Parameter optimization focused on the dihedral terms using a reweighting protocol targeting NMR solution J-coupling data for model glycopeptides. Validation of the model include eight model glycopeptides and four glycoproteins with multiple N- and O-linked glycosylations. The new glycoprotein carbohydrate force field can be used in conjunction with the remainder of Drude polarizable force field through a variety of MD simulation programs including GROMACS, OPENMM, NAMD, and CHARMM and may be accessed through the Drude Prepper module in the CHARMM-GUI.

Carbon Nanodots from an In Silico Perspective

Carbon nanodots (CNDs) are the latest and most shining rising stars among photoluminescent (PL) nanomaterials. These carbon-based surface-passivated nanostructures compete with other related PL materials, including traditional semiconductor quantum dots and organic dyes, with a long list of benefits and emerging applications. Advantages of CNDs include tunable inherent optical properties and high photostability, rich possibilities for surface functionalization and doping, dispersibility, low toxicity, and viable synthesis (top-down and bottom-up) from organic materials. CNDs can be applied to biomedicine including imaging and sensing, drug-delivery, photodynamic therapy, photocatalysis but also to energy harvesting in solar cells and as LEDs. More applications are reported continuously, making this already a research field of its own. Understanding of the properties of CNDs requires one to go to the levels of electrons, atoms, molecules, and nanostructures at different scales using modern molecular modeling and to correlate it tightly with experiments. This review highlights different in silico techniques and studies, from quantum chemistry to the mesoscale, with particular reference to carbon nanodots, carbonaceous nanoparticles whose structural and photophysical properties are not fully elucidated. The role of experimental investigation is also presented. Hereby, we hope to encourage the reader to investigate CNDs and to apply virtual chemistry to obtain further insights needed to customize these amazing systems for novel prospective applications.

Consistent Picture of Phosphate-Divalent Cation Binding from Models with Implicit and Explicit Electronic Polarization

The binding of divalent cations to the ubiquitous phosphate group is essential for a number of key biological processes, such as DNA compaction, RNA folding, or interactions of some proteins with membranes. Yet, probing their binding sites, modes, and associated binding free energy is a challenge for both experiments and simulations. In simulations, standard force fields strongly overestimate the interaction between phosphate groups and divalent cations. Here, we examine how different strategies to include electronic polarization effects in force fields─implicitly, through the use of scaled charges or pair-specific Lennard-Jones parameters, or explicitly, with the polarizable force fields Drude and AMOEBA─capture the interactions of a model phosphate compound, dimethyl phosphate, with calcium and magnesium divalent cations. We show that both implicit and explicit approaches, when carefully parameterized, are successful in capturing the overall binding free energy and that common trends emerge from the comparison of different simulation approaches. Overall, the binding is very moderate, slightly weaker for Ca²⁺ than Mg²⁺, and the solvent-shared ion pair is slightly more stable than the contact monodentate ion pair. The bidentate ion pair is higher in energy (or even fully unstable for Mg²⁺). Our results thus suggest practical ways to capture the divalent cations with biomolecular phosphate groups in complex biochemical systems. In particular, the computational efficiency of implicit models makes them ideally suited for large-scale simulations of biological assemblies, with improved accuracy compared to state-of-the-art fixed-charge force fields.

Improving the Accuracy of Atomistic Simulations of the Electrochemical Interface

Atomistic simulation of the electrochemical double layer is an ambitious undertaking, requiring quantum mechanical description of electrons, phase space sampling of liquid electrolytes, and equilibration of electrolytes over nanosecond time scales. All models of electrochemistry make different trade-offs in the approximation of electrons and atomic configurations, from the extremes of classical molecular dynamics of a complete interface with point-charge atoms to correlated electronic structure methods of a single electrode configuration with no dynamics or electrolyte. Here, we review the spectrum of simulation techniques suitable for electrochemistry, focusing on the key approximations and accuracy considerations for each technique. We discuss promising approaches, such as enhanced sampling techniques for atomic configurations and computationally efficient beyond density functional theory (DFT) electronic methods, that will push electrochemical simulations beyond the present frontier.

Stress tensor and constant pressure simulation for polarizable Gaussian multipole model

Our previous article has established the theory of molecular dynamics (MD) simulations for systems modeled with the polarizable Gaussian multipole (pGM) electrostatics [Wei et al., J. Chem. Phys. 153(11), 114116 (2020)]. Specifically, we proposed the covalent basis vector framework to define the permanent multipoles and derived closed-form energy and force expressions to facilitate an efficient implementation of pGM electrostatics. In this study, we move forward to derive the pGM internal stress tensor for constant pressure MD simulations with the pGM electrostatics. Three different formulations are presented for the flexible, rigid, and short-range screened systems, respectively. The analytical formulations were implemented in the SANDER program in the Amber package and were first validated with the finite-difference method for two different boxes of pGM water molecules. This is followed by a constant temperature and constant pressure MD simulation for a box of 512 pGM water molecules. Our results show that the simulation system stabilized at a physically reasonable state and maintained the balance with the externally applied pressure. In addition, several fundamental differences were observed between the pGM and classic point charge models in terms of the simulation behaviors, indicating more extensive parameterization is necessary to utilize the pGM electrostatics.

Deep Neural Network Model to Predict the Electrostatic Parameters in the Polarizable Classical Drude Oscillator Force Field

The Drude polarizable force field (FF) captures electronic polarization effects via auxiliary Drude particles that are attached to non-hydrogen atoms, distinguishing it from commonly used additive FFs that rely on fixed charges. The Drude FF currently includes parameters for biomolecules such as proteins, nucleic acids, lipids, and carbohydrates and small-molecule representative of those classes of molecules as well as a range of atomic ions. Extension of the Drude FF to novel small druglike molecules is challenging as it requires the assignment of partial charges, atomic polarizabilities, and Thole scaling factors. In the present article, deep neural network (DNN) models are trained on quantum mechanical (QM)-based partial charges and atomic polarizabilities along with Thole scale factors trained to target QM molecular dipole moments and polarizabilities. Training of the DNN model used a collection of 39 421 molecules with molecular weights up to 200 Da and containing H, C, N, O, P, S, F, Cl, Br, or I atoms. The DNN model utilizes bond connectivity, including 1,2, 1,3, 1,4, and 1,5 terms and distances of Drude FF atom types as the feature vector to build the model, allowing it to capture both local and nonlocal effects in the molecules. Novel methods have been developed to determine restrained electrostatic potential (RESP) charges on atoms and external points representing lone pairs and to determine Thole scale factors, which have no QM analogue. A penalty scheme is devised as a performance predictor of the trained model. Validation studies show that these DNN models can precisely predict molecular dipole and polarizabilities of Food and Drug Administration (FDA)-approved drugs compared to reference MP2 calculations. The availability of the DNN model allowing for the rapid estimation of the Drude electrostatic parameters will facilitate its applicability to a wider range of molecular species.

Modelling of the dynamic polarizability of macromolecules for single-molecule optical biosensing

The structural dynamics of macromolecules is important for most microbiological processes, from protein folding to the origins of neurodegenerative disorders. Noninvasive measurements of these dynamics are highly challenging. Recently, optical sensors have been shown to allow noninvasive time-resolved measurements of the dynamic polarizability of single-molecules. Here we introduce a method to efficiently predict the dynamic polarizability from the atomic configuration of a given macromolecule. This provides a means to connect the measured dynamic polarizability to the underlying structure of the molecule, and therefore to connect temporal measurements to structural dynamics. To illustrate the methodology we calculate the change in polarizability as a function of time based on conformations extracted from molecular dynamics simulations and using different conformations of motor proteins solved crystalographically. This allows us to quantify the magnitude of the changes in polarizablity due to thermal and functional motions.

Learning Atomic Multipoles: Prediction of the Electrostatic Potential with Equivariant Graph Neural Networks

The accurate description of electrostatic interactions remains a challenging problem for classical potential-energy functions. The commonly used fixed partial-charge approximation fails to reproduce the electrostatic potential at short range due to its insensitivity to conformational changes and anisotropic effects. At the same time, possibly more accurate machine-learned (ML) potentials struggle with the long-range behavior due to their inherent locality ansatz. Employing a multipole expansion offers in principle an exact treatment of the electrostatic potential such that the long-range and short-range electrostatic interactions can be treated simultaneously with high accuracy. However, such an expansion requires the calculation of the electron density using computationally expensive quantum-mechanical (QM) methods. Here, we introduce an equivariant graph neural network (GNN) to address this issue. The proposed model predicts atomic multipoles up to the quadrupole, circumventing the need for expensive QM computations. By using an equivariant architecture, the model enforces the correct symmetry by design without relying on local reference frames. The GNN reproduces the electrostatic potential of various systems with high fidelity. Possible uses for such an approach include the separate treatment of long-range interactions in ML potentials, the analysis of electrostatic potential surfaces, and static multipoles in polarizable force fields.

CHARMM-GUI Drude prepper for molecular dynamics simulation using the classical Drude polarizable force field

Explicit treatment of electronic polarizability in empirical force fields (FFs) represents an extension over a traditional additive or pairwise FF and provides a more realistic model of the variations in electronic structure in condensed phase, macromolecular simulations. To facilitate utilization of the polarizable FF based on the classical Drude oscillator model, Drude Prepper has been developed in CHARMM-GUI. Drude Prepper ingests additive CHARMM protein structures file (PSF) and pre-equilibrated coordinates in CHARMM, PDB, or NAMD format, from which the molecular components of the system are identified. These include all residues and patches connecting those residues along with water, ions, and other solute molecules. This information is then used to construct the Drude FF-based PSF using molecular generation capabilities in CHARMM, followed by minimization and equilibration. In addition, inputs are generated for molecular dynamics (MD) simulations using CHARMM, GROMACS, NAMD, and OpenMM. Validation of the Drude Prepper protocol and inputs is performed through conversion and MD simulations of various heterogeneous systems that include proteins, nucleic acids, lipids, polysaccharides, and atomic ions using the aforementioned simulation packages. Stable simulations are obtained in all studied systems, including 5 μs simulation of ubiquitin, verifying the integrity of the generated Drude PSFs. In addition, the ability of the Drude FF to model variations in electronic structure is shown through dipole moment analysis in selected systems. The capabilities and availability of Drude Prepper in CHARMM-GUI is anticipated to greatly facilitate the application of the Drude FF to a range of condensed phase, macromolecular systems.

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.

Understanding the avidin-biotin binding based on polarized protein-specific charge

In this study, charge updating schemes based on the local polarized protein-specific charge (LPPC) were introduced to vary the atomic charges of the biotin molecule and the residues in close contact during the simulation of the avidin-biotin complexes. The need of the charge variation of the ligand in response to changes in its surroundings was thoroughly studied. The results show that the calculated binding energy difference between biotin (BTN1) and 2'-iminobiotin (BTN2) and avidin is in excellent agreement with the experimental value, thus verifying the feasibility of updating the atomic charges of ligands during the simulation.

MB-Fit: Software infrastructure for data-driven many-body potential energy functions

Many-body potential energy functions (MB-PEFs), which integrate data-driven representations of many-body short-range quantum mechanical interactions with physics-based representations of many-body polarization and long-range interactions, have recently been shown to provide high accuracy in the description of molecular interactions from the gas to the condensed phase. Here, we present MB-Fit, a software infrastructure for the automated development of MB-PEFs for generic molecules within the TTM-nrg (Thole-type model energy) and MB-nrg (many-body energy) theoretical frameworks. Besides providing all the necessary computational tools for generating TTM-nrg and MB-nrg PEFs, MB-Fit provides a seamless interface with the MBX software, a many-body energy and force calculator for computer simulations. Given the demonstrated accuracy of the MB-PEFs, particularly within the MB-nrg framework, we believe that MB-Fit will enable routine predictive computer simulations of generic (small) molecules in the gas, liquid, and solid phases, including, but not limited to, the modeling of quantum isomeric equilibria in molecular clusters, solvation processes, molecular crystals, and phase diagrams.

Machine Learning Force Fields

In recent years, the use of machine learning (ML) in computational chemistry has enabled numerous advances previously out of reach due to the computational complexity of traditional electronic-structure methods. One of the most promising applications is the construction of ML-based force fields (FFs), with the aim to narrow the gap between the accuracy of ab initio methods and the efficiency of classical FFs. The key idea is to learn the statistical relation between chemical structure and potential energy without relying on a preconceived notion of fixed chemical bonds or knowledge about the relevant interactions. Such universal ML approximations are in principle only limited by the quality and quantity of the reference data used to train them. This review gives an overview of applications of ML-FFs and the chemical insights that can be obtained from them. The core concepts underlying ML-FFs are described in detail, and a step-by-step guide for constructing and testing them from scratch is given. The text concludes with a discussion of the challenges that remain to be overcome by the next generation of ML-FFs.

The phase diagram of carbon dioxide from correlation functions and a many-body potential

The phase stability and equilibria of carbon dioxide are investigated from 125-325 K and 1-10 000 atm using extensive molecular dynamics (MD) simulations and the Two-Phase Thermodynamics (2PT) method. We devise a direct approach for calculating phase diagrams, in general, by considering the separate chemical potentials of the isolated phase at specific points on the P-T diagram. The unique ability of 2PT to accurately and efficiently approximate the entropy and Gibbs energy of liquids allows for assignment of phase boundaries from relatively short (∼100 ps) MD simulations. We validate our approach by calculating the critical properties of the flexible elementary physical model 2, showing good agreement with previous results. We show, however, that the incorrect description of the short-range Pauli force and the lack of molecular charge polarization lead to deviations from experiments at high pressures. We, thus, develop a many-body, fluctuating charge model for CO₂, termed CO₂-Fq, from high level quantum mechanics (QM) calculations that accurately capture the condensed phase vibrational properties of the solid (including the Fermi resonance at 1378 cm^-1) as well as the diffusional properties of the liquid, leading to overall excellent agreement with experiments over the entire phase diagram. This work provides an efficient computational approach for determining phase diagrams of arbitrary systems and underscores the critical role of QM charge reorganization physics in molecular phase stability.

An enhanced sampling QM/AMOEBA approach: The case of the excited state intramolecular proton transfer in solvated 3-hydroxyflavone

We present an extension of the polarizable quantum mechanical (QM)/AMOEBA approach to enhanced sampling techniques. This is achieved by connecting the enhanced sampling PLUMED library to the machinery based on the interface of Gaussian and Tinker to perform QM/AMOEBA molecular dynamics. As an application, we study the excited state intramolecular proton transfer of 3-hydroxyflavone in two solvents: methanol and methylcyclohexane. By using a combination of molecular dynamics and umbrella sampling, we find an ultrafast component of the transfer, which is common to the two solvents, and a much slower component, which is active in the protic solvent only. The mechanisms of the two components are explained in terms of intramolecular vibrational redistribution and intermolecular hydrogen-bonding, respectively. Ground and excited state free energies along an effective reaction coordinate are finally obtained allowing for a detailed analysis of the solvent mediated mechanism.

Efficient formulation of polarizable Gaussian multipole electrostatics for biomolecular simulations

Molecular dynamics simulations of biomolecules have been widely adopted in biomedical studies. As classical point-charge models continue to be used in routine biomolecular applications, there have been growing demands on developing polarizable force fields for handling more complicated biomolecular processes. Here, we focus on a recently proposed polarizable Gaussian Multipole (pGM) model for biomolecular simulations. A key benefit of pGM is its screening of all short-range electrostatic interactions in a physically consistent manner, which is critical for stable charge-fitting and is needed to reproduce molecular anisotropy. Another advantage of pGM is that each atom's multipoles are represented by a single Gaussian function or its derivatives, allowing for more efficient electrostatics than other Gaussian-based models. In this study, we present an efficient formulation for the pGM model defined with respect to a local frame formed with a set of covalent basis vectors. The covalent basis vectors are chosen to be along each atom's covalent bonding directions. The new local frame can better accommodate the fact that permanent dipoles are primarily aligned along covalent bonds due to the differences in electronegativity of bonded atoms. It also allows molecular flexibility during molecular simulations and facilitates an efficient formulation of analytical electrostatic forces without explicit torque computation. Subsequent numerical tests show that analytical atomic forces agree excellently with numerical finite-difference forces for the tested system. Finally, the new pGM electrostatics algorithm is interfaced with the particle mesh Ewald (PME) implementation in Amber for molecular simulations under the periodic boundary conditions. To validate the overall pGM/PME electrostatics, we conducted an NVE simulation for a small water box of 512 water molecules. Our results show that to achieve energy conservation in the polarizable model, it is important to ensure enough accuracy on both PME and induction iteration. It is hoped that the reformulated pGM model will facilitate the development of future force fields based on the pGM electrostatics for applications in biomolecular systems and processes where polarization plays crucial roles.

Multiscale Models for Light-Driven Processes

Multiscale models combining quantum mechanical and classical descriptions are a very popular strategy to simulate properties and processes of complex systems. Many alternative formulations have been developed, and they are now available in all of the most widely used quantum chemistry packages. Their application to the study of light-driven processes, however, is more recent, and some methodological and numerical problems have yet to be solved. This is especially the case for the polarizable formulation of these models, the recent advances in which we review here. Specifically, we identify and describe the most important specificities that the polarizable formulation introduces into both the simulation of excited-state dynamics and the modeling of excitation energy and electron transfer processes.

A practical guide to biologically relevant molecular simulations with charge scaling for electronic polarization

Molecular simulations can elucidate atomistic-level mechanisms of key biological processes, which are often hardly accessible to experiment. However, the results of the simulations can only be as trustworthy as the underlying simulation model. In many of these processes, interactions between charged moieties play a critical role. Current empirical force fields tend to overestimate such interactions, often in a dramatic way, when polyvalent ions are involved. The source of this shortcoming is the missing electronic polarization in these models. Given the importance of such biomolecular systems, there is great interest in fixing this deficiency in a computationally inexpensive way without employing explicitly polarizable force fields. Here, we review the electronic continuum correction approach, which accounts for electronic polarization in a mean-field way, focusing on its charge scaling variant. We show that by pragmatically scaling only the charged molecular groups, we qualitatively improve the charge-charge interactions without extra computational costs and benefit from decades of force field development on biomolecular force fields.

Extended Lagrangian Born-Oppenheimer molecular dynamics for orbital-free density-functional theory and polarizable charge equilibration models

Extended Lagrangian Born-Oppenheimer molecular dynamics (XL-BOMD) [A. M. N. Niklasson, Phys. Rev. Lett. 100, 123004 (2008)] is formulated for orbital-free Hohenberg-Kohn density-functional theory and for charge equilibration and polarizable force-field models that can be derived from the same orbital-free framework. The purpose is to introduce the most recent features of orbital-based XL-BOMD to molecular dynamics simulations based on charge equilibration and polarizable force-field models. These features include a metric tensor generalization of the extended harmonic potential, preconditioners, and the ability to use only a single Coulomb summation to determine the fully equilibrated charges and the interatomic forces in each time step for the shadow Born-Oppenheimer potential energy surface. The orbital-free formulation has a charge-dependent, short-range energy term that is separate from long-range Coulomb interactions. This enables local parameterizations of the short-range energy term, while the long-range electrostatic interactions can be treated separately. The theory is illustrated for molecular dynamics simulations of an atomistic system described by a charge equilibration model with periodic boundary conditions. The system of linear equations that determines the equilibrated charges and the forces is diagonal, and only a single Ewald summation is needed in each time step. The simulations exhibit the same features in accuracy, convergence, and stability as are expected from orbital-based XL-BOMD.

Uncovering a Universal Molecular Mechanism of Salt Ion Adsorption at Solid/Water Interfaces

Solid/water interfaces, in which salt ions come in close proximity to solids, are ubiquitous in nature. Because water is a polar solvent and salt ions are charged, a long-standing puzzle involving solid/water interfaces is how do the electric fields exerted by the salt ions and the interfacial water molecules polarize the charge distribution in the solid and how does this polarization, in turn, influence ion adsorption at any solid/water interface. Here, using state-of-the-art polarizable force fields derived from quantum chemical simulations, we perform all-atomistic molecular dynamics simulations to investigate the adsorption of various ions comprising the well-known Hofmeister series at the graphene/water interface, including comparing with available experimental data. Our findings reveal that, in vacuum, the ionic electric field-induced polarization of graphene results in a significantly large graphene-ion polarization energy, which drives all salt ions to adsorb to graphene. On the contrary, in the presence of water molecules, we show that the ions and the water molecules exert waves of molecular electric fields on graphene which destructively interfere with each other. This remarkable phenomenon is shown to cause a water-mediated screening of more than 85% of the graphene-ion polarization energy. Finally, by investigating superhydrophilic and superhydrophobic model surfaces, we demonstrate that this phenomenon occurs universally at all solid/water interfaces and results in a significant weakening of the ion-solid interactions, such that ion specific effects are governed primarily by a competition between the ion-water and water-water interactions, irrespective of the nature of the solid/water interface.