Project

General

Profile

Bug #14

Rahman Parrinello barostat not accurate for the anisotropic case.

Added by Ramon Garcia about 14 years ago. Updated over 12 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

The implementation of the Rahman Parrinello barostat does not match the one
published in M. Parrinello, A. Rahaman J. Appl. Phys vol 52 n 12 page 7182

For the case where only isotropic pressure is applied, the implementation is
correct. However, when there is an non isotropic stress tensor, Gromacs uses an
incorrect generalization that consists in replacing the pressure by the stress
tensor in the formula:

W box'' = (pressure_observed - pressure) box^{-1} * volume

This seems a natural generalization. In fact first paper anticipating the Rahman
Parrinello thermostat (PRL vol 45, n 14, pg 1196) suggested this extension. But
the paper (earlier mentioned) that developed the method gives a different formula.

A posible solution is to support isotropic pressure only, even when allowing box
deformations. This is done by DL_POLY. In fact, the review article of Nose and
Klein (Molecular Physics, vol 50, n 5 pg 1055) cited by Gromacs does not mention
how to implement a barostat to a stress tensor. Furthermore according to the
paper by Rahman and Parrinello, if a external stress tensor is applied, both the
hydrostatic pressure and the tensor have to be specified. So ref_p for the
anisotropic case would need 7 components, one for the isotropic pressure and 6
for the external stress tensor.

History

#1 Updated by Erik Lindahl about 14 years ago

Hi Ramon,

Thanks for the report. This will probably not be fixed right now, since Michael Shirts is working on a
complete overhaul of the Parinello-Rahman (and Anderson) pressure scaling.

The problem is that to conserve enthalphy we cannot reformulate the equation of motion from scaled to
normal coordinates like we do now. It also requires velocity verlet integration, which in turn means we
need the velocity constraint correction step in Lincs and Settle (Rattle). Too many things depending on
each other...

Nevertheless, I'll try to have a look at the article.

#2 Updated by David van der Spoel over 12 years ago

I'm closing this bug, it has now been listed as a development target at the gromacs wiki:

http://wiki.gromacs.org/index.php/Barostat

Also available in: Atom PDF