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

[Andrew DeYoung]

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.  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.


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!

Andrew DeYoung
Carnegie Mellon University