Nuclear Quantum Effects from the Analysis of Smoothed Trajectories: Pilot Study for Water.

Nuclear quantum effects have significant contributions to thermodynamic quantities and structural properties; furthermore, very expensive methods are necessary for their accurate computation. In most calculations, these effects, for instance, zero-point energies, are simply neglected or only taken into account within the quantum harmonic oscillator approximation. Herein, we present a new method, Generalized Smoothed Trajectory Analysis, to determine nuclear quantum effects from molecular dynamics simulations. The broad applicability is demonstrated with the examples of a harmonic oscillator and different states of water. Ab initio molecular dynamics simulations have been performed for ideal gas up to the temperature of 5000 K. Classical molecular dynamics have been carried out for hexagonal ice, liquid water, and vapor at atmospheric pressure. With respect to the experimental heat capacity, our method outperforms previous calculations in the literature in a wide temperature range at lower computational cost than other alternatives. Dynamic and structural nuclear quantum effects of water are also discussed.

Introduction

Calculations of reaction free energy profiles and activation barriers are routinely performed within the rigid-rotor and harmonic-oscillator approximation;[1] meanwhile, the more accurate computation of thermodynamic quantities or vibrational spectra is still a great challenge.[2−10] The inclusion of nuclear quantum effects (NQEs), such as zero-point energy (ZPE) or tunneling, is even more difficult.[11−13] Path-integral molecular dynamics (PIMD) and path integral Monte Carlo (PIMC) simulations are accurate, yet highly expensive methods to incorporate NQEs.[14,15] The computational cost of PIMD simulations can be significantly reduced by advanced techniques.[16−18] Recently developed algorithms, such as a colored noise thermostat and quantum thermal bath, are more effective to add quantum effects to classical simulations,[19−21] but settings need to be chosen carefully to prevent zero-point energy leakage.[22,23]

When empirical water models were used in PIMD simulations, several properties deviated more from the experiments than in the classical simulations.[24−27] In these quantum simulations, the liquid water becomes less structured and less viscous. This has been explained by double counting of quantum effects: once in the parameter optimization using experimental data, second in the quantum simulations. This is why several water models were reparametrized for accurate PIMD simulations resulting in q-SPC/Fw,[28] q-TIP4P/f,[24] and TIP4PQ/2005 models.[29] Another solution to avoid double counting is the application of PIMD with ab initio methods[9] or force fields trained on ab initio data.[25,27,30] Numerous methodological developments have also been made to calculate quantum free energy values from PIMD simulations.[31−35]

In routine DFT calculations, with the optimized geometry in hand, the free energies are summed for all the different motions such as translation, rotation, and vibration, using the partition functions of the particle in the box, rigid rotor, and harmonic oscillator (RRHO) models.[36] This approach works satisfactorily for small molecules at ambient temperatures, where the normal modes of vibrations can be considered as decoupled harmonic oscillators. For systems, where anharmonicity is significant, and/or the conformational space is extended, the RRHO fails to reproduce the exact thermodynamic quantities. Recognizing the need to address this issue, more sophisticated approaches use slightly modified partition functions on optimized geometries.[37,38] There are a few methods which can estimate quantum corrections from classical MD trajectories, for example, one- and two-phase thermodynamics methods (1PT, 2PT).[39−41] In those cases the vibrational density of states (VDOS) is determined from molecular dynamics by the Fourier transformation of the velocity autocorrelation function. Quantum corrections are computed by the multiplication of VDOS with weight functions derived from the partition functions of motions. In the 1PT model, only vibrational modes are considered as harmonic oscillators; in the 2PT model, gas phase motions such as rotational and translational modes are also taken into account. The 2PT model is an improved method based on the original work of Berens et al. which corresponds to the 1PT method with anharmonic correction (1PT+AC).[39]

The 2PT and 1PT+AC methods were successfully applied for the calculation of thermodynamic properties of several systems such as Lennard-Jones fluids,[40] water,[39,41−47] organic liquids,[48,49] carbon dioxide,[50] adsorbed urea to cellulose,[51] ionic liquids,[52,53] carbohydrates,[54] mixtures,[55] and interfaces.[56] Heat capacity is generally used as a reference property for the benchmark of force fields.[42,48,49,53] The 1PT/2PT methods are still in continuous development in respect to accuracy and applicability.[57−63]

Here, we propose the Generalized Smoothed Trajectory Analysis (GSTA) method, which is numerically beneficial to the 1PT/2PT methods and, moreover, addresses their limitations arising from the used approximations. Our theory is demonstrated on the exact reproduction of heat capacity and internal energy of a harmonic oscillator. We have chosen different states of water as real-world examples as the heat capacity varies widely between phases, and it is still one of the most investigated materials in computations.[64] Beyond thermodynamic properties, structural and dynamic NQEs are also investigated.

Theory

In a molecular dynamics simulations, the velocity autocorrelation function (VACF) can be defined as follows:where v is the velocity as a function of the time (t). Here, we refer to mass (m)-weighted VACF, but we assume identical masses in the derivations for simplicity. The vibrational density of states (VDOS) is the representation of the autocorrelation function (VACF) in the Fourier domain:where ν is the frequency. Since in our calculations real numbers are used, the Fourier cosine transform is applied. Consequently, the VACF function is the inverse Fourier transform of the VDOS function:Using t = 0, we get the norm of the VDOS function:Applying ν = 0 in eq , we get the norm of the VACF function:

Originally, Berens et al. proposed that quantum corrected density of states can be determined by the multiplication of VDOS with an appropriate weight function w:[39]

Corrections of Vibrational Density of States (VDOS)

One-Phase Thermodynamics (1PT)

In the 1PT method, only vibrational motions are considered. The internal energy iswhere U0 is the reference energy and β = (kBT)−1, kB is the Boltzmann constant, and T is the temperature.

The vibrational weight function w for the energy originates from the quantum harmonic oscillator model:[36]where h is the Planck constant and coth denotes the hyperbolic cotangent function.

