diff --git a/examples/rmg/liquid_phase/input.py b/examples/rmg/liquid_phase/input.py new file mode 100644 index 00000000000..236771e8a49 --- /dev/null +++ b/examples/rmg/liquid_phase/input.py @@ -0,0 +1,59 @@ +# Data sources +database( + thermoLibraries = ['primaryThermoLibrary'], + reactionLibraries = [], + seedMechanisms = [], + kineticsDepositories = ['training'], + kineticsFamilies = ['!Intra_Disproportionation'], + kineticsEstimator = 'rate rules', +) + +# List of species +species( + label='octane', + reactive=True, + structure=SMILES("C(CCCCC)CC"), +) + +species( + label='oxygen', + reactive=True, + structure=SMILES("[O][O]"), +) + +# Reaction systems +liquidReactor( + temperature=(500,'K'), + initialConcentrations={ + "octane": (6.154e-3,'mol/cm^3'), + "oxygen": (4.953e-6,'mol/cm^3') + }, + terminationConversion={ + 'octane': 0.9, + }, + terminationTime=(1e6,'s'), +) + +solvation( + solvent='octane' +) + +simulator( + atol=1e-16, + rtol=1e-8, +) + +model( + toleranceKeepInEdge=1E-9, + toleranceMoveToCore=0.001, + toleranceInterruptSimulation=0.1, + maximumEdgeSpecies=100000 +) + +options( + units='si', + saveRestartPeriod=None, + drawMolecules=False, + generatePlots=False, + saveConcentrationProfiles=True, +) diff --git a/examples/rmg/minimal/input.py b/examples/rmg/minimal/input.py index adce7e561a2..10839f4cbad 100644 --- a/examples/rmg/minimal/input.py +++ b/examples/rmg/minimal/input.py @@ -28,6 +28,10 @@ terminationTime=(1e6,'s'), ) +solvation( + solvent='water' +) + simulator( atol=1e-16, rtol=1e-8, diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index acd0b1bafd4..6a0311d6864 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -39,7 +39,7 @@ from rmgpy.data.base import Database, Entry, LogicNode, LogicOr, ForbiddenStructures,\ ForbiddenStructureException, getAllCombinations -from rmgpy.reaction import Reaction, ReactionError +from rmgpy.reaction import Reaction from rmgpy.kinetics import Arrhenius, ArrheniusEP, ThirdBody, Lindemann, Troe, \ PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, \ Chebyshev, KineticsData, PDepKineticsModel diff --git a/rmgpy/data/rmg.py b/rmgpy/data/rmg.py index 80be92a6a4c..c729d6bf17d 100644 --- a/rmgpy/data/rmg.py +++ b/rmgpy/data/rmg.py @@ -39,6 +39,7 @@ from thermo import ThermoDatabase from kinetics import KineticsDatabase from statmech import StatmechDatabase +from solvation import SolvationDatabase # Module-level variable to store the (only) instance of RMGDatabase in use. database = None @@ -55,6 +56,7 @@ def __init__(self): self.forbiddenStructures = None self.kinetics = None self.statmech = None + self.solvation = None # Store the newly created database in the module. global database @@ -69,7 +71,8 @@ def load(self, kineticsFamilies=None, kineticsDepositories=None, statmechLibraries=None, - depository=True + depository=True, + solvation=True, ): """ Load the RMG database from the given `path` on disk, where `path` @@ -87,6 +90,9 @@ def load(self, kineticsDepositories ) self.loadStatmech(os.path.join(path, 'statmech'), statmechLibraries, depository) + + if solvation: + self.loadSolvation(os.path.join(path, 'solvation')) def loadThermo(self, path, thermoLibraries=None, depository=True): """ @@ -133,6 +139,14 @@ def loadKinetics(self, depositories=kineticsDepositories ) + def loadSolvation(self, path): + """ + Load the RMG solvation database from the given `path` on disk, where + `path` points to the top-level folder of the RMG solvation database. + """ + self.solvation = SolvationDatabase() + self.solvation.load(path) + def loadStatmech(self, path, statmechLibraries=None, depository=True): """ Load the RMG statmech database from the given `path` on disk, where diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py new file mode 100644 index 00000000000..2e40bd06ae4 --- /dev/null +++ b/rmgpy/data/solvation.py @@ -0,0 +1,765 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +################################################################################ +# +# RMG - Reaction Mechanism Generator +# +# Copyright (c) 2002-2010 Prof. William H. Green (whgreen@mit.edu) and the +# RMG Team (rmg_dev@mit.edu) +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the 'Software'), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# +################################################################################ + +""" + +""" +import rmgpy.quantity as quantity + +import os +import os.path +import math +import logging +import numpy +from copy import copy, deepcopy + +from base import Database, Entry, makeLogicNode + +import rmgpy.constants as constants +from rmgpy.molecule import Molecule, Atom, Bond, Group + +################################################################################ + +def saveEntry(f, entry): + """ + Write a Pythonic string representation of the given `entry` in the thermo + database to the file object `f`. + """ + raise NotImplementedError() + +def generateOldLibraryEntry(data): + """ + Return a list of values used to save entries to the old-style RMG + thermo database based on the thermodynamics object `data`. + """ + raise NotImplementedError() + +def processOldLibraryEntry(data): + """ + Process a list of parameters `data` as read from an old-style RMG + thermo database, returning the corresponding thermodynamics object. + """ + raise NotImplementedError() + + +class SolventData(): + """ + Stores Abraham/Mintz parameters for characterizing a solvent. + """ + def __init__(self, s_h=None, b_h=None, e_h=None, l_h=None, a_h=None, + c_h=None, s_g=None, b_g=None, e_g=None, l_g=None, a_g=None, c_g=None, A=None, B=None, + C=None, D=None, E=None, alpha=None, beta=None, eps=None): + self.s_h = s_h + self.b_h = b_h + self.e_h = e_h + self.l_h = l_h + self.a_h = a_h + self.c_h = c_h + self.s_g = s_g + self.b_g = b_g + self.e_g = e_g + self.l_g = l_g + self.a_g = a_g + self.c_g = c_g + # These are parameters for calculating viscosity + self.A = A + self.B = B + self.C = C + self.D = D + self.E = E + # These are SOLUTE parameters used for intrinsic rate correction in H-abstraction rxns + self.alpha = alpha + self.beta = beta + # This is the dielectric constant + self.eps = eps + + def getIntrinsicCorrection(self): + """ + If solvation is on, this will give the log10 of the ratio of the intrinsic rate + constants log10(k_sol/k_gas) for H-abstraction rxns + """ + return -8.3*self.alpha*self.beta + + def getSolventViscosity(self, T): + """ + Returns the viscosity in Pa s, according to correlation in Perry's Handbook + and coefficients in DIPPR + """ + return math.exp(self.A + (self.B / T) + (self.C*math.log(T)) + (self.D * (T**self.E))) + +class SolvationCorrection(): + """ + Stores corrections for enthalpy, entropy, and Gibbs free energy when a species is solvated. + """ + def __init__(self, enthalpy=None, gibbs=None, entropy=None): + self.enthalpy = enthalpy + self.entropy = entropy + self.gibbs = gibbs + +class SoluteData(): + """ + Stores Abraham parameters to characterize a solute + """ + def __init__(self, S=None, B=None, E=None, L=None, A=None, V=None, comment=""): + self.S = S + self.B = B + self.E = E + self.L = L + self.A = A + self.V = V + self.comment = comment + def __repr__(self): + return "SoluteData(S={0},B={1},E={2},L={3},A={4},comment={5!r})".format(self.S, self.B, self.E, self.L, self.A, self.comment) + + def getStokesDiffusivity(self, T, solventViscosity): + """ + Get diffusivity of solute using the Stokes-Einstein sphere relation. Radius is + found from the McGowan volume. + """ + k_b = 1.3806488e-23 # m2*kg/s2/K + radius = ((75*self.V/3.14159)**(1/3))/100 # in meters + D = k_b*T/6/3.14159/solventViscosity/radius # m2/s + return D + + def setMcGowanVolume(self, species): + """ + Find and store the McGowan's Volume + Returned volumes are in cm^3/mol/100 (see note below) + See Table 2 in Abraham & McGowan, Chromatographia Vol. 23, No. 4, p. 243. April 1987 + doi: 10.1007/BF02311772 + + "V is scaled to have similar values to the other + descriptors by division by 100 and has units of (cm3mol−1/100)." + the contibutions in this function are in cm3/mol, and the division by 100 is done at the very end. + """ + molecule = species.molecule[0] # any will do, use the first. + Vtot = 0 + + for atom in molecule.atoms: + thisV = 0.0 + if atom.isCarbon(): + thisV = 16.35 + elif (atom.element.number == 7): # nitrogen, do this way if we don't have an isElement method + thisV = 14.39 + elif atom.isOxygen(): + thisV = thisV + 12.43 + #else if (element.equals("F")) + #thisV = thisV + 10.48; + elif atom.isHydrogen(): + thisV = thisV + 8.71 + # else if (element.equals("Si")) +# thisV = thisV + 26.83; +# else if (element.equals("P")) +# thisV = thisV + 24.87; +# else if (element.equals("S")) +# thisV = thisV + 22.91; +# else if (element.equals("Cl")) +# thisV = thisV + 20.95; +# else if (element.equals("B")) +# thisV = thisV + 18.32; +# else if (element.equals("Ge")) +# thisV = thisV + 31.02; +# else if (element.equals("As")) +# thisV = thisV + 29.42; +# else if (element.equals("Se")) +# thisV = thisV + 27.81; +# else if (element.equals("Br")) +# thisV = thisV + 26.21; +# else if (element.equals("Sn")) +# thisV = thisV + 39.35; +# else if (element.equals("Te")) +# thisV = thisV + 36.14; +# else if (element.equals("I")) +# thisV = thisV + 34.53; + else: + raise Exception() + Vtot = Vtot + thisV + + for bond in molecule.getBonds(atom): + # divide contribution in half since all bonds would be counted twice this way + Vtot = Vtot - 6.56/2 + + self.V= Vtot / 100; # division by 100 to get units correct. + +################################################################################ + + +################################################################################ + +class SolventLibrary(Database): + """ + A class for working with a RMG solvent library. + """ + def __init__(self, label='', name='', shortDesc='', longDesc=''): + Database.__init__(self, label=label, name=name, shortDesc=shortDesc, longDesc=longDesc) + + def loadEntry(self, + index, + label, + solvent, + reference=None, + referenceType='', + shortDesc='', + longDesc='', + history=None + ): + self.entries[label] = Entry( + index = index, + label = label, + data = solvent, + reference = reference, + referenceType = referenceType, + shortDesc = shortDesc, + longDesc = longDesc.strip(), + history = history or [], + ) + + def load(self, path): + """ + Load the solvent library from the given path + """ + Database.load(self, path, local_context={'SolventData': SolventData}, global_context={}) + + def saveEntry(self, f, entry): + """ + Write the given `entry` in the solute database to the file object `f`. + """ + return saveEntry(f, entry) + + def getSolventData(self, label): + """ + Get a solvent's data from its name + """ + return self.entries[label].data + + +class SoluteLibrary(Database): + """ + A class for working with a RMG solute library. Not currently used. + """ + def __init__(self, label='', name='', shortDesc='', longDesc=''): + Database.__init__(self, label=label, name=name, shortDesc=shortDesc, longDesc=longDesc) + + def loadEntry(self, + index, + label, + molecule, + solute, + reference=None, + referenceType='', + shortDesc='', + longDesc='', + history=None + ): + self.entries[label] = Entry( + index = index, + label = label, + item = Molecule(SMILES=molecule), + data = solute, + reference = reference, + referenceType = referenceType, + shortDesc = shortDesc, + longDesc = longDesc.strip(), + history = history or [], + ) + + def load(self, path): + pass + + def saveEntry(self, f, entry): + """ + Write the given `entry` in the solute database to the file object `f`. + """ + return saveEntry(f, entry) + + def generateOldLibraryEntry(self, data): + """ + Return a list of values used to save entries to the old-style RMG + thermo database based on the thermodynamics object `data`. + """ + return generateOldLibraryEntry(data) + + def processOldLibraryEntry(self, data): + """ + Process a list of parameters `data` as read from an old-style RMG + thermo database, returning the corresponding thermodynamics object. + """ + return processOldLibraryEntry(data) + +################################################################################ + +class SoluteGroups(Database): + """ + A class for working with an RMG solute group additivity database. + """ + + def __init__(self, label='', name='', shortDesc='', longDesc=''): + Database.__init__(self, label=label, name=name, shortDesc=shortDesc, longDesc=longDesc) + + def loadEntry(self, + index, + label, + group, + solute, + reference=None, + referenceType='', + shortDesc='', + longDesc='', + history=None + ): + if group[0:3].upper() == 'OR{' or group[0:4].upper() == 'AND{' or group[0:7].upper() == 'NOT OR{' or group[0:8].upper() == 'NOT AND{': + item = makeLogicNode(group) + else: + item = Group().fromAdjacencyList(group) + self.entries[label] = Entry( + index = index, + label = label, + item = item, + data = solute, + reference = reference, + referenceType = referenceType, + shortDesc = shortDesc, + longDesc = longDesc.strip(), + history = history or [], + ) + + def saveEntry(self, f, entry): + """ + Write the given `entry` in the thermo database to the file object `f`. + """ + return saveEntry(f, entry) + + def generateOldLibraryEntry(self, data): + """ + Return a list of values used to save entries to the old-style RMG + thermo database based on the thermodynamics object `data`. + """ + + return generateOldLibraryEntry(data) + + def processOldLibraryEntry(self, data): + """ + Process a list of parameters `data` as read from an old-style RMG + thermo database, returning the corresponding thermodynamics object. + """ + return processOldLibraryEntry(data) + +################################################################################ + +class SolvationDatabase(object): + """ + A class for working with the RMG solute database. + """ + + def __init__(self): + self.solventLibrary = SolventLibrary() + self.soluteLibrary = SoluteLibrary() + self.groups = {} + self.local_context = { + 'SoluteData': SoluteData, + 'SolventData': SolventData + } + self.global_context = {} + + def __reduce__(self): + """ + A helper function used when pickling a SoluteDatabase object. + """ + d = { + 'libraries': self.libraries, + 'groups': self.groups, + 'libraryOrder': self.libraryOrder, + } + return (SoluteDatabase, (), d) + + def __setstate__(self, d): + """ + A helper function used when unpickling a SoluteDatabase object. + """ + self.libraries = d['libraries'] + self.groups = d['groups'] + self.libraryOrder = d['libraryOrder'] + + def load(self, path, libraries=None, depository=True): + """ + Load the solvation database from the given `path` on disk, where `path` + points to the top-level folder of the solvation database. + + Load the solvent and solute (not used) libraries, then the solute groups. + """ + + self.solventLibrary.load(os.path.join(path,'libraries','solvent.py')) + self.soluteLibrary.load(os.path.join(path,'libraries','solute.py')) + + self.loadGroups(os.path.join(path, 'groups')) + + def getSolventData(self, solvent_name): + return self.solventLibrary.getSolventData(solvent_name) + + + def loadGroups(self, path): + """ + Load the solute database from the given `path` on disk, where `path` + points to the top-level folder of the solute database. + + Two sets of groups for additivity, atom-centered ('abraham') and non atom-centered + ('nonacentered'). + """ + logging.info('Loading Platts additivity group database from {0}...'.format(path)) + self.groups = {} + self.groups['abraham'] = SoluteGroups(label='abraham').load(os.path.join(path, 'abraham.py' ), self.local_context, self.global_context) + self.groups['nonacentered'] = SoluteGroups(label='nonacentered').load(os.path.join(path, 'nonacentered.py' ), self.local_context, self.global_context) + + def save(self, path): + """ + Save the solvation database to the given `path` on disk, where `path` + points to the top-level folder of the solvation database. + """ + path = os.path.abspath(path) + if not os.path.exists(path): os.mkdir(path) + self.saveLibraries(os.path.join(path, 'libraries')) + self.saveGroups(os.path.join(path, 'groups')) + + def saveLibraries(self, path): + """ + Save the solute libraries to the given `path` on disk, where `path` + points to the top-level folder of the solute libraries. + """ + if not os.path.exists(path): os.mkdir(path) + for library in self.libraries.values(): + library.save(os.path.join(path, '{0}.py'.format(library.label))) + + def saveGroups(self, path): + """ + Save the solute groups to the given `path` on disk, where `path` + points to the top-level folder of the solute groups. + """ + if not os.path.exists(path): os.mkdir(path) + self.groups['abraham'].save(os.path.join(path, 'abraham.py')) + self.groups['nonacentered'].save(os.path.join(path, 'nonacentered.py')) + + def loadOld(self, path): + """ + Load the old RMG solute database from the given `path` on disk, where + `path` points to the top-level folder of the old RMG database. + """ + + for (root, dirs, files) in os.walk(os.path.join(path, 'thermo_libraries')): + if os.path.exists(os.path.join(root, 'Dictionary.txt')) and os.path.exists(os.path.join(root, 'Library.txt')): + library = SoluteLibrary(label=os.path.basename(root), name=os.path.basename(root)) + library.loadOld( + dictstr = os.path.join(root, 'Dictionary.txt'), + treestr = '', + libstr = os.path.join(root, 'Library.txt'), + numParameters = 5, + numLabels = 1, + pattern = False, + ) + library.label = os.path.basename(root) + self.libraries[library.label] = library + + self.groups = {} + self.groups['abraham'] = SoluteGroups(label='abraham', name='Platts Group Additivity Values for Abraham Solute Descriptors').loadOld( + dictstr = os.path.join(path, 'thermo_groups', 'Abraham_Dictionary.txt'), + treestr = os.path.join(path, 'thermo_groups', 'Abraham_Tree.txt'), + libstr = os.path.join(path, 'thermo_groups', 'Abraham_Library.txt'), + numParameters = 5, + numLabels = 1, + pattern = True, + ) + + def saveOld(self, path): + """ + Save the old RMG Abraham database to the given `path` on disk, where + `path` points to the top-level folder of the old RMG database. + """ + # Depository not used in old database, so it is not saved + + librariesPath = os.path.join(path, 'thermo_libraries') + if not os.path.exists(librariesPath): os.mkdir(librariesPath) + for library in self.libraries.values(): + libraryPath = os.path.join(librariesPath, library.label) + if not os.path.exists(libraryPath): os.mkdir(libraryPath) + library.saveOld( + dictstr = os.path.join(libraryPath, 'Dictionary.txt'), + treestr = '', + libstr = os.path.join(libraryPath, 'Library.txt'), + ) + + groupsPath = os.path.join(path, 'thermo_groups') + if not os.path.exists(groupsPath): os.mkdir(groupsPath) + self.groups['abraham'].saveOld( + dictstr = os.path.join(groupsPath, 'Abraham_Dictionary.txt'), + treestr = os.path.join(groupsPath, 'Abraham_Tree.txt'), + libstr = os.path.join(groupsPath, 'Abraham_Library.txt'), + ) + + def getSoluteData(self, species): + """ + Return the solute descriptors for a given :class:`Species` + object `species`. This function first searches the loaded libraries + in order, returning the first match found, before falling back to + estimation via Platts group additivity. + """ + soluteData = None + + # Check the library first + soluteData = self.getSoluteDataFromLibrary(species, self.soluteLibrary) + if soluteData is not None: + soluteData[0].comment = 'solute' + else: + # Solute not found in any loaded libraries, so estimate + soluteData = self.getSoluteDataFromGroups(species) + # No Platts group additivity for V, so set using atom sizes + soluteData.setMcGowanVolume(species) + # Return the resulting solute parameters S, B, E, L, A + return soluteData + + def getAllSoluteData(self, species): + """ + Return all possible sets of Abraham solute descriptors for a given + :class:`Species` object `species`. The hits from the library come + first, then the group additivity estimate. This method is useful + for a generic search job. + """ + raise NotImplementedError() + + def getSoluteDataFromLibrary(self, species, library): + """ + Return the set of Abraham solute descriptors corresponding to a given + :class:`Species` object `species` from the specified solute + `library`. If `library` is a string, the list of libraries is searched + for a library with that name. If no match is found in that library, + ``None`` is returned. If no corresponding library is found, a + :class:`DatabaseError` is raised. + """ + for label, entry in library.entries.iteritems(): + for molecule in species.molecule: + if molecule.isIsomorphic(entry.item) and entry.data is not None: + return (deepcopy(entry.data), library, entry) + return None + + def getSoluteDataFromGroups(self, species): + """ + Return the set of Abraham solute parameters corresponding to a given + :class:`Species` object `species` by estimation using the Platts group + additivity method. If no group additivity values are loaded, a + :class:`DatabaseError` is raised. + + It averages (linearly) over the desciptors for each Molecule (resonance isomer) + in the Species. + """ + soluteData = SoluteData(0.0,0.0,0.0,0.0,0.0) + count = 0 + comments = [] + for molecule in species.molecule: + molecule.clearLabeledAtoms() + molecule.updateAtomTypes() + sdata = self.estimateSoluteViaGroupAdditivity(molecule) + soluteData.S += sdata.S + soluteData.B += sdata.B + soluteData.E += sdata.E + soluteData.L += sdata.L + soluteData.A += sdata.A + count += 1 + comments.append(sdata.comment) + + soluteData.S /= count + soluteData.B /= count + soluteData.E /= count + soluteData.L /= count + soluteData.A /= count + + # Print groups that are used for debugging purposes + soluteData.comment = "Average of {0}".format(" and ".join(comments)) + + return soluteData + + def estimateSoluteViaGroupAdditivity(self, molecule): + """ + Return the set of Abraham solute parameters corresponding to a given + :class:`Molecule` object `molecule` by estimation using the Platts' group + additivity method. If no group additivity values are loaded, a + :class:`DatabaseError` is raised. + """ + # For thermo estimation we need the atoms to already be sorted because we + # iterate over them; if the order changes during the iteration then we + # will probably not visit the right atoms, and so will get the thermo wrong + molecule.sortVertices() + + # Create the SoluteData object with the intercepts from the Platts groups + soluteData = SoluteData( + S = 0.277, + B = 0.071, + E = 0.248, + L = 0.13, + A = 0.003 + ) + + if sum([atom.radicalElectrons for atom in molecule.atoms]) > 0: # radical species + + # Make a copy of the structure so we don't change the original + saturatedStruct = molecule.copy(deep=True) + + # Saturate structure by replacing all radicals with bonds to + # hydrogen atoms + added = {} + for atom in saturatedStruct.atoms: + for i in range(atom.radicalElectrons): + H = Atom('H') + bond = Bond(atom, H, 'S') + saturatedStruct.addAtom(H) + saturatedStruct.addBond(bond) + if atom not in added: + added[atom] = [] + added[atom].append([H, bond]) + atom.decrementRadical() + + # Update the atom types of the saturated structure (not sure why + # this is necessary, because saturating with H shouldn't be + # changing atom types, but it doesn't hurt anything and is not + # very expensive, so will do it anyway) + saturatedStruct.updateConnectivityValues() + saturatedStruct.sortVertices() + saturatedStruct.updateAtomTypes() + + # Get solute descriptor estimates for saturated form of structure + soluteData = self.estimateSoluteViaGroupAdditivity(saturatedStruct) + assert soluteData is not None, "Solute data of saturated {0} of molecule {1} is None!".format(saturatedStruct, molecule) + + # For each radical site, get radical correction + # Only one radical site should be considered at a time; all others + # should be saturated with hydrogen atoms + for atom in added: + + # Remove the added hydrogen atoms and bond and restore the radical + for H, bond in added[atom]: + saturatedStruct.removeBond(bond) + saturatedStruct.removeAtom(H) + atom.incrementRadical() + + saturatedStruct.updateConnectivityValues() + + else: # non-radical species + # Generate estimate of solute data + for atom in molecule.atoms: + # Iterate over heavy (non-hydrogen) atoms + if atom.isNonHydrogen(): + # Get initial solute data from main group database. Every atom must + # be found in the main abraham database + try: + self.__addGroupSoluteData(soluteData, self.groups['abraham'], molecule, {'*':atom}) + except KeyError: + logging.error("Couldn't find in main abraham database:") + logging.error(molecule) + logging.error(molecule.toAdjacencyList()) + raise + # Get solute data for non-atom centered groups (being found in this group + # database is optional) + try: + self.__addGroupSoluteData(soluteData, self.groups['nonacentered'], molecule, {'*':atom}) + except KeyError: pass + + return soluteData + + def __addGroupSoluteData(self, soluteData, database, molecule, atom): + """ + Determine the Platts group additivity solute data for the atom `atom` + in the structure `structure`, and add it to the existing solute data + `soluteData`. + """ + + node0 = database.descendTree(molecule, atom, None) + + if node0 is None: + raise KeyError('Node not found in database.') + + # It's possible (and allowed) that items in the tree may not be in the + # library, in which case we need to fall up the tree until we find an + # ancestor that has an entry in the library + node = node0 + + while node is not None and node.data is None: + node = node.parent + if node is None: + raise KeyError('Node has no parent with data in database.') + data = node.data + comment = node.label + while isinstance(data, basestring) and data is not None: + for entry in database.entries.values(): + if entry.label == data: + data = entry.data + comment = entry.label + break + comment = '{0}({1})'.format(database.label, comment) + + # This code prints the hierarchy of the found node; useful for debugging + #result = '' + #while node is not None: + # result = ' -> ' + node + result + # node = database.tree.parent[node] + #print result[4:] + + # Add solute data for each atom to the overall solute data for the molecule. + soluteData.S += data.S + soluteData.B += data.B + soluteData.E += data.E + soluteData.L += data.L + soluteData.A += data.A + soluteData.comment += comment + "+" + + return soluteData + + + def calcH(self, soluteData, solventData): + # Use Mintz parameters for solvents + delH = 1000*((soluteData.S*solventData.s_h)+(soluteData.B*solventData.b_h)+(soluteData.E*solventData.e_h)+(soluteData.L*solventData.l_h)+(soluteData.A*solventData.a_h)+solventData.c_h) + return delH + + def calcG(self, soluteData, solventData): + # Use Abraham parameters for solvents + logK = (soluteData.S*solventData.s_g)+(soluteData.B*solventData.b_g)+(soluteData.E*solventData.e_g)+(soluteData.L*solventData.l_g)+(soluteData.A*solventData.a_g)+solventData.c_g + delG = -8.314*298*2.303*logK + return delG + + def calcS(self, delG, delH): + delS = (delH-delG)/298 + return delS + + def getSolvationCorrection(self, soluteData, solventData): + """ + Given a soluteData and solventData object, calculates the enthalpy, entropy, + and Gibbs free energy of solvation at 298 K + """ + correction = SolvationCorrection(0.0, 0.0, 0.0) + correction.enthalpy = self.calcH(soluteData, solventData) + correction.gibbs = self.calcG(soluteData, solventData) + correction.entropy = self.calcS(correction.gibbs, correction.enthalpy) + return correction \ No newline at end of file diff --git a/rmgpy/data/solvationTest.py b/rmgpy/data/solvationTest.py new file mode 100644 index 00000000000..af8e3157f18 --- /dev/null +++ b/rmgpy/data/solvationTest.py @@ -0,0 +1,109 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import os +from unittest import TestCase + +from rmgpy import settings +from rmgpy.species import Species +from rmgpy.data.solvation import * +from rmgpy.molecule.molecule import Molecule + +################################################### + +class TestSoluteDatabase(TestCase): + + def runTest(self): + pass + + def testSoluteGeneration(self): + + self.database = SolvationDatabase() + self.database.load(os.path.join(settings['database.directory'], 'solvation')) + + self.testCases = [ + # from RMG-Java test runs by RWest (mostly in agreement with Jalan et. al. supplementary data) + # ['octane', 'CCCCCCCC', 0.127, 0.085, 0.04, 3.766, 0.0030 , 1.2358], +# ['water', 'O', 0.524, 0.378, 0.309, 0.802, 0.348, 0.1673], +# ['chrysene', 'C1=CC=CC3=C1C=CC4=C2C=CC=CC2=CC=C34', 1.603, 0.317, 2.864, 10.222, 0.003, None], +# ['anthracene', 'C2=CC=CC3=CC1=CC=CC=C1C=C23', 1.261, 0.257, 2.128, 7.796, 0.003, None], + ['1,2-ethanediol', 'C(CO)O', 0.823, 0.685, 0.327, 2.572, 0.693, None], + # ['1-decanol', 'C(CCCCCCCCC)O', 0.449, 0.385, 0.205, 5.614, 0.348, None], +# ['acetylaldehyde', 'CC=O', 0.622, 0.423, 0.171, 1.415, 0.0030, 0.4061], +# ['acetic acid', 'C(C)(=O)O', 0.508, 0.411, 0.152, 1.873, 0.591, None], +# ['ethyloctadecanoate', 'C(C)OC(CCCCCCCCCCCCCCCCC)=O', 0.558, 0.424, 0.08, 10.344, 0.003, None], +# ['naphthalene', 'C1=CC=CC2=CC=CC=C12', 0.919, 0.197, 1.392, 5.37, 0.003, None], +# ['phenol', 'C1(=CC=CC=C1)O', 0.875, 0.433, 0.829, 3.771, 0.546, None], +# ['m-cresol', 'C1=C(C=CC=C1O)C', 0.851, 0.429, 0.837, 4.247, 0.546, None], +# ['m-hydroxybenzaldehyde', 'OC1=CC=CC(=C1)C=O', 1.346, 0.767, 0.968, 4.89, 0.546, None], +# ['ethylbenzene', 'C(C)C1=CC=CC=C1', 0.553, 0.133, 0.664, 3.919, 0.003, None], +# ['toluene', 'C1(=CC=CC=C1)C', 0.553, 0.133, 0.664, 3.42, 0.003, None], +# ['1-butene', 'C=CCC', 0.167, 0.108, 0.167, 1.663, 0.003, None], +# ['g-decalactone', 'C(C1OC(=O)CC1)CCCCC', 0.669, 0.428, 0.273, 5.482, 0.003, None], + + #['dimethyl ether'] + ] + + for name, smiles, S, B, E, L, A, V in self.testCases: + species = Species(molecule=[Molecule(SMILES=smiles)]) + soluteData = self.database.getSoluteDataFromGroups(Species(molecule=[species.molecule[0]])) + print name, soluteData + print self.assertAlmostEqual(soluteData.S, S) + print self.assertAlmostEqual(soluteData.B, B) + print self.assertAlmostEqual(soluteData.E, E) + print self.assertAlmostEqual(soluteData.L, L) + print self.assertAlmostEqual(soluteData.A, A) + + def testCorrectionGeneration(self): + self.database = SolvationDatabase() + self.database.load(os.path.join(settings['database.directory'], 'solvation')) + self.testCases = [ + # solventName, soluteName, soluteSMILES, Hsolv, Gsolv, Ssolv + ['water', 'methane', 'C', -12000, 2000*4.184, None], + ['water', 'octane', 'CCCCCCCC', -36000, 2890*4.184, None], + ['water', '1,2-ethanediol', 'C(CO)O', -77300, -9300*4.184, None], + ['water', 'acetic acid', 'C(C)(=O)O', -56500, -6700*4.184, None], + ['water', 'naphthalene', 'C1=CC=CC2=CC=CC=C12', -42800, -2390*4.184, None], + ['water', 'm-hydroxybenzaldehyde', 'OC1=CC=CC(=C1)C=O', -70700, -9510*4.184, None], + ['water', 'ethylbenzene', 'C(C)C1=CC=CC=C1', -39400, -800*4.184, None], + ['water', 'toluene', 'C1(=CC=CC=C1)C', -32400, -890*4.184, None], + ['water', 'ethane', 'CC', -17900, 1830*4.184, None], + ['water', 'propane', 'CCC', -20400, 1960*4.184, None], + ['water', 'ethene', 'C=C', -13700, 1270*4.184, None], + ['water', 'propene', 'CC=C', -21600, 1270*4.184, None], + ['water', 'dimethyl ether', 'COC', -34000, -1920*4.184, None], + ['water', 'diethyl ether', 'C(C)OCC', -45300, -1760*4.184, None], + ['water', 'tetrahydrofuran', 'C1CCOC1', -47300, -3470*4.184, None], + ['water', '1,4-dioxane', 'C1COCCO1', -48400, -5050*4.184, None], + ['water', 'methanol', 'CO', -52000, -5110*4.184, None], + ['water', 'ethanol', 'C(C)O', -50600, -5010*4.184, None], + ['water', '1,2 propanediol', 'C(CO)CO', -81100, None, None], + ['1-octanol', 'methane', 'C', -3900, 510*4.184, None], + ['1-octanol', 'octane', 'CCCCCCCC', -40080, -4180*4.184, None], + ['1-octanol', 'ethylbenzene', 'C(C)C1=CC=CC=C1', -39910, -5080*4.184, None], + ['1-octanol', 'toluene', 'C1(=CC=CC=C1)C', -35990, -4550*4.184, None], + ['1-octanol', 'tetrahydrofuran', 'C1CCOC1', -28320, -3930*4.184, None], + ['benzene', 'octane', 'CCCCCCCC', -35480, -5350*4.184, None], + ['benzene', 'toluene', 'C1(=CC=CC=C1)C', -37660, -5320*4.184, None], + ['benzene', '1,4-dioxane', 'C1COCCO1', -39030, -5210*4.184, None], + ['benzene', 'methanol', 'CO', None, -2580*4.184, None], + ['benzene', 'ethanol', 'C(CO)CO', -29620, -3420*4.184, None], + + ] + + for solventName, soluteName, smiles, H, G, S in self.testCases: + species = Species(molecule=[Molecule(SMILES=smiles)]) + soluteData = self.database.getSoluteDataFromGroups(Species(molecule=[species.molecule[0]])) + solventData = self.database.getSolventData(solventName) + solvationCorrection = self.database.getSolvationCorrection(soluteData, solventData) + #print("Enthalpy of solvation for {0} in {1} is {2} J/mol".format(soluteName, solventName, solvationCorrection.enthalpy)) + #print("Enthalpy: {0} {1} {2}".format(soluteName, H, solvationCorrection.enthalpy)) + print("Gibbs: {0} {1} {2}".format(soluteName, G, solvationCorrection.gibbs)) + #self.assertAlmostEqual(solvationCorrection.enthalpy/10000., H/10000.,0) #0 decimal place, in 10kJ. + +##################################################### + +if __name__ == '__main__': + myTest = TestSoluteDatabase() + #myTest.testSoluteGeneration() + myTest.testCorrectionGeneration() \ No newline at end of file diff --git a/rmgpy/kinetics/diffusionLimited.py b/rmgpy/kinetics/diffusionLimited.py new file mode 100644 index 00000000000..04eb9d12f23 --- /dev/null +++ b/rmgpy/kinetics/diffusionLimited.py @@ -0,0 +1,90 @@ +import rmgpy.quantity as quantity +import logging +from rmgpy.species import Species +from rmgpy.data.solvation import SolventData, SoluteData, SoluteGroups, SolvationDatabase +from rmgpy.reaction import Reaction + + +class DiffusionLimited(): + + def __init__(self): + # default is false, enabled if there is a solvent + self.enabled = False + + def enable(self, solventData, solvationDatabase, comment=''): + # diffusionLimiter is enabled if a solvent has been added to the RMG object. + logging.info("Enabling diffusion-limited kinetics...") + diffusionLimiter.enabled = True + diffusionLimiter.database = solvationDatabase + diffusionLimiter.solventData = solventData + + def getSolventViscosity(self, T): + return self.solventData.getSolventViscosity(T) + + def getEffectiveRate(self, reaction, T): + """ + Return the ratio of k_eff to k_intrinsic, which is between 0 and 1. + + It is 1.0 if diffusion has no effect. + + For 1<=>2 reactions, the reverse rate is limited. + For 2<=>2 reactions, the faster direction is limited. + For 2<=>1 or 2<=>3 reactions, the forward rate is limited. + """ + intrinsicKinetics = reaction.kinetics + reactants = len(reaction.reactants) + products = len(reaction.products) + k_forward = intrinsicKinetics.getRateCoefficient(T,P=100e5) + Keq = reaction.getEquilibriumConstant(T) # Kc + k_reverse = k_forward / Keq + k_eff = k_forward + + if reactants == 1: + if products == 1: + k_eff = k_forward + else: # two products; reverse rate is limited + k_diff = self.getDiffusionLimit(T, reaction, forward=False) + k_eff_reverse = k_reverse*k_diff/(k_reverse+k_diff) + k_eff = k_eff_reverse * Keq + else: # 2 reactants + if products == 1 or products == 3: + k_diff = self.getDiffusionLimit(T, reaction, forward=True) + k_eff = k_forward*k_diff/(k_forward+k_diff) + else: # 2 products + if Keq > 1.0: # forward rate is faster and thus limited + k_diff = self.getDiffusionLimit(T, reaction, forward=True) + k_eff = k_forward*k_diff/(k_forward+k_diff) + else: # reverse rate is faster and thus limited + k_diff = self.getDiffusionLimit(T, reaction, forward=False) + k_eff_reverse = k_reverse*k_diff/(k_reverse+k_diff) + k_eff = k_eff_reverse * Keq + return k_eff + + def getDiffusionLimit(self, T, reaction, forward=True): + """ + Return the diffusive limit on the rate coefficient, k_diff. + + This is the upper limit on the rate, in the specified direction. + (ie. forward direction if forward=True [default] or reverse if forward=False) + """ + if forward: + reacting = reaction.reactants + else: + reacting = reaction.products + assert len(reacting)==2, "Can only calculate diffusion limit in a bimolecular direction" + radii = 0.0 + diffusivities = 0.0 + for spec in reacting: + soluteData = self.database.getSoluteData(spec) + # calculate radius with the McGowan volume and assuming sphere + radius = ((75*soluteData.V/3.14159)**(1/3))/100 + diff = soluteData.getStokesDiffusivity(T, self.getSolventViscosity(T)) + radii += radius + diffusivities += diff + N_a = 6.022e23 # Avogadro's Number + k_diff = 4*3.14159*radii*diffusivities*N_a + return k_diff + + +# module level variable. There should only ever be one. It starts off disabled +diffusionLimiter = DiffusionLimited() \ No newline at end of file diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 1eba3e230ad..2a141d95956 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -54,6 +54,8 @@ from rmgpy.kinetics import KineticsData, ArrheniusEP, ThirdBody, Lindemann, Troe, Chebyshev, PDepArrhenius, MultiArrhenius, MultiPDepArrhenius #PyDev: @UnresolvedImport from rmgpy.pdep.reaction import calculateMicrocanonicalRateCoefficient +from rmgpy.kinetics.diffusionLimited import diffusionLimiter + ################################################################################ class ReactionError(Exception): @@ -109,6 +111,9 @@ def __init__(self, self.duplicate = duplicate self.degeneracy = degeneracy self.pairs = pairs + + if diffusionLimiter.enabled: + self.__k_effective_cache = {} def __repr__(self): """ @@ -477,6 +482,8 @@ def getEquilibriumConstant(self, T, type='Kc'): K *= P0 ** (len(self.products) - len(self.reactants)) elif type != 'Ka' and type != '': raise ReactionError('Invalid type "%s" passed to Reaction.getEquilibriumConstant(); should be "Ka", "Kc", or "Kp".') + if K == 0: + raise ReactionError('Got equilibrium constant of 0') return K def getEnthalpiesOfReaction(self, Tlist): @@ -530,8 +537,20 @@ def getRateCoefficient(self, T, P=0): Return the overall rate coefficient for the forward reaction at temperature `T` in K and pressure `P` in Pa, including any reaction path degeneracies. - """ - return self.kinetics.getRateCoefficient(T, P) + + If diffusionLimiter is enabled, the reaction is in the liquid phase and we use + a diffusion limitation to correct the rate. If not, then use the intrinsic rate + coefficient. + """ + if diffusionLimiter.enabled: + try: + k = self.__k_effective_cache[T] + except KeyError: + k = diffusionLimiter.getEffectiveRate(self, T) + self.__k_effective_cache[T] = k + return k + else: + return self.kinetics.getRateCoefficient(T, P) def getRate(self, T, P, conc, totalConc=-1.0): """ @@ -580,7 +599,7 @@ def getRate(self, T, P, conc, totalConc=-1.0): # Return rate return rateConstant * (forward - reverse / equilibriumConstant) - + def fixBarrierHeight(self, forcePositive=False): """ Turns the kinetics into Arrhenius (if they were ArrheniusEP) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 5c50924ab5a..84456727ca6 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -40,6 +40,7 @@ from rmgpy.quantity import Quantity from rmgpy.solver.base import TerminationTime, TerminationConversion from rmgpy.solver.simple import SimpleReactor +from rmgpy.solver.liquid import LiquidReactor from model import CoreEdgeReactionModel @@ -143,9 +144,35 @@ def simpleReactor(temperature, system = SimpleReactor(T, P, initialMoleFractions, termination, sensitivitySpecies, sensitivityThreshold) rmg.reactionSystems.append(system) + +# Reaction systems +def liquidReactor(temperature, initialConcentrations, terminationConversion=None, terminationTime=None): + logging.debug('Found LiquidReactor reaction system') + T = Quantity(temperature) + for spec,conc in initialConcentrations.iteritems(): + concentration = Quantity(conc) + # check the dimensions are ok + # convert to mol/m^3 (or something numerically nice? or must it be SI) + initialConcentrations[spec] = concentration.value_si + termination = [] + if terminationConversion is not None: + for spec, conv in terminationConversion.iteritems(): + termination.append(TerminationConversion(speciesDict[spec], conv)) + if terminationTime is not None: + termination.append(TerminationTime(Quantity(terminationTime))) + if len(termination) == 0: + raise InputError('No termination conditions specified for reaction system #{0}.'.format(len(rmg.reactionSystems)+2)) + system = LiquidReactor(T, initialConcentrations, termination) + rmg.reactionSystems.append(system) + def simulator(atol, rtol): rmg.absoluteTolerance = atol rmg.relativeTolerance = rtol + +def solvation(solvent): + # If solvation module in input file, set the RMG solvent variable + assert isinstance(solvent,str), "solvent should be a string like 'water'" + rmg.solvent = solvent def model(toleranceMoveToCore, toleranceKeepInEdge=0.0, toleranceInterruptSimulation=1.0, maximumEdgeSpecies=None): rmg.fluxToleranceKeepInEdge = toleranceKeepInEdge @@ -278,7 +305,9 @@ def readInputFile(path, rmg0): 'InChI': InChI, 'adjacencyList': adjacencyList, 'simpleReactor': simpleReactor, + 'liquidReactor': liquidReactor, 'simulator': simulator, + 'solvation': solvation, 'model': model, 'quantumMechanics': quantumMechanics, 'pressureDependence': pressureDependence, @@ -295,11 +324,9 @@ def readInputFile(path, rmg0): finally: f.close() + # convert keys from species names into species objects. for reactionSystem in rmg.reactionSystems: - initialMoleFractions = {} - for label, moleFrac in reactionSystem.initialMoleFractions.iteritems(): - initialMoleFractions[speciesDict[label]] = moleFrac - reactionSystem.initialMoleFractions = initialMoleFractions + reactionSystem.convertInitalKeysToSpeciesObjects(speciesDict) logging.info('') @@ -367,7 +394,9 @@ def saveInputFile(path, rmg): f.write(' sensitivityThreshold = {0},\n'.format(system.sensitivity)) f.write(')\n\n') - + + if rmg.solvent: + f.write("solvation(\n solvent = '{0!s}'\n)\n\n".format(solvent)) # Simulator tolerances f.write('simulator(\n') diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 18bf39f3feb..1ef86ab157d 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -52,6 +52,10 @@ from rmgpy.data.rmg import RMGDatabase from rmgpy.data.kinetics import KineticsLibrary, KineticsFamily, LibraryReaction, TemplateReaction +from rmgpy.reaction import Reaction +from rmgpy.kinetics.diffusionLimited import diffusionLimiter +from rmgpy.quantity import Quantity + from model import Species, CoreEdgeReactionModel from pdep import PDepNetwork @@ -76,6 +80,7 @@ class RMG: `kineticsFamilies` The kinetics families to use for reaction generation `kineticsDepositories` The kinetics depositories to use for looking up kinetics in each family `kineticsEstimator` The method to use to estimate kinetics: 'group additivity' or 'rate rules' + `solvent` If solvation estimates are required, the name of the solvent. --------------------------- ------------------------------------------------ `reactionModel` The core-edge reaction model generated by this job `reactionSystems` A list of the reaction systems used in this job @@ -127,6 +132,8 @@ def clear(self): self.kineticsFamilies = None self.kineticsDepositories = None self.kineticsEstimator = 'group additivity' + self.solvent = None + self.diffusionLimiter = None self.reactionModel = None self.reactionSystems = None @@ -271,6 +278,13 @@ def initialize(self, args): # Load databases self.loadDatabase() + + # Do all liquid-phase startup things: + if self.solvent: + Species.solventData = self.database.solvation.getSolventData(self.solvent) + Species.solventName = self.solvent + diffusionLimiter.enable(Species.solventData, self.database.solvation) + logging.info("Setting solvent data for {0}".format(self.solvent)) # Set wall time if args.walltime == '0': diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 199391b9c54..980aa272832 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -48,6 +48,7 @@ from rmgpy.statmech import Conformer from rmgpy.data.thermo import * +from rmgpy.data.solvation import * from rmgpy.data.kinetics import * from rmgpy.data.statmech import * import rmgpy.data.rmg @@ -62,6 +63,10 @@ ################################################################################ class Species(rmgpy.species.Species): + solventName = None + solventData = None + solventViscosity = None + diffusionTemp = None def __init__(self, index=-1, label='', thermo=None, conformer=None, molecule=None, lennardJones=None, molecularWeight=None, @@ -137,9 +142,9 @@ def generateThermoData(self, database, thermoClass=NASA, quantumMechanics=None): if thermo0 is None: thermo0 = database.thermo.getThermoData(self) - return self.processThermoData(thermo0, thermoClass) + return self.processThermoData(database, thermo0, thermoClass) - def processThermoData(self, thermo0, thermoClass=NASA): + def processThermoData(self, database, thermo0, thermoClass=NASA): """ Converts via Wilhoit into required `thermoClass` and sets `E0`. @@ -163,6 +168,20 @@ def processThermoData(self, thermo0, thermoClass=NASA): wilhoit = thermo0.toWilhoit(Cp0=Cp0, CpInf=CpInf) wilhoit.comment = thermo0.comment + + # Add on solvation correction + if Species.solventData: + logging.info("Making solvent correction for {0}".format(Species.solventName)) + soluteData = database.solvation.getSoluteData(self) + solvation_correction = database.solvation.getSolvationCorrection(soluteData, Species.solventData) + # correction is added to the entropy and enthalpy + wilhoit.S0.value_si = (wilhoit.S0.value_si + solvation_correction.entropy) + wilhoit.H0.value_si = (wilhoit.H0.value_si + solvation_correction.enthalpy) + + # Compute E0 by extrapolation to 0 K + if self.conformer is None: + self.conformer = Conformer() + self.conformer.E0 = (wilhoit.getEnthalpy(1.0)*1e-3,"kJ/mol") # Convert to desired thermo class if isinstance(thermo0, thermoClass): diff --git a/rmgpy/rmg/output.py b/rmgpy/rmg/output.py index a05d574fa9d..e5d796612ef 100644 --- a/rmgpy/rmg/output.py +++ b/rmgpy/rmg/output.py @@ -66,7 +66,7 @@ def saveOutputHTML(path, reactionModel): try: import jinja2 except ImportError: - logging.warning("jinja package not found; HTML output will not be saved.") + logging.warning("jinja2 package not found; HTML output will not be saved.") return path = os.path.abspath(path) diff --git a/rmgpy/solver/liquid.pyx b/rmgpy/solver/liquid.pyx new file mode 100644 index 00000000000..715ac2e17e1 --- /dev/null +++ b/rmgpy/solver/liquid.pyx @@ -0,0 +1,696 @@ +################################################################################ +# +# RMG - Reaction Mechanism Generator +# +# Copyright (c) 2002-2010 Prof. William H. Green (whgreen@mit.edu) and the +# RMG Team (rmg_dev@mit.edu) +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the 'Software'), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# +################################################################################ + +""" +Contains the :class:`SimpleReactor` class, providing a reaction system +consisting of a homogeneous, isothermal, isobaric batch reactor. +""" + +import numpy +cimport numpy +from pydas cimport DASSL +from base cimport ReactionSystem +cimport cython + +import rmgpy.constants as constants +cimport rmgpy.constants as constants +from rmgpy.quantity import Quantity +from rmgpy.quantity cimport ScalarQuantity, ArrayQuantity + +cdef class LiquidReactor(ReactionSystem): + """ + A reaction system consisting of a homogeneous, isothermal, constant volume batch + reactor. These assumptions allow for a number of optimizations that enable + this solver to complete very rapidly, even for large kinetic models. + """ + + cdef public ScalarQuantity T + cdef public ScalarQuantity V + cdef public ScalarQuantity P + + cdef public dict initialConcentrations + + cdef public numpy.ndarray reactantIndices + cdef public numpy.ndarray productIndices + cdef public numpy.ndarray networkIndices + cdef public numpy.ndarray forwardRateCoefficients + cdef public numpy.ndarray reverseRateCoefficients + cdef public numpy.ndarray networkLeakCoefficients + cdef public numpy.ndarray jacobianMatrix + + def __init__(self, T, initialConcentrations, termination): + ReactionSystem.__init__(self, termination) + self.T = Quantity(T) + self.P = Quantity(100000.,'kPa') # Arbitrary high pressure (1000 Bar) to get reactions in the high-pressure limit! + self.initialConcentrations = initialConcentrations # should be passed in SI + self.V = None # will be set from initialConcentrations in initializeModel + + # These are helper variables used within the solver + self.reactantIndices = None + self.productIndices = None + self.networkIndices = None + self.forwardRateCoefficients = None + self.reverseRateCoefficients = None + self.jacobianMatrix = None + + def convertInitalKeysToSpeciesObjects(self, speciesDict): + """ + Convert the initialConcentrations dictionary from species names into species objects, + using the given dictionary of species. + """ + initialConcentrations = {} + for label, moleFrac in self.initialConcentrations.iteritems(): + initialConcentrations[speciesDict[label]] = moleFrac + self.initialConcentrations = initialConcentrations + + cpdef initializeModel(self, list coreSpecies, list coreReactions, list edgeSpecies, list edgeReactions, list pdepNetworks=None, atol=1e-16, rtol=1e-8): + """ + Initialize a simulation of the simple reactor using the provided kinetic + model. + """ + + # First call the base class version of the method + # This initializes the attributes declared in the base class + ReactionSystem.initializeModel(self, coreSpecies, coreReactions, edgeSpecies, edgeReactions, pdepNetworks, atol, rtol) + + cdef int numCoreSpecies, numCoreReactions, numEdgeSpecies, numEdgeReactions, numPdepNetworks + cdef int i, j, l, index + cdef double V + cdef dict speciesIndex, reactionIndex + cdef numpy.ndarray[numpy.int_t, ndim=2] reactantIndices, productIndices, networkIndices + cdef numpy.ndarray[numpy.float64_t, ndim=1] forwardRateCoefficients, reverseRateCoefficients, networkLeakCoefficients + + pdepNetworks = pdepNetworks or [] + + numCoreSpecies = len(coreSpecies) + numCoreReactions = len(coreReactions) + numEdgeSpecies = len(edgeSpecies) + numEdgeReactions = len(edgeReactions) + numPdepNetworks = len(pdepNetworks) + + # Assign an index to each species (core first, then edge) + speciesIndex = {} + for index, spec in enumerate(coreSpecies): + speciesIndex[spec] = index + for index, spec in enumerate(edgeSpecies): + speciesIndex[spec] = index + numCoreSpecies + # Assign an index to each reaction (core first, then edge) + reactionIndex = {} + for index, rxn in enumerate(coreReactions): + reactionIndex[rxn] = index + for index, rxn in enumerate(edgeReactions): + reactionIndex[rxn] = index + numCoreReactions + + # Generate reactant and product indices + # Generate forward and reverse rate coefficients k(T,P) + reactantIndices = -numpy.ones((numCoreReactions + numEdgeReactions, 3), numpy.int ) + productIndices = -numpy.ones_like(reactantIndices) + forwardRateCoefficients = numpy.zeros((numCoreReactions + numEdgeReactions), numpy.float64) + reverseRateCoefficients = numpy.zeros_like(forwardRateCoefficients) + for rxnList in [coreReactions, edgeReactions]: + for rxn in rxnList: + j = reactionIndex[rxn] + forwardRateCoefficients[j] = rxn.getRateCoefficient(self.T.value_si, self.P.value_si) + reverseRateCoefficients[j] = forwardRateCoefficients[j] / rxn.getEquilibriumConstant(self.T.value_si) + for l, spec in enumerate(rxn.reactants): + i = speciesIndex[spec] + reactantIndices[j,l] = i + for l, spec in enumerate(rxn.products): + i = speciesIndex[spec] + productIndices[j,l] = i + + networkIndices = -numpy.ones((numPdepNetworks, 3), numpy.int ) + networkLeakCoefficients = numpy.zeros((numPdepNetworks), numpy.float64) + for j, network in enumerate(pdepNetworks): + networkLeakCoefficients[j] = network.getLeakCoefficient(self.T.value_si, self.P.value_si) + for l, spec in enumerate(network.source): + i = speciesIndex[spec] + networkIndices[j,l] = i + + self.reactantIndices = reactantIndices + self.productIndices = productIndices + self.forwardRateCoefficients = forwardRateCoefficients + self.reverseRateCoefficients = reverseRateCoefficients + self.networkIndices = networkIndices + self.networkLeakCoefficients = networkLeakCoefficients + + # Set initial conditions + t0 = 0.0 + y0 = numpy.zeros((numCoreSpecies), numpy.float64) + + for spec, conc in self.initialConcentrations.iteritems(): + self.coreSpeciesConcentrations[speciesIndex[spec]] = conc + V = 1.0 / numpy.sum(self.coreSpeciesConcentrations) + self.V = Quantity(V,'m^3') #: volume (m3) required to contain one mole total of core species at start + for j in range(y0.shape[0]): + y0[j] = self.coreSpeciesConcentrations[j] * V + + # Initialize the model + dydt0 = - self.residual(t0, y0, numpy.zeros((numCoreSpecies), numpy.float64))[0] + DASSL.initialize(self, t0, y0, dydt0, atol, rtol) + + cpdef writeWorksheetHeader(self, worksheet): + """ + Write some descriptive information about the reaction system to the + first two rows of the given `worksheet`. + """ + import xlwt + style0 = xlwt.easyxf('font: bold on') + worksheet.write(0, 0, 'Liquid Reactor', style0) + worksheet.write(1, 0, 'T = {0:g} K'.format(self.T.value_si)) + + @cython.boundscheck(False) + def residual(self, double t, numpy.ndarray[numpy.float64_t, ndim=1] y, numpy.ndarray[numpy.float64_t, ndim=1] dydt): + + """ + Return the residual function for the governing DAE system for the + simple reaction system. + """ + cdef numpy.ndarray[numpy.int_t, ndim=2] ir, ip, inet + cdef numpy.ndarray[numpy.float64_t, ndim=1] res, kf, kr, knet + cdef int numCoreSpecies, numCoreReactions, numEdgeSpecies, numEdgeReactions, numPdepNetworks + cdef int j, first, second, third + cdef double k, V, reactionRate + cdef numpy.ndarray[numpy.float64_t, ndim=1] coreSpeciesConcentrations, coreSpeciesRates, coreReactionRates, edgeSpeciesRates, edgeReactionRates, networkLeakRates + cdef numpy.ndarray[numpy.float64_t, ndim=1] C + + res = numpy.zeros(y.shape[0], numpy.float64) + + ir = self.reactantIndices + ip = self.productIndices + kf = self.forwardRateCoefficients + kr = self.reverseRateCoefficients + inet = self.networkIndices + knet = self.networkLeakCoefficients + + numCoreSpecies = len(self.coreSpeciesRates) + numCoreReactions = len(self.coreReactionRates) + numEdgeSpecies = len(self.edgeSpeciesRates) + numEdgeReactions = len(self.edgeReactionRates) + numPdepNetworks = len(self.networkLeakRates) + + coreSpeciesConcentrations = numpy.zeros_like(self.coreSpeciesConcentrations) + coreSpeciesRates = numpy.zeros_like(self.coreSpeciesRates) + coreReactionRates = numpy.zeros_like(self.coreReactionRates) + edgeSpeciesRates = numpy.zeros_like(self.edgeSpeciesRates) + edgeReactionRates = numpy.zeros_like(self.edgeReactionRates) + networkLeakRates = numpy.zeros_like(self.networkLeakRates) + + C = numpy.zeros_like(self.coreSpeciesConcentrations) + V = self.V.value_si # constant volume reactor + + for j in range(y.shape[0]): + C[j] = y[j] / V + coreSpeciesConcentrations[j] = C[j] + + for j in range(ir.shape[0]): + k = kf[j] + if ir[j,0] >= numCoreSpecies or ir[j,1] >= numCoreSpecies or ir[j,2] >= numCoreSpecies: + reactionRate = 0.0 + elif ir[j,1] == -1: # only one reactant + reactionRate = k * C[ir[j,0]] + elif ir[j,2] == -1: # only two reactants + reactionRate = k * C[ir[j,0]] * C[ir[j,1]] + else: # three reactants!! (really?) + reactionRate = k * C[ir[j,0]] * C[ir[j,1]] * C[ir[j,2]] + k = kr[j] + if ip[j,0] >= numCoreSpecies or ip[j,1] >= numCoreSpecies or ip[j,2] >= numCoreSpecies: + pass + elif ip[j,1] == -1: # only one reactant + reactionRate -= k * C[ip[j,0]] + elif ip[j,2] == -1: # only two reactants + reactionRate -= k * C[ip[j,0]] * C[ip[j,1]] + else: # three reactants!! (really?) + reactionRate -= k * C[ip[j,0]] * C[ip[j,1]] * C[ip[j,2]] + + # Set the reaction and species rates + if j < numCoreReactions: + # The reaction is a core reaction + coreReactionRates[j] = reactionRate + + # Add/substract the total reaction rate from each species rate + # Since it's a core reaction we know that all of its reactants + # and products are core species + first = ir[j,0] + coreSpeciesRates[first] -= reactionRate + second = ir[j,1] + if second != -1: + coreSpeciesRates[second] -= reactionRate + third = ir[j,2] + if third != -1: + coreSpeciesRates[third] -= reactionRate + first = ip[j,0] + coreSpeciesRates[first] += reactionRate + second = ip[j,1] + if second != -1: + coreSpeciesRates[second] += reactionRate + third = ip[j,2] + if third != -1: + coreSpeciesRates[third] += reactionRate + + else: + # The reaction is an edge reaction + edgeReactionRates[j-numCoreReactions] = reactionRate + + # Add/substract the total reaction rate from each species rate + # Since it's an edge reaction its reactants and products could + # be either core or edge species + # We're only interested in the edge species + first = ir[j,0] + if first >= numCoreSpecies: edgeSpeciesRates[first-numCoreSpecies] -= reactionRate + second = ir[j,1] + if second != -1: + if second >= numCoreSpecies: edgeSpeciesRates[second-numCoreSpecies] -= reactionRate + third = ir[j,2] + if third != -1: + if third >= numCoreSpecies: edgeSpeciesRates[third-numCoreSpecies] -= reactionRate + first = ip[j,0] + if first >= numCoreSpecies: edgeSpeciesRates[first-numCoreSpecies] += reactionRate + second = ip[j,1] + if second != -1: + if second >= numCoreSpecies: edgeSpeciesRates[second-numCoreSpecies] += reactionRate + third = ip[j,2] + if third != -1: + if third >= numCoreSpecies: edgeSpeciesRates[third-numCoreSpecies] += reactionRate + + for j in range(inet.shape[0]): + k = knet[j] + if inet[j,1] == -1: # only one reactant + reactionRate = k * C[inet[j,0]] + elif inet[j,2] == -1: # only two reactants + reactionRate = k * C[inet[j,0]] * C[inet[j,1]] + else: # three reactants!! (really?) + reactionRate = k * C[inet[j,0]] * C[inet[j,1]] * C[inet[j,2]] + networkLeakRates[j] = reactionRate + + self.coreSpeciesConcentrations = coreSpeciesConcentrations + self.coreSpeciesRates = coreSpeciesRates + self.coreReactionRates = coreReactionRates + self.edgeSpeciesRates = edgeSpeciesRates + self.edgeReactionRates = edgeReactionRates + self.networkLeakRates = networkLeakRates + + res = coreSpeciesRates * V - dydt + return res, 0 + + @cython.boundscheck(False) + def jacobian(self, double t, numpy.ndarray[numpy.float64_t, ndim=1] y, numpy.ndarray[numpy.float64_t, ndim=1] dydt, double cj): + """ + Return the analytical Jacobian for the reaction system. + """ + cdef numpy.ndarray[numpy.int_t, ndim=2] ir, ip + cdef numpy.ndarray[numpy.float64_t, ndim=1] kf, kr, C + cdef numpy.ndarray[numpy.float64_t, ndim=2] pd + cdef int numCoreReactions, j + cdef double k, deriv + + pd = -cj * numpy.identity(y.shape[0], numpy.float64) + ir = self.reactantIndices + ip = self.productIndices + kf = self.forwardRateCoefficients + kr = self.reverseRateCoefficients + numCoreReactions = len(self.coreReactionRates) + + V = self.V.value_si # constant volume reactor + + C = numpy.zeros_like(self.coreSpeciesConcentrations) + for j in range(y.shape[0]): + C[j] = y[j] / V + + for j in range(numCoreReactions): + + k = kf[j] + if ir[j,1] == -1: # only one reactant + deriv = k + pd[ir[j,0], ir[j,0]] -= deriv + + pd[ip[j,0], ir[j,0]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,0]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,0]] += deriv + + + elif ir[j,2] == -1: # only two reactants + if ir[j,0] == ir[j,1]: + deriv = 2 * k * C[ir[j,0]] + pd[ir[j,0], ir[j,0]] -= 2 * deriv + + pd[ip[j,0], ir[j,0]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,0]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,0]] += deriv + + else: + # Derivative with respect to reactant 1 + deriv = k * C[ir[j, 1]] + pd[ir[j,0], ir[j,0]] -= deriv + pd[ir[j,1], ir[j,0]] -= deriv + + pd[ip[j,0], ir[j,0]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,0]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,0]] += deriv + + # Derivative with respect to reactant 2 + deriv = k * C[ir[j, 0]] + pd[ir[j,0], ir[j,1]] -= deriv + pd[ir[j,1], ir[j,1]] -= deriv + + pd[ip[j,0], ir[j,1]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,1]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,1]] += deriv + + + else: # three reactants!! (really?) + if (ir[j,0] == ir[j,1] & ir[j,0] == ir[j,2]): + deriv = 3 * k * C[ir[j,0]] * C[ir[j,0]] + pd[ir[j,0], ir[j,0]] -= 3 * deriv + + pd[ip[j,0], ir[j,0]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,0]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,0]] += deriv + + elif ir[j,0] == ir[j,1]: + # derivative with respect to reactant 1 + deriv = 2 * k * C[ir[j,0]] * C[ir[j,2]] + pd[ir[j,0], ir[j,0]] -= 2 * deriv + pd[ir[j,2], ir[j,0]] -= deriv + + pd[ip[j,0], ir[j,0]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,0]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,0]] += deriv + + # derivative with respect to reactant 3 + deriv = k * C[ir[j,0]] * C[ir[j,0]] + pd[ir[j,0], ir[j,2]] -= 2 * deriv + pd[ir[j,2], ir[j,2]] -= deriv + + pd[ip[j,0], ir[j,2]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,2]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,2]] += deriv + + + elif ir[j,1] == ir[j,2]: + # derivative with respect to reactant 1 + deriv = k * C[ir[j,1]] * C[ir[j,1]] + pd[ir[j,0], ir[j,0]] -= deriv + pd[ir[j,1], ir[j,0]] -= 2 * deriv + + pd[ip[j,0], ir[j,0]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,0]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,0]] += deriv + + # derivative with respect to reactant 2 + deriv = 2 * k * C[ir[j,0]] * C[ir[j,1]] + pd[ir[j,0], ir[j,1]] -= deriv + pd[ir[j,1], ir[j,1]] -= 2 * deriv + + pd[ip[j,0], ir[j,1]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,1]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,1]] += deriv + + else: + # derivative with respect to reactant 1 + deriv = k * C[ir[j,1]] * C[ir[j,2]] + pd[ir[j,0], ir[j,0]] -= deriv + pd[ir[j,1], ir[j,0]] -= deriv + pd[ir[j,2], ir[j,0]] -= deriv + + pd[ip[j,0], ir[j,0]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,0]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,0]] += deriv + + # derivative with respect to reactant 2 + deriv = k * C[ir[j,0]] * C[ir[j,2]] + pd[ir[j,0], ir[j,1]] -= deriv + pd[ir[j,1], ir[j,1]] -= deriv + pd[ir[j,2], ir[j,1]] -= deriv + + pd[ip[j,0], ir[j,1]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,1]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,1]] += deriv + + # derivative with respect to reactant 3 + deriv = k * C[ir[j,0]] * C[ir[j,1]] + pd[ir[j,0], ir[j,2]] -= deriv + pd[ir[j,1], ir[j,2]] -= deriv + pd[ir[j,2], ir[j,2]] -= deriv + + pd[ip[j,0], ir[j,2]] += deriv + if ip[j,1] != -1: + pd[ip[j,1], ir[j,2]] += deriv + if ip[j,2] != -1: + pd[ip[j,2], ir[j,2]] += deriv + + + + k = kr[j] + if ip[j,1] == -1: # only one reactant + deriv = k + pd[ip[j,0], ip[j,0]] -= deriv + + pd[ir[j,0], ip[j,0]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,0]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,0]] += deriv + + + elif ip[j,2] == -1: # only two reactants + if ip[j,0] == ip[j,1]: + deriv = 2 * k * C[ip[j,0]] + pd[ip[j,0], ip[j,0]] -= 2 * deriv + + pd[ir[j,0], ip[j,0]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,0]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,0]] += deriv + + else: + # Derivative with respect to reactant 1 + deriv = k * C[ip[j, 1]] + pd[ip[j,0], ip[j,0]] -= deriv + pd[ip[j,1], ip[j,0]] -= deriv + + pd[ir[j,0], ip[j,0]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,0]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,0]] += deriv + + # Derivative with respect to reactant 2 + deriv = k * C[ip[j, 0]] + pd[ip[j,0], ip[j,1]] -= deriv + pd[ip[j,1], ip[j,1]] -= deriv + + pd[ir[j,0], ip[j,1]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,1]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,1]] += deriv + + + else: # three reactants!! (really?) + if (ip[j,0] == ip[j,1] & ip[j,0] == ip[j,2]): + deriv = 3 * k * C[ip[j,0]] * C[ip[j,0]] + pd[ip[j,0], ip[j,0]] -= 3 * deriv + + pd[ir[j,0], ip[j,0]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,0]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,0]] += deriv + + elif ip[j,0] == ip[j,1]: + # derivative with respect to reactant 1 + deriv = 2 * k * C[ip[j,0]] * C[ip[j,2]] + pd[ip[j,0], ip[j,0]] -= 2 * deriv + pd[ip[j,2], ip[j,0]] -= deriv + + pd[ir[j,0], ip[j,0]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,0]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,0]] += deriv + + # derivative with respect to reactant 3 + deriv = k * C[ip[j,0]] * C[ip[j,0]] + pd[ip[j,0], ip[j,2]] -= 2 * deriv + pd[ip[j,2], ip[j,2]] -= deriv + + pd[ir[j,0], ip[j,2]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,2]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,2]] += deriv + + + elif ip[j,1] == ip[j,2]: + # derivative with respect to reactant 1 + deriv = k * C[ip[j,1]] * C[ip[j,1]] + pd[ip[j,0], ip[j,0]] -= deriv + pd[ip[j,1], ip[j,0]] -= 2 * deriv + + pd[ir[j,0], ip[j,0]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,0]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,0]] += deriv + + # derivative with respect to reactant 2 + deriv = 2 * k * C[ip[j,0]] * C[ip[j,1]] + pd[ip[j,0], ip[j,1]] -= deriv + pd[ip[j,1], ip[j,1]] -= 2 * deriv + + pd[ir[j,0], ip[j,1]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,1]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,1]] += deriv + + else: + # derivative with respect to reactant 1 + deriv = k * C[ip[j,1]] * C[ip[j,2]] + pd[ip[j,0], ip[j,0]] -= deriv + pd[ip[j,1], ip[j,0]] -= deriv + pd[ip[j,2], ip[j,0]] -= deriv + + pd[ir[j,0], ip[j,0]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,0]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,0]] += deriv + + # derivative with respect to reactant 2 + deriv = k * C[ip[j,0]] * C[ip[j,2]] + pd[ip[j,0], ip[j,1]] -= deriv + pd[ip[j,1], ip[j,1]] -= deriv + pd[ip[j,2], ip[j,1]] -= deriv + + pd[ir[j,0], ip[j,1]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,1]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,1]] += deriv + + # derivative with respect to reactant 3 + deriv = k * C[ip[j,0]] * C[ip[j,1]] + pd[ip[j,0], ip[j,2]] -= deriv + pd[ip[j,1], ip[j,2]] -= deriv + pd[ip[j,2], ip[j,2]] -= deriv + + pd[ir[j,0], ip[j,2]] += deriv + if ir[j,1] != -1: + pd[ir[j,1], ip[j,2]] += deriv + if ir[j,2] != -1: + pd[ir[j,2], ip[j,2]] += deriv + + self.jacobianMatrix = pd + cj * numpy.identity(y.shape[0], numpy.float64) + return pd + + @cython.boundscheck(False) + def computeRateDerivative(self): + """ + Returns derivative vector df/dk_j where dc/dt = f(c, t, k) and + k_j is the rate parameter for the jth core reaction. + """ + cdef numpy.ndarray[numpy.int_t, ndim=2] ir, ip + cdef numpy.ndarray[numpy.float64_t, ndim=1] y, kf, kr + cdef numpy.ndarray[numpy.float64_t, ndim=2] rateDeriv + cdef double fderiv, rderiv, flux + cdef int j, numCoreReactions + + ir = self.reactantIndices + ip = self.productIndices + + kf = self.forwardRateCoefficients + kr = self.reverseRateCoefficients + y = self.coreSpeciesConcentrations + + + numCoreReactions = len(self.coreReactionRates) + + rateDeriv = numpy.zeros((y.shape[0],numCoreReactions), numpy.float64) + + for j in range(numCoreReactions): + if ir[j,1] == -1: # only one reactant + fderiv = y[ir[j,0]] + elif ir[j,2] == -1: # only two reactants + fderiv = y[ir[j,0]] * y[ir[j,1]] + else: # three reactants!! (really?) + fderiv = y[ir[j,0]] * y[ir[j,1]] * y[ir[j,2]] + + if ip[j,1] == -1: # only one reactant + rderiv = kr[j] / kf [j] * y[ip[j,0]] + elif ip[j,2] == -1: # only two reactants + rderiv = kr[j] / kf [j] * y[ip[j,0]] * y[ip[j,1]] + else: # three reactants!! (really?) + rderiv = kr[j] / kf [j] * y[ip[j,0]] * y[ip[j,1]] * y[ip[j,2]] + + + flux = fderiv - rderiv + if ir[j,1] == -1: + rateDeriv[ir[j,0], j] -= flux + elif ir[j,2] == -1: + rateDeriv[ir[j,0], j] -= flux + rateDeriv[ir[j,1], j] -= flux + else: + rateDeriv[ir[j,0], j] -= flux + rateDeriv[ir[j,1], j] -= flux + rateDeriv[ir[j,2], j] -= flux + + if ip[j,1] == -1: + rateDeriv[ip[j,0], j] += flux + elif ip[j,2] == -1: + rateDeriv[ip[j,0], j] += flux + rateDeriv[ip[j,1], j] += flux + else: + rateDeriv[ip[j,0], j] += flux + rateDeriv[ip[j,1], j] += flux + rateDeriv[ip[j,2], j] += flux + + return rateDeriv diff --git a/rmgpy/solver/simple.pyx b/rmgpy/solver/simple.pyx index f759e221cd5..18a5f5f038f 100644 --- a/rmgpy/solver/simple.pyx +++ b/rmgpy/solver/simple.pyx @@ -77,6 +77,16 @@ cdef class SimpleReactor(ReactionSystem): self.forwardRateCoefficients = None self.reverseRateCoefficients = None self.jacobianMatrix = None + + def convertInitalKeysToSpeciesObjects(self, speciesDict): + """ + Convert the initialMoleFractions dictionary from species names into species objects, + using the given dictionary of species. + """ + initialMoleFractions = {} + for label, moleFrac in self.initialMoleFractions.iteritems(): + initialMoleFractions[speciesDict[label]] = moleFrac + self.initialMoleFractions = initialMoleFractions cpdef initializeModel(self, list coreSpecies, list coreReactions, list edgeSpecies, list edgeReactions, list pdepNetworks=None, atol=1e-16, rtol=1e-8): """ @@ -825,4 +835,4 @@ cdef class SimpleReactor(ReactionSystem): for j in range(numCoreReactions): if c[i] != 0.0: norm[i,j] = kf[j]/c[i] - return norm \ No newline at end of file + return norm diff --git a/setup.py b/setup.py index d05d3d87d16..25080d943a1 100644 --- a/setup.py +++ b/setup.py @@ -118,6 +118,7 @@ def getSolverExtensionModules(): return [ Extension('rmgpy.solver.base', ['rmgpy/solver/base.pyx'], include_dirs=['.']), Extension('rmgpy.solver.simple', ['rmgpy/solver/simple.pyx'], include_dirs=['.']), + Extension('rmgpy.solver.liquid', ['rmgpy/solver/liquid.pyx'], include_dirs=['.']), ] def getCanthermExtensionModules():