Problem with sllod energies and averages

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

We are having problems with measurements using fix nvt/slllod. The potential energy (and therefore the forces) after the first time step appear incorrect (jumps discontinuously). This can most easily be seen by plotting the potential energy every time step and examining the difference between steps (LJeq.lmp attached). If the time step size is shrunk, the effect continues to be only in the first step (LJ_shorter_timestep attached). Similarly if we run the system for 1 timestep, then run it again for 1 timestep, the potential energy jumps going from 0-1, but then is correct at time 1 for the beginning of the next time step, before jumping again from 1-2 (LJ_run_loop.lmp attached). We have also seen a similar effect in the remapping of the box, which may be associated with the same problem.

We have performed MD simulations using our own code (4th order Runge-Kutta), and it appears to be associated with the boundary conditions; the force profile appears consistent with the boundaries not moving on the first time step, suggesting that the boundaries are continually lagging the system state.

The effect is relatively small in the energy, so I have attached scripts for 128 particle LJ system that will generate the error.

LJeq.lmp

  1. Equilibrates the system.

  2. Does 10 time steps with a strain rate of 0.

  3. Imposes a velocity profile consistent with a strain rate of 0.5 and saves the system state (restart.start.lmp).

  4. Does 10 time steps with a strain rate of 0.5 and dt=0.001.

Note that Delta_PE between timesteps both with a strain rate of zero and of 0.5 is of approximately 0.0013, except for the first time step that is 0.0023.

LJ_shorter_timestep.lmp

  1. Loads system initial system state from LJeq.lmp

  2. Does 10 timesteps with a strain rate of 0.5 and dt=0.0001.

Note that delta PE is now of the order of 0.0001 for all time steps except the first.

LJ_run_loop.lmp

  1. Loads system initial system state from LJeq.lmp

  2. Does a loop of 1 timestep 10 times with a strain rate of 0.5 and dt=0.001.

Note that the PE between the end of the last step and the beginning of the next step changes discontinuously, and the start value looks more consistent with the trend from LJeq.lmp (eg Delta PE for the first time step taken between the start of the second time step and the start of the first is approximately 0.001). Performing simulations with this loop structure matches significantly better with our code.

Any suggestions would be appreciated.

James Reid

LJ_run_loop.lmp (680 Bytes)

LJ_shorter_timestep.lmp (702 Bytes)

LJeq.lmp (1.43 KB)

2

Hi - I’m not clear what issue you are flagging here. In general, there
are two related issues with using fix nvt/sllod and fix deform on
the early timesteps, if you are starting a system from a standstill.

  1. If the initial deformation is just turned on instantaneously,
    it can induce a shock at the left edge. A bigger problem for
    solids, but still can be a (hopefully small) effect in liquids.

  2. Fix nvt/sllod subtracts off a velocity profile that it assumes
    matches the box deformation rate. When the system is
    nicely flowing with the deforming box, that is the case. Initially,
    however the fluid is stationary, so the temperature computed by

nvt/sllod is bogus and thus its thermostatting effect is too.

Both these effects should be transient and disappear as
the deformed system comes to an equilibrium flow. I’m
not clear if what you are describing is related to these

effects and thus you should not be worrying about what
happens in the first few timesteps ??

Steve

3

Hi Steve,

Thanks for the reply, but it’s not due to either of those problems, as it can be observed even after a long run in the steady state. In reverse order:

  1. We explicitly impose the correct velocity profile on the system so that it matches the box deformation rate – this can be seen by looking at the temperature, kinetic energy, pressure, or any quantity that removes the streaming rate which match the values for the system before we impose the velocity profile and the sllod flow. Furthermore, the same effect can be seen by pausing and restarting the sllod flow.

  2. Since we have the correct velocity profile, there shouldn’t be a shock at the left edge as the sllod transformation is volume preserving.

To boil the question down to its most basic, if I do a long sllod run so that the system is at the right system point, and I am outputting any variable based on force/potential, why, if I do

Run 1

Run 1

do I get a different result from the same configuration then if I go

Run 2

This is clearly not just a problem with the thermo output, because the trajectories will diverge with time if I keep incrementing them. We are working on a project and we are not getting the expected results from sllod that we are from other non-equilibrium perturbations in lammps. We’re not sure if this is the problem that is causing the disagreement, but something seems to be wrong here.

James Reid