Extensions and improvements of empirical force fields are discussed in view of applications to computational vibrational spectroscopy and reactive molecular dynamics simulations. Particular focus is on quantitative studies, which make contact with experiments and provide complementary information for a molecular-level understanding of processes in the gas phase and in solution. Methods range from including multipolar charge distributions to reproducing kernel Hilbert space approaches and machine learned energy functions based on neural networks.

Atomistic simulations have become a powerful way to provide molecular-level insights into gas- and condensed-phase properties of physico-chemical systems and to interpret and complement experiments. Much of the progress is directly related to the dramatically increased capabilities of modern computer infrastructure. Nevertheless, the field is still far from solving general problems by running either sufficiently long (quasi-) classical trajectories or large-scale quantum dynamics simulations in order to characterize a particular system. Even for an atom + diatom system, the comprehensive calculation of all state-to-state cross sections is a formidable task, and brute-force sampling is a serious computational undertaking.1,2

For studying the dynamics of complex systems, means to determine their total energies and associated forces are required. Empirical force fields were initially developed to characterize the chemical structure and dynamics of macromolecules, including peptides and proteins.3–10 Their primary applications concerned sampling conformations of large molecules, such as proteins in the gas phase11 and in solution.12 The typical mathematical form of an empirical force field contains bonded and nonbonded terms. The bonded part is

V b o n d = K b ( r r e ) 2 , V a n g l e = K θ ( θ θ e ) 2 , V d i h e = K ϕ ( 1 + cos ( n ϕ δ ) )
(1)

and is geared toward conformational sampling, spectroscopy, or thermodynamics and not, e.g., toward describing chemical reactions for which bonds need to be broken and formed. In Eq. (1), the constants K are the force constants associated with the particular type of interaction, re and θe are equilibrium values, n is the periodicity of the dihedral, and δ is the phase that determines the location of the maximum. The sums run over all respective terms. Nonbonded interactions include electrostatic and van der Waals terms

V e l s t a t = 1 4 π ϵ 0 q i q j r i j , V v d W = ϵ i j R m i n , i j r i j 12 2 R m i n , i j r i j 6 ,
(2)

where the sums run over all nonbonded atom pairs. qi and qj are the partial charges of the atoms i and j, and ϵ0 is the vacuum dielectric constant. For the van der Waals terms, the potential energy is expressed as a Lennard-Jones potential with the well depth ϵ i j = ϵ i ϵ j and range Rmin,ij = (Rmin,i + Rmin,j)/2 at the Lennard-Jones minimum, according to standard mixing rules.13 This interaction captures long range dispersion (∝ − R−6) and exchange repulsion (∝ R−12) where the power of the latter is chosen for convenience. The combination of Eqs. (1) and (2) constitutes a minimal model for an empirical force field.

One of the desirable features of an empirical force field is that it is a computationally efficient model that can be improved in various ways based on physical principles such as to allow direct comparison with experiment and eventually provide a basis for interpretation. Forces can be evaluated analytically, which is an advantage for molecular dynamics (MD) simulations. For example, harmonic bonds can be replaced by Morse oscillators, which are either useful to investigate chemical reactions or accurate vibrational spectroscopy. These two aspects are of particular relevance to the present work.

Using point charges (PCs) for the electrostatic interactions can be considered as the zero-order expansion of a more comprehensive multipolar expansion, which considerably improves the description of the charge distribution around molecular building blocks.14–19 In addition, the standard form of an empirical force field can be complemented by including further terms such as polarization.20 The extensibility and specific ways to improve the quality of an empirical force field make it an ideal starting point for using them in increasingly realistic simulations of molecular materials at an atomistic level.

The breaking and formation of chemical bonds as a function of time lie at the heart of chemistry. However, the most obvious and direct route, evaluating the electronic structure along a molecular dynamics (MD) trajectory or time-dependent quantum dynamics simulation, is often impractical for several reasons. For one, evaluating the total energy and forces by solving the electronic Schrödinger equation is computationally prohibitive. This precludes using high-level quantum chemistry methods and/or sufficiently large basis sets for running a statistically significant number of independent trajectories. For converged observables, averages over initial conditions need to be computed. In addition, some electronic structure calculations can be affected by systematic error, such as basis set superposition errors, and convergence of the Hartree–Fock wavefunction to the desired state can be challenging. This can occur for systems containing metal atoms, higher excited states, or systems with a challenging electronic structure.

In the present perspective, several improvements to empirical energy functions and their application that have been pursued in the more recent past will be highlighted. The applications focus on vibrational spectroscopy, reactive dynamics, and thermodynamic properties, which provide a wide range of challenging problems for understanding and interpreting physico-chemical experiments at a molecular level.

Both bonded and non-bonded parts of an empirical force field can be systematically improved. An instructive example is recent efforts to correctly describe condensed-phase and spectroscopic properties of liquid water. Building on the realization that multipolar interactions are required,17,21–26 the TL4P water model was developed.25 TL4P is a rigid water model with the atoms fixed to the experimental liquid phase geometry. The electrostatic interactions are described by three-point charges (one on each hydrogen atom and one on a massless point M that is displaced from the oxygen atom) and a Gaussian inducible dipole. The van der Waals interactions are represented as a Buckingham potential acting between oxygen atoms. With this model, parameterized to reference density functional theory calculations, a wide range of bulk properties could be described reliably compared with the experiment (see Table 7 of Ref. 25).

When using TL4P for computing the low-frequency part of the infrared and Raman spectra of liquid water, it was found that the shapes of the spectral features did not compare well with the experiment.26 Based on earlier work, it was decided to augment the TL4P water model with some amount of charge transfer along the O–HO intramolecular hydrogen bond and to include anisotropic polarizability.27 The parameters were again determined from reference DFT calculations, and both bulk properties and the optical spectra were computed. With these modifications, the bulk properties were in equally good agreement with experiment as for the original TL4P model,25 whereas the relative intensities and peak positions for both, infrared and Raman, spectra compared well with the experiment. In this way, a physically motivated and computationally efficient model for water was developed by refinement, largely guided by comparison with experiment.

Following very different strategies, other successful force fields for water have been developed in the recent past. One of them is the AMOEBA model, which builds explicitly on atomic multipoles.22,23 Another strategy is followed in the E3B model, which is based on adding explicit three-body terms, akin to a many-body expansion.28 Similarly, the HBB (Huang, Braams, Bowman) force field also uses a many-body expansion.29 As a more recently developed force field, the RexPoN force field is entirely based on reference electronic structure calculations without empirical data.30 Finally, the most comprehensive water force field in various phases is probably the MB-Pol model, which also builds on multipolar interactions and many-body polarization.31 

Multipolar interactions have been considered early on in molecular recognition.32 Different from the spherically symmetric field around a point charge, atomic multipoles can capture anisotropic interactions. A typical example is carbon monoxide, which cannot be realistically modeled with two atom-centered point charges located at each of the two atoms [Eq. (2)] because the total charge (Q = 0) and the total molecular dipole μ = 0.048 ea0 lead to two partial charges close to zero but with opposite sign. In order to describe the substantial molecular quadrupole moment33–35 ( Θ = 1.58 e a 0 2 ), either a third interaction site midway the two atoms can be included36 or the two atoms are described by a distributed multipole expansion.14,21,37,38

More generally, the electrostatic potential (ESP) around a molecule can be represented as a superposition of point charges (PC) and higher multipole (MTP) moments. Capturing strongly anisotropic features, e.g., lone pairs, hydrogen bonding, π-electron density, or sigma holes (halogens),39–41 requires schemes beyond point charge representations. The ESP, Φ r , is related to the charge distribution ρ(r) through

4 π ε 0 Φ ( r ) = d r ρ ( r ) r r = l = 0 m = l l Q l m r l + 1 4 π 2 l + 1 Y l m ( θ , ϕ ) ,
(3)

where r and r′ are spatial variables. For the second equation, 1/|rr′| was expanded in powers of r′/r < 142 to represent the ESP as a sum over spherical harmonics Ylm(θ, ϕ) from which the spherical MTP moments Qlm are defined by

Q l m = d r ρ ( r ) ( r ) l 4 π 2 l + 1 Y l m * ( θ , ϕ ) .
(4)

These can be determined from the density by spatial integration and provide a compact atom-based (“distributed”)14,43 representation of the ESP around a molecule. The interaction energy between two sites i and j, each represented as a distributed multipolar expansion, is then U t u i j ( q ) = Q t i Q u j T t u i j ( q ) , where Q t i and Q u j are the MTP moments of order t and u, respectively, T t u i j ( q ) describes the geometrical dependence of the interactions, and q are the direction cosines for the relative orientations of the multipoles on sites i and j.14,43 As an example, the term with t = 0 and u = 2 is for the charge–quadrupole interaction.

Alternatively, multipoles of a given order can be represented by distributed point charges as is done in the distributed charge model (DCM).18,19 Such a model replaces the evaluation of multipole–multipole interactions at the expense of using a larger number of terms. Reducing the number of interaction sites can be accomplished by finding the smallest number of point charges to accurately represent the ESP using differential evolution.19 

