Extended Lagrangian Born-Oppenheimer molecular dynamics with dissipation

Shadow Molecular Dynamics for a Charge-Potential Equilibration Model

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

Shadow Molecular Dynamics with a Machine Learned Flexible Charge Potential

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

Grassmann Extrapolation for Accelerating Geometry Optimization

This study extends the Grassmann extrapolation (G-Ext) method, which was introduced for Born-Oppenheimer molecular dynamics, to the context of geometry optimization. Using density matrices from previous optimization steps, the G-Ext approach applies a nonlinear, structure-preserving mapping onto the Grassmann manifold to provide an initial guess which accelerates the convergence of the self-consistent field (SCF) procedure. Using the optimal parameters identified by employing various descriptors and computational strategies across a diverse set of molecules, G-Ext shows excellent performance improvements, particularly with large molecular systems.

Initial Conditions for Excited-State Dynamics in Solvated Systems: A Case Study

Simulating excited-state dynamics or computing spectra for molecules in condensed phases requires sampling the ground state to generate initial conditions. Initial conditions (or snapshots for spectra) are typically produced by QM/MM Boltzmann sampling following MM equilibration or optimization. Given the switch from a MM to a QM/MM potential energy surface, one should discard a set period of time (which we call the "healing time") from the beginning of the QM/MM trajectory. Ideally, the healing time is as short as possible (to avoid unnecessary computational effort), but long enough to equilibrate to the QM/MM ground state distribution. Healing times in previous studies range from tens of femtoseconds to tens of picoseconds, suggesting the need for guidelines to choose a healing time. We examine the effect of healing time on the nonadiabatic dynamics and spectrum of a first-generation Donor-Acceptor Stenhouse Adduct in chloroform. Insufficient healing times skew the branching ratio of ground state products and alter the relaxation time for one pathway. The influence of the healing time on the absorption spectrum is less pronounced, warning that the spectrum is not a sensitive indicator for the quality of a set of initial conditions for dynamics. We demonstrate that a reasonable estimate for the healing time can be obtained by monitoring the solute temperature during the healing trajectory. We suggest that this procedure should become standard practice for determining healing times to generate initial conditions for nonadiabatic QM/MM simulations in large molecules and condensed phases.

CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

Multiple-time scale integration method based on an interpolated potential energy surface for ab initio path integral molecular dynamics

A new multiple-time scale integration method is presented that propagates ab initio path integral molecular dynamics (PIMD). This method uses a large time step to generate an approximate geometrical configuration whose energy and gradient are evaluated at the level of an ab initio method, and then, a more precise integration scheme, e.g., the Bulirsch-Stoer method or velocity Verlet integration with a smaller time step, is used to integrate from the previous step using the computationally efficient interpolated potential energy surface constructed from two consecutive points. This method makes the integration of PIMD more efficient and accurate compared with the velocity Verlet integration. A Nosé-Hoover chain thermostat combined with this new multiple-time scale method has good energy conservation even with a large time step, which is usually challenging in velocity Verlet integration for PIMD due to the very small chain mass when a large number of beads are used. The new method is used to calculate infrared spectra and free energy profiles to demonstrate its accuracy and capabilities.

A Quasi Time-Reversible Scheme Based on Density Matrix Extrapolation on the Grassmann Manifold for Born-Oppenheimer Molecular Dynamics

This Letter introduces the so-called Quasi Time-Reversible scheme based on Grassmann extrapolation (QTR G-Ext) of density matrices for an accurate calculation of initial guesses in Born-Oppenheimer Molecular Dynamics (BOMD) simulations. The method shows excellent results on four large molecular systems that are representative of real-life production applications, ranging from 21 to 94 atoms simulated with Kohn-Sham (KS) density functional theory surrounded with a classical environment with 6k to 16k atoms. Namely, it clearly reduces the number of self-consistent field iterations while at the same time achieving energy-conserving simulations, resulting in a considerable speed-up of BOMD simulations even when tight convergence of the KS equations is required.

Semi-Empirical Shadow Molecular Dynamics: A PyTorch Implementation

Extended Lagrangian Born-Oppenheimer molecular dynamics (XL-BOMD) in its most recent shadow potential energy version has been implemented in the semiempirical PyTorch-based software PySeQM. The implementation includes finite electronic temperatures, canonical density matrix perturbation theory, and an adaptive Krylov subspace approximation for the integration of the electronic equations of motion within the XL-BOMB approach (KSA-XL-BOMD). The PyTorch implementation leverages the use of GPU and machine learning hardware accelerators for the simulations. The new XL-BOMD formulation allows studying more challenging chemical systems with charge instabilities and low electronic energy gaps. The current public release of PySeQM continues our development of modular architecture for large-scale simulations employing semi-empirical quantum-mechanical treatment. Applied to molecular dynamics, simulation of 840 carbon atoms, one integration time step executes in 4 s on a single Nvidia RTX A6000 GPU.

