Bug fixes in Monte Carlo functionality 27/4327/4
authorradakb <brian.radak@gmail.com>
Thu, 5 Jul 2018 13:37:21 +0000 (09:37 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Mon, 9 Jul 2018 21:25:42 +0000 (16:25 -0500)
This update is an unfortunate hodge-podge of changes. A significant
error in the algorithm was discovered while introducing new features and
reorganizing the output code.

MAJOR BUG FIXES
---------------
1) cphMaxProposalAttempts is an extremely problematic parameter for
  small systems. This keyword is now deprecated and internally assigned
  to be one. Simulations of small systems with non-unit values almost
  certainly yielded incorrect statistics, but the effect is much smaller
  for large systems with uniform weights (the default).

2) Equivalent states were incorrectly accounted for during the inherent
  pKa proposal step (a mistake from the original paper). The pKas are
  now correctly shifted by the log multiplicity to obtain the
  _microscopic_ pKa when moving in the direction of a multiplicity.

  An unrelated error in the code seems to have bizarrely compensated for
  this problem. Namely, the unnormed pairwise weights were taken with
  respect to the current state such that the unnormed weight of the
  current state is exactly one. However, these weights were computed in
  log-space such that the assignment of one was incorrect (it should be
  zero). This meant that all other weights were reduced by a factor of
  1/e in a manner that did not maintain detailed balance in certain
  instances. Since 1/e ~= 1/2 and 1/3, this seems to have semi-corrected
  residues with those multiplicities.

As severe as the second error sounds, tests indicate that the most
significant problems were actually from the first error. Systematic
errors from the latter are actually probably at or below statistical
precision (+/- 1.0 pKa units).

All of these errors were unhelpfully masked by the published WHAM
procedure.

MAJOR CHANGES
-------------
1) The inherent pKa MC step now precludes transitions between equivalent
  states. This should be more efficient in nearly all instances and was
  a fairly heavily requested feature. The original aim was to be able to
  toggle this behavior on a per residue basis, but this was put on hold
  when bugs starting cropping up.

MINOR CHANGES
-------------
1) Re-organized archiving the titrator to cphrst so that this is more
  flexible.

2) Added preliminary support for proton transfer and co-titration moves,
  including separate enumeration in the move summary and persistence in
  the restart file. However, in light of the MC errors this
  functionality is now suspect and thus deactivated until it can be
  retested.

3) Various bug-fixes in the initialization procedure. These were
  probably not encountered by most users and would have lead to crashes.

4) Disable obsolete building of hydrogen coordinates when not restarting
  from a known state. This was introduced for older versions of the
  build procedure and is catastrophic for certain instances such as
  certain primary amines, since "chirality" was occasionally built
  incorrectly.

Change-Id: If2692eef0110d3ad835090481371a65798fecdf6

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

index 2616f02..32c9856 100644 (file)
@@ -19,7 +19,7 @@ namespace eval ::namdcphwrapper {
             cphRestartFile cphRestartFreq\
             cphForceConstant cphMDBasename cphSwitchBasename\
             cphMaxProposalAttempts cphNumMinSteps cphProposalWeight\
-            cphExcludeResidue cphProposeProtonTransfer\
+            cphExcludeResidue\
             cphNumstepsPerSwitch\
             checkResidueDefinitions testResidueSwitch cphAnalyzeForce
 
index 07b9784..29a0d54 100644 (file)
@@ -138,68 +138,111 @@ proc ::cphTitrator::dictPopOrDefault {myDict key {defaultValue {}}} {
 # the information necessary to perform a switch for the given proposal.
 #
 proc ::cphTitrator::proposeMove {pH} {
-    variable ::cphTitrator::maxAttempts
     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]
-    }
+    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]
+    return [list $accept $numsteps $segresidnameList]
 }
 
+# FOR TESTING ONLY - INDEPENDENCE SAMPLING
+#
+#proc ::cphTitrator::proposeResidueMove {segresidname pH} {
+#    lassign [cphSystem compute inherentWeights $pH $segresidname] states weights
+#    lassign [choice $states $weights] trialState
+#    set currState [cphSystem get state $segresidname]
+#    if {$trialState == $currState} {
+#        return 0
+#    } else {
+#        cphSystem set trialState $segresidname $trialState
+#        return 1
+#    }
+#}
+
+# FOR TESTING ONLY - METROPOLIS MONTE CARLO (NON-ERGODIC IN SOME CASES!)
+#
+#proc ::cphTitrator::proposeResidueMove {segresidname pH} {
+#    cphSystem propose titration $segresidname
+#    set du [cphSystem compute inherent $pH $segresidname]
+#    return [metropolisAcceptance $du]
+#}
+
 # ::cphTitrator::proposeResidueMove
 #
 # Propose a move involving a single residue via Metropolized independence
-# sampling. Return True if such a move was proposed, accepted, and stored.
-#
-# The core idea here is that we want to propose the most probable new state
-# based on the estimated inherent pKa values. For two state systems, this is
-# trivial, it's just the other state. For three or more states, this will
-# naturally prefer tautomerization at extreme pKa values.
+# sampling. Return True if such a move was proposed, accepted, and stored,
+# otherwise return False.
 #
 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]
