Wrong temperature value when simulating rigid linear molecules in conjunction with pair_style zero

Dear LAMMPS users,

(LAMMPS version 3 Mar 2020)

I am writing in the hope that someone can help me understand something strange I observed when running NVT simulations of rigid molecules.

The molecules I’m simulating are composed of three particles: a central particle (type 1) and two interaction sites (type 2), located at a distance 0.5 from the center of the central particle.

The three particles lay on the same line, and therefore constitute a rigid body with 5 degrees of freedom. The model is basically the same as the one employed in this paper: https://doi.org/10.1088/0953-8984/20/15/155101.

I was running some tests in which the 1-1 interaction is purely repulsive (WCA, i.e. Lennard-Jones cut and shifted at its minimum, r=2^(1/6)) whereas the 1-2 and 2-2 interactions are completely turned off (pair_style zero).

I am running langevin dynamics at low temperature, T=0.05, using “fix rigid/nve/small molecule langevin”. The damp coefficient is 0.1, the timestep 0.001. Intra-molecular interactions are neglected (neigh_modify exclude molecule/intra all). The input script is reported below.

What I noticed is the following (see attached plots, plot_1 and plot_2):

-The temperature doesn’t have the desired value (0.05), but to a value ~3% lower (~0.0485)
-There seems to be a (very small) drift in the temperature/kinetic energy

-Some “spikes” can be observed in the temperature/kinetic energy at long simulation times.

What is strange is that the temperature is initialized correctly, but drops almost immediately to the lower value.

Here are some things I have tried (with no results):

-Decreasing the timestep

-Changing the damp parameter

-Using fix rigid/nve instead of fix rigid/nve/small

-Changing the 1-1 interaction from WCA to a tabulated purely repulsive potential U®=r^(-A), with A=100,200

I have also run some NVE simulations of the same system with the same timestep and initial temperature, and the total energy is conserved.
Moreover, the problem seems to be even worse if I use some attractive potential instead of the “zero” pair style for the 2-2 interaction (this is what I would like to do in the end).

Finally, here is what puzzles me the most:

-If instead of disposing the 3 particles on a line I arrange them in a rectangular triangle configuration (90 deg. instead of 180 deg.), creating therefore a rigid body with 6 degrees of freedom instead of 5, the problem is not present anymore (see attached plots plot_1 and plot_2). The problem could therefore be related to the counting of the degrees of freedom of the rigid body (?)

-If I take 3 identical type 1 particles, all interacting with the WCA potential, arranged in the same linear configuration, the problem is not present anymore. It could therefore come from the use of the “zero” pair style in conjunction with fix rigid (?)

I report below my input script, and I attach two data files, one with the particles arranged on a line (180.lammpsdata) or in a rectangular triangle (90.lammpsdata).

Thank you in advance for your help.

Best regards,

Valerio

plot_1.pdf (364 KB)

plot_2.pdf (242 KB)

180.lammpsdata (98 KB)

90.lammpsdata (98 KB)

I cannot give you a complete explanation of what is happening, but only some observations and a speculation.

  • this has nothing to do with pair style zero. you get the same behavior with
    pair_style lj/cut ${lj_cutoff}
    pair_modify shift yes
    pair_coeff 1 1 1.0 1.0
    pair_coeff 1 2 0.0 1.0 1.0
    pair_coeff 2 2 0.0 1.0 1.0

and optionally (to regain the speed lost by looping over atoms with epsilon == 0.0):

neigh_modify exclude type 1 2 exclude type 2 2

  • the thermostatting works as expected with fix rigid/nvt and fix rigid/nvt/small

  • i recall that there were some changes to fix langevin for extended particles (spherical and aspherical), where the gamma parameter had a different prefactor, so I would speculate that something similar would be required here and the langevin thermostat built into fix rigid and its derived classes is missing that feature.

axel.

Dear Axel,

Thank you for your reply.

I cannot give you a complete explanation of what is happening, but only some observations and a speculation.

  • this has nothing to do with pair style zero. you get the same behavior with
    pair_style lj/cut ${lj_cutoff}
    pair_modify shift yes
    pair_coeff 1 1 1.0 1.0
    pair_coeff 1 2 0.0 1.0 1.0
    pair_coeff 2 2 0.0 1.0 1.0

Yes, you’re right. Actually I kinda contradicted myself when saying that pair style zero could be the cause, since I said in the previous message that using an attractive potential for the 2-2 interaction gave me the same problem (I used a tabulated generalized Gaussian potential).

So the problem must be more general.

and optionally (to regain the speed lost by looping over atoms with epsilon == 0.0):

neigh_modify exclude type 1 2 exclude type 2 2

  • the thermostatting works as expected with fix rigid/nvt and fix rigid/nvt/small

  • i recall that there were some changes to fix langevin for extended particles (spherical and aspherical), where the gamma parameter had a different prefactor, so I would speculate that something similar would be required here and the langevin thermostat built into fix rigid and its derived classes is missing that feature.

axel.

I see. Hopefully someone will be able to find the cause or some up with some ideas. I’ll do some more tests in the meantime.

Thanks again.

Best,

Valerio

Dear Axel, dear LAMMPS users:

A quick follow up;
Using the same script as reported in the original message, I get “reasonable” behavior (although the temperature problem is still not solved) using as input the following data file:

The best way to ensure that somebody is looking into this (eventually) is to file a bug report with explanations and a complete and easily runnable input deck attached in the GitHub issue tracker of the LAMMPS project: https://github.com/lammps/lammps/issues

The mailing list is not as closely monitored by all LAMMPS developers as it used to be and - given the volume of messages on the list - bug reports can easily drop out of sight easily, while GitHub issues are more persistent and require an explicit dismissal.

Axel.