Changeset 1098 for branches/sandbox


Ignore:
Timestamp:
Nov 28, 2010 11:04:47 PM (10 years ago)
Author:
toby
Message:

add z-matrix conversion

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/rb.tcl

    r1036 r1098  
    582582}
    583583
     584# zmat2coord converts a z-matrix to a set of cartesian coordinates
     585#   a z-matrix is also known as "internal coordinates" or "torsion space"
     586#   (see Journal of Computational Chemistry, Vol 26, #10, p. 1063–1068, 2005 or
     587#    http://www.cmbi.ru.nl/molden/zmat/zmat.html)
     588# INPUT:
     589#   atmlist is a list of ascii lines where each line contains
     590#     lbl c1 distance c2 angle c3 torsion
     591#   where each atom is computed from the previous where the new atom is:
     592#     distance $distance from atom $c1 (angstrom)
     593#     angle $angle from $c1--$c2 (degrees)
     594#     torsion $torsion from $c1--$c2--$c3 (degrees)
     595# OUTPUT:
     596#  zmat2coord returns a list of atom labels and cartesian coordinates,
     597#  with 4 items in each element (label, x, y, z)
     598# this routine was tested against results from Babel via the web interface at
     599# http://www.shodor.org/chemviz/zmatrices/babel.html and sample input at
     600# http://iopenshell.usc.edu/howto/zmatrix/
     601proc zmat2coord {atmlist} {
     602    set torad [expr {acos(0)/90.}]
     603    set i 0
     604    foreach line $atmlist {
     605        incr i
     606        foreach {lbl c1 dist c2 angle c3 torsion} $line {}
     607        if {$i == 1} {
     608            set atm(1) [list $lbl 0 0 0] ; # 1st atom is at origin
     609        } elseif {$i == 2} {
     610            set dist1 $dist
     611            set atm(2) [list $lbl $dist1 0 0] ; # 2nd atom is along x-axis
     612        } elseif {$i == 3} {
     613            # 3rd atom can be bonded to the 1st or 2nd
     614            if {$c1 == 1} {
     615                set atm(3) [list $lbl \
     616                                [expr {$dist * cos($torad * $angle)}] \
     617                                [expr {$dist * sin($torad * $angle)}] \
     618                                0]
     619            } else {
     620                set atm(3) [list $lbl \
     621                                [expr {$dist1 - $dist * cos($torad * $angle)}] \
     622                                [expr {$dist * sin($torad * $angle)}] \
     623                                0]
     624            }
     625        } else {
     626            set atm($i) [concat $lbl \
     627                             [ahcat "atm" $c1 $dist $c2 $angle $c3 $torsion]]
     628        }
     629    }
     630    set coordlist {}
     631    foreach key [lsort -integer [array names atm]] {
     632        lappend coordlist $atm($key)
     633    }
     634    return $coordlist
     635}
     636# Compute the length of a vector
     637proc vlen {a} {
     638    set sum 0.0
     639    foreach ai $a {
     640        set sum [expr {$sum + $ai*$ai}]
     641    }
     642    return [expr sqrt($sum)]
     643}
     644# compute vector (a + z * b) and optionally normalize to length d
     645proc vadd {a b d z} {
     646    set c {}
     647    foreach ai $a bi $b {
     648        lappend c [expr {$bi + $z * $ai}]
     649    }
     650    set v [vlen $c]
     651    if {$d != 0} {
     652        set r {}
     653        foreach ci $c {
     654            lappend r [expr {$d * $ci / $v}]
     655        }
     656        return [list $v $r]
     657    }
     658    return [list $v $c]
     659}
     660# normalize a vector
     661proc vnrm {x} {
     662    set v [vlen $x]
     663    if {abs($v) < 1e-8} {return [list 0 0 0]}
     664    set y {}
     665    foreach xi $x {
     666        lappend y [expr {$xi / $v}]
     667    }
     668    return $y
     669}
     670# compute the coordinates for an atom that is bonded:
     671#   distance $dist from atom $nc
     672#   angle $bondang from $nc--$nb
     673#   torsion $torsang from $nc--$nb--$na
     674#   coordinates are found in array $atmarr in the calling routine
     675# based on a Fortran routine provided by Peter Zavalij (Thanks Peter!)
     676proc ahcat {atmarr nc dist nb bondang na torsang} {
     677    upvar 1 $atmarr atm
     678    set xa [lrange $atm($na) 1 3]
     679    set xb [lrange $atm($nb) 1 3]
     680    set xc [lrange $atm($nc) 1 3]
     681    set torad [expr {acos(0)/90.}]
     682    # x = unit Vector A-B
     683    foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. -1.] 1] {}
     684    # y = unit Vector C-B
     685    set y [lindex [vadd $xb $xc 1. -1.] 1]
     686    foreach {y1 y2 y3} $y {}
     687    set z1 [expr {$x2*$y3 - $x3*$y2}]
     688    set z2 [expr {$x3*$y1 - $x1*$y3}]
     689    set z3 [expr {$x1*$y2 - $x2*$y1}]
     690    set z [vnrm [list $z1 $z2 $z3]]
     691    set q1 [expr {$y2*$z3 - $y3*$z2}]
     692    set q2 [expr {$y3*$z1 - $y1*$z3}]
     693    set q3 [expr {$y1*$z2 - $y2*$z1}]
     694    set q [vnrm [list $q1 $q2 $q3]]
     695    set th [expr {$bondang * $torad}]
     696    set ph [expr {-1. * $torsang * $torad}]
     697    set cth [expr {cos($th)}]
     698    set sth [expr {sin($th)}]
     699    set cph [expr {cos($ph)}]
     700    set sph [expr {sin($ph)}]
     701    set xh {}
     702    foreach xci $xc xi $q zi $z yi $y {
     703        lappend xh [expr {
     704                          $xci +
     705                          $dist*($sth*($cph*$xi + $sph*$zi)-$cth*$yi)
     706                      }]
     707    }
     708    return $xh
     709}
     710 
     711
     712
    584713#AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
    585714#puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"]
     
    634763
    635764
     765# test zmat2coord
     766set atmlist {
     767    {C1 0 0.0 0 0.0 0 0.0}
     768    {O2 1 1.20 0 0.0 0 0.0}
     769    {H3 1 1.10 2 120.0 0 0.0}
     770    {C4 1 1.50 2 120.0 3 180.0}
     771    {H5 4 1.10 1 110.0 2 0.00}
     772    {H6 4 1.10 1 110.0 2 120.0}
     773    {H7 4 1.10 1 110.0 2 -120.0}
     774}
     775#  C        0.00000        0.00000        0.00000
     776#  O        1.20000        0.00000        0.00000
     777#  H       -0.55000        0.95263        0.00000
     778#  C       -0.75000       -1.29904       -0.00000
     779#  H       -0.04293       -2.14169       -0.00000
     780#  H       -1.38570       -1.36644        0.89518
     781#  H       -1.38570       -1.36644       -0.89518
     782set coordlist [zmat2coord $atmlist]
     783set i 0
     784foreach line $atmlist {
     785    incr i
     786    puts "$i) $line"
     787}
     788foreach line $coordlist {
     789    puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
     790}
Note: See TracChangeset for help on using the changeset viewer.