-   }
-   return $accept
+    lassign [cphSystem compute inherentWeights $pH $segresidname] states weights
+    # Construct the Metropolized proposal weights for i-->j.
+    set equivStates [cphSystem get equivStateList $segresidname]
+    set proposalWeights [list]
+    set stateiProposalNorm 1.0
+    foreach state $states weight $weights {
+        if {$state in $equivStates} {
+            lappend proposalWeights 0.0
+            set stateiProposalNorm [expr {$stateiProposalNorm - $weight}]
+        } else {
+            lappend proposalWeights $weight
+        }
+    }
+    lassign [choice $states $proposalWeights] trialState
+    # Construct the Metropolized proposal weights for j-->i.
+    set equivStates [cphSystem get equivStateList $segresidname $trialState]
+    set proposalWeights [list]
+    set statejProposalNorm 1.0
+    foreach state $states weight $weights {
+        if {$state in $equivStates} {
+            lappend proposalWeights 0.0
+            set statejProposalNorm [expr {$statejProposalNorm - $weight}]
+        } else {
+            lappend proposalWeights $weight
+        }
+    }
+    # NB: Invert log argument bc of negative sign in Metropolis routine.
+    set du [expr {log($statejProposalNorm / $stateiProposalNorm)}]
+    set accept [metropolisAcceptance $du]
+    if {$accept} {
+        cphSystem set trialState $segresidname $trialState
+    }
+    return $accept
 }
 
-# ::cphTitrator::proposeProtonTransferMove
+# ::cphTitrator::proposeTwoResidueMove
+#
+# Propose a move involving two residues. This may take two forms:
+#
+# 1) a proton transfer, in which a proton is moved from one resiue to another
 #
-# 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.
+# 2) a co-titration, in which two residues change concurrently
 #
-# Note that this automatically fails if no proton can be transferred.
+# Return True if either kind of move is proposed, accepted, and stored.
 #
-proc ::cphTitrator::proposeProtonTransferMove {moveLabel pH} {
-    lassign [split $moveLabel "/"] segresidname1 segresidname2
+proc ::cphTitrator::proposeTwoResidueMove {moveLabel pH} {
+    lassign [split $moveLabel "/"] srn1 srn2 ;# srn = segresidname
 
-    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}]]
+    set doPT [cphTitrator get protonTransfer $moveLabel]
+    set doCT [cphTitrator get cotitration $moveLabel]
+    if {$doPT} {
+        set doPT [expr {![cphSystem propose protonTransfer $srn1 $srn2]}]
     }
-    return $accept
+    # We don't need to bother with co-titration if proton transfer is already
+    # proposed, since that means it isn't even possible. The order of these
+    # could equally be switched, but it is generally faster to check that
+    # proton transfer is possible, so this is done first.
+    #
+    if {!$doPT && $doCT} {
+        set doCT [expr {![cphSystem propose cotitration $srn1 $srn2]}]
+    }
+    if {$doPT || $doCT} {
+        set du1 [cphSystem compute inherent $pH $srn1]
+        set du2 [cphSystem compute inherent $pH $srn2]
+        return [metropolisAcceptance [expr {$du1 + $du2}]]
+    }
+    return 0
 }
 
 # =============================================================================
