[gmx-developers] virial

[Erik Lindahl]

Hi,

On Nov 22, 2006, at 12:26 AM, Jason de Joannis wrote:

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