Klauda Lab Simulation Wiki

Authors: Xiaohong Zhuang, Viviana Monje, Dr. Jeffery B. Klauda Revised by Eric Zhang (2026): handles asymmetric bilayers and complex lipid types

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

Symmetric Scripts: 2d-rdf.tar.gz
Asymmetric Scripts: 2d-rdf-asym.tar.gz

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

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

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

# ./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

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:

./trim.scr

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

# ./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

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:

extract.scr

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.

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.

# 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

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

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

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.

This workflow uses the asymmetric package instead of the original scripts, and separately averages top and bottom leaflet RDFs. Many of the steps in the original script were combined/automated, but files for the most part can still be run individually following the symmetric workflow to get results. If any steps aren't clear, reach out to ezhang32@umd.edu.

Open system_setup.py and find the marker_atoms dictionary:

marker_atoms = {
    'pe': 'C1', 'pg': 'C1', 'pc': 'C1', 'pi': 'C1', 'pa': 'C1', 'ps': 'C1', 
    'lipa': 'OA6', 'ol': 'C1', 'chl': 'O3', 'sm': 'C1S', 'cer': 'C1S'
}

The keys of this dictionary correspond to the following table and define the atoms used to represent the lipids in the rdf calculations.

If you have lipids that don't belong in any of these categories, you can add your own entry to this dictionary as I have done with ornithine lipids <'ol': C1>. Just confirm that the atom you choose exists in the psf/toppar files. If it's standard diacyl lipid, try to choose a backbone carbon atom at the same level as the C1 atoms of standard phospholipids. If your lipid is bulky or special, you can try other atoms. The idea is to avoid atoms that might overlap each other (since the RDF math only looks at XY distances), which can lead to garbage results.

Next, scroll down to the config dictionary and populate the ntype, resnames, top_counts, and bot_counts sections.
- ntype: Number of resnames/residues in lipid category.
- resnames: The exact resi definition of your lipids found in toppar files (case-sensitive).
- top_counts: The exact number of lipids for each subtype located in the top leaflet.
- bot_counts: The exact number of lipids for each subtype located in the bottom leaflet.

Example: Highly asymmetric membrane with a lipid A outer membrane and phospholipid inner membrane

config = {
    'ntype': {'pe': 1, 'pg': 1, 'lipa': 3},
    'resnames': {
        'pe': ['POPE'],
        'pg': ['POPG'],
        'lipa': ['LIPA4','LIPA5','LIPA6']
    },
    'top_counts': {
        'pe': [0],
        'pg': [0],
        'lipa': [10, 10, 10]
    },
    'bot_counts': {
        'pe': [80],
        'pg': [20],
        'lipa': [0, 0, 0]
    }
}

Note that you still have to write 0 for lipids that don't appear in both leaflets. Ensure the order of the counts mirrors the order of the resnames.

Open system_setup.scr. Populate the lipids array with the lipid types from system_setup.py. If you want the pairwise 2d rdf of every unique pair of lipid types, set pair_mode to “auto”. If you don't want that, then set pair_mode to “manual” and populate the manual_pairs array. They should all follow lip1-lip2 syntax exactly.

Example 1: A membrane that has the same components in both leaflets (uses “auto” pair mode)

set lipids = (pe pg ol pc cer sm)
set pair_mode = "auto"
set manual_pairs = ()
sbatch -A energybio-eng system_setup.csh "$pair_mode" "$lipids" "$manual_pairs"

Example 2: A membrane where some components only exist in one leaflet (uses “manual” pair mode)

set lipids = (pe pg lipa)
set pair_mode = "manual"
set manual_pairs = (pe-pg)
sbatch -A energybio-eng system_setup.csh "$pair_mode" "$lipids" "$manual_pairs"

# The logic behind the manual setting is that you can't compute the rdf of pe-lipa or pg-lip
# if the only outer membrane component is lipid A and the inner membrane is pe and pg for example

Open system_setup.scr. Set the namd and project paths to relative paths that point to the namd directory and project directory with all of your charmm setup files (relative to the current 2d-rdf-asym directory). Next, ensure recenter is set to “yes” for your first run. If you already have an established equilibrated range dcds from a previous run, you can set it to “no” to save time/resource by not recentering again. Similarly, set toppar to “yes” to copy your toppar folder to the right directories, but set it to “no” for future runs. Next, set dcdnum to the number of dcds you want to process. The script uses the next.seqno found in your namd folder to automatically compute the most recent dcd range of size dcdnum.

Example:

set MAN_CHARMM_DIR = "../../"
set MAN_NAMD_DIR   = "../"

# if you are running this for the first time, set this to yes to copy toppar folders to appropriate directories
# if you have already run this once already, set this to "no" to save computational resources
set toppar = "yes"

# Leave this as "yes" for the first run. If you're still waiting on more dcds to finish running. Leave it on "yes"
# Only set to "no" once your namd simulation has finished and you've finalized your dcd range.
set recenter = "yes"

# will compute rdfs for the 20 most recently finished dcds.
set dcdnum = 20

All the hard work is done! Now run ./system_setup.scr, ./run_rdf.scr, ./run_rdf2.scr, and plot_rdf.scr in order once the prior slurm job has finished. ./system_setup.scr creates directories for each pairwise RDF and generates customized input files. ./run_rdf.scr recenters dcds, automates steps 1-5 of the symmetric workflow, and submits coord.csh jobs to Slurm. Once that's done, ./run_rdf2.scr reorders xyz coordinates and submits the rdf.csh jobs. Finally, plot_rdf.scr calls a customizable Python script to separate and plot the top and bottom leaflets, but you can plot however you want using Excel or another program by extracting values from rdf.dat in each lipid directory.

Here are the estimated times for each script.

For those interested in block rdfs, you can check the comments of block.scr to see how smaller ranges of dcds behave. If you have large spikes in your rdfs and it's struggling to converge to 1 for large r, I've also added block_stats.py which analyses the behavior of block rdfs for large r to see if it's just statistical noise or something real.

Example of block.scr

# Say I have 20 dcds (37-56) and I want to split it into 5 blocks
# Use this syntax: ./block.csh <first dcd> <last dcd> <block number (always starts at 1)>
./block.csh 37 40 1
./block.csh 41 44 2
./block.csh 45 48 3
./block.csh 49 52 4
./block.csh 53 56 5

Run block_stats.py as python block_stats.py avg-rdf-p <number of blocks> avg-rdf-p-stats.dat

Example

# ran this: python block_stats.py avg-rdf-p 5 avg-rdf-p-stats.dat
# Sample Output:
--- Tail diagnostic (r > 30 A) ---
Top:  mean deviation from 1 = -0.1609   avg SEM = 0.0632
Bot:  mean deviation from 1 = -0.0520   avg SEM = 0.0629
If |mean deviation| is within ~1-2x the avg SEM, the tail drift
is consistent with statistical noise. If it's much larger than
the SEM, it's a reproducible effect (real depletion/clustering).

# This is mainly diagnostic, and doesn't mean much if your system isn't equilibrated
# The important features of the RDF are that it stays at 0 in the beginning and converges to 1 in the long run.