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