Quantum nuclear effects in water using centroid molecular dynamics

The quantum nuclear effects are studied in water using the method of centroid molecular dynamics (CMD). The aim is the calibration of CMD implementation in LAMMPS. The calculated intramolecular energy, atoms gyration radii and radial distribution functions are shown in comparison with previous works. The work is assumed to be the step toward to solution of the discrepancy between the simulation results and the experimental data of liquid n-alkane properties in our previous works.


Introduction
Classical molecular dynamics (MD) [1][2][3] methods allow to calculate properties of systems in which quantum effects could be neglected [4][5][6][7][8][9]. The interatomic potentials for such simulations can be obtained from experimental data or ab initio density functional theory calculations [10]. The last method gives the average forces between atoms based on electronic density distributions of the different microstates of current system.
The problems with classical MD appear when the effects of zero point energy and tunnelling should be taken into account [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27]. It becomes especially important when the system consists of molecules with light atoms that oscillate with high frequencies ( ω/(k B T ) > 1). In these cases, atomic nuclei should be treated nonclassically. The examples of such systems are liquid water [13,15,17,28] and liquid hydrogen [14,18,19,21]. In the case of condensed helium quantum nuclear effects also should be considered [29]. Authors of this paper faced difficulties in calculating liquid n-alkane properties in [30] using classical MD. The reasons of such discrepancies lie in quantum impact on molecular ordering and motion in liquid which give an impact on equation of state and transport properties. The properties of hydrogen in gas hydrates also can depend on such effects [31].
Feynman's formalism of quantum mechanics is a clue to solving the problems of quantum nuclei in systems [32]. In this theory, quantum particle has an infinite amount of possible paths from one point to another [2]. The application of the idea in practice is a ring polymer that consists of P beads connected by quantum springs. It describes the Trotter decomposition of quantum Hamiltonian. The technique for improving of convergence is presented in [33] The pioneer applications of this theory in practice are path integral Monte Carlo (PIMC) and molecular dynamics (PIMD) methods firstly developed and applied in [11,12] for liquid H 2 O and D 2 O water. The most recent PIMC researches can be found in [22][23][24] correspondingly. PIMC allows to calculate static properties, but for dynamic correlation functions PIMD should be used.
In addition to classical PIMD, there are two realizations of PIMD: centroid molecular dynamics (CMD) [34] and ring polymer molecular dynamics (RPMD) [16]. The first one is applied for calculation of water dynamic properties in [13]. The CMD study of quantum effects in light and heavy water are performed in [15]. The modification of water force field for quantum simulations and simulation of nuclear effects using CMD are done in [20]. The quantum properties of liquid para-hydrogen and liquid ortho-deuterium are studied using CMD in [14], the following comparison between CMD and RPMD methods is done in [19] based on reproducing quantum time correlation functions. Also, these two methods are compared in [18] for para-hydrogen. RPMD with modified TIP4P potential is used in [28] for calculation of water properties. The research of hydrocarbons spectra is performed using RPMD and CMD methods in [35]. A lot of information about the current situation in the area of studying nuclear quantum effects can be found in reviews [25] and [27].
Most of groups use their own codes for PIMD simulations. Here, we present the available ways to carry out PIMD simulation in open-source packages. Well-known MD package LAMMPS [36] has realizations of classical PIMD and CMD within "fix" command in USER-MISC package. Another way to perform PIMD simulation is to use Python wrapper i-PI [26] which takes forces from various MD or ab initio codes and performs PIMD calculation (it also works with LAMMPS). i-PI is mostly aimed at ab initio calculations because Python can become a bottleneck when forces are calculated quickly in MD systems. Also, one can found the following open-source projects: PIMD [37] which is used by Japan Atomic Energy Agency (e.g. [38]) and RPMDrate package [39] for RPMD simulations of bimolecular reactions (e.g. [40]).
The aim of this paper is the calibration of the CMD method implemented in LAMMPS. Section 2 is dedicated to the PIMD methods that are available, the potential used and the computational details. In section 3, the calculated intramolecular energy, atoms gyration radii and radial distribution functions in comparison with previous works are presented. The results are assumed to be the step toward supplementing for the difficulties in calculating liquid n-alkane properties in [30,41].

