Active sparse Bayesian committee machine potential for isothermal-isobaric molecular dynamics simulations.

Recent advancements in machine learning potentials (MLPs) have significantly impacted the fields of chemistry, physics, and biology by enabling large-scale first-principles simulations. Among different machine learning approaches, kernel-based MLPs distinguish themselves through their ability to handle small datasets, quantify uncertainties, and minimize over-fitting. Nevertheless, their extensive computational requirements present considerable challenges. To alleviate these, sparsification methods have been developed, aiming to reduce computational scaling without compromising accuracy. In the context of isothermal and isobaric ML molecular dynamics (MD) simulations, achieving precise pressure estimation is crucial for reproducing reliable system behavior under constant pressure. Despite progress, sparse kernel MLPs struggle with precise pressure prediction. Here, we introduce a virial kernel function that significantly enhances the pressure estimation accuracy of MLPs. Additionally, we propose the active sparse Bayesian committee machine (BCM) potential, an on-the-fly MLP architecture that aggregates local sparse Gaussian process regression (SGPR) MLPs. The sparse BCM potential overcomes the steep computational scaling with the kernel size, and a predefined restriction on the size of kernel allows for fast and efficient on-the-fly training. Our advancements facilitate accurate and computationally efficient machine learning-enhanced MD (MLMD) simulations across diverse systems, including ice-liquid coexisting phases, Li10Ge(PS6)2 lithium solid electrolyte, and high-pressure liquid boron nitride.

Learned mappings for targeted free energy perturbation between peptide conformations.

Targeted free energy perturbation uses an invertible mapping to promote configuration space overlap and the convergence of free energy estimates. However, developing suitable mappings can be challenging. Wirnsberger et al. [J. Chem. Phys. 153, 144112 (2020)] demonstrated the use of machine learning to train deep neural networks that map between Boltzmann distributions for different thermodynamic states. Here, we adapt their approach to the free energy differences of a flexible bonded molecule, deca-alanine, with harmonic biases and different spring centers. When the neural network is trained until "early stopping"-when the loss value of the test set increases-we calculate accurate free energy differences between thermodynamic states with spring centers separated by 1 Å and sometimes 2 Å. For more distant thermodynamic states, the mapping does not produce structures representative of the target state, and the method does not reproduce reference calculations.

The Back Door to the Surface Hydrated Electron.

We use a Mg+ metal to extend the size regime of aqueous clusters to extrapolate to the bulk limit of the vertical detachment energy (VDE) of the solvated electron to >3,200, a value between 1 to over 2 orders of magnitude larger than the one previously measured experimentally or computed theoretically. We relate the VDE to the energy difference between the Mg+(H2O)n and Mg2+(H2O)n systems and the metal's second ionization potential. The extrapolated bulk VDEs of the localized surface electron, which moves away from the metal as n increases, are 1.89 ± 0.01 eV for semiempirical (n ∼ 3,200; PM6-D3H4) and 1.73 ± 0.03 eV (n ∼ 150; HF) and 1.83 ± 0.02 eV (n ∼ 150; MP2) for ab initio, in excellent agreement with the 1.6-1.8 eV range of experimental results. The VDEs converge from above (larger values) to the bulk limit, in a manner that is qualitatively opposite from previous studies and experiments lacking a charged metal, a fact justifying the "back door" approach to the solvated electron.

Identification of the riboflavin cofactor-binding site in the Vibrio cholerae ion-pumping NQR complex: A novel structural motif in redox enzymes.

The ion-pumping NQR complex is an essential respiratory enzyme in the physiology of many pathogenic bacteria. This enzyme transfers electrons from NADH to ubiquinone through several cofactors, including riboflavin (vitamin B2). NQR is the only enzyme reported that is able to use riboflavin as a cofactor. Moreover, the riboflavin molecule is found as a stable neutral semiquinone radical. The otherwise highly reactive unpaired electron is stabilized via an unknown mechanism. Crystallographic data suggested that riboflavin might be found in a superficially located site in the interface of NQR subunits B and E. However, this location is highly problematic, as the site does not have the expected physiochemical properties. In this work, we have located the riboflavin-binding site in an amphipathic pocket in subunit B, previously proposed to be the entry site of sodium. Here, we show that this site contains absolutely conserved residues, including N200, N203, and D346. Mutations of these residues decrease enzymatic activity and specifically block the ability of NQR to bind riboflavin. Docking analysis and molecular dynamics simulations indicate that these residues participate directly in riboflavin binding, establishing hydrogen bonds that stabilize the cofactor in the site. We conclude that riboflavin is likely bound in the proposed pocket, which is consistent with enzymatic characterizations, thermodynamic studies, and distance between cofactors.

