The 50 most-cited articles in The Journal of Chemical Physics, 2022

We present our latest advancements of machine-learned potentials (MLPs) based on the neuroevolution potential (NEP) framework introduced in Fan et al. [Phys. Rev. B 104, 104309 (2021)] and their implementation in the open-source package gpumd. We increase the accuracy of NEP models both by improving the radial functions in the atomic-environment descriptor using a linear combination of Chebyshev basis functions and by extending the angular descriptor with some four-body and five-body contributions as in the atomic cluster expansion approach. We also detail our efficient implementation of the NEP approach in graphics processing units as well as our workflow for the construction of NEP models and demonstrate their application in large-scale atomistic simulations. By comparing to state-of-the-art MLPs, we show that the NEP approach not only achieves above-average accuracy but also is far more computationally efficient. These results demonstrate that the gpumd package is a promising tool for solving challenging problems requiring highly accurate, large-scale atomistic simulations. To…

Machine learning models for the potential energy of multi-atomic systems, such as the deep potential (DP) model, make molecular simulations with the accuracy of quantum mechanical density functional theory possible at a cost only moderately higher than that of empirical force fields. However, the majority of these models lack explicit long-range interactions and fail to describe properties that derive from the Coulombic tail of the forces. To overcome this limitation, we extend the DP model by approximating the long-range electrostatic interaction between ions (nuclei + core electrons) and valence electrons with that of distributions of spherical Gaussian charges located at ionic and electronic sites. The latter are rigorously defined in terms of the centers of the maximally localized Wannier distributions, whose dependence on the local atomic environment is modeled accurately by a deep neural network. In the DP long-range (DPLR) model, the electrostatic energy of the Gaussian charge system is added to short-range interactions that are represented as in the standard DP model. The…

A new composite density functional theory (DFT) method is presented. It is based on ωB97X-V as one of the best-performing density functionals for the GMTKN55 thermochemistry database and completes the family of "3c" methods toward range-separated hybrid DFT. This method is consistently available for all elements up to Rn (Z = 1-86). Its further key ingredients are a polarized valence double-ζ (vDZP) Gaussian basis set, which was fully optimized in molecular DFT calculations, in combination with large-core effective core potentials and a specially adapted D4 dispersion correction. Unlike most existing double-ζ atomic orbital sets, vDZP shows only small basis set superposition errors (BSSEs) and can compete with standard sets of triple-ζ quality. Small residual BSSE effects are efficiently absorbed by the D4 damping scheme, which overall eliminates the need for an explicit treatment or empirical corrections for BSSE. Thorough tests on a variety of thermochemistry benchmark sets show that the new composite method, dubbed ωB97X-3c, is on par with or even outperforms standard hybrid DFT…

Crystal structure prediction has been a subject of topical interest but remains a substantial challenge especially for complex structures as it deals with the global minimization of the extremely rugged high-dimensional potential energy surface. In this paper, a symmetry-orientated divide-and-conquer scheme was proposed to construct a symmetry tree graph, where the entire search space is decomposed into a finite number of symmetry dependent subspaces. An artificial intelligence-based symmetry selection strategy was subsequently devised to select the low-lying subspaces with high symmetries for global exploration and in-depth exploitation. Our approach can significantly simplify the problem of crystal structure prediction by avoiding exploration of the most complex P1 subspace on the entire search space and has the advantage of preserving the crystal symmetry during structure evolution, making it well suitable for predicting the complex crystal structures. The effectiveness of the method has been validated by successful prediction of the candidate structures of binary Lennard-Jones…

Kinetic Monte Carlo (KMC) simulations in combination with first-principles (1p)-based calculations are rapidly becoming the gold-standard computational framework for bridging the gap between the wide range of length scales and time scales over which heterogeneous catalysis unfolds. 1p-KMC simulations provide accurate insights into reactions over surfaces, a vital step toward the rational design of novel catalysts. In this Perspective, we briefly outline basic principles, computational challenges, successful applications, as well as future directions and opportunities of this promising and ever more popular kinetic modeling approach.

. These ions were the remaining alkaline and halogen ions, not previously considered in the Madrid-2019 force field. The force field, denoted as Madrid-2019-Extended, does not include polarizability and uses the TIP4P/2005 model of water and scaled charges for the ions. A charge of ±0.85e is assigned to monovalent ions. The force field developed provides an accurate description of aqueous solution densities over a wide range of concentrations up to the solubility limit of each salt studied. Good predictions of viscosity and diffusion coefficients are obtained for concentrations below 2 m. Structural properties obtained with this force field are also in reasonable agreement with the experiment. The number of contact ion pairs has been controlled to be low so as to avoid precipitation of the system at concentrations close to the experimental solubility limit. A comprehensive comparison of the performance for aqueous solutions of alkaline halides of force fields of electrolytes using scaled and integer charges is now possible. This comparison will help in the future to learn about the…

