Changeset 1107 for branches/sandbox


Ignore:
Timestamp:
Jan 17, 2011 11:05:57 AM (10 years ago)
Author:
toby
Message:

fix rb bugs: ReadRigidBody? needs nested list & SetRigidBody? writes wrong # of translations; Add plot RB capability

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/rb.tcl

    r1106 r1107  
    4040#  not used in GSAS.
    4141proc ReadRigidBody {rbnum} {
    42     if {[RigidBodyCount] > $rbnum} {
     42    if {[RigidBodyCount] < $rbnum} {
    4343        return ""
    4444    }
     
    4949    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
    5050    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
    51     set out $used
     51    set out {}
    5252    for {set i 1} {$i <= $nmult} {incr i} {
    5353        set line [readexp "${key}${i}PARM"]
     
    6464            set z [string trim [string range $line 20 29]]
    6565            set lbl [string trim [string range $line 30 39]]
    66             lappend coordlist [list $x $y $x $lbl]
     66            lappend coordlist [list $x $y $z $lbl]
    6767        }
    6868        lappend out [list $mult $damp $var $coordlist]
    6969    }
    70     return $out
     70    return [list $used $out]
    7171}
    7272
     
    380380    setexp "$key NBDS" $value 1 5
    381381    # number of coordinate matrices
    382     set value [llength multlist]
     382    set value [llength $multlist]
    383383    validint value 5
    384384    makeexprec "$key NSMP"
     
    396396            validint value 3
    397397            makeexprec "${key}${i}SC$value"
    398             setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
     398            if {[llength $item] == 4} {
     399                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
     400            } elseif {[llength $item] == 3} {
     401                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f" $item] 1 30
     402            } else {
     403                return -code 3 "Invalid number of coordinates"
     404            }
    399405        }
    400406    }
     
    800806        set dvs 0
    801807        foreach var {x y z} v1 $frac v2 $coord abc [lrange $cell 0 2] {
    802             set dv [expr {($v2 - $v1)}]
     808            set dv [expr {($v2 - $v1)*$abc}]
    803809            set dvs [expr {$dvs + $dv*$dv}]
    804810            set sumdvs [expr {$sumdvs + $dv*$dv}]
     
    904910            set rmsprev $rms
    905911        }
    906     }   
     912    } 
    907913    #proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin}
    908914    #return "$neworig $rms $fracout $rmsout"
     
    10411047}
    10421048 
    1043 
     1049# convert the rigid body representation reported as the 2nd element in ReadRigidBody
     1050# into cartesian coordinates
     1051#   rblist: a list containing an element for each scaling factor
     1052# in each element there are four items:
     1053#    a multiplier value for the rigid body coordinates
     1054#    a damping value (0-9) for the refinement of the multiplier (not used)
     1055#    a variable number if the multiplier will be refined (not used)
     1056#    a list of cartesian coordinates coordinates
     1057# each cartesian coordinate contains 4 items: x,y,z and a label
     1058# returns a list of coordinate triplets
     1059proc RB2cart {rblist} {
     1060    foreach item $rblist {
     1061        foreach {mult damp ref coords} $item {}
     1062        set i 0
     1063        foreach xyz $coords {
     1064            foreach {x y z} [lrange $xyz 0 2] {}
     1065            foreach val [lrange $xyz 0 2] var {X Y Z} {
     1066                if {[array names $var $i] == ""} {
     1067                    set ${var}($i) [expr {$mult * $val}]
     1068                } else {
     1069                    set ${var}($i) [expr {[set ${var}($i)] + $mult * $val}]
     1070                }
     1071            }
     1072            incr i
     1073        }
     1074    }
     1075    set out ""
     1076    foreach i [lsort -integer [array names X]] {
     1077        lappend out [list $X($i) $Y($i) $Z($i)]
     1078    }
     1079    return $out
     1080}
     1081
     1082proc GetDRAWxtlApp {} {
     1083    # is DRAWxtl installed?
     1084    set app {}
     1085    if {![catch {set fp [open [file join $::env(HOME) .drawxtlrc] r]}]} {
     1086        # line 12 is name of executable
     1087        set i 0
     1088        while {$i < 12} {
     1089            incr i
     1090            gets $fp appname
     1091        }
     1092        close $fp
     1093        set app [auto_execok $appname]
     1094    }
     1095    if {$app != ""} {
     1096        return $appname
     1097    }
     1098    return ""
     1099}
     1100
     1101proc DRAWxtlPlotOrtho {filename title coords bondlist} {
     1102    foreach {xmin ymin zmin} {"" "" ""} {}
     1103    foreach {xmax ymax zmax} {"" "" ""} {}
     1104    foreach xyz $coords {
     1105        foreach {x y z} $xyz {}
     1106        foreach s {x y z} {
     1107            foreach t {min max} {
     1108                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
     1109            }
     1110            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
     1111            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
     1112        }
     1113    }
     1114    puts "$xmin $xmax $ymin $ymax $zmin $zmax"
     1115    set max $xmin
     1116    foreach val "$xmin $xmax $ymin $ymax $zmin $zmax" {
     1117        if {$max < abs($val)} {set max $val}
     1118    }
     1119    set scale [expr {4.*$max}]
     1120    set a 10.
     1121    lappend range [expr -0.01+($xmin/$scale)] [expr 0.01+($xmax/$scale)] \
     1122        [expr -0.01+($ymin/$scale)] [expr 0.01+($ymax/$scale)] \
     1123        [expr -0.01+($zmin/$scale)] [expr 0.01+($zmax/$scale)]
     1124    set fp [open $filename w]
     1125    puts $fp "title $fp"
     1126    puts $fp "box  0.000 Black"
     1127    puts $fp "background White"
     1128    puts $fp "nolabels"
     1129    puts $fp "cell $a $a $a 90 90 90"
     1130    puts $fp "spgr P 1"
     1131    puts $fp "pack $range"
     1132    set i 0
     1133    foreach xyz $coords {
     1134        foreach {x y z} $xyz {}
     1135        incr i
     1136        puts $fp "atom c $i [expr {$x/$scale}] [expr {$y/$scale}] [expr {$z/$scale}]"
     1137    }
     1138    puts $fp "sphere c  [expr 0.100*($a/$scale)] Red"
     1139    puts $fp "finish   0.70   0.30   0.08   0.01"
     1140    foreach bondpair $bondlist {
     1141        foreach {b1 b2 color} $bondpair {}
     1142        if {$color == ""} {set color Red}
     1143        puts $fp "bond c c [expr {0.01*$a/$scale}] [expr {$b1*$a/$scale}] [expr {$b2*$a/$scale}] $color"
     1144    }
     1145    puts $fp "frame"
     1146    set range {}
     1147    lappend range -0.01 [expr 0.01+(0.1*$a/$scale)] \
     1148        -0.01 [expr 0.01+(0.1*$a/$scale)] \
     1149        -0.01 [expr 0.01+(0.1*$a/$scale)]
     1150    puts $fp "cell $a $a $a 90 90 90"
     1151    puts $fp "spgr P 1"
     1152    puts $fp "pack $range"
     1153    puts $fp "atom o 1 0 0 0"
     1154    puts $fp "atom o 2 [expr {0.1*$a/$scale}] 0 0"
     1155    puts $fp "atom o 3 0 [expr {0.1*$a/$scale}] 0"
     1156    puts $fp "atom o 4 0 0 [expr {0.1*$a/$scale}]"
     1157    puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
     1158    puts $fp "sphere o [expr {0.02*$a/$scale}] Blue"
     1159    puts $fp "origin   .0 .0 .0"
     1160    puts $fp "end"
     1161    close $fp
     1162}
     1163
     1164proc PlotRBtype {rbtype "bondlist {}" "file {}"} {
     1165    set app [GetDRAWxtlApp]
     1166    if {$app == ""} {
     1167        MyMessageBox -parent . -title "No DRAWxtl" \
     1168                -message "Sorry, DRAWxtl is not installed no phases are present to write" \
     1169                -icon warning
     1170        return
     1171    }
     1172    if {$::tcl_platform(platform) == "windows" && $file == ""} {
     1173        set file [file join [pwd] rbplot.str]
     1174    } else {
     1175        set file "/tmp/rbplot.str"
     1176    }
     1177    set coords [RB2cart [lindex [ReadRigidBody $rbtype] 1]]
     1178puts $coords
     1179    DRAWxtlPlotOrtho $file "" $coords $bondlist
     1180    if {$app != ""} [exec $app $file]
     1181}
    10441182
    10451183#AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
     
    10631201#RunRecalcRBCoords
    10641202
    1065 puts "Test FitBody"
     1203#puts "Test FitBody"
    10661204set fraclist {
    10671205    { 0.5483305238484277 0.4887545024531055 0.6167996784631056 }
     
    10841222    {     -1.000000   -1.000000   -1.000000}
    10851223}
     1224# test code, generates DRAWxtl imput file from orthogonal coordinate list
     1225# with bonds of ~2, 2.8 and 3.4 A
     1226#DRAWxtlPlotOrtho test4.str "test file" $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
     1227
     1228# test code, plots rigid body type #2 with bonds drawn at ~1.3 & 2 A
     1229#PlotRBtype 2 {{1.9 2.1} {1.28 1.32}}
    10861230
    10871231set useflag {1 1 1 1 1 1 1 1}
     
    10941238
    10951239#puts [La::show $xform]
    1096 puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
     1240#puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
    10971241
    10981242
     
    11141258#  H       -1.38570       -1.36644        0.89518
    11151259#  H       -1.38570       -1.36644       -0.89518
    1116 set coordlist [zmat2coord $atmlist]
    1117 set i 0
    1118 puts "\nZmatrix in"
    1119 foreach line $atmlist {
    1120     incr i
    1121     puts "$i) $line"
    1122 }
    1123 puts "Cartesian out"
    1124 foreach line $coordlist {
    1125     puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
    1126 }
     1260# set coordlist [zmat2coord $atmlist]
     1261 set i 0
     1262# puts "\nZmatrix in"
     1263# foreach line $atmlist {
     1264#     incr i
     1265#     puts "$i) $line"
     1266# }
     1267# puts "Cartesian out"
     1268# foreach line $coordlist {
     1269#     puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
     1270# }
     1271
     1272# AddRigidBody {1 0.75} {
     1273#     {
     1274#       {1 1 1 c}
     1275#       {-1 1 1 c}
     1276#       {      1.000000   -1.000000    1.000000 c}
     1277#       {     -1.000000   -1.000000    1.000000 c}
     1278#       {      1.000000    1.000000   -1.000000 c}
     1279#       {     -1.000000    1.000000   -1.000000 c}
     1280#       {      1.000000   -1.000000   -1.000000 c}
     1281#       {     -1.000000   -1.000000   -1.000000 c}
     1282#       {1 1 1 h}
     1283#       {1 -1 -1 h}
     1284#       {-1 1 -1 h}
     1285#       {-1 -1 1 h}
     1286#     } {
     1287#       {0 0 0 c }
     1288#       {0 0 0 c}
     1289#       {0 0 0 c}
     1290#       {0 0 0 c}
     1291#       {0 0 0 c}
     1292#       {0 0 0 c}
     1293#       {0 0 0 c}
     1294#       {0 0 0 c}
     1295#       {1 1 1 h}
     1296#       {1 -1 -1 h}
     1297#       {-1 1 -1 h}
     1298#       {-1 -1 1 h}
     1299#     }
     1300# }
     1301# MapRigidBody 2 2 1 {0 0 0} {10 15 20}
Note: See TracChangeset for help on using the changeset viewer.