Permit negative pH and pKai
[namd.git] / lib / namdcph / namdcph / namdcph.core.tcl
1
2 # Core utilities for setting up and running constant pH simulations in NAMD.
3 #
4 package require Tcl 8.5
5
6 source [file join [file dirname [info script]] "namdtcl.tcl"]
7 source [file join [file dirname [info script]] "namdmcmc.tcl"]
8 source [file join [file dirname [info script]] "cphtitrator.tcl"]
9 source [file join [file dirname [info script]] "json.tcl"]
10
11 namespace eval ::namdcph {
12     namespace import ::cphTitrator::*
13
14     variable TITLE "namdcph)"
15     variable SystempH
16     variable configFilenames [list]
17     variable restartFilename ""
18     variable restartFreq 0
19     variable outFile ""
20     variable numMinSteps 0
21     variable excludeList [list]
22     variable alchFrcCons 100.0
23     variable useHeat 0
24     variable MDBasename namdcph.md
25     variable SWBasename namdcph.sw
26     variable stateInfo [dict create] ;# cphSystem data and defaults
27     variable moveInfo [dict create] ;# cphTitrator data and defaults
28
29     namespace export *
30 }
31
32 # =============================================================================
33 # Main NAMD Routines
34 # =============================================================================
35 # ::namdcph::cphRun
36 #
37 # Run a set number of neMD/MC cycles.
38 #
39 proc ::namdcph::cphRun {numsteps {numcycles 1}} {
40     variable ::namdcph::SystempH
41
42     # Initialize NAMD and build a constant pH enabled PSF.
43     initialize
44     finalize
45     if {$::namdcph::numMinSteps > 0} {
46         set storedFirstTimeStep [firstTimeStep]
47         minimize $::namdcph::numMinSteps
48         if {[isset temperature]} {
49           reinitvels [temperature]
50         } else {
51           reinitvels [$::thermostatTempCmd]
52         }
53         firstTimeStep $storedFirstTimeStep
54     }
55     set cphlog [openCpHLog]
56     set firstCycle 1
57     set lastCycle [expr {$firstCycle + $numcycles - 1}] 
58     set runArgs [list $numsteps]
59     for {set cycle $firstCycle} {$cycle <= $lastCycle} {incr cycle} { 
60         # (1) Perform whatever equilibrium sampling is desired.
61         #
62         runMD {*}$runArgs
63         # (2) Propose a move from the full move set.
64         #
65         lassign [cphTitrator propose $SystempH] accept swNumsteps\
66                 segresidnameList
67         # (3) At this point we have either selected a switch or rejected a
68         # whole bunch of moves. If the former, perform the switch.
69         #
70         if {!$accept} {
71             cphPrint "Proposal rejected."
72             set runArgs [list norepeat $numsteps]
73         } else {
74             cphPrint "Proposal accepted, attemping a switch."
75             set accept [runSwitch $swNumsteps $segresidnameList]
76             set runArgs [list $numsteps]
77             # Only accumulate statistics for attempted switches.
78             set moveLabel [join $segresidnameList "/"]
79             cphTitrator accumulateStats $accept $moveLabel
80         }
81         cphSystem update $accept $segresidnameList
82         # (4) Output cycle information and a checkpoint file if necessary.
83         #
84         puts $cphlog "[format "%6d" $cycle] [cphSystem get occupancy]"
85         flush $cphlog
86         writeRestart "[outputname].cphrst" $cycle
87     }
88     writeRestart force "[outputname].cphrst" $cycle
89     # Cleanup temporary files
90     file delete {*}[glob [getSWBasename].*]
91     file delete {*}[glob [getMDBasename].*]
92     close $cphlog
93     printTitratorReport
94     return
95 }
96
97 #
98 # FOR ADVANCED USE ONLY!!!
99 #
100 # ::namdcph::checkResidueDefinitions
101 #
102 # Check that appropriate residue definitions are defined for this residue.
103 #
104 # This essentially just scans all possible states for each titratable residue,
105 # of a given type in the system. THIS OFFERS NO VALIDATION AS TO WHETHER OR NOT
106 # THE DEFINITIONS MAKE SENSE, it simply checks that the exist.
107 #
108 proc ::namdcph::checkResidueDefinitions {resnames} {
109     cphWarnExpt
110     pH 7.0
111     cphForceConstant 0.0
112     # Initialize NAMD and build a constant pH enabled PSF.
113     initialize
114     finalize
115     $::thermostatTempCmd 0.0
116     outputEnergies 1
117     alchLambdaFreq 0
118     foreach segresidname [cphSystem get segresidnames] {
119         lassign [split $segresidname ":"] segid resid resname
120         if {$resname ni $resnames} {
121             continue
122         } 
123         set states [cphSystem get stateList $segresidname]
124         foreach state0 $states {
125             foreach state1 $states {
126                 if {[string match $state0 $state1]} {
127                     continue
128                 }
129                 # Force the initial state as state0
130                 output [getMDBasename]
131                 psfgenRead [getMDBasename]
132                 cphSystem initializeState $segresidname $state0
133                 alchemify $segresidname
134                 dealchemify $segresidname
135                 cphSystem update 1 $segresidname
136                 # Run an instantaneous switch
137                 runTestSwitch $segresidname $state1
138             }
139         }
140     }
141     return 
142 }
143
144 #
145 # FOR ADVANCED USE ONLY!!!
146 #
147 # ::namdcph::testResidueSwitch
148 #
149 # Run an instantenous switch for the given residue between the given states.
150 #
151 # This is meant as an interface function for checking energies against
152 # non-constant pH calculations, but cannot perform that comparison on its own.
153 #
154 proc ::namdcph::testResidueSwitch {segresidname state0 state1} {
155     cphWarnExpt
156     pH 7.0
157     cphNumstepsPerSwitch 20
158     cphForceConstant 0.0
159     # Initialize NAMD and build a constant pH enabled PSF.
160     cphSetResidueState $segresidname $state0
161     initialize
162     finalize
163     $::thermostatTempCmd 0.0
164     outputEnergies 1
165     alchLambdaFreq 0
166     return [runTestSwitch $segresidname $state1]
167 }
168
169 #
170 # FOR ADVANCED USE ONLY!!!
171 #
172 # ::namdcph::runTestSwitch
173 #
174 # This is just an internal convenience function for running instantaneous test
175 # switches.
176 #
177 # See checkResidueDefinitions and testResidueSwitch
178 #
179 proc ::namdcph::runTestSwitch {segresidname state1} {
180     # You can't make lists of arrays or arrays of arrays, so the return type
181     # has to be a dict of arrays (callback only returns arrays).
182     set retVals [dict create]
183     run 0
184     storeEnergies
185     dict set retVals MDEnergies0 [array get ::energyArray]
186     cphSystem set trialState $segresidname $state1
187     alchemify $segresidname
188     alchLambda 0.0
189     run 0
190     storeEnergies
191     dict set retVals SWEnergies0 [array get ::energyArray]
192     alchLambda 1.0
193     run 0
194     storeEnergies
195     dict set retVals SWEnergies1 [array get ::energyArray]
196     dealchemify $segresidname
197     cphSystem update 1 $segresidname
198     run 0
199     storeEnergies
200     dict set retVals MDEnergies1 [array get ::energyArray]
201     # Cleanup temporary files
202     file delete {*}[glob [getSWBasename].*]
203     file delete {*}[glob [getMDBasename].*]
204     return $retVals
205 }
206
207 #
208 # FOR ADVANCED USE ONLY!!!
209 #
210 # ::namdcph::cphAnalyzeForce
211 #
212 # Test one or more residue definitions.
213 #
214 proc ::namdcph::cphAnalyzeForce {dcdfilename segresidname state0 state1} {
215     cphWarnExpt
216     pH 7.0
217
218     cphSetResidueState $segresidname $state0
219     initialize
220     set initialPSF [structure]
221     set initialPDB [coordinates]
222     finalize
223
224     set numsteps 500
225     outputEnergies [expr {$numsteps == 0 ? 1 : int($numsteps)}]
226     alchLambdaFreq [expr {$numsteps == 0 ? 0 : 1}]
227
228     set nframes -1
229     coorfile open dcd $dcdfilename
230     while {![coorfile read]} {
231         incr nframes
232
233         # We have to do this so that inputs can be correctly loaded...
234         set basename [format "%s.%d" [getMDBasename] $nframes]
235         output $basename
236         file copy -force $initialPSF "$basename.psf"
237         file copy -force $initialPDB "$basename.pdb"
238         reloadAndReinit $basename false
239         # Assign the correct state and build protons or dummy atoms
240         cphSystem initializeState $segresidname $state0
241         alchemify $segresidname
242         dealchemify $segresidname
243         cphSystem update 1 $segresidname
244         # Now build the alchemical atoms
245         cphSystem set trialState $segresidname $state1
246         alchemify $segresidname
247         firsttimestep 0
248         run $numsteps
249         if {$numsteps} {
250             storeEnergies
251             set DeltaE [cphSystem compute switch $segresidname]
252             set Work [expr {$::energyArray(CUMALCHWORK) + $DeltaE}]
253             set ReducedWork [expr {$Work / [::kBT]}]
254             set tmp [printProposalSummary $segresidname]
255             cphPrint [format "%s WorkCorr % 10.4f CorrWork % 10.4f"\
256                     [join $tmp "/"] $DeltaE $Work]
257         }
258         dealchemify $segresidname
259         file delete {*}[glob $basename.*]
260     }
261     file delete {*}[glob [getMDBasename].*]
262     file delete {*}[glob [getSWBasename].*]
263     coorfile close
264     return
265 }
266
267 # =============================================================================
268 # "Setter" Routines - used as new keywords in NAMD
269 #
270 #   All of these procedures take the desired value as an argument and return 
271 # that value. For default values see the decalaration in the namespace 
272 # header.
273 # =============================================================================
274 # -----------------
275 # Required Keywords
276 # -----------------
277 # ::namdcph::pH
278 #
279 # pH value for the simulation
280 #
281 proc ::namdcph::pH {pHValue} {
282     checkIsNumeric "pH" $pHValue
283     variable ::namdcph::SystempH $pHValue 
284     return
285 }
286
287 # ::namdcph::cphConfigFile
288 #
289 # Configuration file for constant pH residues.
290 #
291 proc ::namdcph::cphConfigFile {filename} {
292     variable ::namdcph::configFilenames
293     lappend configFilenames [string trim $filename]
294     return
295 }
296
297 # ::namdcph::cphNumstepsPerSwitch
298 #
299 # For odd arguments, the first argument is assumed to be a default switch time.
300 # All remaining arguments are presumed to be label/numsteps pairs for specific
301 # moves.
302 #
303 proc ::namdcph::cphNumstepsPerSwitch {args} {
304     variable ::namdcph::moveInfo
305     if {[expr {[llength $args] % 2}]} {
306         set numsteps [lindex $args 0]
307         checkIsNotNegative numSwitchSteps $numsteps
308         dict set moveInfo default numsteps [expr {int($numsteps)}]
309         set args [lrange $args 1 end]
310     }
311     checkArglistIsMultiple $args 2
312     for {set i 0} {$i < [llength $args]} {incr i 2} {
313         lassign [lrange $args $i [expr {$i+1}]] moveLabel numsteps
314         checkIsNotNegative numSwitchSteps $numsteps
315         dict set moveInfo $moveLabel numsteps [expr {int($numsteps)}]
316     }
317     return
318 }
319
320 # ---------------------
321 # Commonly Used Options
322 # ---------------------
323 # ::namdcph::cphSetResidueState
324 #
325 # Set the state of one or more residues using psfgen syntax.
326 #
327 proc ::namdcph::cphSetResidueState {args} {
328     checkArglistIsMultiple $args 2
329     variable ::namdcph::stateInfo
330     for {set i 0} {$i < [llength $args]} {incr i 2} {
331         lassign [lrange $args $i [expr {$i+1}]] segresidname state
332         dict set stateInfo $segresidname state $state
333     }
334     return
335 }
336
337 # ::namdcph::cphSetResiduepKai
338 #
339 # Set the inherent pKa of one or more residues using psfgen syntax.
340 #
341 proc ::namdcph::cphSetResiduepKai {args} {
342     checkArglistIsMultiple $args 2
343     variable ::namdcph::stateInfo
344     for {set i 0} {$i < [llength $args]} {incr i 2} {
345         #NB pKai may be a list of values - check individually.
346         lassign [lrange $args $i [expr {$i+1}]] segresidname pKai
347         foreach pKa $pKai {
348             checkIsNumeric pKai $pKa
349         }
350         dict set stateInfo $segresidname pKai $pKai
351     }
352     return
353 }
354
355 # ::namdcph::cphRestartFile 
356 #
357 # Restart a constant pH run from a restart file.
358 #
359 proc ::namdcph::cphRestartFile {filename} {
360     variable ::namdcph::restartFilename $filename
361     return
362 }
363
364 # ::namdcph::cphRestartFreq
365 #
366 # Frequency (in neMD/MC cycles) at which to save constant pH restart files
367 #
368 proc ::namdcph::cphRestartFreq {frequency} {
369     checkIsNotNegative cphRestartFreq $frequency
370     variable ::namdcph::restartFreq $frequency
371     return
372 }
373
374 # ::namdcph::cphOutFile
375 #
376 # Name for constant pH output file - default is [outputname].cphlog
377 #
378 proc ::namdcph::cphOutFile {filename} {
379     variable ::namdcph::outFile $filename
380     return
381 }
382
383 # ::namdcph::cphProposalWeight
384 #
385 # The (unnormalized) proposal weight assigned to each move. The default is for
386 # all such weights to be equal (uniformally distributed moves).
387 #
388 proc ::namdcph::cphProposalWeight {args} {
389     variable ::namdcph::moveInfo
390     checkArglistIsMultiple $args 2
391     for {set i 0} {$i < [llength $args]} {incr i 2} {
392         lassign [lrange $args $i [expr {$i+1}]] moveLabel weight
393         checkIsNotNegative proposalWeight $weight
394         dict set moveInfo $moveLabel weight [expr {1.*$weight}]
395     }
396     return
397 }
398
399 # ::namdcph::cphProposeProtonTransfer
400 #
401 # Add additional "proton transfer" moves to the move set. These new moves are
402 # only attempted if a proton can be swapped between the given residues. The
403 # two residues are given by their normal <segid>:<resid>:<resname> designation
404 # except that they are further separated by a "/".
405 #
406 proc ::namdcph::cphProposeProtonTransfer {args} {
407     variable ::namdcph::moveInfo
408     foreach moveLabel $args {
409         dict set moveInfo $moveLabel protonTransfer 1
410     }
411     return
412 }
413
414 # ::namdcph::cphProposeCotitration
415 #
416 # Add additional "co-titration" moves to the move set. These new moves are only
417 # attempted if two residues are both protonated or deprotonated. They will then
418 # attempt to change states concurrently. The two residues are given by their
419 # normal <segid>:<resid>:<resname> designation except that they are further
420 # separated by a "/".
421 #
422 proc ::namdcph::cphProposeCotitration {args} {
423     variable ::namdcph::moveInfo
424     foreach moveLabel $args {
425         dict set moveInfo $moveLabel cotitration 1
426     }
427     return
428 }
429
430 # ::namdcph::cphMaxProposalAttempts
431 #
432 # Number of attempted MC proposals from each move set before giving up.
433 #
434 # Values less than 1 default to the number of residues in the system.
435 # NB: This is not the necessarily the same thing as attempting a move once for
436 #     each residue.
437 #
438 proc ::namdcph::cphMaxProposalAttempts {maxAttempts} {
439     checkIsNumeric cphMaxProposalAttempts $maxAttempts
440     variable ::namdcph::moveInfo
441     dict set moveInfo maxProposalAttempts [expr {int($maxAttempts)}]
442     cphWarn "cphMaxProposalAttempts is now deprecated and set to 1!"
443     return
444 }
445
446 # ::namdcph::cphNumMinSteps
447 #
448 # Number of minimization steps to perform before dynamics.
449 #
450 # This is especially useful when the initial states are randomized according to
451 # the pH.
452 #
453 proc ::namdcph::cphNumMinSteps {numsteps} {
454     checkIsNotNegative cphNumMinSteps $numsteps
455     variable ::namdcph::numMinSteps $numsteps
456     return
457 }
458
459 proc ::namdcph::cphExcludeResidue {args} {
460     variable ::namdcph::excludeList
461     foreach segresidname $args {
462          lappend excludeList $segresidname
463     }
464     return
465 }
466
467 # -------------------
468 # Specialized Options
469 # -------------------
470 # ::namdcph::cphForceConstant
471 #
472 # Force constant for zero-length bonds between alchemical atoms in kcal/mol-A^2
473 #
474 proc ::namdcph::cphForceConstant {forceConstant} {
475     checkIsNotNegative cphForceConstant $forceConstant
476     variable ::namdcph::alchFrcCons $forceConstant
477     return
478 }
479
480 # ::namdcph::cphMDBasename
481 #
482 # Basename for (temporary) constant pH MD input files - the only output of
483 # interest to the user should be in the usual places
484 #
485 proc ::namdcph::cphMDBasename {basename} {
486     variable ::namdcph::MDBasename $basename
487     return
488 }
489
490 # ::namdcph::cphSwitchBasename
491 #
492 # Basename for (temporary) constant pH switch trajectory output files
493 #
494 proc ::namdcph::cphSwitchBasename {basename} {
495     variable ::namdcph::SWBasename $basename
496     return
497 }
498
499 # =============================================================================
500 # Nonequilibrium Switching Routines
501 # =============================================================================
502 # ::namdcph::runMD
503 proc ::namdcph::runMD {args} {
504    run {*}$args 
505 }
506
507 # ::namdcph::runSwitch
508 #
509 # Run a neMD/MC switch on the specified residue(s).
510 #
511 # This includes multiple steps:
512 # (1) Modify I/O settings to differentiate from regular MD
513 # (2) Build alchemically enabled side chains for the desired residues
514 # (3) Perform a switching trajectory with appropriate momentum reversals
515 # (4) Perform Metropolis MC using an appropriate work quantity
516 # (5) Reinitialize MD based on the MC result
517 #
518 # Arguments:
519 # ----------
520 # numsteps : int
521 #   Number of steps in the switch
522 # segresidnameList : list 
523 #   One or more "<segid>:<resid>" specifications - this is the same syntax as 
524 #   for the regular psfgen patch command.
525 #
526 # Returns:
527 # --------
528 # accept : boolean
529 #   Result of MC accept/reject test 
530 #
531 proc ::namdcph::runSwitch {numsteps segresidnameList} {
532     variable ::namdcph::useHeat
533
534     # (1) Checkpoint and modify output parameters. 
535     #
536     storeEnergies
537     set storedOutputEnergies [outputEnergies]
538     set storedDCDFreq [dcdFreq]
539     set storedTimestep $::energyArray(TS)
540     outputEnergies $numsteps
541     dcdFreq 0
542     firstTimestep 0
543     if {$::barostatIsSet && [$::barostatCmd]} {
544         $::barostatCmd off
545     }
546     # (2) Build the alchemical switch inputs.
547     #
548     alchemify $segresidnameList
549     # (3) Run the switch trajectory with momentum reversal.
550     #
551     runprpswitch $numsteps
552     outputEnergies $storedOutputEnergies
553     dcdFreq $storedDCDFreq
554     # (4) Compute the work with state dependent energy shifts.
555     #
556     storeEnergies
557     set DeltaE [cphSystem compute switch $segresidnameList]
558     if {$useHeat} {
559         # Computed as totalEnergy - totalEnergy0 - heat.
560         set Work $::energyArray(WORK)
561     } else {
562         # Computed as alchemical protocol work.
563         set Work $::energyArray(CUMALCHWORK)
564     }
565     set ReducedWork [expr {($Work + $DeltaE) / [::kBT]}]
566     set accept [metropolisAcceptance $ReducedWork]
567     set tmp [printProposalSummary $segresidnameList]
568     cphPrint [format "%s WorkCorr % 10.4f CorrWork % 10.4f"\
569             [join $tmp "/"] $DeltaE $Work]
570     # (5) Reinitialize for the next MD step based on accept/reject.
571     #
572     outputEnergies $storedOutputEnergies
573     dcdFreq $storedDCDFreq
574     firstTimestep $storedTimestep
575     if {$::barostatIsSet && ![$::barostatCmd]} {
576         $::barostatCmd on
577     }
578
579     if {$accept} {
580         cphPrint "Switch accepted!"
581         dealchemify $segresidnameList
582     } else {
583         cphPrint "Switch rejected!"
584         alch off
585         reloadAndReinit [getMDBasename] false
586     }
587     return $accept
588 }
589
590 # ::namdcph::alchemify
591 #
592 # Reinitialize the system with alchemical sidechains for the given residues.
593 #
594 # This includes multiple steps:
595 # (1) The NAMD state is written to disk and re-read by PSFGEN.
596 # (2) Appropriate alchemical patches are applied and alchemical atom 
597 #     coordinates and velocities are sampled.
598 # (3) New NAMD inputs are written and then re-read. This includes new 
599 #     extraBonds for the alchemical sidechain.
600 #
601 # Arguments:
602 # ----------
603 # segresidnameList : list of strings
604 #   One or more "<segid>:<resid>" specifications - this is the same syntax as 
605 #   for the regular psfgen patch command.
606 #
607 # Returns:
608 # --------
609 # None
610 #
611 proc ::namdcph::alchemify {segresidnameList} {  
612     variable ::namdcph::alchFrcCons
613     set T [$::thermostatTempCmd]
614
615     alch on
616     # (1) Read in the nonalchemical PSF and apply patches.
617     #
618     output [getMDBasename]
619     psfgenRead [getMDBasename]
620     # (2) Apply patches and build coordinates and velocities.
621     #
622     foreach segresidname $segresidnameList {
623         cphSystem alchemifypsf $segresidname $alchFrcCons $T
624     }
625     regenerate angles dihedrals
626     # (3) Write a new set of inputs and reinitialize.
627     #
628     psfgenWrite [getSWBasename] [mdFilename xsc]
629     # Dummy atoms have been built and a PDB has been written. We can now query
630     # atom indices and build extraBonds. If psfgenWrite were _not_ called, then
631     # the atomid queries would all return zero (that would be bad).
632     #
633     set ExtraBondsFile [open [swFilename extrabonds] "w"]
634     foreach segresidname $segresidnameList {
635         puts $ExtraBondsFile\
636                 [cphSystem get alchBonds $segresidname $alchFrcCons]
637     }
638     close $ExtraBondsFile
639     reloadAndReinit [getSWBasename] true
640     return
641 }
642
643 # ::namdcph::dealchemify
644 #
645 # Remove alchemical side
646 #
647 # Arguments:
648 # ----------
649 # segresidnameList : list of strings
650 #   One or more "<segid>:<resid>" specifications - this is the same syntax as 
651 #   for the regular psfgen patch command.
652 #
653 # Returns:
654 # --------
655 # None
656 #
657 proc ::namdcph::dealchemify {segresidnameList} {
658     output [getSWBasename]
659     psfgenRead [getSWBasename]
660     foreach segresidname $segresidnameList {
661         cphSystem dealchemifypsf $segresidname
662     }
663     psfgenWrite [getMDBasename] [swFilename xsc]
664     alch off
665     reloadAndReinit [getMDBasename] false
666     return
667 }
668
669 # =============================================================================
670 # Constant pH specific I/O
671 # =============================================================================
672 # ::namdcph::readConfig
673 #
674 # Read one or more constant pH config files and return template definitions.
675 #
676 # Arguments:
677 # ----------
678 # None
679 #
680 # Returns:
681 # --------
682 # templateDefs : dict
683 #   A dict conversion of the config file contents from JSON
684 #
685 proc ::namdcph::readConfig {} {
686     variable ::namdcph::configFilenames
687
688     set templateDefs [dict create]
689     foreach configFilename $configFilenames {
690         set ConfigFile [open $configFilename "r"]
691         set tmpDict [json::json2dict [read $ConfigFile]]
692         # Warn the user about re-definitions, but proceed anyway. Note that
693         # dict merge overwrites values right to left.
694         set oldKeys [dict keys $templateDefs]
695         foreach newKey [dict keys $tmpDict] {
696             if {$newKey in $oldKeys} {
697                 cphWarn "Reading and using new definition for residue $newKey"
698             }
699         } 
700         set templateDefs [dict merge $templateDefs $tmpDict] 
701         close $ConfigFile
702     }
703     return $templateDefs
704 }
705
706 # ::namdcph::readRestart
707 #
708 # Read in a constant pH state from file and _merge_ the data into the namespace
709 # variables. That is, give precedence to keyword specifications.
710 #
711 # NB! This also requires that cphSystem has been built so that the system and
712 #     state info can be compared.
713 #
714 proc ::namdcph::readRestart {} {
715     variable ::namdcph::restartFilename
716     variable ::namdcph::SystempH
717     variable ::namdcph::moveInfo
718
719     set RestartFile [open $restartFilename "r"]
720     set Restart [json::json2dict [read $RestartFile]]
721     close $RestartFile
722
723     if {![info exists SystempH] && [dict exists $Restart pH]} {
724         pH [dict get $Restart pH]
725     }
726
727     if {[dict exists $Restart exclude]} {
728         cphExcludeResidue {*}[dict get $Restart exclude]
729     }
730
731     # Read state parameters, if present. 
732     #
733     set stateList {}
734     if {[dict exists $Restart states]} {
735         set stateList [dict get $Restart states]
736     }
737     set pKaiList {}
738     if {[dict exists $Restart pKais]} {
739         set pKaiList [dict get $Restart pKais]
740     }
741
742     # Read MC move parameters, if present. Only reset those values which have
743     # not been set in via keywords.
744     #
745     if {[dict exists $Restart MCmoves]} {
746         set rstrtMoveInfo [dict get $Restart MCmoves]    
747     }
748     if {[info exists rstrtMoveInfo]} {
749         # This is a bit kludgy - the global data is not a nested dict, so it
750         # will crash in the dict for loop. Treat these specially and then unset
751         # them...
752         if {[dict exists $rstrtMoveInfo maxProposalAttempts]} {
753             if {![dict exists $moveInfo maxProposalAttempts]} {
754                 dict set moveInfo maxProposalAttempts\
755                     [dict get $rstrtMoveInfo maxProposalAttempts]
756             }
757             dict unset rstrtMoveInfo maxProposalAttempts
758         }
759         dict for {moveLabel data} $rstrtMoveInfo {
760             dict for {key value} $data {
761                  if {![dict exists $moveInfo $moveLabel $key]} {
762                      dict set moveInfo $moveLabel $key $value
763                  }
764             }
765         }
766     }
767
768     return [list $stateList $pKaiList]
769 }
770
771 # ::namdcph::writeRestart
772 #
773 # Write the current constant pH state information to restartFilename based on 
774 # the restart frequency. If the "force" keyword precedes the arguments, ignore 
775 # the restart frequency and write no matter what.
776 #
777 # Arguments:
778 # ----------
779 # restartFilename : string
780 #   Name of (JSON format) restart file to write
781 # cycle : int
782 #   Last neMD/MC cycle index
783 #
784 # Returns:
785 # --------
786 # None
787 #
788 proc ::namdcph::writeRestart {args} { 
789     variable ::namdcph::restartFreq
790     variable ::namdcph::SystempH
791     variable ::namdcph::excludeList
792
793     if {[string match [lindex $args 0] force]} {
794         set restartFilename [lindex $args 1]
795         set cycle [expr {int([lindex $args 2])}]
796     } else {
797         set restartFilename [lindex $args 0]
798         set cycle [expr {int([lindex $args 1])}]
799         if {!($restartFreq) || !($cycle) || [expr {$cycle % $restartFreq}]} {
800             return
801         }
802     }
803
804     # Here we write the restart as a dict in JSON format, however, this is NOT
805     # the same as converting a Tcl dict to JSON, because knowledge of types is
806     # lost in the conversion. The easiest (and maybe faster?) way is to build a
807     # list that contains all of the dict key/value pairs and then emulate the
808     # JSON format with commas and curly braces.
809     #
810     set rstrtList [list]
811     lappend rstrtList "\"cycle\":$cycle"
812     lappend rstrtList "\"pH\":$SystempH"
813     if {[llength $excludeList]} {
814         lappend rstrtList "\"exclude\":[json::list2json $excludeList 1]"
815     }
816     lappend rstrtList "\"states\":[json::list2json [cphSystem get state] 1]"
817     lappend rstrtList "\"pKais\":[json::list2json [cphSystem get pKai]]"
818     lappend rstrtList "\"MCmoves\":[json::dict2json [cphTitrator get archive]]"
819
820     # Write everything to file.
821     namdFileBackup $restartFilename
822     set RestartFile [open $restartFilename "w"]
823     puts $RestartFile "\{[join $rstrtList ,]\}"
824     close $RestartFile
825
826     # Always checkpoint the PSF and PDB when a restart is written.
827     file copy -force [mdFilename psf] "[outputName].psf"
828     file copy -force [mdFilename pdb] "[outputName].pdb"
829     return
830 }
831
832 # ::namdcph::openCpHLog
833 #
834 # Open a new constant pH log for proton occupancies and return the file object.
835 #
836 proc ::namdcph::openCpHLog {} {
837     variable ::namdcph::outFile
838     variable ::namdcph::SystempH
839
840     if {[string length $outFile] > 0} {
841         set logFilename $outFile
842     } else {
843         set logFilename "[outputname].cphlog"
844     }
845     namdFileBackup $logFilename
846     set cphlog [open $logFilename "w"]
847     puts $cphlog "#pH $SystempH"
848     puts $cphlog "#[join [cphSystem get segresidnames] " "]"
849     return $cphlog 
850 }
851
852 # =============================================================================
853 # Convenience Routines 
854 # =============================================================================
855 # ::namdcph::swFilename
856 #
857 # Get an appropriate filename for a temporary file used during switches
858 #
859 proc ::namdcph::swFilename {ext} {
860     return "[getSWBasename].$ext"
861 }
862
863 # ::namdcph::mdFilename
864 #
865 # Get an appropriate filename for a temporary file used to launch regular MD
866 #
867 proc ::namdcph::mdFilename {ext} {
868     return "[getMDBasename].$ext"
869 }
870
871 # ::namdcph::clearExtraBonds
872 #  
873 # Clear all bonds in the extraBondsFile.
874 #
875 proc ::namdcph::clearExtraBonds {} {
876     set ExtraBondsFilename [swFilename extrabonds]
877     cphPrint "clearing extraBondsFile $ExtraBondsFilename"
878     set ExtraBondsFile [open $ExtraBondsFilename "w"]
879     puts $ExtraBondsFile ""
880     close $ExtraBondsFile
881     return
882 }
883
884 # ::namdcph::reloadAndReinit
885 #
886 # Read in a new PSF/PDB pair and then reinitialize atoms using a basename read.
887 #
888 proc ::namdcph::reloadAndReinit {basename keepExtraBonds} {
889     if {!$keepExtraBonds} {
890         clearExtraBonds
891     }
892     structure "$basename.psf" pdb "$basename.pdb"
893     reinitatoms $basename
894     return
895 }
896
897 proc ::namdcph::cphPrint {msg} {
898     print "$::namdcph::TITLE $msg"
899 }
900
901 proc ::namdcph::cphWarn {msg} {
902     print "$::namdcph::TITLE WARNING! $msg"
903 }
904
905 proc ::namdcph::cphAbort {{msg ""}} {
906     abort "$::namdcph::TITLE $msg"
907 }
908
909 proc ::namdcph::cphWarnExpt {} {
910     cphWarn "THIS FEATURE IS EXPERIMENTAL!"
911     cphPrint "RESULTS ARE NOT GUARANTEEED - USE AT YOUR OWN RISK"
912 }
913
914 # =============================================================================
915 # Setup Routines
916 # =============================================================================
917 # ::namdcph::initialize
918 #
919 # Initialize the system for constant pH. This requires two main things to
920 # happen:
921 #
922 # 1) nonequilibrium alchemical transformations must be enabled
923 # 2) the PSF/PDB must be rebuilt to include dummy atoms and possibly reassigned
924 #    protonation states
925 #
926 proc ::namdcph::initialize {} {
927     variable ::namdcph::restartFilename
928     variable ::namdcph::configFilenames
929     variable ::namdcph::stateInfo
930     variable ::namdcph::moveInfo
931     variable ::namdcph::excludeList
932     variable ::namdcph::SystempH
933     variable ::namdcph::useHeat
934
935     callback energyCallback
936     # 1) Set up alchemical keywords and run startup. 
937     #
938     initializeAlch
939     getThermostat
940     getBarostat
941     # Note that checking the stored name avoids the need to catch in case we
942     # have an older version of NAMD where stochRescale is unavailable.
943     if {[string match $::thermostatName "stochastic-rescaling"]} {
944         if {[isset stochRescaleHeat] && [stochRescaleHeat]} {
945             set useHeat 1
946         }
947     }
948
949     # 2) Rebuild the PSF with dummy protons and modify protonation states as 
950     # needed. Build the residue definitions, assign states to each residue, and 
951     # rebuild the topology to reflect those states.
952     #
953     cphPrint "initializing constant pH PSF..."
954     if {![llength configFilenames]} {
955         cphAbort "At least one constant pH configuration file is required."
956     }
957
958     if {[string length $restartFilename]} {
959         lassign [readRestart] stateList pKaiList
960         lassign [list 0.0 false] temp buildH
961     } else {
962         lassign [list [$::thermostatTempCmd] false] temp buildH
963     }
964     cphSystem build [readConfig] $excludeList
965
966     # Now that we've built the system, we can allocate the state and pKa
967     # information from the restart (if present).
968     if {[info exists stateList] && [llength $stateList] &&
969         [info exists pKaiList] && [llength $pKaiList]} {
970         set segresidnameList [cphSystem get segresidnames]
971
972         if {[llength $stateList] != [llength $pKaiList]} {
973             cphAbort "mismatch in states/pKais in $restartFilename"
974         }
975         if {[llength $stateList] != [llength $segresidnameList]} {
976             cphAbort "Too few/many state/pKai definitions in $restartFilename"
977         }
978
979         foreach segresidname $segresidnameList\
980                 state $stateList pKai $pKaiList {
981             if {![dict exists $stateInfo $segresidname]} {
982                 dict set stateInfo $segresidname state $state
983                 dict set stateInfo $segresidname pKai $pKai
984             } else {
985                 if {![dict exists $stateInfo $segresidname state]} {
986                     dict set stateInfo $segresidname state $state
987                 }
988                 if {![dict exists $stateInfo $segresidname pKai]} {
989                     dict set stateInfo $segresidname pKai $pKai
990                 }
991             }
992         }
993     }
994
995     if {![info exists SystempH]} {
996         cphAbort "A pH value is required."
997     }
998
999     # All residue state info from all sources is now in stateInfo.
1000     # Check that all of the keys are valid.
1001     foreach segresidname [dict keys $stateInfo] {
1002         if {[cphSystem validate $segresidname]} {
1003             cphAbort
1004         }
1005     }
1006     cphSystem initialize $SystempH $temp $buildH $stateInfo
1007
1008     # 3) Build the MC move set (the "titrator").
1009     #
1010     # All neMD/MC move info from all sources is now in moveInfo.
1011     cphTitrator build $moveInfo
1012
1013     # 4) Report to stdout.
1014     printSettingsSummary
1015     printSystemSummary
1016     printTitratorSummary
1017     return
1018 }
1019
1020 proc ::namdcph::finalize {} {
1021     # 5) Write to disk and prepare for MD.
1022     if {[isset extendedSystem]} {
1023         set inputXSCName [extendedSystem]
1024     } else {
1025         set inputXSCName [mdFilename xsc]
1026         set inputXSC [open $inputXSCName "w"]
1027         puts $inputXSC $::dummyXSC 
1028         close $inputXSC
1029     }
1030     psfgenWrite [getMDBasename] $inputXSCName 
1031     reloadAndReinit [getMDBasename] false
1032     if {![isset binVelocities] && ![isset velocities]} {
1033         reinitvels [temperature]
1034     }
1035     return
1036 }
1037
1038 # ::namdcph::initializeAlch
1039 #
1040 #   Initialize the alchemical settings needed for nonequilibrium alchemical
1041 # switching of protonation states. This activates the appropriate nonbonded
1042 # kernels and stores SimParameters settings that will be needed during 
1043 # switches.
1044 #
1045 proc ::namdcph::initializeAlch {} {
1046     if {[isset alch]} {
1047         cphAbort "Constant pH is currently incompatible with alchemical"\
1048                  "transformations. Remove all alch settings."
1049     }
1050     cphPrint "Setting up nonequilibrium alchemical switching."
1051     # Here we require alchemical force computations and complete decoupling of
1052     # _all_ interactions. Unstaged linear coupling without softcore shifting is
1053     # currently assumed/enforced, but this might be ok to change.
1054     #
1055     alch on
1056     alchType TI
1057     alchDecouple off
1058     alchElecLambdaStart 0.0
1059     alchVdwLambdaEnd 1.0
1060     alchVdwShiftCoeff 0.0
1061     alchBondLambdaEnd 1.0
1062     alchBondDecouple on
1063     alchLambda 0.0
1064     alchLambda2 1.0
1065     alchOutFreq 0 ;# Suppress output - this would just be a nightmare.
1066     # Bonds between corresponding alchemical atoms are built dynamically when
1067     # a switch is started. In principle, this does not conflict with any other
1068     # extraBonds settings. However, indices in the extraBondsFile are likely
1069     # corrupted and would lead to wildly unexpected behavior - disable this for
1070     # now.
1071     #
1072     # TODO: Read and modify these files on they fly?
1073     if {[isset extraBonds]} {
1074         cphAbort "Constant pH is currently incompatible with extraBonds."
1075     }
1076     extraBonds on
1077     extraBondsFile [swFilename extrabonds]
1078     clearExtraBonds
1079     startup
1080     alchLambdaFreq 1
1081     alch off
1082     return
1083 }
1084
1085 # ::namdcph::printSettingsSummary
1086 #
1087 # Print a summary of the constant pH settings.
1088 #
1089 proc ::namdcph::printSettingsSummary {} {
1090     variable ::namdcph::configFilenames
1091     variable ::namdcph::restartFilename
1092     variable ::namdcph::SystempH
1093     variable ::namdcph::alchFrcCons
1094     variable ::namdcph::useHeat
1095
1096     set StarBar "***************************************"
1097     cphPrint $StarBar
1098     cphPrint "CONSTANT pH MD ACTIVE"
1099     if {[string length $restartFilename] > 0} {
1100         cphPrint "RESTART FILENAME $restartFilename"
1101     }
1102     cphPrint "SYSTEM pH $SystempH"
1103     cphPrint "CONSTANT pH CONFIGURATION FILE(S)"
1104     foreach configFilename $configFilenames {
1105         cphPrint "$configFilename"
1106     }
1107     cphPrint "NONEQUILIBRIUM SWITCH PARAMETERS:"
1108     cphPrint "ALCHEMICAL FORCE CONSTANT $alchFrcCons kcal/mol-A^2"
1109     cphPrint "neMD/MC CRITERION TEMPERATURE [$::thermostatTempCmd]"
1110     if {$useHeat} {
1111         cphPrint "WORK COMPUTATION: HEAT EXCHANGE W/THERMOSTAT"
1112         cphWarnExpt
1113     } else {
1114         cphPrint "WORK COMPUTATION: ALCHEMICAL PROTOCOL WORK"
1115     }
1116     cphPrint "TEMPORARY FILE INFORMATION:"
1117     cphPrint "cpH TOPOLOGY FILE BASENAME [getMDBasename]"
1118     cphPrint "neMD/MC TRAJECTORY BASENAME [getSWBasename]"
1119     cphPrint $StarBar
1120     return
1121 }
1122
1123 # ::namdcph::printSystemSummary
1124 #
1125 # Print a summary of the titratable system. 
1126 #
1127 proc ::namdcph::printSystemSummary {} {
1128     set StarBar "***************************************"
1129     cphPrint $StarBar
1130     cphPrint "TITRATABLE RESIDUE DEFINITIONS:"
1131     cphPrint "[join [cphSystem get resdefs] " "]"
1132     cphPrint "TITRATABLE SYSTEM SUMMARY:"
1133     cphPrint "[cphSystem get numresidues] RESIDUE(S)"
1134     cphPrint "[llength [cphSystem get occupancy]] PROTONATION SITE(S)"
1135     cphPrint $StarBar
1136     cphPrint [format "%-19s : %5s : %s" segid:resid:resname state pKai]
1137     foreach segresidname [cphSystem get segresidnames]\
1138             state [cphSystem get state]\
1139             pKaiList [cphSystem get pKai] {
1140         cphPrint [format "%-19s : % 5s : %-s" $segresidname $state $pKaiList]
1141     }
1142     cphPrint $StarBar
1143     return
1144 }
1145
1146 # ::namdcph::printTitratorSummary
1147 #
1148 # Print a summary of the MC moves.
1149 #
1150 proc ::namdcph::printTitratorSummary {} {
1151     set StarBar "*************************************************"
1152     set moveLabels [cphTitrator get moveLabels]
1153     set numMoves [llength $moveLabels]
1154
1155     cphPrint $StarBar
1156     cphPrint "CONSTANT pH neMD/MC MOVES:"
1157     cphPrint "MAX. ATTEMPTS PER CYCLE: [cphTitrator get maxAttempts]"
1158
1159     cphPrint "ONE RESIDUE MOVES"
1160     cphPrint [format "%-14s : %8s %8s" "move label" numsteps weight]
1161     foreach moveLabel $moveLabels {
1162         if {[llength [split $moveLabel "/"]] != 1} continue
1163         incr numEntries
1164         set numsteps [cphTitrator get numsteps $moveLabel]
1165         set weight [cphTitrator get weight $moveLabel]
1166         cphPrint [format "%-14s : % 8d % 8.2f" $moveLabel $numsteps $weight]
1167     }
1168     if {$numEntries == $numMoves} {
1169         cphPrint $StarBar
1170         return
1171     }
1172
1173     cphPrint "TWO RESIDUE MOVES : proton transfer (PT) : co-titration (CT)"
1174     cphPrint [format "%-29s : %8s %8s %2s %2s"\
1175             "move label" numsteps weight PT CT]
1176     foreach moveLabel $moveLabels {
1177         if {[llength [split $moveLabel "/"]] != 2} continue
1178         set numsteps [cphTitrator get numsteps $moveLabel]
1179         set weight [cphTitrator get weight $moveLabel]
1180         set PT [cphTitrator get protonTransfer $moveLabel]
1181         set CT [cphTitrator get cotitration $moveLabel]
1182         cphPrint [format "%-29s : % 8d % 8.2f %2d %2d"\
1183                 $moveLabel $numsteps $weight $PT $CT]
1184     }
1185     cphPrint $StarBar
1186     return
1187 }
1188
1189 # ::namdcph::printProposalSummary
1190 proc ::namdcph::printProposalSummary {segresidnameList} {
1191     set retList [list]
1192     foreach segresidname $segresidnameList {
1193         set state [cphSystem get state $segresidname]
1194         set trialState [cphSystem get trialState $segresidname]
1195         lappend retList "$segresidname:$state:$trialState"
1196     }
1197     return $retList
1198 }
1199
1200 # ::namdcph::printTitratorReport
1201 #
1202 # Print a report of the titration MC statistics.
1203 #
1204 proc ::namdcph::printTitratorReport {} {
1205     set StarBar "*************************************************"
1206     cphPrint $StarBar
1207     cphPrint "CONSTANT pH MD STATISTICS:"
1208     cphPrint [format "%-25s : %-8s %-12s" "move label" attempts "accept. rate"]
1209     foreach moveLabel [cphTitrator get moveLabels] {
1210         set attempts [cphTitrator get attempts $moveLabel]
1211         set successes [cphTitrator get successes $moveLabel]
1212         if {$successes && $attempts} {
1213             set rate [expr {100.*$successes/$attempts}]
1214         } else {
1215             set rate 0.0
1216         }
1217         cphPrint [format "%-25s : %8d %12.2f" $moveLabel $attempts $rate] 
1218     }
1219     cphPrint $StarBar
1220     return
1221 }
1222
1223 # =============================================================================
1224 # Getter Routines
1225 #
1226 # These are largely unnecessary, but cut down on "variable" declarations.
1227 # =============================================================================
1228 proc ::namdcph::getMDBasename {} {
1229     return $::namdcph::MDBasename
1230 }
1231
1232 proc ::namdcph::getSWBasename {} {
1233     return $::namdcph::SWBasename
1234 }