AMBER Archive (2009)

Subject: Re: [AMBER] Amber: water-water hydrogen bonding analysis on a subset of waters in a trajectory

From: Thomas Cheatham III (
Date: Wed Mar 04 2009 - 13:59:30 CST

> My current simulation contains a single mononucleotide, 5 Na+ and 5
> Cl- ions, 20 urea molecules and 1093 TIP3P waters. My goal is to
> perform a hydrogen bond analysis between water molecules within 6
> angstroms of the mononucleotide AND beyond 6 angstroms of the
> mononucleotide to characterize water-water hydrogen bonds in a "local"
> domain and a "bulk" domain.

This sounds tricky, but I think you are on the right track assuming the
mask selections work as advertized. Verify this by creating pdbs and
looking at them.

> Can ptraj perform such an analysis? So far, I have thought of two
> ways to attack this problem.

Not automatically.

> 1) Using ptraj, I have tried hbond. The ptraj input script I have
> used to analyze just donor water-water hydrogen bonds within 6
> angstroms of my mononucleotide (residues :1-2; :1 represents the
> phosphate, :2 the sugar and base) is:
> trajin AMP2md.mdcrd.gz
> center :1-2 mass origin
> image origin center familiar
> rms first mass :1-2
> #
> # water donor
> #
> donor mask "(:1-2 < @6) & :WAT_at_O"
> hbond distance 3.5 angle 120 solvent nieghbor 6 solventacceptor WAT O
> H1 solventacceptor WAT O H2 time 1.0 series hbww out
> hbond_wat_watdonor.out

Neighbor is misspelled, and the command is solventneighbor; also you need
continuation lines on the above:

 hbond distance 3.5 angle 120 solventneighbor 6 solventacceptor WAT O \
  H1 solventacceptor WAT O H2 time 1.0 series hbww out \

> This script generates a segmentation fault. I have tried 3.5 angstroms

Yes, not a very useful piece of info.

> as the distance from the mononucleotide to test with fewer waters, but I
> still get the segmentation error. I suppose I could be running out of
> memory even using 3.5 angstroms. My script might also have an error.
> Regardless, I guess that the donor mask does not reselect waters every
> frame. Is this correct? In that case, I suppose I could write a script
> to determine the hydrogen bonds for each frame.

Ah yes, that is true. To get around this I would use strip and/or closest

  strip (:1-2 < @6) & :WAT_at_O


  strip (:1-2 > @6) & :WAT_at_O

Make sure you have imaged the water and placed :1-2 at the center prior

  center :1-2 mass origin
  image origin center familiar

However, strip will still be broken as it does not check every step.
(i.e. the same waters would be stripped). I will have to think about how
to approach this.

So, in the meantime I would do the closest or do it frame by frame as you
suggested in a Perl script or some such via strip. It may be possible to
make a furthest command simply by changing a few less-thans; I will look
at this too...


[note also there have been reports that the selection mechanism with < or
> is hinky and I'll look into that too. ...and note that as mentioned in
the reflector closest is broken in AMBER 10 but I have the fix and will
soon update the WWW sites.]

AMBER mailing list