Computation of MSD

LAMMPS LAMMPS General Discussion
1

Hello, I am currently working on polymer system where I am using “compute MSD” module using the dump file ehich was generated from the simulation run of 100 ns, here is my dump file generation script:-

4. PRODUCTION RUN (50 ns @ 2 fs)

NVT at volume from NPT

==========================================

reset_timestep 0
timestep 2.0 # 2 fs, as requested (watch stability!)

Use NVT for production at fixed volume from NPT

fix myprod all nvt temp 298.0 298.0 100.0

Dumps and restarts

dump myDump all atom 10000 dump.lammpstrj
dump_modify myDump sort id

Write restart files so you can resume if something happens

restart 5000000 restart1.lmp restart2.lmp

thermo 5000

50 ns = 50,000 ps = 50,000,000 fs

100 ns = 100,000 ps = 100,000,000 fs

At 2 fs/step → 25,000,000 steps

At 2 fs/step → 50,000,000 steps

run 50000000 # 100 ns @ 2 fs

After this I created another script of MSD where i used the dump.lammpstrj to compute for msd where the part of script is :-
compute 1 li msd com yes average yes
compute 2 emim msd com yes average yes
compute 3 tfsi msd com yes average yes
fix msdout all ave/time 1 1 10000 c_1[1] file msd_li.dat
fix msdout all ave/time 2 2 10000 c_2[1] file msd_emim.dat
fix msdout all ave/time 3 3 10000 c_3[1] file msd_tfsi.dat
rerun dump.lammpstrj first dump x y z box yes

for the reference the dump.lammpstrj file looks like this (first few lines):-
ITEM: TIMESTEP
419000
ITEM: NUMBER OF ATOMS
10300
ITEM: BOX BOUNDS pp pp pp
1.1403265252281745e+01 7.3894734747722907e+01
1.1390764252281752e+01 7.3882233747722921e+01
1.1244764252281744e+01 7.3736233747722906e+01
ITEM: ATOMS id type xs ys zs
1 6 0.775772 0.86388 0.272189
2 14 0.777518 0.855162 0.287326
3 14 0.785218 0.855428 0.261174
4 14 0.782322 0.880022 0.274059
5 3 0.752195 0.865121 0.264957
6 11 0.746612 0.850849 0.258994
7 11 0.741713 0.871879 0.276264
8 16 0.751356 0.880404 0.24674
9 7 0.762333 0.879665 0.227716
10 15 0.773626 0.866266 0.224697

But my problem is MSD plot which I am getting has too much fluctuations and i am not able to get linear fit , so i also went to search for it in other sources ; that it is wrapped trajectories,
so can anyone help me understand this part and help me in this problem?

Warm regards.

2

Hey,

Well, I do think you are gonna have a problem with the coordinates being wrapped back to the simulation domain. When saving the trajectory you should either save the coordinates in unwrapped style or save the image flags - neither of these two things are done by using the “atom” keyword.

Another thing to consider is the group of atoms you are calculating the msd of. I am not sure what atoms are included in the group “li”, “emim” and “tfsi”, but my guess is that you should probably be computing the msd based on the center of mass of polymer chains instead of individual atoms (see compute msd/chunk command — LAMMPS documentation). This is because part of the movement of the atoms in a chain may not reflect the transnational movement of the chain as a whole: imagine for example changes in atomic positions that underlie simply a change in conformation.

A last thing to consider is… are you really sure this is 50 ns is sufficient to observe significant movement in polymers? My guess is that your polymer is solid and not molten?

Hope it helps !

2 Likes
3

no I did it for 100 ns for production run so this is my full job script :-

Polymer–electrolyte system @ 298 K

Minimization → NVT → NPT → 100 ns prod

==========================================

units real
atom_style full
dimension 3
boundary p p p

---------- Force field styles ----------

bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style cvff

special_bonds lj/coul 0.0 0.0 0.5 # OPLS 1–4 scaling
newton off # matches your previous choice

---------- Nonbonded and electrostatics ----------

pair_style lj/cut/coul/long 12.0 12.0
kspace_style pppm 1.0e-5

neighbor 4.0 bin
neigh_modify every 1 delay 0 check yes one 8000

---------- Read system ----------

read_data mylammps.data
include parameters.inc

---------- Thermo output ----------

group polymer type 1 8 10
group emim type 3 4 5 6 7 11 12 13 14 15 17
group tfsi type 2 9 18 19 20
group li type 16

compute T_poly polymer temp
compute T_emim emim temp
compute T_tfsi tfsi temp
compute T_li li temp

