[[Analysis]]
=====Brief instructions to obtain H-bond by CHARMM script for a Multi Lipid Bilayer=====
Original Authors: Xiaohong Zhuang and Dr. Jeffery B. Klauda\\
Refactored by Eric Zhang (2026) to handle any membrane system\\
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 used to run all H-bond scripts should be at least c41b2; earlier versions may fail to produce the correct number of hydrogen bonds.
Legacy Scripts (Only designed for symmetric membranes with glycerophospholipids and sphingolipids)\\
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 |}}\\
Python Package (Recommended)\\
Download Script: {{ ::hb_asym_lipids_v2.tar.gz |}}
======Legacy Workflow=====
=====Initial Steps=====
1. Copy the script folder provided to your path.\\
2. Go to your path that you would like to work on H-bond analysis.\\
3. Use the following command to copy the script folder to your path
cp -r /afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/hb_chm_mult_lip ./
For multiple lipids system, the intralipid Hbonds and all Hbonds are calculated in section 1 and 2 separately and the interlipid Hbonds are calculated in Section 3. The example system provided has 1 sterol/cholesterol lipid: CHL1, 2 glycerol lipids: DPPC, DOPE, and 1 sphingo lipid: NSM.
=====Section 1 (Intralipid Hbonds)=====
1.1 Edits to: **1_calc_hb_run.csh**
1.1.1 Update the path of the file based on your system\\
''set sdir = "/lustre/mda003/hlms_310.15/pe_pc_nsm_chol/a_traj"''
1.1.2 Update the equilibrated dcd file range based on your system
@ f = 101 # first dcd file to read
@ l = 150 # last dcd file to read
1.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
# 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 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
1.2 Update your **toppar.str** file
Update the full path in your toppar.str so that you don’t need to copy the toppar folder and it allows to call your topology files correctly, i.e. change each short path in your **toppar.str** from
open read card unit 10 name toppar/top_all36_prot.rtf
read rtf card unit 10
To the full path as following but with your path
open read card unit 10 name /lustre/xzhuang/simulation/toppar/top_all36_prot.rtf
read rtf card unit 10
If you have not done this, it is recommended that you can make one **toppar.str** and toppar folder in your path for all your analyses and update the topology files in that particular toppar folder when necessary, so that you don’t need to update the path.
1.3 In **1_calc_hb_intra.inp**, you may need to update the following when necessary
1.3.1. After you have the full path, change the **toppar.str** path to your **toppar.str** path.
''stream "/lustre/xzhuang/simulation/toppar.str"''
1.3.2 The path of trajectory file is based on common CHARMM-GUI package setup, if your trajectory files are in a different location, please update the path
''open read unit @u file name @P/namd/dyn/dyn@F.dcd''
1.3.3 Update data point per dyn file when necessary
''calc ntot = 1000 * @N''
Number of data points in the DYN files = Total time steps divided by dcdfreqency. It can be found by reading your **dyn-2.inp** file to find dcd frequency and total time steps. It is usually 1000 (total times steps 1000000 per dyn / dcd frequency 1000).
1.4 Submit the calculation job by ''1_calc_hb_run.scr''.
After the run is completed, you will see multiple *.txt, which shows the number of intra-lipid H-bonds of each H donor-Acceptor pair of the lipid of each frame. The last column is the total H-bond of each frame.
1.5 Calculate the statistics of the Hbond by command ''2_calc_avg_percent.scr''
The results will be shown in two files: **avg-intra.txt** and **avg-perc-intra.txt**.
In **avg-intra.txt**, as shown in example output below, the 1st columns is lipid name, the 2nd column is number of each lipid, the 3rd column is the H donor and acceptor pair group number (the final row is the Hbond of all the pairs). The 4th and 5th columns show the block averaged number of Hbond and standard error of the block average, respectively. The block size used here are each dyn. You may update the **avg_percent.py** if you want to use different block size. The 6th and 7th columns are the block averaged number of Hbond per lipid and standard error of it.
If you are interested in detail of each pair, you may refer to the files **def_dn_ac.str**, and **def_hb_intra.str** to see the atoms in each pair of your lipid.
Lip Nlip Pair N_avg N_se Na/lip Nse/lip
DOPE 12 1 0.5637 0.0165 0.0470 0.0014
DOPE 12 2 0.1809 0.0169 0.0151 0.0014
DOPE 12 3 0.0000 0.0000 0.0000 0.0000
DOPE 12 4 0.0000 0.0000 0.0000 0.0000
DOPE 12 5 0.0000 0.0000 0.0000 0.0000
DOPE 12 6 0.0000 0.0000 0.0000 0.0000
DOPE 12 TotG 0.7446 0.0223 0.0620 0.0019
…
In **avg-perc-intra.txt**, as shown below, the total donor-acceptor pairs of the block averaged number of Hbond and standard error of the block average are shown in 3rd and 4th columns, respectively. The 6th and 7th columns are the block averaged number of Hbond per lipid and standard error of it.
Lip Nlip N_avg N_se Na/lip Nse/lip
DOPE 12 0.7446 0.0223 0.0620 0.0019
NSM 12 6.5965 0.0231 0.5497 0.0019
TotL 24
=====Section 2 (All-lipid Hbonds)=====
2.1 Edits in **1_calc_hb_run.csh**
2.1.1 Update the path of the file:
''set sdir = "/lustre/mda003/hlms_310.15/pe_pc_nsm_chol/a_traj"''
2.1.2 Update the equilibrated dcd file range
@ f = 101 # first dcd file to read
@ l = 150 # last dcd file to read
2.1.3 Update lipid type, residue name and headgroup type (//It is different from intra lipid HB as all the H donor lipid, including sterol lipids should be listed here//.)
# 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 ;
# 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.
# 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
# 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=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
2.2 Update your **toppar.str** file (Please skip to 2.3 if you already done this in Section 1.2). See details of this in Section 1.2
2.3 Edits for **1_calc_hb_all.inp**
2.3.1 Update total number of H acceptors [includes glycerol and shingo lipid but not Sterol]
set gtyp 2
set styp 1
2.3.2 List glycerols and sphingo H acceptor lipids below [No Sterol lipids, as it has no H aceptor], should be consistent as above
! 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
set rng1 DPPC
set hg1 PC
set iscrb1 Y
set rng2 DOPE
set hg2 PE
set iscrb2 Y
set rns1 NSM
set hs1 SM
2.3.3. After you have the full path, change the **toppar.str** path to your **toppar.str** path:
''stream "/lustre/xzhuang/simulation/toppar.str"''
2.3.4 The path of trajectory file is based on common CHARMM-GUI package setup, if your trajectory files are in different location, please update the path of it:
''open read unit @u file name @P/namd/dyn/dyn@F.dcd''
2.3.5 Update data point per dyn file when necessary
''calc ntot = 1000 * @N''
Number of data points in the DYN files = Total time steps divided by dcdfreqency. It can be found by reading your **dyn-2.inp** file to find dcd frequency and total time steps. It is usually 1000 (total times steps 1000000 per dyn / dcd frequency 1000).
2.4 Submit the calculation job by ''1_calc_hb_run.scr''
After the run completed, you will see multiple **Dlip1-Alip.txt**, which shows the number of all H-bonds of each H donor-Acceptor pair of the lipid of each frame. The last column is the total H-bond of each frame.
2.5 Update **2_calc_avg_percent.scr**
# Update H donor [Dlip] and acceptor lipid [Alip] name
# Sterol lipids are H donor only, not H acceptor, so it is in Dlip list , not in Alip
# PC glycerol lipids are not H donor, but a H acceptor only, so it is not in Dlip, but in Alip list.
foreach Dlip ( CHL1 DOPE NSM )
foreach Alip ( DPPC DOPE NSM )
echo $Dlip-$Alip >> lippair.txt
end
end
2.6 Calculate the statistics of the Hbond by command ''2_calc_avg_percent.scr''
The results will be shown in two files: **avg-all.txt** and **avg-perc-all.txt**.
In **avg-all.txt**, as shown in example output below, the 1st columns is Donor lipid – Acceptor lipid names, the 2nd column is the H donor and acceptor pair group number (the final row is the Hbond of all the pairs). The 3rd and 4th columns show the block averaged number of Hbond and standard error of the block average, respectively. The block size used here are each dyn. You may update the avg_percent.py if you want to use different block size.
If you are interested in details of each pair, you may refer to the files **def-donor.str**, **def-acceptor-g.str**, **def-acceptor-s.str**, and **def_hb_all.str** to see the atoms in each pair of your lipid.
Dlip-Alip pair N_avg N_se
chl1-dppc 1 2.0728 0.1151
chl1-dppc 2 1.4138 0.0714
chl1-dppc 3 0.0000 0.0000
chl1-dppc TotG 3.4866 0.1316
…
In **avg-perc-all.txt**, as shown below, the total lipid number of lipid is shown in 2nd column, the total donor-acceptor pairs of the block averaged number of Hbond and standard error of the block average are shown in 3rd and 4th columns, respectively. The 6th and 7th columns are the block averaged number of Hbond per lipid and standard error of it. The total number of lipids are shown in the last row.
Lip Nlip N_avg N_se Na/lip Nse/lip
CHL1 64 9.6250 0.3658 0.1504 0.0057
DOPE 12 7.8334 0.3803 0.6528 0.0317
DPPC 12 5.5743 0.2490 0.4645 0.0208
NSM 12 12.4328 0.3289 1.0361 0.0274
TotL 100
=====Section 3 (Inter-Lipid HBond)=====
''Ninter_lipid = Ninter_plus_intra - Nintra_lipid''
After finish the 1-intra and 2-all calculation (complete all the steps in them including the average calculation), go to 3-inter folder, calculate inter-lipid Hbond by command ''3_run_calc_inter.scr''.
The file **avg-perc-inter.txt** will be obtained:
Lip Nlip N_avg N_se Na/lip Nse/lip
CHL1 64 9.6250 0.3658 0.1504 0.0057
DPPC 12 5.5743 0.2490 0.4645 0.0208
DOPE 12 7.0888 0.3580 0.5907 0.0298
NSM 12 5.8363 0.3058 0.4864 0.0255
======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.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.
# 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',),],
},
},
}
1.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:
# 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',)],
}
1.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=====
2.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.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!
2.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.