The heat capacity is the temperature derivative of the internal energy:The last term in the integral is the weight function for the heat capacity:[36]The weight functions are shown in Figure . Using this weight function w(ν), the heat capacity can be obtained directly from the classical VDOS by integrating over the frequency domain:In the classical limit (h → 0), the 1PT method always gives kB for the heat capacity, which corresponds to the classical harmonic value, so it cannot model anharmonicity.

Figure 1

Weight functions for VDOS.

Two-Phase Thermodynamics (2PT)

Gas-like motions such as translations and rotations are separated from vibrations in the 2PT method. The total VDOS is decomposed into vibrations, translations, and rotations:

Translation and rotation are determined for the center of mass and principal axes of the molecules, respectively. Different weight functions are used for the different motions in the calculation of the thermodynamic properties:

The weight function of translation and rotation is 1/2 for the heat capacity within the 2PT model.[41]

In the classical limit (h → 0), the 2PT heat capacity can vary between kB/2 and kB per atom, not incorporating any anharmonicity, so the 2PT method cannot describe cases where the classical heat capacity is above kB. Another limitation is that the molecular topology needs to remain the same during the simulation, so bond breaking/formation, thus chemical reactions, cannot be modeled. Moreover, intramolecular rotations cannot be described properly with 2PT.

1PT with Anharmonic Correction (1PT+AC)

Berens et al. proposed a quantum correction added to the classical heat capacities:[39]The quantum correction can be determined from the quantum harmonic weight function cΔ given byIf the integral terms are partitioned differently then we get:where the second term is the anharmonic correction. The 1PT+AC internal energy can be gained analogously as the heat capacity:where the anharmonic correction is the deviation of the classical internal energy (Ucl) from the ideal harmonic case:The 1PT+AC model always satisfies the correspondence principle in contrast with the 1PT or 2PT methods:This also implies that the technique is able to describe anharmonic motions to the extent of the method used to obtain the original classical trajectory.

Generalized Smoothed Trajectory Analysis (GSTA)

In the previous sections, we briefly introduced the relevant methods from the literature about the correction of the VDOS function to get quantized thermodynamic properties. In this section, we present a derivation to show that a similar correction can be performed in the time domain on the VACF function and on the coordinates.

Correction of Velocity Autocorrelation Function

Convolution of two functions f and g is defined asis frequently used in digital processing for the smoothing of signals or the filtering of high frequency noises.[65] On the basis of the convolution theorem, the multiplication of the VDOS function with a weight function w is equivalent to a convolution of the Fourier transforms of the two functions. The quantum corrected VACF function is the inverse Fourier transform of the quantum corrected VDOS function:The Fourier transform of the weight function can be used for the quantum correction of the VACF function:From the corrected VACF function, the thermodynamic properties can also be calculated. For instance, the 1PT+AC heat capacity:where γ is the Fourier transform of the weight function in eq according to eq .where csch is the hyperbolic cosecant function. The 1PT+AC internal energy takes the formand the weight function for VACF with the internal energy is given by

The integral in eq fails to converge as the w weight function increases monotonously. Fortunately, in real systems, the maximum frequency in the VDOS is always finite, so the weight function can be cut at a high finite frequency with the application of an exponentially decreasing function:where variable b controls the exponential decay. The decay is faster as b increases. This function becomes unity as b approaches infinity:The γ function can be determined as a limit of an integral:The integral can be evaluated in eq , and finally, the γ function can be expressed as a limit of a piecewise function:where δ denotes the dirac delta function. The norm of the function in eq is 1. The coth function is the primitive integral of the csch2 function. Since in the simulations the data are represented at discrete time intervals Δt, it is useful to write the γ function at discrete time steps in such a way that it is directly applicable for integration, i.e.,where n is the index of the time step. When n = 0, the γ function takes the simple form:The weight functions of VACF are shown in Figure .

Figure 2

Weight functions for VACF.

Correction of Trajectory

The velocity autocorrelation function is actually a convolution function of the velocity with itself, i.e., f = g = v in eq :According to eqs and 2, the VDOS can be written in the form of:Similarly, the quantum corrected counterpart can be written aswhere ṽ is a modified velocity which satisfies eq . In the following steps, we determine ṽ.

Substituting eqs and 35 into eq :Using the convolution theorem:Assuming that w(ν) is a non-negative real-valued function we getIn the time domain, the multiplication by w(ν) is replaced by convolution.Thus, we arrived to a g(t) function whereby convoluting the velocities one can directly obtain ṽ and the quantum corrected vibrational density of states in eq . If one wants to use a general function, which can be applied any atomic velocity function, then the vibrational weight function (eq ) can be chosen:where sech means the hyperbolic secant function.

In the determination of a kernel function, the weight function of the internal energy can also be used:

Similarly to the determination of the γ function in eq , the integral does not converge here either. We could not derive an analytic form for the Fourier transform in eq . In order to perform this Fourier transform numerically in a practical way, the weight function of the internal energy can be split into two parts:The Fourier transform of the second term in eq can be readily evaluated numerically, while the Fourier transform of the first part can be determined analytically in a similar way as we did with the γ function:

The result can be given as a limit of a piecewise function:

So the g function can be calculated:

When the convolution is performed at discrete time step, Δt, the value of the kernel function in the nth step isAt zero timeThe kernels of g and g are shown in Figure . g is represented in discrete points with Δt = 0.5 fs.

Figure 3

Kernels for trajectory filtration

Utilizing a kernel function defined in eq , one can obtain filtered time-dependent variables such as coordinates, velocities, and forces (x̃, ṽ, F̃):Having these variables smoothed, we are able to derive smoothed energy components. The smoothed kinetic energy (Ẽkin) can be obtained straightforwardly from smoothed velocities.

Correction of Potential Energy

For the definition of the smoothed potential energy (Ẽpot), it can be expected to fulfill the following condition:

