source: branches/sandbox/rb.tcl @ 1102

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

make rb tcl8.4 compatible

File size: 32.3 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 (to move into readexp.tcl)
7#============================================================================
8# returns the number of defined rigid bodies
9proc RigidBodyCount {} {
10    set n [string trim [readexp "RGBD  NRBDS"]]
11    if {$n == ""} {
12        set n 0
13    }
14    return $n
15}
16
17# returns two items:
18#   the number of times the rigid body is used
19#   a list containing an element for each multiplier in rigid body #rbnum.
20# in each element there are four items:
21#    a multiplier value for the rigid body coordinates
22#    a damping value (0-9) for the refinement of the multiplier
23#    a variable number if the multiplier will be refined
24#    a list of cartesian coordinates coordinates
25# each cartesian coordinate contains 4 items: x,y,z and a label
26#  note that the label is present only when the RB is created in EXPGUI and is
27#  not used in GSAS.
28proc ReadRigidBody {rbnum} {
29    if {[RigidBodyCount] > $rbnum} {
30        return ""
31    }
32    set value $rbnum
33    validint value 2
34    set key "RGBD${value}"
35    set n [string trim [string range [readexp "$key NATR"] 0 4]]
36    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
37    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
38    set out $used
39    for {set i 1} {$i <= $nmult} {incr i} {
40        set line [readexp "${key}${i}PARM"]
41        set mult [string trim [string range $line 0 9]]
42        set var [string trim [string range $line 10 14]]
43        set damp [string trim [string range $line 15 19]]
44        set coordlist {}
45        for {set j 1} {$j <= $n} {incr j} {
46            set value $j
47            validint value 3
48            set line [readexp "${key}${i}SC$value"]
49            set x [string trim [string range $line 0 9]]
50            set y [string trim [string range $line 10 19]]
51            set z [string trim [string range $line 20 29]]
52            set lbl [string trim [string range $line 30 39]]
53            lappend coordlist [list $x $y $x $lbl]
54        }
55        lappend out [list $mult $damp $var $coordlist]
56    }
57    return $out
58}
59
60proc RigidBodyMappingCount {phase bodytyp} {
61    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
62    set n [string trim [readexp "$key  NBDS"]]
63    if {$n == ""} {
64        set n 0
65    }
66    return $n
67}
68
69# reads rigid body mapping for phase, body # and instance number
70# returns a list of items (most lists) as follows:
71#   1) sequence # of first atom in body
72#   2) origin of body in fractional coordinates (3 elements)
73#   3) Euler angles as 6 pairs of numbers (see below)
74#   4) variable numbers for the 9 position variables (origin followed by rotations)
75#   5) damping vals for the 9 position variables (origin followed by rotations)
76#   6) the TLS values, in order below (empty list if TLS is not in use)
77#   7) the variable numbers for each TLS values, in order below (or empty)
78#   8) three damping values for the T, L and S terms.
79# returns an empty list if no such body exists.
80#
81# Euler angles are a list of axes and angles to rotate:
82#   { {axis1 angle1} {axis2 angle2} ...}
83# where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
84#
85# The 20 TLS terms are ordered:
86#    T11, T22, T33, T12, T13, T23
87#    L11, L22, L33, L12, L13, L23
88#    S12, S13, S21, S23, S31, S32, SAA, SBB
89#
90proc ReadRigidBodyMapping {phase bodytyp num} {
91    if {[RigidBodyMappingCount $phase $bodytyp] < $num} {
92        return ""
93    }
94    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
95    set first [string trim [string range [readexp "$key  NDA"] 0 4]]
96    set line [readexp "$key BDFL"]
97    set varlist {}
98    set damplist {}
99    foreach i {0 1 2 3 4 5 6 7 8} {
100        lappend varlist [string trim [string range $line [expr {5*$i}] [expr {4 + 5*$i}] ]]
101        lappend damplist [string trim [string range $line [expr {45 + $i}] [expr {45 + $i}] ]]
102    }
103    set TLSdamplist {}
104    foreach i {54 55 56} {
105        lappend TLSdamplist [string trim [string range $line $i $i ]]
106    }
107    set line [readexp "${key} BDLC"]
108    set x [string trim [string range $line 0 9]]
109    set y [string trim [string range $line 10 19]]
110    set z [string trim [string range $line 20 29]]
111    set origin [list $x $y $z]
112    set line [readexp "${key} BDOR"]
113    set rotations {}
114    foreach i {0 10 20 30 40 50} {
115        set angle [string trim [string range $line $i [expr {$i+7}]]]
116        set axis [string trim [string range $line [expr {$i+8}] [expr {$i+9}]]]
117        lappend rotations [list $angle $axis]
118    }
119    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
120    set tlsvars {}
121    set tlsvals {}
122    if {$TLS != 0} {
123        set line [readexp "${key}TLSF1"]
124        for {set j 0} {$j < 20} {incr j} {
125            lappend tlsvars [string trim [string range $line [expr {3*$j}] [expr {3*$j+2}]]]
126        }
127        for {set j 0} {$j < 20} {incr j} {
128            set i 0
129            if {$j == 0} {
130                set i 1
131            } elseif {$j == 8} {
132                set i 2
133            } elseif {$j == 16} {
134                set i 3
135            }
136            if {$i != 0} {
137                set line [readexp "${key}TLSP$i"]
138                set i 0
139                set j1 0
140                set j2 7
141            } else {
142                incr j1 8
143                incr j2 8
144            }
145            lappend tlsvals [string trim [string range $line $j1 $j2]]
146        }
147    }
148    return [list $first $origin $rotations $varlist $damplist $tlsvals $tlsvars $TLSdamplist]
149}
150
151#---------------------------
152# routines to write
153# setTLS on/off
154# vary body flags
155# vary TLS flags
156# edit TLS values
157# set body damping
158# set TLS damping
159#---------------------------
160
161
162# AddRigidBody: add a new rigid body definition into the .EXP file
163# arguments are:
164#   multlist: defines a list of multipliers for each set of coordinates. In the
165#             simplest case this will be {1}
166#   coordlist: a nested list of coordinates such as { { {0 0 0} {.1 .1 .1} {.2 .2 .2} } }
167# note that when the length of multlist > 1 then coordlist must have the same length.
168# for input where
169#     multlist = {s1 s2} and
170#     coordlist = { { {0 0 0} {1 1 0} {.0 .0 .0} ...}
171#                     {0 0 0} {1 1 0} {2 1 2} ...}
172#                 }   
173# the cartesian coordinates are defined from the input as
174#    atom 1 = s1 * (0,0,0) + s2*(0,0,0) [ = (0,0,0)]
175#    atom 2 = s1 * (1,1,0) + s2*(1,1,0) [ = (s1+s2) * (1,1,0)]
176#    atom 3 = s1 * (0,0,0) + s2*(2,1,2) [ = s2 * (2,1,2)]
177#    ...
178# Returns the number of the rigid body that has been created
179proc AddRigidBody {multlist coordlist} {
180    #
181    set n [string trim [readexp "RGBD  NRBDS"]]
182    if {$n == ""} {
183        makeexprec "RGBD  NRBDS"
184        set n 0
185    }
186    incr n
187    validint n 5
188    setexp "RGBD  NRBDS" $n 1 5
189    SetRigidBody $n $multlist $coordlist
190    return $n
191}
192
193# ReplaceRigidBody: replace all the information for rigid body #rbnum
194# Works the sames as AddRigidBody (see above) except that the rigid body is replaced rather
195# than added.
196# Note that count of the # of times the body is used is preserved
197proc ReplaceRigidBody {rbnum multlist coordlist} {
198    set value $rbnum
199    validint value 2
200    set key "RGBD${value}"
201    set line [readexp "$key NBDS"]
202    foreach key [array names ::exparray "${key}*"] {
203        #puts $key
204        delexp $key
205    }
206    SetRigidBody $rbnum $multlist $coordlist
207    setexp "$key NBDS" $line 1 68
208}
209
210# Edit the parameters for rigid body #rbnum
211# (normally called from ReplaceRigidBody or AddRigidBody)
212proc SetRigidBody {rbnum multlist coordlist} {
213    set value $rbnum
214    validint value 2
215    set key "RGBD${value}"
216    # number of atoms
217    set value [llength [lindex $coordlist 0]]
218    validint value 5
219    makeexprec "$key NATR"
220    setexp "$key NATR" $value 1 5
221    # number of times used
222    set value 0
223    validint value 5
224    makeexprec "$key NBDS"
225    setexp "$key NBDS" $value 1 5
226    # number of coordinate matrices
227    set value [llength multlist]
228    validint value 5
229    makeexprec "$key NSMP"
230    setexp "$key NSMP" $value 1 5
231    set i 0
232    foreach mult $multlist coords $coordlist {
233        incr i
234        makeexprec "${key}${i}PARM"
235        setexp "${key}${i}PARM" [format "%10.5f%5d%5d" $mult 0 0] 1 20
236        set j 0
237        foreach item $coords {
238            #puts $item
239            incr j
240            set value $j
241            validint value 3
242            makeexprec "${key}${i}SC$value"
243            setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
244        }
245    }
246}
247
248# convert a decimal to the GSAS hex encoding with a field $digits long.
249proc ToHex {num digits} {
250    return [string toupper [format "%${digits}x" $num]]
251}
252
253# convert a GSAS hex encoding to a decimal integer
254proc FromHex {hex} {
255    return [scan $hex "%x"]
256}
257
258# ApplyRigidBody: define an "instance" of a rigid body: meaning that the coordinates
259# (and optionally U values) for a set of atoms will be generated from the rigid body
260# arguments:
261#   phase: phase number (1-9)
262#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
263#   firstatom: sequence number of the first atom in phase (note that atoms may
264#              not be numbered sequentially)
265#   position: list of three fractional coordinates for the origin of the rigid body coordinates
266#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
267#           cartesian system before the body is translated to position.
268# returns the instance # (number of times body $bodytyp has been used in phase $phase)
269proc ApplyRigidBody {phase bodytyp firstatom position angles} {
270    # number of bodies of this type in this phase
271    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
272    set n [string trim [readexp "$key  NBDS"]]
273    if {$n == ""} {
274        makeexprec "$key  NBDS"
275        set n 0
276    }
277    incr n
278    set value $n
279    validint value 5
280    setexp "$key  NBDS" $value 1 5
281    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $n 1]"
282    set value $firstatom
283    validint value 5
284    makeexprec "$key  NDA"
285    setexp "$key  NDA" $value 1 5
286    set l1 {}
287    set l2 {}
288    for {set i 0} {$i < 9} {incr i} {
289        append l1 [format %5d 0]
290        append l2 [format %1d 0]
291    }
292    makeexprec "$key BDFL"
293    setexp "$key BDFL" $l1$l2 1 54
294    makeexprec "${key} BDLC"
295    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
296    makeexprec "${key} BDOR"
297    set l1 {}
298    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
299        append l1 [format "%8.2f%2d" $val $dir]
300    }
301    setexp "${key} BDOR" $l1 1 60
302    makeexprec "${key} LSTF"
303    setexp "${key} LSTF" [format "%5d" 0] 1 5
304    return $n
305}
306
307# EditRigidBody: edit parameters that define an "instance" of a rigid body (see ApplyRigidBody)
308# arguments:
309#   phase: phase number (1-9)
310#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
311#   bodynum: instance number, as returned by ApplyRigidBody
312#   position: list of three fractional coordinates for the origin of the rigid body coordinates
313#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
314#           cartesian system before the body is translated to position.
315#
316proc EditRigidBody {phase bodytyp bodynum position angles} {
317    # number of bodies of this type in this phase
318    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $bodynum 1]"
319    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
320    set l1 {}
321    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
322        append l1 [format "%8.2f%2d" $val $dir]
323    }
324    setexp "${key} BDOR" $l1 1 60
325}
326#============================================================================
327#============================================================================
328
329# Use the GSAS geometry program to compute a set of cartesian coordinates for a
330# set of atoms in a .EXP file and provide the origin shift and Euler angles needed to
331# place the cartesian system into the crystal coordinates. Used for setting up a rigid body.
332# Returns a nested list of lists:
333#   element 0: a list of the origin location {x y z} in fraction coordinates
334#   element 1: a list of three rotation angles in form {a1 a2 a3}
335#              where a1, a2 and a3 are rotations around the cartesian x, y and z axes
336#   element 2: a list of $natom cartesian coordinate triples {{x1 y1 z1} {x2 y2 z2}...}
337# arguments:
338    # phase: phase #
339    # natom: number of atoms in group
340    # firstatom: sequence # in phase (may be > than number of the atom)
341    # originlist: atoms to define origin (where 1 is first atom in group; <= natom)
342    # vector1: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)
343    # vector2: list of 3 values with X, Y or Z, atom #a and #b (number as in origin)
344    # note that vector2 must define a different axis than vector1
345    # also and vector1 and vector2 cannot use the same atom pair
346proc ExtractRigidBody {phase natom firstatom originlist vector1 vector2} {
347    global expgui
348    set fp [open "geom.inp" "w"]
349    puts $fp "N"
350    puts $fp "N"
351    puts $fp $phase
352    puts $fp "N"
353
354    puts $fp "R"
355    puts $fp "$natom"
356    puts $fp "$firstatom"
357    puts $fp [llength $originlist]
358    foreach i $originlist {
359        puts $fp $i
360    }
361    foreach i [concat $vector1 $vector2] {
362        puts $fp $i
363    }
364    puts $fp "0"
365    puts $fp "X"
366    close $fp
367    catch {
368        exec [file join $expgui(gsasexe) geometry] $expgui(expfile) < geom.inp > geom.out
369    }
370    file delete geom.inp
371    set fp [open geom.out r]
372    while {[gets $fp line] >= 0} {
373        if {[string first "Cell coordinates of origin" $line] != -1} {
374            set origin [lrange [string range $line [string first "are" $line] end] 1 3]
375        }
376        if {[string first "Rotation angles" $line] != -1} {
377            set Euler {}
378            foreach i [lrange [split $line "="] 1 3] {
379                lappend Euler [lindex $i 0]
380            }
381            #puts $line
382            #puts $Euler
383        }
384        if {[string first "Atom   Orthon" $line] != -1} {
385            set coordlist {}
386            for {set i 1} {$i <= $natom} {incr i} {
387                gets $fp line
388                set coord {}
389                lappend coord [string trim [string range $line 9 15]]
390                lappend coord [string trim [string range $line 16 22]]
391                lappend coord [string trim [string range $line 23 29]]
392                lappend coord [string trim [string range $line 0 8]]
393                #puts $line
394                #puts $coord
395                lappend coordlist $coord
396            }
397            #puts $coordlist
398        }
399    }
400    #file delete geom.out
401    return [list $origin $Euler $coordlist]
402}
403
404# updates the coordinates in a .EXP file after a rigid body instance has been added/changed
405proc RunRecalcRBCoords { } {
406    global expgui tcl_platform
407    set input [open resetmult.inp w]
408    puts $input "Y"
409    puts $input "l b"
410    puts $input "n"
411    puts $input "x x x"
412    puts $input "x"
413    close $input
414    # Save the current exp file
415    savearchiveexp
416    # disable the file changed monitor
417    set expgui(expModifiedLast) 0
418    set expnam [file root [file tail $expgui(expfile)]]
419    set err [catch {
420        if {$tcl_platform(platform) == "windows"} {
421            exec [file join $expgui(gsasexe) expedt.exe] $expnam < resetmult.inp >& resetmult.out
422        } else {
423            exec [file join $expgui(gsasexe) expedt] $expnam < resetmult.inp >& resetmult.out
424        }
425    } errmsg]
426    loadexp $expgui(expfile)
427    set fp [open resetmult.out r]
428    set out [read $fp]
429    close $fp
430    set expgui(exptoolout) $out
431    catch {file delete resetmult.inp resetmult.out}
432    if {$err} {
433        return $errmsg
434    } else {
435        return ""
436    }
437}
438
439
440# compute a rotation matrix for orthogonal coordinates (based on MAKMATD in GSAS)
441# rotate angle degrees around axis (1, 2 or 3) for (x, y, or z)
442# returns a list that can be used as a matrix in the La package
443proc RotMatrix {axis angle} {
444    set ang [expr {$angle * acos(0) / 90.}]
445    set mat "1 0 0 0 1 0 0 0 1"
446    if {$axis == 1}  {
447        set i1 1
448        set i2 2
449    } elseif {$axis == 2}  {
450        set i1 2
451        set i2 0
452    } else {
453        set i1 0
454        set i2 1
455    }
456    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
457    foreach item {
458        {$i1 $i1 [expr {cos($ang)}]}
459        {$i2 $i2 [expr {cos($ang)}]}
460        {$i1 $i2 [expr {-sin($ang)}]}
461        {$i2 $i1 [expr {sin($ang)}]}
462    } {
463        foreach {c r val} [subst $item] {}
464        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
465    }
466    return "2 3 3 $mat"
467}
468
469# compute the derivative of the rotation matrix with respect to the angle, see RotMatrix
470# (based on MAKMATD in GSAS)
471# returns a list that can be used as a matrix in the La package
472proc DerivRotMatrix {axis angle} {
473    set ang [expr {$angle * acos(0) / 90.}]
474    set mat "0 0 0 0 0 0 0 0 0"
475    if {$axis == 1}  {
476        set i1 1
477        set i2 2
478    } elseif {$axis == 2}  {
479        set i1 2
480        set i2 0
481    } else {
482        set i1 0
483        set i2 1
484    }
485    proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]}
486    foreach item {
487        {$i1 $i1 [expr {-sin($ang) * acos(0) / 90.}]}
488        {$i2 $i2 [expr {-sin($ang) * acos(0) / 90.}]}
489        {$i1 $i2 [expr {-cos($ang) * acos(0) / 90.}]}
490        {$i2 $i1 [expr {cos($ang) * acos(0) / 90.}]}
491    } {
492        foreach {c r val} [subst $item] {}
493        set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] 
494    }
495    return "2 3 3 $mat"
496}
497
498# compute an orthogonalization matrix from cell parameters (based on AMATRX in GSAS)
499# returns a list that can be used as a matrix in the La package
500proc OrthoMatrix {a b c alpha beta gamma} {
501    set CA [expr {cos($alpha * acos(0) / 90.)}]
502    set CB [expr {cos($beta * acos(0) / 90.)}]
503    set CG [expr {cos($gamma * acos(0) / 90.)}]
504    set SA [expr {sin($alpha * acos(0) / 90.)}]
505    set SB [expr {sin($beta * acos(0) / 90.)}]
506    set SC [expr {sin($gamma * acos(0) / 90.)}]
507    set CASTAR [expr { ($CB*$CG-$CA)/($SB*$SC) }]    ;#! cos(Alpha*)
508    set CBSTAR [expr { ($CA*$CG-$CB)/($SA*$SC) }]    ;#! cos(Beta*)
509    set CCSTAR [expr { ($CA*$CB-$CG)/($SA*$SB) }]    ;#! cos(Gamma*)
510    set SASTAR [expr { sqrt(1.0-($CASTAR*$CASTAR*2)) }]    ;#! sin(Alpha*)
511    set SBSTAR [expr { sqrt(1.0-($CBSTAR*$CBSTAR*2)) }]    ;#! sin(Beta*)
512    set SCSTAR [expr { sqrt(1.0-($CCSTAR*$CCSTAR*2)) }]    ;#! sin(Gamma*)
513
514    set A  "2 3 3      $a 0 0    0 $b 0    0 0 $c"
515    set A1 "2 3 3     1.0 0 0    $CG [expr {$SASTAR*$SC}] [expr {-$CASTAR*$SC}]    $CB 0.0 $SB"
516    #!This matrix is
517    #!   (1.0      0.0            0.0             )
518    #!   (cos(G)  sin(A*)*sin(G) -cos(A*)*sin(G)  )
519    #!   (cos(B)   0.0            sin(B)          )
520    return [La::transpose [La::mmult $A $A1]]
521}
522
523# compute the transformation matrix that converts a rigid body coordinates into fractional
524# coordinates
525# arguments:
526#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
527#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
528#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
529# returns a list that can be used as a matrix in the La package
530proc CalcXformMatrix {rotations cellprms} {
531    set prod {}
532    foreach item $rotations {
533        #puts $item
534        set mat [eval RotMatrix $item]
535        if {$prod == ""} {
536            set prod $mat
537        } else {
538            set prod [La::mmult $prod $mat]
539        }
540    }
541    #puts "--- rotation product ---"
542    #puts [La::show $prod]
543
544    set ortho [eval OrthoMatrix $cellprms]
545    #puts "--- ortho ---"
546    #puts [La::show $ortho]
547    set deortho [La::msolve $ortho [La::mident 3] ]
548    #puts "--- deortho ---"
549    #puts [La::show $deortho]
550    #puts "--- xform ---"
551    set xform [La::mmult $deortho $prod]
552    return $xform
553}
554
555# transforms a triplet of orthogonal coordinates into fractional ones using
556# arguments:
557#    xform: a transformation matrix from CalcXformMatrix
558#    origin: a list of fraction coordinates {x y z} describing the location of the
559#            origin of the orthogonal coordinates in the crystal system
560#    ortho: a triplet of othogonal coordinates
561# returns a triplet of fractional coordinates
562proc Ortho2Xtal {xform origin ortho} {
563    set ocv "2 3 0 $ortho"
564    set frac [La::mmult $xform $ocv]
565    #puts [La::show [La::transpose $frac]]
566    #puts $frac
567    set frac [La::madd $frac "[lrange $frac 0 2] $origin"]
568    #puts [La::show [La::transpose $frac]]
569    return $frac
570}
571
572# compute the derivative of the transformation matrix (see CalcXformMatrix)
573# with respect to every rotation in the $rotations list
574# arguments:
575#   rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
576#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
577#   cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees
578# returns a list of matrices where each matrix is a list that can be used as a
579# matrix in the La package
580proc CalcDerivMatrix {rotations cellprms} {
581    set ortho [eval OrthoMatrix $cellprms]
582    set deortho [La::msolve $ortho [La::mident 3] ]
583    set derivlist {}
584
585    foreach item $rotations {
586        #puts $item
587        set mat [eval DerivRotMatrix $item]
588        #puts $item
589        #puts [La::show $mat]
590        set xform [La::mmult $deortho $mat]
591        lappend derivlist $xform
592    }
593    return $derivlist
594}
595
596# fit a rigid body's origin
597# arguments:
598#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
599#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
600#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
601#  ortholist: list containing triplets with orthogonal coordinates
602#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
603#  fraclist: list containing triplets with fractional coordinates 
604#  origin: a list of fraction coordinates {x y z} describing the location of the
605#            origin of the orthogonal coordinates in the crystal system
606#     note that the length of ortholist, useflag and fraclist should be the same
607# Returns a list with the following elements
608#   0: a list with three new origin values
609#   1: the root-mean square difference between the fraclist coordinates and those computed
610#      using the input values for those atoms where $use is one (in Angstroms)
611#   2: the computed fractional coordinates for all atoms
612#   3: individual rms values for all atoms (in Angstroms)
613# note that items 1-3 are computed with the imput origin, not the revised one
614proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} {
615    set xform [CalcXformMatrix $Euler $cell]
616    foreach var {x y z} {set sum($var) 0.0}
617
618    set i 0
619    set sumdvs 0
620    set fracout {}
621    set rmsout {}
622    foreach oc $ortholist use $useflag coord $fraclist {
623        #puts "ortho: $oc"
624        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
625        lappend fracout $frac
626        if {$use} {incr i}
627        set dvs 0
628        foreach var {x y z} v1 $frac v2 $coord abc [lrange $cell 0 2] {
629            set dv [expr {($v2 - $v1)}]
630            set dvs [expr {$dvs + $dv*$dv}]
631            set sumdvs [expr {$sumdvs + $dv*$dv}]
632            if {$use} {set sum($var) [expr {$sum($var) + $dv}]}
633        }
634        lappend rmsout [expr {sqrt($dvs)}]
635    }
636    set rms 0
637    if {$i > 1} {set rms [expr {sqrt($sumdvs)/$i}]}
638    set neworig {}
639    foreach var {x y z} o $origin {
640        lappend neworig [expr {$o + ($sum($var)/$i)}]
641    }
642    return [list $neworig $rms $fracout $rmsout]
643}
644
645# fit a rigid body's Euler angles using least-squares
646# arguments:
647#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
648#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
649#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
650#  ortholist: list containing triplets with orthogonal coordinates
651#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
652#  fraclist: list containing triplets with fractional coordinates 
653#  origin: a list of fraction coordinates {x y z} describing the location of the
654#            origin of the orthogonal coordinates in the crystal system
655#     note that the length of ortholist, useflag and fraclist should be the same
656# Returns a list of new Euler angles
657proc FitBodyRot {Euler cell ortholist useflag fraclist origin} {
658    set xform [CalcXformMatrix $Euler  $cell]
659    set derivlist [CalcDerivMatrix $Euler  $cell]
660    set A "2 [expr 3*[llength $ortholist]] 3"
661    foreach oc $ortholist use $useflag coord $fraclist {
662        if {! $use} continue
663        foreach deriv $derivlist {
664            foreach xyz [lrange [Ortho2Xtal $deriv "0 0 0" $oc] 3 end] {
665                lappend A $xyz
666            }
667        }
668    }
669    #puts "A: [La::show $A]"
670    set y "2 [expr 3*[llength $ortholist]] 1"
671    foreach oc $ortholist use $useflag coord $fraclist {
672        if {! $use} continue
673        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
674        foreach xyz $coord XYZ $frac {
675            lappend y [expr {$XYZ - $xyz}]
676        }
677    }
678
679    set AtA [La::mmult [La::transpose $A] $A]
680    set Aty [La::mmult [La::transpose $A] $y]
681   
682    set l {}
683    #set shifts {}
684    foreach delta [lrange [La::msolve $AtA $Aty] 3 end] item $Euler {
685        #lappend shifts $delta
686        lappend l "[lindex $item 0] [expr {$delta + [lindex $item 1]}]"
687    }
688    #puts "shifts = $shifts"
689    return $l
690}
691
692# fit a rigid body's Origin and Euler angles
693# arguments:
694#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
695#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
696#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
697#  ortholist: list containing triplets with orthogonal coordinates
698#  useflag: list of flags to indicate if an atom should be used (1) or ignored (0)
699#  fraclist: list containing triplets with fractional coordinates 
700#  origin: a list of fraction coordinates {x y z} describing the location of the
701#            origin of the orthogonal coordinates in the crystal system
702#     note that the length of ortholist, useflag and fraclist should be the same
703# Returns a list containing
704#   new origin
705#   new Euler angles
706#   total rms
707#   fractional coordinates
708#   rms deviation in fractional coordinates of new Euler angles
709proc FitBody {Euler cell ortholist useflag fraclist origin "ncycle 5"} {
710    puts "start origin = $origin"
711    foreach {
712        origin 
713        startrms
714        fracout
715        rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
716    puts "start rms = $startrms"
717    set rmsprev $startrms
718    puts "new origin = $origin"
719    for {set i 0} {$i < $ncycle} {incr i} {
720        set Eulerprev $Euler
721        set Euler [FitBodyRot $Euler $cell $ortholist $useflag $fraclist $origin]
722        puts "New Euler $Euler"
723        puts "after fit" 
724        foreach {
725            origin 
726            rms
727            fracout
728            rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
729        if {$rms > (1.1 * $rmsprev) + 0.01} {
730            puts "rms = $rms, new origin = $origin"
731            set rmsprev $rms
732        }
733    }   
734    #proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} 
735    #return "$neworig $rms $fracout $rmsout"
736    set fmt  {"%8.5f %8.5f %8.5f     %8.5f %8.5f %8.5f   %6.3f"}
737    foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
738        puts "[eval format $fmt $fracin $fraccalc $rmsi]"
739    }
740    return [list $origin $Euler $rms $fracout $rmsout]
741}
742
743# zmat2coord converts a z-matrix to a set of cartesian coordinates
744#   a z-matrix is also known as "internal coordinates" or "torsion space"
745#   (see Journal of Computational Chemistry, Vol 26, #10, p. 1063–1068, 2005 or
746#    http://www.cmbi.ru.nl/molden/zmat/zmat.html)
747# INPUT:
748#   atmlist is a list of ascii lines where each line contains
749#     lbl c1 distance c2 angle c3 torsion
750#   where each atom is computed from the previous where the new atom is:
751#     distance $distance from atom $c1 (angstrom)
752#     angle $angle from $c1--$c2 (degrees)
753#     torsion $torsion from $c1--$c2--$c3 (degrees)
754# OUTPUT:
755#  zmat2coord returns a list of atom labels and cartesian coordinates,
756#  with 4 items in each element (label, x, y, z)
757# this routine was tested against results from Babel via the web interface at
758# http://www.shodor.org/chemviz/zmatrices/babel.html and sample input at
759# http://iopenshell.usc.edu/howto/zmatrix/
760proc zmat2coord {atmlist} { 
761    set torad [expr {acos(0)/90.}]
762    set i 0
763    foreach line $atmlist {
764        incr i
765        foreach {lbl c1 dist c2 angle c3 torsion} $line {} 
766        if {$i == 1} {
767            set atm(1) [list $lbl 0 0 0] ; # 1st atom is at origin
768        } elseif {$i == 2} {
769            set dist1 $dist
770            set atm(2) [list $lbl $dist1 0 0] ; # 2nd atom is along x-axis
771        } elseif {$i == 3} {
772            # 3rd atom can be bonded to the 1st or 2nd
773            if {$c1 == 1} {
774                set atm(3) [list $lbl \
775                                [expr {$dist * cos($torad * $angle)}] \
776                                [expr {$dist * sin($torad * $angle)}] \
777                                0]
778            } else {
779                set atm(3) [list $lbl \
780                                [expr {$dist1 - $dist * cos($torad * $angle)}] \
781                                [expr {$dist * sin($torad * $angle)}] \
782                                0]
783            }
784        } else {
785            set atm($i) [concat $lbl \
786                             [ahcat "atm" $c1 $dist $c2 $angle $c3 $torsion]]
787        }
788    }
789    set coordlist {}
790    foreach key [lsort -integer [array names atm]] {
791        lappend coordlist $atm($key)
792    }
793    return $coordlist
794}
795# Compute the length of a vector
796proc vlen {a} {
797    set sum 0.0
798    foreach ai $a {
799        set sum [expr {$sum + $ai*$ai}]
800    }
801    return [expr sqrt($sum)]
802}
803# compute vector (a + z * b) and optionally normalize to length d
804proc vadd {a b d z} {
805    set c {}
806    foreach ai $a bi $b {
807        lappend c [expr {$bi + $z * $ai}]
808    }
809    set v [vlen $c]
810    if {$d != 0} {
811        set r {}
812        foreach ci $c {
813            lappend r [expr {$d * $ci / $v}]
814        }
815        return [list $v $r]
816    }
817    return [list $v $c]
818}
819# normalize a vector
820proc vnrm {x} {
821    set v [vlen $x]
822    if {abs($v) < 1e-8} {return [list 0 0 0]}
823    set y {}
824    foreach xi $x {
825        lappend y [expr {$xi / $v}]
826    }
827    return $y
828}
829# compute the coordinates for an atom that is bonded:
830#   distance $dist from atom $nc
831#   angle $bondang from $nc--$nb
832#   torsion $torsang from $nc--$nb--$na
833#   coordinates are found in array $atmarr in the calling routine
834# based on a Fortran routine provided by Peter Zavalij (Thanks Peter!)
835proc ahcat {atmarr nc dist nb bondang na torsang} {
836    upvar 1 $atmarr atm
837    set xa [lrange $atm($na) 1 3]
838    set xb [lrange $atm($nb) 1 3]
839    set xc [lrange $atm($nc) 1 3]
840    set torad [expr {acos(0)/90.}]
841    # x = unit Vector A-B
842    foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. -1.] 1] {}
843    # y = unit Vector C-B
844    set y [lindex [vadd $xb $xc 1. -1.] 1]
845    foreach {y1 y2 y3} $y {}
846    set z1 [expr {$x2*$y3 - $x3*$y2}]
847    set z2 [expr {$x3*$y1 - $x1*$y3}]
848    set z3 [expr {$x1*$y2 - $x2*$y1}]
849    set z [vnrm [list $z1 $z2 $z3]]
850    set q1 [expr {$y2*$z3 - $y3*$z2}]
851    set q2 [expr {$y3*$z1 - $y1*$z3}]
852    set q3 [expr {$y1*$z2 - $y2*$z1}]
853    set q [vnrm [list $q1 $q2 $q3]]
854    set th [expr {$bondang * $torad}]
855    set ph [expr {-1. * $torsang * $torad}]
856    set cth [expr {cos($th)}]
857    set sth [expr {sin($th)}]
858    set cph [expr {cos($ph)}]
859    set sph [expr {sin($ph)}]
860    set xh {}
861    foreach xci $xc xi $q zi $z yi $y {
862        lappend xh [expr {
863                          $xci +
864                          $dist*($sth*($cph*$xi + $sph*$zi)-$cth*$yi)
865                      }]
866    }
867    return $xh
868}
869 
870
871
872#AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
873#puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"]
874#puts [GetRB 1 4 8 "1" "X 1 2" "Z 3 4"]
875#ApplyRigidBody 1 1 7 ".11 .22 .33" "11 12 13"
876
877
878#AddRigidBody {1} { {
879#    {1 1 1 o} {-1 1 1 o} {1 -1 1 o} {-1 -1 1 o}
880#    {1 1 -1 o} {-1 1 -1 o} {1 -1 -1 o} {-1 -1 -1 o}
881#} }
882#set n [ApplyRigidBody 1 1 1 ".2 .3 .4" "13 17 19"]
883#puts "body $n created"
884#incr expgui(changed)
885#RunRecalcRBCoords
886#puts "press Enter to continue"
887#gets stdin line
888#ApplyRigidBody 1 1 $n ".5 .5 .5" "0 0 0"
889#incr expgui(changed)
890#RunRecalcRBCoords
891
892puts "Test FitBody"
893set fraclist {
894    { 0.5483305238484277 0.4887545024531055 0.6167996784631056 }
895    { 0.1036801409356145 0.5954016321779562 0.5129448102437683 }
896    { 0.26404665760133855 0.09455414439078394 0.612655365147539 }
897    { -0.18060372531147473 0.20120127411563465 0.5088004969282018 }
898    { 0.5806037253114747 0.3987987258843653 0.2911995030717982 }
899    { 0.13595334239866147 0.5054458556092161 0.18734463485246095 }
900    { 0.2963198590643855 0.004598367822043814 0.2870551897562318 }
901    { -0.1483305238484277 0.1112454975468945 0.1832003215368945 }
902}
903set ortholist {
904    {1 1 1} 
905    {-1 1 1}
906    {      1.000000   -1.000000    1.000000}
907    {     -1.000000   -1.000000    1.000000}
908    {      1.000000    1.000000   -1.000000}
909    {     -1.000000    1.000000   -1.000000}
910    {      1.000000   -1.000000   -1.000000}
911    {     -1.000000   -1.000000   -1.000000} 
912}
913
914set useflag {1 1 1 1 1 1 1 1}
915set cell {4.  5695100105.}
916#set origin ".20 .30 .40"
917set origin ".0 .0 .0"
918#set Euler  {{1 13} {2 17} {3 19}}
919#set Euler  {{1 0} {2 180} {3 0}}
920set Euler  {{1 0} {2 0} {3 0}}
921
922#puts [La::show $xform]
923puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
924
925
926# test zmat2coord
927set atmlist {
928    {C1 0 0.0 0 0.0 0 0.0}
929    {O2 1 1.20 0 0.0 0 0.0}
930    {H3 1 1.10 2 120.0 0 0.0}
931    {C4 1 1.50 2 120.0 3 180.0}
932    {H5 4 1.10 1 110.0 2 0.00}
933    {H6 4 1.10 1 110.0 2 120.0}
934    {H7 4 1.10 1 110.0 2 -120.0}
935}
936#  C        0.00000        0.00000        0.00000
937#  O        1.20000        0.00000        0.00000
938#  H       -0.55000        0.95263        0.00000
939#  C       -0.75000       -1.29904       -0.00000
940#  H       -0.04293       -2.14169       -0.00000
941#  H       -1.38570       -1.36644        0.89518
942#  H       -1.38570       -1.36644       -0.89518
943set coordlist [zmat2coord $atmlist]
944set i 0
945puts "\nZmatrix in"
946foreach line $atmlist {
947    incr i
948    puts "$i) $line"
949}
950puts "Cartesian out"
951foreach line $coordlist {
952    puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
953}
Note: See TracBrowser for help on using the repository browser.