Seamless integration of GEM, a density based-force field, for QM/MM simulations via LICHEM, Psi4, and Tinker-HP.

Hybrid quantum mechanics/molecular mechanics (QM/MM) simulations have become an essential tool in computational chemistry, particularly for analyzing complex biological and condensed phase systems. Building on this foundation, our work presents a novel implementation of the Gaussian Electrostatic Model (GEM), a polarizable density-based force field, within the QM/MM framework. This advancement provides seamless integration, enabling efficient and optimized QM/GEM calculations in a single step using the LICHEM Code. We have successfully applied our implementation to water dimers and hexamers, demonstrating the ability to handle water systems with varying numbers of water molecules. Moreover, we have extended the application to describe the double proton transfer of the aspartic acid dimer in a box of water, which highlights the method's proficiency in investigating heterogeneous systems. Our implementation offers the flexibility to perform on-the-fly density fitting or to utilize pre-fitted coefficients to estimate exchange and Coulomb contributions. This flexibility enhances efficiency and accuracy in modeling molecular interactions, especially in systems where polarization effects are significant.

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.

The Extended Eighth-Shell method for periodic boundary conditions with rotational symmetry.

The Eighth-Shell method for parallelization of molecular dynamics simulations has previously been shown to be the most optimal for efficiency at large process counts. However, in its current formulation only the P1 space group is supported for periodic boundary conditions (PBC) and thus reflection and/or rotational crystal symmetries are not supported. In this work, we outline the development and implementation of the Extended Eighth-Shell (EES) method that allows rotational symmetry by using an extended import region compared to the ES method. It simulates only the asymmetric unit and communicates coordinates and forces with images that correspond to P21 PBC. The P21 PBC has application in lipid bilayer simulations as it can be used to allow lipids to switch leaftlets, thus rapidly balancing the chemical potential difference between the two layers. Our results show that the EES method scales efficiently over large number of processes and can be used for simulations with P21 symmetry in an orthorhombic crystal.

Analytical Hessians for Ewald and particle mesh Ewald electrostatics.

The particle mesh Ewald (PME) method has become ubiquitous in the molecular simulation community due to its ability to deliver long range electrostatics accurately with ON ⁡log(N) complexity. Despite this widespread use, spanning more than two decades, second derivatives (Hessians) have not been available. In this work, we describe the theory and implementation of PME Hessians, which have applications in normal mode analysis, characterization of stationary points, phonon dispersion curve calculation, crystal structure prediction, and efficient geometry optimization. We outline an exact strategy that requires O(1) effort for each Hessian element; after discussing the excessive memory requirements of such an approach, we develop an accurate, efficient approximation that is far more tractable on commodity hardware.

A compression strategy for particle mesh Ewald theory.

Particle Mesh Ewald (PME) has become a standard method for treating long-range electrostatics in molecular simulations. Although the method has inferior asymptotic computational complexity to its linear scaling competitors, it remains enormously popular due to its high efficiency, which stems from the use of fast Fourier transforms (FFTs). This use of FFTs provides great challenges for scaling the method up to massively parallel systems, in large part because of the need to transfer large amounts of data. In this work, we demonstrate that this data transfer volume can be greatly reduced as a natural consequence of the structure of the PME equations. We also suggest an alternative algorithm that supplants the FFT with a linear algebra approach, which further decreases communication costs at the expense of increased asymptotic computational complexity. This linear algebra based approach is demonstrated to have great potential for latency hiding by interleaving communication and computation steps of the short- and long-range electrostatic terms.

A simplified charge projection scheme for long-range electrostatics in ab initio QM/MM calculations.

In a previous work [Pan et al., Molecules 23, 2500 (2018)], a charge projection scheme was reported, where outer molecular mechanical (MM) charges [>10 Å from the quantum mechanical (QM) region] were projected onto the electrostatic potential (ESP) grid of the QM region to accurately and efficiently capture long-range electrostatics in ab initio QM/MM calculations. Here, a further simplification to the model is proposed, where the outer MM charges are projected onto inner MM atom positions (instead of ESP grid positions). This enables a representation of the long-range MM electrostatic potential via augmentary charges (AC) on inner MM atoms. Combined with the long-range electrostatic correction function from Cisneros et al. [J. Chem. Phys. 143, 044103 (2015)] to smoothly switch between inner and outer MM regions, this new QM/MM-AC electrostatic model yields accurate and continuous ab initio QM/MM electrostatic energies with a 10 Å cutoff between inner and outer MM regions. This model enables efficient QM/MM cluster calculations with a large number of MM atoms as well as QM/MM calculations with periodic boundary conditions.

Quantum Chemistry Common Driver and Databases (QCDB) and Quantum Chemistry Engine (QCEngine): Automation and interoperability among computational chemistry programs.

