Klauda Lab Simulation Wiki

1.1 Update the rtfprm.str Update the path for the topology stream file and the psf and pdb files

e.g. 1.1 In rftprm.str, specify the path of the toppar.str file (it is toppar_less_1.str for hypo system), psf and crd file

1.2 In the toppar_less_1.str file, Replace the short path: toppar/ by the full path to your CHARMM toppar folder e.g. /lustre/xzhuang/soybean/hypocotyl-b1/toppar/

Note: “Replace” command in WinSCP would save time to replace the short paths by full path. Since the toppar.str is streamed in the current path, the system will try to find the toppar folder from the current path. So in order to let system to use the correct path of the *rtf and *prm files, the full path is provided in toppar.str file.

2. Update the set-range.str

set top1 = 145 ! first resid of lipid of interest on the top leaflet
set top2 = 165 ! last resid of lipid of interest on the top leaflet
set bot1 = 166 ! first resid of lipid of interest on the bottom leaflet
set bot2 = 186 ! last resid of lipid of interest on the bottom leaflet
set lip = 21   ! total number of lipids per leaflet
set ti = 41     ! first trajectory of interest
set tf = 75    ! last trajectory of interest
set tra = 49    ! first trajectory ti plus 8 to define the first unit# 

Note: a. The resid range can be found in the step5_assembly.crd file, usually the first half resid is on the top leaflet, and rest half on bottom leaflet in hypo system. You may confirm the top and bottom leaflet by VMD

b. The lipid of interest for soybean system is beta-sitosterol (SITO) and stigmasterol (STIG) for hypo system. The above setup is for SITO, which has 21 SITO in each leaflet. Assume system reach equilibrium by dyn 41, so dyn range 41 to 75 are used to calculate average angle distribution.
c. We use @tra=@ti+8 because in tilt.inp, the first unit number is cacluated by calc u = 8 + @i
…
traj firstu @tra nunit @N

3. Update the tilt.inp

1.1 Update the data points per DYN/dcd file

calc mxs = 10 + 16*100 ! TIME SERIES PER MOLECULE +10  
calc mxt = 1000 * @N  ! update first number to data points per DYN file

Note: For mxs, it is just used to save space, it depends on the lines of calculation performed (in 1.2 below). It does not need to be exact. Can increase the number if the error occurs because the space (time seriers) limitation.

For mxt, 1000 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. 1000000/1000=1000)

1.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.

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):
35000 ! # of total data points/frames (= # frames per dcd times # dyn, e.g. 1000 X 35)
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:
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 depends on you run the analyses on DT2. Use # to comment out the one you don’t use

# In DT2
set chm = /lustre/xzhuang/software/charmm/gnu/gfort-xxlg

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/xzhuang/soybean/all_analyses_summary/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 35 21 1 hypo-b1

Note:

As shown in the get.scr file, the four variables are
75: Total number of dcd files, which will be value of N in gettilt.csh
21: last tilt file # which is the # of lipid of interest,
1: first tilt file #
hypo-b1: name to copy final dat file to common directory

After you run get.scr command $chm N=$1 < tilt.inp > tilt.out, the files t1.dat …t#.dat, b1.dat…b#.dat will be created. Each t#.dat is a tilt file, the # of tilt files depends on the number of lipid of interest. Therefore, last tilt file # is the # of lipid of interest.

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

Additonal note: (This section only written for helps understanding the code, so it is fine if you don’t understand since no modification is required)

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 21 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 35 21 1 hypo-b1

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