Differences

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

--- hydrogen_bonds_for_multi-component_bilayers [2020/06/10 12:02] – [Section 2 (All-lipid Hbonds)] edit
+++ hydrogen_bonds_for_multi-component_bilayers [2026/10/11 23:01] (current) – Large update to the modern workflow edit
@@ Line 1: / Line 1: @@
 [[Analysis]]
-=====Brief instruction to obtain H-bond by CHARMM script for a Single Lipid Bilayer=====
+=====Brief instructions to obtain H-bond by CHARMM script for a Multi Lipid Bilayer=====
-Authors: Xiaohong Zhuang and Dr. Jeffery B. Klauda
+Original Authors: Xiaohong Zhuang and Dr. Jeffery B. Klauda\\
+Refactored by Eric Zhang (2026) to handle any membrane system\\
-This instruction is explain the procedure to calculate the number of H-bonding of a multiple lipid bilayer. For the H-bond of single lipid system, please use the instruction: [[hydrogen_bonds_for_single-component_bilayers|]].
+This instruction explains the procedure for calculating the number of H-bonds in a multiple-lipid bilayer. For the H-bond of a single lipid system, you may use the instruction: [[hydrogen_bonds_for_single-component_bilayers|]]. However, it's highly recommended that users use the modern Python package, as it can handle any symmetric/asymmetric, lipid/protein membrane system while providing more detailed statistics. For now, the instructions to run the old, legacy scripts are left in place. To see instructions for the Python package, scroll down to the section titled modern workflow after downloading the package.
-! Special Note: the CHARMM version to run Xiao's H-bond script should be newer or equal to c41b2, versions before this may fail to produce the correct number of hydrogen bonds.
+! Special Note: the CHARMM version used to run all H-bond scripts should be at least c41b2; earlier versions may fail to produce the correct number of hydrogen bonds.
-Script in DT2: **/lustre/jbklauda/scripts/xiao/lipids_analyses/hb_chm_mult_lip**\\
+Legacy Scripts (Only designed for symmetric membranes with glycerophospholipids and sphingolipids)\\
-Example in DT2: **/lustre/jbklauda/scripts/xiao/lipids_analyses/hb_chm_mult_lip_eg**
+Script in ZT1: **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/hb_chm_mult_lip**\\
+Example in ZT1: **/afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/hb_chm_mult_lip_eg**\\
+Download Script: {{ ::hb_mult.tar.gz |}}\\
-Download Script: {{ ::hb_mult.tar.gz |}}
+Python Package (Recommended)\\
+Download Script: {{ ::hb_asym_lipids_v2.tar.gz |}}
+======Legacy Workflow=====
 =====Initial Steps=====
@@ Line 20: / Line 26: @@
 <code commend>
-cp -r /lustre/jbklauda/scripts/xiao/lipids_analyses/hb_chm_mult_lip ./
+cp -r /afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/hb_chm_mult_lip ./
 </code>
@@ Line 119: / Line 125: @@
 <code commend>
-Lip	Nlip	N_avg	N_se	Na/lip	Nse/lip\\
+Lip	Nlip	N_avg	N_se	Na/lip	Nse/lip
-DOPE	12	0.7446	0.0223	0.0620	0.0019\\
+DOPE	12	0.7446	0.0223	0.0620	0.0019
-NSM	12	6.5965	0.0231	0.5497	0.0019\\
+NSM	12	6.5965	0.0231	0.5497	0.0019
 TotL	24
@@ Line 186: / Line 192: @@
 set rns1 NSM
 set hs1 SM
+</code>
 .3.3. After you have the full path, change the **toppar.str** path to your **toppar.str** path:
@@ Line 266: / Line 273: @@
 NSM	12	5.8363	0.3058	0.4864	0.0255
 </code>