Shadow energy functionals and potentials in Born-Oppenheimer molecular dynamics

In Born-Oppenheimer molecular dynamics (BOMD) simulations based on the density functional theory (DFT), the potential energy and the interatomic forces are calculated from an electronic ground state density that is determined by an iterative self-consistent field optimization procedure, which, in practice, never is fully converged. The calculated energies and forces are, therefore, only approximate, which may lead to an unphysical energy drift and instabilities. Here, we discuss an alternative shadow BOMD approach that is based on backward error analysis. Instead of calculating approximate solutions for an underlying exact regular Born-Oppenheimer potential, we do the opposite. Instead, we calculate the exact electron density, energies, and forces, but for an underlying approximate shadow Born-Oppenheimer potential energy surface. In this way, the calculated forces are conservative with respect to the approximate shadow potential and generate accurate molecular trajectories with long-term energy stabilities. We show how such shadow Born-Oppenheimer potentials can be constructed at different levels of accuracy as a function of the integration time step, δt, from the constrained minimization of a sequence of systematically improvable, but approximate, shadow energy density functionals. For each energy functional, there is a corresponding ground state Born-Oppenheimer potential. These pairs of shadow energy functionals and potentials are higher-level generalizations of the original "zeroth-level" shadow energy functionals and potentials used in extended Lagrangian BOMD [Niklasson, Eur. Phys. J. B 94, 164 (2021)]. The proposed shadow energy functionals and potentials are useful only within this extended dynamical framework, where also the electronic degrees of freedom are propagated as dynamical field variables together with the atomic positions and velocities. The theory is quite general and can be applied to MD simulations using approximate DFT, Hartree-Fock, or semi-empirical methods, as well as to coarse-grained flexible charge models.

Graph-based quantum response theory and shadow Born-Oppenheimer molecular dynamics

Graph-based linear scaling electronic structure theory for quantum-mechanical molecular dynamics simulations [A. M. N. Niklasson et al., J. Chem. Phys. 144, 234101 (2016)] is adapted to the most recent shadow potential formulations of extended Lagrangian Born-Oppenheimer molecular dynamics, including fractional molecular-orbital occupation numbers [A. M. N. Niklasson, J. Chem. Phys. 152, 104103 (2020) and A. M. N. Niklasson, Eur. Phys. J. B 94, 164 (2021)], which enables stable simulations of sensitive complex chemical systems with unsteady charge solutions. The proposed formulation includes a preconditioned Krylov subspace approximation for the integration of the extended electronic degrees of freedom, which requires quantum response calculations for electronic states with fractional occupation numbers. For the response calculations, we introduce a graph-based canonical quantum perturbation theory that can be performed with the same natural parallelism and linear scaling complexity as the graph-based electronic structure calculations for the unperturbed ground state. The proposed techniques are particularly well-suited for semi-empirical electronic structure theory, and the methods are demonstrated using self-consistent charge density-functional tight-binding theory both for the acceleration of self-consistent field calculations and for quantum-mechanical molecular dynamics simulations. Graph-based techniques combined with the semi-empirical theory enable stable simulations of large, complex chemical systems, including tens-of-thousands of atoms.

Fast Method for Excited-State Dynamics in Complex Systems and Its Application to the Photoactivation of a Blue Light Using Flavin Photoreceptor

The excited-state dynamics of molecules embedded in complex (bio)matrices is still a challenging goal for quantum chemical models. Hybrid QM/MM models have proven to be an effective strategy, but an optimal combination of accuracy and computational cost still has to be found. Here, we present a method which combines the accuracy of a polarizable embedding QM/MM approach with the computational efficiency of an excited-state self-consistent field method. The newly implemented method is applied to the photoactivation of the blue-light-using flavin (BLUF) domain of the AppA protein. We show that the proton-coupled electron transfer (PCET) process suggested for other BLUF proteins is still valid also for AppA.

Identifying the Molecular Properties that Drive Explosive Sensitivity in a Series of Nitrate Esters

Energetic materials undergo hundreds of chemical reactions during exothermic runaway, generally beginning with the breaking of the weakest chemical bond, the "trigger linkage." Herein we report the syntheses of a series of pentaerythritol tetranitrate (PETN) derivatives in which the energetic nitrate ester groups are systematically substituted by hydroxyl groups. Because all the PETN derivatives have the same nitrate ester-based trigger linkages, quantum molecular dynamics (QMD) simulations show very similar Arrhenius kinetics for the first reactions. However, handling sensitivity testing conducted using drop weight impact indicates that sensitivity decreases precipitously as nitrate esters are replaced by hydroxyl groups. These experimental results are supported by QMD simulations that show systematic decreases in the final temperatures of the products and the energy release as the nitrate ester functional groups are removed. To better interpret these results, we derive a simple model based only on the specific enthalpy of explosion and the kinetics of trigger linkage rupture that accounts qualitatively for the decrease in sensitivity as nitrate ester groups are removed.