Electrostatics and water occlusion regulate covalently-bound flavin mononucleotide cofactors of Vibrio cholerae respiratory complex NQR.

Proteins like NADH:ubiquinone oxidoreductase (NQR), an essential enzyme and ion pump in the physiology of several pathogenic bacteria, tightly regulate the redox properties of their cofactors. Although flavin mononucleotide (FMN) is fully reduced in aqueous solution, FMN in subunits B and C of NQR exclusively undergo one-electron transitions during its catalytic cycle. Here, we perform ab initio calculations and molecular dynamics simulations to elucidate the mechanisms that regulate the redox state of FMN in NQR. QM/MM calculations show that binding site electrostatics disfavor anionic forms of FMNH2 , but permit a neutral form of the fully reduced flavin. The potential energy surface is unaffected by covalent bonding between FMN and threonine. Molecular dynamics simulations show that the FMN binding sites are inaccessible by water, suggesting that further reductions of the cofactors are limited or prohibited by the availability of water and other proton donors. These findings provide a deeper understanding of the mechanisms used by NQR to regulate electron transfer through the cofactors and perform its physiologic role. They also provide the first, to our knowledge, evidence of the simple concept that proteins regulate flavin redox states via water occlusion.

On the polarization of ligands by proteins.

Although ligand-binding sites in many proteins contain a high number density of charged side chains that can polarize small organic molecules and influence binding, the magnitude of this effect has not been studied in many systems. Here, we use a quantum mechanics/molecular mechanics (QM/MM) approach, in which the ligand is the QM region, to compute the ligand polarization energy of 286 protein-ligand complexes from the PDBBind Core Set (release 2016). Calculations were performed both with and without implicit solvent based on the domain decomposition Conductor-like Screening Model. We observe that the ligand polarization energy is linearly correlated with the magnitude of the electric field acting on the ligand, the magnitude of the induced dipole moment, and the classical polarization energy. The influence of protein and cation charges on the ligand polarization diminishes with the distance and is below 2 kcal mol-1 at 9 Å and 1 kcal mol-1 at 12 Å. Compared to these embedding field charges, implicit solvent has a relatively minor effect on ligand polarization. Considering both polarization and solvation appears essential to computing negative binding energies in some crystallographic complexes. Solvation, but not polarization, is essential for achieving moderate correlation with experimental binding free energies.

Molecular-Level Insight of the Effect of Hofmeister Anions on the Interfacial Surface Tension of a Model Protein.

The effect of the Hofmeister anions on the precipitation of proteins is often discussed using liquid-vapor coexisting systems with the assumption that the liquid-vapor interface mimics the liquid-protein interface. Solvated proteins, however, have both hydrophobic and hydrophilic regions on their surfaces rather than just a pure hydrophobic one. Using a solvated parallel ß-sheet layer consisting of both hydrophobic and positively charged hydrophilic surfaces, we investigated the adsorption of kosmotropic (SO42-) and chaotropic (ClO4-) anions toward the protein's hydrophobic and hydrophilic surfaces via Born-Oppenheimer molecular dynamics simulations using the BLYP density functional theory. It was found that both anions prefer to reside on the hydrophilic surface. Furthermore, kosmotropic anions, like SO42-, enhance the interfacial surface tension of the protein and stabilize the protein, whereas, in contrast, chaotropic anions, like ClO4-, weaken the interfacial surface tension of the protein and allow water molecules to penetrate toward the peptide bonds to form water-peptide hydrogen bonds, thus destabilizing the protein.

Ice Ih anomalies: Thermal contraction, anomalous volume isotope effect, and pressure-induced amorphization.