The regularized and restored semilocal meta-generalized gradient approximation (meta-GGA) exchange–correlation functional r2SCAN [Furness et al., J. Phys. Chem. Lett. 11, 8208–8215 (2020)] is used to create three global hybrid functionals with varying admixtures of Hartree–Fock “exact” exchange (HFX). The resulting functionals r2SCANh (10% HFX), r2SCAN0 (25% HFX), and r2SCAN50 (50% HFX) are combined with the semi-classical D4 London dispersion correction. The new functionals are assessed for the calculation of molecular geometries, main-group, and metalorganic thermochemistry at 26 comprehensive benchmark sets. These include the extensive GMTKN55 database, ROST61, and IONPI19 sets. It is shown that a moderate admixture of HFX leads to relative improvements of the mean absolute deviations for thermochemistry of 11% (r2SCANh-D4), 16% (r2SCAN0-D4), and 1% (r2SCAN50-D4) compared to the parental semi-local meta-GGA. For organometallic reaction energies and barriers, r2SCAN0-D4 yields an even larger mean improvement of 35%. The computation of structural parameters (geometry optimization)…

This Perspective provides a brief introduction into the theoretical complexity of polaritonic chemistry, which emerges from the hybrid nature of strongly coupled light-matter states. To tackle this complexity, the importance of ab initio methods is highlighted. Based on those, novel ideas and research avenues are developed with respect to quantum collectivity, as well as for resonance phenomena immanent in reaction rates under vibrational strong coupling. Indeed, fundamental theoretical questions arise about the mesoscopic scale of quantum-collectively coupled molecules when considering the depolarization shift in the interpretation of experimental data. Furthermore, to rationalize recent findings based on quantum electrodynamical density-functional theory (QEDFT), a simple, but computationally efficient, Langevin framework is proposed based on well-established methods from molecular dynamics. It suggests the emergence of cavity-induced non-equilibrium nuclear dynamics, where thermal (stochastic) resonance phenomena could emerge in the absence of external periodic driving. Overall,…

Constant potential methods (CPMs) enable computationally efficient simulations of the solid-liquid interface at conducting electrodes in molecular dynamics. They have been successfully used, for example, to realistically model the behavior of ionic liquids or water-in-salt electrolytes in supercapacitors and batteries. CPMs model conductive electrodes by updating charges of individual electrode atoms according to the applied electric potential and the (time-dependent) local electrolyte structure. Here, we present a feature-rich CPM implementation, called ELECTRODE, for the Large-scale Atomic/Molecular Massively Parallel Simulator, which includes a constrained charge method and a thermo-potentiostat. The ELECTRODE package also contains a finite-field approach, multiple corrections for nonperiodic boundary conditions of the particle-particle particle-mesh solver, and a Thomas-Fermi model for using nonideal metals as electrodes. We demonstrate the capabilities of this implementation for a parallel-plate electrical double-layer capacitor, for which we have investigated the charging…

There is a perceived dichotomy between structure-based and descriptor-based molecular representations used for predictive chemistry tasks. Here, we study the performance, generalizability, and explainability of the quantum mechanics-augmented graph neural network (ml-QM-GNN) architecture as applied to the prediction of regioselectivity (classification) and of activation energies (regression). In our hybrid QM-augmented model architecture, structure-based representations are first used to predict a set of atom- and bond-level reactivity descriptors derived from density functional theory calculations. These estimated reactivity descriptors are combined with the original structure-based representation to make the final reactivity prediction. We demonstrate that our model architecture leads to significant improvements over structure-based GNNs in not only overall accuracy but also in generalization to unseen compounds. Even when provided training sets of only a couple hundred labeled data points, the ml-QM-GNN outperforms other state-of-the-art structure-based architectures that have…

In this work, we present a general purpose deep neural network package for representing energies, forces, dipole moments, and polarizabilities of atomistic systems. This so-called recursively embedded atom neural network model takes advantages of both the physically inspired atomic descriptor based neural networks and the message-passing based neural networks. Implemented in the PyTorch framework, the training process is parallelized on both the central processing unit and the graphics processing unit with high efficiency and low memory in which all hyperparameters can be optimized automatically. We demonstrate the state-of-the-art accuracy, high efficiency, scalability, and universality of this package by learning not only energies (with or without forces) but also dipole moment vectors and polarizability tensors in various molecular, reactive, and periodic systems. An interface between a trained model and LAMMPs is provided for large scale molecular dynamics simulations. We hope that this open-source toolbox will allow for future method development and applications of machine…

Salt aqueous solutions are relevant in many fields, ranging from biological systems to seawater. Thus, the availability of a force-field that is able to reproduce the thermodynamic and dynamic behavior of salt aqueous solutions would be of great interest. Unfortunately, this has been proven challenging, and most of the existing force-fields fail to reproduce much of their behavior. In particular, the diffusion of water or the salt solubility are often not well reproduced by most of the existing force-fields. Recently, the Madrid-2019 model was proposed, and it was shown that this force-field, which uses the TIP4P/2005 model for water and non-integer charges for the ions, provides a good description of a large number of properties, including the solution densities, viscosities, and the diffusion of water. In this work, we assess the performance of this force-field on the evaluation of the freezing point depression. Although the freezing point depression is a colligative property that at low salt concentrations depends solely on properties of pure water, a good model for the…

