2.7. MD simulations of liquid water under periodic boundary conditions

MD simulations of liquid water are commonly performed using periodic spatial boundary conditions and keeping the temperature and pressure constant. The algorithm and parameters of the temperature and pressure coupling were mentioned above. However, when testing total energy conservation, the system is not coupled to a heat bath and a pressure bath, and the center of mass motion of the system is not removed during the simulation (water_liq_ener). The nonbonded interactions were calculated using a single cut‐off radius of 1.4 nm. Minimum‐image periodic boundary conditions were applied. Configurations were saved about every 0.02 ps. These 32 simulations were performed for 20 ps (Table ​(Table3).3). The simulations applying the cf‐MTS algorithm were not stable, so their results were omitted from Table ​Table33.

Energy conservation in 20 ps MD simulations of liquid water, 1000 molecules in a cubic periodic box using the minimum‐image convention to calculate the forces, without coupling to temperature and pressure baths, and without removal of center of mass translation, using no MTS algorithm or the if‐MTS algorithm for the flexible SPC/F model or three distance constraints for the rigid SPC model, as function of the MD time step Δt and the number of time steps n _hf at which only the bond‐stretching and bond‐angle forces (SPC/F) are evaluated and integrated

Note: Relative geometric precision of the distance constraints (SPC): tol _DC = 10⁻⁶. Nonbonded interaction cut‐off radius R _cp = R _cl = 1.4 nm. Outside a sphere of radius R _RF = R _cl, a homogeneous continuum dielectric with ε _RF = 78.5, κ _RF = 0 and ε _cs = 1 is assumed to be present. E _int: Bond‐stretching and bond‐angle bending energy. ΔE _int: Fluctuation of E _int. E _tot: Total energy. ΔE _tot: Fluctuation of E _tot. E _tot ^drift: Total energy drift. ΔE _tot ^drift: Fluctuation around E _tot ^drift . E _kin: Kinetic energy. ΔE _kin: Fluctuation of E _kin. E _kin ^drift: Kinetic energy drift. ΔE _kin ^drift: Fluctuation around E _kin ^drift. All energies given per molecule. T _tr: Temperature of the translational degrees of freedom. T _ir: Temperature of the internal and rotational degrees of freedom. All values are averages calculated from trajectory structures separated by approximately 0.02 ps. Energies in kJ mol⁻¹. Energy drifts in kJ mol⁻¹ ps⁻¹.

In the simulations to evaluate the properties of liquid water (water_liq_prop) the nonbonded interactions were also calculated using a single cut‐off radius of 1.4 nm. Outside the cut‐off radius a reaction‐field approximation 40 , 41 with a relative dielectric permittivity of 78.5 was used. Minimum‐image periodic boundary conditions were applied. The temperature and pressure were kept constant using the weak‐coupling algorithm. In the simulations using the flexible SPC/F water model, the motions of the translational and rotational plus internal degrees of freedom of the molecules were separately coupled to the heat bath. These simulations were performed for 1 ns. Time steps Δt of 0.1, 0.2, 0.5, 1.0, and 2.0 fs were tested using n _hf values of 1, 3 and 5, so for Δt′ values of 0.3, 0.6, 1.5, 3.0, and 6.0 fs in case n _hf = 3, and for Δt′ values of 0.5, 1.0, 2.5, 5.0 and 10.0 fs in case n _hf = 5. In the simulations using constraints, longer time steps Δt = 5.0 and 10.0 fs were also used. When constraints were applied, tol _DC = 10⁻⁶ was used. Translational motion of the center of mass of the system was removed every 2 ps. Configurations were saved about every 0.1 ps and were used to analyze various properties of liquid water as function of the type of MTS algorithm, time‐step size and whether constraints are applied. In total 32 simulations were performed (Table ​(Table44).

Selected properties of liquid water from 1 ns MD simulations of 1000 molecules in a cubic periodic box using the minimum‐image convention to calculate the forces, with coupling to temperature (τ _T = 0.1 ps, T _ref = 298.15 K) and pressure (τ _p = 0.5 ps, p _ref = 1 atm) baths, using no MTS algorithm or two different (if and cf) MTS algorithms for the flexible SPC/F model or three distance constraints for the rigid SPC model, as function of the MD time step Δt and the number of time steps n _hf at which only the bond‐stretching and bond‐angle bending forces (SPC/F) are evaluated and integrated

Note: Removal of center of mass translation every 2 ps. Relative geometric precision of the constraints (SPC): tol _DC = 10⁻⁶. Nonbonded interaction cut‐off radius R _cp = R _cl = 1.4 nm. Outside a sphere of radius R _RF = R _cl, a homogeneous continuum dielectric with ε _RF = 78.5, κ _RF = 0 and ε _cs = 1 is assumed to be present. T _tr: Temperature of the translational degrees of freedom. T _ir: Temperature of the internal and rotational degrees of freedom. ρ: Mass density. ΔH _vap: Heat of vaporization. E _pot: Potential energy per molecule. E _int: Bond‐stretching and bond‐angle bending energy per molecule. <D>: Average molecular diffusion coefficient. <τ ₂ ^OH>: Average molecular rotational correlation time of an O‐H bond. <τ ^H‐bond>: Average lifetime of water–water hydrogen bonds. All values are averages calculated from trajectory structures separated by approximately 0.1 ps.