Changeset 1108 for branches/sandbox


Ignore:
Timestamp:
Jan 18, 2011 3:43:29 PM (10 years ago)
Author:
toby
Message:

relabel xform atoms; more RB stuff

Location:
branches/sandbox
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/expgui

    r1104 r1108  
    996996
    997997    set crsPhase {}
    998     $expgui(atomxform) config -text "Xform Atoms" -state disabled
     998    $expgui(atomxform) config -text "Modify Atoms" -state disabled
    999999    foreach n $expmap(phaselist) type $expmap(phasetype) {
    10001000        if {$n == $num} {
     
    13821382    DisplayRefFlags $atomnum $p
    13831383    $expgui(EditingAtoms) config -text "Editing atom #$atomnum -- [atominfo $p $atomnum label]"
    1384     $expgui(atomxform) config -text "Xform Atom" -state normal
     1384    $expgui(atomxform) config -text "Modify Atom" -state normal
    13851385}
    13861386
     
    14071407    # this needs to track by phase
    14081408    DisplayRefFlags $numberList $p
    1409     $expgui(atomxform) config -text "Xform Atoms" -state normal
     1409    $expgui(atomxform) config -text "Modify Atoms" -state normal
    14101410}
    14111411
  • branches/sandbox/rb.tcl

    r1107 r1108  
    1818# MapRigidBody -- map a rigid body type into a phase
    1919# EditRigidBodyMapping -- change the parameters in a rigid body mapping
     20#  need to unmap a rigid body
     21#  delete rigid body
    2022#============================================================================
    2123# returns the number of defined rigid bodies
     
    7476proc RigidBodyMappingCount {phase bodytyp} {
    7577    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
     78    if {! [existsexp "$key  NBDS"]} {return 0}
    7679    set n [string trim [readexp "$key  NBDS"]]
    7780    if {$n == ""} {
     
    773776}
    774777
     778#calc rigid body fractional coordinates
     779# arguments:
     780#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
     781#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
     782#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
     783#  ortholist: list containing triplets with orthogonal coordinates
     784#  origin: a list of fraction coordinates {x y z} describing the location of the
     785#            origin of the orthogonal coordinates in the crystal system
     786#     note that the length of ortholist, useflag and fraclist should be the same
     787# Returns a list with the computed fractional coordinates for all atoms
     788proc CalcBody {Euler cell ortholist origin} {
     789    set xform [CalcXformMatrix $Euler $cell]
     790    set i 0
     791    set sumdvs 0
     792    set fracout {}
     793    set rmsout {}
     794    foreach oc $ortholist {
     795        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
     796        lappend fracout $frac
     797    }
     798    return $fracout
     799}
     800
     801
    775802# fit a rigid body's origin
    776803# arguments:
     
    809836            set dvs [expr {$dvs + $dv*$dv}]
    810837            set sumdvs [expr {$sumdvs + $dv*$dv}]
    811             if {$use} {set sum($var) [expr {$sum($var) + $dv}]}
     838            if {$use} {set sum($var) [expr {$sum($var) + $dv/$abc}]}
    812839        }
    813840        lappend rmsout [expr {sqrt($dvs)}]
     
    887914#   rms deviation in fractional coordinates of new Euler angles
    888915proc FitBody {Euler cell ortholist useflag fraclist origin "ncycle 5"} {
    889     puts "start origin = $origin"
     916    #puts "start origin = $origin"
    890917    foreach {
    891918        origin
     
    893920        fracout
    894921        rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
    895     puts "start rms = $startrms"
     922    #puts "start rms = $startrms"
    896923    set rmsprev $startrms
    897     puts "new origin = $origin"
     924    #puts "new origin = $origin"
    898925    for {set i 0} {$i < $ncycle} {incr i} {
    899926        set Eulerprev $Euler
    900927        set Euler [FitBodyRot $Euler $cell $ortholist $useflag $fraclist $origin]
    901         puts "New Euler $Euler"
    902         puts "after fit"
     928        #puts "New Euler $Euler"
     929        #puts "after fit"
    903930        foreach {
    904931            origin
     
    907934            rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
    908935        if {$rms > (1.1 * $rmsprev) + 0.01} {
    909             puts "rms = $rms, new origin = $origin"
     936            #puts "rms = $rms, new origin = $origin"
    910937            set rmsprev $rms
    911938        }
     
    914941    #return "$neworig $rms $fracout $rmsout"
    915942    set fmt  {"%8.5f %8.5f %8.5f     %8.5f %8.5f %8.5f   %6.3f"}
    916     foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
    917         puts "[eval format $fmt $fracin $fraccalc $rmsi]"
    918     }
     943    #foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
     944        #puts "[eval format $fmt $fracin $fraccalc $rmsi]"
     945    #}
    919946    return [list $origin $Euler $rms $fracout $rmsout]
    920947}
     
    10801107}
    10811108
     1109# get the name of the DRAWxtl application, if installed
    10821110proc GetDRAWxtlApp {} {
    10831111    # is DRAWxtl installed?
     
    10991127}
    11001128
     1129# plot orthogonal coordinates in DRAWxtl
     1130# input:
     1131#  filename: file name for the .str file to create
     1132#  title: string for title in .str file
     1133#  coords: cartesian coordinates
     1134#  bondlist: list of bonds to draw as min, max length (A) and
     1135#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
    11011136proc DRAWxtlPlotOrtho {filename title coords bondlist} {
    11021137    foreach {xmin ymin zmin} {"" "" ""} {}
     
    11121147        }
    11131148    }
    1114     puts "$xmin $xmax $ymin $ymax $zmin $zmax"
     1149    #puts "$xmin $xmax $ymin $ymax $zmin $zmax"
    11151150    set max $xmin
    11161151    foreach val "$xmin $xmax $ymin $ymax $zmin $zmax" {
     
    11231158        [expr -0.01+($zmin/$scale)] [expr 0.01+($zmax/$scale)]
    11241159    set fp [open $filename w]
    1125     puts $fp "title $fp"
     1160    puts $fp "title $title"
    11261161    puts $fp "box  0.000 Black"
    11271162    puts $fp "background White"
     
    11621197}
    11631198
     1199# plot a rigid body in DRAWxtl
     1200# input:
     1201#  rbtype: # of rigid body
     1202#  bondlist: list of bonds to draw as min, max length (A) and
     1203#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
     1204#  file: file name for the .str file to create
    11641205proc PlotRBtype {rbtype "bondlist {}" "file {}"} {
    11651206    set app [GetDRAWxtlApp]
     
    11761217    }
    11771218    set coords [RB2cart [lindex [ReadRigidBody $rbtype] 1]]
    1178 puts $coords
    11791219    DRAWxtlPlotOrtho $file "" $coords $bondlist
    1180     if {$app != ""} [exec $app $file]
    1181 }
     1220    if {$app != ""} {exec $app $file &}
     1221}
     1222
     1223# plot orthogonal coordinates in DRAWxtl
     1224# input:
     1225#  coords: cartesian coordinates
     1226#  bondlist: list of bonds to draw as min, max length (A) and
     1227#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
     1228#  file: file name for the .str file to create
     1229proc PlotRBcoords {coords "bondlist {}" "file {}"} {
     1230    set app [GetDRAWxtlApp]
     1231    if {$app == ""} {
     1232        MyMessageBox -parent . -title "No DRAWxtl" \
     1233                -message "Sorry, DRAWxtl is not installed no phases are present to write" \
     1234                -icon warning
     1235        return
     1236    }
     1237    if {$::tcl_platform(platform) == "windows" && $file == ""} {
     1238        set file [file join [pwd] rbplot.str]
     1239    } else {
     1240        set file "/tmp/rbplot.str"
     1241    }
     1242    DRAWxtlPlotOrtho $file "" $coords $bondlist
     1243    if {$app != ""} {exec $app $file &}
     1244}
     1245
     1246# plot a rigid body superimposed on a structure
     1247# input:
     1248#  RBcoords: fractional coordinates for rigid body
     1249#  phase:# of phase to plot
     1250#  firstatom: # of 1st atom in structure matching rigid body
     1251#  allatoms: 0 to plot only atoms in phase that are in the rigid body,
     1252#      otherwise plot all atoms
     1253#  bondlist: list of bonds to draw for the phase as min, max length (A) and
     1254#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
     1255#  rbbondlist: list of bonds to draw for the phase as min, max length (A) and
     1256#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
     1257#  file: optional file name for the .str file to create
     1258proc DRAWxtlPlotRBFit {RBcoords phase firstatom "allatoms 0" \
     1259                           "bondlist {}" "rbbondlist {}" "file {}"} {
     1260    set natom [llength $RBcoords]
     1261    set app [GetDRAWxtlApp]
     1262    if {$app == ""} {
     1263        MyMessageBox -parent . -title "No DRAWxtl" \
     1264                -message "Sorry, DRAWxtl is not installed no phases are present to write" \
     1265                -icon warning
     1266        return
     1267    }
     1268    if {$::tcl_platform(platform) == "windows" && $file == ""} {
     1269        set file [file join [pwd] rbplot.str]
     1270    } else {
     1271        set file "/tmp/rbfit.str"
     1272    }
     1273
     1274    # get rigid body coordinate range
     1275    foreach {xmin ymin zmin} {"" "" ""} {}
     1276    foreach {xmax ymax zmax} {"" "" ""} {}
     1277    foreach xyz $RBcoords {
     1278        foreach {x y z} $xyz {}
     1279        foreach s {x y z} {
     1280            foreach t {min max} {
     1281                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
     1282            }
     1283            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
     1284            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
     1285        }
     1286    }
     1287    set rbrange {}
     1288    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
     1289                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
     1290                     [expr -0.01+$zmin] [expr 0.01+$zmax] ] {
     1291        append rbrange [format " %8.4f" $val]
     1292    }
     1293    set rbcenter [list [expr {($xmin+$xmax)/2}] \
     1294                      [expr {($ymin+$ymax)/2}] \
     1295                      [expr {($zmin+$zmax)/2}] ]
     1296    # get matching atoms coordinate range
     1297    foreach atom [lrange \
     1298                      [lrange $::expmap(atomlist_$phase) [expr {$firstatom-1}] end] \
     1299                      0 [expr {$natom-1}]] {
     1300        foreach s {x y z} {
     1301            set $s [atominfo $phase $atom $s]
     1302            foreach t {min max} {
     1303                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
     1304            }
     1305            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
     1306            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
     1307        }
     1308    }
     1309    # expand to cover at least one unit cell
     1310    foreach var {xmin ymin zmin} {
     1311        if {[set $var] > 0.0} {set $var 0.0}
     1312    }
     1313    foreach var {xmax ymax zmax} {
     1314        if {[set $var] < 1.} {set $var 1.}
     1315    }
     1316    set range {}
     1317    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
     1318                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
     1319                     [expr -0.01+$zmin] [expr 0.01+$zmax]] {
     1320        append range [format " %8.4f" $val]
     1321    }
     1322
     1323    set fp [open $file w]
     1324    puts $fp "title structure/rigid-body fit plot"
     1325    # plot the structure
     1326    puts -nonewline $fp "cell"
     1327    foreach p {a b c alpha beta gamma} {
     1328        puts -nonewline $fp " [phaseinfo $phase $p]"
     1329    }
     1330    puts $fp ""
     1331    puts $fp "spgp [phaseinfo $phase spacegroup]"
     1332    puts $fp "pack $range"
     1333    if {$allatoms != 0} {
     1334        set atoms $::expmap(atomlist_$phase)
     1335    } else {
     1336        set atoms [lrange \
     1337                       [lrange $::expmap(atomlist_$phase) [expr {$firstatom-1}] end] \
     1338                       0 [expr {$natom-1}]]
     1339    }
     1340
     1341    # set origin at center of rigid body
     1342    puts $fp "origin  $rbcenter"
     1343    # now loop over atoms
     1344    foreach atom $atoms {
     1345        set type [atominfo $phase $atom type]
     1346        set typelist($type) 1
     1347        puts -nonewline $fp "atom $type $atom "
     1348        foreach v {x y z} {
     1349            puts -nonewline $fp "[atominfo $phase $atom $v] "
     1350        }
     1351        puts $fp ""
     1352               
     1353        set uiso [atominfo $phase $atom Uiso]
     1354        # are there anisotropic atoms? If so convert them to Uequiv
     1355        if {[atominfo $phase $atom temptype] == "A"} {
     1356            puts -nonewline $fp "Uij [atominfo $phase $atom type] $atom "
     1357            foreach v {U11 U22 U33 U12 U13 U23} {
     1358                puts -nonewline $fp "[atominfo $phase $atom $v] "
     1359            }
     1360            puts $fp ""
     1361        }
     1362    }
     1363    foreach type [array names typelist] color {Green Blue Magenta Cyan} {
     1364        if {$type == ""} break
     1365        puts $fp "sphere $type 0.1 $color"
     1366    }
     1367    foreach type [array names typelist] color1 {Green Blue Magenta Cyan} {
     1368        foreach bondpair $bondlist {
     1369            foreach {b1 b2 color} $bondpair {}
     1370            if {$color == ""} {set color $color1}
     1371            puts $fp "bond $type $type 0.02 $b1 $b2 $color"
     1372        }
     1373        foreach type1 [array names typelist] {
     1374            if {$type1 == $type} break
     1375            foreach bondpair $bondlist {
     1376                foreach {b1 b2 color} $bondpair {}
     1377                if {$color == ""} {set color $color1}
     1378                puts $fp "bond $type $type1 0.02 $b1 $b2 $color"
     1379            }
     1380        }
     1381    }
     1382    # plot the rigid body
     1383    puts $fp "frame"
     1384    puts -nonewline $fp "cell"
     1385    foreach p {a b c alpha beta gamma} {
     1386        puts -nonewline $fp " [phaseinfo $phase $p]"
     1387    }
     1388    puts $fp ""
     1389    puts $fp "background White"
     1390    puts $fp "nolabels"
     1391    puts $fp "spgr P 1"
     1392    puts $fp "pack $rbrange"
     1393    set i 0
     1394    foreach xyz $RBcoords {
     1395        foreach {x y z} $xyz {}
     1396        incr i
     1397        puts $fp "atom c $i $x $y $z"
     1398    }
     1399    foreach bondpair $rbbondlist {
     1400        foreach {b1 b2 color} $bondpair {}
     1401        if {$color == ""} {set color Red}
     1402        puts $fp "bond c c 0.02 $b1 $b2 $color"
     1403    }
     1404
     1405    puts $fp "sphere c 0.05 Red"
     1406    puts $fp "finish   0.70   0.30   0.08   0.01"
     1407    puts $fp "end"
     1408
     1409    #puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
     1410    close $fp
     1411    if {$app != ""} {exec $app $file &}
     1412}
     1413
    11821414
    11831415#AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
     
    12281460# test code, plots rigid body type #2 with bonds drawn at ~1.3 & 2 A
    12291461#PlotRBtype 2 {{1.9 2.1} {1.28 1.32}}
     1462
     1463# test code, plots rigid body coords in ortholist with bonds @  ~2, 2.8 and 3.4 A
     1464#PlotRBcoords $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
     1465
    12301466
    12311467set useflag {1 1 1 1 1 1 1 1}
Note: See TracChangeset for help on using the changeset viewer.