Calculation of diffusion coefficient based on msd

[lauraf244]

Hello!
I'm trying to reproduce a diffusion coefficient based on msd of the center of mass of my solvent but i think i have problems with wrapped coordinates. The positions of the COM oscillates with time going up and down even though i'm using a lmp_unwrapped.dcd file from LAMMPS calculations. How can i solve this iussue?
Thank you

[orbeckst]

I don't know how LAMMPS does the "unwrapping".

nojump transformation

In MDAnalysis (and GROMACS) the required transformation is called "nojump". There's a trajectory transformation transformations.NoJump that could do it:

import MDAnalysis as mda
import MDAnalysis.transformations
from MDAnalysis.analysis import msd

# load the data with bond information
u = mda.Universe(PSF, "lmp.dcd")

# add the NoJump transformation: for each timestep, coordinates will be changed
# so that atom trajectories are continuous
u.trajectory.add_transformations(MDAnalysis.transformations.nojump.NoJump())

# optional: write nojump trajectory for visualization or later use
u.atoms.write("nojump.dcd", frames="all")

Then do the MSD analysis with the transformed trajectory (eg for water oxygen atoms)

MSD = msd.EinsteinMSD(u, select='resname SOL and name OW', msd_type='xyz', fft=True)
MSD.run()

unwrapping (making molecules whole)

Note that the molecules must be whole to remain whole, which we can do with atoms.unwrap() if you want to do this manually for a frame

# make molecules (=fragments, i.e., bonded atoms) whole in the first frame;
# (reference=None is more robust than "com" or "cog" in my experience, and it should
# not matter for the use case of analyzing diffusion)
u.atoms.unwrap(compound="fragments", reference=None)

# save first frame
u.atoms.write("unwrapped.pdb")

or if you need to do this every frame, you can add the transformations.unwrap() transformation together with nojump:

transformations = [MDAnalysis.transformations.unwrap(u.atoms), 
                   MDAnalysis.transformations.NoJump()]
u.trajectory.add_transformations(*transformations)

# write trajectory (rather slow, sorry)
u.atoms.write("nojump.dcd", frames="all")

The unwrap() code is very slow so at the moment I'd recommend using another tool such as GROMACS trjconv (for any format that GROMACS can read).