Error in Running LAMMPS in Loops

Dear LAMMPS Community,

I am simulating a single particle using the Langevin equation and have encountered an issue I hope you can help with. My workflow involves running LAMMPS (version 29 Aug 2024) in a loop, where I,
1. Run a LAMMPS simulation for a certain number of timesteps.
2. Save the restart file.
3. Use this restart file to continue the simulation for the next loop iteration.
4. Repeat this process for 100,000 loops.

When I plot the mean squared displacement (MSD) as a function of time (t), I observe superdiffusive behavior at longer time, which is unexpected for this setup—I should be seeing normal diffusive behavior.

Interestingly, when I run LAMMPS continuously (without using loops and restart files), I observe the expected diffusive behavior.

I have checked my simulation parameters carefully and ensured consistency between the looped and continuous setups, but the discrepancy persists. Could the use of restart files in a looped setup introduce subtle errors, or is there something else I might be overlooking?

Thank you in advance for your help!

Here is the input script (LAMMPS version 29 Aug 2024):
############################################################
############################################################

Box and units

units lj
atom_style atomic
boundary p p p

READ “start” data file

read_restart restart.restart

log log.txt append

pair_style lj/cut 1.12246152962189
pair_coeff 1 1 1.0 1.0 1.12246152962189

variable seed equal 70

fix 1 all nve # NVE integrator
fix 2 all langevin 1.0 1.0 1.0 ${seed} # langevin thermostat

Output thermodynamic info (temperature, energy, pressure, etc.)##########################

thermo 10000
thermo_style custom step temp etotal pe ke epair emol press vol
#############################################################################################

set timestep of integrator

timestep 0.1

dump 1 all custom 1 dump.position.txt id type xu yu zu

run 1
write_restart restart.restart

#######################

It is standard for particle mean square displacements to display ballistic behaviour over very short time scales.

I cannot determine on your behalf if that explains your observations – but that is the first thing I would check.

Thank you for your response.

You are correct that MSD typically displays ballistic behavior at very short timescales. I observe this in both cases. However, the issue arises at longer timescales:

  1. Single Continuous Run: The MSD shows ballistic behavior initially, followed by diffusive behavior at longer timescales, as expected.
  2. Looped Runs with Restart Files: The MSD transitions from ballistic to diffusive for some time but then becomes superdiffusive at longer timescales (as shown in fig below).

Are you trying to use that approach in which you collect several different (time, MSD) data sets from a single simulation by “restarting” the t=0 reference?

I dont know why your method of saving the restart files and doing new simulations isnt working (to be seen by someone else in the continuation of your thread), but I can propose a much more computational friendly way for that, which you can also use to compare the results you are getting with. The approach is:

Do a simulation like the one you are doing there but, instead of saving restart files every once in a while, save a trajectory. Afterwards, you can create a new simulation that contemplates a loop where this trajectory is read N times. For each of these N times, a different portion of the trajectory is read and a single (t, MSD) data is computed and output. The idea is that you are going to define some variables that, at each iteration of the loop, assumes the value of tstart and tend corresponding to the portion of the trajectory you intend to read and consider in the MSD calculation.

Ultimately, you will have N sets of (t, MSD) data out of the NVT simulation you initially ran and got the trajectory from. From there, suffices to create a python code to average them all and get the MSD curve you seem to want :slight_smile:

I will leave some commands that are relevant for the approach I am proposing. Tomorrow I can try to help in a more straightforward way if you need.

rerun command — LAMMPS documentation,
jump command — LAMMPS documentation,
variable command — LAMMPS documentation

2 Likes

Thank you for the suggestion! :blush:

I tried your method of using the trajectory and looping through it, and now I’m getting the expected results.
I’m still looking into why my original method with restart files isn’t working as expected.

Thanks again for your help!

1 Like

When you said “restart” again it finally made sense: fix langevin does not store its random generator state in a restart file (which is documented behaviour).

Thus, by starting a series of 1-step “runs” from restarts, the “random” Langevin force experienced by each particle is identical over all timesteps (still different for each particle) as long as no particles move between processors. This is a recipe for overcorrelated motion, which shows up as superdiffusion in the short and long timescales (while still diffusive on the timescale of particle collisions).

Thank you for pointing that out!

However, in my loop, when I restart the simulation using the restart file, I also open the LAMMPS input script and change the seed value for each subsequent run. Because of this, I don’t think the identical random Langevin force is causing the superdiffusive motion.

Well, you could have mentioned right at the start that you’re choosing a new Langevin seed each run. Without any details of how you are doing so it is impossible to determine whether your random Langevin forces are in fact correlated.