1 | #============================================================================ |
---|
2 | # rigid body EXP editing (to move into readexp.tcl) |
---|
3 | #============================================================================ |
---|
4 | # returns the number of defined rigid bodies |
---|
5 | proc 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. |
---|
24 | proc 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 | |
---|
56 | proc 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 | # |
---|
86 | proc 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 |
---|
175 | proc 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 |
---|
193 | proc 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) |
---|
208 | proc 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. |
---|
245 | proc ToHex {num digits} { |
---|
246 | return [string toupper [format "%${digits}x" $num]] |
---|
247 | } |
---|
248 | |
---|
249 | # convert a GSAS hex encoding to a decimal integer |
---|
250 | proc 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) |
---|
265 | proc 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 | # |
---|
312 | proc 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 | |
---|
325 | package 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 |
---|
344 | proc 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 |
---|
403 | proc 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 |
---|
441 | proc 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 |
---|
470 | proc 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 |
---|
498 | proc 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 |
---|
528 | proc 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 |
---|
560 | proc 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 |
---|
578 | proc 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 |
---|
612 | proc 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 |
---|
655 | proc 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 |
---|
707 | proc 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/ |
---|
758 | proc 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 |
---|
794 | proc 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 |
---|
802 | proc 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 |
---|
818 | proc 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!) |
---|
833 | proc 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 | |
---|
890 | puts "Test FitBody" |
---|
891 | set 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 | } |
---|
901 | set 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 | |
---|
912 | set useflag {1 1 1 1 1 1 1 1} |
---|
913 | set cell {4. 5. 6. 95. 100. 105.} |
---|
914 | #set origin ".20 .30 .40" |
---|
915 | set origin ".0 .0 .0" |
---|
916 | #set Euler {{1 13} {2 17} {3 19}} |
---|
917 | #set Euler {{1 0} {2 180} {3 0}} |
---|
918 | set Euler {{1 0} {2 0} {3 0}} |
---|
919 | |
---|
920 | #puts [La::show $xform] |
---|
921 | puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]" |
---|
922 | |
---|
923 | |
---|
924 | # test zmat2coord |
---|
925 | set 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 |
---|
941 | set coordlist [zmat2coord $atmlist] |
---|
942 | set i 0 |
---|
943 | puts "\nZmatrix in" |
---|
944 | foreach line $atmlist { |
---|
945 | incr i |
---|
946 | puts "$i) $line" |
---|
947 | } |
---|
948 | puts "Cartesian out" |
---|
949 | foreach line $coordlist { |
---|
950 | puts [eval format "%-4s%10.5f%10.5f%10.5f" $line] |
---|
951 | } |
---|