Constant-pH MD: Revisions and new functionality 94/3794/3
authorradakb <brian.radak@gmail.com>
Thu, 1 Mar 2018 21:56:16 +0000 (16:56 -0500)
committerJim Phillips <jim@ks.uiuc.edu>
Fri, 2 Mar 2018 18:26:35 +0000 (12:26 -0600)
Minor changes:
 - Clean up of documentation
 - Variable name segresid changed to segresidname throughout (see below)
 - Complete re-hab or residue testing facility. The old system was based
   on the old, directional pattern of patch definition. This is no
   longer necessary and way more robust. Testing that the correct
   patches are defined is still very easy, but a new convenience
   function is added for testing the correctness of the patches is also
   added. Scripting a test with this is very involved...
 - Simplification and performance enhancements for the experimental
   cphAnalyzeForce command.
 - Using minimization via cphNumMinSteps now results in reallocation of
   velocities when done. Either the temperature or thermostat
   temperature are used.
 - The generic namdtcl interface now includes a "::kBT" proc that
   returns the value of kT in the current thermostat.

Major changes:
 - A residue, in the topology sense, can now be divided into multiple
   titratable residues. Any titratable residue which does not have the
   same name as the topological residue is considered a "sub-residue"
   and requires special definition in the configuration file. Note that
   some residues may ONLY contain sub-residues, such as a C-terminal
   alanine or the phosphates in a lipid head group. In order to
   differentiate sub-residues from the regular residues, we expand the
   <segid>:<resid> selection to <segid>:<resid>:<resname>, where
   <resname> may be a sub-residue name. This new label (the
   "segresidname") is the only guaranteed unique selector for a
   titratable residue.
   IMPORTANT! Even if a residue has no sub-residues, the new
   segresidname syntax MUST be used for any keyword that requires a
   residue selection. THIS BREAKS OLD NAMD CONFIGURATION FILES THAT
   SELECT RESIDUES!
 - The initial "run 0" command to kick-start the simulation has been
   replaced by the "startup" command and thus eliminates the spurious
   ENERGY: entry when a simulation begins. As a consequence, the old
   neMD/MC cycle sequence of MC/MD has now been flipped to MD/MC. This
   should not cause any change to statistics, but output will look
   slightly different.

Change-Id: Iad97be3aa5cfd1ee58228703f1c92b78987886b0

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

