VMD-L Mailing List

From: Ban Arn (ban.arn_at_gmail.com)
Date: Mon Apr 16 2012 - 05:16:31 CDT

Dear VMD Users

I'm using the following script to calculate lipid order parameter.

mol new lipid_initial.pdb
mol addfile lipid_traj.trr waitfor all

proc orderparam-c3 { result { seltext all } } {
  upvar $result arr
  set n [molinfo top get numframes]
  for { set i 2 } { $i <= 20 } { incr i } {
# puts $i
    set cp [atomselect top "($seltext) and lipid and name C3$i"]
    if { [$cp num] == 0 } {
# puts "skipping $i"
      continue
    }
    set hx [atomselect top "($seltext) and lipid and name H${i}X"]
    set hy [atomselect top "($seltext) and lipid and name H${i}Y"]
    set hz [atomselect top "($seltext) and lipid and name H${i}Z"]

    set sum 0.0
    set nh 0
    set nres [$cp num]
    for { set frame 500 } { $frame < $n } { incr frame } {
      $cp frame $frame
      $hx frame $frame
      $hy frame $frame
      $hz frame $frame
      set cpx [$cp get x]
      set cpy [$cp get y]
      set cpz [$cp get z]

      set hxx [vecsub $cpx [$hx get x]]
      set hxy [vecsub $cpy [$hx get y]]
      set hxz [vecsub $cpz [$hx get z]]
      foreach dx $hxx dy $hxy dz $hxz {
        set norm2 [expr {$dx*$dx + $dy*$dy + $dz*$dz}]
        set sum [expr {$sum + $dz*$dz/$norm2}]
      }
      incr nh $nres

      if { [$hy num] != 0 } {
        set hyx [vecsub $cpx [$hy get x]]
        set hyy [vecsub $cpy [$hy get y]]
        set hyz [vecsub $cpz [$hy get z]]
        foreach dx $hyx dy $hyy dz $hyz {
          set norm2 [expr {$dx*$dx + $dy*$dy + $dz*$dz}]
          set sum [expr {$sum + $dz*$dz/$norm2}]
        }
        incr nh $nres
      }

      if { [$hz num] != 0 } {
        set hzx [vecsub $cpx [$hz get x]]
        set hzy [vecsub $cpy [$hz get y]]
        set hzz [vecsub $cpz [$hz get z]]
        foreach dx $hzx dy $hzy dz $hzz {
          set norm2 [expr {$dx*$dx + $dy*$dy + $dz*$dz}]
          set sum [expr {$sum + $dz*$dz/$norm2}]
        }
        incr nh $nres
      }

    }
    set arr($i) [expr {-1.5*$sum/$nh + 0.5}]
   puts [expr {-1.5*$sum/$nh + 0.5}]
  }
}

The script works fine for 30 frames and later shows error that the two
vectors dont have the same size.

I couldn't figure out. Kindly advice.

Many Thanks
Balaji