Calculating hydrogen bonds life time.

I was asked to do such an analysis that is calculating average lifetime of hydrogen bonds formed between certain atom and surrounding water molecules . And I found that AMBERTOOLS had such a tool called “lifetime” but the calculation seems wrong when you set up certain values of “fuzz”. And I think it’s also easy to do this by myself, so I used following method to do the calculation.

  1. You need a trajectory file. I choose a 5ns trajectory, with 1ps per frame.
  2. Find surrounding water molecules of certain atom.
  3. Calculate hydrogen bonds formed between these surrounding water molecules and the target atom.
  4. Calculate the average lifetime of all hydrogen bonds formed between the target atom and surrounding water molecules in the 5ns.

The following is the script for Step 2 and 3. O1 on RES is the target atom.

rm -f *_temp.dat
for res in RES; do
 for atom in O1; do
  rm -f ${res}_${atom}.dat
cpptraj ../*.prmtop << EOF
trajin *.dcd
mask "(:${res}@${atom}<:3.6)&:WAT" maskout mask_temp.dat
run
exit
EOF
  water_res=`cat mask_temp.dat | awk '{print $4}' | sort -n | uniq -c | sed /Res/d | awk '{print $2}' | xargs echo`

   for water in $water_res; do
cpptraj ../*prmtop << EOF
trajin *.dcd
hbond A acceptormask :${res}@${atom} donormask :${water}@O= dist 3.5 angle 130.0 out hb_temp.dat
hbond B acceptormask :${water}@O= donormask :${res}@${atom} dist 3.5 angle 130.0 out hb_temp.dat
run
C = A[UU] + B[UU]
writedata ${res}_${atom}_temp.dat C
quit
EOF
     cat ${res}_${atom}_temp.dat >> ${res}_${atom}.dat
    done

 done
done

rm -f *_temp.dat

 

The following is the script for Step 4. All result will output file “lifetime.dat”. The “fuzz” here is the same with that in AMBERTOOLS.

rm -f lifetime.dat
((fuzz=0))
for res in RES; do
 for atom in O1; do
  echo ${res}_${atom}
  count=(`cat ${res}_${atom}.dat | awk '{print $2}' | cut -f1 -d"." | uniq -c | awk '{print $1}'`)
  hb=(`cat ${res}_${atom}.dat | awk '{print $2}' | cut -f1 -d"." | uniq -c | awk '{print $2}'`)
  for ((i=0;i<${#count[*]};i++)); do
   if [ ${hb[$i]} == 0 ] && [ ${count[$i]} == $fuzz ];then hb[$i]=1; fi
   if [ ${hb[$i]} == C ];then hb[$i]=0; fi
   ((temp[$i] = ( ${count[$i]} * ${hb[$i]} )))
  done
  hb_num=`echo ${temp[*]} | tr ' ' '\n' | uniq -c | awk '{print $2}' | grep -vw 0 | wc -l`
  hb_sum=`echo ${temp[*]} | tr ' ' '\n' | awk '{s+=$1} END {print s}'`
  hb_avg=`echo "$hb_sum / $hb_num * 1.0" | bc -l | cut -c -5`
  echo ${res}_${atom} $hb_avg >> lifetime.dat
  unset temp count hb
 done
done

 

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s