Page 1 of 2 1 2 >
Topic Options
#7464 - 07/12/05 09:10 PM bilayer chain tilt
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
The following script uses a vector from C3 to C13 of each lipid chain to estimate the overall chain tilt.


* DPPC 80 lipid gel; chain tilt time series
*
bomblev -2

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

! COMPUTE TIME SERIES LIMITS; 1 NS SAVED AT 1 PS INCR IN EACH FILE
calc mxt = 1000 * @N
! FOUR VECTORS PER LIPID + TEN MORE
calc mxs = 10 + ( 4 * 80 )

correl maxt @MXT maxs @MXS maxa ?NATOM noupdate
! DECLARE ACCUMULATION SERIES (FOR AVGS) FIRST
enter aav zero
enter bav zero
! LOOP OVER LIPIDS; DEFINE VECTOR LENGTH (R) AND Z COMPONENT
set k 1
label elp
enter r@K vect r L @K C23 L @K C213 ! BETA CHAIN
enter s@K vect r L @K C33 L @K C313 ! GAMMA CHAIN
enter z@K vect z L @K C23 L @K C213
enter w@K vect z L @K C33 L @K C313
incr k by 1
if k le 80 goto elp

! OPEN FILES (COPIED TO /scratch FOR EFFICIENCY)
set j 1
label olp
calc u = 7 + @J
open unit @U file name /scratch/dyn@J.trj
incr j by 1
if j le @N goto olp
traj first 8 nunit @N

set k 1
label tlp
! NORMALIZE Z COMPONENT (DIVIDE BY R)
mantim z@K ratio r@K
mantim w@K ratio s@K
! IGNORE SIGN; NEEDED TO AVG OVER BOTH LEAFLETS
mantim z@K abs
mantim w@K abs
! SUM OVER CHAINS, EACH SEPARATELY
mantim aav add z@K
mantim bav add w@K
incr k by 1
if k le 80 goto tlp

! DIVIDE BY N; RESULT IS AVG Z PROJECTION, WHICH IS COS(THETA)
mantim aav divi 80.
mantim bav divi 80.
! CONVERT COS TO AN ANGLE
mantim aav acos
mantim bav acos

! OUTPUT FILE WITH (TIME, BETA-TILT, GAMMA-TILT) FOR EACH TIME POINT
edit aav veccod 2
open unit 2 write card name tiltavg.dat
write aav unit 2 dumb time
* dummy
*

end
stop
_________________________
Rick Venable
computational chemist


Top
#7465 - 01/18/06 11:02 AM Re: bilayer chain tilt [Re: rmv]
beginner Offline
Forum Member

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

I would like to calculate the tilt angles in the upper/lower leaflet separately. Could you please guide me how to do that?

Thank you very much.

Top
#7466 - 01/18/06 05:49 PM Re: bilayer chain tilt [Re: beginner]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
Most of the script would be the same, only how the individual chain vectors are averaged changes. The loop with the label 'tlp' becomes

set k 1
label tlp
mantim z@K ratio r@K
mantim w@K ratio s@K
mantim z@K abs
mantim w@K abs
if k gt 40 goto skp1
mantim aav add z@K
mantim aav add w@K
goto skp2
label skp1
mantim bav add z@K
mantim bav add w@K
label skp2
incr k by 1
if k le 80 goto tlp

In this case, it is assumed that lipids 1 thru 40 are in one leaflet, and lipids 41 thru 80 are in the other.
_________________________
Rick Venable
computational chemist


Top
#7467 - 01/18/06 09:50 PM Re: bilayer chain tilt [Re: rmv]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
A final note on this script for chain tilt-- the value is likely only meaningful for systems with considerable collective tilt, usually with mostly trans alkane chains packed together, as in a gel phase. For the more disordered fluid phase, this value may not be very useful at all.

Top
#7468 - 01/19/06 08:11 AM Re: bilayer chain tilt [Re: rmv]
beginner Offline
Forum Member

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

