Oxidation of Pyrrole to Polypyrrole (Ammonium Persulphate; SDBS)

pyrrole_conversion.py

import argparse

from matplotlib import pyplot, ticker
from scipy import integrate


def f(time, concs, k1, k2, chain_length):
    dPY = dPYdt(concs, k1, k2)
    return [dPY] + [
        x(dPY, chain_length) for x in [dAMPSdt, dSDBSdt, dPPYdt, dAMSdt, dSAdt]
    ]


def dPYdt(concs, k1, k2):
    """ Rate of change of pyrrole concentration with respect to time. """
    return -k1 * concs[1] * concs[0] - k2 * concs[0] * concs[3]


def dAMPSdt(pyrrole_change, chain_length):
    """ Rate of change of ammonium persulphate concentration with respect to time. """
    return (chain_length - 1) / chain_length * (pyrrole_change)


def dSDBSdt(pyrrole_change, chain_length):
    """ Rate of change of SDBS concentration with respect to time. """
    return 1 / chain_length * (pyrrole_change)


def dPPYdt(pyrrole_change, chain_length):
    """ Rate of change of PPY concentration with respect to time. """
    return -1 / chain_length * (pyrrole_change)


def dAMSdt(pyrrole_change, chain_length):
    """ Rate of change of ammonium sulphate concentration with respect to time. """
    return -(chain_length - 1) / chain_length * (pyrrole_change)


def dSAdt(pyrrole_change, chain_length):
    """ Rate of change of sulphuric acid concentration with respect to time. """
    return -(chain_length - 1) / chain_length * (pyrrole_change)


if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description="Polymerisation of Pyrrole (Ammonium Persulphate/SDBS)"
    )
    parser.add_argument(
        "-k1",
        help="Rate constant with respect to the oxidant and monomer [M/min]",
        default=0.355,
        type=float,
    )
    parser.add_argument(
        "-k2",
        help="Rate constant with respect to the polymer and monomer [M/min]",
        default=0.965,
        type=float,
    )
    parser.add_argument(
        "-c",
        "-chain-length",
        help="Number of repeating monomer units",
        default=30,
        type=int,
    )
    parser.add_argument(
        "-t", "-time", help="Reaction time [min]", default=90, type=float
    )
    parser.add_argument(
        "-s", "-step", help="Integration step size [min]", default=0.1, type=float
    )
    parser.add_argument(
        "-pyrrole", help="Initial pyrrole concentration [M]", default=0.1413, type=float
    )
    parser.add_argument(
        "-ammonium-persulphate",
        help="Initial ammonium persulphate concentration [M]",
        default=0.1365,
        type=float,
    )
    parser.add_argument(
        "-sdbs", help="Initial sdbs concentration [M]", default=0.0040, type=float
    )
    parser.add_argument(
        "-polypyrrole",
        help="Initial polypyrrole concentration [M]",
        default=0.0,
        type=float,
    )
    parser.add_argument(
        "-ammonium-sulphate",
        help="Initial ammonium sulphate concentration [M]",
        default=0.0,
        type=float,
    )
    parser.add_argument(
        "-sulphuric-acid",
        help="Initial sulphuric acid concentration [M]",
        default=0.0,
        type=float,
    )
    args = vars(parser.parse_args())

    labels = [
        "Pyrrole",
        "Ammonium Persulphate",
        "SDBS",
        "Polypyrrole",
        "Ammonium Sulphate",
        "Sulphuric Acid",
    ]
    concs = [args[l.lower().replace(" ", "_")] for l in labels]
    k1, k2, chain_length, time, step = (args[i] for i in ["k1", "k2", "c", "t", "s"])
    res = integrate.solve_ivp(
        fun=f,
        t_span=[0, time],
        y0=concs,
        max_step=step,
        args=(k1, k2, chain_length),
    )

    figure = pyplot.figure()

    conversion_data = [1 - c / concs[0] for c in res.y[0]]
    conversion = figure.add_subplot(1, 2, 1)
    conversion.plot(res.t, conversion_data)
    conversion.set_xlabel("Time (min)")
    conversion.set_ylabel("Conversion (%)")
    conversion.yaxis.set_major_formatter(ticker.PercentFormatter(1.0))
    conversion.set_title(f"Pyrrole Conversion (Units: {chain_length})")
    conversion.grid()
    conversion.set_ylim(0, 1)

    concentrations = figure.add_subplot(1, 2, 2)
    for i, label in enumerate(labels):
        concentrations.plot(res.t, res.y[i], label=label)
    concentrations.set_xlabel("Time (min)")
    concentrations.set_ylabel("Concentration (M)")
    concentrations.legend()
    concentrations.set_title("Changes in Concentration with Time")
    concentrations.grid()
    pyplot.show()