07b9784eea60619022ddaf944340e5ce9946c866
[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     variable ::cphTitrator::maxAttempts
142     set weights [cphTitrator get weight]
143
144     set accept 0
145     set attemptsThisCycle 0
146     while {!$accept && $attemptsThisCycle < $maxAttempts} {
147         incr attemptsThisCycle
148         # Clear the previous trial if it was rejected.
149         if {$attemptsThisCycle > 1} {
150              cphSystem update $accept $segresidnameList
151         }
152
153         lassign [choice [cphTitrator get moveLabels] $weights] moveLabel
154         set proposalCmd [cphTitrator get proposalCmd $moveLabel]
155         set accept [eval $proposalCmd $pH]
156         set segresidnameList [cphTitrator get segresidnameList $moveLabel]
157     }
158     set numsteps [cphTitrator get numsteps $moveLabel]
159     return [list $accept $numsteps $segresidnameList $attemptsThisCycle]
160 }
161
162 # ::cphTitrator::proposeResidueMove
163 #
164 # Propose a move involving a single residue via Metropolized independence
165 # sampling. Return True if such a move was proposed, accepted, and stored.
166 #
167 # The core idea here is that we want to propose the most probable new state
168 # based on the estimated inherent pKa values. For two state systems, this is
169 # trivial, it's just the other state. For three or more states, this will
170 # naturally prefer tautomerization at extreme pKa values.
171 #
172 proc ::cphTitrator::proposeResidueMove {segresidname pH} {
173    lassign [cphSystem compute inherentWeights $pH $segresidname] weights states
174
175    set pi [lindex $weights 0]
176    lassign [choice [lrange $weights 1 end] [lrange $weights 1 end]] pj j
177    set du [expr {log((1. - $pj) / (1. - $pi))}]
178    set accept [metropolisAcceptance $du]
179    if {$accept} {
180        cphSystem set trialState $segresidname [lindex $states $j]
181    }
182    return $accept
183 }
184
185 # ::cphTitrator::proposeProtonTransferMove
186 #
187 # Propose a move involving two residues in which a proton is moved from one to
188 # the other. Return True if such a move was proposed, accepted, and stored.
189 #
190 # Note that this automatically fails if no proton can be transferred.
191 #
192 proc ::cphTitrator::proposeProtonTransferMove {moveLabel pH} {
193     lassign [split $moveLabel "/"] segresidname1 segresidname2
194
195     if {[cphSystem propose protonTransfer $segresidname1 $segresidname2]} {
196         set accept 0
197     } else {
198         set du1 [cphSystem compute inherent $pH $segresidname1]
199         set du2 [cphSystem compute inherent $pH $segresidname2]
200         set accept [metropolisAcceptance [expr {$du1 + $du2}]]
201     }
202     return $accept
203 }
204
205 # =============================================================================
206 # "Constructor" Routines 
207 # =============================================================================
208 # ::cphTitrator::buildTitrator
209 #
210 # Construct the titrator given a set of MC move information.
211 #
212 proc ::cphTitrator::buildTitrator {moveInfo} {
213     variable ::cphTitrator::moveDict
214     variable ::cphTitrator::defaultMoveParams
215     variable ::cphTitrator::maxAttempts
216
217     # 1) Pop off global and default move parameters.
218     #
219     set maxAttempts [dictPopOrDefault moveInfo maxProposalAttempts 0]
220     if {$maxAttempts <= 0} {
221         # This is a reasonable default for most systems.
222         set maxAttempts [cphSystem get numresidues]
223     }
224     set defaultMoveParams [dict merge $defaultMoveParams\
225                                       [dictPopOrDefault moveInfo default]]
226     if {![dict exists $defaultMoveParams default weight]} {
227         cphTitrator set weight default 1.0
228     }
229     if {![dict exists $defaultMoveParams default protonTransfer]} {
230         cphTitrator set protonTransfer default 0
231     }
232
233     # 2) Validate all remaining keys - these better be move labels!
234     #
235     dict for {moveLabel data} $moveInfo {
236         if {[validateMoveLabel $moveLabel]} {
237             abort
238         }
239     }
240
241     # The default move set is to titrate each residue independently.
242     foreach segresidname [cphSystem get segresidnames] {
243         set moveLabel $segresidname
244         cphTitrator set proposalCmd $moveLabel\
245                 "proposeResidueMove $segresidname"
246         if {![dict exists $moveInfo $moveLabel]} continue
247         dict for {attr value} [dict get $moveInfo $moveLabel] {
248             cphTitrator set $attr $moveLabel $value
249         }
250     }
251     # Build proton transfer moves, if present.
252     dict for {moveLabel data} $moveInfo {
253         if {![dict exists $data protonTransfer]
254             || ![dict get $data protonTransfer]} {
255             continue
256         }
257         if {[llength [split $moveLabel "/"]] != 2} {
258             abort "Proton transfer moves must have exactly 2 segresidnames!"
259         }
260         cphTitrator set proposalCmd $moveLabel\
261                 "proposeProtonTransferMove $moveLabel"
262         if {![dict exists $moveInfo $moveLabel]} continue
263         dict for {attr value} [dict get $moveInfo $moveLabel] {
264             cphTitrator set $attr $moveLabel $value
265         }
266     }
267
268     # Initialize statistics for any moves that are new.
269     dict for {moveLabel data} $moveDict {
270         if {![dict exists $data attempts]} {
271             cphTitrator set successes $moveLabel 0
272             cphTitrator set attempts $moveLabel 0
273         }
274     }
275
276     return
277 }
278
279 # =============================================================================
280 # Getter Routines
281 # =============================================================================
282 # ::Titrator::cphTitratorGet
283 #
284 # Getter for move attributes, called as:
285 #
286 #   <attribute> [<moveLabel>]
287 #
288 # With some specialized exceptions, <attribute> is either the name of a system
289 # attribute (usually a list) or else a move attribute.
290 #
291 # system attributes  description
292 # -----------------  -----------
293 # moveLabels         list of all the moveLabels
294 # maxAttempts        max # of attempts at move proposals
295 #
296 # move attributes  description
297 # ---------------  -----------
298 # numsteps         number of steps in a switch after successful proposal
299 # weight           weight for consideration in move selection
300 # protonTransfer   if 2 residues are involved, can they proton transfer?
301 # successes        the number of successful _switches_ for this move
302 # attempts         the number of attempted _switches_ for this move
303 # proposalCmd      the command needed to set up trial states
304 # segresidnameList list of the residues involved (NB: this is just shorthand
305 #                  for [split $moveLabel "/"])
306 #
307 # specialized calls  description
308 # -----------------  -----------
309 # nodefaults         Return a dictionary of all current move information, but
310 #                    suppress entries where a move has the default value. This
311 #                    gives a minimal, but full, description of the move set.
312 #
313 # For now this is _much_ simpler than the analogous cphSystemGet. In the future
314 # it may be useful to generalize this though.
315 #
316 proc ::cphTitrator::cphTitratorGet {attr {moveLabel {}}} {
317     variable ::cphTitrator::moveDict    
318     variable ::cphTitrator::maxAttempts
319
320     set getAll [expr {![llength $moveLabel]}]
321     if {!$getAll && ![dict exists $moveDict $moveLabel]
322         && ![string match -nocase $moveLabel default]} {
323         abort "cphTitratorGet: Invalid moveLabel $moveLabel"
324     }
325
326     return [switch -nocase -- $attr {
327         moveLabels {
328             dict keys $moveDict
329         }
330         maxAttempts {
331             expr {$maxAttempts}
332         }
333         numsteps -
334         weight -
335         protonTransfer -
336         successes -
337         attempts -
338         proposalCmd {
339             if {$getAll} {
340                 getAllMoveAttr $attr 
341             } else {
342                 getMoveAttr $attr $moveLabel
343             }
344         }
345         segresidnameList {
346             if {$getAll} {
347                 getAllSegresidnameLists
348             } else {
349                 split $moveLabel "/"
350             }
351         }
352         nodefaults {
353             if {$getAll} {
354                 cannotGetAll $attr
355             } else {
356                 getMoveNoDefaults $moveLabel
357             }
358         }
359         default {
360             abort "cphTitratorGet: Invalid attribute $attr"
361         }
362     }]
363 }
364
365 proc ::cphTitrator::cannotGetAll {attr} {
366     abort "cphTitratorGet: Cannot get all $attr - must select a moveLabel"
367     return -1
368 }
369
370 # ---------------------------
371 # Getters for move attributes
372 # ---------------------------
373 # Return the default value if no specific value was set.
374 proc ::cphTitrator::getMoveAttr {attr moveLabel} {
375     variable ::cphTitrator::moveDict
376     variable ::cphTitrator::defaultMoveParams
377
378     if {[dict exists $moveDict $moveLabel $attr]} {
379         return [dict get $moveDict $moveLabel $attr]
380     } elseif {[dict exists $defaultMoveParams $attr]} {
381         return [dict get $defaultMoveParams $attr]
382     }
383     abort "cphTitratorGet: Error getting $attr for move $moveLabel"
384     return -1
385 }
386
387 proc ::cphTitrator::getAllMoveAttr {attr} {
388     set retList [list]
389     foreach moveLabel [cphTitrator get moveLabels] {
390         lappend retList [cphTitrator get $attr $moveLabel]
391     }
392     return $retList
393 }
394
395 proc ::cphTitrator::getAllSegresidnameLists {} {
396     set retList [list]
397     foreach moveLabel [cphTitrator get moveLabels] {
398         lappend retList [split $moveLabel "/"]
399     }
400     return $retList
401 }
402
403 # -----------------------------
404 # Special getters for archiving
405 # -----------------------------
406 # Return the dict for a given move label, but strip out default values.
407 proc ::cphTitrator::getMoveNoDefaults {moveLabel} {
408     set defaultNumsteps [cphTitrator get numsteps default]
409     set defaultWeight [cphTitrator get weight default]
410     set defaultPT [cphTitrator get protonTransfer default]
411     set numsteps [cphTitrator get numsteps $moveLabel]
412     set weight [cphTitrator get weight $moveLabel]
413     set PT [cphTitrator get protonTransfer $moveLabel]
414
415     set retDict [dict create]
416     if {$numsteps != $defaultNumsteps} {
417         dict set retDict numsteps $numsteps
418     }
419     if {$weight != $defaultWeight} {
420         dict set retDict weight $weight
421     }
422     if {$PT != $defaultPT} {
423         dict set retDict protonTransfer $PT
424     }
425     dict set retDict attempts [cphTitrator get attempts $moveLabel]
426     dict set retDict successes [cphTitrator get successes $moveLabel]
427
428     return $retDict
429 }
430
431 # =============================================================================
432 # Setter Routines
433 # =============================================================================
434 # ::cphTitrator::cphTitratorSet
435 #
436 # Setters for move attributes, called as:
437 #
438 #  <attribute> <moveLabel> <value>
439 #
440 # <attribute> is the name of a move attribute.
441 #
442 # move attributes  description
443 # ---------------  -----------
444 # numsteps         number of steps in a switch after successful proposal
445 # weight           weight for consideration in move selection
446 # protonTransfer   if 2 residues are involved, can they proton transfer?
447 # successes        the number of successful _switches_ for this move
448 # attempts         the number of attempted _switches_ for this move
449 # proposalCmd      the command needed to set up trial states
450 #
451 proc ::cphTitrator::cphTitratorSet {attr moveLabel value} {
452     variable ::cphTitrator::moveDict
453     variable ::cphTitrator::defaultMoveParams
454
455     if {[string match -nocase $moveLabel default]} {
456         set setDefault 1
457     } else {
458         if {[validateMoveLabel $moveLabel]} {
459             abort
460         }
461         set setDefault 0
462     }
463
464     return [switch -nocase -- $attr {
465         numsteps -
466         successes -
467         attempts -
468         protonTransfer { ;# Require integer or boolean argument.
469             set value [expr {int($value)}]
470             if {$setDefault} {
471                 dict set defaultMoveParams $attr $value
472             } else {
473                 dict set moveDict $moveLabel $attr $value
474             }
475             expr {$value}
476         }
477         weight { ;# Require float argument.
478             set value [expr {1.*$value}]
479             if {$setDefault} {
480                 dict set defaultMoveParams $attr $value
481             } else {
482                 dict set moveDict $moveLabel $attr $value
483             }
484             expr {$value}
485         }
486         proposalCmd { ;# Argument is string or pure list.
487             dict set moveDict $moveLabel $attr $value
488             expr {$value}
489         }
490         default {
491             abort "cphTitratorSet: Invalid attribute $attr"
492             expr {-1}
493         }
494     }]
495 }
496
497 # =============================================================================
498 # Routines for tracking and reporting MC statistics
499 # =============================================================================
500 # ::cphTitrator::accumulateAcceptanceRate
501 #
502 # Accumulate statistics for the given move.
503 #
504 proc ::cphTitrator::accumulateAcceptanceRate {accept moveLabel} {
505     # Alas, dict incr does not support nested keys.
506     set successes [cphTitrator get successes $moveLabel]
507     set attempts [cphTitrator get attempts $moveLabel]
508     incr successes $accept
509     incr attempts
510     cphTitrator set successes $moveLabel $successes
511     cphTitrator set attempts $moveLabel $attempts
512     return
513 }
514