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.
{{:sample_rdf_plots.png?200|}}[[Analysis]]

=====Brief instructions to obtain 2D radial distribution function (RDF)=====
Authors: Xiaohong Zhuang, Viviana Monje, Dr. Jeffery B. Klauda\\
Refactored by Eric Zhang (2026): handles any membrane system

These are the brief instructions to calculate the radial distribution function (RDF), also known as a pair correlation function. New users should use the modern Python scripts. they have been backtested and validated against the legacy scripts, and run 95% faster. It's also easier to use. 

Legacy Scripts: {{ :2d-rdf.tar.gz |}}\\
**Note**: The example is made in DT2. It is provided for reference purpose. 

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

Modern Python Scripts: {{ :2d-rdf-asym.tar.gz |}} (Recommended)
=====Legacy Bilayer 2D RDF Workflow=====

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

=====Modern Python 2D RDF Workflow=====
This workflow automates the entire process and runs much faster. If any steps aren't clear, reach out to ezhang32@umd.edu.

====1) Edit Necessary Files====
a) Open ''config.py'' and read the comments I left throughout. Then, adjust the fields to work for your membrane system. I would not recommend you change num_blocks (which is set at 1 by default) because it is already highly optimized and can process 20 dcds in under a minute. There's no point in trying to parallelize the workflow further. Also, custom_grouping_dict should only be modified if you have non segid MEMB lipids (check psf). 
<code>
# Example config_dict: a complex sphingolipid membrane with some of my own custom patched NAAA lipids
config_dict = {
    'MEMB':['AIPE','AMPG','DPPC','CER160','CER221','LSM','BSM'],
    'NAAA':['LLL',],
}

# corresponding custom_grouping_dict which handles non-memb lipids like LLL, where LLL is grouped under an arbitrary label called OL
custom_grouping_dict = {
    'OL':['LLL'],
}

# Example of a simple POPE/POPG membrane which have segid MEMB because they are native lipids
config_dict = {
    'MEMB':['POPE','POPG',],
}

# corresponding custom_grouping_dict which is empty since you only have normal memb lipids
custom_grouping_dict = {
}
</code>

====2) Run Scripts====
a) Run ''./0_setup.py''. If you set recenter_dir = None in config.py, you should see a recenter dir/scripts get generated. Everyone should see rdf.csh and system_info.txt be generated. If there are errors, simply fix config.py and rerun (I added error handling, so you should see some error messages in the terminal for common mistakes).

b) Run ''./1_rdf.scr''. The first time you run it, it will take a minute to set up the Python virtual environment and install MDAnalysis. Then, it will execute the RDF Slurm job. The entire run should take less than a minute, or longer for really big systems.

c) Once the slurm job finishes, run ''./2_analysis.py'' and ''./3_plot_rdf.py''. You should now see .dat files and a combined figure stored inside the results directory. Feel free to adjust the plot settings if you're familiar with matplotlib, or take the .dat files into Excel or any other software and plot them differently. It also averages the top and bottom leaflets and plots the averages automatically if you specified symmetric = True in config.py.