Minor updates and fixes for constant-pH MD 93/3893/3
authorradakb <brian.radak@gmail.com>
Fri, 23 Mar 2018 17:41:09 +0000 (13:41 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Fri, 23 Mar 2018 20:19:18 +0000 (15:19 -0500)
The changes might look much more comprehensive than they actually
are...(this one kind of got away from me).

New features/fixes:
-------------------

- Official support for diprotic phosphate monoesters
  (parameterization pending, publish soon?)

- The pH is now stored in the restart file and therefore does NOT
  need to be set as a keyword during a restart. This is in
  preparation for replica exchange and/or other simulation modes
  that will change the pH value

- Multiple config files can now be read. These will overwrite (and
  warn the user about) duplicate titratable residue definitions,
  with preference for the last loaded.

- MC statistics are now stored in the restart file and accumulated
  over mutiple restarts. These may not make sense if the inherent
  pKa is updated between runs and this might get updated in the
  future

- Fix expected behavior when changing weights and switch times
  after a restart. Previously the restart file was overwriting
  keyword specifications, which are supposed to always have
  preference.

- See previous. The restart format is now a bit different, so OLD
  RESTART FILES MAY NOT HAVE THE BEHAVIOR ADVERTISED HERE! The
  JSON format dict is now minimal and entries that match the
  default are removed. This decreases the file size by about half
  (though it was only order of KB before).

- Hopefully more robust error checking and responses when the user
  gives bad input specs. Things now crash when you probably want
  them to rather than just ignoring the input.

Code clean-up and other backend things the user doesn't see:
------------------------------------------------------------

- Harmonized the Tcl switch idiom across all files - this was
  actually a surpising thing that I didn't know. Tcl switches do
  not adhere to "break" statements and _return_ the final
  expression within the given selection. Thus you can wrap a
  switch statement in []s and get the value back (cool!).

- Change lsearch statements to more readable "in" and "ni"
  statements (this might only work for Tcl 8.5 and later)

- simplified numtcl interface:
    - choice now subsumes weightedChoice (which no longer exists)
      and accepts optional weights. The result now includes the
      list element AND index

    - lincr now subsumes llincr (which no longer exists).

    - lmultiply and lincr now have broadcast behavior similiar to
      Numpy ndarrays. That is, you can add or multiply a list
      in-place by either a value or a list of identical size (in an
      element by element fashion). Both functions follow the
      indexing behavior or lindex, lset, etc. when using nested
      lists

- Massive overhaul of cphTitrator to be analogous to cphSystem
    - The namespace now exports a _single_ function, cphtitrator
      which accepts actions (e.g. get, set, propose) as the first
      argument

    - fuction documentation! hey, how about that?

    - internal data structures are almost totally different

- Clean-up/re-write of matrix reading for residue definitions. This
  started getting really ugly for large residues, so a whole bunch
  of shorthand convenience functions were introduced to make it
  bearable to add new residue patterns.

Change-Id: Ib041f5e70cd25aa14aabe9a2a8ae5f7384b211aa

lib/namdcph/namdcph/cphtitrator.tcl
lib/namdcph/namdcph/cphtoppar.tcl
lib/namdcph/namdcph/namdcph.core.tcl
lib/namdcph/namdcph/numtcl.tcl

index a0cc9c5..c9b9154 100644 (file)
@@ -6,8 +6,7 @@
 #
 #   If cphSystem and cphTitrator were true objects, this would essentially be
 # an object composition. That is, a cphTitrator "contains" a cphSystem and
-# yields information pertaining to it depending on a context (most often, the
-# pH).
+# yields information pertaining to it depending on a context (e.g. the pH).
 #
 package require Tcl 8.5
 
@@ -16,107 +15,151 @@ source [file join [file dirname [info script]] "numtcl.tcl"]
 source [file join [file dirname [info script]] "namdmcmc.tcl"]
 
 namespace eval ::cphTitrator {
+    #   The core information of the cphTitrator is the moveDict, which stores
+    # all possible neMD/MC moves involving one or more titratable residues. The
+    # keys are of the form:
+    #
+    # segresidname[/segresidname..[/segresidname]]
+    #
+    # in analogy to the keys of ::cphSystem::residueDict. Note that a "/" is
+    # used as a syntactic divider between residues. Only _specific_ move
+    # information is kept in moveDict. Another dict, defaultMoveParams is used
+    # to store default information (similar to ::cphSystem::residueDefDict) and
+    # is used to avoid saving copies of redundant information.
+    #
     namespace import ::cphSystem::*
-    variable pH
-    variable MoveSet [dict create]
-    variable WeightSum 0.0
-    variable maxAttempts 0
+
+    variable moveDict [dict create]
+    variable defaultMoveParams [dict create]
+    variable maxAttempts
 
     namespace export cphSystem
-    # These are defined here
-    namespace export buildTitrator proposeMove getMoveSet\
-           accumulateAcceptanceRate
+    namespace export cphTitrator
 }
 
 # =============================================================================
-# "Constructor" Routines 
+# cphTitrator Interface function
 # =============================================================================
-proc ::cphTitrator::buildTitrator {systempH moveInfo} {
-    variable ::cphTitrator::pH $systempH
-    variable ::cphTitrator::maxAttempts
-    variable ::cphTitrator::MoveSet
-    set maxProposalAttempts [dict get $moveInfo maxProposalAttempts]
-    dict unset moveInfo maxProposalAttempts
-    if {$maxProposalAttempts > 0} {
-        set maxAttempts [expr {int($maxProposalAttempts)}]
-    } else {
-        set maxAttempts [cphSystem get numresidues]
+# ::cphTitrator::cphTitrator
+#
+# This is the only exported function from the cphTitrator namespace and
+# provides a complete interface. The first argument is always an action (e.g.
+# "get" or "set") and the rest are variable arguments.
+#
+# action           description
+# ---------------- -----------
+# get              see cphTitratorGet 
+# set              see cphTitratorSet
+# build            see buildTitrator
+# propose          select and propose a move, return switch information
+# accumulateStats  accumulate mean acceptance rate for the given move
+#
+# Some notes about Tcl switch statements for non-guru types like myself:
+#
+#   These can appear a bit delicate to those who are not familiar with their
+# quirks (compared to C). For one, badly placed comments can break them. Also
+# note that a switch returns the expression it first matches and does NOT use
+# break statements.
+#
+proc ::cphTitrator::cphTitrator {action args} {
+    return [switch -nocase -- $action {
+        get {
+            cphTitratorGet {*}$args
+        }
+        set {
+            cphTitratorSet {*}$args
+        }
+        build {
+            buildTitrator {*}$args 
+        }
+        propose { ;# sample the move set and return the selected move info.
+            proposeMove {*}$args
+        }
+        accumulateStats { ;# accumulate MC statistics
+            accumulateAcceptanceRate {*}$args
+        }
+        default {
+            abort "Invalid cphTitrator action $action"
+        }
+    }]
+}
+
+# ::cphTitrator::validateMoveLabel
+#
+# Check that a given move label is valid. That is:
+#
+# 1) is it in the correct segresidname[/segresidname ...] format?
+# 2) are the component segresidnames valid?
+#
+# A non-zero error code and msg are returned for each condition. Zero is
+# returned for a valid segresidname.
+#
+# NB! This should only be called _after_ buildSystem.
+#
+proc ::cphTitrator::validateMoveLabel {moveLabel} {
+    set segresidnameList [split $moveLabel "/"]
+    if {![llength $segresidnameList]} {
+        print "Move labels should be of the form segresidname/segresidname ...]"
+        return -1
     }
-    #   The default move set is to switch each residue independently (i.e. the
-    # move label is just the segresidname for that residue). Settings from the
-    # restart and config files are discarded as they get used so that only
-    # unused options remain in the moveInfo dict when this is finished (also
-    # the default options, which must also be discarded.
-    #
-    foreach segresidname [cphSystem get segresidnames] {
-        dict set MoveSet $segresidname proposalCmd\
-                "proposeResidueMove $segresidname"
-        dict set MoveSet $segresidname segresidnameList $segresidname
-        dict set MoveSet $segresidname weight\
-                [dictPopOrDefault moveInfo $segresidname weight]
-        dict set MoveSet $segresidname numsteps\
-                [dictPopOrDefault moveInfo $segresidname numsteps]
-        dict set MoveSet $segresidname successes 0
-        dict set MoveSet $segresidname attempts 0
-        if {[dict exists $moveInfo $segresidname]\
-            && ![dict size [dict get $moveInfo $segresidname]]} {
-            dict unset moveInfo $segresidname
-        }
-    }
-    # Build proton transfer moves
-    foreach moveLabel [dict get $moveInfo ptransfer] {
-        lassign [split $moveLabel "/"] segresidname1 segresidname2
-        dict set MoveSet $moveLabel proposalCmd\
-                "proposeProtonTransferMove $moveLabel" 
-        dict set MoveSet $moveLabel segresidnameList\
-                [list $segresidname1 $segresidname2]
-        dict set MoveSet $moveLabel weight\
-                [dictPopOrDefault moveInfo $moveLabel weight]
-        dict set MoveSet $moveLabel numsteps\
-                [dictPopOrDefault moveInfo $moveLabel numsteps]
-        dict set MoveSet $moveLabel successes 0
-        dict set MoveSet $moveLabel attempts 0
-        if {[dict exists $moveInfo $moveLabel]\
-            && ![dict size [dict get $moveInfo $moveLabel]]} {
-            dict unset moveInfo $moveLabel
-        }
-    }
-    dict unset moveInfo ptransfer
-    ###
-    computeWeightSum
-    return $moveInfo
-}
-
-# This is a kludgy shorthand for using dicts that contain setting values. If
-# the dict contains the move label, then any specific settings for that move
-# are assumed to be stored as a nested dict with that label as key. All values
-# are removed after being returned. If no such label or setting is present,
-# then a default setting is assumed to exist with just the setting as the key.
-#
-proc ::cphTitrator::dictPopOrDefault {settingsDict label setting} {
-    upvar 1 $settingsDict SettingsDict
-    if {[dict exists $SettingsDict $label $setting]} {
-        set value [dict get $SettingsDict $label $setting]
-        dict unset SettingsDict $label $setting
-    } else {
-       set value [dict get $SettingsDict default $setting] 
+    foreach segresidname $segresidnameList { 
+        if {[cphSystem validate $segresidname]} {
+            print "Bad segresidname $segresidname in move label $moveLabel!"
+            return -2
+        }
     }
-    return $value
+    return 0
 }
 
-# ::cphTitrator::computeWeightSum
-#
-# Compute, store, and return the sum of weights for the MC move set.
+# Pop the value from a dict with the corresponding key. If no key exists,
+# return the the default value instead.
 #
-proc ::cphTitrator::computeWeightSum {} {
-    lassign [getMoveWeightsAndSum] weights tmp
-    set ::cphTitrator::WeightSum [expr {lsum($weights)}]
-    return
+proc ::cphTitrator::dictPopOrDefault {myDict key {defaultValue {}}} {
+    upvar 1 $myDict MyDict
+    if {[dict exists $MyDict $key]} {
+        set value [dict get $MyDict $key]
+        dict unset MyDict $key
+    }
+    if {![info exists value]} {
+        set value $defaultValue
+    }
+    return $value
 }
 
 # =============================================================================
 # Proposal Routines
 # =============================================================================
+# ::cphTitrator::proposeMove
+#
+# Propose a move from all possible moves in the set. This is done by directly
+# sampling the probability mass function of the weighted moves. Once this is
+# done, test the inherent probability of that move given the imposed pH. If
+# rejected, try again up to the given number of "maxAttempts". Always return
+# the information necessary to perform a switch for the given proposal.
+#
+proc ::cphTitrator::proposeMove {pH} {
+    variable ::cphTitrator::maxAttempts
+    variable ::cphTitrator::moveDict
+    set weights [cphTitrator get weight]
+
+    set accept 0
+    set attemptsThisCycle 0
+    while {!$accept && $attemptsThisCycle < $maxAttempts} {
+        incr attemptsThisCycle
+        # Clear the previous trial if it was rejected.
+        if {$attemptsThisCycle > 1} {
+             cphSystem update $accept $segresidnameList
+        }
+
+        lassign [choice [cphTitrator get moveLabels] $weights] moveLabel
+        set proposalCmd [cphTitrator get proposalCmd $moveLabel]
+        set accept [eval $proposalCmd $pH]
+        set segresidnameList [cphTitrator get segresidnameList $moveLabel]
+    }
+    set numsteps [cphTitrator get numsteps $moveLabel]
+    return [list $accept $numsteps $segresidnameList $attemptsThisCycle]
+}
+
 # ::cphTitrator::proposeResidueMove
 #
 # Propose a move involving a single residue via Metropolized independence
@@ -127,17 +170,12 @@ proc ::cphTitrator::computeWeightSum {} {
 # trivial, it's just the other state. For three or more states, this will
 # naturally prefer tautomerization at extreme pKa values.
 #
-proc ::cphTitrator::proposeResidueMove {segresidname} {
-   variable ::cphTitrator::pH
-   lassign [cphSystem compute inherentWeights $pH $segresidname] weights states 
-   # NB: The remaining weights normalize to pic = (1 - pi), NOT 1.
-   #     If you want to see the weights as probabilities, each weight should
-   #     be divided by pic (this is implicit inside weightedChoice).
-   #
-   set pic [expr {1. - [lindex $weights 0]}]
-   set j [weightedChoice [lrange $weights 1 end] $pic]
-   set pj [lindex $weights [expr {$j+1}]] ;# note the index frame due to lrange
-   set du [expr {log((1. - $pj) / $pic)}]
+proc ::cphTitrator::proposeResidueMove {segresidname pH} {
+   lassign [cphSystem compute inherentWeights $pH $segresidname] weights states
+
+   set pi [lindex $weights 0]
+   lassign [choice [lrange $weights 1 end] [lrange $weights 1 end]] pj j
+   set du [expr {log((1. - $pj) / (1. - $pi))}]
    set accept [metropolisAcceptance $du]
    if {$accept} {
        cphSystem set trialState $segresidname [lindex $states $j]
@@ -145,84 +183,329 @@ proc ::cphTitrator::proposeResidueMove {segresidname} {
    return $accept
 }
 
-proc ::cphTitrator::proposeProtonTransferMove {moveLabel} {
-    variable ::cphTitrator::pH
+# ::cphTitrator::proposeProtonTransferMove
+#
+# Propose a move involving two residues in which a proton is moved from one to
+# the other. Return True if such a move was proposed, accepted, and stored.
+#
+# Note that this automatically fails if no proton can be transferred.
+#
+proc ::cphTitrator::proposeProtonTransferMove {moveLabel pH} {
     lassign [split $moveLabel "/"] segresidname1 segresidname2
-    set errCode\
-        [::cphSystem::proposeProtonTransfer $segresidname1 $segresidname2]
-    if {$errCode > 0} {
+
+    if {[cphSystem propose protonTransfer $segresidname1 $segresidname2]} {
+        set accept 0
+    } else {
         set du1 [cphSystem compute inherent $pH $segresidname1]
         set du2 [cphSystem compute inherent $pH $segresidname2]
         set accept [metropolisAcceptance [expr {$du1 + $du2}]]
-    } else {
-        set accept 0
     }
     return $accept
 }
 
-proc ::cphTitrator::proposeMove {} {
+# =============================================================================
+# "Constructor" Routines 
+# =============================================================================
+# ::cphTitrator::buildTitrator
+#
+# Construct the titrator given a set of MC move information.
+#
+proc ::cphTitrator::buildTitrator {moveInfo} {
+    variable ::cphTitrator::moveDict
+    variable ::cphTitrator::defaultMoveParams
     variable ::cphTitrator::maxAttempts
-    variable ::cphTitrator::MoveSet
-    lassign [getMoveWeightsAndSum] weights weightSum
 
-    set accept 0
-    set nattempts 0
-    while {!$accept && $nattempts < $maxAttempts} {
-        incr nattempts
-        # Clear the previous trial if it was rejected.
-        if {$nattempts > 1} {
-             cphSystem update $accept $segresidnameList
+    # 1) Pop off global and default move parameters.
+    #
+    set maxAttempts [dictPopOrDefault moveInfo maxProposalAttempts 0]
+    if {$maxAttempts <= 0} {
+        # This is a reasonable default for most systems.
+        set maxAttempts [cphSystem get numresidues]
+    }
+    set defaultMoveParams [dict merge $defaultMoveParams\
+                                      [dictPopOrDefault moveInfo default]]
+    if {![dict exists $defaultMoveParams default weight]} {
+        cphTitrator set weight default 1.0
+    }
+    if {![dict exists $defaultMoveParams default protonTransfer]} {
+        cphTitrator set protonTransfer default 0
+    }
+
+    # 2) Validate all remaining keys - these better be move labels!
+    #
+    dict for {moveLabel data} $moveInfo {
+        if {[validateMoveLabel $moveLabel]} {
+            abort
+        }
+        dict for {attr value} $data {
+            cphTitrator set $attr $moveLabel $value
         }
-        set index [weightedChoice $weights $weightSum]
-        set moveLabel [lindex [dict keys $MoveSet] $index]
-        set proposalCmd [dict get $MoveSet $moveLabel proposalCmd]
-        set accept [eval $proposalCmd]
-        set numsteps [dict get $MoveSet $moveLabel numsteps]
-        set segresidnameList [dict get $MoveSet $moveLabel segresidnameList]
     }
-    return [list $accept $numsteps $segresidnameList $nattempts]
+    set moveDict [dict merge $moveDict $moveInfo]
+
+    # The default move set is to titrate each residue independently.
+    foreach segresidname [cphSystem get segresidnames] {
+        set moveLabel $segresidname
+        cphTitrator set proposalCmd $moveLabel\
+                "proposeResidueMove $segresidname"
+    }
+    # Build proton transfer moves, if present.
+    dict for {moveLabel data} $moveInfo {
+        if {![dict exists $data protonTransfer]
+            || ![dict get $data protonTransfer]} {
+            continue
+        }
+        if {[llength [split $moveLabel "/"]] != 2} {
+            abort "Proton transfer moves must have exactly 2 segresidnames!"
+        }
+        cphTitrator set proposalCmd $moveLabel\
+                "proposeProtonTransferMove $moveLabel"
+    }
+
+    # Initialize statistics for any moves that are new.
+    dict for {moveLabel data} $moveDict {
+        if {![dict exists $data attempts]} {
+            cphTitrator set successes $moveLabel 0
+            cphTitrator set attempts $moveLabel 0
+        }
+    }
+
+    return
 }
 
 # =============================================================================
-# Routines for tracking and reporting MC statistics
+# Getter Routines
 # =============================================================================
-# ::cphTitrator::accumulateAcceptanceRate
+# ::Titrator::cphTitratorGet
 #
-# Accumulate statistics for the given move.
+# Getter for move attributes, called as:
 #
-proc ::cphTitrator::accumulateAcceptanceRate {accept moveLabel} {
-    variable ::cphTitrator::MoveSet
-    set moveLabel [join $moveLabel "/"]
-    # Alas, dict incr does not support nested keys.
-    set successes [dict get $MoveSet $moveLabel successes]
-    set attempts [dict get $MoveSet $moveLabel attempts]
-    incr successes $accept
-    incr attempts
-    dict set MoveSet $moveLabel successes $successes
-    dict set MoveSet $moveLabel attempts $attempts 
-    return
+#   <attribute> [<moveLabel>]
+#
+# With some specialized exceptions, <attribute> is either the name of a system
+# attribute (usually a list) or else a move attribute.
+#
+# system attributes  description
+# -----------------  -----------
+# moveLabels         list of all the moveLabels
+# maxAttempts        max # of attempts at move proposals
+#
+# move attributes  description
+# ---------------  -----------
+# numsteps         number of steps in a switch after successful proposal
+# weight           weight for consideration in move selection
+# protonTransfer   if 2 residues are involved, can they proton transfer?
+# successes        the number of successful _switches_ for this move
+# attempts         the number of attempted _switches_ for this move
+# proposalCmd      the command needed to set up trial states
+# segresidnameList list of the residues involved (NB: this is just shorthand
+#                  for [split $moveLabel "/"])
+#
+# specialized calls  description
+# -----------------  -----------
+# nodefaults         Return a dictionary of all current move information, but
+#                    suppress entries where a move has the default value. This
+#                    gives a minimal, but full, description of the move set.
+#
+# For now this is _much_ simpler than the analogous cphSystemGet. In the future
+# it may be useful to generalize this though.
+#
+proc ::cphTitrator::cphTitratorGet {attr {moveLabel {}}} {
+    variable ::cphTitrator::moveDict    
+    variable ::cphTitrator::maxAttempts
+
+    set getAll [expr {![llength $moveLabel]}]
+    if {!$getAll && ![dict exists $moveDict $moveLabel]
+        && ![string match -nocase $moveLabel default]} {
+        abort "cphTitratorGet: Invalid moveLabel $moveLabel"
+    }
+
+    return [switch -nocase -- $attr {
+        moveLabels {
+            dict keys $moveDict
+        }
+        maxAttempts {
+            expr {$maxAttempts}
+        }
+        numsteps -
+        weight -
+        protonTransfer -
+        successes -
+        attempts -
+        proposalCmd {
+            if {$getAll} {
+                getAllMoveAttr $attr 
+            } else {
+                getMoveAttr $attr $moveLabel
+            }
+        }
+        segresidnameList {
+            if {$getAll} {
+                getAllSegresidnameLists
+            } else {
+                split $moveLabel "/"
+            }
+        }
+        nodefaults {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                getMoveNoDefaults $moveLabel
+            }
+        }
+        default {
+            abort "cphTitratorGet: Invalid attribute $attr"
+        }
+    }]
+}
+
+proc ::cphTitrator::cannotGetAll {attr} {
+    abort "cphTitratorGet: Cannot get all $attr - must select a moveLabel"
+    return -1
+}
+
+# ---------------------------
+# Getters for move attributes
+# ---------------------------
+# Return the default value if no specific value was set.
+proc ::cphTitrator::getMoveAttr {attr moveLabel} {
+    variable ::cphTitrator::moveDict
+    variable ::cphTitrator::defaultMoveParams
+
+    if {[dict exists $moveDict $moveLabel $attr]} {
+        return [dict get $moveDict $moveLabel $attr]
+    } elseif {[dict exists $defaultMoveParams $attr]} {
+        return [dict get $defaultMoveParams $attr]
+    }
+    abort "cphTitratorGet: Error getting $attr for move $moveLabel"
+    return -1
+}
+
+proc ::cphTitrator::getAllMoveAttr {attr} {
+    set retList [list]
+    foreach moveLabel [cphTitrator get moveLabels] {
+        lappend retList [cphTitrator get $attr $moveLabel]
+    }
+    return $retList
+}
+
+proc ::cphTitrator::getAllSegresidnameLists {} {
+    set retList [list]
+    foreach moveLabel [cphTitrator get moveLabels] {
+        lappend retList [split $moveLabel "/"]
+    }
+    return $retList
+}
+
+# -----------------------------
+# Special getters for archiving
+# -----------------------------
+# Return the dict for a given move label, but strip out default values.
+proc ::cphTitrator::getMoveNoDefaults {moveLabel} {
+    set defaultNumsteps [cphTitrator get numsteps default]
+    set defaultWeight [cphTitrator get weight default]
+    set defaultPT [cphTitrator get protonTransfer default]
+    set numsteps [cphTitrator get numsteps $moveLabel]
+    set weight [cphTitrator get weight $moveLabel]
+    set PT [cphTitrator get protonTransfer $moveLabel]
+
+    set retDict [dict create]
+    if {$numsteps != $defaultNumsteps} {
+        dict set retDict numsteps $numsteps
+    }
+    if {$weight != $defaultWeight} {
+        dict set retDict weight $weight
+    }
+    if {$PT != $defaultPT} {
+        dict set retDict protonTransfer $PT
+    }
+    dict set retDict attempts [cphTitrator get attempts $moveLabel]
+    dict set retDict successes [cphTitrator get successes $moveLabel]
+
+    return $retDict
 }
 
 # =============================================================================
-# Titrator Getter Routines
-#
-# These all return a list of values describing each move in the titrator.
+# Setter Routines
 # =============================================================================
-proc ::cphTitrator::getMoveSet {} {
-    return $::cphTitrator::MoveSet
+# ::cphTitrator::cphTitratorSet
+#
+# Setters for move attributes, called as:
+#
+#  <attribute> <moveLabel> <value>
+#
+# <attribute> is the name of a move attribute.
+#
+# move attributes  description
+# ---------------  -----------
+# numsteps         number of steps in a switch after successful proposal
+# weight           weight for consideration in move selection
+# protonTransfer   if 2 residues are involved, can they proton transfer?
+# successes        the number of successful _switches_ for this move
+# attempts         the number of attempted _switches_ for this move
+# proposalCmd      the command needed to set up trial states
+#
+proc ::cphTitrator::cphTitratorSet {attr moveLabel value} {
+    variable ::cphTitrator::moveDict
+    variable ::cphTitrator::defaultMoveParams
+
+    if {[string match -nocase $moveLabel default]} {
+        set setDefault 1
+    } else {
+        if {[validateMoveLabel $moveLabel]} {
+            abort
+        }
+        set setDefault 0
+    }
+
+    return [switch -nocase -- $attr {
+        numsteps -
+        successes -
+        attempts -
+        protonTransfer { ;# Require integer or boolean argument.
+            set value [expr {int($value)}]
+            if {$setDefault} {
+                dict set defaultMoveParams $attr $value
+            } else {
+                dict set moveDict $moveLabel $attr $value
+            }
+            expr {$value}
+        }
+        weight { ;# Require float argument.
+            set value [expr {1.*$value}]
+            if {$setDefault} {
+                dict set defaultMoveParams $attr $value
+            } else {
+                dict set moveDict $moveLabel $attr $value
+            }
+            expr {$value}
+        }
+        proposalCmd { ;# Argument is string or pure list.
+            dict set moveDict $moveLabel $attr $value
+            expr {$value}
+        }
+        default {
+            abort "cphTitratorSet: Invalid attribute $attr"
+            expr {-1}
+        }
+    }]
 }
 
-# ::cphTitrator::getMoveWeightsAndSum
+# =============================================================================
+# Routines for tracking and reporting MC statistics
+# =============================================================================
+# ::cphTitrator::accumulateAcceptanceRate
 #
-# Return the list of MC move weights along with their (precomputed) sum.
+# Accumulate statistics for the given move.
 #
-proc ::cphTitrator::getMoveWeightsAndSum {} {
-    variable ::cphTitrator::MoveSet
-    variable ::cphTitrator::WeightSum
-    set weights [list]
-    dict for {moveLabel data} $MoveSet {
-        lappend weights [dict get $data weight]
-    }
-    return [list $weights $WeightSum]
+proc ::cphTitrator::accumulateAcceptanceRate {accept moveLabel} {
+    # Alas, dict incr does not support nested keys.
+    set successes [cphTitrator get successes $moveLabel]
+    set attempts [cphTitrator get attempts $moveLabel]
+    incr successes $accept
+    incr attempts
+    cphTitrator set successes $moveLabel $successes
+    cphTitrator set attempts $moveLabel $attempts
+    return
 }
 
index c438b63..9d79187 100644 (file)
@@ -56,63 +56,116 @@ namespace eval ::cphSystem {
 # dealchemifypsf   remove patches so the system is no longer alchemical
 # initializeState  set the initial state of a residue
 #
+# Some notes about Tcl switch statements for non-guru types like myself:
+#
+#   These can appear a bit delicate to those who are not familiar with their
+# quirks (compared to C). For one, badly placed comments can break them. Also
+# note that a switch returns the expression it first matches and does NOT use
+# break statements.
+#
 proc ::cphSystem::cphSystem {action args} {
-    if {[string match -nocase $action get]} {
-        # Getters
-        return [cphSystemGet {*}$args] 
-    } elseif {[string match -nocase $action set]} {
-        # Setters
-        return [cphSystemSet {*}$args]
-    } elseif {[string match -nocase $action build]} {
-        # System building - determine system composition and definitions
-        return [buildSystem {*}$args]
-    } elseif {[string match -nocase $action initialize]} {
-        # System initialization - establish initial states and parameters
-        return [initializeSystem {*}$args]
-    } elseif {[string match -nocase $action propose]} {
-        # Proposal - Set new trial state - return 0 if none is found
-        set type [lindex $args 0]
-        set newArgs [lrange $args 1 end]
-        if {[string match -nocase $type titration]} {
-            return [proposeResidueTitration {*}$newArgs]
-        } elseif {[string match -nocase $type tautomerization]} {
-            return [proposeResidueTautomerization {*}$newArgs]
-        } else {
-            abort "Invalid proposal type $type"
+    return [switch -nocase -- $action {
+        get {
+            cphSystemGet {*}$args
         }
-    } elseif {[string match -nocase $action compute]} {
-        # Compute acceptance energy (inherent or switch correction)
-        set type [lindex $args 0]
-        set newArgs [lrange $args 1 end]
-        if {[string match -nocase $type inherent]} {
-            return [computeInherentAcceptance {*}$newArgs]
-        } elseif {[string match -nocase $type inherentWeights]} {
-            return [computeInherentNormedWeights {*}$newArgs]
-        } elseif {[string match -nocase $type switch]} {
-            return [computeSwitchAcceptance {*}$newArgs]
-        } else {
-            abort "Invalid energy type $type"
+        set {
+            cphSystemSet {*}$args
+        }
+        build {
+            buildSystem {*}$args
+        }
+        initialize {
+            initializeSystem {*}$args
+        }
+        propose { ;# return 0 if no valid proposal is found
+            set type [lindex $args 0]
+            set newArgs [lrange $args 1 end]
+            switch -nocase -- $type {
+                titration {
+                    proposeResidueTitration {*}$newArgs
+                }
+                tautomerization {
+                    proposeResidueTautomerization {*}$newArgs
+                }
+                protonTransfer {
+                    proposeProtonTransfer {*}$newArgs
+                }
+                default {
+                    abort "Invalid proposal type $type"
+                }
+            }
+        }
+        compute { ;# acceptance energy (inherent or switch correction)
+            set type [lindex $args 0]
+            set newArgs [lrange $args 1 end]
+            switch -nocase -- $type {
+                inherent {
+                    computeInherentAcceptance {*}$newArgs
+                }
+                inherentWeights {
+                    computeInherentNormedWeights {*}$newArgs
+                }
+                switch {
+                    computeSwitchAcceptance {*}$newArgs
+                }
+                default {
+                    abort "Invalid energy type $type"
+                }
+            }
+        }
+        update { ;# Update states based on an MC result
+            updateStates {*}$args
         }
-    } elseif {[string match -nocase $action update]} {
-       # Update states based on an MC result
-       return [updateStates {*}$args]
-    } elseif {[string match -nocase $action alchemifypsf]} {
-       # Alchemify the PSF (in memory) based on trial states
-       return [alchemifyResidue {*}$args]
-    } elseif {[string match -nocase $action dealchemifypsf]} {
-       # Dealchemify the PSF (in memory)
-       return [dealchemifyResidue {*}$args]
-    } elseif {[string match -nocase $action initializeState]} {
-       # Assign a state
-       set method [lindex $args 0]
-       if {[string match -nocase $method random]} {
-           return [randomizeState {*}[lrange $args 1 end]]
-       } else {
-           return [assignState {*}$args]
+        alchemifypsf { ;#Alchemify the PSF (in memory) based on trial states
+            alchemifyResidue {*}$args
+        }
+        dealchemifypsf { ;# Dealchemify the PSF (in memory)
+            dealchemifyResidue {*}$args
+        }
+        initializeState { ;# Assign a state
+            set method [lindex $args 0]
+            switch -nocase -- $method {
+                random {
+                    randomizeState {*}[lrange $args 1 end]
+                }
+                default {
+                    assignState {*}$args
+                }
+            }
        }
-    } else {
-        abort "Invalid cphSystem action $action."
+       validate { ;# Validate a segresidname
+           validateSegresidname {*}$args
+       }
+       default {
+           abort "Invalid cphSystem action $action."
+       }
+    }]
+}
+
+# ::cphSystem::validateSegresidname
+#
+# Check that a given segresidname is valid. That is:
+#
+# 1) is it in the correct <segid>:<resid>:<resname> format?
+# 2) does it correspond to a known titratable residue?
+#
+# A non-zero error code is returned for each condition.
+#
+# NB! This should only be called _after_ buildSystem.
+#
+proc ::cphSystem::validateSegresidname {segresidname} {
+    lassign [split $segresidname ":"] segid resid resname
+    if {![info exists segid] || ![info exists resid]
+        || ![info exists resname]} {
+        print "segresidname selections must be of form\
+                 <segid>:<resid>:<resname>!"
+        return -1
     }
+    if {$segresidname ni [cphSystem get segresidnames]} {
+        print "Invalid segresidname $segresidname!"
+        return -2
+    }
+    return 0
 }
 
 # =============================================================================
@@ -153,7 +206,7 @@ proc ::cphSystem::proposeResidueTitration {segresidname} {
     set resname [cphSystem get resname $segresidname]
     set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
     while {true} {
-        set state [choice $possibleStates]
+        lassign [choice $possibleStates] state
         set occ [dict get $resDefDict $resname states $state]
         set numProtonsTrial [expr {lsum($occ)}]
         if {$numProtons != $numProtonsTrial} {
@@ -170,11 +223,11 @@ proc ::cphSystem::proposeResidueTitration {segresidname} {
 # Propose a new trial state requiring a proton transfer - i.e. a movement of a
 # proton from residue to another.
 #
-# If both residues in the pair have/do not currently have a proton, then a
-# transfer cannot occur and this returns false.
+# This returns true if no proton transfer or else a negative error code.
 #
 proc ::cphSystem::proposeProtonTransfer {segresidname1 segresidname2} {
     variable ::cphSystem::resDefDict
+
     set numProtons1 [expr {lsum([cphSystem get occupancy $segresidname1])}]
     set possibleStates1 [cphSystem get trialStateList $segresidname1]
     set resname1 [cphSystem get resname $segresidname1]
@@ -185,27 +238,26 @@ proc ::cphSystem::proposeProtonTransfer {segresidname1 segresidname2} {
 
     if {$dn == 0.0} {     
         # No proton transfer is possible.
-        return 0
+        return 1
     } 
     # transfer from 1 to 2
     while {true} {
-        set state1 [choice $possibleStates1]
+        lassign [choice $possibleStates1] state1
         set occ1 [dict get $resDefDict $resname1 states $state1]
         set numProtonsTrial1 [expr {lsum($occ1)}]
-        set state2 [choice $possibleStates2]
+        lassign [choice $possibleStates2] state2
         set occ2 [dict get $resDefDict $resname2 states $state2]
         set numProtonsTrial2 [expr {lsum($occ2)}]
         if {$dn == [expr {$numProtonsTrial2 - $numProtonsTrial1}]} {
             cphSystem set trialState $segresidname1 $state1
             cphSystem set trialState $segresidname2 $state2
-            return 1
+            return 0
         }
     }
     # This is an error and should never happen.
     return -1
 }
 
-
 # ::cphSystem::proposeResidueTautomerization
 #
 # Propose a new trial state requiring a tautomerization - i.e. the number of 
@@ -239,7 +291,7 @@ proc ::cphSystem::proposeResidueTautomerization {segresidname {maxAttempts 10}}
     set resname [cphSystem get resname $segresidname]
     set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
     while {true} {
-        set state [choice $possibleStates]
+        lassign [choice $possibleStates] state
         set occ [dict get $resDefDict $resname states $state]
         set numProtonsTrial [expr {lsum($occ)}]
         if {$numProtonsTrial == $numProtons} {
@@ -531,26 +583,18 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
             set segresidname [format "%s:%s" $segresid $resname]
             # Perform aliasing to fix residues with multiple names.
             dict for {realName altNames} $resAliases {
-                if {[lsearch -nocase $altNames $resname] < 0} {
+                if {$resname ni $altNames} {
                     continue
                 }
                 print "aliasing $segresidname to $realName"
                 psfset resname $segid $resid $realName
                 set resname $realName
-                set segresid [format "%s:%s" $segresid $resname]
-            }
-            set resIsDefined 0
-            if {[lsearch -nocase $definedResidues $resname] >= 0} {
-                set resIsDefined 1
-            }
-            set resIsSubtyped 0
-            if {[lsearch -nocase $definedSubResidues $resname] >= 0} {
-                set resIsSubtyped 1
-            }
-            set resIsExcluded 0
-            if {[lsearch -nocase $segresExcls $segresidname] >= 0} {
-                set resIsExcluded 1
+                set segresidname [format "%s:%s" $segresid $resname]
             }
+
+            set resIsDefined [expr {$resname in $definedResidues}]          
+            set resIsSubtyped [expr {$resname in $definedSubResidues}]
+            set resIsExcluded [expr {$segresidname in $segresExcls}]
             # Break here if nothing to do.
             if {(!$resIsDefined && !$resIsSubtyped) || $resIsExcluded} {
                 continue
@@ -562,7 +606,6 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
             }
 
             if {!$resIsSubtyped} continue
-
             # The conditions for subtyping may still not be met...
             set atoms [segment atoms $segid $resid]
             foreach subresname [dict get $subresDefDict $resname] {
@@ -578,14 +621,14 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
 
                 set allSubatomsExist 1
                 foreach atom $subatoms {
-                    if {[lsearch -nocase $atoms $atom] < 0} {
+                    if {$atom ni $atoms} {
                         set allSubatomsExist 0
                         break
                     }
                 }
                 set allNotSubatomsDoNotExist 1
-                foreach atom $notsubatoms { 
-                    if {[lsearch -nocase $atoms $atom] >= 0} {
+                foreach atom $notsubatoms {
+                    if {$atom in $atoms} {
                         set allNotSubatomsDoNotExist 0
                         break
                     }
@@ -596,7 +639,6 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
                     dict set resDict $segresidname resname $subresname
                 }
             }
-
         }
     }
     return $resDict
@@ -654,10 +696,10 @@ proc ::cphSystem::checkResDefs {resDefs} {
             #
             set statePatch [format "%s%s" $resname $state]
             set hybridPatch [format "%sH%s" $resname $state]
-            if {[lsearch -nocase [topology patches] $statePatch] < 0} {
+            if {$statePatch ni [topology patches]} {
                 abort "No patch definition in RTFs for $statePatch!"
             }
-            if {[lsearch -nocase [topology patches] $hybridPatch] < 0} {
+            if {$hybridPatch ni [topology patches]} {
                 abort "No patch definition in RTFs for $hybridPatch!"
             }
         }
@@ -723,7 +765,7 @@ proc ::cphSystem::assignState {segresidname state} {
 proc ::cphSystem::randomizeState {segresidname pH} {
     while {true} {
         set states [cphSystem get stateList $segresidname]
-        cphSystem set state $segresidname [choice $states]
+        cphSystem set state $segresidname [lindex [choice $states] 0]
         cphSystem propose titration $segresidname
         set du [cphSystem compute inherent $pH $segresidname]
         if {[metropolisAcceptance $du]} { 
@@ -786,124 +828,138 @@ proc ::cphSystem::cphSystemGet {attr {segresidname {}} args} {
     if {!$getAll && ![dict exists $resDict $segresidname]} {
         abort "cphSystemGet: Invalid segresidname $segresidname"
     }
-    # System attributes - any selection is invalid.
-    #
-    if {[string match -nocase $attr segresidnames]} {
-        return [dict keys $resDict]
-    } elseif {[string match -nocase $attr numresidues]} {
-        return [dict size $resDict]
-    } elseif {[string match -nocase $attr resdefs]} {
-        return [dict keys $resDefDict]
-    } elseif {[string match -nocase $attr numdefs]} {
-        return [dict size $resDefDict]
-    }
-    # Residue attributes - some of these need to be specially computed.
-    #
-    if {[lsearch -nocase {resname state trialState pKai} $attr] > -1} {
-        if {$getAll} {
-            return [getAllResAttr $attr]
-        } else {
-            return [getResAttr $attr $segresidname]
+
+    return [switch -nocase -- $attr {
+        segresidnames {
+            dict keys $resDict
         }
-    } 
-    if {[lsearch -nocase {dG pKa Tref} $attr] > -1} {
-        if {$getAll} {
-            return [getAllResDefAttr $attr]
-        } else {
-            return [getResDefAttr $attr $segresidname]
+        numresidues {
+            dict size $resDict
         }
-    } 
-    if {[string match -nocase $attr occupancy]} {
-        if {$getAll} {
-            return [getAllOccupancy]
-        } else {
-            return [getOccupancy $segresidname]
+        resdefs {
+            dict keys $resDefDict
         }
-    } 
-    if {[string match -nocase $attr trialOccupancy]} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            return [getTrialOccupancy $segresidname]
+        numdefs {
+            dict size $resDefDict
         }
-    } 
-    if {[string match -nocase $attr stateList]} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            set resname [getResAttr resname $segresidname]
-            return [dict keys [dict get $resDefDict $resname states]] 
+        resname -
+        state -
+        trialState -
+        pKai {
+            if {$getAll} {
+                getAllResAttr $attr
+            } else {
+                getResAttr $attr $segresidname
+            }
         }
-    }
-    if {[string match -nocase $attr trialStateList]} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            set resname [getResAttr resname $segresidname]
-            set state [getResAttr state $segresidname]
-            set states [dict keys [dict get $resDefDict $resname states]]
-            return [lsearch -all -inline -not -nocase $states $state]
+        dG -
+        pKa -
+        Tref {
+            if {$getAll} {
+                getAllResDefAttr $attr
+            } else {
+                getResDefAttr $attr $segresidname
+            }
         }
-    }
-    if {[lsearch -nocase {dGPair pKaPair pKaiPair} $attr] > -1} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            if {[lsearch -nocase {dGPair pKaPair} $attr] > -1} {
-                return [getResDefAttr $attr $segresidname]
+        occupancy {
+            if {$getAll} {
+                getAllOccupancy
             } else {
-                return [getResAttr $attr $segresidname]
+                getOccupancy $segresidname
             }
         }
-    }
-    if {[string match -nocase $attr statePatch]} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            set resname [getResAttr resname $segresidname]
-            set state [getResAttr state $segresidname]
-            return [format "%s%s" $resname $state]
+        trialOccupancy {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                getTrialOccupancy $segresidname
+            }
         }
-    }
-    if {[string match -nocase $attr hybridPatch]} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            set resname [getResAttr resname $segresidname]
-            set trialState [getResAttr trialState $segresidname]
-            return [format "%sH%s" $resname $trialState]
+        stateList {
+            if {$getAll} {
+                cannotGetAll $attr 
+            } else {
+                set resname [getResAttr resname $segresidname]
+                dict keys [dict get $resDefDict $resname states]
+            }
         }
-    }
-    if {[string match -nocase $attr alchAtomLists]} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            return [list [getResDefAttr l0atoms $segresidname]\
-                         [getResDefAttr l1atoms $segresidname]]
+        trialStateList {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                set resname [getResAttr resname $segresidname]
+                set state [getResAttr state $segresidname]
+                set states [dict keys [dict get $resDefDict $resname states]]
+                lsearch -all -inline -not -nocase $states $state
+            }
         }
-    }
-    if {[string match -nocase $attr alchBonds]} {
-        if {$getAll} {
-            cannotGetAll $attr
-        } else {
-            lassign $args k
-            set bondEntries [list]
-            foreach l0atom [getResDefAttr l0atoms $segresidname]\
-                    l1atom [getResDefAttr l1atoms $segresidname] {
-                lassign [split $segresidname ":"] segid resid
-                # Note that atomid indices start at one, not zero!
-                set i [expr {[segment atomid $segid $resid $l0atom] - 1}]
-                set j [expr {[segment atomid $segid $resid $l1atom] - 1}]
-                lappend bondEntries [format "bond %d %d %f %f" $i $j $k 0]
+        dGPair -
+        pKaPair {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                getResDefAttr $attr $segresidname
             }
-            return [join $bondEntries "\n"]
         }
-    }
-    abort "cphSystemGet: Invalid attribute $attr"
+        pKaiPair {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                getResAttr $attr $segresidname
+            }
+        }
+        statePatch {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                set resname [getResAttr resname $segresidname]
+                set state [getResAttr state $segresidname]
+                format "%s%s" $resname $state
+            }
+        }
+        hybridPatch {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                set resname [getResAttr resname $segresidname]
+                set trialState [getResAttr trialState $segresidname]
+                format "%sH%s" $resname $trialState
+            }
+        }
+        alchAtomLists {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                list [getResDefAttr l0atoms $segresidname]\
+                     [getResDefAttr l1atoms $segresidname]
+            }
+        }
+        alchBonds {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                lassign $args k
+                set bondEntries [list]
+                foreach l0atom [getResDefAttr l0atoms $segresidname]\
+                        l1atom [getResDefAttr l1atoms $segresidname] {
+                    lassign [split $segresidname ":"] segid resid
+                    # Note that atomid indices start at one, not zero!
+                    set i [expr {[segment atomid $segid $resid $l0atom] - 1}]
+                    set j [expr {[segment atomid $segid $resid $l1atom] - 1}]
+                    lappend bondEntries [format "bond %d %d %f %f" $i $j $k 0]
+                }
+                join $bondEntries "\n"
+            }
+        }
+        default {
+            abort "cphSystemGet: Invalid attribute $attr"
+        }
+    }]
 }
 
 proc ::cphSystem::cannotGetAll {attr} {
     abort "cphSystemGet: Cannot get all $attr - must select a segresidname"
+    return -1
 }
 
 # ------------------------------
@@ -1023,6 +1079,28 @@ proc ::cphSystem::index2flatindex {i j} {
     return [expr {($I*($I - 1) + 2*$J) / 2}]
 }
 
+# Convenience function for assignment of flat, off-diagonal, upper-triangular
+# matrices using 2d indices. Note that not all index combinations are valid!
+#
+# Example:
+#
+# % set myMatrix [lrepeat 3 0.0]
+# 0.0 0.0 0.0
+# % mset myMatrix 0 1 10.0
+# 10.0 0.0 0.0
+#
+proc ::cphSystem::mset {matrix i j value} {
+    upvar 1 $matrix Matrix
+    lset Matrix [index2flatindex $i $j] $value
+    return $Matrix
+}
+
+# Convenience function - inverse of mset.
+#
+proc ::cphSystem::mindex {matrix i j} {
+    return [lindex $matrix [index2flatindex $i $j]]
+}
+
 # =============================================================================
 # Setter Routines
 # =============================================================================
@@ -1043,29 +1121,39 @@ proc ::cphSystem::index2flatindex {i j} {
 proc ::cphSystem::cphSystemSet {attr segresidname value} {
     variable ::cphSystem::resDict
     variable ::cphSystem::resDefDict
+
     if {![dict exists $resDict $segresidname]} {
         abort "cphSystemSet: Invalid segresidname $segresidname"
     }
-    if {[lsearch -nocase {state trialState} $attr] > -1} {
-        set states [cphSystem get stateList $segresidname]
-        if {[lsearch -nocase $states $value] < 0} {
-            if {[string match -nocase $attr trialState] && ![llength $value]} {
-            
-            } else {
-                abort "Invalid state assignment $value for residue $segresidname"
+
+    return [switch -nocase -- $attr {
+        state -
+        trialState {
+            set states [cphSystem get stateList $segresidname]
+            if {$value ni $states} {
+                if {[string match -nocase $attr trialState]
+                    && ![llength $value]} {
+                } else {
+                    abort "Invalid state assignment $value for residue"\
+                            "$segresidname"
+                }
             }
+            dict set resDict $segresidname $attr $value
+            expr {$value}
         }
-        dict set resDict $segresidname $attr $value
-        return $value
-    }
-    if {[string match $attr pKai]} {
-        set resname [cphSystem get resname $segresidname]
-        set resDef [dict get $resDefDict $resname]
-        set pKaiMatrix [resDef2Matrix $resDef $value]
-        dict set resDict $segresidname pKai $value
-        dict set resDict $segresidname pKaiPair $pKaiMatrix 
-        return $pKaiMatrix
-    }
+        pKai {
+            set resname [cphSystem get resname $segresidname]
+            set resDef [dict get $resDefDict $resname]
+            set pKaiMatrix [resDef2Matrix $resDef $value]
+            dict set resDict $segresidname pKai $value
+            dict set resDict $segresidname pKaiPair $pKaiMatrix
+            expr {$pKaiMatrix}
+        }
+        default {
+            abort "cphSystemSet: Invalid attribute $attr"
+            expr {-1}
+        }
+    }]
 }
 
 # ::cphSystem::resDef2Matrix
@@ -1079,8 +1167,10 @@ proc ::cphSystem::cphSystemSet {attr segresidname value} {
 #
 # Arguments:
 # ----------
-# matrix : list
-#   
+# resDef : dict
+#   The residue definition to be assigned
+# data : list
+#   The _minimial_ attribute data needed to fill in the matrix
 #
 # Returns
 # -------
@@ -1109,7 +1199,7 @@ proc ::cphSystem::resDef2Matrix {resDef data} {
     # full set of states (i.e. the size is computed with maxStates, not
     # numStates). The row/column indices are computed as a left-right binary
     # conversion of the occupancy vector to a scalar (see occupancy2Index).
-    #  Just to make things more complicated, many systems have equivalent
+    #   Just to make things more complicated, many systems have equivalent
     # sites such that many fewer pairwise values are needed. For example, a
     # carboxylate formally has two sites and four states, one of which is
     # neglected. The two sites are also equivalent, so this system has, at max
@@ -1140,169 +1230,246 @@ proc ::cphSystem::resDef2Matrix {resDef data} {
     set maxStates [expr {int(pow(2, $numSites))}]
     set numPairs [expr {int($maxStates*($maxStates - 1) / 2)}]
     # Are any protonation counts missing?
-    set nprotonsExists [lrepeat [expr {$numSites+1}] 0]
+    set protCntExists [lrepeat [expr {$numSites+1}] 0]
     dict for {state occ} $states {
-        set nprotons [expr {lsum($occ)}]
-        for {set i 0} {$i <= $numSites} {incr i} {
-            if {$i == $nprotons} {
-                lset nprotonsExists $i 1
-            }
+        lset protCntExists [expr {int(lsum($occ))}] 1
+    }
+    set missingProtCnts [list]
+    set i 0
+    foreach exists $protCntExists {
+        if {!$exists} {
+            lappend missingProtCnts $i
         }
+        incr i
     }
 
     set Matrix [lrepeat $numPairs 0.0]
-    set numValues [llength $data]
+    set numValues [llength $data] 
     set errorCode 0
-    switch -- $numSites 1 {
-        switch -- $numValues 1 {
-            lset Matrix [index2flatindex 1 0] $data
-        } default {
-            set errorCode -1 
-        }
-    } 2 {
-        switch -- $numStates 3 {
-            if {![lindex $nprotonsExists 0]} {;# state (0,0) is deleted
-                switch -- $numValues 1 {
-                    lassign $data attr32
-                    set attr31 $attr32
-                } 2 {;# Ex. HIS
-                    lassign $data attr32 attr31
-                } default {
-                    set errorCode -1
-                }
-                lset Matrix [index2flatindex 3 2] $attr32
-                lset Matrix [index2flatindex 3 1] $attr31
-                lset Matrix [index2flatindex 2 1] [expr {$attr31 - $attr32}]
-            } elseif {![lindex $nprotonsExists 2]} {;# state (1,1) is deleted
-                switch -- $numValues 1 {;# Ex. carboxylates (ASP, GLU)
-                    lassign $data attr20
-                    set attr10 $attr20
-                } 2 {
-                    lassign $data attr20 attr10
-                } default {
-                    set errorCode -1
-                }
-                lset Matrix [index2flatindex 2 1] [expr {$attr20 - $attr10}]
-                lset Matrix [index2flatindex 2 0] $attr20
-                lset Matrix [index2flatindex 1 0] $attr10
-            } else {
-                # Why would state (0,1) or (1,0) be missing?
-                set errorCode -1
-            }
-        } 4 {
-            switch -- $numValues 4 {
-                lassign $data attr32 attr31 attr20 attr10
-            } default {
-                set errorCode -1
-            }
-            set attr30 [expr {$attr31 + $attr10}]
-            if {$attr30 != [expr {$attr32 + $attr20}]} {
-                set errorCode -2
-            }
-            lset Matrix [index2flatindex 3 2] $attr32
-            lset Matrix [index2flatindex 3 1] $attr31
-            lset Matrix [index2flatindex 3 0] $attr30
-            lset Matrix [index2flatindex 2 1] [expr {$attr20 - $attr10}]
-            lset Matrix [index2flatindex 2 0] $attr20
-            lset Matrix [index2flatindex 1 0] $attr10
-        } default {
-            set errorCode -1
-        }
-    } 3 {
-        switch -- $numStates 4 {
-            if {![lindex $nprotonsExists 0] && ![lindex $nprotonsExists 1]} {
-                # state (0,0,0) and 1 proton states missing
-                switch -- $numValues 1 {;# Ex. LYS
-                    lassign $data attr76
-                    set attr75 $attr76
-                    set attr73 $attr76
-                } default {
-                    set errorCode -1
-                }
-                set attr53 [expr {$attr73 - $attr75}]
-                set attr63 [expr {$attr73 - $attr76}]
-                set attr65 [expr {$attr75 - $attr76}]
-                if {$attr53 != [expr {$attr63 - $attr65}]
-                    || $attr63 != [expr {$attr53 - $attr65}]
-                    || $attr65 != [expr {$attr63 - $attr53}]} {
-                    set errorCode -2
-                }
-                lset Matrix [index2flatindex 7 6] $attr76
-                lset Matrix [index2flatindex 7 5] $attr75
-                lset Matrix [index2flatindex 7 3] $attr73
-                lset Matrix [index2flatindex 6 5] $attr65
-                lset Matrix [index2flatindex 6 3] $attr63
-                lset Matrix [index2flatindex 5 3] $attr53
-            } elseif {![lindex $nprotonsExists 2]
-                      && ![lindex $nprotonsExists 3]} {
-                # state (1, 1, 1) and 2 proton states missing
-                switch -- $numValues 1 {;# Ex. phosphate monoesters
-                    lassign $data attr10
-                    set attr20 $attr10
-                    set attr40 $attr10
-                } default {
+    switch -- $numSites {
+        1 {
+            set errorCode [expr {$numStates == 2 ? $errorCode : -1}]
+            set errorCode [expr {$numValues == 1 ? $errorCode : -1}]
+            lassign $data attr10
+            mset Matrix 1 0 $attr10
+        } ;# END 1 site, 1 state
+        2 {
+            switch -- $numStates {
+                3 {
+                    if {0 in $missingProtCnts} {
+                        # state (0,0) is deleted
+                        switch -- $numValues {
+                            1 { ;# Ex. H2O/H3O+
+                                lassign $data attr32
+                                mset Matrix 3 2 $attr32
+                                mset Matrix 3 1 $attr32
+                            }
+                            2 { ;# Ex. HIS
+                                lassign $data attr32 attr31
+                                mset Matrix 3 2 $attr32
+                                mset Matrix 3 1 $attr31
+                            }
+                            default {
+                                set errorCode -1
+                            }
+                        }
+                        if {$errorCode} break
+
+                        mset Matrix 2 1 [expr {[mindex $Matrix 3 1]\
+                                               - [mindex $Matrix 3 2]}]
+                        # END 2 sites, 3 states, no state 0 
+                    } elseif {2 in $missingProtCnts} {
+                        # state (1,1) is deleted
+                        switch -- $numValues {
+                            1 { ;# Ex. carboxylates (ASP, GLU)
+                                lassign $data attr20
+                                mset Matrix 2 0 $attr20
+                                mset Matrix 1 0 $attr20
+                            }
+                            2 { ;# Ex. asymmetric carboxylate?
+                                lassign $data attr20 attr10
+                                mset Matrix 2 0 $attr20
+                                mset Matrix 1 0 $attr10
+                            }
+                            default {
+                                set errorCode -1
+                            }
+                        }
+                        if {$errorCode} break
+
+                        mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
+                                               - [mindex $Matrix 1 0]}]
+                        # END 2 sites, 3 states, no state 3
+                    } else {
+                        # Why would state (0,1) or (1,0) be missing?
+                        set errorCode -1
+                    }
+                } ;# END 2 sites, 3 states
+                4 {
+                    switch -- $numValues {
+                        4 {
+                            lassign $data attr32 attr31 attr20 attr10
+                            mset Matrix 3 2 $attr32
+                            mset Matrix 3 1 $attr31
+                            mset Matrix 2 0 $attr20
+                            mset Matrix 1 0 $attr10
+                        }
+                        default {
+                            set errorCode -1
+                        }
+                    }
+                    if {$errorCode} break
+
+                    mset Matrix 3 0 [expr {[mindex $Matrix 3 1]\
+                                           + [mindex $Matrix 1 0]}]
+                    mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
+                                           - [mindex $Matrix 1 0]}]
+
+                    set alt30 [expr {[mindex $Matrix 3 2]\
+                                     + [mindex $Matrix 2 0]}]
+                    if {[mindex $Matrix 3 0] != $alt30} {
+                        set errorCode -2
+                    }
+                    # END 2 sites, 4 states, no states missing
+                } ;# END 2 sites, 4 states
+                default {
                     set errorCode -1
                 }
-                lset Matrix [index2flatindex 1 0] $attr10
-                lset Matrix [index2flatindex 2 0] $attr20
-                lset Matrix [index2flatindex 4 0] $attr40
-            } else {
-                set errorCode -1
             }
-        } 7 {
-            if {![lindex $nprotonsExists 3]} {
-                switch -- $numValues 2 {;# Ex. PHOsphate w/o H3PO4
-                    lassign $data attr64 attr40
-                    set attr62 $attr64
-                    set attr61 $attr64
-                    set attr54 $attr64
-                    set attr52 $attr64
-                    set attr51 $attr64
-                    set attr34 $attr64
-                    set attr32 $attr64
-                    set attr31 $attr64 
-                    set attr20 $attr40
-                    set attr10 $attr40
+        } ;# END 2 sites
+        3 {
+            switch -- $numStates {
+                4 {
+                    if {0 in $missingProtCnts && 1 in $missingProtCnts} {
+                        # state (0,0,0) and 1 proton states missing
+                        switch -- $numValues {
+                            1 { ;# Ex. LYS
+                                lassign $data attr76
+                                mset Matrix 7 6 $attr76
+                                mset Matrix 7 5 $attr76
+                                mset Matrix 7 3 $attr76
+                            }
+                            default {
+                                set errorCode -1
+                            }
+                        }
+                        if {$errorCode} break
 
-                    set attr60 [expr {$attr64 + $attr40}]
-                    set attr50 [expr {$attr54 + $attr40}]
-                    set attr30 [expr {$attr32 + $attr20}]
-                    # All other pairs are zero
+                        mset Matrix 6 5 [expr {[mindex $Matrix 7 5]\
+                                               - [mindex $Matrix 7 6]}]
+                        mset Matrix 6 3 [expr {[mindex $Matrix 7 3]\
+                                               - [mindex $Matrix 7 6]}]
+                        mset Matrix 5 3 [expr {[mindex $Matrix 7 3]\
+                                               - [mindex $Matrix 7 5]}]
 
-                    lset Matrix [index2flatindex 6 4] $attr64
-                    lset Matrix [index2flatindex 6 2] $attr62
-                    lset Matrix [index2flatindex 6 1] $attr61
-                    lset Matrix [index2flatindex 5 4] $attr54
-                    lset Matrix [index2flatindex 5 2] $attr52
-                    lset Matrix [index2flatindex 5 1] $attr51
-                    lset Matrix [index2flatindex 3 4] $attr34
-                    lset Matrix [index2flatindex 3 2] $attr32
-                    lset Matrix [index2flatindex 3 1] $attr31
-                    lset Matrix [index2flatindex 4 0] $attr40
-                    lset Matrix [index2flatindex 2 0] $attr20
-                    lset Matrix [index2flatindex 1 0] $attr10
-                    lset Matrix [index2flatindex 6 0] $attr60
-                    lset Matrix [index2flatindex 5 0] $attr50
-                    lset Matrix [index2flatindex 3 0] $attr30
-                } default {
+                        set alt53 [expr {[mindex $Matrix 6 3]\
+                                         - [mindex $Matrix 6 5]}]
+                        set alt63 [expr {[mindex $Matrix 5 3]\
+                                         - [mindex $Matrix 6 5]}]
+                        set alt65 [expr {[mindex $Matrix 6 3]\
+                                         - [mindex $Matrix 5 3]}]
+                        if {[mindex $Matrix 5 3] != $alt53
+                            || [mindex $Matrix 6 3] != $alt63
+                            || [mindex $Matrix 6 5] != $alt65} {
+                            set errorCode -2
+                        }
+                        # END 3 sites, 4 states, no state 0, 1, 2, 4
+                    } elseif {2 in $missingProtCnts && 3 in $missingProtCnts} {
+                        # state (1, 1, 1) and 2 proton states missing
+                        switch -- $numValues {
+                            1 { ;# Ex. phosphate monoesters
+                                lassign $data attr40
+                                mset Matrix 4 0 $attr40
+                                mset Matrix 2 0 $attr40
+                                mset Matrix 1 0 $attr40 
+                            } 
+                            default {
+                                set errorCode -1
+                            }
+                        }
+                        if {$errorCode} break
+                        
+                        mset Matrix 4 2 [expr {[mindex $Matrix 4 0]\
+                                               - [mindex $Matrix 2 0]}]
+                        mset Matrix 4 1 [expr {[mindex $Matrix 4 0]\
+                                               - [mindex $Matrix 1 0]}]
+                        mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
+                                               - [mindex $Matrix 1 0]}]
+                        # END 3 sites, 4 states, no state 3, 5, 6, 7
+                    } else {
+                        set errorCode -1
+                    }
+                } ;# END 3 sites, 4 states
+                7 {
+                    if {3 in $missingProtCnts} {
+                        # state (1, 1, 1) missing
+                        switch -- $numValues {
+                            2 { ;# Ex. PHOsphate w/o H3PO4
+                                lassign $data attr64 attr40
+                                mset Matrix 6 4 $attr64
+                                mset Matrix 6 2 $attr64
+                                mset Matrix 6 1 $attr64
+                                mset Matrix 5 4 $attr64
+                                mset Matrix 5 2 $attr64
+                                mset Matrix 5 1 $attr64
+                                mset Matrix 3 4 $attr64
+                                mset Matrix 3 2 $attr64
+                                mset Matrix 3 1 $attr64
+
+                                mset Matrix 4 0 $attr40
+                                mset Matrix 2 0 $attr40
+                                mset Matrix 1 0 $attr40
+                            }
+                            default {
+                                set errorCode -1
+                            }
+                        }
+                        if {$errorCode} break
+
+                        mset Matrix 6 0 [expr {[mindex $Matrix 6 4]\
+                                               + [mindex $Matrix 4 0]}]
+                        mset Matrix 5 0 [expr {[mindex $Matrix 5 4]\
+                                               + [mindex $Matrix 4 0]}]
+                        mset Matrix 3 0 [expr {[mindex $Matrix 3 1]\
+                                               + [mindex $Matrix 1 0]}]
+                        mset Matrix 6 5 [expr {[mindex $Matrix 6 4]\
+                                               - [mindex $Matrix 5 4]}]
+                        mset Matrix 6 3 [expr {[mindex $Matrix 6 2]\
+                                               - [mindex $Matrix 3 2]}]
+                        mset Matrix 5 3 [expr {[mindex $Matrix 5 1]\
+                                               - [mindex $Matrix 3 1]}]
+                        mset Matrix 4 2 [expr {[mindex $Matrix 4 0]\
+                                               - [mindex $Matrix 2 0]}]
+                        mset Matrix 4 1 [expr {[mindex $Matrix 4 0]\
+                                               - [mindex $Matrix 1 0]}]
+                        mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
+                                               - [mindex $Matrix 1 0]}]
+                        # END 3 sites, 7 states, no state 7
+                    } else {
+                        set errorCode -1
+                    }
+                } ;# END 3 sites, 7 states
+                default {
                     set errorCode -1
                 }
-            } else {
-                set errorCode -1
             }
-        # TODO: Add H3PO4 with 8 states total 
-        default {
+        } ;# END 3 sites
+        default {
             set errorCode -1
         }
-    } default {
-        set errorCode -1
     }
-    switch -- $errorCode -1 {
-        abort "Bad or unimplemented specification of $numValues values for\
-                $numSites site residue with $numStates states"
-    } -2 {
-        abort "Error in thermodynamic cycle"
+
+    # Catch errors here.
+    switch -- $errorCode {
+        -1 {
+            abort "Bad or unimplemented specification of $numValues values for\
+                    $numSites site residue with $numStates states"
+        }
+        -2 {
+            abort "Error in thermodynamic cycle"
+        }
     }
+
     return $Matrix
 }
 
index a191c76..4d32a93 100644 (file)
@@ -1,4 +1,3 @@
-# namdcph.core.tcl
 # 
 # Core utilities for setting up and running constant pH simulations in NAMD.
 #
@@ -14,23 +13,17 @@ namespace eval ::namdcph {
 
     variable TITLE "namdcph)"
     variable SystempH
-    variable configFilename ""
+    variable configFilenames [list]
     variable restartFilename ""
     variable restartFreq 0
     variable outFile ""
     variable numMinSteps 0
     variable excludeList [list]
-    variable AlchFrcCons 100.0
+    variable alchFrcCons 100.0
     variable MDBasename namdcph.md
     variable SWBasename namdcph.sw
-    # cphSystem data and defaults
-    variable stateInfo [dict create]
-    # cphTitrator data and defaults
-    variable moveInfo [dict create]
-    dict set moveInfo maxProposalAttempts 0
-    dict set moveInfo default numsteps 20 
-    dict set moveInfo default weight 1.0
-    dict set moveInfo ptransfer [list]
+    variable stateInfo [dict create] ;# cphSystem data and defaults
+    variable moveInfo [dict create] ;# cphTitrator data and defaults
 
     variable residueAliases [dict create]
     dict set residueAliases HIS {HSD HSE HSP}
@@ -46,18 +39,20 @@ namespace eval ::namdcph {
 # Run a set number of neMD/MC cycles.
 #
 proc ::namdcph::cphRun {numsteps {numcycles 1}} {
+    variable ::namdcph::SystempH
+
     # Initialize NAMD and build a constant pH enabled PSF.
     initialize
     finalize
     if {$::namdcph::numMinSteps > 0} {
-        set tmp [firstTimeStep]
+        set storedFirstTimeStep [firstTimeStep]
         minimize $::namdcph::numMinSteps
         if {[isset temperature]} {
           reinitvels [temperature]
         } else {
           reinitvels [$::thermostatTempCmd]
         }
-        firstTimeStep $tmp
+        firstTimeStep $storedFirstTimeStep
     }
     set cphlog [openCpHLog]
     set firstCycle 1
@@ -69,7 +64,8 @@ proc ::namdcph::cphRun {numsteps {numcycles 1}} {
         runMD {*}$runArgs
         # (2) Propose a move from the full move set.
         #
-        lassign [proposeMove] accept swNumsteps segresidnameList nattempts
+        lassign [cphTitrator propose $SystempH] accept swNumsteps\
+                segresidnameList nattempts
         # (3) At this point we have either selected a switch or rejected a
         # whole bunch of moves. If the former, perform the switch.
         #
@@ -81,7 +77,8 @@ proc ::namdcph::cphRun {numsteps {numcycles 1}} {
             set accept [runSwitch $swNumsteps $segresidnameList]
             set runArgs [list $numsteps]
             # Only accumulate statistics for attempted switches.
-            accumulateAcceptanceRate $accept $segresidnameList
+            set moveLabel [join $segresidnameList "/"]
+            cphTitrator accumulateStats $accept $moveLabel
         }
         cphSystem update $accept $segresidnameList
         # (4) Output cycle information and a checkpoint file if necessary.
@@ -122,7 +119,7 @@ proc ::namdcph::checkResidueDefinitions {resnames} {
     alchLambdaFreq 0
     foreach segresidname [cphSystem get segresidnames] {
         lassign [split $segresidname ":"] segid resid resname
-        if {[lsearch $resnames $resname] < 0} {
+        if {$resname ni $resnames} {
             continue
         } 
         set states [cphSystem get stateList $segresidname]
@@ -293,7 +290,8 @@ proc ::namdcph::pH {pHValue} {
 # Configuration file for constant pH residues.
 #
 proc ::namdcph::cphConfigFile {filename} {
-    variable ::namdcph::configFilename [string trim $filename]
+    variable ::namdcph::configFilenames
+    lappend configFilenames [string trim $filename]
     return
 }
 
@@ -402,7 +400,7 @@ proc ::namdcph::cphProposalWeight {args} {
 proc ::namdcph::cphProposeProtonTransfer {args} {
     variable ::namdcph::moveInfo
     foreach moveLabel $args {
-        dict lappend moveInfo ptransfer $moveLabel
+        dict set moveInfo $moveLabel protonTransfer 1
     }
     return
 }
@@ -452,7 +450,7 @@ proc ::namdcph::cphExcludeResidue {args} {
 #
 proc ::namdcph::cphForceConstant {forceConstant} {
     checkIsNotNegative cphForceConstant $forceConstant
-    variable ::namdcph::AlchFrcCons $forceConstant
+    variable ::namdcph::alchFrcCons $forceConstant
     return
 }
 
@@ -580,8 +578,9 @@ proc ::namdcph::runSwitch {numsteps segresidnameList} {
 # None
 #
 proc ::namdcph::alchemify {segresidnameList} {  
+    variable ::namdcph::alchFrcCons
     set T [$::thermostatTempCmd]
-    set FrcCons [getAlchFrcCons]
+
     alch on
     # (1) Read in the nonalchemical PSF and apply patches.
     #
@@ -590,7 +589,7 @@ proc ::namdcph::alchemify {segresidnameList} {
     # (2) Apply patches and build coordinates and velocities.
     #
     foreach segresidname $segresidnameList {
-        cphSystem alchemifypsf $segresidname $FrcCons $T
+        cphSystem alchemifypsf $segresidname $alchFrcCons $T
     }
     regenerate angles dihedrals
     # (3) Write a new set of inputs and reinitialize.
@@ -602,7 +601,8 @@ proc ::namdcph::alchemify {segresidnameList} {
     #
     set ExtraBondsFile [open [swFilename extrabonds] "w"]
     foreach segresidname $segresidnameList {
-        puts $ExtraBondsFile [cphSystem get alchBonds $segresidname $FrcCons]
+        puts $ExtraBondsFile\
+                [cphSystem get alchBonds $segresidname $alchFrcCons]
     }
     close $ExtraBondsFile
     reloadAndReinit [getSWBasename] true
@@ -640,7 +640,7 @@ proc ::namdcph::dealchemify {segresidnameList} {
 # =============================================================================
 # ::namdcph::readConfig
 #
-# Read a constant pH config file and return template definitions.
+# Read one or more constant pH config files and return template definitions.
 #
 # Arguments:
 # ----------
@@ -652,55 +652,109 @@ proc ::namdcph::dealchemify {segresidnameList} {
 #   A dict conversion of the config file contents from JSON
 #
 proc ::namdcph::readConfig {} {
-    #TODO: permit reading of multiple config files
-    variable ::namdcph::configFilename
-    set ConfigFile [open $configFilename "r"]
-    set templateDefs [json::json2dict [read $ConfigFile]]
-    close $ConfigFile
+    variable ::namdcph::configFilenames
+
+    set templateDefs [dict create]
+    foreach configFilename $configFilenames {
+        set ConfigFile [open $configFilename "r"]
+        set tmpDict [json::json2dict [read $ConfigFile]]
+        # Warn the user about re-definitions, but proceed anyway. Note that
+        # dict merge overwrites values right to left.
+        set oldKeys [dict keys $templateDefs]
+        foreach newKey [dict keys $tmpDict] {
+            if {$newKey in $oldKeys} {
+                cphWarn "Reading and using new definition for residue $newKey"
+            }
+        } 
+        set templateDefs [dict merge $templateDefs $tmpDict] 
+        close $ConfigFile
+    }
     return $templateDefs
 }
 
 # ::namdcph::readRestart
 #
-# Read in a constant pH state from file.
+# Read in a constant pH state from file and _merge_ the data into the namespace
+# variables. That is, give precedence to keyword specifications.
 #
-# Arguments:
-# ----------
-# restartFilename : string
-#   Name of (JSON format) restart file to read
-#
-# Returns:
-# --------
-# stateList : list
-#   State codes for each titratable residue
-# pKaiList : list
-#   Latest estimate of the inherent pKa
+# NB! This also requires that cphSystem has been built so that the system and
+#     state info can be compared.
 #
-proc ::namdcph::readRestart {restartFilename} {
+proc ::namdcph::readRestart {} {
+    variable ::namdcph::restartFilename
+    variable ::namdcph::SystempH
+    variable ::namdcph::stateInfo
+    variable ::namdcph::moveInfo
+
     set RestartFile [open $restartFilename "r"]
     set Restart [json::json2dict [read $RestartFile]]
     close $RestartFile
+
+    if {![info exists SystempH] && [dict exists $Restart pH]} {
+        pH [dict get $Restart pH]
+    }
+
+    # Read state parameters, if present. 
+    #
     if {[dict exists $Restart states]} {
         set stateList [dict get $Restart states]
-    } else {
-        set stateList {}
     }
     if {[dict exists $Restart pKais]} {
         set pKaiList [dict get $Restart pKais]
-    } else {
-        set pKaiList {}
     }
-    if {[llength $stateList] > 0 && [llength $pKaiList] > 0} {
+    if {[info exists stateList] && [info exists pKaiList]} {
+        set segresidnameList [cphSystem get segresidnames]
+
         if {[llength $stateList] != [llength $pKaiList]} {
             cphAbort "mismatch in states/pKais in $restartFilename"
         }
+        if {[llength $stateList] != [llength $segresidnameList]} {
+            cphAbort "Too few/many state/pKai definitions in $restartFilename"
+        }
+
+        foreach segresidname $segresidnameList\
+                state $stateList pKai $pKaiList {
+            if {![dict exists $stateInfo $segresidname]} {
+                dict set stateInfo $segresidname state $state
+                dict set stateInfo $segresidname pKai $pKai
+            } else {
+                if {![dict exists $stateInfo $segresidname state]} {
+                    dict set stateInfo $segresidname state $state
+                }
+                if {![dict exists $stateInfo $segresidname pKai]} {
+                    dict set stateInfo $segresidname pKai $pKai
+                }
+            }
+        }
     }
+
+    # Read MC move parameters, if present. Only reset those values which have
+    # not been set in via keywords.
+    #
     if {[dict exists $Restart MCmoves]} {
-        set MCmoves [dict get $Restart MCmoves]    
-    } else {
-        set MCmoves [dict create]
+        set rstrtMoveInfo [dict get $Restart MCmoves]    
+    }
+    if {[info exists rstrtMoveInfo]} {
+        # This is a bit kludgy - the global data is not a nested dict, so it
+        # will crash in the dict for loop. Treat these specially and then unset
+        # them...
+        if {[dict exists $rstrtMoveInfo maxProposalAttempts]} {
+            if {![dict exists $moveInfo maxProposalAttempts]} {
+                dict set moveInfo maxProposalAttempts\
+                    [dict get $rstrtMoveInfo maxProposalAttempts]
+            }
+            dict unset rstrtMoveInfo maxProposalAttempts
+        }
+        dict for {moveLabel data} $rstrtMoveInfo {
+            dict for {key value} $data {
+                 if {![dict exists $moveInfo $moveLabel $key]} {
+                     dict set moveInfo $moveLabel $key $value
+                 }
+            }
+        }
     }
-    return [list $stateList $pKaiList $MCmoves]
+
+    return
 }
 
 # ::namdcph::writeRestart
@@ -722,6 +776,8 @@ proc ::namdcph::readRestart {restartFilename} {
 #
 proc ::namdcph::writeRestart {args} { 
     variable ::namdcph::restartFreq
+    variable ::namdcph::SystempH
+
     if {[string match [lindex $args 0] force]} {
         set restartFilename [lindex $args 1]
         set cycle [expr {int([lindex $args 2])}]
@@ -737,11 +793,14 @@ proc ::namdcph::writeRestart {args} {
     # TODO: This is god-awful. Is there not a useable json encoder that can do
     # this for us with some guarantee of accuracy?
     set cycleStr "\"cycle\":$cycle"
+    set pHStr "\"pH\":$SystempH"
+
     set stateStr "\"states\":\["
     foreach state [cphSystem get state] {
         set stateStr "$stateStr\"$state\","
     }
     set stateStr "[string trimright $stateStr ","]\]"
+
     set pKaiStr "\"pKais\":\["
     foreach pKaiList [cphSystem get pKai] {
         set pKaiStr "$pKaiStr\["
@@ -751,15 +810,26 @@ proc ::namdcph::writeRestart {args} {
         set pKaiStr "[string trimright $pKaiStr ","]\],"
     }
     set pKaiStr "[string trimright $pKaiStr ","]\]"
-    set MCStr "\"MCmoves\":\{"
-    dict for {moveLabel data} [getMoveSet] {
-        set numsteps [dict get $data numsteps]
-        set weight [dict get $data weight]
-        set MCStr "$MCStr\"$moveLabel\":\{\"numsteps\":$numsteps,\"weight\":$weight\},"
+
+    set maxAttempts [cphTitrator get maxAttempts]
+    set MCStr "\"MCmoves\":\{\"maxProposalAttempts\":$maxAttempts,"
+
+    set defaultNumsteps [cphTitrator get numsteps default]
+    set defaultWeight [cphTitrator get weight default]
+    set MCStr "$MCStr\"default\":\{\"numsteps\":$defaultNumsteps,\"weight\":$defaultWeight\},"
+
+    foreach moveLabel [cphTitrator get moveLabels] {
+        set thisMoveInfo [cphTitrator get nodefaults $moveLabel]
+        if {![dict size $thisMoveInfo]} continue
+        set MCStr "$MCStr\"$moveLabel\":\{"
+        dict for {key value} $thisMoveInfo {
+            set MCStr "$MCStr\"$key\":$value,"
+        }
+        set MCStr "[string trimright $MCStr ","]\},"
     }
     set MCStr "[string trimright $MCStr ","]\}"
 
-    puts $RestartFile "\{$cycleStr,$stateStr,$pKaiStr,$MCStr\}"
+    puts $RestartFile "\{$cycleStr,$pHStr,$stateStr,$pKaiStr,$MCStr\}"
     close $RestartFile
     # Always checkpoint the PSF and PDB when a restart is written.
     file copy -force [mdFilename psf] "[outputName].psf"
@@ -773,6 +843,8 @@ proc ::namdcph::writeRestart {args} {
 #
 proc ::namdcph::openCpHLog {} {
     variable ::namdcph::outFile
+    variable ::namdcph::SystempH
+
     if {[string length $outFile] > 0} {
         set logFilename $outFile
     } else {
@@ -780,7 +852,7 @@ proc ::namdcph::openCpHLog {} {
     }
     namdFileBackup $logFilename
     set cphlog [open $logFilename "w"]
-    puts $cphlog "#pH [getpH]"
+    puts $cphlog "#pH $SystempH"
     puts $cphlog "#[join [cphSystem get segresidnames] " "]"
     return $cphlog 
 }
@@ -830,8 +902,12 @@ proc ::namdcph::reloadAndReinit {basename keepExtraBonds} {
     return
 }
 
-proc ::namdcph::cphPrint {text} {
-    print "$::namdcph::TITLE $text"
+proc ::namdcph::cphPrint {msg} {
+    print "$::namdcph::TITLE $msg"
+}
+
+proc ::namdcph::cphWarn {msg} {
+    print "$::namdcph::TITLE WARNING! $msg"
 }
 
 proc ::namdcph::cphAbort {msg} {
@@ -839,33 +915,13 @@ proc ::namdcph::cphAbort {msg} {
 }
 
 proc ::namdcph::cphWarnExpt {} {
-    cphPrint "WARNING! THIS FEATURE IS EXPERIMENTAL!"
+    cphWarn "THIS FEATURE IS EXPERIMENTAL!"
     cphPrint "RESULTS ARE NOT GUARANTEEED - USE AT YOUR OWN RISK"
 }
 
 # =============================================================================
 # Setup Routines
 # =============================================================================
-# ::namdcph::checkSettings
-#
-# Check the settings from the configuration file and verify that they are
-# compatible with constant pH.
-#
-proc ::namdcph::checkSettings {} {
-    variable ::namdcph::configFilename
-    # Check constant pH specific settings. 
-    if {[string length configFilename] <= 0} {
-        cphAbort "A constant pH configuration file is required."
-    }
-    if {![info exists ::namdcph::SystempH]} {
-        cphAbort "A pH value is required."
-    }
-
-    getThermostat
-    getBarostat
-    return
-}
-
 # ::namdcph::initialize
 #
 # Initialize the system for constant pH. This requires two main things to
@@ -877,61 +933,60 @@ proc ::namdcph::checkSettings {} {
 #
 proc ::namdcph::initialize {} {
     variable ::namdcph::restartFilename
+    variable ::namdcph::configFilenames
     variable ::namdcph::stateInfo
     variable ::namdcph::moveInfo
     variable ::namdcph::residueAliases
     variable ::namdcph::excludeList
-    checkSettings
-    callback energyCallback    
-    # 1) Set up alchemical keywords and run an energy evaluation.
+    variable ::namdcph::SystempH
+
+    getThermostat
+    getBarostat
+    callback energyCallback
+    # 1) Set up alchemical keywords and run startup. 
     #
     initializeAlch
+
     # 2) Rebuild the PSF with dummy protons and modify protonation states as 
     # needed. Build the residue definitions, assign states to each residue, and 
     # rebuild the topology to reflect those states.
     #
     cphPrint "initializing constant pH PSF..."
+    # TODO: integrate aliases into the configfile standard...
+    if {![llength configFilenames]} {
+        cphAbort "At least one constant pH configuration file is required."
+    }
     cphSystem build [readConfig] $residueAliases $excludeList
 
-    if {[string length $restartFilename] > 0} {
-        lassign [readRestart $restartFilename] states pKais MCmoves
+    if {[string length $restartFilename]} {
+        readRestart
         lassign [list 0.0 false] temp buildH
     } else {
-        lassign [list {} {} [dict create]] states pKais MCmoves
         lassign [list [$::thermostatTempCmd] true] temp buildH
     }
-    # Use the system topology to make assignments from the restart
-    if {[llength $states] && [llength $pKais]} { 
-        foreach segresidname [cphSystem get segresidnames]\
-                state $states pKai $pKais {
-            if {![dict exists $stateInfo $segresidname]} {
-                dict set stateInfo $segresidname state $state
-                dict set stateinfo $segresidname pKai $pKai
-            } else {
-                if {![dict exists $stateInfo $segresidname state]} {
-                    dict set stateInfo $segresidname state $state
-                }
-                if {![dict exists $stateInfo $segresidname pKai]} {
-                    dict set stateInfo $segresidname pKai $pKai
-                }
-            }
+
+    if {![info exists SystempH]} {
+        cphAbort "A pH value is required."
+    }
+
+    # All residue state info from all sources is now in stateInfo.
+    # Check that all of the keys are valid.
+    foreach segresidname [dict keys $stateInfo] {
+        if {[cphSystem validate $segresidname]} {
+            cphAbort
         }
     }
-    cphSystem initialize [getpH] $temp $buildH $stateInfo
+    cphSystem initialize $SystempH $temp $buildH $stateInfo
+
     # 3) Build the MC move set (the "titrator").
-    set moveInfo [dict merge $MCmoves $moveInfo]
-    set unusedInfo [buildTitrator [getpH] $moveInfo]
+    #
+    # All neMD/MC move info from all sources is now in moveInfo.
+    cphTitrator build $moveInfo
+
     # 4) Report to stdout.
     printSettingsSummary
     printSystemSummary
     printTitratorSummary
-    if {[dict size [dict remove $unusedInfo default]]} {
-        cphPrint "WARNING! Unused/unrecognized MC move info!"
-        cphPrint [format "%-14s : %-s" "Move Label" "Info"]
-        dict for {moveLabel data} [dict remove $unusedInfo default] {
-            cphPrint [format "%-14s : %-s" $moveLabel $data]
-        }
-    }
     return
 }
 
@@ -1005,18 +1060,24 @@ proc ::namdcph::initializeAlch {} {
 # Print a summary of the constant pH settings.
 #
 proc ::namdcph::printSettingsSummary {} {
-    variable ::namdcph::configFilename
+    variable ::namdcph::configFilenames
     variable ::namdcph::restartFilename
+    variable ::namdcph::SystempH
+    variable ::namdcph::alchFrcCons
+
     set StarBar "***************************************"
     cphPrint $StarBar
     cphPrint "CONSTANT pH MD ACTIVE"
     if {[string length $restartFilename] > 0} {
         cphPrint "RESTART FILENAME $restartFilename"
     }
-    cphPrint "SYSTEM pH [getpH]"
-    cphPrint "CONSTANT pH CONFIGURATION FILE $configFilename"
+    cphPrint "SYSTEM pH $SystempH"
+    cphPrint "CONSTANT pH CONFIGURATION FILE(S)"
+    foreach configFilename $configFilenames {
+        cphPrint "$configFilename"
+    }
     cphPrint "NONEQUILIBRIUM SWITCH PARAMETERS:"
-    cphPrint "ALCHEMICAL FORCE CONSTANT [getAlchFrcCons] kcal/mol-A^2"
+    cphPrint "ALCHEMICAL FORCE CONSTANT $alchFrcCons kcal/mol-A^2"
     cphPrint "neMD/MC CRITERION TEMPERATURE [$::thermostatTempCmd]"
     cphPrint "TEMPORARY FILE INFORMATION:"
     cphPrint "cpH TOPOLOGY FILE BASENAME [getMDBasename]"
@@ -1053,14 +1114,15 @@ proc ::namdcph::printSystemSummary {} {
 # Print a summary of the MC moves.
 #
 proc ::namdcph::printTitratorSummary {} {
-    set StarBar "*******************************************"
+    set StarBar "*************************************************"
     cphPrint $StarBar
     cphPrint "CONSTANT pH neMD/MC MOVES:"
-    cphPrint [format "%-19s : %8s %12s" "move label" numsteps weight]
-    dict for {moveLabel data} [getMoveSet] {
-        set numsteps [dict get $data numsteps]
-        set weight [dict get $data weight]
-        cphPrint [format "%-19s : % 8d % 12.2f" $moveLabel $numsteps $weight]
+    cphPrint "MAX. ATTEMPTS PER CYCLE: [cphTitrator get maxAttempts]"
+    cphPrint [format "%-25s : %8s %12s" "move label" numsteps weight]
+    foreach moveLabel [cphTitrator get moveLabels] {
+        set numsteps [cphTitrator get numsteps $moveLabel]
+        set weight [cphTitrator get weight $moveLabel]
+        cphPrint [format "%-25s : % 8d % 12.2f" $moveLabel $numsteps $weight]
     }
     cphPrint $StarBar
     return
@@ -1082,19 +1144,19 @@ proc ::namdcph::printProposalSummary {segresidnameList} {
 # Print a report of the titration MC statistics.
 #
 proc ::namdcph::printTitratorReport {} {
-    set StarBar "*******************************************"
+    set StarBar "*************************************************"
     cphPrint $StarBar
     cphPrint "CONSTANT pH MD STATISTICS:"
-    cphPrint [format "%-19s : %-8s %-12s" "move label" attempts "accept. rate"]
-    dict for {moveLabel data} [getMoveSet] {
-        set attempts [dict get $data attempts]
-        set successes [dict get $data successes]
+    cphPrint [format "%-25s : %-8s %-12s" "move label" attempts "accept. rate"]
+    foreach moveLabel [cphTitrator get moveLabels] {
+        set attempts [cphTitrator get attempts $moveLabel]
+        set successes [cphTitrator get successes $moveLabel]
         if {$successes && $attempts} {
             set rate [expr {100.*$successes/$attempts}]
         } else {
             set rate 0.0
         }
-        cphPrint [format "%-19s : %8d %12.2f" $moveLabel $attempts $rate] 
+        cphPrint [format "%-25s : %8d %12.2f" $moveLabel $attempts $rate] 
     }
     cphPrint $StarBar
     return
@@ -1105,14 +1167,6 @@ proc ::namdcph::printTitratorReport {} {
 #
 # These are largely unnecessary, but cut down on "variable" declarations.
 # =============================================================================
-proc ::namdcph::getpH {} {
-    return $::namdcph::SystempH
-}
-
-proc ::namdcph::getAlchFrcCons {} {
-    return $::namdcph::AlchFrcCons
-}
-
 proc ::namdcph::getMDBasename {} {
     return $::namdcph::MDBasename
 }
index 24f3eda..166bc12 100644 (file)
@@ -32,47 +32,116 @@ proc tcl::mathfunc::lmean {list} {
 
 # lincr
 #
-# Increment an element of a list by the given value.
-# Indexing uses the standard syntax for nested lists, similar to lindex.
-#
+# Increment a list (in place) by the given value. The modified list is also
+# returned.
+#
+# - This supports the standard lindex syntax for nested lists (ala lindex).
+# - NumPy style "broadcasting" is also supported, whereby a list plus a value
+#   results in each element of the list being incremented and a list plus a
+#   list results in element by element addition.
+#   NB: This only works left to right (ala +=) and modifies the left list.
+#
+# Example:
+#
+# % set foo {0 2 1 7 3}
+# % lincr foo 2
+# 2 4 3 9 5
+# % lincr foo {0 1 2 3 4}
+# 2 5 5 12 9
+#
+# % set bar {0 1 {2 3 4} 5}
+# % lincr bar 2 1
+# 0 1 {3 4 5} 5
+# % lincr bar 2 1 1
+# 0 1 {3 5 5} 5
+# 
 proc lincr {list args} {
     upvar 1 $list List
     set value [lindex $args end]
     set index [lrange $args 0 end-1]
-    lset List {*}$index [expr {[lindex $List {*}$index] + $value}]
-    return
-}
+    
+    set target [lindex $List {*}$index]
+    set vlen [llength $value]
+    set tlen [llength $target]
 
-# llincr
-#
-# Increment each element of list1 by the elements in list2.
-#
-proc llincr {list1 args} {
-    upvar 1 $list1 List1
-    set list2 [lindex $args end]
-    set index [lrange $args 0 end-1]
-    set i 0
-    foreach item1 [lindex $List1 {*}$index] item2 $list2 {
-        lset List1 {*}$index $i [expr {$item1 + $item2}]
-        incr i
+    if {$vlen == 1} {
+        if {$tlen == 1} {
+            # increment a single element by a single value
+            lset List {*}$index [expr {$target + $value}]
+        } elseif {$tlen > 1} {
+            # increment a each element in a list by a single value
+            set i 0 
+            foreach item $target {
+                if {[llength $item] != 1} {
+                    error "lincr encountered a nested list during list + value\
+                            operation"
+                }
+                lset List {*}$index $i [expr {$item + $value}]
+                incr i
+            }
+        }
+    } elseif {$vlen > 1} {
+        if {$tlen != $vlen} {
+            error "lincr cannot perform list + list operation with different\
+                    lengths ($tlen != $vlen)"
+        }
+        set i 0
+        foreach item1 $target item2 $value {
+            lset List {*}$index $i [expr {$item1 + $item2}]
+            incr i
+        }
     }
-    return
+    return $List
 }
 
 # lmultiply
 #
-# Multiple each element of a list by the given value.
+# Multiple a list (in place) by the given value. The modified list is also
+# returned.
+#
+# - This supports the standard lindex syntax for nested lists (ala lindex).
+# - NumPy style "broadcasting" is also supported, whereby a list times a value
+#   results in each element of the list being multiplied and a list times a
+#   list results in element by element multiplication.
+#   NB: This only works left to right (ala *=) and modifies the left list.
 #
 proc lmultiply {list args} {
     upvar 1 $list List
     set value [lindex $args end]
     set index [lrange $args 0 end-1]
-    set i 0
-    foreach item [lindex $List {*}$index] {
-        lset List {*}$index $i [expr {$item*$value}]
-        incr i
+
+    set target [lindex $List {*}$index]
+    set vlen [llength $value]
+    set tlen [llength $target]
+
+    if {$vlen == 1} {
+        if {$tlen == 1} {
+            # increment a single element by a single value
+            lset List {*}$index [expr {$target*$value}]
+        } elseif {$tlen > 1} {
+            # increment a each element in a list by a single value
+            set i 0
+            foreach item $target {
+                if {[llength $item] != 1} {
+                    error "lmultiply encountered a nested list during list *\
+                            value operation"
+                }
+                lset List {*}$index $i [expr {$item*$value}]
+                incr i
+            }
+        }
+    } elseif {$vlen > 1} {
+        if {$tlen != $vlen} {
+            error "lmultiply cannot perform list * list operation with\
+                    different lengths ($tlen != $vlen)"
+        }
+        set i 0
+        foreach item1 $target item2 $value {
+            lset List {*}$index $i [expr {$item1*$item2}]
+            incr i
+        }
     }
-    return
+    return $List
 }
 
 # =============================================================================
@@ -96,46 +165,47 @@ proc tcl::mathfunc::normal {{mu 0.0} {sigma 1.0}} {
 
 # choice
 #
-# Choose a random element from a list with equal probability.
+# Choose a random element from a list. The index of the element is also
+# returned.
+#
+# The default is to uniformly weight the elements, but (unnormalized)
+# non-uniform weights may also be provided.
 #
 # Arguments:
 # ----------
 # list : list
 #    list to choose an element from
-#
-proc choice {list} {
-    return [lindex $list [expr {int(rand()*[llength $list])}]]
-}
-
-# weightedChoice
-#
-# Given a set of probability weights, return the index of a given outcome with 
-# the correct probability.
-#
-# Arguments:
-# ----------
-# weights : list of floats
+# weights : list of floats (optional)
 #   The probability weight of each choice. These need not be normalized.
-# weightSum : float (optional, calculated if not specified)
-#   The sum of all entries in weights. This will not function correctly if this
-#   is incorrectly specified!
-#
-proc weightedChoice {weights {weightSum {}}} {
-    if {[llength $weightSum]} {
-        set WeightSum $weightSum
-    } else {
-        set WeightSum [expr {lsum($weights)}]
+#
+# Returns:
+# --------
+# element
+#    the selected list element
+# index : int
+#    the index of the selected element
+#
+proc choice {list {weights {}}} {
+    # Uniform weights is really easy.
+    if {![llength $weights]} {
+        set j [expr {int(rand()*[llength $list])}]
+        return [list [lindex $list $j] $j]
+    }
+    # Non-uniform weighting is more involved.
+    if {[llength $list] != [llength $weights]} {
+        error "Mismatch in list/weights in choice!"
     }
+    set WeightSum [expr {lsum($weights)}]
     set Rand [expr {$WeightSum*rand()}]
     set j 0
     foreach Weight $weights {
         set Rand [expr {$Rand - $Weight}]
         if {$Rand <= 0.0} {
-            return $j
+            return [list [lindex $list $j] $j]
         }
         incr j
     }
     # This should never be reached.
-    error "Something bad happened in weightedChoice!"
+    error "Something bad happened in choice!"
 }