The design of reactive force fields has followed different routes. Early work by Penney44 and Pauling45 related the bond order, bond length, and bond strength. The finding that the bond order and bond length follow a near linear relation45 and a log–log plot of dissociation energies vs bond order is also almost linear was used to develop “bond energy bond order” (BEBO) potentials.46 One of the essential assumptions underlying this approach is that—at least for hydrogen-atom transfer reactions—the sum of the bond orders of the broken bond, n1, and of the newly formed bond, n2, fulfills n1 + n2 = 1. In other words, “at all stages of the reaction, the formation of the second bond must be “paying for” the breaking of the first bond.”46 

An extension of the concept underlying BEBO is ReaxFF.47 Here, nonbonded terms are included from the beginning, and the dissociation and reaction curves are derived from fitting to electronic structure calculations.47,48 Central to ReaxFF is that the bond order is related to the distance between two atoms. As an example, for two carbon atoms, there can be anything between “no bond” (bond order = 0) to a triple bond. In ReaxFF, the total energy consists of a sum of individual terms such as Ebond, a correction for over-coordination, a penalty term Eover for under-coordinated atoms reflecting resonance energies between π−electrons, a conjugation energy, and others. Contrary to empirical force fields,7–9 ReaxFF does not build on the concept of atom types. While the functional form of ReaxFF is universal, its application to concrete problems always involves more or less extensive fitting to reference calculations.49 

Chemical reactivity can be encapsulated in a natural fashion within the framework of valence bond theory. For the ionic bond cleavage AB → A + B+, three chemically relevant states are considered within empirical valence bond (EVB)50 theory: ψ1 = AB, ψ2 = AB+, and ψ3 = A+B. If A is more electronegative than B (i.e., formation B is unlikely in the presence of A), the resonance structure ψ3 is largely irrelevant, and the process can be described by ψ1 and ψ2 alone. For a more general chemical reaction, one must choose, “… on the basis of experience and intuition, a set of bonding arrangements or “resonance structures”[…]” as was noted in the original EVB publication.50 The choice of the states in EVB is not always obvious a priori but may also need to be based, e.g., on the requirement to best reproduce the reference electronic structure calculations.51 

The total n state EVB Hamiltonian Hnn contains diagonal terms (force fields for each of the n states) and the coupling between them (Hij). Representation and parameterization of the off-diagonal terms has been a source of considerable discussion;52–54 in particular, the assumption that upon transfer of the reaction from the gas phase to the solution phase, these elements do not change significantly.55 Applications of EVB include enzymatic reactions (for which it was originally developed56), proton transfer processes, and the autodissociation of water.57 More recently, this has been extended to other types of reactions, including bimolecular reactions58,59 or the association and dissociation of CH3 from diamond surfaces.60 Furthermore, several extensions have been suggested to the original EVB method, allowing its application to a wider class of problems.61–63 

Time resolved experiments have contributed considerably to our understanding of chemical reactivity over the past three decades. Short laser pulses (“femtochemistry”) allow us to follow chemical transformations on time scales relevant to the actual chemical step (bond breaking or bond formation). Typical examples are the time resolved studies of ligand (re)binding in myoglobin (Mb)64–67 or vibrationally induced reactivity.68–70 This prompted the development of time-based reactive MD71,72 whereby a chemical reaction is followed in time, which is the natural progression coordinate in such an experiment. For the example of two states (“R” and “P”), the two PESs are mixed according to V e ff ( x ) = f ( t ) V R ( x ) + ( 1 f ( t ) ) V P ( x ) , where f(t) is a sigmoid switching function, which changes from 1 to 0 between t = −ts/2 and t = ts/2. At the beginning of the mixing, the system is fully in the R-state [f(t = −ts/2) = 1], while at the end, it is fully in the P-state [f(t > ts/2) = 1]. Here, the only free variable is the switching time ts. The algorithm of time-based ARMD is schematically shown for a collinear atom transfer reaction in Fig. 1.

FIG. 1.

Schematic illustration of the ARMD method for a collinear reaction, where atom B is transferred from donor atom A to acceptor atom C. During crossing, the surfaces are switched in time, and the Morse bond AB is replaced by van der Waals (vdW) interactions and vice versa.

FIG. 1.

Schematic illustration of the ARMD method for a collinear reaction, where atom B is transferred from donor atom A to acceptor atom C. During crossing, the surfaces are switched in time, and the Morse bond AB is replaced by van der Waals (vdW) interactions and vice versa.

Close modal

During the mixing, the system is propagated under the influence of a time-dependent Hamiltonian, which does not strictly conserve total energy.73,74 The magnitude of this effect can be shown to scale with 1/m, where m is the reduced mass involved in the reaction,74 which is inconsequential for heavy systems such as NO rebinding to Mb but affects the dynamics in the product channel for proton transfer processes. This led to the development of multi state adiabatic reactive MD (MS-ARMD), which mixes the PESs with energy-dependent weights and strictly conserves energy.74 

In MS-ARMD, the PESs are mixed in energy-space according to

V M S A R M D ( x ) = i = 1 n w i ( x ) V i ( x ) .
(5)

The weights (see Fig. 2) w i ( x ) are obtained by renormalizing the raw weights w i , 0 ( x ) ,

w i ( x ) = w i , 0 ( x ) i = 1 n w i , 0 ( x )    w h e r e    w i , 0 ( x ) = exp ( V i ( x ) V m i n ( x ) ) Δ V ,
(6)

where V m i n ( x ) is the minimal energy for a given configuration x and ΔV is a parameter. The raw weights [Eq. (6)] depend exponentially on the energy difference between surface i and the minimum energy surface V m i n ( x ) over a characteristic energy scale ΔV (switching parameter). Only those surfaces will have significant weights, whose energy is within a few times of ΔV from the lowest-energy surface V m i n ( x ) for instantaneous configuration x . The performance of MS-ARMD is demonstrated for crossings of 1D and 2D surfaces in Fig. 2.

FIG. 2.

Schematic for the MS-ARMD method applied in one (left panel) and two (right panel) dimensions to 3 and 2 surfaces (V1,2,3), respectively. The effective surface (VMS−ARMD) is always close to the lowest-energy surface for a particular configuration x, except for regions where other surfaces are within a few times ΔV (here = 0.5 units) in energy. In these regions, the algorithm switches smoothly between the PESs by varying the weights ( w 1 , 2 , 3 ; lower left panel).

FIG. 2.

Schematic for the MS-ARMD method applied in one (left panel) and two (right panel) dimensions to 3 and 2 surfaces (V1,2,3), respectively. The effective surface (VMS−ARMD) is always close to the lowest-energy surface for a particular configuration x, except for regions where other surfaces are within a few times ΔV (here = 0.5 units) in energy. In these regions, the algorithm switches smoothly between the PESs by varying the weights ( w 1 , 2 , 3 ; lower left panel).

Close modal

Finally, ARMD with energy-dependent weights mixes the different PESs by using Gaussian and polynomial functions (GAPOs) in the neighborhood of crossing points between the states. The parameters of these GAPOs need to be determined through fitting to the reference data [e.g., from a calculation along the intrinsic reaction coordinate (IRC)]. This step is the most demanding part in MS-ARMD.74 A smooth global surface is obtained everywhere, even for energies where more than two surfaces approach one another. Because the mixed PES V M S A R M D ( x ) depends on the energies of the different states through the weights w i , which, in turn, are analytical functions of the coordinates x , the derivatives can be readily determined, which leads to energy conservation in MD simulations.

A recent extension75 of MS-ARMD is its combination with VALBOND, a force field that allows us to describe the geometries and dynamics of metal complexes.76–78 Here, the formulation is reminiscent of EVB whereby the diagonal terms are VALBOND descriptions of the states involved, and the off-diagonal elements describe the orbital overlap. Furthermore, MS-ARMD can also be combined with molecular mechanics with proton transfer (MMPT)79 to follow proton transfer in the gas- and condensed phase.80–83 

Vibrational spectroscopy is a sensitive tool for quantitatively characterizing the structure of and the structural dynamics around a solute embedded in an environment.84 Such studies can also provide important insight into the dynamical coupling between the solute and the solvent. One example is the possibility to infer ligand binding strengths from infrared spectroscopic measurements.85 A quantitative study for benzonitrile in lysozyme has been carried out recently using atomistic simulations with non-conventional force fields.86 Another interesting future experimental application concerns the elucidation of structural dynamics and biological function, which, however, requires sufficient sensitivity and selectivity of the spectroscopic probe.87 In addition, in the realm of materials design, infrared spectroscopy can be potentially useful to characterize structural properties of strongly aggregating species such as ionic liquids.88 For all these applications, sensitive and selective vibrational probes are required.

