Differences

This shows you the differences between two versions of the page.

--- 2d_rdfs [2026/07/04 23:02] – [1) Edit Necessary Files] edit
+++ 2d_rdfs [2026/09/24 00:24] (current) – [1) Edit Necessary Files] edit
@@ Line 1: / Line 1: @@
-[[Analysis]]
+{{:sample_rdf_plots.png?200|}}[[Analysis]]
 =====Brief instructions to obtain 2D radial distribution function (RDF)=====
-Authors: Xiaohong Zhuang, Viviana Monje, Dr. Jeffery B. Klauda
+Authors: Xiaohong Zhuang, Viviana Monje, Dr. Jeffery B. Klauda\\
-Revised by Eric Zhang (2026): handles asymmetric bilayers and complex lipid types (In Progress)
+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.
+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.
-Symmetric Scripts: {{ :2d-rdf.tar.gz |}}\\
-Asymmetric Scripts: To be added
+Legacy Scripts: {{ :2d-rdf.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**\\
+The legacy scripts are also on ZT1 path:  **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/2d-rdf**\\
-=====Symmetric Bilayer 2D RDF Workflow=====
+Modern Python Scripts: {{ :2d-rdf-asym.tar.gz |}} (Recommended)
+=====Legacy Bilayer 2D RDF Workflow=====
 === Part I. ===
@@ Line 116: / Line 115: @@
 .	The final averaged RDF data are in rdf.dat file, Plot the data with RDF vs. radial distance r.
-=====Asymmetric Bilayer 2D RDF Workflow (In Progress)=====
+=====Modern Python 2D RDF Workflow=====
-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.
+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) Edit ''system_setup.py''===
+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).
-Open ''system_setup.py'' and find the marker_atoms dictionary:
 <code>
-marker_atoms = {
+# Example config_dict: a complex sphingolipid membrane with some of my own custom patched NAAA lipids
-    'pe': 'C1', 'pg': 'C1', 'pc': 'C1', 'pi': 'C1', 'pa': 'C1', 'ps': 'C1',
+config_dict = {
-    'lipa': 'OA6', 'ol': 'C1', 'chl': 'O3', 'sm': 'C1S', 'cer': 'C1S'
+    'MEMB':['AIPE','AMPG','DPPC','CER160','CER221','LSM','BSM'],
+    'NAAA':['LLL',],
 }
-</code>
-The keys of this dictionary correspond to the following table and define the atoms used to represent the lipids in the rdf calculations.
-^ Lipid Key ^ Marker Atom ^ Typical Lipid Structural Group ^
-| pe | C1 | Phosphatidylethanolamine |
-| pg | C1 | Phosphatidylglycerol |
-| pc | C1 | Phosphatidylcholine |
-| pi | C1 | Phosphatidylinositol |
-| pa | C1 | Phosphatidic acid |
-| ps | C1 | Phosphatidylserine |
-| lipa | OA6 | Lipid A |
-| ol | C1 | Ornithine (custom lipid type) |
-| chl | O3 | Cholesterol / Sterols |
-| sm | C1S | Sphingomyelin |
-| cer | C1S | Ceramide |
-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.
+# 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'],
+}
-Next, scroll down to the config dictionary and populate the ntype, resnames, top_counts, and bot_counts sections.\\
+# Example of a simple POPE/POPG membrane which have segid MEMB because they are native lipids
-- ntype: Number of resnames/residues in lipid category.\\
+config_dict = {
-- resnames: The exact resi definition of your lipids found in toppar files (case-sensitive).\\
+    'MEMB':['POPE','POPG',],
-- 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
+# corresponding custom_grouping_dict which is empty since you only have normal memb lipids
-<code>
+custom_grouping_dict = {
-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]
-    }
 }
-</code>
-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.
-===b) Edit ''system_setup.scr''===
-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)
-<code>
-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"
-</code>
-Example 2: A membrane where some components only exist in one leaflet (uses "manual" pair mode)
-<code>
-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
-</code>
-===c) Edit ''run_rdf.scr''===
-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:
-<code>
-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
 </code>
 ====2) Run Scripts====
-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.
+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).
-Here are the estimated times for each script.
-^ Script/Jobs ^ Time ^
-|''./system_setup.scr'' | <5 mins |
-|''./run_rdf.scr'' | 10-15 mins |
-|''./coord.csh'' jobs | 0.5-2 hrs |
-|''./run_rdf2.scr'' | <30 secs |
-|''./rdf.csh'' jobs | 0.5-1 hrs |
-|''plot_rdf.scr'' | <30 secs |
-===(Optional)===
+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.
-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''
+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.
-<code>
-# 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
-</code>
-Run ''block_stats.py'' as ''python block_stats.py avg-rdf-p <number of blocks> avg-rdf-p-stats.dat''
-Example
-<code>
-# 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.
-</code>