> I want to know how many water molecules are around a given amide
> hydrogen at a given distance. And do this throughout a trajectory for
> each residue??
The hydrogen bond facility of ptraj can be useful for this purpose,
however it doesn't directly report the number of atoms within a given
distance but the average occupancy (which relates to the number of atoms).
The occupancy is the percent over the whole trajectory that the distance and
angle is satisfied. The command does not (currently) give the frame by
frame data, however it could be modified for this purpose (i.e. if you
wanted to know the instantaneous occupancy).
The documentation is a bit sparse so I will elaborate a bit (and hopefully
answer earlier e-mails on this question to the reflector at the same
time).
Basically, this utility is a way to keep track of a list of interactions,
defined as electron pair donors (a single atom, i.e. O) and electron pair
acceptors (two atoms, i.e. N-H) subject to a distance and angle cutoff.
The distance is the heavy atom distance and the angle is the angle off the
plane (so zero represents a linear bond)
N---H-------
\ angle
\
O
The code however is general in that the acceptor can be a single atom
(i.e. if you list the same atom twice in the specification, i.e. Na+ Na+)
to allow for tracking atom-atom interactions (i.e. ion binding to donors).
To use it, first you set up a list of possible "donor" and
"acceptor" groups with the "donor" and "acceptor" commands in ptraj.
| --> means OR
<filename> --> means input the name of the file
filename --> means a command to ptraj called filename
donor usage:
donor clear
...clear the current list of hbond donors
donor print
...print the current list of hbond donors
donor <resname> <atomname> | mask <mask>
i.e. donor ADE N7
...means set all atoms named N7 in residue ADE as potential donors.
donor mask :10.N7
...means set the atom named N7 in residue 10 as a potential donor.
acceptor usage:
acceptor clear
acceptor print
acceptor <resname> <atomname> <atomnameH> | mask <mask> <maskH>
i.e. acceptor ADE N6 H61
...means set the acceptor to be for residues ADE heavy atom N6 and
light atom H61. Note that <atomname> and <atomnameH> can be the same
atom (for example to look at donor -- ion interactions).
The above only sets up the list of donor and acceptors. To actually do
something, use the hbond command. All options are optional.
distance <value> ...set the cutoff distance
angle <value> ...set the angle cutoff
solventneighbor <value> ...how many "solvents" to keep track of
solventdonor <donor-spec> ...what is a solvent donor, syntax as above
hbond solventacceptor <acceptor-spec> ...what is a solvent acceptor
nosort ...don't sort by occupancy
time <value> ...what is the time between frames
print <value> ...print the results
series <name> ...create a time series (calc. occ/lifetimes)
The solvent stuff needs a little more explaining. The idea is that not
only might we want to keep track of a "static" list of possible h-bond
interactions, we may want to keep track of a dynamic list. An example is
with water since for example, we don't care which water is in contact with
a donor in the minor groove, just that water is there. Since a given
donor/acceptor may be able to interact with more than 1 water (for example
a Na+ ion which can interact with 6 in the first solvation shell), we have
to specify how many (within our cutoff) we may contact. This is
"solventneighbor". For the lifetimes/occupancy, exchanges of the water
are properly factored in...
So, a general input for DNA:
trajin ../charmm/t4cga4_rhdo_cornell_1.crd 1 10 1
trajout ctraj.strip PDB
center :1-20 mass origin
image origin center
prnlev 2
donor mask :1-20.O3',O4',O5',O1P,O2P
donor mask :ADE.N1,N3,N7
donor mask :THY.O2,O4
donor mask :GUA.O6,N3
donor mask :CYT.N1,N3,O2
acceptor mask :1.O5' :1.H5T
acceptor mask :11.O5' :11.H5T
acceptor mask :10.O3' :10.H3T
acceptor mask :20.O3' :20.H3T
acceptor mask :ADE.N6 :ADE.H61
acceptor mask :ADE.N6 :ADE.H62
acceptor mask :THY.N3 :THY.H3
acceptor mask :GUA.N2 :GUA.H21
acceptor mask :GUA.N2 :GUA.H22
acceptor mask :GUA.N1 :GUA.H1
acceptor mask :CYT.N4 :CYT.H41
acceptor mask :CYT.N4 :CYT.H42
prnlev 0
hbond series h1 time 1.0 distance 3.0 solventdonor WAT OH2 \
solventacceptor WAT OH2 H1 solventacceptor WAT OH2 H2
go
Note that you need to store
(donor+solventneighbor)x(acceptor+solventneighbor)
potential interactions so this can get memory intensive. Therefore, you
may want to limit (and break it up into multiple runs) the donor/acceptor
lists.
The output lists the atoms involved in the interaction, the %occupied
value (i.e. how often over the trajectory the interaction was satisfied),
the average distance (and std dev) when the distances was satisfied and
the average angle off the plane (and std dev). If the "series" keyword is
specified, in addition the lifetime and maxocc is printed along with a
time series schematic of the interaction over the trajectory. The
lifetime in units of time (by default 1 ps between frames; this can be
changed with the time keyword) is the average time the interaction was
satisified. The maxocc is the maximum number of frames in a row the
interaction was satisfied; it is NOT in units of time.
To clarify the meaning of the lifetime, basically if the "bond" is
satisfied, you increment a counter. If it breaks, you start the count
again. The sum of the counters / the number of occurences equals the
lifetimes. This description is a little bit simplified for the water case
since there we allow "gaps" without breaking the lifetime, i.e. if the
bond breaks for only for 1 frame, we do not count it as a break...
Example: x's below means the bond is formed, spaces mean it isn't
xxxxx xxxx xxx xx X
The lifetime would be 5+4+3+2+1 / 5 = 3
The () mean standard deviations. If the lifetime = the length of the
trajectory, the stdev would be zero.
Good luck and Happy Holidays,
\ Thomas E. Cheatham, III (Assistant Professor) College of Pharmacy, Depts of
| Medicinal Chemistry and of Pharmaceutics and Pharmaceutical Chemistry
| Adjunct Asst Prof of Bioengineering; Center for High Performance Computing
| University of Utah, 30 South 2000 East, Skaggs 201, Salt Lake City, UT 84112
|
| tec3.utah.edu (801) 587-9652; FAX: (801) 585-9119
\ BPRP295A / INSCC 418
http://www.chpc.utah.edu/~cheatham
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Jan 14 2004 - 15:53:09 PST
