
BarostatTable of contentsNo headers(This text was copied from bugzilla 14) Ramon Garcia Fernandez: 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 possible 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. Erik Lindahl: 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...
m_i d^2s_i/dt^2 = h^(1) f_i  m_i G^(1) dG/dt d s_i/dt s = coordinates in the box frame h = box matrix, transposed of Gromacs box G = h^t t Looking at Gromacs code in do_update_md and parrinellorahman_pcoupl and comparing to the equation above we see that the velocity is defined as: v_i = h ds_i/dt therefore, the updating of the coordinates should be dx_i/dt = d(h s_i)/dt = dh/dt * s_i + h * d s_i/dt = = dh/dt * h^(1) * x_i + v_i = = M * x_i + v_i (the last step assumes that dh/dt * h^(1) is symmetric) Michael Shirts: I've got a working version of the Anderson barostat (ParrinelloRahman is the generalization that allow nonisotropic changes), but it will require velocity verlet to be implemented in them main code to work correctly, as Erik mentioned. 