Tinker-HP: Accelerating Molecular Dynamics Simulations of Large Complex Systems with Advanced Point Dipole Polarizable Force Fields Using GPUs and Multi-GPU Systems

Lambda-ABF: Simplified, Portable, Accurate, and Cost-Effective Alchemical Free-Energy Computation

We introduce the lambda-Adaptive Biasing Force (lambda-ABF) method for the computation of alchemical free-energy differences. We propose a software implementation and showcase it on biomolecular systems. The method arises from coupling multiple-walker adaptive biasing force with λ-dynamics. The sampling of the alchemical variable is continuous and converges toward a uniform distribution, making manual optimization of the λ schedule unnecessary. Contrary to most other approaches, alchemical free-energy estimates are obtained immediately without any postprocessing. Free diffusion of λ improves orthogonal relaxation compared to fixed-λ thermodynamic integration or free-energy perturbation. Furthermore, multiple walkers provide generic orthogonal space coverage with minimal user input and negligible computational overhead. We show that our high-performance implementations coupling the Colvars library with NAMD and Tinker-HP can address real-world cases including ligand-receptor binding with both fixed-charge and polarizable models, with a demonstrably richer sampling than fixed-λ methods. The implementation is fully open-source, publicly available, and readily usable by practitioners of current alchemical methods. Thanks to the portable Colvars library, lambda-ABF presents a unified user interface regardless of the back-end (NAMD, Tinker-HP, or any software to be interfaced in the future), sparing users the effort of learning multiple interfaces. Finally, the Colvars Dashboard extension of the visual molecular dynamics (VMD) software provides an interactive monitoring and diagnostic tool for lambda-ABF simulations.

Incorporating Neural Networks into the AMOEBA Polarizable Force Field

Neural network potentials (NNPs) offer significant promise to bridge the gap between the accuracy of quantum mechanics and the efficiency of molecular mechanics in molecular simulation. Most NNPs rely on the locality assumption that ensures the model's transferability and scalability and thus lack the treatment of long-range interactions, which are essential for molecular systems in the condensed phase. Here we present an integrated hybrid model, AMOEBA+NN, which combines the AMOEBA potential for the short- and long-range noncovalent atomic interactions and an NNP to capture the remaining local covalent contributions. The AMOEBA+NN model was trained on the conformational energy of the ANI-1x data set and tested on several external data sets ranging from small molecules to tetrapeptides. The hybrid model demonstrated substantial improvements over the baseline models in term of accuracy as the molecule size increased, suggesting its potential as a next-generation approach for chemically accurate molecular simulations.

Critical Evaluation of Polarizable and Nonpolarizable Force Fields for Proteins Using Experimentally Derived Nitrile Electric Fields

Molecular dynamics (MD) simulations are frequently carried out for proteins to investigate the role of electrostatics in their biological function. The choice of force field (FF) can significantly alter the MD results, as the simulated local electrostatic interactions lack benchmarking in the absence of appropriate experimental methods. We recently reported that the transition dipole moment (TDM) of the popular nitrile vibrational probe varies linearly with the environmental electric field, overcoming well-known hydrogen bonding (H-bonding) issues for the nitrile frequency and, thus, enabling the unambiguous measurement of electric fields in proteins (J. Am. Chem. Soc. 2022, 144 (17), 7562-7567). Herein, we utilize this new strategy to enable comparisons of experimental and simulated electric fields in protein environments. Specifically, previously determined TDM electric fields exerted onto nitrile-containing o-cyanophenylalanine residues in photoactive yellow protein are compared with MD electric fields from the fixed-charge AMBER FF and the polarizable AMOEBA FF. We observe that the electric field distributions for H-bonding nitriles are substantially affected by the choice of FF. As such, AMBER underestimates electric fields for nitriles experiencing moderate field strengths; in contrast, AMOEBA robustly recapitulates the TDM electric fields. The FF dependence of the electric fields can be partly explained by the presence of additional negative charge density along the nitrile bond axis in AMOEBA, which is due to the inclusion of higher-order multipole parameters; this, in turn, begets more head-on nitrile H-bonds. We conclude by discussing the implications of the FF dependence for the simulation of nitriles and proteins in general.

Enforcing Local DNA Kinks by Sequence-Selective Trisintercalating Oligopeptides of a Tricationic Porphyrin: A Polarizable Molecular Dynamics Study

