Changeset 1108

Timestamp:

Jan 18, 2011 3:43:29 PM (13 years ago)

Author:

toby

Message:

relabel xform atoms; more RB stuff

Location:

branches/sandbox

Files:

• expgui (modified) (3 diffs)
• rb.tcl (modified) (15 diffs)

branches/sandbox/expgui

@@ -996 +996 @@
 
     set crsPhase {}
-    $expgui(atomxform) config -text "Xform Atoms" -state disabled
+    $expgui(atomxform) config -text "Modify Atoms" -state disabled
     foreach n $expmap(phaselist) type $expmap(phasetype) {
         if {$n == $num} {
@@ -1382 +1382 @@
     DisplayRefFlags $atomnum $p
     $expgui(EditingAtoms) config -text "Editing atom #$atomnum -- [atominfo $p $atomnum label]"
-    $expgui(atomxform) config -text "Xform Atom" -state normal
+    $expgui(atomxform) config -text "Modify Atom" -state normal
 }
 
@@ -1407 +1407 @@
     # this needs to track by phase
     DisplayRefFlags $numberList $p
-    $expgui(atomxform) config -text "Xform Atoms" -state normal
+    $expgui(atomxform) config -text "Modify Atoms" -state normal
 }
 

branches/sandbox/rb.tcl

@@ -18 +18 @@
 # MapRigidBody -- map a rigid body type into a phase
 # EditRigidBodyMapping -- change the parameters in a rigid body mapping
+#  need to unmap a rigid body
+#  delete rigid body
 #============================================================================
 # returns the number of defined rigid bodies 
@@ -74 +76 @@
 proc RigidBodyMappingCount {phase bodytyp} {
     set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
+    if {! [existsexp "$key  NBDS"]} {return 0}
     set n [string trim [readexp "$key  NBDS"]]
     if {$n == ""} {
@@ -773 +776 @@
 }
 
+#calc rigid body fractional coordinates
+# arguments: 
+#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...} 
+#              where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
+#  cell: a list with "a b c alpha beta gamma" in Angstroms and degrees
+#  ortholist: list containing triplets with orthogonal coordinates 
+#  origin: a list of fraction coordinates {x y z} describing the location of the 
+#            origin of the orthogonal coordinates in the crystal system
+#     note that the length of ortholist, useflag and fraclist should be the same
+# Returns a list with the computed fractional coordinates for all atoms
+proc CalcBody {Euler cell ortholist origin} {
+    set xform [CalcXformMatrix $Euler $cell]
+    set i 0
+    set sumdvs 0
+    set fracout {}
+    set rmsout {}
+    foreach oc $ortholist {
+        set frac [lrange [Ortho2Xtal $xform $origin $oc] 3 end]
+        lappend fracout $frac
+    }
+    return $fracout
+}
+
+
 # fit a rigid body's origin 
 # arguments: 
@@ -809 +836 @@
             set dvs [expr {$dvs + $dv*$dv}]
             set sumdvs [expr {$sumdvs + $dv*$dv}]
-            if {$use} {set sum($var) [expr {$sum($var) + $dv}]}
+            if {$use} {set sum($var) [expr {$sum($var) + $dv/$abc}]}
         }
         lappend rmsout [expr {sqrt($dvs)}]
@@ -887 +914 @@
 #   rms deviation in fractional coordinates of new Euler angles 
 proc FitBody {Euler cell ortholist useflag fraclist origin "ncycle 5"} {
-    puts "start origin = $origin"
+    #puts "start origin = $origin"
     foreach {
         origin 
@@ -893 +920 @@
         fracout 
         rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
-    puts "start rms = $startrms"
+    #puts "start rms = $startrms"
     set rmsprev $startrms
-    puts "new origin = $origin"
+    #puts "new origin = $origin"
     for {set i 0} {$i < $ncycle} {incr i} {
         set Eulerprev $Euler
         set Euler [FitBodyRot $Euler $cell $ortholist $useflag $fraclist $origin]
-        puts "New Euler $Euler"
-        puts "after fit" 
+        #puts "New Euler $Euler"
+        #puts "after fit" 
         foreach {
             origin 
@@ -907 +934 @@
             rmsout } [FitBodyOrigin $Euler $cell $ortholist $useflag $fraclist $origin] {}
         if {$rms > (1.1 * $rmsprev) + 0.01} {
-            puts "rms = $rms, new origin = $origin"
+            #puts "rms = $rms, new origin = $origin"
             set rmsprev $rms
         }
@@ -914 +941 @@
     #return "$neworig $rms $fracout $rmsout"
     set fmt  {"%8.5f %8.5f %8.5f     %8.5f %8.5f %8.5f   %6.3f"}
-    foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
-        puts "[eval format $fmt $fracin $fraccalc $rmsi]"
-    }
+    #foreach fracin $fraclist fraccalc $fracout rmsi $rmsout {
+        #puts "[eval format $fmt $fracin $fraccalc $rmsi]"
+    #}
     return [list $origin $Euler $rms $fracout $rmsout]
 }
@@ -1080 +1107 @@
 }
 