@@ -217,17 +260,14 @@ proc ::cphTitrator::buildTitrator {moveInfo} {
     # 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 maxAttempts 1
     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
+    foreach attr {weight protonTransfer cotitration}\
+            hardcodedDefault {1.0 0 0} {
+        if {![dict exists $defaultMoveParams $attr]} {
+            cphTitrator set $attr default $hardcodedDefault
+        }
     }
 
     # 2) Validate all remaining keys - these better be move labels!
@@ -244,27 +284,41 @@ proc ::cphTitrator::buildTitrator {moveInfo} {
         cphTitrator set proposalCmd $moveLabel\
                 "proposeResidueMove $segresidname"
         if {![dict exists $moveInfo $moveLabel]} continue
+
         dict for {attr value} [dict get $moveInfo $moveLabel] {
             cphTitrator set $attr $moveLabel $value
         }
+        dict unset moveInfo $moveLabel
     }
-    # Build proton transfer moves, if present.
+    # Build proton transfer and co-titration moves, if present.
     dict for {moveLabel data} $moveInfo {
-        if {![dict exists $data protonTransfer]
-            || ![dict get $data protonTransfer]} {
-            continue
+        set isTwoResidueMove 0
+        foreach type {protonTransfer cotitration} {
+            if {[dict exists $data $type] && [dict get $data $type]} {
+                set isTwoResidueMove 1
+                break
+            }
         }
+        if {!$isTwoResidueMove} continue
         if {[llength [split $moveLabel "/"]] != 2} {
-            abort "Proton transfer moves must have exactly 2 segresidnames!"
+            abort "Proton transfer and co-titration moves must have exactly 2\
+                   segresidnames!"
         }
         cphTitrator set proposalCmd $moveLabel\
-                "proposeProtonTransferMove $moveLabel"
-        if {![dict exists $moveInfo $moveLabel]} continue
-        dict for {attr value} [dict get $moveInfo $moveLabel] {
+                "proposeTwoResidueMove $moveLabel"
+        dict for {attr value} $data {
             cphTitrator set $attr $moveLabel $value
         }
+        dict unset moveInfo $moveLabel
+    }
+    # Abort if extra invalid move information was provided.
+    if {[dict size $moveInfo]} {
+        print "ERROR! Invalid move specifications:"
+        dict for {moveLabel data} $moveInfo {
+            print "$moveLabel $data"
+        }
+        abort
     }
-
     # Initialize statistics for any moves that are new.
     dict for {moveLabel data} $moveDict {
         if {![dict exists $data attempts]} {
@@ -298,6 +352,7 @@ proc ::cphTitrator::buildTitrator {moveInfo} {
 # 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?
+# cotitration      if 2 residues are involved, can they cotitrate?
 # 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
@@ -333,6 +388,7 @@ proc ::cphTitrator::cphTitratorGet {attr {moveLabel {}}} {
         numsteps -
         weight -
         protonTransfer -
+        cotitration -
         successes -
         attempts -
         proposalCmd {
@@ -349,12 +405,8 @@ proc ::cphTitrator::cphTitratorGet {attr {moveLabel {}}} {
                 split $moveLabel "/"
             }
         }
-        nodefaults {
-            if {$getAll} {
-                cannotGetAll $attr
-            } else {
-                getMoveNoDefaults $moveLabel
-            }
+        archive {
+            getMinimalMoveDict
         }
         default {
             abort "cphTitratorGet: Invalid attribute $attr"
@@ -400,32 +452,33 @@ proc ::cphTitrator::getAllSegresidnameLists {} {
     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]
+# ----------------------------
+# Special getter for archiving
+# ----------------------------
+# Return a minimal version of moveDict. This removes inferred values such as
+# command names as well as specific settings that now match the defaults
+# (since these can change after a restart).
+#
+proc ::cphTitrator::getMinimalMoveDict {} {
+    variable ::cphTitrator::defaultMoveParams
 
-    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
+    set minMoveDict [dict create]
+    dict set minMoveDict maxProposalAttempts [cphTitrator get maxAttempts]
+    dict set minMoveDict default $defaultMoveParams
+    foreach moveLabel [cphTitrator get moveLabels] {
+        foreach attr [dict keys $defaultMoveParams] {
+            set defaultValue [cphTitrator get $attr default]
+            set thisValue [cphTitrator get $attr $moveLabel]
+            if {$thisValue != $defaultValue} {
+                dict set minMoveDict $moveLabel $attr $thisValue
+            }
+        }
+        dict set minMoveDict $moveLabel attempts\
+                [cphTitrator get attempts $moveLabel]
+        dict set minMoveDict $moveLabel successes\
+                [cphTitrator get successes $moveLabel]
     }
-    dict set retDict attempts [cphTitrator get attempts $moveLabel]
-    dict set retDict successes [cphTitrator get successes $moveLabel]
-
-    return $retDict
+    return $minMoveDict
 }
 
 # =============================================================================
@@ -444,6 +497,7 @@ proc ::cphTitrator::getMoveNoDefaults {moveLabel} {
 # 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?
+# cotitration      if 2 residues are involved, can they cotitrate?
 # 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
@@ -465,7 +519,8 @@ proc ::cphTitrator::cphTitratorSet {attr moveLabel value} {
         numsteps -
         successes -
         attempts -
-        protonTransfer { ;# Require integer or boolean argument.
+        protonTransfer -
+        cotitration { ;# Require integer or boolean argument.
             set value [expr {int($value)}]
             if {$setDefault} {
                 dict set defaultMoveParams $attr $value
index 68606da..83ba5f7 100644 (file)
@@ -10,7 +10,7 @@ source [file join [file dirname [info script]] "cphpsfgen.tcl"]
 source [file join [file dirname [info script]] "namdmcmc.tcl"]
 
 namespace eval ::cphSystem {
-    #   The core information of the cphSystem is the residueDict, which stores
+    #   The core information of the cphSystem is the resDict, which stores
     # the unique and mutable information for each instance of a titratable
     # residue. The keys are "segresidnames" of the form:
     #
@@ -18,18 +18,18 @@ namespace eval ::cphSystem {
     #
     # Note the <resname> is a _titratable residue_ name and may or may not be
     # the same as a recognizable residue from the force field. For example,
-    # amino acid termini are titratable residues. The residueDefDict stores
+    # amino acid termini are titratable residues. The resDefDict stores
     # non-unique and immutable information for each residue definition. The
     # keys are residue names (e.g. ASP, HIS, etc.). Generally, when looking up
     # information for a specific residue, one simply needs look up the name in
-    # residueDict and then look up the info in residueDefDict.
+    # resDict and then look up the info in resDefDict.
     #
     # Note that a <segid>:<resid> combination can correspond to _multiple_
     # residue names. This permits non-overlapping "sub-residues," such as amino
     # acid termini, which titrate independently of the sidechain.
     #
-    variable residueDict [dict create]
-    variable residueDefDict [dict create]
+    variable resDict [dict create]
+    variable resDefDict [dict create]
 
     namespace export cphSystem
 }
@@ -90,6 +90,9 @@ proc ::cphSystem::cphSystem {action args} {
                 protonTransfer {
                     proposeProtonTransfer {*}$newArgs
                 }
+                cotitration {
+                    proposeCotitration {*}$newArgs
+                }
                 default {
                     abort "Invalid proposal type $type"
                 }
@@ -201,13 +204,11 @@ proc ::cphSystem::updateStates {accept segresidnameList} {
 # to confirm that a titration was in fact found.
 #
 proc ::cphSystem::proposeResidueTitration {segresidname} {
-    variable ::cphSystem::resDefDict
     set possibleStates [cphSystem get trialStateList $segresidname]
-    set resname [cphSystem get resname $segresidname]
     set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
     while {true} {
         lassign [choice $possibleStates] state
-        set occ [dict get $resDefDict $resname states $state]
+        set occ [cphSystem get occupancy $segresidname $state]
         set numProtonsTrial [expr {lsum($occ)}]
         if {$numProtons != $numProtonsTrial} {
             cphSystem set trialState $segresidname $state
@@ -226,16 +227,11 @@ proc ::cphSystem::proposeResidueTitration {segresidname} {
 # 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]
     set numProtons2 [expr {lsum([cphSystem get occupancy $segresidname2])}]
+    set possibleStates1 [cphSystem get trialStateList $segresidname1]
     set possibleStates2 [cphSystem get trialStateList $segresidname2]
-    set resname2 [cphSystem get resname $segresidname2]
     set dn [expr {$numProtons1 - $numProtons2}]
-
     if {$dn == 0.0} {     
         # No proton transfer is possible.
         return 1
@@ -243,12 +239,13 @@ proc ::cphSystem::proposeProtonTransfer {segresidname1 segresidname2} {
     # transfer from 1 to 2
     while {true} {
         lassign [choice $possibleStates1] state1
-        set occ1 [dict get $resDefDict $resname1 states $state1]
+        set occ1 [cphSystem get occupancy $segresidname1 $state1]
         set numProtonsTrial1 [expr {lsum($occ1)}]
         lassign [choice $possibleStates2] state2
-        set occ2 [dict get $resDefDict $resname2 states $state2]
+        set occ2 [cphSystem get occupancy $segresidname2 $state2]
         set numProtonsTrial2 [expr {lsum($occ2)}]
-        if {$dn == [expr {$numProtonsTrial2 - $numProtonsTrial1}]} {
+        set dnTrial [expr {$numProtonsTrial2 - $numProtonsTrial1}]
+        if {$dn == $dnTrial} {
             cphSystem set trialState $segresidname1 $state1
             cphSystem set trialState $segresidname2 $state2
             return 0
@@ -258,6 +255,44 @@ proc ::cphSystem::proposeProtonTransfer {segresidname1 segresidname2} {
     return -1
 }
 
+# ::cphSystem::proposeCotitration
+#
+# Propose a new trial state requiring cotitration - i.e. two residues
+# concurrently changing from protonated to deprotonated or vice versa.
+#
+# This returns true if no cotitration or else a negative error code.
+#
+proc ::cphSystem::proposeCotitration {segresidname1 segresidname2} {
+    set maxAttempts 10
+
+    set numProtons1 [expr {lsum([cphSystem get occupancy $segresidname1])}]
+    set numProtons2 [expr {lsum([cphSystem get occupancy $segresidname2])}]
+    set possibleStates1 [cphSystem get trialStateList $segresidname1]
+    set possibleStates2 [cphSystem get trialStateList $segresidname2]
+    while {true} {
+        lassign [choice $possibleStates1] state1
+        set occ1 [cphSystem get occupancy $segresidname1 $state1]
+        set numProtonsTrial1 [expr {lsum($occ1)}]
+        set dn1 [expr {$numProtonsTrial1 - $numProtons1}]
+        lassign [choice $possibleStates2] state2
+        set occ2 [cphSystem get occupancy $segresidname2 $state2]
+        set numProtonsTrial2 [expr {lsum($occ2)}]
+        set dn2 [expr {$numProtonsTrial2 - $numProtons2}]
+        if {$dn1 == $dn2} {
+            cphSystem set trialState $segresidname1 $state1
+            cphSystem set trialState $segresidname2 $state2
+            return 0
+        }
+        incr attempts
+        if {$attempts >= $maxAttempts} {
+            # This probably implies that no cotitration exists.
+            return 1
+        }
+    }
+    # 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 
@@ -285,14 +320,14 @@ proc ::cphSystem::proposeProtonTransfer {segresidname1 segresidname2} {
 #
 # For P(k>0, N) = 0.999, this gives N ~= 10 (the default).
 #
-proc ::cphSystem::proposeResidueTautomerization {segresidname {maxAttempts 10}} {
-    variable ::cphSystem::resDefDict
+proc ::cphSystem::proposeResidueTautomerization {segresidname} {
+    set maxAttempts 10
+
     set possibleStates [cphSystem get trialStateList $segresidname]
-    set resname [cphSystem get resname $segresidname]
     set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
     while {true} {
         lassign [choice $possibleStates] state
-        set occ [dict get $resDefDict $resname states $state]
+        set occ [cphSystem get occupancy $segresidname $state]
         set numProtonsTrial [expr {lsum($occ)}]
         if {$numProtonsTrial == $numProtons} {
             cphSystem set trialState $segresidname $state
@@ -308,99 +343,127 @@ proc ::cphSystem::proposeResidueTautomerization {segresidname {maxAttempts 10}}
     return -1
 }
 
-# ::cphSystem::computeInherentAcceptance
+# ::cphSystem::computeInherentLogRatio
 #
-# Compute the (reduced) energy difference for a Monte Carlo move based on the 
-# given segresidname, its current and trial state, and the given pH.
+# Compute the (common) log ratio* of inherent (i.e. pH-dependent) weights of
+# the given states within the given residue.
 #
-# The proposal energy is based exclusively on the "intrinsic" pKa and the
-# change in the protonation vector. There are two cases: 1) tautomerizations
-# and 2) titrations, the former of which is not a proper pKa as the move does
-# not depend on pH (since no protons are entering or leaving the bath).
+# Because pKa is strictly pairwise, the relative unnormed weights (Q) are most
+# clearly computed as ratios (and as log ratios, for numerical stability).
+# There are two cases: 1) tautomerizations and 2) titrations, the former is not
+# a proper pKa since the pH dependence cancels (i.e. no protons are entering or
+# leaving the bath).
 #
 # case 1, n = n':
 #
-#     P(s --> s')
-#     ----------- = 10^[sgn(s' - s) pKa_i(s, s')]
-#     P(s' --> s)
+#     P(l')
+#     ----- = 10^[sgn(l, l') pKa_i(l, l')]
+#     P(l)
 #
 # case 2, n != n':
 #     
-#     P(s --> s')
-#     ----------- = 10^[sgn(n' - n) pKa_i(s, s') - (n' - n)pH]
-#     P(s' --> s)
-#
-# where s and s' are the current and trial state indices, respectively, with
-# number of protons (i.e. sum of the occupation vector elements) n and n',
-# respectively. By convention, pKa_i(s, s') = pKa_i(s', s) and the antisymmetry
-# of adding vs deleting protons is accounted for by the sgn function. Note
-# that for tautomers, the sgn is computed by the _state index_, not the number
-# of protons - this is also an arbitrary internal convention.
-#
-# The above ratios are correctly sampled by a Metropolis sampling of the form:
+#     P(l')
+#     ----- = 10^{-[sgn(n' - n) pKa_i(l, l') - (n' - n)pH]}
+#     P(l)
 #
-#     P(s --> s') = min{1, e^[-du(s, s')]},
+# where l and l' are occupation vectors for the given states and n and n' are
+# the number of protons in each state (i.e. the sums of the occupation vector
+# elements). By convention, pKa_i(l, l') = pKa_i(l', l) and the antisymmetry
+# of adding vs deleting protons is accounted for by the sgn function. Note 
+# that, for tautomers, l and l' are converted to indices and the sgn is
+# computed by an arbitrary convention.
 #
-# where du is either of the the exponents above times -ln(10).
+# Arguments:
+# ----------
+# pH : float
+#   The pH value to evaluate the ratio at.
+# segresidname : string
+#   Residue specification as "<segid>:<resid>:<resname>" 
+# statei : string
+#   The state whose weight is in the numerator
+# statej : string
+#   The state whose weight is in the denomenator
 #
-proc ::cphSystem::computeInherentAcceptance {pH segresidname} {
-    set l [cphSystem get occupancy $segresidname]
-    set lp [cphSystem get trialOccupancy $segresidname]
-    set dn [expr {lsum($lp) - lsum($l)}]
-    set s [occupancy2Index $l]
-    set sp [occupancy2Index $lp]
-    set ssp [index2flatindex $s $sp]
-    set pKai [lindex [cphSystem get pKaiPair $segresidname] $ssp]
+# Returns:
+# --------
+# logRatio : float
+#   The natural log of the ratio Qj/Qi, i.e. ln(Qj/Qi)
+#
+proc ::cphSystem::computeInherentLogRatio {pH segresidname statei statej} {
+    set multij [cphSystem get multiplicity $segresidname $statej]
+    if {$statei == $statej} {
+        # dn = pKai = 0 exactly, but lookup will fail.
+        return [expr {-log($multij)}]
+    }
+
+    set li [cphSystem get occupancy $segresidname $statei]
+    set lj [cphSystem get occupancy $segresidname $statej]
+    set dn [expr {lsum($lj) - lsum($li)}]
+    set si [occupancy2Index $li]
+    set sj [occupancy2Index $lj]
+    set sij [index2flatindex $si $sj]
+    set pKai [lindex [cphSystem get pKaiPair $segresidname] $sij]
     if {$dn == 0.0} {
         # tautomerization: "pKa" is positive in direction of lower index
-        set sgn [expr {$sp > $s} ? 1.0 : -1.0]
+        set sgn [expr {$sj > $si} ? 1.0 : -1.0]
     } else {
         # titration: pKa is positive in direction of fewer protons
         set sgn [expr {$dn > 0} ? 1.0 : -1.0]
     }
-    return [expr {-$::LN10*($sgn*$pKai - $dn*$pH)}] 
+    return [expr {$::LN10*($sgn*$pKai - $dn*$pH) - log($multij)}]
 }
 
-# ::cphSystem::computeInherentNormedWeights
+# ::cphSystem::computeInherentAcceptance
+#
+# Compute the (reduced) energy difference for a Monte Carlo move based on the 
+# given segresidname, its current and trial state, and the given pH.
+#
+# The Metropolis criterion can be expressed as:
 #
-# Compute the normalized inherent pKa weights of all states.
+#     P(l --> l') = min{1, Q'/Q} = min{1, e^[-du(l, l')]},
+#
+# where
+#
+#     du = -ln(Q'/Q)
+#
+proc ::cphSystem::computeInherentAcceptance {pH segresidname} {
+    set statei [cphSystem get state $segresidname]
+    set statej [cphSystem get trialState $segresidname]
+    set du [computeInherentLogRatio $pH $segresidname $statei $statej]
+    set multij [cphSystem get multiplicity $segresidname $statej]
+    # Note the sign flip for compatibility with Metropolis acceptance!
+    return [expr {-($du + log($multij))}]
+}
+
+# ::cphSystem::computeInherentNormedWeights
 #
-# This is a bit tricky, bc only the relative unnormed weights are known. This
-# is solved by taking all weights relative to the current state and assigning
-# that state an unnormed weight of one.
+# Compute the normalized inherent pKa weights (i.e. probability mass function,
+# or PMF) of all states within the given residue at the given pH.
 #
-# NB: The order of states is randomly shuffled here, because the current state
-#   is omitted from the stateList and given the index 0 in the list of weights
-#   (hence the weight list is one element longer!).
+# The resulting PMF can be directly sampled in order to choose a new state ala
+# independence sampling. However, this is often not an efficent scheme for
+# Markov chain Monte Carlo and the output can also be used for other approaches
+# such as Metropolized independence sampling.
 #
 proc ::cphSystem::computeInherentNormedWeights {pH segresidname} {
-    set l [cphSystem get occupancy $segresidname]
-    set s [occupancy2Index $l]
-    set resname [cphSystem get resname $segresidname]
-    set stateList [cphSystem get trialStateList $segresidname]
-    # We implicitly reference against the current state, so its unnormed weight
-    # is exactly one.
-    set logQs [list 1.0]
-    set logQMax 1.0
+    set currState [cphSystem get state $segresidname]
+    set stateList [cphSystem get stateList $segresidname] 
+
+    # Note that all ratios are taken wrt the current state and are thus
+    # equivalent to unnormed weights. This would NOT be correct if a different
+    # state were in the denomenator each time.
+    #
+    set logQs [list]
+    set logQMax 0.0
     foreach state $stateList {
-        set lp [state2Occupancy $resname $state]
-        set dn [expr {lsum($lp) - lsum($l)}]
-        set sp [occupancy2Index $lp]
-        set ssp [index2flatindex $s $sp]
-        set pKai [lindex [cphSystem get pKaiPair $segresidname] $ssp]
-        if {$dn == 0.0} {
-            # tautomerization: "pKa" is positive in direction of lower index
-            set sgn [expr {$sp > $s} ? 1.0 : -1.0]
-        } else {
-            # titration: pKa is positive in direction of fewer protons
-            set sgn [expr {$dn > 0} ? 1.0 : -1.0]
-        }
-        lappend logQs [expr {$::LN10*($sgn*$pKai - $dn*$pH)}]
-        if {[lindex $logQs end] > $logQMax} {
-            set logQMax [lindex $logQs end]
+        set logQ [computeInherentLogRatio $pH $segresidname $currState $state]
+        lappend logQs $logQ
+        if {$logQ > $logQMax} {
+            set logQMax $logQ
         }
     }
-    return [list [normalizeLogWeights $logQs $logQMax] $stateList]
+    set normedWeights [normalizeLogWeights $logQs $logQMax]
+    return [list $stateList $normedWeights]
 }
 
 # ::cphSystem::computeSwitchAcceptance
@@ -436,6 +499,7 @@ proc ::cphSystem::computeSwitchAcceptance {segresidnameList} {
         set dG [lindex [cphSystem get dGPair $segresidname] $ssp]
         set pKa [lindex [cphSystem get pKaPair $segresidname] $ssp]
         set pKai [lindex [cphSystem get pKaiPair $segresidname] $ssp]
+        set trialState [cphSystem get trialState $segresidname]
         if {$dn == 0.0} {
             # tautomerization: "pKa" is positive in direction of lower index
             set sgn [expr {$sp > $s} ? 1.0 : -1.0]
@@ -566,7 +630,7 @@ proc ::cphSystem::buildSystem {resDefs excludeList} {
     # Check for residue aliases (i.e. titratable residues whose states may span
     # multiple residue definitions, such as HIS).
     set resAliases [dict create]
-    foreach {resname resDef} $resDefDict {
+    dict for {resname resDef} $resDefDict {
         if {![dict exists $resDef aliases]} continue
         dict set resAliases $resname [dict get $resDef aliases]
     }
@@ -758,8 +822,7 @@ proc ::cphSystem::checkResDefs {resDefs} {
 # Arguments:
 # ----------
 # segresidname : string
-#   residue specification as "<segid>:<resid>" - this is the same syntax as for
-#   the regular psfgen patch command.
+#   residue specification as "<segid>:<resid>:<resname>"
 # state : string
 #   state to assign 
 #
@@ -777,39 +840,34 @@ proc ::cphSystem::assignState {segresidname state} {
 
 # ::cphSystem::randomizeState
 #
-# Assign a protonation state by randomly choosing a state and performing MC 
-# moves based on the pH until a new state is accepted.
-#
-#   This is a bit sneaky. Implementing a proper independence sampling of all of
-# the states seems overkill. Rather, an originating state is chosen uniformly 
-# from all possible states and the state to be assigned is chosen by a pairwise
-# pH/pKai based Metropolis criterion. If this is rejected, then a new 
-# originating state is chosen. This should be equivalent to independence
-# sampling for many such trials. 
+# Assign a protonation state by independence sampling at the given pH.
 #
 # Arguments:
 # ----------
 # segresidname : string
-#   residue specification as "<segid>:<resid>" - this is the same syntax as for
-#   the regular psfgen patch command.
+#   residue specification as "<segid>:<resid>:<resname>"
 # pH : float
-#    pH value at which to assign protonation states
+#   pH value at which to assign protonation states
 #
 # Returns:
 # --------
 # None 
 #
 proc ::cphSystem::randomizeState {segresidname pH} {
-    while {true} {
-        set states [cphSystem get stateList $segresidname]
-        cphSystem set state $segresidname [lindex [choice $states] 0]
-        cphSystem propose titration $segresidname
-        set du [cphSystem compute inherent $pH $segresidname]
-        if {[metropolisAcceptance $du]} { 
-            return 0
-        }
-    }
-    return -1
+    # Note that computeInherentNormedWeights requires a state to be set for
+    # normalization purposes. Therefore choose a random state first and then
+    # rigorously sample the distribution.
+    set states [cphSystem get stateList $segresidname]
+    lassign [choice $states] state
+    cphSystem set state $segresidname $state
+    lassign [computeInherentNormedWeights $pH $segresidname] states weights
+    lassign [choice $states $weights] state
+    cphSystem set state $segresidname $state
+    # Use the same hack as in assignState.
+    cphSystem propose titration $segresidname
+    cphSystem set state $segresidname [cphSystem get trialState $segresidname]
+    cphSystem set trialState $segresidname $state
+    return 0
 }
 
 # =============================================================================
@@ -845,6 +903,7 @@ proc ::cphSystem::randomizeState {segresidname pH} {
 # occupancy          occupancy vector for the current state
 # trialOccupancy     occupancy vector for the trial state
 # stateList          all possible states
+# equivStateList     all states equivalent to the given state
 # trialStateList     all possible trial states (not the current state)
 # dGPair             pair dG for current/trial states
 # pKaPair            pair pKa for current/trial states
@@ -902,7 +961,12 @@ proc ::cphSystem::cphSystemGet {attr {segresidname {}} args} {
             if {$getAll} {
                 getAllOccupancy
             } else {
-                getOccupancy $segresidname
+                if {[llength $args]} {
+                    lassign $args state
+                    getOccupancy $segresidname $state
+                } else {
+                    getCurrOccupancy $segresidname
+                }
             }
         }
         trialOccupancy {
@@ -920,6 +984,32 @@ proc ::cphSystem::cphSystemGet {attr {segresidname {}} args} {
                 dict keys [dict get $resDefDict $resname states]
             }
         }
+        equivStateList {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                if {[llength $args]} {
+                    lassign $args state
+                } else {
+                    # Default to current state.
+                    set state [getResAttr state $segresidname]
+                }
+                getEquivStates $segresidname $state
+            }
+        }
+        multiplicity {
+            if {$getAll} {
+                cannotGetAll $attr
+            } else {
+                if {[llength $args]} {
+                    lassign $args state
+                } else {
+                    # Default to current state.
+                    set state [getResAttr state $segresidname]
+                }
+                llength [getEquivStates $segresidname $state]
+            }
+        }
         trialStateList {
             if {$getAll} {
                 cannotGetAll $attr
@@ -1044,20 +1134,22 @@ proc ::cphSystem::state2Occupancy {resname state} {
     return [dict get $::cphSystem::resDefDict $resname states $state]
 }
 
-# NB: When returning all occupancies, only the current state can be used and
-#   the list is flattened.
-proc ::cphSystem::getOccupancy {segresidname} {
+proc ::cphSystem::getOccupancy {segresidname state} {
     variable ::cphSystem::resDict
     set resname [dict get $resDict $segresidname resname]
-    set state [dict get $resDict $segresidname state]
     return [state2Occupancy $resname $state]
 }
 
+proc ::cphSystem::getCurrOccupancy {segresidname} {
+    variable ::cphSystem::resDict
+    set state [dict get $resDict $segresidname state]
+    return [getOccupancy $segresidname $state]
+}
+
 proc ::cphSystem::getTrialOccupancy {segresidname} {
     variable ::cphSystem::resDict
-    set resname [dict get $resDict $segresidname resname]
     set state [dict get $resDict $segresidname trialState]
-    return [state2Occupancy $resname $state]
+    return [getOccupancy $segresidname $state] 
 }
 
 proc ::cphSystem::getAllOccupancy {} {
@@ -1072,6 +1164,34 @@ proc ::cphSystem::getAllOccupancy {} {
     return $retList
 }
 
+# ::cphSystem::getEquivStates
+#
+# Return the list of states equivalent to the given state of the given residue. # This includes the given state such that the length of the list is exactly the
+# degeneracy.
+#
+proc ::cphSystem::getEquivStates {segresidname testState} {
+    set l [cphSystem get occupancy $segresidname $testState]
+    set s [occupancy2Index $l]
+
+    set equivStateList [list]
+    foreach state [cphSystem get stateList $segresidname] {
+        if {$state == $testState} {
+            lappend equivStateList $state
+            continue
+        }
+
+        set lp [cphSystem get occupancy $segresidname $state]
+        set dn [expr {lsum($lp) - lsum($l)}]
+        set sp [occupancy2Index $lp]
+        set ssp [index2flatindex $s $sp]
+        set pKa [lindex [cphSystem get pKaPair $segresidname] $ssp]
+        if {$dn == 0.0 && $pKa == 0.0} {
+            lappend equivStateList $state
+        }
+    }
+    return $equivStateList
+}
+
 # ----------------------
 # Other helper functions
 # ----------------------
index 6f707f2..ce0e21d 100644 (file)
@@ -62,15 +62,15 @@ proc ::namdcph::cphRun {numsteps {numcycles 1}} {
         # (2) Propose a move from the full move set.
         #
         lassign [cphTitrator propose $SystempH] accept swNumsteps\
-                segresidnameList nattempts
+                segresidnameList
         # (3) At this point we have either selected a switch or rejected a
         # whole bunch of moves. If the former, perform the switch.
         #
         if {!$accept} {
-            cphPrint "All proposals rejected ($nattempts total)."
+            cphPrint "Proposal rejected."
             set runArgs [list norepeat $numsteps]
         } else {
-            cphPrint "Proposal accepted ($nattempts), attemping a switch."
+            cphPrint "Proposal accepted, attemping a switch."
             set accept [runSwitch $swNumsteps $segresidnameList]
             set runArgs [list $numsteps]
             # Only accumulate statistics for attempted switches.
@@ -153,6 +153,7 @@ proc ::namdcph::checkResidueDefinitions {resnames} {
 proc ::namdcph::testResidueSwitch {segresidname state0 state1} {
     cphWarnExpt
     pH 7.0
+    cphNumstepsPerSwitch 20
     cphForceConstant 0.0
     # Initialize NAMD and build a constant pH enabled PSF.
     cphSetResidueState $segresidname $state0
@@ -394,6 +395,13 @@ proc ::namdcph::cphProposalWeight {args} {
     return
 }
 
+# ::namdcph::cphProposeProtonTransfer
+#
+# Add additional "proton transfer" moves to the move set. These new moves are
+# only attempted if a proton can be swapped between the given residues. The
+# two residues are given by their normal <segid>:<resid>:<resname> designation
+# except that they are further separated by a "/".
+#
 proc ::namdcph::cphProposeProtonTransfer {args} {
     variable ::namdcph::moveInfo
     foreach moveLabel $args {
@@ -402,6 +410,22 @@ proc ::namdcph::cphProposeProtonTransfer {args} {
     return
 }
 
+# ::namdcph::cphProposeCotitration
+#
+# Add additional "co-titration" moves to the move set. These new moves are only
+# attempted if two residues are both protonated or deprotonated. They will then
+# attempt to change states concurrently. The two residues are given by their
+# normal <segid>:<resid>:<resname> designation except that they are further
+# separated by a "/".
+#
+proc ::namdcph::cphProposeCotitration {args} {
+    variable ::namdcph::moveInfo
+    foreach moveLabel $args {
+        dict set moveInfo $moveLabel cotitration 1
+    }
+    return
+}
+
 # ::namdcph::cphMaxProposalAttempts
 #
 # Number of attempted MC proposals from each move set before giving up.
@@ -414,6 +438,7 @@ proc ::namdcph::cphMaxProposalAttempts {maxAttempts} {
     checkIsNumeric cphMaxProposalAttempts $maxAttempts
     variable ::namdcph::moveInfo
     dict set moveInfo maxProposalAttempts [expr {int($maxAttempts)}]
+    cphWarn "cphMaxProposalAttempts is now deprecated and set to 1!"
     return
 }
 
@@ -528,7 +553,7 @@ proc ::namdcph::runSwitch {numsteps segresidnameList} {
     storeEnergies
     set DeltaE [cphSystem compute switch $segresidnameList]
     set Work [expr {$::energyArray(CUMALCHWORK) + $DeltaE}]
-    set ReducedWork [expr {$Work / [::kBT]}] 
+    set ReducedWork [expr {$Work / [::kBT]}]
     set accept [metropolisAcceptance $ReducedWork]
     set tmp [printProposalSummary $segresidnameList]
     cphPrint [format "%s WorkCorr % 10.4f CorrWork % 10.4f"\
@@ -781,21 +806,7 @@ proc ::namdcph::writeRestart {args} {
     }
     lappend rstrtList "\"states\":[json::list2json [cphSystem get state] 1]"
     lappend rstrtList "\"pKais\":[json::list2json [cphSystem get pKai]]"
-
-    set minMCDict [dict create]
-    dict set minMCDict MCmoves maxProposalAttempts\
-            [cphTitrator get maxAttempts]
-    dict set minMCDict MCmoves default numsteps\
-            [cphTitrator get numsteps default]
-    dict set minMCDict MCmoves default weight [cphTitrator get weight default]
-    foreach moveLabel [cphTitrator get moveLabels] {
-        set thisMoveInfo [cphTitrator get nodefaults $moveLabel]
-        if {![dict size $thisMoveInfo]} continue
-        dict set minMCDict MCmoves $moveLabel $thisMoveInfo
-    }
-    # Note that dict2json returns curly braces on either end which need to be
-    # discarded since we are hacking the format into another dict.
-    lappend rstrtList [string range [json::dict2json $minMCDict] 1 end-1] 
+    lappend rstrtList "\"MCmoves\":[json::dict2json [cphTitrator get archive]]"
 
     # Write everything to file.
     namdFileBackup $restartFilename
@@ -911,12 +922,12 @@ proc ::namdcph::initialize {} {
     variable ::namdcph::excludeList
     variable ::namdcph::SystempH
 
-    getThermostat
-    getBarostat
     callback energyCallback
     # 1) Set up alchemical keywords and run startup. 
     #
     initializeAlch
+    getThermostat
+    getBarostat
 
     # 2) Rebuild the PSF with dummy protons and modify protonation states as 
     # needed. Build the residue definitions, assign states to each residue, and 
@@ -931,13 +942,14 @@ proc ::namdcph::initialize {} {
         lassign [readRestart] stateList pKaiList
         lassign [list 0.0 false] temp buildH
     } else {
-        lassign [list [$::thermostatTempCmd] true] temp buildH
+        lassign [list [$::thermostatTempCmd] false] temp buildH
     }
     cphSystem build [readConfig] $excludeList
 
     # Now that we've built the system, we can allocate the state and pKa
     # information from the restart (if present).
-    if {[info exists stateList] && [info exists pKaiList]} {
+    if {[info exists stateList] && [llength $stateList] &&
+        [info exists pKaiList] && [llength $pKaiList]} {
         set segresidnameList [cphSystem get segresidnames]
 
         if {[llength $stateList] != [llength $pKaiList]} {
@@ -1113,14 +1125,38 @@ proc ::namdcph::printSystemSummary {} {
 #
 proc ::namdcph::printTitratorSummary {} {
     set StarBar "*************************************************"
+    set moveLabels [cphTitrator get moveLabels]
+    set numMoves [llength $moveLabels]
+
     cphPrint $StarBar
     cphPrint "CONSTANT pH neMD/MC MOVES:"
     cphPrint "MAX. ATTEMPTS PER CYCLE: [cphTitrator get maxAttempts]"
-    cphPrint [format "%-25s : %8s %12s" "move label" numsteps weight]
-    foreach moveLabel [cphTitrator get moveLabels] {
+
+    cphPrint "ONE RESIDUE MOVES"
+    cphPrint [format "%-14s : %8s %8s" "move label" numsteps weight]
+    foreach moveLabel $moveLabels {
+        if {[llength [split $moveLabel "/"]] != 1} continue
+        incr numEntries
+        set numsteps [cphTitrator get numsteps $moveLabel]
+        set weight [cphTitrator get weight $moveLabel]
+        cphPrint [format "%-14s : % 8d % 8.2f" $moveLabel $numsteps $weight]
+    }
+    if {$numEntries == $numMoves} {
+        cphPrint $StarBar
+        return
+    }
+
+    cphPrint "TWO RESIDUE MOVES : proton transfer (PT) : co-titration (CT)"
+    cphPrint [format "%-29s : %8s %8s %2s %2s"\
+            "move label" numsteps weight PT CT]
+    foreach moveLabel $moveLabels {
+        if {[llength [split $moveLabel "/"]] != 2} continue
         set numsteps [cphTitrator get numsteps $moveLabel]
         set weight [cphTitrator get weight $moveLabel]
-        cphPrint [format "%-25s : % 8d % 12.2f" $moveLabel $numsteps $weight]
+        set PT [cphTitrator get protonTransfer $moveLabel]
+        set CT [cphTitrator get cotitration $moveLabel]
+        cphPrint [format "%-29s : % 8d % 8.2f %2d %2d"\
+                $moveLabel $numsteps $weight $PT $CT]
     }
     cphPrint $StarBar
     return