Bisacridinyl-bisarginyl porphyrin (BABAP) is a trisintercalating derivative of a tricationic porphyrin, formerly designed and synthesized in order to selectively target and photosensitize the ten-base pair palindromic sequence d(CGGGCGCCCG)2 . We resorted to the previously derived (Far et al., 2004) lowest energy-minimized (EM) structure of the BABAP complex with this sequence as a starting point. We performed polarizable molecular dynamics (MD) on this complex. It showed, over a 150 ns duration, the persistent binding of the Arg side-chain on each BABAP arm to the two G bases upstream from the central porphyrin intercalation site. We subsequently performed progressive shortenings of the connector chain linking the Arg-Gly backbone to the acridine, from n=6 methylenes to 4, followed by removal of the Gly backbone and further connector shortenings, from n=4 to n=1. These resulted into progressive deformations ('kinks') of the DNA backbone. In its most accented kinked structure, the DNA backbone was found to have a close overlap with that of DNA bound to Cre recombinase, with, at the level of one acridine intercalation site, negative roll and positive tilt values consistent with those experimentally found for this DNA at its own kinked dinucleotide sequence. Thus, in addition to their photosensitizing properties, some BABAP derivatives could induce sequence-selective, controlled DNA deformations, which are targets for cleavage by endonucleases or for repair enzymes.

Force-field-enhanced neural network interactions: from local equivariant embedding to atom-in-molecule properties and long-range effects

We introduce FENNIX (Force-Field-Enhanced Neural Network InteraXions), a hybrid approach between machine-learning and force-fields. We leverage state-of-the-art equivariant neural networks to predict local energy contributions and multiple atom-in-molecule properties that are then used as geometry-dependent parameters for physically-motivated energy terms which account for long-range electrostatics and dispersion. Using high-accuracy ab initio data (small organic molecules/dimers), we trained a first version of the model. Exhibiting accurate gas-phase energy predictions, FENNIX is transferable to the condensed phase. It is able to produce stable Molecular Dynamics simulations, including nuclear quantum effects, for water predicting accurate liquid properties. The extrapolating power of the hybrid physically-driven machine learning FENNIX approach is exemplified by computing: (i) the solvated alanine dipeptide free energy landscape; (ii) the reactive dissociation of small molecules.

Smooth particle mesh Ewald-integrated stochastic Lanczos many-body dispersion algorithm

We derive and implement an alternative formulation of the Stochastic Lanczos algorithm to be employed in connection with the Many-Body Dispersion model (MBD). Indeed, this formulation, which is only possible due to the Stochastic Lanczos' reliance on matrix-vector products, introduces generalized dipoles and fields. These key quantities allow for a state-of-the-art treatment of periodic boundary conditions via the O(Nlog(N)) Smooth Particle Mesh Ewald (SPME) approach which uses efficient fast Fourier transforms. This SPME-Lanczos algorithm drastically outperforms the standard replica method which is affected by a slow and conditionally convergence rate that limits an efficient and reliable inclusion of long-range periodic boundary conditions interactions in many-body dispersion modelling. The proposed algorithm inherits the embarrassingly parallelism of the original Stochastic Lanczos scheme, thus opening up for a fully converged and efficient periodic boundary conditions treatment of MBD approaches.

Crystallization of FAPbI3: Polytypes and stacking faults

Molecular dynamics simulations are performed to study the crystallization of formamidinium lead iodide. From all-atom simulations of the crystal growth process and the δ-α-phase transitions, we try to reveal the formation of various stack-faulted intermediate defected structures and report various polytypes of formamidinium lead iodide that are observed from simulations.

Enhanced Configurational Sampling Approaches to Alchemical Ligand Binding Free Energy Simulations: Current Status and Challenges

Ligand binding free energy simulations (LB-FES) have been routine tasks in modern drug discovery campaign. A long-standing challenge for LB-FES is the difficulty in adequately sampling nontrivial environmental reorganizations in response to ligand binding. Therefore, various enhanced configurational sampling (ECS) approaches were devised to speed up fluctuations of relevant slow degrees of freedom (SDOF) and ensure simulation convergence. However, in contrast to the achievements in parametrization, software performance, and workflow automation, efficient ECS methodology suitable for high throughput screening remains in an early stage of development. Here, a review of ECS developments with LB-FES is presented, revisiting current approaches and underlining the major technical pitfalls and challenges. This Perspective focuses on alchemical LB-FES on account of their predominant role in high throughput drug screening as well as the established partnership with ECS. The critical aspects of designing ECS approaches, from both theoretical and applied perspectives, are described. This work is intended to provide a contemporary review of the scientific, technical, and practical issues associated with the accelerating convergence of alchemical LB-FES.

AMOEBA Force Field Predicts Accurate Hydrogen Bond Counts of Nitriles in SNase by Revealing Water-Protein Interaction in Vibrational Absorption Frequencies

