Kinetic pressure calculation

[mariohm1311]

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.

[tmpchem]

Hi Mario,

Thanks for pointing out that discrepancy. I would certainly treat this feature as experimental and not draw any important conclusions from it's value. As you indicated, the implementation appears to be incorrect.

This is somewhat legacy code as I spent some time testing out this feature for dynamic pressure-based boundary condition volume resizing, but it never really worked well so I didn't really discuss it in the MD chapter. Looks like I should have taken it out completely, as is a good software engineering practice for untested alpha features.

-Trent

On Thu, Feb 15, 2018 at 11:07 AM, Mario Hernandez notifications@github.com wrote:

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: image] https://user-images.githubusercontent.com/13200607/36275441-22d5cc68-128b-11e8-85dc-44d8b660bb16.png

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 http://url Page 1: https://spiral.imperial.ac.uk/bitstream/10044/1/262/1/ edm2006jcp.pdf http://url Page 5: http://phys.ubbcluj.ro/~tbeu/MD/C2_for.pdf http://url

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.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/tmpchem/computational_chemistry/issues/2, or mute the thread https://github.com/notifications/unsubscribe-auth/AGMmkuE7-rESixJe2U7_d3M1sIlu4GNGks5tVICPgaJpZM4SHWce .

[tmpchem]

Fixed in above linked commit.