The smoothed total energy is simply the sum of the smoothed kinetic and potential energies:

The kernel function can depend on several parameters like the temperature or the Planck constant h. In order to connect the classical systems with the quantum systems continuously, a fictitious η variable is introduced. η = 0 corresponds to the classical picture, and η = h is the quantum one. With η, we can perform an integration from the classical state to the quantum state.

It can be shown that the total energy remains conservative upon smoothing:The total differential of the smoothed total energy isThe negative of the first term can be called as work of smoothing:

The smoothed potential energy is defined at a specific time (t0) as a correction on the original potential energy (Epot) with the work of smoothing:Integrating this equation according to time yields the expectation value of smoothed potential energy:

The average change of kinetic energy is equal to the average change of the potential energy (more details in Appendix A):So the mean smoothed total energy isNote that eq corresponds to U1PT+AC in eq . So we arrived to the same result as the 1PT+AC method, which is true for the heat capacity as well:

Quantum-corrected state functions can be determined from the presented smoothed quantities, as shown by the example of heat capacity below.

Computational Details

Calculation of Heat Capacity

According to the Theory section, several estimators can be designed to determine the heat capacity depending on what is corrected: VDOS, VACF, or the trajectory. The quantum correction can be introduced with different functions which correspond to the heat capacity or the internal energy, but the resulting thermodynamic functions can be transformed into each other by integration/differentiation. We used the kernel function of the heat capacity (eq ) for the filtration since it is more convenient; i.e., its analytical form is known.

We performed several (50–120) independent NVE simulations around the T target temperature. The classical temperature (Tcl) was determined simply from the average classical kinetic energy for a particular trajectory:

The isochoric heat capacity can be determined from the slope of the mean total energies vs Tcl functions:The classic temperature originates from the classical normalization factor in eq . The isobaric heat capacity can be determined from the slope of the H vs Tcl (or H̃ vs Tcl) function:

For condensed phases, p = 1 atm was applied instead of the calculated average pressure because its fluctuation was larger than several hundred atm. We used linear regression to determine the heat capacities and their errors. The uncertainties of our calculations are given at the 95% confidence level. Representative fittings are shown in the Supporting Information.

Born–Oppenheimer Molecular Dynamics Simulations

We carried out classical microcanonical normal mode sampling[66] with Gaussian 09 Rev. E for the vibrations of one water molecule using B3PW91/6-311G(d,p) level of theory.[67,68] This functional was chosen since it reproduced the exact frequencies of a water molecule the most accurately with an average error of 6.1 cm–1.[69] Initial energies were distributed equally between modes according to the equipartition theorem. The equations of motion for the nuclear evolution were integrated employing the velocity Verlet algorithm with a 0.1 fs time step.[70] Then, 50 independent 1 ps long trajectories were generated around the desired temperature, and the total energies represented a χ2 distribution.

Classical Molecular Dynamics Simulations

Classical molecular dynamics simulations were performed with the Tinker molecular modeling software.[71] The simulation box always contained 432 water molecules. For the liquid and vapor phases, a cubic box was applied. For Ih ice, water molecules were arranged according to the Bernal–Fowler ice rules[72] having a vanishing net dipole moment.

First, we performed 10 ps long NpT simulations starting from previously equilibrated structures followed by 11 ps long NVE simulations. The last 1 ps data were used for the trajectory analysis. Here, 120 independent trajectories were generated to determine the heat capacity at a given temperature. In order to make the linear fitting more effective, the temperature of the thermostat was varied around the target temperature according to a χ2 distribution. The equations of motion were integrated employing the Bussi–Parrinello algorithm[73] for the NpT ensemble and a modified Beeman algorithm[74] for the NVE ensemble with a 0.5 fs time step. In condensed phase simulations, the particle mesh Ewald (PME) method was applied for the long-range electrostatic interaction with 9.0 Å cutoff distance.[75] For the vapor phase, a larger cutoff of 113.0 Å was used with a 30 × 30 × 30 grid size for the Ewald summation.

We performed simulations for an NVT ensemble to determine the structural properties of the SPC/Fw water model at 298.15 K, with a simulation box size of 23.439 Å. Then, 500 configurations were collected, and after 11 ps long NVE simulations, the last 1 ps trajectory was filtered with the g function with Δt = 0.1 fs to get one filtered structure from each trajectory. These 500 independent structures were used to calculate the distributions of the intramolecular distances as well as the radial distribution functions.

For the calculation of the IR absorption spectrum of the AMOEBA14 water model, 120 independent NVE trajectories were used. All simulations were 20 ps long, and they were equilibrated at 298.15 K before. We used a four term Blackman–Harris window before the Fourier transform of the dipole autocorrelation function.[76]

Path Integral Molecular Dynamics Calculations

PIMD simulations were carried out with AMBER12[77] in a canonical ensemble for 216 water molecules using the SPC/Fw model. The settings were taken from ref (28), but 32 beads were used instead of 24. The length of the cubic simulation box was 18.68 Å, according to the equilibrium density. After 1 ns long equilibration, 1000 structures were collected for each bead in an additional 1 ns long simulation. For the calculation of the isochoric heat capacity, we performed canonical simulations at 288.15 and 308.15 K as well.

In order to determine the isobaric heat capacity for the liquid phase, NpT simulations were also performed at atmospheric pressure in the temperature range from 260.65 to 385.65 K.

Results and Discussion

Harmonic Oscillator Model

Here, we show the effect of two filters on the sum of two noncoupled oscillators at 298.15 K. One oscillator has high frequency (3000 cm–1), while the other has low frequency (100 cm–1). The analytic curves are shown in Figure . Here, g smooths the vibration with high frequency, while g enhances the high frequency motion. The filtered function with g corresponds to the quantum fluctuation.

Figure 4

Filtration of the sum of two harmonic functions.

In the following section, we determine the fluctuation of the coordinates for a single harmonic oscillator. The time evolution of the position:where X is the amplitude in a particular trajectory.