Community efforts in the computational molecular sciences (CMS) are evolving toward modular, open, and interoperable interfaces that work with existing community codes to provide more functionality and composability than could be achieved with a single program. The Quantum Chemistry Common Driver and Databases (QCDB) project provides such capability through an application programming interface (API) that facilitates interoperability across multiple quantum chemistry software packages. In tandem with the Molecular Sciences Software Institute and their Quantum Chemistry Archive ecosystem, the unique functionalities of several CMS programs are integrated, including CFOUR, GAMESS, NWChem, OpenMM, Psi4, Qcore, TeraChem, and Turbomole, to provide common computational functions, i.e., energy, gradient, and Hessian computations as well as molecular properties such as atomic charges and vibrational frequency analysis. Both standard users and power users benefit from adopting these APIs as they lower the language barrier of input styles and enable a standard layout of variables and data. These designs allow end-to-end interoperable programming of complex computations and provide best practices options by default.

Psi4 1.4: Open-source software for high-throughput quantum chemistry.

PSI4 is a free and open-source ab initio electronic structure program providing implementations of Hartree-Fock, density functional theory, many-body perturbation theory, configuration interaction, density cumulant theory, symmetry-adapted perturbation theory, and coupled-cluster theory. Most of the methods are quite efficient, thanks to density fitting and multi-core parallelism. The program is a hybrid of C++ and Python, and calculations may be run with very simple text files or using the Python API, facilitating post-processing and complex workflows; method developers also have access to most of PSI4's core functionalities via Python. Job specification may be passed using The Molecular Sciences Software Institute (MolSSI) QCSCHEMA data format, facilitating interoperability. A rewrite of our top-level computation driver, and concomitant adoption of the MolSSI QCARCHIVE INFRASTRUCTURE project, makes the latest version of PSI4 well suited to distributed computation of large numbers of independent tasks. The project has fostered the development of independent software components that may be reused in other quantum chemistry programs.

OpenMM 7: Rapid development of high performance algorithms for molecular dynamics.

OpenMM is a molecular dynamics simulation toolkit with a unique focus on extensibility. It allows users to easily add new features, including forces with novel functional forms, new integration algorithms, and new simulation protocols. Those features automatically work on all supported hardware types (including both CPUs and GPUs) and perform well on all of them. In many cases they require minimal coding, just a mathematical description of the desired function. They also require no modification to OpenMM itself and can be distributed independently of OpenMM. This makes it an ideal tool for researchers developing new simulation methods, and also allows those new methods to be immediately available to the larger community.

Mapping the Drude polarizable force field onto a multipole and induced dipole model.

The induced dipole and the classical Drude oscillator represent two major approaches for the explicit inclusion of electronic polarizability into force field-based molecular modeling and simulations. In this work, we explore the equivalency of these two models by comparing condensed phase properties computed using the Drude force field and a multipole and induced dipole (MPID) model. Presented is an approach to map the electrostatic model optimized in the context of the Drude force field onto the MPID model. Condensed phase simulations on water and 15 small model compounds show that without any reparametrization, the MPID model yields properties similar to the Drude force field with both models yielding satisfactory reproduction of a range of experimental values and quantum mechanical data. Our results illustrate that the Drude oscillator model and the point induced dipole model are different representations of essentially the same physical model. However, results indicate the presence of small differences between the use of atomic multipoles and off-center charge sites. Additionally, results on the use of dispersion particle mesh Ewald further support its utility for treating long-range Lennard Jones dispersion contributions in the context of polarizable force fields. The main motivation in demonstrating the transferability of parameters between the Drude and MPID models is that the more than 15 years of development of the Drude polarizable force field can now be used with MPID formalism without the need for dual-thermostat integrators nor self-consistent iterations. This opens up a wide range of new methodological opportunities for polarizable models.

Blind prediction of distribution in the SAMPL5 challenge with QM based protomer and pK a corrections.

The computation of distribution coefficients between polar and apolar phases requires both an accurate characterization of transfer free energies between phases and proper accounting of ionization and protomerization. We present a protocol for accurately predicting partition coefficients between two immiscible phases, and then apply it to 53 drug-like molecules in the SAMPL5 blind prediction challenge. Our results combine implicit solvent QM calculations with classical MD simulations using the non-Boltzmann Bennett free energy estimator. The OLYP/DZP/SMD method yields predictions that have a small deviation from experiment (RMSD = 2.3 [Formula: see text] D units), relative to other participants in the challenge. Our free energy corrections based on QM protomer and [Formula: see text] calculations increase the correlation between predicted and experimental distribution coefficients, for all methods used. Unfortunately, these corrections are overly hydrophilic, and fail to account for additional effects such as aggregation, water dragging and the presence of polar impurities in the apolar phase. We show that, although expensive, QM-NBB free energy calculations offer an accurate and robust method that is superior to standard MM and QM techniques alone.

