Fix HMR bug in constant-pH MD 63/4363/1
authorradakb <brian.radak@gmail.com>
Tue, 17 Jul 2018 14:59:01 +0000 (10:59 -0400)
committerradakb <brian.radak@gmail.com>
Tue, 17 Jul 2018 14:59:01 +0000 (10:59 -0400)
This patch makes sure that constant-pH MD does not undo mass
repartitioning within titratable residues.

Also, a warning is now issued whenever a proton is automagically
inserted during setup. This is because it is impossible to
on-the-fly repartition mass on those protons and thus the insertions
will not match any pre-existing HMR.

Unfortunately, this means that HMR constant-pH MD will be an
advanced feature for now, since one cannot rely on the auto-setup
within constant pH.

Change-Id: I0243f2e9a781b923a9b5bd83b68b867a0a6fe701

lib/namdcph/namdcph/cphtoppar.tcl

index 83ba5f7..a1b627c 100644 (file)
@@ -576,7 +576,28 @@ proc ::cphSystem::initializeSystem {pH temperature buildH stateInfo} {
     # Final pass - Apply the patches
     foreach segresidname [cphSystem get segresidnames] {
         lassign [split $segresidname ":"] segid resid
+        lassign [cphSystem get alchAtomLists $segresidname] atoms
+        set masses [list]
+        foreach atom $atoms {
+            if {[catch {eval segment mass $segid $resid $atom}]} {
+                lappend masses -1.0
+            } else {
+                lappend masses [segment mass $segid $resid $atom]
+            }
+        }
+        # The patch statement reassigns all masses based on RTF definitions.
+        # This must be undone in the case of HMR.
         patch [cphSystem get statePatch $segresidname] "$segid:$resid"
+        foreach atom $atoms mass $masses {
+            if {$mass > 0.0} {
+                psfset mass $segid $resid $atom $mass
+            } else {
+                print "WARNING! Building new atom $segid:$resid:$atom with"\
+                      "unknown mass."
+                print "         This may cause issues, for example, with"\
+                      "mass repartitioning."
+            }
+        }
     }
     guesscoord
     foreach segresidname [cphSystem get segresidnames] {