Bulk water#

The NMR relaxation time \(T_1\) of water is calculated

Water molecules simulated with lammps - NMR relaxation time calculation Water molecules simulated with lammps - NMR relaxation time calculation

In this tutorial, the NMR relaxation time \(T_1\) of water is measured using NMRforMD. MDAnalysis, numpy, and matplotlib and NMRforMD must be installed.

The system is made of 398 TIP4P water molecules simulated in the NPT ensemble with LAMMPS at a temperature of 300 K. The total duration of the simulation is 10,ns, and the timestep is 2 fs. You can access the input files in this repository, which you can use to create larger system or longer trajectory. If you are not familiar with LAMMPS, you can find tutorials here.

File preparation#

Either download the files from the Github repository, or clone the NMRforMD repository:

git clone git@github.com:simongravelle/nmrformd.git

The datasets are located in ‘examples/raw-data/bulk-water/N398/’.

Open a Python script or a Jupyter notebook, and start by defining a path to the data files. In my case, since I am working from ‘examples/analyse-scripts/bulk-water/’, its simply:

datapath = "../../raw-data/bulk-water/N398/"

Import the libraries#

Import numpy, MDAnalysis, and NMRforMD:

import numpy as np
import MDAnalysis as mda
import nmrformd as nmrmd

Create a MDAnalysis universe#

Import the configuration file and the trajectory:

u = mda.Universe(datapath+"topology.data", datapath+"traj.xtc")
u.transfer_to_memory(stop=501)

The u.transfer_to_memory(stop=501), is optional, it only serve to reduce the number of frames, and therefore reduce the duration of the calculation. Feel free to remove it, or change its value.

Note : the figures here have been generated using the full trajectory (i.e. without u.transfer_to_memory(stop=501)), but its takes a few minutes to complete.

The MDAnalysis universe u contains both topology (atoms types, masses, etc.) and trajectory (atom positions at every frame).

Let us extract a few information from the universe, such as number of molecules, timestep, and total duration:

n_molecules = u.atoms.n_residues
print(f"The number of water molecules is {n_molecules}")

>> The number of water molecules is 398

timestep = np.int32(u.trajectory.dt)
print(f"The timestep is {timestep} ps")

>> The timestep is 1 ps

total_time = np.int32(u.trajectory.totaltime)
print(f"The total simulation time is {total_time} ps")

>> The total simulation time is 500 ps

Run NMRforMD#

Let us isolate a group of atoms containing all the hydrogen atoms (i.e. atoms of type 2) of the system:

group_i = u.select_atoms("type 2")

Then, let us run NMRforMD, using the same group as i and j types:

nmr_result = nmrmd.NMR(u, group_i, number_i=40)

With ‘number_i = 40’, only 40 randomly selected atoms within ‘group_i’ are considered for the calculation. Increase this number for better resolution. Use ‘number_i = 0’ to consider all the atoms.

Extract results#

Let us access the calculated value of the NMR relaxation time T1:

T1 = np.round(nmr_result.T1,2)
print(f"NMR relaxation time T1 = {T1} s")

>> NMR relaxation time T1 = 3.08 s

The value you get may vary a little, depending on which hydrogen atoms were randomly selected by NMRforMD.

The T1 spectrum can be extracted as 1/nmr_result.R1 (i.e. the invert of R1), and the corresponding frequency is given by nmr_result.f. Let up plot T1 as a function of f:

NMR results obtained from the LAMMPS simulation of water NMR results obtained from the LAMMPS simulation of water

The correlation function Gij can be accessed from nmr_result.gij[0], and the time from nmr_result.t. Let us plot Gij as a function of t:

NMR results obtained from the LAMMPS simulation of water NMR results obtained from the LAMMPS simulation of water

Intra vs inter-molecular#

Let us calculate the intra-molecular and inter-molecular contributions to the NMR spectrum R1 by calling NMRforMD twice:

nmr_result_intra = nmrmd.NMR(u, group_i, type_analysis="intra_molecular", number_i=40)
nmr_result_inter = nmrmd.NMR(u, group_i, type_analysis="inter_molecular", number_i=20)

Note that the intra_molecular contribution is always noisier than the inter_molecular, which is why more atoms were included in the analysis. We can plot both intra-molecular and inter-molecular contributions separately:

NMR results obtained from the LAMMPS simulation of water NMR results obtained from the LAMMPS simulation of water

It appears that the intra-molecular contribution is the dominant one, which is expected for bulk water. We can also look at the correlation functions:

NMR results obtained from the LAMMPS simulation of water NMR results obtained from the LAMMPS simulation of water

Another expected result: the inter-molecular contribution, which is typically associated with the translational motion of the molecules, has longer characteristic times than the intra-contribution, which is typically associated to the rotation of the molecules.