Bug fixes in Monte Carlo functionality
[namd.git] / lib / namdcph / namdcph / cphtitrator.tcl
1 # cphtitrator.tcl
2 #
3 #   This file provides the ::cphTitrator namespace, which effectively emulates
4 # an object containing all constant-pH information and routines relevant to
5 # proposing and performing Monte Carlo titrations.
6 #
7 #   If cphSystem and cphTitrator were true objects, this would essentially be
8 # an object composition. That is, a cphTitrator "contains" a cphSystem and
9 # yields information pertaining to it depending on a context (e.g. the pH).
10 #
11 package require Tcl 8.5
12
13 source [file join [file dirname [info script]] "cphtoppar.tcl"]
14 source [file join [file dirname [info script]] "numtcl.tcl"]
15 source [file join [file dirname [info script]] "namdmcmc.tcl"]
16
17 namespace eval ::cphTitrator {
18     #   The core information of the cphTitrator is the moveDict, which stores
19     # all possible neMD/MC moves involving one or more titratable residues. The
20     # keys are of the form:
21     #
22     # segresidname[/segresidname..[/segresidname]]
23     #
24     # in analogy to the keys of ::cphSystem::residueDict. Note that a "/" is
25     # used as a syntactic divider between residues. Only _specific_ move
26     # information is kept in moveDict. Another dict, defaultMoveParams is used
27     # to store default information (similar to ::cphSystem::residueDefDict) and
28     # is used to avoid saving copies of redundant information.
29     #
30     namespace import ::cphSystem::*
31
32     variable moveDict [dict create]
33     variable defaultMoveParams [dict create]
34     variable maxAttempts
35
36     namespace export cphSystem
37     namespace export cphTitrator
38 }
39
40 # =============================================================================
41 # cphTitrator Interface function
42 # =============================================================================
43 # ::cphTitrator::cphTitrator
44 #
45 # This is the only exported function from the cphTitrator namespace and
46 # provides a complete interface. The first argument is always an action (e.g.
47 # "get" or "set") and the rest are variable arguments.
48 #
49 # action           description
50 # ---------------- -----------
51 # get              see cphTitratorGet 
52 # set              see cphTitratorSet
53 # build            see buildTitrator
54 # propose          select and propose a move, return switch information
55 # accumulateStats  accumulate mean acceptance rate for the given move
56 #
57 # Some notes about Tcl switch statements for non-guru types like myself:
58 #
59 #   These can appear a bit delicate to those who are not familiar with their
60 # quirks (compared to C). For one, badly placed comments can break them. Also
61 # note that a switch returns the expression it first matches and does NOT use
62 # break statements.
63 #
64 proc ::cphTitrator::cphTitrator {action args} {
65     return [switch -nocase -- $action {
66         get {
67             cphTitratorGet {*}$args
68         }
69         set {
70             cphTitratorSet {*}$args
71         }
72         build {
73             buildTitrator {*}$args 
74         }
75         propose { ;# sample the move set and return the selected move info.
76             proposeMove {*}$args
77         }
78         accumulateStats { ;# accumulate MC statistics
79             accumulateAcceptanceRate {*}$args
80         }
81         default {
82             abort "Invalid cphTitrator action $action"
83         }
84     }]
85 }
86
87 # ::cphTitrator::validateMoveLabel
88 #
89 # Check that a given move label is valid. That is:
90 #
91 # 1) is it in the correct segresidname[/segresidname ...] format?
92 # 2) are the component segresidnames valid?
93 #
94 # A non-zero error code and msg are returned for each condition. Zero is
95 # returned for a valid segresidname.
96 #
97 # NB! This should only be called _after_ buildSystem.
98 #
99 proc ::cphTitrator::validateMoveLabel {moveLabel} {
100     set segresidnameList [split $moveLabel "/"]
101     if {![llength $segresidnameList]} {
102         print "Move labels should be of the form segresidname/segresidname ...]"
103         return -1
104     }
105     foreach segresidname $segresidnameList { 
106         if {[cphSystem validate $segresidname]} {
107             print "Bad segresidname $segresidname in move label $moveLabel!"
108             return -2
109         }
110     }
111     return 0
112 }
113
114 # Pop the value from a dict with the corresponding key. If no key exists,
115 # return the the default value instead.
116 #
117 proc ::cphTitrator::dictPopOrDefault {myDict key {defaultValue {}}} {
118     upvar 1 $myDict MyDict
119     if {[dict exists $MyDict $key]} {
120         set value [dict get $MyDict $key]
121         dict unset MyDict $key
122     }
123     if {![info exists value]} {
124         set value $defaultValue
125     }
126     return $value
127 }
128
129 # =============================================================================
130 # Proposal Routines
131 # =============================================================================
132 # ::cphTitrator::proposeMove
133 #
134 # Propose a move from all possible moves in the set. This is done by directly
135 # sampling the probability mass function of the weighted moves. Once this is
136 # done, test the inherent probability of that move given the imposed pH. If
137 # rejected, try again up to the given number of "maxAttempts". Always return
138 # the information necessary to perform a switch for the given proposal.
139 #
140 proc ::cphTitrator::proposeMove {pH} {
141     set weights [cphTitrator get weight]
142     lassign [choice [cphTitrator get moveLabels] $weights] moveLabel
143     set proposalCmd [cphTitrator get proposalCmd $moveLabel]
144     set accept [eval $proposalCmd $pH]
145     set segresidnameList [cphTitrator get segresidnameList $moveLabel]
146     set numsteps [cphTitrator get numsteps $moveLabel]
147     return [list $accept $numsteps $segresidnameList]
148 }
149
150 # FOR TESTING ONLY - INDEPENDENCE SAMPLING
151 #
152 #proc ::cphTitrator::proposeResidueMove {segresidname pH} {
153 #    lassign [cphSystem compute inherentWeights $pH $segresidname] states weights
154 #    lassign [choice $states $weights] trialState
155 #    set currState [cphSystem get state $segresidname]
156 #    if {$trialState == $currState} {
157 #        return 0
158 #    } else {
159 #        cphSystem set trialState $segresidname $trialState
160 #        return 1
161 #    }
162 #}
163
164 # FOR TESTING ONLY - METROPOLIS MONTE CARLO (NON-ERGODIC IN SOME CASES!)
165 #
166 #proc ::cphTitrator::proposeResidueMove {segresidname pH} {
167 #    cphSystem propose titration $segresidname
168 #    set du [cphSystem compute inherent $pH $segresidname]
169 #    return [metropolisAcceptance $du]
170 #}
171
172 # ::cphTitrator::proposeResidueMove
173 #
174 # Propose a move involving a single residue via Metropolized independence
175 # sampling. Return True if such a move was proposed, accepted, and stored,
176 # otherwise return False.
177 #
178 proc ::cphTitrator::proposeResidueMove {segresidname pH} {
179     lassign [cphSystem compute inherentWeights $pH $segresidname] states weights
180     # Construct the Metropolized proposal weights for i-->j.
181     set equivStates [cphSystem get equivStateList $segresidname]
182     set proposalWeights [list]
183     set stateiProposalNorm 1.0
184     foreach state $states weight $weights {
185         if {$state in $equivStates} {
186             lappend proposalWeights 0.0
187             set stateiProposalNorm [expr {$stateiProposalNorm - $weight}]
188         } else {
189             lappend proposalWeights $weight
190         }
191     }
192     lassign [choice $states $proposalWeights] trialState
193     # Construct the Metropolized proposal weights for j-->i.
194     set equivStates [cphSystem get equivStateList $segresidname $trialState]
195     set proposalWeights [list]
196     set statejProposalNorm 1.0
197     foreach state $states weight $weights {
198         if {$state in $equivStates} {
199             lappend proposalWeights 0.0
200             set statejProposalNorm [expr {$statejProposalNorm - $weight}]
201         } else {
202             lappend proposalWeights $weight
203         }
204     }
205     # NB: Invert log argument bc of negative sign in Metropolis routine.
206     set du [expr {log($statejProposalNorm / $stateiProposalNorm)}]
207     set accept [metropolisAcceptance $du]
208     if {$accept} {
209         cphSystem set trialState $segresidname $trialState
210     }
211     return $accept
212 }
213
214 # ::cphTitrator::proposeTwoResidueMove
215 #
216 # Propose a move involving two residues. This may take two forms:
217 #
218 # 1) a proton transfer, in which a proton is moved from one resiue to another
219 #
220 # 2) a co-titration, in which two residues change concurrently
221 #
222 # Return True if either kind of move is proposed, accepted, and stored.
223 #
224 proc ::cphTitrator::proposeTwoResidueMove {moveLabel pH} {
225     lassign [split $moveLabel "/"] srn1 srn2 ;# srn = segresidname
226
227     set doPT [cphTitrator get protonTransfer $moveLabel]
228     set doCT [cphTitrator get cotitration $moveLabel]
229     if {$doPT} {
230         set doPT [expr {![cphSystem propose protonTransfer $srn1 $srn2]}]
231     }
232     # We don't need to bother with co-titration if proton transfer is already
233     # proposed, since that means it isn't even possible. The order of these
234     # could equally be switched, but it is generally faster to check that
235     # proton transfer is possible, so this is done first.
236     #
237     if {!$doPT && $doCT} {
238         set doCT [expr {![cphSystem propose cotitration $srn1 $srn2]}]
239     }
240     if {$doPT || $doCT} {
241         set du1 [cphSystem compute inherent $pH $srn1]
242         set du2 [cphSystem compute inherent $pH $srn2]
243         return [metropolisAcceptance [expr {$du1 + $du2}]]
244     }
245     return 0
246 }
247
248 # =============================================================================
249 # "Constructor" Routines 
250 # =============================================================================
251 # ::cphTitrator::buildTitrator
252 #
253 # Construct the titrator given a set of MC move information.
254 #
255 proc ::cphTitrator::buildTitrator {moveInfo} {
256     variable ::cphTitrator::moveDict
257     variable ::cphTitrator::defaultMoveParams
258     variable ::cphTitrator::maxAttempts
259
260     # 1) Pop off global and default move parameters.
261     #
262     set maxAttempts [dictPopOrDefault moveInfo maxProposalAttempts 0]
263     set maxAttempts 1
264     set defaultMoveParams [dict merge $defaultMoveParams\
265                                       [dictPopOrDefault moveInfo default]]
266     foreach attr {weight protonTransfer cotitration}\
267             hardcodedDefault {1.0 0 0} {
268         if {![dict exists $defaultMoveParams $attr]} {
269             cphTitrator set $attr default $hardcodedDefault
270         }
271     }
272
273     # 2) Validate all remaining keys - these better be move labels!
274     #
275     dict for {moveLabel data} $moveInfo {
276         if {[validateMoveLabel $moveLabel]} {
277             abort
278         }
279     }
280
281     # The default move set is to titrate each residue independently.
282     foreach segresidname [cphSystem get segresidnames] {
283         set moveLabel $segresidname
284         cphTitrator set proposalCmd $moveLabel\
285                 "proposeResidueMove $segresidname"
286         if {![dict exists $moveInfo $moveLabel]} continue
287
288         dict for {attr value} [dict get $moveInfo $moveLabel] {
289             cphTitrator set $attr $moveLabel $value
290         }
291         dict unset moveInfo $moveLabel
292     }
293     # Build proton transfer and co-titration moves, if present.
294     dict for {moveLabel data} $moveInfo {
295         set isTwoResidueMove 0
296         foreach type {protonTransfer cotitration} {
297             if {[dict exists $data $type] && [dict get $data $type]} {
298                 set isTwoResidueMove 1
299                 break
300             }
301         }
302         if {!$isTwoResidueMove} continue
303         if {[llength [split $moveLabel "/"]] != 2} {
304             abort "Proton transfer and co-titration moves must have exactly 2\
305                    segresidnames!"
306         }
307         cphTitrator set proposalCmd $moveLabel\
308                 "proposeTwoResidueMove $moveLabel"
309         dict for {attr value} $data {
310             cphTitrator set $attr $moveLabel $value
311         }
312         dict unset moveInfo $moveLabel
313     }
314     # Abort if extra invalid move information was provided.
315     if {[dict size $moveInfo]} {
316         print "ERROR! Invalid move specifications:"
317         dict for {moveLabel data} $moveInfo {
318             print "$moveLabel $data"
319         }
320         abort
321     }
322     # Initialize statistics for any moves that are new.
323     dict for {moveLabel data} $moveDict {
324         if {![dict exists $data attempts]} {
325             cphTitrator set successes $moveLabel 0
326             cphTitrator set attempts $moveLabel 0
327         }
328     }
329
330     return
331 }
332
333 # =============================================================================
334 # Getter Routines
335 # =============================================================================
336 # ::Titrator::cphTitratorGet
337 #
338 # Getter for move attributes, called as:
339 #
340 #   <attribute> [<moveLabel>]
341 #
342 # With some specialized exceptions, <attribute> is either the name of a system
343 # attribute (usually a list) or else a move attribute.
344 #
345 # system attributes  description
346 # -----------------  -----------
347 # moveLabels         list of all the moveLabels
348 # maxAttempts        max # of attempts at move proposals
349 #
350 # move attributes  description
351 # ---------------  -----------
352 # numsteps         number of steps in a switch after successful proposal
353 # weight           weight for consideration in move selection
354 # protonTransfer   if 2 residues are involved, can they proton transfer?
355 # cotitration      if 2 residues are involved, can they cotitrate?
356 # successes        the number of successful _switches_ for this move
357 # attempts         the number of attempted _switches_ for this move
358 # proposalCmd      the command needed to set up trial states
359 # segresidnameList list of the residues involved (NB: this is just shorthand
360 #                  for [split $moveLabel "/"])
361 #
362 # specialized calls  description
363 # -----------------  -----------
364 # nodefaults         Return a dictionary of all current move information, but
365 #                    suppress entries where a move has the default value. This
366 #                    gives a minimal, but full, description of the move set.
367 #
368 # For now this is _much_ simpler than the analogous cphSystemGet. In the future
369 # it may be useful to generalize this though.
370 #
371 proc ::cphTitrator::cphTitratorGet {attr {moveLabel {}}} {
372     variable ::cphTitrator::moveDict    
373     variable ::cphTitrator::maxAttempts
374
375     set getAll [expr {![llength $moveLabel]}]
376     if {!$getAll && ![dict exists $moveDict $moveLabel]
377         && ![string match -nocase $moveLabel default]} {
378         abort "cphTitratorGet: Invalid moveLabel $moveLabel"
379     }
380
381     return [switch -nocase -- $attr {
382         moveLabels {
383             dict keys $moveDict
384         }
385         maxAttempts {
386             expr {$maxAttempts}
387         }
388         numsteps -
389         weight -
390         protonTransfer -
391         cotitration -
392         successes -
393         attempts -
394         proposalCmd {
395             if {$getAll} {
396                 getAllMoveAttr $attr 
397             } else {
398                 getMoveAttr $attr $moveLabel
399             }
400         }
401         segresidnameList {
402             if {$getAll} {
403                 getAllSegresidnameLists
404             } else {
405                 split $moveLabel "/"
406             }
407         }
408         archive {
409             getMinimalMoveDict
410         }
411         default {
412             abort "cphTitratorGet: Invalid attribute $attr"
413         }
414     }]
415 }
416
417 proc ::cphTitrator::cannotGetAll {attr} {
418     abort "cphTitratorGet: Cannot get all $attr - must select a moveLabel"
419     return -1
420 }
421
422 # ---------------------------
423 # Getters for move attributes
424 # ---------------------------
425 # Return the default value if no specific value was set.
426 proc ::cphTitrator::getMoveAttr {attr moveLabel} {
427     variable ::cphTitrator::moveDict
428     variable ::cphTitrator::defaultMoveParams
429
430     if {[dict exists $moveDict $moveLabel $attr]} {
431         return [dict get $moveDict $moveLabel $attr]
432     } elseif {[dict exists $defaultMoveParams $attr]} {
433         return [dict get $defaultMoveParams $attr]
434     }
435     abort "cphTitratorGet: Error getting $attr for move $moveLabel"
436     return -1
437 }
438
439 proc ::cphTitrator::getAllMoveAttr {attr} {
440     set retList [list]
441     foreach moveLabel [cphTitrator get moveLabels] {
442         lappend retList [cphTitrator get $attr $moveLabel]
443     }
444     return $retList
445 }
446
447 proc ::cphTitrator::getAllSegresidnameLists {} {
448     set retList [list]
449     foreach moveLabel [cphTitrator get moveLabels] {
450         lappend retList [split $moveLabel "/"]
451     }
452     return $retList
453 }
454
455 # ----------------------------
456 # Special getter for archiving
457 # ----------------------------
458 # Return a minimal version of moveDict. This removes inferred values such as
459 # command names as well as specific settings that now match the defaults
460 # (since these can change after a restart).
461 #
462 proc ::cphTitrator::getMinimalMoveDict {} {
463     variable ::cphTitrator::defaultMoveParams
464
465     set minMoveDict [dict create]
466     dict set minMoveDict maxProposalAttempts [cphTitrator get maxAttempts]
467     dict set minMoveDict default $defaultMoveParams
468     foreach moveLabel [cphTitrator get moveLabels] {
469         foreach attr [dict keys $defaultMoveParams] {
470             set defaultValue [cphTitrator get $attr default]
471             set thisValue [cphTitrator get $attr $moveLabel]
472             if {$thisValue != $defaultValue} {
473                 dict set minMoveDict $moveLabel $attr $thisValue
474             }
475         }
476         dict set minMoveDict $moveLabel attempts\
477                 [cphTitrator get attempts $moveLabel]
478         dict set minMoveDict $moveLabel successes\
479                 [cphTitrator get successes $moveLabel]
480     }
481     return $minMoveDict
482 }
483
484 # =============================================================================
485 # Setter Routines
486 # =============================================================================
487 # ::cphTitrator::cphTitratorSet
488 #
489 # Setters for move attributes, called as:
490 #
491 #  <attribute> <moveLabel> <value>
492 #
493 # <attribute> is the name of a move attribute.
494 #
495 # move attributes  description
496 # ---------------  -----------
497 # numsteps         number of steps in a switch after successful proposal
498 # weight           weight for consideration in move selection
499 # protonTransfer   if 2 residues are involved, can they proton transfer?
500 # cotitration      if 2 residues are involved, can they cotitrate?
501 # successes        the number of successful _switches_ for this move
502 # attempts         the number of attempted _switches_ for this move
503 # proposalCmd      the command needed to set up trial states
504 #
505 proc ::cphTitrator::cphTitratorSet {attr moveLabel value} {
506     variable ::cphTitrator::moveDict
507     variable ::cphTitrator::defaultMoveParams
508
509     if {[string match -nocase $moveLabel default]} {
510         set setDefault 1
511     } else {
512         if {[validateMoveLabel $moveLabel]} {
513             abort
514         }
515         set setDefault 0
516     }
517
518     return [switch -nocase -- $attr {
519         numsteps -
520         successes -
521         attempts -
522         protonTransfer -
523         cotitration { ;# Require integer or boolean argument.
524             set value [expr {int($value)}]
525             if {$setDefault} {
526                 dict set defaultMoveParams $attr $value
527             } else {
528                 dict set moveDict $moveLabel $attr $value
529             }
530             expr {$value}
531         }
532         weight { ;# Require float argument.
533             set value [expr {1.*$value}]
534             if {$setDefault} {
535                 dict set defaultMoveParams $attr $value
536             } else {
537                 dict set moveDict $moveLabel $attr $value
538             }
539             expr {$value}
540         }
541         proposalCmd { ;# Argument is string or pure list.
542             dict set moveDict $moveLabel $attr $value
543             expr {$value}
544         }
545         default {
546             abort "cphTitratorSet: Invalid attribute $attr"
547             expr {-1}
548         }
549     }]
550 }
551
552 # =============================================================================
553 # Routines for tracking and reporting MC statistics
554 # =============================================================================
555 # ::cphTitrator::accumulateAcceptanceRate
556 #
557 # Accumulate statistics for the given move.
558 #
559 proc ::cphTitrator::accumulateAcceptanceRate {accept moveLabel} {
560     # Alas, dict incr does not support nested keys.
561     set successes [cphTitrator get successes $moveLabel]
562     set attempts [cphTitrator get attempts $moveLabel]
563     incr successes $accept
564     incr attempts
565     cphTitrator set successes $moveLabel $successes
566     cphTitrator set attempts $moveLabel $attempts
567     return
568 }
569