Simulation technique 2.1. PIMD methods
The essence of the PIMD methods are the application of Feynman's formalism to molecular dynamics. This approach has been used in a lot of previous works, so the theoretical justification can be found in [2]. Here we only present the difference between the classical molecular dynamics and briefly give key ideas of different PIMD approaches.
In classical molecular dynamics a Hamiltonian H of N particles can be written as follows: where p j and r j are j-th particle linear momentum and radii, V (r 1 , . . . , r N ) is a potential that acts between the particles. In the PIMD case, a Hamiltonian become quantized which is presented via polymer rings with classical interaction. Each polymer ring represents a single atom which becomes distributed via this ring. The interatomic potential V is realized in such way that the each bead of polymer ring acts only with the corresponding bead of another polymer ring. The beads interaction inside the polymer is presented via harmonic potential that acts only between nearest neighbours. The dynamics of polymer rings system where rings oscillations imitate quantum delocalization of atoms allows to calculate quantum statistical properties of such liquids. The quantum Hamiltonian has a form where P is a number of beads (Trotter slices), p i and m ′ j are a fictitious momentum and mass correspondingly, is a potential that acts between beads with the same i of two different atoms.
There are three different approaches: classical PIMD, CMD and RPMD. All of them allow to calculate dynamical properties of systems with such Hamiltonian H q . The difference lies mostly in choice of fictitious masses m ′ j and the use of thermostats. As a result, PIMD, CMD and RPMD have different techniques for averaging of observables and obtaining autocorrelation functions for dynamic properties (can be found in [19,21,42]).
In the classical PIMD, the value used for a beads fictitious mass is arbitrary. A common choice is to use m ′ j = m j /P , which results in the mass of the entire ring-polymer being equal to the real quantum particle. The efficiency problem of simulating system with Hamiltonian H q lies in the size of MD timestep that should be less than the highest frequency in considering system. Due to this fact, very long runs should be performed to treat low-frequency modes. The solution lies in transformation of coordinates to normal mode coordinates that diagonalizes the harmonic term [2]. In this case, the fictitious masses become equal to λ j m j where λ j are eigenvalues of transformation matrix. All the modes should be thermostated in this method.
In CMD, the fictitious masses are the same as in PIMD. It acts in normal modes by default. The first (centroid) mode which is a centre of quantum polymer is uncoupled from thermostat. Such technique allows to obtain the correlation functions of quantum system from centroid mode only. The additional adiabatic parameter γ 2 is used [14,19] for fictitious masses (which become to γ 2 λ j m j , γ < 1) to separate non centroid modes from the first one. It is needed for obtaining more accurate mean forces that act on centroid and for shifting quantum frequencies off the system physical spectra [35]. Also, it allows to integrate non-centroid modes separately.
In contrast to PIMD and CMD, RPMD uses physical masses as fictitious m ′ j = m j . This choice leads to physical dynamics of all beads and no thermostats should be used. The correlation functions in this case should be averaged over all beads. The original formulation of RPMD can be found in [25,43] which is shown to be equivalent to described here [42].

Model
The authors of [28] show that interatomic potential should be fitted for quantum simulations with H q . They use TIP4P/2005 [44] water force field with anharmonic O-H bonds and fitted the parameters for reproducing equation of state. We use the designed in [28] q-TIP4P/F potential that has the following form: where intramolecular part consists of 2 anharmonic bond and harmonic angle interactions where D r , α r and k θ are bond and angle constants, r eq and θ eq are the equilibrium values of bond distance and angle.

Intermolecular part is presented via Lennard-Jones and Coulomb interactions correspondingly
where ε and σ are Lennard-Jones constants, r ij is the distance between the oxygen atoms, r mn is the distance between partial charge sites in molecules i and j (for more details see [44]). All the parameters are taken from [28].

Computational details
The starting classical configuration is a gas of the 216 water molecules. The distance between them is larger than the force field cutoff radius (9Å). The integration timestep is 0.5 fs. At the first stage, the temperature is set to 298 K and the system is compressed to the experimental density 0.997 g/cm 3 [28] for 100 ps. The second stage is a relaxation in the NVT ensemble for 100 ps (T = 298 K). The CMD simulations for the system described above are carried out in the LAMMPS package using fix CMD command from USER-MISC package. The number of Nose-Hoover thermostat chains [45] is 4 for all runs. The adiabatic parameter γ 2 is chosen to be 1/P depending on the number of beads as in [19]. The integration timesteps are 0.01 fs for P = 4, P = 8, 0.008 fs for P = 16 and 0.005 fs for P = 24, P = 32 which found to be satisfactory for water (found in [20]). The trajectory length is 40 ps for each P . For P = 32, 6 runs with different initial velocity distributions are performed to obtain better statistics for the radial distribution functions.

