source: branches/sandbox/rb.tcl @ 1172

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

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