Thank you very much for your advice.

Unfortunately, I got this following error message
fmt: end of file
apparent state: unit 5 (unnamed)
last format: (A)
lately reading sequential formatted external IO
Abort

I doubt it came from reading trajectories. However, it went through when I computed both leaflet at the same time. How could I solve this problem?

Thank you very much.

Top
#7469 - 01/19/06 02:54 PM Re: bilayer chain tilt [Re: beginner]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
It's hard to say-- although UNIT 5 is standard input, which indicates a problem reading the input itself, and not some other file. Note that because of the loops, there may be problems using an MPI-based version for this; two solutions are: [1] use a non-parallel version of CHARMM, or [2] use the STREAM command to read the input.
_________________________
Rick Venable
computational chemist


Top
#7470 - 03/27/07 10:45 AM Re: bilayer chain tilt [Re: rmv]
mzzt Offline
Forum Member

Registered: 01/03/05
Posts: 86
Dear Charmmers,

I am interested in calculating the angle between the helical axis of peptide and the bilayer normal.
If I replace all the DMPCs with the peptide in which I am interested, I get the warning mesg. Am I going in the right direction of using this above script for calculating angle between helix axis and membrane normal?

regards,
MzZt.

Top
#7471 - 03/27/07 12:53 PM Re: bilayer chain tilt [Re: mzzt]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
Probably not-- the above is for averaging many vectors, and a helical axis is only a single vector. A time series of the principal axis of the helix backbone atoms from an INERTIA time series would be more what you need; the angle can be computed from the unit vector Z projection.

Top
#7472 - 11/21/08 09:17 AM Re: bilayer chain tilt [Re: rmv]
beginner Offline
Forum Member

Registered: 07/19/04
Posts: 165
dear rick,

as I would like to focus on the distribution of that angle,
let say, I would like to get the histogram from 0-180 degree on the x-axis and probability on the y-axis.

and I did it via following script (based on your script), is it correct?

enter aav zero
edit aav total 180 delta 1. offset 0 skip 1
enter bav dupl aav

! LOOP OVER LIPIDS; DEFINE VECTOR LENGTH (R) AND Z COMPONENT
set k 1
label elp
enter r@K vect r L @K C23 L @K C213 ! BETA CHAIN
enter s@K vect r L @K C33 L @K C313 ! GAMMA CHAIN
enter z@K vect z L @K C23 L @K C213
enter w@K vect z L @K C33 L @K C313
incr k by 1
if k le 224 goto elp
! OPEN FILES (COPIED TO /scratch FOR EFFICIENCY)
open unit 8 file name lipid.dcd
traj first 8 nunit 1 skip 2500
set k 1
label tlp
! NORMALIZE Z COMPONENT (DIVIDE BY R)
mantim z@K ratio r@K
mantim w@K ratio s@K
! IGNORE SIGN; NEEDED TO AVG OVER BOTH LEAFLETS
mantim z@K acos
mantim w@K acos
mantim z@K hist 0 180 180
mantim w@K hist 0 180 180
mantime aav add z@K
mantime bav add w@K
incr k by 1
if k le 224 goto tlp
mantim aav divi 224.
mantim bav divi 224.
!mantim aav acos
!mantim bav acos
! OUTPUT FILE WITH (TIME, BETA-TILT, GAMMA-TILT) FOR EACH TIME POINT
edit aav total 180 delta 1 offset 0 skip 1 veccod 2
open unit 2 write card name cahin_distribu.dat
write aav unit 2 dumb time
* dummy
*


Edited by beginner (11/21/08 09:39 AM)

Top
#7473 - 11/21/08 01:12 PM Re: bilayer chain tilt [Re: beginner]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8373
Loc: 39 03 48 N, 77 06 54 W
I don't see any obvious errors, but I usually allow the output to guide me.

Top
Page 1 of 2 1 2 >

Moderator:  chmgr, John Legato, petrella