The probability distribution of the position for the harmonic oscillator in canonical ensemble:where κ is the force constant of the harmonic potential. The filtered position isThe probability distribution of the filtered position is given by

This corresponds to the exact quantum fluctuation of the position for the harmonic oscillator. So the GSTA method is exact for the harmonic model for not only in the thermodynamic properties but for the coordinate distribution as well. This is a clear advance of GSTA compared to the 2PT method, since structural quantum effects cannot be investigated with 2PT.

BOMD Simulations of One Water Molecule

Apart from cases where analytical solution exists, like a harmonic oscillator, we have chosen to test the proposed method on a single water molecule as a more realistic system. Born–Oppenheimer molecular dynamics (BOMD) simulations were carried out without rotation in the temperature range from 25 to 5000 K using the B3PW91 functional.[67] The vibrational part of the isochoric heat capacity (cvib) was obtained as the slope of the internal energy vs temperature function. For the sake of comparison, we estimated the isobaric heat capacity (including rotation and translation) based on these data as c = cvib + 1.5R + 1.5R + R = cvib + 4R.

We compare GSTA results with heat capacities determined using several different methods in Figure . The graph also includes the ideal classical value, 7R = 58.2 J/K/mol (blue line) as well as the values determined with the application of the quantum harmonic oscillator model to the optimized geometry (green circles). At low temperature, all methods give very similar results, while at the higher temperature only the trajectory smoothing (red crosses) is close to the experimental values (black triangles).[78] This illustrates that anharmonicity can be described properly by GSTA in contrast with the quantum harmonic oscillator model or the 1PT method (purple squares).

Figure 5

Isobaric heat capacities of a single water molecule from BOMD simulations.

Properties of Bulk Water

In order to demonstrate that our approach can work in very different circumstances, bulk water was considered as a test case in a wide temperature range. The quantum effects are decreasing while anharmonicity is increasing with temperature, and the heat capacity of water changes significantly from zero to 76.0 J/mol/K between different phases at atmospheric pressure.[79,80] The intramolecular vibrations remain in the ground state at ambient temperature; therefore, many rigid water models are successful in reproducing the experimental values.[64,81] For the sake of generality, a flexible model is desired with which similar approaches can be examined even at high temperatures where excited vibrations can be populated as well. Three-site models are beneficial because analytical Hessian can be calculated for the optimized structures, which is not possible for polarizable water models. Considering these requirements, we have chosen the recently developed SPC/Fw[82] water model which performs well on several physical properties.[81] The parameters of the SPC/Fw water model were developed for classical molecular dynamics simulations to reproduce experimental data such as self-diffusion coefficient, dielectric constant, vibrational frequencies, oxygen–oxygen radial pair distribution function, and heat of evaporation.[82]

Static Properties

Classical simulations cannot distinguish between the static properties of light and heavy water. In principle, exactly the same thermodynamic and structural properties should be determined for H2O and D2O since classical statistical thermodynamics predicts that equilibrium properties (heat capacity, molar density, surface tension, etc.) are independent of the masses of atoms. To reveal the different NQEs of classical trajectories, we need post processing methods like 1PT, 2PT, 1PT+AC, or GSTA. An important issue is to decide the reference which the computed quantities are compared to. Converged PIMD simulations provide exact static properties within the limit of the particular water model. However, empirical potentials are developed to reproduce various experimental data with classical simulations. Thus, our strategy is to compare the different approximate values to both PIMD results and experiments as well.

Heat Capacity

At a given temperature, we determined the isobaric heat capacity as the derivative of the enthalpy vs temperature function from 120 independent 1 ps microcanonical trajectories. Three states such as Ih ice, liquid water, and vapor were simulated at temperatures from 25 to 1000 K at atmospheric pressure (Figure ). The classical ideal values are also shown with blue lines for the condensed and vapor phases in the figure (9R = 74.8 J/K/mol, 7R = 58.2 J/K/mol, respectively).

Figure 6

Isobaric heat capacities of bulk water in different phases.

The nuclear motions can be accurately described by independent harmonic oscillators at low temperature; therefore, the heat capacity can be estimated quite well with all methods in the case of hexagonal ice (black triangles in Figure ). Our method (denoted with red crosses) also gives values close to the experimental ones.[79] The maximum deviation is 6.0 J/mol/K at T = 200 K, probably due to the fact that the SPC/Fw water model was developed for liquid water and not for ice. A typical indication of this is that the melting point of the SPC/Fw water model is similar to that of the SPC/E water model (215 K), well below the experimental value (273.15 K).[83]

For liquid water, the harmonic oscillator model is qualitatively wrong as it fails to describe anharmonicity captured already by classical simulations (yellow diamonds in Figure ). The smoothing correction successfully takes the coupled motions of the molecules into account, thus performing outstandingly when compared to the experimental values. The isobaric heat capacity calculated by PIMD is somewhat higher (green circles in Figure ). Most importantly, GSTA performs much better than the classical model or the 1PT method no matter what reference we use, i.e., experimental value or PIMD results. Intramolecular vibrations are filtered out in the g-smoothed trajectories. In this aspect, water molecules behave rigidly in these liquid simulations (see the animation in the Supporting Information). This is in line with the experience that the overall performance of rigid water models is comparable to that of the flexible ones at room temperature.[64,81]

For the vapor phase, we used the same water model and trajectory analysis as in condensed phases (Figure ). GSTA gives reasonable results but overestimates the experimental heat capacities by about 3.2 J/mol/K. This is probably due to the fact that the dipole moment of the SPC/Fw water model is adjusted to liquid properties and considering that it is high compared to the dipole moment of a single water molecule in the gas phase (2.39 vs 1.85 D, respectively). The 1PT method overestimates the experimental values with 16 J/mol/K as a consequence that it considers only harmonic vibrations.

