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 instructions to obtain 2D radial distribution function (RDF)**\\
Authors: Xiaohong Zhuang, Viviana Monje, Dr. Jeffery B. Klauda

This is the brief instructions to calculate the radial distribution function (RDF), also known as a pair correlation function.

Downloaded Scripts: {{ :2d-rdf.tar.gz |}}

**Note**: The example is made in DT2. It is provided for reference purpose. 

The scripts are on ZT1 path:  **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/2d-rdf**\\

=== Part I. ===
First of all, please check the dcdfreq and xstfreq in your //dyn-2.inp// file, if the dcdfreq is equal to xstfreq (which is the case almost all the time), then do **A** and skip **B**;  if the dcdfreq is not equal to xstfreq, then skip **A** and do **B** only. The **A** and **B** session both generate the //*-cut.dat//. (Since the coord of each lipid is based on each frame, so we want to have the number of x-dimension data point agree with number of frame. You may know this better after you finish this rdf analysis.) 


**A**. Update the //trim.scr//

<code commend>
# ./trim.csh (first dcd) (last dcd) (path of xst file)

# UPDATE the three variables mentioned above
set f = 10
set l = 11
set p = /lustre/jbklauda/lipid/analysis-test/gl-chl1
</code>

**Note**: The p is the charmm-gui path that has the //dyn*.xst.gz// files and f and l are the first and last dyn*.dcd file.

Run //trim.scr// by command:
<code commend>./trim.scr</code>
This will get the x-dimension data for each frame. (If you do **A**, you don’t need the extract folder, please ignore it.) The //dyn*-cut.dat// will be generated in rdf folder after run.

**B**.  Make sure you have copied the extract folder and then go to extract folder.

Update the //extract.scr//

<code commend>
# ./extract.csh (first dcd) (last dcd) (path of psf file)

# UPDATE the three variables mentioned above
set f = 1
set l = 75
set p = /lustre/aou1/soybean/hypocotyl-b1
</code>

**Note**: The p is the charmm-gui path that has the //step5_assembly.psf// file and f and l are the first and last dyn*.dcd file.

Update the path of topology files in //gen_crd.inp// and //extract.inp// if needed.

Run //extract.scr// by command:
<code commend>extract.scr</code>

This does two steps,  1) extract the crd file (The extract crd is used as reference crd file instead of //step5_assembly.crd//, just to make the image of lipid agree in all dcd.) 2) x-dimension data for each frame is extracted.  The //dyn*-cut.dat// will be generated in rdf folder after run.


===Part II.===

1.	Copy the toppar.str and toppar subdirectory from the CHARMM-GUI folder (CHARMM input files) to the 2d-rdf directory. Also, update the psf and crd file in //psfcrd.str//. This must be an absolute path so do not use relative paths for the psf/crd files.

2.	In 2d-rdf folder, make a directory called "inp" if you don’t have it (to store the input file that will be read into the *.exe program when calculating the RDF between lipids).

3.	Update the size and data point information in //gen-inp.csh// if necessary.
Obtain the x dimension value of the last value from //dyn.xsc// (in your dyn folder), first number is about 2/3 of x-dim in simulation box e.g. if the value is 72.1711, then 72*2/3=48≈50, then 50 is used as the first number.

<code commend>
# first number is ~2/3 of x-dim in simulation box 
# second number is number of points in g(r)
# ~0.1 incremets across (r) so funct looks smooth
  echo "50 500" >> inp/inp$k
</code>

In the example, assure two third of your system size is about 50 Å (the first value above), then 50/0.1=500, so you use 500 points (the second value above), if it is 60 Å, then use 600 points, etc.

4.	Update the dyn file range in gen-inp.scr and then run command //./gen-inp.scr// to generate input file for RDF calculation.

You will have many pairs of heagroup that you want to calculate based on your system. Please work on the folders with same lipid headgroup first, after they are all finished, then work on the different (mixed) lipids because we need to find the coordinate of each lipid first, and then the mixed lipids pair just using the coordinated of lipid head group already calculated.  In the example, after finish the calculation in same lipid pair, e.g. pc,pe,pi,sterols, then we start to work on pc-pe, pc-pi,pc-sterol, etc.

The scripts for same lipid pair for glycerol (please use eg_same-glycerol as script template) and sterol lipids (eg_same-sterol), and also the mixed lipids pairs are slightly different (eg_diff). So for each lipid headgroup pair, please copy the correct folder/ scripts as template and rename the folder based on the headgroup names in your system (e.g. use following command to copy the glycerol script folder for your pc lipid group:  ''scp -r eg_same-glycerol pc'' ), and then do the following steps: 

5.	**In the each of same (single) lipid folder (eg. pc,pe,pi,sterol etc.)**, Update the residue names and number of types of each headgroup in //coor-*.inp//

<code commend>
set rn1 = PLPE                                           ! Update the resname
set rn2 = LLPE                                           ! Update the resname
set rn3 = L2PE                                           ! Update the resname

set Ntyp = 3   ! number of types of lipids, which should be consistent as resnames types above
</code>

5.1	Update the job name in //coord.csh// and //rdf.csh//  (This just to distinguish to make debugging easier if error occurs.) Also, copy the toppar.str and toppar subdirectory into the respective eg* directories.

5.2	Updating //coord.scr// by update the range and path of DCD files to be read.  If needed, update the charmm version used in coor.csh.

5.3     Run command ''./coord.scr'' 

After run //coor.scr//, the x,y, and z coordinate of specified atom (P or C2 for glycerol lipids (selecting baesd on need), O3 for sterols, etc.) of each lipid in each frame are in xyz-*.dat.

5.4	Update the dcd range in //rdf.scr// as you will compute the RDF only for the equilibrated portion of the trajectory //rdf.csh// (first dcd) (last dcd) (first equilibrated dcd). The (first equilibrated dcd) can be the same as the (first dcd), but if you wan the rdfs for all dcds but only do averages for a equilibrated section then (first dcd) would be 1 and (irst equilibrated dcd) would be the first dcd in the equilbrated section. Once the values are entered run ''./rdf.scr''.


6.	**After the same headgroup calculation is finished, then we can start to work for the different (mixed) lipid folder(eg. pc-pe,pc-pi,pc-sterol etc.),** 
 
6.1	Update the lipid head group (folder) name in //gen-in.scr//, edit the range of r as step 3 of Part II, and then run ''gen-in.scr''\\
6.2	Update the dcd range in //rdf.scr// as you will compute the RDF only for the equilibrated portion of the trajectory rdf.csh (first dcd) (last dcd) (first equilibrated dcd). And then run ''rdf.scr''.

Note: You don’t have //coor.csh// nor //coor.scr// since the all the coordinate file that we need are already calculated in single lipid folder. All we need to do is to generate the new input files that select the specified coordinate files (in order to calculate 2d rdf in the next step).  The //gen-inp.csh// for same lipid (one directory out) which select only one //xyz-p*.dat// file, which is different with //gen-inp.csh// in mixed lipid folder, which select two //xyz-p*.dat// files of different lipids.

7.	 If you are interested in the block 2d rdf instead of each dcd, update the dcd range in each block, and run ''block.scr'', the block values will be shown in //avg-rdf-p.dat//. The number of block can be selected based on your interest. (e.g. 5 to 10, or any others)

8.	The final averaged RDF data are in rdf.dat file, Plot the data with RDF vs. radial distance r.