Npt/asphere pressure calculation

LAMMPS LAMMPS General Discussion
1

Dear all,

I am simulating a system of ellipsoids using Gay-Berne potential using an npt/asphere fix.

Checking through the code I noticed npt/asphere calculates pressure according to the total kinetic energy (translational + rotational). This yields a kinetic part of the pressure equal to ndof/3 * kbT/V, which is equal to NkbT/V for point particles (where N is the number of atoms), while it is equal to 5/3 N kbT/V for uniaxial LCs (or ellipsoids) and to 2 NkbT/V for biaxial nematics.

First I would like to ask if the paragraph above is mathematically sound?

I am aware that npt/asphere will calculate everything correctly given the selected keywords. The only problem is your desired pressure (the numbers after iso keyword) will not be the same as the “real” pressure (which in my opinion is calculated using only the translational degrees of freedom).

My way around it was to use the following commands: (lammps version 22 Oct 2020)

group spheroid type 1
variable dof equal count(spheroid)

fix 1 all npt/asphere temp 1.4 1.4 0.1 iso 4.0 4.0 0.5

compute_modify 1_temp extra/dof ${dof}
uncompute 1_press
compute nonrot all temp
compute 1_press all pressure nonrot

So I first uncomputed the pressure compute created by the fix npt/asphere and then computed pressure
again with the temperature compute for point-like atoms which include only the translational degrees of freedom.

Now the kinetic part of the pressure, which I defined as
compute mypress all pressure nonrot ke

c_mypress outputs in equilibrium approximately NkbT/V, which is the correct value.

So my question is if the definition for the pressure in npt/asphere is standard for MD of aspherical particles, or is my way physically the correct way?

Best regards

2

I made the same assumption when developing the MOLC model, i.e. to use the pressure computed with the translational temperature instead of the default pressure for the NPT/NPH time integrators. I would say that the rotational component of the temperature should not contribute to the pressure, but I, too, wondered why it was implemented in this way in the ASPHERE package.

These are the settings I currently use:

atom_style hybrid molecular ellipsoid
units real
# ...
compute temp_trasl     all temp
compute temp_rototrasl all temp/asphere dof all
compute press_trasl    all pressure temp_trasl
# ...
thermo_modify temp temp_rototrasl press press_trasl flush yes
# ...
timestep 10.
fix 1 all npt/asphere temp 300. 300. 1000. iso 1. 1. 10000.
fix_modify 1 press press_trasl
fix_modify 1 temp  temp_rototrasl

Here, I use the full virial and the translational temperature to compute the pressure. I haven’t tested your settings with the ke keyword, but I think this is incorrect as it will not include the force field contributions to the virial.

I hope this helps.
Otello

3

Dear hothello,

thanks for the answer and the link to the paper.

Yes, I think rotational part should never be included in the pressure (be it for point particles, spherical particles or aspherical particles). Although, I have seen 1 paper on LCs that takes (for some reason) the rotational part of the kinetic energy as well. When one derives pressure from virial theorem, one takes only translational kinetic energy.

Many books only state “total kinetic energy”, which is probably the source of all confusion, as people then think they should include the rotational kinetic energy as well. The term “total” means in this case the sum over all particles. Virial theorem is nicely explained in the book “Theory of molecular fluids” by K.E. Gubbins (Appendix E, page 602-615).

To answer you comment on the “ke” keyword. I only used this at the end in an “auxiliary” compute to output the value in the “thermo custom” command. I never actually use it. In the npt/asphere I use the “full” pressure (kinetic part + virial part).

Thanks again! Now I can sleep at night.

Best regards,
Entropic

4

I performed an NpT ensemble simulation of uniaxial Gay-Berne (3, 5, 2, 1) ellipsoidal molecules in a 3D system using LAMMPS (version 29Aug2024). The system parameters are:
N = 1000
T = 2.0
P = 5.0

The ellipsoids are defined and initialized as follows:

set          type 1 shape 3 1 1
pair_style   gayberne 1 1 2 4
pair_coeff   * * 1 1 0.2 1 1.0 0 0 0 4

To understand how pressure is being computed in the system, I included the following compute statements in the script. These were informed by comments from hothelo here.

compute temp_trasl         all temp
compute temp_rototrasl     all temp/asphere dof all

compute press_trasl        all pressure temp_trasl
compute press_rototrasl    all pressure temp_rototrasl

compute press_virial       all pressure temp_trasl virial

compute press_ke_trasl     all pressure temp_trasl ke
compute press_ke_rototrasl all pressure temp_rototrasl ke

compute ke_trasl           all ke

The following commands were used to average and output the computed properties:

fix AVE all ave/time 10 50 500 c_temp_trasl c_temp_rototrasl c_press_trasl c_press_rototrasl press_virial press_ke_trasl press_ke_rototrasl mode scalar ave running

thermo_style  custom step temp c_temp_rototrasl epair etotal press c_press_trasl c_press_rototrasl c_press_virial c_press_ke_trasl c_press_ke_rototrasl c_ke_trasl vol density 
thermo_modify temp temp_rototrasl press press_trasl flush yes