Overall, the smoothed SPC/Fw simulations are able to reproduce the heat capacity with an average error of 3.1 J/mol/K in an extended temperature range at atmospheric pressure. We have found two computational studies in the literature which investigated the heat capacities of the three most important phases of water. Yeh et al. used the 2PT method with different rigid water models, and all calculated heat capacities deviated more than 6.7 J/mol/K from the experiments.[45] Shinoda and Shiga performed PIMD simulations for the three phases using the flexible SPC/F water model.[84] They reproduced the experimental values excellently for ice and gas phase, but the heat capacity of liquid water was significantly underestimated by 16.6 J/mol/K.

We have investigated other water models to examine their performance for the calculated heat capacities of liquid water (at T = 298.15 K, p = 1 atm). Besides the discussed SPC/Fw model, other potentials were also tested including rigid, polarizable, and four-site models. Table includes results from previous works as well for comparison. Regarding the classical heat capacities (ccl), we have reproduced the literature data with small differences, which validates our simulations. GSTA-corrected heat capacities vary between 72.8 and 80.9 J/mol/K with the different water models around the experimental value. The 1PT and 2PT methods give values from 52.2 to 86.6 J/mol/K, and generally, rigid models perform better. The PIMD heat capacities range from 58.2 to 92.8 J/mol/K. The PIMD values can be considered as the exact heat capacity for a particular water model. Previously, the PIMD heat capacity of the SPC/Fw water model was found to be 88.7 J/mol/K, significantly higher than that of the experimental (75.4 J/mol/K) or the GSTA values (72.8 J/mol/K). In our PIMD simulations on an NpT ensemble, the isobaric heat capacity is 79.6 J/mol/K. In order to resolve this contradiction, we performed PIMD simulations on a canonical ensemble with the same settings as in the literature.[28] Using the same number of beads, 24, we obtain 79.8 J/mol/K for the heat capacity, but with 32 beads it decreases to 76.6 J/mol/K. This indicates that there is a discrepancy between the value obtained here and those from earlier simulations. Still, we note that our PIMD heat capacity is in reasonable agreement with the experimental and the GSTA values. Our results indicate that the SPC/Fw water model reproduces the experimental heat capacity sufficiently. On the basis of the comparisons discussed here, we conclude that GSTA is an accurate and robust method in the sense that the calculated heat capacities are less sensitive to the chosen potential.

Table 1

Calculated Isobaric Heat Capacities of Liquid Watera

potential	flexibility	cp	cpcl
experiment		75.4	–
1PT or 1PT+AC
AMOEBA14	flexible	85.7b	110.6b
F3C	flexible	84.5c,r	109.6c,r
SPC/E	flexible	86.6d,r	123.7d,r
WATTS	flexible	72.5e,r	107.3e,r
SPC/E	rigid	73.0d,r	87.4d,r
2PT
AIMD–PBE	flexible	52.2f,r	90.1f,r
F3C	flexible	67.4f,r	104.6f,r
ReaxFF	flexible	73.0f,r	108.4f,r
SPC-Fw	flexible	69.3f,r	109.5f,r
SPC/E	rigid	76.7f,r; 68g	87.3f,r; 87g
TIP3P	rigid	58.3f,r; 60g	67.5f,r; 86g
TIP4P/2005	rigid	77.5f,r; 68g	85.5f,r; 89g
TIP4P-Ew	rigid	69.0f,r; 69g	80.8f,r; 86g
PIMD
SPC/F	flexible	59.0h,r	117.2h,r
SPC/F2	flexible	58.2i,r	117.2i,r
SPC-Fw	flexible	88.7j,r; 79.6 ± 0.5	113.8j,r
q-SPC/Fw	flexible	92.8j,r	–
q-TIP4P/F	flexible	66.7k,r	–
TIP4PQ/2005	rigid	75.3l	–
GSTA
AMOEBA14b	flexible	80.9 ± 0.7	113.5 ± 0.6
F3Cc	flexible	76.3 ± 0.7	108.4 ± 0.5
SPC-Fwm	flexible	72.8 ± 0.6	109.3 ± 0.4
TIP3Fn	flexible	72.8 ± 0.7	106.9 ± 0.6
TIP4P/2005Fo	flexible	80.1 ± 0.9	113.4 ± 0.7
TIP3Pp	rigid	75.0 ± 1.0	81.2 ± 0.9
TIP4P/2005q	rigid	79.0 ± 0.7	87.1 ± 0.6

T = 298.15 K, p = 1 atm; all values are given in J/mol/K.

ref (42).

ref (43).

ref (44).

ref (39); 2PT.

ref (46).

ref (45); PIMD.

ref (84).

ref (85).

ref (28).

ref (26).

ref (86); GSTA: this work.

ref (82).

ref (87).

ref (88).

ref (89).

ref (90).

c was estimated from the given c as c = c + 0.8 J/mol/K.

Structure of Liquid Water

A quantum-corrected structure can be obtained after the filtration of the coordinates with the kernel function of the energy (g). The largest NQE in the structural properties of water is the quantum fluctuation of the hydrogen atoms. This effect is illustrated in the distribution of the intramolecular distances in Figure and also in the animation (Supporting Information). The distribution determined from PIMD simulations can be considered as the exact reference. The distributions of the intramolecular distances become broader compared to the classic results according to the Heisenberg uncertainty principle. The classical distributions are too narrow, but the average distances are almost identical (Table ). The GSTA method reproduces almost perfectly the exact distributions with a slight shift of 0.01 Å for the maxima.

Figure 7

Probability distribution functions of intramolecular distances.

Table 2

Intramolecular Distances of Bulk Water

OH distance/Å
	classic	GSTA	PIMD
mean	1.03	1.04	1.03
standard deviation	0.03	0.07	0.07

The classical and the filtered radial distribution functions are shown in Figure together with the experimental data.[91] The positions of the classical peaks are not altered significantly by the filtration, but they become broadened and get closer to the experimental curves. Note that the 1PT and 2PT methods do not correct the structural properties.