Intramolecular energy
The number of discretizations P should be enough to achieve the convergence of system thermodynamic properties. In figure 1, we present the dependence of intramolecular energy E intra on the inverse number of polymer beads P . The numbers show the corresponding amount of the discretizations P .
The intramolecular energy per molecule is calculated as E intra = P i=1 V intra /P and then divided by the number of molecules. It begins to converge after P = 16. This result sets the minimum number of P required for the quantum simulations. This dependence matches qualitatively the results calculated in [13] using CMD in the SPC/F water model. Totally, the intramolecular energy per molecule increases by a factor of 5.

Gyration radii
In theory, a free quantum particle is localized in the area determined by the size of de Broglie wavelength λ(T ) = h/ √ 2πmk B T . The gyration radii r 2 G = P i=1 (r i − r c ) 2 /P for free particle reproduced by polymer with P beads can be estimated as: [46]. These two formulas for λ(T ) and r 2 G provide the following relation between the gyration radii and the size of de Broglie wavelength r G = λ(T )(1 − 1/P 2 )/ √ 8π. The gyration radii for the hydrogen and oxygen atoms are obtained by averaging over atoms in the system: r 2 G = N n=1 r 2 G /N . The results are shown in figure 2. The blue, black and red open squares correspond to the different numbers of beads. The black and blue lines are theoretical estimates for r 2 G (P → ∞) and r 2 G (P = 4) respectively. The calculated gyration radii become constant value after short time of initialization and stay the same during the whole calculation. After P = 4, obtained r 2 G lie at the same distance from the predicted one. We also see that the size of quantum polymer does not depend on the number of beads. It can be explained by the fact that beads lie closer to each other due to the higher quantum spring constants which are equal to m 2 j P/(β ) 2 with the increase of P .   The theoretical estimations for free particle predict the values of λ(T )/ √ 8π that lie close to the data obtained with q-TIP4P/F potential. This result verifies the used model.

Radial distribution functions
The radial distribution functions g(r) are calculated for centroids to show the difference between classical and quantum approaches (figure 3). The black and red curves correspond to classical (CMD, P = 1) and quantum (CMD, P = 32) simulations. The red squares show data from [28].
In the quantum system, the atoms become more delocalized in the space than in the classical system. The g H-H (r) and g O-H (r) peaks are much thinner in the classical system which means that the atoms are better structured. The widths of quantum g H-H (r) and g O-H (r) first peaks match to the value of hydrogen gyration radii value (≈ 0.2Å) presented above. Both quantum and classical simulations give almost similar g O-O (r). It can be explained by the smaller r 2 G of the oxygen polymer due to its mass which leads to very similar structure in both cases.
The results from [28] obtained using classical PIMD lie on the g(r) curves calculated using CMD in this work. This coincidence indicates the fact that both classical PIMD and CMD methods give the same result.

Diffusion
The consideration of quantum nuclear effects influences upon the calculation of self-diffusion coefficient D. In the cases of liquid para-hydrogen and ortho-deuterium [14,18,19,21], it is shown that the values of D become higher for H 2 and lower for D 2 than in the classical simulations. In the quantum simulations diffusivity values are closer to the experimental data. The results calculated using CMD and RPMD techniques can differ within 15% [18,19,21].
In water, these effects have the same influence as in para-hydrogen: the diffusivity becomes higher in the quantum model [13,15,17,28]. Authors of [13,15] show about 1.4-1.5 increase of the diffusion constants for light [13] and heavy water [15] using CMD method. RPMD gives similar results in [17]. In [28], the diffusivity calculation is done using RPMD technique in q-TIP4P/F potential which is shown to give better agreement with the experimental data.
We also see the same behavior for D, but the results require additional computations for better statistics.

Conclusion
In this work, the centroid molecular dynamics (CMD) quantum simulations of liquid water are performed using q-TIP4P/F potential. The static properties such as intramolecular energy, gyration atoms radii and radial distribution functions are calculated.
The intramolecular energy dependence on the number of quantum discretizations P begin to converge after P = 16. This result sets the minimum number of P required for the quantum simulations. The same trend is shown in work [13].
The calculated quantum polymer gyration radii almost coincide with the theoretical estimates based on the size of de Brogile wavelength. The quantum delocalization of hydrogen atoms leads to less atomic ordering in the simulations which is shown in calculated radial distribution functions g H-H (r) and g O-H (r). All the g(r) curves match the results obtained in [28].
The results of this work additionally verify the CMD method implemented in LAMMPS.