[gmx-users] Virial from non-bonded forces

[Omololu Akin-Ojo]

Dear David,

Thanks for your response. It was very valuable. I looked in the codes to locate the place(s) where GROMACS obtains/accumulates the shifted forces and as far as I could see, GROMACS does this correctly for the short-range forces and re-distributes the forces on the virtual sites unto the real/atomic  sites before calculating the virial.
For the long-range forces (Ewald, SPME), however, one has to use all the possible (infinite number of) images but GROMACS does not do this; personally, I think it is impossible to do this.

There is, however, this important point: if the position of the dummy scales the same way as those of the real sites, then one does not need to re-distribute the forces in order to calculate the virial -- simply treat all sites in the same way. Thus, for dummy types 1,2,3, the virial calculated (even with Ewald) will be correct since r_d(L*r_1, ..., L*r_n) = L*r_d(r_1, ..., r_n) where r_d is the dummy site and r_1, ..., r_n are the real/atomic sites. For types 3fd,  3fad, 3out, and 4fd for which the position of the virtual site does not scale like those of the real sites, the virial calculated in GROMACS will be correct only for the short-range forces (since this is done correctly in the code) but not for the long-range ones. Note, however, that for rigid molecules, since the intramolecular bondlengths are constant, one can replace certain magnitudes |( )|  by constants in  3fd,  3fad, and 4fd. [In 3out replace by c/sqrt(|rij| * |rik|) ] so that r_d will
 scale just like any r_i or any r_ij. Thus, for rigid monomers (eg TIP5P), one again gets the correct contribution to the virial from long-range forces.
Fortunately, I do not know any flexible models out there that use 3fd,  3fad, 3out, and 4fd for the virtual sites.

Thanks for your time!

o.




________________________________
From: David van der Spoel <spoel at xray.bmc.uu.se>
To: Discussion list for GROMACS users <gmx-users at gromacs.org>
Sent: Tuesday, February 10, 2009 12:52:06 PM
Subject: Re: [gmx-users] Virial from non-bonded forces

Omololu Akin-Ojo wrote:
> Hi,
>  I have two main issues:
>  1. In Appendix B of the GROMACS user manual, there is a discussion of the calculation of the virial. The top of the second page of the appendix (just before B.1.2) states, "In a triclinic system there are 27 possible images of i, ...".
> My question is: In this case, there should also be 27 possible forces on particle  j due to i. Granted, only the force corresponding to the shortest i-j distance will be non-zero but when the forces are multiplied by \delta_i, (in Eq. B.11) the correct force (out of all the 27 possible ones) should be multiplied by the corresponding \delta_i.
> In other words, if we label each of the 27 possible \delta_i by n ((\delta_i)^n), we should do the same with the forces F_i  ((F_i)^n) and the sum in B.11 should also contain a sum over n.
> For ver. 3.3.1, I see that in src/mdlib/calcvir.c each of the different (\delta_i)^n is used but for fixed i, the same force F_i is used irrespective of the value of n.
> (See also, J. Chem. Phys. vol. 117, p2449 (2002), the 3rd line from the bottom of the 2nd column on the first page.)
The force on the non+central boxes are also stored separately in an arrya fshift. The displacement part of the position is thus taken care of separately.

>  2. How does GROMACS calculate the virial for charged virtual sites using the Ewald summation (or SPME)? One needs to redistribute the forces before calculating the virial but if one uses this approach one will need not only 27 possible \delta_i but, technically, an infinite number.
>  
The order of redistributing forces and virial calculation is important, and one should indeed do it in the order you describe. This is not completely obvious, since the position used for computing the forces are those of the virtual sites. I do not see why one would need to compute multiple deltas, but it could be that the contribution to fshift (see above) is computed before redistributing the forces, which would not be consistent in that case. GROMACS does reproduce e.g. the density of TIP4P and TIP5P water correctly, so it seems that there is not a large error if there is one at all.


> Thanks for your input.
>  o.
> 
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php


-- David van der Spoel, Ph.D., Professor of Biology
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:    +46184714205. Fax: +4618511755.
spoel at xray.bmc.uu.se    spoel at gromacs.org  http://folding.bmc.uu.se
_______________________________________________
gmx-users mailing list    gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php



      
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090213/872a4c5a/attachment.html>