Ice Ih displays several anomalous thermodynamic properties such as thermal contraction at low temperatures, an anomalous volume isotope effect (VIE) rendering the volume of D2O ice greater than that of H2O ice, and a pressure-induced transition to the high-density amorphous (HDA) phase. Furthermore, the anomalous VIE increases with temperature, despite its quantum-mechanical origin. Here, embedded-fragment ab initio second-order many-body perturbation (MP2) theory in the quasiharmonic approximation (QHA) is applied to the Gibbs energy of an infinite, proton-disordered crystal of ice Ih at wide ranges of temperatures and pressures. The quantum effect of nuclei moving in anharmonic potentials is taken into account from first principles without any empirical or nonsystematic approximation to either the electronic or vibrational Hamiltonian. MP2 predicts quantitatively correctly the thermal contraction at low temperatures, which is confirmed to originate from the volume-contracting hydrogen-bond bending modes (acoustic phonons). It qualitatively reproduces (but underestimates) the thermal expansion at higher temperatures, caused by the volume-expanding hydrogen-bond stretching (and to a lesser extent librational) modes. The anomalous VIE is found to be the result of subtle cancellations among closely competing isotope effects on volume from all modes. Consequently, even ab initio MP2 with the aug-cc-pVDZ and aug-cc-pVTZ basis sets has difficulty reproducing this anomaly, yielding qualitatively varied predictions of the sign of the VIE depending on such computational details as the choice of the embedding field. However, the temperature growth of the anomalous VIE is reproduced robustly and is ascribed to the librational modes. These solid-state MP2 calculations, as well as MP2 Born-Oppenheimer molecular dynamics, find a volume collapse and a loss of symmetry and long-range order in ice Ih upon pressure loading of 2.35 GPa or higher. Concomitantly, rapid softening of acoustic phonons is observed starting around 2 GPa. They constitute a computational detection of a mechanical instability in ice Ih and the resulting pressure-induced amorphization to HDA.

Why Is MP2-Water "Cooler" and "Denser" than DFT-Water?

Density functional theory (DFT) with a dispersionless generalized gradient approximation (GGA) needs much higher temperature and pressure than the ambient conditions to maintain water in the liquid phase at the correct (1 g/cm(3)) density during first-principles simulations. Conversely, ab initio second-order many-body perturbation (MP2) calculations of liquid water require lower temperature and pressure than DFT/GGA to keep water liquid. Here we present a unifying explanation of these trends derived from classical water simulations using a polarizable force field with different sets of parameters. We show that the different temperatures and pressures between DFT/GGA and MP2 at which the simulated water displays the experimentally observed liquid structure under the ambient conditions can be largely explained by their differences in polarizability and dispersion interaction, respectively. In DFT/GGA, the polarizability and thus the induced dipole moments and the hydrogen-bond strength are all overestimated. This hinders the rotational motion of molecules and requires a higher temperature for DFT-water to be liquid. MP2 gives a stronger dispersion interaction and thus shorter intermolecular distances than dispersionless DFT/GGA, which is why MP2-water is denser than DFT-water under the same external pressure.

Ab initio molecular dynamics of liquid water using embedded-fragment second-order many-body perturbation theory towards its accurate property prediction.

A direct, simultaneous calculation of properties of a liquid using an ab initio electron-correlated theory has long been unthinkable. Here we present structural, dynamical, and response properties of liquid water calculated by ab initio molecular dynamics using the embedded-fragment spin-component-scaled second-order many-body perturbation method with the aug-cc-pVDZ basis set. This level of theory is chosen as it accurately and inexpensively reproduces the water dimer potential energy surface from the coupled-cluster singles, doubles, and noniterative triples with the aug-cc-pVQZ basis set, which is nearly exact. The calculated radial distribution function, self-diffusion coefficient, coordinate number, and dipole moment, as well as the infrared and Raman spectra are in excellent agreement with experimental results. The shapes and widths of the OH stretching bands in the infrared and Raman spectra and their isotropic-anisotropic Raman noncoincidence, which reflect the diverse local hydrogen-bond environment, are also reproduced computationally. The simulation also reveals intriguing dynamic features of the environment, which are difficult to probe experimentally, such as a surprisingly large fluctuation in the coordination number and the detailed mechanism by which the hydrogen donating water molecules move across the first and second shells, thereby causing this fluctuation.

