AMBER Archive (2009)

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

From: Jeff Schwinefus (schwinef_at_stolaf.edu)
Date: Wed Mar 04 2009 - 23:22:34 CST

Thomas

Thank you for your input.

On Wed, Mar 4, 2009 at 1:59 PM, Thomas Cheatham III <tec3_at_utah.edu> wrote:
>
>> 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 \
>  hbond_wat_watdonor.out

Sorry about the sloppiness of my first script. I did find those
errors. My script for a single frame in the trajectory now is:

trajin AMP2md.mdcrd.gz 1 1

center :1-2 mass origin
image origin center familiar
rms first mass :1-2

#
# water donor
#
donor mask "(:1-2 <@3.5) & :WAT_at_O"

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

This script either generates a segmentation fault or I receive the
message (when it gets to donor mask):

PTRAJ: donor mask "(:1-2 < @3.5) & :WAT_at_O"
*** glibc detected *** ptraj: free(): invalid pointer: 0x00007fff2cf5eda0 ***

Unless I have an error in my donor mask, there might be something
going on with the <@ operator.

>
>> 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
> waters.
>
>  strip (:1-2 < @6) & :WAT_at_O
>
> or
>
>  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

Some might have spotted this earlier than me, but there is a problem
using closest. I specified the donor mask as water oxygens and
acceptor mask as water hydrogens and hbond seemed to work fine after
closest (selecting 32 nearest WAT based on a previous watershell test
run with a distance of 3.5 angstroms) over 500 frames. However, after
summing water hydrogen bonds for a specific oxygen atom, it seemed
that the %occupancy was too small (~60%). Even if everything is
working correctly, I realize now that by using closest, the waters on
the outermost solvent "shell" furthest from the mononucleotide do not
have hydrogen bonding partners because they were removed with closest.

At this point, I am going to have to spend some time thinking about this.

Thanks again.

Jeff

> at this too...
>
> --tec3
>
> [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
> AMBER_at_ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
> 

-- 
Jeff Schwinefus
Associate Professor of Chemistry
St. Olaf College
Northfield, MN 55057
(507) 786-3105
FAX (507) 786-3968

_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber