To make it easier to get started with ChemApp for Python, we developed Equi2Py, a module that bridges the gap between FactSage’s Equilib calculations and modern Python workflows. ChemApp for Python provides the functionality to perform equilibrium calculations in a Python environment, and Equi2Py streamlines the process by converting FactSage Equilib files directly into ready-to-use Python scripts or Jupyter Notebooks.
By automatically generating a fully functional script, Equi2Py removes the need for manual setup and lowers the entry barrier for process modeling and materials informatics. The output is clean, transparent, and easy to adapt, making it straightforward to extend calculations, test variations, or integrate results into broader data pipelines.
Many features of Equilib are supported. The following figure shows the list of supported and not supported features. Note that there are some combinations of features that are not supported.
The ChemApp for Python script can either be in the form of a python file (.py) or in the form of a Jupyter Notebook (.ipynb). The following example shows the resulting Jupyter Notebook that is created when you export the example Equilib file (Ex_CH4-O2) that shows how to calculate the adiabatic flame temperature for methane combustion:
Initialization and general setup¶
This file is automatically generated using equi2py. Created on 2023-09-25 15:35:16
import os
from pathlib import Path
import numpy as np
from collections import OrderedDict
import chemapp as ca
from chemapp.friendly import (
Info,
StreamCalculation as casc,
ThermochemicalSystem as cats,
EquilibriumCalculation as caec,
Units as cau,
TemperatureUnit,
PressureUnit,
VolumeUnit,
AmountUnit,
EnergyUnit,
Status,
)
from chemapp.core import (
LimitVariable,
ChemAppError,
)
Load data file¶
cst_file = Path("pfadtest.cst")
# load the database file
cats.load(cst_file)
Info about the CST file thermodynamic system¶
print(cats.get_cst_info_str())
# system components
print(cats.get_str_scs())
# phases and constituents
print(cats.get_str_phs_pcs())
# only phases: cats.get_str_phs()
Transparent Data File Header Information
========================================
dtc File Creation Timestamp 2023-09-25 15:34:51
dte File Expiration Timestamp 2025-11-01 00:00:00
rem Remarks
id Permitted User IDs 5001
usr Permitted User Names GTT - Technologies
nrp Permitted Reading Program IDs ????
vnr Minimum Reading Program Version -1.-1.-1
ver Header Format Version 1
nwp Source Program Name Equilib
vnw Source Program Version 8.3.0
SYSTEM COMPONENTS
INDEX NAME
0 O
1 C
2 H
PHASES AND CONSTITUENTS
INDEX PHASE NAME STATUS INDEX CONSTITUENT NAME STATUS
0 gas_ideal ENTERED 0 H ENTERED
1 H2 ENTERED
2 C ENTERED
3 C2 ENTERED
4 C3 ENTERED
5 C4 ENTERED
6 C5 ENTERED
7 CH ENTERED
8 CH2 ENTERED
9 CH3 ENTERED
10 CH4 ENTERED
11 C2H ENTERED
12 C2H2 ENTERED
13 C2H3 ENTERED
14 C2H4 ENTERED
15 C2H6 ENTERED
16 O ENTERED
17 O2 ENTERED
18 O3 ENTERED
19 OH ENTERED
20 H2O ENTERED
21 HOO ENTERED
22 HOOH ENTERED
23 CO ENTERED
24 C2O ENTERED
25 CO2 ENTERED
26 C3O2 ENTERED
27 HCO ENTERED
28 H2CO ENTERED
29 CH3OH ENTERED
30 CH2CO ENTERED
31 C2H4O(g) ENTERED
32 C2H4O(g2) ENTERED
33 CH3CH2OH(g) ENTERED
34 CH3CH2OH(g2) ENTERED
35 HCOOH ENTERED
36 CH3COOH ENTERED
1 C_Graphite(s) ENTERED
2 C_diamond(s2) ENTERED
3 H2O_Ice(s) ENTERED
Units setup¶
It is possible to change the units before retrieving results, or in between making multiple inputs. However, mixed unit input is currently not generated automatically.
cau.set(
T = TemperatureUnit.C,
P = PressureUnit.atm,
V = VolumeUnit.dm3,
A = AmountUnit.mol,
E = EnergyUnit.J
)
Fixed conditions¶
casc.set_eq_P(1.0)
casc.set_eq_H(0.0)
Stream setup¶
Create streams. The names are defaults of “#1”, “#2”, … Feel free to make these more reasonable, but take care to edit everywhere they are used
casc.create_st(name="#1", T=25, P=1)
Variable conditions¶
define variables which are supposedly modified in iterations.
A = 0.0
# original condition for A was 1.0, which would not include 1.0.
# Increasing it by 0.05 yields the inclusive upper boundary.
A_range = np.arange(0.0, 1.05, 0.05)
# we keep the results in an OrderedDict, to access by variable value
results = OrderedDict()
Calculation¶
# iter_ic_tg_1var
for A in A_range:
# set incoming amounts
casc.set_IA_pc("#1", "gas_ideal", "CH4", 1-A )
casc.set_IA_pc("#1", "gas_ideal", "O2", A )
try:
# calculate
casc.calculate_eq_T(T_guess=1000.0)
#fetch results
results[A] = casc.get_result_object()
except ChemAppError as c_err:
print(f"[<A> = {A}] ChemApp Error {c_err.errno}:")
print(c_err)
# change this ...!
pass
[<A> = 0.0] ChemApp Error 712: Target calculation aborted; the value of the target variable is less than th lowest permitted The error appeared after a call of "TQCE" [<A> = 1.0] ChemApp Error 711: Target calculation aborted; no solution is found within the permitted interval The error appeared after a call of "TQCE"
Postprocessing¶
import pandas as pd
Temperatures along the axis of calculation¶
TvsA = pd.DataFrame(
data=[(A, x.T) for A, x in results.items()],
columns=["A", "T"],
).set_index("A")
# For more elaborate plotting, use e.g. matplotlib.
# If you work with a Jupyter notebooks, you can remove the part
# .figure.savefig("TvsA.png")
# to avoid saving the figure.
TvsA.plot(grid=True, xlabel='Alpha', ylabel=(f'T in {cau.get_T_unit().name}')).figure.savefig("TvsA.png")
Oxygen activity over the iterating variable¶
import numpy as np
# Check if oxygen is in the dataset
if ('gas_ideal', 'O2') in cats.get_name_ph_pcs():
ActivityO2vsA = pd.DataFrame(
data=[(A, r.phs["gas_ideal"].pcs["O2"].AC) for A, r in results.items()],
columns=["A", "Oxygen activity"],
).set_index("A")
ActivityO2vsA["Oxygen activity (log10)"] = ActivityO2vsA["Oxygen activity"].apply(np.log10)
# For more elaborate plotting, use e.g. matplotlib.
# If you work with a Jupyter notebooks, you can remove the part
# .figure.savefig("ActivityO2vsA.png")
# to avoid saving the figure.
ActivityO2vsA["Oxygen activity (log10)"].plot(
grid=True,
xlabel='Alpha',
ylabel=f'log(p(O2)), p(O2) in {cau.get_P_unit().name}',
).figure.savefig("ActivityO2vsA.png")
Phase constituent fractions¶
# Setting a threshold to filter
threshold = 0.001
# change this line for the phase you want
# we simply pick index 0 (typically the gas phase)
# you can also just set the name, e.g.
# phase = "gas_real"
phase = cats.get_name_phs()[0]
# Set up DataFrame, empty
constituent_amounts = pd.DataFrame(
columns=["A"] + cats.get_name_pcs_in_ph(phase),
).set_index("A")
# Iterate over each result object
for A, res in results.items():
# go through each phase constituent object
for pcname, pc in res.phs[phase].pcs.items():
# check for the threshold value
if pc.A > threshold:
# Add the datapoint value -> .A is property of PhaseConstituentState class
constituent_amounts.loc[A, pcname] = pc.A
# Plot the DataFrame (after dropping those never above the threshold)
constituent_amounts.dropna(how="all", axis=1).plot(grid=True)
<Axes: xlabel='A'>
In case of questions, please contact us via chemapp@factsage.com!
Interested in news about ChemApp for Python (new versions with new features, training courses)? Please register for our newsletter here!