Klauda Lab Simulation Wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
[[Analysis]]

**Brief instruction to obtain tilt angle distribution** \\
Authors: Xiaohong Zhuang and Viviana Monje\\
Updated: Dr. Jeff Klauda (2021)


This brief instruction is used to calculate tilt angle distribution of sterols/cholestrols. The scripts that you need are in ZT1 path: **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/tilt
**
or using zip file: {{ :tilt.gz |}}

In each sbatch file, the ZT1 charmm are provided. 

Copy the scripts folder “tilt” from path highlighted above to your directory that have all the dyn files.

1. Update the //rtfprm.str//, if needed and copy toppar files.

Copy the toppar.str and the toppar subdirectory from the CHARMM-GUI build directory (for the CHARMM files not NAMD). Then, update the PSF and CRD file in the rtfprm.str associated with your files, if needed.


2.	Update the //set-range.str // 
<code commend>
set rnc = CHL1   ! Update the resname of sterol
set ti = 10      ! Update the dyn number first trajectory
set tf = 11      ! Update the dyn number of last trajectory
set tra = 18     ! Update the value of (first trajectory plus 8) to define unit# to be read
</code>

Note:\\
a. Choose ti and tf based on the equilibrated data. So if SA/lipid shows that you have equilbration at dyn10 and higher then use ti=10 and tf is the number of DCD files.\\
b. In this example the sterol of interest is CHL1. Other sterols can be used and rnc is changed based on their resname.\\
c. We use @tra=@ti+8 because in tilt.inp, the first unit number is calculated by calc u = 8 + @i \\
… \\
''traj firstu @tra nunit @N''\\


3.	Update the  //tilt.inp// 

3.1	Update the data points per DYN/dcd file\\
<code commend>
calc mxs = 10 + 16*100 ! TIME SERIES PER MOLECULE +10  
calc mxt = 500 * @N  ! update first number to data points per DYN file
</code>


Note:\\
For mxs, it depends on the lines of calculation performed (in 1.2 below) and for most cases does not need to change. You can increase the value if the error occurs because the space (time series) limitation.
 

For mxt, 500 is # of data points/frames per DYN file (calculated by # of total time steps divided by dcd frequency, both can be obtained from dyn2-inp, e.g. numsteps/dcdfreq=500)


3.2 Update the reference atoms to calculate the tilt angles 


<code commend>
enter r@y vect r MEMB @t C3 MEMB @t C17
enter t@y vect z MEMB @t C3 MEMB @t C17
</code>


Note:\\
MEMB is segid @t is resid (which is based on the  resid range of  top1,top2,bot1 and bot 2 set up in set_range.str), C3 (first C close to the OH group near head end) or C17 (carbon at the "base" of the last rings near tail end) are atom type names selected (which is more clear when you look at the topology file). The reference vectors are created based on connecting the reference atoms C3 and C17. Then, the angle between reference vector and the bilayer normal are calculated. 

3.3  Update location of DCD file if needed.

<code commend>
calc u = 8 + @i
open unit @u file read name ../dyn@i.dcd
incr i by 1
</code>

Update the path for the DCD file if needed.

4.	Update the //top.inp// and //bot.inp//


The top/bot.inp are the top and bottom leaflet variable input data files for dist.exe. It is executed in the //dist.csh// file (More detail about //dist.csh// will be explained in later section) by the following command lines: \\
<code commend>
dist.exe < top.inp
dist.exe < bot.inp
</code>

