Hi everybody,
We have been trying to run Langevin dynamics simulations of WCA spheres with shear deformation. Here is a minimal script that we use for that purpose:
units lj
atom_style atomic
lattice sc 0.97
region myBox prism 0 10 0 10 0 10 0 0 0
create_box 1 myBox
create_atoms 1 region myBox
pair_style lj/cut 1.1225
pair_modify shift yes
pair_coeff 1 1 1.0 1.0 1.1225
mass 1 1.0
timestep 0.002
fix dfmFix all deform 1 xy erate 0.2 remap v
compute dfmTemp all temp/deform
fix langFix all langevin 0.5 0.5 1.0 141 zero yes
fix_modify langFix temp dfmTemp
fix ensFix all nve
thermo 50
thermo_style custom step c_dfmTemp pe etotal
log log-lammps-shear0.2_timestep0.002-tD1.0.txt
run 500000
The temperature histogram that we get from the dfmTmp compute (made italic above), however, deviates significantly from the target temperature of 0.5. The deviation is absent in conventional Langevin dynamics and is only observed when fix deform is utilized. It also becomes worse when the shear rate increases. I tried this with multiple LAMMPS versions (including 12Dec18 and 29Sep21) and with time steps as small as 0.0002 (10 times smaller than what I have above) and always observe this behavior. Indeed the temperature histogram is independent from time step and only depends on the shear rate. I would also like to note that this effect is not a consequence of what one would typically consider as “heating up”, since the temperature does not tend to drift upwards along a trajectory but rather fluctuates around the wrong value for extended periods of time. For this particular script, for instance, the average temperature [the output of c_dfmTemp] is around 0.58.
I looked at the LAMMPS source code, and as I understand this should not happen as the random and drag forces are applied to unbiased velocities. Am I missing something in setting up this script? Or might there be a bug in the thermostat that I am not easily able to note?
Regards,
Amir