In this Communication, we propose the time-domain Prony fitting decomposition (t-PFD) as an accurate and efficient exponential series method, applicable to arbitrary bath correlation functions. The resulting numerical efficiency of hierarchical equations of motion (HEOM) formalism is greatly optimized, especially in low temperature regimes that would be inaccessible with other methods. For demonstration, we calibrate the present t-PFD against the celebrated Padé spectrum decomposition method, followed by converged HEOM evaluations on the single-impurity Anderson model system.

Based on self-consistent field (SCF) atomic mean-field (amf) quantities, we present two simple yet computationally efficient and numerically accurate matrix-algebraic approaches to correct both scalar-relativistic and spin–orbit two-electron picture-change effects (PCEs) arising within an exact two-component (X2C) Hamiltonian framework. Both approaches, dubbed amfX2C and e(xtended)amfX2C, allow us to uniquely tailor PCE corrections to mean-field models, viz. Hartree–Fock or Kohn–Sham DFT, in the latter case also avoiding the need for a point-wise calculation of exchange–correlation PCE corrections. We assess the numerical performance of these PCE correction models on spinor energies of group 18 (closed-shell) and group 16 (open-shell) diatomic molecules, achieving a consistent ≈10−5 Hartree accuracy compared to reference four-component data. Additional tests include SCF calculations of molecular properties such as absolute contact density and contact density shifts in copernicium fluoride compounds (CnFn, n = 2,4,6), as well as equation-of-motion coupled-cluster calculations of…

Cross-linking is known to play a pivotal role in the relaxation dynamics and mechanical properties of thermoset polymers, which are commonly used in structural applications because of their light weight and inherently strong nature. Here, we employ a coarse-grained (CG) polymer model to systematically explore the effect of cross-link density on basic thermodynamic properties as well as corresponding changes in the segmental dynamics and elastic properties of these network materials upon approaching their glass transition temperatures (Tg). Increasing the cross-link density unsurprisingly leads to a significant slowing down of the segmental dynamics, and the fragility K of glass formation shifts in lockstep with Tg, as often found in linear polymer melts when the polymer mass is varied. As a consequence, the segmental relaxation time τα becomes almost a universal function of reduced temperature, (T − Tg)/Tg, a phenomenon that underlies the applicability of the “universal” Williams–Landel–Ferry (WLF) relation to many polymer materials. We also test a mathematical model of the…

Electrochemistry is central to many applications, ranging from biology to energy science. Studies now involve a wide range of techniques, both experimental and theoretical. Modeling and simulations methods, such as density functional theory or molecular dynamics, provide key information on the structural and dynamic properties of the systems. Of particular importance are polarization effects of the electrode/electrolyte interface, which are difficult to simulate accurately. Here, we show how these electrostatic interactions are taken into account in the framework of the Ewald summation method. We discuss, in particular, the formal setup for calculations that enforce periodic boundary conditions in two directions, a geometry that more closely reflects the characteristics of typical electrolyte/electrode systems and presents some differences with respect to the more common case of periodic boundary conditions in three dimensions. These formal developments are implemented and tested in MetalWalls, a molecular dynamics software that captures the polarization of the electrolyte and…

A water microdroplet in air or oil typically possesses an electric double layer (EDL) from the preferential adsorption of surface-bound ions at the periphery. We present the calculations of the ion gradients within a microdroplet at equilibrium, including systems containing buffers and water autoionization. These ion gradients are used to calculate the potential energy stored within the microdroplet. We consider how this stored potential energy can be utilized to drive chemical reactions, much like an electrochemical cell. Effective voltages as high as 111 mV are found for microdroplets having a low surface charge density (0.01 ions per nm2). Two sources of potential energy are investigated: (1) the electrostatic energy of the EDL of the microdroplet and (2) shifts in other chemical equilibria coupled to the main reaction through the EDL. A particularly important example of the latter is water autoionization, wherein the reaction of interest causes a flattening of the [H+] gradient within the EDL, resulting in a net recombination of H+ and OH− throughout the microdroplet. Numerical…

We combine ab initio molecular electronic Hamiltonians with a cavity quantum electrodynamics model for dissipative photonic modes and apply mean-field theories to the ground- and excited-states of resulting polaritonic systems. In particular, we develop a non-Hermitian configuration interaction singles theory for mean-field ground- and excited-states of the molecular system strongly interacting with a photonic mode and apply these methods to elucidating the phenomenology of paradigmatic polaritonic systems. We leverage the Psi4Numpy framework to yield open-source and accessible reference implementations of these methods.