Refined multipolar force fields have shown to allow quantitative comparisons with experiments and their interpretation.21,38 One of the noticeable examples is the infrared spectrum of photodissociated carbon monoxide (CO) in myoglobin (Mb).89 The strong (43 MV/cm)90 inhomogeneous electric field in the heme pocket leads to characteristic shifting and splitting of the spectral lines due to the Stark effect. Several attempts were made91–93 to correctly interpret the experimental infrared spectrum using computational methods. Although some of them were capable of correctly modeling the width of the spectrum, they usually were unable to find the characteristic splitting of the CO spectrum (i.e., ≈10 cm−1). A first successful attempt used a fluctuating point-charge model21 based on an earlier three-point model for CO.36 This was later generalized to a rigorous fluctuating MTP model, which reproduced most features of the spectrum known from experiments.38 In particular, the splitting, width, and relative intensities of the computed spectrum agreed favourably with the experimentally known properties. Based on this, it was then also possible to assign the two spectroscopic signatures to distinct conformational substates. These agreed with the results from more indirect experiments (based on mutations in the active site)94 or more difficult-to-converge mixed QM/MM simulations.95 

Understanding structural, spectroscopic, and dynamical properties of water in its different phases provides a rich variety for developing and improving special-purpose, tailored, and physics-based force fields. Among the most successful for spectroscopic applications are MB-pol31 and WHBB.96 Other similarly successful models include the explicit 3-body (E3B) parametrization,28 the (i) AMOEBA family of force fields,22,23 the TLnP family,25 and a recent modification of it.27 

The ManyBody-polarizable (MB-pol) model is based on two- and three-body water interactions calculated at the CCSD(T) level of theory. It is entirely developed from first-principles.31,97 MB-pol explicitly treats the one-body (monomer distortion energy) term and the short-ranged two- and three-body terms akin to a polarizable potential supplemented by short-range two- and three-body terms that effectively represent quantum–mechanical interactions arising from the overlap of the monomer electron densities. Specifically, at all separations, the total MB-pol two-body term includes (damped) dispersion forces derived from ab initio computed asymptotic expansions of the dispersion energy along with electrostatic contributions due to the interactions between the molecular permanent and induced moments. Similarly, the MB-pol three-body term includes a three-body polarization term at all separations, which is supplemented by a short-range fourth-degree permutationally invariant polynomial that effectively corrects for the deficiencies of a purely classical representation of the three-body interactions in regions where the electron densities of the three monomers overlap. This short-range three-body contribution is smoothly switched off once the oxygen–oxygen separation between any water molecule and the other two water molecules of a trimer reaches a value of 4.5 Å. In MB-pol, all induced interactions are described through many-body polarization. MB-pol thus contains many-body effects at all monomer separations, as well as at all orders, in an explicit way up to the third order and in a mean-field fashion at all higher orders. Using MD simulations, the computed infrared and Raman spectra for liquid water can be realistically modeled. Finally, molecular dynamics (MB-MD) simulations carried out with MB-pol describe the infrared (IR) and Raman spectra of liquid water well compared with experiments.98 

The WHBB model96 is the sum of a monomer and two-body and three-body, full-dimensional potentials, with the option of including 4 and higher–body interactions from other polarizable potentials. The monomer PES is spectroscopically accurate and the two- and three-body potentials are mathematical fits to tens of thousand of ab initio electronic energies [CCSD(T)/aug-cc-pVTZ (aVTZ) and MP2/aVTZ energies] using permutationally invariant polynomials.99 The WHBB model has been applied to the infrared spectroscopy of liquid water and showed to perform well compared with experiment.100 

Multipolar force fields have recently been used for a number of 1D- and 2D-spectroscopy studies including N-methylacetamide,101 fluoro-acetonitrile,102 cyanide,103 insulin in water,104 and azide in the gas phase and in water.105 As an example for using such simulations in a concrete context, the 1d-infrared spectroscopy of wild type insulin monomer (WT-MO, chains A and B, see Fig. 3) and dimer (WT-DI, chains A to D) in solution is discussed here. For this, the x-ray crystal structure of the WT insulin dimer resolved at 1.5 Å [protein data bank (PDB106,107) code: 4INS]108 was solvated in a cubic box (753 Å3) of TIP3P109 water molecules. For all molecular dynamics (MD) simulations, the CHARMM110 package together with the CHARMM36111 force fields was used. The systems were minimized and equilibrated for 20 ps in the NVT ensemble. Production runs, 1 ns in length, were carried out in the NPT ensemble, with coordinates saved every 10 fs for subsequent analysis. A velocity Verlet112 integrator and Nosé–Hoover thermostat113,114 were employed in the NVT simulations. For the NPT simulations, an Andersen and Nosé–Hoover constant pressure and temperature algorithm114–116 was used together with a leapfrog integrator.117 

FIG. 3.

(a) Structure of the WT insulin dimer with the CO labels considered specifically. (b) Comparison between the maxima of the frequency distribution for WT-MO and WT-DI. Residues of chain A/B and C/D are shown in green and red, respectively.

FIG. 3.

(a) Structure of the WT insulin dimer with the CO labels considered specifically. (b) Comparison between the maxima of the frequency distribution for WT-MO and WT-DI. Residues of chain A/B and C/D are shown in green and red, respectively.

Close modal

Time-ordered snapshots were extracted from the MD simulation every 10 fs, and then, the anharmonic frequencies were computed. For this purpose, one-dimensional potential energy curves were calculated along the –CO normal modes for each snapshot. The anharmonic transition frequencies ( v = 0 → 1) were then calculated by solving the 1D time-independent Schrödinger equation to obtain the frequency trajectory. Figure 3(b) compares the maximum of the frequency distribution for 105 snapshots between the monomer and the dimer for several residues (6, 19, and 20 from chains A/C and 6 and 19–26 from chains B/D) for the WT protein. Residues 24–26 of the B/D chains are involved in hydrogen bonding at the dimerization interface. As is shown in Fig. 3(b), the maximum frequencies for residues 24 and 26 of chain B/D are both blue shifted when going from the monomer to the dimer, while for residue 25B, a small red shift is observed and clearly blue shifted for residue 25D. This is also consistent with previous results,104 which showed that the IR spectra of the symmetry-related pairs of chain B/D do not superimpose. This also suggests that the structural dynamics and spectroscopy of the CO probes involved in the protein–protein contact (24–26 B/D) represent characteristic changes, and this feature can be used to distinguish the monomeric from dimeric insulin in solution.

As a second example, the infrared spectroscopy of N-methyl acetamide (NMA) is considered. For this case, the particular focus is on the influence of how the instantaneous frequency of the oscillator in question is determined as there are multiple possibilities. Contrary to the case of CN in solution for which the oscillator is one-dimensional,103 obtaining an instantaneous frequency for a given molecular configuration for polyatomic molecules can follow different protocols. For one, a harmonic vibrational analysis (instantaneous normal mode)118 can be carried out. If anharmonicities should be correctly captured, scans along the potential energy surface can be carried out to map out an instantaneous energy landscape from which a transition frequency can be calculated by solving the nuclear Schrödinger equation.105 Such a scan can be carried out in different ways, and two of them together with their influence on the results are explored in the following.

Simulations of N-methyl acetamide (NMA) were carried out in which the entire NMA molecule was treated with MTP electrostatics.25,26,101,119 For the bonded terms, a reproducing kernel Hilbert space (RKHS-) based representation of the PES has been constructed for the amide-I stretch (CONH) fragment. For convenience and illustration, this was based on density functional theory calculations with the PBE (Perdew, Burke, Ernzerhof) functional120,121 for correlation and exchange using the cc-pVTZ basis set for the reference calculations although the approach can also be applied to the reference data calculated at much higher levels of theory. In the simulations, the electrostatic interactions are treated using Particle–Mesh Ewald122 (PME) with a grid size spacing of 1 Å, characteristic reciprocal length κ = 0.32 Å−1, and interpolation order 4. All bonds involving hydrogen atoms were constrained using SHAKE123 except for the NH bond in the amide-I mode.

To construct the RKHS PES (for details on RKHS, see Refs. 124 and 125), the NMA structure was optimized first, and then, the normal modes are calculated by fixing all hydrogen atoms of the methyl groups. Next, gas phase energies and forces were calculated using the Gaussian09 software126 for 2400 structures sampled randomly along the 12 normal modes. In order to describe the PES for the NMA fragment [X(CONH)Y where “X” and “Y” are the two methyl groups as united atoms], a multidimensional reproducing kernel based function is considered, which is the sum of 15 four-body interactions (6C4 = 15, where nCk is the binomial coefficient). Each of the four-body interactions is represented by a six-dimensional tensor product kernel for six radial dimensions (the interatomic distances in tetra-atomic systems). A reciprocal power decay kernel of type (n = 3, m = 0) was used for all the bond distances.124,125 The kernel coefficients were calculated using least square fitting. Finally, the RKHS PES has been constructed for the amide-I stretch mode. This RKHS PES describes all relevant bonded terms for the CONH motif and also the bonded interaction in H3C–C and H3C–N. The harmonic frequency calculated for the amide-I stretch from the RKHS PES compares well with the one obtained from the DFT calculations (1696 cm−1 vs 1698 cm−1), which serves as a validation of the present approach.

