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