Changeset 1098
 Timestamp:
 Nov 28, 2010 11:04:47 PM (11 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/sandbox/rb.tcl
r1036 r1098 582 582 } 583 583 584 # zmat2coord converts a zmatrix to a set of cartesian coordinates 585 # a zmatrix 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/ 601 proc 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 xaxis 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 637 proc 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 645 proc 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 661 proc vnrm {x} { 662 set v [vlen $x] 663 if {abs($v) < 1e8} {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!) 676 proc 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 AB 683 foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. 1.] 1] {} 684 # y = unit Vector CB 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 584 713 #AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} } 585 714 #puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"] … … 634 763 635 764 765 # test zmat2coord 766 set 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 782 set coordlist [zmat2coord $atmlist] 783 set i 0 784 foreach line $atmlist { 785 incr i 786 puts "$i) $line" 787 } 788 foreach 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.