[gmx-users] Bugs in calculation of virial

[Berk Hess]

>From: Vincent Ballenegger <vcb25 at cam.ac.uk>
>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>To: gmx-users at gromacs.org
>Subject: [gmx-users] Bugs in calculation of virial
>Date: Tue, 20 Jul 2004 10:40:01 +0100
>
>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.
>

I noticed and fixed these virial problems (the errors were not
my work though) one month ago in CVS.

I have rewritten and moved the whole LR correction section
in the manual.

In the code there is only a problem with the virial correction,
the pressure correction is correct. The total error in
the virial correction due to both problems is a factor of 6.
This virial correction has no influence on the dynamics,
only the pressure correction is used. So all runs are
correct, only the reported virial is incorrect.
I have fixed this in cvs.

>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.

I made this update and now I looked at it again is is clearly wrong.
The virial call only calculates the shift part of the virial and
calc_f_el only adds to the force array. So the electric field force
is still part of the virial.
This can be solved by moving f_calc_el call back to before
the virial call and subtracting the electric force from the shift
force of the central cell.
Should I send you corrected code?

Thanks for reporting these bugs,

Berk.

_________________________________________________________________
Hotmail en Messenger on the move 
http://www.msn.nl/communicatie/smsdiensten/hotmailsmsv2/