[[Analysis]]

**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

<code commend>
cp -r /lustre/xzhuang/simulation/script_folders/lipids_analyses/hb_chm_1_lip ./
</code>


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
<code commend>
@ f = 36       # first dcd file to read
@ l = 50       # last dcd file to read
</code>

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**

<code commend>
# Update the lipid that have both H donor and acceptor for intra-lipid H-bonding.  
set type = G
set rn   = DMPI
set h    = PI
</code>


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
<code commend>
open read card unit 10 name toppar/top_all36_prot.rtf
read  rtf card unit 10
</code>

To the full path as following
<code commend>
open read card unit 10 name /lustre/xzhuang/simulation/toppar/top_all36_prot.rtf
read  rtf card unit 10
</code>
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
<code commend>
open read unit @u file name @P/namd/dyn/dyn@F.dcd  
</code>
4.2 Update data point per dyn file when necessary
<code commend>
calc ntot = 1000 * @N
</code>
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//:
<code commend>
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
</code>
9.	It is recommend to run //4-gzip-output.scr// to gzip the output file.