[gmx-users] Virial from non-bonded forces

[David van der Spoel]

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