Understanding Explosive Sensitivity with Effective Trigger Linkage Kinetics

We present a simple linear model for ranking the drop weight impact sensitivity of organic explosives that is based explicitly on chemical kinetics. The model is parameterized to specific heats of explosion, Q, and Arrhenius kinetics for the onset of chemical reactions that are obtained from gas-phase Born-Oppenheimer molecular dynamics simulations for a chemically diverse set of 24 molecules. Reactive molecular dynamics simulations sample all possible decomposition pathways of the molecules with the appropriate probabilities to provide an effective reaction barrier. In addition, the calculations of effective trigger linkage kinetics can be accomplished without prior physical intuition of the most likely decomposition pathways. We found that the specific heat of explosion tends to reduce the effective barrier for decomposition in accordance with the Bell-Evans-Polanyi principle, which accounts naturally for the well-known correlations between explosive performance and sensitivity. Our model indicates that sensitive explosives derive their properties from a combination of weak trigger linkages that react at relatively low temperatures and large specific heats of explosion that further reduce the effective activation energy.

Accelerating Hybrid Density Functional Theory Molecular Dynamics Simulations by Seminumerical Integration, Resolution-of-the-Identity Approximation, and Graphics Processing Units

The computationally very demanding evaluation of the 4-center-2-electron (4c2e) integrals and their respective integral derivatives typically represents the major bottleneck within hybrid Kohn-Sham density functional theory molecular dynamics simulations. Building upon our previous works on seminumerical exact-exchange (sn-LinK) [Laqua, H., Thompsons, T. H., Kussmann, J., Ochsenfeld, C., J. Chem. Theory Comput. 2020, 16, 1465] and resolution-of-the-identity Coulomb (RI-J) [Kussmann, J., Laqua, H., Ochsenfeld, C., J. Chem. Theory Comput. 2021, 17, 1512], the expensive 4c2e integral evaluation can be avoided entirely, resulting in a highly efficient electronic structure theory method, allowing for fast ab initio molecular dynamics (AIMD) simulations even with large basis sets. Moreover, we propose to combine the final self-consistent field (SCF) step with the subsequent nuclear forces evaluation, providing the forces at virtually no additional cost after a converged SCF calculation, reducing the total runtime of an AIMD simulation by about another 25%. In addition, multiple independent MD trajectories can be computed concurrently on a single node, leading to a greatly increased utilization of the available hardware─especially when combined with graphics processing unit acceleration─improving the overall throughput by up to another 5 times in this way. With all of those optimizations combined, our proposed method provides nearly 3 orders of magnitude faster execution times than traditional 4c2e integral-based methods. To demonstrate the practical utility of the approach, quantum-mechanical/molecular-mechanical dynamics simulations on double-stranded DNA were performed, investigating the relative hydrogen bond strength between adenine-thymine and guanine-cytosine base pairs. In addition, this illustrative application also contains a general accuracy assessment of the introduced approximations (integration grids, resolution-of-the-identity) within AIMD simulations, serving as a protocol on how to apply these new methods to practical problems.

Polymer degradation through chemical change: a quantum-based test of inferred reactions in irradiated polydimethylsiloxane

Chemical reaction schemes are key conceptual tools for interpreting the results of experiments and simulations, but often carry implicit assumptions that remain largely unverified for complicated systems. Established schemes for chemical damage through crosslinking in irradiated silicone polymers comprised of polydimethylsiloxane (PDMS) date to the 1950's and correlate small-molecule off-gassing with specific crosslink features. In this regard, we use a somewhat reductionist model to develop a general conditional probability and correlation analysis approach that tests these types of causal connections between proposed experimental observables to reexamine this chemistry through quantum-based molecular dynamics (QMD) simulations. Analysis of the QMD simulations suggests that the established reaction schemes are qualitatively reasonable, but lack strong causal connections under a broad set of conditions that would enable making direct quantitative connections between off-gassing and crosslinking. Further assessment of the QMD data uncovers a strong (but nonideal) quantitative connection between exceptionally hard-to-measure chain scission events and the formation of silanol (Si-OH) groups. Our analysis indicates that conventional notions of radiation damage to PDMS should be further qualified and not necessarily used ad hoc. In addition, our efforts enable independent quantum-based tests that can inform confidence in assumed connections between experimental observables without the burden of fully elucidating entire reaction networks.

Exponential averaging versus umbrella sampling for computing the QM/MM free energy barrier of the initial step of the desuccinylation reaction catalyzed by sirtuin 5

