Use cases for pyrrole

pyrrole_use_cases.py

#!/usr/bin/python3

import pandas as pd
import pyrrole as pyr
from matplotlib import pyplot as plt

# if you want a single molecule

# pyr.Atoms(series)
# pyr.Atoms(dataframe["A"])
# pyr.atoms.from_cclib(ccopen("A.out"))  # Atoms.data index will be jobfilename

# if you want all your molecules

# df = ccframe(...)
# molecules = pyr.atoms.from_dataframe(dataframe)  # returns a list of pyr.Atoms, one for each row in df

# if you want an equation

# df = pd.DataFrame([m.to_series() for m in molecules])
# equation = pyr.ChemicalEquation("A -> B", data=dataframe)
# equation.calculate("freeenergy")  # things are always calculated on demand such that equations are always stateless

# if you want a system

# system = pyr.ChemicalSystem(..., data=dataframe)

def main():
    system = pyr.ChemicalSystem(["BH <- AH <=> A- + H+",
                                 "A- -> B-",
                                 "BH <=> B- + H+"])

    # $ cat reactions.txt
    # BH <- AH <=> A- + H+
    # A- -> B-
    # BH <=> B- + H+

    # with open("reactions.txt") as f:
    #     system = pyr.ChemicalSystem(f.readlines())  # if you have a file for the model

    # $ cat data.tsv
    # label freeenergy
    # A-    0.0
    # B-    1.0
    # AH    2.0
    # BH    3.0
    # H+    4.0

    # internal data is almost aways pd.Series or pd.DataFrame and we say shamelessly so:
    # pass things that can be at least converted!
    
    df = pd.read_table("data.tsv").set_index("label")
    
    # if species is a dict of labels to paths, i.e., {"A": "A.out", ...}
    
    # from cclib import ccopen
    # df = pd.DataFrame({label: pyr.atoms.from_cclib(ccopen(path)).to_series()
    #                    for label, path in species.items()})

    # in the code, any attribute named "data" will be a species (row) by property (column) dataframe

    # only accept dataframes, can also be passed to constructor as ChemicalSystem(..., data=df),
    # delegates data to each ChemicalEquation.data
    system.data = df
    
    # https://en.wikipedia.org/wiki/Energy_profile_(chemistry)
    # calls system.all_simple_paths(source, target), which in turn uses networkx,
    # and assumes energies are as in cclib ("freeenergy", ..., and in hartree)
    system.energy_profile("AH", "BH")
    plt.show()

    print(system.to_latex())  # calls system.to_dataframe(), which calls ChemicalEquation.to_series() to each equation

    # x, y = system.spectrum("ir")


if __name__ == "__main__":
    main()