Page 1 of 2 1 2 >
Topic Options
#7271 - 06/24/05 11:56 AM lipid (or other alkane) chain torsion histograms
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
The input script in the next post extracts torsion time series for 2 torsions of both chains of a lipid; for a simulation with 72 lipids, this is 72*4 time series, which are defined (via ENTER) in a loop, the trajectory is processed, and each time series is converted to a normalized histogram in another loop. The script expects 3 arguments: N (number of .trj files), and T1 and T2, which are the higher-numbered C atom of the central bond of the torsion, i.e.

charmm N:25 T1:3 T2:4 < chnhst.inp >& chnhst.out

will compute histograms from 25 .trj files for all instances of

C21-C22-C23-C24
C31-C32-C33-C34

and

C22-C23-C24-C25
C32-C33-C34-C35
_________________________
Rick Venable
computational chemist


Top
#7272 - 06/24/05 12:10 PM Re: lipid (or other alkane) chain torsion histograms [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W

* lipid; extract chain torsion histograms @T1 @T2 from @N files
*

! READ RTF, PARAM, PSF, COOR files
stream rtfprm.str
stream psfcrd.str

! NO. OF LIPIDS
set nlpd 72

! COMPUTE REMAINING C ATOM INDICES
calc a1 @T1 - 2
calc a2 @T1 - 1
set a3 @T1
calc a4 @T1 + 1
calc c1 @T2 - 2
calc c2 @T2 - 1
set c3 @T2
calc c4 @T2 + 1

calc mxt 100 * @N
calc mxs 10 + ( 4 * @NLPD )
correl maxt @MXT maxs @MXS maxa ?NATOM noupdate
! FOUR TIME SERIES FOR FINAL RESULT (AVERAGED HISTOGRAMS)
enter hb zero
edit hb total 360 delta 1. offset -180. skip 1
enter hg dupl hb
enter hd dupl hb
enter hq dupl hb

! BETA CHAIN ATTACHED TO GLYCEROL C2, CHAIN C ATOMS START WITH C2
set r 1
label elp
enter b@r dihe L @R C2@A1 L @R C2@A2 L @R C2@A3 L @R C2@A4
enter g@r dihe L @R C3@A1 L @R C3@A2 L @R C3@A3 L @R C3@A4
enter d@r dihe L @R C2@C1 L @R C2@C2 L @R C2@C3 L @R C2@C4
enter q@r dihe L @R C3@C1 L @R C3@C2 L @R C3@C3 L @R C3@C4
incr r by 1
if r le @NLPD goto elp

! OPEN AND PROCESS .trj FILES; START FROM UNIT 8, UP TO 90 FILES
set k 1
label olp
calc u @K + 7
open unit @U file read name dyn@K.trj
incr k by 1
if k le @N goto olp
traj first 8 nunit @N

! CONVERT TO HISTOGRAMS, ACCUMLATE
set k 1
label klp
mantime g@k hist -180. 180. 360
mantime hg add g@k
mantime b@k hist -180. 180. 360
mantime hb add b@k
mantime q@k hist -180. 180. 360
mantime hq add q@k
mantime d@k hist -180. 180. 360
mantime hd add d@k
incr k by 1
if k .le. @NLPD goto klp
! NORMALIZE FOR NUMBER OF LIPIDS
mantime hb divi @NLPD
mantime hg divi @NLPD
mantime hd divi @NLPD
mantime hq divi @NLPD

! FIVE COLUMNS; TIME T1-beta T1-gamma T2-beta T2-gamma
edit hb total 361 delta 1. offset -180. skip 1 veccod 4
open unit 2 write card name tor/chst@T1@T2.dat
write hb unit 2 dumb time
* dumb
*

end
stop

Top
#7273 - 06/26/05 10:53 AM Re: lipid (or other alkane) chain torsion histograms [Re: rmv]
mzzt Offline
Forum Member

Registered: 01/03/05
Posts: 86
Dear Prof. Rick,

" edit hb total 360 delta 1. offset -180. skip 1 "

Is it correct to say that the advantage of redefining the time series with total as 360 etc is because the script deals with the dihedral angle so a range upto 360 degrees is considered?

Also for chst34.dat my output is:

-180.000000 0.016250 0.021667 0.021667 0.021354
-179.000000 0.017292 0.021354 0.020833 0.019792
-178.000000 0.015313 0.021042 0.022188 0.021563
-177.000000 0.014479 0.020208 0.022083 0.021979
-176.000000 0.013958 0.021458 0.023333 0.019896
-175.000000 0.014271 0.019479 0.021250 0.019792
-174.000000 0.011250 0.019167 0.019167 0.020729
-173.000000 0.012083 0.017813 0.017188 0.017292
-172.000000 0.012813 0.018646 0.021146 0.018125
-171.000000 0.011354 0.015000 0.018542 0.018542
-170.000000 0.009896 0.016771 0.017708 0.016563
-169.000000 0.010521 0.010625 0.016771 0.016042
-168.000000 0.008333 0.013438 0.015521 0.013854
-167.000000 0.008333 0.010938 0.014479 0.011667
-166.000000 0.009688 0.012708 0.011875 0.010521
-165.000000 0.006458 0.007813 0.010521 0.010313
-164.000000 0.005208 0.010833 0.008854 0.008438
-163.000000 0.006458 0.007188 0.010000 0.009375
-162.000000 0.004271 0.006354 0.007813 0.007708
-161.000000 0.003750 0.005313 0.008021 0.007813
-160.000000 0.003750 0.005000 0.008021 0.006146
-159.000000 0.003750 0.003542 0.005938 0.006042
-158.000000 0.002396 0.005000 0.006563 0.005729
-157.000000 0.002292 0.003333 0.005521 0.004792
.......
......
.......

Does this give the population of the torsion angle i.e Occurrence vs the torsion angle?

The output is not clear to me, specifically the use of the
mantime g@k hist -180. 180. 360.

Eagerly waiting for your reply,
Regards,
Priyanka.

Top
#7274 - 06/26/05 01:08 PM Re: lipid (or other alkane) chain torsion histograms [Re: mzzt]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
! FIVE COLUMNS; TIME T1-beta T1-gamma T2-beta T2-gamma
edit hb total 360 delta 1. offset -180. skip 1 veccod 4


The posted script has a minor error, as indicated above; the total in the EDIT command should be 360 (not 361) to be consistent with the ENTER declarations and the use of MANTime name HIST. The main purpose of the EDIT command is to allow writing a single file instead of 4 separate files, via VECCOD 4; without that, one could write each of the four averaged histograms (hb hg hd hq) to its own file.


From correl.doc description of MANTime operations:

Code:
 
HIST min max nbins
Q(ibin) = Fraction of Q(t) values within ibin
This command replaces a time series with a
histogram of the time series divided into "nbins" with
a range from "min" to "max". The histogram values sum to 1.



The result in this case is the fractional occupancy of each 1 deg bin, averaged over all of the coordinate sets in the trajectories processed. The HIST processing converts each time series to a histogram of 360 bins from -180. to 180. (1 deg wide); the script sums these over all lipids (into hb, etc.), and then divides by the number of lipids. Since the data sums to 1.0, it is the probability of finding a conformation within a 1 deg torsion bin from the simulation. The populations of the g-, trans and g+ wells can be obtained by integration; I've found it convenient to transform to 0-360 first, because there's usually not much (if any) population in the fully eclipsed state.
_________________________
Rick Venable
computational chemist


Top
#7275 - 08/04/05 07:09 AM Re: lipid (or other alkane) chain torsion histograms [Re: rmv]
mzzt Offline
Forum Member

Registered: 01/03/05
Posts: 86
Dear Prof. Rick,

I am unable to understand one thing in chain torsion histograms script.
Why is it that offset = -180 calculates the fractional occupancy for the range -180 to +180 whereas when offset = 0 gives the values from 0 to +180 only and not upto 360?

Also, I am doing the analysis of a neat bilayer system for which I want to calculate the trans/gauche percentage.

The way I have done this using the torsion histograms script is:
1. calculate the torsions for -180 to +180 for each acyl chain C keeping offset = -180

2. calculate the torsions for -120 to +120 by making offset = -120.

3. Summing the values obtained from the first two points above and then evaluating the fraction.

Is this correct method of calculating trans/gauche percentage or is it that gauche means +60 or -60 strictly and nothing above and nothing less?

