source: branches/sandbox/rb.tcl @ 1106

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

more fixes to validreal for small nums; add more RB control routines

File size: 38.6 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 $used
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 $x $lbl]
67        }
68        lappend out [list $mult $damp $var $coordlist]
69    }
70    return $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            setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
399        }
400    }
401}
402
403# convert a decimal to the GSAS hex encoding with a field $digits long.
404proc ToHex {num digits} {
405    return [string toupper [format "%${digits}x" $num]]
406}
407
408# convert a GSAS hex encoding to a decimal integer
409proc FromHex {hex} {
410    return [scan $hex "%x"]
411}
412
413# MapRigidBody: define an "instance" of a rigid body: meaning that the coordinates
414# (and optionally U values) for a set of atoms will be generated from the rigid body
415# arguments:
416#   phase: phase number (1-9)
417#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
418#   firstatom: sequence number of the first atom in phase (note that atoms may
419#              not be numbered sequentially)
420#   position: list of three fractional coordinates for the origin of the rigid body coordinates
421#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
422#           cartesian system before the body is translated to position.
423# returns the instance # (number of times body $bodytyp has been used in phase $phase)
424proc MapRigidBody {phase bodytyp firstatom position angles} {
425    # number of bodies of this type in this phase
426    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
427    set n [string trim [readexp "$key  NBDS"]]
428    if {$n == ""} {
429        makeexprec "$key  NBDS"
430        set n 0
431    }
432    incr n
433    set value $n
434    validint value 5
435    setexp "$key  NBDS" $value 1 5
436    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $n 1]"
437    set value $firstatom
438    validint value 5
439    makeexprec "$key  NDA"
440    setexp "$key  NDA" $value 1 5
441    set l1 {}
442    set l2 {}
443    for {set i 0} {$i < 9} {incr i} {
444        append l1 [format %5d 0]
445        append l2 [format %1d 0]
446    }
447    makeexprec "$key BDFL"
448    setexp "$key BDFL" $l1$l2 1 54
449    makeexprec "${key} BDLC"
450    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
451    makeexprec "${key} BDOR"
452    set l1 {}
453    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
454        append l1 [format "%8.2f%2d" $val $dir]
455    }
456    setexp "${key} BDOR" $l1 1 60
457    makeexprec "${key} LSTF"
458    setexp "${key} LSTF" [format "%5d" 0] 1 5
459    return $n
460}
461
462# EditRigidBodyMapping: edit parameters that define an "instance" of a rigid body (see MapRigidBody)
463# arguments:
464#   phase: phase number (1-9)
465#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
466#   bodynum: instance number, as returned by MapRigidBody
467#   position: list of three fractional coordinates for the origin of the rigid body coordinates
468#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
469#           cartesian system before the body is translated to position.
470#
471proc EditRigidBodyMapping {phase bodytyp bodynum position angles} {
472    # number of bodies of this type in this phase
473    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $bodynum 1]"
474    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
475    set l1 {}
476    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
477        append l1 [format "%8.2f%2d" $val $dir]
478    }
479    setexp "${key} BDOR" $l1 1 60
480}
481#============================================================================
482#============================================================================
483# Returns a list of the variable numbers used already for rigid body variable parameters
484proc RigidBodyGetVarNums {} {
485    set varlist {}
486    set bodies [RigidBodyCount]
487    for {set type 1} {$type <= $bodies} {incr type} {
488        foreach phase $::expmap(phaselist) {
489            set count [RigidBodyMappingCount $phase $type]
490            for {set i 1} {$i <= $bodies} {incr i} {
491                set items [ReadRigidBodyMapping $phase $type $i]
492                set varlist [concat $varlist [lindex $items 3]]
493                if {[llength [lindex $items 6]] > 0} {
494                    set varlist [concat $varlist [lindex $items 6]]
495                }
496            }
497        }
498    }
499    return [lsort -unique $varlist]
500}
501
502# Use the GSAS geometry program to compute a set of cartesian coordinates for a
503# set of atoms in a .EXP file and provide the origin shift and Euler angles needed to
504# place the cartesian system into the crystal coordinates. Used for setting up a rigid body.
505# Returns a nested list of lists:
506#   element 0: a list of the origin location {x y z} in fraction coordinates
507#   element 1: a list of three rotation angles in form {a1 a2 a3}
508#              where a1, a2 and a3 are rotations around the cartesian x, y and z axes
509#   element 2: a list of $natom cartesian coordinate triples {{x1 y1 z1} {x2 y2 z2}...}
510# arguments:
511    # phase: phase #
512    # natom: number of atoms in group
513    # firstatom: sequence # in phase (may be > than number of the atom)
514    # originlist: atoms to define origin (where 1 is first atom in group; <= natom)
515    # vector1: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)
516    # vector2: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)
517    # note that vector2 must define a different axis than vector1
518    # also and vector1 and vector2 cannot use the same atom pair
519proc ExtractRigidBody {phase natom firstatom originlist vector1 vector2} {
520    global expgui
521    set fp [open "geom.inp" "w"]
522    puts $fp "N"
523    puts $fp "N"
524    puts $fp $phase
525    puts $fp "N"
526
527    puts $fp "R"
528    puts $fp "$natom"
529    puts $fp "$firstatom"
530    puts $fp [llength $originlist]
531    foreach i $originlist {
532        puts $fp $i
533    }
534    foreach i [concat $vector1 $vector2] {
535        puts $fp $i
536    }
537    puts $fp "0"
538    puts $fp "X"
539    close $fp
540    catch {
541        exec [file join $expgui(gsasexe) geometry] $expgui(expfile) < geom.inp > geom.out
542    }
543    file delete geom.inp
544    set fp [open geom.out r]
545    while {[gets $fp line] >= 0} {
546        if {[string first "Cell coordinates of origin" $line] != -1} {
547            set origin [lrange [string range $line [string first "are" $line] end] 1 3]
548        }
549        if {[string first "Rotation angles" $line] != -1} {
550            set Euler {}
551            foreach i [lrange [split $line "="] 1 3] {
552                lappend Euler [lindex $i 0]
553            }
554            #puts $line
555            #puts $Euler
556        }
557        if {[string first "Atom   Orthon" $line] != -1} {
558            set coordlist {}
559            for {set i 1} {$i <= $natom} {incr i} {
560                gets $fp line
561                set coord {}
562                lappend coord [string trim [string range $line 9 15]]
563                lappend coord [string trim [string range $line 16 22]]
564                lappend coord [string trim [string range $line 23 29]]
565                lappend coord [string trim [string range $line 0 8]]
566                #puts $line
567                #puts $coord
568                lappend coordlist $coord
569            }
570            #puts $coordlist
571        }
572    }
573    #file delete geom.out
574    return [list $origin $Euler $coordlist]
575}
576
577# updates the coordinates in a .EXP file after a rigid body instance has been added/changed
578proc RunRecalcRBCoords { } {
579    global expgui tcl_platform
580    set input [open resetmult.inp w]
581    puts $input "Y"
582    puts $input "l b"
583    puts $input "n"
584    puts $input "x x x"
585    puts $input "x"
586    close $input
587    # Save the current exp file
588    savearchiveexp
589    # disable the file changed monitor
590    set expgui(expModifiedLast) 0
591    set expnam [file root [file tail $expgui(expfile)]]
592    set err [catch {
593        if {$tcl_platform(platform) == "windows"} {
594            exec [file join $expgui(gsasexe) expedt.exe] $expnam < resetmult.inp >& resetmult.out
595        } else {
596            exec [file join $expgui(gsasexe) expedt] $expnam < resetmult.inp >& resetmult.out
597        }
598    } errmsg]
599    loadexp $expgui(expfile)
600    set fp [open resetmult.out r]
601    set out [read $fp]
602    close $fp
603    set expgui(exptoolout) $out
604    catch {file delete resetmult.inp resetmult.out}
605    if {$err} {
606        return $errmsg
607    } else {
608        return ""
609    }
610}
611
612
613# compute a rotation matrix for orthogonal coordinates (based on MAKMATD in GSAS)
614# rotate angle degrees around axis (1, 2 or 3) for (x, y, or z)
615# returns a list that can be used as a matrix in the La package
616proc RotMatrix {axis angle} {
617    set ang [expr {$angle * acos(0) / 90.}]
618    set mat "1 0 0 0 1 0 0 0 1"
619    if {$axis == 1}  {
620        set i1 1
621        set i2 2
622    } elseif {$axis == 2}  {
623        set i1 2
624        set i2 0
625    } else {
626        set i1 0
627        set i2 1
628    }
629    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
630    foreach item {
631        {$i1 $i1 [expr {cos($ang)}]}
632        {$i2 $i2 [expr {cos($ang)}]}
633        {$i1 $i2 [expr {-sin($ang)}]}
634        {$i2 $i1 [expr {sin($ang)}]}
635    } {
636        foreach {c r val} [subst $item] {}
637        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
638    }
639    return "2 3 3 $mat"
640}
641
642# compute the derivative of the rotation matrix with respect to the angle, see RotMatrix
643# (based on MAKMATD in GSAS)
644# returns a list that can be used as a matrix in the La package
645proc DerivRotMatrix {axis angle} {
646    set ang [expr {$angle * acos(0) / 90.}]
647    set mat "0 0 0 0 0 0 0 0 0"
648    if {$axis == 1}  {
649        set i1 1
650        set i2 2
651    } elseif {$axis == 2}  {
652        set i1 2
653        set i2 0
654    } else {
655        set i1 0
656        set i2 1
657    }
658    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
659    foreach item {
660        {$i1 $i1 [expr {-sin($ang) * acos(0) / 90.}]}
661        {$i2 $i2 [expr {-sin($ang) * acos(0) / 90.}]}
662        {$i1 $i2 [expr {-cos($ang) * acos(0) / 90.}]}
663        {$i2 $i1 [expr {cos($ang) * acos(0) / 90.}]}
664    } {
665        foreach {c r val} [subst $item] {}
666        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
667    }
668    return "2 3 3 $mat"
669}
670
671# compute an orthogonalization matrix from cell parameters (based on AMATRX in GSAS)
672# returns a list that can be used as a matrix in the La package
673proc OrthoMatrix {a b c alpha beta gamma} {
674    set CA [expr {cos($alpha * acos(0) / 90.)}]
675    set CB [expr {cos($beta * acos(0) / 90.)}]
676    set CG [expr {cos($gamma * acos(0) / 90.)}]
677    set SA [expr {sin($alpha * acos(0) / 90.)}]
678    set SB [expr {sin($beta * acos(0) / 90.)}]
679    set SC [expr {sin($gamma * acos(0) / 90.)}]
680    set CASTAR [expr { ($CB*$CG-$CA)/($SB*$SC) }]    ;#! cos(Alpha*)
681    set CBSTAR [expr { ($CA*$CG-$CB)/($SA*$SC) }]    ;#! cos(Beta*)
682    set CCSTAR [expr { ($CA*$CB-$CG)/($SA*$SB) }]    ;#! cos(Gamma*)
683    set SASTAR [expr { sqrt(1.0-($CASTAR*$CASTAR*2)) }]    ;#! sin(Alpha*)
684    set SBSTAR [expr { sqrt(1.0-($CBSTAR*$CBSTAR*2)) }]    ;#! sin(Beta*)
685    set SCSTAR [expr { sqrt(1.0-($CCSTAR*$CCSTAR*2)) }]    ;#! sin(Gamma*)
686
687    set A  "2 3 3      $a 0 0    0 $b 0    0 0 $c"
688    set A1 "2 3 3     1.0 0 0    $CG [expr {$SASTAR*$SC}] [expr {-$CASTAR*$SC}]    $CB 0.0 $SB"
689    #!This matrix is
690    #!   (1.0      0.0            0.0             )
691    #!   (cos(G)  sin(A*)*sin(G) -cos(A*)*sin(G)  )
692    #!   (cos(B)   0.0            sin(B)          )
693    return [La::transpose [La::mmult $A $A1]]
694}
695
696# compute the transformation matrix that converts a rigid body coordinates into fractional
697# coordinates
698# arguments:
699#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
700#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
701#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
702# returns a list that can be used as a matrix in the La package
703proc CalcXformMatrix {rotations cellprms} {
704    set prod {}
705    foreach item $rotations {
706        #puts $item
707        set mat [eval RotMatrix $item]
708        if {$prod == ""} {
709            set prod $mat
710        } else {
711            set prod [La::mmult $prod $mat]
712        }
713    }
714    #puts "--- rotation product ---"
715    #puts [La::show $prod]
716
717    set ortho [eval OrthoMatrix $cellprms]
718    #puts "--- ortho ---"
719    #puts [La::show $ortho]
720    set deortho [La::msolve $ortho [La::mident 3] ]
721    #puts "--- deortho ---"
722    #puts [La::show $deortho]
723    #puts "--- xform ---"
724    set xform [La::mmult $deortho $prod]
725    return $xform
726}
727
728# transforms a triplet of orthogonal coordinates into fractional ones using
729# arguments:
730#    xform: a transformation matrix from CalcXformMatrix
731#    origin: a list of fraction coordinates {x y z} describing the location of the
732#            origin of the orthogonal coordinates in the crystal system
733#    ortho: a triplet of othogonal coordinates
734# returns a triplet of fractional coordinates
735proc Ortho2Xtal {xform origin ortho} {
736    set ocv "2 3 0 $ortho"
737    set frac [La::mmult $xform $ocv]
738    #puts [La::show [La::transpose $frac]]
739    #puts $frac
740    set frac [La::madd $frac "[lrange $frac 0 2] $origin"]
741    #puts [La::show [La::transpose $frac]]
742    return $frac
743}
744
745# compute the derivative of the transformation matrix (see CalcXformMatrix)
746# with respect to every rotation in the $rotations list
747# arguments:
748#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
749#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
750#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
751# returns a list of matrices where each matrix is a list that can be used as a
752# matrix in the La package
753proc CalcDerivMatrix {rotations cellprms} {
754    set ortho [eval OrthoMatrix $cellprms]
755    set deortho [La::msolve $ortho [La::mident 3] ]
756    set derivlist {}
757
758    foreach item $rotations {
759        #puts $item
760        set mat [eval DerivRotMatrix $item]
761        #puts $item
762        #puts [La::show $mat]
763        set xform [La::mmult $deortho $mat]
764        lappend derivlist $xform
765    }
766    return $derivlist
767}
768
769# fit a rigid body's origin
770# arguments:
771#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
772#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
773#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
774#  ortholist: list containing triplets with orthogonal coordinates
775#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
776#  fraclist: list containing triplets with fractional coordinates 
777#  origin: a list of fraction coordinates {x y z} describing the location of the
778#            origin of the orthogonal coordinates in the crystal system
779#     note that the length of ortholist, useflag and fraclist should be the same
780# Returns a list with the following elements
781#   0: a list with three new origin values
782#   1: the root-mean square difference between the fraclist coordinates and those computed
783#      using the input values for those atoms where $use is one (in Angstroms)
784#   2: the computed fractional coordinates for all atoms
785#   3: individual rms values for all atoms (in Angstroms)
786# note that items 1-3 are computed with the imput origin, not the revised one
787proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} {
788    set xform [CalcXformMatrix $Euler $cell]
789    foreach var {x y z} {set sum($var) 0.0}
790
791    set i 0
792    set sumdvs 0
793    set fracout {}
794    set rmsout {}
795    foreach oc $ortholist use $useflag coord $fraclist {
796        #puts "ortho: $oc"
797        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
798        lappend fracout $frac
799        if {$use} {incr i}
800        set dvs 0
801        foreach var {x y z} v1 $frac v2 $coord abc [lrange $cell 0 2] {
802            set dv [expr {($v2 - $v1)}]
803            set dvs [expr {$dvs + $dv*$dv}]
804            set sumdvs [expr {$sumdvs + $dv*$dv}]
805            if {$use} {set sum($var) [expr {$sum($var) + $dv}]}
806        }
807        lappend rmsout [expr {sqrt($dvs)}]
808    }
809    set rms 0
810    if {$i > 1} {set rms [expr {sqrt($sumdvs)/$i}]}
811    set neworig {}
812    foreach var {x y z} o $origin {
813        lappend neworig [expr {$o + ($sum($var)/$i)}]
814    }
815    return [list $neworig $rms $fracout $rmsout]
816}
817
818# fit a rigid body's Euler angles using least-squares
819# arguments:
820#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
821#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
822#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
823#  ortholist: list containing triplets with orthogonal coordinates
824#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
825#  fraclist: list containing triplets with fractional coordinates 
826#  origin: a list of fraction coordinates {x y z} describing the location of the
827#            origin of the orthogonal coordinates in the crystal system
828#     note that the length of ortholist, useflag and fraclist should be the same
829# Returns a list of new Euler angles
830proc FitBodyRot {Euler cell ortholist useflag fraclist origin} {
831    set xform [CalcXformMatrix $Euler  $cell]
832    set derivlist [CalcDerivMatrix $Euler  $cell]
833    set A "2 [expr 3*[llength $ortholist]] 3"
834    foreach oc $ortholist use $useflag coord $fraclist {
835        if {! $use} continue
836        foreach deriv $derivlist {
837            foreach xyz [lrange [Ortho2Xtal $deriv "0 0 0" $oc] 3 end] {
838                lappend A $xyz
839            }
840        }
841    }
842    #puts "A: [La::show $A]"
843    set y "2 [expr 3*[llength $ortholist]] 1"
844    foreach oc $ortholist use $useflag coord $fraclist {
845        if {! $use} continue
846        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
847        foreach xyz $coord XYZ $frac {
848            lappend y [expr {$XYZ - $xyz}]
849        }
850    }
851
852    set AtA [La::mmult [La::transpose $A] $A]
853    set Aty [La::mmult [La::transpose $A] $y]
854   
855    set l {}
856    #set shifts {}
857    foreach delta [lrange [La::msolve $AtA $Aty] 3 end] item $Euler {
858        #lappend shifts $delta
859        lappend l "[lindex $item 0] [expr {$delta + [lindex $item 1]}]"
860    }
861    #puts "shifts = $shifts"
862    return $l
863}
864
865# fit a rigid body's Origin and Euler angles
866# arguments:
867#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
868#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
869#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
870#  ortholist: list containing triplets with orthogonal coordinates
871#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
872#  fraclist: list containing triplets with fractional coordinates 
873#  origin: a list of fraction coordinates {x y z} describing the location of the
874#            origin of the orthogonal coordinates in the crystal system
875#     note that the length of ortholist, useflag and fraclist should be the same
876# Returns a list containing
877#   new origin
878#   new Euler angles
879#   total rms
880#   fractional coordinates
881#   rms deviation in fractional coordinates of new Euler angles
882proc FitBody {Euler cell ortholist useflag fraclist origin "ncycle 5"} {
883    puts "start origin = $origin"
884    foreach {
885        origin 
886        startrms
887        fracout
888        rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
889    puts "start rms = $startrms"
890    set rmsprev $startrms
891    puts "new origin = $origin"
892    for {set i 0} {$i < $ncycle} {incr i} {
893        set Eulerprev $Euler
894        set Euler [FitBodyRot $Euler $cell $ortholist $useflag $fraclist $origin]
895        puts "New Euler $Euler"
896        puts "after fit" 
897        foreach {
898            origin 
899            rms
900            fracout
901            rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
902        if {$rms > (1.1 * $rmsprev) + 0.01} {
903            puts "rms = $rms, new origin = $origin"
904            set rmsprev $rms
905        }
906    }   
907    #proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} 
908    #return "$neworig $rms $fracout $rmsout"
909    set fmt  {"%8.5f %8.5f %8.5f     %8.5f %8.5f %8.5f   %6.3f"}
910    foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
911        puts "[eval format $fmt $fracin $fraccalc $rmsi]"
912    }
913    return [list $origin $Euler $rms $fracout $rmsout]
914}
915
916# zmat2coord converts a z-matrix to a set of cartesian coordinates
917#   a z-matrix is also known as "internal coordinates" or "torsion space"
918#   (see Journal of Computational Chemistry, Vol 26, #10, p. 1063–1068, 2005 or
919#    http://www.cmbi.ru.nl/molden/zmat/zmat.html)
920# INPUT:
921#   atmlist is a list of ascii lines where each line contains
922#     lbl c1 distance c2 angle c3 torsion
923#   where each atom is computed from the previous where the new atom is:
924#     distance $distance from atom $c1 (angstrom)
925#     angle $angle from $c1--$c2 (degrees)
926#     torsion $torsion from $c1--$c2--$c3 (degrees)
927# OUTPUT:
928#  zmat2coord returns a list of atom labels and cartesian coordinates,
929#  with 4 items in each element (label, x, y, z)
930# this routine was tested against results from Babel via the web interface at
931# http://www.shodor.org/chemviz/zmatrices/babel.html and sample input at
932# http://iopenshell.usc.edu/howto/zmatrix/
933proc zmat2coord {atmlist} { 
934    set torad [expr {acos(0)/90.}]
935    set i 0
936    foreach line $atmlist {
937        incr i
938        foreach {lbl c1 dist c2 angle c3 torsion} $line {} 
939        if {$i == 1} {
940            set atm(1) [list $lbl 0 0 0] ; # 1st atom is at origin
941        } elseif {$i == 2} {
942            set dist1 $dist
943            set atm(2) [list $lbl $dist1 0 0] ; # 2nd atom is along x-axis
944        } elseif {$i == 3} {
945            # 3rd atom can be bonded to the 1st or 2nd
946            if {$c1 == 1} {
947                set atm(3) [list $lbl \
948                                [expr {$dist * cos($torad * $angle)}] \
949                                [expr {$dist * sin($torad * $angle)}] \
950                                0]
951            } else {
952                set atm(3) [list $lbl \
953                                [expr {$dist1 - $dist * cos($torad * $angle)}] \
954                                [expr {$dist * sin($torad * $angle)}] \
955                                0]
956            }
957        } else {
958            set atm($i) [concat $lbl \
959                             [ahcat "atm" $c1 $dist $c2 $angle $c3 $torsion]]
960        }
961    }
962    set coordlist {}
963    foreach key [lsort -integer [array names atm]] {
964        lappend coordlist $atm($key)
965    }
966    return $coordlist
967}
968# Compute the length of a vector
969proc vlen {a} {
970    set sum 0.0
971    foreach ai $a {
972        set sum [expr {$sum + $ai*$ai}]
973    }
974    return [expr sqrt($sum)]
975}
976# compute vector (a + z * b) and optionally normalize to length d
977proc vadd {a b d z} {
978    set c {}
979    foreach ai $a bi $b {
980        lappend c [expr {$bi + $z * $ai}]
981    }
982    set v [vlen $c]
983    if {$d != 0} {
984        set r {}
985        foreach ci $c {
986            lappend r [expr {$d * $ci / $v}]
987        }
988        return [list $v $r]
989    }
990    return [list $v $c]
991}
992# normalize a vector
993proc vnrm {x} {
994    set v [vlen $x]
995    if {abs($v) < 1e-8} {return [list 0 0 0]}
996    set y {}
997    foreach xi $x {
998        lappend y [expr {$xi / $v}]
999    }
1000    return $y
1001}
1002# compute the coordinates for an atom that is bonded:
1003#   distance $dist from atom $nc
1004#   angle $bondang from $nc--$nb
1005#   torsion $torsang from $nc--$nb--$na
1006#   coordinates are found in array $atmarr in the calling routine
1007# based on a Fortran routine provided by Peter Zavalij (Thanks Peter!)
1008proc ahcat {atmarr nc dist nb bondang na torsang} {
1009    upvar 1 $atmarr atm
1010    set xa [lrange $atm($na) 1 3]
1011    set xb [lrange $atm($nb) 1 3]
1012    set xc [lrange $atm($nc) 1 3]
1013    set torad [expr {acos(0)/90.}]
1014    # x = unit Vector A-B
1015    foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. -1.] 1] {}
1016    # y = unit Vector C-B
1017    set y [lindex [vadd $xb $xc 1. -1.] 1]
1018    foreach {y1 y2 y3} $y {}
1019    set z1 [expr {$x2*$y3 - $x3*$y2}]
1020    set z2 [expr {$x3*$y1 - $x1*$y3}]
1021    set z3 [expr {$x1*$y2 - $x2*$y1}]
1022    set z [vnrm [list $z1 $z2 $z3]]
1023    set q1 [expr {$y2*$z3 - $y3*$z2}]
1024    set q2 [expr {$y3*$z1 - $y1*$z3}]
1025    set q3 [expr {$y1*$z2 - $y2*$z1}]
1026    set q [vnrm [list $q1 $q2 $q3]]
1027    set th [expr {$bondang * $torad}]
1028    set ph [expr {-1. * $torsang * $torad}]
1029    set cth [expr {cos($th)}]
1030    set sth [expr {sin($th)}]
1031    set cph [expr {cos($ph)}]
1032    set sph [expr {sin($ph)}]
1033    set xh {}
1034    foreach xci $xc xi $q zi $z yi $y {
1035        lappend xh [expr {
1036                          $xci +
1037                          $dist*($sth*($cph*$xi + $sph*$zi)-$cth*$yi)
1038                      }]
1039    }
1040    return $xh
1041}
1042 
1043
1044
1045#AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
1046#puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"]
1047#puts [GetRB 1 4 8 "1" "X 1 2" "Z 3 4"]
1048#MapRigidBody 1 1 7 ".11 .22 .33" "11 12 13"
1049
1050
1051#AddRigidBody {1} { {
1052#    {1 1 1 o} {-1 1 1 o} {1 -1 1 o} {-1 -1 1 o}
1053#    {1 1 -1 o} {-1 1 -1 o} {1 -1 -1 o} {-1 -1 -1 o}
1054#} }
1055#set n [MapRigidBody 1 1 1 ".2 .3 .4" "13 17 19"]
1056#puts "body $n created"
1057#incr expgui(changed)
1058#RunRecalcRBCoords
1059#puts "press Enter to continue"
1060#gets stdin line
1061#MapRigidBody 1 1 $n ".5 .5 .5" "0 0 0"
1062#incr expgui(changed)
1063#RunRecalcRBCoords
1064
1065puts "Test FitBody"
1066set fraclist {
1067    { 0.5483305238484277 0.4887545024531055 0.6167996784631056 }
1068    { 0.1036801409356145 0.5954016321779562 0.5129448102437683 }
1069    { 0.26404665760133855 0.09455414439078394 0.612655365147539 }
1070    { -0.18060372531147473 0.20120127411563465 0.5088004969282018 }
1071    { 0.5806037253114747 0.3987987258843653 0.2911995030717982 }
1072    { 0.13595334239866147 0.5054458556092161 0.18734463485246095 }
1073    { 0.2963198590643855 0.004598367822043814 0.2870551897562318 }
1074    { -0.1483305238484277 0.1112454975468945 0.1832003215368945 }
1075}
1076set ortholist {
1077    {1 1 1} 
1078    {-1 1 1}
1079    {      1.000000   -1.000000    1.000000}
1080    {     -1.000000   -1.000000    1.000000}
1081    {      1.000000    1.000000   -1.000000}
1082    {     -1.000000    1.000000   -1.000000}
1083    {      1.000000   -1.000000   -1.000000}
1084    {     -1.000000   -1.000000   -1.000000} 
1085}
1086
1087set useflag {1 1 1 1 1 1 1 1}
1088set cell {4.  5695100105.}
1089#set origin ".20 .30 .40"
1090set origin ".0 .0 .0"
1091#set Euler  {{1 13} {2 17} {3 19}}
1092#set Euler  {{1 0} {2 180} {3 0}}
1093set Euler  {{1 0} {2 0} {3 0}}
1094
1095#puts [La::show $xform]
1096puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
1097
1098
1099# test zmat2coord
1100set atmlist {
1101    {C1 0 0.0 0 0.0 0 0.0}
1102    {O2 1 1.20 0 0.0 0 0.0}
1103    {H3 1 1.10 2 120.0 0 0.0}
1104    {C4 1 1.50 2 120.0 3 180.0}
1105    {H5 4 1.10 1 110.0 2 0.00}
1106    {H6 4 1.10 1 110.0 2 120.0}
1107    {H7 4 1.10 1 110.0 2 -120.0}
1108}
1109#  C        0.00000        0.00000        0.00000
1110#  O        1.20000        0.00000        0.00000
1111#  H       -0.55000        0.95263        0.00000
1112#  C       -0.75000       -1.29904       -0.00000
1113#  H       -0.04293       -2.14169       -0.00000
1114#  H       -1.38570       -1.36644        0.89518
1115#  H       -1.38570       -1.36644       -0.89518
1116set coordlist [zmat2coord $atmlist]
1117set i 0
1118puts "\nZmatrix in"
1119foreach line $atmlist {
1120    incr i
1121    puts "$i) $line"
1122}
1123puts "Cartesian out"
1124foreach line $coordlist {
1125    puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
1126}
Note: See TracBrowser for help on using the repository browser.