To simulate the dynamics of NMA, the RKHS-based representation for the bonded terms is implemented in CHARMM and supplemented with conventional force field parameters to describe the interactions between the methyl hydrogen atoms and all other atoms. The harmonic frequency of the amide-I stretch for NMA is at 1714 cm−1, slightly higher than that computed using the stand-alone RKHS PES (1696 cm−1) from above. This is due to coupling between the CONH motif and the rest of the NMA molecule. Next, MD simulations are run for the solvated system and the one-dimensional potential energy is scanned along two different directions, the CO and the CONH directions, see Fig. 5, to investigate the differences depending on the approach chosen.

The anharmonic frequencies for the fundamental and overtones of a particular mode and corresponding eigenfunctions can be calculated by solving the 1D Schrödinger equation using a discrete variable representation (DVR) approach.127 For the amide-I mode in NMA—which can also be applied to the protein backbone, see above—two different schemes are followed: (a) scanning along the CO stretch and (b) displacing the CONH group along the amide-I normal mode (see Fig. 4). The potential energies are again represented by a 1D-reproducing kernel, and the Schrödinger equation is then solved using a DVR approach for − < r < , where r is the displacement from the equilibrium geometry. The anharmonic stretching frequencies in the gas phase are computed as 1671 cm−1 and 1695 cm−1 using the potential obtained following schemes (a) and (b), respectively. The reduced mass is 1 amu for both cases since unnormalized displacement coordinates are used to distort the molecule along the normal modes.

FIG. 4.

Displacement vector for two different scanning to construct a 1D potential. (a) Distortion along the CO stretching and (b) movement along the CONH stretching.

FIG. 4.

Displacement vector for two different scanning to construct a 1D potential. (a) Distortion along the CO stretching and (b) movement along the CONH stretching.

Close modal

Using the RKHS-based PES described above, a 4 ns trajectory of solvated NMA (303 Å3 TIP3P water109) in the NVE ensemble was run, following 20 ps of heating and 1.2 ns of equilibration at 300 K. The trajectory was saved every 5 fs, and the anharmonic frequencies were determined from scanning the RKHS-potential along the CO and the CONH directions (see Fig. 4) for 8 × 105 snapshots using MTP electrostatics. The MTPs mainly account for the NMA-water electrostatic interactions, while the gas phase RKHS potential accurately represents the amide-I mode within NMA. Figure 5 shows the decay of the frequency fluctuation correlation function (FFCF) and frequency distribution (see inset in Fig. 5) for NMA in water with anharmonic frequency calculated via scanning along the two directions. The average frequencies for scanning along the CO (a) and CONH mode (b) of NMA in water are at 1657 cm−1 and 1671 cm−1, respectively. A tri-exponential fit to the frequency fluctuation correlation functions yields three timescales, which are close to one another from approaches (a) and (b), and compare quite well with experiment128,129 and previous simulations,101 see Table I.

FIG. 5.

Frequency fluctuation correlation functions (FFCF) of NMA in water with the anharmonic frequency calculated via scanning along the CONH mode (red) and CO stretching mode (blue) using multipolar force fields and gas phase RKHS potential of NMA described above. The black lines refer to the tri-exponential fit to the FFCF. Inset shows the corresponding frequency distributions.

FIG. 5.

Frequency fluctuation correlation functions (FFCF) of NMA in water with the anharmonic frequency calculated via scanning along the CONH mode (red) and CO stretching mode (blue) using multipolar force fields and gas phase RKHS potential of NMA described above. The black lines refer to the tri-exponential fit to the FFCF. Inset shows the corresponding frequency distributions.

Close modal
TABLE I.

The FFCF decay times and their amplitudes for scanning along the CO and CONH normal modes for solvated NMA obtained from a tri-exponential fit, see Fig. 5. The frequency calculations used the RKHS PES and MTP electrostatics. For comparison, vibrational decay times from experiment128 and simulations129 are also given.

Method a1 (ps−2) t1 (ps) a2 (ps−2) t2 (ps) a3 (ps−2) t3(ps)
CO scan  0.71  0.021  0.302  0.208  0.103  1.0 
CONH scan  0.70  0.020  0.25  0.20  0.106  0.81 
Sim.129     0.06        0.66 
Exp.128     (0.05–0.1)        1.6 
Exp.129     0.01        1.0 
Method a1 (ps−2) t1 (ps) a2 (ps−2) t2 (ps) a3 (ps−2) t3(ps)
CO scan  0.71  0.021  0.302  0.208  0.103  1.0 
CONH scan  0.70  0.020  0.25  0.20  0.106  0.81 
Sim.129     0.06        0.66 
Exp.128     (0.05–0.1)        1.6 
Exp.129     0.01        1.0 

As a final example for a spectroscopic application, a recent model for azide in the gas phase and in solution is discussed. The bonded terms were again based on RKHS, which exactly reproduces the energies from the reference ab initio calculations. For N 3 , a three-dimensional PES based in multireference CI (MRCI) calculations was represented as an RKHS, and quantum bound states were determined.105 The agreement with gas phase experiments was found to be good, typically within ∼10 cm−1. Simulations in solution provided the correct blue shift, consistent with experiment. This opens now the possibility to use the new parameterization of N 3 for simulations of proteins in solution and to determine the site-specific dynamics depending on the position of the label.

It is known that the charge distribution ρ(x) around a molecule and therefore the charges or multipoles used to represent ρ(x) depend on the geometry. Fluctuating point charge models based on the charge equilibration scheme have been developed to capture this.130,131 However, they have been mainly applied to simulating the structure and diffusivity of water, ions in water, or amides in water and not for infrared spectroscopy. Furthermore, generalizations to simulating larger molecules are difficult, and for multipoles, no such extensions are available although their conformational dependence has been investigated.132 For CO in myoglobin38 and CN in water,103 such conformationally dependent MTPs have demonstrated to perform very well. Future developments for computational infrared spectroscopy in solution will benefit from a rigorous treatment to capture such effects. A first step in this direction has been made by machine learning the point charges for an ensemble of structures of the same molecule, e.g., within the PhysNet neural network approach.133 

The investigation of chemical reactions using atomistic simulations dates back to the 1960s when the first studies of the H + H2 reactions were carried out.134 One unexpected finding was that such quasi-classical trajectory (QCT) simulations agree quite well with quantum simulations135 despite the fact that the H + H2 system is particularly susceptible to quantum effects including zero point energy and tunneling. In the following, QCT simulations became an established way to investigate small- and large-molecule reactive processes.

The field of small molecule reactions was heavily driven by the advent of refined experimental techniques, such as molecular beam studies or photodissociation spectroscopy, which provide energy-, state-, and conformationally-resolved information about reactive molecular collisions. Initially, empirical energy functions, such as the London–Eyring–Polanyi–Sato (LEPS) surface,136–138 were used to carry out QCT simulations. Once electronic structure calculations became sufficiently accurate and the configurational space could be covered by computing energies for many geometries, such empirical PESs were largely superseded. The primary problem then shifted to representing the computed points such that the total potential energy can be evaluated with comparable accuracy as the underlying quantum chemical calculations.

For triatomic molecules, a convenient representation was the one that used an expansion into Legendre polynomials Pλ(cos θ) and corresponding radial strength functions Vλ(R, r), i.e., V ( R , r , θ ) = λ = 0 λ m a x V λ ( R , r ) P λ ( cos θ ) . Such expansions have been particularly useful for studying van der Waals complexes.139–142 Such representations have been used either together with explicit fits to experimental data or by fitting them to reference electronic structure calculations. An alternative is to separate the total interaction energy into long- and short-range parts and to represent them separately. Often, the long range part in such complexes can be described very accurately by resorting to classical electrostatics based on atomic and molecular multipoles and polarizabilities and hyperpolarizabilities.143,144

An alternative is an expansion into permutationally invariant polynomials (PIPs),99 which has, for example, been applied to N4 for studying reactive collisions for N2 + N2 → N2 + 2N and N2 + N2 → 4N.145 PIPs employ a basis of exponential or Morse-type functions to expand the PES and fit products of such basis functions to the reference data. The permutational symmetry of chemically identical atoms is preserved in such an approach. Similarly, a study of N 4 + using an RKHS representation also incorporated permutational invariance explicitly.146 

Another method is a representation based on a Reproducing Kernel Hilbert Space (RKHS).124,125,147 Such an approach exactly reproduces the input data, which are usually results from electronic structure calculations on a predefined grid. One of the advantageous features of an RKHS is that the physical long range decay for a radial coordinate → can be encoded in the functional form of the kernel. Such RKHS-based representations have been used for entire PESs148 or in a QM/MM-type approach to treat part of an extended system with higher accuracy.149–151 

Recently, RKHS-based reactive PESs have been successfully used for a range of atom + diatom reactions. For these systems, the focus was primarily on thermal reaction rates, final state distributions, and vibrational relaxation times, which are properties of particular interest to hypersonics and atmospheric re-entry.152–154 Reference data for the systems of interest ([NOO], i.e., N + O2 and O + NO; [NNO], and [CNO]) can be determined at high levels of theory (MRCI) with good-quality basis sets (aug-cc-pVTZ or larger).148,155–157 A typical PES for one electronic state is based on 104 energies, and the RKHS-interpolated surfaces reproduce the reference data to within much better than chemical accuracy, typically within a few cm−1. Extensive QCT simulations on these PESs now provide a comprehensive and consistent set of thermal rates and vibrational relaxation times for the most relevant systems involved in hypersonics.158 

