source: branches/sandbox/rb.tcl @ 1152

Last change on this file since 1152 was 1152, checked in by toby, 12 years ago

a bit of cleanup
File size: 39.4 KB
Line 
1#============================================================================
2# Rigid body utility routines
3#============================================================================
4# RigidBodyGetVarNums: Returns a list of the variable numbers in use
5#       for rigid body variable parameters.
6# RigidBodyAtomNums: returns a list of atom numbers that are mapped to
7#       rigid bodies in a selected phase
8# RigidStartAtoms: returns a list of atoms that are allowed for creation of RB
9# ExtractRigidBody: Use the GSAS geometry program to cartesian coordinates &
10#       setting info for a RB from fractional coordinates for atoms in a phase
11# RunRecalcRBCoords: updates the coordinates in all phases after changes have
12#       been made to rigid parameters.
13# CalcBody: Convert ortho to fractional coordinates using RB parameters
14# FitBody: Optimize the origin and Euler angles to match a rigid body to a
15#       set of fractional coordinates
16# zmat2coord: convert a z-matrix to a set of cartesian coordinates
17# RB2cart: convert the representation used for rigid bodies into
18#       cartesian coordinates
19# PlotRBtype: plot a rigid body with DRAWxtl
20# PlotRBcoords: plot orthogonal coordinates with DRAWxtl
21# DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed
22#       on a structure read from a phase with DRAWxtl
23#============================================================================
24#============================================================================
25# RigidBodyGetVarNums: Returns a list of the variable numbers used already
26# for rigid body variable parameters
27proc RigidBodyGetVarNums {} {
28    set varlist {}
29    foreach type [RigidBodyList] {
30        set typelist [lindex [ReadRigidBody $type] 1]
31        foreach item $typelist {
32            lappend varlist [lindex $item 2]
33        }
34        foreach phase $::expmap(phaselist) {
35            foreach i [RigidBodyMappingList $phase $type] {
36                set items [ReadRigidBodyMapping $phase $type $i]
37                set varlist [concat $varlist [lindex $items 3]]
38                if {[llength [lindex $items 6]] > 0} {
39                    set varlist [concat $varlist [lindex $items 6]]
40                }
41            }
42        }
43    }
44    return [lsort -integer -unique $varlist]
45}
46
47# RigidBodyAtomNums: Returns a list of the atoms mapped to rigid bodies in
48# phase $phase
49proc RigidBodyAtomNums {phase} {
50    if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""}
51    set allatoms $::expmap(atomlist_$phase)
52    # get matching atoms coordinate range
53    set mappedlist {}
54    foreach type [RigidBodyList] {
55        foreach i [RigidBodyMappingList $phase $type] {
56            # get the number of atoms in this type of body
57            set natoms [llength [lindex [lindex [lindex [ReadRigidBody $type] 1] 0] 3]]
58            set natom1 [expr {$natoms - 1}]
59            set items [ReadRigidBodyMapping $phase $type $i]
60            set firstatom [lindex $items 0]
61            set firstind [lsearch $allatoms $firstatom]
62            set mappedlist [concat $mappedlist \
63                                [lrange \
64                                     [lrange $allatoms $firstind end] \
65                                     0 $natom1] \
66                               ]
67        }
68    }
69    return [lsort -integer $mappedlist]
70}
71
72# RigidStartAtoms: Find allowed starting atoms for a rigid body in a phase
73# Input:
74#   phase is the phase number
75#   natoms is the number of atoms in the RB to be mapped
76# Returns a list of valid "start" atoms.
77# Example: if the atom numbers in the phase are {2 4 5 6 7 8} and no rigid bodies
78# are mapped, then a 4-atom body can be mapped starting with atom 2, 4 or 5 only,
79# so {2 4 5} is returned
80# If atoms 2-6 were already mapped, then this routine would return an empty
81# list, as it is not possible to map the body.
82proc RigidStartAtoms {phase natoms} {
83    if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""}
84    set allatoms $::expmap(atomlist_$phase)
85    set usedatoms [RigidBodyAtomNums $phase]
86    set startatomlist {}
87    for {set i 0} {$i < [llength $allatoms]} {incr i} {
88        set al [lrange $allatoms $i [expr {$i+$natoms-1}]]
89        if {[llength $al] < $natoms} break
90        set ok 1
91        foreach atom $al {
92            if {[lsearch $usedatoms $atom] != -1} {
93                set ok 0
94                break
95            }
96        }
97        if $ok {lappend startatomlist [lindex $al 0]}
98    }
99    return $startatomlist
100}
101
102# ExtractRigidBody: Use the GSAS geometry program to compute a set of cartesian coordinates for a
103# set of atoms in a .EXP file and provide the origin shift and Euler angles needed to
104# place the cartesian system into the crystal coordinates. Used for setting up a rigid body.
105# Returns a nested list of lists:
106#   element 0: a list of the origin location {x y z} in fraction coordinates
107#   element 1: a list of three rotation angles in form {a1 a2 a3}
108#              where a1, a2 and a3 are rotations around the cartesian x, y and z axes
109#   element 2: a list of $natom cartesian coordinate triples {{x1 y1 z1} {x2 y2 z2}...}
110# arguments:
111    # phase: phase #
112    # natom: number of atoms in group
113    # firstatom: sequence # in phase (may be > than number of the atom)
114    # originlist: atoms to define origin (where 1 is first atom in group; <= natom)
115    # vector1: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)  (for example {X 1 3})
116    # vector2: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)
117    # note that vector2 must define a different axis than vector1
118    # also and vector1 and vector2 cannot use the same atom pair
119proc ExtractRigidBody {phase natom firstatom originlist vector1 vector2} {
120    global expgui
121    set fp [open "geom.inp" "w"]
122    puts $fp "N"
123    if {[llength ::expmap(phaselist)] > 1} {
124       # select phase
125       puts $fp "N"
126       puts $fp $phase
127       puts $fp "N"
128    }
129    puts $fp "R"
130    puts $fp "$natom"
131    puts $fp "$firstatom"
132    puts $fp [llength $originlist]
133    foreach i $originlist {
134        puts $fp $i
135    }
136    foreach i [concat $vector1 $vector2] {
137        puts $fp $i
138    }
139    puts $fp "0"
140    puts $fp "X"
141    close $fp
142    #puts "[file join $expgui(gsasexe) geometry] $expgui(expfile) < geom.inp > geom.out"
143    catch {
144        exec [file join $expgui(gsasexe) geometry] $expgui(expfile) < geom.inp > geom.out
145    } err
146    #puts $err
147    file delete geom.inp
148    set fp [open geom.out r]
149    set origin {}
150    set Euler {}
151    set coordlist {}
152    while {[gets $fp line] >= 0} {
153        if {[string first "Cell coordinates of origin" $line] != -1} {
154            set origin [lrange [string range $line [string first "are" $line] end] 1 3]
155            #puts "origin in rb = $origin"
156        }
157        if {[string first "Rotation angles" $line] != -1} {
158            set Euler {}
159            foreach i [lrange [split $line "="] 1 3] {
160                lappend Euler [lindex $i 0]
161            }
162            #puts $line
163            #puts $Euler
164        }
165        if {[string first "Atom   Orthon" $line] != -1} {
166            set coordlist {}
167            for {set i 1} {$i <= $natom} {incr i} {
168                gets $fp line
169                set coord {}
170                lappend coord [string trim [string range $line 9 15]]
171                lappend coord [string trim [string range $line 16 22]]
172                lappend coord [string trim [string range $line 23 29]]
173                lappend coord [string trim [string range $line 0 8]]
174                #puts $line
175                #puts $coord
176                lappend coordlist $coord
177            }
178            #puts $coordlist
179        }
180    }
181    #file delete geom.out
182    if {[llength $origin] == 0 || [llength $Euler] == 0 || [llength $coordlist] == 0} {
183       puts "Error: run of GEOMETRY failed"
184    }
185    return [list $origin $Euler $coordlist]
186}
187
188# RunRecalcRBCoords: updates the coordinates in a .EXP file after a rigid
189# body has been changed, mapped or the setting info is changed
190proc RunRecalcRBCoords { } {
191    global expgui tcl_platform
192    set input [open resetmult.inp w]
193    puts $input "Y"
194    puts $input "l b"
195    puts $input "n"
196    puts $input "x x x"
197    puts $input "x"
198    close $input
199    # Save the current exp file
200    savearchiveexp
201    # disable the file changed monitor
202    set expgui(expModifiedLast) 0
203    set expnam [file root [file tail $expgui(expfile)]]
204    set err [catch {
205        if {$tcl_platform(platform) == "windows"} {
206            exec [file join $expgui(gsasexe) expedt.exe] $expnam < resetmult.inp >& resetmult.out
207        } else {
208            exec [file join $expgui(gsasexe) expedt] $expnam < resetmult.inp >& resetmult.out
209        }
210    } errmsg]
211    loadexp $expgui(expfile)
212    set fp [open resetmult.out r]
213    set out [read $fp]
214    close $fp
215    set expgui(exptoolout) $out
216    catch {file delete resetmult.inp resetmult.out}
217    if {$err} {
218        return $errmsg
219    } else {
220        return ""
221    }
222}
223
224
225# compute a rotation matrix for orthogonal coordinates (based on MAKMATD in GSAS)
226# rotate angle degrees around axis (1, 2 or 3) for (x, y, or z)
227# returns a list that can be used as a matrix in the La package
228proc RotMatrix {axis angle} {
229    set ang [expr {$angle * acos(0) / 90.}]
230    set mat "1 0 0 0 1 0 0 0 1"
231    if {$axis == 1}  {
232        set i1 1
233        set i2 2
234    } elseif {$axis == 2}  {
235        set i1 2
236        set i2 0
237    } else {
238        set i1 0
239        set i2 1
240    }
241    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
242    foreach item {
243        {$i1 $i1 [expr {cos($ang)}]}
244        {$i2 $i2 [expr {cos($ang)}]}
245        {$i1 $i2 [expr {-sin($ang)}]}
246        {$i2 $i1 [expr {sin($ang)}]}
247    } {
248        foreach {c r val} [subst $item] {}
249        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
250    }
251    return "2 3 3 $mat"
252}
253
254# compute the derivative of the rotation matrix with respect to the angle, see RotMatrix
255# (based on MAKMATD in GSAS)
256# returns a list that can be used as a matrix in the La package
257proc DerivRotMatrix {axis angle} {
258    set ang [expr {$angle * acos(0) / 90.}]
259    set mat "0 0 0 0 0 0 0 0 0"
260    if {$axis == 1}  {
261        set i1 1
262        set i2 2
263    } elseif {$axis == 2}  {
264        set i1 2
265        set i2 0
266    } else {
267        set i1 0
268        set i2 1
269    }
270    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
271    foreach item {
272        {$i1 $i1 [expr {-sin($ang) * acos(0) / 90.}]}
273        {$i2 $i2 [expr {-sin($ang) * acos(0) / 90.}]}
274        {$i1 $i2 [expr {-cos($ang) * acos(0) / 90.}]}
275        {$i2 $i1 [expr {cos($ang) * acos(0) / 90.}]}
276    } {
277        foreach {c r val} [subst $item] {}
278        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
279    }
280    return "2 3 3 $mat"
281}
282
283# compute an orthogonalization matrix from cell parameters (based on AMATRX in GSAS)
284# returns a list that can be used as a matrix in the La package
285proc OrthoMatrix {a b c alpha beta gamma} {
286    set CA [expr {cos($alpha * acos(0) / 90.)}]
287    set CB [expr {cos($beta * acos(0) / 90.)}]
288    set CG [expr {cos($gamma * acos(0) / 90.)}]
289    set SA [expr {sin($alpha * acos(0) / 90.)}]
290    set SB [expr {sin($beta * acos(0) / 90.)}]
291    set SC [expr {sin($gamma * acos(0) / 90.)}]
292    set CASTAR [expr { ($CB*$CG-$CA)/($SB*$SC) }]    ;#! cos(Alpha*)
293    set CBSTAR [expr { ($CA*$CG-$CB)/($SA*$SC) }]    ;#! cos(Beta*)
294    set CCSTAR [expr { ($CA*$CB-$CG)/($SA*$SB) }]    ;#! cos(Gamma*)
295    set SASTAR [expr { sqrt(1.0-($CASTAR*$CASTAR*2)) }]    ;#! sin(Alpha*)
296    set SBSTAR [expr { sqrt(1.0-($CBSTAR*$CBSTAR*2)) }]    ;#! sin(Beta*)
297    set SCSTAR [expr { sqrt(1.0-($CCSTAR*$CCSTAR*2)) }]    ;#! sin(Gamma*)
298
299    set A  "2 3 3      $a 0 0    0 $b 0    0 0 $c"
300    set A1 "2 3 3     1.0 0 0    $CG [expr {$SASTAR*$SC}] [expr {-$CASTAR*$SC}]    $CB 0.0 $SB"
301    #!This matrix is
302    #!   (1.0      0.0            0.0             )
303    #!   (cos(G)  sin(A*)*sin(G) -cos(A*)*sin(G)  )
304    #!   (cos(B)   0.0            sin(B)          )
305    return [La::transpose [La::mmult $A $A1]]
306}
307
308# compute the transformation matrix that converts a rigid body coordinates into fractional
309# coordinates
310# arguments:
311#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
312#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
313#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
314# returns a list that can be used as a matrix in the La package
315proc CalcXformMatrix {rotations cellprms} {
316    set prod {}
317    foreach item $rotations {
318        #puts $item
319        set mat [eval RotMatrix $item]
320        if {$prod == ""} {
321            set prod $mat
322        } else {
323            set prod [La::mmult $prod $mat]
324        }
325    }
326    #puts "--- rotation product ---"
327    #puts [La::show $prod]
328
329    set ortho [eval OrthoMatrix $cellprms]
330    #puts "--- ortho ---"
331    #puts [La::show $ortho]
332    set deortho [La::msolve $ortho [La::mident 3] ]
333    #puts "--- deortho ---"
334    #puts [La::show $deortho]
335    #puts "--- xform ---"
336    set xform [La::mmult $deortho $prod]
337    return $xform
338}
339
340# transforms a triplet of orthogonal coordinates into fractional ones using
341# arguments:
342#    xform: a transformation matrix from CalcXformMatrix
343#    origin: a list of fraction coordinates {x y z} describing the location of the
344#            origin of the orthogonal coordinates in the crystal system
345#    ortho: a triplet of othogonal coordinates
346# returns a triplet of fractional coordinates
347proc Ortho2Xtal {xform origin ortho} {
348    set ocv "2 3 0 $ortho"
349    set frac [La::mmult $xform $ocv]
350    #puts [La::show [La::transpose $frac]]
351    #puts $frac
352    set frac [La::madd $frac "[lrange $frac 0 2] $origin"]
353    #puts [La::show [La::transpose $frac]]
354    return $frac
355}
356
357# compute the derivative of the transformation matrix (see CalcXformMatrix)
358# with respect to every rotation in the $rotations list
359# arguments:
360#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
361#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
362#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
363# returns a list of matrices where each matrix is a list that can be used as a
364# matrix in the La package
365proc CalcDerivMatrix {rotations cellprms} {
366    set ortho [eval OrthoMatrix $cellprms]
367    set deortho [La::msolve $ortho [La::mident 3] ]
368    set derivlist {}
369
370    foreach item $rotations {
371        #puts $item
372        set mat [eval DerivRotMatrix $item]
373        #puts $item
374        #puts [La::show $mat]
375        set xform [La::mmult $deortho $mat]
376        lappend derivlist $xform
377    }
378    return $derivlist
379}
380
381# CalcBody: Calculate fractional coordinates using rigid body setting parameters
382# arguments:
383#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
384#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
385#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
386#  ortholist: list containing triplets with orthogonal coordinates
387#  origin: a list of fraction coordinates {x y z} describing the location of the
388#            origin of the orthogonal coordinates in the crystal system
389#     note that the length of ortholist, useflag and fraclist should be the same
390# Returns a list with the computed fractional coordinates for all atoms
391proc CalcBody {Euler cell ortholist origin} {
392    set xform [CalcXformMatrix $Euler $cell]
393    set i 0
394    set sumdvs 0
395    set fracout {}
396    set rmsout {}
397    foreach oc $ortholist {
398        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
399        lappend fracout $frac
400    }
401    return $fracout
402}
403
404
405# fit a rigid body's origin
406# arguments:
407#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
408#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
409#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
410#  ortholist: list containing triplets with orthogonal coordinates
411#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
412#  fraclist: list containing triplets with fractional coordinates 
413#  origin: a list of fraction coordinates {x y z} describing the location of the
414#            origin of the orthogonal coordinates in the crystal system
415#     note that the length of ortholist, useflag and fraclist should be the same
416# Returns a list with the following elements
417#   0: a list with three new origin values
418#   1: the root-mean square difference between the fraclist coordinates and those computed
419#      using the input values for those atoms where $use is one (in Angstroms)
420#   2: the computed fractional coordinates for all atoms
421#   3: individual rms values for all atoms (in Angstroms)
422# note that items 1-3 are computed with the imput origin, not the revised one
423proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} {
424puts $fraclist
425    set xform [CalcXformMatrix $Euler $cell]
426    #puts "entering FitBodyOrigin"
427    foreach var {x y z} {set sum($var) 0.0}
428    set i 0
429    set sumdvs 0
430    set fracout {}
431    set rmsout {}
432    foreach oc $ortholist use $useflag coord $fraclist {
433        #puts "ortho: $oc"
434        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
435        lappend fracout $frac
436        if {$use} {incr i}
437        set dvs 0
438        foreach var {x y z} v1 $frac v2 $coord abc [lrange $cell 0 2] {
439            #puts "v2 = $v2"
440            #puts "v1 = $v1"
441            #puts "abc = $abc"
442            set dv [expr {($v2 - $v1)*$abc}]
443            set dvs [expr {$dvs + $dv*$dv}]
444            set sumdvs [expr {$sumdvs + $dv*$dv}]
445            if {$use} {set sum($var) [expr {$sum($var) + $dv/$abc}]}
446            #puts "round and round"
447        }
448        lappend rmsout [expr {sqrt($dvs)}]
449    }
450    set rms 0
451    if {$i > 1} {set rms [expr {sqrt($sumdvs)/$i}]}
452    set neworig {}
453    foreach var {x y z} o $origin {
454        lappend neworig [expr {$o + ($sum($var)/$i)}]
455    }
456    return [list $neworig $rms $fracout $rmsout]
457}
458
459# fit a rigid body's Euler angles using least-squares
460# arguments:
461#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
462#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
463#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
464#  ortholist: list containing triplets with orthogonal coordinates
465#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
466#  fraclist: list containing triplets with fractional coordinates 
467#  origin: a list of fraction coordinates {x y z} describing the location of the
468#            origin of the orthogonal coordinates in the crystal system
469#     note that the length of ortholist, useflag and fraclist should be the same
470# Returns a list of new Euler angles
471proc FitBodyRot {Euler cell ortholist useflag fraclist origin} {
472    set xform [CalcXformMatrix $Euler  $cell]
473    set derivlist [CalcDerivMatrix $Euler  $cell]
474    set A "2 [expr 3*[llength $ortholist]] 3"
475    foreach oc $ortholist use $useflag coord $fraclist {
476        if {! $use} continue
477        foreach deriv $derivlist {
478            foreach xyz [lrange [Ortho2Xtal $deriv "0 0 0" $oc] 3 end] {
479                lappend A $xyz
480            }
481        }
482    }
483    #puts "A: [La::show $A]"
484    set y "2 [expr 3*[llength $ortholist]] 1"
485    foreach oc $ortholist use $useflag coord $fraclist {
486        if {! $use} continue
487        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
488        foreach xyz $coord XYZ $frac {
489            lappend y [expr {$XYZ - $xyz}]
490        }
491    }
492
493    set AtA [La::mmult [La::transpose $A] $A]
494    set Aty [La::mmult [La::transpose $A] $y]
495   
496    set l {}
497    #set shifts {}
498    foreach delta [lrange [La::msolve $AtA $Aty] 3 end] item $Euler {
499        #lappend shifts $delta
500        lappend l "[lindex $item 0] [expr {$delta + [lindex $item 1]}]"
501    }
502    #puts "shifts = $shifts"
503    return $l
504}
505
506# FitBody: fit a rigid body's Origin and Euler angles
507# arguments:
508#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
509#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
510#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
511#  ortholist: list containing triplets with orthogonal coordinates
512#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
513#  fraclist: list containing triplets with fractional coordinates 
514#  origin: a list of fraction coordinates {x y z} describing the location of the
515#            origin of the orthogonal coordinates in the crystal system
516#     note that the length of ortholist, useflag and fraclist should be the same
517# Returns a list containing
518#   new origin
519#   new Euler angles
520#   total rms
521#   fractional coordinates
522#   rms deviation in fractional coordinates of new Euler angles
523proc FitBody {Euler cell ortholist useflag fraclist origin "ncycle 5"} {
524    #puts "start origin = $origin"
525    foreach {
526        origin 
527        startrms
528        fracout
529        rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
530    #puts "start rms = $startrms"
531    set rmsprev $startrms
532    #puts "new origin = $origin"
533    for {set i 0} {$i < $ncycle} {incr i} {
534        set Eulerprev $Euler
535        set Euler [FitBodyRot $Euler $cell $ortholist $useflag $fraclist $origin]
536        #puts "New Euler $Euler"
537        #puts "after fit"
538        foreach {
539            origin 
540            rms
541            fracout
542            rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
543        if {$rms > (1.1 * $rmsprev) + 0.01} {
544            #puts "rms = $rms, new origin = $origin"
545            set rmsprev $rms
546        }
547    } 
548    #proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} 
549    #return "$neworig $rms $fracout $rmsout"
550    set fmt  {"%8.5f %8.5f %8.5f     %8.5f %8.5f %8.5f   %6.3f"}
551    #foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
552        #puts "[eval format $fmt $fracin $fraccalc $rmsi]"
553    #}
554    return [list $origin $Euler $rms $fracout $rmsout]
555}
556
557# zmat2coord: convert a z-matrix to a set of cartesian coordinates
558#   a z-matrix is also known as "internal coordinates" or "torsion space"
559#   (see Journal of Computational Chemistry, Vol 26, #10, p. 1063–1068, 2005 or
560#    http://www.cmbi.ru.nl/molden/zmat/zmat.html)
561# INPUT:
562#   atmlist is a list of ascii lines where each line contains
563#     lbl c1 distance c2 angle c3 torsion
564#   where each atom is computed from the previous where the new atom is:
565#     distance $distance from atom $c1 (angstrom)
566#     angle $angle from $c1--$c2 (degrees)
567#     torsion $torsion from $c1--$c2--$c3 (degrees)
568# OUTPUT:
569#  zmat2coord returns a list of atom labels and cartesian coordinates,
570#  with 4 items in each element (label, x, y, z)
571# this routine was tested against results from Babel via the web interface at
572# http://www.shodor.org/chemviz/zmatrices/babel.html and sample input at
573# http://iopenshell.usc.edu/howto/zmatrix/
574proc zmat2coord {atmlist} { 
575    set torad [expr {acos(0)/90.}]
576    set i 0
577    foreach line $atmlist {
578        incr i
579        foreach {lbl c1 dist c2 angle c3 torsion} $line {} 
580        if {$i == 1} {
581            set atm(1) [list $lbl 0 0 0] ; # 1st atom is at origin
582        } elseif {$i == 2} {
583            set dist1 $dist
584            set atm(2) [list $lbl $dist1 0 0] ; # 2nd atom is along x-axis
585        } elseif {$i == 3} {
586            # 3rd atom can be bonded to the 1st or 2nd
587            if {$c1 == 1} {
588                set atm(3) [list $lbl \
589                                [expr {$dist * cos($torad * $angle)}] \
590                                [expr {$dist * sin($torad * $angle)}] \
591                                0]
592            } else {
593                set atm(3) [list $lbl \
594                                [expr {$dist1 - $dist * cos($torad * $angle)}] \
595                                [expr {$dist * sin($torad * $angle)}] \
596                                0]
597            }
598        } else {
599            set atm($i) [concat $lbl \
600                             [ahcat "atm" $c1 $dist $c2 $angle $c3 $torsion]]
601        }
602    }
603    set coordlist {}
604    foreach key [lsort -integer [array names atm]] {
605        lappend coordlist $atm($key)
606    }
607    return $coordlist
608}
609# Compute the length of a vector
610proc vlen {a} {
611    set sum 0.0
612    foreach ai $a {
613        set sum [expr {$sum + $ai*$ai}]
614    }
615    return [expr sqrt($sum)]
616}
617# compute vector (a + z * b) and optionally normalize to length d
618proc vadd {a b d z} {
619    set c {}
620    foreach ai $a bi $b {
621        lappend c [expr {$bi + $z * $ai}]
622    }
623    set v [vlen $c]
624    if {$d != 0} {
625        set r {}
626        foreach ci $c {
627            lappend r [expr {$d * $ci / $v}]
628        }
629        return [list $v $r]
630    }
631    return [list $v $c]
632}
633# normalize a vector
634proc vnrm {x} {
635    set v [vlen $x]
636    if {abs($v) < 1e-8} {return [list 0 0 0]}
637    set y {}
638    foreach xi $x {
639        lappend y [expr {$xi / $v}]
640    }
641    return $y
642}
643# compute the coordinates for an atom that is bonded:
644#   distance $dist from atom $nc
645#   angle $bondang from $nc--$nb
646#   torsion $torsang from $nc--$nb--$na
647#   coordinates are found in array $atmarr in the calling routine
648# based on a Fortran routine provided by Peter Zavalij (Thanks Peter!)
649proc ahcat {atmarr nc dist nb bondang na torsang} {
650    upvar 1 $atmarr atm
651    set xa [lrange $atm($na) 1 3]
652    set xb [lrange $atm($nb) 1 3]
653    set xc [lrange $atm($nc) 1 3]
654    set torad [expr {acos(0)/90.}]
655    # x = unit Vector A-B
656    foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. -1.] 1] {}
657    # y = unit Vector C-B
658    set y [lindex [vadd $xb $xc 1. -1.] 1]
659    foreach {y1 y2 y3} $y {}
660    set z1 [expr {$x2*$y3 - $x3*$y2}]
661    set z2 [expr {$x3*$y1 - $x1*$y3}]
662    set z3 [expr {$x1*$y2 - $x2*$y1}]
663    set z [vnrm [list $z1 $z2 $z3]]
664    set q1 [expr {$y2*$z3 - $y3*$z2}]
665    set q2 [expr {$y3*$z1 - $y1*$z3}]
666    set q3 [expr {$y1*$z2 - $y2*$z1}]
667    set q [vnrm [list $q1 $q2 $q3]]
668    set th [expr {$bondang * $torad}]
669    set ph [expr {-1. * $torsang * $torad}]
670    set cth [expr {cos($th)}]
671    set sth [expr {sin($th)}]
672    set cph [expr {cos($ph)}]
673    set sph [expr {sin($ph)}]
674    set xh {}
675    foreach xci $xc xi $q zi $z yi $y {
676        lappend xh [expr {
677                          $xci +
678                          $dist*($sth*($cph*$xi + $sph*$zi)-$cth*$yi)
679                      }]
680    }
681    return $xh
682}
683 
684# RB2cart: convert the rigid body representation reported as the 2nd element
685# in ReadRigidBody into cartesian coordinates
686#   rblist: a list containing an element for each scaling factor
687# in each element there are four items:
688#    a multiplier value for the rigid body coordinates
689#    a damping value (0-9) for the refinement of the multiplier (not used)
690#    a variable number if the multiplier will be refined (not used)
691#    a list of cartesian coordinates coordinates
692# each cartesian coordinate contains 4 items: x,y,z and a label
693# returns a list of coordinate triplets
694proc RB2cart {rblist} {
695    foreach item $rblist {
696        foreach {mult damp ref coords} $item {}
697        set i 0
698        foreach xyz $coords {
699            foreach {x y z} [lrange $xyz 0 2] {}
700            foreach val [lrange $xyz 0 2] var {X Y Z} {
701                if {[array names $var $i] == ""} {
702                    set ${var}($i) [expr {$mult * $val}]
703                } else {
704                    set ${var}($i) [expr {[set ${var}($i)] + $mult * $val}]
705                }
706            }
707            incr i
708        }
709    }
710    set out ""
711    foreach i [lsort -integer [array names X]] {
712        lappend out [list $X($i) $Y($i) $Z($i)]
713    }
714    return $out
715}
716
717# get the name of the DRAWxtl application, if installed
718proc GetDRAWxtlApp {} {
719    # is DRAWxtl installed?
720    set app {}
721    if {![catch {set fp [open [file join $::env(HOME) .drawxtlrc] r]}]} {
722        # line 12 is name of executable
723        set i 0
724        while {$i < 12} {
725            incr i
726            gets $fp appname
727        }
728        close $fp
729        set app [auto_execok $appname]
730    }
731    if {$app != ""} {
732        return $appname
733    }
734    return ""
735}
736
737# DRAWxtlPlotOrtho: plot orthogonal coordinates in DRAWxtl
738# input:
739#  filename: file name for the .str file to create
740#  title: string for title in .str file
741#  coords: cartesian coordinates
742#  bondlist: list of bonds to draw as min, max length (A) and
743#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
744proc DRAWxtlPlotOrtho {filename title coords bondlist} {
745    foreach {xmin ymin zmin} {"" "" ""} {}
746    foreach {xmax ymax zmax} {"" "" ""} {}
747    foreach xyz $coords {
748        foreach {x y z} $xyz {}
749        foreach s {x y z} {
750            foreach t {min max} {
751                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
752            } 
753            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
754            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
755        }
756    }
757    #puts "$xmin $xmax $ymin $ymax $zmin $zmax"
758    set max $xmin
759    foreach val "$xmin $xmax $ymin $ymax $zmin $zmax" {
760        if {$max < abs($val)} {set max $val}
761    }
762    set scale [expr {4.*$max}]
763    set a 10.
764    lappend range [expr -0.01+($xmin/$scale)] [expr 0.01+($xmax/$scale)] \
765        [expr -0.01+($ymin/$scale)] [expr 0.01+($ymax/$scale)] \
766        [expr -0.01+($zmin/$scale)] [expr 0.01+($zmax/$scale)]
767    set fp [open $filename w]
768    puts $fp "title $title"
769    puts $fp "box  0.000 Black"
770    puts $fp "background White"
771    #puts $fp "nolabels"
772    puts $fp "cell $a $a $a 90 90 90"
773    puts $fp "spgr P 1"
774    puts $fp "pack $range"
775    set i 0
776    foreach xyz $coords {
777        foreach {x y z} $xyz {}
778        incr i
779        puts $fp "atom c $i [expr {$x/$scale}] [expr {$y/$scale}] [expr {$z/$scale}]"
780        puts $fp "labeltext [expr {0.02 + $x/$scale}] [expr {0.01 + $y/$scale}] [expr {0.01 + $z/$scale}] $i"
781    }
782    puts $fp "sphere c  [expr 0.100*($a/$scale)] Red"
783    puts $fp "finish   0.70   0.30   0.08   0.01"
784    foreach bondpair $bondlist {
785        foreach {b1 b2 color} $bondpair {}
786        if {$color == ""} {set color Red}
787        puts $fp "bond c c [expr {0.01*$a/$scale}] [expr {$b1*$a/$scale}] [expr {$b2*$a/$scale}] $color"
788    }
789    puts $fp "frame"
790    set range {}
791    lappend range -0.01 [expr 0.01+(0.1*$a/$scale)] \
792        -0.01 [expr 0.01+(0.1*$a/$scale)] \
793        -0.01 [expr 0.01+(0.1*$a/$scale)]
794    puts $fp "cell $a $a $a 90 90 90"
795    puts $fp "spgr P 1"
796    puts $fp "pack $range"
797    puts $fp "atom o 1 0 0 0"
798    puts $fp "atom o 2 [expr {0.1*$a/$scale}] 0 0"
799    puts $fp "atom o 3 0 [expr {0.1*$a/$scale}] 0"
800    puts $fp "atom o 4 0 0 [expr {0.1*$a/$scale}]"
801    puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
802    puts $fp "labelscale 0.5"
803    puts $fp "labeltext [expr {0.11*$a/$scale}] 0 0 x"
804    puts $fp "labeltext 0 [expr {0.11*$a/$scale}] 0 y"
805    puts $fp "labeltext 0 0 [expr {0.11*$a/$scale}] z"
806    puts $fp "sphere o [expr {0.02*$a/$scale}] Blue"
807    puts $fp "origin   .0 .0 .0"
808    puts $fp "end"
809    close $fp
810}
811
812# PlotRBtype: plot a rigid body in DRAWxtl
813# input:
814#  rbtype: # of rigid body
815#  bondlist: list of bonds to draw as min, max length (A) and
816#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
817#  file: file name for the .str file to create
818proc PlotRBtype {rbtype "bondlist {}" "file {}"} {
819    set app [GetDRAWxtlApp]
820    if {$app == ""} {
821        MyMessageBox -parent . -title "No DRAWxtl" \
822                -message "Sorry, DRAWxtl is not installed" \
823                -icon warning
824        return
825    }
826    if {$::tcl_platform(platform) == "windows" && $file == ""} {
827        set file [file join [pwd] rbplot.str]
828    } else {
829        set file "/tmp/rbplot.str"
830    }
831    set coords [RB2cart [lindex [ReadRigidBody $rbtype] 1]]
832    DRAWxtlPlotOrtho $file "" $coords $bondlist
833    if {$app != ""} {exec $app $file &}
834}
835
836# PlotRBcoords: plot orthogonal coordinates in DRAWxtl
837# input:
838#  coords: cartesian coordinates
839#  bondlist: list of bonds to draw as min, max length (A) and
840#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
841#  file: file name for the .str file to create
842proc PlotRBcoords {coords "bondlist {}" "file {}"} {
843    set app [GetDRAWxtlApp]
844    if {$app == ""} {
845        MyMessageBox -parent . -title "No DRAWxtl" \
846                -message "Sorry, DRAWxtl is not installed" \
847                -icon warning
848        return
849    }
850    if {$::tcl_platform(platform) == "windows" && $file == ""} {
851        set file [file join [pwd] rbplot.str]
852    } else {
853        set file "/tmp/rbplot.str"
854    }
855    DRAWxtlPlotOrtho $file "" $coords $bondlist
856    if {$app != ""} {exec $app $file &}
857}
858
859# DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed
860# on a structure read from a phase
861# input:
862#  RBcoords: fractional coordinates for rigid body
863#  phase:# of phase to plot
864#  firstatom: seq # of 1st atom in structure to be mapped to rigid body
865#  allatoms: 0 to plot only atoms in phase that are in the rigid body,
866#      otherwise plot all atoms
867#  bondlist: list of bonds to draw for the phase as min, max length (A) and
868#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
869#  rbbondlist: list of bonds to draw for the phase as min, max length (A) and
870#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
871#  file: optional file name for the .str file to create
872proc DRAWxtlPlotRBFit {RBcoords phase firstatom "allatoms 0" \
873                           "bondlist {}" "rbbondlist {}" "file {}"} {
874    set natom [llength $RBcoords]
875    set app [GetDRAWxtlApp]
876    if {$app == ""} {
877        MyMessageBox -parent . -title "No DRAWxtl" \
878                -message "Sorry, DRAWxtl is not installed" \
879                -icon warning
880        return
881    }
882    if {$::tcl_platform(platform) == "windows" && $file == ""} {
883        set file [file join [pwd] rbplot.str]
884    } else {
885        set file "/tmp/rbfit.str"
886    }
887
888    # get rigid body coordinate range
889    foreach {xmin ymin zmin} {"" "" ""} {}
890    foreach {xmax ymax zmax} {"" "" ""} {}
891    foreach xyz $RBcoords {
892        foreach {x y z} $xyz {}
893        foreach s {x y z} {
894            foreach t {min max} {
895                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
896            } 
897            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
898            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
899        }
900    }
901    set rbrange {}
902    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
903                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
904                     [expr -0.01+$zmin] [expr 0.01+$zmax] ] {
905        append rbrange [format " %8.4f" $val]
906    }
907    set rbcenter [list [expr {($xmin+$xmax)/2}] \
908                      [expr {($ymin+$ymax)/2}] \
909                      [expr {($zmin+$zmax)/2}] ]
910    # get matching atoms coordinate range
911    set firstind [lsearch $::expmap(atomlist_$phase) $firstatom]
912    set matchedatomlist [lrange \
913                      [lrange $::expmap(atomlist_$phase) $firstind end] \
914                      0 [expr {$natom-1}]]
915    foreach atom $matchedatomlist {
916        foreach s {x y z} {
917            set $s [atominfo $phase $atom $s]
918            foreach t {min max} {
919                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
920            } 
921            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
922            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
923        }
924    }
925    # expand to cover at least one unit cell
926    foreach var {xmin ymin zmin} { 
927        if {[set $var] > 0.0} {set $var 0.0}
928    }
929    foreach var {xmax ymax zmax} { 
930        if {[set $var] < 1.} {set $var 1.}
931    }
932    set range {}
933    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
934                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
935                     [expr -0.01+$zmin] [expr 0.01+$zmax]] {
936        append range [format " %8.4f" $val]
937    }
938
939    set fp [open $file w]
940    puts $fp "title structure/rigid-body fit plot"
941    # plot the structure
942    puts -nonewline $fp "cell"
943    foreach p {a b c alpha beta gamma} {
944        puts -nonewline $fp " [phaseinfo $phase $p]"
945    }
946    puts $fp ""
947    puts $fp "spgp [phaseinfo $phase spacegroup]"
948    puts $fp "pack $range"
949    if {$allatoms != 0} {
950        set atoms $::expmap(atomlist_$phase)
951    } else {
952        set firstind [lsearch $::expmap(atomlist_$phase) $firstatom]
953        set atoms [lrange \
954                       [lrange $::expmap(atomlist_$phase) $firstind end] \
955                       0 [expr {$natom-1}]]
956    }
957
958    # set origin at center of rigid body
959    puts $fp "origin  $rbcenter"
960    # now loop over atoms
961    foreach atom $atoms {
962        set type [atominfo $phase $atom type]
963        set typelist($type) 1
964        set xyz ""
965        foreach v {x y z} {
966            append xyz "[atominfo $phase $atom $v] "
967        }
968        puts $fp "atom $type $atom $xyz"
969        if {[lsearch $matchedatomlist $atom] != -1} {
970            puts $fp "labeltext $xyz $atom"
971        }
972               
973        set uiso [atominfo $phase $atom Uiso]
974        # are there anisotropic atoms? If so convert them to Uequiv
975        if {[atominfo $phase $atom temptype] == "A"} {
976            puts -nonewline $fp "Uij [atominfo $phase $atom type] $atom "
977            foreach v {U11 U22 U33 U12 U13 U23} {
978                puts -nonewline $fp "[atominfo $phase $atom $v] "
979            }
980            puts $fp ""
981        }
982    }
983
984    foreach type [array names typelist] color {Green Blue Magenta Cyan} {
985        if {$type == ""} break
986        puts $fp "sphere $type 0.1 $color"
987    }
988    foreach type [array names typelist] color1 {Green Blue Magenta Cyan} {
989        foreach bondpair $bondlist {
990            foreach {b1 b2 color} $bondpair {}
991            if {$color == ""} {set color $color1}
992            puts $fp "bond $type $type 0.02 $b1 $b2 $color"
993        }
994        foreach type1 [array names typelist] {
995            if {$type1 == $type} break
996            foreach bondpair $bondlist {
997                foreach {b1 b2 color} $bondpair {}
998                if {$color == ""} {set color $color1}
999                puts $fp "bond $type $type1 0.02 $b1 $b2 $color"
1000            }
1001        }
1002    }
1003    # plot the rigid body
1004    puts $fp "frame"
1005    puts -nonewline $fp "cell"
1006    foreach p {a b c alpha beta gamma} {
1007        puts -nonewline $fp " [phaseinfo $phase $p]"
1008    }
1009    puts $fp ""
1010    puts $fp "background White"
1011    #puts $fp "nolabels"
1012    puts $fp "labelscale 0.5"
1013    puts $fp "spgr P 1"
1014    puts $fp "pack $rbrange"
1015    set i 0
1016    foreach xyz $RBcoords {
1017        foreach {x y z} $xyz {}
1018        incr i
1019        puts $fp "atom c $i $x $y $z"
1020        puts $fp "labeltext $x $y $z r$i"
1021    }
1022    foreach bondpair $rbbondlist {
1023        foreach {b1 b2 color} $bondpair {}
1024        if {$color == ""} {set color Red}
1025        puts $fp "bond c c 0.02 $b1 $b2 $color"
1026    }
1027
1028    puts $fp "sphere c 0.05 Red"
1029    puts $fp "finish   0.70   0.30   0.08   0.01"
1030    puts $fp "end"
1031
1032    #puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
1033    close $fp
1034    MyMessageBox -parent . -title "Info" \
1035                -message "Note that the phase is drawn in green, blue, cyan & magenta and the rigid body in red."
1036    if {$app != ""} {exec $app $file &}
1037}
1038
1039
1040#AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
1041#puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"]
1042#puts [GetRB 1 4 8 "1" "X 1 2" "Z 3 4"]
1043#MapRigidBody 1 1 7 ".11 .22 .33" "11 12 13"
1044
1045
1046#AddRigidBody {1} { {
1047#    {1 1 1 o} {-1 1 1 o} {1 -1 1 o} {-1 -1 1 o}
1048#    {1 1 -1 o} {-1 1 -1 o} {1 -1 -1 o} {-1 -1 -1 o}
1049#} }
1050#set n [MapRigidBody 1 1 1 ".2 .3 .4" "13 17 19"]
1051#puts "body $n created"
1052#incr expgui(changed)
1053#RunRecalcRBCoords
1054#puts "press Enter to continue"
1055#gets stdin line
1056#MapRigidBody 1 1 $n ".5 .5 .5" "0 0 0"
1057#incr expgui(changed)
1058#RunRecalcRBCoords
1059
1060#puts "Test FitBody"
1061set fraclist {
1062    { 0.5483305238484277 0.4887545024531055 0.6167996784631056 }
1063    { 0.1036801409356145 0.5954016321779562 0.5129448102437683 }
1064    { 0.26404665760133855 0.09455414439078394 0.612655365147539 }
1065    { -0.18060372531147473 0.20120127411563465 0.5088004969282018 }
1066    { 0.5806037253114747 0.3987987258843653 0.2911995030717982 }
1067    { 0.13595334239866147 0.5054458556092161 0.18734463485246095 }
1068    { 0.2963198590643855 0.004598367822043814 0.2870551897562318 }
1069    { -0.1483305238484277 0.1112454975468945 0.1832003215368945 }
1070}
1071set ortholist {
1072    {1 1 1} 
1073    {-1 1 1}
1074    {      1.000000   -1.000000    1.000000}
1075    {     -1.000000   -1.000000    1.000000}
1076    {      1.000000    1.000000   -1.000000}
1077    {     -1.000000    1.000000   -1.000000}
1078    {      1.000000   -1.000000   -1.000000}
1079    {     -1.000000   -1.000000   -1.000000} 
1080}
1081# test code, generates DRAWxtl imput file from orthogonal coordinate list
1082# with bonds of ~2, 2.8 and 3.4 A
1083#DRAWxtlPlotOrtho test4.str "test file" $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
1084
1085# test code, plots rigid body type #2 with bonds drawn at ~1.3 & 2 A
1086#PlotRBtype 2 {{1.9 2.1} {1.28 1.32}}
1087
1088# test code, plots rigid body coords in ortholist with bonds @  ~2, 2.8 and 3.4 A
1089#PlotRBcoords $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
1090
1091
1092set useflag {1 1 1 1 1 1 1 1}
1093set cell {4.  5.  6.  95.  100.  105.}
1094#set origin ".20 .30 .40"
1095set origin ".0 .0 .0"
1096#set Euler  {{1 13} {2 17} {3 19}}
1097#set Euler  {{1 0} {2 180} {3 0}}
1098set Euler  {{1 0} {2 0} {3 0}}
1099
1100#puts [La::show $xform]
1101#puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
1102
1103
1104# test zmat2coord
1105set atmlist {
1106    {C1 0 0.0 0 0.0 0 0.0}
1107    {O2 1 1.20 0 0.0 0 0.0}
1108    {H3 1 1.10 2 120.0 0 0.0}
1109    {C4 1 1.50 2 120.0 3 180.0}
1110    {H5 4 1.10 1 110.0 2 0.00}
1111    {H6 4 1.10 1 110.0 2 120.0}
1112    {H7 4 1.10 1 110.0 2 -120.0}
1113}
1114#  C        0.00000        0.00000        0.00000
1115#  O        1.20000        0.00000        0.00000
1116#  H       -0.55000        0.95263        0.00000
1117#  C       -0.75000       -1.29904       -0.00000
1118#  H       -0.04293       -2.14169       -0.00000
1119#  H       -1.38570       -1.36644        0.89518
1120#  H       -1.38570       -1.36644       -0.89518
1121# set coordlist [zmat2coord $atmlist]
1122 set i 0
1123# puts "\nZmatrix in"
1124# foreach line $atmlist {
1125#     incr i
1126#     puts "$i) $line"
1127# }
1128# puts "Cartesian out"
1129# foreach line $coordlist {
1130#     puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
1131# }
1132
1133# AddRigidBody {1 0.75} {
1134#     {
1135#       {1 1 1 c}
1136#       {-1 1 1 c}
1137#       {      1.000000   -1.000000    1.000000 c}
1138#       {     -1.000000   -1.000000    1.000000 c}
1139#       {      1.000000    1.000000   -1.000000 c}
1140#       {     -1.000000    1.000000   -1.000000 c}
1141#       {      1.000000   -1.000000   -1.000000 c}
1142#       {     -1.000000   -1.000000   -1.000000 c}
1143#       {1 1 1 h}
1144#       {1 -1 -1 h}
1145#       {-1 1 -1 h}
1146#       {-1 -1 1 h}
1147#     } {
1148#       {0 0 0 c }
1149#       {0 0 0 c}
1150#       {0 0 0 c}
1151#       {0 0 0 c}
1152#       {0 0 0 c}
1153#       {0 0 0 c}
1154#       {0 0 0 c}
1155#       {0 0 0 c}
1156#       {1 1 1 h}
1157#       {1 -1 -1 h}
1158#       {-1 1 -1 h}
1159#       {-1 -1 1 h}
1160#     }
1161# }
1162# MapRigidBody 2 2 1 {0 0 0} {10 15 20}
Note: See TracBrowser for help on using the repository browser.