Calculating distribution coefficients based on multi-scale free energy simulations: an evaluation of MM and QM/MM explicit solvent simulations of water-cyclohexane transfer in the SAMPL5 challenge.

One of the central aspects of biomolecular recognition is the hydrophobic effect, which is experimentally evaluated by measuring the distribution coefficients of compounds between polar and apolar phases. We use our predictions of the distribution coefficients between water and cyclohexane from the SAMPL5 challenge to estimate the hydrophobicity of different explicit solvent simulation techniques. Based on molecular dynamics trajectories with the CHARMM General Force Field, we compare pure molecular mechanics (MM) with quantum-mechanical (QM) calculations based on QM/MM schemes that treat the solvent at the MM level. We perform QM/MM with both density functional theory (BLYP) and semi-empirical methods (OM1, OM2, OM3, PM3). The calculations also serve to test the sensitivity of partition coefficients to solute polarizability as well as the interplay of the quantum-mechanical region with the fixed-charge molecular mechanics environment. Our results indicate that QM/MM with both BLYP and OM2 outperforms pure MM. However, this observation is limited to a subset of cases where convergence of the free energy can be achieved.

An efficient protocol for obtaining accurate hydration free energies using quantum chemistry and reweighting from molecular dynamics simulations.

The non-Boltzmann Bennett (NBB) free energy estimator method is applied to 21 molecules from the blind subset of the SAMPL4 challenge. When NBB is applied with the SMD implicit solvent model, and the OLYP/DZP level of quantum chemistry, highly accurate hydration free energy calculations are obtained with respect to experiment (RMSD=0.89kcal·mol-1). Other quantum chemical methods are also tested, and the effects of solvent model, density functional, basis set are explored in this benchmarking study, providing a framework for improvements in calculating hydration free energies. We provide a practical guide for using the best QM-NBB protocols that are consistently more accurate than either pure QM or pure MM alone. In situations where high accuracy hydration free energy predictions are needed, the QM-NBB method with SMD implicit solvent should be the first choice of quantum chemists.

An empirical extrapolation scheme for efficient treatment of induced dipoles.

Many cutting edge force fields include polarization, to enhance their accuracy and range of applicability. In this work, we develop efficient strategies for the induced dipole polarization method. By fitting various orders of perturbation theory (PT) dipoles to a diverse training set, we arrive at a family of fully analytic methods - whose nth order is referred to OPTn - that span the full spectrum of polarization methods from the fast zeroth-order approach that neglects mutual dipole coupling, approaching the fully variational approach at high order. Our training set contains many difficult cases where the PT series diverges, and we demonstrate that our OPTn methods still deliver excellent results in these cases. Our tests show that the OPTn methods exhibit rapid convergence towards the exact answer with each increasing PT order. The fourth order OPT4 method, whose costs are commensurate with three iterations of the leading conjugate gradient method, is a particularly promising candidate to be used as a drop-in replacement for existing solvers without further parameterization.

Numerical study on the partitioning of the molecular polarizability into fluctuating charge and induced atomic dipole contributions.

In order to carry out a detailed analysis of the molecular static polarizability, which is the response of the molecule to a uniform external electric field, the molecular polarizability was computed using the finite-difference method for 21 small molecules, using density functional theory. Within nine charge population schemes (Löwdin, Mulliken, Becke, Hirshfeld, CM5, Hirshfeld-I, NPA, CHELPG, MK-ESP) in common use, the charge fluctuation contribution is found to dominate the molecular polarizability, with its ratio ranging from 59.9% with the Hirshfeld or CM5 scheme to 96.2% with the Mulliken scheme. The Hirshfeld-I scheme is also used to compute the other contribution to the molecular polarizability coming from the induced atomic dipoles, and the atomic polarizabilities in eight small molecules and water pentamer are found to be highly anisotropic for most atoms. Overall, the results suggest that (a) more emphasis probably should be placed on the charge fluctuation terms in future polarizable force field development and (b) an anisotropic polarizability might be more suitable than an isotropic one in polarizable force fields based entirely or partially on the induced atomic dipoles.

Efficient treatment of induced dipoles.

