Refinement of thermostated molecular dynamics using backward error analysis

MD Simulation of Water Using a Rigid Body Description Requires a Small Time Step to Ensure Equipartition

In simulations of aqueous systems, it is common to freeze the bond vibration and angle bending modes in water to allow for a longer time step δt for integrating the equations of motion. Thus, δt = 2 fs is often used in simulating rigid models of water. We simulate the SPC/E model of water using δt from 0.5 to 3.0 fs and up to 4 fs using hydrogen mass repartitioning. In these simulations, we find that for all but δt = 0.5 fs, equipartition is not obtained between translational and rotational modes, with the rotational modes exhibiting a lower temperature than the translation modes. To probe the reasons for the lack of equipartition, we study the autocorrelation of the translational velocity of the center of mass and the angular velocity of the rigid water molecule, respectively. We find that the rotational relaxation occurs on a timescale comparable to vibrational periods, calling into question the original motivations for freezing the vibrations. Furthermore, a time step with δt ≥ 1 fs is not able to capture accurately the fast rotational relaxation, which reveals its impact as an effective slowing-down of rotational relaxation. The fluctuation-dissipation relation then leads to the conclusion that the rotational temperature should be cooler for δt greater than the reference value of 0.5 fs. Consideration of fluctuation-dissipation in equilibrium molecular dynamics simulations also emphasizes the need to capture the temporal evolution of fluctuations with fidelity and the role of δt in this regard. The time step also influences the solution thermodynamic properties: both the mean system potential energies and the excess entropy of hydration of a soft repulsive cavity are sensitive to δt.

Molecular Dynamics with Very Large Time Steps for the Calculation of Solvation Free Energies

In multiple time scale molecular dynamics, the use of isokinetic constraints along with massive thermostatting has enabled the adoption of very large integration steps, well beyond the limits imposed by resonance artifacts in standard algorithms. In this work, we present two new contributions to this topic. First, we investigate the velocity distribution and the temperature-kinetic energy relationship associated with the isokinetic Nosé-Hoover family of methods, showing how they depend on the number of thermostats attached to each atomic degree of freedom. Second, we investigate the performance of these methods in the calculation of solvation free energies, the determination of which is often key for understanding the partition of a chemical species among distinct environments. We show how one can extract this property from canonical (constant-NVT) simulations and compare the result to experimental data obtained at a specific pressure. Finally, we demonstrate that large time steps can, in fact, be used to improve the efficiency of these calculations and that attaching multiple thermostats per degree of freedom is beneficial for effectively exploring the configurational space of a molecular system.