One of the most intriguing phenomena in active matter has been the gas-liquid-like motility-induced phase separation (MIPS) observed in repulsive active particles. However, experimentally, no particle can be a perfect sphere, and the asymmetric shape, mass distribution, or catalysis coating can induce an active torque on the particle, which makes it a chiral active particle. Here, using computer simulations and dynamic mean-field theory, we demonstrate that the large enough torque of circle active Brownian particles in two dimensions generates a dynamical clustering state interrupting the conventional MIPS. Multiple clusters arise from the combination of the conventional MIPS cohesion, and the circulating current caused disintegration. The nonvanishing current in non-equilibrium steady states microscopically originates from the motility "relieved" by automatic rotation, which breaks the detailed balance at the continuum level. This suggests that no equilibrium-like phase separation theory can be constructed for chiral active colloids even with tiny active torque, in which no visible…

Understanding the thermodynamic stability and metastability of materials can help us to, for example, gauge whether crystalline polymorphs in pharmaceutical formulations are likely to be durable. It can also help us to design experimental routes to novel phases with potentially interesting properties. In this Perspective, we provide an overview of how thermodynamic phase behavior can be quantified both in computer simulations and machine-learning approaches to determine phase diagrams, as well as combinations of the two. We review the basic workflow of free-energy computations for condensed phases, including some practical implementation advice, ranging from the Frenkel-Ladd approach to thermodynamic integration and to direct-coexistence simulations. We illustrate the applications of such methods on a range of systems from materials chemistry to biological phase separation. Finally, we outline some challenges, questions, and practical applications of phase-diagram determination which we believe are likely to be possible to address in the near future using such state-of-the-art…

The function of light-harvesting complexes is determined by a complex network of dynamic interactions among all the different components: the aggregate of pigments, the protein, and the surrounding environment. Complete and reliable predictions on these types of composite systems can be only achieved with an atomistic description. In the last few decades, there have been important advances in the atomistic modeling of light-harvesting complexes. These advances have involved both the completeness of the physical models and the accuracy and effectiveness of the computational protocols. In this Perspective, we present an overview of the main theoretical and computational breakthroughs attained so far in the field, with particular focus on the important role played by the protein and its dynamics. We then discuss the open problems in their accurate modeling that still need to be addressed. To illustrate an effective computational workflow for the modeling of light harvesting complexes, we take as an example the plant antenna complex CP29 and its H111N mutant.

Molecular Dynamics (MD) simulation is a powerful tool for understanding the dynamics and structure of matter. Since the resolution of MD is atomic-scale, achieving long timescale simulations with femtosecond integration is very expensive. In each MD step, numerous iterative computations are performed to calculate energy based on different types of interaction and their corresponding spatial gradients. These repetitive computations can be learned and surrogated by a deep learning model, such as a Graph Neural Network (GNN). In this work, we developed a GNN Accelerated MD (GAMD) model that directly predicts forces, given the state of the system (atom positions, atom types), bypassing the evaluation of potential energy. By training the GNN on a variety of data sources (simulation data derived from classical MD and density functional theory), we show that GAMD can predict the dynamics of two typical molecular systems, Lennard-Jones system and water system, in the NVT ensemble with velocities regulated by a thermostat. We further show that GAMD's learning and inference are agnostic to…

We report the implementation of effective quantum electrodynamics (QED) potentials for all-electron four-component relativistic molecular calculations using the DIRAC code. The potentials are also available for two-component calculations, being properly picture-change transformed. The latter point is important; we demonstrate through atomic calculations that picture-change errors are sizable. Specifically, we have implemented the Uehling potential [E. A. Uehling, Phys. Rev. 48, 55 (1935)] for vacuum polarization and two effective potentials [P. Pyykkö and L.-B. Zhao, J. Phys. B: At., Mol. Opt. Phys. 36, 1469 (2003) and V. V. Flambaum and J. S. M. Ginges, Phys. Rev. A 72, 052115 (2005)] for electron self-energy. We provide extensive theoretical background for these potentials, hopefully reaching an audience beyond QED specialists. We report the following sample applications: (i) We first confirm the conjecture of P. Pyykkö that QED effects are observable for the AuCN molecule by directly calculating ground-state rotational constants B0 of the three isotopomers studied by microwave…

Kohn-Sham density-functional theory (DFT), the predominant framework for electronic structure computations in chemistry today, has undergone considerable evolution in the past few decades. The earliest DFT approximations were based on uniform electron gas models completely free of empirical parameters. Tremendous improvements were made by incorporating density gradients and a small number of parameters, typically one or two, obtained from fits to atomic data. Incorporation of exact exchange and fitting to molecular data, such as experimental heats of formation, allowed even further improvements. This, however, opened a Pandora's Box of fitting possibilities, given the limitless choices of chemical reactions that can be fit. The result is a recent explosion of DFT approximations empirically fit to hundreds, or thousands, of chemical reference data. These fitted density functionals may contain several dozen empirical parameters. What has been lost in this fitting trend is physical modeling based on theory. In this work, we present a density functional comprising our best efforts to…

