source: branches/sandbox/rb.tcl @ 1107

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

fix rb bugs: ReadRigidBody? needs nested list & SetRigidBody? writes wrong # of translations; Add plot RB capability

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