Ignore:
Timestamp:
Jan 21, 2011 4:43:16 PM (10 years ago)
Author:
toby
Message:

fix save as bug, expand rb.tcl to delete & handle gaps in body lists; likewise in rigid.tcl

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/rb.tcl

    r1108 r1112  
    55#============================================================================
    66# rigid body EXP editing routines (to move into readexp.tcl)
    7 # RigidBodyCount -- returns the number of defined rigid bodies (body types)
    8 # ReadRigidBody  -- # of times a body is mapped & scaling factors
    9 # RigidBodyMappingCount -- # of times a rigid body is mapped in phase
     7# RigidBodyList -- returns a list of the defined rigid body types
     8# ReadRigidBody  -- # of times a body is mapped, scaling factors, var #s & coordinates
     9# RigidBodyMappingList - return a list instances where a RB is mapped in phase
    1010# RigidBodyEnableTLS -- Enable or Disable TLS use for a rigid body mapping
    1111# RigidBodySetTLS  -- change the TLS values for a rigid body mapping
     
    1414# RigidBodyTLSVary -- set TLS refinement variable nums for a rigid body mapping
    1515# AddRigidBody -- defines a new rigid body type
     16# DeleteRigidBody -- remove a rigid body definition
    1617# ReplaceRigidBody -- replaces a previous rigid body type
    1718# ReadRigidBodyMapping  -- get parameters for a rigid body mapping
    1819# MapRigidBody -- map a rigid body type into a phase
    1920# EditRigidBodyMapping -- change the parameters in a rigid body mapping
    20 #  need to unmap a rigid body
    21 #  delete rigid body
     21# UnMapRigidBody --remove a rigid body constraint by removing a RB "instance"
     22#----- note that these older routines should not be used ------
     23# RigidBodyCount -- returns the number of defined rigid bodies (body types)
     24#    use RigidBodyList instead
     25# RigidBodyMappingCount -- # of times a rigid body is mapped in phase
     26#    use RigidBodyMappingList instead
    2227#============================================================================
    2328# returns the number of defined rigid bodies
     
    2833    }
    2934    return $n
     35}
     36
     37# returns a list of the defined rigid body types
     38proc RigidBodyList {} {
     39    set n [string trim [readexp "RGBD  NRBDS"]]
     40    if {$n == ""} {
     41        set n 0
     42    }
     43    set rblist {}
     44    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
     45        set value $rbnum
     46        validint value 2
     47        set key "RGBD${value}"
     48        if {[existsexp "$key NATR "]} {
     49            lappend rblist $rbnum
     50        }
     51        if {[llength $rblist] == $n} break
     52    }
     53    return $rblist
    3054}
    3155
     
    4266#  not used in GSAS.
    4367proc ReadRigidBody {rbnum} {
    44     if {[RigidBodyCount] < $rbnum} {
     68    if {[lsearch [RigidBodyList] $rbnum] == -1} {
    4569        return ""
    4670    }
     
    83107    return $n
    84108}
     109# return a list of the instances where rigid body $bodytyp is mapped in phase $phase
     110proc RigidBodyMappingList {phase bodytyp} {
     111    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
     112    if {! [existsexp "$key  NBDS"]} {return {}}
     113    set n [string trim [readexp "$key  NBDS"]]
     114    if {$n == ""} {
     115        set n 0
     116    }
     117    set rblist {}
     118    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15} {
     119        set value $rbnum
     120        validint value 2
     121        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
     122        if {[existsexp "$key  NDA"]} {
     123            lappend rblist $rbnum
     124        }
     125        if {[llength $rblist] == $n} break
     126    }
     127    return $rblist
     128}
     129
     130
    85131
    86132# reads rigid body mapping parameters for phase ($phase), body type # ($bodytyp) and instance # ($num)
     
    106152#
    107153proc ReadRigidBodyMapping {phase bodytyp num} {
    108     if {[RigidBodyMappingCount $phase $bodytyp] < $num} {
     154    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
    109155        return ""
    110156    }
     
    170216# Enable TLS use if TLS is non-zero (true). Disable if zero
    171217proc RigidBodyEnableTLS {phase bodytyp num TLS} {
    172     if {[RigidBodyMappingCount $phase $bodytyp] < $num} {
     218    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
    173219        return ""
    174220    }
     
    220266#    S12, S13, S21, S23, S31, S32, SAA, SBB
    221267proc RigidBodySetTLS {phase bodytyp num TLSvals} {
    222     if {[RigidBodyMappingCount $phase $bodytyp] < $num} {
     268    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
    223269        return ""
    224270    }
     
    255301# three damping terms (for all T terms; for all L terms and for all S terms)
    256302proc RigidBodySetDamp {phase bodytyp num RBdamp "TLSdamp {}"} {
    257     if {[RigidBodyMappingCount $phase $bodytyp] < $num} {
     303    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
    258304        return ""
    259305    }
     
    283329# note that the variable values should be unique integers
    284330proc RigidBodyVary {phase bodytyp num RBvar} {
    285     if {[RigidBodyMappingCount $phase $bodytyp] < $num} {
     331    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
    286332        return ""
    287333    }
     
    304350# note that the variable values should be unique integers
    305351proc RigidBodyTLSVary {phase bodytyp num TLSvar} {
    306     if {[RigidBodyMappingCount $phase $bodytyp] < $num} {
     352    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $num] == -1} {
    307353        return ""
    308354    }
     
    336382# Returns the number of the rigid body that has been created
    337383proc AddRigidBody {multlist coordlist} {
    338     #
     384    # find the first unused body #
     385    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
     386        set value $rbnum
     387        validint value 2
     388        set key "RGBD${value}"
     389        if {! [existsexp "$key NATR "]} {break}
     390    }
     391    # did we go too far?
     392    if {$rbnum == 16} {return ""}
     393    # increment the RB counter
    339394    set n [string trim [readexp "RGBD  NRBDS"]]
    340395    if {$n == ""} {
     
    345400    validint n 5
    346401    setexp "RGBD  NRBDS" $n 1 5
    347     SetRigidBody $n $multlist $coordlist
    348     return $n
     402    SetRigidBody $rbnum $multlist $coordlist
     403    return $rbnum
     404}
     405
     406# DeleteRigidBody: remove a rigid body definition from the .EXP file
     407# The body may not be mapped. I am not sure if GSAS allows more than 9 bodies,
     408# but if it does, the simplifed approach used here will fail, so this
     409# is not allowed.
     410# Input:
     411#   Rigid body number
     412# Returns:
     413#   1 on success
     414#   -1 if the body number is 11 or greater
     415#   -2 if the body is mapped
     416#   -3 if the body is not defined
     417proc DeleteRigidBody {rbnum} {
     418    # can't delete bodies with numbers higher than 10, since the key prefix
     419    # (RGBD11... will overlap with rigid body instance records, which would be
     420    # deleted below
     421    if {$rbnum > 10} {
     422        return -1
     423    }
     424    set value $rbnum
     425    validint value 2
     426    set key "RGBD${value}"
     427    if {![existsexp "$key NATR "]} {
     428        return -2
     429    }
     430    # make sure the body is not mapped
     431    if {[string trim [string range [readexp "$key NBDS"] 0 4]] != 0} {
     432        return -3
     433    }
     434    # delete the records starting with "RGBD x" or "RGBD10"
     435    foreach key [array names ::exparray "${key}*"] {
     436        #puts $key
     437        delexp $key
     438    }
     439    # decrement the RB counter
     440    set n [string trim [readexp "RGBD  NRBDS"]]
     441    if {$n == ""} {
     442        set n 0
     443    }
     444    incr n -1
     445    validint n 5
     446    if {$n > 0} {
     447        setexp "RGBD  NRBDS" $n 1 5
     448    } else {
     449        delexp "RGBD  NRBDS"
     450    }
     451    return 1
    349452}
    350453
     
    432535# returns the instance # (number of times body $bodytyp has been used in phase $phase)
    433536proc MapRigidBody {phase bodytyp firstatom position angles} {
    434     # number of bodies of this type in this phase
     537    # find the first unused body # for this phase & type
     538    foreach rbnum {1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16} {
     539        set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
     540        if {! [existsexp "$key  NDA"]} {break}
     541    }
     542    # did we go too far?
     543    if {$rbnum == 16} {return ""}
     544    # increment number of mapped bodies of this type overall
     545    set value $bodytyp
     546    validint value 2
     547    set key "RGBD${value}"
     548    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
     549    incr used
     550    set value $used
     551    validint value 5
     552    setexp "$key NBDS" $value 1 5
     553    # increment number of mapped bodies of this type in this phase
    435554    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
    436     set n [string trim [readexp "$key  NBDS"]]
    437     if {$n == ""} {
     555    if {[existsexp "$key  NBDS"]} {
     556        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
     557    } else {
    438558        makeexprec "$key  NBDS"
    439         set n 0
    440     }
    441     incr n
    442     set value $n
     559        set used 0
     560    }
     561    incr used
     562    set value $used
    443563    validint value 5
    444564    setexp "$key  NBDS" $value 1 5
    445     set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $n 1]"
     565    # now write the mapping parameters
     566    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $rbnum 1]"
    446567    set value $firstatom
    447568    validint value 5
     
    466587    makeexprec "${key} LSTF"
    467588    setexp "${key} LSTF" [format "%5d" 0] 1 5
    468     return $n
     589    return $rbnum
    469590}
    470591
     
    488609    setexp "${key} BDOR" $l1 1 60
    489610}
     611
     612# UnMapRigidBody: remove a rigid body constraint by removing a RB "instance"
     613# (undoes MapRigidBody)
     614# arguments:
     615#   phase: phase number (1-9)
     616#   bodytyp: number of rigid body (1-15) as returned from AddRigidBody
     617#   bodynum: instance number, as returned by MapRigidBody
     618proc UnMapRigidBody {phase bodytyp mapnum} {
     619    if {[lsearch [RigidBodyMappingList $phase $bodytyp] $mapnum] == -1} {
     620        return ""
     621    }
     622    # decrement number of mapped bodies of this type overall
     623    set value $bodytyp
     624    validint value 2
     625    set key "RGBD${value}"
     626    set used [string trim [string range [readexp "$key NBDS"] 0 4]]
     627    incr used -1
     628    set value $used
     629    validint value 5
     630    setexp "$key NBDS" $value 1 5
     631    # decrement number of mapped bodies of this type in this phase
     632    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1]"
     633    if {[existsexp "$key  NBDS"]} {
     634        set used [string trim [string range [readexp "$key  NBDS"] 0 4]]
     635    } else {
     636        set used 0
     637    }
     638    incr used -1
     639    set value $used
     640    validint value 5
     641    if {$used > 0} {
     642        setexp "$key  NBDS" $value 1 5
     643    } else {
     644        delexp "$key  NBDS"
     645    }
     646    # now delete the mapping parameter records
     647    set key "RGBD[ToHex $phase 1][ToHex $bodytyp 1][ToHex $mapnum 1]"
     648    foreach key [array names ::exparray "${key}*"] {
     649        delexp $key
     650    }
     651    return $used
     652}
     653
     654#============================================================================
     655# Rigid body utility routines
     656#============================================================================
     657# RigidBodyGetVarNums: Returns a list of the variable numbers in use
     658#       for rigid body variable parameters.
     659# RigidBodyAtomNums: returns a list of atom numbers that are mapped to
     660#       rigid bodies in a selected phase
     661# RigidStartAtoms: returns a list of atoms that are allowed for creation of RB
     662# ExtractRigidBody: Use the GSAS geometry program to cartesian coordinates &
     663#       setting info for a RB from fractional coordinates for atoms in a phase
     664# RunRecalcRBCoords: updates the coordinates in all phases after changes have
     665#       been made to rigid parameters.
     666# CalcBody: Convert ortho to fractional coordinates using RB parameters
     667# FitBody: Optimize the origin and Euler angles to match a rigid body to a
     668#       set of fractional coordinates
     669# zmat2coord: convert a z-matrix to a set of cartesian coordinates
     670# RB2cart: convert the representation used for rigid bodies into
     671#       cartesian coordinates
     672# PlotRBtype: plot a rigid body with DRAWxtl
     673# PlotRBcoords: plot orthogonal coordinates with DRAWxtl
     674# DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed
     675#       on a structure read from a phase with DRAWxtl
    490676#============================================================================
    491677#============================================================================
    492 # Returns a list of the variable numbers used already for rigid body variable parameters
     678# RigidBodyGetVarNums: Returns a list of the variable numbers used already
     679# for rigid body variable parameters
    493680proc RigidBodyGetVarNums {} {
    494681    set varlist {}
    495     set bodies [RigidBodyCount]
    496     for {set type 1} {$type <= $bodies} {incr type} {
     682    foreach type [RigidBodyList] {
    497683        foreach phase $::expmap(phaselist) {
    498             set count [RigidBodyMappingCount $phase $type]
    499             for {set i 1} {$i <= $bodies} {incr i} {
     684            foreach i [RigidBodyMappingList $phase $type] {
    500685                set items [ReadRigidBodyMapping $phase $type $i]
    501686                set varlist [concat $varlist [lindex $items 3]]
     
    506691        }
    507692    }
    508     return [lsort -unique $varlist]
    509 }
    510 
    511 # Use the GSAS geometry program to compute a set of cartesian coordinates for a
     693    return [lsort -integer -unique $varlist]
     694}
     695
     696# RigidBodyAtomNums: Returns a list of the atoms mapped to rigid bodies in
     697# phase $phase
     698proc RigidBodyAtomNums {phase} {
     699    if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""}
     700    set allatoms $::expmap(atomlist_$phase)
     701    # get matching atoms coordinate range
     702    set mappedlist {}
     703    foreach type [RigidBodyList] {
     704        foreach i [RigidBodyMappingList $phase $type] {
     705            # get the number of atoms in this type of body
     706            set natoms [llength [lindex [lindex [lindex [ReadRigidBody $type] 1] 0] 3]]
     707            set natom1 [expr {$natoms - 1}]
     708            set items [ReadRigidBodyMapping $phase $type $i]
     709            set firstatom [lindex $items 0]
     710            set firstind [lsearch $allatoms $firstatom]
     711            set mappedlist [concat $mappedlist \
     712                                [lrange \
     713                                     [lrange $allatoms $firstind end] \
     714                                     0 $natom1] \
     715                               ]
     716        }
     717    }
     718    return [lsort -integer $mappedlist]
     719}
     720
     721# RigidStartAtoms: Find allowed starting atoms for a rigid body in a phase
     722# Input:
     723#   phase is the phase number
     724#   natoms is the number of atoms in the RB to be mapped
     725# Returns a list of valid "start" atoms.
     726# Example: if the atom numbers in the phase are {2 4 5 6 7 8} and no rigid bodies
     727# are mapped, then a 4-atom body can be mapped starting with atom 2, 4 or 5 only,
     728# so {2 4 5} is returned
     729# If atoms 2-6 were already mapped, then this routine would return an empty
     730# list, as it is not possible to map the body.
     731proc RigidStartAtoms {phase natoms} {
     732    if {[lsearch $::expmap(phaselist) $phase] == -1} {return ""}
     733    set allatoms $::expmap(atomlist_$phase)
     734    set usedatoms [RigidBodyAtomNums $phase]
     735    set startatomlist {}
     736    for {set i 0} {$i < [llength $allatoms]} {incr i} {
     737        set al [lrange $allatoms $i [expr {$i+$natoms-1}]]
     738        if {[llength $al] < $natoms} break
     739        set ok 1
     740        foreach atom $al {
     741            if {[lsearch $usedatoms $atom] != -1} {
     742                set ok 0
     743                break
     744            }
     745        }
     746        if $ok {lappend startatomlist [lindex $al 0]}
     747    }
     748    return $startatomlist
     749}
     750
     751# ExtractRigidBody: Use the GSAS geometry program to compute a set of cartesian coordinates for a
    512752# set of atoms in a .EXP file and provide the origin shift and Euler angles needed to
    513753# place the cartesian system into the crystal coordinates. Used for setting up a rigid body.
     
    584824}
    585825
    586 # updates the coordinates in a .EXP file after a rigid body instance has been added/changed
     826# RunRecalcRBCoords: updates the coordinates in a .EXP file after a rigid
     827# body has been changed, mapped or the setting info is changed
    587828proc RunRecalcRBCoords { } {
    588829    global expgui tcl_platform
     
    7761017}
    7771018
    778 #calc rigid body fractional coordinates
     1019# CalcBody: Calculate fractional coordinates using rigid body setting parameters
    7791020# arguments:
    7801021#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
     
    8961137}
    8971138
    898 # fit a rigid body's Origin and Euler angles
     1139# FitBody: fit a rigid body's Origin and Euler angles
    8991140# arguments:
    9001141#  Euler: a list of axes and angles to rotate: { {axis1 angle1} {axis2 angle2} ...}
     
    9471188}
    9481189
    949 # zmat2coord converts a z-matrix to a set of cartesian coordinates
     1190# zmat2coord: convert a z-matrix to a set of cartesian coordinates
    9501191#   a z-matrix is also known as "internal coordinates" or "torsion space"
    9511192#   (see Journal of Computational Chemistry, Vol 26, #10, p. 1063–1068, 2005 or
     
    10741315}
    10751316 
    1076 # convert the rigid body representation reported as the 2nd element in ReadRigidBody
    1077 # into cartesian coordinates
     1317# RB2cart: convert the rigid body representation reported as the 2nd element
     1318# in ReadRigidBody into cartesian coordinates
    10781319#   rblist: a list containing an element for each scaling factor
    10791320# in each element there are four items:
     
    11271368}
    11281369
    1129 # plot orthogonal coordinates in DRAWxtl
     1370# DRAWxtlPlotOrtho: plot orthogonal coordinates in DRAWxtl
    11301371# input:
    11311372#  filename: file name for the .str file to create
     
    11971438}
    11981439
    1199 # plot a rigid body in DRAWxtl
     1440# PlotRBtype: plot a rigid body in DRAWxtl
    12001441# input:
    12011442#  rbtype: # of rigid body
     
    12211462}
    12221463
    1223 # plot orthogonal coordinates in DRAWxtl
     1464# PlotRBcoords: plot orthogonal coordinates in DRAWxtl
    12241465# input:
    12251466#  coords: cartesian coordinates
     
    12441485}
    12451486
    1246 # plot a rigid body superimposed on a structure
     1487# DRAWxtlPlotRBFit: plot a set of fraction coordinates superimposed
     1488# on a structure read from a phase
    12471489# input:
    12481490#  RBcoords: fractional coordinates for rigid body
    12491491#  phase:# of phase to plot
    1250 #  firstatom: # of 1st atom in structure matching rigid body
     1492#  firstatom: seq # of 1st atom in structure to be mapped to rigid body
    12511493#  allatoms: 0 to plot only atoms in phase that are in the rigid body,
    12521494#      otherwise plot all atoms
     
    12951537                      [expr {($zmin+$zmax)/2}] ]
    12961538    # get matching atoms coordinate range
     1539    set firstind [lsearch $::expmap(atomlist_$phase) $firstatom]
    12971540    foreach atom [lrange \
    1298                       [lrange $::expmap(atomlist_$phase) [expr {$firstatom-1}] end] \
     1541                      [lrange $::expmap(atomlist_$phase) $firstind end] \
    12991542                      0 [expr {$natom-1}]] {
    13001543        foreach s {x y z} {
     
    13341577        set atoms $::expmap(atomlist_$phase)
    13351578    } else {
     1579        set firstind [lsearch $::expmap(atomlist_$phase) $firstatom]
    13361580        set atoms [lrange \
    1337                        [lrange $::expmap(atomlist_$phase) [expr {$firstatom-1}] end] \
     1581                       [lrange $::expmap(atomlist_$phase) $firstind end] \
    13381582                       0 [expr {$natom-1}]]
    13391583    }
Note: See TracChangeset for help on using the changeset viewer.