Changeset 1107

Timestamp:

Jan 17, 2011 11:05:57 AM (13 years ago)

Author:

toby

Message:

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

File:

• branches/sandbox/rb.tcl (modified) (12 diffs)

branches/sandbox/rb.tcl

@@ -40 +40 @@
 #  not used in GSAS.
 proc ReadRigidBody {rbnum} {
-    if {[RigidBodyCount] > $rbnum} {
+    if {[RigidBodyCount] < $rbnum} {
         return ""
     }
@@ -49 +49 @@
     set used [string trim [string range [readexp "$key NBDS"] 0 4]]
     set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
-    set out $used
+    set out {}
     for {set i 1} {$i <= $nmult} {incr i} {
         set line [readexp "${key}${i}PARM"]
@@ -64 +64 @@
             set z [string trim [string range $line 20 29]]
             set lbl [string trim [string range $line 30 39]]
-            lappend coordlist [list $x $y $x $lbl]
+            lappend coordlist [list $x $y $z $lbl]
         }
         lappend out [list $mult $damp $var $coordlist]
     }
-    return $out
+    return [list $used $out]
 }
 
@@ -380 +380 @@
     setexp "$key NBDS" $value 1 5
     # number of coordinate matrices
-    set value [llength multlist]
+    set value [llength $multlist]
     validint value 5
     makeexprec "$key NSMP"
@@ -396 +396 @@
             validint value 3
             makeexprec "${key}${i}SC$value"
-            setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
+            if {[llength $item] == 4} {
+                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
+            } elseif {[llength $item] == 3} {
+                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f" $item] 1 30
+            } else {
+                return -code 3 "Invalid number of coordinates"
+            }
         }
     }
@@ -800 +806 @@
         set dvs 0
         foreach var {x y z} v1 $frac v2 $coord abc [lrange $cell 0 2] {
-            set dv [expr {($v2 - $v1)}]
+            set dv [expr {($v2 - $v1)*$abc}]
             set dvs [expr {$dvs + $dv*$dv}]
             set sumdvs [expr {$sumdvs + $dv*$dv}]
@@ -904 +910 @@
             set rmsprev $rms
         }
-    }   
+    } 
     #proc FitBodyOrigin {Euler cell ortholist useflag fraclist origin} 
     #return "$neworig $rms $fracout $rmsout"
@@ -1041 +1047 @@
 }
  
-
+# convert the rigid body representation reported as the 2nd element in ReadRigidBody
+# into cartesian coordinates
+#   rblist: a list containing an element for each scaling factor
+# in each element there are four items: 
+#    a multiplier value for the rigid body coordinates
+#    a damping value (0-9) for the refinement of the multiplier (not used)
+#    a variable number if the multiplier will be refined (not used)
+#    a list of cartesian coordinates coordinates 
+# each cartesian coordinate contains 4 items: x,y,z and a label
+# returns a list of coordinate triplets
+proc RB2cart {rblist} {
+    foreach item $rblist {
+        foreach {mult damp ref coords} $item {}
+        set i 0
+        foreach xyz $coords {
+            foreach {x y z} [lrange $xyz 0 2] {}
+            foreach val [lrange $xyz 0 2] var {X Y Z} {
+                if {[array names $var $i] == ""} {
+                    set ${var}($i) [expr {$mult * $val}]
+                } else {
+                    set ${var}($i) [expr {[set ${var}($i)] + $mult * $val}]
+                }
+            }
+            incr i
+        }
+    }
+    set out ""
+    foreach i [lsort -integer [array names X]] {
+        lappend out [list $X($i) $Y($i) $Z($i)]
+    }
+    return $out
+}
+
+proc GetDRAWxtlApp {} {
+    # is DRAWxtl installed?
+    set app {}
+    if {![catch {set fp [open [file join $::env(HOME) .drawxtlrc] r]}]} {
+        # line 12 is name of executable
+        set i 0
+        while {$i < 12} {
+            incr i
+            gets $fp appname
+        }
+        close $fp
+        set app [auto_execok $appname]
+    }
+    if {$app != ""} {
+        return $appname
+    }
+    return ""
+}
+
+proc DRAWxtlPlotOrtho {filename title coords bondlist} {
+    foreach {xmin ymin zmin} {"" "" ""} {}
+    foreach {xmax ymax zmax} {"" "" ""} {}
+    foreach xyz $coords {
+        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]}
+        }
+    }
+    puts "$xmin $xmax $ymin $ymax $zmin $zmax"
+    set max $xmin
+    foreach val "$xmin $xmax $ymin $ymax $zmin $zmax" {
+        if {$max < abs($val)} {set max $val}
+    }
+    set scale [expr {4.*$max}]
+    set a 10.
+    lappend range [expr -0.01+($xmin/$scale)] [expr 0.01+($xmax/$scale)] \
+        [expr -0.01+($ymin/$scale)] [expr 0.01+($ymax/$scale)] \
+        [expr -0.01+($zmin/$scale)] [expr 0.01+($zmax/$scale)]
+    set fp [open $filename w]
+    puts $fp "title $fp"
+    puts $fp "box  0.000 Black"
+    puts $fp "background White"
+    puts $fp "nolabels"
+    puts $fp "cell $a $a $a 90 90 90"
+    puts $fp "spgr P 1"
+    puts $fp "pack $range"
+    set i 0
+    foreach xyz $coords {
+        foreach {x y z} $xyz {}
+        incr i
+        puts $fp "atom c $i [expr {$x/$scale}] [expr {$y/$scale}] [expr {$z/$scale}]"
+    }
+    puts $fp "sphere c  [expr 0.100*($a/$scale)] Red"
+    puts $fp "finish   0.70   0.30   0.08   0.01"
+    foreach bondpair $bondlist {
+        foreach {b1 b2 color} $bondpair {}
+        if {$color == ""} {set color Red}
+        puts $fp "bond c c [expr {0.01*$a/$scale}] [expr {$b1*$a/$scale}] [expr {$b2*$a/$scale}] $color"
+    }
+    puts $fp "frame"
+    set range {}
+    lappend range -0.01 [expr 0.01+(0.1*$a/$scale)] \
+        -0.01 [expr 0.01+(0.1*$a/$scale)] \
+        -0.01 [expr 0.01+(0.1*$a/$scale)]
+    puts $fp "cell $a $a $a 90 90 90"
+    puts $fp "spgr P 1"
+    puts $fp "pack $range"
+    puts $fp "atom o 1 0 0 0"
+    puts $fp "atom o 2 [expr {0.1*$a/$scale}] 0 0"
+    puts $fp "atom o 3 0 [expr {0.1*$a/$scale}] 0"
+    puts $fp "atom o 4 0 0 [expr {0.1*$a/$scale}]"
+    puts $fp "bond o o [expr {0.01*$a/$scale}] [expr {-0.1 + $a/$scale}] [expr {0.1 + $a/$scale}] Black"
+    puts $fp "sphere o [expr {0.02*$a/$scale}] Blue"
+    puts $fp "origin   .0 .0 .0"
+    puts $fp "end"
+    close $fp
+}
+
+proc PlotRBtype {rbtype "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"
+    }
+    set coords [RB2cart [lindex [ReadRigidBody $rbtype] 1]]
+puts $coords
+    DRAWxtlPlotOrtho $file "" $coords $bondlist
+    if {$app != ""} [exec $app $file]
+}
 
 #AddRigidBody {1} { {{0 0 0 xe} {1 1 1 o} {2 2 2 si+4}} }