The application of such methods to bimolecular reactive collisions involving larger gas phase systems, in particular CH4, has recently been reviewed.159,160 The methods can also be applied to SN2 reactions in the gas phase161 to explore alternative reaction mechanisms and in solution162 to quantify the differences between reactions in the gas- and condensed phase.

For larger molecules in the gas phase, MS-ARMD has also been used together with established empirical force fields. This has been applied to the study of specific reaction channels such as proton transfer within H2SO4 or photodissociation of H2SO4 → H2O + SO3163 and other atmospherically relevant molecules following vibrational excitation of the OH stretch and164,165 the Claisen rearrangement reaction,166 or to investigate Diels–Alder reactions.167 Such studies provide molecular-level details into the reaction mechanisms and relevant coordinates driving the reaction based on dynamics studies, which goes beyond scanning and relaxing the energy function along a minimum energy path. As an example, for the Diels–Alder reaction between 2,3-dibromo-1,3-butadiene and maleic anhydride such reactive dynamics simulations emphasized the important role played by rotations or the two reactants to reach the transition state.167 

Finally, biological systems can also be studied with a combination of RKHS-based PES and an empirical force field.150,151 Such studies allowed a structural interpretation of metastable states in MbNO and a molecularly refined understanding of ligand exchange (NO vs O2) at the heme-iron in truncated hemoglobin.

In this outlook, recent developments in the energy function design for molecular simulations are discussed. These developments concern improvements in the accuracy of the methods and their generalizations to larger molecules or to applications in condensed phase simulations.

Very recently, accurate representations of multidimensional PESs using PIPs168 or machine learning within the context of kernel ridge regression169 or neural networks133,170,171 have been presented. Common to all these approaches is the fact that they are based on extensive reference datasets. In addition, the evaluation time of such representations is computationally considerably more expensive than that of a parameterized empirical function. For the Diels–Alder reaction between 2,3-dibromo-1,3-butadiene and maleic anhydride, evaluating the MS-ARMD and NN-based PESs differs by a factor of ∼200.167 Such an increase in computing time is problematic for gas-phase studies but may be less detrimental when using such NN-based energy functions in a mixed quantum/molecular mechanics context. For gas-phase studies, the number of trajectories to be run for converged observables can easily reach 106 or more. On the other hand, once the environment (treated with molecular mechanics) becomes sufficiently large, the evaluation time for the molecular mechanics part becomes comparable to that for the ML (i.e., “quantum”) part. Under such circumstances, the added benefit afforded by the much increased accuracy of a machine learned energy function will outweigh the increased computational cost.

As an example for future applications and improvement, the possibility to run MD simulations of molecules with machine learned energy functions and charges in solution is mentioned. PhysNet133 can be used to train energies, forces, and partial charges, which yields accurate representations of high-dimensional, reactive PESs.167,172,173 As has been noted above, accurate representations of the ESP are required for quantitative investigations of the vibrational spectroscopy in the condensed phase. As PhysNet provides geometry-dependent charges in a natural way, simulations of solutes in aqueous and other environments will benefit greatly from such more realistic representations. This also brings atomistic simulations based on empirical energy functions closer to more rigorous mixed quantum/classical simulations. First steps in this direction have been undertaken for protonated water in helium as the solvent.174 

Alternatively, the “bonded terms” [see Eq. (1)] can be represented as a multidimensional tensor product over one-dimensional reproducing kernels, which has been recently done for molecules with up to ten atoms.175 For CH2O, such an approach yields a globally valid PES (errors of ∼0.5 kcal/mol for energies ∼100 kcal/mol above the global minimum), which can be used in molecular dynamics simulations. Combining this with advanced representations for the electrostatics (MTPs,25,26 polarizable DCM,18,19 or electrostatics plus polarization) is also expected to provide considerably enhanced performance for condensed phase simulations.

Another aspect of the force field design that deserves further attention concerns van der Waals interactions, which are often treated rather empirically in simulations using empirical force fields. While for the long range part reliable models based on dispersion coefficients (C6) including their angular anisotropies are available14 and considerable progress has been made in particular for applications in solids,176 far less work has been devoted to the short range, repulsive part. In many widely used empirical force fields, the “mixing rules” (or combination rules) for dissimilar atoms are a weakness177 that can be potentially addressed by using extensions of the techniques discussed in the present work.

It is expected that with the increased performance of computer architecture and the recent advances in machine learning and representing high-dimensional potential energy surfaces, the quality of simulations involving small molecules will further increase. This will bring experiment, simulation, and theory closer together and provide new opportunities to characterize, analyze, and predict molecular properties and processes at an atomic scale. A next step in this endeavor is to transfer and apply this technology to larger systems, including proteins and materials.

D.K., S.M.S., and P.M. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work was supported by the Swiss National Science Foundation Grant Nos. 200021-117810 and 200020-188724, the NCCR MUST, the AFOSR, and the University of Basel, which is gratefully acknowledged.

1.
D.
Koner
,
O. T.
Unke
,
K.
Boe
,
R. J.
Bemish
, and
M.
Meuwly
,
J. Chem. Phys.
150
,
211101
(
2019
).
2.
M. S.
Grover
,
E.
Torres
, and
T. E.
Schwartzentruber
,
Phys. Fluids
31
,
076107
(
2019
).
3.
S.
Lifson
and
A.
Warshel
,
J. Chem. Phys.
49
,
5116
(
1968
).
4.
M.
Levitt
and
S.
Lifson
,
J. Mol. Biol.
46
,
269
(
1969
).
5.
M. J.
Hwang
,
T. P.
Stockfisch
, and
A. T.
Hagler
,
J. Am. Chem. Soc.
116
,
2515
(
1994
).
6.
J. R.
Maple
,
M.-J.
Hwang
,
T. P.
Stockfisch
,
U.
Dinur
,
M.
Waldman
,
C. S.
Ewig
, and
A. T.
Hagler
,
J. Comput. Chem.
15
,
162
(
1994
).
7.
B.
Brooks
,
R.
Bruccoleri
,
B.
Olafson
,
D.
States
,
S.
Swaminathan
, and
M.
Karplus
,
J. Comput. Chem.
4
,
187
(
1983
).
8.
S. J.
Weiner
,
P. A.
Kollman
,
D. A.
Case
,
U.
Singh
,
C.
Ghio
,
G.
Alagona
,
S.
Profeta
, Jr.
, and
P.
Weiner
,
J. Am. Chem. Soc.
106
,
765
(
1984
).
9.
W. L.
Jorgensen
and
J.
Tirado-Rives
,
J. Am. Chem. Soc.
110
,
1657
(
1988
).
10.
J.
Hermans
,
H. J. C.
Berendsen
,
W. F.
van Gunsteren
, and
J. P. M.
Postma
,
Biopolymers
23
,
1513
(
1984
).
11.
J.
McCammon
,
B.
Gelin
, and
M.
Karplus
,
Nature
267
,
585
(
1977
).
12.
P.
Steinbach
and
B.
Brooks
,
Proc. Natl. Acad. Sci. U. S. A.
90
,
9135
(
1993
).
13.
M. S.
Al-Manthari
,
M. A.
Al-Wadhahi
,
K.
Nasrifar
, and
G. R.
Vakili-Nezhaad
,
Int. J. Thermodyn.
22
,
107
(
2019
).
14.
A. J.
Stone
,
The Theory of Intermolecular Forces
(
Clarendon Press Oxford
,
1996
), Vol. 32.
15.
C.
Kramer
,
P.
Gedeck
, and
M.
Meuwly
,
J. Comput. Chem.
33
,
1673
(
2012
).
16.
T.
Bereau
,
C.
Kramer
, and
M.
Meuwly
,
J. Chem. Theory Comput.
9
,
5450
(
2013
).
17.
C.
Kramer
,
T.
Bereau
,
A.
Spinn
,
K. R.
Liedl
,
P.
Gedeck
, and
M.
Meuwly
, “
Deriving static atomic multipoles from the electrostatic potential
,”
J. Chem. Inf. Model.
53
,
3410
(
2013
).
18.
M.
Devereux
,
S.
Raghunathan
,
D. G.
Fedorov
, and
M.
Meuwly
,
J. Chem. Theory Comput.
10
,
4229
(
2014
).
19.
O. T.
Unke
,
M.
Devereux
, and
M.
Meuwly
,
J. Chem. Phys.
147
,
161712
(
2017
).
20.
K.
Vanommeslaeghe
and
A. D.
MacKerell
, Jr.
,
Biochim. Biophys. Acta, Gen. Subj.
1850
,
861
(
2015
).
21.
D.
Nutt
and
M.
Meuwly
,
Biophys. J.
85
,
3612
(
2003
).
22.
P.
Ren
and
J.
Ponder
,
J. Phys. Chem. B
107
,
5933
(
2003
).
23.
L.-P.
Wang
,
T.
Head-Gordon
,
J. W.
Ponder
,
P.
Ren
,
J. D.
Chodera
,
P. K.
Eastman
,
T. J.
Martinez
, and
V. S.
Pande
,
J. Phys. Chem. B
117
,
9956
(
2013
).
24.
A. J.
Stone
,
J. Chem. Theory Comput.
1
,
1128
(
2005
).
25.
P.
Troester
,
K.
Lorenzen
,
M.
Schwoerer
, and
P.
Tavan
,
J. Phys. Chem. B
117
,
9486
(
2013
).
26.
J.
Bertie
and
Z.
Lan
,
Appl. Spectr.
50
,
1047
(
1996
).
27.
D.
Sidler
,
M.
Meuwly
, and
P.
Hamm
,
J. Chem. Phys.
148
,
244504
(
2018
).
28.
R.
Kumar
and
J. L.
Skinner
,
J. Phys. Chem. B
112
,
8311
(
2008
).
29.
X.
Huang
,
B.
Braams
, and
J.
Bowman
,
J. Phys. Chem. A
110
,
445
(
2006
).
30.
S.
Naserifar
and
W. A.
Goddard
 III
