Differences

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

--- comp-sa [2026/06/08 23:13] – [Multi-lipid SA Average] edit
+++ comp-sa [2026/09/03 09:56] (current) – [Step 1: edit config.py] edit
@@ Line 7: / Line 7: @@
 Updated Authors: Yalun Yu and Jeffery Klauda\\
 Updated rtfpsf folder: Joshua Lucker \\
-Updated to support asymmetric bilayers: Eric Zhang
+Completely refactored, supports any membrane: Eric Zhang
-Systems involving lipid A should use the asymmetric bilayer package, even if the lipid A is symmetric. If there are any questions or issues with the code, reach out to ezhang32@umd.edu (may need fine-tuning).
+All new users should use the modern Python workflow. It is highly customizable, works for any membrane system, and runs 95% faster than the legacy scripts. It has been backtested with data and reproduces numbers to 4 decimal places of accuracy compared to the legacy scripts. If there are any questions or issues with the code, reach out to ezhang32@umd.edu (may need fine-tuning).
-The script is in ZT1 path: **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/area_mult_lipids** \\
+Legacy scripts are in ZT1 path: **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/area_mult_lipids** \\
 or using zip file: {{ area_mult_lipids.gz |}}
-The script for asymmetric systems is in the zip file : {{ area_asym_lipids.tar.gz |}}
+Modern Python scripts can be found here: {{ area_asym_lipids.tar.gz |}}
-**How to Edit and Run the Scripts (Symmetric Bilayer Only):**
+**How to Edit and Run the Scripts (Legacy Instructions, Discouraged):**
 **Edit the scripts:** \\
@@ Line 167: / Line 167: @@
 Check that the system go to equilibrium by comparing the distribution data, and adjust the dyn file range when necessary. Also check the distribution, the edge is not 0, increase the area range for broader distribution.
