source: branches/sandbox/rb.tcl @ 1101

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

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