,
J. Chem. Phys.
149
,
174502
(
2018
).
31.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
,
J. Chem. Theory Comput.
9
,
5395
(
2013
).
32.
R.
Rein
,
Adv. Quantum Chem.
7
,
335
(
1973
).
33.
G.
Maroulis
,
J. Phys. Chem.
100
,
13466
(
1996
).
34.
J. M. M.
Roco
,
A. C.
Hernández
, and
S.
Velasco
,
J. Chem. Phys.
103
,
9161
(
1995
).
35.
J. M. M.
Roco
,
A.
Medina
,
A. C.
Hernández
, and
S.
Velasco
,
J. Chem. Phys.
103
,
9175
(
1995
).
36.
J. E.
Straub
and
M.
Karplus
,
Chem. Phys.
158
,
221
(
1991
).
37.
D.
Nutt
and
M.
Meuwly
,
Proc. Natl. Acad. Sci. U. S. A.
101
,
5998
(
2004
).
38.
N.
Plattner
and
M.
Meuwly
,
Biophys. J.
94
,
2505
(
2008
).
39.
T.
Clark
,
M.
Hennemann
,
J.
Murray
, and
P.
Politzer
,
J. Mol. Model.
13
,
291
(
2007
).
40.
J.
Murray
,
P.
Lane
, and
P.
Politzer
,
J. Mol. Model.
15
,
723
(
2009
).
41.
P.
Politzer
,
J. S.
Murray
, and
T.
Clark
,
Phys. Chem. Chem. Phys.
12
,
7748
(
2010
).
42.
J. D.
Jackson
,
Classical Electrodynamics
(
John Wiley & Sons
,
New York
,
1998
).
43.
A.
Stone
,
Chem. Phys. Lett.
83
,
233
(
1981
).
44.
W.
Penney
,
Proc. R. Soc. London, Ser. A
158
,
0306
(
1937
).
45.
L.
Pauling
,
J. Am. Chem. Soc.
69
,
542
(
1947
).
46.
H. S.
Johnston
and
C.
Parr
,
J. Am. Chem. Soc.
85
,
2544
(
1963
).
47.
A. C. T.
van Duin
,
S.
Dasgupta
,
F.
Lorant
, and
W. A.
Goddard
 III
,
J. Phys. Chem. A
105
,
9396
(
2001
).
48.
K.
Chenoweth
,
A. C. T.
van Duin
, and
W. A.
Goddard
 III
