Changeset 1166 for trunk/readexp.tcl


Ignore:
Timestamp:
Aug 17, 2011 6:17:04 PM (9 years ago)
Author:
toby
Message:

bring sandbox changes over to main release

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/readexp.tcl

    r1033 r1166  
    99# returns 1 if the file is 80 char/line + cr/lf
    1010# returns 2 if the file is sequential but not fixed-record length
    11 proc expload {expfile} {
    12     global exparray tcl_platform
    13     # $expfile is the path to the data file.
    14 
     11proc expload {expfile "ns {}"} {
     12    # expfile is the path to the data file.
     13    # ns is the namespace to place the output array (default is global)
     14    if {$ns != ""} {
     15        namespace eval $ns {}
     16    }
    1517    if [catch {set fil [open "$expfile" r]}] {
    1618        tk_dialog .expFileErrorMsg "File Open Error" \
     
    2729    }
    2830    catch {
    29         unset exparray
     31        unset ${ns}::exparray
    3032    }
    3133    if {$len > 160} {
     
    3941            incr i2 80
    4042            set key [string range $nline 0 11]
    41             set exparray($key) [string range $nline 12 end]
     43            set ${ns}::exparray($key) [string range $nline 12 end]
    4244        }
    4345    } else {
     
    4547        while {$len > 0} {
    4648            set key [string range $line 0 11]
    47             set exparray($key) [string range $line 12 79]
     49            set ${ns}::exparray($key) [string range $line 12 79]
    4850            if {$len != 81 || [string range $line end end] != "\r"} {set fmt 2}
    4951            set len [gets $fil line]
     
    187189        }
    188190    }
     191    # load the constrained parameters
     192    atom_constraint_load
    189193    set expgui(mapstat) 1
    190194}
     
    294298            set digits 1
    295299        }
     300        if {$digits <= 0} {set digits 1}
    296301        if {$digits + $sign >= $length} {
    297302            # the number is much too big -- use exponential notation
     
    307312            set decimal [expr {$length - $digits - $sign - 1}]
    308313            set tmp [format "%#.${decimal}f" $value]
    309         } elseif {abs($value) < pow(10,2-$decimal) && $length > 6} {
     314        } elseif {abs($value) < pow(10,2-$decimal) && $length > 6 && $value != 0} {
    310315            # for small values, switch to exponential notation (2-$decimal -> three sig figs.)
    311316            set decimal [expr {$length - 6 - $sign}]
     
    13881393#     absdamp -- damping for absorption refinement (*)
    13891394#     absref -- refinement damping for absorption refinement (*)
    1390 #     proftbl -- returns number of profile table termns
     1395#     proftbl -- returns number of profile table terms
     1396#     anomff -- returns a list of elements, f' and f"
    13911397#   parameters transferred from the instrument parameter file:
    13921398#     ITYP    -- returns the contents of the ITYP record
     
    19771983                return $val
    19781984            }
     1985            anomff-get {
     1986                set l {}
     1987                foreach i {1 2 3 4 5 6 7 8 9} {
     1988                    if {![existsexp "${key}FFANS$i"]} continue
     1989                    set line [readexp "${key}FFANS$i"]
     1990                    set elem [string trim [string range $line 2 9]]
     1991                    set fp [string trim [string range $line 10 19]]
     1992                    set fpp [string trim [string range $line 20 29]]
     1993                    lappend l [list $elem $fp $fpp]
     1994                }
     1995                return $l
     1996            }
     1997            anomff-set {
     1998                # match up input against elements in list.
     1999                # change elements included, return any elements that are
     2000                # not found.
     2001                set errorlist {}
     2002                foreach triplet $value {
     2003                    foreach {e fp fpp} $triplet {}               
     2004                    foreach i {1 2 3 4 5 6 7 8 9} {
     2005                        if {![existsexp "${key}FFANS$i"]} continue
     2006                        # note that the element name is not used or validated
     2007                        set elem [string trim [string range \
     2008                                                   [readexp "${key}FFANS$i"] 2 9]]
     2009                        if {[string match -nocase $e $elem]} {
     2010                            if ![validreal fp 10 3] {return 0}
     2011                            setexp "${key}FFANS$i" $fp 11 10
     2012                            if ![validreal fpp 10 3] {return 0}
     2013                            setexp "${key}FFANS$i" $fpp 21 10
     2014                            set e {}
     2015                            break
     2016                        }
     2017                    }
     2018                    if {$e != ""} {
     2019                        # oops, no match
     2020                        lappend errorlist $e
     2021                    }
     2022                }
     2023                if {$errorlist != ""} {return [list 0 $errorlist]}
     2024            }
    19792025            default {
    19802026                set msg "Unsupported histinfo access: parm=$parm action=$action"
     
    21732219    return 1
    21742220}
     2221
     2222#  read fixed constraints
     2223
     2224proc atom_constraint_read {phase} {
     2225    set fix_list ""
     2226    foreach k {1 2 3 4 5 6 7 8 9} {
     2227        set key [format "LEQV HOLD%1d%2d" $phase $k]
     2228        set line [readexp $key]
     2229        foreach j {2 10 18 26 34 42 50 58} {
     2230            set fix_param [string range $line $j [expr $j+7]]
     2231            if {[string trim $fix_param] == ""} {return $fix_list}
     2232            lappend fix_list $fix_param
     2233        }
     2234    }
     2235}
     2236
     2237# load all atom constraints into global array fix_param
     2238proc atom_constraint_load { } {
     2239    catch {unset ::fix_param}
     2240    foreach i $::expmap(phaselist) {
     2241        set temp [atom_constraint_read $i]
     2242        foreach j $temp {
     2243            set atomnum [string trim [string range $j 2 3]]
     2244            set param [string trim [string range $j 4 6]]
     2245            set ::fix_param($i,$atomnum,$param) 1   
     2246        }
     2247    }
     2248}
     2249
     2250proc atom_constraint_write {phase fix_list} {
     2251    foreach key [array names ::exparray "LEQV HOLD$phase*"] {
     2252        delexp $key
     2253    }
     2254    set k 0
     2255    set j 1
     2256    set line ""
     2257    foreach fix $fix_list {
     2258        incr k
     2259        append line $fix
     2260        if {$k == 8} {
     2261            set key [format "LEQV HOLD%1d%2d" $phase $j]
     2262            makeexprec $key
     2263            setexp $key $line 3 [expr ($k * 8) + 2]
     2264            set k 0
     2265            incr j
     2266            set line ""
     2267        }
     2268    }
     2269    if {$line != ""} {
     2270        set key [format "LEQV HOLD%1d%2d" $phase $j]
     2271        makeexprec $key
     2272        setexp $key $line 3 [expr ($k * 8) + 2]
     2273    }   
     2274}
     2275
    21752276
    21762277#  get a logical constraint
     
    30993200        }
    31003201    }
    3101     if {$HST == ""} {return "1"}
    31023202    if {$HST <=9} {
    31033203        set key "HST  $HST"
     
    31193219            makeexprec $key
    31203220        }
    3121         setexp $key "RSN" [expr 3 + 5*(($HST-1) % 12)] 5
     3221        setexp $key "RSN " [expr 3 + 5*(($HST-1) % 12)] 5
    31223222        # create other HST  xx recs
    31233223        if {$HST <=9} {
     
    31303230        makeexprec "$key FACTR"
    31313231        makeexprec "$key NBNDS"
     3232        mapexp
     3233    } elseif {$HST == ""} {
     3234        if $::expgui(debug) {puts "no restraints"}
     3235        return "1"
    31323236    }
    31333237
     
    34493553    }
    34503554}
     3555#============================================================================
     3556# rigid body EXP editing routines (to move into readexp.tcl)
     3557# RigidBodyList -- returns a list of the defined rigid body types
     3558# SetRigidBodyVar -- set variables and damping for rigid body type multipliers
     3559# ReadRigidBody  -- # of times a body is mapped, scaling factors, var #s & coordinates
     3560# RigidBodyMappingList - return a list instances where a RB is mapped in phase
     3561# RigidBodyEnableTLS -- Enable or Disable TLS use for a rigid body mapping
     3562# RigidBodySetTLS  -- change the TLS values for a rigid body mapping
     3563# RigidBodySetDamp -- change the damping values for a rigid body mapping
     3564# RigidBodyVary    -- set refinement variable numbers for a rigid body mapping
     3565# RigidBodyTLSVary -- set TLS refinement variable nums for a rigid body mapping
     3566# AddRigidBody -- defines a new rigid body type
     3567# DeleteRigidBody -- remove a rigid body definition
     3568# ReplaceRigidBody -- replaces a previous rigid body type
     3569# ReadRigidBodyMapping  -- get parameters for a rigid body mapping
     3570# MapRigidBody -- map a rigid body type into a phase
     3571# EditRigidBodyMapping -- change the parameters in a rigid body mapping
     3572# UnMapRigidBody --remove a rigid body constraint by removing a RB "instance"
     3573#----- note that these older routines should not be used ------
     3574# RigidBodyCount -- returns the number of defined rigid bodies (body types)
     3575#    use RigidBodyList instead
     3576# RigidBodyMappingCount -- # of times a rigid body is mapped in phase
     3577#    use RigidBodyMappingList instead
     3578#============================================================================
     3579# returns the number of defined rigid bodies
     3580proc RigidBodyCount {} {
     3581    set n [string trim [readexp "RGBD  NRBDS"]]
     3582    if {$n == ""} {
     3583        set n 0
     3584    }
     3585    return $n
     3586}
     3587
     3588# returns a list of the defined rigid body types
     3589proc RigidBodyList {} {
     3590    set n [string trim [readexp "RGBD  NRBDS"]]
     3591    if {$n == ""} {
     3592        set n 0
     3593    }
     3594    set rblist {}
     3595    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
     3596        set value $rbnum
     3597        validint value 2
     3598        set key "RGBD${value}"
     3599        if {[existsexp "$key NATR "]} {
     3600            lappend rblist $rbnum
     3601        }
     3602        if {[llength $rblist] == $n} break
     3603    }
     3604    return $rblist
     3605}
     3606
     3607# ReadRigidBody provides all information associated with a rigid body type
     3608#  rbnum is the rigid body type number
     3609# it returns two items:
     3610#   the number of times the rigid body is mapped
     3611#   a list containing an element for each scaling factor in rigid body #rbnum.
     3612# in each element there are four items:
     3613#    a multiplier value for the rigid body coordinates
     3614#    a damping value (0-9) for the refinement of the multiplier
     3615#    a variable number if the multiplier will be refined
     3616#    a list of cartesian coordinates coordinates
     3617# each cartesian coordinate contains 4 items: x,y,z and a label
     3618#  note that the label is present only when the RB is created in EXPGUI and is
     3619#  not used in GSAS.
     3620proc ReadRigidBody {rbnum} {
     3621    if {[lsearch [RigidBodyList] $rbnum] == -1} {
     3622        return ""
     3623    }
     3624    set value $rbnum
     3625    validint value 2
     3626    set key "RGBD${value}"
     3627    set n [string trim [string range [readexp "$key NATR"] 0 4]]
     3628    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
     3629    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
     3630    set out {}
     3631    for {set i 1} {$i <= $nmult} {incr i} {
     3632        set line [readexp "${key}${i}PARM"]
     3633        set mult [string trim [string range $line 0 9]]
     3634        set var [string trim [string range $line 10 14]]
     3635        set damp [string trim [string range $line 15 19]]
     3636        set coordlist {}
     3637        for {set j 1} {$j <= $n} {incr j} {
     3638            set value $j
     3639            validint value 3
     3640            set line [readexp "${key}${i}SC$value"]
     3641            set x [string trim [string range $line 0 9]]
     3642            set y [string trim [string range $line 10 19]]
     3643            set z [string trim [string range $line 20 29]]
     3644            set lbl [string trim [string range $line 30 39]]
     3645            lappend coordlist [list $x $y $z $lbl]
     3646        }
     3647        lappend out [list $mult $damp $var $coordlist]
     3648    }
     3649    return [list $used $out]
     3650}
     3651
     3652# SetRigidBodyVar
     3653#   rbnum is the rigid body type number
     3654#   varnumlist is a list of variable numbers
     3655#      note that if this list is shorter than the number of actual multipliers
     3656#      for the body, the unspecified variable will not be changed
     3657#   damplist   is a list of damping values (0-9)
     3658#      note that if the damplist is shorter than the number of actual multipliers
     3659#      the unspecified values are not changed
     3660#  SetRigidBodVar 2 {1 2 3} {}
     3661#       will vary the (first 3) translations in body #3 and will not change the
     3662#       damping values
     3663#  SetRigidBodVar 3 {} {0 0 0}
     3664#       will not change variable settings but will change the (first 3) damping values
     3665#  SetRigidBodVar 4 {11 11} {8 8}
     3666#      changes both variable numbers and damping at the same time
     3667# Nothing is returned
     3668proc SetRigidBodyVar {rbnum varnumlist damplist} {
     3669    if {[lsearch [RigidBodyList] $rbnum] == -1} {
     3670        return ""
     3671    }
     3672    set value $rbnum
     3673    validint value 2
     3674    set key "RGBD${value}"
     3675    set nmult [string trim [string range [readexp "$key NSMP"] 0 4]]
     3676    for {set i 1} {$i <= $nmult} {incr i} {
     3677        set j $i
     3678        incr j -1
     3679        set var [lindex $varnumlist $j]
     3680        if {$var != ""} {
     3681            validint var 5
     3682            setexp "${key}${i}PARM" $var 11 15
     3683        }
     3684        set damp [lindex $damplist $j]
     3685        if {$damp != ""} {
     3686            if {$damp > 9} {set damp 9}
     3687            if {$damp < 0} {set damp 0}
     3688            validint damp 5
     3689        }
     3690        setexp "${key}${i}PARM" $damp 16 20
     3691    }
     3692}
     3693
     3694
     3695# return the number of times rigid body $bodytyp is mapped in phase $phase
     3696proc RigidBodyMappingCount {phase bodytyp} {
     3697    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
     3698    if {! [existsexp "$key  NBDS"]} {return 0}
     3699    set n [string trim [readexp "$key  NBDS"]]
     3700    if {$n == ""} {
     3701        set n 0
     3702    }
     3703    return $n
     3704}
     3705# return a list of the instances where rigid body $bodytyp is mapped in phase $phase
     3706proc RigidBodyMappingList {phase bodytyp} {
     3707    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
     3708    if {! [existsexp "$key  NBDS"]} {return {}}
     3709    set n [string trim [readexp "$key  NBDS"]]
     3710    if {$n == ""} {
     3711        set n 0
     3712    }
     3713    set rblist {}
     3714    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
     3715        set value $rbnum
     3716        validint value 2
     3717        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
     3718        if {[existsexp "$key  NDA"]} {
     3719            lappend rblist $rbnum
     3720        }
     3721        if {[llength $rblist] == $n} break
     3722    }
     3723    return $rblist
     3724}
     3725
     3726
     3727
     3728# reads rigid body mapping parameters for phase ($phase), body type # ($bodytyp) and instance # ($num)
     3729# returns a list of items (most lists) as follows:
     3730#   1) sequence # of first atom in body
     3731#   2) origin of body in fractional coordinates (3 elements)
     3732#   3) Euler angles as 6 pairs of numbers (see below)
     3733#   4) variable numbers for the 9 position variables (origin followed by rotations)
     3734#   5) damping vals for the 9 position variables (origin followed by rotations)
     3735#   6) the TLS values, in order below (empty list if TLS is not in use)
     3736#   7) the variable numbers for each TLS values, in order below (or empty)
     3737#   8) three damping values for the T, L and S terms.
     3738# returns an empty list if no such body exists.
     3739#
     3740# Euler angles are a list of axes and angles to rotate:
     3741#   { {axis1 angle1} {axis2 angle2} ...}
     3742# where axis1,... can be 1, 2 or 3 corresponding to the cartesian X, Y or Z axes
     3743#
     3744# The 20 TLS terms are ordered:
     3745#    T11, T22, T33, T12, T13, T23
     3746#    L11, L22, L33, L12, L13, L23
     3747#    S12, S13, S21, S23, S31, S32, SAA, SBB
     3748#
     3749proc ReadRigidBodyMapping {phase bodytyp num} {
     3750    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
     3751        return ""
     3752    }
     3753    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
     3754    set first [string trim [string range [readexp "$key  NDA"] 0 4]]
     3755    set line [readexp "$key BDFL"]
     3756    set varlist {}
     3757    set damplist {}
     3758    foreach i {0 1 2 3 4 5 6 7 8} {
     3759        lappend varlist [string trim [string range $line [expr {5*$i}] [expr {4 + 5*$i}] ]]
     3760        lappend damplist [string trim [string range $line [expr {45 + $i}] [expr {45 + $i}] ]]
     3761    }
     3762    set TLSdamplist {}
     3763    foreach i {54 55 56} {
     3764        lappend TLSdamplist [string trim [string range $line $i $i ]]
     3765    }
     3766    set line [readexp "${key} BDLC"]
     3767    set x [string trim [string range $line 0 9]]
     3768    set y [string trim [string range $line 10 19]]
     3769    set z [string trim [string range $line 20 29]]
     3770    set origin [list $x $y $z]
     3771    set line [readexp "${key} BDOR"]
     3772    set rotations {}
     3773    foreach i {0 10 20 30 40 50} {
     3774        set angle [string trim [string range $line $i [expr {$i+7}]]]
     3775        set axis [string trim [string range $line [expr {$i+8}] [expr {$i+9}]]]
     3776        lappend rotations [list $angle $axis]
     3777    }
     3778    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
     3779    set tlsvars {}
     3780    set tlsvals {}
     3781    if {$TLS != 0} {
     3782        set line [readexp "${key}TLSF1"]
     3783        for {set j 0} {$j < 20} {incr j} {
     3784            set var [string trim [string range $line [expr {3*$j}] [expr {3*$j+2}]]]
     3785            if {$var == ""} {set var 0}
     3786            lappend tlsvars $var
     3787        }
     3788        for {set j 0} {$j < 20} {incr j} {
     3789            set i 0
     3790            if {$j == 0} {
     3791                set i 1
     3792            } elseif {$j == 8} {
     3793                set i 2
     3794            } elseif {$j == 16} {
     3795                set i 3
     3796            }
     3797            if {$i != 0} {
     3798                set line [readexp "${key}TLSP$i"]
     3799                set i 0
     3800                set j1 0
     3801                set j2 7
     3802            } else {
     3803                incr j1 8
     3804                incr j2 8
     3805            }
     3806            set val [string trim [string range $line $j1 $j2]]
     3807            if {$val == ""} {set val 0}
     3808            lappend tlsvals $val
     3809        }
     3810    }
     3811    return [list $first $origin $rotations $varlist $damplist $tlsvals $tlsvars $TLSdamplist]
     3812}
     3813
     3814# Control TLS representation for phase, body # and instance number of a Rigid body mapping
     3815#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
     3816# Enable TLS use if TLS is non-zero (true). Disable if zero
     3817proc RigidBodyEnableTLS {phase bodytyp num TLS} {
     3818    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
     3819        return ""
     3820    }
     3821    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
     3822    if {$TLS} {
     3823        setexp "${key} LSTF" [format "%5d" 1] 1 5
     3824        if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
     3825        if {![existsexp "${key}TLSP1"]} {
     3826            makeexprec "${key}TLSP1"
     3827            set str {}
     3828            foreach v {.01 .01 .01 0 0 0 0 0} d {4 4 4 4 4 4 2 2} {
     3829                validreal v 8 $d
     3830                append str $v
     3831            }
     3832            setexp "${key}TLSP1" $str 1 64
     3833        }
     3834        if {![existsexp "${key}TLSP2"]} {
     3835            makeexprec "${key}TLSP2"
     3836            set str {}
     3837            set v 0
     3838            foreach d {2 2 2 2 4 4 4 4} {
     3839                validreal v 8 $d
     3840                append str $v
     3841            }
     3842            setexp "${key}TLSP2" $str 1 64
     3843        }
     3844        if {![existsexp "${key}TLSP3"]} {
     3845            makeexprec "${key}TLSP3"
     3846            set str {}
     3847            set v 0
     3848            foreach d {4 4 4 4} {
     3849                validreal v 8 $d
     3850                append str $v
     3851            }
     3852            setexp "${key}TLSP3" $str 1 64
     3853        }
     3854    } else {
     3855        setexp "${key} LSTF" [format "%5d" 0] 1 5
     3856    }
     3857    return 1
     3858}
     3859
     3860# Control the TLS values for Rigid body mapping for mapping with
     3861#    phase ($phase), body type # ($bodytyp) and instance # ($num)
     3862# set the 20 TLS values to the values in TLSvals
     3863# There must be exactly 20 TLS terms, which are ordered:
     3864#    T11, T22, T33, T12, T13, T23
     3865#    L11, L22, L33, L12, L13, L23
     3866#    S12, S13, S21, S23, S31, S32, SAA, SBB
     3867proc RigidBodySetTLS {phase bodytyp num TLSvals} {
     3868    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
     3869        return ""
     3870    }
     3871    if {[llength $TLSvals] != 20} {return ""}
     3872    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
     3873    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
     3874    if {$TLS == 0} {return ""}
     3875    if {![existsexp "${key}TLSF1"]} {makeexprec "${key}TLSF1"}
     3876    foreach n {1 2 3} {
     3877        if {![existsexp "${key}TLSP$n"]} {makeexprec "${key}TLSP$n"}
     3878    }
     3879    set str {}
     3880    set n 1
     3881    set i 0
     3882    foreach v $TLSvals d {4 4 4 4 4 4 2 2 2 2 2 2 4 4 4 4 4 4 4 4} {
     3883        incr i
     3884        validreal v 8 $d
     3885        append str $v
     3886        if {$i == 8} {
     3887            set i 0
     3888            setexp "${key}TLSP$n" $str 1 64
     3889            incr n
     3890            set str {}
     3891        }
     3892    }
     3893    setexp "${key}TLSP$n" $str 1 64
     3894    return 1
     3895}
     3896
     3897# set damping values for a Rigid body mapping
     3898#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
     3899# there must be 9 damping values in RBdamp for the 9 position variables (origin followed by rotations)
     3900# Use of TLSdamp is optional, but to be used, TLS representation must be enabled and there must be
     3901# three damping terms (for all T terms; for all L terms and for all S terms)
     3902proc RigidBodySetDamp {phase bodytyp num RBdamp "TLSdamp {}"} {
     3903    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
     3904        return ""
     3905    }
     3906    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
     3907    if {[llength $RBdamp] != 9} {return ""}
     3908    set str {}
     3909    foreach v $RBdamp {
     3910        if {[validint v 1] != 1} {set v " "}
     3911        append str $v
     3912    }
     3913    setexp "$key BDFL" $str 46 9
     3914    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
     3915    if {$TLS != 0 &&  [llength $TLSdamp] == 3} {
     3916        set str {}
     3917        foreach v $TLSdamp {
     3918        if {[validint v 1] != 1} {set v " "}
     3919            append str $v
     3920        }
     3921        setexp "$key BDFL" $str 55 3
     3922    }
     3923    return 1
     3924}
     3925
     3926# set refinement variable numbers for a Rigid body mapping
     3927#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
     3928# there must be 9 variable values in RBvar for the 9 position variables (origin followed by rotations)
     3929# note that the variable values should be unique integers
     3930proc RigidBodyVary {phase bodytyp num RBvar} {
     3931    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
     3932        return ""
     3933    }
     3934    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
     3935    if {[llength $RBvar] != 9} {return ""}
     3936    set str {}
     3937    foreach v $RBvar {
     3938        if {[validint v 5] != 1} {set v " "}
     3939        append str $v
     3940    }
     3941    setexp "$key BDFL" $str 1 45   
     3942}
     3943
     3944# set TLS refinement variable numbers for a Rigid body mapping
     3945#   for mapping with phase ($phase), body type # ($bodytyp) and instance # ($num)
     3946# there must be 20 variable values in TLSvar for the 20 parameters:
     3947#    T11, T22, T33, T12, T13, T23
     3948#    L11, L22, L33, L12, L13, L23
     3949#    S12, S13, S21, S23, S31, S32, SAA, SBB
     3950# note that the variable values should be unique integers
     3951proc RigidBodyTLSVary {phase bodytyp num TLSvar} {
     3952    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
     3953        return ""
     3954    }
     3955    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $num 1]"
     3956    if {[llength $TLSvar] != 20} {return ""}
     3957    set TLS [string trim [string range [readexp "${key} LSTF"] 0 4]]
     3958    if {$TLS == 0} {return ""}
     3959    set str {}
     3960    foreach v $TLSvar {
     3961        if {[validint v 3] != 1} {set v " "}
     3962        append str $v
     3963    }
     3964    setexp "${key}TLSF1" $str 1 60
     3965
     3966# AddRigidBody: add a new rigid body definition into the .EXP file
     3967# arguments are:
     3968#   multlist: defines a list of multipliers for each set of coordinates. In the
     3969#             simplest case this will be {1}
     3970#   coordlist: a nested list of coordinates such as { { {0 0 0} {.1 .1 .1} {.2 .2 .2} } }
     3971# note that when the length of multlist > 1 then coordlist must have the same length.
     3972# for input where
     3973#     multlist = {s1 s2} and
     3974#     coordlist = { { {0 0 0} {1 1 0} {.0 .0 .0} ...}
     3975#                     {0 0 0} {1 1 0} {2 1 2} ...}
     3976#                 }
     3977# the cartesian coordinates are defined from the input as
     3978#    atom 1 = s1 * (0,0,0) + s2*(0,0,0) [ = (0,0,0)]
     3979#    atom 2 = s1 * (1,1,0) + s2*(1,1,0) [ = (s1+s2) * (1,1,0)]
     3980#    atom 3 = s1 * (0,0,0) + s2*(2,1,2) [ = s2 * (2,1,2)]
     3981#    ...
     3982# Returns the number of the rigid body that has been created
     3983proc AddRigidBody {multlist coordlist} {
     3984    # find the first unused body #
     3985    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
     3986        set value $rbnum
     3987        validint value 2
     3988        set key "RGBD${value}"
     3989        if {! [existsexp "$key NATR "]} {break}
     3990    }
     3991    # did we go too far?
     3992    if {$rbnum == 16} {return ""}
     3993    # increment the RB counter
     3994    set n [string trim [readexp "RGBD  NRBDS"]]
     3995    if {$n == ""} {
     3996        makeexprec "RGBD  NRBDS"
     3997        set n 0
     3998    }
     3999    incr n
     4000    validint n 5
     4001    setexp "RGBD  NRBDS" $n 1 5
     4002    SetRigidBody $rbnum $multlist $coordlist
     4003    return $rbnum
     4004}
     4005
     4006# DeleteRigidBody: remove a rigid body definition from the .EXP file
     4007# The body may not be mapped. I am not sure if GSAS allows more than 9 bodies,
     4008# but if it does, the simplifed approach used here will fail, so this
     4009# is not allowed.
     4010# Input:
     4011#   Rigid body number
     4012# Returns:
     4013#   1 on success
     4014#   -1 if the body number is 11 or greater
     4015#   -2 if the body is mapped
     4016#   -3 if the body is not defined
     4017proc DeleteRigidBody {rbnum} {
     4018    # can't delete bodies with numbers higher than 10, since the key prefix
     4019    # (RGBD11... will overlap with rigid body instance records, which would be
     4020    # deleted below
     4021    if {$rbnum > 10} {
     4022        return -1
     4023    }
     4024    set value $rbnum
     4025    validint value 2
     4026    set key "RGBD${value}"
     4027    if {![existsexp "$key NATR "]} {
     4028        return -2
     4029    }
     4030    # make sure the body is not mapped
     4031    if {[string trim [string range [readexp "$key NBDS"] 0 4]] != 0} {
     4032        return -3
     4033    }
     4034    # delete the records starting with "RGBD x" or "RGBD10"
     4035    foreach key [array names ::exparray "${key}*"] {
     4036        #puts $key
     4037        delexp $key
     4038    }
     4039    # decrement the RB counter
     4040    set n [string trim [readexp "RGBD  NRBDS"]]
     4041    if {$n == ""} {
     4042        set n 0
     4043    }
     4044    incr n -1
     4045    validint n 5
     4046    if {$n > 0} {
     4047        setexp "RGBD  NRBDS" $n 1 5
     4048    } else {
     4049        delexp "RGBD  NRBDS"
     4050    }
     4051    return 1
     4052}
     4053
     4054# ReplaceRigidBody: replace all the information for rigid body #rbnum
     4055# Works the sames as AddRigidBody (see above) except that the rigid body is replaced rather
     4056# than added.
     4057# Note that count of the # of times the body is used is preserved
     4058proc ReplaceRigidBody {rbnum multlist coordlist {varlist ""} {damplist ""}} {
     4059    set value $rbnum
     4060    validint value 2
     4061    set key "RGBD${value}"
     4062    set line [readexp "$key NBDS"]
     4063    foreach key [array names ::exparray "${key}*"] {
     4064        #puts $key
     4065        delexp $key
     4066    }
     4067    SetRigidBody $rbnum $multlist $coordlist $varlist $damplist
     4068    setexp "$key NBDS" $line 1 68
     4069}
     4070
     4071# Edit the parameters for rigid body #rbnum
     4072# (normally called from ReplaceRigidBody or AddRigidBody)
     4073proc SetRigidBody {rbnum multlist coordlist {varlist ""} {damplist ""}} {
     4074    set value $rbnum
     4075    validint value 2
     4076    set key "RGBD${value}"
     4077    # number of atoms
     4078    set value [llength [lindex $coordlist 0]]
     4079    validint value 5
     4080    makeexprec "$key NATR"
     4081    setexp "$key NATR" $value 1 5
     4082    # number of times used
     4083    set value 0
     4084    validint value 5
     4085    makeexprec "$key NBDS"
     4086    setexp "$key NBDS" $value 1 5
     4087    # number of coordinate matrices
     4088    set value [llength $multlist]
     4089    validint value 5
     4090    makeexprec "$key NSMP"
     4091    setexp "$key NSMP" $value 1 5
     4092    set i 0
     4093    foreach mult $multlist coords $coordlist {
     4094        set var [lindex $varlist $i]
     4095        if {$var == ""} {set var 0}
     4096        set damp [lindex $damplist $i]
     4097        if {$damp == ""} {set damp 0}
     4098        incr i
     4099        makeexprec "${key}${i}PARM"
     4100        setexp "${key}${i}PARM" [format "%10.5f%5d%5d" $mult $var $damp] 1 20
     4101        set j 0
     4102        foreach item $coords {
     4103            #puts $item
     4104            incr j
     4105            set value $j
     4106            validint value 3
     4107            makeexprec "${key}${i}SC$value"
     4108            if {[llength $item] == 4} {
     4109                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f%10s" $item] 1 40
     4110            } elseif {[llength $item] == 3} {
     4111                setexp "${key}${i}SC$value" [eval format "%10.6f%10.6f%10.6f" $item] 1 30
     4112            } else {
     4113                return -code 3 "Invalid number of coordinates"
     4114            }
     4115        }
     4116    }
     4117}
     4118
     4119# convert a decimal to the GSAS hex encoding with a field $digits long.
     4120proc ToHex {num digits} {
     4121    return [string toupper [format "%${digits}x" $num]]
     4122}
     4123
     4124# convert a GSAS hex encoding to a decimal integer
     4125proc FromHex {hex} {
     4126    return [scan $hex "%x"]
     4127}
     4128
     4129# MapRigidBody: define an "instance" of a rigid body: meaning that the coordinates
     4130# (and optionally U values) for a set of atoms will be generated from the rigid body
     4131# arguments:
     4132#   phase: phase number (1-9)
     4133#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
     4134#   firstatom: sequence number of the first atom in phase (note that atoms may
     4135#              not be numbered sequentially)
     4136#   position: list of three fractional coordinates for the origin of the rigid body coordinates
     4137#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
     4138#           cartesian system before the body is translated to position.
     4139# returns the instance # (number of times body $bodytyp has been used in phase $phase)
     4140proc MapRigidBody {phase bodytyp firstatom position angles} {
     4141    # find the first unused body # for this phase & type
     4142    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
     4143        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
     4144        if {! [existsexp "$key  NDA"]} {break}
     4145    }
     4146    # did we go too far?
     4147    if {$rbnum == 16} {return ""}
     4148    # increment number of mapped bodies of this type overall
     4149    set value $bodytyp
     4150    validint value 2
     4151    set key "RGBD${value}"
     4152    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
     4153    incr used
     4154    set value $used
     4155    validint value 5
     4156    setexp "$key NBDS" $value 1 5
     4157    # increment number of mapped bodies of this type in this phase
     4158    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
     4159    if {[existsexp "$key  NBDS"]} {
     4160        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
     4161    } else {
     4162        makeexprec "$key  NBDS"
     4163        set used 0
     4164    }
     4165    incr used
     4166    set value $used
     4167    validint value 5
     4168    setexp "$key  NBDS" $value 1 5
     4169    # now write the mapping parameters
     4170    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
     4171    set value $firstatom
     4172    validint value 5
     4173    makeexprec "$key  NDA"
     4174    setexp "$key  NDA" $value 1 5
     4175    set l1 {}
     4176    set l2 {}
     4177    for {set i 0} {$i < 9} {incr i} {
     4178        append l1 [format %5d 0]
     4179        append l2 [format %1d 0]
     4180    }
     4181    makeexprec "$key BDFL"
     4182    setexp "$key BDFL" $l1$l2 1 54
     4183    makeexprec "${key} BDLC"
     4184    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
     4185    makeexprec "${key} BDOR"
     4186    set l1 {}
     4187    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
     4188        append l1 [format "%8.2f%2d" $val $dir]
     4189    }
     4190    setexp "${key} BDOR" $l1 1 60
     4191    makeexprec "${key} LSTF"
     4192    setexp "${key} LSTF" [format "%5d" 0] 1 5
     4193    return $rbnum
     4194}
     4195
     4196# EditRigidBodyMapping: edit parameters that define an "instance" of a rigid body (see MapRigidBody)
     4197# arguments:
     4198#   phase: phase number (1-9)
     4199#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
     4200#   bodynum: instance number, as returned by MapRigidBody
     4201#   position: list of three fractional coordinates for the origin of the rigid body coordinates
     4202#   angles: list of 3 angles to rotate the rigid body coordinates around x, y, z of the
     4203#           cartesian system before the body is translated to position.
     4204#
     4205proc EditRigidBodyMapping {phase bodytyp bodynum position angles} {
     4206    # number of bodies of this type in this phase
     4207    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $bodynum 1]"
     4208    setexp "${key} BDLC" [eval format "%10.6f%10.6f%10.6f" $position] 1 30
     4209    set l1 {}
     4210    foreach val "$angles 0 0 0" dir "1 2 3 1 1 1" {
     4211        append l1 [format "%8.2f%2d" $val $dir]
     4212    }
     4213    setexp "${key} BDOR" $l1 1 60
     4214}
     4215
     4216# UnMapRigidBody: remove a rigid body constraint by removing a RB "instance"
     4217# (undoes MapRigidBody)
     4218# arguments:
     4219#   phase: phase number (1-9)
     4220#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
     4221#   bodynum: instance number, as returned by MapRigidBody
     4222proc UnMapRigidBody {phase bodytyp mapnum} {
     4223    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $mapnum] == -1} {
     4224        return ""
     4225    }
     4226    # decrement number of mapped bodies of this type overall
     4227    set value $bodytyp
     4228    validint value 2
     4229    set key "RGBD${value}"
     4230    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
     4231    incr used -1
     4232    set value $used
     4233    validint value 5
     4234    setexp "$key NBDS" $value 1 5
     4235    # decrement number of mapped bodies of this type in this phase
     4236    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
     4237    if {[existsexp "$key  NBDS"]} {
     4238        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
     4239    } else {
     4240        set used 0
     4241    }
     4242    incr used -1
     4243    set value $used
     4244    validint value 5
     4245    if {$used > 0} {
     4246        setexp "$key  NBDS" $value 1 5
     4247    } else {
     4248        delexp "$key  NBDS"
     4249    }
     4250    # now delete the mapping parameter records
     4251    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $mapnum 1]"
     4252    foreach key [array names ::exparray "${key}*"] {
     4253        delexp $key
     4254    }
     4255    return $used
     4256}
     4257
Note: See TracChangeset for help on using the changeset viewer.