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_mult_lip.tar.gz

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.

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

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

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

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 check out DEVELOPERS_NOTES.md for a very detailed explanation of the Python pipeline (viewable via any Markdown extension or simply convert it to html/pdf). Any questions can be directed to ezhang32@umd.edu. As a heads up, the terms “lipid type” and “lipid headgroup” will be used interchangeably in the instructions below.

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. The top_counts and bot_counts subdicts use the same keys as the resnames subdicts, but contain non-negative integers instead of resnames to specify the quantity of each lipid in the top and bottom leaflets of your membrane. 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)
config_dict = {
    'resnames': {
        'PE': ['POPE','AIPE'],
        'PG': ['AMPG'],
        'NAAA':[('ORN', 'LLL')],
        'LPSA':[('LIPA4',)],
        'LPSB':[('LIPA5',)],
        'LPSC':[('LIPA6',)],
    },
    'top_counts': {
        'PE': [0, 0],
        'PG': [0],
        'NAAA': [0],
        'LPSA':[10],
        'LPSB':[10],
        'LPSC':[10],
    },
    'bot_counts': {
        'PE': [30, 30],
        'PG': [20],
        'NAAA': [20],
        'LPSA':[0],
        'LPSB':[0],
        'LPSC':[0],
    },
}

1.3 Stay in config.py and specify the dcd_range, charmm_dir, dyn_dir, asymmetric, num_blocks, output_csv, and compress_data fields as outlined in the comments.

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.

2.2 Run ./1_calc_hb_run.scr to submit your slurm scripts. They can take from 1-6 hours depending on the size of your system. By default, the Slurm jobs are set with an 8-hour time limit, which should be enough, but feel free to manually edit 1_calc_hb_run.csh which is located in your freshly generated directories, to change Slurm specifics such as increasing the job time or running on the scavenger queue instead of the standard queue.

2.3 Once all Slurm jobs finish running, run 2_analysis.scr. This merges data, runs statistics, and compresses data if specified. If you navigate to the *-hbonds directories, you should *.txt files with your hbond data, heatmaps for all and pairwise lipid h-bond interactions, and CSV files if specified earlier in config.py. If you want to see which donors and acceptors are bonding the most, look at the pair number for all-hbonds or the DN-AC pairing in all/inter-hbond results and determine the corresponding group based on your def_dn_ac.str file. It's worth noting that whenever it says normalized, it means the raw h-bond data was divided by the donor lipid count in your system. Here's some sample asymmetric inter-hbond output. The results for a symmetric bilayer are similar (just averaged across leaflets):

--- Top Leaflet Aggregate Pairwise Inter-lipid Hbond Interactions ---
DLip     ALip     N_avg    N_se     N_avg/lip N_se/lip
AIPE     AIPE     40.5140  0.6165   0.3585    0.0055  
AIPE     AMPG     9.7620   0.2851   0.0864    0.0025  
AIPE     ORN      34.9095  0.5685   0.3089    0.0050  
AMPG     AIPE     4.5250   0.1493   0.1810    0.0060  
AMPG     AMPG     1.2820   0.0803   0.0513    0.0032  
AMPG     ORN      1.8625   0.1378   0.0745    0.0055  
ORN      AIPE     34.3005  0.5171   0.5532    0.0083  
ORN      AMPG     9.6245   0.2501   0.1552    0.0040  
ORN      ORN      25.2025  0.5834   0.4065    0.0094  

--- Bottom Leaflet Aggregate Pairwise Inter-lipid Hbond Interactions ---
DLip     ALip     N_avg    N_se     N_avg/lip N_se/lip
AIPE     AIPE     35.2030  0.4736   0.3353    0.0045  
AIPE     AMPG     12.4245  0.1884   0.1183    0.0018  
AIPE     ORN      30.8595  0.6351   0.2939    0.0060  
AMPG     AIPE     6.8320   0.1324   0.1846    0.0036  
AMPG     AMPG     2.3940   0.1175   0.0647    0.0032  
AMPG     ORN      2.7050   0.1255   0.0731    0.0034  
ORN      AIPE     30.5265  0.7185   0.5263    0.0124  
ORN      AMPG     13.9795  0.3651   0.2410    0.0063  
ORN      ORN      20.5845  0.4912   0.3549    0.0085  

--- Top Leaflet Normalized Inter-lipid Hbond Interactions ---
Lip      Nlip   DN_avg/lip DN_se/lip AC_avg/lip AC_se/lip Tot_avg/lip Tot_se/lip
AIPE     113    0.7539     0.0057    0.3436     0.0042    1.0974      0.0073    
AMPG     25     0.3068     0.0093    0.7755     0.0127    1.0822      0.0176    
ORN      62     1.1150     0.0112    0.5931     0.0089    1.7081      0.0137    

--- Top Leaflet Raw Inter-lipid Hbond Interactions ---
Lip      Nlip   DN_avg   DN_se    AC_avg   AC_se    Tot_avg  Tot_se  
AIPE     113    85.1855  0.6454   38.8255  0.4782   124.0110 0.8248  
AMPG     25     7.6695   0.2329   19.3865  0.3178   27.0560  0.4409  
ORN      62     69.1275  0.6927   36.7720  0.5525   105.8995 0.8501  

--- Bottom Leaflet Normalized Inter-lipid Hbond Interactions ---
Lip      Nlip   DN_avg/lip DN_se/lip AC_avg/lip AC_se/lip Tot_avg/lip Tot_se/lip
AIPE     105    0.7475     0.0064    0.3558     0.0067    1.1033      0.0108    
AMPG     37     0.3225     0.0041    0.7136     0.0096    1.0361      0.0116    
ORN      58     1.1222     0.0068    0.5787     0.0109    1.7009      0.0137    

--- Bottom Leaflet Raw Inter-lipid Hbond Interactions ---
Lip      Nlip   DN_avg   DN_se    AC_avg   AC_se    Tot_avg  Tot_se  
AIPE     105    78.4870  0.6699   37.3585  0.6986   115.8455 1.1339  
AMPG     37     11.9310  0.1530   26.4040  0.3535   38.3350  0.4298  
ORN      58     65.0905  0.3944   33.5645  0.6329   98.6550  0.7965