Conformational dynamics of biomolecules are of fundamental importance for their function. Single-molecule studies of Förster Resonance Energy Transfer (smFRET) between a tethered donor and acceptor dye pair are a powerful tool to investigate the structure and dynamics of labeled molecules. However, capturing and quantifying conformational dynamics in intensity-based smFRET experiments remains challenging when the dynamics occur on the sub-millisecond timescale. The method of multiparameter fluorescence detection addresses this challenge by simultaneously registering fluorescence intensities and lifetimes of the donor and acceptor. Together, two FRET observables, the donor fluorescence lifetime τD and the intensity-based FRET efficiency E, inform on the width of the FRET efficiency distribution as a characteristic fingerprint for conformational dynamics. We present a general framework for analyzing dynamics that relates average fluorescence lifetimes and intensities in two-dimensional burst frequency histograms. We present parametric relations of these observables for interpreting…

We propose a new overarching model for self-propelled particles that flexibly generates a full family of "descendants." The general dynamics introduced in this paper, which we denote as the "parental" active model (PAM), unifies two special cases commonly used to describe active matter, namely, active Brownian particles (ABPs) and active Ornstein-Uhlenbeck particles (AOUPs). We thereby document the existence of a deep and close stochastic relationship between them, resulting in the subtle balance between fluctuations in the magnitude and direction of the self-propulsion velocity. Besides illustrating the relation between these two common models, the PAM can generate additional offsprings, interpolating between ABP and AOUP dynamics, that could provide more suitable models for a large class of living and inanimate active matter systems, possessing characteristic distributions of their self-propulsion velocity. Our general model is evaluated in the presence of a harmonic external confinement. For this reference example, we present a two-state phase diagram that sheds light on the…

We develop a quartic-scaling implementation of coupled-cluster singles and doubles (CCSD) based on low-rank tensor hypercontraction (THC) factorizations of both the electron repulsion integrals (ERIs) and the doubles amplitudes. This extends our rank-reduced (RR) coupled-cluster method to incorporate higher-order tensor factorizations. The THC factorization of the doubles amplitudes accounts for most of the gain in computational efficiency as it is sufficient, in conjunction with a Cholesky decomposition of the ERIs, to reduce the computational complexity of most contributions to the CCSD amplitude equations. Further THC factorization of the ERIs reduces the complexity of certain terms arising from nested commutators between the doubles excitation operator and the two-electron operator. We implement this new algorithm using graphical processing units and demonstrate that it enables CCSD calculations for molecules with 250 atoms and 2500 basis functions using a single computer node. Furthermore, we show that the new method computes correlation energies with comparable accuracy to the…

Generalized Langevin equations with non-linear forces and position-dependent linear friction memory kernels, such as commonly used to describe the effective dynamics of coarse-grained variables in molecular dynamics, are rigorously derived within the Mori-Zwanzig formalism. A fluctuation-dissipation theorem relating the properties of the noise to the memory kernel is shown. The derivation also yields Volterra-type equations for the kernel, which can be used for a numerical parametrization of the model from all-atom simulations.

Molecular dynamics (MD) simulations of complex electrochemical systems, such as ionic liquid supercapacitors, are increasingly including the constant potential method (CPM) to model conductive electrodes at a specified potential difference, but the inclusion of CPM can be computationally expensive. We demonstrate the computational savings available in CPM MD simulations of ionic liquid supercapacitors when the usual non-periodic slab geometry is replaced with fully periodic boundary conditions. We show how a doubled cell approach, previously used in non-CPM MD simulations of charged interfaces, can be used to enable fully periodic CPM MD simulations. Using either a doubled cell approach or a finite field approach previously reported by others, fully periodic CPM MD simulations produce comparable results to the traditional slab geometry simulations with a nearly double speedup in computational time. Indeed, these savings can offset the additional cost of the CPM algorithm, resulting in periodic CPM MD simulations that are computationally competitive with the non-periodic, fixed…

Automatic differentiation represents a paradigm shift in scientific programming, where evaluating both functions and their derivatives is required for most applications. By removing the need to explicitly derive expressions for gradients, development times can be shortened and calculations can be simplified. For these reasons, automatic differentiation has fueled the rapid growth of a variety of sophisticated machine learning techniques over the past decade, but is now also increasingly showing its value to support ab initio simulations of quantum systems and enhance computational quantum chemistry. Here, we present an open-source differentiable quantum chemistry simulation code and explore applications facilitated by automatic differentiation: (1) calculating molecular perturbation properties, (2) reoptimizing a basis set for hydrocarbons, (3) checking the stability of self-consistent field wave functions, and (4) predicting molecular properties via alchemical perturbations.

