Gromacs

Barostat

    Table of contents
    No 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...


    Ramon: Although there is some ongoing work to enhance the barostat with Gromacs, the current barostat can be made more accurate with a painless change. In do_update_md, modify the updating of the coordinates to include the change of the box. In fact, I believe that not doing so was a misunderstanding of the equations. The differential equation of the coordinates, in the box reference frame is:

     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 (Parrinello-Rahman 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.

    Page last modified 13:49, 14 Sep 2009 by rossen