-** Modified Workflow for Asymmetric Bilayers **\\
+=====Modern Python Workflow=====
+This workflow is designed to work for any symmetric, asymmetric membrane system. The current copy available for download has predefined inputs for Charmm's native glycerolipids, sphingolipids, sterols, and 4-, 5-, and 6-chain lipid A. If you have a question, please reach out to ezhang32@umd.edu.
-**1. edit system_setup.py **\\
+Note: anytime you run into "no permission" issues, just run ''chmod u+x *''
-''system_setup.py'' auto-generates ''def_nlip_(top/bot).str'', ''*type_(top/bot).str'', and ''lipid_(top/bot).str'' which are essential CHARMM stream and data files needed for Comp-SA analysis. To customize the membrane composition, edit the top_config and bot_config blocks near the bottom of the file. I'll summarize the essentials, but feel free to read the comments in the file for a more detailed breakdown.
-Let’s start with the top leaflet. Identify how many lipid classes it has. Each lipid class has a resname and type prefix.
+====Step 1: edit config.py====
-^ Class ^ Resn Prefix ^ Type Prefix ^ Typical Lipids ^
+Open ''config.py''. Read the comments carefully and fill in the specified fields. As mentioned in the file, you probably don't need to touch fields in the last section, but you have to change config_dict, dyn_dir, charmm_dir, and dcd_range. The dcd_range should be the equilibrated dcd range you found from the geteq analysis. You also do not need to specify how many lipids are in each leaflet; that is automatically determined at runtime. 90% of users don't need to worry about this, but for those working with patched lipids, add only the subresidues with marker atoms.
-| Sterol | ''rnc'' | ''ctype'' | CHOL, ERG |
-| Glycerophospholipid | ''rng'' | ''gtype'' | DPPC, POPE, POPG |
-| Ornithine lipid | ''rno'' | ''otype'' | - |
-| Sphingolipid | ''rns'' | ''stype'' | CER, LSM, BSM |
-| 4-chained lipid A | ''rnl4'' | ''l4type'' | - |
-| 5-chained lipid A | ''rnl5'' | ''l5type'' | - |
-| 6-chained lipid A | ''rnl6'' | ''l6type'' | - |
-| 7-chained lipid A | ''rnl7'' | ''l7type'' | - |
-**(i) Edit the counts dictionary:** Edit this in the config blocks. This holds the various types/class of lipids you have and the number of different lipids of each type.
+I embedded a lot of error handling in the scripts, so don't worry too much about precision if something isn't clear.
-Example:
+Example Membrane with just POPE/POPG lipids (which have segid MEMB since they are part of CHARMM's native lipid library).
 <code>
- 'counts': {
+config_dict = {
-        'ctype': 0, 'gtype': 3, 'otype': 1, 'stype': 4, 'l4type': 0, 'l5type': 0, 'l6type': 0, 'l7type': 0
+    'MEMB':['POPE','POPG',],
-    },
+}
+# ignore the current config_dict in the files with lipids like 'AIPE' and 'LLL'; that's specific to the membranes I was working with at the time when developing the scripts.
 </code>
-**(ii) Edit the resnames dictionary:** Specify the resnames for each lipid per lipid type.
-Example:
+====Step 2: run scripts====
-<code>
+First, run ''0_setup.py''. If you see errors, read the terminal output and fix ''config.py'' before rerunning.
-   'resnames': {
-        'rng': ['AIPE', 'AMPG', 'DPPC'],
-        'rno': ['AIOL'],
-        'rns': ['CER160', 'CER221', 'LSM', 'BSM'],
-    },
-</code>
- **Note:** The order of the actual resnames in the lists doesn't matter, but the lipid type/resname prefixes serving as dictionary keys must follow the order listed in the table from top to bottom.
+Next, if you left "recenter_dir = None" indicating you need to recenter your dcds, then ''recenter.scr'' will be created and you should run it using ''./recenter.scr''. This should take < 5 minutes. Otherwise, skip to the next step.
-**(iii) Edit the mol_counts dictionary:** Specify how many of each lipid is in the configuration (per top or bottom leaflet since each one gets its own configuration). This should mirror the order the resnames were defined in step (ii).
+Next, run ''./1_comp_sa.scr''. This sets up a virtual Python environment and installs the MDAnalysis package before submitting the Slurm jobs. This entire process can take 1-10 minutes.
-Example:
+Finally, run ''./2_analysis.py''. This computes statistics and writes txt and csv results to the ''results'' directory, which was created automatically. Here's some sample data:
 <code>
-'mol_counts': {
+leaflet resn  N_lip  Avg_SA  Std_SA  Stderr_SA
-        'rng': [35, 27, 51],
+    Top AIPE    113 58.3686  0.4927     0.1102
-        'rno': [11],
+    Top AMPG     25 61.0204  1.3520     0.3023
-        'rns': [18, 6, 10, 12],
+    Top  LLL     62 52.8182  0.5267     0.1178
-    }
+    Bot AIPE    105 58.2105  0.4785     0.1070
+    Bot AMPG     37 59.9090  0.8467     0.1893
+    Bot  LLL     58 52.8821  0.4872     0.1089
+N_blocks: 20
+N_comp_top: 3
+N_lip_top: 200
+N_comp_bot: 3
+N_lip_bot: 200
+Comp_SA_sum_top: 11395.8933
+Comp_SA_sum_bot: 11395.8933
+Total_bilayer_SA: 22791.7866
+Weighted average across leaflets (by N_lip):
+resn  N_lip_total  Avg_SA  Stderr_SA
+AIPE          218 58.2925     0.0769
+AMPG           62 60.3571     0.1662
+ LLL          120 52.8491     0.0805
 </code>
-**(iv) Finalize system setup:** Repeat steps 1-3 for the bottom leaflet and run ''./system_setup.scr''. This will display in the terminal the specified leaflets for your system and should look something like this:
+====Additional Info====
+At a high level, here's what the scripts are doing. Also, I highly recommend you watch a video on Voronoi tessellation so you can visually understand how the areas are computed.
-<code>
+**config.py**
-TOP leaflet (NTOT=336):
+System definition: segid/resname mapping, dcd range, file paths, and ''marker_map'' (atom selections per lipid class for Voronoi markers).
-  AIPE: 7 lipids
-  LIPA4: 21 lipids
-  LIPA5: 21 lipids
-  LIPA6: 21 lipids
-BOT leaflet (NTOT=510):
+**0_setup.py**
-  AIPE: 35 lipids
+Reads config.py, validates required CHARMM/dyn files, resolves marker atoms per lipid from marker_map, splits the dcd range into blocks, writes system_info.txt. Generates comp_sa.csh from src/tmpls/comp_sa.csh.tmpl.
-  AMPG: 27 lipids
-  DPPC: 51 lipids
+**system_info.txt**
-  LLL: 11 lipids
+Generated configuration file consumed by all downstream steps: dcd range, marker dictionary, block assignments, paths.
-  CER160: 18 lipids
-  CER221: 6 lipids
+**CHARMM recentering**
-  LSM: 10 lipids
+Recenters raw dcds so the bilayer midplane sits at z=0; output written to recenter/. Runs outside this pipeline; all downstream steps assume it has already been done.
-  BSM: 12 lipids
-</code>
+**1_comp_sa.scr**
+SLURM submission wrapper. Reads block boundaries from system_info.txt, submits one sbatch job per block via comp_sa.csh.
+**comp_sa.csh**
+SLURM job script. Sets up the run environment (venv/module) and executes src/comp_sa.py for its assigned dcd range.
+**src/comp_sa.py**
+Core computation. Per dcd: loads marker atoms, splits into top/bot leaflets, runs periodic 2D Voronoi tessellation, sums polygon areas per lipid, writes raw_data/raw_voronoi_{dcd}.dat (per-frame, per-leaflet, per-component SA/lipid).
-Make adjustments as necessary. **NTOT** is auto-calculated here and can be directly used in the following steps.
+**2_analysis.py**
+Final averaging step. Reads all raw_voronoi_*.dat files, block-averages each dcd, computes mean/std/stderr across blocks per leaflet+component and a leaflet-weighted whole-bilayer average per component. Writes results/comp_sa_results.txt and comp_sa_results.csv, comp_sa_weighted_avg.csv, comp_sa_totals_unweighted.csv, comp_sa_totals_weighted.csv.
-** 2. Run Symmetric Workflow Steps **
-Do steps 2, 4, 7, 8, and 10 from the workflow for symmetric membranes. Everything else can be skipped (handled automatically).
-** 3. Execute Scripts **
-Run script by command: ''./1_test.scr'' as mentioned above. Once finished, then run ''./2_run-dist-avg.scr''.
-**Note:** ''./2_run-nochol-dist-avg.scr'' is no longer used.