[gmx-developers] virial

[Jason de Joannis]

I am still baffled, can someone help? I have even gone as far as to try 
computing the fshifts myself, but again I am getting different numbers.

I don't know enough assembly to decipher what is Gromacs is doing
differently than I am.

/Jason

Quoting Jason de Joannis <jdejoan at emory.edu>:

> Let me revive an old topic concerning the fshift correction
> to the virial. This term is added up in the calc_vir function and
> I believe it is equivalent to the second term in Eq(2) of the 2005 
> Gromacs paper or Eq. B.11 of the manual. Since that term is
> rather simple it should be easy to check with the following snippet of code:
>
>  clear_mat(virtest);
>  for (n=start; n<start+homenr; n++)
>    for (i=0;i<DIM;i++) for (j=0;j<DIM;j++)
>      virtest[i][j] -= 0.25*box[i][i]*graph->ishift[n][i]*f[n][j];
>
> Unfortunately I am not getting the same number for the trace of the
> virial tensor. Any suggestions would be appreciated.
>
> /Jason
>
>
>> Quoting Erik Lindahl <lindahl at cbr.su.se>:
>>
>>> Hi Jason,
>>>
>> Hi.
>>>> On Nov 21, 2006, at 9:41 PM, Jason de Joannis wrote:
>>>>
>>>> From the manual and the code, it looks like Gromacs uses the atomic
>>>> virial. I think that is besides the point.
>>>>
>>>> Why does Gromacs use f[i] which (I think) contains all bonded  and 
>>>>  nonbonded forces? After all the derivative resp. volume  scaling  
>>>> vanishes for some of the bonded forces.
>>>
>>> Forces are forces are forces, which give rise to pressure. There  
>>> is  no fundamental physical difference between bonded and  
>>> nonbonded  interactions.
>> Hmm. I generally think from the partition function rather than  
>> mechanically. Anyways, why do you have Appendix B.1.4 which treats
>> the virial for harmonic bonds explicitely?
>
> That could be an old remnant - currently all interactions are treated 
>  the same way.
>
> There is one slight technical difference in how we treat bonded vs.  
> nonbonded interactions, though, and this is the reason for the fshift 
>  term. For bonded interactions we make the molecule continuous and  
> just calculate all forces without having to worry about periodicity.  
> For nonbonded interactions the "shift" depends on the particle we're  
> interacting with, so that's a bit more complicated. If you look at  
> the JCC gromacs paper in 2005 we have a small figure there explaning  
> it better than I can do in words here!
>
>>>
>>> Gromos (and some other codes) use whole-molecule-virials to  reduce 
>>>  the effect of some numerical errors, but there are other  
>>> drawbacks  with this - primarily performance.
>>>
>>>>
>>>> Also, the manual describes the SHAKE virial but what about LINCS?
>>>
>>> We use "Shake virial" as a generic name for virial contributions  
>>> from  constraints, so it's identical for SHAKE and LINCS.
>> Okay, thats true. I can see real numbers when I print out  shake_vir. That
>> means I will have to go inside LINCS and somehow zero the unwanted
>> virial terms.
>>>
>>>>
>>>> Also, where can I find a little more info on the fr->fshift   
>>>> contribution.
>>>> On the mailing list archive someone explained that this comes  
>>>> from  forces on the 27 unit cells. Is that explained in the  
>>>> manual or  somewhere?
>>>>
>>>
>>> It's mentioned in the manual where we derive the expression for   
>>> single-sum-virial. The performance advantage of the single sum is   
>>> also the reason why we prefer the atomic virial, btw!
>> Well I must need a stronger prescription on my contact lenses because
>> I have been looking closely at that section. It seems to me that
>> Eq. B.11 is all calculated within the normal virial routine.
>
> The first half of Eq. B11 is the "normal" virial sum, while the  
> second is the fshift contribution.
>
> Normally, when you calculate nonbonded interactions you have to check 
>  the x/y/z periodicity and select the minimum image every time you  
> calculate the distance between the two atoms. This is extremely bad  
> to have in the innermost kernels since modern CPUs just hate  
> conditionals (it breaks the instruction pipeline).
>
> In Gromacs we take care of this during neighborsearching. Instead of  
> having one big list with e.g. 100 neighbors of "particle i", we split 
>  the list and essentially store it like:
>
> (i particle translated one box to the left): 10 neighbors
> (i particle not translated): 30 neighbors
> (i particle translated one box to the right): 0 neighbors
> ... etc ...
>
> The "translations" here are the shift indices, in this case -1,0,1 in 
>  the x dimension.
>
> Now, in the nonbonded kernel the i particle coordinates are loaded in 
>  the _outer_ loop, where we correct it for the current shift. After  
> this we know that the periodicity is correct, so we don't have to  
> include any periodic boundary condition checks in the innermost loop, 
>  which is the time-critical part.
>
> If we combined this with a molecular virial (which would need to be  
> calculated in the innermost loop) all would be fine.
>
> However, when we want to use the atomic virial and apply the single  
> sum of Appendix B1 it turns out that we are making an error, since we 
>  modify the particle coordinates due to the shift in the nonbonded  
> kernel (different for different nlists), but the single sum always  
> use the coordinate in the central box. This error is the second term  
> of Eq. B11.
>
> Fortunately, it turns out to be easy and computationally cheap to  
> correct for. That term can be calculated by simply summing all  
> interaction forces for each shift vector (3*3*3=27), which we can do  
> entirely in the outer nonbonded loop - after we have added the force  
> to the i particle force array we also increment the fshift array  
> corresponding to the shift used for this i particle.
>
> It's perhaps not entirely obvious without a figure, but pretty neat:  
> it completely removes all virial and PBC calculations from them  
> nonbonded innerloop, and we only do PBC checks on neighborsearching.
>
> Cheers,
>
> Erik
>
> PS: Unfortunately we also use the term "molecular shift" in Gromacs,  
> when we create a graph of all bonded interactions in a molecule to  
> move it between  the states where all atoms are in the central box or 
>  where the molecule is whole.  This is really entirely unrelated to  
> the nonbonded shifts, apart from the fact that both deal with  
> periodic boundary conditions. Sorry about the confusing naming ;-)
>
>
>
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-developers-request at gromacs.org.
>
>
>



-- 
\ Jason de Joannis  \
\ Emory University  \
  \ jdejoan at emory.edu \
   \ 404-402-1332      \