Eagerly waiting for your reply,
warm regards,
Priyanka

Top
#7276 - 08/04/05 01:03 PM Re: lipid (or other alkane) chain torsion histograms [Re: mzzt]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
Why is it that offset = -180 calculates the fractional occupancy for the range -180 to +180 whereas when offset = 0 gives the values from 0 to +180 only and not upto 360?

The data limits for "MANTIM series-name HIST" must match what's in the actual time series; changing them won't transform the data. Dihedral time series are computed over the range -180 to 180, so those are the data limits that must be used. If I want data from 0 to 360, I use an external program to transform the histogram data file (by adding 360 to negative torsion values).

1. calculate the torsions for -180 to +180 for each acyl chain C keeping offset = -180
2. calculate the torsions for -120 to +120 by making offset = -120
3. Summing the values obtained from the first two points above and then evaluating the fraction.

I can't quite follow this, but I'm not sure it's right. I usually use the integral of the complete histogram to graphically evaluate the populations for a given torsion. The histograms have been normalized, so that they represent p(phi), and should sum to 1, so that the integral from -120 to 120 would be the fraction gauche.
_________________________
Rick Venable
computational chemist


Top
#7277 - 03/09/06 11:43 AM Re: lipid (or other alkane) chain torsion histograms [Re: rmv]
beginner Offline
Forum Member

Registered: 07/19/04
Posts: 165
Dear Prof Venable,

I did analyze the conformation of alkyl chain in gel state lipid bilayer, so I expected the reporting of trans conformation.

Please find the attachment. The script went through without warning/errors, but I got very low amount of both gauche and trans conformation.

Have I did something nonsense?

Thank you very much for your suggestion!!


Attachments
9568-chst78.txt (480 downloads)


Top
#7278 - 03/09/06 11:46 AM Re: lipid (or other alkane) chain torsion histograms [Re: beginner]
beginner Offline
Forum Member

Registered: 07/19/04
Posts: 165
My output!!


Attachments
9569-torsion78.txt (501 downloads)


Top
#7279 - 03/09/06 01:45 PM Re: lipid (or other alkane) chain torsion histograms [Re: beginner]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
The output indicates almost all trans for the initial time series averages (ca. 180 and -180). If you want the exact populations, you'd need to integrate the histograms produced in the output .hst file. The CHARMM output does not report the populations directly. One could add additional MANTINE commands to integrate the 4 histograms, or do it externally.
_________________________
Rick Venable
computational chemist


Top
#7280 - 03/09/06 06:14 PM Re: lipid (or other alkane) chain torsion histograms [Re: rmv]
beginner Offline
Forum Member

Registered: 07/19/04
Posts: 165
Dear Prof. Venable,

Thank you very much.

I am not clear; the following commands produced histogram, didnít they?

! CONVERT TO HISTOGRAMS, ACCUMLATE
set k 1
label klp
mantime g@k hist -180. 180. 360
mantime hg add g@k
mantime b@k hist -180. 180. 360
mantime hb add b@k
mantime q@k hist -180. 180. 360
mantime hq add q@k
mantime d@k hist -180. 180. 360
mantime hd add d@k
incr k by 1
if k .le. @NLPD goto klp
! NORMALIZE FOR NUMBER OF LIPIDS
mantime hb divi @NLPD
mantime hg divi @NLPD
mantime hd divi @NLPD
mantime hq divi @NLPD

As you can see from my fist attachment, I got very low population of trans conformation, even though I made integration. What Did I Do Wrong?

Top
#7281 - 03/09/06 06:42 PM Re: lipid (or other alkane) chain torsion histograms [Re: beginner]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8379
Loc: 39 03 48 N, 77 06 54 W
The problem is your short time interval, and the way the script computes the @MXT value to set the time series length. As a consequence, you've only allocated space for 100 points in a time series, and then proceed to declare time series of 360 points for the histograms. Either analyze more data, or add the following line (in green) to the script:

calc mxt 100 * @N
if mxt lt 360 set mxt 360
calc mxs 10 + ( 4 * @NLPD )
correl maxt @MXT maxs @MXS maxa ?NATOM noupdate
_________________________
Rick Venable
computational chemist


Top
Page 1 of 2 1 2 >

Moderator:  chmgr, John Legato, petrella