Jason de Joannis
jdejoan at emory.edu
Thu Jan 11 01:11:50 CET 2007
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.
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:
> 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.
>> Quoting Erik Lindahl <lindahl at cbr.su.se>:
>>> Hi Jason,
>>>> 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
>>>> 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.
> 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
> 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 \
More information about the gromacs.org_gmx-developers