[gmx-users] Bugs in calculation of virial

[Vincent Ballenegger]

Dear all,

I think I found three bugs in the way the virial (and hence the 
pressure) is calculated in gromacs.

1) Bug in long range correction term:
The virial in gromacs is defined as (eq. 3.18 of manual)
    tensor(Xi) = - (1/2) sum_{i<j} vect{r_ij} x vect{F_ij}.
The factor (1/2) has then be forgotten in equations (C.8) and (C.10) of 
the manual, as well as in the code in sim_util.c that computes this LR 
correction term (C.10).

2) Bug in long range correction term (again)
The correction (C.10) is to be applied to the scalar virial Xi. Gromacs 
adds this term to each element on the diagonal of the virial tensor, so 
a factor 1/3 is missing, since Xi = Trace(tensor(Xi)).
This bug was probably induced by the fact that the scalar pressure is 
defined as *one third of* the trace of the pressure tensor, so it's 
right to add the LR correction term (C.11) to each diagonal element of 
the pressure tensor.

3) virial and external electric field
When a uniform external electric field is applied on the system, the 
pressure should be computed from the *internal* virial, which excludes 
the contribution from the forces due to the external electric field. As 
far as I can see in the code, gromacs does not exclude this contribution 
(which can be big if the applied field is strong).
In the revision notes, it is written that the virial contribution of 
forces from an electric field has been removed in version 3.0.2. I 
compiled however versions 3.0.1 and 3.0.2, and both give identical 
results for the virial of SPC water in a strong electric field, so it 
seems that this bug has not been corrected in version 3.0.2 for some 
strange reason.

Regards,
Vincent Ballenegger