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
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
