global pressure with rigid and non-rigid bodies

LAMMPS LAMMPS Mailing List Mirror
1

Hello lammps users,

I have a system with rigid and non-rigid bodies that I am attempting to barostat and thermostat. I have implemented npt for the non-rigid bodies and rigid/nvt for the rigid bodies as suggested in the fix rigid documentation:

fix 2 filler rigid/nvt molecule temp 1.0 1.0 2.0
fix 3 polymer npt temp 1.0 1.0 2 iso 0 0 2 dilate all

The thermostat seems to be working just fine, but the global pressure is not going to zero. Rather, it is consistently positive. Even though all atoms in the system are re-scaled with the changing box size, I am wondering if the large volume fraction of rigid bodies (0.33) is somehow causing issues with the barostat.

I am using lammps-31Oct15. Thank you for any feedback.

2

Hello lammps users,

I have a system with rigid and non-rigid bodies that I am attempting to
barostat and thermostat. I have implemented npt for the non-rigid bodies and
rigid/nvt for the rigid bodies as suggested in the fix rigid documentation:

fix 2 filler rigid/nvt molecule temp 1.0 1.0 2.0
fix 3 polymer npt temp 1.0 1.0 2 iso 0 0 2 dilate all

The thermostat seems to be working just fine, but the global pressure is not
going to zero. Rather, it is consistently positive. Even though all atoms in

*how* positive?

the system are re-scaled with the changing box size, I am wondering if the
large volume fraction of rigid bodies (0.33) is somehow causing issues with
the barostat.

check out using "neigh_modify exclude molecule" to avoid computing
force internal to the rigid bodies, which can cause bogus
contributions to the total stress and thus pressure.

in addition, i would use the npt integrator on the rigid bodies, as
that guarantees that the cell changes will displace the rigid bodies
and not individual atoms, but that

axel.

3

I am using LJ units and the pressure is around 0.06. I tried your suggestion of using rigid/npt on the rigid bodies, but the results are the same as from the rigid/nvt integration (positive pressure). I have been using “neigh_modify exclude molecule” in all runs.

4

I am using LJ units and the pressure is around 0.06. I tried your suggestion
of using rigid/npt on the rigid bodies, but the results are the same as from
the rigid/nvt integration (positive pressure). I have been using
"neigh_modify exclude molecule" in all runs.

without more specific information about your simulations (i.e. a
complete input deck), it is difficult to make any specific statements.
the only detail that still stands out, is that you use isotropic cell
changes. your system might prefer a different box shape due to packing
effects. but as mentioned before, that is pure speculation.

axel.

5

Please find my attached input script and data file. My system contains 5000 polymer chains that are 20 beads in length, 2500 single beads that function as crosslinkers in later phases of my simulation, and 198 rigid bodies, which contain 255 beads each. The input script can be broken down into five phases: 1) Apply a box deformation (compression) in all 3 dimensions to only the non-rigid elements to get the simulation box close to the zero pressure equilibrium size; do not time integrate the rigid bodies, as their data file coordinates are already positioned within the dimensions of a box at equilibrium zero pressure 2) Apply a soft potential to the rigid bodies in case of overlap; use LJ/cut for the non-rigid and cross interactions 3) Switch back to LJ/cut for all nonbonded interactions; time integrate the non-rigid elements with NPT and apply positive pressure to ensure the box does not blow up; time integrate the filler with rigid/nvt 4) Ramp the pressure to zero 5) run long enough to see the pressure stabilize and take note of the value (should be zero).

When I reach the final phase of this run, the pressure is not going to zero–it stabilizes around 0.06-0.08. I have also tried the suggestion of using “aniso” on the box pressure, which results in the same outcome.

bulk_PS00007.data (7.79 MB)

combined_input_PS00007.inp (3.07 KB)

6

Please find my attached input script and data file. My system contains 5000
polymer chains that are 20 beads in length, 2500 single beads that function
as crosslinkers in later phases of my simulation, and 198 rigid bodies,
which contain 255 beads each. The input script can be broken down into five
phases: 1) Apply a box deformation (compression) in all 3 dimensions to only
the non-rigid elements to get the simulation box close to the zero pressure
equilibrium size; do not time integrate the rigid bodies, as their data file
coordinates are already positioned within the dimensions of a box at
equilibrium zero pressure 2) Apply a soft potential to the rigid bodies in
case of overlap; use LJ/cut for the non-rigid and cross interactions 3)
Switch back to LJ/cut for all nonbonded interactions; time integrate the
non-rigid elements with NPT and apply positive pressure to ensure the box
does not blow up; time integrate the filler with rigid/nvt 4) Ramp the
pressure to zero 5) run long enough to see the pressure stabilize and take
note of the value (should be zero).