The computational characterization of enzymatic reactions poses a great challenge which arises from the high dimensional and often rough potential energy surfaces commonly explored by static QM/MM methods such as adiabatic mapping (AM). The present study highlights the difficulties in estimating free energy barriers via exponential averaging over AM pathways. Based on our previous study [von der Esch et al., J. Chem. Theory Comput., 2019, 15, 6660-6667], where we analyzed the first reaction step of the desuccinylation reaction catalyzed by human sirtuin 5 (SIRT5) by means of QM/MM adiabatic mapping and machine learning, we use, here, umbrella sampling to compute the free energy profile of the initial reaction step. The computational investigations show that the initial step of the desuccinylation reaction proceeds via an S_N2-type reaction mechanism in SIRT5, suggesting that the first step of the deacylation reactions catalyzed by sirtuins is highly conserved. In addition, the direct comparison of the extrapolated free energy barrier from minimal energy paths and the computed free energy path from umbrella sampling further underlines the importance of extensive sampling.

Flexible boundary layer using exchange for embedding theories. II. QM/MM dynamics of the hydrated electron

The FlexiBLE embedding method introduced in Paper I [Z. Shen and W. J. Glover, J. Chem. Phys. 155, 224112 (2021)] is applied to explore the structure and dynamics of the aqueous solvated electron at an all-electron density functional theory Quantum Mechanics/Molecular Mechanics level. Compared to a one-electron mixed quantum/classical description, we find the dynamics of the many-electron model of the hydrated electron exhibits enhanced coupling to water OH stretch modes. Natural bond orbital analysis reveals this coupling is due to significant population of water OH σ^* orbitals, reaching 20%. Based on this, we develop a minimal frontier orbital picture of the hydrated electron involving a cavity orbital and important coupling to 4-5 coordinating OH σ^* orbitals. Implications for the interpretation of the spectroscopy of this interesting species are discussed.

Grassmann Extrapolation of Density Matrices for Born-Oppenheimer Molecular Dynamics

Born-Oppenheimer molecular dynamics (BOMD) is a powerful but expensive technique. The main bottleneck in a density functional theory BOMD calculation is the solution to the Kohn-Sham (KS) equations that requires an iterative procedure that starts from a guess for the density matrix. Converged densities from previous points in the trajectory can be used to extrapolate a new guess; however, the nonlinear constraint that an idempotent density needs to satisfy makes the direct use of standard linear extrapolation techniques not possible. In this contribution, we introduce a locally bijective map between the manifold where the density is defined and its tangent space so that linear extrapolation can be performed in a vector space while, at the same time, retaining the correct physical properties of the extrapolated density using molecular descriptors. We apply the method to real-life, multiscale, polarizable QM/MM BOMD simulations, showing that sizeable performance gains can be achieved, especially when a tighter convergence to the KS equations is required.

Quantum-Based Molecular Dynamics Simulations Using Tensor Cores

Tensor cores, along with tensor processing units, represent a new form of hardware acceleration specifically designed for deep neural network calculations in artificial intelligence applications. Tensor cores provide extraordinary computational speed and energy efficiency but with the caveat that they were designed for tensor contractions (matrix-matrix multiplications) using only low-precision floating-point operations. Despite this perceived limitation, we demonstrate how tensor cores can be applied with high efficiency to the challenging and numerically sensitive problem of quantum-based Born-Oppenheimer molecular dynamics, which requires highly accurate electronic structure optimizations and conservative force evaluations. The interatomic forces are calculated on-the-fly from an electronic structure that is obtained from a generalized deep neural network, where the computational structure naturally takes advantage of the exceptional processing power of the tensor cores and allows for high performance in excess of 100 Tflops on a single Nvidia A100 GPU. Stable molecular dynamics trajectories are generated using the framework of extended Lagrangian Born-Oppenheimer molecular dynamics, which combines computational efficiency with long-term stability, even when using approximate charge relaxations and force evaluations that are limited in accuracy by the numerically noisy conditions caused by the low-precision tensor core floating-point operations. A canonical ensemble simulation scheme is also presented, where the additional numerical noise in the calculated forces is absorbed into a Langevin-like dynamics.

Resolving the ultrafast dynamics of the anionic green fluorescent protein chromophore in water