Precisely quantifying the magnitude and direction of electric fields in proteins has long been an outstanding challenge in understanding biological functions. Nitrile vibrational Stark effect probes have been shown to be minimally disruptive to the protein structure and can be better direct reporters of local electrostatic field in the native state of a protein than other measures such as pKa shifts of titratable residues. However, interpretations of the connection between measured vibrational energy and electric field rely on the accurate molecular understanding of interactions of the nitrile group and its environment, particularly from hydrogen bonding. In this work, we compared the extent of hydrogen bonding calculated in two common force fields, the fixed charge force field Amber03 and polarizable force field AMOEBA, at 10 locations of cyanocysteine (CNC) in staphylococcal nuclease (SNase) against the experimental nitrile absorption frequency in terms of full width at half-maximum (FWHM) and frequency temperature line slope (FTLS). We observed that the number of hydrogen bonds correlated well in AMOEBA trajectories with respect to both the FWHM (r = 0.88) and the FTLS (r = -0.85), whereas the correlation of Amber03 trajectories was less reliable because the Amber03 force field predicted more hydrogen bonds in some mutants. Moreover, we demonstrated that contributions from the interactions between CNC and nearby water molecules were significant in AMOEBA trajectories but were not predicted by Amber03. We conclude that although the nitrile absorption peak shape could be qualitatively predicted by the fixed charge Amber03 force field, the detailed electrostatic environment measured by the nitrile probe in terms of the extent of hydrogen bonding could only be accurately observed in the AMOEBA trajectories, where the permanent dipole, quadrupole, and dipole-induced-dipole polarizable interactions were all taken into account. The significance of this finding to the goal of accurately predicting electric fields in complex biomolecular environments is discussed.

Scalable hybrid deep neural networks/polarizable potentials biomolecular simulations including long-range effects

Deep-HP is a scalable extension of the Tinker-HP multi-GPU molecular dynamics (MD) package enabling the use of Pytorch/TensorFlow Deep Neural Network (DNN) models. Deep-HP increases DNNs' MD capabilities by orders of magnitude offering access to ns simulations for 100k-atom biosystems while offering the possibility of coupling DNNs to any classical (FFs) and many-body polarizable (PFFs) force fields. It allows therefore the introduction of the ANI-2X/AMOEBA hybrid polarizable potential designed for ligand binding studies where solvent-solvent and solvent-solute interactions are computed with the AMOEBA PFF while solute-solute ones are computed by the ANI-2X DNN. ANI-2X/AMOEBA explicitly includes AMOEBA's physical long-range interactions via an efficient Particle Mesh Ewald implementation while preserving ANI-2X's solute short-range quantum mechanical accuracy. The DNN/PFF partition can be user-defined allowing for hybrid simulations to include key ingredients of biosimulation such as polarizable solvents, polarizable counter ions, etc.… ANI-2X/AMOEBA is accelerated using a multiple-timestep strategy focusing on the model's contributions to low-frequency modes of nuclear forces. It primarily evaluates AMOEBA forces while including ANI-2X ones only via correction-steps resulting in an order of magnitude acceleration over standard Velocity Verlet integration. Simulating more than 10 μs, we compute charged/uncharged ligand solvation free energies in 4 solvents, and absolute binding free energies of host-guest complexes from SAMPL challenges. ANI-2X/AMOEBA average errors are discussed in terms of statistical uncertainty and appear in the range of chemical accuracy compared to experiment. The availability of the Deep-HP computational platform opens the path towards large-scale hybrid DNN simulations, at force-field cost, in biophysics and drug discovery.

Routine Molecular Dynamics Simulations Including Nuclear Quantum Effects: From Force Fields to Machine Learning Potentials

We report the implementation of a multi-CPU and multi-GPU massively parallel platform dedicated to the explicit inclusion of nuclear quantum effects (NQEs) in the Tinker-HP molecular dynamics (MD) package. The platform, denoted Quantum-HP, exploits two simulation strategies: the Ring-Polymer Molecular Dynamics (RPMD) that provides exact structural properties at the cost of a MD simulation in an extended space of multiple replicas and the adaptive Quantum Thermal Bath (adQTB) that imposes the quantum distribution of energy on a classical system via a generalized Langevin thermostat and provides computationally affordable and accurate (though approximate) NQEs. We discuss some implementation details, efficient numerical schemes, and parallelization strategies and quickly review the GPU acceleration of our code. Our implementation allows an efficient inclusion of NQEs in MD simulations for very large systems, as demonstrated by scaling tests on water boxes with more than 200,000 atoms (simulated using the AMOEBA polarizable force field). We test the compatibility of the approach with Tinker-HP's recently introduced Deep-HP machine learning potentials module by computing water properties using the DeePMD potential with adQTB thermostatting. Finally, we show that the platform is also compatible with the alchemical free energy estimation capabilities of Tinker-HP and fast enough to perform simulations. Therefore, we study how NQEs affect the hydration free energy of small molecules solvated with the recently developed Q-AMOEBA water force field. Overall, the Quantum-HP platform allows users to perform routine quantum MD simulations of large condensed-phase systems and will help to shed new light on the quantum nature of important interactions in biological matter.

