From: Ioana Cozmuta (ioana_at_nas.nasa.gov)
Date: Thu Apr 29 2004 - 23:09:49 CDT

Hi VMD-users,

I am using the script below to calculate the average distance between
bases in a DNA strand during a MD simulation.
I am trying to make this script work in a pbs shell via
vmd -dispdev text -e script.vmd

1. At the moment I am defining the output file inside the procedure
myBaseDist. I've tried to define the name of the output file in the main
body of the script
set fa [open name.file w]
$fa global

and then inside the procedure use
global $fa

but vmd does not seem to recognize the name of the file in the procedure
this way.
What is the correct way to make this variable global?

2. Although my test trajectory file has 10 frames, the output contains
only 3 frames.

Could anyone please help me figure out why (without getting any errors)
the bigdcd procedure only processes 3 frames out of 10?

Thank you,
Ioana

************************************************

#This script calculates the distance between bases of a DNA strand

proc bigdcd { script args } {
  global bigdcd_frame bigdcd_proc bigdcd_firstframe vmd_frame

  set bigdcd_frame 0
  set bigdcd_firstframe [molinfo top get numframes]
  set bigdcd_proc $script

  uplevel #0 trace variable vmd_frame w bigdcd_callback
  foreach dcd $args {
    animate read dcd $dcd waitfor 0
  }
}

proc bigdcd_callback { name1 name2 op } {
  global bigdcd_frame bigdcd_proc bigdcd_firstframe vmd_frame

  # If we're out of frames, we're also done
  set thisframe $vmd_frame($name2)
  if { $thisframe < $bigdcd_firstframe } {
    bigdcd_done
    return
  }

  incr bigdcd_frame
  if { [catch {uplevel #0 $bigdcd_proc $bigdcd_frame} msg] } {
    puts stderr "bigdcd aborting at frame $bigdcd_frame\n$msg"
    bigdcd_done
    return
  }
  animate delete beg $thisframe end $thisframe
  return $msg
}

proc bigdcd_done { } {
  puts "bigdcd_done"
  after idle uplevel #0 trace vdelete vmd_frame w bigdcd_callback
}

proc dist_CM {select1 select2} {

if {[$select1 num] <= 0 || [$select2 num] <= 0 } {
error "center_of_mass: needs a selection with atoms"
 }

 set cm1 [veczero]
 set mass1 0
 set cm2 [veczero]
 set mass2 0

foreach crd1 [$select1 get {x y z}] m1 [$select1 get mass] {
 set mass1 [expr $mass1 + $m1]
 set cm1 [vecadd $cm1 [vecscale $m1 $crd1]]
 }

foreach crd2 [$select2 get {x y z}] m2 [$select2 get mass] {
 set mass2 [expr $mass2 + $m2]
 set cm2 [vecadd $cm2 [vecscale $m2 $crd2]]
}

 if {$mass1 == 0 || $mass2 == 0} {
  error "center_of_mass: total mass is zero"
 }

return [vecdist [vecscale [expr 1.0/$mass1] $cm1] [vecscale [expr 1.0/$mass2] $cm2]]

}

proc myBaseDist { frame } {
#global fa
set fa [open DistBDNA.dat a]
set residues [lsort -integer -unique [[atomselect top "resid 1 to 19 or resid 21 to 39"] get resid]]
set avg 0
set i 1
foreach residue $residues {
        set sel1 [atomselect top "resid $residue and not backbone"]
        set next [expr $residue + 1]
        set sel2 [atomselect top "resid $next and not backbone"]
        set dist [dist_CM $sel1 $sel2]
        set avg [expr $avg +$dist]
        set i [expr $i +1]
        }
set N $i
set avg [expr $avg/$N]
puts $fa "$frame: $avg"
close $fa
}

mol new {./m20AT/Eq20ATmod.restrt.coor} type {pdb}
# set fa [open DistBDNA.dat a]
# $fa global
bigdcd myBaseDist test.dcd

wait 1200
exit