source: branches/sandbox/rb.tcl @ 1124

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

remove debug code

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