Long-time scale simulations of virus-like particles from three human-norovirus strains

The dynamics of the virus like particles (VLPs) corresponding to the GII.4 Houston, GII.2 SMV, and GI.1 Norwalk strains of human noroviruses (HuNoV) that cause gastroenteritis was investigated by means of long-time (about 30 μs in the laboratory timescale) molecular dynamics simulations with the coarse-grained UNRES force field. The main motion of VLP units turned out to be the bending at the junction between the P1 subdomain (that sits in the VLP shell) and the P2 subdomain (that protrudes outside) of the major VP1 protein, this resulting in a correlated wagging motion of the P2 subdomains with respect to the VLP surface. The fluctuations of the P2 subdomain were found to be more pronounced and the P2 domain made a greater angle with the normal to the VLP surface for the GII.2 strain, which could explain the inability of this strain to bind the histo-blood group antigens (HBGAs).

Predicting efficacy of drug-carrier nanoparticle designs for cancer treatment: a machine learning-based solution

Molecular Dynamic (MD) simulations are very effective in the discovery of nanomedicines for treating cancer, but these are computationally expensive and time-consuming. Existing studies integrating machine learning (ML) into MD simulation to enhance the process and enable efficient analysis cannot provide direct insights without the complete simulation. In this study, we present an ML-based approach for predicting the solvent accessible surface area (SASA) of a nanoparticle (NP), denoting its efficacy, from a fraction of the MD simulations data. The proposed framework uses a time series model for simulating the MD, resulting in an intermediate state, and a second model to calculate the SASA in that state. Empirically, the solution can predict the SASA value 260 timesteps ahead 7.5 times faster with a very low average error of 1956.93. We also introduce the use of an explainability technique to validate the predictions. This work can reduce the computational expense of both processing and data size greatly while providing reliable solutions for the nanomedicine design process.

Δ-Machine Learned Potential Energy Surfaces and Force Fields

There has been great progress in developing machine-learned potential energy surfaces (PESs) for molecules and clusters with more than 10 atoms. Unfortunately, this number of atoms generally limits the level of electronic structure theory to less than the "gold standard" CCSD(T) level. Indeed, for the well-known MD17 dataset for molecules with 9-20 atoms, all of the energies and forces were obtained with DFT calculations (PBE). This Perspective is focused on a Δ-machine learning method that we recently proposed and applied to bring DFT-based PESs to close to CCSD(T) accuracy. This is demonstrated for hydronium, N-methylacetamide, acetyl acetone, and ethanol. For 15-atom tropolone, it appears that special approaches (e.g., molecular tailoring, local CCSD(T)) are needed to obtain the CCSD(T) energies. A new aspect of this approach is the extension of Δ-machine learning to force fields. The approach is based on many-body corrections to polarizable force field potentials. This is examined in detail using the TTM2.1 water potential. The corrections make use of our recent CCSD(T) datasets for 2-b, 3-b, and 4-b interactions for water. These datasets were used to develop a new fully ab initio potential for water, termed q-AQUA.

Predicting molecular properties of α-synuclein using force fields for intrinsically disordered proteins

Independent force field validation is an essential practice to keep track of developments and for performing meaningful Molecular Dynamics simulations. In this work, atomistic force fields for intrinsically disordered proteins (IDP) are tested by simulating the archetypical IDP α-synuclein in solution for 2.5 μs. Four combinations of protein and water force fields were tested: ff19SB/OPC, ff19SB/TIP4P-D, ff03CMAP/TIP4P-D, and a99SB-disp/TIP4P-disp, with four independent repeat simulations for each combination. We compare our simulations to the results of a 73 μs simulation using the a99SB-disp/TIP4P-disp combination, provided by D. E. Shaw Research. From the trajectories, we predict a range of experimental observations of α-synuclein and compare them to literature data. This includes protein radius of gyration and hydration, intramolecular distances, NMR chemical shifts, and ³ J-couplings. Both ff19SB/TIP4P-D and a99SB-disp/TIP4P-disp produce extended conformational ensembles of α-synuclein that agree well with experimental radius of gyration and intramolecular distances while a99SB-disp/TIP4P-disp reproduces a balanced α-synuclein secondary structure content. It was found that ff19SB/OPC and ff03CMAP/TIP4P-D produce overly compact conformational ensembles and show discrepancies in the secondary structure content compared to the experimental data.

Targeting the Major Groove of the Palindromic d(GGCGCC)2 Sequence by Oligopeptide Derivatives of Anthraquinone Intercalators