,
J. Phys. Chem. A
112
,
1040
(
2008
).
49.
T. P.
Senftle
,
S.
Hong
,
M. M.
Islam
,
S. B.
Kylasa
,
Y.
Zheng
,
Y. K.
Shin
,
C.
Junkermeier
,
R.
Engel-Herbert
,
M. J.
Janik
,
H. M.
Aktulga
,
T.
Verstraelen
,
A.
Grama
, and
A. C. T.
van Duin
,
npj Comput. Mater.
2
,
15011
(
2016
).
50.
A.
Warshel
and
R. M.
Weiss
,
J. Am. Chem. Soc.
102
,
6218
(
1980
).
51.
D. R.
Glowacki
,
A. J.
Orr-Ewing
, and
J. N.
Harvey
,
J. Chem. Phys.
143
,
044120
(
2015
).
52.
R.
Valero
,
L.
Song
,
J.
Gao
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
5
,
1
(
2009
).
53.
R.
Valero
,
L.
Song
,
J.
Gao
, and
D. G.
Truhlar
,
J. Chem. Theory Comput.
5
,
2191
(
2009
).
54.
S. C. L.
Kamerlin
,
J.
Cao
,
E.
Rosta
, and
A.
Warshel
,
J. Phys. Chem. B
113
,
10905
(
2009
).
55.
G.
Hong
,
E.
Rosta
, and
A.
Warshel
,
J. Phys. Chem. B
110
,
19570
(
2006
).
56.
A.
Warshel
,
Proc. Natl. Acad. Sci. U. S. A.
81
,
444
(
1984
).
57.
M.
Strajbl
,
G.
Hong
, and
A.
Warshel
,
J. Phys. Chem. B
106
,
13333
(
2002
).
58.
S. J.
Greaves
,
R. A.
Rose
,
T. A. A.
Oliver
,
D. R.
Glowacki
,
M. N. R.
Ashfold
,
J. N.
Harvey
,
I. P.
Clark
,
G. M.
Greetham
,
A. W.
Parker
,
M.
Towrie
, and
A. J.
Orr-Ewing
,
Science
331
,
1423
(
2011
).
59.
D. R.
Glowacki
,
R. A.
Rose
,
S. J.
Greaves
,
A. J.
Orr-Ewing
, and
J. N.
Harvey
,
Nat. Chem.
3
,
850
(
2011
).
60.
D. R.
Glowacki
,
W. J.
Rodgers
,
R.
Shannon
,
S. H.
Robertson
, and
J. N.
Harvey
,
Proc. R. Soc. London, Ser. A
375
,
20160206
(
2017
).
61.
Y. T.
Chang
and
W. H.
Miller
,
J. Phys. Chem.
94
,
5884
(
1990
).
62.
U.
Schmitt
and
G.
Voth
,
J. Phys. Chem. B
102
,
5547
(
1998
).
63.
H. B.
Schlegel
and
J. L.
Sonnenberg
,
J. Chem. Theory Comput.
2
,
905
(
2006
).
64.
J. W.
Petrich
,
J. C.
Lambry
,
K.
Kuczera
,
M.
Karplus
,
C.
Poyart
, and
J. L.
Martin
,
Biochemistry
30
,
3975
(
1991
).
65.
S. G.
Kruglik
,
B.-K.
Yoo
,
S.
Franzen
,
M. H.
Vos
,
J.-L.
Martin
, and
M.
Negrerie
,
Proc. Natl. Acad. Sci. U. S. A.
107
,
13678
(
2010
).
66.
J.
Kim
,
J.
Park
,
T.
Lee
, and
M.
Lim
,
J. Phys. Chem. B
116
,
13663
(
2012
).
67.
B.-K.
Yoo
,
S. G.
Kruglik
,
I.
Lamarre
,
J.-L.
Martin
, and
M.
Negrerie
,
J. Phys. Chem. B
116
,
4106
(
2012
).
68.
V.
Vaida
,
H. G.
Kjaergaard
,
P. E.
Hintze
, and
D. J.
Donaldson
,
Science
299
,
1566
(
2003
).
69.
A. C.
Crowther
,
S. L.
Carrier
,
T. J.
Preston
, and
F. F.
Crim
,
J. Phys. Chem. A
113
,
3758
(
2009
).
70.
P. M.
Hundt
,
B.
Jiang
,
M. E.
van Reijzen
,
H.
Guo
, and
R. D.
Beck
,
Science
344
,
504
(
2014
).
71.
D. R.
Nutt
and
M.
Meuwly
,
Biophys. J.
90
,
1191
(
2006
).
72.
J.
Danielsson
and
M.
Meuwly
,
J. Chem. Theory Comput.
4
,
1083
(
2008
).
73.
J.
Yosa
and
M.
Meuwly
,
J. Phys. Chem. A
115
,
14350
(
2011
).
74.
T.
Nagy
,
J.
Yosa Reyes
, and
M.
Meuwly
,
J. Chem. Theory Comput.
10
,
1366
(
2014
).
75.
M. H.
Schmid
,
A. K.
Das
,
C. R.
Landis
, and
M.
Meuwly
,
J. Chem. Theory Comput.
14
,
3565
(
2018
).
76.
C. R.
Landis
,
T.
Cleveland
, and
T. K.
Firman
,
J. Am. Chem. Soc.
120
,
2641
(
1998
).
77.
T. K.
Firman
and
C. R.
Landis
,
J. Am. Chem. Soc.
123
,
11728
(
2001
).
78.
I.
Tubert-Brohman
,
M.
Schmid
, and
M.
Meuwly
,
J. Chem. Theory Comput.
5
,
530
(
2009
).
79.
S.
Lammers
,
S.
Lutz
, and
M.
Meuwly
,
J. Comput. Chem.
29
,
1048
(
2008
).
80.
K.
Mackeprang
,
Z.-H.
Xu
,
Z.
Maroun
,
M.
Meuwly
, and
H. G.
Kjaergaard
,
Phys. Chem. Chem. Phys.
18
,
24654
(
2016
).
81.
Z.-H.
Xu
and
M.
Meuwly
,
J. Phys. Chem. A
121
,
5389
(
2017
).
82.
K.
Karandashev
,
Z.-H.
Xu
,
M.
Meuwly
,
J.
Vaníček
, and
J. O.
Richardson
,
Struct. Dyn.
4
,
061501
(
2017
).
83.
Z.-H.
Xu
and
M.
Meuwly
,
J. Phys. Chem. B
123
,
9846
(
2019
).
84.
P.
Hamm
,
M.
Lim
,
W.
DeGrado
, and
R.
Hochstrasser
,
Proc. Natl. Acad. Sci. U. S. A.
96
,
2036
(
1999
).
85.
I. T.
Suydam
,
C. D.
Snow
,
V. S.
Pande
, and
S. G.
Boxer
,
Science
313
,
200
(
2006
).
86.
P.
Mondal
and
M.
Meuwly
,
Phys. Chem. Chem. Phys.
19
,
16131
(
2017
).
87.
K. L.
Koziol
,
P. J. M.
Johnson
,
B.
Stucki-Buchli
,
S. A.
Waldauer
, and
P.
Hamm
,
Curr. Opin. Struct. Biol.
34
,
1
(
2015
).
88.
H.-O.
Hamaguchi
and
R.
Ozawa
, in
Advances in Chemical Physics
, edited by
S. A.
Rice
(
2005
), Vol. 131, pp.
85
104
.
89.
M.
Lim
,
T. A.
Jackson
, and
P. A.
Anfinrud
,
J. Chem. Phys.
102
,
4355
(
1995
).
90.
E. S.
Park
,
S. S.
Andrews
,
R. B.
Hu
, and
S. G.
Boxer
,
J. Phys. Chem. B
103
,
9813
(
1999
).
91.
J.
Ma
,
S.
Huo
, and
J.
Straub
,
J. Am. Chem. Soc.
119
,
2541
(
1997
).
92.
J.
Meller
and
R.
Elber
,
Biophys. J.
74
,
789
(
1998
).
93.
M.
Anselmi
,
M.
Aschi
,
A.
Di Nola
, and
A.
Amadei
,
Biophys. J.
92
,
3442
(
2007
).
94.
K.
Nienhaus
,
J.
Olson
,
S.
Franzen
, and
G.
Nienhaus
,
J. Am. Chem. Soc.
127
,
40
(
2005
).
95.
M.
Meuwly
,
ChemPhysChem
7
,
2061
(
2006
).
96.
Y.
Wang
,
X.
Huang
,
B. C.
Shepler
,
B. J.
Braams
, and
J. M.
Bowman
,
J. Chem. Phys.
134
,
094509
(
2011
).
97.
G. R.
Medders
,
V.
Babin
, and
F.
Paesani
,
J. Chem. Theory Comput.
10
,
2906
(
2014
).
98.
G. R.
Medders
and
F.
Paesani
,
J. Chem. Theory Comput.
11
,
1145
(
2015
).
99.
B. J.
Braams
and
J. M.
Bowman
,
Int. Rev. Phys. Chem.
28
,
577
(
2009
).
100.
H.
Liu
,
Y.
Wang
, and
J. M.
Bowman
,
J. Chem. Phys.
142
,
194502
(
2015
).
101.
P.-A.
Cazade
,
T.
Bereau
, and
M.
Meuwly
,
J. Phys. Chem. B
118
,
8135
(
2014
).
102.
P.-A.
Cazade
,
H.
Tran
,
T.
Bereau
,
A. K.
Das
,
F.
Kläsi
,
P.
Hamm
, and
M.
Meuwly
,
J. Chem. Phys.
142
,
212415
(
2015
).
103.
M. W.
Lee
,
J. K.
Carr
,
M.
Göllner
,
P.
Hamm
, and
M.
Meuwly
,
J. Chem. Phys.
139
,
054506
(
2013
).
104.
J.
Desmond
,
D.
Koner
, and
M.
Meuwly
,
J. Phys. Chem. B
123
,
6588
(
2019
).
105.
S. M.
Salehi
,
D.
Koner
, and
M.
Meuwly
,
J. Phys. Chem. B
123
,
3282
(
2019
).
106.
F.
Bernstein
,
T.
Koetzle
,
G. J.
Williams
,
E. F.
Meyer
, Jr.
,
M. D.
Brice
,
J.
Rodgers
,
O.
Kennard
,
T.
Shimanouchi
, and
M.
Tasumi
,
J. Mol. Biol.
112
,
535
(
1977
).
107.
H. M.
Berman
,
T.
Battistuz
,
T.
Bhat
,
W. F.
Bluhm
,
P. E.
Bourne
,
K.
Burkhardt
,
Z.
Feng
,
G. L.
Gilliland
,
L.
Iype
,
S.
Jain
,
P.
Fagan
,
J.
Marvin
,
D.
Padilla
,
V.
Ravichandran
,
B.
Schneider
,
N.
Thanki
,
H.
Weissig
,
J. D.
Westbrook
, and
C.
Zardecki
,
Biol. Crystallogr.
58
,
899
(
2002
).
108.
E. N.
Baker
,
T. L.
Blundell
,
J. F.
Cutfield
,
E. J.
Dodson
,
G. G.
Dodson
,
D. M. C.
Hodgkin
,
R. E.
Hubbard
,
N. W.
Isaacs
,
C. D.
Reynolds
,
K.
Sakabe
,
N.
Sakabe
, and
N. M.
Vijayan
,
Philos. Trans. R. Soc. London, Ser. B
319
,
369
(
1988
).
109.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
,
J. Chem. Phys.
79
,
926
(
1983
).
110.
B. R.
Brooks
,
C. L.
Brooks
 III
