Visualize connected voids in the unit cell

In the previous post, I have shown how to find the void in VMD. But if you want to visualize the connected voids, it will become a little complicated. At the beginning, I tried to use watershed algorithm to do that, during which I had an idea to do that in a much simpler way. I will introduce below.

  1. Use the following Tcl script to find all voids and put a Si atom in each void. And the end of the script, a XYZ file is generated to use Si to represent the voids in the unit cell.
  2. Use VMD to open the XYZ file.
  3. In the representation setting, you can select “fragment 0” to “fragment N” to see the connected voids. VMD can figure out the connected voids and treat them as different fragments. With Selection tool, you can even move the connected voids to anywhere you want for better representation.
  4. You can also use color method “fragment” to use different colors to highlight different fragments, but the default color scheme is not good.


#Change background color to white
color Display Background white

#set sphere radius
set probe 0.8

#set the box side-length, this is a cube
set length 39.91

#the side-length is divided into N grid, 50 is enough to visualize
set step [expr { $length/1.4 }]

#delete all graphics
draw delete all

#draw box, wrap the system
pbc box
pbc wrap
#draw color green
#recalculate the bond formation
mol bondsrecalc 0

set outfile [open "" w+]

set i 0
#find void
for {set x 0.5} {$x <= $length} { set x [expr {$x + $length/$step}]} {
        for {set y 0.5} {$y <= $length} { set y [expr {$y + $length/$step}]} {
                for {set z 0.5} {$z <= $length} { set z [expr {$z + $length/$step}]} {
                        set a [[atomselect top "sqrt( sqr(x-$x) +sqr(y-$y) + sqr(z-$z) )  < 2"] num]
                        if { $a == 0 } {
#                               draw sphere "$x $y $z" radius $probe
                                incr i
                                set X($i) $x; set Y($i) $y; set Z($i) $z


puts $outfile $i
puts $outfile "Connected Pores"
for {set m 1} {$m <= $i } {incr m} {
        puts $outfile "Si $X($m)  $Y($m)  $Z($m)"
close $outfile
mol new



Steered MD simulation in AMBER using Torsion as the reaction coordination.

I need to generate a set of trajectory of the rotation of one dihedral in a polymer. AMBERTOOLS14 has a command “dihedralscan” to do that, but it only works for dihedrals of proteins. It’s not allowed to define your own dihedral. So I chose SMD to do this.

It’s not complicated to set up the system. The only confusing part is the value of coordinate of torsion is rad not degree.

The input file of NPT part is listed below in case you might need. The initial torsion was about -1rad. To rotate 360degree. I set the reaction coordinate change from the initial value to 5.3. The diehdral is on residue 10 ( atom  100, 196, 155, 175). The rest parts are restrained with 500kcal/mol/angstrom^2.

  imin = 0, irest = 1, ntx = 7,
  ntb = 2, pres0 = 1.0, ntp = 1,
  taup = 2.0,
  cut = 10.0, ntr = 1,
  ntc = 1, ntf = 1,
  tempi = 300.0, temp0 = 300.0,
  ntt = 3, gamma_ln = 1.0,
  nstlim = 100000, dt = 0.001,
  ntpr = 1000, ntwx = 1000, ntwr = 1000, ig = -1,
Weak Hold
RES 1 9
RES 11 20

 output_file = 'smd.txt'
 output_freq = 1000
  type = TORSION
  i = (100, 196, 155, 175)
  path = (X, 5.3) path_mode = LINES
  harm = (100.0)
 end variable
end ncsu_smd


DSSP refers to Dictionary of Secondary Structures of Proteins. It’s a rough method to analyze protein secondary structures along with time in MD simulations.

The following is an example. Cited from

There are many softwares or MD simulation package can do such an analysis. I prefer to use AMBERTOOLS.

The script is below.

parm test.prmtop
trajin test.mdcrd
secstruct :1-100 out dssp.gnu

The GNU file has some problem and you can not get a picture from it.

I did the following modification to generate PNG file or PDF file.

  1. At the end of GNU file change pause -1 to set output.
  2. Delete line 1 to 12. The error part is at line 5, which is the setting of cbtics.
  3. Change line 13 splot “-” with pm3d title “dssp.gnu” to splot “-” u ($1/10):2:3 w pm3d title “dssp.gnu”. You may change $1/10 according to your relationship between frame and time.

The head of GNU file can be changed to the following content. 

unset key
set grid
set terminal pngcairo
set size 0.96,1
set encoding iso_8859_1
set xtics nomirror
set ytics nomirror
set pm3d map corners2color c1
set cbrange [  -0.500:   7.500]
set cbtics    0.000,1,7
set palette maxcolors 8
set palette defined (0 "white",1 "#0000FF",4 "#00FF00", 5 "grey", 7 "#FF0000")
set cbtics("None"    0.000,"Para"    1.000,"Anti"    2.000,"3-10"    3.000,"Alpha"    4.000,"Pi"    5.000,"Turn"    6.000,"Bend"    7.000)
set xlabel "Time (ns)"
set ylabel "Residue Number"
set yrange [   0.000: 600]
set xrange [   0.000: 84]
set title "DSSP of this protein"
set output "dssp.png"

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



Systems with multiple segments in CHARMM.

If you read through my blog, you will see I did post one blog about setting up PEG in a TIP3P water box with CHARMM ether force field.


I got new tasks from my dear advisor. Two PEGs in one TIP3P water box. In AMBER, this is pretty simple, I just need to prepare a PDB files contains two PEG and add TER between the two molecules. Never mind, my story in this blog is using CHARMM ether force field for PEG. By the way, PEG with GAFF seemed not that rational compared with that using CHARMM ether force field.


As a beginner of CHARMM, the learning process is really painful. I got help from I used a multiple chains protein in PDB as input files, and look into the input scripts. I find the clue that CHARMM will generate multiple PSF files for each segments at the beginning, and then attached all PSF files together as a whole.


Based on the script from the website, I tried several times and got the answer. The following procedure is what I used for the system.

  1. I got the PDB files of one PEG from my previous simulation.
  2. Use VMD to prepare two PEG PDB files with proper coordinates.
  3. Save the two PEG PDB files separately as 1.pdb and 2.pdb.
  4. Use the following 1.inp to generate 1.psf, and change 1 to 2 and generate 2.psf.
  5. Use combine.inp to generate two PEG vacuum PSF file.
  6. Use VMD to generate PSF and PDB of the water box.



* Run Segment Through CHARMM

dimension chsize 1500000

prnlev 5
wrnlev 5
bomlev -1

! read topology and parameter files

wrnlev 1
wrnlev 5

read rtf card name ../top_all35_ethers-oh.rtf
read param card name ../par_all35_ethers-oh.prm

read sequence pdb name 1.pdb
generate a-good setu first GCL0 last GCL3 noangle nodihedral

read coor pdb name 1.pdb


! write out the protein structure file (psf) and
! the coordinate file in pdb and crd format.

write psf card name 1.psf

write coor card name 1.crd
* Coords

write coor pdb name 1_charmm.pdb
* PDB for JSmol visualization





* Append multiple segments into a final structure

dimension chsize 1500000

prnlev 5
wrnlev 5
bomlev -1

! read topology and parameter files
read rtf card name ../top_all35_ethers-oh.rtf

wrnlev 1
read param card name ../par_all35_ethers-oh.prm
wrnlev 5

! Read PSF from file
read psf card name 1.psf
read psf card append name 2.psf

autogenerate angles dihedrals

! read coordinates
read coor card name 1.crd
read coor card append name 2.crd

! orient structure
!! BTM, don't do this, it might muck with re-assembling the
!! structure for redox calculations
!coor orient

! calculate energy

ioform extended

write psf card name twoPEGvac.psf
* twoPEGvac.psf

write coor card name twoPEGvac.crd
* twoPEGvac.crd

write coor pdb name twoPEGvac.pdb
* twoPEGvac.pdb