GC-rich sequences are recurring motifs in oncogenes and retroviruses and could be targeted by noncovalent major-groove therapeutic ligands. We considered the palindromic sequence d(G₁G₂C₃G₄C₅C₆)₂, and designed several oligopeptide derivatives of the anticancer intercalator mitoxantrone. The stability of their complexes with an 18-mer oligonucleotide encompassing this sequence in its center was validated using polarizable molecular dynamics. We report the most salient structural features of two novel compounds, having a dialkylammonium group as a side chain on both arms. The anthraquinone ring is intercalated in the central d(CpG)₂ sequence with its long axis perpendicular to that of the two base pairs. On each strand, this enables each ammonium group to bind in-register to O₆/N₇ of the two facing G bases upstream. We subsequently designed tris-intercalating derivatives, each dialkylammonium substituted with a connector to an N₉-aminoacridine intercalator extending our target range from a six- to a ten-base-pair palindromic sequence, d(C₁G₂G₃G₄C₅G₆C₇C₈C₉G₁₀)₂. The structural features of the complex of the most promising derivative are reported. The present design strategy paves the way for designing intercalator-oligopeptide derivatives with even higher selectivity, targeting an increased number of DNA bases, going beyond ten.

Improving Condensed-Phase Water Dynamics with Explicit Nuclear Quantum Effects: The Polarizable Q-AMOEBA Force Field

We introduce a new parametrization of the AMOEBA polarizable force field for water denoted Q-AMOEBA, for use in simulations that explicitly account for nuclear quantum effects (NQEs). This study is made possible thanks to the recently introduced adaptive Quantum Thermal Bath (adQTB) simulation technique which computational cost is comparable to classical molecular dynamics. The flexible Q-AMOEBA model conserves the initial AMOEBA functional form, with an intermolecular potential including an atomic multipole description of electrostatic interactions (up to quadrupole), a polarization contribution based on the Thole interaction model and a buffered 14-7 potential to model van der Waals interactions. It has been obtained by using a ForceBalance fitting strategy including high-level quantum chemistry reference energies and selected condensed-phase properties targets. The final Q-AMOEBA model is shown to accurately reproduce both gas-phase and condensed-phase properties, notably improving the original AMOEBA water model. This development allows the fine study of NQEs on water liquid phase properties such as the average H-O-H angle compared to its gas-phase equilibrium value, isotope effects, and so on. Q-AMOEBA also provides improved infrared spectroscopy prediction capabilities compared to AMOEBA03. Overall, we show that the impact of NQEs depends on the underlying model functional form and on the associated strength of hydrogen bonds. Since adQTB simulations can be performed at near classical computational cost using the Tinker-HP package, Q-AMOEBA can be extended to organic molecules, proteins, and nucleic acids opening the possibility for the large-scale study of the importance of NQEs in biophysics.

Flexible multipole moments in smooth particle mesh Ewald

The smooth particle mesh Ewald sum is extended with additional force terms that arise from the so-called flexible multipole moments. These are multipole moments (of any rank) that depend explicitly on atomic positions in some local environment that can be made arbitrarily large. By introducing explicit dependence on atomic positions, flexible multipole moments are polarized by their local environment, allowing both intramolecular and intermolecular polarizations to be captured. Multipolar torques are discussed in detail, and it is shown that they arise naturally in the presented framework. Furthermore, we give details of how we validated our implementation of the flexible smooth particle mesh Ewald sum by considering two mathematical limits of the smooth particle mesh Ewald summation.

Development of the Quantum-Inspired SIBFA Many-Body Polarizable Force Field: Enabling Condensed-Phase Molecular Dynamics Simulations

We present the extension of the Sum of Interactions Between Fragments Ab initio Computed (SIBFA) many-body polarizable force field to condensed-phase molecular dynamics (MD) simulations. The quantum-inspired SIBFA procedure is grounded on simplified integrals obtained from localized molecular orbital theory and achieves full separability of its intermolecular potential. It embodies long-range multipolar electrostatics (up to quadrupole) coupled to a short-range penetration correction (up to charge-quadrupole), exchange repulsion, many-body polarization, many-body charge transfer/delocalization, exchange dispersion, and dispersion (up to C10). This enables the reproduction of all energy contributions of ab initio symmetry-adapted perturbation theory (SAPT(DFT)) gas-phase reference computations. The SIBFA approach has been integrated within the Tinker-HP massively parallel MD package. To do so, all SIBFA energy gradients have been derived and the approach has been extended to enable periodic boundary conditions simulations using smooth particle mesh Ewald. This novel implementation also notably includes a computationally tractable simplification of the many-body charge transfer/delocalization contribution. As a proof of concept, we perform a first computational experiment defining a water model fitted on a limited set of SAPT(DFT) data. SIBFA is shown to enable a satisfactory reproduction of both gas-phase energetic contributions and condensed-phase properties highlighting the importance of its physically motivated functional form.