Tensor network decompositions of path integrals for simulating open quantum systems have recently been proven to be useful. However, these methods scale exponentially with the system size. This makes it challenging to simulate the non-equilibrium dynamics of extended quantum systems coupled with local dissipative environments. In this work, we extend the tensor network path integral (TNPI) framework to efficiently simulate such extended systems. The Feynman-Vernon influence functional is a popular approach used to account for the effect of environments on the dynamics of the system. In order to facilitate the incorporation of the influence functional into a multisite framework (MS-TNPI), we combine a matrix product state (MPS) decomposition of the reduced density tensor of the system along the sites with a corresponding tensor network representation of the time axis to construct an efficient 2D tensor network. The 2D MS-TNPI network, when contracted, yields the time-dependent reduced density tensor of the extended system as an MPS. The algorithm presented is independent of the…

Finding a low dimensional representation of data from long-timescale trajectories of biomolecular processes, such as protein folding or ligand-receptor binding, is of fundamental importance, and kinetic models, such as Markov modeling, have proven useful in describing the kinetics of these systems. Recently, an unsupervised machine learning technique called VAMPNet was introduced to learn the low dimensional representation and the linear dynamical model in an end-to-end manner. VAMPNet is based on the variational approach for Markov processes and relies on neural networks to learn the coarse-grained dynamics. In this paper, we combine VAMPNet and graph neural networks to generate an end-to-end framework to efficiently learn high-level dynamics and metastable states from the long-timescale molecular dynamics trajectories. This method bears the advantages of graph representation learning and uses graph message passing operations to generate an embedding for each datapoint, which is used in the VAMPNet to generate a coarse-grained dynamical model. This type of molecular representation…

a Hamiltonian vibration vs a bath vibration on the 2D electronic spectra of a homodimer

The hypothesis that the anomalous behavior of liquid water is related to the existence of a second critical point in deeply supercooled states has long been the subject of intense debate. Recent, sophisticated experiments designed to observe the transformation between the two subcritical liquids on nano- and microsecond time scales, along with demanding numerical simulations based on classical (rigid) models parameterized to reproduce thermodynamic properties of water, have provided support to this hypothesis. A stronger numerical proof requires demonstrating that the critical point, which occurs at temperatures and pressures far from those at which the models were optimized, is robust with respect to model parameterization, specifically with respect to incorporating additional physical effects. Here, we show that a liquid-liquid critical point can be rigorously located also in the WAIL model of water [Pinnick et al., J. Chem. Phys. 137, 014510 (2012)], a model parameterized using ab initio calculations only. The model incorporates two features not present in many previously studied…

The electron attachment variant of equation-of-motion coupled-cluster theory (EOM-EA-CC) is generalized to the case of strong light-matter coupling within the framework of cavity quantum electrodynamics (QED). The resulting EOM-EA-QED-CC formalism provides an ab initio, correlated, and non-perturbative description of cavity-induced effects in many-electron systems that complements other recently proposed cavity-QED-based extensions of CC theory. Importantly, this work demonstrates that QED generalizations of EOM-CC theory are useful frameworks for exploring particle-non-conserving sectors of Fock space, thereby establishing a path forward for the simultaneous description of both strong electron-electron and electron-photon correlation effects.

Permutationally invariant polynomial (PIP) regression has been used to obtain machine-learned potential energy surfaces, including analytical gradients, for many molecules and chemical reactions. Recently, the approach has been extended to moderate size molecules with up to 15 atoms. The algorithm, including "purification of the basis," is computationally efficient for energies; however, we found that the recent extension to obtain analytical gradients, despite being a remarkable advance over previous methods, could be further improved. Here, we report developments to further compact a purified basis and, more significantly, to use the reverse differentiation approach to greatly speed up gradient evaluation. We demonstrate this for our recent four-body water interaction potential. Comparisons of training and testing precision on the MD17 database of energies and gradients (forces) for ethanol against numerous machine-learning methods, which were recently assessed by Dral and co-workers, are given. The PIP fits are as precise as those using these methods, but the PIP computation time…

Coarse-grained (CG) models facilitate an efficient exploration of complex systems by reducing the unnecessary degrees of freedom of the fine-grained (FG) system while recapitulating major structural correlations. Unlike structural properties, assessing dynamic properties in CG modeling is often unfeasible due to the accelerated dynamics of the CG models, which allows for more efficient structural sampling. Therefore, the ultimate goal of the present series of articles is to establish a better correspondence between the FG and CG dynamics. To assess and compare dynamical properties in the FG and the corresponding CG models, we utilize the excess entropy scaling relationship. For Paper I of this series, we provide evidence that the FG and the corresponding CG counterpart follow the same universal scaling relationship. By carefully reviewing and examining the literature, we develop a new theory to calculate excess entropies for the FG and CG systems while accounting for entropy representability. We demonstrate that the excess entropy scaling idea can be readily applied to liquid water…

