Updates and fixes to constant pH 68/3768/1
authorradakb <brian.radak@gmail.com>
Fri, 23 Feb 2018 19:13:21 +0000 (14:13 -0500)
committerradakb <brian.radak@gmail.com>
Fri, 23 Feb 2018 19:13:21 +0000 (14:13 -0500)
- Add basic support for titratable phosphate monoesters
- Fix potential bug when using mass repartitioning
- Remove extraneous code:
  o checkpoint/revert is redundant with reinitatoms
  o bfactors are already zeroed by psfgen
  o clean up old-style proposal summaries

Change-Id: I670aab3eb3630edf9ee8d622707f83951590881e

lib/namdcph/namdcph/cphpsfgen.tcl
lib/namdcph/namdcph/cphtoppar.tcl
lib/namdcph/namdcph/namdcph.core.tcl

index c5ac655..c093270 100644 (file)
@@ -89,6 +89,7 @@ proc alchPatch {patch segresid l0atomList l1atomList
     set NSkipH 0
     foreach l0atom $l0atomList l1atom $l1atomList {
         set Elem [string index $l0atom 0]
+        set Mass [segment mass $segid $resid $l0atom]
         # Coordinates
         if {($buildH) && ($Elem == "H") 
             && (($PrevElem == "H") || ($PrevElem == "C") 
@@ -102,7 +103,6 @@ proc alchPatch {patch segresid l0atomList l1atomList
         set PrevElem $Elem
         # Velocities
         if {$temp > 0} {
-            set Mass [segment mass $segid $resid $l0atom]
             set SigmaV [expr {($Mass > 0) ? [expr {sqrt($kT/$Mass)}] : 0.0}]
             set v {0.0 0.0 0.0}
             AddNoiseToVector v 0.0 $SigmaV
@@ -110,6 +110,7 @@ proc alchPatch {patch segresid l0atomList l1atomList
             set v [segment velocities $segid $resid $l0atom]
         }
         psfset vel $segid $resid $l1atom $v
+        psfset mass $segid $resid $l1atom $Mass
         # Alchemical Groups
         psfset beta $segid $resid $l0atom -1.0
         psfset beta $segid $resid $l1atom 1.0
index 282f56f..f3d73ab 100644 (file)
@@ -467,11 +467,6 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
                 psfset resname $segid $resid $realName
                 set resname $realName
             }
-            # Make sure we have no stray nonzero B-factors
-            # TODO: Check that this is really ok and warn the user?
-            foreach atom [segment atoms $segid $resid] {
-                psfset beta $segid $resid $atom 0.0
-            }
             # 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
@@ -1140,6 +1135,19 @@ proc ::cphSystem::resDef2Matrix {resDef data} {
                 lset Matrix [index2flatindex 6 5] $attr65
                 lset Matrix [index2flatindex 6 3] $attr63
                 lset Matrix [index2flatindex 5 3] $attr53
+            } elseif {![lindex $nprotonsExists 2]
+                      && ![lindex $nprotonsExists 3]} {
+                # state (1, 1, 1) and 2 proton states missing
+                switch -- $numValues 1 {;# Ex. phosphate monoesters
+                    lassign $data attr10
+                    set attr20 $attr10
+                    set attr40 $attr10
+                } default {
+                    set errorCode -1
+                }
+                lset Matrix [index2flatindex 1 0] $attr10
+                lset Matrix [index2flatindex 2 0] $attr20
+                lset Matrix [index2flatindex 4 0] $attr40
             } else {
                 set errorCode -1
             }
index e4d52aa..be7bcd5 100644 (file)
@@ -68,7 +68,6 @@ proc ::namdcph::cphRun {numsteps {numcycles 1}} {
             cphPrint "All proposals rejected ($nattempts total)."
             set runArgs [list norepeat $numsteps]
         } else {
-#            printProposalSummary $segresidList
             cphPrint "Proposal accepted ($nattempts), attemping a switch."
             set accept [runSwitch $swNumsteps $segresidList]
             set runArgs [list $numsteps]
@@ -501,7 +500,6 @@ proc ::namdcph::runMD {args} {
 proc ::namdcph::runSwitch {numsteps segresidList} {
     # (1) Checkpoint and modify output parameters. 
     #
-    checkpoint
     storeEnergies
     set storedOutputEnergies [outputEnergies]
     set storedDCDFreq [dcdFreq]
@@ -527,7 +525,6 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
     set Work [expr {$::energyArray(CUMALCHWORK) + $DeltaE}]
     set ReducedWork [expr {$Work / ($::BOLTZMANN*[$::thermostatTempCmd])}]
     set accept [metropolisAcceptance $ReducedWork]
-#    printProposalSummary $segresidList
     set tmp [printProposalSummary $segresidList]
     cphPrint [format "%s WorkCorr % 10.4f CorrWork % 10.4f"\
             [join $tmp "/"] $DeltaE $Work]
@@ -546,7 +543,6 @@ proc ::namdcph::runSwitch {numsteps segresidList} {
     } else {
         cphPrint "Switch rejected!"
         alch off
-        revert
         reloadAndReinit [getMDBasename] false
     }
     return $accept
@@ -1070,7 +1066,6 @@ proc ::namdcph::printProposalSummary {segresidList} {
         set state [cphSystem get state $segresid]
         set trialState [cphSystem get trialState $segresid]
         lappend retList "$reslabel:$state:$trialState"
-#        cphPrint "Switch $reslabel $state --> $trialState"
     }
     return $retList
 }