# Changeset 1098 for branches

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

add z-matrix conversion

File:
1 edited

Unmodified
Added
Removed
• ## branches/sandbox/rb.tcl

 r1036 } # 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"] # 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] }
Note: See TracChangeset for help on using the changeset viewer.