Figure 8

Radial distribution functions of atom pairs in liquid water.

Static Dielectric Constant

The static dielectric constant is another equilibrium property which exhibits nuclear quantum effects. Although, for liquid water the effect of deuterium substitution is small, the experimental dielectric constants are 78.39 vs 78.06 for H2O and D2O;[92] in ice, the isotope effect is the opposite and more enhanced. The static dielectric constants of ordinary and heavy ice are 110 and 124 parallel to the c axis of the crystal at −25 °C, and the difference is even larger at lower temperature.

We have calculated the dielectric constant for both the classical and filtered structures at 298.15 K using the following relationship:where ϵ0 is the vacuum permittivity, V is the volume of the simulation box, and M is the total dipole moment of the system. We computed the dielectric constant for SPC/Fw and AMOEBA14 water models. In Table , we collected our GSTA results along with previous classic and PIMD values from the literature. Our classical values are in good agreement with the literature values. GSTA slightly increases the dielectric constant with about 0.5, which is smaller than its uncertainty. The differences between the classical and filtered trajectories can be calculated pairwise, resulting in a more accurate estimation 0.47 ± 0.03 and 0.6 ± 0.3 for AMOEABA14 and SPC/Fw, respectively. The good agreement with the experimental value is due to the fact that the parameters of the SPC/Fw and AMOEBA14 water models were optimized to reproduce the dielectric constant. The filtered dielectric constant is only slightly larger than the classical value, which is in line with the small difference between the dielectric constants of H2O and D2O (78.39 vs 78.06, respectively).[92]

Table 3

Static Dielectric Constant of Liquid Water

model	classic	GSTA	PIMD
SPC/Fw	80 ± 2a; 79.1 ± 6	79.7 ± 6	64 ± 4b
AMOEBA14	79.4 ± 1.6c; 83.6 ± 7	84.1 ± 7
SPC/F	80 ± 8d		67 ± 9d
SPC/F2	94 ± 10d		74 ± 11d
TTM2.1-F	67 ± 2e		74 ± 2e
TTM3-F	94 ± 1f		81 ± 2g
q-SPC/Fw			86 ± 4b; 90 ± 3h
q-TIP4P/F	58 ± 1g		60 ± 3h
fm-TIP4P/F-TPSS-D3	48.38i		45.69i
experimentj	D2O: 78.06; H2O: 78.39

ref (82).

ref (28).

ref (42).

ref (93).

ref (94).

ref (95).

ref (96).

ref (24).

ref (25).

ref (92).

For the SPC/Fw water model, PIMD gives a smaller dielectric constant by 16 than the classical value. Thus, GSTA predicts a significantly lower quantum effect for the dielectric constant than PIMD (+0.6 vs −16), so even the effect is the opposite. The quantum effect determined by PIMD strongly depends on the water model. In most cases, the dielectric constant decreases drastically, but in some cases, it increases. This implies that the overall effect comes from competing quantum effects, as Habershon et al. proposed for the self-diffusion coefficient of water.[24] They also mentioned that this may affect equilibrium properties such as the melting point of ice.

In order to investigate the possible source of these competing effects, we calculated the molecular dipole moment μ and the Kirkwood g-factor GK for the SPC/Fw water model. (For the AMOEBA14 model, it is not so straightforward to calculate the molecular dipole moment because it is a polarizable water model.) The molecular dipole moment is calculated as a root-mean-square of the individual dipole vectors:where Nmol is the number of the molecules. The Kirkwood g-factor is defined asThe dielectric constant can be expressed from the molecular dipole moment and the GK factor.The molecular dipole moments and the Kirkwood g-factor are collected in Table as well as the static dielectric constants.

Table 4

Components of Static Dielectric Constant of Water

	model	classic	GSTA	difference
ϵs	AMOEBA14	83.6 ± 7	84.1 ± 7	0.47 ± 0.03
ϵs	SPC/Fw	79.1 ± 6	79.7 ± 6	0.58 ± 0.33
μ/D	SPC/Fw	2.3727 ± 0.0003	2.3993 ± 0.0007	0.0266 ± 0.0005
GK	SPC/Fw	4.06 ± 0.31	4.00 ± 0.31	–0.059 ± 0.017

The overall effect is a product of both the molecular dipole moment and the Kirkwood g-factor. The sign of these effects can be positive or negative depending on the water model, but in experiments, the overall effect is small. Most of the PIMD simulations indicate large quantum effects on the static dielectric constant which implies a large effect on the molecular dipole moment and/or on the Kirkwood g-factor as well. This does not mean that PIMD would be inaccurate, just indicates that these force fields were not designed to reproduce the experimental NQE on dielectric constant with PIMD simulations. GSTA predicts small effects on these properties, which is in line with the experiments. The q-TIP4P/F force field predicts the smallest NQE in PIMD simulations, and the difference between the dielectric constants is smaller than the uncertainty (Table ).

Dynamical Properties

The investigation of NQEs on dynamical properties is less straightforward than that on static properties. Even classical simulations show drastic differences for H2O and D2O simply because the deuterium moves slower than protium due to the mass difference. Standard PIMD calculations cannot model time-dependent processes, but there are several approximate PIMD-based methods such as RPMD, CMD, or LSC-IVR simulations which can imitate quantum dynamics.[13,18,97] Hence, it is challenging to separate quantum and classical effects in isotope substitutions, and this makes it difficult to validate new methods like GSTA on dynamical NQEs.

Self-Diffusion Coefficient