The chromophore of the green fluorescent protein (GFP) is critical for probing environmental influences on fluorescent protein behavior. Using the aqueous system as a bridge between the unconfined vacuum system and a constricting protein scaffold, we investigate the steric and electronic effects of the environment on the photodynamical behavior of the chromophore. Specifically, we apply ab initio multiple spawning to simulate five picoseconds of nonadiabatic dynamics after photoexcitation, resolving the excited-state pathways responsible for internal conversion in the aqueous chromophore. We identify an ultrafast pathway that proceeds through a short-lived (sub-picosecond) imidazolinone-twisted (I-twisted) species and a slower (several picoseconds) channel that proceeds through a long-lived phenolate-twisted (P-twisted) intermediate. The molecule navigates the non-equilibrium energy landscape via an aborted hula-twist-like motion toward the one-bond-flip dominated conical intersection seams, as opposed to following the pure one-bond-flip paths proposed by the excited-state equilibrium picture. We interpret our simulations in the context of time-resolved fluorescence experiments, which use short- and long-time components to describe the fluorescence decay of the aqueous GFP chromophore. Our results suggest that the longer time component is caused by an energetically uphill approach to the P-twisted intersection seam rather than an excited-state barrier to reach the twisted intramolecular charge-transfer species. Irrespective of the location of the nonadiabatic population events, the twisted intersection seams are inefficient at facilitating isomerization in aqueous solution. The disordered and homogeneous nature of the aqueous solvent environment facilitates non-selective stabilization with respect to I- and P-twisted species, offering an important foundation for understanding the consequences of selective stabilization in heterogeneous and rigid protein environments.

Predicted Reaction Mechanisms, Product Speciation, Kinetics, and Detonation Properties of the Insensitive Explosive 2,6-Diamino-3,5-dinitropyrazine-1-oxide (LLM-105)

2,6-Diamino-3,5-dinitropyrazine-1-oxide (LLM-105) is a relatively new and promising insensitive high-explosive (IHE) material that remains only partially characterized. IHEs are of interest for a range of applications and from a fundamental science standpoint, as the root causes behind insensitivity are poorly understood. We adopt a multitheory approach based on reactive molecular dynamic simulations performed with density functional theory, density functional tight-binding, and reactive force fields to characterize the reaction pathways, product speciation, reaction kinetics, and detonation performance of LLM-105. We compare and contrast these predictions to 1,3,5-triamino-2,4,6-trinitrobenzene (TATB), a prototypical IHE, and 1,3,5,7-tetranitro-1,3,5,7-tetrazoctane (HMX), a more sensitive and higher performance material. The combination of different predictive models allows access to processes operative on progressively longer timescales while providing benchmarks for assessing uncertainties in the predictions. We find that the early reaction pathways of LLM-105 decomposition are extremely similar to TATB; they involve intra- and intermolecular hydrogen transfer. Additionally, the detonation performance of LLM-105 falls between that of TATB and HMX. We find agreement between predictive models for first-step reaction pathways but significant differences in final product formations. Predictions of detonation performance result in a wide range of values, and one-step kinetic parameters show the similar reaction rates at high temperatures for three out of four models considered.

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

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

Quantitative Comparison of Experimental and Computed IR-Spectra Extracted from Ab Initio Molecular Dynamics

Experimentally measured infrared spectra are often compared to their computed equivalents. However, the accordance is typically characterized by visual inspection, which is prone to subjective judgment. The primary challenge for a similarity-based analysis is that the artifacts introduced by each approach are very different and, therefore, may require preprocessing steps to determine and correct impeding irregularities. To allow for automated objective assessment, we propose a practical and comprehensive workflow involving scaling factors, a novel baseline correction scheme, and peak smoothing. The resulting spectra can then easily be compared quantitatively using similarity measures, for which we found the Pearson correlation coefficient to be the most suitable. The proposed procedure is then applied to compare the agreement of the experimental infrared spectra from the NIST Chemistry Web book with the calculated spectra using standard harmonic frequency analysis and spectra extracted from ab initio molecular dynamics simulations at different levels of theory. We conclude that the direct, quantitative comparison of calculated and measured IR spectra might become a novel, sophisticated approach to benchmark quantum-chemical methods. In the present benchmark, simulated spectra based on ab initio molecular dynamics show in general better agreement with the experiment than static calculations.

A Quantum-Based Approach to Predict Primary Radiation Damage in Polymeric Networks

Initial atomistic-level radiation damage in chemically reactive materials is thought to induce reaction cascades that can result in undesirable degradation of macroscale properties. Ensembles of quantum-based molecular dynamics (QMD) simulations can accurately predict these cascades, but extracting chemical insights from the many underlying trajectories is a labor-intensive process that can require substantial a priori intuition. We develop here a general and automated graph-based approach to extract all chemically distinct structures sampled in QMD simulations and apply our approach to predict primary radiation damage of polydimethylsiloxane (PDMS), the main constituent of silicones. A postprocessing protocol is developed to identify underlying polymer backbone structures as connected components in QMD trajectories. These backbones form a repository of radiation-damaged structures. A scheme for extracting and updating a library of isomorphically distinct structures is proposed to identify the spanning set and aid chemical interpretation of the repository. The analyses are applied to ensembles of cascade QMD simulations in which the four element types in PDMS are selectively excited in primary knock-on atom events. Our approach reveals a much higher degree of combinatorial complexity in this system than was inferred through radiolysis experiments. Probabilities are extracted for radiation-induced network changes including formation of branch points, carbon linkages, cycles, bond scissions, and carbon uptake into the Si-O siloxane backbone network. The general analysis framework presented here is readily extendable to modeling chemical degradation of other polymers and molecular materials and provides a basis for future quantum-informed multiscale modeling of radiation damage.

