Enable proton transfer moves 34/3134/3
authorBrian Radak <brian.radak@gmail.com>
Tue, 17 Oct 2017 18:56:30 +0000 (13:56 -0500)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Sat, 4 Nov 2017 15:10:13 +0000 (10:10 -0500)
This patch adds the (experimental, for now) ability to add proton
transfer neMD/MC moves. These do not exist by default and must be
added with the new (possibly fragile) keyword:

cphProposeProtonTransfer <segid:resid>/<segid:resid> [..[..]]

which takes psfgen style segid:resid combinations _in pairs_ and
separated by a "/". A new move is then added to the pool of MC
moves (the default move weight is always 1.0) which will propose
a proton transfer between the two residues (i.e. one residue will
lose a proton and the other will gain a proton). If no transfer is
possible, then the move is rejected.

Preliminary tests indicate that the switches for these moves need
to be longer than for single residue titrations, but validation
tests are still being run.

Change-Id: I9bbcd564504f17606908b84cf870b531ea77a7fb

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

index 1c19af5..9354cb0 100644 (file)
@@ -19,6 +19,7 @@ namespace eval ::namdcphwrapper {
             cphRestartFile cphRestartFreq\
             cphForceConstant cphMDBasename cphSwitchBasename\
             cphMaxProposalAttempts cphNumMinSteps cphProposalWeight\
+            cphProposeProtonTransfer\
             cphNumstepsPerSwitch\
             testResidue cphAnalyzeForce cphExcludeResidue
 
index 45874a4..9b5978d 100644 (file)
@@ -62,6 +62,25 @@ proc ::cphTitrator::buildTitrator {systempH moveInfo} {
             dict unset moveInfo $segresid
         }
     }
+    # Build proton transfer moves
+    foreach moveLabel [dict get $moveInfo ptransfer] {
+        lassign [split $moveLabel "/"] segresid1 segresid2
+        dict set MoveSet $moveLabel proposalCmd\
+                "proposeProtonTransferMove $moveLabel" 
+        dict set MoveSet $moveLabel segresidList [list $segresid1 $segresid2]
+        dict set MoveSet $moveLabel weight\
+                [dictPopOrDefault moveInfo $moveLabel weight]
+        dict set MoveSet $moveLabel numsteps\
+                [dictPopOrDefault moveInfo $moveLabel numsteps]
+        dict set MoveSet $moveLabel successes 0
+        dict set MoveSet $moveLabel attempts 0
+        if {[dict exists $moveInfo $moveLabel]\
+            && ![dict size [dict get $moveInfo $moveLabel]]} {
+            dict unset moveInfo $moveLabel
+        }
+    }
+    dict unset moveInfo ptransfer
+    ###
     computeWeightSum
     return $moveInfo
 }
@@ -124,6 +143,20 @@ proc ::cphTitrator::proposeResidueMove {segresid} {
    return $accept
 }
 
+proc ::cphTitrator::proposeProtonTransferMove {moveLabel} {
+    variable ::cphTitrator::pH
+    lassign [split $moveLabel "/"] segresid1 segresid2
+    set errCode [::cphSystem::proposeProtonTransfer $segresid1 $segresid2]
+    if {$errCode > 0} {
+        set du1 [cphSystem compute inherent $pH $segresid1]
+        set du2 [cphSystem compute inherent $pH $segresid2]
+        set accept [metropolisAcceptance [expr {$du1 + $du2}]]
+    } else {
+        set accept 0
+    }
+    return $accept
+}
+
 proc ::cphTitrator::proposeMove {} {
     variable ::cphTitrator::maxAttempts
     variable ::cphTitrator::MoveSet
index a4a0ab7..282f56f 100644 (file)
@@ -134,13 +134,46 @@ proc ::cphSystem::proposeResidueTitration {segresid} {
         set numProtonsTrial [expr {lsum($occ)}]
         if {$numProtons != $numProtonsTrial} {
             cphSystem set trialState $segresid $state
-            return 0
+            return 1
         }
     }
     # This is an error and should never happen.
     return -1
 }
 
+proc ::cphSystem::proposeProtonTransfer {segresid1 segresid2} {
+    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 dn [expr {$numProtons1 - $numProtons2}]
+
+    if {$dn == 0.0} {     
+        # No proton transfer is possible.
+        return 0
+    } 
+    # transfer from 1 to 2
+    while {true} {
+        set state1 [choice $possibleStates1]
+        set occ1 [dict get $resDefDict $resname1 states $state1]
+        set numProtonsTrial1 [expr {lsum($occ1)}]
+        set state2 [choice $possibleStates2]
+        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
+            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 
@@ -307,25 +340,29 @@ proc ::cphSystem::computeInherentNormedWeights {pH segresid} {
 # computed, not the temperature of the simulation. See 
 # computeInherentAcceptance for definition of the remaining notation.
 #
-proc ::cphSystem::computeSwitchAcceptance {segresid} {
-    set l [cphSystem get occupancy $segresid]
-    set lp [cphSystem get trialOccupancy $segresid]
-    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]
-    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]
+proc ::cphSystem::computeSwitchAcceptance {segresidList} {
+    set du 0.0
+    foreach segresid $segresidList {
+        set l [cphSystem get occupancy $segresid]
+        set lp [cphSystem get trialOccupancy $segresid]
+        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]
+        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]
+        }
+        set kT [expr {$::BOLTZMANN*[cphSystem get Tref $segresid]}]
+        set du [expr {$du + $sgn*($dG - $kT*$::LN10*($pKa - $pKai))}]
     }
-    set kT [expr {$::BOLTZMANN*[cphSystem get Tref $segresid]}]
-    return [expr {$sgn*($dG - $kT*$::LN10*($pKa - $pKai))}]
+    return $du
 }
 
 # ---------------
index 516d0c1..e4d52aa 100644 (file)
@@ -30,6 +30,7 @@ namespace eval ::namdcph {
     dict set moveInfo maxProposalAttempts 0
     dict set moveInfo default numsteps 20 
     dict set moveInfo default weight 1.0
+    dict set moveInfo ptransfer [list]
 
     variable residueAliases [dict create]
     dict set residueAliases HIS {HSD HSE HSP}
@@ -389,6 +390,14 @@ proc ::namdcph::cphProposalWeight {args} {
     return
 }
 
+proc ::namdcph::cphProposeProtonTransfer {args} {
+    variable ::namdcph::moveInfo
+    foreach moveLabel $args {
+        dict lappend moveInfo ptransfer $moveLabel
+    }
+    return
+}
+
 # ::namdcph::cphMaxProposalAttempts
 #
 # Number of attempted MC proposals from each move set before giving up.
@@ -517,7 +526,7 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
     set DeltaE [cphSystem compute switch $segresidList]
     set Work [expr {$::energyArray(CUMALCHWORK) + $DeltaE}]
     set ReducedWork [expr {$Work / ($::BOLTZMANN*[$::thermostatTempCmd])}]
-    set accept [metropolisAcceptance $ReducedWork] 
+    set accept [metropolisAcceptance $ReducedWork]
 #    printProposalSummary $segresidList
     set tmp [printProposalSummary $segresidList]
     cphPrint [format "%s WorkCorr % 10.4f CorrWork % 10.4f"\