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,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
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,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
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,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
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,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
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,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
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,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
I don't see any obvious errors, but I usually allow the output to guide me.

Re: bilayer chain tilt
rmv #27023 04/04/11 06:23 PM
Joined: Feb 2011
Posts: 31
J
Forum Member
Offline
Forum Member
J
Joined: Feb 2011
Posts: 31
Hi Dr. Venerable,

why not just use coor helix?

Re: bilayer chain tilt
jeremy adler #27024 04/04/11 06:55 PM
Joined: Sep 2003
Posts: 8,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
COOR HELIX does not process trajectory files; in general, an analysis using CORREL will be much more efficient than a loop over frames.

There often can be a couple different ways to do a particular analysis, and sometimes it's useful to try more than one and compare the results.


Rick Venable
computational chemist

Re: bilayer chain tilt
rmv #27026 04/04/11 07:07 PM
Joined: Sep 2003
Posts: 4,793
Likes: 2
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,793
Likes: 2
There is a helix timeseries in correl, but this is (just as COOR HELIx) intended for helical structures, and may not be appropriate for a fatty acid hydrocarbon tail.


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Re: bilayer chain tilt
lennart #27030 04/04/11 07:46 PM
Joined: Sep 2003
Posts: 8,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
In either case, it should be noted that the tilt calculation is only meaningful when the chains are fairly ordered, such as more gel-like states.


Rick Venable
computational chemist

Re: bilayer chain tilt
rmv #34530 11/02/14 06:08 PM
Joined: Jan 2014
Posts: 8
S
Forum Member
Offline
Forum Member
S
Joined: Jan 2014
Posts: 8
Dear rick,

Hello, I have a question:
Actually you know, DPPC lipid has 16 carbons in each of chain.
But, script in this archive can calculate tilt angle from C3 to C13.
Are they have specific reason?
Why you don't calculate from C1 to C16?

Re: bilayer chain tilt
rmv #34531 11/02/14 06:32 PM
Joined: Sep 2003
Posts: 8,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
The free ends of the chains are highly disordered, and the upper ends near the glycerol are influenced by the ester binding geometry and headgroup packing.


Rick Venable
computational chemist

Re: bilayer chain tilt
rmv #34600 11/28/14 05:36 AM
Joined: Jan 2014
Posts: 8
S
Forum Member
Offline
Forum Member
S
Joined: Jan 2014
Posts: 8
Dear Rick,

Thanks, I understood that in case of C3.
However, I think that setting the end to end vector (C3-C16) is more reasonable to check the tilt angle.
Is it does not matter if all values are calculated from C3 to C13?

Last edited by shanKim; 11/28/14 05:37 AM.
Re: bilayer chain tilt
rmv #34601 11/28/14 05:12 PM
Joined: Sep 2003
Posts: 8,498
rmv Online Content OP
Forum Member
OP Online Content
Forum Member
Joined: Sep 2003
Posts: 8,498
What matters most is that you are consistent in your methods, and clearly state what you've done in reports and publications.

Because of the flexibility of the free chain ends, I do not agree that using the terminal methyl is more reasonable; if you are concerned, do the calculation both ways and see how much difference it makes.


Rick Venable
computational chemist

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.015s Queries: 50 (0.006s) Memory: 1.0651 MB (Peak: 1.2633 MB) Data Comp: Off Server Time: 2020-09-24 23:41:28 UTC
Valid HTML 5 and Valid CSS