Brief instruction to obtain H-bond by CHARMM script for a Single Lipid Bilayer
Authors: Xiaohong Zhuang and Dr. Jeffery B. Klauda
This instruction is to explain the procedure to calculate the number of H-bonding in single lipid bilayer including the PBC.
! 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.
Script in ZT1: /afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/hb_chm_1_lip
Example in ZT1: /afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/xiao/lipids_analyses/hb_chm_1_lip_eg
Download Scripts: hb_chm_1_lip.tar.gz
Steps
1. Copy the script folder provided to your path.
1.1 Go to your path that you would like to work on H-bond analysis
1.2 For example, use the following command to copy the script folder to your path
cp -r /lustre/xzhuang/simulation/script_folders/lipids_analyses/hb_chm_1_lip ./
2. Updates to the 1_calc_hb_run.csh file
2.1 Update the path of your input file
set sdir = “/lustre/xzhuang/simulation/dmpi/charmm-gui”
2.2 Update equilibrated dcd file range
@ f = 36 # first dcd file to read @ l = 50 # last dcd file to read
2.3 Update lipid information including lipid type, lipid name, and headgroup name
2.3.1 Lipid type [type] for Glycerol lipids:G ; Sphingo lipids:S ;
2.3.2 Update lipid residue name [rn]
2.3.3. Update headgroups [h] for Glycerol:PA/PC/PE/PG/PI/PS; Sphingo: SM/CER ;
DON'T includes the sterol,cholestrol lipids here as it has no H acceptor. DON'T includes the PC glycerol lipids here as it has no H donor.
Example for a DMPI membrane
# Update the lipid that have both H donor and acceptor for intra-lipid H-bonding. set type = G set rn = DMPI set h = PI
3. Updater 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.
For example change each short path in your toppar.str from the following
open read card unit 10 name toppar/top_all36_prot.rtf read rtf card unit 10
To the full path as following
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.
4. In 1_calc_hb_all.inp and 1_calc_hb_intra.inp, you may need to update the following when necessary
4.1 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
open read unit @u file name @P/namd/dyn/dyn@F.dcd
4.2 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).
5. Submit the calculation job by 1_calc_hb_run.scr.
6. After the run completed, you will see two output files, hb_all.txt and hb_intra.txt, which shows the number of all H-bonds and 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. If you are interested in detail of each pair, you may refer to the files def_don_acc.str, def_hb_all.str, and def_hb_intra.str to see the atoms in each pair of your lipid.
7. Calculate the inter-lipid Hbond by command 2_calc_inter_hb.scr
Ninter_lipid = Ninter_plus_intra - Nintra_lipid
The file hb_inter.txt will be obtained.
8. Calculate the statistics of the Hbond by command 3_calc_avg_percent.scr
The results will be shown in two files. avg-perc-intra.txt and avg-perc-inter.txt.
8.1 As shown in example output below, the 1st column is the H donor and acceptor pair group number (the final row is the Hbond of all the pairs).
8.2 The 2nd and 3rd 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.
8.3 The 4th and 5th columns are the block averaged number of Hbond per lipid and standard error of it.
avg-perc-intra.txt:
Grp N_avg N_se N/lip Nse/lip 0 15.4661 0.5143 0.2148 0.0071 1 2.4287 0.1527 0.0337 0.0021 2 0.0000 0.0000 0.0000 0.0000 3 0.0000 0.0000 0.0000 0.0000 Tot 17.8948 0.5387 0.2485 0.0075
9. It is recommend to run 4-gzip-output.scr to gzip the output file.