A disconnect between "compute pressure" and "press/berendsen" when bpm package is used in 2023 feature lammps

I have identified another issue that when i used bpm package to create clumps for simulations in this latest version.
the mean stress of the system calculated from “compute pressure” is different compared to my desired stress used in “press/berendsen”. it only happen for clumps created in bpm package in LAMMPS Feature Release 28 Mar 2023.

i got same mean stress of system between “compute pressure” and “press/berendsen” for pure sand (only from granular package) in both LAMMPS Stable Release 23 Jun 2022 and LAMMPS Feature Release 28 Mar 2023.

i also got same mean stress of system between “compute pressure” and “press/berendsen” for pure clumps created by bpm package in LAMMPS Stable Release 23 Jun 2022, but not in this latest version.

for details,

i have tried to use bpm package to create clumps and use press/berendsen to obtain a desired stress like

fix         9 all press/berendsen x 30000 30000 10 y 30000 30000 10 z 30000 30000 10
fix	        3 all nve/bpm/sphere

it worked perfectly in 23 Jun 2022 stable lammps, I used “0 30000 1000” at very begining, until i got the close stress value (checked by compute pressure), i would restart using “30000 30000” 10 so that i got almost mean effective stress of system is just 30000. it works very well in 23 Jun 2022 stable lammps.

but in LAMMPS Feature Release 28 Mar 2023. i always got a big gap. for instance, in my case, the desired stress is 30000, however, the simulation will always give value of 22800. After a long time simulation. it seems that in current code, they believe that they reached the target stress, but not really reach the target illusrated by press/berendse (30000)

Since it only happen for this 2023 feature lammps and only happen for clump created by bpm. could this a potential logical problem?
is this anything i made mistake? its weired, i use 100% same input and coordinates but got two distinct things from two verision of lammps.
i have attached related files for your reference
@jtclemm is this right way to let the developer identify this potential problem?
rubber_sand100small12.data (8.2 MB)
check.in (2.8 KB)
check_mean_stresses_iso.txt (2.8 KB)
the restart is attached in google drive

fix nph/sphere also has the similar problem. it can not be used to reach the targe stress when created bpm is used

This is probably just due to changes/bug fixes to the BPM package between these two releases. June 2022 was right after the BPM package was added to LAMMPS, missing any subsequent patches and edits. I recall that bond style bpm/rotational did not originally include tangential forces in the stress calculation. This may just be a consequence.

Hi jtclemm
many thanks for your response and help. So do you means that we have fixed bugs for this bpm package? but how could we deal with such a big gap between the pressure calculation between between “compute pressure” and “press/berendsen”?

There is aother issue for pair_style granular, which may show some problems.
since in lammps, we can always use the force of each contact point and its corresponding branch vector to calculate the mean pressure of the system.
in x direction, the mean pressure can be calculated by the sum (fx * lx) ./ volume of system. I tried to use this method to calculate the mean pressure of the system based on the dumped contacts force files. But it always give the wrong result for the estimation of system pressure.

compute     2 all pair/local dist dx dy dz force fx fy fz p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 cutoff radius

This method works well in the 2022 version and a modified lammps version. I used 100% same input script, coordinates and matlab code, everything is same. but cannot got right answer in 2023 version. i may assume that the dumped forces and branch vector may have some potential problems in this 2023 features. unfortunately, i am not good at c++. I may try to look the code and maybe try to see the reason.

hope you have a nice day!!

best regards

2023.contact (7.7 MB)
2022.contact (3.7 MB)
i attached these two dumped contact files for your reference. the 8-10 column are fx, fy and fz. 18-20 are lx ly and lz.

Hazarding a guess, you might have different pressure computes with different contributions to the stress tensor. Note that fix press/berendsen does not directly calculate pressures. I suggest reading the documentation page for fix press/berendsen in particular about issuing compute pressure commands as that may clarify your question. Hopefully that helps resolve this issue.

Thanks for being careful and double checking everything. Just to confirm, does your pressure calculation include any other contributions such as the virial contribution from bond forces or the kinetic energy which could account for an inconsistency? If you can create a minimal system (ideally with only a few contacts and no bond forces etc.) that demonstrates an inconsistency I can take a look at it.

