Best practices

Accurate NMR relaxation calculations from molecular dynamics simulations require careful attention to both the simulation protocol and the subsequent analysis. Because relaxation rates depend on molecular structure and dynamics over a broad range of timescales, they are sensitive to simulation parameters such as the force field, trajectory length, sampling frequency, simulation box size, and analysis settings. Here, some of the main factors that influence the accuracy of NMR relaxation calculations from molecular dynamics are discussed, and practical recommendations for obtaining reliable and reproducible results are provided.

Force field

The agreement between experiments and simulations is limited by the quality of the chosen force field. While some force fields show excellent agreement with NMR experimental data, for instance in simulations of water, hydrocarbons, or polymer melts [5, 10, 20], it is important to recognize that they are generally not parametrized specifically to reproduce NMR relaxation observables. Instead, they are typically optimized for selected thermodynamic, structural, and, in some cases, dynamical properties [25, 26]. These target properties may include densities, heats of vaporization, phase equilibria, solvation energies, radial distribution functions, or dynamical observables. Since NMR relaxation rates are governed by time correlation functions that reflect both equilibrium structure and molecular dynamics, the suitability of a force field for relaxation studies depends on its ability to capture the relevant motions on the corresponding timescales.

Impact of the water model

As an example, the NMR relaxation properties of bulk water were calculated from molecular dynamics simulations for three water models, \(\text{TIP4P/2005}\) [27], \(\text{SPC/E}\) [28], and \(\text{TIP3P}\) [29].

Correlation functions extracted for all three models at \(T = 300\,\text{K}\) highlight the differences between them, with \(\text{TIP3P}\) showing much faster decorrelation compared to \(\text{SPC/E}\), which itself decorrelates slightly faster than \(\text{TIP4P-2005}\) (Fig. 1, panel A). This ordering of the molecular dynamics timescales is consistent with the relative viscosities reported for these models [30] at \(T = 298\,\text{K}\): \(0.321 \, \text{mPa s}\) for \(\text{TIP3P}\), \(0.729 \, \text{mPa s}\) for \(\text{SPC/E}\), and \(0.855 \, \text{mPa s}\) for \(\text{TIP4P/2005}\). Among the models considered, the \(\text{TIP4P/2005}\) model has a viscosity closest to the experimental value of \(0.896 \, \text{mPa s}\) [31].

The NMR relaxation rate, \(R_1\), was extracted as a function of temperature. The \(\text{TIP4P/2005}\) model shows excellent agreement with experimental measurements reported by Krynicki [32] and Hindman et al. [33]. By contrast, both \(\text{SPC/E}\) and \(\text{TIP3P}\) underestimate the relaxation rate, consistent with their larger deviations in viscosity and with previous observations by Calero et al. [4].

Figure 1: A) Correlation functions extracted from molecular dynamics simulations for three water models: \(\text{TIP4P/2005}\) (blue disks), \(\text{SPC/E}\) (cyan squares), and \(\text{TIP3P}\) (green pentagons) at \(T = 300\,\text{K}\). B) NMR relaxation rate, \(R_1\), as a function of the scaled inverse temperature, \(1000 [\text{K}]/T\), for bulk water obtained from molecular dynamics simulations using the \(\text{TIP4P/2005}\), \(\text{SPC/E}\), and \(\text{TIP3P}\) models. Simulation results are compared with experimental measurements reported by Krynicki [32] and Hindman et al. [33].

The water example illustrates a more general principle: NMR relaxation provides a demanding test of molecular dynamics force fields because it is sensitive to both equilibrium structure and molecular motions over a wide range of timescales. Agreement with thermodynamic or structural properties alone does not guarantee accurate relaxation rates. Consequently, NMR relaxation calculations from MD simulations can be used not only to interpret experimental measurements, but also to assess and compare force fields based on their ability to reproduce dynamical observables.

Simulation protocol

NMR relaxation calculations are sensitive to both thermodynamic and dynamical properties. To ensure accurate simulations, the simulation protocol must be carefully chosen alongside the force field discussed above. Important aspects of the simulation protocol include the integration timestep, cutoff distances, thermostat and barostat settings, and the equilibration procedure [34, 35]. An excessively large integration timestep introduces errors in the equations of motion, insufficient sampling can bias calculated properties, and inappropriate thermostat or barostat coupling parameters can artificially affect dynamical properties. These effects can alter the structural and dynamical properties of the system, which are directly reflected in the time correlation functions used to calculate NMR relaxation rates and can therefore lead to inaccurate relaxation predictions.

Cutoff

The treatment of non-bonded interactions is an important consideration for NMR relaxation calculations, as these interactions influence both the equilibrium structure and molecular dynamics. These properties influence the time correlation functions underlying relaxation rates. To quantify this effect, the NMR relaxation rate \(R_1\) of bulk water was calculated for different Lennard-Jones cutoff distances \(r_\text{LJ}\) using the \(\text{TIP4P/2005}\) model.

The intermolecular characteristic time, \(\tau_\text{T}\), shows a clear dependence on the cutoff distance. For the smallest cutoff considered (\(r_\text{LJ} = 0.6\,\text{nm}\)), \(\tau_\text{T}\) increases by about \(10\,\%\) compared to the converged value obtained for the largest cutoff distances (Fig. 2, panel A). This reflects the weakening of intermolecular interactions caused by the truncation of the Lennard-Jones potential, which modifies the liquid structure and the molecular motions contributing to the intermolecular correlation functions. For a cutoff of \(1\,\text{nm}\), which is a commonly used value in molecular dynamics simulations, the deviation is reduced to approximately \(2\,\%\).

The variation of \(\tau_\text{T}\) affects the calculated intermolecular contribution to the relaxation rate, \(R_\text{1, T}\), which is overestimated for the smallest cutoffs (Fig. 2, panel B). These observations are consistent with previous measurements, see for instance Ref. [10].

Figure 2: a) Inter-molecular characteristic time, \(\tau_\text{T}\), as a function of the LJ cutoff, \(r_\text{LJ}\). The dashed line is a guide to the eye, indicating the value \(\tau_\text{T} = 3.78 \, \text{ps}\) obtained for the largest cutoffs. b) Inter-molecular NMR relaxation rate, \(R_\text{1, T}\), as a function of \(r_\text{LJ}\). The dashed line is \(R_\text{1, T} = 0.125 \, \text{s}^{-1}\).

Integration timestep

The integration timestep determines the numerical accuracy of the molecular dynamics trajectory. If the timestep is too large, the equations of motion are not accurately integrated, resulting in systematic errors in both structural and dynamical properties [36]. Since NMR relaxation rates are directly related to molecular motions, these integration errors can propagate into the calculated correlation functions and relaxation rates. The timestep should therefore be chosen according to established best practices for the selected force field and simulation conditions, and its adequacy should be verified through convergence testing when high accuracy is required.

For rigid water models such as \(\text{TIP4P/2005}\) combined with the SHAKE algorithm [37], a timestep of \(2\,\text{fs}\) is commonly used.

Thermostat

The thermostat controls the temperature of the simulated system [38, 39, 40, 41], but it can also influence molecular dynamics if applied too aggressively [34, 42, 43]. Strong coupling or inappropriate thermostat parameters may artificially damp or modify translational and rotational motions, leading to biased time-correlation functions and relaxation rates. For NMR relaxation calculations, it is therefore important to employ a thermostat that preserves realistic dynamics and to use coupling parameters that minimally perturb the natural motion of the system.

Simulation length

The required simulation duration depends on the quantity of interest. For the zero-frequency relaxation rate, \(R_1(0)\), the trajectory must be sufficiently long for the correlation function \(G(t)\) to decay to zero. In practice, the total simulation time should significantly exceed the longest correlation time, \(\tau_c\), of the system so that all relevant molecular motions contributing to the relaxation process are adequately sampled. For water at ambient temperature, the longest correlation time is the intermolecular time, \(\tau_\text{T} \approx 4\,\text{ps}\). In systems with slower dynamics, such as polymer melts or highly viscous liquids, \(\tau_c\) can be orders of magnitude larger, requiring substantially longer simulations to obtain a converged estimate of \(R_1(0)\).

A different consideration applies when evaluating the frequency-dependent relaxation rate, \(R_1(f)\). Because the trajectory has a finite duration \(T_\text{sim}\), its frequency resolution is approximately \(1/T_\text{sim}\). Consequently, frequencies smaller than \(1/T_\text{sim}\) cannot be resolved, irrespective of the sampling interval or analysis method. Extending the simulation is therefore the only way to improve the accessible low-frequency range of \(R_1(f)\).

Box size

NMR relaxation calculations can be sensitive to finite-size effects arising from the limited size of the simulation box [2]. To assess this effect, NMR relaxation properties were calculated for water boxes containing different numbers of molecules, \(N \in [25,\,6100]\), corresponding to equilibrium box lengths of \(L \in [0.9,\,5.6]\,\text{nm}\).

The influence of the box size is not immediately apparent from the intermolecular correlation function, \(G_{ij, \text{T}}(t)\), over the main decay region. However, finite-size effects become significant at long times, where \(G_{ij, \text{T}}(t)\) follows the hydrodynamic decay \(G_{ij, \text{T}}(t) \sim t^{-3/2}\) (Fig. 3, panel A). In small simulation boxes, this long-time contribution is artificially truncated, affecting the integral of the correlation function and therefore the intermolecular contribution to the relaxation rate, \(R_{1,\text{T}}\).

The effect on \(R_{1,\text{T}}\) remains significant even for the largest boxes considered in this study. Our results indicate that even for the largest boxes with \(N = 6100\), \(R_{1,\text{T}}\) may not be fully converged (Fig. 3, panel B). It should be noted that the intramolecular contribution, \(R_{1,\text{T}}\), which dominates the total relaxation rate of water at ambient temperature, is only weakly affected by the box size.

Figure 3: A) Inter-molecular correlation function \(G_{ij, \text{T}}\) for different numbers of molecules (see the legend). B) Inter-molecular NMR relaxation rate, \(R_\text{1, T}\), as a function of the number of molecules \(N\) for a bulk water system. For the smallest systems, results were averaged from up to 20 independent simulations and the error bar is calculated from the standard deviation. The dashed line is \(R_\text{1, T} = 0.125 \, \text{s}^{-1}\).

Trajectory output frequency

The trajectory output frequency determines the temporal resolution of the analysis and therefore the fastest molecular motions that can be resolved. The sampling interval, \(\Delta t\), must be significantly smaller than the shortest relevant correlation time of the system. Otherwise, the rapid decay of the correlation function \(G_{ij}(t)\) is undersampled (Fig. 4, panel A), leading to an inaccurate estimate of the characteristic times, \(\tau\), and consequently of the relaxation rates. When the relevant correlation times are not known a priori, the appropriate value of \(\Delta t\) can be determined through convergence testing. A smaller sampling interval, however, increases the trajectory size and the computational cost of the subsequent analysis.

As an illustration, the intermolecular NMR relaxation rate, \(R_{1,T}\), of bulk water was calculated for sampling intervals ranging from \(\Delta t = 20\,\text{fs}\) to \(1.2\,\text{ps}\). Sampling intervals larger than approximately \(100\,\text{fs}\) lead to an overestimation of \(R_{1,T}\). This behavior is consistent with an insufficient temporal resolution of the fast molecular motions contributing to the intermolecular correlation function.

Figure 4: A) Intermolecular correlation functions extracted from molecular dynamics simulations using two sampling intervals, \(\Delta t = 20\,\text{fs}\) and \(\Delta t = 1.3\,\text{ps}\). The characteristic time \(\tau_\text{T} \approx 3.8 ~\text{ps}\) is indicated by the vertical dashed line. B) Intermolecular NMR relaxation rate, \(R_{1,T}\), as a function of the sampling interval \(\Delta t\). The dashed line is \(R_\text{1, T} = 0.125 \, \text{s}^{-1}\).

Analysis parameters

The accuracy of NMR relaxation calculations depends not only on the quality of the molecular dynamics simulation, but also on the parameters used during the analysis. While these parameters do not alter the underlying trajectory, they can affect the statistical uncertainty of the results or determine which physical contributions are included in the calculation. Their influence should therefore be assessed through appropriate convergence tests.

Number of reference atoms

The parameter number_i controls how many reference atoms are randomly sampled during the calculation. Because the selection is stochastic, results will vary slightly between runs when number_i > 0. The statistical uncertainty decreases as number_i increases, and setting number_i = 0 includes all eligible atoms, providing the most accurate result at the highest computational cost. In practice, convergence should be verified by repeating the calculation with increasing values of number_i until the relaxation rates stabilize.

Cross-species interactions

When neighbor_group is not specified, intermolecular contributions are computed only between atoms belonging to the same chemical species. In a mixture such as polymer–water, this means that water–polymer cross-interactions are excluded from the relaxation calculation. If cross-species contributions are expected to be significant, the appropriate neighbor_group must be set explicitly. Neglecting these contributions may lead to an underestimation of the total intermolecular relaxation rate.

Table 1 Summary of the main factors affecting the accuracy of NMR relaxation calculations from molecular dynamics.

Parameter	If not properly chosen	Consequence
Force field	Incorrect description of molecular structure or dynamics	Systematic deviations in relaxation rates due to inaccurate structural properties, molecular motions, or correlation times.
Simulation protocol	Inappropriate timestep, cutoff distances, thermostat/barostat settings, or equilibration procedure	Errors in structural and dynamical properties, leading to inaccurate correlation functions and relaxation rates.
Simulation length	Trajectory duration insufficient for correlation functions to fully decay	Incomplete convergence of relaxation rates, especially in the low-frequency limit.
Box size	Simulation box too small, causing finite-size effects	Truncation of long-time intermolecular correlations and underestimation of intermolecular relaxation contributions.
Trajectory output frequency	Sampling interval too large compared with relevant correlation times	Fast molecular motions are undersampled, causing errors in correlation functions and relaxation rates.
Number of reference atoms	Too few atoms sampled during the calculation	Increased statistical uncertainty and reduced precision of the calculated relaxation rates.
Cross-species interactions	Interactions between different chemical species are not included	Underestimation of intermolecular relaxation contributions in mixtures.