energy drift when using fix rigid in 2d simulations

Dear LAMMPS users,

(LAMMPS version 7 Aug 2019)

I am running some test NVE simulations of 2d rigid molecules, and I have noticed a systematic drift of the kinetic energy over time. A sample input file is reported below.

In the simulations, each molecule is composed by 10 atoms arranged in a pentagonal shape, with 2 atoms per edge (see attached data file “data”). The 2 atoms in one “edge” of the pentagon overlap strongly; however, this should not play a role in the dynamics since I turn off intra-molecular forces by using the command “neigh_modify exclude molecule/intra all”.

For long enough simulations, the kinetic energy drifts markedly (see attached plot ke_all.pdf).

Some things I’ve tried:

-Using a smaller time step, down to 1e-4. The drift is still present, and the slope (dE/dt) is basically unchanged. This is rather surprising to me since I would expect the slope of the drift to become smaller if the time step is smaller.

-Using different fixes. I’ve tried fix (a) rigid, (b) rigid/nve, © rigid/small and (d) rigid/nve/small. For (a-c), I observe the same qualitative behavior, as shown in the attached plot. For (d), I systematically get a segmentation fault error (which could be related to the energy drift problem).

-Using different interaction potentials. In particular, I’ve tried pair_style soft, which seems to worsen the problem.

I have checked that the initial configuration does not contain inter-molecular overlaps.

I have also tried to run the same simulations with similar pentagonal 2d molecules, but with each molecule composed by only 5 atoms arranged on the vertices of a regular pentagon (see attached file “data_5”). When simulating this kind of molecules, the energy drift seems to be absent over comparable time scales.
It is therefore possible that the drift is related to the strong overlaps present in the molecules, although, as mentioned above, I make sure to turn off intra-molecule interactions.

Thanks for your patience. Any suggestion is greatly appreciated.

Best regards,

data (16.6 KB)

ke_all.pdf (25.7 KB)

data_5 (8.34 KB)

PS: I should have specified this before, however with this choice of interactions the potential energy flucyuates close to zero and doesn’t show an appreciable drift, but rather increasingly strong fluctuations. The total energy is approximately equal to the kinetic energy (and therefore shows the same drift).

I am CCing Trung Nguyen who may have ideas about this.


Dear LAMMPS users and dear Steve,

Thanks for considering my question.

I was finally able, with some external help, to solve the problem!

