Constant-pH - more config and restart format updates 20/3920/1
authorradakb <brian.radak@gmail.com>
Thu, 29 Mar 2018 19:27:34 +0000 (15:27 -0400)
committerradakb <brian.radak@gmail.com>
Thu, 29 Mar 2018 19:27:34 +0000 (15:27 -0400)
- Residue aliases are now part of the config file definition,
  rather than being hardcoded ad hoc. THIS MEANS YOU WILL NEED TO
  UPDATE YOUR JSON FILES FROM HEREON.

- Residue exclusions are now stored in the restart file, which
  means you no longer have to keep using cphExcludeResidue commands
  during subsequent runs (it would have crashed before if you
  didn't).

TODO: All of the changes to the restart format are starting to
  make the writer a bit heavy and over complicated. It looks like
  I'll finally have to sit down and write some general JSON
  encoding procs...

Change-Id: Ia8868b0cec9da0348bfa31f9f2a82346b1f7b6d6

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

index 9d79187..68606da 100644 (file)
@@ -531,7 +531,7 @@ proc ::cphSystem::initializeSystem {pH temperature buildH stateInfo} {
 # Build the residue definitions and residue objects for the system based on 
 # the NAMD inputs.
 #
-proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
+proc ::cphSystem::buildSystem {resDefs excludeList} {
     variable ::cphSystem::resDict
     variable ::cphSystem::resDefDict [checkResDefs $resDefs]
 
@@ -563,6 +563,14 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
         }
     }
 
+    # 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 {
+        if {![dict exists $resDef aliases]} continue
+        dict set resAliases $resname [dict get $resDef aliases]
+    }
+
     # Read in whatever files were specified to NAMD.
     set Args [list [structure] pdb [coordinates]]
     if {[isset binCoordinates]} {
@@ -594,53 +602,82 @@ proc ::cphSystem::buildSystem {resDefs resAliases segresExcls} {
 
             set resIsDefined [expr {$resname in $definedResidues}]          
             set resIsSubtyped [expr {$resname in $definedSubResidues}]
-            set resIsExcluded [expr {$segresidname in $segresExcls}]
-            # Break here if nothing to do.
+            set resIsExcluded [expr {$segresidname in $excludeList}]
+            # Break here if nothing to do. Be sure to remove successful
+            # exclusions from the list so that we can identify bad selections
+            # later.
+            if {$resIsExcluded} {
+                set excludeList [lsearch -inline -all -not $excludeList\
+                        $segresidname]
+            }
             if {(!$resIsDefined && !$resIsSubtyped) || $resIsExcluded} {
                 continue
             }
 
-            if {$resIsDefined} { 
+            # Add the residue to the system.
+            if {$resIsDefined} {
                 dict set resDict $segresidname [dict create]
                 dict set resDict $segresidname resname $resname
             }
 
+            # If the residue is subtype-able, check that subresidues exist.
+            # Note that this might be true even if the residue itself is not
+            # titratable.
             if {!$resIsSubtyped} continue
-            # The conditions for subtyping may still not be met...
+
             set atoms [segment atoms $segid $resid]
             foreach subresname [dict get $subresDefDict $resname] {
+                # Atoms that MUST exist.
+                set allSubatomsExist 1
                 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 {
+                foreach atom $subatoms { 
                     if {$atom ni $atoms} {
                         set allSubatomsExist 0
                         break
                     }
                 }
+
+                # Atoms that must NOT exist.
                 set allNotSubatomsDoNotExist 1
+                set notsubatoms [list]
+                if {[dict exists $resDefDict $subresname notsubatoms]} {
+                    set notsubatoms\
+                            [dict get $resDefDict $subresname notsubatoms]
+                }
                 foreach atom $notsubatoms {
                     if {$atom in $atoms} {
                         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
+
+                # Break here if anything doesn't match. 
+                if {!$allSubatomsExist || !$allNotSubatomsDoNotExist} continue
+
+                set segresidname [format "%s:%s" $segresid $subresname]
+                # After all of that, the residue might be excluded!
+                set resIsExcluded [expr {$segresidname in $excludeList}]
+                if {$resIsExcluded} {
+                    set excludeList [lsearch -inline -all -not $excludeList\
+                            $segresidname]
+                    continue
                 }
+
+                # Add the sub-residue to the system.
+                dict set resDict $segresidname [dict create]
+                dict set resDict $segresidname resname $subresname
             }
         }
     }
+    # This is a bit of a hack way to check that all of the selected exclusions
+    # were valid, but it works.
+    if {[llength $excludeList]} {
+        abort "Cannot exclude non-existant residue(s): $excludeList"
+    }
+
     return $resDict
 }
 
index 4d32a93..d0fb046 100644 (file)
@@ -25,9 +25,6 @@ namespace eval ::namdcph {
     variable stateInfo [dict create] ;# cphSystem data and defaults
     variable moveInfo [dict create] ;# cphTitrator data and defaults
 
-    variable residueAliases [dict create]
-    dict set residueAliases HIS {HSD HSE HSP}
-
     namespace export *
 }
 
@@ -683,7 +680,6 @@ proc ::namdcph::readConfig {} {
 proc ::namdcph::readRestart {} {
     variable ::namdcph::restartFilename
     variable ::namdcph::SystempH
-    variable ::namdcph::stateInfo
     variable ::namdcph::moveInfo
 
     set RestartFile [open $restartFilename "r"]
@@ -694,39 +690,20 @@ proc ::namdcph::readRestart {} {
         pH [dict get $Restart pH]
     }
 
+    if {[dict exists $Restart exclude]} {
+        cphExcludeResidue {*}[dict get $Restart exclude]
+    }
+
     # Read state parameters, if present. 
     #
+    set stateList {}
     if {[dict exists $Restart states]} {
         set stateList [dict get $Restart states]
     }
+    set pKaiList {}
     if {[dict exists $Restart pKais]} {
         set pKaiList [dict get $Restart pKais]
     }
-    if {[info exists stateList] && [info exists pKaiList]} {
-        set segresidnameList [cphSystem get segresidnames]
-
-        if {[llength $stateList] != [llength $pKaiList]} {
-            cphAbort "mismatch in states/pKais in $restartFilename"
-        }
-        if {[llength $stateList] != [llength $segresidnameList]} {
-            cphAbort "Too few/many state/pKai definitions in $restartFilename"
-        }
-
-        foreach segresidname $segresidnameList\
-                state $stateList pKai $pKaiList {
-            if {![dict exists $stateInfo $segresidname]} {
-                dict set stateInfo $segresidname state $state
-                dict set stateInfo $segresidname pKai $pKai
-            } else {
-                if {![dict exists $stateInfo $segresidname state]} {
-                    dict set stateInfo $segresidname state $state
-                }
-                if {![dict exists $stateInfo $segresidname pKai]} {
-                    dict set stateInfo $segresidname pKai $pKai
-                }
-            }
-        }
-    }
 
     # Read MC move parameters, if present. Only reset those values which have
     # not been set in via keywords.
@@ -754,7 +731,7 @@ proc ::namdcph::readRestart {} {
         }
     }
 
-    return
+    return [list $stateList $pKaiList]
 }
 
 # ::namdcph::writeRestart
@@ -777,6 +754,7 @@ proc ::namdcph::readRestart {} {
 proc ::namdcph::writeRestart {args} { 
     variable ::namdcph::restartFreq
     variable ::namdcph::SystempH
+    variable ::namdcph::excludeList
 
     if {[string match [lindex $args 0] force]} {
         set restartFilename [lindex $args 1]
@@ -795,6 +773,15 @@ proc ::namdcph::writeRestart {args} {
     set cycleStr "\"cycle\":$cycle"
     set pHStr "\"pH\":$SystempH"
 
+    set exclusionStr ""
+    if {[llength $excludeList]} {
+        set exclusionStr "\"exclude\":\["
+        foreach segresidname $excludeList {
+            set exclusionStr "$exclusionStr\"$segresidname\","
+        }
+        set exclusionStr "[string trimright $exclusionStr ","]\]"
+    }
+
     set stateStr "\"states\":\["
     foreach state [cphSystem get state] {
         set stateStr "$stateStr\"$state\","
@@ -829,7 +816,7 @@ proc ::namdcph::writeRestart {args} {
     }
     set MCStr "[string trimright $MCStr ","]\}"
 
-    puts $RestartFile "\{$cycleStr,$pHStr,$stateStr,$pKaiStr,$MCStr\}"
+    puts $RestartFile "\{$cycleStr,$pHStr,$exclusionStr,$stateStr,$pKaiStr,$MCStr\}"
     close $RestartFile
     # Always checkpoint the PSF and PDB when a restart is written.
     file copy -force [mdFilename psf] "[outputName].psf"
@@ -936,7 +923,6 @@ proc ::namdcph::initialize {} {
     variable ::namdcph::configFilenames
     variable ::namdcph::stateInfo
     variable ::namdcph::moveInfo
-    variable ::namdcph::residueAliases
     variable ::namdcph::excludeList
     variable ::namdcph::SystempH
 
@@ -952,18 +938,45 @@ proc ::namdcph::initialize {} {
     # rebuild the topology to reflect those states.
     #
     cphPrint "initializing constant pH PSF..."
-    # TODO: integrate aliases into the configfile standard...
     if {![llength configFilenames]} {
         cphAbort "At least one constant pH configuration file is required."
     }
-    cphSystem build [readConfig] $residueAliases $excludeList
 
     if {[string length $restartFilename]} {
-        readRestart
+        lassign [readRestart] stateList pKaiList
         lassign [list 0.0 false] temp buildH
     } else {
         lassign [list [$::thermostatTempCmd] true] 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]} {
+        set segresidnameList [cphSystem get segresidnames]
+
+        if {[llength $stateList] != [llength $pKaiList]} {
+            cphAbort "mismatch in states/pKais in $restartFilename"
+        }
+        if {[llength $stateList] != [llength $segresidnameList]} {
+            cphAbort "Too few/many state/pKai definitions in $restartFilename"
+        }
+
+        foreach segresidname $segresidnameList\
+                state $stateList pKai $pKaiList {
+            if {![dict exists $stateInfo $segresidname]} {
+                dict set stateInfo $segresidname state $state
+                dict set stateInfo $segresidname pKai $pKai
+            } else {
+                if {![dict exists $stateInfo $segresidname state]} {
+                    dict set stateInfo $segresidname state $state
+                }
+                if {![dict exists $stateInfo $segresidname pKai]} {
+                    dict set stateInfo $segresidname pKai $pKai
+                }
+            }
+        }
+    }
 
     if {![info exists SystempH]} {
         cphAbort "A pH value is required."