VMD-L Mailing List

From: Thomas Schmidt (schmidt.nedlitz_at_hotmail.de)
Date: Wed Jun 12 2013 - 15:37:37 CDT

Hi Dear VMD Community

My name is Thomas, I am currently working in a NMR lab at USC (California). I performed a MD simulation on a membrane protein and would like to compare the order parameter (S^2) to the measured order parameter (NMR). I am trying to determine the S^2 value in VMD, for that I utilize a script which Justin Gullingsrud has written previously (lipid_order.tcl; see below #my adopted script is "Find order param for the aa of Protein" ) for the C2 and C3 tail in Lipids. In his previous script he utilized the C-H vectors of individual carbons in the lipid tail .... I am hoping to rewrite the script utilizing the N-H vector of the protein backbone. I was hoping you might have additional insights in what have to be changed in his previous script, in order to adopt it to proteins .... or you might have knowledge of a method developed by you or someone else which would calculate the order parameter of a protein.

Thank you in advance for your help,
Thomas

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is a set of scripts for finding the lipid order parameter from
# coordinate data. The lipid order parameter is defined as
# S_CD = < 3/2 cos^2 w - 1/2 >
# where w is the angle made by the bond between an aliphatic carbon
# and a hydrogen.

# Author: Justin Gullingsrud
# Email: justinrocks_at_gmail.com

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Find order param for the aa of Protein.
proc orderparam-protein { result { seltext all } } {
  upvar $result arr
  set n [molinfo top get numframes]
  for { set i 964 } { $i <= 990 } { incr i } {
    puts $i
    set cp [atomselect top "($seltext) and resid $i and name N"]
    if { [$cp num] == 0 } {
      puts "skipping $i"
      continue
    }
    set hx [atomselect top "($seltext) and resid $i and name H"]

    set sum 0.0
    set nh 0
    set nres [$cp num]
    for { set frame 0 } { $frame < $n } { incr frame } {
      $cp frame $frame
      $hx 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

    

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

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Find order param for the c3 tail.
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 0 } { $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}]
  }
}

# c2 tail
proc orderparam-c2 { 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 C2$i"]
    if { [$cp num] == 0 } {
      puts "skipping $i"
    }
    set hx [atomselect top "($seltext) and lipid and name H${i}R"]
    set hy [atomselect top "($seltext) and lipid and name H${i}S"]
    set hz [atomselect top "($seltext) and lipid and name H${i}T"]
    set h9 [atomselect top "($seltext) and lipid and name H${i}1"]
    set sum 0.0
    set nh 0
    set nres [$cp num]
    for { set frame 0 } { $frame < $n } { incr frame } {
      $cp frame $frame
      $hx frame $frame
      $hy frame $frame
$hz frame $frame
$h9 frame $frame
      set cpx [$cp get x]
      set cpy [$cp get y]
      set cpz [$cp get z]

      if { [$hx num] != 0 } {
        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
      }

      if { [$h9 num] != 0 } {
        set h9x [vecsub $cpx [$h9 get x]]
        set h9y [vecsub $cpy [$h9 get y]]
        set h9z [vecsub $cpz [$h9 get z]]
        foreach dx $h9x dy $h9y dz $h9z {
          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}]"
  }
}