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
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
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
calc mxs = 10 + 16*100 ! TIME SERIES PER MOLECULE +10 calc mxt = 500 * @N ! update first number to data points per DYN file
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
enter r@y vect r MEMB @t C3 MEMB @t C17 enter t@y vect z MEMB @t C3 MEMB @t C17
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.
calc u = 8 + @i open unit @u file read name ../dyn@i.dcd incr i by 1
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:
dist.exe < top.inp dist.exe < bot.inp
The variable values in top.inp are the following (which can be found in dist.c file):
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
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:
# In DT2\\ set chm = **/homes/jbklauda/charmm/c43b1/exec/gnu/charmm** \\
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:
# INVOKE CHARMM $''chm N=$1 < tilt.inp > tilt.out'' # GET DISTRIBUTION # especify last and first tilt data file ''./dist.csh $2 $3''
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 10 1 chl1
Note:
As shown in the 1_get.scr file, the four variables are
2: Total number of dcd files, 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
./combine.csh t $last $first # (leader) (last #) (first #) ./combine.csh b $last $first
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:
set last = $1 # last tilt data file set first = $2 # first tilt data file
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