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

@jtclemm has already given some of the answer above, but to be explicit, you need to add

compute temp_sphere all temp/sphere
compute press_sphere all pressure temp_sphere
fix_modify 9 temp temp_sphere press press_sphere

to your script and see if that gives you the results you want.

By the way, when you make statements like

please be explicit about how you have measured pressure. In any given simulation you will have both the thermo pressure compute and the barostat’s pressure compute, which may have slightly different settings, and you may also have set up your own different pressure compute to report pressure. You may have told us earlier how you are measuring pressure but it’s worth repeating – especially since you will be rapidly adjusting your scripts as you try different things.

Hi srtee
many thanks for your respones and help. my apologies for any confusions.
many thanks for your suggestions for the commands. althought i am not familar about these commands, i have tried to add this to my script and see what would happen. unforunately, i got:
ERROR: Could not find fix_modify pressure compute ID: press_sphere (…/fix_press_berendsen.cpp:494)
Last command: fix_modify 9 temp temp_sphere press press_sphere
I think i may need to study these compute temp_sphere and fix_modify firstly. i can see temp is related to temperature. if the simulation has nothing to do with temperature, should we still have this?

the way i measured pressure is:
compute 7 all pressure NULL pair
based on this command, compute pressure command — LAMMPS documentation

the equation is illusrated from documentation. and the pressure calculated from “compute_pressure” is what i want in my simulation. This pressure can also be calculated by sum(fx*lx)/v (i.e. for x direction) based on contact force dump files. they are same in my script and simulation.

But actually, what i am really care about is how can we achieve the target pressure (i.e. from compute_pressure and also can be calculated from contact dump files). since press/berendsen cannot really help if bpm package is included.

So in other words, if we assume that pressure from “fix press/berendsen” is different compared to that computed from “compute pressure”. then how can we achieve and hold a desired system pressure (i.e. taht from compute pressure) if we do not use “fix press/berendsen”

This error could only happen if you haven’t defined compute press_sphere all pressure .... If you have defined a compute with that name, LAMMPS should find it.

Your simulation always has a measurable temperature, even if you are not controlling the temperature. (If there is zero kinetic energy, the temperature is 0 K.) As the documentation states, LAMMPS’s standard calculation of pressure includes this kinetic, or ke, component.

As the documentation states, a barostat will create its own pressure compute, and that pressure compute will in turn create its own temperature compute, even if there is no system thermostat. This all happens unless fix_modify is used to “replace” the default pressure compute with a user-specified alternative.

To be explicit (and again this is in the documentation!) the barostat fix 9 all press/berendsen will create the following computes:

compute 9_temp all temp
compute 9_press all pressure 9_temp

The barostat fix 9 will then implement equations of motion such that the pressure compute 9_press will have the target value (say 30kPa). Now, if you have another pressure compute

compute myPress all pressure NULL pair

I hope you can see why myPress will be systematically lower than the target pressure. Furthermore, since the barostat fix 10 all nph/sphere will create the following computes instead:

compute 10_temp all temp/sphere
compute 10_press all pressure 10_temp

I hope you can see that in such a simulation both computes 9_press and myPress would be systematically lower than the target pressure.

And, on the other hand, a simulation with the barostat and barometer

fix myBeren all press/berendsen ...
compute myPress all pressure NULL pair
fix_modify myBeren press myPress

will be one where the value of myPress will match its target on average. In that simulation, any pressure compute which includes the kinetic contribution will systematically read a higher pressure than the target.

Note that in all of these posts (I think) I have not used the phrase “the pressure”, because LAMMPS doesn’t know what “the pressure” is either – it simply adjusts a fix so that a compute has the specified average value. It is up to you to know which pressure compute setting is physically correct, and which is not.

Hi Srtee
Many thanks to your patient reponse and help. and i started to understand your comments. the problem could be that temp and related kinetic energy. in my case, i can see that Temp is 2.32e13, which is really high. in contrast, for my compute pressure command, i used NULL, which suggests that the kinetic energy is not included in the pressure. This may lead to the inconsistency. but it cannot really let me see why i do not have this problem for pure sand and in previous version of lammps.

and you tried to teach me to use fix_modify, which replace the default pressure compute with my specific “compute mypress all pressure NULL pair”. i have tried this, but always got the error:
ERROR: Could not find fix_modify pressure compute ID: mypress (…/fix_press_berendsen.cpp:494)
but its weired since i have defined the pressure compute ID.