The variable values in top.inp are the following (which can be found in dist.c file):
<code commend>
1000        ! # of total data points/frames (= # frames per dcd times # dyn, e.g. 500 X 2)
90.10         ! Maximum (value, which is the one step larger than the real maximum (90+0.10)
0                    ! Minimum (angle) value
t-all-nc1.dat  ! input file name
t-tilt.dat        ! output file name with the distribution
</code>


Note:\\
The only change you likely need is the fist line based on the number of frames per DCD and number of DCD file. The other lines are described here. //t-all-nc1.dat// is same as //t-all.dat// but with the first time step columns removed. Since //dist.c// calculate distribution (with bin size 0.1) from the first column, the input file is //t-all-nc1.dat// (you won’t see it as it is deleted after the calculation) is used instead of t-all.dat which having the first time steps columns.


//bot.inp// is similar to above,but with b-all-nc1.dat and b-tilt.dat instead of 
//t-all-nc1.dat// and //t-tilt.dat//


5.	 Update //gettilt.csh//
 

5.1 Update the charm location if needed but the current working code is typically in this file:\\

<code commend>
# In DT2\\
set chm = **/homes/jbklauda/charmm/c43b1/exec/gnu/charmm** \\
</code>

5.2 Update the path to where you want to save your final tilt distributions file, it can be a common folder than you want to save all your tilt angle data, or it can be current folder path.


''set all = "/lustre/jbklauda/lipid/analysis-test/gl-chl1/tilt"''


It is recommend to save a copy of final result to a common directory, which may help to organize/plot the data later if you analyze a lot of systems)


Note:\\
In the following commands:

<code commend>
# INVOKE CHARMM
$''chm N=$1 < tilt.inp > tilt.out'' 
# GET DISTRIBUTION 
# especify last and first tilt data file 
''./dist.csh $2 $3''
</code>


The variables $1, $2, $3 $4 \\
$1 in the line: ''$chm N=$1 < tilt.inp > tilt.out ''\\
$2 $3 in the line: ''./dist.csh $2 $3 ''\\
$4 in the line: ''cp -f avg-tilt.dat $all/$4.dat ''\\
 are read from the variable specification when execute the //gettilt.csh// by //get.scr// which is explained in the following section.


6.	Update //1_get.scr//


As mentioned above, the number and file name following //gettilt.csh// are the variables used in //gettilt.csh//
 

''sbatch -A energybio-hi gettilt.csh 2 20 1 chl1'' \\


Note: 


As shown in the 1_get.scr file, the four variables are \\
2: Total number of dcd files to be read, which will be value of N in //gettilt.csh// \\
20: last tilt file number, which is the number of lipids of interest per leaflet (not the total), \\
1: first tilt file number \\
chl1: name to copy final dat file to common directory \\


7.  Submit the job by command: ''1_get.scr'' 


//T1.dat …t#.dat and b1.dat … b#.dat// will be obtained, and the final top and bottom leaflet averaged tilt angle distribution is provided in //avg-tilt-nh.dat//. The first column is angle, and the second column is normalized probability in each angle.


As the probability is already normalized, the sum of probability should equal to 1. So to double check the distribution, check the file 1_sump.dat ( the three columns are: angle, probability, accumulative sum of probabilities) to make sure the end of third column equal to 1.


8.Plot the //avg-tilt-nh.dat// to visualize the distribution in a program of your choice. 


//**Some additional details of the codes here:**//


This section is to aid in the understanding the codes for this tilt calculation. 


1.	//combine.csh// is used to combine the data of  all the lipids, which is executed by the following command lines in //dist.csh//

<code commend>
./combine.csh t $last $first # (leader) (last #) (first #)
./combine.csh b $last $first
</code>

2.	The numbers of 20 and 1 in //get.scr// are defined as the variables
 $2 $3 in the command line in the //gettilt.csh// file, 
 ''./dist.csh $2 $3 ''\\

but they are defined as $1 $2 in the command lines in the //dist.csh//:
<code commend>
set last = $1 # last tilt data file
set first = $2 # first tilt data file
</code>


because the variable number are defined based on the order of specification when a script is executed. Since 4 variables are specified when gettilt.csh is excecuted by command  line in //1_get.scr// \\
''sbatch -A energybio-hi gettilt.csh 2 20 1 chl1''


While only 2 variables are specified when //dist.csh// is executed by command line in //gettilt.csh// \\
''./dist.csh $2 $3''