,
A. D.
MacKerell
, Jr.
,
L.
Nilsson
,
R. J.
Petrella
,
B.
Roux
,
Y.
Won
,
G.
Archontis
,
C.
Bartels
,
S.
Boresch
,
A.
Caflisch
,
L.
Caves
,
Q.
Cui
,
A. R.
Dinner
,
M.
Feig
,
S.
Fischer
,
J.
Gao
,
M.
Hodoscek
,
W.
Im
,
K.
Kuczera
,
T.
Lazaridis
,
J.
Ma
,
V.
Ovchinnikov
,
E.
Paci
,
R. W.
Pastor
,
C. B.
Post
,
J. Z.
Pu
,
M.
Schaefer
,
B.
Tidor
,
R. M.
Venable
,
H. L.
Woodcock
,
X.
Wu
,
W.
Yang
,
D. M.
York
, and
M.
Karplus
,
J. Comput. Chem.
30
,
1545
(
2009
).
111.
J. B.
Klauda
,
R. M.
Venable
,
J. A.
Freites
,
J. W.
O’Connor
,
D. J.
Tobias
,
C.
Mondragon-Ramirez
,
I.
Vorobyov
,
A. D.
MacKerell
, and
R. W.
Pastor
,
J. Phys. Chem. B
114
,
7830
(
2010
).
112.
W. C.
Swope
,
H. C.
Andersen
,
P. H.
Berens
, and
K. R.
Wilson
,
J. Chem. Phys.
76
,
637
(
1982
).
113.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
114.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
115.
H. C.
Andersen
,
J. Chem. Phys.
72
,
2384
(
1980
).
116.
S.
Nosé
and
M. L.
Klein
,
Mol. Phys.
50
,
1055
(
1983
).
117.
E.
Hairer
,
C.
Lubich
, and
G.
Wanner
,
Acta Numer.
12
,
399
(
2003
).
118.
R.
Stratt
,
Acc. Chem. Res.
28
,
201
(
1995
).
119.
P.-A.
Cazade
,
F.
Hédin
,
Z.-H.
Xu
, and
M.
Meuwly
,
J. Phys. Chem. B
119
,
3112
(
2015
).
120.
J.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
121.
J.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
78
,
1396
(
1997
).
122.
T.
Darden
,
D.
York
, and
L.
Pedersen
,
J. Chem. Phys.
98
,
10089
(
1993
).
123.
W. V.
Gunsteren
and
H.
Berendsen
,
Mol. Phys.
34
,
1311
(
1997
).
124.
T. S.
Ho
and
H.
Rabitz
,
J. Chem. Phys.
104
,
2584
(
1996
).
125.
O. T.
Unke
and
M.
Meuwly
,
J. Chem. Inf. Model.
57
,
1923
(
2017
).
126.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
B.
Mennucci
,
G. A.
Petersson
,
H.
Nakatsuji
,
M.
Caricato
,
X.
Li
,
H. P.
Hratchian
,
A. F.
Izmaylov
,
J.
Bloino
,
G.
Zheng
,
J. L.
Sonnenberg
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M.
Bearpark
,
J. J.
Heyd
,
E.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
N.
Rega
,
J. M.
Millam
,
M.
Klene
,
J. E.
Knox
,
J. B.
Cross
,
V.
Bakken
,
C.
Adamo
,
J.
Jaramillo
,
R.
Gomperts
,
R. E.
Stratmann
,
O.
Yazyev
,
A. J.
Austin
,
R.
Cammi
,
C.
Pomelli
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
V. G.
Zakrzewski
,
G. A.
Voth
,
P.
Salvador
,
J. J.
Dannenberg
,
S.
Dapprich
,
A. D.
Daniels
,
Ö.
Farkas
,
J. B.
Foresman
,
J. V.
Ortiz
,
J.
Cioslowski
, and
D. J.
Fox
, Gaussian09 revision d.01,
Gaussian, Inc.
,
Wallingford CT
,
2009
.
127.
D. T.
Colbert
and
W. H.
Miller
,
J. Chem. Phys.
96
,
1982
(
1992
).
128.
S.
Woutersen
,
R.
Pfister
,
P.
Hamm
,
Y.
Mu
,
D. S.
Kosov
, and
G.
Stock
,
J. Chem. Phys.
117
,
6833
(
2002
).
129.
M. F.
Decamp
,
L.
Deflores
,
J. M.
Mccracken
,
A.
Tokmakoff
,
K.
Kwac
, and
M.
Cho
,
J. Phys. Chem. B
109
,
11016
(
2005
).
130.
S.
Rick
,
S.
Stuart
, and
B.
Berne
,
J. Chem. Phys.
101
,
6141
(
1994
).
131.
P. E. M.
Lopes
,
B.
Roux
, and
A. D.
MacKerell
, Jr.
,
Theor. Chem. Acc.
124
,
11
(
2009
).
132.
U.
Koch
,
P.
Popelier
, and
A.
Stone
,
Chem. Phys. Lett.
238
,
253
(
1995
).
133.
O. T.
Unke
and
M.
Meuwly
,
J. Chem. Theory Comput.
15
,
3678
(
2019
).
134.
M.
Karplus
,
R. N.
Porter
, and
R. D.
Sharma
,
J. Chem. Phys.
40
,
2033
(
1964
).
135.
G.
Schatz
and
A.
Kuppermann
,
J. Chem. Phys.
65
,
4668
(
1976
).
136.
F.
London
,
Z. Elektrochem.
35
,
552
(
1929
).
137.
H.
Eyring
and
M.
Polanyi
,
Z. Phys. Chem., Abt. B
12
,
279
(
1931
).
138.
S.
Sato
,
J. Chem. Phys.
23
,
2465
(
1955
).
139.
J. M.
Hutson
,
Annu. Rev. Phys. Chem.
41
,
123
(
1990
).
140.
J. M.
Hutson
,
Adv. Mol. Vib. Collision Dyn.
1A
,
1
(
1991
).
141.
M.
Meuwly
and
J.
Hutson
,
J. Chem. Phys.
110
,
3418
(
1999
).
142.
M.
Meuwly
and
J. M.
Hutson
,
J. Chem. Phys.
110
,
8338
(
1999
).
143.
T.
Karman
,
A.
van der Avoird
, and
G. C.
Groenenboom
,
J. Chem. Phys.
147
,
084306
(
2017
).
144.
D.
Koner
,
J. C. S. V.
Veliz
,
A.
van der Avoird
, and
M.
Meuwly
,
Phys. Chem. Chem. Phys.
21
,
24976
(
2019
).
145.
J. D.
Bender
,
S.
Doraiswamy
,
D. G.
Truhlar
, and
G. V.
Candler
,
J. Chem. Phys.
140
,
054302
(
2014
).
146.
X.
Tong
,
T.
Nagy
,
J. Y.
Reyes
,
M.
Germann
,
M.
Meuwly
, and
S.
Willitsch
,
Chem. Phys. Lett.
547
,
1
(
2012
).
147.
T.
Hollebeek
,
T.-S.
Ho
, and
H.
Rabitz
,
Annu. Rev. Phys. Chem.
50
,
537
(
1999
).
148.
D.
Koner
,
R. J.
Bemish
, and
M.
Meuwly
,
J. Chem. Phys.
149
,
094305
(
2018
).
149.
M.
Soloviov
and
M.
Meuwly
,
J. Chem. Phys.
140
,
145101
(
2014
).
150.
M.
Soloviov
,
A. K.
Das
, and
M.
Meuwly
,
Angew. Chem., Int. Ed.
55
,
10126
(
2016
).
151.
A. K.
Das
and
M.
Meuwly
,
Angew. Chem., Int. Ed.
57
,
3509
(
2018
).
152.
G.
Sarma
,
Prog. Aerosp. Sci.
36
,
281
(
2000
).
153.
J.
Bertin
and
R.
Cummings
,
Prog. Aerosp. Sci.
39
,
511
(
2003
).
154.
I. D.
Boyd
and
T. E.
Schwartzentruber
,
Nonequilibrium Gas Dynamics and Molecular Simulation
(
Cambridge University Press
,
New York
,
2017
).
155.
J. C. S. V.
Veliz
,
D.
Koner
,
M.
Schwilk
,
R. J.
Bemish
, and
M.
Meuwly
,
Phys. Chem. Chem. Phys.
22
,
3927
(
2020
).
156.
O.
Denis-Alpizar
,
R. J.
Bemish
, and
M.
Meuwly
,
Phys. Chem. Chem. Phys.
19
,
2392
(
2017
).
157.
D.
Koner
,
J. C. S. V.
Veliz
,
R. J.
Bemish
, and
M.
Meuwly
, arXiv:2002.02310 (
2020
).
158.
D.
Koner
,
R. J.
Bemish
, and
M.
Meuwly
,
J. Phys. Chem. A
(in press) (
2020
), arXiv:2002.05087.
159.
G.
Czakó
and
J. M.
Bowman
,
J. Phys. Chem. A
118
,
2839
(
2014
).
160.
B.
Fu
,
X.
Shan
,
D. H.
Zhang
, and
D. C.
Clary
,
Chem. Soc. Rev.
46
,
7625
(
2017
).
161.
I.
Szabo
and
G.
Czako
,
J. Phys. Chem. A
121
,
9005
(
2017
).
162.
S.
Brickel
,
A. K.
Das
,
O. T.
Unke
,
H. T.
Turan
, and
M.
Meuwly
,
Electron. Struct.
1
,
024002
(
2019
).
163.
J.
Yosa Reyes
,
T.
Nagy
, and
M.
Meuwly
,
Phys. Chem. Chem. Phys.
16
,
18533
(
2014
).
164.
J.
Yosa Reyes
,
S.
Brickel
,
O. T.
Unke
,
T.
Nagy
, and
M.
Meuwly
,
Phys. Chem. Chem. Phys.
18
,
6780
(
2016
).
165.
S.
Brickel
and
M.
Meuwly
,
J. Phys. Chem. A
121
,
5079
(
2017
).
166.
S.
Brickel
and
M.
Meuwly
,
J. Phys. Chem. B
123
,
448
(
2019
).
167.
U.
Rivero
,
O. T.
Unke
,
M.
Meuwly
, and
S.
Willitsch
,
J. Chem. Phys.
151
,
104301
(
2019
).
168.
A.
Nandi
,
C.
Qu
, and
J. M.
Bowman
,
J. Chem. Phys.
151
,
084306
(
2019
).
169.
F. A.
Faber
,
A. S.
Christensen
,
B.
Huang
, and
O. A.
von Lilienfeld
,
J. Chem. Phys.
148
,
241717
(
2018
).
170.
J. S.
Smith
,
O.
Isayev
, and
A. E.
Roitberg
,
Chem. Sci.
8
,
3192
(
2017
).
171.
K. T.
Schuett
,
H. E.
Sauceda
,
P. J.
Kindermans
,
A.
Tkatchenko
, and
K. R.
Mueller
,
J. Chem. Phys.
148
,
241722
(
2018
).
172.
S.
Käser
,
O. T.
Unke
, and
M.
Meuwly
,
New J. Phys.
22
,
055002
(
2020
).
173.
S.
Käser
,
O. T.
Unke
, and
M.
Meuwly
,
J. Chem. Phys.
152
,
214304
(
2020
).
174.
C.
Schran
,
F.
Uhl
,
J.
Behler
, and
D.
Marx
,
J. Chem. Phys.
148
,
102310
(
2018
).
175.
D.
Koner
and
M.
Meuwly
, arXiv:2005.04667 (
2020
).
176.
A.
Tkatchenko
and
M.
Scheffler
,
Phys. Rev. Lett.
102
,
073005
(
2009
).
177.
W.
Song
,
P.
Rossky
, and
M.
Maroncelli
,
J. Chem. Phys.
119
,
9145
(
2003
).