[gmx-users] Calculating the average separation between two multi-atom groups

[Mark Abraham]

On 31/05/2012 11:24 AM, Andrew DeYoung wrote:
> Hi,
>
> I have a system in a slab geometry.  A surface exists at z = 0; many
> hydrogens protrude from the surface, and these hydrogens are mostly (but not
> precisely) immobile.  Above the surface, there is liquid, including the
> anion BF4- (tetrahedral arrangement of fluorines around a central boron).
> The liquid region is large, extending far from the surface, (far above the
> surface, the liquid behaves like liquid in the bulk).
>
> There is hydrogen-like bonding between the hydrogens protruding from the
> surface and the fluorines in BF4-.  I would like to calculate the "average"
> H-F distance, where H atoms protrude from the stationary surface and F atoms
> exist on the BF4- ions.  But saying that I want to compute the "average" H-F
> distance is very vague.  I can think of at least two possible, hopefully
> reasonable, ways to formulate the problem:
>
> (i) For the purpose of calculating the average H-F separation, only consider
> fluorines on BF4- ions which are within a certain perpendicular distance z0
> from the surface.  In other words, consider the BF4- ions which lie in the
> region 0<  z<  z0 (where z0 is positive and very small compared to the z
> dimension of the simulation box).  Then, using those BF4- ions, I calculate
> the (time-averaged) H-F separation.
>
> (ii) For the purpose of calculating the average H-F separation, only
> consider fluorines when they are a certain small distance from any hydrogen.
>
> Are (i) or (ii) these feasible?
>
> For (i), I can think about using g_select to select BF4- ions which are a
> distance of z0 or less from the surface at z = 0.  Maybe I would use a
> selection like 'res_com of resname BF4 and z<  10' (where z0 = 10).  The
> problem with this is that, I think, I would obtain an index file for each
> simulation timestep.

Yes, you'd need to do some iteration over g_select and then another 
tool. GROMACS 5 will likely do that for you, but nirvana is some way off...

>    So, I guess then if I have 200,000 simulation
> timesteps, I would have to run g_bond 200,000 times!  (Or would g_dist be
> appropriate here?)  Also, even my formulation in (i) is a little awkward;
> fluorines at one edge of the xy dimension would be far from hydrogens
> immobilized at the other side of the xy dimension, so I would get artifacts.

These tools *should* work correctly with PBC, but I would advise 
checking that they do, e.g. by making some very small subsets that you 
know cross boundaries and checking their stats agree with subsets that 
you know do not cross boundaries.

> For (ii), it seems that g_hbond might be useful.  However, it does not seem
> that fluorine is currently implemented as a hydrogen bond acceptor for use
> in g_hbond.  I have never attempted to modify the Gromacs code and I am not
> sure how easy this would be.  But if H-F is a hydrogen-like bond, then
> (average) H-F bond length is what I am going after, I guess.
>
> Do you know of any recipes with which to do this, or do you have any
> suggestions?  Thanks so much!

You can match the trajectory to some .tpr that didn't create it. So if 
you generate a new topology that uses PO4 instead of BF4, while 
preserving the atom ordering over the whole system, you can fool g_hbond 
with a new .tpr without needing to touch any code. You don't even need 
sensible parameters for PO4, of course.

Alternatively, adding

((*top->atoms.atomname[n])[0] == 'F') ||

to the conditional that resembles

             if ((bContact ||
                  (((*top->atoms.atomname[n])[0] == 'O') ||
                   (bNitAcc && ((*top->atoms.atomname[n])[0] == 'N')))) &&
                 ISINGRP(datable[n])) {
                 datable[n] |= ACC; /* set the atom's acceptor flag in 
datable. */
                 add_acc(a,n,grp);
             }

in search_acceptors() of src/gmxlib/gmx_hbond.c and recompiling is 
faster still - assuming you have no other atoms whose names start with 
"F". For safety, name the executable something that will remind you that 
it is not the standard one :-)

Mark