Differences

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

--- comp-sa [2021/06/09 08:44] – admin
+++ comp-sa [2026/09/03 09:56] (current) – [Step 1: edit config.py] edit
@@ Line 5: / Line 5: @@
 Instruction to obtain surface area per lipid of **__multicomponent system__** \\
 Original Authors: Xiaohong Zhuang and Viviana Monje\\
-Updated Authors: Yalun Yu and Jeffery Klauda
+Updated Authors: Yalun Yu and Jeffery Klauda\\
+Updated rtfpsf folder: Joshua Lucker \\
+Completely refactored, supports any membrane: Eric Zhang
-The script is in DT2 path: **/lustre/jbklauda/scripts/area_mult_lipids** \\
+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).
+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 |}}
-**How to Edit and Run the Scripts:**
+Modern Python scripts can be found here: {{ area_asym_lipids.tar.gz |}}
+**How to Edit and Run the Scripts (Legacy Instructions, Discouraged):**
 **Edit the scripts:** \\
@@ Line 112: / Line 118: @@
 **>> area-a**:  area-a is input file for top leaflet, area-b is for bottom leaflet
-**>> 2**: The lipid types of 2 are used instead of 8. (1-atom represented choline/sterol, and other 3-atoms represented lipids, e.g. glycerol lipids, sphingo lipids). So if you have sterols set this two 2 , if none set to 1.
+**>> 2**: This is the dimension of the tessellation, which is always 2.
 **>> 2196**: "image atoms" = 9*(#chol+3*#lip)=9*(28+3*72)
@@ Line 123: / Line 129: @@
 </code>
-**9. calc_avg_vertical.py & calc_avg_final.py:** \\
+**9a. If you have cholesterol: ** calc_avg_vertical.py & calc_avg_final.py: \\
-Update the number of each lipid component per leaflet in the same order as in lipid.dat but remove the sterol flag.
+Update the number of each lipid component per leaflet in the same order as in lipid.dat.
 <code>
 # Update the number of each lipid component per leaflet in the same order as in area-a.inp
 Nlist = [28,12,13,9,5,22,7,4]
 </code>
+**9b. If you lack cholesterol: ** calc_avg_vertical_nochol.py & calc_avg_final_nochol.py:
+Update the number of each lipid component per leaflet in the same order as in lipid.dat and no zeros for cholesterol.
+**10. rtfpsf.str:** \\
+The file names in the rtfpsf.str file will most likely need to be updated.
 **Run the scripts:** \\
@@ Line 138: / Line 151: @@
-.	Get the averages and standard errors by command: ''./2_run-dist-avg.scr''. You will see 2_avg_std_final.txt which contains the component area for each lipid type following the order in def_nlip.str.
+a.	Get the averages and standard errors if you have sterols by command: ''./2_run-dist-avg.scr''. You will see 2_avg_std_final.txt which contains the component area for each lipid type following the order in def_nlip.str.
+b.     Get the averages and standard errors if you DO NOT have sterols by command: ''./2_run-nochol-dist-avg.scr''. You will see 2_avg_std_final.txt which contains the component area for each lipid type following the order in def_nlip.str.
 ** Ignore the rest of this wiki for now...Yalun will updated as needed **
@@ Line 150: / Line 166: @@
       Simulations of multiple replicates (e.g. hypocotyl and root, each has 3 replicates) are performed, the REAL final average are taken by the average value of all the replicates. And the standard error is calculated by standard deviation of the **average SA/lipid** divided by square root of the replicate number. E.g. for hypocotyl system of 3 replicate, it is divided by sqrt(3).
 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=====
+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.