Most existing treatments of induced dipoles in polarizable molecular mechanics force field calculations use either the self-consistent variational method, which is solved iteratively, or the "direct" approximation that is non-iterative as a result of neglecting coupling between induced dipoles. The variational method is usually implemented using assumptions that are only strictly valid under tight convergence of the induced dipoles, which can be computationally demanding to enforce. In this work, we discuss the nature of the errors that result from insufficient convergence and suggest a strategy that avoids such problems. Using perturbation theory to reintroduce the mutual coupling into the direct algorithm, we present a computationally efficient method that combines the precision of the direct approach with the accuracy of the variational approach. By analyzing the convergence of this perturbation series, we derive a simple extrapolation formula that delivers a very accurate approximation to the infinite order solution at the cost of only a few iterations. We refer to the new method as extrapolated perturbation theory. Finally, we draw connections to our previously published permanent multipole algorithm to develop an efficient implementation of the electric field and Thole terms and also derive some necessary, but not sufficient, criteria that force field parameters must obey.

Capture of hydroxymethylene and its fast disappearance through tunnelling.

Singlet carbenes exhibit a divalent carbon atom whose valence shell contains only six electrons, four involved in bonding to two other atoms and the remaining two forming a non-bonding electron pair. These features render singlet carbenes so reactive that they were long considered too short-lived for isolation and direct characterization. This view changed when it was found that attaching the divalent carbon atom to substituents that are bulky and/or able to donate electrons produces carbenes that can be isolated and stored. N-heterocyclic carbenes are such compounds now in wide use, for example as ligands in metathesis catalysis. In contrast, oxygen-donor-substituted carbenes are inherently less stable and have been less studied. The pre-eminent case is hydroxymethylene, H-C-OH; although it is the key intermediate in the high-energy chemistry of its tautomer formaldehyde, has been implicated since 1921 in the photocatalytic formation of carbohydrates, and is the parent of alkoxycarbenes that lie at the heart of transition-metal carbene chemistry, all attempts to observe this species or other alkoxycarbenes have failed. However, theoretical considerations indicate that hydroxymethylene should be isolatable. Here we report the synthesis of hydroxymethylene and its capture by matrix isolation. We unexpectedly find that H-C-OH rearranges to formaldehyde with a half-life of only 2 h at 11 K by pure hydrogen tunnelling through a large energy barrier in excess of 30 kcal mol(-1).

An efficient algorithm for multipole energies and derivatives based on spherical harmonics and extensions to particle mesh Ewald.

Next-generation molecular force fields deliver accurate descriptions of non-covalent interactions by employing more elaborate functional forms than their predecessors. Much work has been dedicated to improving the description of the electrostatic potential (ESP) generated by these force fields. A common approach to improving the ESP is by augmenting the point charges on each center with higher-order multipole moments. The resulting anisotropy greatly improves the directionality of the non-covalent bonding, with a concomitant increase in computational cost. In this work, we develop an efficient strategy for enumerating multipole interactions, by casting an efficient spherical harmonic based approach within a particle mesh Ewald (PME) framework. Although the derivation involves lengthy algebra, the final expressions are relatively compact, yielding an approach that can efficiently handle both finite and periodic systems without imposing any approximations beyond PME. Forces and torques are readily obtained, making our method well suited to modern molecular dynamics simulations.

OpenMM 8: Molecular Dynamics Simulation with Machine Learning Potentials.

Machine learning plays an important and growing role in molecular simulation. The newest version of the OpenMM molecular dynamics toolkit introduces new features to support the use of machine learning potentials. Arbitrary PyTorch models can be added to a simulation and used to compute forces and energy. A higher-level interface allows users to easily model their molecules of interest with general purpose, pretrained potential functions. A collection of optimized CUDA kernels and custom PyTorch operations greatly improves the speed of simulations. We demonstrate these features in simulations of cyclin-dependent kinase 8 (CDK8) and the green fluorescent protein chromophore in water. Taken together, these features make it practical to use machine learning to improve the accuracy of simulations with only a modest increase in cost.

New potential energy surface features for the Li + HF → LiF + H reaction.

The existing potential energy surfaces for the Li + HF system have been challenged by the experiments of Loesch, Stienkemeier, and co-workers. Here a very accurate potential energy surface has been obtained with rather rigorous theoretical methods. Methods up to full CCSDT have been pursued with basis sets as large as core correlated quintuple ζ. Reported here are the reactants, products, two transition states, and three intermediate complexes for this reaction. These reveal one previously undiscovered equilibrium geometry. The stationary point relative energies are very sensitive to level of theory. The reaction has a classical endothermicity of 2.6 kcal mol(-1). The complex Li···HF in the entrance valley lies 6.1 kcal/mol below the reactants. The expected transition state Li···H···F is bent with an angle of 72.2° and lies 4.5 kcal/mol above the reactants. The latter predicted classical barrier should be no more than one kcal/mol above the exact barrier. Not one but two product complexes lie 1.6 and 2.2 kcal/mol above reactants, respectively. Between the two product complexes, a second transition state, very broad, is found. The vibrational frequencies and zero-point vibrational energies (ZPVE) of all stationary points are reported, and significantly affect the relative energies.