Neural network (NN) potentials are a natural choice for coarse-grained (CG) models. Their many-body capacity allows highly accurate approximations of the potential of mean force, promising CG simulations of unprecedented accuracy. CG NN potentials trained bottom-up via force matching (FM), however, suffer from finite data effects: They rely on prior potentials for physically sound predictions outside the training data domain, and the corresponding free energy surface is sensitive to errors in the transition regions. The standard alternative to FM for classical potentials is relative entropy (RE) minimization, which has not yet been applied to NN potentials. In this work, we demonstrate, for benchmark problems of liquid water and alanine dipeptide, that RE training is more data efficient, due to accessing the CG distribution during training, resulting in improved free energy surfaces and reduced sensitivity to prior potentials. In addition, RE learns to correct time integration errors, allowing larger time steps in CG molecular dynamics simulation, while maintaining accuracy. Thus,…

Since Arrhenius first proposed an equation to account for the behavior of thermally activated reactions in 1889, significant progress has been made in our understanding of chemical reactivity. A number of capture theory models have been developed over the past several decades to predict the rate coefficients for reactions between ions and molecules-ranging from the Langevin equation (for reactions between ions and non-polar molecules) to more recent fully quantum theories (for reactions at ultracold temperatures). A number of different capture theory methods are discussed, with the key assumptions underpinning each approach clearly set out. The strengths and limitations of these capture theory methods are examined through detailed comparisons between low-temperature experimental measurements and capture theory predictions. Guidance is provided on the selection of an appropriate capture theory method for a given class of ion-molecule reaction and set of experimental conditions-identifying when a capture-based model is likely to provide an accurate prediction. Finally, the impact of…

We extend the twin-space formulation of the hierarchical equations of motion approach in combination with the matrix product state representation [R. Borrelli, J. Chem. Phys. 150, 234102 (2019)] to nonequilibrium scenarios where the open quantum system is coupled to a hybrid fermionic and bosonic environment. The key ideas used in the extension are a reformulation of the hierarchical equations of motion for the auxiliary density matrices into a time-dependent Schrödinger-like equation for an augmented multi-dimensional wave function as well as a tensor decomposition into a product of low-rank matrices. The new approach facilitates accurate simulations of non-equilibrium quantum dynamics in larger and more complex open quantum systems. The performance of the method is demonstrated for a model of a molecular junction exhibiting current-induced mode-selective vibrational excitation.

The kinetics of a dynamical system dominated by two metastable states is examined from the perspective of the activated-dynamics reactive flux formalism, Markov state eigenvalue spectral decomposition, and committor-based transition path theory. Analysis shows that the different theoretical formulations are consistent, clarifying the significance of the inherent microscopic lag-times that are implicated, and that the most meaningful one-dimensional reaction coordinate in the region of the transition state is along the gradient of the committor in the multidimensional subspace of collective variables. It is shown that the familiar reactive flux activated dynamics formalism provides an effective route to calculate the transition rate in the case of a narrow sharp barrier but much less so in the case of a broad flat barrier. In this case, the standard reactive flux correlation function decays very slowly to the plateau value that corresponds to the transmission coefficient. Treating the committor function as a reaction coordinate does not alleviate all issues caused by the slow…

The vibrational spectra of condensed and gas-phase systems are influenced by thequantum-mechanical behavior of light nuclei. Full-dimensional simulations of approximate quantum dynamics are possible thanks to the imaginary time path-integral (PI) formulation of quantum statistical mechanics, albeit at a high computational cost which increases sharply with decreasing temperature. By leveraging advances in machine-learned coarse-graining, we develop a PI method with the reduced computational cost of a classical simulation. We also propose a simple temperature elevation scheme to significantly attenuate the artifacts of standard PI approaches as well as eliminate the unfavorable temperature scaling of the computational cost. We illustrate the approach, by calculating vibrational spectra using standard models of water molecules and bulk water, demonstrating significant computational savings and dramatically improved accuracy compared to more expensive reference approaches. Our simple, efficient, and accurate method has prospects for routine calculations of vibrational spectra for a wide…

