AMBER Archive (2002)

Subject: Re: ion distances

From: Thomas E. Cheatham, III (cheatham_at_chpc.utah.edu)
Date: Thu Nov 21 2002 - 17:54:46 CST

[p.s. if people have better solutions, I'd love to hear them...]

> how can I get the distances from the ions to the DNA,
> with either carnal or ptraj, is there a way other
> than typing explicitly all the atoms?

If you want all the distances, there is no simple way in ptraj or carnal
(to my knowledge) to extract this information without listing all possible
distances (which can be done with the distance command in ptraj to
generate files for each distance and then a perl script or some such to
process these files).

If however you are simply interested in finding out if ions are "binding"
to your DNA molecule (i.e. you do not care about long distance
interactions), this can be done with the hydrogen bond capability
in ptraj quite simply. You can do this either for each and every ion
individually, or likely better, assume that the ions are indistinguishable
and therefore as long as any ion is close, keep track of it. This is done
with reference to a particular atom (or set of atoms) in the DNA (i.e. the
program monitors interactions between sets of donor and sets of acceptors
subject to a distance/angle criteria).

However note that ion binding to DNA in MD simulation appears to be rather
sensitive to the initial conditions, i.e. if an ion is initially bound, it
may stay there for quite some time! For this reason, I advocate
randomizing the initial ion positions (by swapping with water) such that
ions are not initially close to the DNA (or another ion). This can be
done with the randomizeion command in newer versions of ptraj as has been
discussed in earlier postings to the reflector.

Using the hbond commands in ptraj:

To use it, first you set up a list of possible "donor" and "acceptor"
groups with the "donor" and "acceptor" commands. In my terminology,
currently a donor means an electron pair donor (O, N, ...) and an acceptor
means an electron pair acceptor (N-H, O-H). However, if you simple want
to understand contact pairs, you can define donors and acceptors to be any
atom you want. If the acceptor atom specification is repeated twice, no
angle is calculated, i.e.

  donor GLY O
  acceptor SER O O

will keep track of all GLY-O --- O-SER interactions.

In the specification below, | means OR.

SPECIFY DONORS, "donor" command 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 select all atoms named N7 in residue ADE as potential donors.

      donor mask :10_at_N7

   ...means select the atom named N7 in residue 10 as a potential donor.

SPECIFY ACCEPTORS:

   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!

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)
        out <name> ...output the results to a file

So, to keep track of all interactions within 3.0 angstroms with a 120
degree angle cutoff,

  hbond distance 3.5 angle 120.0 series out hbondinfo.dat

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 (of
indistinguishable interactions). An example is with water; we often do
not 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". Note that in this
manner the occupancies can be > 1 (100%).

So, a input for DNA might be:

   trajin ../charmm/t4cga4_rhdo_cornell_1.crd 1 10 1
   trajout ctraj.strip PDB
   center :1-20 mass origin
   image origin center familiar

   prnlev 2

   donor mask :1-20_at_O3',O4',O5',O1P,O2P
   donor mask :ADE_at_N1,N3,N7
   donor mask :THY_at_O2,O4
   donor mask :GUA_at_O6,N3
   donor mask :CYT_at_N1,N3,O2

   acceptor mask :1_at_O5' :1_at_H5T
   acceptor mask :11_at_O5' :11_at_H5T
   acceptor mask :10_at_O3' :10_at_H3T
   acceptor mask :20_at_O3' :20_at_H3T
   acceptor mask :ADE_at_N6 :ADE_at_H61
   acceptor mask :ADE_at_N6 :ADE_at_H62
   acceptor mask :THY_at_N3 :THY_at_H3
   acceptor mask :GUA_at_N2 :GUA_at_H21
   acceptor mask :GUA_at_N2 :GUA_at_H22
   acceptor mask :GUA_at_N1 :GUA_at_H1
   acceptor mask :CYT_at_N4 :CYT_at_H41
   acceptor mask :CYT_at_N4 :CYT_at_H42

   prnlev 0

   hbond series h1 time 1.0 distance 3.0 out hbondinfo.dat \
      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.

Similarly, if you were interested in tracking ion interactions, you could
either list the ions as donor/acceptors or list as
solventdonors/solventacceptors (better); for an Na+ ion (in a residue
named Na+) this would look like:

   hbond series h1 time 1.0 distance 3.0 out hbond_ions.dat \
      solventdonor Na+ Na+ \
      solventacceptor Na+ Na+ Na+

CAVEAT: there are subtleties with the "occupancy" related to how one
defines the distance/angle cutoff, i.e. lifetimes and occupancy critically
depend on your definition of the interaction length.

Here is the partial output from a long DNA trajectory of mine where the
"solventacceptor" (and donors) is Na+ ions.

  HBOND SUMMARY:

  Dumping schematic of time series after each h-bond, key follows:
   | . - o x * @ |
      0-5% 5-20% 20-40% 40-60% 60-80% 80-95% 95-100% occupancy

        DONOR ACCEPTORH ACCEPTOR
  atom# :res_at_atom atom# :res_at_atom atom# :res_at_atom %occupied distance angle lifetime maxocc
| 1234 :39_at_O2 | solvent acceptor | 28.55 2.423 ( 0.13) 180.00 ( 0.00) 2727.5 (2589.) 7054 |@@* |
| 215 :7_at_O2 | solvent acceptor | 28.39 2.420 ( 0.12) 180.00 ( 0.00) 3615.7 (5019.) 10713 |@@* |
| 1128 :36_at_N7 | solvent acceptor | 11.41 2.598 ( 0.20) 180.00 ( 0.00) 40.0 ( 71.6) 529 |-.... . -|
| 1160 :37_at_N7 | solvent acceptor | 10.10 2.601 ( 0.21) 180.00 ( 0.00) 31.1 ( 53.9) 272 | ........|
| 979 :31_at_O2 | solvent acceptor | 9.16 2.469 ( 0.16) 180.00 ( 0.00) 194.5 (299.9) 1002 | -..-. |
| 809 :26_at_N7 | solvent acceptor | 7.19 2.595 ( 0.20) 180.00 ( 0.00) 37.1 ( 83.5) 637 | . .. -|
| 460 :15_at_N7 | solvent acceptor | 6.76 2.622 ( 0.22) 180.00 ( 0.00) 23.9 ( 54.1) 293 | ... . -|
| 469 :15_at_N3 | solvent acceptor | 6.76 2.671 ( 0.20) 180.00 ( 0.00) 80.7 (118.1) 523 | - -. |
| 1192 :38_at_N7 | solvent acceptor | 6.43 2.590 ( 0.18) 180.00 ( 0.00) 54.6 (108.6) 621 | .-. |
| 492 :16_at_N7 | solvent acceptor | 6.33 2.592 ( 0.18) 180.00 ( 0.00) 71.1 (102.6) 542 | . .. -|

Please let me know if you have any problems with the code or in your
understanding. Sometime I will attempt to update the DNA tutorial to show
an example of this in common usage.

Good luck,

\ Thomas E. Cheatham, III (Assistant Professor) College of Pharmacy
| Departments of Medicinal Chemistry and of University of Utah
| Pharmaceutics and Pharmaceutical Chemistry 30 South 2000 East, Room 201
| & Center for High Performance Computing Salt Lake City, Utah 84112
|
| e-mail: tec3_at_utah.edu phone: (801) 587-9652 FAX: (801) 585-9119
\ http://www.chpc.utah.edu/~cheatham Offices: BPRP295A / INSCC 418