@@ -1063 +1201 @@
 #RunRecalcRBCoords
 
-puts "Test FitBody"
+#puts "Test FitBody"
 set fraclist {
     { 0.5483305238484277 0.4887545024531055 0.6167996784631056 }
@@ -1084 +1222 @@
     {     -1.000000   -1.000000   -1.000000} 
 }
+# test code, generates DRAWxtl imput file from orthogonal coordinate list
+# with bonds of ~2, 2.8 and 3.4 A
+#DRAWxtlPlotOrtho test4.str "test file" $ortholist {{1.9 2.1} {3.4 3.5 Blue} {2.8 2.83 Green} }
+
+# 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}}
 
 set useflag {1 1 1 1 1 1 1 1}
@@ -1094 +1238 @@
 
 #puts [La::show $xform]
-puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
+#puts "out: [FitBody $Euler $cell $ortholist $useflag $fraclist $origin 30]"
 
 
@@ -1114 +1258 @@
 #  H       -1.38570       -1.36644        0.89518
 #  H       -1.38570       -1.36644       -0.89518
-set coordlist [zmat2coord $atmlist]
-set i 0
-puts "\nZmatrix in"
-foreach line $atmlist {
-    incr i
-    puts "$i) $line"
-}
-puts "Cartesian out"
-foreach line $coordlist {
-    puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
-}
+# set coordlist [zmat2coord $atmlist]
+ set i 0
+# puts "\nZmatrix in"
+# foreach line $atmlist {
+#     incr i
+#     puts "$i) $line"
+# }
+# puts "Cartesian out"
+# foreach line $coordlist {
+#     puts [eval format "%-4s%10.5f%10.5f%10.5f" $line]
+# }
+
+# AddRigidBody {1 0.75} { 
+#     {
+#       {1 1 1 c} 
+#       {-1 1 1 c}
+#       {      1.000000   -1.000000    1.000000 c}
+#       {     -1.000000   -1.000000    1.000000 c}
+#       {      1.000000    1.000000   -1.000000 c}
+#       {     -1.000000    1.000000   -1.000000 c}
+#       {      1.000000   -1.000000   -1.000000 c}
+#       {     -1.000000   -1.000000   -1.000000 c} 
+#       {1 1 1 h} 
+#       {1 -1 -1 h} 
+#       {-1 1 -1 h} 
+#       {-1 -1 1 h} 
+#     } {
+#       {0 0 0 c }
+#       {0 0 0 c}
+#       {0 0 0 c}
+#       {0 0 0 c}
+#       {0 0 0 c}
+#       {0 0 0 c}
+#       {0 0 0 c}
+#       {0 0 0 c}
+#       {1 1 1 h} 
+#       {1 -1 -1 h} 
+#       {-1 1 -1 h} 
+#       {-1 -1 1 h} 
+#     }
+# }
+# MapRigidBody 2 2 1 {0 0 0} {10 15 20}