Npt/asphere pressure calculation

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

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

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