source: branches/sandbox/rb.tcl @ 1127

Last change on this file since 1127 was 1127, checked in by toby, 10 years ago

fix TLS reporting

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