[gmx-developers] Bond-angle coupling

[Bogdan Costescu]

On Thu, Oct 4, 2012 at 9:10 PM, MURAT OZTURK <murozturk at ku.edu.tr> wrote:
> ... Comments
> being sparse as they are (not complaining, just observing) I thought I would
> try and steal some more of your time.

Indeed, comments are quite sparse. But where they are, they are good :)

> I have simplified my hack even further. I will hard-code a threshold into
> tab_angles() function so that it will be turned off when the distance
> between atoms ai and aj is above the threshold.

OK, sounds simple enough.
But what about the distance aj-ak ? That's another bond... You should
not make assumptions about the order of the atoms, so very likely
you'll have to test for both these distances.

> I am trying to figure out how the distance is calculated in bonds(), but it
> is surprisingly fancy :
>
> ki = pbc_rvec_sub(pbc,x[ai],x[aj],dx);
> dr2 = iprod(dx,dx);
> dr = dr2*gmx_invsqrt(dr2);
>
> I have trouble understanding what is what here. pbc_rvec_sub takes care of
> pbc, but how?

You should follow the code. Usually the .h files which carry the
functions' prototypes have a comment describing the functions,
otherwise look up the actual implementation.

> dx is defined, but never initialized before, so what is it?

dx is the result, it's initialized as result of this call. It's the
vector distance.

> I can't find any comment as to what iprod() or gmx_invsqrt() does.

Look them up, vec.h might be the place :)
iprod() does inner product of 2 vectors such that dr2 is the square of
the inter-atom distance, gmx_invsqrt() does an optimized 1/sqrt().
You'll see this kind of sequence quite often in the code, it's the
fastest way to calculate the distance.

> But never mind the details, I assume in the end dr is the distance between
> atoms ai and aj. Is this correct?

Yes, it's correct.

> Then we skip to tab_angles() :
>
> *dvdlambda += bonded_tab(blabla , &va, &dVdt);
>
> So we get the potential and force from the table and put them into va and
> dVdT.
>
> What if I calculate the distance just before this line, and put this line
> within an if clause. Never call bonded_tab() if distance is above threshold,
> and assign va=0 and dvDt=0. Is it that simple?

It could be, if that's what you want. But you should decide also
whether it makes sense at all to execute the rest of the loop...
meaning that you can put the if around the whole code between the
bonded_tab() call and the end of the loop.

I would also suggest to do the distance calculation slightly more
efficient, by using a modified version of bond_angle(). It already
calculates the vector distances which are needed again for the
distance calculations... Leave the original function in its place as
it's called from other functions and create another one which
calculates both the angle and the distances and call it only from
tab_angles().

> And lastly, as a bonus question, what does dvdlambda keep? Do I even want to
> know?

Do you do free energy calculations ? If not, then you don't need to know :)
If not doing free energy calculations, forceparams[type].tab.kA is
equal to forceparams[type].tab.kB so dvdlambda will be zero.

Again, good luck!
Bogdan