The reason for this behavior was the data file: indeed, by looking at the z coordinates of the particles in “data”, one can see that some atoms have z coordinates very close to zlo=0, and other atoms have z coordinate very close to zhi=0.1 (I am using a finite but narrow zlo-zhi range, as suggested on the manual The latter particles also have a z flag value of -1.

It was surprising to me to discover this problem with the data file, and I’ll explain why below (although probably I should do so in a separate post…)

When I start the simulation, I assign to each atom a z coordinate z=0.05=(zhi-zlo)/2 (again, as suggested on the manual: “For each atom in the file, assign a z coordinate so it falls inside the z-boundaries of the box”).
The initial coordinates of the molecules (pentagons) are random, so that initially overlap are present.
In order to remove these overlaps, I use the script below, in which two consecutive runs are performed with a soft potential (pair_style soft).
The strength of the potential (“A” constant) is ramped from 0 to 20 during the first (shorter) run, then from 20 to 100 during the second (longer) run.

The data file written at the end of this simulation (through the write_data command) is then used as input for the script reported in the original post.

When running this script on a data file in which all the atoms have z=0.05 (see attached file “data_initial”), the result is the change of value of the z coordinates described above.

However, it is sufficient to set z=0 for all particles for this problem to disappear.

If this is done, the energy drift problem described in my initial email also disappears!

It is not clear to me what causes the script below to give different results depending on whether the atoms initially have z=0.05 or z=0.

I attach two data files, one in which all the atoms have z=0 and one in which they have z=0.05. In both cases, zlo=0 and zhi=0.1.

Below is the script used to remove the overlaps.

You can verify that the data file written at the end of the simulation is different depending on which file is used as input.

Thanks again for your help.

Best regards,


PS: Another interesting (?) observation: when simulating the pentagons composed by 5 atoms (in which no intra-molecular overlaps are present – file data_5_initial attached below), the z coordinate of all the atoms is set to 0 after the script to remove the overlap (below) is run (!). This behavior is different from the one found for the pentagons composed by 10 atoms, which I described above.

I meant to dig a bit more into this problem, to see whether it was possible to reproduce the same behaviors with simpler shapes instead of pentagons, but I still haven’t had the time…


variable tpushoff_1 equal 1e4
variable tpushoff_2 equal 1e6
variable mytemp equal 1.0
variable thermodump equal ((v_tpushoff_1+v_tpushoff_2)/1e3) variable trajdump equal ((v_tpushoff_1+v_tpushoff_2)/20)
variable tstep equal 0.0001
variable damp equal 1.0
variable soft_cutoff equal 1.1225
variable soft_energy_max_1 equal 20.0
variable soft_energy_max_2 equal 100.0
variable lseed equal 34314142
variable fname string data
variable thermofile string thermo.dat
variable simname string push

units lj
boundary p p p
atom_style bond
dimension 2

read_data ${fname}

pair_style soft {soft_cutoff} pair_coeff * * 0.0 {soft_cutoff}
variable soft_energy equal ramp(0,${soft_energy_max_1})

neighbor 0.5 bin
neigh_modify every 5 delay 0 check yes
neigh_modify exclude molecule/intra all

variable step equal step
variable temp equal temp
variable etot equal etotal
variable ke equal ke
variable pe equal pe
variable press equal press

thermo_style custom step etotal ke pe temp press
thermo ${thermodump}

timestep ${tstep}

fix fix_print all print {thermodump} "{step} {etot} {ke} {pe} {temp} {press}" file {thermofile} screen no title “#1:step 2:etot 3:ke 4:pe 5:temp 6:press”

fix fix_rigid_damp all rigid molecule langevin {mytemp} {mytemp} {damp} {lseed}
fix fix_zero_linmom all momentum 10 linear 1 1 1
fix fix_zero_angmom all momentum 10 angular
fix fix_push_1 all adapt 1 pair soft a * * v_soft_energy
fix fix_2d_1 all enforce2d

run ${tpushoff_1}

unfix fix_push_1
unfix fix_2d_1

variable soft_energy equal ramp({soft_energy_max_1},{soft_energy_max_2})

fix fix_push_2 all adapt 1 pair soft a * * v_soft_energy
fix fix_2d_2 all enforce2d

run ${tpushoff_2}

write_data data.{simname} write_restart restart/restart.{simname}

data_initial_z0 (12.4 KB)

data_initial (12.4 KB)

data_5_initial (6.21 KB)

I think we should add a test for valid imageflags to the data file processing.
It should refuse to run with imz != 0 and also use utils::inumeric() to converting to detect non-integer data.
steve, i am proposing such a change in pull request #1945.


Hi Valerio,

when I set the image flags in z for the particles to be zero the energy drift seems to go away, so it seems to be related to the issue and solution Axel mentioned.

You can try the attached data file.


data2 (16.5 KB)

Dear Trung,

Thanks for your reply.

I see that in your version of the data file, in addition to set z_flag=0 for all particles, you set the z coordinates which are close to z_hi=0.1 to a very small value (order 10^-12), and change the z range from [0, 0.1] to [-5, 0.5]. This indeed seem to solve the problem.

As explained in my previous email, also setting z=0 for all particles while leaving the z range unchanged solves the problem (independently of the values of the z flags…!).

What I still don’t understand is why the z coordinates end up having these values, since in the data file I start from -before removing the overlaps- all the particles have z=(zhi-zlo)/2 (as explained in my previous email), and why the z coordinates are unchanged when the 5-atoms pentagons are simulated (also discussed in the previous email).

This, however, should probably be discussed in a separate thread.

First, I want to check the script which removes the overlaps more carefully and to get repeatable results.