Changeset 1098

Timestamp:

Nov 28, 2010 11:04:47 PM (10 years ago)

Author:

toby

Message:

add z-matrix conversion

File:

• branches/sandbox/rb.tcl (modified) (2 diffs)

branches/sandbox/rb.tcl

@@ -582 +582 @@
 }
 
+# zmat2coord converts a z-matrix to a set of cartesian coordinates
+#   a z-matrix is also known as "internal coordinates" or "torsion space"
+#   (see Journal of Computational Chemistry, Vol 26, #10, p. 1063–1068, 2005 or
+#    http://www.cmbi.ru.nl/molden/zmat/zmat.html)
+# INPUT:
+#   atmlist is a list of ascii lines where each line contains
+#     lbl c1 distance c2 angle c3 torsion
+#   where each atom is computed from the previous where the new atom is: 
+#     distance $distance from atom $c1 (angstrom)
+#     angle $angle from $c1--$c2 (degrees)
+#     torsion $torsion from $c1--$c2--$c3 (degrees)
+# OUTPUT:
+#  zmat2coord returns a list of atom labels and cartesian coordinates, 
+#  with 4 items in each element (label, x, y, z)
+# this routine was tested against results from Babel via the web interface at 
+# http://www.shodor.org/chemviz/zmatrices/babel.html and sample input at 
+# http://iopenshell.usc.edu/howto/zmatrix/
+proc zmat2coord {atmlist} { 
+    set torad [expr {acos(0)/90.}]
+    set i 0
+    foreach line $atmlist {
+        incr i
+        foreach {lbl c1 dist c2 angle c3 torsion} $line {} 
+        if {$i == 1} {
+            set atm(1) [list $lbl 0 0 0] ; # 1st atom is at origin
+        } elseif {$i == 2} {
+            set dist1 $dist
+            set atm(2) [list $lbl $dist1 0 0] ; # 2nd atom is along x-axis
+        } elseif {$i == 3} {
+            # 3rd atom can be bonded to the 1st or 2nd
+            if {$c1 == 1} {
+                set atm(3) [list $lbl \
+                                [expr {$dist * cos($torad * $angle)}] \
+                                [expr {$dist * sin($torad * $angle)}] \
+                                0]
+            } else {
+                set atm(3) [list $lbl \
+                                [expr {$dist1 - $dist * cos($torad * $angle)}] \
+                                [expr {$dist * sin($torad * $angle)}] \
+                                0]
+            }
+        } else {
+            set atm($i) [concat $lbl \
+                             [ahcat "atm" $c1 $dist $c2 $angle $c3 $torsion]]
+        }
+    }
+    set coordlist {}
+    foreach key [lsort -integer [array names atm]] {
+        lappend coordlist $atm($key)
+    }
+    return $coordlist
+}
+# Compute the length of a vector
+proc vlen {a} {
+    set sum 0.0
+    foreach ai $a {
+        set sum [expr {$sum + $ai*$ai}]
+    }
+    return [expr sqrt($sum)]
+}
+# compute vector (a + z * b) and optionally normalize to length d
+proc vadd {a b d z} {
+    set c {}
+    foreach ai $a bi $b {
+        lappend c [expr {$bi + $z * $ai}]
+    }
+    set v [vlen $c]
+    if {$d != 0} {
+        set r {}
+        foreach ci $c {
+            lappend r [expr {$d * $ci / $v}]
+        }
+        return [list $v $r]
+    }
+    return [list $v $c]
+}
+# normalize a vector
+proc vnrm {x} {
+    set v [vlen $x]
+    if {abs($v) < 1e-8} {return [list 0 0 0]}
+    set y {}
+    foreach xi $x {
+        lappend y [expr {$xi / $v}]
+    }
+    return $y
+}
+# compute the coordinates for an atom that is bonded:
+#   distance $dist from atom $nc
+#   angle $bondang from $nc--$nb
+#   torsion $torsang from $nc--$nb--$na
+#   coordinates are found in array $atmarr in the calling routine
+# based on a Fortran routine provided by Peter Zavalij (Thanks Peter!)
+proc ahcat {atmarr nc dist nb bondang na torsang} {
+    upvar 1 $atmarr atm
+    set xa [lrange $atm($na) 1 3]
+    set xb [lrange $atm($nb) 1 3]
+    set xc [lrange $atm($nc) 1 3]
+    set torad [expr {acos(0)/90.}]
+    # x = unit Vector A-B
+    foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. -1.] 1] {}
+    # y = unit Vector C-B
+    set y [lindex [vadd $xb $xc 1. -1.] 1]
+    foreach {y1 y2 y3} $y {}
+    set z1 [expr {$x2*$y3 - $x3*$y2}]
+    set z2 [expr {$x3*$y1 - $x1*$y3}]
+    set z3 [expr {$x1*$y2 - $x2*$y1}]
+    set z [vnrm [list $z1 $z2 $z3]]
+    set q1 [expr {$y2*$z3 - $y3*$z2}]
+    set q2 [expr {$y3*$z1 - $y1*$z3}]
+    set q3 [expr {$y1*$z2 - $y2*$z1}]
+    set q [vnrm [list $q1 $q2 $q3]]
+    set th [expr {$bondang * $torad}]
+    set ph [expr {-1. * $torsang * $torad}]
+    set cth [expr {cos($th)}]
+    set sth [expr {sin($th)}]
+    set cph [expr {cos($ph)}]
+    set sph [expr {sin($ph)}]
+    set xh {}
+    foreach xci $xc xi $q zi $z yi $y {
+        lappend xh [expr {
+                          $xci +
+                          $dist*($sth*($cph*$xi + $sph*$zi)-$cth*$yi)
+                      }]
+    }
+    return $xh
+}
+ 
+
+
 #AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
 #puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"]
@@ -634 +763 @@
 
 
+# test zmat2coord
+set atmlist {
+    {C1 0 0.0 0 0.0 0 0.0}
+    {O2 1 1.20 0 0.0 0 0.0}
+    {H3 1 1.10 2 120.0 0 0.0}
+    {C4 1 1.50 2 120.0 3 180.0}
+    {H5 4 1.10 1 110.0 2 0.00}
+    {H6 4 1.10 1 110.0 2 120.0}
+    {H7 4 1.10 1 110.0 2 -120.0}
+}
+#  C        0.00000        0.00000        0.00000
+#  O        1.20000        0.00000        0.00000
+#  H       -0.55000        0.95263        0.00000
+#  C       -0.75000       -1.29904       -0.00000
+#  H       -0.04293       -2.14169       -0.00000
+#  H       -1.38570       -1.36644        0.89518
+#  H       -1.38570       -1.36644       -0.89518
+set coordlist [zmat2coord $atmlist]
+set i 0
+foreach line $atmlist {
+    incr i
+    puts "$i) $line"
+}
+foreach line $coordlist {
+    puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
+}