[lammps-users] Thermostating issues with Langevin dynamics alongside shear deformation

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 (highlighted in light blue 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 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

Using fix temp/deform is only appropriate if you have a liquid that is flowing consistent with the
shear imposed by fix deform. Is that the case with your model? You can test that by plotting
the velocity profile of your liquid in the flow direction and seeing if is at equilibrium with the
deformation rate. If is not (yet) then the temperature calculated by compute temp/deform
will be off, b/c it subtracts that ideal velocity profile from the liquid before computing the temp.

Steve

Thanks Steve for your response. We actually re-checked the velocity profiles and they match perfectly with the fully developed profiles expected from the applied shear rate. Moreover, the obtained profiles are independent of the sampling time, while the apparent temperature increases with sampling time. So this does not seem to arise from insufficient sampling.

Then I’m not sure. Here are a couple more ideas.

(1) Shearing adds energy to the system. Shearing fast adds it fast.
If your Langevin time constant is too long it may not be able to remove the energy
fast enough and thus T will rise

(2) Two alternatives to using Langevin + compute temp/deform. The first is
to use Langevin with compute temp/partial and only thermostat the 2 non-shearing dims.
Y and Z in your case. Make sure you are also monitoring the temp only in those 2 dims.
If your x-velocity profiles are good, that takes care of the x-direction temperature, so you
should have a well-equilibrated system if the YZ temp is stable.

(3) The other alternative is to use fix nve/sllod as the thermostat. It incudes the flow to match
the box deformation and it does things “more correctly” in the macroscopic sense than Langevin by adding a
special term to the global thermostat.

Steve