Hello lammps users,
I am running into an issue of the barostat
not achieving the
target pressure (zero) in a system with rigid and non-rigid bodies. This
issue was discussed previously on the forum, with the subject titled
"global pressure with rigid and non-rigid bodies." I was unsuccessful in
replying to the old thread, so I have started a new one. From the old
thread, I was not able to figure out why the barostat (nose-hoover) was
not achieving the target pressure, and so I eventually moved forward
with a Berendsen barostat because it was reaching the target pressure as
desired. However, my latest simulations require calculation of a bulk
modulus, and Berendsen does not capture the correct physics for a bulk
modulus determination. I recently tried going back to the Nose-Hoover
barostat, and not surprisingly, the pressure issue has returned.
I
believe there is a bug in the way the
pressure is handled for systems with rigid and non-rigid bodies when
using a Nose-Hoover barostat, and I have provided supporting evidence
below.
In my system, I use fix
rigid/nvt on the rigid bodies and npt on the non-rigid. I have attached a
sample input file and data file for a system with 16 rigid bodies in a
polymer melt; I apply a small amount of positive pressure
initially (ramp the LJ pressure from 0.01 to 0 to prevent the system
from
vaporizing) and then switch the target pressure to zero and monitor the
pressure. To reduce noise in the output, I have also printed a
pressure that is averaged over 1000 timesteps using the fix ave/time
command. In my experience, this output reveals a positive pressure that
scales linearly in the number of rigid bodies included in the simulation
(although the target pressure is zero in each case), as shown in the
attachment.
I am running the March 2017 version of lammps. The
following list includes troubleshooting techniques that I have employed
in trying to understand what is going on with the system:
- Using "neigh_modify exclude molecule/intra" to prevent bogus pressure contributions due to intraparticle forces
- reducing the timestep from 0.005 tau LJ to 0.001 tau LJ (no change in the results)
- changing the pressure damping parameter from 2 tau to 20 tau (no change in the results)
- Use
"compute pressure" to get the pressure that the barostat uses/computes
to make box volume adjustments (this number is consistently positive and
fluctuates around a value greater than zero; I am not sure why the
barostat does not increase the box volume to reduce the pressure to the
target value of zero)
- invoking the barostat on the rigid body
time integration: fix rigid/npt for rigid body, fix nvt for non-rigid
(this throws many FENE bond warnings, and the pressure still does not go
to the target value of zero)
- using "aniso" instead of "iso" for altering the box dimensions due to pressure fluctuations (no change)
- I am allowing the rigid body COM's to dilate with the changing box size ("dilate all"; this command is a default anyways)
- Reducing the volume fraction of rigid bodies to rule out possible jamming effects
I
want to emphasize that this is not an issue of fluctuations in the
pressure; I am seeing fluctuations around a mean pressure value that is
not the target value. I would greatly appreciate if someone could
confirm this result by running the attached sample system and observing
the pressure magnitude in the thermo output. Thank you for any
feedback.
Scott