MDA diffusisivty vs. LAMMPS diffusivity

[Margaret Brooks]

Hello,

I'm working on a computational research project where we need to find the diffusivity of our LAMMPS-modeled protein; first we used LAMMPS to find to find the MSD's and calculated the diffusivity from there, but we felt that MDAnalysis would yield more accurate results. However, when I wrote the code to find the diffusivity using MDAnalysis it did not line up well with the values we got from LAMMPS, and for simulations of more polymers in one box the number went up instead of down (as we expected it to and as it did in the case with LAMMPS). Even if the numbers don't completely match up, I don't understand why they would trend in different directions. 

Notably, we want the diffusivity based on the center of mass of our polymer, so I'm passing the LAMMPS output to MDA, using that find the center of mass, then using that data to build a new MDA universe that holds only the centers of mass of each polymer. 

Any ideas would be helpful! I'm feeling pretty stuck on this problem.   

[Micaela Matta]

Hi there,

I'd like a few clarifications from you regarding your simulation:

- what trajectory format are you using? if you are running the diffusivity calculation in lammps on the fly you are probably using higher precision than what you are printing in the trajectory file you then feed to MDAnalysis.

- how is periodic box information accounted for? when finding the centre of mass in MDAnalysis what options did you use?

Thanks,

Micaela

[Margaret Brooks]

HI, sorry for the delay; I thought I had replied but it appears it didn't go through. I'll answer your questions to the best of my ability:

- We use a custom LAMMPS dump that outputs atom id, type, and the unwrapped coordinates of each particle

- Like I mentioned above, we're using unwrapped coordinates. Our box is 400x400x400 nm. I found the center of mass by averaging the positions of each bead in the polymer as they are all identical. When I stopped running that bit of the code, the numbers began trending the way we expected them to, so the problem is definitely with how I'm calculating center of mass. 

Here's the bit of my code that's finding CoM from the first MDA universe and then passing it the a second: