[gmx-developers] virial

[Jason de Joannis]

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 ;-)