Previous Thread
Next Thread
Print Thread
Page 1 of 2 1 2
bilayer chain tilt
#7464 07/13/05 01:10 AM
Joined: Sep 2003
Posts: 8,507
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,507
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

Re: bilayer chain tilt
rmv #7465 01/18/06 04:02 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
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.

Re: bilayer chain tilt
beginner #7466 01/18/06 10:49 PM
Joined: Sep 2003
Posts: 8,507
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,507
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

Re: bilayer chain tilt
rmv #7467 01/19/06 02:50 AM
Joined: Sep 2003
Posts: 8,507
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,507
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.

Re: bilayer chain tilt
rmv #7468 01/19/06 01:11 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
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.

Re: bilayer chain tilt
beginner #7469 01/19/06 07:54 PM
Joined: Sep 2003
Posts: 8,507
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,507
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

Re: bilayer chain tilt
rmv #7470 03/27/07 02:45 PM
Joined: Jan 2005
Posts: 86
M
Forum Member
Offline
Forum Member
M
Joined: Jan 2005
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.

Re: bilayer chain tilt
mzzt #7471 03/27/07 04:53 PM
Joined: Sep 2003
Posts: 8,507
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,507
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.

Re: bilayer chain tilt
rmv #7472 11/21/08 02:17 PM
Joined: Jul 2004
Posts: 165
B
Forum Member
Offline
Forum Member
B
Joined: Jul 2004
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
*

Last edited by beginner; 11/21/08 02:39 PM.
Re: bilayer chain tilt
beginner #7473 11/21/08 06:12 PM
Joined: Sep 2003
Posts: 8,507
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,507
I don't see any obvious errors, but I usually allow the output to guide me.

Page 1 of 2 1 2

Moderated by  lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.4
(Release build 20200307)
Responsive Width:

PHP: 5.6.33-0+deb8u1 Page Time: 0.017s Queries: 35 (0.008s) Memory: 0.9870 MB (Peak: 1.1203 MB) Data Comp: Off Server Time: 2020-10-31 09:29:58 UTC
Valid HTML 5 and Valid CSS