i attached the error log, we can identify that red lines, i have defined the compute “mypress” . but i still got the error when i tried to use fix_modify. do i make any mistakes for that?

many thanks for your response. it is really helpful.
i will also try to control the tempature as well. try to be similar to that of normal temp and see what will happen.

Yup, there was a bug in the fix_modify code for fix press/berendsen. See this PR: debug fix_modify press for press/berendsen by srtee · Pull Request #3826 · lammps/lammps · GitHub it should be extremely easy for you to fix this in your own source code and recompile. Thanks for your detailed documentation.

Hi srtee
Many thanks for your help for that. i have revised the source code based on your modifcation. and this time, fix_modify work perfectly. i got 30 kPa for fix press/barendesn, compute pressure NULL and also the pressure calculated from contact force dump files. everything looks wonderful now. thank you very much for that. it is really helpful.
i am planing to control the temp as well. i hope to identify the exact reason why i used to get disconnect before and will update this if i identify any.
hope you have a nice day

best regards

Hi Jtclemm
my apologies for bothering you.
it just there is one thing, which potentially could be a bug and i therefore want to illusrate this.
for the latest develop branch, saying 15 Jun 2023
i got the error
ERROR: illegal bond bpm command, invalid argument break/no.
when i give the bpm/rotational break/no

I can see break/no was added in the latest feature lammps version.So do we delete this break/no or something wrong about this?


please find attached error log.

many thanks for your time and help.

I think srtee hit a home run with the advice he gave. Only one additional comment:

Seems odd not to include the contribution from bonds when using the BPM package.

The syntax was unified in a recent pull request, Misc minor patches/features in BPM/Granular packages by jtclemm · Pull Request #3807 · lammps/lammps · GitHub. Now the option is selected as “break no” or “break yes”. Hope this helps.

2 Likes

Hi jtclemm
many thanks for your response and help.
and also thanks for your suggestions for that contribution from bonds. based on your suggestions, i
used new compute pressure with extra keywords bond.
compute mypress all pressure NULL pair bond
and now i can see you are definitely right. The fact that i see different values of the pressure is due to the bond (inside a “clump”) interaction, which can also contribute to the virial.
But it seems this consideration may be really weired to my study aims, my priority is to create clump and consider this clump as single particle. its hard to consider that interaction within one grain can contribute to the pressure of system, is this really odd not include the contribution from bonds to the system pressure. Will this virtual contribution significantly affect the behaviour of system?

many thanks for your message for that break no but. i still got the error when i used break no.

please find attached image for your reference. does i make any mistake for that?

hope you have a nice day

best regards

Good, glad to know it works.

Nope, it looks like that was my mistake. I forgot to increment the argument parser twice - once for each keyword. The patch is in this pull request: Bugfix and documentation corrections for BPM+Granular packages by jtclemm · Pull Request #3827 · lammps/lammps · GitHub. Thanks for catching it.

Hi jtclemm
my apologies for bothering you again, since you are the expert for the Granular package, if you do not mind, may i get some idea from you about this?
i currently got weird results for compute stress/atom results and I wrote a new topic for that. see

is this a potential bug or I made something wrong?

many thanks for your help. hope you have a nice day

best regards
deyun

Hi jtclemm
Please accept my apology if this message bothers you.
if you do not mind, may i ask some questions related to the BPM package since you are definitely expert in this field.

as you mentioned before,
the contribution from bonds will contribute to the “pressure of the system” as well as the “particle stress of atom” .

I am therefore wondering whether we can directly turn off or ignore this virtual contribution of bonds to the “pressure of the system” and the “particle stress of atom”.

Since i am trying to use bpm package to create clump, and consider this clump as a “single particle”, i am trying to look the interaction between this clump and other particles. it is therefore not suitable to consider this "internal stress contribution " from the generated clump.

thanks to the help of srtee. i am able to use fix_modify to control system pressure and ignore the virtual pressure contribution from the bonds.

but now i am sure that the atom stress computed from " compute stress/atom command" also have the contribution from the bonds, which is not what i want in my case.
So is that possible that we just ignore the virtual stress/pressure contribution from the bonds in a bpm package. or is this any way to fix_modify this “compute stress/atom”?

hope you have a nice day
best regards

I’m not sure I’d call it a virtual pressure, contributions from bonds represent real stresses due to elastic deformation. But the doc pages for compute stress/atom describes how to control which virial components contribute to the stress.

Hi Jtclemm

many thanks for your reponse and explanation for this question, it is really helpful. i will be careful when i used this bpm package. hope you have a nice day

best regards