Differences

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

--- comp-sa [2026/08/23 17:23] – 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 \\
-Completely refactored, supports any membrane: Eric Zhang (IN PROGRESS)
+Completely refactored, supports any membrane: Eric Zhang
 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).
@@ 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.
-** Modern Python Workflow (IN PROGRESS) **\\
+=====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.
+Note: anytime you run into "no permission" issues, just run ''chmod u+x *''
+====Step 1: edit config.py====
+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.
+I embedded a lot of error handling in the scripts, so don't worry too much about precision if something isn't clear.
+Example Membrane with just POPE/POPG lipids (which have segid MEMB since they are part of CHARMM's native lipid library).
+<code>
+config_dict = {
+    '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>
+====Step 2: run scripts====
+First, run ''0_setup.py''. If you see errors, read the terminal output and fix ''config.py'' before rerunning.
+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.
+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.
+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>
+leaflet resn  N_lip  Avg_SA  Std_SA  Stderr_SA
+    Top AIPE    113 58.3686  0.4927     0.1102
+    Top AMPG     25 61.0204  1.3520     0.3023
+    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>
+====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.
+**config.py**
+System definition: segid/resname mapping, dcd range, file paths, and ''marker_map'' (atom selections per lipid class for Voronoi markers).
+**0_setup.py**
+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.
+**system_info.txt**
+Generated configuration file consumed by all downstream steps: dcd range, marker dictionary, block assignments, paths.
+**CHARMM recentering**
+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.
+**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).
+**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.