Computationally driven discovery of SARS-CoV-2 Mpro inhibitors: from design to experimental validation

We report a fast-track computationally driven discovery of new SARS-CoV-2 main protease (M^pro) inhibitors whose potency ranges from mM for the initial non-covalent ligands to sub-μM for the final covalent compound (IC₅₀ = 830 ± 50 nM). The project extensively relied on high-resolution all-atom molecular dynamics simulations and absolute binding free energy calculations performed using the polarizable AMOEBA force field. The study is complemented by extensive adaptive sampling simulations that are used to rationalize the different ligand binding poses through the explicit reconstruction of the ligand–protein conformation space. Machine learning predictions are also performed to predict selected compound properties. While simulations extensively use high performance computing to strongly reduce the time-to-solution, they were systematically coupled to nuclear magnetic resonance experiments to drive synthesis and for in vitro characterization of compounds. Such a study highlights the power of in silico strategies that rely on structure-based approaches for drug design and allows the protein conformational multiplicity problem to be addressed. The proposed fluorinated tetrahydroquinolines open routes for further optimization of M^pro inhibitors towards low nM affinities.

The dominant binding mode of the QUB-00006-Int-07 main protease inhibitor during absolute binding free energy simulations.

Electronic Polarization Is Essential for the Stabilization and Dynamics of Buried Ion Pairs in Staphylococcal Nuclease Mutants

Buried charged residues play important roles in the modulation of protein stabilities and conformational dynamics and make crucial contributions to protein functions. Considering the generally nonpolar nature of protein interior, a key question concerns the contribution of electronic polarization to the stabilization and properties of buried charges. We answer this question by conducting free energy simulations using the latest polarizable CHARMM force field based on Drude oscillators for a series of Staphylococcal nuclease mutants that involve a buried Glu-Lys pair in different titration states and orientations. While a nonpolarizable model suggests that the ionized form of the buried Glu-Lys pair is more than 40 kcal/mol less stable than the charge-neutral form, the two titration states are comparable in stability when electronic polarization is included explicitly, a result better reconcilable with available experimental data. Analysis of free energy components suggests that additional stabilization of the ionized Glu-Lys pair has contributions from both the enhanced salt-bridge strength and stronger interaction between the ion-pair and surrounding protein residues and penetrated water. Despite the stronger direct interaction between Glu and Lys, the ion-pair exhibits considerably larger and faster structural fluctuations when polarization is included, due to compensation of interactions in the cavity. Collectively, observations from this work provide compelling evidence that electronic polarization is essential to the stability, hydration, dynamics, and therefore function of buried charges in proteins. Therefore, our study advocates for the explicit consideration of electronic polarization for mechanistic and engineering studies that implicate buried charged residues, such as enzymes and ion transporters.

O(N) Stochastic Evaluation of Many-Body van der Waals Energies in Large Complex Systems

We propose a new strategy to solve the key equations of the many-body dispersion (MBD) model by Tkatchenko, DiStasio Jr., and Ambrosetti. Our approach overcomes the original O(N3) computational complexity that limits its applicability to large molecular systems within the context of O(N) density functional theory. First, to generate the required frequency-dependent screened polarizabilities, we introduce an efficient solution to the Dyson-like self-consistent screening equations. The scheme reduces the number of variables and, coupled to a direct inversion of the iterative subspace extrapolation, exhibits linear-scaling performances. Second, we apply a stochastic Lanczos trace estimator resolution to the equations evaluating the many-body interaction energy of coupled quantum harmonic oscillators. While scaling linearly, it also enables communication-free pleasingly parallel implementations. As the resulting O(N) stochastic massively parallel MBD approach is found to exhibit minimal memory requirements, it opens up the possibility of computing accurate many-body van der Waals interactions of millions-atoms' complex materials and solvated biosystems with computational times in the range of minutes.

An Efficient Gaussian-Accelerated Molecular Dynamics (GaMD) Multilevel Enhanced Sampling Strategy: Application to Polarizable Force Fields Simulations of Large Biological Systems

We introduce a novel multilevel enhanced sampling strategy grounded on Gaussian-accelerated Molecular Dynamics (GaMD). First, we propose a GaMD multi-GPUs-accelerated implementation within the Tinker-HP molecular dynamics package. We introduce the new "dual-water" mode and its use with the flexible AMOEBA polarizable force field. By adding harmonic boosts to the water stretching and bonding terms, it accelerates the solvent-solute interactions while enabling speedups, thanks to the use of fast multiple-time step integrators. To further reduce the time-to-solution, we couple GaMD to Umbrella Sampling (US). The GaMD─US/dual-water approach is tested on the 1D Potential of Mean Force (PMF) of the solvated CD2-CD58 system (168 000 atoms), allowing the AMOEBA PMF to converge within 1 kcal/mol of the experimental value. Finally, Adaptive Sampling (AS) is added, enabling AS-GaMD capabilities but also the introduction of the new Adaptive Sampling-US-GaMD (ASUS-GaMD) scheme. The highly parallel ASUS-GaMD setup decreases time to convergence by, respectively, 10 and 20 times, compared to GaMD-US and US. Overall, beside the acceleration of PMF computations, Tinker-HP now allows for the simultaneous use of Adaptive Sampling and GaMD-"dual water" enhanced sampling approaches increasing the applicability of polarizable force fields to large-scale simulations of biological systems.

