[gmx-users] Bugs in calculation of virial
Berk Hess
gmx3 at hotmail.com
Tue Jul 20 12:20:14 CEST 2004
>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/
More information about the gromacs.org_gmx-users
mailing list