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