tmpchem/computational_chemistry

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:

image

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.