Dear Daniel,
I was away and now I am back.
Thanks for your reply.
Since I am working on disaccharide glycolipids, I have 11 acceptors
(O11,O12,O13,O14,O15,O16,O22,O23,O24,O25,O26), and 7 hydroxyl groups
(O12-H12, O13-H13, O16-H16, O22-H22, O23-H23, O24-H24, O26-H26).
My main target was to know, in average, how long each acceptor (O) of a
glycolipid involve in hydrogen bonding interaction no matter which donor is
involved, in the glycolipid bilayer assembly.
So to achieve this, I selected columns only with "O22" (for example:
BMR_42.O22-BMR, BMR_49.O22-BMR, BMR_59.O22-BMR, BMR_23.O22-BMR,
BMR_10.O22-BMR BMR_61.O22-BMR) and sum up all the "1s" and divided with the
total number of frames involved. I repeated this procedure for each
acceptors.
If I am not mistaken, I guess what you have suggested is what I want
actually. But I did not use *readdata* and *lifeltime* commands to get the
final value. I used an awk script to do the selections column wise and
average over all frames.
Welcome any feedback from your side to enhance my understanding.
Thank you.
Vijay
On Sun, Mar 22, 2015 at 12:06 AM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
> Hi,
>
> You still haven't said whether I was right or wrong about you wanting
> to calculate the lifetime of a certain solute atom (or residue) being
> involved in *any* hydrogen bond. Since you haven't said no, I will
> proceed on the assumption that this is what you want. In that case
> this is the procedure I would use.
>
> 1) Calculate the hydrogen bond time series data (you have done this).
>
> 2) Use a script to sum up all columns involving your atom/residue of
> interest. This could be as simple as visually identifying which
> columns contain the data you need (e.g. columns 2, 4, 5, and 7) and
> using something like awk:
>
> awk '{print $2 + $4 + $5 + $7}' solutehb.dat > sum.dat
>
> 3) Read the summed data back into cpptraj with readdata:
>
> readdata sum.dat
>
> 4) Perform lifetime analysis on the summed data. The default settings
> for cutoff etc will still work since anything greater than 0 means a
> hydrogen bond is present:
>
> lifetime sum.dat out life.sum.dat ...
>
> -Dan
>
> On Mon, Mar 16, 2015 at 10:50 AM, Vijay Achari <glycoamber.gmail.com>
> wrote:
> > Dear Dan,
> >
> > Below is the example you gave to explain on how the calculation of HB can
> > done with raw data.
> >
> >
> > HB1-1 HB1-2
> > 1 0
> > 0 1
> > 1 0
> > 1 1
> >
> > Using the method I illustrated before you'd come up with an average of
> > (5 / 4 = 1.25). However, the actual answer of how often is HB1
> > involved in a hydrogen bond is clearly all 4 frames, since when HB1-1
> > is broken, HB1-2 is formed, so its always involved in some kind of
> > hydrogen bond. So to do what you want you will need to write out the
> > raw time series data, sum up the columns corresponding to the hydrogen
> > bonds you are interested in, then run 'lifetime' analysis on that data
> > set. So using the above sets as an example, the set I would actually
> > run lifetime analysis on would look like:
> >
> >
> > Relating to the above example, I understand that we need to add the "1"s
> > and divide by total number of lines.
> >
> > (total of all "1"s) / (total rows with "1"s)
> >
> >
> > So how with the example below? (I modified the above sample).
> >
> > HB1-1 HB1-2
> > 1 0
> > 0 1
> > 1 0
> > 1 1
> > 0 1
> > 0 0
> > 0 0
> >
> > For the above example is the the average is 6/7 or 6/5. The number 7
> > stand for total number of rows and the later stand for rows with the
> > presence of "1" (at least once).
> >
> > Means do I need to count the total rows or the rows with the occurrence
> of
> > "1" only?
> >
> >
> >
> > Your explanation would help me to write the script.
> >
> > Thanks.
> >
> >
> >
> > On Thu, Mar 12, 2015 at 4:53 PM, Vijay Achari <glycoamber.gmail.com>
> wrote:
> >
> >> Dear Dan,
> >>
> >> Following your explanation, I would like to verify few things.
> >>
> >> #1) The file that contains the raw data; is this one "solutehb.dat"
> >>
> >> if yes,
> >>
> >> #2) The format of the data in the file is as below?
> >>
> >> #Frame *BMR_42.O22-BMR_1.O23-H23* BMR_49.O13-BMR_1.O13-H13
> >> BMR_59.O13-BMR_2.O22-H22 BMR_10.O13-BMR_2.O26-H26 BMR_61.O26-BMR_3
> .O22-H22
> >> BMR_23.
> >> O25-BMR_4.O12-H12 BMR_23.O26-BMR_4.O13-H13 BMR_23.O25-BMR_4.O13-H13
> >> BMR_20.O13-BMR_5.O22-H22 BMR_10.O24-BMR_6.O13-H13 BMR_41.O23-BMR_7
> .O23-H2
> >> 3 BMR_41.O24-BMR_7.O23-H23 ....
> >>
> >>
> >> * 1 * 1
> >> 1 1 1
> >> 1 1 1
> >> 1 1 1
> >> 1 1 1
> >> 1 1 1
> >> 1 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1
> >> 1 1 1 . . .
> >>
> >>
> >> #3) So need I operate on this file?
> >>
> >> #4) Is the first data "BMR_42.O22-BMR_1.O23-H23" is corresponding to
> the
> >> value "1" (below #2)
> >>
> >> Thank you.
> >>
> >>
> >>
> >>
> >> On Wed, Mar 11, 2015 at 11:11 PM, Daniel Roe <daniel.r.roe.gmail.com>
> >> wrote:
> >>
> >>> Hi,
> >>>
> >>> Sorry, it's not clear to me what you're trying to calculate. You said
> >>> "I want to know the lifetime of HB for O22 (ACCEPTOR) atom only"; do
> >>> you want the average lifetime of any given O22 involved in a hydrogen
> >>> bond? In that case you will probably need to write your own script for
> >>> it for the time being. I'll try to illustrate why with an example. Say
> >>> I want to ask how often is BMR_49.O22 involved in any hydrogen bond as
> >>> an acceptor. I could naively sum up the total number of frames the
> >>> bond is present (TotFrames) for each instance of this acceptor atom
> >>> and then divide by the total number of lifetimes:
> >>>
> >>> 4 20 5 1.8500 37
> >>> BMR_49.O22-BMR_1.O26-H26
> >>> 7 28 2 1.0357 29
> >>> BMR_49.O22-BMR_1.O13-H13
> >>>
> >>> the total number frames there is a hydrogen bond is 66, and the total
> >>> number of lifetimes is 48, so the average lifetime for BMR_49.O22
> >>> involved in any hydrogen bond (given the data here) is 1.375 frames.
> >>> However, this is only true if the data were sequential, because it
> >>> doesn't take into account the times when one is present and the other
> >>> isn't. For example, given this hydrogen bond data (where HB1
> >>> represents a single acceptor):
> >>>
> >>> HB1-1 HB1-2
> >>> 1 0
> >>> 0 1
> >>> 1 0
> >>> 1 1
> >>>
> >>> Using the method I illustrated before you'd come up with an average of
> >>> (5 / 4 = 1.25). However, the actual answer of how often is HB1
> >>> involved in a hydrogen bond is clearly all 4 frames, since when HB1-1
> >>> is broken, HB1-2 is formed, so its always involved in some kind of
> >>> hydrogen bond. So to do what you want you will need to write out the
> >>> raw time series data, sum up the columns corresponding to the hydrogen
> >>> bonds you are interested in, then run 'lifetime' analysis on that data
> >>> set. So using the above sets as an example, the set I would actually
> >>> run lifetime analysis on would look like:
> >>>
> >>> (HB1-1)+(HB1-2)
> >>> 1
> >>> 1
> >>> 1
> >>> 2
> >>>
> >>> Hope this helps,
> >>>
> >>> -Dan
> >>>
> >>> On Wed, Mar 11, 2015 at 12:57 AM, Vijay Achari <glycoamber.gmail.com>
> >>> wrote:
> >>> > Dear Dan,
> >>> >
> >>> > In my case the time gap in between two frames (taken for analysis) is
> >>> 5ps.
> >>> >
> >>> > #Set Nlifetimes MaxLT AvgLT TotFrames SetName
> >>> > 0 17 215 55.7059 947
> >>> > BMR_42.O22-BMR_1.O23-H23
> >>> > 1 1 1 1.0000 1
> >>> > BMR_42.O23-BMR_1.O23-H23
> >>> > 2 31 3 1.0968 34
> >>> > BMR_31.O22-BMR_1.O26-H26
> >>> > 3 11 2 1.0909 12
> >>> > BMR_31.O13-BMR_1.O26-H26
> >>> > 4 20 5 1.8500 37
> >>> > BMR_49.O22-BMR_1.O26-H26
> >>> > 5 24 4 1.7917 43
> >>> > BMR_49.O23-BMR_1.O26-H26
> >>> > 6 41 2 1.0488 43
> >>> > BMR_49.O12-BMR_1.O12-H12
> >>> > 7 28 2 1.0357 29
> >>> > BMR_49.O22-BMR_1.O13-H13
> >>> > 8 1 1000 1000.0000 1000
> >>> > BMR_49.O13-BMR_1.O13-H13
> >>> > 9 1 1 1.0000 1
> >>> > BMR_22.O12-BMR_2.O22-H22
> >>> > 10 50 62 18.9400 947
> >>> > BMR_59.O13-BMR_2.O22-H22
> >>> > 11 1 1 1.0000 1
> >>> > BMR_59.O25-BMR_2.O23-H23
> >>> > 12 34 5 1.1471 39
> >>> > BMR_59.O13-BMR_2.O23-H23
> >>> > 13 7 6 2.1429 15
> >>> > BMR_59.O26-BMR_2.O24-H24
> >>> > 14 6 8 3.6667 22
> >>> > BMR_59.O25-BMR_2.O24-H24
> >>> > 15 28 2 1.1071 31
> >>> > BMR_10.O26-BMR_2.O26-H26
> >>> > 16 52 121 15.1346 787
> >>> > BMR_10.O25-BMR_2.O26-H26
> >>> > 17 65 9 1.9231 125
> >>> > BMR_10.O13-BMR_2.O26-H26
> >>> > 18 4 4 2.2500 9
> >>> > BMR_59.O26-BMR_2.O26-H26
> >>> > 19 9 1 1.0000 9
> >>> > BMR_22.O15-BMR_2.O13-H13
> >>> > 20 1 1 1.0000 1
> >>> > BMR_10.O11-BMR_2.O16-H16
> >>> > 21 7 2 1.1429 8
> >>> > BMR_53.O13-BMR_2.O16-H16
> >>> > 22 2 3 2.5000 5
> >>> > BMR_33.O13-BMR_3.O22-H22
> >>> > 23 97 54 6.8247 662
> >>> > BMR_61.O26-BMR_3.O22-H22
> >>> > 24 1 1 1.0000 1
> >>> > BMR_29.O26-BMR_3.O23-H23
> >>> > 25 31 36 3.3226 103
> >>> > BMR_29.O16-BMR_3.O23-H23
> >>> > 26 22 2 1.2273 27
> >>> > BMR_30.O22-BMR_3.O23-H23
> >>> > 27 4 1 1.0000 4
> >>> > BMR_33.O22-BMR_3.O23-H23
> >>> > 28 38 2 1.1316 43
> >>> > BMR_33.O13-BMR_3.O23-H23
> >>> > (there are more than 2000 lines, but I work on only 28 lines to get
> some
> >>> > understanding)
> >>> >
> >>> >
> >>> > For clarity, I shall show you how I worked on the given results
> above.
> >>> >
> >>> > I want to know the lifetime of HB for O22 (ACCEPTOR) atom only. So, I
> >>> did
> >>> > in this way,
> >>> >
> >>> > 1) sum the values from column AvgLT for the only occurrences of O22.
> >>> > 2) than find the average of that, where the denominator would be the
> >>> number
> >>> > of occurrences of O22.
> >>> >
> >>> > Based on the above steps, the results are :
> >>> >
> >>> > O22 occur 6 times,
> >>> > sum of O22 is 61.916, and
> >>> > average of O22 is 10.319.
> >>> >
> >>> > So, I figured out the lifetime of O22 would be 10.319 x 5ps = 51.595
> >>> ps.
> >>> >
> >>> > Is this correct? Did I choose the correct column "AvgLT"?
> >>> >
> >>> > I hope to get some feedback if this is correct way to do it.
> >>> >
> >>> > Many thanks in advance.
> >>> >
> >>> > Vijay
> >>> >
> >>> >
> >>> > On Wed, Mar 11, 2015 at 11:50 AM, Daniel Roe <daniel.r.roe.gmail.com
> >
> >>> wrote:
> >>> >
> >>> >> On Tue, Mar 10, 2015 at 9:24 PM, Vijay Achari <glycoamber.gmail.com
> >
> >>> >> wrote:
> >>> >> > Could you explain on how to get the lifetime value in pico second
> >>> (ps)?
> >>> >>
> >>> >> This depends on how often you recorded your coordinate trajectory.
> For
> >>> >> example, say you ran a simulation with a timestep of 2 fs (dt=0.002
> >>> >> ps) and you recorded a trajectory frame every 500 steps (ntwx=500).
> >>> >> This means that each frame in your trajectory has been recorded at 1
> >>> >> ps intervals. However, say you recorded your trajectory every 5000
> >>> >> steps instead - your trajectory then will have been recorded at 10
> ps
> >>> >> intervals. Since lifetimes are always given in frames, it should be
> >>> >> easy to convert to ps based on how often your coordinate trajectory
> >>> >> was written to (e.g. in the latter case a max lifetime of 1 frame
> >>> >> would mean 10 ps).
> >>> >>
> >>> >> Hope this helps,
> >>> >>
> >>> >> -Dan
> >>> >>
> >>> >> >
> >>> >> > Thanks in advance.
> >>> >> > Vijay
> >>> >> >
> >>> >> >
> >>> >> > On Tue, Mar 10, 2015 at 10:15 PM, Daniel Roe <
> daniel.r.roe.gmail.com
> >>> >
> >>> >> wrote:
> >>> >> >
> >>> >> >> Hi,
> >>> >> >>
> >>> >> >> On Tue, Mar 10, 2015 at 2:53 AM, Vijay Achari <
> glycoamber.gmail.com
> >>> >
> >>> >> >> wrote:
> >>> >> >> > generate two files with lifetime information. The
> >>> >> *solute.lifetime.dat*
> >>> >> >> contain
> >>> >> >> > info like:
> >>> >> >> >
> >>> >> >> > #Set Nlifetimes MaxLT AvgLT TotFrames
> SetName
> >>> >> >> > 0 22 1 1.0000
> 22
> >>> >> >> BMR_3.O16-BMR_1.O22-H22
> >>> >> >> > 1 296 346 15.6959 4646
> >>> >> >> BMR_42.O22-BMR_1.O22-H22
> >>> >> >> > 2 992 12 1.4688 1457
> >>> >> >> BMR_42.O14-BMR_1.O22-H22
> >>> >> >> > 3 1 1 1.0000
> 1
> >>> >> >> BMR_42.O13-BMR_1.O22-H22
> >>> >> >> > 4 189 12 1.1429 216
> >>> >> >> BMR_57.O25-BMR_1.O22-H22
> >>> >> >> > 5 462 410 12.3074 5686
> >>> >> >> BMR_57.O16-BMR_1.O22-H22
> >>> >> >> >
> >>> >> >> > I would like to know how I can go from here to calculate the
> >>> >> hb-lifetime
> >>> >> >> between
> >>> >> >> > solute and solute?
> >>> >> >>
> >>> >> >> I'm not really sure I understand your question. The data output
> you
> >>> >> >> posted is exactly the lifetime calculation. For example, the
> second
> >>> >> >> set (1) contains lifetime information for the hydrogen bond
> between
> >>> >> >> residue 42, atom O22 and residue 1, atoms O22-H22; there were 296
> >>> >> >> individual lifetimes (i.e. the hbond formed 296 times), the max
> of
> >>> >> >> which lasted 346 frames, the average lifetime is ~15.7 frames.
> Let
> >>> me
> >>> >> >> know if I'm not understanding you or if I can explain more.
> >>> >> >>
> >>> >> >> -Dan
> >>> >> >>
> >>> >> >> >
> >>> >> >> > I have read the pages 556-557 from AMBER 14 manual, but I find
> it
> >>> >> >> difficult
> >>> >> >> > to see how one should start processing and getting the lifetime
> >>> value.
> >>> >> >> >
> >>> >> >> > I think simple example would help me much in this case.
> >>> >> >> >
> >>> >> >> > Could you give me some example on how this can be obtained?
> >>> >> >> >
> >>> >> >> > Your help is much appreciated.
> >>> >> >> >
> >>> >> >> > Thank you.
> >>> >> >> > Vijay
> >>> >> >> > _______________________________________________
> >>> >> >> > AMBER mailing list
> >>> >> >> > AMBER.ambermd.org
> >>> >> >> > http://lists.ambermd.org/mailman/listinfo/amber
> >>> >> >>
> >>> >> >>
> >>> >> >>
> >>> >> >> --
> >>> >> >> -------------------------
> >>> >> >> Daniel R. Roe, PhD
> >>> >> >> Department of Medicinal Chemistry
> >>> >> >> University of Utah
> >>> >> >> 30 South 2000 East, Room 307
> >>> >> >> Salt Lake City, UT 84112-5820
> >>> >> >> http://home.chpc.utah.edu/~cheatham/
> >>> >> >> (801) 587-9652
> >>> >> >> (801) 585-6208 (Fax)
> >>> >> >>
> >>> >> >> _______________________________________________
> >>> >> >> AMBER mailing list
> >>> >> >> AMBER.ambermd.org
> >>> >> >> http://lists.ambermd.org/mailman/listinfo/amber
> >>> >> >>
> >>> >> > _______________________________________________
> >>> >> > AMBER mailing list
> >>> >> > AMBER.ambermd.org
> >>> >> > http://lists.ambermd.org/mailman/listinfo/amber
> >>> >>
> >>> >>
> >>> >>
> >>> >> --
> >>> >> -------------------------
> >>> >> Daniel R. Roe, PhD
> >>> >> Department of Medicinal Chemistry
> >>> >> University of Utah
> >>> >> 30 South 2000 East, Room 307
> >>> >> Salt Lake City, UT 84112-5820
> >>> >> http://home.chpc.utah.edu/~cheatham/
> >>> >> (801) 587-9652
> >>> >> (801) 585-6208 (Fax)
> >>> >>
> >>> >> _______________________________________________
> >>> >> AMBER mailing list
> >>> >> AMBER.ambermd.org
> >>> >> http://lists.ambermd.org/mailman/listinfo/amber
> >>> >>
> >>> > _______________________________________________
> >>> > AMBER mailing list
> >>> > AMBER.ambermd.org
> >>> > http://lists.ambermd.org/mailman/listinfo/amber
> >>>
> >>>
> >>>
> >>> --
> >>> -------------------------
> >>> Daniel R. Roe, PhD
> >>> Department of Medicinal Chemistry
> >>> University of Utah
> >>> 30 South 2000 East, Room 307
> >>> Salt Lake City, UT 84112-5820
> >>> http://home.chpc.utah.edu/~cheatham/
> >>> (801) 587-9652
> >>> (801) 585-6208 (Fax)
> >>>
> >>> _______________________________________________
> >>> AMBER mailing list
> >>> AMBER.ambermd.org
> >>> http://lists.ambermd.org/mailman/listinfo/amber
> >>>
> >>
> >>
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
>
>
>
> --
> -------------------------
> Daniel R. Roe, PhD
> Department of Medicinal Chemistry
> University of Utah
> 30 South 2000 East, Room 307
> Salt Lake City, UT 84112-5820
> http://home.chpc.utah.edu/~cheatham/
> (801) 587-9652
> (801) 585-6208 (Fax)
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Sat Mar 21 2015 - 22:30:02 PDT
