Pressure Calculation in npt/asphere fix of Uniaxial Gay-Berne Molecules

LAMMPS LAMMPS General Discussion
1

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 @hothello 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   1.9855451      1.9855451     -1.5981413      4.3555156      4.679503       4.679503       4.679503       3.6518803      0.61270336     1.0276228      3.5497711      3862.4141      0.25890544

Corresponding fix AVE output:

Temp_trasl         2.3966076465001076023
Temp_rototrasl     1.9995185764328464728
Pressure_trasl     4.999155980778550834
Pressure_rototrasl 4.999155980778550834
Pressure_virial    3.9593455670215988995    # virial component of the pressure
ke_trasl           0.62284063344093165693
ke_rototrasl       1.0398104137569483818

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.5981406   -1.5981406    0.5178109    3.6518812    3.6518812    4.1696920    4.1696921 3862.4140770    0.2589054

It can be seen that the per-particle pair potential (-1.5981406) and virial component (3.6518812) of the pressure matches very well with the corresponding values obtained in LAMMPS thermo output. However, the total pressure (kinetic term + virial = 4.1696920) 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 makes the corresponding components of thermo output different, while all other aspect of the simulation (output) remains exactly 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.368883       1.9855451     -1.5981413      1.9516298      4.2645836      4.679503       4.679503       3.6518803      0.61270336     1.0276228      3.5497711      3862.4141      0.25890544 

Temp_trasl         2.3966076465001076023
Temp_rototrasl     1.9995185764328464728
Pressure_trasl     4.999155980778550834
Pressure_rototrasl 4.999155980778550834
Pressure_virial    3.9593455670215988995    # virial component of the pressure
ke_trasl           0.62284063344093165693
ke_rototrasl       1.0398104137569483818

Here are the input script, log files and command used to run the simulation:

mpirun -np 4 lmp_gpu -in in.test.set.quat -var seed_vel 22222 -var P 5 -var T 2.0 -log log.npt.annealing -nocite -sf gpu -pk gpu 1 -var seed_quat 11111

log.npt.annealing (161 Bytes)
log.npt.P.5.T.2.00 (30.9 KB)
in.test.set.quat (4.8 KB)

[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).

2

Hi Deb,

Most of your results are consistent – the one thing I am scratching my head about is how the “translational pressure” and “rototranslational pressure” has come out the same.

Other than that, here is the accounting for the system:

For 1000 particles in 3D, there are 3 translational DOFs and 3 rotational DOFs per particle resulting in:

3000 translational DOFs and
3000 rotational DOFs in the box.

In LJ units, equipartition says that at a temperature of 2 units, each DOF should hold 1 unit of kinetic energy on average, resulting in 3000 units of translational KE and 3000 units of rotational KE, and thus in total 6000 units of KE.

Presumably you do actually have 6000 units of KE in total which results in temp_rototrasl reading a temperature of 2. However, your thermo output shows an average of 3.6 units of translational KE per particle (in LJ, thermo readouts are “intensive” or normalised per particle by default), which corresponds to a total of 3600 units translational KE and thus a temp_trasl of 2.4 (3600 times 2, then divided by 3000 translational DOFs).

The rotational DOFs are thus underthermalised – since I don’t usually do extended-particle MD I wouldn’t know why.

Now, the general kinetic contribution to pressure is 2KE/3V. The “total kinetic pressure” is 2*6000 units of KE/3*(just under) 4000 units of volume, which comes out to (just over) 1 unit of pressure. Similarly, the “translational kinetic temperature” of 3600 units of translational KE comes out as 0.6 units of pressure.

So the only mystery is how press_rototrasl and press_trasl are both returning the same pressure. That is extremely unlikely if they both have the same virial and different kinetic pressures. If you can print the virial contribution to both pressures (not just press_rototrasl) we can get to the bottom of that.

1 Like
3

Thank you for putting the number crunching so clearly.

Here is the virial contribution to both the pressures, namely press_trasl and press_rototrasl. And yet again the virial contribution to both the pressure computes are same, indicating the mystery of how the kinetic contribution to the pressures are calculated.

Temp_trasl                2.3966076465001076023
Temp_rototrasl            1.9995185764328464728
Pressure_trasl            4.999155980778550834
Pressure_rototrasl        4.999155980778550834
Pressure_trasl_virial     3.9593455670215988995
Pressure_rototrasl_virial 3.9593455670215988995
Press_ke_trasl            0.62284063344093165693
Press_ke_rototrasl        1.0398104137569483818
ke_trasl                  3591.3165582804113001

The two computes for the virial calculation are set as:

compute press_trasl_virial       all pressure temp_trasl virial
compute press_rototrasl_virial   all pressure temp_rototrasl virial

The mystery: If trasl_virial and rototrasl_virial are identical and temp_trasl and temp_rototrasl are different, then how come press_trasl and press_rototrasl are still same!

4

I’ve figured it out: that behaviour results from the following line in your script:

fix           2 all npt/asphere ...
fix_modify    2 press press_trasl temp temp_rototrasl

This does not merely set the NPT integrator’s compute temp to temp_rototrasl – it also sets the pressure compute’s compute temp to temp_rototrasl.

This is documented behaviour, and I hope you will agree that it is a sensible default for ensuring that the barostat follows a consistent temperature with the thermostat (although interestingly if you had happened to put the temp keyword before the press keyword in fix_modify, you would have made more sense of the results immediately):

If both the temp and press keywords are used in a single thermo_modify command (or in two separate commands), then the order in which the keywords are specified is important. Note that a pressure compute defines its own temperature compute as an argument when it is specified. The temp keyword will override this (for the pressure compute being used by fix npt), but only if the temp keyword comes after the press keyword. If the temp keyword comes before the press keyword, then the new pressure compute specified by the press keyword will be unaffected by the temp setting. (emphasis added)

Thus, after this command, the particular pressure compute press_trasl uses the “rototranslational KE” in its pressure calculation – this makes it identical to press_rototransl as you have observed. If you issue a new, different pressure compute coupled to just the translational KE you should observe the expected results.

4 Likes
5

Dear Shern,

Thank you for your brilliant analysis, I was actually not aware of the behaviour of fix_modify.

The point that I have raised in my previous topics on this subject is that the correct target pressure to use in fix npt/asphere should be the one computed only on the traslational DOF, overriding the default behaviour. In my scripts with finite-size (a)spherical particles, I only set fix_modify NPT press press_trasl, as the default temperature computed by the fix is usually the correct one[1].

There is only one small caveat to your analysis.

For uniaxial ellipsoids (i.e. \sigma_x=\sigma_y), the Gay-Berne potential does not produce any torque on the XY plane, therefore the rotational DOF per particle are 2. With the following modifications, the kinetic energy is correctly thermalised to an average value of 5 LJ units per particle:

compute temp_trasl         all temp
compute temp_rototrasl     all temp/asphere dof all
compute_modify temp_rototrasl extra/dof $(atoms) # Subtract 1 DOF per particle.
compute press_trasl        all pressure temp_trasl
compute ke_trasl           all ke
compute ke_rot             all erotate/asphere
variable ke_tot equal c_ke_trasl+c_ke_rot

thermo_style  custom step epair etotal temp c_temp_trasl press c_press_rototrasl ke 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 temp temp_rototrasl press press_trasl # Both temperature and pressure are different from the default ones.

# Average T and P.
fix AVE all ave/time 10 50 500 c_temp_trasl c_temp_rototrasl c_press_trasl c_press_rototrasl c_press_virial  c_ke_trasl v_ke_tot mode scalar ave running

Which yields the following values:

Temp_trasl         2.0006614894864211962
Temp_rototrasl     2.0000857802874287827
Pressure_trasl     5.0000040701043726443
Pressure_rototrasl 5.3733741644110919466
Pressure_virial    4.4409510439032864326
KE_trasl           2997.9912419954016514
KE_total           5000.2144507185712428

PS I really like visualising thermodynamic phases of prolate ellipsoids :slight_smile:

image
image800Ă—831 96.9 KB

  1. but not in this case, see later. ↩︎

3 Likes
6

Dear @srtee and @hothello

Thank you so much for time and effort that you have put on the issue. You guys deserve a treat. The issue is absolutely solved now. Each and every value of different measurements are accounted for.

One great outcome of this discussion is @srtee 's find related to the importance of the order of temp and press compute in the fix_modify of npt/asphere.

@hothello Of course the reduction in the DoF for the Uni-axial GB particles is the cherry on the top.

compute_modify temp_rototrasl extra/dof $(atoms)

Extra Note: I was actually working in a system of GB ellipses in two dimensions (2D) and hence I just just want to add here (for those who would be interested in 2D) that: The DoF for the 2D case is well accounted for in the npt/asphere fix and temp/asphere compute. So, no need to add any compute_modify statement for 2D simulations.

I have found couple of more things about the 2D simulations using LAMMPS that I will post later, again for those who would be interested.

1 Like
7

Here is the post related to simulation of GB ellipses in 2D

I have found couple of more things about the 2D simulations using LAMMPS that I will post later, again for those who would be interested.