The self-diffusion coefficient Ds can be determined from the zero frequency value of VDOS:Since the weight function w(ν) is always 1 at zero frequency, the self-diffusion coefficient does not change with the application of the GSTA method (see eqs and 77). The self-diffusion coefficient can be calculated from the mean square displacement of the atoms using the Einstein equation:From this equation it is also evident that one obtains the same self-diffusion coefficient after filtering the classical trajectory, since GSTA perturbs the classical coordinates with a few tenths of Å which becomes negligible compared to the movement of the atoms at a sufficiently long time. Thus, previous classical self-diffusion coefficients are collected in Table which are identical with the GSTA values as well. PIMD-based approximate quantum dynamics methods like CMD, RPMD, or LSC-IVR change the classical value of self-diffusion coefficient. The quantum effect (relative change to the classical value) varies between 180% and −30% depending on the water model and quantum dynamics method.

Table 5

Calculated Self-Diffusion Coefficients of Liquid Water in Å2/ps

	classic/GSTA		quantum
model	H2O	D2O	method	H2O	D2O	quantum effect	isotope effect
SPC/Fw	0.232a		CMD	0.32b		38%
q-SPC/Fw	0.18b		CMD	0.24b		39%
q-SPC/Fw	0.18b		LSC-IVR	0.50c		178%
SPC/F	0.30d		CMD	0.42d		40%
SPC/F2	0.22d		CMD	0.38d		73%
q-TIP4P/F	0.192e		RPMD	0.221e	0.172e	15%	–29%
MB-pol	0.23f		CMD	0.22g		–4%
revPBE-D3	0.222h		TRPMD	0.16h		–28%
revPBE0-D3	0.267h		TRPMD	0.229h		–14%
fm-TIP4P/F-TPSS-D3	0.3i		RPMD	0.288i		–4%
TTM2.1-F	0.150j	0.140j	CMD	0.225j	0.185j	50%; 32%	–7%; −18%
TTM3F	0.27k		CMD	0.30k		11%
vdW-cx	0.23k		CMD	0.34k		48%
MCY	0.25l,m	0.198m					–21%
SPC rigid	0.44n	0.33n					–25%
SPC flexible	0.61n	0.59n					–3%
SPC/E	0.224o	0.201o					–10%
SPC-MPG	0.275p	0.244p					–11%
TIP4P/2005F	0.22p	0.182p					–17%
experiment				0.230q	0.177r		–23%

ref (82).

ref (28).

ref (97).

ref (93).

ref (24).

ref (98).

ref (30).

ref (9).

ref (25).

ref (94).

ref (27).

ref (99).

ref (100).

ref (101).

ref (102).

ref (103).

ref (104).

ref (105).

Although GSTA does not show any nuclear quantum effect on the self-diffusion coefficient, classical simulations can also reproduce the experimental isotope effect. The experimental self-diffusion of heavy water is 23% slower than that of the light water, and in classical simulations, this varies between 3% and 25%. We have found only two quantum simulations on heavy water with q-TIP4P/F and TTM2.1-F water models, and the calculated isotope effects (−29% and −18%) are in good agreement with the experiment as well as the absolute values.

Infrared Absorption Spectrum

In the classical simulation, the infrared spectrum of the SPC/Fw water model failed to reproduce the experimental spectrum qualitatively because of the harmonic model and the fixed point charges.[28,106] This is why we calculated the infrared absorption spectrum of the AMOEBA14 water model which is polarizable, and the OH bonds are anharmonic. We determined the infrared absorption spectra both from the classical and filtered trajectory at different levels of theory (see details in Appendix B).

For classical simulation, the absorption cross section isUsing the harmonic quantum correction factor on the classical line shape function we obtain[107,108]The absorption cross section is calculated from the filtered trajectories by the expression

Interestingly, almost identical spectra were obtained from the filtered trajectories and from the classical simulation applying the harmonic quantum correction factor (see orange and blue lines in Figure ). This implies that the dipole moment changes linearly with the coordinates during the filtration. The absorption spectrum of the classical trajectory without quantum correction shows much lower intensities (green line). The relative intensities are also different, but the positions of the peaks are the same as in the quantum corrected spectra. Comparison with the experimental spectrum (black line) shows that the quantum corrected peaks deviate with ±100 cm–1, and the relative intensities are reproduced qualitatively.[109] The peaks of the symmetric and antisymmetric OH stretches are partially overlapped due to the anharmonicity of the AMOEBA14 water model. The calculated spectrum also reproduces the smaller peaks at around 200 and 2200 cm–1 which are completely missing from the spectrum of nonpolarizable water models.

Figure 9

Infrared absorption spectrum of AMOEBA14 water model at 298.15 K. Calculated intensities are scaled by a factor of 0.417.

In previous dynamic simulations, the quantum spectra deviated significantly from the ones computed from classical simulations. The classical peaks are always blue-shifted compared to the peaks of quantum simulations, the largest difference is in the OH stretching being shifted by ∼100 cm–1.[9,24,28,64,97,110−114] Assuming that GSTA does not change the position of the peaks, the spectra of the filtered trajectories are always blue-shifted relative to the spectra of quantum simulations independently from the water model. The justification is that in classic simulation the atoms vibrate at the bottom of the potential well and experience less anharmonicity than in real quantum dynamics where the atoms reach higher potential energies because of the zero-point energy. As mentioned in the Self-Diffusion Coefficient section, filtration does not change the frequencies of the vibrations but their intensities. Since GSTA is based on a harmonic approximation, it cannot reproduce the exact position of the vibrational frequencies as Marx already showed for quantum harmonic correction.[107]

In previous simulations, imaginary frequencies were also identified in the vibrational spectrum of water.[112,115] GSTA cannot reproduce these imaginary frequencies, because in classical simulations, the VACF functions are even. Therefore, only real frequencies appear in VDOS after Fourier transform.

Real time autocorrelation functions can also be determined by the maximum entropy analytic continuation (MEAC) method for the calculation of quantum properties.[116−119] In MEAC, the imaginary time correlation function is converted into a real time correlation function, and the former is generated from PIMD simulations. MEAC needs an approximate real time function so-called prior which can be a general function (i.e., flat prior) or a more specific one which comes from a CMD, RPMD, or LSC-IVR calculation. The input of GSTA is only a classical trajectory in the traditional sense that it propagates in real time as a single bead. MEAC is useful well below room temperature; therefore, it is not used for liquid water. Typically, it is applied for liquid para-hydrogen.[117]

