Python code to create random kinetic mechanisms

File contents

from sympy import var
from random import random
from KinMech import *

def randomrxn(reagentlist, maxreactantstoich,
              maxproductstoich):  # generates a random reaction from reagentlist, where the order of the reaction
    # is controlled by maxreactantstoich and maxproductstoich.
    rxn = []
    reagent = []
    phaselist = ['l', 's', 'g', 'aq']
    puritylist = ['solution', 'solvent', 'pure']
    maxnumberofreactants = min(int(round(1 + (maxreactantstoich - 1) * random())), (len(reagentlist) - 1))
    totalstoich = 0
    # print ('maxnumberofreactants: '+str(maxnumberofreactants))
    for i in range(0, (maxnumberofreactants)):
        stoich = int(round(1 + (maxreactantstoich - 1) * random()))  # 1..maxreactantstoich
        # print ('reactants: '+str(i)+', stoich: '+str(stoich))
        if (maxreactantstoich - stoich - totalstoich) < 0:
            stoich = max(0, (maxreactantstoich - totalstoich))
        if stoich > 0:
            phase = phaselist[int(round((len(phaselist) - 1) * random()))]
            purity = puritylist[int(round((len(puritylist) - 1) * random()))]
            if phase == 'aq':
                purity = 'solution'
            if phase == 's':
                purity = 'pure'
            reagent = [stoich, reagentlist[i], 'reactant', phase, purity]
            totalstoich = totalstoich + stoich
            productlist = []
            for j in range(0, (len(reagentlist))):
                productlist.append(reagentlist[j])
            # print i
            for j in range(0, (i + 1)):
                if j <= i:
                    del productlist[0]
                # print productlist
            rxn.append(reagent)
    # print ('reagentlist: '+str(reagentlist)+ ' productlist: '+str(productlist))
    totalstoich = 0
    maxnumberofproducts = min(int(round(1 + (maxproductstoich - 1) * random())), (len(productlist)))
    for i in range(0, (maxnumberofproducts)):
        stoich = int(1 + (maxproductstoich - 1) * random())  # 1..maxproductstoich
        if (maxproductstoich - stoich - totalstoich) < 0:
            stoich = max(0, (maxproductstoich - totalstoich))
        if stoich > 0:
            phase = phaselist[int(round((len(phaselist) - 1) * random()))]
            purity = puritylist[int(round((len(puritylist) - 1) * random()))]
            if phase == 'aq':
                purity = 'solution'
            if phase == 's':
                purity = 'pure'
            if purity == 'solvent':
                purity = 'solution'
            reagent = [stoich, productlist[i], 'product', phase, purity]
            totalstoich = totalstoich + stoich
            rxn.append(reagent)
    return (rxn)


def randommech(maxrxn,
               maxreagents):  # generates a random mechanism with at most maxrxn reversible reactions (2*maxrxn reactions).
    # the mechanism will have at most maxreagents different species in it. Maxreagents > maxrxn.
    # This will occasionally generate mechanisms which are unreasonable, but not often.  Even the
    # unreasonable ones are good for practice with mechanisms.
    maxreagents = max(maxreagents, (maxrxn + 1))  # insisting that the number of reagents is high enough
    # print 'maxreagents: '+str(maxreagents)
    mech = []
    possiblereagents = 'ABCDEFGHJKLMNPRSTUVWXYZ'
    reagentslist = []
    rateconstantindex = 1
    for i in range(0, min((maxreagents), (len(possiblereagents)))):
        reagentslist.append(possiblereagents[i])
        # print 'reagentslist: '+str(reagentslist)
    for i in range(0, max(1, (maxrxn - 1))):
        rateconstantstr = 'k_' + str(rateconstantindex)
        rateconstant = var(rateconstantstr)
        rxn = randomrxn(reagentslist, 2, 2)
        # set everything to aqueous solutions
        for j in range(0, len(rxn)):
            rxn[j][3] = 'aq'
            rxn[j][4] = 'solution'
        appendrxntomech(mech, rxn, rateconstant)
        rateconstantindex = rateconstantindex + 1
        if (round(random()) == 1.0):  # reversible rxn, so add the reaction in reverse
            rateconstantstr = 'k_' + str(rateconstantindex)
            rateconstant = var(rateconstantstr)
            reverserxn = []
            for j in range(0, len(rxn)):
                reverseindex = len(rxn) - j - 1
                if (rxn[reverseindex][2] == 'product'):
                    reverserxn.append([rxn[reverseindex][0], rxn[reverseindex][1], 'reactant', rxn[reverseindex][3],
                                       rxn[reverseindex][4]])
                else:
                    reverserxn.append([rxn[reverseindex][0], rxn[reverseindex][1], 'product', rxn[reverseindex][3],
                                       rxn[reverseindex][4]])
            appendrxntomech(mech, reverserxn, rateconstant)
            rateconstantindex = rateconstantindex + 1
        # now readjust the reagentlist so that reactions will not be repeated.  The dumb way is to just remove
        # the first reagent in the list.
        del reagentslist[0]
    # print 'reagentslist: '+str(reagentslist)
    return (mech)