# [gmx-developers] virial

Jason de Joannis jdejoan at emory.edu
Mon Jan 8 22:32:20 CET 2007

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

```