Temperature of ellipsoidal fluid growing in NVE ensemble despite significant NVT equilibration

Dear LAMMPS team,

I would appreciate your insight into why I’ve been running into this runaway temperature problem while using LAMMPS’s aspherical methods. In this simulation, I use the settings in the example ‘ellipse’ in LAMMPS package (in.ellipse.gayberne, version 15Apr2020); my full script is at the end of this message. I set only one type of ellipsoid in the system and keep the mass and volume of the ellipsoid to be the same as in the original example. I varied the aspect ratio of the ellipses between 1 and 3, which interestingly led to very different behaviors in the NVE ensemble.

I ran the simulation in three stages: 100000 steps of nve/limit to eliminate/mitigate overlaps in case any were present, 100000 steps of nvt/asphere to equilibrate the temperature, and then 1000000 steps of nve/asphere (during which I hoped the temperature would be stable).

Unfortunately this is not what I found. I plot below the temperature in NVT and NVE. When the atoms are perfectly spherical (aspect ratio of 1), the temperature remains stable during NVE, with fluctuations that are roughly on the order I would expect from stat mech. However, when I changed the aspect ratio (x) of the ellipsoid to be 3:1, the temperature continuously grows in the NVE ensemble (left side of the green line is NVT ensemble and right side is NVE).

图片.png

The growing trend can also be observed when we change the original non-dimensional temperature to be 0.5 or 5. In addition, I’ve tried changing the length of the NVT run as well as changing the damping parameter (0.1 through 10), but the trend of growing temperature still persists. As an additional note, I’ve tried some extremely long NVE runs, and the temperature increase persists for the entire duration after the switch from NVT to NVE.

图片.png图片.png图片.png

As a sanity check, I’ve plotted the energy change during the NVE run, and confirmed that indeed the kinetic energy (and thus the total energy) keeps increasing in the NVE ensembles. Thus I suspect that I have either misused the aspherical integrators or, perhaps, that there is a bug in these integrators.

图片.png

In summary, the problem we find is that when we use the NVE ensemble to collect data for aspherical atoms, we observe a stark trend of temperature increase, which cannot be solved by changing equilibration time, damping parameter, or original temperature. I would greatly appreciate any ideas that anybody has on this matter.

Sincerely,

Shuyuan Wang

The full script of the simulation:

# GayBerne ellipsoids in LJ background fluid
units       lj
atom_style   ellipsoid
dimension    2

lattice     sq 0.02
region      box block 0 20 0 20 -0.5 0.5
create_box   1 box
create_atoms 1 box

set         group all type/fraction 1 0.1 95392
set         type 1 mass 1.0
set         type 1 shape 2.08 0.693 0.693
set         group all quat/random 18238

compute      rot all temp/asphere dof all
velocity     all create 1 187287 loop geom

pair_style   gayberne 1.0 3.0 1.0 4.0
pair_coeff   1 1 1.0 1.0 1 1 1 1 1 1 2.5

neighbor     0.8 bin
thermo_style custom step c_rot epair ke pe etotal press vol
thermo      100
timestep     0.002

compute      orient all property/atom quatw quati quatj quatk
compute shape all property/atom shapex shapey shapez

fix 1 all nve/limit 0.01
run 100000

reset_timestep 0
unfix 1
fix temp all ave/time 1 1 100 c_rot file temp.profile
fix 1 all nvt/asphere temp 1 1 1
run 100000

unfix       1
fix         1 all nve/asphere
dump        1 all custom 100 dump.ellipse.gayberne.lammpstrj &
        id type x y z c_orient[1] c_orient[2] c_orient[3] c_orient[4] c_shape[*]
run         1000000

i am missing any mention of reducing the time step. please note that examples in the LAMMPS distribution are meant to serve as examples on how to use certain commands (i.e. demonstrate the syntax and combinations of commands that are useful) usually in short runs, NOT as a guideline for how to do good simulations.

i am not very surprised that you see better energy conservation for spherical particles and see a degradation with a significant degree of aspherical nature. i have very little practical experience with these kind of systems, but am familiar with a similar behavior for rigid objects, where you - similarly - have separate time integration of the center of mass translation and of the rotation around the center of mass. the integration of the rotational degrees of freedom is much more fragile and sensitive to the choice of timestep.

axel.

图片.png

图片.png

图片.png

图片.png

图片.png