index 9354cb0..2616f02 100644 (file)
@@ -19,9 +19,9 @@ namespace eval ::namdcphwrapper {
             cphRestartFile cphRestartFreq\
             cphForceConstant cphMDBasename cphSwitchBasename\
             cphMaxProposalAttempts cphNumMinSteps cphProposalWeight\
-            cphProposeProtonTransfer\
+            cphExcludeResidue cphProposeProtonTransfer\
             cphNumstepsPerSwitch\
-            testResidue cphAnalyzeForce cphExcludeResidue
+            checkResidueDefinitions testResidueSwitch cphAnalyzeForce
 
     # Old command name for legacy reasons.
     proc runcph {args} {
index 9b5978d..a0cc9c5 100644 (file)
@@ -43,31 +43,33 @@ proc ::cphTitrator::buildTitrator {systempH moveInfo} {
         set maxAttempts [cphSystem get numresidues]
     }
     #   The default move set is to switch each residue independently (i.e. the
-    # move label is just the segresid for that residue). Settings from 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 segresid [cphSystem get segresids] {
-        dict set MoveSet $segresid proposalCmd "proposeResidueMove $segresid"
-        dict set MoveSet $segresid segresidList $segresid
-        dict set MoveSet $segresid weight\
-                [dictPopOrDefault moveInfo $segresid weight]
-        dict set MoveSet $segresid numsteps\
-                [dictPopOrDefault moveInfo $segresid numsteps]
-        dict set MoveSet $segresid successes 0
-        dict set MoveSet $segresid attempts 0
-        if {[dict exists $moveInfo $segresid]\
-            && ![dict size [dict get $moveInfo $segresid]]} {
-            dict unset moveInfo $segresid
+    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 "/"] segresid1 segresid2
+        lassign [split $moveLabel "/"] segresidname1 segresidname2
         dict set MoveSet $moveLabel proposalCmd\
                 "proposeProtonTransferMove $moveLabel" 
-        dict set MoveSet $moveLabel segresidList [list $segresid1 $segresid2]
+        dict set MoveSet $moveLabel segresidnameList\
+                [list $segresidname1 $segresidname2]
         dict set MoveSet $moveLabel weight\
                 [dictPopOrDefault moveInfo $moveLabel weight]
         dict set MoveSet $moveLabel numsteps\
@@ -117,7 +119,7 @@ proc ::cphTitrator::computeWeightSum {} {
 # =============================================================================
 # ::cphTitrator::proposeResidueMove
 #
-# Propose a move invovling a single residue via Metropolized independence
+# 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
@@ -125,9 +127,9 @@ 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 {segresid} {
+proc ::cphTitrator::proposeResidueMove {segresidname} {
    variable ::cphTitrator::pH
-   lassign [cphSystem compute inherentWeights $pH $segresid] weights states 
+   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).
@@ -138,18 +140,19 @@ proc ::cphTitrator::proposeResidueMove {segresid} {
    set du [expr {log((1. - $pj) / $pic)}]
    set accept [metropolisAcceptance $du]
    if {$accept} {
-       cphSystem set trialState $segresid [lindex $states $j]
+       cphSystem set trialState $segresidname [lindex $states $j]
    }
    return $accept
 }
 
 proc ::cphTitrator::proposeProtonTransferMove {moveLabel} {
     variable ::cphTitrator::pH
-    lassign [split $moveLabel "/"] segresid1 segresid2
-    set errCode [::cphSystem::proposeProtonTransfer $segresid1 $segresid2]
+    lassign [split $moveLabel "/"] segresidname1 segresidname2
+    set errCode\
+        [::cphSystem::proposeProtonTransfer $segresidname1 $segresidname2]
     if {$errCode > 0} {
-        set du1 [cphSystem compute inherent $pH $segresid1]
-        set du2 [cphSystem compute inherent $pH $segresid2]
+        set du1 [cphSystem compute inherent $pH $segresidname1]
+        set du2 [cphSystem compute inherent $pH $segresidname2]
         set accept [metropolisAcceptance [expr {$du1 + $du2}]]
     } else {
         set accept 0
@@ -168,16 +171,16 @@ proc ::cphTitrator::proposeMove {} {
         incr nattempts
         # Clear the previous trial if it was rejected.
         if {$nattempts > 1} {
-             cphSystem update $accept $segresidList
+             cphSystem update $accept $segresidnameList
         }
         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 segresidList [dict get $MoveSet $moveLabel segresidList]
+        set segresidnameList [dict get $MoveSet $moveLabel segresidnameList]
     }
-    return [list $accept $numsteps $segresidList $nattempts]
+    return [list $accept $numsteps $segresidnameList $nattempts]
 }
 
 # =============================================================================
index f3d73ab..450d879 100644 (file)
@@ -12,13 +12,22 @@ source [file join [file dirname [info script]] "namdmcmc.tcl"]
 namespace eval ::cphSystem {
     #   The core information of the cphSystem is the residueDict, which stores
     # the unique and mutable information for each instance of a titratable
-    # residue. The keys are "segresids" of the form <segid>:<resid>, just as in
-    # the psfgen patch command. The residueDefDict stores non-unique and
-    # immutable information for each residue definition (or type). The keys are
-    # residue names (e.g. ASP, HIS, etc.). Generally, when looking up 
+    # residue. The keys are "segresidnames" of the form:
+    #
+    # <segid>:<resid>:<resname>.
+    #
+    # 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
+    # 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.
     #
+    # 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]
 
@@ -31,7 +40,21 @@ namespace eval ::cphSystem {
 # ::cphSystem::cphSystem
 #
 # This is the only exported function from the cphSystem namespace and provides
-# a complete interface to
+# a complete interface to the cphSystem namespace. The first argument is always
+# an action (e.g. "get" or "set") and the rest are variable arguments.
+#
+# action           description
+# ---------------- -----------
+# get              see cphSystemGet 
+# set              see cphSystemSet
+# build            see buildSystem
+# initialize       see initializeSystem
+# propose          propose a move of the given type
+# compute          compute various quantities needed for MC
+# update           accept/reject the proposed changes
+# alchemifypsf     apply patches that make the system alchemical
+# dealchemifypsf   remove patches so the system is no longer alchemical
+# initializeState  set the initial state of a residue
 #
 proc ::cphSystem::cphSystem {action args} {
     if {[string match -nocase $action get]} {
@@ -97,19 +120,20 @@ proc ::cphSystem::cphSystem {action args} {
 # =============================================================================
 # proc ::cphSystem::updateStates
 #
-# Update the state of one or more residues with the given <segid>:<resid> 
-# specifications based on acceptance/rejection of the trial states. Reset the 
-# trial states to a null value.
+# Update the state of one or more residues with the given
+# <segid>:<resid>:<resname> specifications based on acceptance/rejection of the
+# trial states. Reset the trial states to a null value.
 #
-proc ::cphSystem::updateStates {accept segresidList} {
+proc ::cphSystem::updateStates {accept segresidnameList} {
     if {$accept} {
-        foreach segresid $segresidList {
-            cphSystem set state $segresid [cphSystem get trialState $segresid]
-            cphSystem set trialState $segresid {}
+        foreach segresidname $segresidnameList {
+            cphSystem set state $segresidname\
+                    [cphSystem get trialState $segresidname]
+            cphSystem set trialState $segresidname {}
         }
     } else {
-        foreach segresid $segresidList {
-            cphSystem set trialState $segresid {}
+        foreach segresidname $segresidnameList {
+            cphSystem set trialState $segresidname {}
         }
     }
     return 0
@@ -123,17 +147,17 @@ proc ::cphSystem::updateStates {accept segresidList} {
 # For consistency with tautomers (see below), this always returns true in order
 # to confirm that a titration was in fact found.
 #
-proc ::cphSystem::proposeResidueTitration {segresid} {
+proc ::cphSystem::proposeResidueTitration {segresidname} {
     variable ::cphSystem::resDefDict
-    set possibleStates [cphSystem get trialStateList $segresid]
-    set resname [cphSystem get resname $segresid]
-    set numProtons [expr {lsum([cphSystem get occupancy $segresid])}]
+    set possibleStates [cphSystem get trialStateList $segresidname]
+    set resname [cphSystem get resname $segresidname]
+    set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
     while {true} {
         set state [choice $possibleStates]
         set occ [dict get $resDefDict $resname states $state]
         set numProtonsTrial [expr {lsum($occ)}]
         if {$numProtons != $numProtonsTrial} {
-            cphSystem set trialState $segresid $state
+            cphSystem set trialState $segresidname $state
             return 1
         }
     }
@@ -141,14 +165,22 @@ proc ::cphSystem::proposeResidueTitration {segresid} {
     return -1
 }
 
-proc ::cphSystem::proposeProtonTransfer {segresid1 segresid2} {
+# ::cphSystem::proposeProtonTransfer
+#
+# 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.
+#
+proc ::cphSystem::proposeProtonTransfer {segresidname1 segresidname2} {
     variable ::cphSystem::resDefDict
-    set numProtons1 [expr {lsum([cphSystem get occupancy $segresid1])}]
-    set possibleStates1 [cphSystem get trialStateList $segresid1]
-    set resname1 [cphSystem get resname $segresid1]
-    set numProtons2 [expr {lsum([cphSystem get occupancy $segresid2])}]
-    set possibleStates2 [cphSystem get trialStateList $segresid2]
-    set resname2 [cphSystem get resname $segresid2]
+    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 possibleStates2 [cphSystem get trialStateList $segresidname2]
+    set resname2 [cphSystem get resname $segresidname2]
     set dn [expr {$numProtons1 - $numProtons2}]
 
     if {$dn == 0.0} {     
@@ -164,8 +196,8 @@ proc ::cphSystem::proposeProtonTransfer {segresid1 segresid2} {
         set occ2 [dict get $resDefDict $resname2 states $state2]
         set numProtonsTrial2 [expr {lsum($occ2)}]
         if {$dn == [expr {$numProtonsTrial2 - $numProtonsTrial1}]} {
-            cphSystem set trialState $segresid1 $state1
-            cphSystem set trialState $segresid2 $state2
+            cphSystem set trialState $segresidname1 $state1
+            cphSystem set trialState $segresidname2 $state2
             return 1
         }
     }
@@ -201,17 +233,17 @@ proc ::cphSystem::proposeProtonTransfer {segresid1 segresid2} {
 #
 # For P(k>0, N) = 0.999, this gives N ~= 10 (the default).
 #
-proc ::cphSystem::proposeResidueTautomerization {segresid {maxAttempts 10}} {
+proc ::cphSystem::proposeResidueTautomerization {segresidname {maxAttempts 10}} {
     variable ::cphSystem::resDefDict
-    set possibleStates [cphSystem get trialStateList $segresid]
-    set resname [cphSystem get resname $segresid]
-    set numProtons [expr {lsum([cphSystem get occupancy $segresid])}]
+    set possibleStates [cphSystem get trialStateList $segresidname]
+    set resname [cphSystem get resname $segresidname]
+    set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
     while {true} {
         set state [choice $possibleStates]
         set occ [dict get $resDefDict $resname states $state]
         set numProtonsTrial [expr {lsum($occ)}]
         if {$numProtonsTrial == $numProtons} {
-            cphSystem set trialState $segresid $state
+            cphSystem set trialState $segresidname $state
             return 1
         }
         incr attempts
@@ -227,7 +259,7 @@ proc ::cphSystem::proposeResidueTautomerization {segresid {maxAttempts 10}} {
 # ::cphSystem::computeInherentAcceptance
 #
 # Compute the (reduced) energy difference for a Monte Carlo move based on the 
-# given <segid>:<resid>, its current and trial state, and the given pH.
+# given segresidname, its current and trial state, and the given pH.
 #
 # The proposal energy is based exclusively on the "intrinsic" pKa and the
 # change in the protonation vector. There are two cases: 1) tautomerizations
@@ -247,7 +279,7 @@ proc ::cphSystem::proposeResidueTautomerization {segresid {maxAttempts 10}} {
 #     P(s' --> s)
 #
 # where s and s' are the current and trial state indices, respectively, with
-# number of protons (i.e. magnitude of the occupation vector) n and n',
+# 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
@@ -259,14 +291,14 @@ proc ::cphSystem::proposeResidueTautomerization {segresid {maxAttempts 10}} {
 #
 # where du is either of the the exponents above times -ln(10).
 #
-proc ::cphSystem::computeInherentAcceptance {pH segresid} {
-    set l [cphSystem get occupancy $segresid]
-    set lp [cphSystem get trialOccupancy $segresid]
+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 $segresid] $ssp]
+    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]
@@ -289,11 +321,11 @@ proc ::cphSystem::computeInherentAcceptance {pH segresid} {
 #   is omitted from the stateList and given the index 0 in the list of weights
 #   (hence the weight list is one element longer!).
 #
-proc ::cphSystem::computeInherentNormedWeights {pH segresid} {
-    set l [cphSystem get occupancy $segresid]
+proc ::cphSystem::computeInherentNormedWeights {pH segresidname} {
+    set l [cphSystem get occupancy $segresidname]
     set s [occupancy2Index $l]
-    set resname [cphSystem get resname $segresid]
-    set stateList [cphSystem get trialStateList $segresid]
+    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]
@@ -303,7 +335,7 @@ proc ::cphSystem::computeInherentNormedWeights {pH segresid} {
         set dn [expr {lsum($lp) - lsum($l)}]
         set sp [occupancy2Index $lp]
         set ssp [index2flatindex $s $sp]
-        set pKai [lindex [cphSystem get pKaiPair $segresid] $ssp]
+        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]
@@ -340,18 +372,18 @@ proc ::cphSystem::computeInherentNormedWeights {pH segresid} {
 # computed, not the temperature of the simulation. See 
 # computeInherentAcceptance for definition of the remaining notation.
 #
-proc ::cphSystem::computeSwitchAcceptance {segresidList} {
+proc ::cphSystem::computeSwitchAcceptance {segresidnameList} {
     set du 0.0
-    foreach segresid $segresidList {
-        set l [cphSystem get occupancy $segresid]
-        set lp [cphSystem get trialOccupancy $segresid]
+    foreach segresidname $segresidnameList {
+        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 dG [lindex [cphSystem get dGPair $segresid] $ssp]
-        set pKa [lindex [cphSystem get pKaPair $segresid] $ssp]
-        set pKai [lindex [cphSystem get pKaiPair $segresid] $ssp]
+        set dG [lindex [cphSystem get dGPair $segresidname] $ssp]
+        set pKa [lindex [cphSystem get pKaPair $segresidname] $ssp]
+        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]
@@ -359,7 +391,7 @@ proc ::cphSystem::computeSwitchAcceptance {segresidList} {
             # titration: pKa is positive in direction of fewer protons
             set sgn [expr {$dn > 0} ? 1.0 : -1.0]
         }
-        set kT [expr {$::BOLTZMANN*[cphSystem get Tref $segresid]}]
+        set kT [expr {$::BOLTZMANN*[cphSystem get Tref $segresidname]}]
         set du [expr {$du + $sgn*($dG - $kT*$::LN10*($pKa - $pKai))}]
     }
     return $du
@@ -372,10 +404,12 @@ proc ::cphSystem::computeSwitchAcceptance {segresidList} {
 #
 # Apply a trial alchemical patch to a residue.
 #
-proc ::cphSystem::alchemifyResidue {segresid frcCons temp {buildH false}} {
-    lassign [cphSystem get alchAtomLists $segresid] l0atoms l1atoms
-    set patch [cphSystem get hybridPatch $segresid]
-    alchPatch $patch $segresid $l0atoms $l1atoms $frcCons $temp $buildH
+proc ::cphSystem::alchemifyResidue {segresidname frcCons temp {buildH false}} {
+    lassign [cphSystem get alchAtomLists $segresidname] l0atoms l1atoms
+    set patch [cphSystem get hybridPatch $segresidname]
+    lassign [split $segresidname ":"] segid resid
+    # NB: alchPatch uses psfgen style selections!
+    alchPatch $patch "$segid:$resid" $l0atoms $l1atoms $frcCons $temp $buildH
     return 0
 }
 
@@ -383,9 +417,11 @@ proc ::cphSystem::alchemifyResidue {segresid frcCons temp {buildH false}} {
 #
 # Remove an alchemical patch from a residue.
 #
-proc ::cphSystem::dealchemifyResidue {segresid} {
-    lassign [cphSystem get alchAtomLists $segresid] l0atoms l1atoms
-    alchUnpatch $segresid $l0atoms $l1atoms
+proc ::cphSystem::dealchemifyResidue {segresidname} {
+    lassign [cphSystem get alchAtomLists $segresidname] l0atoms l1atoms
+    lassign [split $segresidname ":"] segid resid 
+    # NB: alchUnpatch uses psfgen style selections!
+    alchUnpatch "$segid:$resid" $l0atoms $l1atoms
     return 0
 }
 
@@ -404,31 +440,33 @@ proc ::cphSystem::dealchemifyResidue {segresid} {
 #
 proc ::cphSystem::initializeSystem {pH temperature buildH stateInfo} {
     variable ::cphSystem::resDict
-    dict for {segresid resData} $resDict {
+    dict for {segresidname resData} $resDict {
         # Assign inherent pKa values.
-        if {[dict exists $stateInfo $segresid pKai]} {
-            cphSystem set pKai $segresid [dict get $stateInfo $segresid pKai]
+        if {[dict exists $stateInfo $segresidname pKai]} {
+            cphSystem set pKai $segresidname\
+                    [dict get $stateInfo $segresidname pKai]
         } else { ;# Default to reference pKa.
-            cphSystem set pKai $segresid [cphSystem get pKa $segresid]
+            cphSystem set pKai $segresidname [cphSystem get pKa $segresidname]
         }
 
         # Assign states.
-        if {[dict exists $stateInfo $segresid state]} {
-            set state [dict get $stateInfo $segresid state]
-            cphSystem initializeState $segresid $state 
+        if {[dict exists $stateInfo $segresidname state]} {
+            set state [dict get $stateInfo $segresidname state]
+            cphSystem initializeState $segresidname $state 
         } else { ;# default randomization based on pKai and pH
-            cphSystem initializeState random $segresid $pH
+            cphSystem initializeState random $segresidname $pH
         }
     }
     # Final pass - Apply the patches
-    foreach segresid [cphSystem get segresids] {
-        patch [cphSystem get statePatch $segresid] $segresid
+    foreach segresidname [cphSystem get segresidnames] {
+        lassign [split $segresidname ":"] segid resid
+        patch [cphSystem get statePatch $segresidname] "$segid:$resid"
     }
     guesscoord
-    foreach segresid [cphSystem get segresids] {
-        cphSystem alchemifypsf $segresid 0.0 $temperature $buildH
-        cphSystem dealchemifypsf $segresid
-        cphSystem update 1 $segresid
+    foreach segresidname [cphSystem get segresidnames] {
+        cphSystem alchemifypsf $segresidname 0.0 $temperature $buildH
+        cphSystem dealchemifypsf $segresidname
+        cphSystem update 1 $segresidname
     }
     regenerate angles dihedrals
     # NB - These changes are only reflected in _memory_ for the cphSystem and
@@ -444,6 +482,35 @@ proc ::cphSystem::initializeSystem {pH temperature buildH stateInfo} {
 proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
     variable ::cphSystem::resDict
     variable ::cphSystem::resDefDict [checkResDefs $resDefs]
+
+    # Check for special sub-type residues.
+    # In general, THESE WILL NOT BE MATCHED BY NAME, but rather are sub-typed
+    # within one or more other residue definitions. Therefore we keep a
+    # separate dict of residue names that have a sub-type (these may not even
+    # be titratable otherwise, such as terminal amino acids). The keys are
+    # the actual residue names and the values are a list of the possible
+    # sub-types. However, even if a sub-typed residue is present, that doesn't
+    # mean the sub-residue exists. There are two (optional) additional checks
+    # in the sub-residue definition:
+    #
+    # subatoms - a list of atom names that MUST exist
+    # notsubatoms - a list of atom names that MUST NOT exist
+    #
+    # For example, alanine is not titratable, but may contain a titratable
+    # C-terminus. The terminal oxygen atoms must exist, but are only titratable
+    # if they are not blocked by amination or methylation.
+    #
+    set subresDefDict [dict create]
+    dict for {subresname resDef} $resDefDict {
+        if {![dict exists $resDef subtypeof]} continue
+        foreach subtypedRes [dict get $resDef subtypeof] {
+            if {![dict exists $subresDefDict $subtypedRes]} {
+                dict set subresDefDict $subtypedRes [list]
+            }
+            dict lappend subresDefDict $subtypedRes $subresname
+        }
+    }
+
     # Read in whatever files were specified to NAMD.
     set Args [list [structure] pdb [coordinates]]
     if {[isset binCoordinates]} {
@@ -454,27 +521,81 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
     }
     resetpsf
     readpsf {*}$Args
+
+    set definedResidues [dict keys $resDefDict]
+    set definedSubResidues [dict keys $subresDefDict]
     foreach segid [segment segids] {
         foreach resid [segment resids $segid] {
             set resname [segment residue $segid $resid]
             set segresid [format "%s:%s" $segid $resid]
+            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} {
                     continue
                 }
-                print "aliasing $segid:$resid:$resname to $realName"
+                print "aliasing $segresidname to $realName"
                 psfset resname $segid $resid $realName
                 set resname $realName
             }
-            # Bail here if the residue name does not match any of the
-            # definitions or the segresid is explicitly excluded.
-            if {[lsearch -nocase [dict keys $resDefDict] $resname] < 0
-                || [lsearch -nocase $segresExcls $segresid] >= 0} {
+            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
+            }
+            # Break here if nothing to do.
+            if {(!$resIsDefined && !$resIsSubtyped) || $resIsExcluded} {
                 continue
             }
-            dict set resDict $segresid [dict create]
-            dict set resDict $segresid resname $resname
+
+            if {$resIsDefined} { 
+                dict set resDict $segresidname [dict create]
+                dict set resDict $segresidname resname $resname
+            }
+
+            if {!$resIsSubtyped} continue
+
+            # The conditions for subtyping may still not be met...
+            set atoms [segment atoms $segid $resid]
+            foreach subresname [dict get $subresDefDict $resname] {
+                set subatoms [list]
+                if {[dict exists $resDefDict $subresname subatoms]} {
+                    set subatoms [dict get $resDefDict $subresname subatoms]
+                }
+                set notsubatoms [list]
+                if {[dict exists $resDefDict $subresname notsubatoms]} {
+                    set notsubatoms\
+                            [dict get $resDefDict $subresname notsubatoms]
+                }
+
+                set allSubatomsExist 1
+                foreach atom $subatoms {
+                    if {[lsearch -nocase $atoms $atom] < 0} {
+                        set allSubatomsExist 0
+                        break
+                    }
+                }
+                set allNotSubatomsDoNotExist 1
+                foreach atom $notsubatoms { 
+                    if {[lsearch -nocase $atoms $atom] >= 0} {
+                        set allNotSubatomsDoNotExist 0
+                        break
+                    }
+                }
+                if {$allSubatomsExist && $allNotSubatomsDoNotExist} {
+                    set segresidname [format "%s:%s" $segresid $subresname]
+                    dict set resDict $segresidname [dict create]
+                    dict set resDict $segresidname resname $subresname
+                }
+            }
+
         }
     }
     return $resDict
@@ -509,7 +630,6 @@ proc ::cphSystem::checkResDefs {resDefs} {
             set l1atoms [dict get $resDefs $resname l1atoms]
         }
 
-
         if {[llength $l0atoms] != [llength $l1atoms]} { 
             abort "Mismatch in atom definitions for residue $resname"
         }
@@ -557,7 +677,7 @@ proc ::cphSystem::checkResDefs {resDefs} {
 #
 # Arguments:
 # ----------
-# segresid : string
+# segresidname : string
 #   residue specification as "<segid>:<resid>" - this is the same syntax as for
 #   the regular psfgen patch command.
 # state : string
@@ -567,12 +687,11 @@ proc ::cphSystem::checkResDefs {resDefs} {
 # --------
 # None
 #
-proc ::cphSystem::assignState {segresid state} {
-    cphSystem set state $segresid $state
-    cphSystem propose titration $segresid
-    set trialState [cphSystem get state $segresid]
-    cphSystem set state $segresid [cphSystem get trialState $segresid]
-    cphSystem set trialState $segresid $trialState
+proc ::cphSystem::assignState {segresidname state} {
+    cphSystem set state $segresidname $state
+    cphSystem propose titration $segresidname
+    cphSystem set state $segresidname [cphSystem get trialState $segresidname]
+    cphSystem set trialState $segresidname $state
     return 0
 }
 
@@ -590,7 +709,7 @@ proc ::cphSystem::assignState {segresid state} {
 #
 # Arguments:
 # ----------
-# segresid : string
+# segresidname : string
 #   residue specification as "<segid>:<resid>" - this is the same syntax as for
 #   the regular psfgen patch command.
 # pH : float
@@ -600,12 +719,12 @@ proc ::cphSystem::assignState {segresid state} {
 # --------
 # None 
 #
-proc ::cphSystem::randomizeState {segresid pH} {
+proc ::cphSystem::randomizeState {segresidname pH} {
     while {true} {
-        set states [cphSystem get stateList $segresid]
-        cphSystem set state $segresid [choice $states]
-        cphSystem propose titration $segresid
-        set du [cphSystem compute inherent $pH $segresid]
+        set states [cphSystem get stateList $segresidname]
+        cphSystem set state $segresidname [choice $states]
+        cphSystem propose titration $segresidname
+        set du [cphSystem compute inherent $pH $segresidname]
         if {[metropolisAcceptance $du]} { 
             return 0
         }
@@ -620,16 +739,16 @@ proc ::cphSystem::randomizeState {segresid pH} {
 #
 # Getters for system and residue attributes, called as:
 #
-#   <attribute> [<segresid> [<args>]]
+#   <attribute> [<segresidname> [<args>]]
 #
-# <attribute> is the name of either a system attribute (segresid selections are
-# invalid) or else a residue attribute. For SOME residue attributes, not
-# specifying a segresid will return a list for all residues. Some attributes
-# also require some additional arguments (see below)
+# <attribute> is the name of either a system attribute (segresidname selections
+# are invalid) or else a residue attribute. For SOME residue attributes, not
+# specifying a segresidname will return a list for all residues. Some
+# attributes also require some additional arguments (see below)
 #
 # system attributes  description
 # -----------------  -----------
-# segresids          list of all segresids
+# segresidnames      list of all segresidnames
 # numresidues        number of residues in the system
 # resdefs            list of defined resnames
 # numdefs            number of defined resnames
@@ -645,7 +764,6 @@ proc ::cphSystem::randomizeState {segresid pH} {
 # Tref               reference temperature for pKa
 # occupancy          occupancy vector for the current state
 # trialOccupancy     occupancy vector for the trial state
-# reslabel           residue name with prepended segresid
 # stateList          all possible states
 # trialStateList     all possible trial states (not the current state)
 # dGPair             pair dG for current/trial states
@@ -656,20 +774,20 @@ proc ::cphSystem::randomizeState {segresid pH} {
 # alchAtomLists*     lists of alchemical atoms at 0/1
 # alchBonds*^        extraBonds entries
 #
-# * - segresid selection is required
+# * - segresidname selection is required
 # ^ - requires additional arguments
 #
-proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
+proc ::cphSystem::cphSystemGet {attr {segresidname {}} args} {
     variable ::cphSystem::resDict
     variable ::cphSystem::resDefDict
 
-    set getAll [expr {![llength $segresid]}]
-    if {!$getAll && ![dict exists $resDict $segresid]} {
-        abort "cphSystemGet: Invalid segresid $segresid"
+    set getAll [expr {![llength $segresidname]}]
+    if {!$getAll && ![dict exists $resDict $segresidname]} {
+        abort "cphSystemGet: Invalid segresidname $segresidname"
     }
     # System attributes - any selection is invalid.
     #
-    if {[string match -nocase $attr segresids]} {
+    if {[string match -nocase $attr segresidnames]} {
         return [dict keys $resDict]
     } elseif {[string match -nocase $attr numresidues]} {
         return [dict size $resDict]
@@ -684,42 +802,35 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
         if {$getAll} {
             return [getAllResAttr $attr]
         } else {
-            return [getResAttr $attr $segresid]
+            return [getResAttr $attr $segresidname]
         }
     } 
     if {[lsearch -nocase {dG pKa Tref} $attr] > -1} {
         if {$getAll} {
             return [getAllResDefAttr $attr]
         } else {
-            return [getResDefAttr $attr $segresid]
+            return [getResDefAttr $attr $segresidname]
         }
     } 
     if {[string match -nocase $attr occupancy]} {
         if {$getAll} {
             return [getAllOccupancy]
         } else {
-            return [getOccupancy $segresid]
+            return [getOccupancy $segresidname]
         }
     } 
     if {[string match -nocase $attr trialOccupancy]} {
         if {$getAll} {
             cannotGetAll $attr
         } else {
-            return [getTrialOccupancy $segresid]
+            return [getTrialOccupancy $segresidname]
         }
     } 
-    if {[string match -nocase $attr reslabel]} {
-        if {$getAll} {
-            return [getAllReslabel]
-        } else {
-            return [getReslabel $segresid]
-        }
-    }
     if {[string match -nocase $attr stateList]} {
         if {$getAll} {
             cannotGetAll $attr
         } else {
-            set resname [getResAttr resname $segresid]
+            set resname [getResAttr resname $segresidname]
             return [dict keys [dict get $resDefDict $resname states]] 
         }
     }
@@ -727,8 +838,8 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
         if {$getAll} {
             cannotGetAll $attr
         } else {
-            set resname [getResAttr resname $segresid]
-            set state [getResAttr state $segresid]
+            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]
         }
@@ -738,9 +849,9 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
             cannotGetAll $attr
         } else {
             if {[lsearch -nocase {dGPair pKaPair} $attr] > -1} {
-                return [getResDefAttr $attr $segresid]
+                return [getResDefAttr $attr $segresidname]
             } else {
-                return [getResAttr $attr $segresid]
+                return [getResAttr $attr $segresidname]
             }
         }
     }
@@ -748,8 +859,8 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
         if {$getAll} {
             cannotGetAll $attr
         } else {
-            set resname [getResAttr resname $segresid]
-            set state [getResAttr state $segresid]
+            set resname [getResAttr resname $segresidname]
+            set state [getResAttr state $segresidname]
             return [format "%s%s" $resname $state]
         }
     }
@@ -757,8 +868,8 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
         if {$getAll} {
             cannotGetAll $attr
         } else {
-            set resname [getResAttr resname $segresid]
-            set trialState [getResAttr trialState $segresid]
+            set resname [getResAttr resname $segresidname]
+            set trialState [getResAttr trialState $segresidname]
             return [format "%sH%s" $resname $trialState]
         }
     }
@@ -766,8 +877,8 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
         if {$getAll} {
             cannotGetAll $attr
         } else {
-            return [list [getResDefAttr l0atoms $segresid]\
-                         [getResDefAttr l1atoms $segresid]]
+            return [list [getResDefAttr l0atoms $segresidname]\
+                         [getResDefAttr l1atoms $segresidname]]
         }
     }
     if {[string match -nocase $attr alchBonds]} {
@@ -776,9 +887,9 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
         } else {
             lassign $args k
             set bondEntries [list]
-            foreach l0atom [getResDefAttr l0atoms $segresid]\
-                    l1atom [getResDefAttr l1atoms $segresid] {
-                lassign [split $segresid ":"] segid resid
+            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}]
@@ -791,15 +902,15 @@ proc ::cphSystem::cphSystemGet {attr {segresid {}} args} {
 }
 
 proc ::cphSystem::cannotGetAll {attr} {
-    abort "cphSystemGet: Cannot get all $attr - must select a segresid"
+    abort "cphSystemGet: Cannot get all $attr - must select a segresidname"
 }
 
 # ------------------------------
 # Getters for residue attributes
 # ------------------------------
-proc ::cphSystem::getResAttr {attr segresid} {
+proc ::cphSystem::getResAttr {attr segresidname} {
     variable ::cphSystem::resDict
-    return [dict get $resDict $segresid $attr]
+    return [dict get $resDict $segresidname $attr]
 }
 
 proc ::cphSystem::getAllResAttr {attr} {
@@ -814,10 +925,10 @@ proc ::cphSystem::getAllResAttr {attr} {
 # -----------------------------------------
 # Getters for residue definition attributes
 # -----------------------------------------
-proc ::cphSystem::getResDefAttr {attr segresid} {
+proc ::cphSystem::getResDefAttr {attr segresidname} {
     variable ::cphSystem::resDict
     variable ::cphSystem::resDefDict
-    set resname [dict get $resDict $segresid resname]
+    set resname [dict get $resDict $segresidname resname]
     return [dict get $resDefDict $resname $attr]
 }
 
@@ -841,17 +952,17 @@ proc ::cphSystem::state2Occupancy {resname state} {
 
 # NB: When returning all occupancies, only the current state can be used and
 #   the list is flattened.
-proc ::cphSystem::getOccupancy {segresid} {
+proc ::cphSystem::getOccupancy {segresidname} {
     variable ::cphSystem::resDict
-    set resname [dict get $resDict $segresid resname]
-    set state [dict get $resDict $segresid state]
+    set resname [dict get $resDict $segresidname resname]
+    set state [dict get $resDict $segresidname state]
     return [state2Occupancy $resname $state]
 }
 
-proc ::cphSystem::getTrialOccupancy {segresid} {
+proc ::cphSystem::getTrialOccupancy {segresidname} {
     variable ::cphSystem::resDict
-    set resname [dict get $resDict $segresid resname]
-    set state [dict get $resDict $segresid trialState]
+    set resname [dict get $resDict $segresidname resname]
+    set state [dict get $resDict $segresidname trialState]
     return [state2Occupancy $resname $state]
 }
 
@@ -867,20 +978,6 @@ proc ::cphSystem::getAllOccupancy {} {
     return $retList
 }
 
-proc ::cphSystem::getReslabel {segresid} {
-    variable ::cphSystem::resDict
-    return "$segresid:[dict get $resDict $segresid resname]"
-}
-
-proc ::cphSystem::getAllReslabel {} {
-    variable ::cphSystem::resDict
-    set retList [list]
-    dict for {segresid resData} $resDict {
-        lappend retList "$segresid:[dict get $resData resname]"
-    }
-    return $retList
-}
-
 # ----------------------
 # Other helper functions
 # ----------------------
@@ -932,7 +1029,7 @@ proc ::cphSystem::index2flatindex {i j} {
 #
 # Setters for residue attributes, called as:
 #
-#   <attribute> <segresid> <value>
+#   <attribute> <segresidname> <value>
 #
 # <attribute> is the name of a residue attribute.
 #
@@ -942,30 +1039,30 @@ proc ::cphSystem::index2flatindex {i j} {
 # trialState         proposed trial state
 # pKai               minimal pKai list for this residue
 #
-proc ::cphSystem::cphSystemSet {attr segresid value} {
+proc ::cphSystem::cphSystemSet {attr segresidname value} {
     variable ::cphSystem::resDict
     variable ::cphSystem::resDefDict
-    if {![dict exists $resDict $segresid]} {
-        abort "cphSystemSet: Invalid segresid $segresid"
+    if {![dict exists $resDict $segresidname]} {
+        abort "cphSystemSet: Invalid segresidname $segresidname"
     }
     if {[lsearch -nocase {state trialState} $attr] > -1} {
-        set states [cphSystem get stateList $segresid]
+        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 $segresid"
+                abort "Invalid state assignment $value for residue $segresidname"
             }
         }
-        dict set resDict $segresid $attr $value
+        dict set resDict $segresidname $attr $value
         return $value
     }
     if {[string match $attr pKai]} {
-        set resname [cphSystem get resname $segresid]
+        set resname [cphSystem get resname $segresidname]
         set resDef [dict get $resDefDict $resname]
         set pKaiMatrix [resDef2Matrix $resDef $value]
-        dict set resDict $segresid pKai $value
-        dict set resDict $segresid pKaiPair $pKaiMatrix 
+        dict set resDict $segresidname pKai $value
+        dict set resDict $segresidname pKaiPair $pKaiMatrix 
         return $pKaiMatrix
     }
 }
index be7bcd5..a191c76 100644 (file)
@@ -52,16 +52,25 @@ proc ::namdcph::cphRun {numsteps {numcycles 1}} {
     if {$::namdcph::numMinSteps > 0} {
         set tmp [firstTimeStep]
         minimize $::namdcph::numMinSteps
+        if {[isset temperature]} {
+          reinitvels [temperature]
+        } else {
+          reinitvels [$::thermostatTempCmd]
+        }
         firstTimeStep $tmp
     }
     set cphlog [openCpHLog]
     set firstCycle 1
     set lastCycle [expr {$firstCycle + $numcycles - 1}] 
+    set runArgs [list $numsteps]
     for {set cycle $firstCycle} {$cycle <= $lastCycle} {incr cycle} { 
-        # (1) Propose a move from the full move set.
+        # (1) Perform whatever equilibrium sampling is desired.
         #
-        lassign [proposeMove] accept swNumsteps segresidList nattempts
-        # (2) At this point we have either selected a switch or rejected a
+        runMD {*}$runArgs
+        # (2) Propose a move from the full move set.
+        #
+        lassign [proposeMove] 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.
         #
         if {!$accept} {
@@ -69,15 +78,12 @@ proc ::namdcph::cphRun {numsteps {numcycles 1}} {
             set runArgs [list norepeat $numsteps]
         } else {
             cphPrint "Proposal accepted ($nattempts), attemping a switch."
-            set accept [runSwitch $swNumsteps $segresidList]
+            set accept [runSwitch $swNumsteps $segresidnameList]
             set runArgs [list $numsteps]
             # Only accumulate statistics for attempted switches.
-            accumulateAcceptanceRate $accept $segresidList
+            accumulateAcceptanceRate $accept $segresidnameList
         }
-        cphSystem update $accept $segresidList
-        # (3) Perform whatever equilibrium sampling is desired.
-        #
-        runMD {*}$runArgs
+        cphSystem update $accept $segresidnameList
         # (4) Output cycle information and a checkpoint file if necessary.
         #
         puts $cphlog "[format "%6d" $cycle] [cphSystem get occupancy]"
@@ -96,11 +102,15 @@ proc ::namdcph::cphRun {numsteps {numcycles 1}} {
 #
 # FOR ADVANCED USE ONLY!!!
 #
-# ::namdcph::testResidue
+# ::namdcph::checkResidueDefinitions
 #
-# Test one or more residue definitions.
+# Check that appropriate residue definitions are defined for this residue.
 #
-proc ::namdcph::testResidue {resnames {verbose 0}} {
+# This essentially just scans all possible states for each titratable residue,
+# of a given type in the system. THIS OFFERS NO VALIDATION AS TO WHETHER OR NOT
+# THE DEFINITIONS MAKE SENSE, it simply checks that the exist.
+#
+proc ::namdcph::checkResidueDefinitions {resnames} {
     cphWarnExpt
     pH 7.0
     cphForceConstant 0.0
@@ -110,15 +120,12 @@ proc ::namdcph::testResidue {resnames {verbose 0}} {
     $::thermostatTempCmd 0.0
     outputEnergies 1
     alchLambdaFreq 0
-    set report [dict create]
-    set terms [list BOND ANGLE DIHED IMPRP ELECT VDW POTENTIAL]
-    foreach resLabel [cphSystem get reslabel] {
-        lassign [split $resLabel ":"] segid resid residue
-        if {[lsearch $resnames $residue] < 0} {
+    foreach segresidname [cphSystem get segresidnames] {
+        lassign [split $segresidname ":"] segid resid resname
+        if {[lsearch $resnames $resname] < 0} {
             continue
         } 
-        set segresid "$segid:$resid"
-        set states [cphSystem get stateList $segresid]
+        set states [cphSystem get stateList $segresidname]
         foreach state0 $states {
             foreach state1 $states {
                 if {[string match $state0 $state1]} {
@@ -127,83 +134,78 @@ proc ::namdcph::testResidue {resnames {verbose 0}} {
                 # Force the initial state as state0
                 output [getMDBasename]
                 psfgenRead [getMDBasename]
-                cphSystem initializeState $segresid $state0
-                cphSystem alchemifypsf $segresid 0.0 0.0
-                cphSystem dealchemifypsf $segresid
-                regenerate angles dihedrals 
-                cphSystem update 1 $segresid
-                psfgenWrite [getMDBasename] 
-                reloadAndReinit [getMDBasename] false
-                # Now perform an alchemifying and dealchemifying transition.
-                run 0
-                storeEnergies
-                array set MDEnergies0 [array get ::energyArray]
-                cphSystem set trialState $segresid $state1
-                alchemify $segresid
-                reloadAndReinit [getSWBasename] true
-                alchLambda 0.0
-                run 0
-                storeEnergies
-                array set SWEnergies0 [array get ::energyArray]
-                alchLambda 1.0
-                run 0
-                storeEnergies
-                array set SWEnergies1 [array get ::energyArray]
-                dealchemify $segresid
-                cphSystem update 1 $segresid
-                run 0
-                storeEnergies
-                array set MDEnergies1 [array get ::energyArray]
-                # Record the differences in energy.
-                set label [format "%s%sH%s" $residue $state0 $state1]
-                dict set report $label [dict create]
-                dict set report $label 0 [dict create]
-                dict set report $label 1 [dict create]
-                foreach term $terms {
-                    set diff0 [expr {$MDEnergies0($term)-$SWEnergies0($term)}]
-                    set diff1 [expr {$MDEnergies1($term)-$SWEnergies1($term)}]
-                    if {[expr {abs($diff0)}] < 1e-5} {
-                        set diff0 0.0
-                    }
-                    if {[expr {abs($diff1)}] < 1e-5} {
-                        set diff0 0.0
-                    }
-                    dict set report $label 0 $term $diff0
-                    dict set report $label 1 $term $diff1
-                }
-            }
-        }
-    }
-    dict for {label data} $report {
-        cphPrint [format "%-14s  %14s" $label "RMS-error"]
-        foreach lambda {0 1} label {alchemify dealchemify} {
-            set i 1
-            set termStr " "
-            set diffStr " "
-            set rmse 0.0
-            dict for {term diff} [dict get $data 0] {
-                set termStr [format "%s % 14s" $termStr $term]
-                set diffStr [format "%s % 14.4f" $diffStr $diff]
-                set rmse [expr {$rmse + $diff*$diff}]
-                incr i
-                if {$i == 5} {
-                    set termStr "$termStr     "
-                    set diffStr "$diffStr     "
-                    set i 0
-                }
-            }
-            set rmse [expr {2*$rmse / [dict size $data]}]
-            cphPrint [format "%-14s: % 14.4f" $label $rmse]
-            if {$verbose} {
-                cphPrint $termStr
-                cphPrint $diffStr
+                cphSystem initializeState $segresidname $state0
+                alchemify $segresidname
+                dealchemify $segresidname
+                cphSystem update 1 $segresidname
+                # Run an instantaneous switch
+                runTestSwitch $segresidname $state1
             }
         }
     }
+    return 
+}
+
+#
+# FOR ADVANCED USE ONLY!!!
+#
+# ::namdcph::testResidueSwitch
+#
+# Run an instantenous switch for the given residue between the given states.
+#
+# This is meant as an interface function for checking energies against
+# non-constant pH calculations, but cannot perform that comparison on its own.
+#
+proc ::namdcph::testResidueSwitch {segresidname state0 state1} {
+    cphWarnExpt
+    pH 7.0
+    cphForceConstant 0.0
+    # Initialize NAMD and build a constant pH enabled PSF.
+    cphSetResidueState $segresidname $state0
+    initialize
+    finalize
+    $::thermostatTempCmd 0.0
+    outputEnergies 1
+    alchLambdaFreq 0
+    return [runTestSwitch $segresidname $state1]
+}
+
+#
+# FOR ADVANCED USE ONLY!!!
+#
+# ::namdcph::runTestSwitch
+#
+# This is just an internal convenience function for running instantaneous test
+# switches.
+#
+# See checkResidueDefinitions and testResidueSwitch
+#
+proc ::namdcph::runTestSwitch {segresidname state1} {
+    # You can't make lists of arrays or arrays of arrays, so the return type
+    # has to be a dict of arrays (callback only returns arrays).
+    set retVals [dict create]
+    run 0
+    storeEnergies
+    dict set retVals MDEnergies0 [array get ::energyArray]
+    cphSystem set trialState $segresidname $state1
+    alchemify $segresidname
+    alchLambda 0.0
+    run 0
+    storeEnergies
+    dict set retVals SWEnergies0 [array get ::energyArray]
+    alchLambda 1.0
+    run 0
+    storeEnergies
+    dict set retVals SWEnergies1 [array get ::energyArray]
+    dealchemify $segresidname
+    cphSystem update 1 $segresidname
+    run 0
+    storeEnergies
+    dict set retVals MDEnergies1 [array get ::energyArray]
     # Cleanup temporary files
     file delete {*}[glob [getSWBasename].*]
     file delete {*}[glob [getMDBasename].*]
-    return 
+    return $retVals
 }
 
 #
@@ -213,48 +215,56 @@ proc ::namdcph::testResidue {resnames {verbose 0}} {
 #
 # Test one or more residue definitions.
 #
-proc ::namdcph::cphAnalyzeForce {dcdfilename segresid stateList} {
+proc ::namdcph::cphAnalyzeForce {dcdfilename segresidname state0 state1} {
     cphWarnExpt
     pH 7.0
-    lassign $stateList state0 state1
 
-    cphSetResidueState $segresid $state0
+    cphSetResidueState $segresidname $state0
     initialize
-    set nframes 0
-    coorfile open dcd $dcdfilename
-    while {![coorfile read]} {
-        output [format "%s.%d" [getMDBasename] $nframes]
-        incr nframes
-    }
-    coorfile close
     set initialPSF [structure]
     set initialPDB [coordinates]
     finalize
 
-    for {set n 0} {$n < $nframes} {incr n} {
+    set numsteps 500
+    outputEnergies [expr {$numsteps == 0 ? 1 : int($numsteps)}]
+    alchLambdaFreq [expr {$numsteps == 0 ? 0 : 1}]
+
+    set nframes -1
+    coorfile open dcd $dcdfilename
+    while {![coorfile read]} {
+        incr nframes
+
         # We have to do this so that inputs can be correctly loaded...
-        set basename [format "%s.%d" [getMDBasename] $n]
+        set basename [format "%s.%d" [getMDBasename] $nframes]
+        output $basename
         file copy -force $initialPSF "$basename.psf"
         file copy -force $initialPDB "$basename.pdb"
         reloadAndReinit $basename false
         # Assign the correct state and build protons or dummy atoms
-        cphSystem initializeState $segresid $state0
-        cphSystem alchemifypsf $segresid 0.0 0.0
-        cphSystem dealchemifypsf $segresid
-        regenerate angles dihedrals
-        cphSystem update 1 $segresid
-        psfgenWrite [getMDBasename]
-        reloadAndReinit [getMDBasename] false
+        cphSystem initializeState $segresidname $state0
+        alchemify $segresidname
+        dealchemify $segresidname
+        cphSystem update 1 $segresidname
         # Now build the alchemical atoms
-        cphSystem set trialState $segresid $state1
-        alchemify $segresid
-        reloadAndReinit [getSWBasename] true
-        firsttimestep $n
-        run 0
-        dealchemify $segresid 
+        cphSystem set trialState $segresidname $state1
+        alchemify $segresidname
+        firsttimestep 0
+        run $numsteps
+        if {$numsteps} {
+            storeEnergies
+            set DeltaE [cphSystem compute switch $segresidname]
+            set Work [expr {$::energyArray(CUMALCHWORK) + $DeltaE}]
+            set ReducedWork [expr {$Work / [::kBT]}]
+            set tmp [printProposalSummary $segresidname]
+            cphPrint [format "%s WorkCorr % 10.4f CorrWork % 10.4f"\
+                    [join $tmp "/"] $DeltaE $Work]
+        }
+        dealchemify $segresidname
+        file delete {*}[glob $basename.*]
     }
     file delete {*}[glob [getMDBasename].*]
     file delete {*}[glob [getSWBasename].*]
+    coorfile close
     return
 }
 
@@ -321,8 +331,8 @@ proc ::namdcph::cphSetResidueState {args} {
     checkArglistIsMultiple $args 2
     variable ::namdcph::stateInfo
     for {set i 0} {$i < [llength $args]} {incr i 2} {
-        lassign [lrange $args $i [expr {$i+1}]] segresid state
-        dict set stateInfo $segresid state $state
+        lassign [lrange $args $i [expr {$i+1}]] segresidname state
+        dict set stateInfo $segresidname state $state
     }
     return
 }
@@ -336,11 +346,11 @@ proc ::namdcph::cphSetResiduepKai {args} {
     variable ::namdcph::stateInfo
     for {set i 0} {$i < [llength $args]} {incr i 2} {
         #NB pKai may be a list of values - check individually.
-        lassign [lrange $args $i [expr {$i+1}]] segresid pKai
+        lassign [lrange $args $i [expr {$i+1}]] segresidname pKai
         foreach pKa $pKai {
             checkIsPositive pKai $pKa
         }
-        dict set stateInfo $segresid pKai $pKai
+        dict set stateInfo $segresidname pKai $pKai
     }
     return
 }
@@ -427,8 +437,8 @@ proc ::namdcph::cphNumMinSteps {numsteps} {
 
 proc ::namdcph::cphExcludeResidue {args} {
     variable ::namdcph::excludeList
-    foreach segresid $args {
-         lappend excludeList $segresid
+    foreach segresidname $args {
+         lappend excludeList $segresidname
     }
     return
 }
@@ -488,7 +498,7 @@ proc ::namdcph::runMD {args} {
 # ----------
 # numsteps : int
 #   Number of steps in the switch
-# segresidList : list 
+# segresidnameList : list 
 #   One or more "<segid>:<resid>" specifications - this is the same syntax as 
 #   for the regular psfgen patch command.
 #
@@ -497,7 +507,7 @@ proc ::namdcph::runMD {args} {
 # accept : boolean
 #   Result of MC accept/reject test 
 #
-proc ::namdcph::runSwitch {numsteps segresidList} {
+proc ::namdcph::runSwitch {numsteps segresidnameList} {
     # (1) Checkpoint and modify output parameters. 
     #
     storeEnergies
@@ -512,7 +522,7 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
     }
     # (2) Build the alchemical switch inputs.
     #
-    alchemify $segresidList
+    alchemify $segresidnameList
     # (3) Run the switch trajectory with momentum reversal.
     #
     runprpswitch $numsteps
@@ -521,11 +531,11 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
     # (4) Compute the work with state dependent energy shifts.
     #
     storeEnergies
-    set DeltaE [cphSystem compute switch $segresidList]
+    set DeltaE [cphSystem compute switch $segresidnameList]
     set Work [expr {$::energyArray(CUMALCHWORK) + $DeltaE}]
-    set ReducedWork [expr {$Work / ($::BOLTZMANN*[$::thermostatTempCmd])}]
+    set ReducedWork [expr {$Work / [::kBT]}] 
     set accept [metropolisAcceptance $ReducedWork]
-    set tmp [printProposalSummary $segresidList]
+    set tmp [printProposalSummary $segresidnameList]
     cphPrint [format "%s WorkCorr % 10.4f CorrWork % 10.4f"\
             [join $tmp "/"] $DeltaE $Work]
     # (5) Reinitialize for the next MD step based on accept/reject.
@@ -539,7 +549,7 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
 
     if {$accept} {
         cphPrint "Switch accepted!"
-        dealchemify $segresidList
+        dealchemify $segresidnameList
     } else {
         cphPrint "Switch rejected!"
         alch off
@@ -561,7 +571,7 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
 #
 # Arguments:
 # ----------
-# segresidList : list of strings
+# segresidnameList : list of strings
 #   One or more "<segid>:<resid>" specifications - this is the same syntax as 
 #   for the regular psfgen patch command.
 #
@@ -569,7 +579,7 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
 # --------
 # None
 #
-proc ::namdcph::alchemify {segresidList} {  
+proc ::namdcph::alchemify {segresidnameList} {  
     set T [$::thermostatTempCmd]
     set FrcCons [getAlchFrcCons]
     alch on
@@ -579,8 +589,8 @@ proc ::namdcph::alchemify {segresidList} {
     psfgenRead [getMDBasename]
     # (2) Apply patches and build coordinates and velocities.
     #
-    foreach segresid $segresidList {
-        cphSystem alchemifypsf $segresid $FrcCons $T
+    foreach segresidname $segresidnameList {
+        cphSystem alchemifypsf $segresidname $FrcCons $T
     }
     regenerate angles dihedrals
     # (3) Write a new set of inputs and reinitialize.
@@ -591,8 +601,8 @@ proc ::namdcph::alchemify {segresidList} {
     # the atomid queries would all return zero (that would be bad).
     #
     set ExtraBondsFile [open [swFilename extrabonds] "w"]
-    foreach segresid $segresidList {
-        puts $ExtraBondsFile [cphSystem get alchBonds $segresid $FrcCons]
+    foreach segresidname $segresidnameList {
+        puts $ExtraBondsFile [cphSystem get alchBonds $segresidname $FrcCons]
     }
     close $ExtraBondsFile
     reloadAndReinit [getSWBasename] true
@@ -605,7 +615,7 @@ proc ::namdcph::alchemify {segresidList} {
 #
 # Arguments:
 # ----------
-# segresidList : list of strings
+# segresidnameList : list of strings
 #   One or more "<segid>:<resid>" specifications - this is the same syntax as 
 #   for the regular psfgen patch command.
 #
@@ -613,11 +623,11 @@ proc ::namdcph::alchemify {segresidList} {
 # --------
 # None
 #
-proc ::namdcph::dealchemify {segresidList} {
+proc ::namdcph::dealchemify {segresidnameList} {
     output [getSWBasename]
     psfgenRead [getSWBasename]
-    foreach segresid $segresidList {
-        cphSystem dealchemifypsf $segresid
+    foreach segresidname $segresidnameList {
+        cphSystem dealchemifypsf $segresidname
     }
     psfgenWrite [getMDBasename] [swFilename xsc]
     alch off
@@ -771,7 +781,7 @@ proc ::namdcph::openCpHLog {} {
     namdFileBackup $logFilename
     set cphlog [open $logFilename "w"]
     puts $cphlog "#pH [getpH]"
-    puts $cphlog "#[join [cphSystem get reslabel] " "]"
+    puts $cphlog "#[join [cphSystem get segresidnames] " "]"
     return $cphlog 
 }
 
@@ -892,16 +902,17 @@ proc ::namdcph::initialize {} {
     }
     # Use the system topology to make assignments from the restart
     if {[llength $states] && [llength $pKais]} { 
-        foreach segresid [cphSystem get segresids] state $states pKai $pKais {
-            if {![dict exists $stateInfo $segresid]} {
-                dict set stateInfo $segresid state $state
-                dict set stateinfo $segresid pKai $pKai
+        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 $segresid state]} {
-                    dict set stateInfo $segresid state $state
+                if {![dict exists $stateInfo $segresidname state]} {
+                    dict set stateInfo $segresidname state $state
                 }
-                if {![dict exists $stateInfo $segresid pKai]} {
-                    dict set stateInfo $segresid pKai $pKai
+                if {![dict exists $stateInfo $segresidname pKai]} {
+                    dict set stateInfo $segresidname pKai $pKai
                 }
             }
         }
@@ -983,11 +994,7 @@ proc ::namdcph::initializeAlch {} {
     extraBonds on
     extraBondsFile [swFilename extrabonds]
     clearExtraBonds
-    run 0
-    # NB: Because alchLambda is incremented _before_ dynamics steps, if 
-    # alchLambdaFreq were nonzero when "run 0" is called we would get a 
-    # confusing energy log with a fractional alchLambda.
-    #
+    startup
     alchLambdaFreq 1
     alch off
     return
@@ -1032,9 +1039,10 @@ proc ::namdcph::printSystemSummary {} {
     cphPrint "[llength [cphSystem get occupancy]] PROTONATION SITE(S)"
     cphPrint $StarBar
     cphPrint [format "%-19s : %5s : %s" segid:resid:resname state pKai]
-    foreach resLabel [cphSystem get reslabel] state [cphSystem get state]\
+    foreach segresidname [cphSystem get segresidnames]\
+            state [cphSystem get state]\
             pKaiList [cphSystem get pKai] {
-        cphPrint [format "%-19s : % 5s : %-s" $resLabel $state $pKaiList]
+        cphPrint [format "%-19s : % 5s : %-s" $segresidname $state $pKaiList]
     }
     cphPrint $StarBar
     return
@@ -1059,13 +1067,12 @@ proc ::namdcph::printTitratorSummary {} {
 }
 
 # ::namdcph::printProposalSummary
-proc ::namdcph::printProposalSummary {segresidList} {
+proc ::namdcph::printProposalSummary {segresidnameList} {
     set retList [list]
-    foreach segresid $segresidList {
-        set reslabel [cphSystem get reslabel $segresid]
-        set state [cphSystem get state $segresid]
-        set trialState [cphSystem get trialState $segresid]
-        lappend retList "$reslabel:$state:$trialState"
+    foreach segresidname $segresidnameList {
+        set state [cphSystem get state $segresidname]
+        set trialState [cphSystem get trialState $segresidname]
+        lappend retList "$segresidname:$state:$trialState"
     }
     return $retList
 }
index e07c147..42f140e 100644 (file)
@@ -199,6 +199,18 @@ proc getThermostat {{forbidBerendsen true}} {
     return 1 
 }
 
+# kBT
+#
+# Return kBT using the current thermostat temperature. If no thermostat is set
+# then return 0.
+#
+proc kBT {} {
+    if {$::thermostatIsSet} {
+        return [expr {$::BOLTZMANN*[$::thermostatTempCmd]}]
+    }
+    return 0.0
+}
+
 # getBarostat
 #
 # This only requires that the barostat itself has been invoked, not that a