Bug fixes in Monte Carlo functionality
[namd.git] / lib / namdcph / namdcph / cphtoppar.tcl
1 # cphtoppar.tcl
2 #
3 #   This file provides the ::cphSystem namespace, which effectively emulates an
4 # object containing all constant-pH specific topology and parameter (toppar)
5 # information.
6 #
7 package require Tcl 8.5
8
9 source [file join [file dirname [info script]] "cphpsfgen.tcl"]
10 source [file join [file dirname [info script]] "namdmcmc.tcl"]
11
12 namespace eval ::cphSystem {
13     #   The core information of the cphSystem is the resDict, which stores
14     # the unique and mutable information for each instance of a titratable
15     # residue. The keys are "segresidnames" of the form:
16     #
17     # <segid>:<resid>:<resname>.
18     #
19     # Note the <resname> is a _titratable residue_ name and may or may not be
20     # the same as a recognizable residue from the force field. For example,
21     # amino acid termini are titratable residues. The resDefDict stores
22     # non-unique and immutable information for each residue definition. The
23     # keys are residue names (e.g. ASP, HIS, etc.). Generally, when looking up
24     # information for a specific residue, one simply needs look up the name in
25     # resDict and then look up the info in resDefDict.
26     #
27     # Note that a <segid>:<resid> combination can correspond to _multiple_
28     # residue names. This permits non-overlapping "sub-residues," such as amino
29     # acid termini, which titrate independently of the sidechain.
30     #
31     variable resDict [dict create]
32     variable resDefDict [dict create]
33
34     namespace export cphSystem
35 }
36
37 # =============================================================================
38 # cphSystem Interface function
39 # =============================================================================
40 # ::cphSystem::cphSystem
41 #
42 # This is the only exported function from the cphSystem namespace and provides
43 # a complete interface to the cphSystem namespace. The first argument is always
44 # an action (e.g. "get" or "set") and the rest are variable arguments.
45 #
46 # action           description
47 # ---------------- -----------
48 # get              see cphSystemGet 
49 # set              see cphSystemSet
50 # build            see buildSystem
51 # initialize       see initializeSystem
52 # propose          propose a move of the given type
53 # compute          compute various quantities needed for MC
54 # update           accept/reject the proposed changes
55 # alchemifypsf     apply patches that make the system alchemical
56 # dealchemifypsf   remove patches so the system is no longer alchemical
57 # initializeState  set the initial state of a residue
58 #
59 # Some notes about Tcl switch statements for non-guru types like myself:
60 #
61 #   These can appear a bit delicate to those who are not familiar with their
62 # quirks (compared to C). For one, badly placed comments can break them. Also
63 # note that a switch returns the expression it first matches and does NOT use
64 # break statements.
65 #
66 proc ::cphSystem::cphSystem {action args} {
67     return [switch -nocase -- $action {
68         get {
69             cphSystemGet {*}$args
70         }
71         set {
72             cphSystemSet {*}$args
73         }
74         build {
75             buildSystem {*}$args
76         }
77         initialize {
78             initializeSystem {*}$args
79         }
80         propose { ;# return 0 if no valid proposal is found
81             set type [lindex $args 0]
82             set newArgs [lrange $args 1 end]
83             switch -nocase -- $type {
84                 titration {
85                     proposeResidueTitration {*}$newArgs
86                 }
87                 tautomerization {
88                     proposeResidueTautomerization {*}$newArgs
89                 }
90                 protonTransfer {
91                     proposeProtonTransfer {*}$newArgs
92                 }
93                 cotitration {
94                     proposeCotitration {*}$newArgs
95                 }
96                 default {
97                     abort "Invalid proposal type $type"
98                 }
99             }
100         }
101         compute { ;# acceptance energy (inherent or switch correction)
102             set type [lindex $args 0]
103             set newArgs [lrange $args 1 end]
104             switch -nocase -- $type {
105                 inherent {
106                     computeInherentAcceptance {*}$newArgs
107                 }
108                 inherentWeights {
109                     computeInherentNormedWeights {*}$newArgs
110                 }
111                 switch {
112                     computeSwitchAcceptance {*}$newArgs
113                 }
114                 default {
115                     abort "Invalid energy type $type"
116                 }
117             }
118         }
119         update { ;# Update states based on an MC result
120             updateStates {*}$args
121         }
122         alchemifypsf { ;#Alchemify the PSF (in memory) based on trial states
123             alchemifyResidue {*}$args
124         }
125         dealchemifypsf { ;# Dealchemify the PSF (in memory)
126             dealchemifyResidue {*}$args
127         }
128         initializeState { ;# Assign a state
129             set method [lindex $args 0]
130             switch -nocase -- $method {
131                 random {
132                     randomizeState {*}[lrange $args 1 end]
133                 }
134                 default {
135                     assignState {*}$args
136                 }
137             }
138        }
139        validate { ;# Validate a segresidname
140            validateSegresidname {*}$args
141        }
142        default {
143            abort "Invalid cphSystem action $action."
144        }
145     }]
146 }
147
148 # ::cphSystem::validateSegresidname
149 #
150 # Check that a given segresidname is valid. That is:
151 #
152 # 1) is it in the correct <segid>:<resid>:<resname> format?
153 # 2) does it correspond to a known titratable residue?
154 #
155 # A non-zero error code is returned for each condition.
156 #
157 # NB! This should only be called _after_ buildSystem.
158 #
159 proc ::cphSystem::validateSegresidname {segresidname} {
160     lassign [split $segresidname ":"] segid resid resname
161     if {![info exists segid] || ![info exists resid]
162         || ![info exists resname]} {
163         print "segresidname selections must be of form\
164                  <segid>:<resid>:<resname>!"
165         return -1
166     }
167     if {$segresidname ni [cphSystem get segresidnames]} {
168         print "Invalid segresidname $segresidname!"
169         return -2
170     }
171     return 0
172 }
173
174 # =============================================================================
175 # Proposal Routines
176 # =============================================================================
177 # proc ::cphSystem::updateStates
178 #
179 # Update the state of one or more residues with the given
180 # <segid>:<resid>:<resname> specifications based on acceptance/rejection of the
181 # trial states. Reset the trial states to a null value.
182 #
183 proc ::cphSystem::updateStates {accept segresidnameList} {
184     if {$accept} {
185         foreach segresidname $segresidnameList {
186             cphSystem set state $segresidname\
187                     [cphSystem get trialState $segresidname]
188             cphSystem set trialState $segresidname {}
189         }
190     } else {
191         foreach segresidname $segresidnameList {
192             cphSystem set trialState $segresidname {}
193         }
194     }
195     return 0
196 }
197
198 # ::cphSystem::proposeResidueTitration
199 #
200 # Propose a new trial state requiring a titration - i.e. a net change in the 
201 # number of protons.
202 #
203 # For consistency with tautomers (see below), this always returns true in order
204 # to confirm that a titration was in fact found.
205 #
206 proc ::cphSystem::proposeResidueTitration {segresidname} {
207     set possibleStates [cphSystem get trialStateList $segresidname]
208     set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
209     while {true} {
210         lassign [choice $possibleStates] state
211         set occ [cphSystem get occupancy $segresidname $state]
212         set numProtonsTrial [expr {lsum($occ)}]
213         if {$numProtons != $numProtonsTrial} {
214             cphSystem set trialState $segresidname $state
215             return 1
216         }
217     }
218     # This is an error and should never happen.
219     return -1
220 }
221
222 # ::cphSystem::proposeProtonTransfer
223 #
224 # Propose a new trial state requiring a proton transfer - i.e. a movement of a
225 # proton from residue to another.
226 #
227 # This returns true if no proton transfer or else a negative error code.
228 #
229 proc ::cphSystem::proposeProtonTransfer {segresidname1 segresidname2} {
230     set numProtons1 [expr {lsum([cphSystem get occupancy $segresidname1])}]
231     set numProtons2 [expr {lsum([cphSystem get occupancy $segresidname2])}]
232     set possibleStates1 [cphSystem get trialStateList $segresidname1]
233     set possibleStates2 [cphSystem get trialStateList $segresidname2]
234     set dn [expr {$numProtons1 - $numProtons2}]
235     if {$dn == 0.0} {     
236         # No proton transfer is possible.
237         return 1
238     } 
239     # transfer from 1 to 2
240     while {true} {
241         lassign [choice $possibleStates1] state1
242         set occ1 [cphSystem get occupancy $segresidname1 $state1]
243         set numProtonsTrial1 [expr {lsum($occ1)}]
244         lassign [choice $possibleStates2] state2
245         set occ2 [cphSystem get occupancy $segresidname2 $state2]
246         set numProtonsTrial2 [expr {lsum($occ2)}]
247         set dnTrial [expr {$numProtonsTrial2 - $numProtonsTrial1}]
248         if {$dn == $dnTrial} {
249             cphSystem set trialState $segresidname1 $state1
250             cphSystem set trialState $segresidname2 $state2
251             return 0
252         }
253     }
254     # This is an error and should never happen.
255     return -1
256 }
257
258 # ::cphSystem::proposeCotitration
259 #
260 # Propose a new trial state requiring cotitration - i.e. two residues
261 # concurrently changing from protonated to deprotonated or vice versa.
262 #
263 # This returns true if no cotitration or else a negative error code.
264 #
265 proc ::cphSystem::proposeCotitration {segresidname1 segresidname2} {
266     set maxAttempts 10
267
268     set numProtons1 [expr {lsum([cphSystem get occupancy $segresidname1])}]
269     set numProtons2 [expr {lsum([cphSystem get occupancy $segresidname2])}]
270     set possibleStates1 [cphSystem get trialStateList $segresidname1]
271     set possibleStates2 [cphSystem get trialStateList $segresidname2]
272     while {true} {
273         lassign [choice $possibleStates1] state1
274         set occ1 [cphSystem get occupancy $segresidname1 $state1]
275         set numProtonsTrial1 [expr {lsum($occ1)}]
276         set dn1 [expr {$numProtonsTrial1 - $numProtons1}]
277         lassign [choice $possibleStates2] state2
278         set occ2 [cphSystem get occupancy $segresidname2 $state2]
279         set numProtonsTrial2 [expr {lsum($occ2)}]
280         set dn2 [expr {$numProtonsTrial2 - $numProtons2}]
281         if {$dn1 == $dn2} {
282             cphSystem set trialState $segresidname1 $state1
283             cphSystem set trialState $segresidname2 $state2
284             return 0
285         }
286         incr attempts
287         if {$attempts >= $maxAttempts} {
288             # This probably implies that no cotitration exists.
289             return 1
290         }
291     }
292     # This is an error and should never happen.
293     return -1
294 }
295
296 # ::cphSystem::proposeResidueTautomerization
297 #
298 # Propose a new trial state requiring a tautomerization - i.e. the number of 
299 # protons remains unchanged. 
300
301 # Unlike a state titration, a state tautomerization is not guaranteed to exist 
302 # for all residues. Return true if one is found, else return false.
303 #
304 # In order to ensure that a tautomer is found (if it exists), but also not
305 # cause an infinite loop when no tautomer exists, the number of state proposals
306 # is capped at maxAttempts. For a typical residue with a tautomer, there are
307 # three states, two of which can interconvert via tautomerization. The
308 # probability of selecting the tautomeric state k times in N trials is thus
309 # binomially distributed with p = (1 - p) = 0.5:
310 #
311 # P(k, N) = [N! / (k!(N-k)!)] 0.5^N
312 #
313 # The probability of picking that correct state at least once is thus:
314 #
315 # P(k>0, N) = 1 - P(0, N) = 1 - 0.5^N 
316 #
317 # or
318 #
319 # N = log[1 - P(k>0, N)] / log(0.5)
320 #
321 # For P(k>0, N) = 0.999, this gives N ~= 10 (the default).
322 #
323 proc ::cphSystem::proposeResidueTautomerization {segresidname} {
324     set maxAttempts 10
325
326     set possibleStates [cphSystem get trialStateList $segresidname]
327     set numProtons [expr {lsum([cphSystem get occupancy $segresidname])}]
328     while {true} {
329         lassign [choice $possibleStates] state
330         set occ [cphSystem get occupancy $segresidname $state]
331         set numProtonsTrial [expr {lsum($occ)}]
332         if {$numProtonsTrial == $numProtons} {
333             cphSystem set trialState $segresidname $state
334             return 1
335         }
336         incr attempts
337         if {$attempts >= $maxAttempts} {
338             # This probably implies that no tautomer exists. 
339             return 0 
340         }
341     }
342     # This is an error and should never happen.
343     return -1
344 }
345
346 # ::cphSystem::computeInherentLogRatio
347 #
348 # Compute the (common) log ratio* of inherent (i.e. pH-dependent) weights of
349 # the given states within the given residue.
350 #
351 # Because pKa is strictly pairwise, the relative unnormed weights (Q) are most
352 # clearly computed as ratios (and as log ratios, for numerical stability).
353 # There are two cases: 1) tautomerizations and 2) titrations, the former is not
354 # a proper pKa since the pH dependence cancels (i.e. no protons are entering or
355 # leaving the bath).
356 #
357 # case 1, n = n':
358 #
359 #     P(l')
360 #     ----- = 10^[sgn(l, l') pKa_i(l, l')]
361 #     P(l)
362 #
363 # case 2, n != n':
364 #     
365 #     P(l')
366 #     ----- = 10^{-[sgn(n' - n) pKa_i(l, l') - (n' - n)pH]}
367 #     P(l)
368 #
369 # where l and l' are occupation vectors for the given states and n and n' are
370 # the number of protons in each state (i.e. the sums of the occupation vector
371 # elements). By convention, pKa_i(l, l') = pKa_i(l', l) and the antisymmetry
372 # of adding vs deleting protons is accounted for by the sgn function. Note 
373 # that, for tautomers, l and l' are converted to indices and the sgn is
374 # computed by an arbitrary convention.
375 #
376 # Arguments:
377 # ----------
378 # pH : float
379 #   The pH value to evaluate the ratio at.
380 # segresidname : string
381 #   Residue specification as "<segid>:<resid>:<resname>" 
382 # statei : string
383 #   The state whose weight is in the numerator
384 # statej : string
385 #   The state whose weight is in the denomenator
386 #
387 # Returns:
388 # --------
389 # logRatio : float
390 #   The natural log of the ratio Qj/Qi, i.e. ln(Qj/Qi)
391 #
392 proc ::cphSystem::computeInherentLogRatio {pH segresidname statei statej} {
393     set multij [cphSystem get multiplicity $segresidname $statej]
394     if {$statei == $statej} {
395         # dn = pKai = 0 exactly, but lookup will fail.
396         return [expr {-log($multij)}]
397     }
398
399     set li [cphSystem get occupancy $segresidname $statei]
400     set lj [cphSystem get occupancy $segresidname $statej]
401     set dn [expr {lsum($lj) - lsum($li)}]
402     set si [occupancy2Index $li]
403     set sj [occupancy2Index $lj]
404     set sij [index2flatindex $si $sj]
405     set pKai [lindex [cphSystem get pKaiPair $segresidname] $sij]
406     if {$dn == 0.0} {
407         # tautomerization: "pKa" is positive in direction of lower index
408         set sgn [expr {$sj > $si} ? 1.0 : -1.0]
409     } else {
410         # titration: pKa is positive in direction of fewer protons
411         set sgn [expr {$dn > 0} ? 1.0 : -1.0]
412     }
413     return [expr {$::LN10*($sgn*$pKai - $dn*$pH) - log($multij)}]
414 }
415
416 # ::cphSystem::computeInherentAcceptance
417 #
418 # Compute the (reduced) energy difference for a Monte Carlo move based on the 
419 # given segresidname, its current and trial state, and the given pH.
420 #
421 # The Metropolis criterion can be expressed as:
422 #
423 #     P(l --> l') = min{1, Q'/Q} = min{1, e^[-du(l, l')]},
424 #
425 # where
426 #
427 #     du = -ln(Q'/Q)
428 #
429 proc ::cphSystem::computeInherentAcceptance {pH segresidname} {
430     set statei [cphSystem get state $segresidname]
431     set statej [cphSystem get trialState $segresidname]
432     set du [computeInherentLogRatio $pH $segresidname $statei $statej]
433     set multij [cphSystem get multiplicity $segresidname $statej]
434     # Note the sign flip for compatibility with Metropolis acceptance!
435     return [expr {-($du + log($multij))}]
436 }
437
438 # ::cphSystem::computeInherentNormedWeights
439 #
440 # Compute the normalized inherent pKa weights (i.e. probability mass function,
441 # or PMF) of all states within the given residue at the given pH.
442 #
443 # The resulting PMF can be directly sampled in order to choose a new state ala
444 # independence sampling. However, this is often not an efficent scheme for
445 # Markov chain Monte Carlo and the output can also be used for other approaches
446 # such as Metropolized independence sampling.
447 #
448 proc ::cphSystem::computeInherentNormedWeights {pH segresidname} {
449     set currState [cphSystem get state $segresidname]
450     set stateList [cphSystem get stateList $segresidname] 
451
452     # Note that all ratios are taken wrt the current state and are thus
453     # equivalent to unnormed weights. This would NOT be correct if a different
454     # state were in the denomenator each time.
455     #
456     set logQs [list]
457     set logQMax 0.0
458     foreach state $stateList {
459         set logQ [computeInherentLogRatio $pH $segresidname $currState $state]
460         lappend logQs $logQ
461         if {$logQ > $logQMax} {
462             set logQMax $logQ
463         }
464     }
465     set normedWeights [normalizeLogWeights $logQs $logQMax]
466     return [list $stateList $normedWeights]
467 }
468
469 # ::cphSystem::computeSwitchAcceptance
470 #
471 # Compute the pairwise energy difference used to shift the work value in the
472 # Monte Carlo acceptance criteria.
473 #
474 #   This is only the configuration _independent_ portion of that difference 
475 # (i.e. it only depends on the current and trial protonation states). The
476 # configuration _dependent_ portion is the total work W applied to effect the
477 # switch. As for the proposal energy, the sign of the correction depends on
478 # whether this is tautomerization (1) or titration (2) move:
479 #
480 # dE(s, s') = c*{-dG(s, s') + kT ln(10) [pKa(s, s') - pKa_i(s, s')]}
481
482 # case 1, n = n':  c = sgn(s' - s)
483 #
484 # case 2, n != n': c = sgn(n' - n)
485 #
486 # Here T is the _reference_ temperature at which pKa is measured and dG is
487 # computed, not the temperature of the simulation. See 
488 # computeInherentAcceptance for definition of the remaining notation.
489 #
490 proc ::cphSystem::computeSwitchAcceptance {segresidnameList} {
491     set du 0.0
492     foreach segresidname $segresidnameList {
493         set l [cphSystem get occupancy $segresidname]
494         set lp [cphSystem get trialOccupancy $segresidname]
495         set dn [expr {lsum($lp) - lsum($l)}]
496         set s [occupancy2Index $l]
497         set sp [occupancy2Index $lp]
498         set ssp [index2flatindex $s $sp]
499         set dG [lindex [cphSystem get dGPair $segresidname] $ssp]
500         set pKa [lindex [cphSystem get pKaPair $segresidname] $ssp]
501         set pKai [lindex [cphSystem get pKaiPair $segresidname] $ssp]
502         set trialState [cphSystem get trialState $segresidname]
503         if {$dn == 0.0} {
504             # tautomerization: "pKa" is positive in direction of lower index
505             set sgn [expr {$sp > $s} ? 1.0 : -1.0]
506         } else {
507             # titration: pKa is positive in direction of fewer protons
508             set sgn [expr {$dn > 0} ? 1.0 : -1.0]
509         }
510         set kT [expr {$::BOLTZMANN*[cphSystem get Tref $segresidname]}]
511         set du [expr {$du + $sgn*($dG - $kT*$::LN10*($pKa - $pKai))}]
512     }
513     return $du
514 }
515
516 # ---------------
517 # Psfgen Routines
518 # ---------------
519 # ::cphSystem::alchemifyResidue
520 #
521 # Apply a trial alchemical patch to a residue.
522 #
523 proc ::cphSystem::alchemifyResidue {segresidname frcCons temp {buildH false}} {
524     lassign [cphSystem get alchAtomLists $segresidname] l0atoms l1atoms
525     set patch [cphSystem get hybridPatch $segresidname]
526     lassign [split $segresidname ":"] segid resid
527     # NB: alchPatch uses psfgen style selections!
528     alchPatch $patch "$segid:$resid" $l0atoms $l1atoms $frcCons $temp $buildH
529     return 0
530 }
531
532 # ::cphSystem::dealchemifyResidue
533 #
534 # Remove an alchemical patch from a residue.
535 #
536 proc ::cphSystem::dealchemifyResidue {segresidname} {
537     lassign [cphSystem get alchAtomLists $segresidname] l0atoms l1atoms
538     lassign [split $segresidname ":"] segid resid 
539     # NB: alchUnpatch uses psfgen style selections!
540     alchUnpatch "$segid:$resid" $l0atoms $l1atoms
541     return 0
542 }
543
544 # =============================================================================
545 # "Constructor" Routines
546 # =============================================================================
547 # ::cphSystem::initializeSystem
548 #
549 # Initialize the state of the system using:
550 #
551 # 1) input data from the user (usually from a restart file)
552 #
553 # AND/OR
554 #
555 # 2) the ensemble information (i.e. the pH and temperature)
556 #
557 proc ::cphSystem::initializeSystem {pH temperature buildH stateInfo} {
558     variable ::cphSystem::resDict
559     dict for {segresidname resData} $resDict {
560         # Assign inherent pKa values.
561         if {[dict exists $stateInfo $segresidname pKai]} {
562             cphSystem set pKai $segresidname\
563                     [dict get $stateInfo $segresidname pKai]
564         } else { ;# Default to reference pKa.
565             cphSystem set pKai $segresidname [cphSystem get pKa $segresidname]
566         }
567
568         # Assign states.
569         if {[dict exists $stateInfo $segresidname state]} {
570             set state [dict get $stateInfo $segresidname state]
571             cphSystem initializeState $segresidname $state 
572         } else { ;# default randomization based on pKai and pH
573             cphSystem initializeState random $segresidname $pH
574         }
575     }
576     # Final pass - Apply the patches
577     foreach segresidname [cphSystem get segresidnames] {
578         lassign [split $segresidname ":"] segid resid
579         patch [cphSystem get statePatch $segresidname] "$segid:$resid"
580     }
581     guesscoord
582     foreach segresidname [cphSystem get segresidnames] {
583         cphSystem alchemifypsf $segresidname 0.0 $temperature $buildH
584         cphSystem dealchemifypsf $segresidname
585         cphSystem update 1 $segresidname
586     }
587     regenerate angles dihedrals
588     # NB - These changes are only reflected in _memory_ for the cphSystem and
589     # psfgen. Nothing has happened to the NAMD PSF/PDB files.
590     return
591 }
592
593 # ::cphSystem::buildSystem
594 #
595 # Build the residue definitions and residue objects for the system based on 
596 # the NAMD inputs.
597 #
598 proc ::cphSystem::buildSystem {resDefs excludeList} {
599     variable ::cphSystem::resDict
600     variable ::cphSystem::resDefDict [checkResDefs $resDefs]
601
602     # Check for special sub-type residues.
603     # In general, THESE WILL NOT BE MATCHED BY NAME, but rather are sub-typed
604     # within one or more other residue definitions. Therefore we keep a
605     # separate dict of residue names that have a sub-type (these may not even
606     # be titratable otherwise, such as terminal amino acids). The keys are
607     # the actual residue names and the values are a list of the possible
608     # sub-types. However, even if a sub-typed residue is present, that doesn't
609     # mean the sub-residue exists. There are two (optional) additional checks
610     # in the sub-residue definition:
611     #
612     # subatoms - a list of atom names that MUST exist
613     # notsubatoms - a list of atom names that MUST NOT exist
614     #
615     # For example, alanine is not titratable, but may contain a titratable
616     # C-terminus. The terminal oxygen atoms must exist, but are only titratable
617     # if they are not blocked by amination or methylation.
618     #
619     set subresDefDict [dict create]
620     dict for {subresname resDef} $resDefDict {
621         if {![dict exists $resDef subtypeof]} continue
622         foreach subtypedRes [dict get $resDef subtypeof] {
623             if {![dict exists $subresDefDict $subtypedRes]} {
624                 dict set subresDefDict $subtypedRes [list]
625             }
626             dict lappend subresDefDict $subtypedRes $subresname
627         }
628     }
629
630     # Check for residue aliases (i.e. titratable residues whose states may span
631     # multiple residue definitions, such as HIS).
632     set resAliases [dict create]
633     dict for {resname resDef} $resDefDict {
634         if {![dict exists $resDef aliases]} continue
635         dict set resAliases $resname [dict get $resDef aliases]
636     }
637
638     # Read in whatever files were specified to NAMD.
639     set Args [list [structure] pdb [coordinates]]
640     if {[isset binCoordinates]} {
641         lappend Args namdbin [binCoordinates]
642         if {[isset binVelocities]} {
643             lappend Args velnamdbin [binVelocities]
644         }
645     }
646     resetpsf
647     readpsf {*}$Args
648
649     set definedResidues [dict keys $resDefDict]
650     set definedSubResidues [dict keys $subresDefDict]
651     foreach segid [segment segids] {
652         foreach resid [segment resids $segid] {
653             set resname [segment residue $segid $resid]
654             set segresid [format "%s:%s" $segid $resid]
655             set segresidname [format "%s:%s" $segresid $resname]
656             # Perform aliasing to fix residues with multiple names.
657             dict for {realName altNames} $resAliases {
658                 if {$resname ni $altNames} {
659                     continue
660                 }
661                 print "aliasing $segresidname to $realName"
662                 psfset resname $segid $resid $realName
663                 set resname $realName
664                 set segresidname [format "%s:%s" $segresid $resname]
665             }
666
667             set resIsDefined [expr {$resname in $definedResidues}]          
668             set resIsSubtyped [expr {$resname in $definedSubResidues}]
669             set resIsExcluded [expr {$segresidname in $excludeList}]
670             # Break here if nothing to do. Be sure to remove successful
671             # exclusions from the list so that we can identify bad selections
672             # later.
673             if {$resIsExcluded} {
674                 set excludeList [lsearch -inline -all -not $excludeList\
675                         $segresidname]
676             }
677             if {(!$resIsDefined && !$resIsSubtyped) || $resIsExcluded} {
678                 continue
679             }
680
681             # Add the residue to the system.
682             if {$resIsDefined} {
683                 dict set resDict $segresidname [dict create]
684                 dict set resDict $segresidname resname $resname
685             }
686
687             # If the residue is subtype-able, check that subresidues exist.
688             # Note that this might be true even if the residue itself is not
689             # titratable.
690             if {!$resIsSubtyped} continue
691
692             set atoms [segment atoms $segid $resid]
693             foreach subresname [dict get $subresDefDict $resname] {
694                 # Atoms that MUST exist.
695                 set allSubatomsExist 1
696                 set subatoms [list]
697                 if {[dict exists $resDefDict $subresname subatoms]} {
698                     set subatoms [dict get $resDefDict $subresname subatoms]
699                 }
700                 foreach atom $subatoms { 
701                     if {$atom ni $atoms} {
702                         set allSubatomsExist 0
703                         break
704                     }
705                 }
706
707                 # Atoms that must NOT exist.
708                 set allNotSubatomsDoNotExist 1
709                 set notsubatoms [list]
710                 if {[dict exists $resDefDict $subresname notsubatoms]} {
711                     set notsubatoms\
712                             [dict get $resDefDict $subresname notsubatoms]
713                 }
714                 foreach atom $notsubatoms {
715                     if {$atom in $atoms} {
716                         set allNotSubatomsDoNotExist 0
717                         break
718                     }
719                 }
720
721                 # Break here if anything doesn't match. 
722                 if {!$allSubatomsExist || !$allNotSubatomsDoNotExist} continue
723
724                 set segresidname [format "%s:%s" $segresid $subresname]
725                 # After all of that, the residue might be excluded!
726                 set resIsExcluded [expr {$segresidname in $excludeList}]
727                 if {$resIsExcluded} {
728                     set excludeList [lsearch -inline -all -not $excludeList\
729                             $segresidname]
730                     continue
731                 }
732
733                 # Add the sub-residue to the system.
734                 dict set resDict $segresidname [dict create]
735                 dict set resDict $segresidname resname $subresname
736             }
737         }
738     }
739     # This is a bit of a hack way to check that all of the selected exclusions
740     # were valid, but it works.
741     if {[llength $excludeList]} {
742         abort "Cannot exclude non-existant residue(s): $excludeList"
743     }
744
745     return $resDict
746 }
747
748 # ::cphSystem::checkResDefs
749 #
750 # Check that a dictionary of residue definitions is valid and consistent with
751 # the topology (RTF) files currently loaded by psfgen.
752 #
753 proc ::cphSystem::checkResDefs {resDefs} {
754     dict for {resname resDef} $resDefs {
755         # Check that all required fields are present.
756         foreach dtype [list dG pKa states l0atoms] {
757             if {![dict exists $resDef $dtype]} {
758                 abort "No $dtype entry for residue $resname!"
759             }
760         }
761         set pKa [dict get $resDef pKa]
762         set dG [dict get $resDef dG]
763         set states [dict get $resDef states]
764         set l0atoms [dict get $resDef l0atoms]
765         # If not specified, use lazy method of appending a "1". If an l1atoms
766         # field is present, check that it matches l0atoms.
767         if {![dict exists $resDef l1atoms]} {
768             set l1atoms [list]
769             foreach atom $l0atoms { 
770                 lappend l1atoms [format "%s1" $atom]
771             }
772             dict set resDefs $resname l1atoms $l1atoms
773         } else {
774             set l1atoms [dict get $resDefs $resname l1atoms]
775         }
776
777         if {[llength $l0atoms] != [llength $l1atoms]} { 
778             abort "Mismatch in atom definitions for residue $resname"
779         }
780         # Check that the definitions for pKa and dG are consistent.
781         if {[llength $pKa] != [llength $dG]} {
782             abort "Mismatch in dG/pKa definitions for residue $resname."
783         }
784         # Check that the state and occupancy definitions are consistent.
785         set numSites [llength [dict get $states [lindex $states 0]]]
786         set maxStates [expr {int(pow(2, $numSites))}]
787         if {[dict size $states] > $maxStates} {
788             abort "Too many states defined for residue $resname!"
789         }
790         dict for {state occ} $states {
791             if {[llength $occ] != $numSites} {
792                 abort "Bad occupancy definition for $resname state $state."
793             }
794             # Check that the RTFs define two patches for each state:
795             # 1) state patches, with the naming convention: <resname><state>
796             # 2) hybrid patches, with the naming convention: <resname>H<state>
797             #
798             set statePatch [format "%s%s" $resname $state]
799             set hybridPatch [format "%sH%s" $resname $state]
800             if {$statePatch ni [topology patches]} {
801                 abort "No patch definition in RTFs for $statePatch!"
802             }
803             if {$hybridPatch ni [topology patches]} {
804                 abort "No patch definition in RTFs for $hybridPatch!"
805             }
806         }
807         # Build the pairwise parameters for this definition.
808         dict set resDefs $resname dGPair [resDef2Matrix $resDef $dG]
809         dict set resDefs $resname pKaPair [resDef2Matrix $resDef $pKa]
810     }
811     return $resDefs
812 }
813
814 # ::cphSystem::assignState
815 #
816 # Assign a protonation state by fiat.
817 #
818 #   This is a bit of a hack. The state is assigned and then an arbitrary trial
819 # state is chosen by rejectionless MC. The states are then swapped and the move
820 # is accepted as if this had happened in reverse.
821 #
822 # Arguments:
823 # ----------
824 # segresidname : string
825 #   residue specification as "<segid>:<resid>:<resname>"
826 # state : string
827 #   state to assign 
828 #
829 # Returns:
830 # --------
831 # None
832 #
833 proc ::cphSystem::assignState {segresidname state} {
834     cphSystem set state $segresidname $state
835     cphSystem propose titration $segresidname
836     cphSystem set state $segresidname [cphSystem get trialState $segresidname]
837     cphSystem set trialState $segresidname $state
838     return 0
839 }
840
841 # ::cphSystem::randomizeState
842 #
843 # Assign a protonation state by independence sampling at the given pH.
844 #
845 # Arguments:
846 # ----------
847 # segresidname : string
848 #   residue specification as "<segid>:<resid>:<resname>"
849 # pH : float
850 #   pH value at which to assign protonation states
851 #
852 # Returns:
853 # --------
854 # None 
855 #
856 proc ::cphSystem::randomizeState {segresidname pH} {
857     # Note that computeInherentNormedWeights requires a state to be set for
858     # normalization purposes. Therefore choose a random state first and then
859     # rigorously sample the distribution.
860     set states [cphSystem get stateList $segresidname]
861     lassign [choice $states] state
862     cphSystem set state $segresidname $state
863     lassign [computeInherentNormedWeights $pH $segresidname] states weights
864     lassign [choice $states $weights] state
865     cphSystem set state $segresidname $state
866     # Use the same hack as in assignState.
867     cphSystem propose titration $segresidname
868     cphSystem set state $segresidname [cphSystem get trialState $segresidname]
869     cphSystem set trialState $segresidname $state
870     return 0
871 }
872
873 # =============================================================================
874 # Getter Routines
875 # =============================================================================
876 # ::cphSystem::cphSystemGet 
877 #
878 # Getters for system and residue attributes, called as:
879 #
880 #   <attribute> [<segresidname> [<args>]]
881 #
882 # <attribute> is the name of either a system attribute (segresidname selections
883 # are invalid) or else a residue attribute. For SOME residue attributes, not
884 # specifying a segresidname will return a list for all residues. Some
885 # attributes also require some additional arguments (see below)
886 #
887 # system attributes  description
888 # -----------------  -----------
889 # segresidnames      list of all segresidnames
890 # numresidues        number of residues in the system
891 # resdefs            list of defined resnames
892 # numdefs            number of defined resnames
893 #
894 # residue attributes description
895 # ------------------ -----------
896 # resname            residue name
897 # state              current state
898 # trialState         proposed trial state
899 # pKai               minimal pKai list for this residue
900 # dG                 minimal dG list for this residue 
901 # pKa                minimal pKa list for this residue
902 # Tref               reference temperature for pKa
903 # occupancy          occupancy vector for the current state
904 # trialOccupancy     occupancy vector for the trial state
905 # stateList          all possible states
906 # equivStateList     all states equivalent to the given state
907 # trialStateList     all possible trial states (not the current state)
908 # dGPair             pair dG for current/trial states
909 # pKaPair            pair pKa for current/trial states
910 # pKaiPair           pair pKai for current/trial states 
911 # statePatch         patch for the current state
912 # hybridPatch        alchemical patch for the trial state
913 # alchAtomLists*     lists of alchemical atoms at 0/1
914 # alchBonds*^        extraBonds entries
915 #
916 # * - segresidname selection is required
917 # ^ - requires additional arguments
918 #
919 proc ::cphSystem::cphSystemGet {attr {segresidname {}} args} {
920     variable ::cphSystem::resDict
921     variable ::cphSystem::resDefDict
922
923     set getAll [expr {![llength $segresidname]}]
924     if {!$getAll && ![dict exists $resDict $segresidname]} {
925         abort "cphSystemGet: Invalid segresidname $segresidname"
926     }
927
928     return [switch -nocase -- $attr {
929         segresidnames {
930             dict keys $resDict
931         }
932         numresidues {
933             dict size $resDict
934         }
935         resdefs {
936             dict keys $resDefDict
937         }
938         numdefs {
939             dict size $resDefDict
940         }
941         resname -
942         state -
943         trialState -
944         pKai {
945             if {$getAll} {
946                 getAllResAttr $attr
947             } else {
948                 getResAttr $attr $segresidname
949             }
950         }
951         dG -
952         pKa -
953         Tref {
954             if {$getAll} {
955                 getAllResDefAttr $attr
956             } else {
957                 getResDefAttr $attr $segresidname
958             }
959         }
960         occupancy {
961             if {$getAll} {
962                 getAllOccupancy
963             } else {
964                 if {[llength $args]} {
965                     lassign $args state
966                     getOccupancy $segresidname $state
967                 } else {
968                     getCurrOccupancy $segresidname
969                 }
970             }
971         }
972         trialOccupancy {
973             if {$getAll} {
974                 cannotGetAll $attr
975             } else {
976                 getTrialOccupancy $segresidname
977             }
978         }
979         stateList {
980             if {$getAll} {
981                 cannotGetAll $attr 
982             } else {
983                 set resname [getResAttr resname $segresidname]
984                 dict keys [dict get $resDefDict $resname states]
985             }
986         }
987         equivStateList {
988             if {$getAll} {
989                 cannotGetAll $attr
990             } else {
991                 if {[llength $args]} {
992                     lassign $args state
993                 } else {
994                     # Default to current state.
995                     set state [getResAttr state $segresidname]
996                 }
997                 getEquivStates $segresidname $state
998             }
999         }
1000         multiplicity {
1001             if {$getAll} {
1002                 cannotGetAll $attr
1003             } else {
1004                 if {[llength $args]} {
1005                     lassign $args state
1006                 } else {
1007                     # Default to current state.
1008                     set state [getResAttr state $segresidname]
1009                 }
1010                 llength [getEquivStates $segresidname $state]
1011             }
1012         }
1013         trialStateList {
1014             if {$getAll} {
1015                 cannotGetAll $attr
1016             } else {
1017                 set resname [getResAttr resname $segresidname]
1018                 set state [getResAttr state $segresidname]
1019                 set states [dict keys [dict get $resDefDict $resname states]]
1020                 lsearch -all -inline -not -nocase $states $state
1021             }
1022         }
1023         dGPair -
1024         pKaPair {
1025             if {$getAll} {
1026                 cannotGetAll $attr
1027             } else {
1028                 getResDefAttr $attr $segresidname
1029             }
1030         }
1031         pKaiPair {
1032             if {$getAll} {
1033                 cannotGetAll $attr
1034             } else {
1035                 getResAttr $attr $segresidname
1036             }
1037         }
1038         statePatch {
1039             if {$getAll} {
1040                 cannotGetAll $attr
1041             } else {
1042                 set resname [getResAttr resname $segresidname]
1043                 set state [getResAttr state $segresidname]
1044                 format "%s%s" $resname $state
1045             }
1046         }
1047         hybridPatch {
1048             if {$getAll} {
1049                 cannotGetAll $attr
1050             } else {
1051                 set resname [getResAttr resname $segresidname]
1052                 set trialState [getResAttr trialState $segresidname]
1053                 format "%sH%s" $resname $trialState
1054             }
1055         }
1056         alchAtomLists {
1057             if {$getAll} {
1058                 cannotGetAll $attr
1059             } else {
1060                 list [getResDefAttr l0atoms $segresidname]\
1061                      [getResDefAttr l1atoms $segresidname]
1062             }
1063         }
1064         alchBonds {
1065             if {$getAll} {
1066                 cannotGetAll $attr
1067             } else {
1068                 lassign $args k
1069                 set bondEntries [list]
1070                 foreach l0atom [getResDefAttr l0atoms $segresidname]\
1071                         l1atom [getResDefAttr l1atoms $segresidname] {
1072                     lassign [split $segresidname ":"] segid resid
1073                     # Note that atomid indices start at one, not zero!
1074                     set i [expr {[segment atomid $segid $resid $l0atom] - 1}]
1075                     set j [expr {[segment atomid $segid $resid $l1atom] - 1}]
1076                     lappend bondEntries [format "bond %d %d %f %f" $i $j $k 0]
1077                 }
1078                 join $bondEntries "\n"
1079             }
1080         }
1081         default {
1082             abort "cphSystemGet: Invalid attribute $attr"
1083         }
1084     }]
1085 }
1086
1087 proc ::cphSystem::cannotGetAll {attr} {
1088     abort "cphSystemGet: Cannot get all $attr - must select a segresidname"
1089     return -1
1090 }
1091
1092 # ------------------------------
1093 # Getters for residue attributes
1094 # ------------------------------
1095 proc ::cphSystem::getResAttr {attr segresidname} {
1096     variable ::cphSystem::resDict
1097     return [dict get $resDict $segresidname $attr]
1098 }
1099
1100 proc ::cphSystem::getAllResAttr {attr} {
1101     variable ::cphSystem::resDict
1102     set retList [list]
1103     foreach resData [dict values $resDict] {
1104         lappend retList [dict get $resData $attr]
1105     }
1106     return $retList
1107 }
1108
1109 # -----------------------------------------
1110 # Getters for residue definition attributes
1111 # -----------------------------------------
1112 proc ::cphSystem::getResDefAttr {attr segresidname} {
1113     variable ::cphSystem::resDict
1114     variable ::cphSystem::resDefDict
1115     set resname [dict get $resDict $segresidname resname]
1116     return [dict get $resDefDict $resname $attr]
1117 }
1118
1119 proc ::cphSystem::getAllResDefAttr {attr} {
1120     variable ::cphSystem::resDict
1121     variable ::cphSystem::resDefDict
1122     set retList [list]
1123     foreach resData [dict values $resDict] { 
1124         set resname [dict get $resData resname]
1125         lappend retList [dict get $resDefDict $resname $attr]
1126     }
1127     return $retList
1128 }
1129
1130 # -----------------------------------------------
1131 # Special getters for more complicated operations
1132 # -----------------------------------------------
1133 proc ::cphSystem::state2Occupancy {resname state} {
1134     return [dict get $::cphSystem::resDefDict $resname states $state]
1135 }
1136
1137 proc ::cphSystem::getOccupancy {segresidname state} {
1138     variable ::cphSystem::resDict
1139     set resname [dict get $resDict $segresidname resname]
1140     return [state2Occupancy $resname $state]
1141 }
1142
1143 proc ::cphSystem::getCurrOccupancy {segresidname} {
1144     variable ::cphSystem::resDict
1145     set state [dict get $resDict $segresidname state]
1146     return [getOccupancy $segresidname $state]
1147 }
1148
1149 proc ::cphSystem::getTrialOccupancy {segresidname} {
1150     variable ::cphSystem::resDict
1151     set state [dict get $resDict $segresidname trialState]
1152     return [getOccupancy $segresidname $state] 
1153 }
1154
1155 proc ::cphSystem::getAllOccupancy {} {
1156     variable ::cphSystem::resDict
1157     variable ::cphSystem::resDefDict
1158     set retList [list]
1159     foreach resData [dict values $resDict] {
1160         set resname [dict get $resData resname]
1161         set state [dict get $resData state]
1162         set retList [list {*}$retList {*}[state2Occupancy $resname $state]]
1163     }
1164     return $retList
1165 }
1166
1167 # ::cphSystem::getEquivStates
1168 #
1169 # Return the list of states equivalent to the given state of the given residue. # This includes the given state such that the length of the list is exactly the
1170 # degeneracy.
1171 #
1172 proc ::cphSystem::getEquivStates {segresidname testState} {
1173     set l [cphSystem get occupancy $segresidname $testState]
1174     set s [occupancy2Index $l]
1175
1176     set equivStateList [list]
1177     foreach state [cphSystem get stateList $segresidname] {
1178         if {$state == $testState} {
1179             lappend equivStateList $state
1180             continue
1181         }
1182
1183         set lp [cphSystem get occupancy $segresidname $state]
1184         set dn [expr {lsum($lp) - lsum($l)}]
1185         set sp [occupancy2Index $lp]
1186         set ssp [index2flatindex $s $sp]
1187         set pKa [lindex [cphSystem get pKaPair $segresidname] $ssp]
1188         if {$dn == 0.0 && $pKa == 0.0} {
1189             lappend equivStateList $state
1190         }
1191     }
1192     return $equivStateList
1193 }
1194
1195 # ----------------------
1196 # Other helper functions
1197 # ----------------------
1198 # Given an occupancy, get the state index.
1199 # This is essentially a left to right conversion from binary to decimal.
1200 # That is, 
1201 #   index = l0*2**0 + l1*2**1 + l2*2**2 + ...
1202 #
1203 proc ::cphSystem::occupancy2Index {occupancy} {
1204     set index 0
1205     for {set i 0} {$i < [llength $occupancy]} {incr i} {
1206         incr index [expr {int([lindex $occupancy $i]*pow(2, $i))}]
1207     }
1208     return $index
1209 }
1210
1211 # Map a pair index to the flat off-diagonal index.
1212 #
1213 # This maps to the _lower_ off-diagonal counting left-right, top-bottom.
1214 #
1215 # Ex. n = 4, 2,0 --> 1 and 3,2 --> 5
1216 #
1217 # |x x x x|
1218 # |0 x x x|
1219 # |1 2 x x|
1220 # |3 4 5 x|
1221 #
1222 # The additional assumption here is that element i,j is the same as j,i.
1223 # Therefore this internally swaps the two indices such that i > j _always_.
1224 #
1225 proc ::cphSystem::index2flatindex {i j} {
1226 #    lassign [list [expr {max($i, $j)}] [expr {min($i, $j)}]] I J
1227     if {$i > $j} {
1228         set I $i
1229         set J $j
1230     } elseif {$i < $j} {
1231         set I $j
1232         set J $i
1233     } else {
1234         abort "invalid transition from state $i to state $i"
1235     }
1236     return [expr {($I*($I - 1) + 2*$J) / 2}]
1237 }
1238
1239 # Convenience function for assignment of flat, off-diagonal, upper-triangular
1240 # matrices using 2d indices. Note that not all index combinations are valid!
1241 #
1242 # Example:
1243 #
1244 # % set myMatrix [lrepeat 3 0.0]
1245 # 0.0 0.0 0.0
1246 # % mset myMatrix 0 1 10.0
1247 # 10.0 0.0 0.0
1248 #
1249 proc ::cphSystem::mset {matrix i j value} {
1250     upvar 1 $matrix Matrix
1251     lset Matrix [index2flatindex $i $j] $value
1252     return $Matrix
1253 }
1254
1255 # Convenience function - inverse of mset.
1256 #
1257 proc ::cphSystem::mindex {matrix i j} {
1258     return [lindex $matrix [index2flatindex $i $j]]
1259 }
1260
1261 # =============================================================================
1262 # Setter Routines
1263 # =============================================================================
1264 # ::cphSystem::cphSystemSet
1265 #
1266 # Setters for residue attributes, called as:
1267 #
1268 #   <attribute> <segresidname> <value>
1269 #
1270 # <attribute> is the name of a residue attribute.
1271 #
1272 # residue attributes description
1273 # ------------------ -----------
1274 # state              current state
1275 # trialState         proposed trial state
1276 # pKai               minimal pKai list for this residue
1277 #
1278 proc ::cphSystem::cphSystemSet {attr segresidname value} {
1279     variable ::cphSystem::resDict
1280     variable ::cphSystem::resDefDict
1281
1282     if {![dict exists $resDict $segresidname]} {
1283         abort "cphSystemSet: Invalid segresidname $segresidname"
1284     }
1285
1286     return [switch -nocase -- $attr {
1287         state -
1288         trialState {
1289             set states [cphSystem get stateList $segresidname]
1290             if {$value ni $states} {
1291                 if {[string match -nocase $attr trialState]
1292                     && ![llength $value]} {
1293                 } else {
1294                     abort "Invalid state assignment $value for residue"\
1295                             "$segresidname"
1296                 }
1297             }
1298             dict set resDict $segresidname $attr $value
1299             expr {$value}
1300         }
1301         pKai {
1302             set resname [cphSystem get resname $segresidname]
1303             set resDef [dict get $resDefDict $resname]
1304             set pKaiMatrix [resDef2Matrix $resDef $value]
1305             dict set resDict $segresidname pKai $value
1306             dict set resDict $segresidname pKaiPair $pKaiMatrix
1307             expr {$pKaiMatrix}
1308         }
1309         default {
1310             abort "cphSystemSet: Invalid attribute $attr"
1311             expr {-1}
1312         }
1313     }]
1314 }
1315
1316 # ::cphSystem::resDef2Matrix
1317 #
1318 # Read a residue definiition for a given thermodynamic attribute (a free energy
1319 # or pKa) and fill in the (flat)matrix elements that are derived from this.
1320
1321 # This permits the configuration files to have "missing" information when it is
1322 # not used. The downside is that essentially every case needs to have its own
1323 # conventions.
1324 #
1325 # Arguments:
1326 # ----------
1327 # resDef : dict
1328 #   The residue definition to be assigned
1329 # data : list
1330 #   The _minimial_ attribute data needed to fill in the matrix
1331 #
1332 # Returns
1333 # -------
1334 #  0 - The matrix was successfully updated.
1335 # -1 - The specific combination of state definitions is not implemented
1336 #      Ex. A state is unexpectedly missing and leaves a gap in the cycle or
1337 #      too many definitions are given
1338 # -2 - There is an inconsistency in the truncated definition
1339 #      Ex. The sum of parameters around the cycle is not 0 - THIS MIGHT SIMPLY
1340 #      IMPLY A NEGATIVE SIGN ERROR!
1341 #
1342 proc ::cphSystem::resDef2Matrix {resDef data} {
1343     #   The logic and mathematics here is a bit dreadful. Basically, we have
1344     # networks of states based on the number of protonation sites (numSites).
1345     # A given network has, at most, maxStates = 2**numSites states, but
1346     # generally many fewer than this, because some of the states don't exist
1347     # for physical reasons (e.g. doubly protonated carboxylates).
1348     #   The problems begin because the residue attributes (dG and pKa) do not
1349     # describe states, but rather transitions between pairs of states - there
1350     # are formally numStates**2 pairs! Fortunately, attributes between
1351     # identical states will always be zero and transitions between pairs of
1352     # states in opposite directions just have flipped signs. There are
1353     # therefore only numStates*(numStates - 1)/2 pairs of states that need to
1354     # be stored (as a flat, off-diagonal matrix, see index2flatindex, note that
1355     # this conversion is symmetric!). Lookup is always done with respect to the
1356     # full set of states (i.e. the size is computed with maxStates, not
1357     # numStates). The row/column indices are computed as a left-right binary
1358     # conversion of the occupancy vector to a scalar (see occupancy2Index).
1359     #   Just to make things more complicated, many systems have equivalent
1360     # sites such that many fewer pairwise values are needed. For example, a
1361     # carboxylate formally has two sites and four states, one of which is
1362     # neglected. The two sites are also equivalent, so this system has, at max
1363     # 2**2*(2**2 - 1)/2 = 6 pairs, but only three are used in practice. Two of
1364     # these pairs are the same, so only a single value is specified.
1365     #
1366     # Long story short: Each case is handled specially and implemented by hand
1367     # (YUCK!). Fortunately, it is the state pattern that is implemented, not
1368     # the residue name or some overly specific nonsense. If a pattern has been
1369     # implemented, then any residue following that pattern works also (e.g.
1370     # ACE, ASP, and GLU are all one pattern).
1371     #
1372     # Final note:
1373     #   Because only half of the diagonal elements are stored, the signs are
1374     # computed on-the-fly. This assumes that all proper pKa values are positive
1375     # - a negative pKa should never be inserted into pKa matrix. For "improper
1376     # pKas" (i.e. tautomerizations), the assumption is that the value is from
1377     # the higher index to the lower index (this may be negative).
1378     #   In the reduced listing all unique transitions from the highest index
1379     # should be specified first before going to the next highest index for
1380     # which unique transitions remain. In practice, that explanation is
1381     # probably pretty useless, so just look at the specific code for the given
1382     # case when making a new definition.
1383     #
1384     set states [dict get $resDef states]
1385     set numSites [llength [dict get $resDef states [lindex $states 0]]]
1386     set numStates [dict size $states]
1387     set maxStates [expr {int(pow(2, $numSites))}]
1388     set numPairs [expr {int($maxStates*($maxStates - 1) / 2)}]
1389     # Are any protonation counts missing?
1390     set protCntExists [lrepeat [expr {$numSites+1}] 0]
1391     dict for {state occ} $states {
1392         lset protCntExists [expr {int(lsum($occ))}] 1
1393     }
1394     set missingProtCnts [list]
1395     set i 0
1396     foreach exists $protCntExists {
1397         if {!$exists} {
1398             lappend missingProtCnts $i
1399         }
1400         incr i
1401     }
1402
1403     set Matrix [lrepeat $numPairs 0.0]
1404     set numValues [llength $data] 
1405     set errorCode 0
1406     switch -- $numSites {
1407         1 {
1408             set errorCode [expr {$numStates == 2 ? $errorCode : -1}]
1409             set errorCode [expr {$numValues == 1 ? $errorCode : -1}]
1410             lassign $data attr10
1411             mset Matrix 1 0 $attr10
1412         } ;# END 1 site, 1 state
1413         2 {
1414             switch -- $numStates {
1415                 3 {
1416                     if {0 in $missingProtCnts} {
1417                         # state (0,0) is deleted
1418                         switch -- $numValues {
1419                             1 { ;# Ex. H2O/H3O+
1420                                 lassign $data attr32
1421                                 mset Matrix 3 2 $attr32
1422                                 mset Matrix 3 1 $attr32
1423                             }
1424                             2 { ;# Ex. HIS
1425                                 lassign $data attr32 attr31
1426                                 mset Matrix 3 2 $attr32
1427                                 mset Matrix 3 1 $attr31
1428                             }
1429                             default {
1430                                 set errorCode -1
1431                             }
1432                         }
1433                         if {$errorCode} break
1434
1435                         mset Matrix 2 1 [expr {[mindex $Matrix 3 1]\
1436                                                - [mindex $Matrix 3 2]}]
1437                         # END 2 sites, 3 states, no state 0 
1438                     } elseif {2 in $missingProtCnts} {
1439                         # state (1,1) is deleted
1440                         switch -- $numValues {
1441                             1 { ;# Ex. carboxylates (ASP, GLU)
1442                                 lassign $data attr20
1443                                 mset Matrix 2 0 $attr20
1444                                 mset Matrix 1 0 $attr20
1445                             }
1446                             2 { ;# Ex. asymmetric carboxylate?
1447                                 lassign $data attr20 attr10
1448                                 mset Matrix 2 0 $attr20
1449                                 mset Matrix 1 0 $attr10
1450                             }
1451                             default {
1452                                 set errorCode -1
1453                             }
1454                         }
1455                         if {$errorCode} break
1456
1457                         mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
1458                                                - [mindex $Matrix 1 0]}]
1459                         # END 2 sites, 3 states, no state 3
1460                     } else {
1461                         # Why would state (0,1) or (1,0) be missing?
1462                         set errorCode -1
1463                     }
1464                 } ;# END 2 sites, 3 states
1465                 4 {
1466                     switch -- $numValues {
1467                         4 {
1468                             lassign $data attr32 attr31 attr20 attr10
1469                             mset Matrix 3 2 $attr32
1470                             mset Matrix 3 1 $attr31
1471                             mset Matrix 2 0 $attr20
1472                             mset Matrix 1 0 $attr10
1473                         }
1474                         default {
1475                             set errorCode -1
1476                         }
1477                     }
1478                     if {$errorCode} break
1479
1480                     mset Matrix 3 0 [expr {[mindex $Matrix 3 1]\
1481                                            + [mindex $Matrix 1 0]}]
1482                     mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
1483                                            - [mindex $Matrix 1 0]}]
1484
1485                     set alt30 [expr {[mindex $Matrix 3 2]\
1486                                      + [mindex $Matrix 2 0]}]
1487                     if {[mindex $Matrix 3 0] != $alt30} {
1488                         set errorCode -2
1489                     }
1490                     # END 2 sites, 4 states, no states missing
1491                 } ;# END 2 sites, 4 states
1492                 default {
1493                     set errorCode -1
1494                 }
1495             }
1496         } ;# END 2 sites
1497         3 {
1498             switch -- $numStates {
1499                 4 {
1500                     if {0 in $missingProtCnts && 1 in $missingProtCnts} {
1501                         # state (0,0,0) and 1 proton states missing
1502                         switch -- $numValues {
1503                             1 { ;# Ex. LYS
1504                                 lassign $data attr76
1505                                 mset Matrix 7 6 $attr76
1506                                 mset Matrix 7 5 $attr76
1507                                 mset Matrix 7 3 $attr76
1508                             }
1509                             default {
1510                                 set errorCode -1
1511                             }
1512                         }
1513                         if {$errorCode} break
1514
1515                         mset Matrix 6 5 [expr {[mindex $Matrix 7 5]\
1516                                                - [mindex $Matrix 7 6]}]
1517                         mset Matrix 6 3 [expr {[mindex $Matrix 7 3]\
1518                                                - [mindex $Matrix 7 6]}]
1519                         mset Matrix 5 3 [expr {[mindex $Matrix 7 3]\
1520                                                - [mindex $Matrix 7 5]}]
1521
1522                         set alt53 [expr {[mindex $Matrix 6 3]\
1523                                          - [mindex $Matrix 6 5]}]
1524                         set alt63 [expr {[mindex $Matrix 5 3]\
1525                                          - [mindex $Matrix 6 5]}]
1526                         set alt65 [expr {[mindex $Matrix 6 3]\
1527                                          - [mindex $Matrix 5 3]}]
1528                         if {[mindex $Matrix 5 3] != $alt53
1529                             || [mindex $Matrix 6 3] != $alt63
1530                             || [mindex $Matrix 6 5] != $alt65} {
1531                             set errorCode -2
1532                         }
1533                         # END 3 sites, 4 states, no state 0, 1, 2, 4
1534                     } elseif {2 in $missingProtCnts && 3 in $missingProtCnts} {
1535                         # state (1, 1, 1) and 2 proton states missing
1536                         switch -- $numValues {
1537                             1 { ;# Ex. phosphate monoesters
1538                                 lassign $data attr40
1539                                 mset Matrix 4 0 $attr40
1540                                 mset Matrix 2 0 $attr40
1541                                 mset Matrix 1 0 $attr40 
1542                             } 
1543                             default {
1544                                 set errorCode -1
1545                             }
1546                         }
1547                         if {$errorCode} break
1548                         
1549                         mset Matrix 4 2 [expr {[mindex $Matrix 4 0]\
1550                                                - [mindex $Matrix 2 0]}]
1551                         mset Matrix 4 1 [expr {[mindex $Matrix 4 0]\
1552                                                - [mindex $Matrix 1 0]}]
1553                         mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
1554                                                - [mindex $Matrix 1 0]}]
1555                         # END 3 sites, 4 states, no state 3, 5, 6, 7
1556                     } else {
1557                         set errorCode -1
1558                     }
1559                 } ;# END 3 sites, 4 states
1560                 7 {
1561                     if {3 in $missingProtCnts} {
1562                         # state (1, 1, 1) missing
1563                         switch -- $numValues {
1564                             2 { ;# Ex. PHOsphate w/o H3PO4
1565                                 lassign $data attr64 attr40
1566                                 mset Matrix 6 4 $attr64
1567                                 mset Matrix 6 2 $attr64
1568                                 mset Matrix 6 1 $attr64
1569                                 mset Matrix 5 4 $attr64
1570                                 mset Matrix 5 2 $attr64
1571                                 mset Matrix 5 1 $attr64
1572                                 mset Matrix 3 4 $attr64
1573                                 mset Matrix 3 2 $attr64
1574                                 mset Matrix 3 1 $attr64
1575
1576                                 mset Matrix 4 0 $attr40
1577                                 mset Matrix 2 0 $attr40
1578                                 mset Matrix 1 0 $attr40
1579                             }
1580                             default {
1581                                 set errorCode -1
1582                             }
1583                         }
1584                         if {$errorCode} break
1585
1586                         mset Matrix 6 0 [expr {[mindex $Matrix 6 4]\
1587                                                + [mindex $Matrix 4 0]}]
1588                         mset Matrix 5 0 [expr {[mindex $Matrix 5 4]\
1589                                                + [mindex $Matrix 4 0]}]
1590                         mset Matrix 3 0 [expr {[mindex $Matrix 3 1]\
1591                                                + [mindex $Matrix 1 0]}]
1592                         mset Matrix 6 5 [expr {[mindex $Matrix 6 4]\
1593                                                - [mindex $Matrix 5 4]}]
1594                         mset Matrix 6 3 [expr {[mindex $Matrix 6 2]\
1595                                                - [mindex $Matrix 3 2]}]
1596                         mset Matrix 5 3 [expr {[mindex $Matrix 5 1]\
1597                                                - [mindex $Matrix 3 1]}]
1598                         mset Matrix 4 2 [expr {[mindex $Matrix 4 0]\
1599                                                - [mindex $Matrix 2 0]}]
1600                         mset Matrix 4 1 [expr {[mindex $Matrix 4 0]\
1601                                                - [mindex $Matrix 1 0]}]
1602                         mset Matrix 2 1 [expr {[mindex $Matrix 2 0]\
1603                                                - [mindex $Matrix 1 0]}]
1604                         # END 3 sites, 7 states, no state 7
1605                     } else {
1606                         set errorCode -1
1607                     }
1608                 } ;# END 3 sites, 7 states
1609                 default {
1610                     set errorCode -1
1611                 }
1612             }
1613         } ;# END 3 sites
1614         default {
1615             set errorCode -1
1616         }
1617     }
1618
1619     # Catch errors here.
1620     switch -- $errorCode {
1621         -1 {
1622             abort "Bad or unimplemented specification of $numValues values for\
1623                     $numSites site residue with $numStates states"
1624         }
1625         -2 {
1626             abort "Error in thermodynamic cycle"
1627         }
1628     }
1629
1630     return $Matrix
1631 }
1632