A refined ring polymer contraction scheme for systems with electrostatic interactions
Graphical abstract
Splitting the Coulomb potential into short- and long-range parts helps to speed up path integral simulations of liquid water.
Introduction
There is considerable current interest in the inclusion of quantum mechanical (zero-point energy and tunnelling) effects in the simulation of condensed phase systems. One of the most appealing approaches to this problem is based on imaginary time path integration [1], and exploits the isomorphism between the quantum statistical mechanics of a single particle and the classical statistical mechanics of a fictitious ring polymer [2]. Static equilibrium properties can be calculated exactly from this isomorphism, using either path integral Monte Carlo (PIMC) [3] or path integral molecular dynamics (PIMD) [4] techniques. Dynamical properties can also be obtained from the isomorphism, within the centroid molecular dynamics (CMD) [5], [6] and ring polymer molecular dynamics (RPMD) [7], [8] approximations. The computational effort required by these methods scales linearly with the number of beads (n) in the ring polymer, where n > ℏωmax/kBT is determined by the highest frequency (ωmax) present in the problem. When high frequencies and/or low temperatures are involved this makes the quantum calculation considerably more demanding than a classical molecular dynamics simulation.
In a recent paper we have shown that, for systems in which the potential can be decomposed into a sum of intramolecular and intermolecular contributions, the computational effort required by path integral methods can be significantly reduced [9]. This is because the long-range intermolecular forces, which dominate the computational effort for large systems, typically vary only slightly over the length scale of a ring polymer, and they may therefore be evaluated on a contracted ring polymer with fewer than the full n beads [9] (or equivalently on a lower-order Fourier decomposition of the imaginary time path [10]). The results obtained with this scheme were found to be quite encouraging, with errors of at most a few percent over a wide range of static and dynamic properties when compared to fully converged path integral simulations of a flexible water model [9]. However, as we pointed out in our paper, the separation of the potential into intramolecular and intermolecular contributions is somewhat naive, as intermolecular (Lennard-Jones and electrostatic) potentials also contain some rapidly-varying short-range components.
The purpose of the present article is to show that it is possible to refine the ring polymer contraction scheme for systems with electrostatic interactions by separating the rapidly- and slowly-varying parts of the Coulomb potential. The importance of long-range electrostatic interactions is well established for a wide variety of complex systems ranging from bulk water to protease structures and DNA [11], [12], [13]. Neglecting these interactions by short-range truncation of the Coulomb potential has been shown in many cases to give spurious results and noticeable deviations from the correct behaviour [14], [15]. Although a variety of approaches to treating long-range electrostatic interactions have been developed [16], one of the most commonly used is Ewald summation, in which the electrostatic interactions are split into a short-range part, summed in real space, and a long-range part which can be evaluated efficiently in reciprocal space [16]. The conventional Ewald approach scales as O(N3/2) while particle mesh Ewald (PME) [17] and its variants scale as O(N log N) with respect to the number of charges (N) present in the system. Hence, in practice, even for moderate system sizes, the long-range Coulomb interactions dominate the calculation compared with the short-range interactions which scale as O(N).
In this Letter we show that, rather than performing a separate Ewald sum for each ring polymer bead in a path integral simulation, one can evaluate a single Ewald sum at the centroid of the ring polymer and cheaply evaluate a short-range correction to the forces acting on each bead. This idea is developed and then illustrated with an application to the same flexible simple point charge (SPC/F) model [18] of liquid water that we considered in our earlier study [9]. The error in the forces obtained from the method is well controlled at a pre-determined level and both the static and dynamic properties converge quickly to their correct values as the length scale that is used to separate the long and short-range forces is increased. For large systems in which the calculation of the long-range electrostatic interactions dominates, this enables path integral simulations to be performed with less than twice the computational effort of classical molecular dynamics simulations.
Section snippets
Theory
The Coulomb interaction energy between two charges, qi and qj, located at positions ri and rj is given in atomic units bywhere rij = ∣ri − rj∣. In a path integral simulation the classical particles i and j are replaced by n bead ring polymers and the interaction is summed over the beads [2], [3], [4],where . This makes the calculation of the potential energy and forces n times more expensive than in a classical molecular dynamics
Application to liquid water
In order to test the above theory, we have performed some PIMD and RPMD simulations of the same flexible simple point charge (SPC/F) water model [18] that we used to test our original ring polymer contraction scheme [9]. These simulations were performed with 216 water molecules in a periodically replicated simulation box at a density of 0.997 g cm−3 and a temperature of 298 K. The intramolecular and oxygen–oxygen Lennard-Jones potentials were each evaluated using n = 32 ring polymer beads. The
Results and discussion
As we discussed in our earlier paper [9], the average intermolecular potential and (virial) kinetic energies 〈Vinter〉 and 〈Tinter〉 provide a particularly severe test of the ring polymer contraction scheme. 〈Vinter〉 provides a direct probe of the accuracy of the approximation to the intermolecular potential and 〈Tinter〉 provides a direct probe of the approximation to the intermolecular forces. Table 1 lists the values of these energies computed using the present contraction scheme with various
Conclusion
In this Letter, we have presented a method to reduce the computational costs involved in evaluating long-range electrostatic interactions in path integral simulations. We have shown that evaluating a single Ewald sum at the ring polymer centroid and a short-range correction between the beads within the first coordination shell (σ = 3 Å) is all that is needed to obtain the vast majority of the quantum mechanical effects in the static and dynamic properties of the SPC/F water model. Including the
Acknowledgements
This work was supported by the US Office of Naval Research under Contract No. N000140510460 and by the UK Engineering and Physical Sciences Research Council under Grant No. E01741X.
References (28)
- et al.
Chem. Phys. Lett.
(1991) - et al.
Biophys. J.
(2000) - et al.
Quantum Mechanics and Path Integrals
(1965) - et al.
J. Chem. Phys.
(1981) J. Chem. Phys.
(1979)- et al.
J. Chem. Phys.
(1984) - et al.
J. Chem. Phys.
(1993) - et al.
J. Chem. Phys.
(1999) - et al.
J. Chem. Phys.
(2004) - et al.
J. Chem. Phys.
(2006)
J. Chem. Phys.
J. Chem. Phys.
J. Chem. Phys.
J. Chem. Phys.
Cited by (94)
i-PI 2.0: A universal force engine for advanced molecular simulations
2019, Computer Physics CommunicationsSimulations of disordered matter in 3D with the morphological autoregressive protocol (MAP) and convolutional neural networks
2024, Journal of Chemical PhysicsModeling nuclear quantum effects on long-range electrostatics in nonuniform fluids
2023, Journal of Chemical PhysicsStatistical mechanics: Theory and molecular simulation
2023, Statistical Mechanics: Theory and Molecular SimulationReal-Time Dynamics and Detailed Balance in Ring Polymer Surface Hopping: The Impact of Frustrated Hops
2023, Journal of Physical Chemistry Letters