+# get the name of the DRAWxtl application, if installed
 proc GetDRAWxtlApp {} {
     # is DRAWxtl installed?
@@ -1099 +1127 @@
 }
 
+# plot orthogonal coordinates in DRAWxtl
+# input:
+#  filename: file name for the .str file to create
+#  title: string for title in .str file
+#  coords: cartesian coordinates
+#  bondlist: list of bonds to draw as min, max length (A) and 
+#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
 proc DRAWxtlPlotOrtho {filename title coords bondlist} {
     foreach {xmin ymin zmin} {"" "" ""} {}
@@ -1112 +1147 @@
         }
     }
-    puts "$xmin $xmax $ymin $ymax $zmin $zmax"
+    #puts "$xmin $xmax $ymin $ymax $zmin $zmax"
     set max $xmin
     foreach val "$xmin $xmax $ymin $ymax $zmin $zmax" {
@@ -1123 +1158 @@
         [expr -0.01+($zmin/$scale)] [expr 0.01+($zmax/$scale)]
     set fp [open $filename w]
-    puts $fp "title $fp"
+    puts $fp "title $title"
     puts $fp "box  0.000 Black"
     puts $fp "background White"
@@ -1162 +1197 @@
 }
 
+# plot a rigid body in DRAWxtl
+# input:
+#  rbtype: # of rigid body
+#  bondlist: list of bonds to draw as min, max length (A) and 
+#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
+#  file: file name for the .str file to create
 proc PlotRBtype {rbtype "bondlist {}" "file {}"} {
     set app [GetDRAWxtlApp]
@@ -1176 +1217 @@
     }
     set coords [RB2cart [lindex [ReadRigidBody $rbtype] 1]]
-puts $coords
     DRAWxtlPlotOrtho $file "" $coords $bondlist
