{ "cells": [ { "cell_type": "markdown", "id": "b640c74d", "metadata": {}, "source": [ "# Performing calculation of reaction path using Reaktoro\n", "\n", "This tutorial demonstrates how to calculate a reaction path between two different chemical states in\n", "equilibrium referred to as the *initial state* and *final state*.\n", "These states can have different temperatures, pressures, and/or molar amounts of elements. If we gradually adjust\n", "temperature, pressure, and elemental amounts to bring the initial state to the final state, slowly\n", "enough so that **every intermediate state is in equilibrium**, the system would trace a co-called *reaction path*.\n", "\n", "Let the initial state have 1 g of calcite (CaCO3) mixed with 1 kg of water. We want to see how the\n", "addition of\n", "hydrochloric acid (HCl), up to 1 mmol, contributes to the dissolution of calcite. Thus, our initial and final\n", "states for a reaction path calculation can be described as follows:\n", "\n", "| Initial state | Final state |\n", "|----------------|----------------|\n", "| 1 kg of H2O | 1 kg of H2O |\n", "| 1 g of CaCO3 | 1 g of CaCO3 |\n", "|

| 1 mmol of HCl |\n", "\n", "As usual, we start by importing the `reaktoro` package:" ] }, { "cell_type": "code", "execution_count": null, "id": "560ee179", "metadata": {}, "outputs": [], "source": [ "from reaktoro import *" ] }, { "cell_type": "markdown", "id": "b540f550", "metadata": {}, "source": [ "Note that the object `editor` from class\n", "[ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) was not initialized with a given\n", "[Database](https://reaktoro.org/cpp/classReaktoro_1_1Database.html) object. Instead, it is initialized using the\n", "default built-in database file\n", "[supcrt98.xml](https://github.com/reaktoro/reaktoro/blob/master/databases/supcrt/supcrt98.xml)." ] }, { "cell_type": "code", "execution_count": null, "id": "6451c195", "metadata": {}, "outputs": [], "source": [ "editor = ChemicalEditor()" ] }, { "cell_type": "markdown", "id": "cfd53325", "metadata": {}, "source": [ "For the aqueous phases, we list the chemical elements composing the phase instead of the species' exact names.\n", "Class [ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) searches for all species in\n", "the database that those elements can form. Only species corresponding to the phase-type are selected\n", "(e.g., only aqueous species are searched in the current case)." ] }, { "cell_type": "code", "execution_count": null, "id": "88a62b24", "metadata": {}, "outputs": [], "source": [ "editor.addAqueousPhaseWithElements(\"H O Ca C Cl\")\n", "editor.addMineralPhase(\"Calcite\")\n", "\n", "system = ChemicalSystem(editor)" ] }, { "cell_type": "markdown", "id": "a6dc48c1", "metadata": {}, "source": [ "In the code below, two instances of the class\n", "[EquilibriumProblem](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html) are created:\n", "`initial_problem` describes the initial state, and `final_problem` corresponds to the final state." ] }, { "cell_type": "code", "execution_count": null, "id": "16e7ab8d", "metadata": {}, "outputs": [], "source": [ "initial_problem = EquilibriumProblem(system)\n", "initial_problem.setTemperature(30.0, \"celsius\")\n", "initial_problem.setPressure(1.0, \"bar\")\n", "initial_problem.add(\"H2O\", 1, \"kg\")\n", "initial_problem.add(\"CaCO3\", 1, \"g\")\n", "\n", "final_problem = EquilibriumProblem(system)\n", "final_problem.setTemperature(30.0, \"celsius\")\n", "final_problem.setPressure(1.0, \"bar\")\n", "final_problem.add(\"H2O\", 1, \"kg\")\n", "final_problem.add(\"CaCO3\", 1, \"g\")\n", "final_problem.add(\"HCl\", 1, \"mmol\")" ] }, { "cell_type": "markdown", "id": "50570523", "metadata": {}, "source": [ "Two instances of the class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) are created\n", "to store the initial and final equilibrium states calculated by the method\n", "[equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e)." ] }, { "cell_type": "code", "execution_count": null, "id": "6c38afc7", "metadata": {}, "outputs": [], "source": [ "initial_state = equilibrate(initial_problem)\n", "final_state = equilibrate(final_problem)" ] }, { "cell_type": "markdown", "id": "692181ad", "metadata": {}, "source": [ "Once the initial and final equilibrium states have been calculated, it is now time to trace the reaction path\n", "between them, with each intermediate state in chemical equilibrium. For this, we use the class\n", "[EquilibriumPath](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumPath.html). Note that its initialization\n", "requires a [ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) instance:" ] }, { "cell_type": "code", "execution_count": null, "id": "bf4161f4", "metadata": {}, "outputs": [], "source": [ "path = EquilibriumPath(system)" ] }, { "cell_type": "markdown", "id": "79824c29", "metadata": {}, "source": [ "Before calling the method\n", "[EquilibriumPath::solve](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumPath.html#a008b74301618ed186caa95ec059eb204)\n", ", one can configure output-file to be generated during the calculation.\n", "To output quantities to a file or terminal during the calculation, use method\n", "[EquilibriumPath::output](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumPath.html#ac700ac6f939acbd4a8524e1346a1e588),\n", "which returns an instance of class [ChemicalOutput](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalOutput.html):" ] }, { "cell_type": "code", "execution_count": null, "id": "af4da9ca", "metadata": {}, "outputs": [], "source": [ "output = path.output()\n", "output.filename(\"result.txt\")\n", "output.add(\"speciesAmount(Cl- units=mmol)\", \"Cl- [mmol]\")\n", "output.add(\"speciesMolality(CO2(aq) units=mmolal)\", \"CO2(aq) [mmol]\")\n", "output.add(\"speciesMolality(CO3-- units=mmolal)\", \"CO3-- [mmol]\")\n", "output.add(\"speciestMolality(Ca++ units=mmolal)\", \"Ca++ [mmolal]\")\n", "output.add(\"pH\")\n", "output.add(\"speciesMass(Calcite units=g)\")" ] }, { "cell_type": "markdown", "id": "c6cd3950", "metadata": {}, "source": [ "The method [ChemicalOutput::add](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalOutput.html#af3b5a7d6b0fbbc870664d6ad100b10dd)\n", "adds a quantity, which we want to be output to the file `result.txt`. The latter filename is specified in the call of\n", "the method [ChemicalOutput::filename](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalOutput.html#ac5cc9d0f90cfe5c6e0972a55b7f7bf5d).\n", "Each call to [ChemicalOutput::add](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalOutput.html#af3b5a7d6b0fbbc870664d6ad100b10dd)\n", "results in a new column of data in the output file.\n", "\n", "> **Note**: When two arguments are provided to the method\n", "[ChemicalOutput::add](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalOutput.html#a54b0e4fd28823c4d1d1884c32eed1cf3),\n", "the first one is the name of the quantity to be output (e.g.,\n", "`time`, `elementAmount(Cl)`, `ionicStrength`). The second one is a label used as the heading of the column of data\n", "in the output file. When only one argument is provided, this single argument is both the label and the quantity name.\n", "\n", "Finally, after output files have been configured, the equilibrium path can be calculated using:" ] }, { "cell_type": "code", "execution_count": null, "id": "238c14b5", "metadata": {}, "outputs": [], "source": [ "path = path.solve(initial_state, final_state)" ] }, { "cell_type": "markdown", "id": "8271e27b", "metadata": {}, "source": [ "### Plotting the results of equilibrium path calculation\n", "\n", "We now use [bokeh](https://docs.bokeh.org/en/latest/docs/gallery.html#standalone-examples)\n", "to do the plotting." ] }, { "cell_type": "code", "execution_count": null, "id": "2a75913e", "metadata": { "lines_to_next_cell": 1 }, "outputs": [], "source": [ "from bokeh.plotting import figure, show\n", "from bokeh.io import output_notebook\n", "output_notebook()" ] }, { "cell_type": "markdown", "id": "3d017f27", "metadata": {}, "source": [ "Besides, we define a custom function that would generate figure of a size 600 x 300:" ] }, { "cell_type": "code", "execution_count": null, "id": "0973b962", "metadata": {}, "outputs": [], "source": [ "def custom_figure(x_axis_label, y_axis_label):\n", " return figure(plot_width=600, plot_height=300,\n", " x_axis_label=x_axis_label,\n", " y_axis_label=y_axis_label)" ] }, { "cell_type": "markdown", "id": "42b4ce12", "metadata": {}, "source": [ "To load results from the outputfile, we use `loadtxt` function provided by the `numpy` package:" ] }, { "cell_type": "code", "execution_count": null, "id": "048a4546", "metadata": {}, "outputs": [], "source": [ "filearray = numpy.loadtxt(\"result.txt\", skiprows=1)\n", "data = filearray.T\n", "[cl_indx, co2aq_indx, co3_indx, ca_indx, ph_indx, calcite_indx] = numpy.arange(6)" ] }, { "cell_type": "markdown", "id": "530ea100", "metadata": {}, "source": [ "The first plot depicts the change of pH with the addition of HCl into the system:" ] }, { "cell_type": "code", "execution_count": null, "id": "13ef7861", "metadata": {}, "outputs": [], "source": [ "fig1 = custom_figure(x_axis_label=\"HCl [mmol]\", y_axis_label=\"pH [-]\")\n", "fig1.line(data[cl_indx], data[ph_indx], line_width=4, color=\"darkviolet\")\n", "show(fig1)" ] }, { "cell_type": "markdown", "id": "b3750d56", "metadata": {}, "source": [ "Adding HCl to the initial state contributes the same amount of Cl-. Therefore, we can study the\n", "concentrations of aqueous species Ca2-. Its growth with respect to the growing amount of Cl-\n", "indicates that more calcite we can get dissolved in the brine. Or in the other words, the solubility of calcite grows\n", "with the addition of HCl." ] }, { "cell_type": "code", "execution_count": null, "id": "10321da6", "metadata": {}, "outputs": [], "source": [ "fig2 = custom_figure(x_axis_label=\"Amount of Cl- [mmol]\", y_axis_label=\"Concentration of Ca++ [mmolal]\")\n", "fig2.line(data[cl_indx], data[ca_indx], line_width=4, color=\"darkgreen\")\n", "show(fig2)" ] }, { "cell_type": "markdown", "id": "0b551392", "metadata": {}, "source": [ "Simultaneously, we plot the concentration of species CO32-, which suppose to grow,\n", "while calcite is dissolving." ] }, { "cell_type": "code", "execution_count": null, "id": "41e6d765", "metadata": {}, "outputs": [], "source": [ "fig3 = custom_figure(x_axis_label=\"pH\", y_axis_label=\"Concentration of CO3--[mmolal]\")\n", "fig3.line(data[ph_indx], data[co3_indx], line_width=4, color=\"orange\")\n", "fig3.legend.location = \"top_left\"\n", "show(fig3)" ] }, { "cell_type": "markdown", "id": "af4847ab", "metadata": {}, "source": [ "The fourth and last figure plots how the mass of calcite (or calcium carbonate) changes with the addition of\n", "HCl in the system. We see that the remaining mass of the mineral after equilibration goes down since its solubility\n", "increases with added hydrogen chloride and increasing pH." ] }, { "cell_type": "code", "execution_count": null, "id": "41c9d9d1", "metadata": {}, "outputs": [], "source": [ "fig4 = custom_figure(x_axis_label=\"HCl [mmol]\", y_axis_label=\"Mass Calcite [g]\")\n", "fig4.line(data[cl_indx], data[calcite_indx], line_width=4, color=\"darkviolet\")\n", "show(fig4)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 5 }