source: branches/sandbox/rb.tcl @ 1108

Last change on this file since 1108 was 1108, checked in by toby, 11 years ago

relabel xform atoms; more RB stuff

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