-    if {$app != ""} [exec $app $file]
-}
+    if {$app != ""} {exec $app $file &}
+}
+
+# plot orthogonal coordinates in DRAWxtl
+# input:
+#  coords: cartesian coordinates
+#  bondlist: list of bonds to draw as min, max length (A) and 
+#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
+#  file: file name for the .str file to create
+proc PlotRBcoords {coords "bondlist {}" "file {}"} {
+    set app [GetDRAWxtlApp]
+    if {$app == ""} {
+        MyMessageBox -parent . -title "No DRAWxtl" \
+                -message "Sorry, DRAWxtl is not installed no phases are present to write" \
+                -icon warning
+        return
+    }
+    if {$::tcl_platform(platform) == "windows" && $file == ""} {
+        set file [file join [pwd] rbplot.str]
+    } else {
+        set file "/tmp/rbplot.str"
+    }
+    DRAWxtlPlotOrtho $file "" $coords $bondlist
+    if {$app != ""} {exec $app $file &}
+}
+
+# plot a rigid body superimposed on a structure
+# input:
+#  RBcoords: fractional coordinates for rigid body
+#  phase:# of phase to plot
+#  firstatom: # of 1st atom in structure matching rigid body
+#  allatoms: 0 to plot only atoms in phase that are in the rigid body, 
+#      otherwise plot all atoms
+#  bondlist: list of bonds to draw for the phase as min, max length (A) and 
+#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
+#  rbbondlist: list of bonds to draw for the phase as min, max length (A) and 
+#      an optional color; for example: {{1.4 1.6} {1.2 1.3 Red}}
+#  file: optional file name for the .str file to create
+proc DRAWxtlPlotRBFit {RBcoords phase firstatom "allatoms 0" \
+                           "bondlist {}" "rbbondlist {}" "file {}"} {
+    set natom [llength $RBcoords]
+    set app [GetDRAWxtlApp]
+    if {$app == ""} {
+        MyMessageBox -parent . -title "No DRAWxtl" \
+                -message "Sorry, DRAWxtl is not installed no phases are present to write" \
+                -icon warning
+        return
+    }
+    if {$::tcl_platform(platform) == "windows" && $file == ""} {
+        set file [file join [pwd] rbplot.str]
+    } else {
+        set file "/tmp/rbfit.str"
+    }
+
+    # get rigid body coordinate range
+    foreach {xmin ymin zmin} {"" "" ""} {}
+    foreach {xmax ymax zmax} {"" "" ""} {}
+    foreach xyz $RBcoords {
+        foreach {x y z} $xyz {}
+        foreach s {x y z} {
+            foreach t {min max} {
+                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
+            } 
+            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
+            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
+        }
+    }
+    set rbrange {}
+    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
+                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
+                     [expr -0.01+$zmin] [expr 0.01+$zmax] ] {
+        append rbrange [format " %8.4f" $val]
+    }
+    set rbcenter [list [expr {($xmin+$xmax)/2}] \
+                      [expr {($ymin+$ymax)/2}] \
+                      [expr {($zmin+$zmax)/2}] ]
+    # get matching atoms coordinate range
+    foreach atom [lrange \
+                      [lrange $::expmap(atomlist_$phase) [expr {$firstatom-1}] end] \
+                      0 [expr {$natom-1}]] {
+        foreach s {x y z} {
+            set $s [atominfo $phase $atom $s]
+            foreach t {min max} {
+                if {[set ${s}${t}] == ""} {set ${s}${t} [set $s]}
+            } 
+            if {[set ${s}min] > [set $s]} {set ${s}min [set $s]}
+            if {[set ${s}max] < [set $s]} {set ${s}max [set $s]}
+        }
+    }
+    # expand to cover at least one unit cell
+    foreach var {xmin ymin zmin} { 
+        if {[set $var] > 0.0} {set $var 0.0}
+    }
+    foreach var {xmax ymax zmax} { 
+        if {[set $var] < 1.} {set $var 1.}
+    }
+    set range {}
+    foreach val [list [expr -0.01+$xmin] [expr 0.01+$xmax] \
+                     [expr -0.01+$ymin] [expr 0.01+$ymax] \
+                     [expr -0.01+$zmin] [expr 0.01+$zmax]] {
+        append range [format " %8.4f" $val]
+    }
+
+    set fp [open $file w]
+    puts $fp "title structure/rigid-body fit plot"
+    # plot the structure
+    puts -nonewline $fp "cell"
+    foreach p {a b c alpha beta gamma} {
+        puts -nonewline $fp " [phaseinfo $phase $p]"
+    }
+    puts $fp ""
+    puts $fp "spgp [phaseinfo $phase spacegroup]"
+    puts $fp "pack $range"
+    if {$allatoms != 0} {
+        set atoms $::expmap(atomlist_$phase)
+    } else {
+        set atoms [lrange \
+                       [lrange $::expmap(atomlist_$phase) [expr {$firstatom-1}] end] \
+                       0 [expr {$natom-1}]]
+    }
+
+    # set origin at center of rigid body
+    puts $fp "origin  $rbcenter"
+    # now loop over atoms
+    foreach atom $atoms {
+        set type [atominfo $phase $atom type]
+        set typelist($type) 1
+        puts -nonewline $fp "atom $type $atom "
+        foreach v {x y z} {
+            puts -nonewline $fp "[atominfo $phase $atom $v] "
+        }
+        puts $fp ""
+                
+        set uiso [atominfo $phase $atom Uiso]
+        # are there anisotropic atoms? If so convert them to Uequiv
+        if {[atominfo $phase $atom temptype] == "A"} {
+            puts -nonewline $fp "Uij [atominfo $phase $atom type] $atom "
+            foreach v {U11 U22 U33 U12 U13 U23} {
+                puts -nonewline $fp "[atominfo $phase $atom $v] "
+            }
+            puts $fp ""
+        }
+    }
+    foreach type [array names typelist] color {Green Blue Magenta Cyan} {
+        if {$type == ""} break
+        puts $fp "sphere $type 0.1 $color"
+    }
+    foreach type [array names typelist] color1 {Green Blue Magenta Cyan} {
+        foreach bondpair $bondlist {
+            foreach {b1 b2 color} $bondpair {}
+            if {$color == ""} {set color $color1}
+            puts $fp "bond $type $type 0.02 $b1 $b2 $color"
+        }
+        foreach type1 [array names typelist] {
+            if {$type1 == $type} break
+            foreach bondpair $bondlist {
+                foreach {b1 b2 color} $bondpair {}
+                if {$color == ""} {set color $color1}
+                puts $fp "bond $type $type1 0.02 $b1 $b2 $color"
+            }
+        }
+    }
+    # plot the rigid body
+    puts $fp "frame"
+    puts -nonewline $fp "cell"
+    foreach p {a b c alpha beta gamma} {
+        puts -nonewline $fp " [phaseinfo $phase $p]"
+    }
+    puts $fp ""
+    puts $fp "background White"
+    puts $fp "nolabels"
+    puts $fp "spgr P 1"
+    puts $fp "pack $rbrange"
+    set i 0
+    foreach xyz $RBcoords {
+        foreach {x y z} $xyz {}
+        incr i
+        puts $fp "atom c $i $x $y $z"
+    }
+    foreach bondpair $rbbondlist {
+        foreach {b1 b2 color} $bondpair {}
+        if {$color == ""} {set color Red}
+        puts $fp "bond c c 0.02 $b1 $b2 $color"
+    }
+
+    puts $fp "sphere c 0.05 Red"
+    puts $fp "finish   0.70   0.30   0.08   0.01"
+    puts $fp "end"
+
+    #puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
+    close $fp
+    if {$app != ""} {exec $app $file &}
+}
+
 
 #AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
@@ -1228 +1460 @@
 # test code, plots rigid body type #2 with bonds drawn at ~1.3 & 2 A
 #PlotRBtype 2 {{1.9 2.1} {1.28 1.32}}
+
+# test code, plots rigid body coords in ortholist with bonds @  ~2, 2.8 and 3.4 A
+#PlotRBcoords $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
+
 
 set useflag {1 1 1 1 1 1 1 1}