Operando Characterization of Organic Mixed Ionic/Electronic Conducting Materials

Operando characterization plays an important role in revealing the structure-property relationships of organic mixed ionic/electronic conductors (OMIECs), enabling the direct observation of dynamic changes during device operation and thus guiding the development of new materials. This review focuses on the application of different operando characterization techniques in the study of OMIECs, highlighting the time-dependent and bias-dependent structure, composition, and morphology information extracted from these techniques. We first illustrate the needs, requirements, and challenges of operando characterization then provide an overview of relevant experimental techniques, including spectroscopy, scattering, microbalance, microprobe, and electron microscopy. We also compare different in silico methods and discuss the interplay of these computational methods with experimental techniques. Finally, we provide an outlook on the future development of operando for OMIEC-based devices and look toward multimodal operando techniques for more comprehensive and accurate description of OMIECs.

Water Nanoconfined in a Hydrophobic Pore: Molecular Dynamics Simulations of Transmembrane Protein 175 and the Influence of Water Models

Water molecules within biological ion channels are in a nanoconfined environment and therefore exhibit behaviors which differ from that of bulk water. Here, we investigate the phenomenon of hydrophobic gating, the process by which a nanopore may spontaneously dewet to form a "vapor lock" if the pore is sufficiently hydrophobic and/or narrow. This occurs without steric occlusion of the pore. Using molecular dynamics simulations with both rigid fixed-charge and polarizable (AMOEBA) force fields, we investigate this wetting/dewetting behavior in the transmembrane protein 175 ion channel. We examine how a range of rigid fixed-charge and polarizable water models affect wetting/dewetting in both the wild-type structure and in mutants chosen to cover a range of nanopore radii and pore-lining hydrophobicities. Crucially, we find that the rigid fixed-charge water models lead to similar wetting/dewetting behaviors, but that the polarizable water model resulted in an increased wettability of the hydrophobic gating region of the pore. This has significant implications for molecular simulations of nanoconfined water, as it implies that polarizability may need to be included if we are to gain detailed mechanistic insights into wetting/dewetting processes. These findings are of importance for the design of functionalized biomimetic nanopores (e.g., sensing or desalination) as well as for furthering our understanding of the mechanistic processes underlying biological ion channel function.

Nuclear Quantum Effects in Liquid Water at Near Classical Computational Cost Using the Adaptive Quantum Thermal Bath

We demonstrate the accuracy and efficiency of a recently introduced approach to account for nuclear quantum effects (NQEs) in molecular simulations: the adaptive quantum thermal bath (adQTB). In this method, zero-point energy is introduced through a generalized Langevin thermostat designed to precisely enforce the quantum fluctuation-dissipation theorem. We propose a refined adQTB algorithm with improved accuracy and report adQTB simulations of liquid water. Through extensive comparison with reference path integral calculations, we demonstrate that it provides excellent accuracy for a broad range of structural and thermodynamic observables as well as infrared vibrational spectra. The adQTB has a computational cost comparable to that of classical molecular dynamics, enabling simulations of up to millions of degrees of freedom.

AENET-LAMMPS and AENET-TINKER: Interfaces for accurate and efficient molecular dynamics simulations with machine learning potentials

Machine-learning potentials (MLPs) trained on data from quantum-mechanics based first-principles methods can approach the accuracy of the reference method at a fraction of the computational cost. To facilitate efficient MLP-based molecular dynamics and Monte Carlo simulations, an integration of the MLPs with sampling software is needed. Here, we develop two interfaces that link the atomic energy network (ænet) MLP package with the popular sampling packages TINKER and LAMMPS. The three packages, ænet, TINKER, and LAMMPS, are free and open-source software that enable, in combination, accurate simulations of large and complex systems with low computational cost that scales linearly with the number of atoms. Scaling tests show that the parallel efficiency of the ænet-TINKER interface is nearly optimal but is limited to shared-memory systems. The ænet-LAMMPS interface achieves excellent parallel efficiency on highly parallel distributed-memory systems and benefits from the highly optimized neighbor list implemented in LAMMPS. We demonstrate the utility of the two MLP interfaces for two relevant example applications: the investigation of diffusion phenomena in liquid water and the equilibration of nanostructured amorphous battery materials.