The ab initio GW plus Bethe–Salpeter equation (GW-BSE, where G is the one particle Green's function and W is the screened Coulomb interaction) approach has emerged as a leading method for predicting excitations in both solids and molecules with a predictive power contingent upon several factors. Among these factors are the (1) generalized Kohn–Sham eigensystem used to construct the GW self-energy and to solve the BSE and (2) the efficacy and suitability of the Tamm–Dancoff approximation. Here, we present a detailed benchmark study of low-lying singlet excitations from a generalized Kohn–Sham (gKS) starting point based on an optimally tuned range-separated hybrid (OTRSH) functional. We show that the use of this gKS starting point with one-shot G0W0 and G0W0-BSE leads to the lowest mean absolute errors (MAEs) and mean signed errors (MSEs), with respect to high-accuracy reference values, demonstrated in the literature thus far for the ionization potentials of the GW100 benchmark set and for low-lying neutral excitations of Thiel’s set molecules in the gas phase, without the need for…

By using the direct coexistence method, we have calculated the melting points of ice Ih at normal pressure for three recently proposed water models, namely, TIP3P-FB, TIP4P-FB, and TIP4P-D. We obtained Tm = 216 K for TIP3P-FB, Tm = 242 K for TIP4P-FB, and Tm = 247 K for TIP4P-D. We revisited the melting point of TIP4P/2005 and TIP5P obtaining Tm = 250 and 274 K, respectively. We summarize the current situation of the melting point of ice Ih for a number of water models and conclude that no model is yet able to simultaneously reproduce the melting temperature of ice Ih and the temperature of the maximum in density at room pressure. This probably points toward our both still incomplete knowledge of the potential energy surface of water and the necessity of incorporating nuclear quantum effects to describe both properties simultaneously.

Quantum chemical calculations using density functional theory are reported for the diatomic molecules LiF, BeO, and BN. The nature of the interatomic interactions is analyzed with the Energy Decomposition Analysis–Natural Orbitals of Chemical Valence (EDA-NOCV) method, and the results are critically discussed and compared with data from Quantum Theory of Atoms in Molecules, Natural Bond Orbital, and Mayer approaches. Polar bonds, like nonpolar bonds, are caused by the interference of wave functions, which lead to an accumulation of electronic charge in the bonding region. Polar bonds generally have a larger percentage of electrostatic bonding to the total attraction, but nonpolar bonds may also possess large contributions from Coulombic interaction. The term “ionic contribution” refers to valence bond structures and is misleading because it refers to separate fragments with negligible overlap that occur only in the solid state and in solution, not in a molecule. The EDA-NOCV method gives detailed information about the individual orbital contributions, which can be identified by…

A method for obtaining appropriate reaction coordinates is required to identify transition states distinguishing the product and reactant in complex molecular systems. Recently, abundant research has been devoted to obtaining reaction coordinates using artificial neural networks from deep learning literature, where many collective variables are typically utilized in the input layer. However, it is difficult to explain the details of which collective variables contribute to the predicted reaction coordinates owing to the complexity of the nonlinear functions in deep neural networks. To overcome this limitation, we used Explainable Artificial Intelligence (XAI) methods of the Local Interpretable Model-agnostic Explanation (LIME) and the game theory-based framework known as Shapley Additive exPlanations (SHAP). We demonstrated that XAI enables us to obtain the degree of contribution of each collective variable to reaction coordinates that is determined by nonlinear regressions with deep learning for the committor of the alanine dipeptide isomerization in vacuum. In particular, both…

Data-driven schemes that associate molecular and crystal structures with their microscopic properties share the need for a concise, effective description of the arrangement of their atomic constituents. Many types of models rely on descriptions of atom-centered environments, which are associated with an atomic property or with an atomic contribution to an extensive macroscopic quantity. Frameworks in this class can be understood in terms of atom-centered density correlations (ACDC), which are used as a basis for a body-ordered, symmetry-adapted expansion of the targets. Several other schemes that gather information on the relationship between neighboring atoms using "message-passing" ideas cannot be directly mapped to correlations centered around a single atom. We generalize the ACDC framework to include multi-centered information, generating representations that provide a complete linear basis to regress symmetric functions of atomic coordinates, and provide a coherent foundation to systematize our understanding of both atom-centered and message-passing and invariant and…

We introduce an extension to the PySCF package, which makes it automatically differentiable. The implementation strategy is discussed, and example applications are presented to demonstrate the automatic differentiation framework for quantum chemistry methodology development. These include orbital optimization, properties, excited-state energies, and derivative couplings, at the mean-field level and beyond, in both molecules and solids. We also discuss some current limitations and directions for future work.

Most computational studies in chemistry and materials science are based on the use of density functional theory. Although the exact density functional is unknown, several density functional approximations (DFAs) offer a good balance of affordable computational cost and semi-quantitative accuracy for applications. The development of DFAs still continues on many fronts, and several new DFAs aiming for improved accuracy are published every year. However, the numerical behavior of these DFAs is an often-overlooked problem. In this work, we look at all 592 DFAs for three-dimensional systems available in Libxc 5.2.2 and examine the convergence of the density functional total energy based on tabulated atomic Hartree–Fock wave functions. We show that several recent DFAs, including the celebrated SCAN family of functionals, show impractically slow convergence with typically used numerical quadrature schemes, making these functionals unsuitable both for routine applications and high-precision studies, as thousands of radial quadrature points may be required to achieve sub-μEh accurate total…

By generalizing the recently developed path integral molecular dynamics for identical bosons and fermions, we consider the finite-temperature thermodynamic properties of fictitious identical particles with a real parameter ξ interpolating continuously between bosons (ξ = 1) and fermions (ξ = -1). Through general analysis and numerical experiments, we find that the average energy may have good analytical properties as a function of this real parameter ξ, which provides the chance to calculate the thermodynamical properties of identical fermions by extrapolation with a simple polynomial function after accurately calculating the thermodynamic properties of the fictitious particles for ξ ≥ 0. Using several examples, it is shown that our method can efficiently give accurate energy values for finite-temperature fermionic systems. Our work provides a chance to circumvent the fermion sign problem for some quantum systems.