Feature: Enable solid-liquid equilibrium via phreeqpython

[rkingsbury]

phreeqpython contains an equalize method demonstrated in this tutorial that exposes PHREEQC's solid-liquid equilibrium calculations.

A method to access this could be added to Solution. The best way to do so is probably as a keyword argument to equilibrate. For example, one could write

from pyEQL import Solution
s1 = Solution({"Ca+2": "10 mM", "Cl-": "20 mM"})
s1.equilibrate(solids=['Calcite'])

Using the same paradigm, it would be easy to add equilibration with atmospheric gases as another kwarg, e.g.

s1.equilibrate(solids=['Calcite'], gases=['CO2'])

Inside equlibrate, the contents of solids and gases would simply be combined into a single list and passed to equalize of the underlying phreeqpython object (accessible via the s1.engine.pp attribute)

[DhruvDuseja]

Hi @rkingsbury! Hope you are doing well. I am interested in making this enhancement.

If I've understood this correctly, the equilibrate method in the Solution class would be modified to be something like below?

def equilibrate(self, solids=[None], gases=[None], **kwargs) -> None:
    """
    Update the composition of the Solution using the thermodynamic engine.

    Any kwargs specified are passed through to self.engine.equilibrate()

    Returns:
        Nothing. The .components attribute of the Solution is updated.
    """
    phases = solids + gases
    self.engine.pp.equalize(phases)
    self.engine.equilibrate(self, **kwargs)

[rkingsbury]

Hi @rkingsbury! Hope you are doing well. I am interested in making this enhancement.
If I've understood this correctly, the equilibrate method in the Solution class would be modified to be something like below?

Hi @DhruvDuseja , yes, something along those lines is exactly what I had in mind. Just a few initial comments:

• In general when writing python functions, you should not use a list as a default value (because lists are "mutable"), so instead of [None] you'd want to use None for the defaults
• It might be nice to add a convenience kwarg atmosphere: bool = False that, if set to True, automatically equilibrates the solution with atmospheric gases (i.e., CO2 and O2). This would keep a user from having to look up the correct partial pressures of these gases. If we implement this, you'd want to put in a check so that atmosphere=True and gases cannot both be supplied at the same time (i.e., gases will overrride atmosphere)

[DhruvDuseja]

Thanks for the feedback.

So, if atmospheric is set to True, it would set the gases list to [CO2, O2] (if the user passes in no gases) and then combine that with the solids, yes?
And if it's False, try to combine the lists.

Is there a possibility that either of the lists is empty and we pass it to the equalize anyway?

[rkingsbury]

Note on call signature

I realize that the syntax for gases (and maybe solids) needs to be slightly different than what I had brainstormed above. None of this is very well documented, so let me walk you through it. If you examine this example from phreeqpython

solution2 = pp.add_solution({}).equalize(['Calcite', 'CO2(g)'], [0, -3.5])

You see that the first argument to equalize is a list of phases, and the second argument is a list of saturation indices, which are defined as follows:

saturation index --Target saturation index for the pure phase in the aqueous phase (line 1a); for gases, this number is the log of the partial pressure (line 1b). The target saturation index (partial pressure) may not be attained if the amount of the phase in the assemblage is insufficient. Default is 0.0.

I figured this out be examining the definition of equalize, which calls equalize_solution, which sets the EQUILIBRIUM PHASES block of a PHREEQC input file.

So for the pyEQL implementation, we'll have to decide how to specify the saturation indices. I think the best way would be to use dict instead of list for both solids and gases, where the keys are phases and the values are the saturation index or log partial pressure, respectively:

s1.equilibrate(solids={'Calcite': 0}, gases={'CO2':-3.5})

This format would parallel the way that solutes are entered (e.g., {'Na+': '1 mol/L'}). An alternative would be to just copy the way it's done in phreeqpython, where you have one list for all phases (not separating liquids and solids), and one list of saturation indices, but I think that's a little less clear.

Your Questions

So, if atmospheric is set to True, it would set the gases list to [CO2, O2] (if the user passes in no gases) and then combine that with the solids, yes? And if it's False, try to combine the lists.

I was thinking the following: If atmospheric is True, then atmospheric gases and partial pressures ({"CO2": -3.5, "O2": -0.6778}) would be added to whatever else the user has supplied in gases. If the user has supplied CO2 or O2, their input should be overwritten with the values listed here, and a warning should be raised.

Is there a possibility that either of the lists is empty and we pass it to the equalize anyway?

Yes, by default (when both lists are set to None), then you would not pass anything to equalize.

[DhruvDuseja]

So, this is what I have so far based on what I've understood. Can you please take a look and see if I'm headed in the right direction? Thanks!

    def equilibrate(
        self, solids: Optional[dict] = None, gases: Optional[dict] = None, atmospheric: bool = False, **kwargs
    ) -> None:
 
        atmospheric_dict = {"CO2": -3.5, "O2": -0.6778}
        phases_list = []
        si_list = []
        if atmospheric:
            if any(gas in gases for gas in ["CO2", "O2"]):
                raise Warning("Updating CO2 & O2 partial pressures with default values.")
            gases.update(atmospheric_dict)
        phases = {**solids, **gases}
        for phase in phases:
            phases_list.append(phase)
            si_list.append(phases[phase])
        if phases_list and si_list:
            self.engine.pp.equalize(phases_list, si_list)
        self.engine.equilibrate(self, **kwargs)

[rkingsbury]

Thanks so much @DhruvDuseja . Yes, definitely the right direction. There are some small edits you can make to streamline the code (which I can comment on directly when you open a pull request), but overall looks good. For example, instead of the for loop to populate phases_list and si_list, you could use list comprehensions, e.g.:

phases_list = [k for k,,v in phases.items()]
si_list = [v for k,v in phases.items()]

Also, you don't need to call self.engine.equlibrate at the end, because the code you've written above would actually be placed inside self.engine.equilibrate (as additional, optional kwargs).