MSD is not the same as the one calculated by lammps

[yann claveau]

Hello everyone,

I am still testing MDA on a box filled with 512 TIP4P/2005 water molecules. positions are ouput every 1ps, the run last 200 ps. The density is converged = 0.996, RDF is OK.

I tried to plot the MSD vs time using the two msd classes :

MDAnalysis.analysis.msd.EinsteinMSD

and

MDAnalysis.analysis.waterdynamics.MeanSquareDisplacement

They gave me similar results.

import numpy as np
import MDAnalysis as mda
import MDAnalysis.analysis.msd as msd

dump      ="wat_unwraped.lammpsdump"
lammpsData="wat.data"

u = mda.Universe(lammpsData,dump,dt=1,memory=True)

MSD = msd.EinsteinMSD(u, select='type 2', msd_type='xyz', fft=True)
MSD.run()

gives me a sel-diffusion coef = slope/(2*n) = 130.95 A²/ps

while it should be ~0.23 A²/ps

The MSD output by lammps gives me a correct D~0.22A²/ps

So I tried to calculate the MSD by hand :

O = u.select_atoms("type 2")

n_frames = u.trajectory.n_frames
n_at     = O.n_atoms

#  reference position at t=0
for ts in u.trajectory[0]:
    ref=O.ts.positions 

# displacement vector (dx, dy, dz) = (x,y,z) - (x0,y0,z0)

displacement = np.zeros(((n_frames),(n_at),3))
i = 0
for ts in u.trajectory:
    displacement[i]=O.ts.positions - ref
    i+=1

MSD=displacement*displacement # dx*dx, dy*dy, dz*dz
MSD=np.sum(MSD,axis=2)        # dx² + dy² + dz²
MSD=np.mean(MSD,axis=1)       # average over atoms

tf   = n_frames
time = np.linspace(0, tf, n_frames)

fit = np.polyfit(time, MSD, 1)
y   = np.poly1d(fit)

y_t = y(time)

fig, ax = plt.subplots()
plt.plot(time, MSD)
plt.plot(time, y_t)
ax.set_xlabel("time ns")
ax.set_ylabel("MSD A²")
print("diff coef = slope/(2*n) =", y[1]/6 )

which gives D = 130 A²/ps.

I have to admit that I am at loss...

[Hugo Macdermott-Opeskin]

Hi Yann,

I suspect you have run into a problem that we have with the MSD code (or rather a lack of something in the Transformations code). 
For the MSD to be calculated correctly by the code in MDA, you need to use unwrapped coordinates. That is, any time a molecule passes through the periodic boundary, it must be modified so that it does not wrap onto the other side of the box. In GROMACS this is achieved with `gmx trjconf -pbc nojump`. I am far from a LAMMPs expert, but perhaps something on this page will help you out. There will hopefully be a way for you to post-process your trajectory with LAMMPS to put the coordinates in the right convention.

Without using an unwrapped trajectory, the water molecules appear to the MSD code to rapidly oscillate in displacement if they pass through the PBC eg from [-boxlength, -boxlength, -boxlength] to [+boxlength, +boxlength, +boxlength] in the worst case.

We/I am are working on implementing the correct pbc unwrapping transformation, but it's a bit involved so progress has been slow.

Perhaps I will add a much more noticeable warning on the MSD docs.

Thanks

Hugo.

[yann claveau]

Thanks for your answer Hugo,

Actually I am post processing unwrapped coordinates, so it is not the problem here.

I've found something that can be the beginning of the answer :

I've tried to calculate the MSD using Ovito on the lammps dump. It gives the correct MSD.

So I've done the following :

1 - Read dump with MDA

2 - write coordinates in xyz using MDA

3 - read xyz with ovito

4 - calculate MSD

When I read the lammps dump with ovito I got the correct MSD. But after 1,2,3,4 , I obtain the same wrong MSD with Ovito than when using MDA (130 A²/ps). So I think the problem is the reader of MDA. Is it possible that there is a problem in the atom's indexes ?

Yann

[Hugo Macdermott-Opeskin]

Hi Yann, 

Would you be able to send me the original DCD  and topology I can have a look?

Cheers

[yann claveau]

I've sent you a private message with the files, let me know if it did not work.

[yann claveau]

I think I found the problem. My lammps dump coordinates are absolute coordinates while MDA read as reduced ones.

In the dump file, it is specified with xu, yu, zu for absolute, unwrap coordinate, it is xs, ys, zs for reduced coordinates.

It could be interesting to modify the lammps parser to read this line and to automatically detect the coordinate types.

[Hugo Macdermott-Opeskin]

Thanks Yann for your hard work investigating. 

I'll have a look through your files and confirm. 

I think this may be related to a known problem in the LAMMPSDUMPParser. 

I'll investigate and fix as soon as possible.

Cheers

[Hugo Macdermott-Opeskin]

I can confirm that the unscaled coordinates are not being parsed correctly. 

I will work on a fix ASAP.

[yann claveau]

Hello all,

For now, I have defined a transformation:

def scaled_to_unscaled():
    def wrapped(ts):
        """Handles the actual Timestep"""
        ts.positions = ts.positions/ ts.dimensions[0:3]
        return ts
    return wrapped

u = mda.Universe(lammpsData, dump, dt=1, transformation=scaled_to_unscaled, memory=True)

 It seems to work correctly.

Cheers,

Yann

[Hugo Macdermott-Opeskin]

Hi Yann,

We have a hotfix that should be ready very soon. 

https://github.com/MDAnalysis/mdanalysis/pull/3360

Cheers