Convert Image to Grey Scale Image

convert logot.png -channel RGBA -matte -colorspace gray logotg.png

Cited from


Use Jump Host for RSYNC

I have computer A, storage server B and HPCC.

Previously, I use A to SSH B and RSYNC files from HPCC to B. One day, B can not PING or SSH or RSYNC or SCP to HPCC. But B can PING or SSH or RSYNC or SCP to A, and A can also communicate with HPCC. So I am trying to use A as a jump host to help B to download files from HPCC.


Got help form

The first thing is make sure the SSH KEY works between A and B, and between A and HPCC.

In ~/.ssh directory on B, build a new file config. The content is below.

User remote_server_user_id
ProxyCommand ssh nc 22

Change VMD startup script.

I use iMac. The version of VMD is 1.9.3. There is a .vmdrc file at /Applications/VMD

You can copy this file to your home directory ~ and then modify it.

What I want to modify is turn on all four lights, remove the axes and change the background color to gray.

The following modification is based on the original file from VMD.

  1. line 25: light 2 on
  2. line 26: light 3 on
  3. line 29: axes location off
  4. add color Display Background gray somewhere

Change VMD Movie Maker Default Setting.

I use it every week to make gif animation and show them in the group meeting. But the default setting is not proper for me. I use an iMac. The default working directory is lkjadfljaslkfjlasdfjlsadfdd……. (which is the system variable $TEMPDIR). The default animation format is MPEG-1. The default movement of molecules is Rock and Roll.


What I want is GIF, rotation along y axis and a specific working directory. If you need to generate tons of animation and sometimes VMD crashed, the repeatable setting will make me crazy.


Let’s work on this. The file is /Applications/VMD

The first thing I change is the working directory in Line 63, delete frames and write to your specific directory. And delete line 130 to 147, this part will find a working directory according to different system.


The second thing is to change the movement of molecules. In line 67, change rockandroll to rotation.


The third thing is to change the animation format. In line 80, change ppmtompeg to imgif.


That’s all.



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
  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
C = A[UU] + B[UU]
writedata ${res}_${atom}_temp.dat C
     cat ${res}_${atom}_temp.dat >> ${res}_${atom}.dat


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
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]} )))
  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



HEATMAP for intramolecular hydrogen bond analysis of a 10mer copolymer.

I am working on some simulations of certain copolymer. This copolymer could form intensive intramolecular hydrogen bonds. To observe how each residue interacts with other residues to form intramolecular hydrogen bonds, I choose HEATMAP to present the results.


So the above figure can show the average hydrogen bonds formed among different residue pairs on one copolymer during a period.


I did this in three step.

  1. Calculate intramolecular hydrogen bonds occupancy with VMD.
  2. Find pairs from the output files.
  3. Plot with GNUPLOT (version 5.0 patchlevel 1).


So easy. Mama will never worry about my study anymore ~~~~~~


The VMD script is below. My system only contains water and polymer, so I use “not water” to simplify my selection.

package require hbonds
mol new copolymer.prmtop type parm7 first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
mol addfile trajectory.dcd type dcd first 1 last -1 step 1 filebonds 1 autobonds 1 waitfor all
set sel1 [atomselect top "not water"]
set sel2 [atomselect top "not water"]
hbonds -polar yes -sel1 $sel1 -sel2 $sel2 -dist 3.5 -ang 50 -plot no -writefile yes -type pair -outfile vmd_hbond_chain.dat -detailout vmd_hbond_chain_detail.dat



The BASH to find pairs is below. This is 10mer copolymer. The names of residue in AMBER are MMH1 MMM2 MMM3 … MMM9 MME10. I chose perl to speed up the floating computing. It might save several minutes.

for i in {1..10};do echo \,${i};done | tr '\n' ' ' > hbond_intramolecular_occupancy.dat
echo >> hbond_intramolecular_occupancy.dat
for i in {1..10};do
 if (( $i == 1)); then j=MMH${i};fi
 if (( $i == 10)); then j=MME${i};fi
 for m in {1..10};do
  if (( $m == 1)); then n=MMH${m};fi
  if (( $m == 10)); then n=MME${m};fi
  a+=(`grep -E "${j}-.*${n}-" vmd_hbond_chain_detail.dat | awk '{print $3}' | cut -f1 -d"%" | perl -lne '$x += $_; END { print $x/100; }'`)
 echo ${i} ${a[*]} | tr ' ' ',' >> hbond_intramolecular_occupancy.dat
 unset a
unset a



The GNUPLOT script is below. I got help from And I did some slight modifications to fit my system.

set grid front
set xrange [*:]
set terminal pngcairo
set encoding iso_8859_1
set xl "Acceptor Residue"
set yl "Donor Residue"
# Color runs from white to green
set palette rgbformula -7,2,-7
set cbrange [0:2]
set cblabel "Number of Hydrogen Bonds per Frame"
set cbtics

set output "hbond_intramolecular_occupancy.png"
set title "Intramolecular Hydrogen Bonds"
set datafile separator ","
plot 'hbond_intramolecular_occupancy.dat' matrix rowheaders columnheaders using 1:2:3 with image noti
set datafile separator
set output