Mechanochemical synthesis of glycine oligomers in a virtual rotational diamond anvil cell

Mechanochemistry of glycine under compression and shear at room temperature is predicted using quantum-based molecular dynamics (QMD) and a simulation design based on rotational diamond anvil cell (RDAC) experiments. Ensembles of high throughput semiempirical density functional tight binding (DFTB) simulations are used to identify chemical trends and bounds for glycine chemistry during rapid shear under compressive loads of up to 15.6 GPa. Significant chemistry is found to occur during compressive shear above 10 GPa. Recovered products consist of small molecules such as water, structural analogs to glycine, heterocyclic molecules, large oligomers, and polypeptides including the simplest polypeptide glycylglycine at up to 4% mass fraction. The population and size of oligomers generally increases with pressure. A number of oligomeric polypeptide precursors and intermediates are also identified that consist of two or three glycine monomers linked together through C-C, C-N, and/or C-O bridges. Even larger oligomers also form that contain peptide C-N bonds and exhibit branched structures. Many of the product molecules exhibit one or more chiral centers. Our simulations demonstrate that athermal mechanical compressive shearing of glycine is a plausible prebiotic route to forming polypeptides.

The ONETEP linear-scaling density functional theory program

We present an overview of the onetep program for linear-scaling density functional theory (DFT) calculations with large basis set (plane-wave) accuracy on parallel computers. The DFT energy is computed from the density matrix, which is constructed from spatially localized orbitals we call Non-orthogonal Generalized Wannier Functions (NGWFs), expressed in terms of periodic sinc (psinc) functions. During the calculation, both the density matrix and the NGWFs are optimized with localization constraints. By taking advantage of localization, onetep is able to perform calculations including thousands of atoms with computational effort, which scales linearly with the number or atoms. The code has a large and diverse range of capabilities, explored in this paper, including different boundary conditions, various exchange-correlation functionals (with and without exact exchange), finite electronic temperature methods for metallic systems, methods for strongly correlated systems, molecular dynamics, vibrational calculations, time-dependent DFT, electronic transport, core loss spectroscopy, implicit solvation, quantum mechanical (QM)/molecular mechanical and QM-in-QM embedding, density of states calculations, distributed multipole analysis, and methods for partitioning charges and interactions between fragments. Calculations with onetep provide unique insights into large and complex systems that require an accurate atomic-level description, ranging from biomolecular to chemical, to materials, and to physical problems, as we show with a small selection of illustrative examples. onetep has always aimed to be at the cutting edge of method and software developments, and it serves as a platform for developing new methods of electronic structure simulation. We therefore conclude by describing some of the challenges and directions for its future developments and applications.

Adiabatic motion and statistical mechanics via mass-zero constrained dynamics

In recent work [Coretti et al., J. Chem. Phys., 2018, 149, 191102], a new algorithm to solve numerically the dynamics of the shell model for polarization was presented. The approach, broadly applicable to systems involving adiabatically separated dynamical variables, employs constrained molecular dynamics to strictly enforce the condition that the force on the fast degrees of freedom, modeled as having zero mass, is null at each time step. The algorithm is symplectic and fully time reversible, and results in stable and efficient propagation. In this paper we complete the discussion of the mechanics of mass-zero constrained dynamics by showing how to adapt it to problems where the fast degrees of freedom must satisfy additional conditions. This extension includes, in particular, the important case of first principles molecular dynamics. We then consider the statistical mechanics of the mass-zero constrained dynamical system demonstrating that the marginal probability sampled by the dynamics in the physical phase space recovers the form of the Born-Oppenheimer probability density. The effectiveness of the approach and the favorable scaling of the algorithm with system size are illustrated in test calculations of solid Na via orbital-free density functional dynamics.

Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials

Optically active molecular materials, such as organic conjugated polymers and biological systems, are characterized by strong coupling between electronic and vibrational degrees of freedom. Typically, simulations must go beyond the Born-Oppenheimer approximation to account for non-adiabatic coupling between excited states. Indeed, non-adiabatic dynamics is commonly associated with exciton dynamics and photophysics involving charge and energy transfer, as well as exciton dissociation and charge recombination. Understanding the photoinduced dynamics in such materials is vital to providing an accurate description of exciton formation, evolution, and decay. This interdisciplinary field has matured significantly over the past decades. Formulation of new theoretical frameworks, development of more efficient and accurate computational algorithms, and evolution of high-performance computer hardware has extended these simulations to very large molecular systems with hundreds of atoms, including numerous studies of organic semiconductors and biomolecules. In this Review, we will describe recent theoretical advances including treatment of electronic decoherence in surface-hopping methods, the role of solvent effects, trivial unavoided crossings, analysis of data based on transition densities, and efficient computational implementations of these numerical methods. We also emphasize newly developed semiclassical approaches, based on the Gaussian approximation, which retain phase and width information to account for significant decoherence and interference effects while maintaining the high efficiency of surface-hopping approaches. The above developments have been employed to successfully describe photophysics in a variety of molecular materials.

Nonadiabatic Molecular Dynamics on Graphics Processing Units: Performance and Application to Rotary Molecular Motors

Nonadiabatic molecular dynamics (NAMD) simulations of molecular systems require the efficient evaluation of excited-state properties, such as energies, gradients, and nonadiabatic coupling vectors. Here, we investigate the use of graphics processing units (GPUs) in addition to central processing units (CPUs) to efficiently calculate these properties at the time-dependent density functional theory (TDDFT) level of theory. Our implementation in the FermiONs++ program package uses the J-engine and a preselective screening procedure for the calculation of Coulomb and exchange kernels, respectively. We observe good speed-ups for small and large molecular systems (comparable to those observed in ground-state calculations) and reduced (down to sublinear) scaling behavior with respect to the system size (depending on the spatial locality of the investigated excitation). As a first illustrative application, we present efficient NAMD simulations of a series of newly designed light-driven rotary molecular motors and compare their S1 lifetimes. Although all four rotors show different S1 excitation energies, their ability to rotate upon excitation is conserved, making the series an interesting starting point for rotary molecular motors with tunable excitation energies.

Towards large scale hybrid QM/MM dynamics of complex systems with advanced point dipole polarizable embeddings

In this work, we present a general route to hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) Molecular Dynamics for complex systems using a polarizable embedding. We extend the capabilities of our hybrid framework, combining the Gaussian and Tinker/Tinker-HP packages in the context of the AMOEBA polarizable force field to treat large (bio)systems where the QM and the MM subsystems are covalently bound, adopting pseudopotentials at the boundaries between the two regions. We discuss in detail the implementation and demonstrate the global energy conservation of our QM/MM Born-Oppenheimer molecular dynamics approach using Density Functional Theory. Finally, the approach is assessed on the electronic absorption properties of a 16 500 atom complex encompassing an organic dye embedded in a DNA matrix in solution, extending the hybrid method to a time-dependent Density Functional Theory approach. The results obtained comparing different partitions between the quantum and the classical subsystems also suggest that large QM portions are not necessary if accurate polarizable force fields are used in a variational formulation of the embedding, properly including the QM/MM mutual polarization.

Synthesis of functionalized nitrogen-containing polycyclic aromatic hydrocarbons and other prebiotic compounds in impacting glycine solutions

Proteinogenic amino acids can be produced on or delivered to a planet via impacting abiotic sources and consequently were likely present before the emergence of life on Earth. However, the role that these materials played in prebiotic scenarios remains an open question, in part because little is known about the survivability and reactivity of astrophysical organic compounds upon impact with a planetary surface. To this end, we use a force-matched semi-empirical quantum simulation method to study impacts of aqueous proteinogenic amino acids at conditions reaching 48 GPa and 3000 K. Here, we probe a relatively unstudied mechanism for prebiotic synthesis where sudden heating and pressurization causes condensation of complex carbon-rich structures from mixtures of glycine, the simplest protein-forming amino acid. These carbon-containing clusters are stable on short timescales and undergo a fundamental structural transition upon expansion and cooling from predominantly sp3-bonded tetrahedral-like moieties to those that are more sp2-bonded and planar. The recovered sp2-bonded structures include large nitrogen containing polycyclic aromatic hydrocarbons (NPAHs) with a number of different functional groups and embedded bonded regions akin to oligo-peptides. A number of small organic molecules with prebiotic relevance are also predicted to form. This work presents an alternate route to gas-phase synthesis for the formation of NPAHs of high complexity and highlights the significance of both the thermodynamic path and local chemical self-assembly in forming prebiotic species during shock synthesis. Our results help determine the role of comets and other celestial bodies in both the delivery and synthesis of potentially significant life building compounds on early Earth.

Parallel replica dynamics simulations of reactions in shock compressed liquid benzene

