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