thermo 1000
thermo_style custom step temp c_T_poly c_T_emim c_T_tfsi c_T_li press etotal pe ke vol density lx ly lz

Optional: quick sanity check of energy

run 0

==========================================

1. ENERGY MINIMIZATION (SD → CG)

==========================================

reset_timestep 0
timestep 0.5 # smaller dt during minimization is safer

(a) SD minimization with box relaxation

fix relax all box/relax iso 1.0 vmax 0.001
min_style sd
minimize 1.0e-4 1.0e-6 20000 100000

(b) CG minimization (refine)

min_style cg
minimize 1.0e-6 1.0e-8 50000 200000

unfix relax
write_data min.data

==========================================

2. NVT EQUILIBRATION (200 ps @ 1 fs)

==========================================

reset_timestep 0
timestep 1.0 # 1 fs

Initialize velocities AFTER minimization

velocity all create 298.0 12345 mom yes rot yes dist gaussian

fix mynvt all nvt temp 298.0 298.0 100.0 # Tdamp ~ 100 fs
thermo 1000

200 ps = 200,000 steps at 1 fs

run 200000

unfix mynvt
write_data nvt_200ps.data

==========================================

3. NPT EQUILIBRATION (500 ps @ 1 fs)

==========================================

reset_timestep 0
timestep 1.0 # keep 1 fs

use last configuration from NVT stage

fix mynpt all npt temp 298.0 298.0 100.0 iso 1.0 1.0 1000.0

Pdamp ~ 1000 fs (1 ps): reasonably gentle barostat

thermo 1000

500 ps = 500,000 steps at 1 fs

run 500000

unfix mynpt
write_data npt_500ps.data

==========================================

4. PRODUCTION RUN (100 ns @ 2 fs)

NVT at volume from NPT

==========================================

reset_timestep 0
timestep 2.0 # 2 fs, as requested (watch stability!)

Use NVT for production at fixed volume from NPT

fix myprod all nvt temp 298.0 298.0 100.0

Dumps and restarts

dump myDump all atom 10000 dump.lammpstrj
dump_modify myDump sort id

Write restart files so you can resume if something happens

restart 5000000 restart1.lmp restart2.lmp

thermo 5000

50 ns = 50,000 ps = 50,000,000 fs

100 ns = 100,000 ps = 100,000,000 fs

At 2 fs/step → 25,000,000 steps

At 2 fs/step → 50,000,000 steps

run 50000000 # 100 ns @ 2 fs

unfix myprod

------------- End of script -------------

Then this is my MSD script :-
-------------------- Initialization --------------------

units real
atom_style full

bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style cvff

special_bonds lj/coul 0.0 0.0 0.5 # OPLS 1–4 scaling
newton off # matches your previous choice

pair_style lj/cut/coul/long 12.0 12.0
kspace_style pppm 1.0e-5

neighbor 4.0 bin
neigh_modify every 1 delay 0 check yes one 8000

read_data min.data
include parameters.inc

-------------------- Groups --------------------

group polymer type 1 8 10
group emim type 3 4 5 6 7 11 12 13 14 15 17
group tfsi type 2 9 18 19 20
group li type 16

-------------------- Computes --------------------

Bin Z direction into thicker bins (6 Å)

compute 1 li msd com yes average yes
compute 2 emim msd com yes average yes
compute 3 tfsi msd com yes average yes

-------------------- Fixes --------------------

fix msdout all ave/time 1 1 10000 c_1[1] file msd_li.dat
fix msdout1 all ave/time 1 1 10000 c_2[1] file msd_emim.dat
fix msdout2 all ave/time 1 1 10000 c_3[1] file msd_tfsi.dat

-------------------- Rerun --------------------

Process the entire trajectory

rerun dump.lammpstrj first dump x y z box yes

4

@Shivam_Tripathy it is pointless to flood this discussion with your input files. That does not provide much additional information will not at all address the fundamental flaw in your approach.

If you spend some more time reading in your favorite text book about what is computed and how, when you try to determine diffusivity via MSD, and if you also pay attention to the details of the LAMMPS documentation about it, you will understand that without having unwrapped coordinates, i.e. include the information of atoms moving to the next periodic image (either via storing image flags or via writing unwrapped coordinates), you cannot compute the MSD correctly.

If you wrote your trajectory sufficiently frequently, you may be able to reconstruct the missing information by comparing the positions of the same atom in adjacent frames and use the heuristic that if an atom moves by more than half a box dimension that you have to add to or subtract from the atom coordinate a box length in that direction. Some MD analysis tools have this kind of feature included. Otherwise you have to repeat your simulation and then it is probably best to compute the MSD directly.