The Melting Temperature of Liquid Water with the Effective Fragment Potential.

The direct simulation of the solid-liquid water interface with the effective fragment potential (EFP) via the constant enthalpy and pressure (NPH) ensemble was used to estimate the melting temperature (T(m)) of ice-I(h). Initial configurations and velocities, taken from equilibrated constant pressure and temperature (NPT) simulations at P = 1 atm and T = 305 K, 325 K and 399 K, respectively, yielded corresponding T(m) values of 378 ± 16 K, 382 ± 14 K and 384 ± 15 K. These estimates are consistently higher than experiment, albeit to the same degree as previously reported estimates using density functional theory (DFT)-based Born-Oppenheimer simulations with the Becke-Lee-Yang-Parr functional plus dispersion corrections (BLYP-D).

Stochastic, real-space, imaginary-time evaluation of third-order Feynman-Goldstone diagrams.

A new, alternative set of interpretation rules of Feynman-Goldstone diagrams for many-body perturbation theory is proposed, which translates diagrams into algebraic expressions suitable for direct Monte Carlo integrations. A vertex of a diagram is associated with a Coulomb interaction (rather than a two-electron integral) and an edge with the trace of a Green's function in real space and imaginary time. With these, 12 diagrams of third-order many-body perturbation (MP3) theory are converted into 20-dimensional integrals, which are then evaluated by a Monte Carlo method. It uses redundant walkers for convergence acceleration and a weight function for importance sampling in conjunction with the Metropolis algorithm. The resulting Monte Carlo MP3 method has low-rank polynomial size dependence of the operation cost, a negligible memory cost, and a naturally parallel computational kernel, while reproducing the correct correlation energies of small molecules within a few mEh after 10(6) Monte Carlo steps.

Communication: stochastic evaluation of explicitly correlated second-order many-body perturbation energy.

A stochastic algorithm is proposed that can compute the basis-set-incompleteness correction to the second-order many-body perturbation (MP2) energy of a polyatomic molecule. It evaluates the sum of two-, three-, and four-electron integrals over an explicit function of electron-electron distances by a Monte Carlo (MC) integration at an operation cost per MC step increasing only quadratically with size. The method can reproduce the corrections to the MP2/cc-pVTZ energies of H2O, CH4, and C6H6 within a few mEh after several million MC steps. It circumvents the resolution-of-the-identity approximation to the nonfactorable three-electron integrals usually necessary in the conventional explicitly correlated (R12 or F12) methods.

Stochastic evaluation of second-order Dyson self-energies.

A stochastic method is proposed that evaluates the second-order perturbation corrections to the Dyson self-energies of a molecule (i.e., quasiparticle energies or correlated ionization potentials and electron affinities) directly and not as small differences between two large, noisy quantities. With the aid of a Laplace transform, the usual sum-of-integral expressions of the second-order self-energy in many-body Green's function theory are rewritten into a sum of just four 13-dimensional integrals, 12-dimensional parts of which are evaluated by Monte Carlo integration. Efficient importance sampling is achieved with the Metropolis algorithm and a 12-dimensional weight function that is analytically integrable, is positive everywhere, and cancels all the singularities in the integrands exactly and analytically. The quasiparticle energies of small molecules have been reproduced within a few mEh of the correct values with 10(8) Monte Carlo steps. Linear-to-quadratic scaling of the size dependence of computational cost is demonstrated even for these small molecules.

Molecular Insight into Different Denaturing Efficiency of Urea, Guanidinium, and Methanol: A Comparative Simulation Study.