it will take a massive amount of compute time to redo the calculation.
if you want to have somebody else look into it, you should create a data file
for a smaller system, and also include a data file from a restart
taken during the last section of the overall run.

just taking a quick look at your system, i have one question: could it
be that you have some kind of "jamming" happen for your system?

axel.

7

I have created data files for much smaller systems and modified the LAMMPS script to shorten the computation time needed to observe the pressure discrepancy. The five data files all contain 500 crosslinker beads and 1000 polymer chains of length 20; The data files vary in the number of filler particles (rigid bodies) that are included (0, 2, 4, 8, 16), which corresponds to a volume fraction of 0-0.13. As Axel pointed out, a jamming transition could come into play at higher volume fractions, but I kept the volume fractions relatively low here to rule out the possibility of jamming phenomena. When running the attached lammps input file, I observe a systematic increase in the average pressure (during second half of the run, when barostat is set to 0.00) as the volume fraction of rigid bodies is increased.

As a further test, I omitted the ‘fix rigid’ command and instead used stiff harmonic bonds to maintain the shape of my filler particles. In this test, the pressure stabilized near 0.0 as desired. I suspect that something is going on with the fix rigid command.

generation_test_PS00007.inp (1.33 KB)

bulk_PS00007_no_filler_000.data (1.14 MB)

bulk_PS00007_2_filler_000.data (1.17 MB)

bulk_PS00007_4_filler_000.data (1.19 MB)

bulk_PS00007_8_filler_000.data (1.23 MB)

bulk_PS00007_16_filler_000.data (1.32 MB)

8

Hi Scott,

can you try either

  1. reduce Tdamp of fix rigid/nvt to lower values, such as 0.1, or

  2. switch to fix rigid instead of rigid/nvt

to see if there is any difference?

For the system you are simulating, the number of bodies is rather small and the thermostat coupled with the bodies (internal to rigid/nvt) would lead to big fluctuations in the kinetic energy of the bodies (hence global pressure). fix rigid/nvt (and rigid/npt) are well suited for thermostating a sufficiently large number of rigid bodies, so I would expect pressure get closer to the set value if you have more than 100 filler particles.

Cheers,
-Trung

9

Trung,

Thank you for the suggestions. I reduced Tdamp of fix rigid/nvt to 0.1 and there was no difference–still a consistent positive pressure of the same magnitude. This problem was originally detected in much larger systems, which included >100 rigid bodies; I reduced the system size as per Axel’s request. Although I understand what you are saying about the thermostat working better for systems with a larger number of rigid bodies, deviation from the pressure set value seems to grow with the number of rigid bodies.

Please find the attached spreadsheet, which contains the pressure data from the systems attached in previous emails (0-196 rigid bodies). I have plotted the average pressure in the x direction for systems of varying numbers of rigid bodies. I arbitrarily chose the x direction but pyy and pzz are nearly identical to pxx in all cases.

Is it possible that the barostat is not using the global pressure tensor in making box volume adjustments? It appears that the virial is only being calculated for the non-rigid bodies, which results in over-pressuring the system (and over-pressurizing grows with the number of rigid bodies), as suggested by the global pressure tensor output.

pressure_comparisons_lammps_users.ods (22.1 KB)

10

Is it possible that the barostat is not using the global pressure tensor in making box volume adjustments?

I’ll ask this another way. Is it possible the barostat (fix npt in your original post)

is not using the same pressure you are monitoring in your thermo output?

The fix npt doc page tells you the IDs of the temperature and

pressure compute it creates to calculate pressure used by

the barostat (temp is part of pressure). You can monitor

those temp/pressure IDs in your thermo output, as well as

the global ones that are default created/used by the thermo

output (what you get with the “temp” and “press” keywords).

See the section discussing this on the fix npt doc page.

If they are different, that could be the issue.

Steve

11

scott,

have you checked, whether the issue could be related to the length of
your time step?
your rigid bodies are fairly large and thus even small updates to the
rotational degrees of freedom can lead to unexpectedly large overlaps
with the surrounding atoms.
that would be consistent with your observation that the error grows
with the number of rigid bodies.

axel.

12

