Kai
Töpfer
,
Silvan
Käser
and
Markus
Meuwly
*
Department of Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056 Basel, Switzerland. E-mail: m.meuwly@unibas.ch
First published on 11th May 2022
The double proton transfer (DPT) reaction in the hydrated formic acid dimer (FAD) is investigated at molecular-level detail. For this, a global and reactive machine learned (ML) potential energy surface (PES) is developed to run extensive (more than 100 ns) mixed ML/MM molecular dynamics (MD) simulations in explicit molecular mechanics (MM) solvent at MP2-quality for the solute. Simulations with fixed – as in a conventional empirical force field – and conformationally fluctuating – as available from the ML-based PES – charge models for FAD show a significant impact on the competition between DPT and dissociation of FAD into two formic acid monomers. With increasing temperature the barrier height for DPT in solution changes by about 10% (∼1 kcal mol−1) between 300 K and 600 K. The rate for DPT is largest, ∼1 ns−1, at 350 K and decreases for higher temperatures due to destabilisation and increased probability for dissociation of FAD. The water solvent is found to promote the first proton transfer by exerting a favourable solvent-induced Coulomb force along the O–H⋯O hydrogen bond whereas the second proton transfer is significantly controlled by the O–O separation and other conformational degrees of freedom. Double proton transfer in hydrated FAD is found to involve a subtle interplay and balance between structural and electrostatic factors.
More specifically, polar solvents such as water, can stabilize the transition state (TS) of the reactant.4,13,14 For the Claisen rearrangement the origin of the catalytic effect of the water solvent is the differentially stronger hydrogen bonding to the TS structure (allyl phenyl ether) compared with that to the reactant structure because the TS is more polar which lowers the reaction barrier.15,16 Further stabilization of the TS and lowering of the reaction barrier by the water solvent can be affected by solvent polarizability and dipole–dipole interaction between solute in its TS structure and the solvent.16 The explicit description of solvent effects by full ab initio MD simulations for reactive systems in solution is computationally very demanding and can be typically only followed on the picosecond time scale.17,18 To avoid this, the (reactive) solute can be treated quantum mechanically (QM – ranging from tight-binding to post-Hartree–Fock methods) whereas the remaining degrees of freedom are described by an empirical force field.19–21 Alternatively, reactive force fields have been developed which allow to follow chemical reactions in solution on extended time scales.1,22,23 More recently, another possibility has emerged by developing a machine learned (ML) representation for the solute that can be used together with an implicit24,25 or explicit description of the environment at an empirical level. The present work develops such a scheme and applies it to the formic acid dimer in solution.
The formic acid dimer (FAD) has been the subject of several computational26–41 and experimental studies.36,42–48 In the gas phase, formic acid exists as hydrogen bonded dimers49 making it a prototype for complexes with hydrogen bonds such as enzymes or DNA base pairs.50 Such conformations allow double proton transfer (DPT) to occur which leads to broadening of the line shape for the O–H frequency in IR-spectra.36 The dissociation energy of FAD into two formic acid monomers (FAMs) has been determined from spectroscopy and statistical thermodynamics to be D0 = −59.5(5) kJ mol−1 (∼−14.22 kcal mol−1)47 which compares with −14.23 ± 0.08 kcal mol−1 from recent computations.41
Another quantity that has been accurately measured experimentally is the tunneling splitting for which a recent value of 331.2 MHz is available from microwave spectroscopy.51 From this, a barrier for DPT of 2559 cm−1 (30.6 kJ mol−1 or 7.3 kcal mol−1) from analysis of a 3D model was determined.51. The tunneling splitting from the microwave data compares with values of 0.011 cm−1 (340.8 MHz)52 and 0.016 cm−1 (473.7 MHz)53,54 from infrared spectroscopy. Computationally, a high-quality (CCSD(T)-F12a/haTZ), full-dimensional PES for FAD has been determined and represented as permutationally invariant polynomials which features a DPT barrier of 2853 cm−1 (8.16 kcal mol−1).35 The computed splitting on this PES from using a ring-polymer instanton approach was 0.014 cm−1 (420 MHz).37
Given its thorough characterization in the gas phase, FAD is also an interesting system to study reactions in solution. In particular the balance between DPT and the dissociation of FAD into two monomers provides an attractive aspect of the system to be explored also as a proxy for DNA base-pair dimers.55,56 The homologous series of carboxylic acids (formic, acetic, propionic, butyric) in water have also served as model systems for studying the balance and role of hydrogen bonding, hydrophobicity and entropy changes which is of particular relevance in the context of protein folding and stability.57
For FAD in pure water detailed experimental results are scarce. Early work reported that dimerization is most relevant for the low-concentration range for all carboxylic acids investigated, including formic acid.58 This was subsequently challenged and measurements in 3 m NaCl by potentiometric titrations were interpreted as leading to singly H-bonded, extended formic acid oligomers.59 The hydration of formic acid in pure water was also investigated with neutron diffraction experiments for higher formic acid concentrations.60 At all conditions studied, a pronounced peak for the O⋯O separation with a maximum at 2.7 Å was reported which, together with a peak for the OH⋯O angle at 180°, supports formation of singly H-bonded or cyclic FAD. It is of interest to note that the extended structures proposed by Scheraga are incompatible with the neutron scattering results.60 Surface scattering experiments found that FAD colliding with the liquid water surface almost completely inserts in its cyclic form into bulk water even at low collision energies.61 This indicates that H-bond breaking in FAD is not essential for uptake. Whether or not and for how long FAD remains in its cyclic form inside water was, however, not determined in these studies.61 Accompanying ab initio MD simulations reported that the cyclic dimer has fewer FAD–water hydrogen bonds (i.e. is “more hydrophobic”) than the branched, singly H-bonded FAD isomer.62 Finally, a study considering clathrate formation involving formic acid and water interpreted the measured infrared spectra as partially originating from water-mediated cyclic formic acid dimers in which the two formic acid monomers are H-bonded by a bridging water molecule.63 An independent MD simulation for FAD in water reported57 a negligible stabilization free energy for cyclic FAD of 0.6 kcal mol−1 which is surprising given the substantial gas phase stability of FAD of −16.8 kcal mol−1.33,41 Very recent Raman spectroscopy measurements reported formation of both, singly and cyclic H-bonded FADs.64 Taken together, the available studies suggest that FAD on and in pure water can exist as an equilibrium between cyclic and singly H-bonded dimers. Therefore, the present work will focus on the dynamics of cyclic FAD in pure water.
The aim of the present work is to characterize DPT dynamics in FAD in aqueous solution. First, the methods and the ML/MM embedding are described. Next, the quality of the PES is reported and the structural dynamics, free energies, and rates together with conformational coordinates relevant for DPT of FAD in solution are discussed. This is followed by an analysis of the solvent distributions for reactant and transition state structures, and the effect of the solvent-generated electrical field. Finally, the results are discussed in a broader context.
PhysNet computes the total potential energy E for given coordinates and nuclear charges Z of N atoms from atomic energies Ei and pairwise electrostatic and dispersion interactions according to
![]() | (1) |
For the water solvent the TIP3P model67 as implemented in the GPAW program package75 is used with a cutoff range of 10 Å due to its improved performance compared to that of ASE. The interaction between formic acid and the water solvent includes van der Waals and electrostatic terms which were treated with a cutoff at 14 Å or half the simulation box edge length at most. The van der Waals interactions are calculated using the Lennard-Jones-potential and the parameters from CGenFF76 for formic acid.
The electrostatic interactions between the TIP3P atomic charges of water and the NN-predicted atomic charges of FAD are computed within PhysNet. This allows exploitation of reverse mode automatic differentiation to calculate the necessary derivatives of the fluctuating atomic charges from the atomic coordinates and ensures energy conservation. For the van der Waals and electrostatic interactions between the “high level” treatment (here PhysNet) and MM atoms a switch function s(r)
![]() | (2) |
![]() | (3) |
![]() | (4) |
![]() | (5) |
![]() | (6) |
First, the structure of the system was minimized to release strain. A heating run is performed for 10 ps to reach the respective target temperatures of 300 K, 350 K, 400 K, 450 K, 500 K and 600 K. As no rigorous implementation of a barostat is available in ASE, simulations at constant volume but increased temperature were carried out to qualitatively assess the influence of increased pressure on the DPT dynamics. The system equilibrates within 5 ps during the heating runs. Equilibration is followed by production simulations of 100 ps each, accumulating to a total of 16 ns simulation time at each temperature. Positions and momenta are stored (1) every 10 fs for the complete trajectory and (2) every 1 fs in a time window of Δt = [−20, +50] fs around every attempt of a reaction registered during the simulation (see Analysis).
The implementation is validated by examining energy conservation for a 100 ps NVE MD simulation at 600 K. The total energies at each time step are within a standard deviation of 6.1 × 10−3 kcal mol−1 from the average total energy with a maximum deviation of ±25.2 × 10−3 kcal mol−1 and no drift is found. This confirms that the forces are correctly implemented.
To determine the impact of the fluctuating atomic charges available from PhysNet on the dynamics and energetics, MD simulations were also performed by fixing the atomic charges to those of the TS structure of FAD to retain the symmetry of the system before and after PT. These fixed atomic charges were used for the solute–solvent interaction only. For the internal interaction of the ML part (FAD) fluctuating charges were retained in order to correctly represent the PES for the isolated solute. Finally, as FAD was found to be largely hydrophobic from earlier simulations57,62 and to probe the effect of the electrostatic interactions between solute and solvent on DPT, simulations were also carried out with zeroed charges on all atoms of the solute.
In addition to unbiased simulations, biased simulations were carried out at different temperatures for FAD in the gas phase and in solution. For this, an artificial harmonic potential involving both progression coordinates ξ1 and ξ2 was added to the total potential energy. For FAD in solution the force constant was kbias = 16.1 kcal mol−1 Å−2 and 5 × 106 snapshots were recorded from 1 ns trajectories at each temperature. Both constraining harmonic potentials are centered at the TS coordinate, ξ1,TS = ξ2,TS = 0, and the magnitude of the force constant was chosen to lift the potential of the equilibrium conformation (ξeq = ±0.66 Å) of FAD roughly to the level of the TS conformation. Additional tests were run with lower values of kbias. The resulting potentials of mean force (PMFs) were found to be insensitive of this choice, see Fig. S2 (ESI†).
The MP2 method is a good compromise between speed and accuracy. For example, the barrier for DPT of 6.7 kcal mol−1 compares with a value of 7.2 kcal mol−1 from a morphed MMPT-MP2 force field that is consistent with infrared spectroscopy36 and a barrier height of 7.3 kcal mol−1 from analysis of microwave spectroscopy data.51 Calculations at the CCSD(T) level of theory report somewhat higher barriers of 7.9 kcal mol−1 and 8.2 kcal mol−1, depending on the basis set and additional F12 corrections used.35,36 Furthermore, the dissociation of the dimer into two monomers is of interest. In the gas phase the PhysNet model yields De = −16.8 kcal mol−1 which is also the value from the reference MP2/aug-cc-pVTZ calculations. Ab initio calculations at the higher CCSD(T)/aug-cc-pVTZ level of theory find De = −16.8 kcal mol−133,41 and −16.0 kcal mol−1 at the basis set limit.33
Next, the energetics of increasingly hydrated FAD is determined from the present ML/MM energy function, from the CGenFF76 parametrization of FAD together with the TIP3P water model, and from electronic structure calculations. For this, 50 snapshots from a 2 ns ML/MD simulation with FAD in its dimeric structure were extracted. The 15 water molecules closest to the center of mass of FAD were retained for each snapshot. Electronic structure calculations for this analysis were carried out at the B3LYP+D3/aug-cc-pVDZ level of theory82–84 because MP2/aug-cc-pVTZ calculations for FAD surrounded by up to 15 water molecules are computationally too demanding. The total interaction energies for FAD–(H2O)n (n = 1,…,15) complexes were determined. Water molecules were retained in increasing order of their distance from FAD. The results are reported in Fig. S3 and S4 (ESI†) and the analysis shows that the interaction energies with 15 water molecules from the MM/ML energy function differ on average by less than 1 kcal mol−1 from the reference B3LYP+D3 calculations whereas CGenFF underestimates the reference results by about 4 kcal mol−1.
The 50 snapshots were also analyzed from retaining the 4, 8, and 12 closest water molecules and correlating the MM/ML and CGenFF interaction energies with those from the reference B3LYP+D3 calculations (Fig. S5–S7, ESI†). Mean absolute errors between MM/ML and DFT reference calculations are 1.13, 2.06 and 2.85 kcal mol−1 compared with 2.42, 3.36 and 3.62 kcal mol−1 for CGenFF. It is also worth to point out that CGenFF finds a dimer stabilization energy of only 9 kcal mol−1 in the gas phase compared with high-level electronic structure calculations that yield −16.8 kcal mol−133,41 and −17.4 kcal mol−1 at the B3LYP+D3/aug-cc-pVDZ level which is −0.6 kcal mol−1 lower than for the MP2 and CCSD(T) levels of theory with the aug-cc-pVTZ basis set. Overall, the comparison with the B3LYP+D3 results validates the quality of the ML/MM energy function whereas the CGenFF parametrization is found to considerably underestimate the stability of FAD.
In conventional force fields the partial charges are fixed76 and do not change with conformation. PhysNet provides geometry-dependent charges, see Fig. S8 (ESI†), which fluctuate by ±0.05 e to ±0.1 e around their mean. The average charges for the C, HC and one of the oxygen atoms from PhysNet are similar to those in CGenFF whereas for the HO and the second oxygen they differ by about 0.1 e. It is also found that the charge distributions from PhysNet differ for cyclic and branched FAD. Given these differences it is of interest to determine the differences between the simulations analyzed with conformationally fluctuating charges and those with fixed partial charges on the FAD. The charges assigned to the atoms are those of the TS structure to ensure a symmetric atomic charge distribution before and after successful DPT. The distribution of formed H-bonds for FAD in water with fluctuating and fixed charges in Fig. 2 shows that the propensity for H-bond formation significantly depends on fixing the charges. With the fixed charge model the probability for the doubly H-bonded dimer to exist is larger at all temperatures from 300 K (98.7%) to 600 K (3.3%) compared with the fluctuating charge model (97.0% and 0.4% at 300 K and 600 K, respectively). In addition, the probability for singly H-bonded FAD and two separate monomers at 600 K differs with the fixed charge model with 6.4% and 90.3%, respectively, from the fluctuating charge model with 0.7% and 98.8%. Thus, a fixed charge model with charges from the TS structure predicts a higher propensity for FAD to exist in solution compared to the fluctuating charge model.
The higher propensity for dimer dissociation in solution using fluctuating rather than fixed charges correlates with the differences in the magnitudes of the atomic charges between the two charge models. From the equilibrium conformation of FAD towards separated FAMs the polarization along the atoms in the H-bond increases. The fluctuating charges along the Odon–H⋯Oacc bond change from [−0.39, +0.32, −0.43] e in the equilibrium dimer conformation to [−0.42, +0.36, −0.45] e for separated monomers, respectively, which compares with [−0.40, +0.30, −0.40] e in the fixed charge models with charges taken from the TS structure. Comparable amplitudes for charge transfer (δq = 0.025 e) have been found for the HOdon–H⋯OaccH2 hydrogen bond in bulk water when the energy function was fit to match experimental frequencies and intensities of the infrared and Raman bands.85 The larger absolute values of the atomic charges in the fluctuating charge model (differences of up to 20% and 13% for the hydrogen and oxygen atoms of the H-bonds relative to the charges in the TS structure, respectively) lead to stronger electrostatic interactions with the surrounding solvent. This also impacts on the conformational sampling and the reaction barrier height.
To evaluate the effect of the water solvent on dissociation, NVT simulations for FAD in the gas phase were performed using the same setup as for the simulations in solution. The presence of both H-bonds is a prerequisite for DPT to occur as the top of the barrier becomes inaccessible if even one H-bond is broken. Thus, dissociation of the FAD competes with DPT and therefore governs the DPT rate. In the gas phase it is found that – similar to the situation in solution – the propensity towards broken H-bonds increases for higher temperatures. At 300 K both H-bonds are formed for 99.4% of the propagation time. This fraction decreases at higher temperatures: 27.3% and 2.1% at 500 K and 600 K, respectively.
Simulations were also carried out in the NVE ensemble for FAD in the gas phase and yield a free energy barrier for DPT of 8.25 kcal mol−1 at a temperature of 310 K which was determined from the kinetic energy of all atoms in this biased simulation. Earlier classical ab initio metadynamics simulations with the same reaction coordinates as those used in the present work at the BLYP level with a plane wave basis for FAD in the gas phase reported a free energy barrier height of 4.7 kcal mol−1 at 300 K compared with the potential barrier height of 5.0 kcal mol−1.32 This differs qualitatively from the present results for which the free energy barrier in the gas phase increases from the potential barrier of 6.7 kcal mol−1.
To quantify the electrostatic contribution of the solute/solvent interactions to the free energy barrier, biased simulations of FAD in solution were performed with the atomic charges of FAD set to zero. At 300 K the free energy of activation is 7.46 kcal mol−1 and increases by +0.87 kcal mol−1 to 8.33 kcal mol−1 at 600 K. This compares with 7.45 and 7.53 kcal mol−1 with fluctuating and fixed charges at 300 K, respectively, which increase to 8.13 and 8.39 kcal mol−1 at 600 K. At 300 K the contribution of fluctuating electrostatics is thus −0.01 kcal mol−1 as opposed to + 0.07 kcal mol−1 if the charges are fixed and a reduction by −0.20 kcal mol−1 and an increase of +0.06 kcal mol−1 at 600 K for fluctuating and fixed charges, respectively. This finding suggests that FAD behaves largely as a hydrophobic solute for which direct electrostatic interactions with the solvent are less relevant.
The T-dependence to the free energy (−TΔS) suggests that the change ΔS between the reactant and the TS state is negative for all charge models used but the magnitude of ΔS differs. This difference is accommodated in the solvent ordering around the solute. Linear interpolation of the free energy barriers at the different temperatures yields an activation enthalpy and activation entropy of [ΔU, ΔS] = [6.87 kcal mol−1, −2.1 cal mol−1 K−1] for fluctuating charges, [6.62 kcal mol−1, −2.9 cal mol−1 K−1] for fixed charges, and [6.70 kcal mol−1, −2.7 cal mol−1 K−1] for zero charges.
Because the constraints in the biased simulations favour FAD over two separated monomers (and the dimer has been found to dissociate at higher temperatures for free dynamics) it is mandatory to also run and analyze unbiased simulations that allow DPT to occur in solution. For this, an aggregate of 16 ns simulations was run and analyzed at each temperature. As the barrier for DPT is relatively high, DPT is a rare event. DPT rates in unbiased simulations are highest at 350 K for both fluctuating and fixed charge models, see Fig. 2A and B. With the fluctuating charge model (Fig. 2A) the rate is 0.9 ns−1 at 350 K and decreases to 0.4 ns−1 at 500 K, while the dimer is only present for 82.4% and 6.4% of the propagation time, respectively. Hence, the true rates are ∼6.5 times higher at 500 K compared with those at 350 K. Within transition state theory energy barriers of ∼8 kcal mol−1 from the biased simulations correspond to rates of ∼0.1 ns−1 which is consistent with results from unbiased simulations shown in Fig. 2.
With the fixed charge model a similar increase of the DPT rates is observed. At 300 K, no DPT was observed during 8 ns of propagation time whereas for T = 350 K and 500 K they are 1.0 ns−1 and 0.5 ns−1, respectively. As the simulation time with fixed charges is 8 ns – shorter than simulations with the fluctuating charge model – the rates for DPT are based on fewer events with larger uncertainties which are not quantified in this work.
The first PT leads to a short-lived, metastable ion-pair with one protonated formic acid monomer and the corresponding anion, see Fig. S1 (ESI†). Subsequently, the transferred hydrogen atom returns either to where it came from (attempted DPT) within an average delay time of 4.8 ± 2.9 fs at 350 K or the hydrogen atom along the second H-bond transfers (successful DPT) with an average delay time of 5.9 ± 2.6 fs. At 450 K the analysis of unbiased simulations find delay times of 5.1 ± 3.9 fs and 8.2 ± 7.1 fs for attempted and successful DPT, respectively. This is consistent with results from earlier work that reported an average delay time of ∼8 fs for successful DPT.30
Fig. 4A and D compare the instantaneous Odon–Oacc separations at the time of the second PT. The probability distribution functions P(dOO,A, dOO,B) for the O–O separations in unbiased simulations at 350 K and 400 K peak at (2.75 ± 0.18) Å and (2.78 ± 0.22) Å, respectively. These distances shorten significantly to (2.40 ± 0.05) Å and (2.42 ± 0.06) Å for successful DPT; see Fig. 4. For successful DPT both O–O separations are between 2.4 Å and 2.6 Å whereas for attempted DPT the first O–O separation ranges from 2.3 Å to 2.5 Å but the second oxygen donor-acceptor distance covers a range between 2.4 and 2.9 Å. This suggests that a more symmetrical geometry at the time of the first PT increases the probability for successful DPT during the second step, i.e. if the distances between both oxygen atoms in the H-bonds are simultaneously contracted.
Next, the CO bond lengths are analysed. Fig. 4B, C, E, and F shows the distribution of the COacc double and COdon single bond distances from all MD simulations at 350 K and 400 K. The maximum peak positions are clearly distinguishable and demonstrate the single and double bond character of the CO bonds with BO = 1 and BO = 2, respectively. In addition, the histograms for the CO distances at the second PT of a successful (panels B and E) and attempted (panels C and F) DPTs are reported. For successful DPT the distributions overlap for both temperatures and BO = 1.5 whereas for attempted DPT they are clearly non-overlapping and closer to BO = 1 and BO = 2, respectively.
The correlation between narrowly overlapping distance distributions of the CO single and double bond and the higher probability for successful DPT can be explained by the respective OH bond potential in the protonated formic acid monomer. Fig. 5 shows the PhysNet potential along the progression coordinate ξ1 of a hydrogen atom in the TS of FAD. With all other bond lengths frozen in the TS conformation, a scan along ξ1 is performed for three different distances of the CO bond involved in the first H-bond as in A) the TS conformation, B) a CO single bond and C) CO double bond (see illustration in Fig. 5). The remaining CO bonds are kept frozen at the TS bond length of 1.26 Å. For dTS(CO) = 1.26 Å (TS conformation with BO = 1.5) the potential minimum along the reactive coordinate is at ξ1 = 0 Å, which changes if the CO bond contracts (B) or extends (C). The OH-stretch potential is energetically more attractive in the covalent bonding range if the CO bond has BO = 1 rather than BO = 2. Vice versa, the OH-stretch potential becomes repulsive at shorter CO distances as in CO double bonds. Thus, hydrogen abstraction from an oxygen donor atom involved in a CO double bond is energetically more favourable compared to that involved in a CO single bond.
![]() | ||
Fig. 6 The radial distribution function g(d) for the water–oxygen atoms around hydrogen atoms HA and HB of FAD in the reactant (A) and TS (B) structure from the from ML/MM (blue line) and CGenFF (dashed green line) simulations at 300 K. Panel (A) also shows g(d) from one entire ML/MM simulation at 400 K (orange line) including formic acid dominantly as cyclic dimer and separated monomers (see Fig. 2A). Panels (C and D) report the solvent distribution functions P(x, y) around FAD from ML/MM MD simulations in its reactant (C) and TS (D) structures, respectively, and panel (E) reports the difference ΔP(x, y) = PD − PC between the two. |
The solvent radial distribution functions around the hydrogen atoms of the H-bonds in Fig. 6A and B exhibit a weak double peak between 4 and 6 Å. The maxima indicate the first solvation shell around FAD and the double peak structure is more pronounced for simulations with FAD in its TS structure. In the biased simulations FAD is more rigid which leads to a more structured solvent distribution. The radial distribution also shows a broad maximum around 8 Å that can be associated with a second solvation shell. The solvent radial distribution from simulation at 400 K (orange line in Fig. 6A) features an additional peak around 2 Å that represents H-bond formation between branched, i.e. non-cyclic FAD, separated formic acid monomers and water molecules. Such H-bonds become available in simulation at temperatures higher than 300 K as the probability for almost exclusively cyclic FAD decreases significantly (see Fig. 2A).
For a more comprehensive characterization of the changes in the hydration of FAD between product/reactant geometries and the TS structure, 2-dimensional solvent distributions were determined. Fig. 6C and D show the actual distributions projected onto the xy-plane which is the average plane containing the solute. Fig. 6E reports the solvent density difference around FAD in its reactant and TS conformation, i.e. the difference density from panels C and D in Fig. 6. In the 2d-histograms the first and second solvation shells are clearly visible as was also found from the radial distribution functions (Fig. 6A and B). There is also a dent for the first solvation shell at (−4, 0) Å and (4, 0) Å along the horizontal axis joining the hydrogen atoms along the H-bonds of FAD.
The solvent distribution difference in Fig. 6E indicates a shift of the first solvation shell closer to the side of the CH group of FAD for the biased simulations due to the short monomer–monomer distance in the TS structure compared with the reactant state. Hence, the solvent probability in this region is higher for FAD in its TS than in its reactant state. This is consistent with the shorter monomer–monomer distance of FAD in the TS conformation. The monomer–monomer contraction in the TS also flattens the dent towards the hydrogen atoms in the H-bond as shown by the blue area that indicates a lower solvent distribution density around FAD in the TS than the reactant state conformation. Finally, the solvent distribution difference exhibits the symmetry of FAD which indicates the convergence of the simulations.
In the following the electric field at the position of the transferring hydrogen atoms along the hydrogen bond was analyzed. For this, the force on the transferring hydrogen atom at the time of PT was determined as the linear interpolation of the forces for the 1 fs-frame before and after PT. The “time of PT” was determined by the sign change of the linearly interpolated reactive coordinate of the respective H-bond along which PT occurs. In Fig. 7 and Fig. S10 (ESI†) the first PT – to which hydrogen atom HA is attributed – is associated with the horizontal x-axis. In a successful DPT, the second PT is performed by hydrogen atom HB (vertical y-axis) whereas in an attempted DPT the second PT is also performed by HA (back-transfer).
![]() | ||
Fig. 7 Correlation between the magnitude and direction of the electrostatic component of the solvent-generated Coulomb force ![]() |
Fig. 7 and Fig. S10 (ESI†) report the projection of the solvent-generated (MM) Coulomb (C) force FMMC,A along the Odon–HA⋯Oacc and Odon–HB⋯Oacc hydrogen bonds for successful and attempted DPTs during the first and second PT, respectively. As per definition, a positive force always supports PT of a hydrogen atom whereas a negative force has an inhibitory effect. The analysis was carried out for the second PT independent of whether transfer of HB did (successful DPT) or did not (attempted DPT with HA transferring back) occur.
In the unbiased simulations, the first PT is typically accompanied by a supportive Coulomb force. This effect is more pronounced at higher temperatures (T = 400 K, see Fig. 7B) than at lower ones (see Fig. 7A). Biased simulations increase the probability for DPT and provide improved statistics which also allows to fit 2-dimensional Gauss distributions to the data (Fig. 7C and D). These results demonstrate that for the first PT the field along the HA-bond is typically supportive. Interestingly, for an attempted DPT the field along the first HA-bond at the time of PT is supportive but FMMC,B along the second HB-bond is more inhibiting. This applies to both, unbiased and biased simulations, see red distributions in Fig. 7. Regarding Fig. 7 the gravity center of the average solvent-generated Coulomb force is [MMC,A,
MMC,B] = [0.24, −0.23] and [0.86, −1.50] kcal mol−1 Å−1 for successful and attempted DPT at 350 K, respectively, and shifts to [0.51, −0.56] and [−0.93, −1.65] kcal mol−1 Å−1 at 400 K. In biased simulations the gravity centers are at [0.45, −0.43] and [0.56, −0.78] kcal mol−1 Å−1, [0.56, −0.55] and [−0.51, −0.81] kcal mol−1 Å−1 accordingly. Hence, for successful DPT the opposing electrostatic field along the second H-bond is considerably smaller in magnitude than for attempted DPT in all cases.
For the second PT the solvent-generated Coulomb force along both H-bonds is in general inhibiting (red distributions in Fig. S10, ESI†). In unbiased simulations at 350 K the center of gravity of the average solvent-generated Coulomb force is at [MMC,A,
MMC,B] = [−0.21, −0.22] and [−0.84, −1.52] kcal mol−1 Å−1 for successful and attempted DPT, respectively. This shifts to [−0.52, −0.58] and [−0.92, −1.64] kcal mol−1 Å−1 at 400 K. In biased simulations the gravity centers are at [−0.45, − 0.42] and [−0.56, −0.78] kcal mol−1 Å−1 at 350 K and [−0.56, −0.54] and [−0.51, −0.81] kcal mol−1 Å−1 at 400 K. Except for the sign change in the projection of the Coulomb force on the hydrogen atoms of the HA-bond the center of gravity does not change significantly between the first and second PT in successful and attempted DPT.
Comparison of the solvent-generated Coulomb force for successful and attempted DPT further establishes that the field of the water molecules at the position of the transferring hydrogen atom does not change appreciably within ∼14 fs, which is the maximum time between first and second PT observed in the simulations at 350 K. This is expected as the rotational reorientation time of water in solution is on the order of several picoseconds.88 The time evolution of the Coulomb force along the H-bonds is shown in Fig. S11 and S12 (ESI†) for an example of a successful and attempted DPT, respectively.
In summary, successful DPT for FAD in water is often accompanied by a mildly supporting solvent-generated Coulomb force along the first PT in the HA-bond whereas for the second PT in the HB-bond the force is inhibiting but smaller in magnitude than for the back transfer along the first HA-bond. As the inhibiting force along the HB-bond increases, back-transfer along the first HA-bond becomes favoured and results in attempted DPT with FAD returning to its initial reactant conformation, see Fig. S1 (ESI†). The findings complete the picture that the probability for successful DPT primarily depends on the conformation of FAD itself at the first PT event. At conformation of high symmetry and low reaction barriers for PT events in both H-bonds, the solvent-generated Coulomb field impacts the probability of performing one PT in both H-bonds each (successful DPT) or a forth and back PT along one H-bond (attempted DPT).
Spontaneous DPT is observed on the nanosecond time scale at 300 K to 350 K with a maximum rate of 0.9 ns−1 at 350 K. The total ML/MM energy function used was validated vis-à-vis B3LYP+D3 calculations and agrees to within ∼1 kcal mol−1 for hydrated FAD with up to 15 surrounding water molecules. One particularly interesting aspect of the dynamics of FAD in solution is the fact that dissociation into a branched, singly H-bonded structure or into two separate formic acid monomers competes with DPT. The present results (see Fig. 2) suggest that on the ∼14 ns time scale the probability to form branched, singly H-bonded structures is considerably smaller than for cyclic FAD to exist (which decreases as T increases) or that of forming two separated monomers (which increases as T increases). Simulations at higher temperature allow to qualitatively assess conditions at higher pressure which has been used to induce proton transfer in other systems.89,90
Results from previous simulation studies are somewhat ambiguous. Recent ab initio MD simulations using the BLYP functional including dispersion correction point towards pronounced hydrophobicity of cyclic FAD62 as was also found from earlier simulations using a fully empirical force field.57 Both are consistent with the present work. On the other hand, simulations with this non-reactive, fully empirical force field report facile transitions from FAD to branched structures.57 Such branched structures were also implied from experiments which were, however, carried out in 3 molal NaCl.59 Analysis of the T-dependent free energies for DPT yields a change in solvation entropy of ΔS = −2.1 cal mol−1 K−1. It is of interest to juxtapose this value with earlier work that rather focused on dissociation of FAD into two separated monomers for which ΔSD = −3 cal mol−1 K−1 was reported also using the TIP3P water model.57 However, it is quite likely that the two entropy values which differ by about 30% are not directly comparable as dimerization is driven by the different arrangements of solvent molecules around FAD and the FAMs whereas for DPT the reorganization is less pronounced; see Fig. 6E.
As the calculated pKD underestimates the experimentally measured one by 1 to 2 units (0.25 vs. 1.2 to 2.1),57 it is expected that the energy function used in these earlier free energy simulations for cyclic FAD57 underestimate the enthalpic contribution by about one order of magnitude. The reported stabilization from umbrella sampling simulations using the CHARMM22 force field was ΔG ∼ −0.5 kcal mol−1. Hence, to be consistent with the experimentally determined pKD, the enthalpic contribution to the stabilization of FAD in solution should rather be −5 kcal mol−1 or larger. Adding this to the stabilization energy for FAD in the gas phase using the CGenFF energy function (−9 kcal mol−1) yields an estimated stabilization of −14 kcal mol−1 or larger which is more consistent with results from recent high-level electronic structure calculations which find a gas-phase stability for cyclic FAD of −16 kcal mol−1.33,41 In other words, the unreactive, empirical CGenFF energy function considerably underestimates the stabilization of cyclic FAD in the gas phase and in solution which is also confirmed by the results in Fig. S3–S7 (ESI†) and leads to ready dissociation of FAD into two FAMs from such simulations.57
Successful DPT is found to chiefly depend on the conformation of the cyclic FAD. In particular, increased symmetry with simultaneously contracted O–O distances in both H-bonds and CO bonds with BO ≈ 1.5 lowers the reaction barrier for PT along both H-bonds. Successful DPT requires the solvent induced electrostatic force to be only mildly inhibiting ∼−0.5 kcal mol−1 Å−1 along the second H-bond. If the inhibiting force originating from the solvent is too strong along the second H-bond, attempted DPT dominates. For DPT, the MD simulations show an average delay between two PT events of ∼5 and ∼6 fs for attempted and successful DPT at 350 K, respectively, and the solvent-generated Coulomb field assists the first PT. A time difference of ∼6 fs between the successive PT events for successful DPT compares with a time scale of ∼155 fs for the Odon–Oacc vibration. This suggests that DPT in hydrated FAD is essentially concerted and not stepwise; both PTs occur within one Odon–Oacc vibration period. For a stepwise DPT the intermediate state (both protons on one FAM) would need to stabilize, i.e. have a lifetime, which is, however, not what is found here.
From a chemical perspective, one point of particular note is the successful and correct description of the change between bond orders 1 and 2 for the CO bond in PhysNet. Conventional empirical force fields encode the bond order in the equilibrium separation of the bonded term which does not allow easily for changes in the bond order depending on changes in the chemical environment. Although there are examples for capturing such effects91 making provisions for it in the context of empirical force fields is cumbersome. As Fig. 4B, C, E, and F demonstrate, a NN-trained energy function based on reference electronic structure calculations successfully captures such chemical effects.
In summary, a ML/MM MD scheme was implemented and applied to DPT for FAD in solution. The results from ∼100 ns of ML/MM MD simulations show that DPT and dissociation into two FAMs compete depending on temperature. FAD is predominantly hydrophobic which agrees with earlier findings, and the rate for DPT of ∼0.1 ns−1 is consistent between biased and unbiased simulations. Extending the present work to other reaction types and protein–ligand binding will provide deeper chemical understanding and improved models suitable for statistically significant sampling to give molecular level insight into processes in the condensed phase. Contrary to straight AIMD simulations at the correlated level, an ML/MM MD ansatz is feasible even for nanosecond time scale simulations.
Footnote |
† Electronic supplementary information (ESI) available: It reports additional graphics S1–S14. See DOI: https://doi.org/10.1039/d2cp01583h |
This journal is © the Owner Societies 2022 |