source: branches/sandbox/rb.tcl @ 1138

Last change on this file since 1138 was 1138, checked in by toby, 14 years ago

fixes for SetRigidBody? and SetRigidBodyVar?
File size: 63.3 KB
Line 
1# test code (package already loaded in expgui)
2lappend auto_path [file dirname [info script]]
3package require La
4
5#============================================================================
6# rigid body EXP editing routines (to move into readexp.tcl)
7# RigidBodyList -- returns a list of the defined rigid body types
8# SetRigidBodyVar -- set variables and damping for rigid body type multipliers
9# ReadRigidBody  -- # of times a body is mapped, scaling factors, var #s & coordinates
10# RigidBodyMappingList - return a list instances where a RB is mapped in phase
11# RigidBodyEnableTLS -- Enable or Disable TLS use for a rigid body mapping
12# RigidBodySetTLS  -- change the TLS values for a rigid body mapping
13# RigidBodySetDamp -- change the damping values for a rigid body mapping
14# RigidBodyVary    -- set refinement variable numbers for a rigid body mapping
15# RigidBodyTLSVary -- set TLS refinement variable nums for a rigid body mapping
16# AddRigidBody -- defines a new rigid body type
17# DeleteRigidBody -- remove a rigid body definition
18# ReplaceRigidBody -- replaces a previous rigid body type
19# ReadRigidBodyMapping  -- get parameters for a rigid body mapping
20# MapRigidBody -- map a rigid body type into a phase
21# EditRigidBodyMapping -- change the parameters in a rigid body mapping
22# UnMapRigidBody --remove a rigid body constraint by removing a RB "instance"
23#----- note that these older routines should not be used ------
24# RigidBodyCount -- returns the number of defined rigid bodies (body types)
25#    use RigidBodyList instead
26# RigidBodyMappingCount -- # of times a rigid body is mapped in phase
27#    use RigidBodyMappingList instead
28#============================================================================
29# returns the number of defined rigid bodies
30proc RigidBodyCount {} {
31    set n [string trim [readexp "RGBD  NRBDS"]]
32    if {$n == ""} {
33        set n 0
34    }
35    return $n
36}
37
38# returns a list of the defined rigid body types
39proc RigidBodyList {} {
40    set n [string trim [readexp "RGBD  NRBDS"]]
41    if {$n == ""} {
42        set n 0
43    }
44    set rblist {}
45    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
46        set value $rbnum
47        validint value 2
48        set key "RGBD${value}"
49        if {[existsexp "$key NATR "]} {
50            lappend rblist $rbnum
51        }
52        if {[llength $rblist] == $n} break
53    }
54    return $rblist
55}
56
57# ReadRigidBody provides all information associated with a rigid body type
58#  rbnum is the rigid body type number
59# it returns two items:
60#   the number of times the rigid body is mapped
61#   a list containing an element for each scaling factor in rigid body #rbnum.
62# in each element there are four items:
63#    a multiplier value for the rigid body coordinates
64#    a damping value (0-9) for the refinement of the multiplier
65#    a variable number if the multiplier will be refined
66#    a list of cartesian coordinates coordinates
67# each cartesian coordinate contains 4 items: x,y,z and a label
68#  note that the label is present only when the RB is created in EXPGUI and is
69#  not used in GSAS.
70proc ReadRigidBody {rbnum} {
71    if {[lsearch [RigidBodyList] $rbnum] == -1} {
72        return ""
73    }
74    set value $rbnum
75    validint value 2
76    set key "RGBD${value}"
77    set n [string trim [string range [readexp "$key NATR"] 0 4]]
78    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
79    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
80    set out {}
81    for {set i 1} {$i <= $nmult} {incr i} {
82        set line [readexp "${key}${i}PARM"]
83        set mult [string trim [string range $line 0 9]]
84        set var [string trim [string range $line 10 14]]
85        set damp [string trim [string range $line 15 19]]
86        set coordlist {}
87        for {set j 1} {$j <= $n} {incr j} {
88            set value $j
89            validint value 3
90            set line [readexp "${key}${i}SC$value"]
91            set x [string trim [string range $line 0 9]]
92            set y [string trim [string range $line 10 19]]
93            set z [string trim [string range $line 20 29]]
94            set lbl [string trim [string range $line 30 39]]
95            lappend coordlist [list $x $y $z $lbl]
96        }
97        lappend out [list $mult $damp $var $coordlist]
98    }
99    return [list $used $out]
100}
101
102# SetRigidBodyVar
103#   rbnum is the rigid body type number
104#   varnumlist is a list of variable numbers
105#      note that if this list is shorter than the number of actual multipliers
106#      for the body, the unspecified variable will not be changed
107#   damplist   is a list of damping values (0-9)
108#      note that if the damplist is shorter than the number of actual multipliers
109#      the unspecified values are not changed
110#  SetRigidBodVar 2 {1 2 3} {}
111#       will vary the (first 3) translations in body #3 and will not change the
112#       damping values
113#  SetRigidBodVar 3 {} {0 0 0}
114#       will not change variable settings but will change the (first 3) damping values
115#  SetRigidBodVar 4 {11 11} {8 8}
116#      changes both variable numbers and damping at the same time
117# Nothing is returned
118proc SetRigidBodyVar {rbnum varnumlist damplist} {
119    if {[lsearch [RigidBodyList] $rbnum] == -1} {
120        return ""
121    }
122    set value $rbnum
123    validint value 2
124    set key "RGBD${value}"
125    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
126    for {set i 1} {$i <= $nmult} {incr i} {
127        set j $i
128        incr j -1
129        set var [lindex $varnumlist $j]
130        if {$var != ""} {
131            validint var 5
132            setexp "${key}${i}PARM" $var 11 15
133        }
134        set damp [lindex $damplist $j]
135        if {$damp != ""} {
136            if {$damp > 9} {set damp 9}
137            if {$damp < 0} {set damp 0}
138            validint damp 5
139            setexp "${key}${i}PARM" $damp 16 20
140        }
141    }
142}
143
144
145# return the number of times rigid body $bodytyp is mapped in phase $phase
146proc RigidBodyMappingCount {phase bodytyp} {
147    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
148    if {! [existsexp "$key  NBDS"]} {return 0}
149    set n [string trim [readexp "$key  NBDS"]]
150    if {$n == ""} {
151        set n 0
152    }
153    return $n
154}
155# return a list of the instances where rigid body $bodytyp is mapped in phase $phase
156proc RigidBodyMappingList {phase bodytyp} {
157    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
158    if {! [existsexp "$key  NBDS"]} {return {}}
159    set n [string trim [readexp "$key  NBDS"]]
160    if {$n == ""} {
161        set n 0
162    }
163    set rblist {}
164    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
165        set value $rbnum
166        validint value 2
167        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
168        if {[existsexp "$key  NDA"]} {
169            lappend rblist $rbnum
170        }
171        if {[llength $rblist] == $n} break
172    }
173    return $rblist
174}
175
176
177
178# reads rigid body mapping parameters for phase ($phase), body type # ($bodytyp) and instance # ($num)
179# returns a list of items (most lists) as follows:
180#   1) sequence # of first atom in body
181#   2) origin of body in fractional coordinates (3 elements)
182#   3) Euler angles as 6 pairs of numbers (see below)
183#   4) variable numbers for the 9 position variables (origin followed by rotations)
184#   5) damping vals for the 9 position variables (origin followed by rotations)
185#   6) the TLS values, in order below (empty list if TLS is not in use)
186#   7) the variable numbers for each TLS values, in order below (or empty)
187#   8) three damping values for the T, L and S terms.
188# returns an empty list if no such body exists.
189#
190# Euler angles are a list of axes and angles to rotate:
191#   { {axis1 angle1} {axis2 angle2} ...}
192# where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
193#
194# The 20 TLS terms are ordered:
195#    T11, T22, T33, T12, T13, T23
196#    L11, L22, L33, L12, L13, L23
197#    S12, S13, S21, S23, S31, S32, SAA, SBB
198#
199proc ReadRigidBodyMapping {phase bodytyp num} {
200    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
201        return ""
202    }
203    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
204    set first [string trim [string range [readexp "$key  NDA"] 0 4]]
205    set line [readexp "$key BDFL"]
206    set varlist {}
207    set damplist {}
208    foreach i {0 1 2 3 4 5 6 7 8} {
209        lappend varlist [string trim [string range $line [expr {5*$i}] [expr {4 + 5*$i}] ]]
210        lappend damplist [string trim [string range $line [expr {45 + $i}] [expr {45 + $i}] ]]
211    }
212    set TLSdamplist {}
213    foreach i {54 55 56} {
214        lappend TLSdamplist [string trim [string range $line $i $i ]]
215    }
216    set line [readexp "${key} BDLC"]
217    set x [string trim [string range $line 0 9]]
218    set y [string trim [string range $line 10 19]]
219    set z [string trim [string range $line 20 29]]
220    set origin [list $x $y $z]
221    set line [readexp "${key} BDOR"]
222    set rotations {}
223    foreach i {0 10 20 30 40 50} {
224        set angle [string trim [string range $line $i [expr {$i+7}]]]
225        set axis [string trim [string range $line [expr {$i+8}] [expr {$i+9}]]]
226        lappend rotations [list $angle $axis]
227    }
228    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
229    set tlsvars {}
230    set tlsvals {}
231    if {$TLS != 0} {
232        set line [readexp "${key}TLSF1"]
233        for {set j 0} {$j < 20} {incr j} {
234            set var [string trim [string range $line [expr {3*$j}] [expr {3*$j+2}]]]
235            if {$var == ""} {set var 0}
236            lappend tlsvars $var
237        }
238        for {set j 0} {$j < 20} {incr j} {
239            set i 0
240            if {$j == 0} {
241                set i 1
242            } elseif {$j == 8} {
243                set i 2
244            } elseif {$j == 16} {
245                set i 3
246            }
247            if {$i != 0} {
248                set line [readexp "${key}TLSP$i"]
249                set i 0
250                set j1 0
251                set j2 7
252            } else {
253                incr j1 8
254                incr j2 8
255            }
256            set val [string trim [string range $line $j1 $j2]]
257            if {$val == ""} {set val 0}
258            lappend tlsvals $val
259        }
260    }
261    return [list $first $origin $rotations $varlist $damplist $tlsvals $tlsvars $TLSdamplist]
262}
263
264# Control TLS representation for phase, body # and instance number of a Rigid body mapping
265#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
266# Enable TLS use if TLS is non-zero (true). Disable if zero
267proc RigidBodyEnableTLS {phase bodytyp num TLS} {
268    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
269        return ""
270    }
271    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
272    if {$TLS} {
273        setexp "${key} LSTF" [format "%5d" 1] 1 5
274        if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
275        if {![existsexp "${key}TLSP1"]} {
276            makeexprec "${key}TLSP1"
277            set str {}
278            foreach v {.01 .01 .01 0 0 0 0 0} d {4 4 4 4 4 4 2 2} {
279                validreal v 8 $d
280                append str $v
281            }
282            setexp "${key}TLSP1" $str 1 64
283        }
284        if {![existsexp "${key}TLSP2"]} {
285            makeexprec "${key}TLSP2"
286            set str {}
287            set v 0
288            foreach d {2 2 2 2 4 4 4 4} {
289                validreal v 8 $d
290                append str $v
291            }
292            setexp "${key}TLSP2" $str 1 64
293        }
294        if {![existsexp "${key}TLSP3"]} {
295            makeexprec "${key}TLSP3"
296            set str {}
297            set v 0
298            foreach d {4 4 4 4} {
299                validreal v 8 $d
300                append str $v
301            }
302            setexp "${key}TLSP3" $str 1 64
303        }
304    } else {
305        setexp "${key} LSTF" [format "%5d" 0] 1 5
306    }
307    return 1
308}
309
310# Control the TLS values for Rigid body mapping for mapping with
311#    phase ($phase), body type # ($bodytyp) and instance # ($num)
312# set the 20 TLS values to the values in TLSvals
313# There must be exactly 20 TLS terms, which are ordered:
314#    T11, T22, T33, T12, T13, T23
315#    L11, L22, L33, L12, L13, L23
316#    S12, S13, S21, S23, S31, S32, SAA, SBB
317proc RigidBodySetTLS {phase bodytyp num TLSvals} {
318    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
319        return ""
320    }
321    if {[llength $TLSvals] != 20} {return ""}
322    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
323    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
324    if {$TLS == 0} {return ""}
325    if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
326    foreach n {1 2 3} {
327        if {![existsexp "${key}TLSP$n"]} {makeexprec "${key}TLSP$n"}
328    }
329    set str {}
330    set n 1
331    set i 0
332    foreach v $TLSvals d {4 4 4 4 4 4 2 2 2 2 2 2 4 4 4 4 4 4 4 4} {
333        incr i
334        validreal v 8 $d
335        append str $v
336        if {$i == 8} {
337            set i 0
338            setexp "${key}TLSP$n" $str 1 64
339            incr n
340            set str {}
341        }
342    }
343    setexp "${key}TLSP$n" $str 1 64
344    return 1
345}
346
347# set damping values for a Rigid body mapping
348#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
349# there must be 9 damping values in RBdamp for the 9 position variables (origin followed by rotations)
350# Use of TLSdamp is optional, but to be used, TLS representation must be enabled and there must be
351# three damping terms (for all T terms; for all L terms and for all S terms)
352proc RigidBodySetDamp {phase bodytyp num RBdamp "TLSdamp {}"} {
353    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
354        return ""
355    }
356    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
357    if {[llength $RBdamp] != 9} {return ""}
358    set str {}
359    foreach v $RBdamp {
360        if {[validint v 1] != 1} {set v " "}
361        append str $v
362    }
363    setexp "$key BDFL" $str 46 9
364    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
365    if {$TLS != 0 &&  [llength $TLSdamp] == 3} {
366        set str {}
367        foreach v $TLSdamp {
368        if {[validint v 1] != 1} {set v " "}
369            append str $v
370        }
371        setexp "$key BDFL" $str 55 3
372    }
373    return 1
374}
375
376# set refinement variable numbers for a Rigid body mapping
377#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
378# there must be 9 variable values in RBvar for the 9 position variables (origin followed by rotations)
379# note that the variable values should be unique integers
380proc RigidBodyVary {phase bodytyp num RBvar} {
381    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
382        return ""
383    }
384    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
385    if {[llength $RBvar] != 9} {return ""}
386    set str {}
387    foreach v $RBvar {
388        if {[validint v 5] != 1} {set v " "}
389        append str $v
390    }
391    setexp "$key BDFL" $str 1 45   
392}
393
394# set TLS refinement variable numbers for a Rigid body mapping
395#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
396# there must be 20 variable values in TLSvar for the 20 parameters:
397#    T11, T22, T33, T12, T13, T23
398#    L11, L22, L33, L12, L13, L23
399#    S12, S13, S21, S23, S31, S32, SAA, SBB
400# note that the variable values should be unique integers
401proc RigidBodyTLSVary {phase bodytyp num TLSvar} {
402    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
403        return ""
404    }
405    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
406    if {[llength $TLSvar] != 20} {return ""}
407    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
408    if {$TLS == 0} {return ""}
409    set str {}
410    foreach v $TLSvar {
411        if {[validint v 3] != 1} {set v " "}
412        append str $v
413    }
414    setexp "${key}TLSF1" $str 1 60
415
416# AddRigidBody: add a new rigid body definition into the .EXP file
417# arguments are:
418#   multlist: defines a list of multipliers for each set of coordinates. In the
419#             simplest case this will be {1}
420#   coordlist: a nested list of coordinates such as { { {0 0 0} {.1 .1 .1} {.2 .2 .2} } }
421# note that when the length of multlist > 1 then coordlist must have the same length.
422# for input where
423#     multlist = {s1 s2} and
424#     coordlist = { { {0 0 0} {1 1 0} {.0 .0 .0} ...}
425#                     {0 0 0} {1 1 0} {2 1 2} ...}
426#                 }   
427# the cartesian coordinates are defined from the input as
428#    atom 1 = s1 * (0,0,0) + s2*(0,0,0) [ = (0,0,0)]
429#    atom 2 = s1 * (1,1,0) + s2*(1,1,0) [ = (s1+s2) * (1,1,0)]
430#    atom 3 = s1 * (0,0,0) + s2*(2,1,2) [ = s2 * (2,1,2)]
431#    ...
432# Returns the number of the rigid body that has been created
433proc AddRigidBody {multlist coordlist} {
434    # find the first unused body #
435    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
436        set value $rbnum
437        validint value 2
438        set key "RGBD${value}"
439        if {! [existsexp "$key NATR "]} {break}
440    }
441    # did we go too far?
442    if {$rbnum == 16} {return ""} 
443    # increment the RB counter
444    set n [string trim [readexp "RGBD  NRBDS"]]
445    if {$n == ""} {
446        makeexprec "RGBD  NRBDS"
447        set n 0
448    }
449    incr n
450    validint n 5
451    setexp "RGBD  NRBDS" $n 1 5
452    SetRigidBody $rbnum $multlist $coordlist
453    return $rbnum
454}
455
456# DeleteRigidBody: remove a rigid body definition from the .EXP file
457# The body may not be mapped. I am not sure if GSAS allows more than 9 bodies,
458# but if it does, the simplifed approach used here will fail, so this
459# is not allowed.
460# Input:
461#   Rigid body number
462# Returns:
463#   1 on success
464#   -1 if the body number is 11 or greater
465#   -2 if the body is mapped
466#   -3 if the body is not defined
467proc DeleteRigidBody {rbnum} {
468    # can't delete bodies with numbers higher than 10, since the key prefix
469    # (RGBD11... will overlap with rigid body instance records, which would be
470    # deleted below
471    if {$rbnum > 10} {
472        return -1
473    }
474    set value $rbnum
475    validint value 2
476    set key "RGBD${value}"
477    if {![existsexp "$key NATR "]} {
478        return -2
479    }
480    # make sure the body is not mapped
481    if {[string trim [string range [readexp "$key NBDS"] 0 4]] != 0} {
482        return -3
483    }
484    # delete the records starting with "RGBD x" or "RGBD10"
485    foreach key [array names ::exparray "${key}*"] {
486        #puts $key
487        delexp $key
488    }
489    # decrement the RB counter
490    set n [string trim [readexp "RGBD  NRBDS"]]
491    if {$n == ""} {
492        set n 0
493    }
494    incr n -1
495    validint n 5
496    if {$n > 0} {
497        setexp "RGBD  NRBDS" $n 1 5
498    } else {
499        delexp "RGBD  NRBDS"
500    }
501    return 1
502}
503
504# ReplaceRigidBody: replace all the information for rigid body #rbnum
505# Works the sames as AddRigidBody (see above) except that the rigid body is replaced rather
506# than added.
507# Note that count of the # of times the body is used is preserved
508proc ReplaceRigidBody {rbnum multlist coordlist} {
509    set value $rbnum
510    validint value 2
511    set key "RGBD${value}"
512    set line [readexp "$key NBDS"]
513    foreach key [array names ::exparray "${key}*"] {
514        #puts $key
515        delexp $key
516    }
517    SetRigidBody $rbnum $multlist $coordlist
518    setexp "$key NBDS" $line 1 68
519}
520
521# Edit the parameters for rigid body #rbnum
522# (normally called from ReplaceRigidBody or AddRigidBody)
523proc SetRigidBody {rbnum multlist coordlist} {
524    set value $rbnum
525    validint value 2
526    set key "RGBD${value}"
527    # number of atoms
528    set value [llength [lindex $coordlist 0]]
529    validint value 5
530    if {! [existsexp "$key NATR"]} {
531        makeexprec "$key NATR"
532    }
533    setexp "$key NATR" $value 1 5
534    # number of times used
535    set value 0
536    validint value 5
537    if {! [existsexp "$key NBDS"]} {
538        makeexprec "$key NBDS"
539    }
540    setexp "$key NBDS" $value 1 5
541    # number of coordinate matrices
542    set value [llength $multlist]
543    validint value 5
544    if {! [existsexp "$key NSMP"]} {
545        makeexprec "$key NSMP"
546    }
547    setexp "$key NSMP" $value 1 5
548    set i 0
549    foreach mult $multlist coords $coordlist {
550        incr i
551        if {! [existsexp  "${key}${i}PARM"]} {
552            makeexprec "${key}${i}PARM"
553        }
554        setexp "${key}${i}PARM" [format "%10.5f%5d%5d" $mult 0 0] 1 20
555        set j 0
556        foreach item $coords {
557            #puts $item
558            incr j
559            set value $j
560            validint value 3
561            if {! [existsexp "${key}${i}SC$value"]} {
562                makeexprec "${key}${i}SC$value"
563            }
564            if {[llength $item] == 4} {
565                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
566            } elseif {[llength $item] == 3} {
567                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f" $item] 1 30
568            } else {
569                return -code 3 "Invalid number of coordinates"
570            }
571        }
572    }
573}
574
575# convert a decimal to the GSAS hex encoding with a field $digits long.
576proc ToHex {num digits} {
577    return [string toupper [format "%${digits}x" $num]]
578}
579
580# convert a GSAS hex encoding to a decimal integer
581proc FromHex {hex} {
582    return [scan $hex "%x"]
583}
584
585# MapRigidBody: define an "instance" of a rigid body: meaning that the coordinates
586# (and optionally U values) for a set of atoms will be generated from the rigid body
587# arguments:
588#   phase: phase number (1-9)
589#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
590#   firstatom: sequence number of the first atom in phase (note that atoms may
591#              not be numbered sequentially)
592#   position: list of three fractional coordinates for the origin of the rigid body coordinates
593#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
594#           cartesian system before the body is translated to position.
595# returns the instance # (number of times body $bodytyp has been used in phase $phase)
596proc MapRigidBody {phase bodytyp firstatom position angles} {
597    # find the first unused body # for this phase & type
598    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
599        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
600        if {! [existsexp "$key  NDA"]} {break}
601    }
602    # did we go too far?
603    if {$rbnum == 16} {return ""} 
604    # increment number of mapped bodies of this type overall
605    set value $bodytyp
606    validint value 2
607    set key "RGBD${value}"
608    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
609    incr used
610    set value $used
611    validint value 5
612    setexp "$key NBDS" $value 1 5
613    # increment number of mapped bodies of this type in this phase
614    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
615    if {[existsexp "$key  NBDS"]} {
616        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
617    } else {
618        makeexprec "$key  NBDS"
619        set used 0
620    }
621    incr used
622    set value $used
623    validint value 5
624    setexp "$key  NBDS" $value 1 5
625    # now write the mapping parameters
626    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
627    set value $firstatom
628    validint value 5
629    makeexprec "$key  NDA"
630    setexp "$key  NDA" $value 1 5
631    set l1 {}
632    set l2 {}
633    for {set i 0} {$i < 9} {incr i} {
634        append l1 [format %5d 0]
635        append l2 [format %1d 0]
636    }
637    makeexprec "$key BDFL"
638    setexp "$key BDFL" $l1$l2 1 54
639    makeexprec "${key} BDLC"
640    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
641    makeexprec "${key} BDOR"
642    set l1 {}
643    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
644        append l1 [format "%8.2f%2d" $val $dir]
645    }
646    setexp "${key} BDOR" $l1 1 60
647    makeexprec "${key} LSTF"
648    setexp "${key} LSTF" [format "%5d" 0] 1 5
649    return $rbnum
650}
651
652# EditRigidBodyMapping: edit parameters that define an "instance" of a rigid body (see MapRigidBody)
653# arguments:
654#   phase: phase number (1-9)
655#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
656#   bodynum: instance number, as returned by MapRigidBody
657#   position: list of three fractional coordinates for the origin of the rigid body coordinates
658#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
659#           cartesian system before the body is translated to position.
660#
661proc EditRigidBodyMapping {phase bodytyp bodynum position angles} {
662    # number of bodies of this type in this phase
663    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $bodynum 1]"
664    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
665    set l1 {}
666    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
667        append l1 [format "%8.2f%2d" $val $dir]
668    }
669    setexp "${key} BDOR" $l1 1 60
670}
671
672# UnMapRigidBody: remove a rigid body constraint by removing a RB "instance"
673# (undoes MapRigidBody)
674# arguments:
675#   phase: phase number (1-9)
676#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
677#   bodynum: instance number, as returned by MapRigidBody
678proc UnMapRigidBody {phase bodytyp mapnum} {
679    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $mapnum] == -1} {
680        return ""
681    }
682    # decrement number of mapped bodies of this type overall
683    set value $bodytyp
684    validint value 2
685    set key "RGBD${value}"
686    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
687    incr used -1
688    set value $used
689    validint value 5
690    setexp "$key NBDS" $value 1 5
691    # decrement number of mapped bodies of this type in this phase
692    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
693    if {[existsexp "$key  NBDS"]} {
694        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
695    } else {
696        set used 0
697    }
698    incr used -1
699    set value $used
700    validint value 5
701    if {$used > 0} {
702        setexp "$key  NBDS" $value 1 5
703    } else {
704        delexp "$key  NBDS"
705    }
706    # now delete the mapping parameter records
707    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $mapnum 1]"
708    foreach key [array names ::exparray "${key}*"] {
709        delexp $key
710    }
711    return $used
712}
713
714#============================================================================
715# Rigid body utility routines
716#============================================================================
717# RigidBodyGetVarNums: Returns a list of the variable numbers in use
718#       for rigid body variable parameters.
719# RigidBodyAtomNums: returns a list of atom numbers that are mapped to
720#       rigid bodies in a selected phase
721# RigidStartAtoms: returns a list of atoms that are allowed for creation of RB
722# ExtractRigidBody: Use the GSAS geometry program to cartesian coordinates &
723#       setting info for a RB from fractional coordinates for atoms in a phase
724# RunRecalcRBCoords: updates the coordinates in all phases after changes have
725#       been made to rigid parameters.
726# CalcBody: Convert ortho to fractional coordinates using RB parameters
727# FitBody: Optimize the origin and Euler angles to match a rigid body to a
728#       set of fractional coordinates
729# zmat2coord: convert a z-matrix to a set of cartesian coordinates
730# RB2cart: convert the representation used for rigid bodies into
731#       cartesian coordinates
732# PlotRBtype: plot a rigid body with DRAWxtl
733# PlotRBcoords: plot orthogonal coordinates with DRAWxtl
734# DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed
735#       on a structure read from a phase with DRAWxtl
736#============================================================================
737#============================================================================
738# RigidBodyGetVarNums: Returns a list of the variable numbers used already
739# for rigid body variable parameters
740proc RigidBodyGetVarNums {} {
741    set varlist {}
742    foreach type [RigidBodyList] {
743        set typelist [lindex [ReadRigidBody $type] 1]
744        foreach item $typelist {
745            lappend varlist [lindex $item 2]
746        }
747        foreach phase $::expmap(phaselist) {
748            foreach i [RigidBodyMappingList $phase $type] {
749                set items [ReadRigidBodyMapping $phase $type $i]
750                set varlist [concat $varlist [lindex $items 3]]
751                if {[llength [lindex $items 6]] > 0} {
752                    set varlist [concat $varlist [lindex $items 6]]
753                }
754            }
755        }
756    }
757    return [lsort -integer -unique $varlist]
758}
759
760# RigidBodyAtomNums: Returns a list of the atoms mapped to rigid bodies in
761# phase $phase
762proc RigidBodyAtomNums {phase} {
763    if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""}
764    set allatoms $::expmap(atomlist_$phase)
765    # get matching atoms coordinate range
766    set mappedlist {}
767    foreach type [RigidBodyList] {
768        foreach i [RigidBodyMappingList $phase $type] {
769            # get the number of atoms in this type of body
770            set natoms [llength [lindex [lindex [lindex [ReadRigidBody $type] 1] 0] 3]]
771            set natom1 [expr {$natoms - 1}]
772            set items [ReadRigidBodyMapping $phase $type $i]
773            set firstatom [lindex $items 0]
774            set firstind [lsearch $allatoms $firstatom]
775            set mappedlist [concat $mappedlist \
776                                [lrange \
777                                     [lrange $allatoms $firstind end] \
778                                     0 $natom1] \
779                               ]
780        }
781    }
782    return [lsort -integer $mappedlist]
783}
784
785# RigidStartAtoms: Find allowed starting atoms for a rigid body in a phase
786# Input:
787#   phase is the phase number
788#   natoms is the number of atoms in the RB to be mapped
789# Returns a list of valid "start" atoms.
790# Example: if the atom numbers in the phase are {2 4 5 6 7 8} and no rigid bodies
791# are mapped, then a 4-atom body can be mapped starting with atom 2, 4 or 5 only,
792# so {2 4 5} is returned
793# If atoms 2-6 were already mapped, then this routine would return an empty
794# list, as it is not possible to map the body.
795proc RigidStartAtoms {phase natoms} {
796    if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""}
797    set allatoms $::expmap(atomlist_$phase)
798    set usedatoms [RigidBodyAtomNums $phase]
799    set startatomlist {}
800    for {set i 0} {$i < [llength $allatoms]} {incr i} {
801        set al [lrange $allatoms $i [expr {$i+$natoms-1}]]
802        if {[llength $al] < $natoms} break
803        set ok 1
804        foreach atom $al {
805            if {[lsearch $usedatoms $atom] != -1} {
806                set ok 0
807                break
808            }
809        }
810        if $ok {lappend startatomlist [lindex $al 0]}
811    }
812    return $startatomlist
813}
814
815# ExtractRigidBody: Use the GSAS geometry program to compute a set of cartesian coordinates for a
816# set of atoms in a .EXP file and provide the origin shift and Euler angles needed to
817# place the cartesian system into the crystal coordinates. Used for setting up a rigid body.
818# Returns a nested list of lists:
819#   element 0: a list of the origin location {x y z} in fraction coordinates
820#   element 1: a list of three rotation angles in form {a1 a2 a3}
821#              where a1, a2 and a3 are rotations around the cartesian x, y and z axes
822#   element 2: a list of $natom cartesian coordinate triples {{x1 y1 z1} {x2 y2 z2}...}
823# arguments:
824    # phase: phase #
825    # natom: number of atoms in group
826    # firstatom: sequence # in phase (may be > than number of the atom)
827    # originlist: atoms to define origin (where 1 is first atom in group; <= natom)
828    # vector1: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)  (for example {X 1 3})
829    # vector2: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)
830    # note that vector2 must define a different axis than vector1
831    # also and vector1 and vector2 cannot use the same atom pair
832proc ExtractRigidBody {phase natom firstatom originlist vector1 vector2} {
833    global expgui
834    set fp [open "geom.inp" "w"]
835    puts $fp "N"
836    if {[llength ::expmap(phaselist)] > 1} {
837       # select phase
838       puts $fp "N"
839       puts $fp $phase
840       puts $fp "N"
841    }
842    puts $fp "R"
843    puts $fp "$natom"
844    puts $fp "$firstatom"
845    puts $fp [llength $originlist]
846    foreach i $originlist {
847        puts $fp $i
848    }
849    foreach i [concat $vector1 $vector2] {
850        puts $fp $i
851    }
852    puts $fp "0"
853    puts $fp "X"
854    close $fp
855    #puts "[file join $expgui(gsasexe) geometry] $expgui(expfile) < geom.inp > geom.out"
856    catch {
857        exec [file join $expgui(gsasexe) geometry] $expgui(expfile) < geom.inp > geom.out
858    } err
859    #puts $err
860    file delete geom.inp
861    set fp [open geom.out r]
862    set origin {}
863    set Euler {}
864    set coordlist {}
865    while {[gets $fp line] >= 0} {
866        if {[string first "Cell coordinates of origin" $line] != -1} {
867            set origin [lrange [string range $line [string first "are" $line] end] 1 3]
868            #puts "origin in rb = $origin"
869        }
870        if {[string first "Rotation angles" $line] != -1} {
871            set Euler {}
872            foreach i [lrange [split $line "="] 1 3] {
873                lappend Euler [lindex $i 0]
874            }
875            #puts $line
876            #puts $Euler
877        }
878        if {[string first "Atom   Orthon" $line] != -1} {
879            set coordlist {}
880            for {set i 1} {$i <= $natom} {incr i} {
881                gets $fp line
882                set coord {}
883                lappend coord [string trim [string range $line 9 15]]
884                lappend coord [string trim [string range $line 16 22]]
885                lappend coord [string trim [string range $line 23 29]]
886                lappend coord [string trim [string range $line 0 8]]
887                #puts $line
888                #puts $coord
889                lappend coordlist $coord
890            }
891            #puts $coordlist
892        }
893    }
894    #file delete geom.out
895    if {[llength $origin] == 0 || [llength $Euler] == 0 || [llength $coordlist] == 0} {
896       puts "Error: run of GEOMETRY failed"
897    }
898    return [list $origin $Euler $coordlist]
899}
900
901# RunRecalcRBCoords: updates the coordinates in a .EXP file after a rigid
902# body has been changed, mapped or the setting info is changed
903proc RunRecalcRBCoords { } {
904    global expgui tcl_platform
905    set input [open resetmult.inp w]
906    puts $input "Y"
907    puts $input "l b"
908    puts $input "n"
909    puts $input "x x x"
910    puts $input "x"
911    close $input
912    # Save the current exp file
913    savearchiveexp
914    # disable the file changed monitor
915    set expgui(expModifiedLast) 0
916    set expnam [file root [file tail $expgui(expfile)]]
917    set err [catch {
918        if {$tcl_platform(platform) == "windows"} {
919            exec [file join $expgui(gsasexe) expedt.exe] $expnam < resetmult.inp >& resetmult.out
920        } else {
921            exec [file join $expgui(gsasexe) expedt] $expnam < resetmult.inp >& resetmult.out
922        }
923    } errmsg]
924    loadexp $expgui(expfile)
925    set fp [open resetmult.out r]
926    set out [read $fp]
927    close $fp
928    set expgui(exptoolout) $out
929    catch {file delete resetmult.inp resetmult.out}
930    if {$err} {
931        return $errmsg
932    } else {
933        return ""
934    }
935}
936
937
938# compute a rotation matrix for orthogonal coordinates (based on MAKMATD in GSAS)
939# rotate angle degrees around axis (1, 2 or 3) for (x, y, or z)
940# returns a list that can be used as a matrix in the La package
941proc RotMatrix {axis angle} {
942    set ang [expr {$angle * acos(0) / 90.}]
943    set mat "1 0 0 0 1 0 0 0 1"
944    if {$axis == 1}  {
945        set i1 1
946        set i2 2
947    } elseif {$axis == 2}  {
948        set i1 2
949        set i2 0
950    } else {
951        set i1 0
952        set i2 1
953    }
954    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
955    foreach item {
956        {$i1 $i1 [expr {cos($ang)}]}
957        {$i2 $i2 [expr {cos($ang)}]}
958        {$i1 $i2 [expr {-sin($ang)}]}
959        {$i2 $i1 [expr {sin($ang)}]}
960    } {
961        foreach {c r val} [subst $item] {}
962        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
963    }
964    return "2 3 3 $mat"
965}
966
967# compute the derivative of the rotation matrix with respect to the angle, see RotMatrix
968# (based on MAKMATD in GSAS)
969# returns a list that can be used as a matrix in the La package
970proc DerivRotMatrix {axis angle} {
971    set ang [expr {$angle * acos(0) / 90.}]
972    set mat "0 0 0 0 0 0 0 0 0"
973    if {$axis == 1}  {
974        set i1 1
975        set i2 2
976    } elseif {$axis == 2}  {
977        set i1 2
978        set i2 0
979    } else {
980        set i1 0
981        set i2 1
982    }
983    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
984    foreach item {
985        {$i1 $i1 [expr {-sin($ang) * acos(0) / 90.}]}
986        {$i2 $i2 [expr {-sin($ang) * acos(0) / 90.}]}
987        {$i1 $i2 [expr {-cos($ang) * acos(0) / 90.}]}
988        {$i2 $i1 [expr {cos($ang) * acos(0) / 90.}]}
989    } {
990        foreach {c r val} [subst $item] {}
991        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
992    }
993    return "2 3 3 $mat"
994}
995
996# compute an orthogonalization matrix from cell parameters (based on AMATRX in GSAS)
997# returns a list that can be used as a matrix in the La package
998proc OrthoMatrix {a b c alpha beta gamma} {
999    set CA [expr {cos($alpha * acos(0) / 90.)}]
1000    set CB [expr {cos($beta * acos(0) / 90.)}]
1001    set CG [expr {cos($gamma * acos(0) / 90.)}]
1002    set SA [expr {sin($alpha * acos(0) / 90.)}]
1003    set SB [expr {sin($beta * acos(0) / 90.)}]
1004    set SC [expr {sin($gamma * acos(0) / 90.)}]
1005    set CASTAR [expr { ($CB*$CG-$CA)/($SB*$SC) }]    ;#! cos(Alpha*)
1006    set CBSTAR [expr { ($CA*$CG-$CB)/($SA*$SC) }]    ;#! cos(Beta*)
1007    set CCSTAR [expr { ($CA*$CB-$CG)/($SA*$SB) }]    ;#! cos(Gamma*)
1008    set SASTAR [expr { sqrt(1.0-($CASTAR*$CASTAR*2)) }]    ;#! sin(Alpha*)
1009    set SBSTAR [expr { sqrt(1.0-($CBSTAR*$CBSTAR*2)) }]    ;#! sin(Beta*)
1010    set SCSTAR [expr { sqrt(1.0-($CCSTAR*$CCSTAR*2)) }]    ;#! sin(Gamma*)
1011
1012    set A  "2 3 3      $a 0 0    0 $b 0    0 0 $c"
1013    set A1 "2 3 3     1.0 0 0    $CG [expr {$SASTAR*$SC}] [expr {-$CASTAR*$SC}]    $CB 0.0 $SB"
1014    #!This matrix is
1015    #!   (1.0      0.0            0.0             )
1016    #!   (cos(G)  sin(A*)*sin(G) -cos(A*)*sin(G)  )
1017    #!   (cos(B)   0.0            sin(B)          )
1018    return [La::transpose [La::mmult $A $A1]]
1019}
1020
1021# compute the transformation matrix that converts a rigid body coordinates into fractional
1022# coordinates
1023# arguments:
1024#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
1025#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
1026#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
1027# returns a list that can be used as a matrix in the La package
1028proc CalcXformMatrix {rotations cellprms} {
1029    set prod {}
1030    foreach item $rotations {
1031        #puts $item
1032        set mat [eval RotMatrix $item]
1033        if {$prod == ""} {
1034            set prod $mat
1035        } else {
1036            set prod [La::mmult $prod $mat]
1037        }
1038    }
1039    #puts "--- rotation product ---"
1040    #puts [La::show $prod]
1041
1042    set ortho [eval OrthoMatrix $cellprms]
1043    #puts "--- ortho ---"
1044    #puts [La::show $ortho]
1045    set deortho [La::msolve $ortho [La::mident 3] ]
1046    #puts "--- deortho ---"
1047    #puts [La::show $deortho]
1048    #puts "--- xform ---"
1049    set xform [La::mmult $deortho $prod]
1050    return $xform
1051}
1052
1053# transforms a triplet of orthogonal coordinates into fractional ones using
1054# arguments:
1055#    xform: a transformation matrix from CalcXformMatrix
1056#    origin: a list of fraction coordinates {x y z} describing the location of the
1057#            origin of the orthogonal coordinates in the crystal system
1058#    ortho: a triplet of othogonal coordinates
1059# returns a triplet of fractional coordinates
1060proc Ortho2Xtal {xform origin ortho} {
1061    set ocv "2 3 0 $ortho"
1062    set frac [La::mmult $xform $ocv]
1063    #puts [La::show [La::transpose $frac]]
1064    #puts $frac
1065    set frac [La::madd $frac "[lrange $frac 0 2] $origin"]
1066    #puts [La::show [La::transpose $frac]]
1067    return $frac
1068}
1069
1070# compute the derivative of the transformation matrix (see CalcXformMatrix)
1071# with respect to every rotation in the $rotations list
1072# arguments:
1073#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
1074#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
1075#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
1076# returns a list of matrices where each matrix is a list that can be used as a
1077# matrix in the La package
1078proc CalcDerivMatrix {rotations cellprms} {
1079    set ortho [eval OrthoMatrix $cellprms]
1080    set deortho [La::msolve $ortho [La::mident 3] ]
1081    set derivlist {}
1082
1083    foreach item $rotations {
1084        #puts $item
1085        set mat [eval DerivRotMatrix $item]
1086        #puts $item
1087        #puts [La::show $mat]
1088        set xform [La::mmult $deortho $mat]
1089        lappend derivlist $xform
1090    }
1091    return $derivlist
1092}
1093
1094# CalcBody: Calculate fractional coordinates using rigid body setting parameters
1095# arguments:
1096#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
1097#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
1098#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
1099#  ortholist: list containing triplets with orthogonal coordinates
1100#  origin: a list of fraction coordinates {x y z} describing the location of the
1101#            origin of the orthogonal coordinates in the crystal system
1102#     note that the length of ortholist, useflag and fraclist should be the same
1103# Returns a list with the computed fractional coordinates for all atoms
1104proc CalcBody {Euler cell ortholist origin} {
1105    set xform [CalcXformMatrix $Euler $cell]
1106    set i 0
1107    set sumdvs 0
1108    set fracout {}
1109    set rmsout {}
1110    foreach oc $ortholist {
1111        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
1112        lappend fracout $frac
1113    }
1114    return $fracout
1115}
1116
1117
1118# fit a rigid body's origin
1119# arguments:
1120#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
1121#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
1122#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
1123#  ortholist: list containing triplets with orthogonal coordinates
1124#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
1125#  fraclist: list containing triplets with fractional coordinates 
1126#  origin: a list of fraction coordinates {x y z} describing the location of the
1127#            origin of the orthogonal coordinates in the crystal system
1128#     note that the length of ortholist, useflag and fraclist should be the same
1129# Returns a list with the following elements
1130#   0: a list with three new origin values
1131#   1: the root-mean square difference between the fraclist coordinates and those computed
1132#      using the input values for those atoms where $use is one (in Angstroms)
1133#   2: the computed fractional coordinates for all atoms
1134#   3: individual rms values for all atoms (in Angstroms)
1135# note that items 1-3 are computed with the imput origin, not the revised one
1136proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} {
1137puts $fraclist
1138    set xform [CalcXformMatrix $Euler $cell]
1139    #puts "entering FitBodyOrigin"
1140    foreach var {x y z} {set sum($var) 0.0}
1141    set i 0
1142    set sumdvs 0
1143    set fracout {}
1144    set rmsout {}
1145    foreach oc $ortholist use $useflag coord $fraclist {
1146        #puts "ortho: $oc"
1147        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
1148        lappend fracout $frac
1149        if {$use} {incr i}
1150        set dvs 0
1151        foreach var {x y z} v1 $frac v2 $coord abc [lrange $cell 0 2] {
1152            #puts "v2 = $v2"
1153            #puts "v1 = $v1"
1154            #puts "abc = $abc"
1155            set dv [expr {($v2 - $v1)*$abc}]
1156            set dvs [expr {$dvs + $dv*$dv}]
1157            set sumdvs [expr {$sumdvs + $dv*$dv}]
1158            if {$use} {set sum($var) [expr {$sum($var) + $dv/$abc}]}
1159            #puts "round and round"
1160        }
1161        lappend rmsout [expr {sqrt($dvs)}]
1162    }
1163    set rms 0
1164    if {$i > 1} {set rms [expr {sqrt($sumdvs)/$i}]}
1165    set neworig {}
1166    foreach var {x y z} o $origin {
1167        lappend neworig [expr {$o + ($sum($var)/$i)}]
1168    }
1169    return [list $neworig $rms $fracout $rmsout]
1170}
1171
1172# fit a rigid body's Euler angles using least-squares
1173# arguments:
1174#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
1175#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
1176#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
1177#  ortholist: list containing triplets with orthogonal coordinates
1178#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
1179#  fraclist: list containing triplets with fractional coordinates 
1180#  origin: a list of fraction coordinates {x y z} describing the location of the
1181#            origin of the orthogonal coordinates in the crystal system
1182#     note that the length of ortholist, useflag and fraclist should be the same
1183# Returns a list of new Euler angles
1184proc FitBodyRot {Euler cell ortholist useflag fraclist origin} {
1185    set xform [CalcXformMatrix $Euler  $cell]
1186    set derivlist [CalcDerivMatrix $Euler  $cell]
1187    set A "2 [expr 3*[llength $ortholist]] 3"
1188    foreach oc $ortholist use $useflag coord $fraclist {
1189        if {! $use} continue
1190        foreach deriv $derivlist {
1191            foreach xyz [lrange [Ortho2Xtal $deriv "0 0 0" $oc] 3 end] {
1192                lappend A $xyz
1193            }
1194        }
1195    }
1196    #puts "A: [La::show $A]"
1197    set y "2 [expr 3*[llength $ortholist]] 1"
1198    foreach oc $ortholist use $useflag coord $fraclist {
1199        if {! $use} continue
1200        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
1201        foreach xyz $coord XYZ $frac {
1202            lappend y [expr {$XYZ - $xyz}]
1203        }
1204    }
1205
1206    set AtA [La::mmult [La::transpose $A] $A]
1207    set Aty [La::mmult [La::transpose $A] $y]
1208   
1209    set l {}
1210    #set shifts {}
1211    foreach delta [lrange [La::msolve $AtA $Aty] 3 end] item $Euler {
1212        #lappend shifts $delta
1213        lappend l "[lindex $item 0] [expr {$delta + [lindex $item 1]}]"
1214    }
1215    #puts "shifts = $shifts"
1216    return $l
1217}
1218
1219# FitBody: fit a rigid body's Origin and Euler angles
1220# arguments:
1221#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
1222#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
1223#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
1224#  ortholist: list containing triplets with orthogonal coordinates
1225#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
1226#  fraclist: list containing triplets with fractional coordinates 
1227#  origin: a list of fraction coordinates {x y z} describing the location of the
1228#            origin of the orthogonal coordinates in the crystal system
1229#     note that the length of ortholist, useflag and fraclist should be the same
1230# Returns a list containing
1231#   new origin
1232#   new Euler angles
1233#   total rms
1234#   fractional coordinates
1235#   rms deviation in fractional coordinates of new Euler angles
1236proc FitBody {Euler cell ortholist useflag fraclist origin "ncycle 5"} {
1237    #puts "start origin = $origin"
1238    foreach {
1239        origin 
1240        startrms
1241        fracout
1242        rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
1243    #puts "start rms = $startrms"
1244    set rmsprev $startrms
1245    #puts "new origin = $origin"
1246    for {set i 0} {$i < $ncycle} {incr i} {
1247        set Eulerprev $Euler
1248        set Euler [FitBodyRot $Euler $cell $ortholist $useflag $fraclist $origin]
1249        #puts "New Euler $Euler"
1250        #puts "after fit"
1251        foreach {
1252            origin 
1253            rms
1254            fracout
1255            rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
1256        if {$rms > (1.1 * $rmsprev) + 0.01} {
1257            #puts "rms = $rms, new origin = $origin"
1258            set rmsprev $rms
1259        }
1260    } 
1261    #proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} 
1262    #return "$neworig $rms $fracout $rmsout"
1263    set fmt  {"%8.5f %8.5f %8.5f     %8.5f %8.5f %8.5f   %6.3f"}
1264    #foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
1265        #puts "[eval format $fmt $fracin $fraccalc $rmsi]"
1266    #}
1267    return [list $origin $Euler $rms $fracout $rmsout]
1268}
1269
1270# zmat2coord: convert a z-matrix to a set of cartesian coordinates
1271#   a z-matrix is also known as "internal coordinates" or "torsion space"
1272#   (see Journal of Computational Chemistry, Vol 26, #10, p. 1063–1068, 2005 or
1273#    http://www.cmbi.ru.nl/molden/zmat/zmat.html)
1274# INPUT:
1275#   atmlist is a list of ascii lines where each line contains
1276#     lbl c1 distance c2 angle c3 torsion
1277#   where each atom is computed from the previous where the new atom is:
1278#     distance $distance from atom $c1 (angstrom)
1279#     angle $angle from $c1--$c2 (degrees)
1280#     torsion $torsion from $c1--$c2--$c3 (degrees)
1281# OUTPUT:
1282#  zmat2coord returns a list of atom labels and cartesian coordinates,
1283#  with 4 items in each element (label, x, y, z)
1284# this routine was tested against results from Babel via the web interface at
1285# http://www.shodor.org/chemviz/zmatrices/babel.html and sample input at
1286# http://iopenshell.usc.edu/howto/zmatrix/
1287proc zmat2coord {atmlist} { 
1288    set torad [expr {acos(0)/90.}]
1289    set i 0
1290    foreach line $atmlist {
1291        incr i
1292        foreach {lbl c1 dist c2 angle c3 torsion} $line {} 
1293        if {$i == 1} {
1294            set atm(1) [list $lbl 0 0 0] ; # 1st atom is at origin
1295        } elseif {$i == 2} {
1296            set dist1 $dist
1297            set atm(2) [list $lbl $dist1 0 0] ; # 2nd atom is along x-axis
1298        } elseif {$i == 3} {
1299            # 3rd atom can be bonded to the 1st or 2nd
1300            if {$c1 == 1} {
1301                set atm(3) [list $lbl \
1302                                [expr {$dist * cos($torad * $angle)}] \
1303                                [expr {$dist * sin($torad * $angle)}] \
1304                                0]
1305            } else {
1306                set atm(3) [list $lbl \
1307                                [expr {$dist1 - $dist * cos($torad * $angle)}] \
1308                                [expr {$dist * sin($torad * $angle)}] \
1309                                0]
1310            }
1311        } else {
1312            set atm($i) [concat $lbl \
1313                             [ahcat "atm" $c1 $dist $c2 $angle $c3 $torsion]]
1314        }
1315    }
1316    set coordlist {}
1317    foreach key [lsort -integer [array names atm]] {
1318        lappend coordlist $atm($key)
1319    }
1320    return $coordlist
1321}
1322# Compute the length of a vector
1323proc vlen {a} {
1324    set sum 0.0
1325    foreach ai $a {
1326        set sum [expr {$sum + $ai*$ai}]
1327    }
1328    return [expr sqrt($sum)]
1329}
1330# compute vector (a + z * b) and optionally normalize to length d
1331proc vadd {a b d z} {
1332    set c {}
1333    foreach ai $a bi $b {
1334        lappend c [expr {$bi + $z * $ai}]
1335    }
1336    set v [vlen $c]
1337    if {$d != 0} {
1338        set r {}
1339        foreach ci $c {
1340            lappend r [expr {$d * $ci / $v}]
1341        }
1342        return [list $v $r]
1343    }
1344    return [list $v $c]
1345}
1346# normalize a vector
1347proc vnrm {x} {
1348    set v [vlen $x]
1349    if {abs($v) < 1e-8} {return [list 0 0 0]}
1350    set y {}
1351    foreach xi $x {
1352        lappend y [expr {$xi / $v}]
1353    }
1354    return $y
1355}
1356# compute the coordinates for an atom that is bonded:
1357#   distance $dist from atom $nc
1358#   angle $bondang from $nc--$nb
1359#   torsion $torsang from $nc--$nb--$na
1360#   coordinates are found in array $atmarr in the calling routine
1361# based on a Fortran routine provided by Peter Zavalij (Thanks Peter!)
1362proc ahcat {atmarr nc dist nb bondang na torsang} {
1363    upvar 1 $atmarr atm
1364    set xa [lrange $atm($na) 1 3]
1365    set xb [lrange $atm($nb) 1 3]
1366    set xc [lrange $atm($nc) 1 3]
1367    set torad [expr {acos(0)/90.}]
1368    # x = unit Vector A-B
1369    foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. -1.] 1] {}
1370    # y = unit Vector C-B
1371    set y [lindex [vadd $xb $xc 1. -1.] 1]
1372    foreach {y1 y2 y3} $y {}
1373    set z1 [expr {$x2*$y3 - $x3*$y2}]
1374    set z2 [expr {$x3*$y1 - $x1*$y3}]
1375    set z3 [expr {$x1*$y2 - $x2*$y1}]
1376    set z [vnrm [list $z1 $z2 $z3]]
1377    set q1 [expr {$y2*$z3 - $y3*$z2}]
1378    set q2 [expr {$y3*$z1 - $y1*$z3}]
1379    set q3 [expr {$y1*$z2 - $y2*$z1}]
1380    set q [vnrm [list $q1 $q2 $q3]]
1381    set th [expr {$bondang * $torad}]
1382    set ph [expr {-1. * $torsang * $torad}]
1383    set cth [expr {cos($th)}]
1384    set sth [expr {sin($th)}]
1385    set cph [expr {cos($ph)}]
1386    set sph [expr {sin($ph)}]
1387    set xh {}
1388    foreach xci $xc xi $q zi $z yi $y {
1389        lappend xh [expr {
1390                          $xci +
1391                          $dist*($sth*($cph*$xi + $sph*$zi)-$cth*$yi)
1392                      }]
1393    }
1394    return $xh
1395}
1396 
1397# RB2cart: convert the rigid body representation reported as the 2nd element
1398# in ReadRigidBody into cartesian coordinates
1399#   rblist: a list containing an element for each scaling factor
1400# in each element there are four items:
1401#    a multiplier value for the rigid body coordinates
1402#    a damping value (0-9) for the refinement of the multiplier (not used)
1403#    a variable number if the multiplier will be refined (not used)
1404#    a list of cartesian coordinates coordinates
1405# each cartesian coordinate contains 4 items: x,y,z and a label
1406# returns a list of coordinate triplets
1407proc RB2cart {rblist} {
1408    foreach item $rblist {
1409        foreach {mult damp ref coords} $item {}
1410        set i 0
1411        foreach xyz $coords {
1412            foreach {x y z} [lrange $xyz 0 2] {}
1413            foreach val [lrange $xyz 0 2] var {X Y Z} {
1414                if {[array names $var $i] == ""} {
1415                    set ${var}($i) [expr {$mult * $val}]
1416                } else {
1417                    set ${var}($i) [expr {[set ${var}($i)] + $mult * $val}]
1418                }
1419            }
1420            incr i
1421        }
1422    }
1423    set out ""
1424    foreach i [lsort -integer [array names X]] {
1425        lappend out [list $X($i) $Y($i) $Z($i)]
1426    }
1427    return $out
1428}
1429
1430# get the name of the DRAWxtl application, if installed
1431proc GetDRAWxtlApp {} {
1432    # is DRAWxtl installed?
1433    set app {}
1434    if {![catch {set fp [open [file join $::env(HOME) .drawxtlrc] r]}]} {
1435        # line 12 is name of executable
1436        set i 0
1437        while {$i < 12} {
1438            incr i
1439            gets $fp appname
1440        }
1441        close $fp
1442        set app [auto_execok $appname]
1443    }
1444    if {$app != ""} {
1445        return $appname
1446    }
1447    return ""
1448}
1449
1450# DRAWxtlPlotOrtho: plot orthogonal coordinates in DRAWxtl
1451# input:
1452#  filename: file name for the .str file to create
1453#  title: string for title in .str file
1454#  coords: cartesian coordinates
1455#  bondlist: list of bonds to draw as min, max length (A) and
1456#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
1457proc DRAWxtlPlotOrtho {filename title coords bondlist} {
1458    foreach {xmin ymin zmin} {"" "" ""} {}
1459    foreach {xmax ymax zmax} {"" "" ""} {}
1460    foreach xyz $coords {
1461        foreach {x y z} $xyz {}
1462        foreach s {x y z} {
1463            foreach t {min max} {
1464                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
1465            } 
1466            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
1467            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
1468        }
1469    }
1470    #puts "$xmin $xmax $ymin $ymax $zmin $zmax"
1471    set max $xmin
1472    foreach val "$xmin $xmax $ymin $ymax $zmin $zmax" {
1473        if {$max < abs($val)} {set max $val}
1474    }
1475    set scale [expr {4.*$max}]
1476    set a 10.
1477    lappend range [expr -0.01+($xmin/$scale)] [expr 0.01+($xmax/$scale)] \
1478        [expr -0.01+($ymin/$scale)] [expr 0.01+($ymax/$scale)] \
1479        [expr -0.01+($zmin/$scale)] [expr 0.01+($zmax/$scale)]
1480    set fp [open $filename w]
1481    puts $fp "title $title"
1482    puts $fp "box  0.000 Black"
1483    puts $fp "background White"
1484    #puts $fp "nolabels"
1485    puts $fp "cell $a $a $a 90 90 90"
1486    puts $fp "spgr P 1"
1487    puts $fp "pack $range"
1488    set i 0
1489    foreach xyz $coords {
1490        foreach {x y z} $xyz {}
1491        incr i
1492        puts $fp "atom c $i [expr {$x/$scale}] [expr {$y/$scale}] [expr {$z/$scale}]"
1493        puts $fp "labeltext [expr {0.02 + $x/$scale}] [expr {0.01 + $y/$scale}] [expr {0.01 + $z/$scale}] $i"
1494    }
1495    puts $fp "sphere c  [expr 0.100*($a/$scale)] Red"
1496    puts $fp "finish   0.70   0.30   0.08   0.01"
1497    foreach bondpair $bondlist {
1498        foreach {b1 b2 color} $bondpair {}
1499        if {$color == ""} {set color Red}
1500        puts $fp "bond c c [expr {0.01*$a/$scale}] [expr {$b1*$a/$scale}] [expr {$b2*$a/$scale}] $color"
1501    }
1502    puts $fp "frame"
1503    set range {}
1504    lappend range -0.01 [expr 0.01+(0.1*$a/$scale)] \
1505        -0.01 [expr 0.01+(0.1*$a/$scale)] \
1506        -0.01 [expr 0.01+(0.1*$a/$scale)]
1507    puts $fp "cell $a $a $a 90 90 90"
1508    puts $fp "spgr P 1"
1509    puts $fp "pack $range"
1510    puts $fp "atom o 1 0 0 0"
1511    puts $fp "atom o 2 [expr {0.1*$a/$scale}] 0 0"
1512    puts $fp "atom o 3 0 [expr {0.1*$a/$scale}] 0"
1513    puts $fp "atom o 4 0 0 [expr {0.1*$a/$scale}]"
1514    puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
1515    puts $fp "labelscale 0.5"
1516    puts $fp "labeltext [expr {0.11*$a/$scale}] 0 0 x"
1517    puts $fp "labeltext 0 [expr {0.11*$a/$scale}] 0 y"
1518    puts $fp "labeltext 0 0 [expr {0.11*$a/$scale}] z"
1519    puts $fp "sphere o [expr {0.02*$a/$scale}] Blue"
1520    puts $fp "origin   .0 .0 .0"
1521    puts $fp "end"
1522    close $fp
1523}
1524
1525# PlotRBtype: plot a rigid body in DRAWxtl
1526# input:
1527#  rbtype: # of rigid body
1528#  bondlist: list of bonds to draw as min, max length (A) and
1529#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
1530#  file: file name for the .str file to create
1531proc PlotRBtype {rbtype "bondlist {}" "file {}"} {
1532    set app [GetDRAWxtlApp]
1533    if {$app == ""} {
1534        MyMessageBox -parent . -title "No DRAWxtl" \
1535                -message "Sorry, DRAWxtl is not installed" \
1536                -icon warning
1537        return
1538    }
1539    if {$::tcl_platform(platform) == "windows" && $file == ""} {
1540        set file [file join [pwd] rbplot.str]
1541    } else {
1542        set file "/tmp/rbplot.str"
1543    }
1544    set coords [RB2cart [lindex [ReadRigidBody $rbtype] 1]]
1545    DRAWxtlPlotOrtho $file "" $coords $bondlist
1546    if {$app != ""} {exec $app $file &}
1547}
1548
1549# PlotRBcoords: plot orthogonal coordinates in DRAWxtl
1550# input:
1551#  coords: cartesian coordinates
1552#  bondlist: list of bonds to draw as min, max length (A) and
1553#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
1554#  file: file name for the .str file to create
1555proc PlotRBcoords {coords "bondlist {}" "file {}"} {
1556    set app [GetDRAWxtlApp]
1557    if {$app == ""} {
1558        MyMessageBox -parent . -title "No DRAWxtl" \
1559                -message "Sorry, DRAWxtl is not installed" \
1560                -icon warning
1561        return
1562    }
1563    if {$::tcl_platform(platform) == "windows" && $file == ""} {
1564        set file [file join [pwd] rbplot.str]
1565    } else {
1566        set file "/tmp/rbplot.str"
1567    }
1568    DRAWxtlPlotOrtho $file "" $coords $bondlist
1569    if {$app != ""} {exec $app $file &}
1570}
1571
1572# DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed
1573# on a structure read from a phase
1574# input:
1575#  RBcoords: fractional coordinates for rigid body
1576#  phase:# of phase to plot
1577#  firstatom: seq # of 1st atom in structure to be mapped to rigid body
1578#  allatoms: 0 to plot only atoms in phase that are in the rigid body,
1579#      otherwise plot all atoms
1580#  bondlist: list of bonds to draw for the phase as min, max length (A) and
1581#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
1582#  rbbondlist: list of bonds to draw for the phase as min, max length (A) and
1583#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
1584#  file: optional file name for the .str file to create
1585proc DRAWxtlPlotRBFit {RBcoords phase firstatom "allatoms 0" \
1586                           "bondlist {}" "rbbondlist {}" "file {}"} {
1587    set natom [llength $RBcoords]
1588    set app [GetDRAWxtlApp]
1589    if {$app == ""} {
1590        MyMessageBox -parent . -title "No DRAWxtl" \
1591                -message "Sorry, DRAWxtl is not installed" \
1592                -icon warning
1593        return
1594    }
1595    if {$::tcl_platform(platform) == "windows" && $file == ""} {
1596        set file [file join [pwd] rbplot.str]
1597    } else {
1598        set file "/tmp/rbfit.str"
1599    }
1600
1601    # get rigid body coordinate range
1602    foreach {xmin ymin zmin} {"" "" ""} {}
1603    foreach {xmax ymax zmax} {"" "" ""} {}
1604    foreach xyz $RBcoords {
1605        foreach {x y z} $xyz {}
1606        foreach s {x y z} {
1607            foreach t {min max} {
1608                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
1609            } 
1610            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
1611            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
1612        }
1613    }
1614    set rbrange {}
1615    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
1616                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
1617                     [expr -0.01+$zmin] [expr 0.01+$zmax] ] {
1618        append rbrange [format " %8.4f" $val]
1619    }
1620    set rbcenter [list [expr {($xmin+$xmax)/2}] \
1621                      [expr {($ymin+$ymax)/2}] \
1622                      [expr {($zmin+$zmax)/2}] ]
1623    # get matching atoms coordinate range
1624    set firstind [lsearch $::expmap(atomlist_$phase) $firstatom]
1625    set matchedatomlist [lrange \
1626                      [lrange $::expmap(atomlist_$phase) $firstind end] \
1627                      0 [expr {$natom-1}]]
1628    foreach atom $matchedatomlist {
1629        foreach s {x y z} {
1630            set $s [atominfo $phase $atom $s]
1631            foreach t {min max} {
1632                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
1633            } 
1634            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
1635            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
1636        }
1637    }
1638    # expand to cover at least one unit cell
1639    foreach var {xmin ymin zmin} { 
1640        if {[set $var] > 0.0} {set $var 0.0}
1641    }
1642    foreach var {xmax ymax zmax} { 
1643        if {[set $var] < 1.} {set $var 1.}
1644    }
1645    set range {}
1646    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
1647                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
1648                     [expr -0.01+$zmin] [expr 0.01+$zmax]] {
1649        append range [format " %8.4f" $val]
1650    }
1651
1652    set fp [open $file w]
1653    puts $fp "title structure/rigid-body fit plot"
1654    # plot the structure
1655    puts -nonewline $fp "cell"
1656    foreach p {a b c alpha beta gamma} {
1657        puts -nonewline $fp " [phaseinfo $phase $p]"
1658    }
1659    puts $fp ""
1660    puts $fp "spgp [phaseinfo $phase spacegroup]"
1661    puts $fp "pack $range"
1662    if {$allatoms != 0} {
1663        set atoms $::expmap(atomlist_$phase)
1664    } else {
1665        set firstind [lsearch $::expmap(atomlist_$phase) $firstatom]
1666        set atoms [lrange \
1667                       [lrange $::expmap(atomlist_$phase) $firstind end] \
1668                       0 [expr {$natom-1}]]
1669    }
1670
1671    # set origin at center of rigid body
1672    puts $fp "origin  $rbcenter"
1673    # now loop over atoms
1674    foreach atom $atoms {
1675        set type [atominfo $phase $atom type]
1676        set typelist($type) 1
1677        set xyz ""
1678        foreach v {x y z} {
1679            append xyz "[atominfo $phase $atom $v] "
1680        }
1681        puts $fp "atom $type $atom $xyz"
1682        if {[lsearch $matchedatomlist $atom] != -1} {
1683            puts $fp "labeltext $xyz $atom"
1684        }
1685               
1686        set uiso [atominfo $phase $atom Uiso]
1687        # are there anisotropic atoms? If so convert them to Uequiv
1688        if {[atominfo $phase $atom temptype] == "A"} {
1689            puts -nonewline $fp "Uij [atominfo $phase $atom type] $atom "
1690            foreach v {U11 U22 U33 U12 U13 U23} {
1691                puts -nonewline $fp "[atominfo $phase $atom $v] "
1692            }
1693            puts $fp ""
1694        }
1695    }
1696
1697    foreach type [array names typelist] color {Green Blue Magenta Cyan} {
1698        if {$type == ""} break
1699        puts $fp "sphere $type 0.1 $color"
1700    }
1701    foreach type [array names typelist] color1 {Green Blue Magenta Cyan} {
1702        foreach bondpair $bondlist {
1703            foreach {b1 b2 color} $bondpair {}
1704            if {$color == ""} {set color $color1}
1705            puts $fp "bond $type $type 0.02 $b1 $b2 $color"
1706        }
1707        foreach type1 [array names typelist] {
1708            if {$type1 == $type} break
1709            foreach bondpair $bondlist {
1710                foreach {b1 b2 color} $bondpair {}
1711                if {$color == ""} {set color $color1}
1712                puts $fp "bond $type $type1 0.02 $b1 $b2 $color"
1713            }
1714        }
1715    }
1716    # plot the rigid body
1717    puts $fp "frame"
1718    puts -nonewline $fp "cell"
1719    foreach p {a b c alpha beta gamma} {
1720        puts -nonewline $fp " [phaseinfo $phase $p]"
1721    }
1722    puts $fp ""
1723    puts $fp "background White"
1724    #puts $fp "nolabels"
1725    puts $fp "labelscale 0.5"
1726    puts $fp "spgr P 1"
1727    puts $fp "pack $rbrange"
1728    set i 0
1729    foreach xyz $RBcoords {
1730        foreach {x y z} $xyz {}
1731        incr i
1732        puts $fp "atom c $i $x $y $z"
1733        puts $fp "labeltext $x $y $z r$i"
1734    }
1735    foreach bondpair $rbbondlist {
1736        foreach {b1 b2 color} $bondpair {}
1737        if {$color == ""} {set color Red}
1738        puts $fp "bond c c 0.02 $b1 $b2 $color"
1739    }
1740
1741    puts $fp "sphere c 0.05 Red"
1742    puts $fp "finish   0.70   0.30   0.08   0.01"
1743    puts $fp "end"
1744
1745    #puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
1746    close $fp
1747    MyMessageBox -parent . -title "Info" \
1748                -message "Note that the phase is drawn in green, blue, cyan & magenta and the rigid body in red."
1749    if {$app != ""} {exec $app $file &}
1750}
1751
1752
1753#AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
1754#puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"]
1755#puts [GetRB 1 4 8 "1" "X 1 2" "Z 3 4"]
1756#MapRigidBody 1 1 7 ".11 .22 .33" "11 12 13"
1757
1758
1759#AddRigidBody {1} { {
1760#    {1 1 1 o} {-1 1 1 o} {1 -1 1 o} {-1 -1 1 o}
1761#    {1 1 -1 o} {-1 1 -1 o} {1 -1 -1 o} {-1 -1 -1 o}
1762#} }
1763#set n [MapRigidBody 1 1 1 ".2 .3 .4" "13 17 19"]
1764#puts "body $n created"
1765#incr expgui(changed)
1766#RunRecalcRBCoords
1767#puts "press Enter to continue"
1768#gets stdin line
1769#MapRigidBody 1 1 $n ".5 .5 .5" "0 0 0"
1770#incr expgui(changed)
1771#RunRecalcRBCoords
1772
1773#puts "Test FitBody"
1774set fraclist {
1775    { 0.5483305238484277 0.4887545024531055 0.6167996784631056 }
1776    { 0.1036801409356145 0.5954016321779562 0.5129448102437683 }
1777    { 0.26404665760133855 0.09455414439078394 0.612655365147539 }
1778    { -0.18060372531147473 0.20120127411563465 0.5088004969282018 }
1779    { 0.5806037253114747 0.3987987258843653 0.2911995030717982 }
1780    { 0.13595334239866147 0.5054458556092161 0.18734463485246095 }
1781    { 0.2963198590643855 0.004598367822043814 0.2870551897562318 }
1782    { -0.1483305238484277 0.1112454975468945 0.1832003215368945 }
1783}
1784set ortholist {
1785    {1 1 1} 
1786    {-1 1 1}
1787    {      1.000000   -1.000000    1.000000}
1788    {     -1.000000   -1.000000    1.000000}
1789    {      1.000000    1.000000   -1.000000}
1790    {     -1.000000    1.000000   -1.000000}
1791    {      1.000000   -1.000000   -1.000000}
1792    {     -1.000000   -1.000000   -1.000000} 
1793}
1794# test code, generates DRAWxtl imput file from orthogonal coordinate list
1795# with bonds of ~2, 2.8 and 3.4 A
1796#DRAWxtlPlotOrtho test4.str "test file" $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
1797
1798# test code, plots rigid body type #2 with bonds drawn at ~1.3 & 2 A
1799#PlotRBtype 2 {{1.9 2.1} {1.28 1.32}}
1800
1801# test code, plots rigid body coords in ortholist with bonds @  ~2, 2.8 and 3.4 A
1802#PlotRBcoords $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
1803
1804
1805set useflag {1 1 1 1 1 1 1 1}
1806set cell {4.  5.  6.  95.  100.  105.}
1807#set origin ".20 .30 .40"
1808set origin ".0 .0 .0"
1809#set Euler  {{1 13} {2 17} {3 19}}
1810#set Euler  {{1 0} {2 180} {3 0}}
1811set Euler  {{1 0} {2 0} {3 0}}
1812
1813#puts [La::show $xform]
1814#puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
1815
1816
1817# test zmat2coord
1818set atmlist {
1819    {C1 0 0.0 0 0.0 0 0.0}
1820    {O2 1 1.20 0 0.0 0 0.0}
1821    {H3 1 1.10 2 120.0 0 0.0}
1822    {C4 1 1.50 2 120.0 3 180.0}
1823    {H5 4 1.10 1 110.0 2 0.00}
1824    {H6 4 1.10 1 110.0 2 120.0}
1825    {H7 4 1.10 1 110.0 2 -120.0}
1826}
1827#  C        0.00000        0.00000        0.00000
1828#  O        1.20000        0.00000        0.00000
1829#  H       -0.55000        0.95263        0.00000
1830#  C       -0.75000       -1.29904       -0.00000
1831#  H       -0.04293       -2.14169       -0.00000
1832#  H       -1.38570       -1.36644        0.89518
1833#  H       -1.38570       -1.36644       -0.89518
1834# set coordlist [zmat2coord $atmlist]
1835 set i 0
1836# puts "\nZmatrix in"
1837# foreach line $atmlist {
1838#     incr i
1839#     puts "$i) $line"
1840# }
1841# puts "Cartesian out"
1842# foreach line $coordlist {
1843#     puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
1844# }
1845
1846# AddRigidBody {1 0.75} {
1847#     {
1848#       {1 1 1 c}
1849#       {-1 1 1 c}
1850#       {      1.000000   -1.000000    1.000000 c}
1851#       {     -1.000000   -1.000000    1.000000 c}
1852#       {      1.000000    1.000000   -1.000000 c}
1853#       {     -1.000000    1.000000   -1.000000 c}
1854#       {      1.000000   -1.000000   -1.000000 c}
1855#       {     -1.000000   -1.000000   -1.000000 c}
1856#       {1 1 1 h}
1857#       {1 -1 -1 h}
1858#       {-1 1 -1 h}
1859#       {-1 -1 1 h}
1860#     } {
1861#       {0 0 0 c }
1862#       {0 0 0 c}
1863#       {0 0 0 c}
1864#       {0 0 0 c}
1865#       {0 0 0 c}
1866#       {0 0 0 c}
1867#       {0 0 0 c}
1868#       {0 0 0 c}
1869#       {1 1 1 h}
1870#       {1 -1 -1 h}
1871#       {-1 1 -1 h}
1872#       {-1 -1 1 h}
1873#     }
1874# }
1875# MapRigidBody 2 2 1 {0 0 0} {10 15 20}
Note: See TracBrowser for help on using the repository browser.