source: branches/sandbox/rb.tcl @ 1100

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

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