Kremer-Grest polymer simulation

Dear all,

I’m trying to do the same simulation with this article.
K. Kremer, G. S. Grest, J. Chem. Phys., 92(8), 5057 (1990)

I did it on LAMMPS, and to check the results, I analyzed the static properties, <R^2> and <RG^2>, and the monomer motion, mean-square displacement (MSD).
The <R^2> and <RG^2> result was good agreement with the article.
However, MSD was not.

I attached the figure of MSD of my simulation.
This is the result of ‘the number of monomers N= 200, the number of chains M= 20’.
Initial structure of this simulation was made from ‘’, and relaxed 5 million steps.
Although this is the average of 20 polymers, it is not smooth line as I looked in the article.
The dots are not convergent at all range of time.

I want to know the reason of this and want to know how to get the correct figure.
I also attached my input file.
Please teach me which line is wrong.

Yours sincerely,

Munehiro Hayashi

monomer_motion.doc (23.5 KB) (404 Bytes)

PARM.FILE (241 Bytes)

Could be a number of ways to go wrong. In general,
the question "I did a long simulation and the answer
isn't what I expect. Why?" is not easy to answer.

Did you correct the MSD for drift of the entire system?
That can be a source of systematic error, if the slope
is wrong. You can also get a less noisy result if
you model a bigger system (better stats).


Another issue in diffusion calculations is that a lot of groups will start the calculation from multiple data points within the trajectory to increase the number of data included in the averaging. That will often smooth out the data at short times, but obviously doesn’t help at long times, where there just isn’t a lot of data to use.


Dear Steve and Ahmed

Thank you for replying so quickly.

The points you indicated are exactly right what I missed.

I’ll try to change the way of averaging data.

May I confirm about the drift of the system?

Do I have to correct the moving distance for MSD, when I use ‘fix momentum command’ in this simulation?

I didn’t do this correction.

Because if the system happens to flow or rotate without that command, I think the data of MSD becomes too high.

I’m sorry. I would like to know more.

Thank you

Muehiro Hayashi

Also, be very careful of your box size in diffusion… The box length should be big enough to keep chains from seeing themselves and there will still be some finite size effect. There are analytical estimates for effect of PBC on short chains, but for entangled chains you can check D vs 1/L.


Dear Sir

The box size of this simulation was set to the density= 0.85.

In the case of short chain, the model across some image cell.

Does it have any trouble in the PBC condition?

This article was come out 20 years ago. I think now we can treat bigger sized simulation.

I’ll check the motion of the chains.

Thank you

Munehiro Hayashi

First, I would simply monitor COM drift with the compute com
command. If there is drift, you can use the compute msd com yes
command to subtract it out. Or use the fix momentum command
you mentioned. Either way you should monitor it.