.. _isotropic-label: Isotropic systems ================= .. container:: justify In this tutorial, the NMR relaxation times :math:`T_1` and :math:`T_2` are measured from a bulk polymer-water mixture using |NMRforMD|. To follow the tutorial, |MDAnalysis|, |numpy|, and |matplotlib| must be installed. .. |NMRforMD| raw:: html NMRforMD .. |MDAnalysis| raw:: html MDAnalysis .. |numpy| raw:: html numpy .. |matplotlib| raw:: html matplotlib MD system --------- .. container:: hatnote Measuring the NMR relaxation time from a bulk water-polymer mixture .. image:: ../figures/tutorials/isotropic-systems/snapshot-dark.png :class: only-dark :alt: PEG-water mixture simulated with LAMMPS - Dipolar NMR relaxation time calculation :width: 250 :align: right .. image:: ../figures/tutorials/isotropic-systems/snapshot-light.png :class: only-light :alt: PEG-water mixture simulated with LAMMPS - Dipolar NMR relaxation time calculation :width: 250 :align: right .. container:: justify The system is made of a bulk mixture of 320 :math:`\text{TIP4P}-\epsilon` water molecules and 32 :math:`\text{PEG}300` polymer molecules. The trajectory was recorded during a :math:`10\,\text{ns}` production run performed with the open source code LAMMPS in the NPT ensemble using a timestep of :math:`1\,\text{fs}`. The imposed was temperature :math:`T = 300\,^\circ\text{K}`, and the pressure :math:`p = 1\,\text{atm}`. The positions of the atoms were recorded in the *prod.xtc* file every :math:`1\,\text{ps}`. .. container:: justify If you are not familiar with LAMMPS, you can find |lammps-tutorials| here. .. |lammps-tutorials| raw:: html tutorials File preparation ---------------- .. container:: justify To access all trajectory and input files, download the *polymer-in-water* repository from Github, or simply clone it on your computer using: .. code-block:: bash git clone https://github.com/simongravelle/polymer-in-water.git .. container:: justify The dataset required to follow this tutorial is located in *raw-data/NPEG32/*. Create a MDAnalysis universe ---------------------------- .. container:: justify Open a new Python script or a new notebook, and define the path to the data files: .. code-block:: python datapath = "mypath/polymer-in-water/raw-data/NPEG32/" .. container:: justify Then, import numpy, MDAnalysis, and NMRforMD: .. code-block:: python import numpy as np import MDAnalysis as mda import nmrformd as nmrmd .. container:: justify From the trajectory files, let us create a MDAnalysis universe. Import the configuration file and the trajectory: .. code-block:: python u = mda.Universe(datapath+"init.data", datapath+"prod.xtc") u.transfer_to_memory(stop=501) .. container:: justify 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 increase its value. The figures here have been generated using the full trajectory (i.e. without the *transfer_to_memory* command). .. container:: justify The MDAnalysis universe *u* contains both topology (atoms types, masses, etc.) and trajectory (atom positions at every frame). .. container:: justify Let us extract a few information from the universe, such as number of molecules (water + PEG), timestep, and total duration: .. code-block:: python n_molecules = u.atoms.n_residues print(f"The number of molecules is {n_molecules}") .. code-block:: bw >> The number of molecules is 352 .. code-block:: python timestep = np.int32(u.trajectory.dt) print(f"The timestep is {timestep} ps") .. code-block:: bw >> The timestep is 1 ps .. code-block:: python total_time = np.int32(u.trajectory.totaltime) print(f"The total simulation time is {total_time} ps") .. code-block:: bw >> The total simulation time is 1000 ps .. container:: justify Note that in the context of MDAnalysis, the *timestep* refers to the duration between two recorded frames, which is different from the actual timestep of :math:`1\,\text{fs}` used for the LAMMPS molecular dynamics simulation. Launch the NMR analysis ----------------------- .. container:: justify Let us create 3 atom groups for respectively the hydrogen atoms of the PEG, the hydrogen atoms of the water, and all the hydrogen atoms: .. code-block:: python H_PEG = u.select_atoms("type 3 5") H_H2O = u.select_atoms("type 7") H_ALL = H_PEG + H_H2O .. container:: justify Then, let us run NMRforMD for all the hydrogen atoms: .. code-block:: python nmr_ALL = nmrmd.NMR(u, atom_group = H_ALL, neighbor_group = H_ALL, number_i=40) .. container:: justify With 'number_i = 40', only 40 randomly selected atoms within 'H_ALL' are considered for the calculation. Increase this number for better resolution, and use 'number_i = 0' to consider all the atoms. Extract the NMR spectra ----------------------- .. container:: justify Let us access the calculated value of the NMR relaxation time :math:`T_1`: .. code-block:: python T1 = np.round(nmr_ALL.T1,2) print(f"The total NMR relaxation time is T1 = {T1} s") .. code-block:: bw >> NMR relaxation time T1 = 2.53 s .. container:: justify The value you obtain may vary, depending on which hydrogen atoms were randomly selected for the calculation. .. container:: justify The full :math:`T_1` and :math:`T_2` spectra can be extracted as well as 1/nmr_ALL.R1 and 1/nmr_ALL.R2, respectively, and the corresponding frequency is given by nmr_ALL.f. .. code-block:: python R1_spectrum = nmr_ALL.R1 R2_spectrum = nmr_ALL.R2 T1_spectrum = 1/R1_spectrum T2_spectrum = 1/R2_spectrum f = nmr_ALL.f .. container:: justify The spectra :math:`T_1` and :math:`T_2` can then be plotted as a function of :math:`f` using pyplot. .. code-block:: python from matplotlib import pyplot as plt plt.loglog(f, T1_spectrum, 'o') plt.loglog(f, T2_spectrum, 's') plt.show() .. image:: ../figures/tutorials/isotropic-systems/T1-dark.png :class: only-dark :alt: NMR results obtained from the LAMMPS simulation of water .. image:: ../figures/tutorials/isotropic-systems/T1-light.png :class: only-light :alt: NMR results obtained from the LAMMPS simulation of water .. container:: figurelegend Figure: NMR relaxation times :math:`T_1` (disks) and :math:`T_2` (squares) as a function of the frequency :math:`f` for the :math:`\text{PEG-H}_2\text{O}` bulk mixture. Calculate the intra-molecular NMR relaxation -------------------------------------------- .. container:: justify Let us measuring the intra-molecular contribution to the NMR relaxation time measurement. To make the analysis easier, let us also differentiate between PEG and :math:`\text{H}_2\text{O}` molecule, and perform 2 separate analyses. .. code-block:: python nmr_H2O_intra = nmrmd.NMR(u, atom_group = H_H2O, type_analysis="intra_molecular", number_i=40) nmr_PEG_intra = nmrmd.NMR(u, atom_group = H_PEG, type_analysis="intra_molecular", number_i=40) .. container:: justify The correlation function Gij can be accessed from nmr_H2O_intra.gij[0], and the time from nmr_H2O_intra.t. .. code-block:: python t = nmr_PEG_intra.t G_intra_H2O = nmr_H2O_intra.gij[0] G_intra_PEG = nmr_PEG_intra.gij[0] .. image:: ../figures/tutorials/isotropic-systems/Gintra-dark.png :class: only-dark :alt: NMR results obtained from the LAMMPS simulation of water-PEG .. image:: ../figures/tutorials/isotropic-systems/Gintra-light.png :class: only-light :alt: NMR results obtained from the LAMMPS simulation of water-PEG .. container:: figurelegend Figure: Intra-molecular correlation function :math:`G_\text{R}` for both PEG (squares) and :math:`\text{H}_2\text{O}` (disks). .. container:: justify From the correlation functions, one can obtain the typical rotational time of the molecules. .. code-block:: python tau_rot_H2O = np.round(np.trapz(G_intra_H2O, t)/G_intra_H2O[0],2) tau_rot_PEG = np.round(np.trapz(G_intra_PEG, t)/G_intra_PEG[0],2) print(f"The rotational time of H2O is = {tau_rot_H2O} ps") print(f"The rotational time of PEG is = {tau_rot_PEG} ps") .. code-block:: bw >> The rotational time of H2O is = 6.35 ps >> The rotational time of PEG is = 8.34 ps