Pair granular recently underwent a large overhaul which fixed many bugs but could have certainly introduced a few new ones. We did extensive testing, but I wouldn’t be surprised if a bug slipped through on the single() method (which is used by compute pair/local) since that is a relatively underutilized feature. Since the new code is much cleaner, the hope is that it’ll be easier to find and patch any new bugs going forward.

Barostats work by measuring the system’s pressure, and then adjusting system properties each timestep to bring the measured pressure in line with the target pressure. This is explained in the manual for (for example) fix press/berendsen under the section starting:

This fix computes a temperature and pressure each timestep. To do this, the fix creates its own computes of style “temp” and “pressure” …

So to start understanding your problem you should work out what pressure your barostat is reading – if your barostat cannot measure the system’s pressure correctly then obviously it will not control the system’s pressure correctly. You can either add the Berendsen barostat’s pressure compute to your printout line:

fix 2 print " ... $(c_9_press)" ...

or directly assign your preferred pressure compute to the barostat using a fix_modify command.

Hi jtclemm
many thanks for your response, i will try to check the documentations. so it may due to the problem o fix press/berendsen. so the pressure computed by compute 7 all pressure NULL pair is right one?
i am also not really understand this point, i cite:
I recall that bond style bpm/rotational did not originally include tangential forces in the stress calculation. This may just be a consequence."
to my own understanding, in my case, bonds are used to create clump, the new clump will have interaction with other clump. but all force or stress induced by bonds within a clump should be something like internal force, why this could affect the system pressure?

Hi jtclemm
many thanks for your reponse for this calculation of stress.
this problem is checked by pure sand. i created a toy case with 500 particles.
same input script and coordinates are used in both 22 and 23 version lammps. right in 22 version, wrong in 23 version.
toy case.zip (1.1 MB)

please find attached related documents for your reference.
many thanks for your contribution and hope you have a nice day

best regards

hi srtee
many thanks for your suggestions. it is really helpful. i will have a look at very begining to see what the barostat is reading in my system .

hope you have a nice day

best regards


It sounds like you are asking LAMMPS to calculate two different pressure tensors. I cannot say which is the right definition for your research application. Srtee gave you some good advice for comparing the two tensors.

Bond style bpm/rotational includes normal forces and tangential forces. When the BPM package was first written, only normal forces were included in the stress calculation. It wasn’t until later that methods were added to include contributions from tangential forces. The fact that you see different values of the stress in the two versions may just be this change.

Forces from both pair (between two ‘clumps’) and bond (inside a ‘clump’) interactions can contribute to the virial. The documentation page explains how these two types of forces can be used to compute stresses and pressures.

In addition to a test system, you also need to explain what you calculated, what result LAMMPS gave, and why LAMMPS is incorrect. While there may indeed be a bug, it isn’t clear to me that your issue isn’t just due to the fact that you are using different definitions of the pressure. I also would prefer plain text files.

You may also consider filling out a formal bug report detailing the problem on LAMMPS’ github page.

Hi jtclemm
many thanks for your reponse, and my apologies for confusing statement.
this time i tried to have a much small coordinates and try to explain this step by step:
firstly. i used

compute     2 all pair/local dist dx dy dz force fx fy fz p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 p12 cutoff radius

to obtain the force for and branch vector for each contact, then use

dump	    1 all local  ${interval_dump} psd30_contact/dump*.contact c_3[1] c_3[2] c_2[6] c_2[7] c_2[8] c_2[18] c_2[19] c_2[20] 

to dump what i want, in this case, c_2[6] is the force in x direction for each contact, saying fx.
c_2[18] is the branch vector in x direction for each contact, saying lx. I just the X as a example, everthing is same in Y and Z directions.
Then sum (lx*fx)/ volume of box is the mean pressure of the system in x dirction.
see this image for reference.

this result makes sense, the pressure obtained from compute pressure is same to the pressure calculated by these contact points and also same to the desired pressure (given in press/berendsen). its obtained in 2022 version.

then in 2023 version.
i did the same thing,see

in this case, the same mathod give pressure of 321 Pa, this is not right. we should obtain a pressure of 30000 which is same to the desired pressure (given in press/berendsen) and also same to the pressure obtained from compute pressure.so i assume that compute pair/local may give wrong values.