We have designed various nanoslit systems, whose opposing surfaces can be either hydrophobic, hydrophilic, or simply a water-vapor interface, for the molecular dynamics simulation of confined water with three different protein denaturants, i.e., urea, guanidinium chloride (GdmCl), and methanol, respectively. Particular attention is placed on the preferential adsorption of the denaturant molecules onto the opposing surfaces and associated resident time in the vicinal layer next to the surfaces, as well as their implication in the denaturing efficiency of different denaturant molecules. Our simulation results show that among the three denaturants, the occupancy of methanol in the vicinal layer is the highest while the residence time of Gdm is the longest. Although the occupancy and the residence time of urea in the vicinal layer is less than those of the other two denaturant molecules, urea entails "all-around" properties for being a highly effective denaturant. The distinct characteristics of three denaturants may suggest a different molecular mechanism for the protein denaturation. This comparative simulation by design allows us to gain additional insights, on the molecular level, into the denaturation effect and related hydrophobic effect.

Convergence Acceleration of Parallel Monte Carlo Second-Order Many-Body Perturbation Calculations Using Redundant Walkers.

A Monte Carlo (MC) integration of the second-order many-body perturbation (MP2) corrections to energies and self-energies eliminates the usual computational bottleneck of the MP2 algorithm (i.e., the basis transformation of two-electron integrals), thereby achieving near-linear size dependence of its operation cost, a negligible core and disk memory cost, and a naturally parallel computational kernel. In this method, the correlation correction expressions are recast into high-dimensional integrals, which are approximated as the sums of integrands evaluated at coordinates of four electron random walkers guided by a Metropolis algorithm for importance sampling. The gravest drawback of this method, however, is the inevitable statistical uncertainties in its results, which decay slowly as the inverse square-root of the number of MC steps. We propose an algorithm that can increase the number of MC sampling points in each MC step by many orders of magnitude by having 2m electron walkers (m > 2) and then using m(m - 1)/2 permutations of their coordinates in evaluating the integrands. Hence, this algorithm brings an O(m(2))-fold increase in the number of MC sampling points at a mere O(m) additional cost of propagating redundant walkers, which is a net O(m)-fold enhancement in sampling efficiency. We have demonstrated a large performance increase in the Monte Carlo MP2 calculations for the correlation energies of benzene and benzene dimer as well as for the correlation corrections to the energy, ionization potential, and electron affinity of C60. The calculation for C60 has been performed with a parallel implementation of this method running on up to 400 processors of a supercomputer, yielding an accurate prediction of its large electron affinity, which makes its derivative useful as an electron acceptor in bulk heterojunction organic solar cells.

Stochastic evaluation of second-order many-body perturbation energies.

With the aid of the Laplace transform, the canonical expression of the second-order many-body perturbation correction to an electronic energy is converted into the sum of two 13-dimensional integrals, the 12-dimensional parts of which are evaluated by Monte Carlo integration. Weight functions are identified that are analytically normalizable, are finite and non-negative everywhere, and share the same singularities as the integrands. They thus generate appropriate distributions of four-electron walkers via the Metropolis algorithm, yielding correlation energies of small molecules within a few mE(h) of the correct values after 10(8) Monte Carlo steps. This algorithm does away with the integral transformation as the hotspot of the usual algorithms, has a far superior size dependence of cost, does not suffer from the sign problem of some quantum Monte Carlo methods, and potentially easily parallelizable and extensible to other more complex electron-correlation theories.

NH4(+) Resides Inside the Water 20-mer Cage As Opposed to H3O(+), Which Resides on the Surface: A First Principles Molecular Dynamics Simulation Study.

Experimental vibrational predissociation spectra of the magic NH4(+)(H2O)20 clusters are close to those of the magic H3O(+)(H2O)20 clusters. It has been assumed that the geometric features of NH4(+)(H2O)20 clusters might be close to those of H3O(+)(H2O)20 clusters, in which H3O(+) resides on the surface. Car-Parrinello molecular dynamics simulations in conjunction with density functional theory calculations are performed to generate the infrared spectra of the magic NH4(+)(H2O)20 clusters. In comparison with the experimental vibrational predissociation spectra of NH4(+)(H2O)20, we find that NH4(+) is inside the cage structure of NH4(+)(H2O)20 as opposed to on the surface structure. This shows a clear distinction between the structures of NH4(+)(H2O)20 and H3O(+)(H2O)20 as well as between the hydration phenomena of NH4(+) and H3O(+).