+======Modern Python Workflow======
+This workflow is designed to automate the old workflow by generating optimized CHARMM files based on a few inputs describing your membrane system. This wiki covers what's necessary, but feel free to ask me anything at ezhang32@umd.edu. As a heads up, the terms "lipid type" and "lipid headgroup" will be used interchangeably in the instructions below.
+=====Step 1: Edit Neccesary Files=====
+.1 Add new h-bond definitions by opening ''hbond_defs.py''. Read through the comments to see how hydrogen bond donors and acceptors are defined for standard and patched lipids. In short, we use nested dictionaries to map each lipid type/headgroup to the atoms that make up its donors and acceptors. By default, definitions are already provided for sterols, glycerophospholipids, and sphingolipids. If you're only using these lipid types, you can skip this step entirely.
+Definitions for typical 4/5/6-tailed Lipid A are also included as patched lipids. To use these, update the headgroup key (the segid base in your PSF, without the trailing number) and the resname key nested underneath it (matching the resname in your PSF/toppar files).
+To add h-bond definitions for a new lipid type, add an entry to ''standard_dict'' (if it's uniquely defined by resid) or ''patched_dict'' (if it spans multiple residues and is uniquely defined by segid), mirroring the structure of the existing entries. Bottom line: if your lipid is uniquely identified by its resid, then use the standard_dict. If your psf shows the lipid being uniquely identified by the segid, use the patched_dict. An example is shown below.
+<code>
+# adding a standard lipid (or protein) type called CUSTOM that has an OH and C=O group
+standard_dict = {
+    ....
+    'CUSTOM': {
+        'donor': [('H1',)],
+        'acceptor': [('O1',), ('OC')], #OC is carbonyl oxygen
+    },
+}
+# adding a patched lipid (or protein) with segids CUSTOM* and residues RESI1 (contains OH) and RESI2 (contains C=O) under each segid
+standard_dict = {
+    ....
+    'CUSTOM': {
+        'donor': {
+            RESI1: [('H1',)],
+            RESI2: [],
+        },
+        'acceptor': {
+            RESI1: [('O1',),],
+            RESI2: [('OC',),],
+        },
+    },
+}
+</code>
+.2 Open ''config.py'' to define your system. First, navigate to standard_dict. The resnames subdictionary (subdict) holds lipid types as keys and lists of corresponding resnames as values. Note: for patched lipid types, the assigned value can only be a list with a single tuple entry that holds all the patched/conjugated residues for that segid. You no longer need to specify the lipid counts per leaflet, as this is checked at runtime. Below is an example:
+<code>
+# top leaflet: 10 4-chained lipid A, 10 5-chained lipid A, 10 6-chained lipid A
+# bottom leaflet: 30 POPE, 30 AIPE, 20 POPG, 20 NAAA (patched ornithine lipid)
+resnames = {
+    'PE': ['AIPE'],
+    'PG': ['AMPG'],
+    'NAAA': [('ORN', 'LLL')],
+    'LPSA': [('LIPA4',)],
+    'LPSB': [('LIPA5',)],
+    'LPSC': [('LIPA6',)],
+}
+</code>
+.3 Stay in ''config.py'' and specify the dcd_range, charmm_dir, dyn_dir, recenter_dir, asymmetric, num_blocks, output_csv, and other fields as outlined in the comments.
+=====Step 2: Run Scripts=====
+.1 Run ''./0_setup.scr''. This generates the directories that store your data and CHARMM scripts. I've also coded lots of error handling and validation functions, so it should catch most of the errors in config.py and hbond_defs.py if they exist. Simply re-run the script once you fix all the errors, if any occur. If you set recenter_dir = None, then it should generate ''./recenter.scr''. Run it, and wait for it to finish before continuing.
+.2 Run ''./1_hbonds.scr'' to submit a slurm job. This used to take up to 6-12 hours, but has been optimized to run in under a minute!
+.3 Once it finishes, run ''2_analysis.scr''. This merges data, runs statistics, and plots data if specified. If you navigate to the results directory, you should now see *-hbonds directories; they should contain *.txt files with your hbond data, heatmaps for all and pairwise lipid h-bond interactions, and CSV files if specified earlier in config.py. There's also an hbond_key.txt file which shows which donor/acceptor symbols in the tables correspond to which atoms. It's worth noting that whenever it says normalized, it means the raw h-bond data was divided by the donor lipid count. The pairwise affinity tables/heatmaps provides a better metric since its normalized against both donors and acceptors.