Equilibrium with some concentrations fixed¶

Sometimes we may have a chemical system in which one or more species has their concentration held constant. For example, we could have a system where a chemical species is constantly fed in such that its concentration in the reaction vessel remains constant. As another example, we might have a weak acid in solution and we add a strong acid to control the pH, thereby keeping the concentration of hydrogen ions constant.

The eqtk.final_value_solve() function allows for solving the coupled equilibrium problem when one or more of of chemical species are present in fixed amounts. Its API, as we will demonstrate, is the same as eqtk.solve(), except an additional input, fixed_c, which specifies the fixed concentration is required. However, eqtk.final_value_solve() only accepts input of the stoichiometric matrix and equilibrium constants; input of free energies and a conservation matrix are not yet implemented.

In what follows, we will assume that Numpy, Pandas, and EQTK have been imported respectively as np, pd, and eqtk.

Example problem¶

To demonstrate the use of eqtk.fixed_value_solve(), we will consider the protonation state of aspartic acid. The reactions describing deprotonation are

H₃D⁺ ⇌ H₂D + H⁺

H₂D ⇌ HD⁻ + H⁺

HD⁻ ⇌ D²⁻ + H⁺

The respective pKₐ’s for the reaction are 1.99, 3.9, and 10.002. We additionally need to take into account the dissociation of water,

H₂O ⇌ OH⁻ + H⁺,

which has an equilibrium constant of 10⁻¹⁴ M².

Our goal is to compute the relative abundances of the four protonation states of aspartic acid as a function of pH. Thus, for each pH value, we want to fix [H⁺] = 10⁻ᵖᴴ and compute the equilibrium concentrations of the other species.

The chemical reactions¶

We can conveniently specify the chemical reactions and species using EQTK’s parser for string representations. Recall that solvent is not explicitly included in the chemical reactions.

rxns = """
     <=> OH⁻ + H⁺  ; 1e-14
H₃D⁺ <=> H₂D + H⁺  ; {Ka1}
H₂D  <=> HD⁻ + H⁺  ; {Ka2}
HD⁻  <=> D²⁻ + H⁺  ; {Ka3}
""".format(Ka1=10**(-1.99), Ka2=10**(-3.9), Ka3=10**(-10.002))

N = eqtk.parse_rxns(rxns)

Initial and fixed concentrations¶

The initial concentrations of all species, \(\mathbf{c}^0\), is specified the same as for eqtk.solve(). For this calculation, we will begin with 1 mM H₂D.

c0 = {
    "H⁺":   0,
    "OH⁻":  0,
    "H₃D⁺": 0,
    "H₂D":  0.001,
    "HD⁻":  0,
    "D²⁻":  0,
}

While it is convenient to specify the concentrations of the species as a dictionary, we could also use a Pandas Series or DataFrame. (We could also specify N, K, and c0 as Numpy arrays.) The specification of inputs is the same as for eqtk.solve().

We need to further specify the fixed concentrations. We will consider pH ranging from zero to 14, with 400 points in between. To set up the fixed concentrations, we create a DataFrame with six columns, one for each species, with 400 rows, one for each pH value we want to consider. Entries in the data frame that are either negative or np.nan denote concentrations that are not fixed. So, we initialize the data frame with negative ones. Following that, we fill the column for H⁺ with our desired fixed values.

fixed_c = pd.DataFrame(data=-np.ones((400, 6)), columns=c0.keys())

pH = np.linspace(0, 14, 400)
fixed_c["H⁺"] = 10**(-pH)

The resulting fixed_c DataFrame has the following first five rows.

         H⁺  OH⁻  H₃D⁺  H₂D  HD⁻  D²⁻
0.100000 -1.0  -1.0 -1.0 -1.0 -1.0
0.082386 -1.0  -1.0 -1.0 -1.0 -1.0
0.066777 -1.0  -1.0 -1.0 -1.0 -1.0
0.053177 -1.0  -1.0 -1.0 -1.0 -1.0
0.041543 -1.0  -1.0 -1.0 -1.0 -1.0

Solving for the relative abundance¶

We call eqtk.fixed_value_solve() similarly to eqtk.solve(), except with the additional argument fixed_c.

c = eqtk.fixed_value_solve(c0=c0, fixed_c=fixed_c, N=N, units='M')

Because we supplied our inputs as data frames, the output is also a data frame. The output is as for eqtk.solve(), except there are additional columns with names like '[H⁺]__fixed (M)', denoting concentrations that were fixed in the calculation.

A plot of the relative abundances computed from c is shown below.

Over-constraining a problem¶

In a given calculation, only some species may have fixed concentration, lest the coupled equilibrium problem be over-constrained. EQTK will check for this and raise an exception if the problem becomes over-constrained. For example, we cannot fix both the H⁺ and OH⁻ concentrations.

fixed_c = {
    "H⁺":   0.01,
    "OH⁻":  0.0001,
    "H₃D⁺": -1,
    "H₂D":  -1,
    "HD⁻":  -1,
    "D²⁻":  -1,
}

c = eqtk.fixed_value_solve(c0=c0, fixed_c=fixed_c, N=N, units='M')

This results in a ValueError saying that the stoichiometric matrix is rank deficient.