Changeset 1213 for branches/sandbox


Ignore:
Timestamp:
Aug 16, 2012 11:44:38 PM (8 years ago)
Author:
toby
Message:

routine for chemical restraint added (ChemConst?)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/sandbox/readexp.tcl

    r1211 r1213  
    35893589        default {
    35903590            set msg "Unsupported phaseinfo access: parm=$parm action=$action"
     3591            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
     3592        }
     3593    return 1
     3594    }
     3595}
     3596
     3597# read/edit chemical restraint info
     3598#  parm:
     3599#    weight -- histogram weight (factr) *
     3600#    restraintlist -- list of restraints *
     3601#  action: get (default) or set
     3602#  value: used only with set
     3603#      value is a list of constraints
     3604#      each constrain contains {sum esd cons1 cons2...}
     3605#      each consN contains {phase atomnum multiplier}
     3606#  * =>  read+write supported
     3607# Examples:
     3608#
     3609#ChemConst restraintlist set { {10 1.1 {1 1 2} {2 2 3}} {0 1 {1 1 1} {1 2 -2}} }
     3610#
     3611#ChemConst restraintlist get
     3612#{10.00000 1.10000 {1 1 2.00000} {2 2 3.00000}} {0.00000 1.00000 {1 1 1.00000} {1 2 -2.00000}}
     3613# constraint one 2*(1:1) + 3*(2:2) = 10(1.1)
     3614# constraint two 1*(1:1) - 2*(1:2) = 0(1)
     3615#   where (1:2) is the total number of atoms (multiplicity*occupancy) for atom 2 in phase 1
     3616
     3617proc ChemConst {parm "action get" "value {}"} {
     3618    set HST {}
     3619    # look for CMP record
     3620    set n 0
     3621    for {set i 0} {$i < $::expmap(nhst)} {incr i} {
     3622        set ihist [expr {$i + 1}]
     3623        if {[expr {$i % 12}] == 0} {
     3624            incr n
     3625            set line [readexp " EXPR  HTYP$n"]
     3626            if {$line == ""} {
     3627                set msg "No HTYP$n entry for Histogram $ihist. This is an invalid .EXP file"
     3628                tk_dialog .badexp "Error in readexp" $msg error 0 Exit
     3629            }
     3630            set j 0
     3631        } else {
     3632            incr j
     3633        }
     3634        if {[string range $line [expr 2+5*$j] [expr 5*($j+1)]] == "CMP "} {
     3635            set HST $ihist
     3636        }
     3637    }
     3638    if {$HST <=9} {
     3639        set key "HST  $HST"
     3640    } else {
     3641        set key "HST $HST"
     3642    }
     3643    if {$HST == "" && $action == "set"} {
     3644        # no CMP found need to add the soft constr. histogram
     3645        # increment number of histograms
     3646        set hst [string trim [string range [readexp { EXPR  NHST }] 0 4]]
     3647        incr hst
     3648        set HST $hst
     3649        if ![validint hst 5] {return 0}
     3650        setexp  { EXPR  NHST } $hst 1 5
     3651        # add to EXPR HTYPx rec, creating if needed
     3652        set n [expr { 1+ (($HST - 1) / 12) }]
     3653        set key " EXPR  HTYP$n"
     3654        if {[array names ::exparray $key] == ""} {
     3655            makeexprec $key
     3656        }
     3657        setexp $key "CMP " [expr 3 + 5*(($HST-1) % 12)] 5
     3658        # create other HST  xx recs
     3659        if {$HST <=9} {
     3660            set key "HST  $HST"
     3661        } else {
     3662            set key "HST $HST"
     3663        }
     3664        makeexprec "$key  HNAM"
     3665        setexp "$key  HNAM" "Chemical composition restraints" 3 31
     3666        makeexprec "$key FACTR"
     3667#       makeexprec "$key NBNDS"
     3668        makeexprec "$key NCMPS"
     3669        mapexp
     3670    } elseif {$HST == ""} {
     3671        if $::expgui(debug) {puts "no restraints"}
     3672        return "1"
     3673    }
     3674
     3675    switch -glob ${parm}-$action {
     3676        weight-get {
     3677            return [string trim [string range [readexp "$key FACTR"] 0 14]]
     3678        }
     3679        weight-set {
     3680            # update FACTR
     3681            if ![validreal value 15 6] {return 0}
     3682            setexp "$key FACTR" $value 1 15
     3683        }
     3684        restraintlist-get {
     3685            set ncons [string trim [string range [readexp "$key NCMPS"] 0 4]]
     3686            set conslist {}
     3687            for {set i 1} {$i <= $ncons} {incr i} {
     3688                set const {}
     3689                set line [readexp "${key} CM$i  "]
     3690                # number of terms
     3691                set nterm [string trim [string range $line 0 4]]
     3692                if {$nterm == ""} {set nterm 0}
     3693                # chemical sum and esd
     3694                lappend const [string trim [string range $line 5 14]]
     3695                lappend const [string trim [string range $line 15 24]]
     3696                for {set j 1} {$j <= $nterm} {incr j} {
     3697                    set n [expr {($j + 2)/3}]
     3698                    set o1 [expr {20*(($j-1)%3)}]
     3699                    set o2 [expr {19 + 20*(($j-1)%3)}]
     3700                    validint n 2
     3701                    if {$o1 == 0} {
     3702                        set line [readexp "${key} CM${i}${n}"]
     3703                    }
     3704                    set frag [string range $line $o1 $o2]                   
     3705                    lappend const [list \
     3706                                       [string trim [string range $frag 0 4]] \
     3707                                       [string trim [string range $frag 5 9]] \
     3708                                       [string trim [string range $frag 10 19]] \
     3709                                   ]
     3710                }
     3711                lappend conslist $const
     3712            }
     3713            return $conslist
     3714        }
     3715        restraintlist-set {
     3716            set num [llength $value]
     3717            if ![validint num 5] {return 0}
     3718            setexp "$key NCMPS" $num 1 5
     3719            # delete all old records
     3720            foreach i [array names ::exparray "${key} CM*"] {
     3721                unset ::exparray($i)
     3722            }
     3723            set i 0
     3724            foreach cons $value {
     3725                incr i
     3726                set sum [lindex $cons 0]
     3727                set esd [lindex $cons 1]
     3728                set terms [lrange $cons 2 end]
     3729                set nterms [llength $terms]
     3730                validint nterms 5
     3731                validreal sum 10 5
     3732                validreal esd 10 5
     3733                makeexprec "${key} CM$i  "
     3734                setexp "${key} CM$i  " "${nterms}${sum}${esd}" 1 25
     3735                set j 0
     3736                set str {}
     3737                foreach term $terms {
     3738                    incr j
     3739                    set n [expr {($j + 2)/3}]
     3740                    validint n 2
     3741                    foreach {phase atom mult} $term {}
     3742                    validint phase 5
     3743                    validint atom 5
     3744                    validreal mult 10 5
     3745                    append str "${phase}${atom}${mult}"
     3746                    if {[expr {$j%3}] == 0} {
     3747                        #puts [readexp "${key} CM${i}${n}"]
     3748                        makeexprec "${key} CM${i}${n}"
     3749                        setexp "${key} CM${i}${n}" $str 1 60
     3750                        set str {}
     3751                    }
     3752                }
     3753                if {[string length $str] > 0} {
     3754                    makeexprec "${key} CM${i}${n}"
     3755                    setexp "${key} CM${i}${n}" $str 1 60
     3756                }
     3757            }
     3758        }
     3759        default {
     3760            set msg "Unsupported phaseinfo access: parm=$parm action=$action"
     3761            puts $msg
    35913762            tk_dialog .badexp "Error in readexp" $msg error 0 Exit
    35923763        }
Note: See TracChangeset for help on using the changeset viewer.