Differences

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

--- hydrogen_bonds_for_multi-component_bilayers [2017/02/10 14:36] – admin
+++ 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.
-Script in DT1:  **/lustre/xzhuang/simulation/script_folders/lipids_analyses/hb_chm_mult_lip**\\
+! 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.
-Example in DT1: **/lustre/xzhuang/simulation/script_folders/lipids_analyses/hb_chm_mult_lip_eg**
-Script in DT2: **/lustre/xzhuang/scripts/lipids_analyses/hb_chm_mult_lip**\\
+Legacy Scripts (Only designed for symmetric membranes with glycerophospholipids and sphingolipids)\\
-Example in DT2: **/lustre/xzhuang/scripts/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:
+Python Package (Recommended)\\
+Download Script: {{ ::hb_asym_lipids_v2.tar.gz |}}
+======Legacy Workflow=====
 =====Initial Steps=====
@@ Line 21: / Line 26: @@
 <code commend>
-cp -r /lustre/xzhuang/simulation/script_folders/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 41: / Line 46: @@
 </code>
-.1.3 Update lipid information including lipid type, lipid name, and headgroup name
+.1.3 Update lipid information including lipid type, lipid name, and headgroup name, as well as whether glycerol lipids are lyso lipids and have carbonyl groups
 <code commend>
 # 1. Update lipid type [type] for Glycerol lipids:G ;                Sphingo lipids:S ;
 # 2. Update lipid resdue name [rn]
 # 3. Update headgroups [h] for Glycerol lipids:PA/PC/PE/PG/PI/PS;    Sphingo lipids: SM/CER ;
+# 4. For glycerol lipids, specify the following for the lipid tails:
+#      a) If there is a carbonyl moiety [iscrb] as Y (usually case), otherwise N (if one tail of each, use Y)
+#      b) If there is no tail (just an alcohol group) [islys] as Y, otherwise N (if one tail of each, use Y)
 # DON'T include the sterol,cholestrol lipids here as it has no H acceptor, so no intralipid HB.
 # DON'T include the PC glycerol lipids here as it has no H donor, so no intralipid HB.
+# Exception: lyso PC lipid should be included, as alcohol group provides H donor
 # List lipids in your system that have both H donor and acceptor for intra-lipid H-bonding
-$chm f=$f l=$l p=${sdir} type=G rn=DOPE h=PE < 1_calc_hb_intra.inp >& 1_calc_hb_intra-DOPE.out
+$chm f=$f l=$l p=${sdir} type=G rn=DOPE h=PE iscrb=Y islys=N < 1_calc_hb_intra.inp >& 1_calc_hb_intra-DOPE.out
-$chm f=$f l=$l p=${sdir} type=S  rn=NSM h=SM < 1_calc_hb_intra.inp >& 1_calc_hb_intra-NSM.out
+$chm f=$f l=$l p=${sdir} type=S rn=NSM  h=SM                 < 1_calc_hb_intra.inp >& 1_calc_hb_intra-NSM.out
 </code>
@@ Line 116: / 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 141: / Line 150: @@
 <code commend>
-# Lipid type for Glycerol lipids:G ;          Sphingo lipids:S ;       Sterol/Cholestrol lipids: C ;
+# Lipid type for Glycerol lipids:G ;                Sphingo lipids:S ;           Sterol/Cholestrol lipids: C ;
-# Headgroups for Glycerol lipids:PA/PC/PE/PG/PI/PS; Sphingo lipids: SM/CER ;
+# headgroups for Glycerol lipids:PA/PC/PE/PG/PI/PS; Sphingo lipids: SM/CER ;
-# DON'T put PC lipid here as it has no H donor, but included PC lipids in 1_calc_hb_all.inp
+# DON'T put PC lipid here as it has no H donor, but included PC lipids in 1_calc_hb_all.inp as it is H acceptor.
-    as it is H acceptor.
+# Exception: lyso PC lipid should be included, since alcohol group gives H donor
 # Update the lipids with H donors with lipid_type, lipid_name, and headgroup in each of the following lines
-# Do not need to specify h for Sterol lipids
+# For glycerol lipids: specify [islys], as with part 1 of the script (but [iscrb] not specified here)
+# Donot need to specify headgroup for Sterol lipids
-$chm f=$f l=$l p=${sdir} type=C rn=CHL1            < 1_calc_hb_all.inp >& 1_calc_hb_all_CHL1.out
+$chm f=$f l=$l p=${sdir} type=C rn=CHL1              < 1_calc_hb_all.inp >& 1_calc_hb_all_CHL1.out
-$chm f=$f l=$l p=${sdir} type=G rn=DOPE h=PE < 1_calc_hb_all.inp >& 1_calc_hb_all_DOPE.out
+$chm f=$f l=$l p=${sdir} type=G rn=DOPE h=PE islys=N < 1_calc_hb_all.inp >& 1_calc_hb_all_DOPE.out
-$chm f=$f l=$l p=${sdir} type=S  rn=NSM  h=SM < 1_calc_hb_all.inp >& 1_calc_hb_all_NSM.out
+$chm f=$f l=$l p=${sdir} type=S rn=NSM  h=SM         < 1_calc_hb_all.inp >& 1_calc_hb_all_NSM.out
 </code>
@@ Line 168: / Line 178: @@
 <code commend>
+! List glycerols and sphingo lipids below [No Sterol lipids, as it has no H aceptor]
 ! Use rng and hg for resname and headgroup for glycerol lipids, and rns and hs for sphingo lipids
+! Specify whether the glycerol lipids have a carbonyl moiety [iscrb] as Y (usually case), otherwise N
+! H acceptors, should be consistent as above
-! H acceptors
 set rng1 DPPC
-set hg1  PC
+set hg1 PC
+set iscrb1 Y
 set rng2 DOPE
-set hg2  PE
+set hg2 PE
+set iscrb2 Y
 set rns1 NSM
-set hs1  SM
+set hs1 SM
 </code>
@@ Line 258: / 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.