Thank you, Steve and Axel. I tried a smaller timestep–decreasing from 0.005 to 0.001 tau–and saw no change. I also started writing the output from the temperature and pressure compute that ‘fix npt’ creates for the barostat and thermostat. These temp/pressure values are identical to the global values created by default by the thermo output. The barostat pressure set value is zero yet it is showing a constant positive value in the output, so it seems that something is fighting the barostat.

13

I should also add that I tried the berendsen barostat with a Nose-Hoover thermostat (fix nvt) in place of ‘fix npt’, and the pressure stabilizes near zero:

fix 1 all momentum 100 linear 1 1 1 angular
fix 2 filler rigid/nvt molecule temp 1.0 1.0 2.0
fix 3 polymer nvt temp 1.0 1.0 2.0
fix 4 polymer press/berendsen iso 0.0 0.0 2.0 dilate all

14

I also started writing the output from the temperature and pressure compute that ‘fix npt’ creates for the >barostat and thermostat. These temp/pressure values are identical to the global values created by default >by the thermo output.

good - that’s the right check

The barostat pressure set value is zero yet it is showing a constant positive value in the output,
so it seems that something is fighting the barostat.

Well, the barostat knows nothing about your system,

or rigid bodies, etc. All it knows is there is a target

pressure and a current pressure, which in your case

(iso) are both scalars. It simply adjusts the box

volume as a function of the difference between the

2 scalars. As long as the actual pressure responds

as expected (increases as volume shrinks), it should work.

Things that can mess it up: bad timestep, bad time

constant in fix npt, a system that responds in an unusual

way (e,g. a tiny change in volume produces a huge spike

in pressure, or a non-monotonic response for some reason).

Note that with rigid bodies, what you should see is when

the volume shrinks, the COMs of the bodies get close

together. But the size of the bodies does not shrink

(b/c they are rigid). So if somehow a small volume

shrink is inducing overlap in the outer atoms of 2 nearby

rigid bodies, that could look like a huge pressure spike.

Also, how noisy the pressure is around your target

value can be hard to interpret, and hard for fix npt to track.

But that is just the way pressure works, esp in solids.

Maybe you can monitor the volume and the pressure

on a fine timescale. I can not imagine how the pressure

could consistently stay below the target and the volume

not continually shrink until the pressure went above

the target. Or vice versa.

I’m CCing Aidan (npt author/expert). Maybe he has

other ideas on how this could occur, or ideas on

fix npt options to try (chains, etc). Info on your

units and time damp constant, also a plot of your

pressure vs time with the target pressure as a line,

would be helpful. Are you doing anything else non-equilibrium,

like shearing the system?

Steve

15

Hi Scott,

I would suggest you first equilibrate the system in the NVT ensemble, i.e. using fix rigid/nvt and fix nvt for the polymers, for a certain period of time, e.g. 10000-50000 steps, depending on the system size, to see what the current average pressure is. After the first equilibration stage, you then switch from fix nvt to fix npt for the polymers, ramping the current average pressure to the set pressure, in this case being 0.0, and finally equilibrating at the set pressure.

When I tried that procedure for 4 filler particles, the average pressure during the first equilibration stage is already close to 0.0 (the initial configuration is very dilute), and I can switch right to the second equilibration stage at the set pressure of 0.0.

My guess for the incorrect pressures is that you use fix npt right at the beginning with Pstart = 0.2, which causes the system to be compressed too fast at the early stage- you can check this by writing out the trajectory. And when the number of the filler particles is increased, it would take longer than 500000 steps to really relax from Pstart=0.2 to the target pressure 0.0 and at the set pressure of 0.0 using Pdamp=2.0 time units. You can try to reduce Pdamp further (say, 0.1 and 0.01), as long as the simulation does not crash, to see if pressure would relax to 0.0 faster.

Hope it helps,
-Trung

16

Trung,

I have attempted a routine on the 4-filler system that is similar to what you suggested. I used fix deform to gradually shrink the intitial box volume (packmol generated) to a volume that is still much larger than the equilibrium, zero-pressure volume (dilute). This routine was performed over 5e05 timesteps and the pressure was near zero during the deformation. Then I switched to npt with pressure set to 0.0 and monitored the system for 5 million timesteps, and the system is showing the same issue. I have attached the pressure data.

This problem goes away when I use a Berendsen barostat insead of Nose-Hoover.

pressure_vs_time.ods (55 KB)

17

Hi Scott,

sorry I have no idea of what could be the cause and how to resolve the issue. In my test run with the 4-filler system (attached are the input script, log file and plot) pressure looks as expected.

Best,
-Trung

log.lammps (422 KB)

generation_test_PS00007.inp (1.41 KB)

pressure.pdf (155 KB)