Inter chain hydrogen bond calculation

User discussions
1

Hello everyone,
I need your help to calculate inter-chain hydrogen bonds.
I have three chains: A, B, and C. Their respective indices are as follows:

$comm make_ndx -f "$TPR_FILE" -o "$index_FILE"<< EOF
a 1-860
name 17 chain_a
a 861-1720 
name 18 chain_b
a 1721-2580 
name 19 chain_c
18 | 19
name 20 chain_bc
17 | 19
name 21 chain_ac
17 | 18
name 22 chain_ab

I am using the following code to generate the results.
Here, as for example, when chain_a (17) is a donar then chain_bc (20) is acceptor.Same for the others.

chains=("chaina" "chainb" "chainc")

# Iterate over each chain
for chain in "${chains[@]}"; do
    GRO_FILE="${BASE_PATH}/trimer_traj5_npt_55.tpr"
    XTC_FILE="${BASE_PATH}/trimer_traj5_55.xtc"
    index_file="${BASE_PATH}/analysis/hbond/index.ndx"
    hb_FILE="${BASE_PATH}/analysis/hbond/1_hbonds_${chain}.xvg"
    o_FILE="${BASE_PATH}analysis/hbond/hbond_${chain}.ndx"
    
    # Run GROMACS hbond command
    if [ "$chain" = "chaina" ]; then
        (echo 17; echo 20) | $comm hbond -f "$XTC_FILE" -s "$GRO_FILE" -n "$index_file" -num "$hb_FILE" -cutoff 0.35 -m 
    elif [ "$chain" = "chainb" ]; then
        (echo 18; echo 21) | $comm hbond -f "$XTC_FILE" -s "$GRO_FILE" -n "$index_file" -num "$hb_FILE" -cutoff 0.35 -m 
    elif [ "$chain" = "chainc" ]; then
        (echo 19; echo 22) | $comm hbond -f "$XTC_FILE" -s "$GRO_FILE" -n "$index_file" -num "$hb_FILE" -cutoff 0.35 -m 
    fi
done

However, when I add all three files together, it generates a higher number of hydrogen bonds than expected.
Could you please tell me the right way to calculate inter-chain hbond?
Thanks in advance!

2

Hi!

How many selections do you have in your index file? If you have three chains you should have three selections, one for each chain. You can of course define the chains manually with indices like you have done, but it’s probably easier (and safer) to use “splitch 1” to get the three protein chains (1 should represent the protein). Then, you can just pass this index file to gmx hbond, select the two chains between which you want to calculate the number of hbonds and then add the results together. Note that you can also use gmx select to create the index file, but splitch in make_ndx is probably the easiest option.

Let me know if this helped!

1 Like
4

Hello Cathrine,
Thank you so much. It works!

Below is my bash script:

BASE_PATH="/home/md_simulation"
comm=gmx

$comm make_ndx -f "$TPR_FILE" -o "$index_FILE"<< EOF
a 1-860
name 17 chain_a
a 861-1720 
name 18 chain_b
a 1721-2580 
name 19 chain_c
q
EOF
chains=("chaina" "chainb" "chainc")

# Create the analysis directory if it doesn't exist
mkdir -p "${BASE_PATH}"

# Iterate over each chain
for chain in "${chains[@]}"; do
    GRO_FILE="${BASE_PATH}/trimer_traj5_npt_55.tpr"
    XTC_FILE="${BASE_PATH}/trimer_traj5_55.xtc"
    index_file="${BASE_PATH}/analysis/hbond/index.ndx"
    hb_FILE="${BASE_PATH}/analysis/hbond/1_hbonds_${chain}.xvg"
    o_FILE="${BASE_PATH}/analysis/hbond/hbond_${chain}.ndx"
    
    # Run GROMACS hbond command
    if [ "$chain" = "chaina" ]; then
        $comm hbond -f "$XTC_FILE" -s "$GRO_FILE" -r chain_a -t chain_b -n "$index_file" -num "$hb_FILE" -cutoff 0.35 -m -o "$o_FILE"
    elif [ "$chain" = "chainb" ]; then
        $comm hbond -f "$XTC_FILE" -s "$GRO_FILE" -r chain_b -t chain_c -n "$index_file" -num "$hb_FILE" -cutoff 0.35 -m -o "$o_FILE"
    elif [ "$chain" = "chainc" ]; then
        $comm hbond -f "$XTC_FILE" -s "$GRO_FILE" -r chain_c -t chain_a -n "$index_file" -num "$hb_FILE" -cutoff 0.35 -m -o "$o_FILE"
    fi
done