Kinetic pressure calculation
Closed this issue · 2 comments
I was looking at the calculations for ensemble properties of the system, and I found something confusing in regarding the calculation for kinetic pressure, defined as follows:
def get_virial(mol):
"""Clausius virial function for all atoms, force,s and coordinates.
Args:
mol (mmlib.molecule.Molecule): Molecule object with coordinate
and force data.
"""
mol.virial = 0.0
for i in range(mol.n_atoms):
for j in range(3):
mol.virial += -mol.atoms[i].coords[j] * mol.g_total[i][j]
def get_pressure(mol):
"""Update total pressure of a system of molecules with boundary.
Args:
mol (mmlib.molecule.Molecule): Molecule object with temperature,
volume, and virial data.
"""
get_virial(mol)
pv = mol.n_atoms * energy.kb() * mol.temp
pv += mol.virial / (3.0 * mol.n_atoms)
mol.press = kcalamol2pa() * pv / mol.vol
In equation form:
After spending an afternoon going through as much MD literature I could find, I've yet to see any paper that lists that N dividing the virial when calculating the pressure. Some examples:
Page 2: https://arxiv.org/pdf/1111.2705.pdf
Page 1: https://spiral.imperial.ac.uk/bitstream/10044/1/262/1/edm2006jcp.pdf
Page 5: http://phys.ubbcluj.ro/~tbeu/MD/C2_for.pdf
Not only that, but the literature is consistent as expected, and apparently in disagreement with this implementation. I am probably missing something, since I'm no expert in MD and the notation could be the culprit here. In any case, I'd like to know what's the problem.
Fixed in above linked commit.