check.zip (114.1 KB)
for the github page. i am not really familar about this thing, i will try to learn how to use this one and try to filling out.

thanks for your suggestion and contribution. hope you have a nice day

Thank you very much, this is much easier for me to look at. Notably, I can see it’s not the pressure compute but rather the output of pair/local. For pair 2-18 fx is 2.65565 in the 22 version and 0.041776 in the 23 version. That really helps clarify the problem and should be much easier to chase down. In particular, it looks like the distance r between particles is ~0.01, which makes me think this is another instance of a missing factor of 1/r. I can try and take a look at it tomorrow.

It does look like there was a missing factor of 1/r. Here’s a branch that includes the patch: GitHub - jtclemm/lammps at BPM

If you’re able to download and run it, please let me know whether or not you still see any differences. Thanks again for the thorough tests.

Hi jtclemm

Many thanks for your response and help. i will try to download this new bpm and try to check the consistenceny.
many thanks for your response and contribution.

hope you have a nice day

best regards

hi jtclemm
many thanks for your contribution. this problem seems fixed right now. the pressure calculated from contact forces is now same to that identifed from compute pressure.
Many thanks for your help and reponse for that.
will this also be merged into the develop branch?

best regards


Great, thank you for confirming it works. It’s very helpful to have someone carefully checking the code.

Yes, I’m gradually building up documentation edits and bug fixes in that branch and will probably submit a request to merge it into the develop branch in around a week.

hi jtclemm

my apologies for bothering you, Following the question about the disconnect between “compute pressure” and “press/berendsen”.
i have tried to use the post-contact force data to recalculate the pressure of the system and now confirms that the “compute pressure” is same to the pressure calculated from post-contact force dump files.
So, as you stated before, which could be right, that fix press/berendsen may refer to a different stress tensor. and the same thing can be identified for fix nph/sphere. there is a disconnect between “compute pressure” and “fix nph/sphere”
Then there is a new question? how can we really reach a target pressure if bpm package is applied.
as shown in this documenation, we have two methods. however, both Both fix nph and fix press/berendsen cannot really help to maintain the target stress once the bpm package is used. do you have anything new or different to keep a certain pressure of system?

i do have a idea, i can adjust the values applied in fix press/berendsen and check the pressure computed from “compute pressure” and then error and trial until i got the right answer, does this make sense?

According to the documentation, fix nph/sphere uses:

compute fix-ID_temp all temp/sphere
compute fix-ID_press all pressure fix-ID_temp

while fix press/berendsen uses:

compute fix-ID_temp group-ID temp
compute fix-ID_press group-ID pressure fix-ID_temp

Does this account for the inconsistency you see? Compute temp/sphere includes rotational kinetic energy.

Unfortunately there’s no version of fix nh that supports atom style bpm/sphere. However, fix press/berendsen works since it doesn’t integrate trajectories and can be used with fix nve/bpm/sphere.

Hi jtclemm
many thanks for your reponse and help.
it seems my statement leads to confusion. I would like to show a example to explain my confusion.
in a simulation with bpm package. i want to reach a isotropic compression saying 30000 Pa.
for the latest version of lammps. i used

fix         9 all press/berendsen x 0 30000 1000 y  0 30000 1000 z 0 30000 1000

it worked well before the latest version of lammps. my system reaches the target pressure of 30000 kPa. But currently. this system will reach a pressure of approximately 23000 Pa and then would not increase to 30000. “press/berendsen” believe it has reached the target pressure, but not really.
Pressure of 23000 Pa is confirmed by “compute pressure” and calculation from post-contact force dump files. So this 23000 Pa should be the real pressure of system. “press/berendsen” did not work effectively.

as you stated that the update of bpm package may lead to this issue. but there is new question. in the latest version of lammps. how could we reach a target system pressure if “press/berendsen” cannot reach the target value? fix nph/sphere can also not work. This suggests that currently, we do not have the suitable commands to let a system reach a desired pressure once bpm package is empolyed? could this be a potentially a problem?
by the way, if bpm package is not included, “press/berendsen” works well in the latest version.