The study of the long-term evolution of slow chemical reactions is challenging because quantum-based reactive molecular dynamics simulation times are typically limited to hundreds of picoseconds. Here, the extended Lagrangian Born-Oppenheimer molecular dynamics formalism is used in conjunction with parallel replica dynamics to obtain an accurate tool to describe the long-term chemical dynamics of shock-compressed benzene. Langevin dynamics has been employed at different temperatures to calculate the first reaction times in liquid benzene at pressures and temperatures consistent with its unreacted Hugoniot. Our coupled engine runs for times on the order of nanoseconds (one to two orders of magnitude longer than traditional techniques) and is capable of detecting reactions that are characterized by rates significantly slower than we could study before. At lower pressures and temperatures, we mainly observe Diels-Alder metastable reactions, whereas at higher pressures and temperatures we observe stable polymerization reactions.

Molecular Simulations with in-deMon2k QM/MM, a Tutorial-Review

deMon2k is a readily available program specialized in Density Functional Theory (DFT) simulations within the framework of Auxiliary DFT. This article is intended as a tutorial-review of the capabilities of the program for molecular simulations involving ground and excited electronic states. The program implements an additive QM/MM (quantum mechanics/molecular mechanics) module relying either on non-polarizable or polarizable force fields. QM/MM methodologies available in deMon2k include ground-state geometry optimizations, ground-state Born-Oppenheimer molecular dynamics simulations, Ehrenfest non-adiabatic molecular dynamics simulations, and attosecond electron dynamics. In addition several electric and magnetic properties can be computed with QM/MM. We review the framework implemented in the program, including the most recently implemented options (link atoms, implicit continuum for remote environments, metadynamics, etc.), together with six applicative examples. The applications involve (i) a reactivity study of a cyclic organic molecule in water; (ii) the establishment of free-energy profiles for nucleophilic-substitution reactions by the umbrella sampling method; (iii) the construction of two-dimensional free energy maps by metadynamics simulations; (iv) the simulation of UV-visible absorption spectra of a solvated chromophore molecule; (v) the simulation of a free energy profile for an electron transfer reaction within Marcus theory; and (vi) the simulation of fragmentation of a peptide after collision with a high-energy proton.

Effects of pressure on the structure and lattice dynamics of α-glycine: a combined experimental and theoretical study

This concomitant experimental and theoretical study provides a definitive EOS for α-glycine up to the record pressure of 50 GPa.

Inertial extended-Lagrangian scheme for solving charge equilibration models

The inertial EL/SCF method is developed to solve charge equilibration models for molecular dynamics, reducing the number of SCFs by 50–80% at each time step when compared to a conjugate gradient SCF solver and tested on diverse reactive systems.

Chemical Degradation Pathways in Siloxane Polymers Following Phenyl Excitations

We use ensembles of quantum-based molecular dynamics simulations to predict the chemical reactions that follow radiation-induced excitations of phenyl groups in a model copolymer of polydimethylsiloxane and polydiphenylsiloxane. Our simulations span a wide range of highly porous and condensed phase densities and include both wet and dry conditions. We observe that in the absence of water, excited phenyl groups tend to abstract hydrogen from other methyl or phenyl side groups to produce benzene, with the under-hydrogenated group initiating subsequent intrachain cyclization reactions. These systems also yield minor products of diphenyl moieties formed by the complete abstraction of both phenyl groups from a single polydiphenylsiloxane subunit. In contrast, we find that the presence of water promotes the formation of free benzene and silanol side groups, reduces the likelihood for intrachain cyclization reactions, and completely suppresses the formation of diphenyl species. In addition, we predict that water plays a critical role in chain scission reactions, which indicates a possible synergistic effect between environmental moisture and radiation that could promote alterations of a larger polymer network. These results could have impact in interpreting accelerated aging experiments, where polymer decomposition reactions and network rearrangements are thought to have a significant effect on the ensuing mechanical properties.

Interfacing the Core-Shell or the Drude Polarizable Force Field With Car-Parrinello Molecular Dynamics for QM/MM Simulations

We report a quantum mechanics/polarizable–molecular mechanics (QM/p–MM) potential based molecular dynamics (MD) technique where the core–shell (or the Drude) type polarizable MM force field is interfaced with the plane-wave density functional theory based QM force field which allows Car–Parrinello MD for the QM subsystem. In the QM/p-MM Lagrangian proposed here, the shell (or the Drude) MM variables are treated as extended degrees of freedom along with the Kohn–Sham (KS) orbitals describing the QM wavefunction. The shell and the KS orbital degrees of freedom are then adiabatically decoupled from the nuclear degrees of freedom. In this respect, we also present here the Nosé–Hoover Chain thermostat implementation for the dynamical subsystems. Our approach is then used to investigate the effect of MM polarization on the QM/MM results. Especially, the consequence of MM polarization on reaction free energy barriers, defect formation energy, and structural and dynamical properties are investigated. A low point charge polarizable potential (p–MZHB) for pure siliceous systems is also reported here.