Computational Cost and Accuracy

According to the convolution theorem, the same results can be obtained for the calculation of thermodynamic properties both in the time and frequency domains. The Fourier transform should not affect the results. This is true for continuous periodic signals, but a finite number of data are represented in discrete points in MD simulations. This is why the calculated quantities carry numerical errors, and different estimators of the same property can give different values. Here, we illustrate that the same correction of the VDOS and VACF can give significantly different results.

The VDOS function is generally determined via the Fourier transformation of the VACF function. The Fourier transform can be carried out with different numerical integrators. The most common technique corresponds to the left Riemann integral, which is probably the default algorithm in most of the simulation programs:

Using the trapezoidal rule for discrete Fourier transform we get

The calculated 1PT heat capacities are collected in Table . In the generally applied method, the VDOS correction using the left Riemann sum has a large error. It overestimates the classical 1PT heat capacity 9R = 74.826 J/mol/K by 15%. The corrected 1PT heat capacity is overestimated by 4%, but if it is normalized with the classical value (), then the underestimation is 9%. Using the trapezoidal rule for integration, all methods give the same results within 0.01%.

Table 6

1PT Heat Capacities of Liquid Water in J/mol/K

	VDOS correction		VACF correction
	left Riemann sum	trapezoidal rule	left Riemann sum	trapezoidal rule	trajectory smoothing
classical	86.068	74.821	74.826	74.826	74.826
quantum	39.319	37.759	39.290	37.760	37.757
normalized quantum	34.183	37.762	39.290	37.760	37.757

In order to sample the low frequency motions for the calculation of the VDOS, a long simulation time is necessary, and as a rule of thumb, the length of the VACF should be at least as long as the time period of the lowest frequency motion. As was already shown in the Theory section, the calculation and the correction of the VDOS are performed with a double integral (eqs and 16), which is more CPU-demanding than the correction of the VACF which requires only a single integral (eq ). One can expect that the calculations will be more accurate if all data are used; i.e., the length of the VACF equals the simulation length. The γ weight function and the g kernel function converges exponentially to zero; therefore, the contribution of the long times are negligible to the thermodynamic properties. The question is this: what is the maximum time separation which has a significant contribution to the particular property?

We can give an upper limit for the maximum time separation (tmax) with the calculation of the heat capacity of ideal monatomic gases. The velocities can be considered as constant between two collisions, and its autocorrelation function is 1. In this case, there is no quantum effect, and the exact heat capacity is 3/2kB. The heat capacity with the filtration of the velocity considering the integration up to tmax iswhere tanh denotes the hyperbolic tangent function. For the VACF correction

Both estimators lead to the exact classical result with a long tmax. The convergence of the functions are shown in Figure at T = 298.15 K. The exact result is reached within 50 fs with an error smaller than 0.01%. This means that a 50 fs long VACF is enough for the quantum correction, and it is not necessary to calculate the VACF for several ps long. Calculations of the VACF consume considerable amounts of memory which correlate with tmax, and the number of the mathematical operations correlate with the square of tmax.

Figure 10

Convergence of heat capacity of ideal monatomic gas with maximum time separation.

Conclusions

In summary, we have proposed a new method, GSTA for quantum correction of classical and BOMD simulations. Our qualitative findings are summarized in Table regarding the capability and accuracy of the GSTA method compared to other methods. A clear advance of GSTA compared to 1PT or 2PT methods is that the effect of anharmonicity can be determined rigorously using the work of smoothing defined by eq . Another novelty is that structural NQEs can be investigated with the filtration of the coordinates. In more advanced methods, where anharmonicity can be described, the classical dynamics is modified to incorporate NQEs (e.g., ZPEs are added to the different normal modes of vibrations). The good agreement with the experiments indicates the plausibility of our smoothing technique. Zero-point vibrations are automatically taken into account by the proper enhancement of high frequency motions from the classical trajectories. The necessary simulations are orders of magnitude faster than with the golden standard technique, PIMD.

Table 7

Reliability of Different Methods for Estimation of Nuclear Quantum Effectsa

	static properties			dynamical properties
method	heat capacity	structure	dielectric constant	diffusion	absorption spectrum
classical/BOMD	–	–	–	–	–
1PT	+	–	–	–	–
2PT	+	–	–	–	–
1PT+AC	++	–	–	–	–
GSTA	++	++	+	–	+
PIMD	+++	+++	+++	–	–
CMD, RPMD, LSC-IVR	+++	+++	+++	++	++

–, no effect; +, rough estimation; ++, good approximation; +++, exact for a particular potential.

GSTA reproduced large NQEs for heat capacity and structural properties. In contrast to PIMD computations, GSTA does not change the IR absorption spectrum or the dielectric constant significantly, and the self-diffusion coefficient remains exactly the same as the classical value. The main reason is that the input of GSTA is a classical trajectory. The quantum fluctuations are added after the simulation, so it does not change what the initial and final states were and how long it took to get there. In contrast, in PIMD calculations, these properties deviate notably from the classical values. While PIMD requires reparametrized or ab initio models to avoid double counting, GSTA can be applied on empirically derived force fields or model potentials.

Our method offers an alternative way to estimate NQEs routinely in theoretical investigations. In addition, force fields and water models can be improved using GSTA. The proposed method can be easily combined with molecular modeling programs to perform simulations and analyze the trajectories for which our source code is available at GitHub.[120]

Here, we have determined the heat capacity and some structural and dynamical properties for various systems as a proof of concept, but in subsequent studies, we show that other thermodynamic quantities, such as entropy and free energy, which are strongly related to the heat capacity, can be estimated by GSTA. We also intend to test the applicability and the limitation of our method on other NQEs, like tunneling, and further spectroscopic and dynamic properties.