Harnessing the Power of Multi-GPU Acceleration into the Quantum Interaction Computational Kernel Program

We report a new multi-GPU capable ab initio Hartree-Fock/density functional theory implementation integrated into the open source QUantum Interaction Computational Kernel (QUICK) program. Details on the load balancing algorithms for electron repulsion integrals and exchange correlation quadrature across multiple GPUs are described. Benchmarking studies carried out on up to four GPU nodes, each containing four NVIDIA V100-SXM2 type GPUs demonstrate that our implementation is capable of achieving excellent load balancing and high parallel efficiency. For representative medium to large size protein/organic molecular systems, the observed parallel efficiencies remained above 82% for the Kohn-Sham matrix formation and above 90% for nuclear gradient calculations. The accelerations on NVIDIA A100, P100, and K80 platforms also have realized parallel efficiencies higher than 68% in all tested cases, paving the way for large-scale ab initio electronic structure calculations with QUICK.

Interfacial Water Many-Body Effects Drive Structural Dynamics and Allosteric Interactions in SARS-CoV-2 Main Protease Dimerization Interface

Following our previous work ( Chem. Sci. 2021, 12, 4889-4907), we study the structural dynamics of the SARS-CoV-2 Main Protease dimerization interface (apo dimer) by means of microsecond adaptive sampling molecular dynamics simulations (50 μs) using the AMOEBA polarizable force field (PFF). This interface is structured by a complex H-bond network that is stable only at physiological pH. Structural correlations analysis between its residues and the catalytic site confirms the presence of a buried allosteric site. However, noticeable differences in allosteric connectivity are observed between PFFs and non-PFFs. Interfacial polarizable water molecules are shown to appear at the heart of this discrepancy because they are connected to the global interface H-bond network and able to adapt their dipole moment (and dynamics) to their diverse local physicochemical microenvironments. The water-interface many-body interactions appear to drive the interface volume fluctuations and to therefore mediate the allosteric interactions with the catalytic cavity.

Atomistic Polarizable Embeddings: Energy, Dynamics, Spectroscopy, and Reactivity

The computational modeling of realistic extended systems, relevant in, e.g., Chemistry and Biophysics, is a fundamental problem of paramount importance in contemporary research. Enzymatic catalysis and photoinduced processes in pigment-protein complexes are typical problems targeted by computer-aided approaches, to complement experiments as interpretative tools at a molecular scale. The daunting complexity of this task lies in between the opposite stringent requirements of results' reliability for structural/dynamical properties and related intermolecular interactions, and a mandatory principle of realism in the modeling strategy. Therefore, in practice, a truly realistic computational model of a biologically relevant system can easily fail to meet the accuracy requirement, in order to balance the excessive computational cost necessary to reach the desired precision.To address such an "accuracy vs reality" dualistic requirement, mixed quantum mechanics/classical mechanics approaches within Atomistic (i.e., preserving the discrete particle configuration) Polarizable Embeddings (QM/APEs) methods have been proposed over the years. In this Account, we review recent developments in the design and application of general QM/APE methods, targeting situations where a local intrinsically quantum behavior is coupled to a large molecular system (i.e., an environment), often involving processes with different dynamical time scales, in order to avoid brute-force, unpractical quantum chemistry calculations on the complete system.In the first place, our interest is devoted to the available APEs models presently implemented in computational software, highlighting the quantum chemistry methods that can be used to treat the QM subsystem. We review the coupling strategy between the QM subsystem and the APE, which requires to examine the way the QM/MM mutual interactions are accounted for and how the polarization of the classical environment is considered with respect to (wrt) the quantum variables. Because of the need of reliable molecular and macromolecular structures, a pivotal aspect to address here is the handling of the system dynamics (i.e., gradients wrt nuclear positions are required), especially for large molecular assemblies composed by an overwhelming number of atoms, exploring many conformations on a complex energy landscape.Alongside, we highlight our views on the necessary steps to take toward more accurate general-purposes and transferable explicit embeddings. The main objective to achieve here is to design a more physically grounded multiscale approach. To do so, one should apply advanced new generation classical models to account for refined induction effects that are able to (i) improve the quality of QM/MM interaction energies; (ii) enhance transferability by avoiding the compulsory partial (or total) reparameterization of the classical model. Moreover, the extension of recent developments originating from the field of advanced classical molecular dynamics (MD) to the realm of QM/APE methods is a key direction to improve both speed and efficiency for the phase space exploration of systems of growing size and complexity.Lastly, we point out specific research topics where an advanced QM/APE dynamics can certainly shed some light. For example, we discuss chemical reactions in "harsh" environments and the case of spectroscopic theoretical modeling where the inclusion of refined environment effects is often mandatory.