| 1 | # test code (package already loaded in expgui) |
|---|
| 2 | lappend auto_path [file dirname [info script]] |
|---|
| 3 | package require La |
|---|
| 4 | |
|---|
| 5 | #============================================================================ |
|---|
| 6 | # rigid body EXP editing routines (to move into readexp.tcl) |
|---|
| 7 | # RigidBodyList -- returns a list of the defined rigid body types |
|---|
| 8 | # ReadRigidBody -- # of times a body is mapped, scaling factors, var #s & coordinates |
|---|
| 9 | # RigidBodyMappingList - return a list instances where a RB is mapped in phase |
|---|
| 10 | # RigidBodyEnableTLS -- Enable or Disable TLS use for a rigid body mapping |
|---|
| 11 | # RigidBodySetTLS -- change the TLS values for a rigid body mapping |
|---|
| 12 | # RigidBodySetDamp -- change the damping values for a rigid body mapping |
|---|
| 13 | # RigidBodyVary -- set refinement variable numbers for a rigid body mapping |
|---|
| 14 | # RigidBodyTLSVary -- set TLS refinement variable nums for a rigid body mapping |
|---|
| 15 | # AddRigidBody -- defines a new rigid body type |
|---|
| 16 | # DeleteRigidBody -- remove a rigid body definition |
|---|
| 17 | # ReplaceRigidBody -- replaces a previous rigid body type |
|---|
| 18 | # ReadRigidBodyMapping -- get parameters for a rigid body mapping |
|---|
| 19 | # MapRigidBody -- map a rigid body type into a phase |
|---|
| 20 | # EditRigidBodyMapping -- change the parameters in a rigid body mapping |
|---|
| 21 | # UnMapRigidBody --remove a rigid body constraint by removing a RB "instance" |
|---|
| 22 | #----- note that these older routines should not be used ------ |
|---|
| 23 | # RigidBodyCount -- returns the number of defined rigid bodies (body types) |
|---|
| 24 | # use RigidBodyList instead |
|---|
| 25 | # RigidBodyMappingCount -- # of times a rigid body is mapped in phase |
|---|
| 26 | # use RigidBodyMappingList instead |
|---|
| 27 | #============================================================================ |
|---|
| 28 | # returns the number of defined rigid bodies |
|---|
| 29 | proc RigidBodyCount {} { |
|---|
| 30 | set n [string trim [readexp "RGBD NRBDS"]] |
|---|
| 31 | if {$n == ""} { |
|---|
| 32 | set n 0 |
|---|
| 33 | } |
|---|
| 34 | return $n |
|---|
| 35 | } |
|---|
| 36 | |
|---|
| 37 | # returns a list of the defined rigid body types |
|---|
| 38 | proc RigidBodyList {} { |
|---|
| 39 | set n [string trim [readexp "RGBD NRBDS"]] |
|---|
| 40 | if {$n == ""} { |
|---|
| 41 | set n 0 |
|---|
| 42 | } |
|---|
| 43 | set rblist {} |
|---|
| 44 | foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} { |
|---|
| 45 | set value $rbnum |
|---|
| 46 | validint value 2 |
|---|
| 47 | set key "RGBD${value}" |
|---|
| 48 | if {[existsexp "$key NATR "]} { |
|---|
| 49 | lappend rblist $rbnum |
|---|
| 50 | } |
|---|
| 51 | if {[llength $rblist] == $n} break |
|---|
| 52 | } |
|---|
| 53 | return $rblist |
|---|
| 54 | } |
|---|
| 55 | |
|---|
| 56 | # returns two items: |
|---|
| 57 | # the number of times the rigid body is mapped |
|---|
| 58 | # a list containing an element for each scaling factor in rigid body #rbnum. |
|---|
| 59 | # in each element there are four items: |
|---|
| 60 | # a multiplier value for the rigid body coordinates |
|---|
| 61 | # a damping value (0-9) for the refinement of the multiplier |
|---|
| 62 | # a variable number if the multiplier will be refined |
|---|
| 63 | # a list of cartesian coordinates coordinates |
|---|
| 64 | # each cartesian coordinate contains 4 items: x,y,z and a label |
|---|
| 65 | # note that the label is present only when the RB is created in EXPGUI and is |
|---|
| 66 | # not used in GSAS. |
|---|
| 67 | proc ReadRigidBody {rbnum} { |
|---|
| 68 | if {[lsearch [RigidBodyList] $rbnum] == -1} { |
|---|
| 69 | return "" |
|---|
| 70 | } |
|---|
| 71 | set value $rbnum |
|---|
| 72 | validint value 2 |
|---|
| 73 | set key "RGBD${value}" |
|---|
| 74 | set n [string trim [string range [readexp "$key NATR"] 0 4]] |
|---|
| 75 | set used [string trim [string range [readexp "$key NBDS"] 0 4]] |
|---|
| 76 | set nmult [string trim [string range [readexp "$key NSMP"] 0 4]] |
|---|
| 77 | set out {} |
|---|
| 78 | for {set i 1} {$i <= $nmult} {incr i} { |
|---|
| 79 | set line [readexp "${key}${i}PARM"] |
|---|
| 80 | set mult [string trim [string range $line 0 9]] |
|---|
| 81 | set var [string trim [string range $line 10 14]] |
|---|
| 82 | set damp [string trim [string range $line 15 19]] |
|---|
| 83 | set coordlist {} |
|---|
| 84 | for {set j 1} {$j <= $n} {incr j} { |
|---|
| 85 | set value $j |
|---|
| 86 | validint value 3 |
|---|
| 87 | set line [readexp "${key}${i}SC$value"] |
|---|
| 88 | set x [string trim [string range $line 0 9]] |
|---|
| 89 | set y [string trim [string range $line 10 19]] |
|---|
| 90 | set z [string trim [string range $line 20 29]] |
|---|
| 91 | set lbl [string trim [string range $line 30 39]] |
|---|
| 92 | lappend coordlist [list $x $y $z $lbl] |
|---|
| 93 | } |
|---|
| 94 | lappend out [list $mult $damp $var $coordlist] |
|---|
| 95 | } |
|---|
| 96 | return [list $used $out] |
|---|
| 97 | } |
|---|
| 98 | |
|---|
| 99 | # return the number of times rigid body $bodytyp is mapped in phase $phase |
|---|
| 100 | proc RigidBodyMappingCount {phase bodytyp} { |
|---|
| 101 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]" |
|---|
| 102 | if {! [existsexp "$key NBDS"]} {return 0} |
|---|
| 103 | set n [string trim [readexp "$key NBDS"]] |
|---|
| 104 | if {$n == ""} { |
|---|
| 105 | set n 0 |
|---|
| 106 | } |
|---|
| 107 | return $n |
|---|
| 108 | } |
|---|
| 109 | # return a list of the instances where rigid body $bodytyp is mapped in phase $phase |
|---|
| 110 | proc RigidBodyMappingList {phase bodytyp} { |
|---|
| 111 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]" |
|---|
| 112 | if {! [existsexp "$key NBDS"]} {return {}} |
|---|
| 113 | set n [string trim [readexp "$key NBDS"]] |
|---|
| 114 | if {$n == ""} { |
|---|
| 115 | set n 0 |
|---|
| 116 | } |
|---|
| 117 | set rblist {} |
|---|
| 118 | foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} { |
|---|
| 119 | set value $rbnum |
|---|
| 120 | validint value 2 |
|---|
| 121 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]" |
|---|
| 122 | if {[existsexp "$key NDA"]} { |
|---|
| 123 | lappend rblist $rbnum |
|---|
| 124 | } |
|---|
| 125 | if {[llength $rblist] == $n} break |
|---|
| 126 | } |
|---|
| 127 | return $rblist |
|---|
| 128 | } |
|---|
| 129 | |
|---|
| 130 | |
|---|
| 131 | |
|---|
| 132 | # reads rigid body mapping parameters for phase ($phase), body type # ($bodytyp) and instance # ($num) |
|---|
| 133 | # returns a list of items (most lists) as follows: |
|---|
| 134 | # 1) sequence # of first atom in body |
|---|
| 135 | # 2) origin of body in fractional coordinates (3 elements) |
|---|
| 136 | # 3) Euler angles as 6 pairs of numbers (see below) |
|---|
| 137 | # 4) variable numbers for the 9 position variables (origin followed by rotations) |
|---|
| 138 | # 5) damping vals for the 9 position variables (origin followed by rotations) |
|---|
| 139 | # 6) the TLS values, in order below (empty list if TLS is not in use) |
|---|
| 140 | # 7) the variable numbers for each TLS values, in order below (or empty) |
|---|
| 141 | # 8) three damping values for the T, L and S terms. |
|---|
| 142 | # returns an empty list if no such body exists. |
|---|
| 143 | # |
|---|
| 144 | # Euler angles are a list of axes and angles to rotate: |
|---|
| 145 | # { {axis1 angle1} {axis2 angle2} ...} |
|---|
| 146 | # where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes |
|---|
| 147 | # |
|---|
| 148 | # The 20 TLS terms are ordered: |
|---|
| 149 | # T11, T22, T33, T12, T13, T23 |
|---|
| 150 | # L11, L22, L33, L12, L13, L23 |
|---|
| 151 | # S12, S13, S21, S23, S31, S32, SAA, SBB |
|---|
| 152 | # |
|---|
| 153 | proc ReadRigidBodyMapping {phase bodytyp num} { |
|---|
| 154 | if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} { |
|---|
| 155 | return "" |
|---|
| 156 | } |
|---|
| 157 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]" |
|---|
| 158 | set first [string trim [string range [readexp "$key NDA"] 0 4]] |
|---|
| 159 | set line [readexp "$key BDFL"] |
|---|
| 160 | set varlist {} |
|---|
| 161 | set damplist {} |
|---|
| 162 | foreach i {0 1 2 3 4 5 6 7 8} { |
|---|
| 163 | lappend varlist [string trim [string range $line [expr {5*$i}] [expr {4 + 5*$i}] ]] |
|---|
| 164 | lappend damplist [string trim [string range $line [expr {45 + $i}] [expr {45 + $i}] ]] |
|---|
| 165 | } |
|---|
| 166 | set TLSdamplist {} |
|---|
| 167 | foreach i {54 55 56} { |
|---|
| 168 | lappend TLSdamplist [string trim [string range $line $i $i ]] |
|---|
| 169 | } |
|---|
| 170 | set line [readexp "${key} BDLC"] |
|---|
| 171 | set x [string trim [string range $line 0 9]] |
|---|
| 172 | set y [string trim [string range $line 10 19]] |
|---|
| 173 | set z [string trim [string range $line 20 29]] |
|---|
| 174 | set origin [list $x $y $z] |
|---|
| 175 | set line [readexp "${key} BDOR"] |
|---|
| 176 | set rotations {} |
|---|
| 177 | foreach i {0 10 20 30 40 50} { |
|---|
| 178 | set angle [string trim [string range $line $i [expr {$i+7}]]] |
|---|
| 179 | set axis [string trim [string range $line [expr {$i+8}] [expr {$i+9}]]] |
|---|
| 180 | lappend rotations [list $angle $axis] |
|---|
| 181 | } |
|---|
| 182 | set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]] |
|---|
| 183 | set tlsvars {} |
|---|
| 184 | set tlsvals {} |
|---|
| 185 | if {$TLS != 0} { |
|---|
| 186 | set line [readexp "${key}TLSF1"] |
|---|
| 187 | for {set j 0} {$j < 20} {incr j} { |
|---|
| 188 | lappend tlsvars [string trim [string range $line [expr {3*$j}] [expr {3*$j+2}]]] |
|---|
| 189 | } |
|---|
| 190 | for {set j 0} {$j < 20} {incr j} { |
|---|
| 191 | set i 0 |
|---|
| 192 | if {$j == 0} { |
|---|
| 193 | set i 1 |
|---|
| 194 | } elseif {$j == 8} { |
|---|
| 195 | set i 2 |
|---|
| 196 | } elseif {$j == 16} { |
|---|
| 197 | set i 3 |
|---|
| 198 | } |
|---|
| 199 | if {$i != 0} { |
|---|
| 200 | set line [readexp "${key}TLSP$i"] |
|---|
| 201 | set i 0 |
|---|
| 202 | set j1 0 |
|---|
| 203 | set j2 7 |
|---|
| 204 | } else { |
|---|
| 205 | incr j1 8 |
|---|
| 206 | incr j2 8 |
|---|
| 207 | } |
|---|
| 208 | lappend tlsvals [string trim [string range $line $j1 $j2]] |
|---|
| 209 | } |
|---|
| 210 | } |
|---|
| 211 | return [list $first $origin $rotations $varlist $damplist $tlsvals $tlsvars $TLSdamplist] |
|---|
| 212 | } |
|---|
| 213 | |
|---|
| 214 | # Control TLS representation for phase, body # and instance number of a Rigid body mapping |
|---|
| 215 | # for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num) |
|---|
| 216 | # Enable TLS use if TLS is non-zero (true). Disable if zero |
|---|
| 217 | proc RigidBodyEnableTLS {phase bodytyp num TLS} { |
|---|
| 218 | if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} { |
|---|
| 219 | return "" |
|---|
| 220 | } |
|---|
| 221 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]" |
|---|
| 222 | if {$TLS} { |
|---|
| 223 | setexp "${key} LSTF" [format "%5d" 1] 1 5 |
|---|
| 224 | if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"} |
|---|
| 225 | if {![existsexp "${key}TLSP1"]} { |
|---|
| 226 | makeexprec "${key}TLSP1" |
|---|
| 227 | set str {} |
|---|
| 228 | foreach v {.01 .01 .01 0 0 0 0 0} d {4 4 4 4 4 4 2 2} { |
|---|
| 229 | validreal v 8 $d |
|---|
| 230 | append str $v |
|---|
| 231 | } |
|---|
| 232 | setexp "${key}TLSP1" $str 1 64 |
|---|
| 233 | } |
|---|
| 234 | if {![existsexp "${key}TLSP2"]} { |
|---|
| 235 | makeexprec "${key}TLSP2" |
|---|
| 236 | set str {} |
|---|
| 237 | set v 0 |
|---|
| 238 | foreach d {2 2 2 2 4 4 4 4} { |
|---|
| 239 | validreal v 8 $d |
|---|
| 240 | append str $v |
|---|
| 241 | } |
|---|
| 242 | setexp "${key}TLSP2" $str 1 64 |
|---|
| 243 | } |
|---|
| 244 | if {![existsexp "${key}TLSP3"]} { |
|---|
| 245 | makeexprec "${key}TLSP3" |
|---|
| 246 | set str {} |
|---|
| 247 | set v 0 |
|---|
| 248 | foreach d {4 4 4 4} { |
|---|
| 249 | validreal v 8 $d |
|---|
| 250 | append str $v |
|---|
| 251 | } |
|---|
| 252 | setexp "${key}TLSP3" $str 1 64 |
|---|
| 253 | } |
|---|
| 254 | } else { |
|---|
| 255 | setexp "${key} LSTF" [format "%5d" 0] 1 5 |
|---|
| 256 | } |
|---|
| 257 | return 1 |
|---|
| 258 | } |
|---|
| 259 | |
|---|
| 260 | # Control the TLS values for Rigid body mapping for mapping with |
|---|
| 261 | # phase ($phase), body type # ($bodytyp) and instance # ($num) |
|---|
| 262 | # set the 20 TLS values to the values in TLSvals |
|---|
| 263 | # There must be exactly 20 TLS terms, which are ordered: |
|---|
| 264 | # T11, T22, T33, T12, T13, T23 |
|---|
| 265 | # L11, L22, L33, L12, L13, L23 |
|---|
| 266 | # S12, S13, S21, S23, S31, S32, SAA, SBB |
|---|
| 267 | proc RigidBodySetTLS {phase bodytyp num TLSvals} { |
|---|
| 268 | if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} { |
|---|
| 269 | return "" |
|---|
| 270 | } |
|---|
| 271 | if {[llength $TLSvals] != 20} {return ""} |
|---|
| 272 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]" |
|---|
| 273 | set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]] |
|---|
| 274 | if {$TLS == 0} {return ""} |
|---|
| 275 | if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"} |
|---|
| 276 | foreach n {1 2 3} { |
|---|
| 277 | if {![existsexp "${key}TLSP$n"]} {makeexprec "${key}TLSP$n"} |
|---|
| 278 | } |
|---|
| 279 | set str {} |
|---|
| 280 | set n 1 |
|---|
| 281 | set i 0 |
|---|
| 282 | foreach v $TLSvals d {4 4 4 4 4 4 2 2 2 2 2 2 4 4 4 4 4 4 4 4} { |
|---|
| 283 | incr i |
|---|
| 284 | validreal v 8 $d |
|---|
| 285 | append str $v |
|---|
| 286 | if {$i == 8} { |
|---|
| 287 | set i 0 |
|---|
| 288 | setexp "${key}TLSP$n" $str 1 64 |
|---|
| 289 | incr n |
|---|
| 290 | set str {} |
|---|
| 291 | } |
|---|
| 292 | } |
|---|
| 293 | setexp "${key}TLSP$n" $str 1 64 |
|---|
| 294 | return 1 |
|---|
| 295 | } |
|---|
| 296 | |
|---|
| 297 | # set damping values for a Rigid body mapping |
|---|
| 298 | # for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num) |
|---|
| 299 | # there must be 9 damping values in RBdamp for the 9 position variables (origin followed by rotations) |
|---|
| 300 | # Use of TLSdamp is optional, but to be used, TLS representation must be enabled and there must be |
|---|
| 301 | # three damping terms (for all T terms; for all L terms and for all S terms) |
|---|
| 302 | proc RigidBodySetDamp {phase bodytyp num RBdamp "TLSdamp {}"} { |
|---|
| 303 | if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} { |
|---|
| 304 | return "" |
|---|
| 305 | } |
|---|
| 306 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]" |
|---|
| 307 | if {[llength $RBdamp] != 9} {return ""} |
|---|
| 308 | set str {} |
|---|
| 309 | foreach v $RBdamp { |
|---|
| 310 | if {[validint v 1] != 1} {set v " "} |
|---|
| 311 | append str $v |
|---|
| 312 | } |
|---|
| 313 | setexp "$key BDFL" $str 46 9 |
|---|
| 314 | set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]] |
|---|
| 315 | if {$TLS != 0 && [llength $TLSdamp] == 3} { |
|---|
| 316 | set str {} |
|---|
| 317 | foreach v $TLSdamp { |
|---|
| 318 | if {[validint v 1] != 1} {set v " "} |
|---|
| 319 | append str $v |
|---|
| 320 | } |
|---|
| 321 | setexp "$key BDFL" $str 55 3 |
|---|
| 322 | } |
|---|
| 323 | return 1 |
|---|
| 324 | } |
|---|
| 325 | |
|---|
| 326 | # set refinement variable numbers for a Rigid body mapping |
|---|
| 327 | # for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num) |
|---|
| 328 | # there must be 9 variable values in RBvar for the 9 position variables (origin followed by rotations) |
|---|
| 329 | # note that the variable values should be unique integers |
|---|
| 330 | proc RigidBodyVary {phase bodytyp num RBvar} { |
|---|
| 331 | if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} { |
|---|
| 332 | return "" |
|---|
| 333 | } |
|---|
| 334 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]" |
|---|
| 335 | if {[llength $RBvar] != 9} {return ""} |
|---|
| 336 | set str {} |
|---|
| 337 | foreach v $RBvar { |
|---|
| 338 | if {[validint v 5] != 1} {set v " "} |
|---|
| 339 | append str $v |
|---|
| 340 | } |
|---|
| 341 | setexp "$key BDFL" $str 1 45 |
|---|
| 342 | } |
|---|
| 343 | |
|---|
| 344 | # set TLS refinement variable numbers for a Rigid body mapping |
|---|
| 345 | # for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num) |
|---|
| 346 | # there must be 20 variable values in TLSvar for the 20 parameters: |
|---|
| 347 | # T11, T22, T33, T12, T13, T23 |
|---|
| 348 | # L11, L22, L33, L12, L13, L23 |
|---|
| 349 | # S12, S13, S21, S23, S31, S32, SAA, SBB |
|---|
| 350 | # note that the variable values should be unique integers |
|---|
| 351 | proc RigidBodyTLSVary {phase bodytyp num TLSvar} { |
|---|
| 352 | if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} { |
|---|
| 353 | return "" |
|---|
| 354 | } |
|---|
| 355 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]" |
|---|
| 356 | if {[llength $TLSvar] != 20} {return ""} |
|---|
| 357 | set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]] |
|---|
| 358 | if {$TLS == 0} {return ""} |
|---|
| 359 | set str {} |
|---|
| 360 | foreach v $TLSvar { |
|---|
| 361 | if {[validint v 3] != 1} {set v " "} |
|---|
| 362 | append str $v |
|---|
| 363 | } |
|---|
| 364 | setexp "${key}TLSF1" $str 1 60 |
|---|
| 365 | |
|---|
| 366 | # AddRigidBody: add a new rigid body definition into the .EXP file |
|---|
| 367 | # arguments are: |
|---|
| 368 | # multlist: defines a list of multipliers for each set of coordinates. In the |
|---|
| 369 | # simplest case this will be {1} |
|---|
| 370 | # coordlist: a nested list of coordinates such as { { {0 0 0} {.1 .1 .1} {.2 .2 .2} } } |
|---|
| 371 | # note that when the length of multlist > 1 then coordlist must have the same length. |
|---|
| 372 | # for input where |
|---|
| 373 | # multlist = {s1 s2} and |
|---|
| 374 | # coordlist = { { {0 0 0} {1 1 0} {.0 .0 .0} ...} |
|---|
| 375 | # {0 0 0} {1 1 0} {2 1 2} ...} |
|---|
| 376 | # } |
|---|
| 377 | # the cartesian coordinates are defined from the input as |
|---|
| 378 | # atom 1 = s1 * (0,0,0) + s2*(0,0,0) [ = (0,0,0)] |
|---|
| 379 | # atom 2 = s1 * (1,1,0) + s2*(1,1,0) [ = (s1+s2) * (1,1,0)] |
|---|
| 380 | # atom 3 = s1 * (0,0,0) + s2*(2,1,2) [ = s2 * (2,1,2)] |
|---|
| 381 | # ... |
|---|
| 382 | # Returns the number of the rigid body that has been created |
|---|
| 383 | proc AddRigidBody {multlist coordlist} { |
|---|
| 384 | # find the first unused body # |
|---|
| 385 | foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} { |
|---|
| 386 | set value $rbnum |
|---|
| 387 | validint value 2 |
|---|
| 388 | set key "RGBD${value}" |
|---|
| 389 | if {! [existsexp "$key NATR "]} {break} |
|---|
| 390 | } |
|---|
| 391 | # did we go too far? |
|---|
| 392 | if {$rbnum == 16} {return ""} |
|---|
| 393 | # increment the RB counter |
|---|
| 394 | set n [string trim [readexp "RGBD NRBDS"]] |
|---|
| 395 | if {$n == ""} { |
|---|
| 396 | makeexprec "RGBD NRBDS" |
|---|
| 397 | set n 0 |
|---|
| 398 | } |
|---|
| 399 | incr n |
|---|
| 400 | validint n 5 |
|---|
| 401 | setexp "RGBD NRBDS" $n 1 5 |
|---|
| 402 | SetRigidBody $rbnum $multlist $coordlist |
|---|
| 403 | return $rbnum |
|---|
| 404 | } |
|---|
| 405 | |
|---|
| 406 | # DeleteRigidBody: remove a rigid body definition from the .EXP file |
|---|
| 407 | # The body may not be mapped. I am not sure if GSAS allows more than 9 bodies, |
|---|
| 408 | # but if it does, the simplifed approach used here will fail, so this |
|---|
| 409 | # is not allowed. |
|---|
| 410 | # Input: |
|---|
| 411 | # Rigid body number |
|---|
| 412 | # Returns: |
|---|
| 413 | # 1 on success |
|---|
| 414 | # -1 if the body number is 11 or greater |
|---|
| 415 | # -2 if the body is mapped |
|---|
| 416 | # -3 if the body is not defined |
|---|
| 417 | proc DeleteRigidBody {rbnum} { |
|---|
| 418 | # can't delete bodies with numbers higher than 10, since the key prefix |
|---|
| 419 | # (RGBD11... will overlap with rigid body instance records, which would be |
|---|
| 420 | # deleted below |
|---|
| 421 | if {$rbnum > 10} { |
|---|
| 422 | return -1 |
|---|
| 423 | } |
|---|
| 424 | set value $rbnum |
|---|
| 425 | validint value 2 |
|---|
| 426 | set key "RGBD${value}" |
|---|
| 427 | if {![existsexp "$key NATR "]} { |
|---|
| 428 | return -2 |
|---|
| 429 | } |
|---|
| 430 | # make sure the body is not mapped |
|---|
| 431 | if {[string trim [string range [readexp "$key NBDS"] 0 4]] != 0} { |
|---|
| 432 | return -3 |
|---|
| 433 | } |
|---|
| 434 | # delete the records starting with "RGBD x" or "RGBD10" |
|---|
| 435 | foreach key [array names ::exparray "${key}*"] { |
|---|
| 436 | #puts $key |
|---|
| 437 | delexp $key |
|---|
| 438 | } |
|---|
| 439 | # decrement the RB counter |
|---|
| 440 | set n [string trim [readexp "RGBD NRBDS"]] |
|---|
| 441 | if {$n == ""} { |
|---|
| 442 | set n 0 |
|---|
| 443 | } |
|---|
| 444 | incr n -1 |
|---|
| 445 | validint n 5 |
|---|
| 446 | if {$n > 0} { |
|---|
| 447 | setexp "RGBD NRBDS" $n 1 5 |
|---|
| 448 | } else { |
|---|
| 449 | delexp "RGBD NRBDS" |
|---|
| 450 | } |
|---|
| 451 | return 1 |
|---|
| 452 | } |
|---|
| 453 | |
|---|
| 454 | # ReplaceRigidBody: replace all the information for rigid body #rbnum |
|---|
| 455 | # Works the sames as AddRigidBody (see above) except that the rigid body is replaced rather |
|---|
| 456 | # than added. |
|---|
| 457 | # Note that count of the # of times the body is used is preserved |
|---|
| 458 | proc ReplaceRigidBody {rbnum multlist coordlist} { |
|---|
| 459 | set value $rbnum |
|---|
| 460 | validint value 2 |
|---|
| 461 | set key "RGBD${value}" |
|---|
| 462 | set line [readexp "$key NBDS"] |
|---|
| 463 | foreach key [array names ::exparray "${key}*"] { |
|---|
| 464 | #puts $key |
|---|
| 465 | delexp $key |
|---|
| 466 | } |
|---|
| 467 | SetRigidBody $rbnum $multlist $coordlist |
|---|
| 468 | setexp "$key NBDS" $line 1 68 |
|---|
| 469 | } |
|---|
| 470 | |
|---|
| 471 | # Edit the parameters for rigid body #rbnum |
|---|
| 472 | # (normally called from ReplaceRigidBody or AddRigidBody) |
|---|
| 473 | proc SetRigidBody {rbnum multlist coordlist} { |
|---|
| 474 | set value $rbnum |
|---|
| 475 | validint value 2 |
|---|
| 476 | set key "RGBD${value}" |
|---|
| 477 | # number of atoms |
|---|
| 478 | set value [llength [lindex $coordlist 0]] |
|---|
| 479 | validint value 5 |
|---|
| 480 | makeexprec "$key NATR" |
|---|
| 481 | setexp "$key NATR" $value 1 5 |
|---|
| 482 | # number of times used |
|---|
| 483 | set value 0 |
|---|
| 484 | validint value 5 |
|---|
| 485 | makeexprec "$key NBDS" |
|---|
| 486 | setexp "$key NBDS" $value 1 5 |
|---|
| 487 | # number of coordinate matrices |
|---|
| 488 | set value [llength $multlist] |
|---|
| 489 | validint value 5 |
|---|
| 490 | makeexprec "$key NSMP" |
|---|
| 491 | setexp "$key NSMP" $value 1 5 |
|---|
| 492 | set i 0 |
|---|
| 493 | foreach mult $multlist coords $coordlist { |
|---|
| 494 | incr i |
|---|
| 495 | makeexprec "${key}${i}PARM" |
|---|
| 496 | setexp "${key}${i}PARM" [format "%10.5f%5d%5d" $mult 0 0] 1 20 |
|---|
| 497 | set j 0 |
|---|
| 498 | foreach item $coords { |
|---|
| 499 | #puts $item |
|---|
| 500 | incr j |
|---|
| 501 | set value $j |
|---|
| 502 | validint value 3 |
|---|
| 503 | makeexprec "${key}${i}SC$value" |
|---|
| 504 | if {[llength $item] == 4} { |
|---|
| 505 | setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40 |
|---|
| 506 | } elseif {[llength $item] == 3} { |
|---|
| 507 | setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f" $item] 1 30 |
|---|
| 508 | } else { |
|---|
| 509 | return -code 3 "Invalid number of coordinates" |
|---|
| 510 | } |
|---|
| 511 | } |
|---|
| 512 | } |
|---|
| 513 | } |
|---|
| 514 | |
|---|
| 515 | # convert a decimal to the GSAS hex encoding with a field $digits long. |
|---|
| 516 | proc ToHex {num digits} { |
|---|
| 517 | return [string toupper [format "%${digits}x" $num]] |
|---|
| 518 | } |
|---|
| 519 | |
|---|
| 520 | # convert a GSAS hex encoding to a decimal integer |
|---|
| 521 | proc FromHex {hex} { |
|---|
| 522 | return [scan $hex "%x"] |
|---|
| 523 | } |
|---|
| 524 | |
|---|
| 525 | # MapRigidBody: define an "instance" of a rigid body: meaning that the coordinates |
|---|
| 526 | # (and optionally U values) for a set of atoms will be generated from the rigid body |
|---|
| 527 | # arguments: |
|---|
| 528 | # phase: phase number (1-9) |
|---|
| 529 | # bodytyp: number of rigid body (1-15) as returned from AddRigidBody |
|---|
| 530 | # firstatom: sequence number of the first atom in phase (note that atoms may |
|---|
| 531 | # not be numbered sequentially) |
|---|
| 532 | # position: list of three fractional coordinates for the origin of the rigid body coordinates |
|---|
| 533 | # angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the |
|---|
| 534 | # cartesian system before the body is translated to position. |
|---|
| 535 | # returns the instance # (number of times body $bodytyp has been used in phase $phase) |
|---|
| 536 | proc MapRigidBody {phase bodytyp firstatom position angles} { |
|---|
| 537 | # find the first unused body # for this phase & type |
|---|
| 538 | foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} { |
|---|
| 539 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]" |
|---|
| 540 | if {! [existsexp "$key NDA"]} {break} |
|---|
| 541 | } |
|---|
| 542 | # did we go too far? |
|---|
| 543 | if {$rbnum == 16} {return ""} |
|---|
| 544 | # increment number of mapped bodies of this type overall |
|---|
| 545 | set value $bodytyp |
|---|
| 546 | validint value 2 |
|---|
| 547 | set key "RGBD${value}" |
|---|
| 548 | set used [string trim [string range [readexp "$key NBDS"] 0 4]] |
|---|
| 549 | incr used |
|---|
| 550 | set value $used |
|---|
| 551 | validint value 5 |
|---|
| 552 | setexp "$key NBDS" $value 1 5 |
|---|
| 553 | # increment number of mapped bodies of this type in this phase |
|---|
| 554 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]" |
|---|
| 555 | if {[existsexp "$key NBDS"]} { |
|---|
| 556 | set used [string trim [string range [readexp "$key NBDS"] 0 4]] |
|---|
| 557 | } else { |
|---|
| 558 | makeexprec "$key NBDS" |
|---|
| 559 | set used 0 |
|---|
| 560 | } |
|---|
| 561 | incr used |
|---|
| 562 | set value $used |
|---|
| 563 | validint value 5 |
|---|
| 564 | setexp "$key NBDS" $value 1 5 |
|---|
| 565 | # now write the mapping parameters |
|---|
| 566 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]" |
|---|
| 567 | set value $firstatom |
|---|
| 568 | validint value 5 |
|---|
| 569 | makeexprec "$key NDA" |
|---|
| 570 | setexp "$key NDA" $value 1 5 |
|---|
| 571 | set l1 {} |
|---|
| 572 | set l2 {} |
|---|
| 573 | for {set i 0} {$i < 9} {incr i} { |
|---|
| 574 | append l1 [format %5d 0] |
|---|
| 575 | append l2 [format %1d 0] |
|---|
| 576 | } |
|---|
| 577 | makeexprec "$key BDFL" |
|---|
| 578 | setexp "$key BDFL" $l1$l2 1 54 |
|---|
| 579 | makeexprec "${key} BDLC" |
|---|
| 580 | setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30 |
|---|
| 581 | makeexprec "${key} BDOR" |
|---|
| 582 | set l1 {} |
|---|
| 583 | foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" { |
|---|
| 584 | append l1 [format "%8.2f%2d" $val $dir] |
|---|
| 585 | } |
|---|
| 586 | setexp "${key} BDOR" $l1 1 60 |
|---|
| 587 | makeexprec "${key} LSTF" |
|---|
| 588 | setexp "${key} LSTF" [format "%5d" 0] 1 5 |
|---|
| 589 | return $rbnum |
|---|
| 590 | } |
|---|
| 591 | |
|---|
| 592 | # EditRigidBodyMapping: edit parameters that define an "instance" of a rigid body (see MapRigidBody) |
|---|
| 593 | # arguments: |
|---|
| 594 | # phase: phase number (1-9) |
|---|
| 595 | # bodytyp: number of rigid body (1-15) as returned from AddRigidBody |
|---|
| 596 | # bodynum: instance number, as returned by MapRigidBody |
|---|
| 597 | # position: list of three fractional coordinates for the origin of the rigid body coordinates |
|---|
| 598 | # angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the |
|---|
| 599 | # cartesian system before the body is translated to position. |
|---|
| 600 | # |
|---|
| 601 | proc EditRigidBodyMapping {phase bodytyp bodynum position angles} { |
|---|
| 602 | # number of bodies of this type in this phase |
|---|
| 603 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $bodynum 1]" |
|---|
| 604 | setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30 |
|---|
| 605 | set l1 {} |
|---|
| 606 | foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" { |
|---|
| 607 | append l1 [format "%8.2f%2d" $val $dir] |
|---|
| 608 | } |
|---|
| 609 | setexp "${key} BDOR" $l1 1 60 |
|---|
| 610 | } |
|---|
| 611 | |
|---|
| 612 | # UnMapRigidBody: remove a rigid body constraint by removing a RB "instance" |
|---|
| 613 | # (undoes MapRigidBody) |
|---|
| 614 | # arguments: |
|---|
| 615 | # phase: phase number (1-9) |
|---|
| 616 | # bodytyp: number of rigid body (1-15) as returned from AddRigidBody |
|---|
| 617 | # bodynum: instance number, as returned by MapRigidBody |
|---|
| 618 | proc UnMapRigidBody {phase bodytyp mapnum} { |
|---|
| 619 | if {[lsearch [RigidBodyMappingList $phase $bodytyp] $mapnum] == -1} { |
|---|
| 620 | return "" |
|---|
| 621 | } |
|---|
| 622 | # decrement number of mapped bodies of this type overall |
|---|
| 623 | set value $bodytyp |
|---|
| 624 | validint value 2 |
|---|
| 625 | set key "RGBD${value}" |
|---|
| 626 | set used [string trim [string range [readexp "$key NBDS"] 0 4]] |
|---|
| 627 | incr used -1 |
|---|
| 628 | set value $used |
|---|
| 629 | validint value 5 |
|---|
| 630 | setexp "$key NBDS" $value 1 5 |
|---|
| 631 | # decrement number of mapped bodies of this type in this phase |
|---|
| 632 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]" |
|---|
| 633 | if {[existsexp "$key NBDS"]} { |
|---|
| 634 | set used [string trim [string range [readexp "$key NBDS"] 0 4]] |
|---|
| 635 | } else { |
|---|
| 636 | set used 0 |
|---|
| 637 | } |
|---|
| 638 | incr used -1 |
|---|
| 639 | set value $used |
|---|
| 640 | validint value 5 |
|---|
| 641 | if {$used > 0} { |
|---|
| 642 | setexp "$key NBDS" $value 1 5 |
|---|
| 643 | } else { |
|---|
| 644 | delexp "$key NBDS" |
|---|
| 645 | } |
|---|
| 646 | # now delete the mapping parameter records |
|---|
| 647 | set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $mapnum 1]" |
|---|
| 648 | foreach key [array names ::exparray "${key}*"] { |
|---|
| 649 | delexp $key |
|---|
| 650 | } |
|---|
| 651 | return $used |
|---|
| 652 | } |
|---|
| 653 | |
|---|
| 654 | #============================================================================ |
|---|
| 655 | # Rigid body utility routines |
|---|
| 656 | #============================================================================ |
|---|
| 657 | # RigidBodyGetVarNums: Returns a list of the variable numbers in use |
|---|
| 658 | # for rigid body variable parameters. |
|---|
| 659 | # RigidBodyAtomNums: returns a list of atom numbers that are mapped to |
|---|
| 660 | # rigid bodies in a selected phase |
|---|
| 661 | # RigidStartAtoms: returns a list of atoms that are allowed for creation of RB |
|---|
| 662 | # ExtractRigidBody: Use the GSAS geometry program to cartesian coordinates & |
|---|
| 663 | # setting info for a RB from fractional coordinates for atoms in a phase |
|---|
| 664 | # RunRecalcRBCoords: updates the coordinates in all phases after changes have |
|---|
| 665 | # been made to rigid parameters. |
|---|
| 666 | # CalcBody: Convert ortho to fractional coordinates using RB parameters |
|---|
| 667 | # FitBody: Optimize the origin and Euler angles to match a rigid body to a |
|---|
| 668 | # set of fractional coordinates |
|---|
| 669 | # zmat2coord: convert a z-matrix to a set of cartesian coordinates |
|---|
| 670 | # RB2cart: convert the representation used for rigid bodies into |
|---|
| 671 | # cartesian coordinates |
|---|
| 672 | # PlotRBtype: plot a rigid body with DRAWxtl |
|---|
| 673 | # PlotRBcoords: plot orthogonal coordinates with DRAWxtl |
|---|
| 674 | # DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed |
|---|
| 675 | # on a structure read from a phase with DRAWxtl |
|---|
| 676 | #============================================================================ |
|---|
| 677 | #============================================================================ |
|---|
| 678 | # RigidBodyGetVarNums: Returns a list of the variable numbers used already |
|---|
| 679 | # for rigid body variable parameters |
|---|
| 680 | proc RigidBodyGetVarNums {} { |
|---|
| 681 | set varlist {} |
|---|
| 682 | foreach type [RigidBodyList] { |
|---|
| 683 | foreach phase $::expmap(phaselist) { |
|---|
| 684 | foreach i [RigidBodyMappingList $phase $type] { |
|---|
| 685 | set items [ReadRigidBodyMapping $phase $type $i] |
|---|
| 686 | set varlist [concat $varlist [lindex $items 3]] |
|---|
| 687 | if {[llength [lindex $items 6]] > 0} { |
|---|
| 688 | set varlist [concat $varlist [lindex $items 6]] |
|---|
| 689 | } |
|---|
| 690 | } |
|---|
| 691 | } |
|---|
| 692 | } |
|---|
| 693 | return [lsort -integer -unique $varlist] |
|---|
| 694 | } |
|---|
| 695 | |
|---|
| 696 | # RigidBodyAtomNums: Returns a list of the atoms mapped to rigid bodies in |
|---|
| 697 | # phase $phase |
|---|
| 698 | proc RigidBodyAtomNums {phase} { |
|---|
| 699 | if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""} |
|---|
| 700 | set allatoms $::expmap(atomlist_$phase) |
|---|
| 701 | # get matching atoms coordinate range |
|---|
| 702 | set mappedlist {} |
|---|
| 703 | foreach type [RigidBodyList] { |
|---|
| 704 | foreach i [RigidBodyMappingList $phase $type] { |
|---|
| 705 | # get the number of atoms in this type of body |
|---|
| 706 | set natoms [llength [lindex [lindex [lindex [ReadRigidBody $type] 1] 0] 3]] |
|---|
| 707 | set natom1 [expr {$natoms - 1}] |
|---|
| 708 | set items [ReadRigidBodyMapping $phase $type $i] |
|---|
| 709 | set firstatom [lindex $items 0] |
|---|
| 710 | set firstind [lsearch $allatoms $firstatom] |
|---|
| 711 | set mappedlist [concat $mappedlist \ |
|---|
| 712 | [lrange \ |
|---|
| 713 | [lrange $allatoms $firstind end] \ |
|---|
| 714 | 0 $natom1] \ |
|---|
| 715 | ] |
|---|
| 716 | } |
|---|
| 717 | } |
|---|
| 718 | return [lsort -integer $mappedlist] |
|---|
| 719 | } |
|---|
| 720 | |
|---|
| 721 | # RigidStartAtoms: Find allowed starting atoms for a rigid body in a phase |
|---|
| 722 | # Input: |
|---|
| 723 | # phase is the phase number |
|---|
| 724 | # natoms is the number of atoms in the RB to be mapped |
|---|
| 725 | # Returns a list of valid "start" atoms. |
|---|
| 726 | # Example: if the atom numbers in the phase are {2 4 5 6 7 8} and no rigid bodies |
|---|
| 727 | # are mapped, then a 4-atom body can be mapped starting with atom 2, 4 or 5 only, |
|---|
| 728 | # so {2 4 5} is returned |
|---|
| 729 | # If atoms 2-6 were already mapped, then this routine would return an empty |
|---|
| 730 | # list, as it is not possible to map the body. |
|---|
| 731 | proc RigidStartAtoms {phase natoms} { |
|---|
| 732 | if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""} |
|---|
| 733 | set allatoms $::expmap(atomlist_$phase) |
|---|
| 734 | set usedatoms [RigidBodyAtomNums $phase] |
|---|
| 735 | set startatomlist {} |
|---|
| 736 | for {set i 0} {$i < [llength $allatoms]} {incr i} { |
|---|
| 737 | set al [lrange $allatoms $i [expr {$i+$natoms-1}]] |
|---|
| 738 | if {[llength $al] < $natoms} break |
|---|
| 739 | set ok 1 |
|---|
| 740 | foreach atom $al { |
|---|
| 741 | if {[lsearch $usedatoms $atom] != -1} { |
|---|
| 742 | set ok 0 |
|---|
| 743 | break |
|---|
| 744 | } |
|---|
| 745 | } |
|---|
| 746 | if $ok {lappend startatomlist [lindex $al 0]} |
|---|
| 747 | } |
|---|
| 748 | return $startatomlist |
|---|
| 749 | } |
|---|
| 750 | |
|---|
| 751 | # ExtractRigidBody: Use the GSAS geometry program to compute a set of cartesian coordinates for a |
|---|
| 752 | # set of atoms in a .EXP file and provide the origin shift and Euler angles needed to |
|---|
| 753 | # place the cartesian system into the crystal coordinates. Used for setting up a rigid body. |
|---|
| 754 | # Returns a nested list of lists: |
|---|
| 755 | # element 0: a list of the origin location {x y z} in fraction coordinates |
|---|
| 756 | # element 1: a list of three rotation angles in form {a1 a2 a3} |
|---|
| 757 | # where a1, a2 and a3 are rotations around the cartesian x, y and z axes |
|---|
| 758 | # element 2: a list of $natom cartesian coordinate triples {{x1 y1 z1} {x2 y2 z2}...} |
|---|
| 759 | # arguments: |
|---|
| 760 | # phase: phase # |
|---|
| 761 | # natom: number of atoms in group |
|---|
| 762 | # firstatom: sequence # in phase (may be > than number of the atom) |
|---|
| 763 | # originlist: atoms to define origin (where 1 is first atom in group; <= natom) |
|---|
| 764 | # vector1: list of 3 values with X, Y or Z, atom #a and #b (number as in origin) |
|---|
| 765 | # vector2: list of 3 values with X, Y or Z, atom #a and #b (number as in origin) |
|---|
| 766 | # note that vector2 must define a different axis than vector1 |
|---|
| 767 | # also and vector1 and vector2 cannot use the same atom pair |
|---|
| 768 | proc ExtractRigidBody {phase natom firstatom originlist vector1 vector2} { |
|---|
| 769 | global expgui |
|---|
| 770 | set fp [open "geom.inp" "w"] |
|---|
| 771 | puts $fp "N" |
|---|
| 772 | puts $fp "N" |
|---|
| 773 | puts $fp $phase |
|---|
| 774 | puts $fp "N" |
|---|
| 775 | |
|---|
| 776 | puts $fp "R" |
|---|
| 777 | puts $fp "$natom" |
|---|
| 778 | puts $fp "$firstatom" |
|---|
| 779 | puts $fp [llength $originlist] |
|---|
| 780 | foreach i $originlist { |
|---|
| 781 | puts $fp $i |
|---|
| 782 | } |
|---|
| 783 | foreach i [concat $vector1 $vector2] { |
|---|
| 784 | puts $fp $i |
|---|
| 785 | } |
|---|
| 786 | puts $fp "0" |
|---|
| 787 | puts $fp "X" |
|---|
| 788 | close $fp |
|---|
| 789 | catch { |
|---|
| 790 | exec [file join $expgui(gsasexe) geometry] $expgui(expfile) < geom.inp > geom.out |
|---|
| 791 | } |
|---|
| 792 | file delete geom.inp |
|---|
| 793 | set fp [open geom.out r] |
|---|
| 794 | while {[gets $fp line] >= 0} { |
|---|
| 795 | if {[string first "Cell coordinates of origin" $line] != -1} { |
|---|
| 796 | set origin [lrange [string range $line [string first "are" $line] end] 1 3] |
|---|
| 797 | } |
|---|
| 798 | if {[string first "Rotation angles" $line] != -1} { |
|---|
| 799 | set Euler {} |
|---|
| 800 | foreach i [lrange [split $line "="] 1 3] { |
|---|
| 801 | lappend Euler [lindex $i 0] |
|---|
| 802 | } |
|---|
| 803 | #puts $line |
|---|
| 804 | #puts $Euler |
|---|
| 805 | } |
|---|
| 806 | if {[string first "Atom Orthon" $line] != -1} { |
|---|
| 807 | set coordlist {} |
|---|
| 808 | for {set i 1} {$i <= $natom} {incr i} { |
|---|
| 809 | gets $fp line |
|---|
| 810 | set coord {} |
|---|
| 811 | lappend coord [string trim [string range $line 9 15]] |
|---|
| 812 | lappend coord [string trim [string range $line 16 22]] |
|---|
| 813 | lappend coord [string trim [string range $line 23 29]] |
|---|
| 814 | lappend coord [string trim [string range $line 0 8]] |
|---|
| 815 | #puts $line |
|---|
| 816 | #puts $coord |
|---|
| 817 | lappend coordlist $coord |
|---|
| 818 | } |
|---|
| 819 | #puts $coordlist |
|---|
| 820 | } |
|---|
| 821 | } |
|---|
| 822 | #file delete geom.out |
|---|
| 823 | return [list $origin $Euler $coordlist] |
|---|
| 824 | } |
|---|
| 825 | |
|---|
| 826 | # RunRecalcRBCoords: updates the coordinates in a .EXP file after a rigid |
|---|
| 827 | # body has been changed, mapped or the setting info is changed |
|---|
| 828 | proc RunRecalcRBCoords { } { |
|---|
| 829 | global expgui tcl_platform |
|---|
| 830 | set input [open resetmult.inp w] |
|---|
| 831 | puts $input "Y" |
|---|
| 832 | puts $input "l b" |
|---|
| 833 | puts $input "n" |
|---|
| 834 | puts $input "x x x" |
|---|
| 835 | puts $input "x" |
|---|
| 836 | close $input |
|---|
| 837 | # Save the current exp file |
|---|
| 838 | savearchiveexp |
|---|
| 839 | # disable the file changed monitor |
|---|
| 840 | set expgui(expModifiedLast) 0 |
|---|
| 841 | set expnam [file root [file tail $expgui(expfile)]] |
|---|
| 842 | set err [catch { |
|---|
| 843 | if {$tcl_platform(platform) == "windows"} { |
|---|
| 844 | exec [file join $expgui(gsasexe) expedt.exe] $expnam < resetmult.inp >& resetmult.out |
|---|
| 845 | } else { |
|---|
| 846 | exec [file join $expgui(gsasexe) expedt] $expnam < resetmult.inp >& resetmult.out |
|---|
| 847 | } |
|---|
| 848 | } errmsg] |
|---|
| 849 | loadexp $expgui(expfile) |
|---|
| 850 | set fp [open resetmult.out r] |
|---|
| 851 | set out [read $fp] |
|---|
| 852 | close $fp |
|---|
| 853 | set expgui(exptoolout) $out |
|---|
| 854 | catch {file delete resetmult.inp resetmult.out} |
|---|
| 855 | if {$err} { |
|---|
| 856 | return $errmsg |
|---|
| 857 | } else { |
|---|
| 858 | return "" |
|---|
| 859 | } |
|---|
| 860 | } |
|---|
| 861 | |
|---|
| 862 | |
|---|
| 863 | # compute a rotation matrix for orthogonal coordinates (based on MAKMATD in GSAS) |
|---|
| 864 | # rotate angle degrees around axis (1, 2 or 3) for (x, y, or z) |
|---|
| 865 | # returns a list that can be used as a matrix in the La package |
|---|
| 866 | proc RotMatrix {axis angle} { |
|---|
| 867 | set ang [expr {$angle * acos(0) / 90.}] |
|---|
| 868 | set mat "1 0 0 0 1 0 0 0 1" |
|---|
| 869 | if {$axis == 1} { |
|---|
| 870 | set i1 1 |
|---|
| 871 | set i2 2 |
|---|
| 872 | } elseif {$axis == 2} { |
|---|
| 873 | set i1 2 |
|---|
| 874 | set i2 0 |
|---|
| 875 | } else { |
|---|
| 876 | set i1 0 |
|---|
| 877 | set i2 1 |
|---|
| 878 | } |
|---|
| 879 | proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]} |
|---|
| 880 | foreach item { |
|---|
| 881 | {$i1 $i1 [expr {cos($ang)}]} |
|---|
| 882 | {$i2 $i2 [expr {cos($ang)}]} |
|---|
| 883 | {$i1 $i2 [expr {-sin($ang)}]} |
|---|
| 884 | {$i2 $i1 [expr {sin($ang)}]} |
|---|
| 885 | } { |
|---|
| 886 | foreach {c r val} [subst $item] {} |
|---|
| 887 | set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] |
|---|
| 888 | } |
|---|
| 889 | return "2 3 3 $mat" |
|---|
| 890 | } |
|---|
| 891 | |
|---|
| 892 | # compute the derivative of the rotation matrix with respect to the angle, see RotMatrix |
|---|
| 893 | # (based on MAKMATD in GSAS) |
|---|
| 894 | # returns a list that can be used as a matrix in the La package |
|---|
| 895 | proc DerivRotMatrix {axis angle} { |
|---|
| 896 | set ang [expr {$angle * acos(0) / 90.}] |
|---|
| 897 | set mat "0 0 0 0 0 0 0 0 0" |
|---|
| 898 | if {$axis == 1} { |
|---|
| 899 | set i1 1 |
|---|
| 900 | set i2 2 |
|---|
| 901 | } elseif {$axis == 2} { |
|---|
| 902 | set i1 2 |
|---|
| 903 | set i2 0 |
|---|
| 904 | } else { |
|---|
| 905 | set i1 0 |
|---|
| 906 | set i2 1 |
|---|
| 907 | } |
|---|
| 908 | proc imat {i1 i2} {return [expr {(3*$i2) + $i1}]} |
|---|
| 909 | foreach item { |
|---|
| 910 | {$i1 $i1 [expr {-sin($ang) * acos(0) / 90.}]} |
|---|
| 911 | {$i2 $i2 [expr {-sin($ang) * acos(0) / 90.}]} |
|---|
| 912 | {$i1 $i2 [expr {-cos($ang) * acos(0) / 90.}]} |
|---|
| 913 | {$i2 $i1 [expr {cos($ang) * acos(0) / 90.}]} |
|---|
| 914 | } { |
|---|
| 915 | foreach {c r val} [subst $item] {} |
|---|
| 916 | set mat [lreplace $mat [imat $c $r] [imat $c $r] $val] |
|---|
| 917 | } |
|---|
| 918 | return "2 3 3 $mat" |
|---|
| 919 | } |
|---|
| 920 | |
|---|
| 921 | # compute an orthogonalization matrix from cell parameters (based on AMATRX in GSAS) |
|---|
| 922 | # returns a list that can be used as a matrix in the La package |
|---|
| 923 | proc OrthoMatrix {a b c alpha beta gamma} { |
|---|
| 924 | set CA [expr {cos($alpha * acos(0) / 90.)}] |
|---|
| 925 | set CB [expr {cos($beta * acos(0) / 90.)}] |
|---|
| 926 | set CG [expr {cos($gamma * acos(0) / 90.)}] |
|---|
| 927 | set SA [expr {sin($alpha * acos(0) / 90.)}] |
|---|
| 928 | set SB [expr {sin($beta * acos(0) / 90.)}] |
|---|
| 929 | set SC [expr {sin($gamma * acos(0) / 90.)}] |
|---|
| 930 | set CASTAR [expr { ($CB*$CG-$CA)/($SB*$SC) }] ;#! cos(Alpha*) |
|---|
| 931 | set CBSTAR [expr { ($CA*$CG-$CB)/($SA*$SC) }] ;#! cos(Beta*) |
|---|
| 932 | set CCSTAR [expr { ($CA*$CB-$CG)/($SA*$SB) }] ;#! cos(Gamma*) |
|---|
| 933 | set SASTAR [expr { sqrt(1.0-($CASTAR*$CASTAR*2)) }] ;#! sin(Alpha*) |
|---|
| 934 | set SBSTAR [expr { sqrt(1.0-($CBSTAR*$CBSTAR*2)) }] ;#! sin(Beta*) |
|---|
| 935 | set SCSTAR [expr { sqrt(1.0-($CCSTAR*$CCSTAR*2)) }] ;#! sin(Gamma*) |
|---|
| 936 | |
|---|
| 937 | set A "2 3 3 $a 0 0 0 $b 0 0 0 $c" |
|---|
| 938 | set A1 "2 3 3 1.0 0 0 $CG [expr {$SASTAR*$SC}] [expr {-$CASTAR*$SC}] $CB 0.0 $SB" |
|---|
| 939 | #!This matrix is |
|---|
| 940 | #! (1.0 0.0 0.0 ) |
|---|
| 941 | #! (cos(G) sin(A*)*sin(G) -cos(A*)*sin(G) ) |
|---|
| 942 | #! (cos(B) 0.0 sin(B) ) |
|---|
| 943 | return [La::transpose [La::mmult $A $A1]] |
|---|
| 944 | } |
|---|
| 945 | |
|---|
| 946 | # compute the transformation matrix that converts a rigid body coordinates into fractional |
|---|
| 947 | # coordinates |
|---|
| 948 | # arguments: |
|---|
| 949 | # rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...} |
|---|
| 950 | # where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes |
|---|
| 951 | # cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees |
|---|
| 952 | # returns a list that can be used as a matrix in the La package |
|---|
| 953 | proc CalcXformMatrix {rotations cellprms} { |
|---|
| 954 | set prod {} |
|---|
| 955 | foreach item $rotations { |
|---|
| 956 | #puts $item |
|---|
| 957 | set mat [eval RotMatrix $item] |
|---|
| 958 | if {$prod == ""} { |
|---|
| 959 | set prod $mat |
|---|
| 960 | } else { |
|---|
| 961 | set prod [La::mmult $prod $mat] |
|---|
| 962 | } |
|---|
| 963 | } |
|---|
| 964 | #puts "--- rotation product ---" |
|---|
| 965 | #puts [La::show $prod] |
|---|
| 966 | |
|---|
| 967 | set ortho [eval OrthoMatrix $cellprms] |
|---|
| 968 | #puts "--- ortho ---" |
|---|
| 969 | #puts [La::show $ortho] |
|---|
| 970 | set deortho [La::msolve $ortho [La::mident 3] ] |
|---|
| 971 | #puts "--- deortho ---" |
|---|
| 972 | #puts [La::show $deortho] |
|---|
| 973 | #puts "--- xform ---" |
|---|
| 974 | set xform [La::mmult $deortho $prod] |
|---|
| 975 | return $xform |
|---|
| 976 | } |
|---|
| 977 | |
|---|
| 978 | # transforms a triplet of orthogonal coordinates into fractional ones using |
|---|
| 979 | # arguments: |
|---|
| 980 | # xform: a transformation matrix from CalcXformMatrix |
|---|
| 981 | # origin: a list of fraction coordinates {x y z} describing the location of the |
|---|
| 982 | # origin of the orthogonal coordinates in the crystal system |
|---|
| 983 | # ortho: a triplet of othogonal coordinates |
|---|
| 984 | # returns a triplet of fractional coordinates |
|---|
| 985 | proc Ortho2Xtal {xform origin ortho} { |
|---|
| 986 | set ocv "2 3 0 $ortho" |
|---|
| 987 | set frac [La::mmult $xform $ocv] |
|---|
| 988 | #puts [La::show [La::transpose $frac]] |
|---|
| 989 | #puts $frac |
|---|
| 990 | set frac [La::madd $frac "[lrange $frac 0 2] $origin"] |
|---|
| 991 | #puts [La::show [La::transpose $frac]] |
|---|
| 992 | return $frac |
|---|
| 993 | } |
|---|
| 994 | |
|---|
| 995 | # compute the derivative of the transformation matrix (see CalcXformMatrix) |
|---|
| 996 | # with respect to every rotation in the $rotations list |
|---|
| 997 | # arguments: |
|---|
| 998 | # rotations: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...} |
|---|
| 999 | # where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes |
|---|
| 1000 | # cellprms: a list with "a b c alpha beta gamma" in Angstroms and degrees |
|---|
| 1001 | # returns a list of matrices where each matrix is a list that can be used as a |
|---|
| 1002 | # matrix in the La package |
|---|
| 1003 | proc CalcDerivMatrix {rotations cellprms} { |
|---|
| 1004 | set ortho [eval OrthoMatrix $cellprms] |
|---|
| 1005 | set deortho [La::msolve $ortho [La::mident 3] ] |
|---|
| 1006 | set derivlist {} |
|---|
| 1007 | |
|---|
| 1008 | foreach item $rotations { |
|---|
| 1009 | #puts $item |
|---|
| 1010 | set mat [eval DerivRotMatrix $item] |
|---|
| 1011 | #puts $item |
|---|
| 1012 | #puts [La::show $mat] |
|---|
| 1013 | set xform [La::mmult $deortho $mat] |
|---|
| 1014 | lappend derivlist $xform |
|---|
| 1015 | } |
|---|
| 1016 | return $derivlist |
|---|
| 1017 | } |
|---|
| 1018 | |
|---|
| 1019 | # CalcBody: Calculate fractional coordinates using rigid body setting parameters |
|---|
| 1020 | # arguments: |
|---|
| 1021 | # Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...} |
|---|
| 1022 | # where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes |
|---|
| 1023 | # cell: a list with "a b c alpha beta gamma" in Angstroms and degrees |
|---|
| 1024 | # ortholist: list containing triplets with orthogonal coordinates |
|---|
| 1025 | # origin: a list of fraction coordinates {x y z} describing the location of the |
|---|
| 1026 | # origin of the orthogonal coordinates in the crystal system |
|---|
| 1027 | # note that the length of ortholist, useflag and fraclist should be the same |
|---|
| 1028 | # Returns a list with the computed fractional coordinates for all atoms |
|---|
| 1029 | proc CalcBody {Euler cell ortholist origin} { |
|---|
| 1030 | set xform [CalcXformMatrix $Euler $cell] |
|---|
| 1031 | set i 0 |
|---|
| 1032 | set sumdvs 0 |
|---|
| 1033 | set fracout {} |
|---|
| 1034 | set rmsout {} |
|---|
| 1035 | foreach oc $ortholist { |
|---|
| 1036 | set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end] |
|---|
| 1037 | lappend fracout $frac |
|---|
| 1038 | } |
|---|
| 1039 | return $fracout |
|---|
| 1040 | } |
|---|
| 1041 | |
|---|
| 1042 | |
|---|
| 1043 | # fit a rigid body's origin |
|---|
| 1044 | # arguments: |
|---|
| 1045 | # Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...} |
|---|
| 1046 | # where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes |
|---|
| 1047 | # cell: a list with "a b c alpha beta gamma" in Angstroms and degrees |
|---|
| 1048 | # ortholist: list containing triplets with orthogonal coordinates |
|---|
| 1049 | # useflag: list of flags to indicate if an atom should be used (1) or ignored (0) |
|---|
| 1050 | # fraclist: list containing triplets with fractional coordinates |
|---|
| 1051 | # origin: a list of fraction coordinates {x y z} describing the location of the |
|---|
| 1052 | # origin of the orthogonal coordinates in the crystal system |
|---|
| 1053 | # note that the length of ortholist, useflag and fraclist should be the same |
|---|
| 1054 | # Returns a list with the following elements |
|---|
| 1055 | # 0: a list with three new origin values |
|---|
| 1056 | # 1: the root-mean square difference between the fraclist coordinates and those computed |
|---|
| 1057 | # using the input values for those atoms where $use is one (in Angstroms) |
|---|
| 1058 | # 2: the computed fractional coordinates for all atoms |
|---|
| 1059 | # 3: individual rms values for all atoms (in Angstroms) |
|---|
| 1060 | # note that items 1-3 are computed with the imput origin, not the revised one |
|---|
| 1061 | proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} { |
|---|
| 1062 | set xform [CalcXformMatrix $Euler $cell] |
|---|
| 1063 | foreach var {x y z} {set sum($var) 0.0} |
|---|
| 1064 | |
|---|
| 1065 | set i 0 |
|---|
| 1066 | set sumdvs 0 |
|---|
| 1067 | set fracout {} |
|---|
| 1068 | set rmsout {} |
|---|
| 1069 | foreach oc $ortholist use $useflag coord $fraclist { |
|---|
| 1070 | #puts "ortho: $oc" |
|---|
| 1071 | set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end] |
|---|
| 1072 | lappend fracout $frac |
|---|
| 1073 | if {$use} {incr i} |
|---|
| 1074 | set dvs 0 |
|---|
| 1075 | foreach var {x y z} v1 $frac v2 $coord abc [lrange $cell 0 2] { |
|---|
| 1076 | set dv [expr {($v2 - $v1)*$abc}] |
|---|
| 1077 | set dvs [expr {$dvs + $dv*$dv}] |
|---|
| 1078 | set sumdvs [expr {$sumdvs + $dv*$dv}] |
|---|
| 1079 | if {$use} {set sum($var) [expr {$sum($var) + $dv/$abc}]} |
|---|
| 1080 | } |
|---|
| 1081 | lappend rmsout [expr {sqrt($dvs)}] |
|---|
| 1082 | } |
|---|
| 1083 | set rms 0 |
|---|
| 1084 | if {$i > 1} {set rms [expr {sqrt($sumdvs)/$i}]} |
|---|
| 1085 | set neworig {} |
|---|
| 1086 | foreach var {x y z} o $origin { |
|---|
| 1087 | lappend neworig [expr {$o + ($sum($var)/$i)}] |
|---|
| 1088 | } |
|---|
| 1089 | return [list $neworig $rms $fracout $rmsout] |
|---|
| 1090 | } |
|---|
| 1091 | |
|---|
| 1092 | # fit a rigid body's Euler angles using least-squares |
|---|
| 1093 | # arguments: |
|---|
| 1094 | # Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...} |
|---|
| 1095 | # where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes |
|---|
| 1096 | # cell: a list with "a b c alpha beta gamma" in Angstroms and degrees |
|---|
| 1097 | # ortholist: list containing triplets with orthogonal coordinates |
|---|
| 1098 | # useflag: list of flags to indicate if an atom should be used (1) or ignored (0) |
|---|
| 1099 | # fraclist: list containing triplets with fractional coordinates |
|---|
| 1100 | # origin: a list of fraction coordinates {x y z} describing the location of the |
|---|
| 1101 | # origin of the orthogonal coordinates in the crystal system |
|---|
| 1102 | # note that the length of ortholist, useflag and fraclist should be the same |
|---|
| 1103 | # Returns a list of new Euler angles |
|---|
| 1104 | proc FitBodyRot {Euler cell ortholist useflag fraclist origin} { |
|---|
| 1105 | set xform [CalcXformMatrix $Euler $cell] |
|---|
| 1106 | set derivlist [CalcDerivMatrix $Euler $cell] |
|---|
| 1107 | set A "2 [expr 3*[llength $ortholist]] 3" |
|---|
| 1108 | foreach oc $ortholist use $useflag coord $fraclist { |
|---|
| 1109 | if {! $use} continue |
|---|
| 1110 | foreach deriv $derivlist { |
|---|
| 1111 | foreach xyz [lrange [Ortho2Xtal $deriv "0 0 0" $oc] 3 end] { |
|---|
| 1112 | lappend A $xyz |
|---|
| 1113 | } |
|---|
| 1114 | } |
|---|
| 1115 | } |
|---|
| 1116 | #puts "A: [La::show $A]" |
|---|
| 1117 | set y "2 [expr 3*[llength $ortholist]] 1" |
|---|
| 1118 | foreach oc $ortholist use $useflag coord $fraclist { |
|---|
| 1119 | if {! $use} continue |
|---|
| 1120 | set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end] |
|---|
| 1121 | foreach xyz $coord XYZ $frac { |
|---|
| 1122 | lappend y [expr {$XYZ - $xyz}] |
|---|
| 1123 | } |
|---|
| 1124 | } |
|---|
| 1125 | |
|---|
| 1126 | set AtA [La::mmult [La::transpose $A] $A] |
|---|
| 1127 | set Aty [La::mmult [La::transpose $A] $y] |
|---|
| 1128 | |
|---|
| 1129 | set l {} |
|---|
| 1130 | #set shifts {} |
|---|
| 1131 | foreach delta [lrange [La::msolve $AtA $Aty] 3 end] item $Euler { |
|---|
| 1132 | #lappend shifts $delta |
|---|
| 1133 | lappend l "[lindex $item 0] [expr {$delta + [lindex $item 1]}]" |
|---|
| 1134 | } |
|---|
| 1135 | #puts "shifts = $shifts" |
|---|
| 1136 | return $l |
|---|
| 1137 | } |
|---|
| 1138 | |
|---|
| 1139 | # FitBody: fit a rigid body's Origin and Euler angles |
|---|
| 1140 | # arguments: |
|---|
| 1141 | # Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...} |
|---|
| 1142 | # where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes |
|---|
| 1143 | # cell: a list with "a b c alpha beta gamma" in Angstroms and degrees |
|---|
| 1144 | # ortholist: list containing triplets with orthogonal coordinates |
|---|
| 1145 | # useflag: list of flags to indicate if an atom should be used (1) or ignored (0) |
|---|
| 1146 | # fraclist: list containing triplets with fractional coordinates |
|---|
| 1147 | # origin: a list of fraction coordinates {x y z} describing the location of the |
|---|
| 1148 | # origin of the orthogonal coordinates in the crystal system |
|---|
| 1149 | # note that the length of ortholist, useflag and fraclist should be the same |
|---|
| 1150 | # Returns a list containing |
|---|
| 1151 | # new origin |
|---|
| 1152 | # new Euler angles |
|---|
| 1153 | # total rms |
|---|
| 1154 | # fractional coordinates |
|---|
| 1155 | # rms deviation in fractional coordinates of new Euler angles |
|---|
| 1156 | proc FitBody {Euler cell ortholist useflag fraclist origin "ncycle 5"} { |
|---|
| 1157 | #puts "start origin = $origin" |
|---|
| 1158 | foreach { |
|---|
| 1159 | origin |
|---|
| 1160 | startrms |
|---|
| 1161 | fracout |
|---|
| 1162 | rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {} |
|---|
| 1163 | #puts "start rms = $startrms" |
|---|
| 1164 | set rmsprev $startrms |
|---|
| 1165 | #puts "new origin = $origin" |
|---|
| 1166 | for {set i 0} {$i < $ncycle} {incr i} { |
|---|
| 1167 | set Eulerprev $Euler |
|---|
| 1168 | set Euler [FitBodyRot $Euler $cell $ortholist $useflag $fraclist $origin] |
|---|
| 1169 | #puts "New Euler $Euler" |
|---|
| 1170 | #puts "after fit" |
|---|
| 1171 | foreach { |
|---|
| 1172 | origin |
|---|
| 1173 | rms |
|---|
| 1174 | fracout |
|---|
| 1175 | rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {} |
|---|
| 1176 | if {$rms > (1.1 * $rmsprev) + 0.01} { |
|---|
| 1177 | #puts "rms = $rms, new origin = $origin" |
|---|
| 1178 | set rmsprev $rms |
|---|
| 1179 | } |
|---|
| 1180 | } |
|---|
| 1181 | #proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} |
|---|
| 1182 | #return "$neworig $rms $fracout $rmsout" |
|---|
| 1183 | set fmt {"%8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %6.3f"} |
|---|
| 1184 | #foreach fracin $fraclist fraccalc $fracout rmsi $rmsout { |
|---|
| 1185 | #puts "[eval format $fmt $fracin $fraccalc $rmsi]" |
|---|
| 1186 | #} |
|---|
| 1187 | return [list $origin $Euler $rms $fracout $rmsout] |
|---|
| 1188 | } |
|---|
| 1189 | |
|---|
| 1190 | # zmat2coord: convert a z-matrix to a set of cartesian coordinates |
|---|
| 1191 | # a z-matrix is also known as "internal coordinates" or "torsion space" |
|---|
| 1192 | # (see Journal of Computational Chemistry, Vol 26, #10, p. 1063â1068, 2005 or |
|---|
| 1193 | # http://www.cmbi.ru.nl/molden/zmat/zmat.html) |
|---|
| 1194 | # INPUT: |
|---|
| 1195 | # atmlist is a list of ascii lines where each line contains |
|---|
| 1196 | # lbl c1 distance c2 angle c3 torsion |
|---|
| 1197 | # where each atom is computed from the previous where the new atom is: |
|---|
| 1198 | # distance $distance from atom $c1 (angstrom) |
|---|
| 1199 | # angle $angle from $c1--$c2 (degrees) |
|---|
| 1200 | # torsion $torsion from $c1--$c2--$c3 (degrees) |
|---|
| 1201 | # OUTPUT: |
|---|
| 1202 | # zmat2coord returns a list of atom labels and cartesian coordinates, |
|---|
| 1203 | # with 4 items in each element (label, x, y, z) |
|---|
| 1204 | # this routine was tested against results from Babel via the web interface at |
|---|
| 1205 | # http://www.shodor.org/chemviz/zmatrices/babel.html and sample input at |
|---|
| 1206 | # http://iopenshell.usc.edu/howto/zmatrix/ |
|---|
| 1207 | proc zmat2coord {atmlist} { |
|---|
| 1208 | set torad [expr {acos(0)/90.}] |
|---|
| 1209 | set i 0 |
|---|
| 1210 | foreach line $atmlist { |
|---|
| 1211 | incr i |
|---|
| 1212 | foreach {lbl c1 dist c2 angle c3 torsion} $line {} |
|---|
| 1213 | if {$i == 1} { |
|---|
| 1214 | set atm(1) [list $lbl 0 0 0] ; # 1st atom is at origin |
|---|
| 1215 | } elseif {$i == 2} { |
|---|
| 1216 | set dist1 $dist |
|---|
| 1217 | set atm(2) [list $lbl $dist1 0 0] ; # 2nd atom is along x-axis |
|---|
| 1218 | } elseif {$i == 3} { |
|---|
| 1219 | # 3rd atom can be bonded to the 1st or 2nd |
|---|
| 1220 | if {$c1 == 1} { |
|---|
| 1221 | set atm(3) [list $lbl \ |
|---|
| 1222 | [expr {$dist * cos($torad * $angle)}] \ |
|---|
| 1223 | [expr {$dist * sin($torad * $angle)}] \ |
|---|
| 1224 | 0] |
|---|
| 1225 | } else { |
|---|
| 1226 | set atm(3) [list $lbl \ |
|---|
| 1227 | [expr {$dist1 - $dist * cos($torad * $angle)}] \ |
|---|
| 1228 | [expr {$dist * sin($torad * $angle)}] \ |
|---|
| 1229 | 0] |
|---|
| 1230 | } |
|---|
| 1231 | } else { |
|---|
| 1232 | set atm($i) [concat $lbl \ |
|---|
| 1233 | [ahcat "atm" $c1 $dist $c2 $angle $c3 $torsion]] |
|---|
| 1234 | } |
|---|
| 1235 | } |
|---|
| 1236 | set coordlist {} |
|---|
| 1237 | foreach key [lsort -integer [array names atm]] { |
|---|
| 1238 | lappend coordlist $atm($key) |
|---|
| 1239 | } |
|---|
| 1240 | return $coordlist |
|---|
| 1241 | } |
|---|
| 1242 | # Compute the length of a vector |
|---|
| 1243 | proc vlen {a} { |
|---|
| 1244 | set sum 0.0 |
|---|
| 1245 | foreach ai $a { |
|---|
| 1246 | set sum [expr {$sum + $ai*$ai}] |
|---|
| 1247 | } |
|---|
| 1248 | return [expr sqrt($sum)] |
|---|
| 1249 | } |
|---|
| 1250 | # compute vector (a + z * b) and optionally normalize to length d |
|---|
| 1251 | proc vadd {a b d z} { |
|---|
| 1252 | set c {} |
|---|
| 1253 | foreach ai $a bi $b { |
|---|
| 1254 | lappend c [expr {$bi + $z * $ai}] |
|---|
| 1255 | } |
|---|
| 1256 | set v [vlen $c] |
|---|
| 1257 | if {$d != 0} { |
|---|
| 1258 | set r {} |
|---|
| 1259 | foreach ci $c { |
|---|
| 1260 | lappend r [expr {$d * $ci / $v}] |
|---|
| 1261 | } |
|---|
| 1262 | return [list $v $r] |
|---|
| 1263 | } |
|---|
| 1264 | return [list $v $c] |
|---|
| 1265 | } |
|---|
| 1266 | # normalize a vector |
|---|
| 1267 | proc vnrm {x} { |
|---|
| 1268 | set v [vlen $x] |
|---|
| 1269 | if {abs($v) < 1e-8} {return [list 0 0 0]} |
|---|
| 1270 | set y {} |
|---|
| 1271 | foreach xi $x { |
|---|
| 1272 | lappend y [expr {$xi / $v}] |
|---|
| 1273 | } |
|---|
| 1274 | return $y |
|---|
| 1275 | } |
|---|
| 1276 | # compute the coordinates for an atom that is bonded: |
|---|
| 1277 | # distance $dist from atom $nc |
|---|
| 1278 | # angle $bondang from $nc--$nb |
|---|
| 1279 | # torsion $torsang from $nc--$nb--$na |
|---|
| 1280 | # coordinates are found in array $atmarr in the calling routine |
|---|
| 1281 | # based on a Fortran routine provided by Peter Zavalij (Thanks Peter!) |
|---|
| 1282 | proc ahcat {atmarr nc dist nb bondang na torsang} { |
|---|
| 1283 | upvar 1 $atmarr atm |
|---|
| 1284 | set xa [lrange $atm($na) 1 3] |
|---|
| 1285 | set xb [lrange $atm($nb) 1 3] |
|---|
| 1286 | set xc [lrange $atm($nc) 1 3] |
|---|
| 1287 | set torad [expr {acos(0)/90.}] |
|---|
| 1288 | # x = unit Vector A-B |
|---|
| 1289 | foreach {x1 x2 x3} [lindex [vadd $xb $xa 1. -1.] 1] {} |
|---|
| 1290 | # y = unit Vector C-B |
|---|
| 1291 | set y [lindex [vadd $xb $xc 1. -1.] 1] |
|---|
| 1292 | foreach {y1 y2 y3} $y {} |
|---|
| 1293 | set z1 [expr {$x2*$y3 - $x3*$y2}] |
|---|
| 1294 | set z2 [expr {$x3*$y1 - $x1*$y3}] |
|---|
| 1295 | set z3 [expr {$x1*$y2 - $x2*$y1}] |
|---|
| 1296 | set z [vnrm [list $z1 $z2 $z3]] |
|---|
| 1297 | set q1 [expr {$y2*$z3 - $y3*$z2}] |
|---|
| 1298 | set q2 [expr {$y3*$z1 - $y1*$z3}] |
|---|
| 1299 | set q3 [expr {$y1*$z2 - $y2*$z1}] |
|---|
| 1300 | set q [vnrm [list $q1 $q2 $q3]] |
|---|
| 1301 | set th [expr {$bondang * $torad}] |
|---|
| 1302 | set ph [expr {-1. * $torsang * $torad}] |
|---|
| 1303 | set cth [expr {cos($th)}] |
|---|
| 1304 | set sth [expr {sin($th)}] |
|---|
| 1305 | set cph [expr {cos($ph)}] |
|---|
| 1306 | set sph [expr {sin($ph)}] |
|---|
| 1307 | set xh {} |
|---|
| 1308 | foreach xci $xc xi $q zi $z yi $y { |
|---|
| 1309 | lappend xh [expr { |
|---|
| 1310 | $xci + |
|---|
| 1311 | $dist*($sth*($cph*$xi + $sph*$zi)-$cth*$yi) |
|---|
| 1312 | }] |
|---|
| 1313 | } |
|---|
| 1314 | return $xh |
|---|
| 1315 | } |
|---|
| 1316 | |
|---|
| 1317 | # RB2cart: convert the rigid body representation reported as the 2nd element |
|---|
| 1318 | # in ReadRigidBody into cartesian coordinates |
|---|
| 1319 | # rblist: a list containing an element for each scaling factor |
|---|
| 1320 | # in each element there are four items: |
|---|
| 1321 | # a multiplier value for the rigid body coordinates |
|---|
| 1322 | # a damping value (0-9) for the refinement of the multiplier (not used) |
|---|
| 1323 | # a variable number if the multiplier will be refined (not used) |
|---|
| 1324 | # a list of cartesian coordinates coordinates |
|---|
| 1325 | # each cartesian coordinate contains 4 items: x,y,z and a label |
|---|
| 1326 | # returns a list of coordinate triplets |
|---|
| 1327 | proc RB2cart {rblist} { |
|---|
| 1328 | foreach item $rblist { |
|---|
| 1329 | foreach {mult damp ref coords} $item {} |
|---|
| 1330 | set i 0 |
|---|
| 1331 | foreach xyz $coords { |
|---|
| 1332 | foreach {x y z} [lrange $xyz 0 2] {} |
|---|
| 1333 | foreach val [lrange $xyz 0 2] var {X Y Z} { |
|---|
| 1334 | if {[array names $var $i] == ""} { |
|---|
| 1335 | set ${var}($i) [expr {$mult * $val}] |
|---|
| 1336 | } else { |
|---|
| 1337 | set ${var}($i) [expr {[set ${var}($i)] + $mult * $val}] |
|---|
| 1338 | } |
|---|
| 1339 | } |
|---|
| 1340 | incr i |
|---|
| 1341 | } |
|---|
| 1342 | } |
|---|
| 1343 | set out "" |
|---|
| 1344 | foreach i [lsort -integer [array names X]] { |
|---|
| 1345 | lappend out [list $X($i) $Y($i) $Z($i)] |
|---|
| 1346 | } |
|---|
| 1347 | return $out |
|---|
| 1348 | } |
|---|
| 1349 | |
|---|
| 1350 | # get the name of the DRAWxtl application, if installed |
|---|
| 1351 | proc GetDRAWxtlApp {} { |
|---|
| 1352 | # is DRAWxtl installed? |
|---|
| 1353 | set app {} |
|---|
| 1354 | if {![catch {set fp [open [file join $::env(HOME) .drawxtlrc] r]}]} { |
|---|
| 1355 | # line 12 is name of executable |
|---|
| 1356 | set i 0 |
|---|
| 1357 | while {$i < 12} { |
|---|
| 1358 | incr i |
|---|
| 1359 | gets $fp appname |
|---|
| 1360 | } |
|---|
| 1361 | close $fp |
|---|
| 1362 | set app [auto_execok $appname] |
|---|
| 1363 | } |
|---|
| 1364 | if {$app != ""} { |
|---|
| 1365 | return $appname |
|---|
| 1366 | } |
|---|
| 1367 | return "" |
|---|
| 1368 | } |
|---|
| 1369 | |
|---|
| 1370 | # DRAWxtlPlotOrtho: plot orthogonal coordinates in DRAWxtl |
|---|
| 1371 | # input: |
|---|
| 1372 | # filename: file name for the .str file to create |
|---|
| 1373 | # title: string for title in .str file |
|---|
| 1374 | # coords: cartesian coordinates |
|---|
| 1375 | # bondlist: list of bonds to draw as min, max length (A) and |
|---|
| 1376 | # an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}} |
|---|
| 1377 | proc DRAWxtlPlotOrtho {filename title coords bondlist} { |
|---|
| 1378 | foreach {xmin ymin zmin} {"" "" ""} {} |
|---|
| 1379 | foreach {xmax ymax zmax} {"" "" ""} {} |
|---|
| 1380 | foreach xyz $coords { |
|---|
| 1381 | foreach {x y z} $xyz {} |
|---|
| 1382 | foreach s {x y z} { |
|---|
| 1383 | foreach t {min max} { |
|---|
| 1384 | if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]} |
|---|
| 1385 | } |
|---|
| 1386 | if {[set ${s}min] > [set $s]} {set ${s}min [set $s]} |
|---|
| 1387 | if {[set ${s}max] < [set $s]} {set ${s}max [set $s]} |
|---|
| 1388 | } |
|---|
| 1389 | } |
|---|
| 1390 | #puts "$xmin $xmax $ymin $ymax $zmin $zmax" |
|---|
| 1391 | set max $xmin |
|---|
| 1392 | foreach val "$xmin $xmax $ymin $ymax $zmin $zmax" { |
|---|
| 1393 | if {$max < abs($val)} {set max $val} |
|---|
| 1394 | } |
|---|
| 1395 | set scale [expr {4.*$max}] |
|---|
| 1396 | set a 10. |
|---|
| 1397 | lappend range [expr -0.01+($xmin/$scale)] [expr 0.01+($xmax/$scale)] \ |
|---|
| 1398 | [expr -0.01+($ymin/$scale)] [expr 0.01+($ymax/$scale)] \ |
|---|
| 1399 | [expr -0.01+($zmin/$scale)] [expr 0.01+($zmax/$scale)] |
|---|
| 1400 | set fp [open $filename w] |
|---|
| 1401 | puts $fp "title $title" |
|---|
| 1402 | puts $fp "box 0.000 Black" |
|---|
| 1403 | puts $fp "background White" |
|---|
| 1404 | puts $fp "nolabels" |
|---|
| 1405 | puts $fp "cell $a $a $a 90 90 90" |
|---|
| 1406 | puts $fp "spgr P 1" |
|---|
| 1407 | puts $fp "pack $range" |
|---|
| 1408 | set i 0 |
|---|
| 1409 | foreach xyz $coords { |
|---|
| 1410 | foreach {x y z} $xyz {} |
|---|
| 1411 | incr i |
|---|
| 1412 | puts $fp "atom c $i [expr {$x/$scale}] [expr {$y/$scale}] [expr {$z/$scale}]" |
|---|
| 1413 | } |
|---|
| 1414 | puts $fp "sphere c [expr 0.100*($a/$scale)] Red" |
|---|
| 1415 | puts $fp "finish 0.70 0.30 0.08 0.01" |
|---|
| 1416 | foreach bondpair $bondlist { |
|---|
| 1417 | foreach {b1 b2 color} $bondpair {} |
|---|
| 1418 | if {$color == ""} {set color Red} |
|---|
| 1419 | puts $fp "bond c c [expr {0.01*$a/$scale}] [expr {$b1*$a/$scale}] [expr {$b2*$a/$scale}] $color" |
|---|
| 1420 | } |
|---|
| 1421 | puts $fp "frame" |
|---|
| 1422 | set range {} |
|---|
| 1423 | lappend range -0.01 [expr 0.01+(0.1*$a/$scale)] \ |
|---|
| 1424 | -0.01 [expr 0.01+(0.1*$a/$scale)] \ |
|---|
| 1425 | -0.01 [expr 0.01+(0.1*$a/$scale)] |
|---|
| 1426 | puts $fp "cell $a $a $a 90 90 90" |
|---|
| 1427 | puts $fp "spgr P 1" |
|---|
| 1428 | puts $fp "pack $range" |
|---|
| 1429 | puts $fp "atom o 1 0 0 0" |
|---|
| 1430 | puts $fp "atom o 2 [expr {0.1*$a/$scale}] 0 0" |
|---|
| 1431 | puts $fp "atom o 3 0 [expr {0.1*$a/$scale}] 0" |
|---|
| 1432 | puts $fp "atom o 4 0 0 [expr {0.1*$a/$scale}]" |
|---|
| 1433 | puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black" |
|---|
| 1434 | puts $fp "sphere o [expr {0.02*$a/$scale}] Blue" |
|---|
| 1435 | puts $fp "origin .0 .0 .0" |
|---|
| 1436 | puts $fp "end" |
|---|
| 1437 | close $fp |
|---|
| 1438 | } |
|---|
| 1439 | |
|---|
| 1440 | # PlotRBtype: plot a rigid body in DRAWxtl |
|---|
| 1441 | # input: |
|---|
| 1442 | # rbtype: # of rigid body |
|---|
| 1443 | # bondlist: list of bonds to draw as min, max length (A) and |
|---|
| 1444 | # an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}} |
|---|
| 1445 | # file: file name for the .str file to create |
|---|
| 1446 | proc PlotRBtype {rbtype "bondlist {}" "file {}"} { |
|---|
| 1447 | set app [GetDRAWxtlApp] |
|---|
| 1448 | if {$app == ""} { |
|---|
| 1449 | MyMessageBox -parent . -title "No DRAWxtl" \ |
|---|
| 1450 | -message "Sorry, DRAWxtl is not installed no phases are present to write" \ |
|---|
| 1451 | -icon warning |
|---|
| 1452 | return |
|---|
| 1453 | } |
|---|
| 1454 | if {$::tcl_platform(platform) == "windows" && $file == ""} { |
|---|
| 1455 | set file [file join [pwd] rbplot.str] |
|---|
| 1456 | } else { |
|---|
| 1457 | set file "/tmp/rbplot.str" |
|---|
| 1458 | } |
|---|
| 1459 | set coords [RB2cart [lindex [ReadRigidBody $rbtype] 1]] |
|---|
| 1460 | DRAWxtlPlotOrtho $file "" $coords $bondlist |
|---|
| 1461 | if {$app != ""} {exec $app $file &} |
|---|
| 1462 | } |
|---|
| 1463 | |
|---|
| 1464 | # PlotRBcoords: plot orthogonal coordinates in DRAWxtl |
|---|
| 1465 | # input: |
|---|
| 1466 | # coords: cartesian coordinates |
|---|
| 1467 | # bondlist: list of bonds to draw as min, max length (A) and |
|---|
| 1468 | # an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}} |
|---|
| 1469 | # file: file name for the .str file to create |
|---|
| 1470 | proc PlotRBcoords {coords "bondlist {}" "file {}"} { |
|---|
| 1471 | set app [GetDRAWxtlApp] |
|---|
| 1472 | if {$app == ""} { |
|---|
| 1473 | MyMessageBox -parent . -title "No DRAWxtl" \ |
|---|
| 1474 | -message "Sorry, DRAWxtl is not installed no phases are present to write" \ |
|---|
| 1475 | -icon warning |
|---|
| 1476 | return |
|---|
| 1477 | } |
|---|
| 1478 | if {$::tcl_platform(platform) == "windows" && $file == ""} { |
|---|
| 1479 | set file [file join [pwd] rbplot.str] |
|---|
| 1480 | } else { |
|---|
| 1481 | set file "/tmp/rbplot.str" |
|---|
| 1482 | } |
|---|
| 1483 | DRAWxtlPlotOrtho $file "" $coords $bondlist |
|---|
| 1484 | if {$app != ""} {exec $app $file &} |
|---|
| 1485 | } |
|---|
| 1486 | |
|---|
| 1487 | # DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed |
|---|
| 1488 | # on a structure read from a phase |
|---|
| 1489 | # input: |
|---|
| 1490 | # RBcoords: fractional coordinates for rigid body |
|---|
| 1491 | # phase:# of phase to plot |
|---|
| 1492 | # firstatom: seq # of 1st atom in structure to be mapped to rigid body |
|---|
| 1493 | # allatoms: 0 to plot only atoms in phase that are in the rigid body, |
|---|
| 1494 | # otherwise plot all atoms |
|---|
| 1495 | # bondlist: list of bonds to draw for the phase as min, max length (A) and |
|---|
| 1496 | # an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}} |
|---|
| 1497 | # rbbondlist: list of bonds to draw for the phase as min, max length (A) and |
|---|
| 1498 | # an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}} |
|---|
| 1499 | # file: optional file name for the .str file to create |
|---|
| 1500 | proc DRAWxtlPlotRBFit {RBcoords phase firstatom "allatoms 0" \ |
|---|
| 1501 | "bondlist {}" "rbbondlist {}" "file {}"} { |
|---|
| 1502 | set natom [llength $RBcoords] |
|---|
| 1503 | set app [GetDRAWxtlApp] |
|---|
| 1504 | if {$app == ""} { |
|---|
| 1505 | MyMessageBox -parent . -title "No DRAWxtl" \ |
|---|
| 1506 | -message "Sorry, DRAWxtl is not installed no phases are present to write" \ |
|---|
| 1507 | -icon warning |
|---|
| 1508 | return |
|---|
| 1509 | } |
|---|
| 1510 | if {$::tcl_platform(platform) == "windows" && $file == ""} { |
|---|
| 1511 | set file [file join [pwd] rbplot.str] |
|---|
| 1512 | } else { |
|---|
| 1513 | set file "/tmp/rbfit.str" |
|---|
| 1514 | } |
|---|
| 1515 | |
|---|
| 1516 | # get rigid body coordinate range |
|---|
| 1517 | foreach {xmin ymin zmin} {"" "" ""} {} |
|---|
| 1518 | foreach {xmax ymax zmax} {"" "" ""} {} |
|---|
| 1519 | foreach xyz $RBcoords { |
|---|
| 1520 | foreach {x y z} $xyz {} |
|---|
| 1521 | foreach s {x y z} { |
|---|
| 1522 | foreach t {min max} { |
|---|
| 1523 | if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]} |
|---|
| 1524 | } |
|---|
| 1525 | if {[set ${s}min] > [set $s]} {set ${s}min [set $s]} |
|---|
| 1526 | if {[set ${s}max] < [set $s]} {set ${s}max [set $s]} |
|---|
| 1527 | } |
|---|
| 1528 | } |
|---|
| 1529 | set rbrange {} |
|---|
| 1530 | foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \ |
|---|
| 1531 | [expr -0.01+$ymin] [expr 0.01+$ymax] \ |
|---|
| 1532 | [expr -0.01+$zmin] [expr 0.01+$zmax] ] { |
|---|
| 1533 | append rbrange [format " %8.4f" $val] |
|---|
| 1534 | } |
|---|
| 1535 | set rbcenter [list [expr {($xmin+$xmax)/2}] \ |
|---|
| 1536 | [expr {($ymin+$ymax)/2}] \ |
|---|
| 1537 | [expr {($zmin+$zmax)/2}] ] |
|---|
| 1538 | # get matching atoms coordinate range |
|---|
| 1539 | set firstind [lsearch $::expmap(atomlist_$phase) $firstatom] |
|---|
| 1540 | foreach atom [lrange \ |
|---|
| 1541 | [lrange $::expmap(atomlist_$phase) $firstind end] \ |
|---|
| 1542 | 0 [expr {$natom-1}]] { |
|---|
| 1543 | foreach s {x y z} { |
|---|
| 1544 | set $s [atominfo $phase $atom $s] |
|---|
| 1545 | foreach t {min max} { |
|---|
| 1546 | if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]} |
|---|
| 1547 | } |
|---|
| 1548 | if {[set ${s}min] > [set $s]} {set ${s}min [set $s]} |
|---|
| 1549 | if {[set ${s}max] < [set $s]} {set ${s}max [set $s]} |
|---|
| 1550 | } |
|---|
| 1551 | } |
|---|
| 1552 | # expand to cover at least one unit cell |
|---|
| 1553 | foreach var {xmin ymin zmin} { |
|---|
| 1554 | if {[set $var] > 0.0} {set $var 0.0} |
|---|
| 1555 | } |
|---|
| 1556 | foreach var {xmax ymax zmax} { |
|---|
| 1557 | if {[set $var] < 1.} {set $var 1.} |
|---|
| 1558 | } |
|---|
| 1559 | set range {} |
|---|
| 1560 | foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \ |
|---|
| 1561 | [expr -0.01+$ymin] [expr 0.01+$ymax] \ |
|---|
| 1562 | [expr -0.01+$zmin] [expr 0.01+$zmax]] { |
|---|
| 1563 | append range [format " %8.4f" $val] |
|---|
| 1564 | } |
|---|
| 1565 | |
|---|
| 1566 | set fp [open $file w] |
|---|
| 1567 | puts $fp "title structure/rigid-body fit plot" |
|---|
| 1568 | # plot the structure |
|---|
| 1569 | puts -nonewline $fp "cell" |
|---|
| 1570 | foreach p {a b c alpha beta gamma} { |
|---|
| 1571 | puts -nonewline $fp " [phaseinfo $phase $p]" |
|---|
| 1572 | } |
|---|
| 1573 | puts $fp "" |
|---|
| 1574 | puts $fp "spgp [phaseinfo $phase spacegroup]" |
|---|
| 1575 | puts $fp "pack $range" |
|---|
| 1576 | if {$allatoms != 0} { |
|---|
| 1577 | set atoms $::expmap(atomlist_$phase) |
|---|
| 1578 | } else { |
|---|
| 1579 | set firstind [lsearch $::expmap(atomlist_$phase) $firstatom] |
|---|
| 1580 | set atoms [lrange \ |
|---|
| 1581 | [lrange $::expmap(atomlist_$phase) $firstind end] \ |
|---|
| 1582 | 0 [expr {$natom-1}]] |
|---|
| 1583 | } |
|---|
| 1584 | |
|---|
| 1585 | # set origin at center of rigid body |
|---|
| 1586 | puts $fp "origin $rbcenter" |
|---|
| 1587 | # now loop over atoms |
|---|
| 1588 | foreach atom $atoms { |
|---|
| 1589 | set type [atominfo $phase $atom type] |
|---|
| 1590 | set typelist($type) 1 |
|---|
| 1591 | puts -nonewline $fp "atom $type $atom " |
|---|
| 1592 | foreach v {x y z} { |
|---|
| 1593 | puts -nonewline $fp "[atominfo $phase $atom $v] " |
|---|
| 1594 | } |
|---|
| 1595 | puts $fp "" |
|---|
| 1596 | |
|---|
| 1597 | set uiso [atominfo $phase $atom Uiso] |
|---|
| 1598 | # are there anisotropic atoms? If so convert them to Uequiv |
|---|
| 1599 | if {[atominfo $phase $atom temptype] == "A"} { |
|---|
| 1600 | puts -nonewline $fp "Uij [atominfo $phase $atom type] $atom " |
|---|
| 1601 | foreach v {U11 U22 U33 U12 U13 U23} { |
|---|
| 1602 | puts -nonewline $fp "[atominfo $phase $atom $v] " |
|---|
| 1603 | } |
|---|
| 1604 | puts $fp "" |
|---|
| 1605 | } |
|---|
| 1606 | } |
|---|
| 1607 | foreach type [array names typelist] color {Green Blue Magenta Cyan} { |
|---|
| 1608 | if {$type == ""} break |
|---|
| 1609 | puts $fp "sphere $type 0.1 $color" |
|---|
| 1610 | } |
|---|
| 1611 | foreach type [array names typelist] color1 {Green Blue Magenta Cyan} { |
|---|
| 1612 | foreach bondpair $bondlist { |
|---|
| 1613 | foreach {b1 b2 color} $bondpair {} |
|---|
| 1614 | if {$color == ""} {set color $color1} |
|---|
| 1615 | puts $fp "bond $type $type 0.02 $b1 $b2 $color" |
|---|
| 1616 | } |
|---|
| 1617 | foreach type1 [array names typelist] { |
|---|
| 1618 | if {$type1 == $type} break |
|---|
| 1619 | foreach bondpair $bondlist { |
|---|
| 1620 | foreach {b1 b2 color} $bondpair {} |
|---|
| 1621 | if {$color == ""} {set color $color1} |
|---|
| 1622 | puts $fp "bond $type $type1 0.02 $b1 $b2 $color" |
|---|
| 1623 | } |
|---|
| 1624 | } |
|---|
| 1625 | } |
|---|
| 1626 | # plot the rigid body |
|---|
| 1627 | puts $fp "frame" |
|---|
| 1628 | puts -nonewline $fp "cell" |
|---|
| 1629 | foreach p {a b c alpha beta gamma} { |
|---|
| 1630 | puts -nonewline $fp " [phaseinfo $phase $p]" |
|---|
| 1631 | } |
|---|
| 1632 | puts $fp "" |
|---|
| 1633 | puts $fp "background White" |
|---|
| 1634 | puts $fp "nolabels" |
|---|
| 1635 | puts $fp "spgr P 1" |
|---|
| 1636 | puts $fp "pack $rbrange" |
|---|
| 1637 | set i 0 |
|---|
| 1638 | foreach xyz $RBcoords { |
|---|
| 1639 | foreach {x y z} $xyz {} |
|---|
| 1640 | incr i |
|---|
| 1641 | puts $fp "atom c $i $x $y $z" |
|---|
| 1642 | } |
|---|
| 1643 | foreach bondpair $rbbondlist { |
|---|
| 1644 | foreach {b1 b2 color} $bondpair {} |
|---|
| 1645 | if {$color == ""} {set color Red} |
|---|
| 1646 | puts $fp "bond c c 0.02 $b1 $b2 $color" |
|---|
| 1647 | } |
|---|
| 1648 | |
|---|
| 1649 | puts $fp "sphere c 0.05 Red" |
|---|
| 1650 | puts $fp "finish 0.70 0.30 0.08 0.01" |
|---|
| 1651 | puts $fp "end" |
|---|
| 1652 | |
|---|
| 1653 | #puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black" |
|---|
| 1654 | close $fp |
|---|
| 1655 | if {$app != ""} {exec $app $file &} |
|---|
| 1656 | } |
|---|
| 1657 | |
|---|
| 1658 | |
|---|
| 1659 | #AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} } |
|---|
| 1660 | #puts [GetRB 1 6 8 "1 2" "X 1 2" "Y 1 3"] |
|---|
| 1661 | #puts [GetRB 1 4 8 "1" "X 1 2" "Z 3 4"] |
|---|
| 1662 | #MapRigidBody 1 1 7 ".11 .22 .33" "11 12 13" |
|---|
| 1663 | |
|---|
| 1664 | |
|---|
| 1665 | #AddRigidBody {1} { { |
|---|
| 1666 | # {1 1 1 o} {-1 1 1 o} {1 -1 1 o} {-1 -1 1 o} |
|---|
| 1667 | # {1 1 -1 o} {-1 1 -1 o} {1 -1 -1 o} {-1 -1 -1 o} |
|---|
| 1668 | #} } |
|---|
| 1669 | #set n [MapRigidBody 1 1 1 ".2 .3 .4" "13 17 19"] |
|---|
| 1670 | #puts "body $n created" |
|---|
| 1671 | #incr expgui(changed) |
|---|
| 1672 | #RunRecalcRBCoords |
|---|
| 1673 | #puts "press Enter to continue" |
|---|
| 1674 | #gets stdin line |
|---|
| 1675 | #MapRigidBody 1 1 $n ".5 .5 .5" "0 0 0" |
|---|
| 1676 | #incr expgui(changed) |
|---|
| 1677 | #RunRecalcRBCoords |
|---|
| 1678 | |
|---|
| 1679 | #puts "Test FitBody" |
|---|
| 1680 | set fraclist { |
|---|
| 1681 | { 0.5483305238484277 0.4887545024531055 0.6167996784631056 } |
|---|
| 1682 | { 0.1036801409356145 0.5954016321779562 0.5129448102437683 } |
|---|
| 1683 | { 0.26404665760133855 0.09455414439078394 0.612655365147539 } |
|---|
| 1684 | { -0.18060372531147473 0.20120127411563465 0.5088004969282018 } |
|---|
| 1685 | { 0.5806037253114747 0.3987987258843653 0.2911995030717982 } |
|---|
| 1686 | { 0.13595334239866147 0.5054458556092161 0.18734463485246095 } |
|---|
| 1687 | { 0.2963198590643855 0.004598367822043814 0.2870551897562318 } |
|---|
| 1688 | { -0.1483305238484277 0.1112454975468945 0.1832003215368945 } |
|---|
| 1689 | } |
|---|
| 1690 | set ortholist { |
|---|
| 1691 | {1 1 1} |
|---|
| 1692 | {-1 1 1} |
|---|
| 1693 | { 1.000000 -1.000000 1.000000} |
|---|
| 1694 | { -1.000000 -1.000000 1.000000} |
|---|
| 1695 | { 1.000000 1.000000 -1.000000} |
|---|
| 1696 | { -1.000000 1.000000 -1.000000} |
|---|
| 1697 | { 1.000000 -1.000000 -1.000000} |
|---|
| 1698 | { -1.000000 -1.000000 -1.000000} |
|---|
| 1699 | } |
|---|
| 1700 | # test code, generates DRAWxtl imput file from orthogonal coordinate list |
|---|
| 1701 | # with bonds of ~2, 2.8 and 3.4 A |
|---|
| 1702 | #DRAWxtlPlotOrtho test4.str "test file" $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} } |
|---|
| 1703 | |
|---|
| 1704 | # test code, plots rigid body type #2 with bonds drawn at ~1.3 & 2 A |
|---|
| 1705 | #PlotRBtype 2 {{1.9 2.1} {1.28 1.32}} |
|---|
| 1706 | |
|---|
| 1707 | # test code, plots rigid body coords in ortholist with bonds @ ~2, 2.8 and 3.4 A |
|---|
| 1708 | #PlotRBcoords $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} } |
|---|
| 1709 | |
|---|
| 1710 | |
|---|
| 1711 | set useflag {1 1 1 1 1 1 1 1} |
|---|
| 1712 | set cell {4. 5. 6. 95. 100. 105.} |
|---|
| 1713 | #set origin ".20 .30 .40" |
|---|
| 1714 | set origin ".0 .0 .0" |
|---|
| 1715 | #set Euler {{1 13} {2 17} {3 19}} |
|---|
| 1716 | #set Euler {{1 0} {2 180} {3 0}} |
|---|
| 1717 | set Euler {{1 0} {2 0} {3 0}} |
|---|
| 1718 | |
|---|
| 1719 | #puts [La::show $xform] |
|---|
| 1720 | #puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]" |
|---|
| 1721 | |
|---|
| 1722 | |
|---|
| 1723 | # test zmat2coord |
|---|
| 1724 | set atmlist { |
|---|
| 1725 | {C1 0 0.0 0 0.0 0 0.0} |
|---|
| 1726 | {O2 1 1.20 0 0.0 0 0.0} |
|---|
| 1727 | {H3 1 1.10 2 120.0 0 0.0} |
|---|
| 1728 | {C4 1 1.50 2 120.0 3 180.0} |
|---|
| 1729 | {H5 4 1.10 1 110.0 2 0.00} |
|---|
| 1730 | {H6 4 1.10 1 110.0 2 120.0} |
|---|
| 1731 | {H7 4 1.10 1 110.0 2 -120.0} |
|---|
| 1732 | } |
|---|
| 1733 | # C 0.00000 0.00000 0.00000 |
|---|
| 1734 | # O 1.20000 0.00000 0.00000 |
|---|
| 1735 | # H -0.55000 0.95263 0.00000 |
|---|
| 1736 | # C -0.75000 -1.29904 -0.00000 |
|---|
| 1737 | # H -0.04293 -2.14169 -0.00000 |
|---|
| 1738 | # H -1.38570 -1.36644 0.89518 |
|---|
| 1739 | # H -1.38570 -1.36644 -0.89518 |
|---|
| 1740 | # set coordlist [zmat2coord $atmlist] |
|---|
| 1741 | set i 0 |
|---|
| 1742 | # puts "\nZmatrix in" |
|---|
| 1743 | # foreach line $atmlist { |
|---|
| 1744 | # incr i |
|---|
| 1745 | # puts "$i) $line" |
|---|
| 1746 | # } |
|---|
| 1747 | # puts "Cartesian out" |
|---|
| 1748 | # foreach line $coordlist { |
|---|
| 1749 | # puts [eval format "%-4s%10.5f%10.5f%10.5f" $line] |
|---|
| 1750 | # } |
|---|
| 1751 | |
|---|
| 1752 | # AddRigidBody {1 0.75} { |
|---|
| 1753 | # { |
|---|
| 1754 | # {1 1 1 c} |
|---|
| 1755 | # {-1 1 1 c} |
|---|
| 1756 | # { 1.000000 -1.000000 1.000000 c} |
|---|
| 1757 | # { -1.000000 -1.000000 1.000000 c} |
|---|
| 1758 | # { 1.000000 1.000000 -1.000000 c} |
|---|
| 1759 | # { -1.000000 1.000000 -1.000000 c} |
|---|
| 1760 | # { 1.000000 -1.000000 -1.000000 c} |
|---|
| 1761 | # { -1.000000 -1.000000 -1.000000 c} |
|---|
| 1762 | # {1 1 1 h} |
|---|
| 1763 | # {1 -1 -1 h} |
|---|
| 1764 | # {-1 1 -1 h} |
|---|
| 1765 | # {-1 -1 1 h} |
|---|
| 1766 | # } { |
|---|
| 1767 | # {0 0 0 c } |
|---|
| 1768 | # {0 0 0 c} |
|---|
| 1769 | # {0 0 0 c} |
|---|
| 1770 | # {0 0 0 c} |
|---|
| 1771 | # {0 0 0 c} |
|---|
| 1772 | # {0 0 0 c} |
|---|
| 1773 | # {0 0 0 c} |
|---|
| 1774 | # {0 0 0 c} |
|---|
| 1775 | # {1 1 1 h} |
|---|
| 1776 | # {1 -1 -1 h} |
|---|
| 1777 | # {-1 1 -1 h} |
|---|
| 1778 | # {-1 -1 1 h} |
|---|
| 1779 | # } |
|---|
| 1780 | # } |
|---|
| 1781 | # MapRigidBody 2 2 1 {0 0 0} {10 15 20} |
|---|