fix        2 all npt/asphere temp ${T} ${T} ${Tdamp} iso ${P} ${P} ${Pdamp} mtk yes pchain 2 tchain 4
fix_modify 2 press press_trasl temp temp_rototrasl

A typical thermo output looks like:

Step          Temp      c_temp_rototrasl     E_pair         TotEng         Press      c_press_trasl  c_press_rototrasl c_press_virial c_press_ke_trasl c_press_ke_rototrasl   c_ke_trasl       Volume        Density    
100000   2.0075482      2.0075482     -1.4597026      4.5599308      5.018496       5.018496       5.018496       3.9747296      0.62832055     1.0437664      3.6236645      3844.8151      0.26009053

Corresponding fix AVE output:

Temp_trasl         2.4148500600939311056
Temp_rototrasl     2.0107547671483891882
Pressure_trasl     4.9854189312010142032
Pressure_rototrasl 4.9854189312010142032
Pressure_virial    3.9437766225578632096     # virial component of the pressure
ke_trasl           0.62495933946304771389
ke_rototrasl       1.041642308643152548

Question 1: How come Temp_trasl is grater than Temp_rototrasl? In calculating Temp_rototrasl both the translational and rotational kinetic energy are used, where as in Temp_trasl only the translational kinetic energy is included and hence I think Temp_trasl should be smaller than Temp_rototrasl.

To check this is have used the press_ke_trasl and press_ke_rototrasl computes to calculate the kinetic energies (ke) (used the calulaing the temperature in press_trasl and press_rototrasl computes). However looking at the AVE outputs I am not sure if the kinetic energies (and hence the temperature) are calculated correctly or not.

Additional Consistency Check:
To verify the internal calculations, I implemented my own post-processing scripts based on the methods from:
1. Ricci et al. — referred to as the Matrix Method (MM)
2. Calderón-Alcaraz et al. — referred to as the No-Matrix Method (NM)
Using dump files from LAMMPS, I computed pair energy, virial, and pressure. Sample output from my code is:

#time_step, Upair_NM, Upair_MM, P_kin(=rho*T), P_vir_NM, P_MM_vir, P_NM, P_MM, Vol, Dens
100000   -1.4597023   -1.4597023    0.5201811    3.9747307    3.9747308    4.4949118    4.4949118 3844.8151469    0.2600905

It can be seen that the per-particle pair potential (-1.4597023) and virial component (3.9747308) of the pressure matches very well with the LAMMPS thermo output. However, the total pressure (kinetic term + virial = 4.4949118) does not match with any of the pressure calculated in LAMMPS thermo/AVE. I believe, this mis-match in P is due to the kenetic term (rho * T) used in calculating the total pressure.

Final question: How and what changes I need to make in the LAMMPS script to correctly sample the temperature and pressure?

NOTE: I have tried with the thermo_modify temp temp_rototrasl press press_trasl flush yes statement switched off too. However this only changes the corresponding components of thermo output different, while all other aspect of the simulation (output) remains same.

#Step          Temp      c_temp_rototrasl     E_pair         TotEng         Press      c_press_trasl  c_press_rototrasl c_press_virial c_press_ke_trasl c_press_ke_rototrasl   c_ke_trasl       Volume        Density
100000   2.4181946      2.0075482     -1.4597026      2.1639619      4.6030501      5.018496       5.018496       3.9747296      0.62832055     1.0437664      3.6236645      3844.8151      0.26009053

[1] Ricci, M., Roscioni, O. M., Querciagrossa, L. & Zannoni, C. MOLC. A reversible coarse grained approach using anisotropic beads for the modelling of organic functional materials. Phys. Chem. Chem. Phys. 21, 26195–26211 (2019).
[2] CalderĂłn-Alcaraz, A. et al. A Bidimensional Gay-Berne Calamitic Fluid: Structure and Phase Behavior in Bulk and Strongly Confined Systems. Front. Phys. 8, (2021).
5

Hi,

Please note it is generally very bad practice to post a reply to a very old thread, as it is difficult to keep threads organised under such circumstances.

Please start a new thread with your post instead – all of the LAMMPS calculations are actually internally consistent except for the translational-temperature pressure (you are missing a factor of two when including the rototranslational KE into the corresponding “rototranslational pressure”).

6

In addition to @srtee, please also include the complete set of files to reproduce your results. It looks suspicious that your averaged Pressure_trasl and Pressure_rototrasl are numerically identical.

1 Like
7

Thank you very much for the response.

I have created a new post with the same query.

@hothello Yes indeed Pressure_trasl and Pressure_rototrasl are numerically identical! I am not yet able to figure it out the reason for the same. In fact, I think I am looking for help to understand exactly this. One more thing to notice is that the KE calculated using the computes press_ke_trasl and press_ke_rototrasl are different (can be seen from the thermo output).

I have attached the infile and corresponding log files in the separate post mentioned above.

Thanks again.

1 Like