search

ReactionMechanismGenerator / RMG-Py

Python version of the amazing Reaction Mechanism Generator (RMG).
http://reactionmechanismgenerator.github.io/RMG-Py/
Other
400 stars 228 forks source link

Can't cope with missing bath gas when running in PDep mode. #37

Closed rwest closed 12 years ago

rwest commented 13 years ago

I am running with:

# Data sources
database(
    '../RMG-database/input',
    thermoLibraries = ['primaryThermoLibrary','DFT_QCI_thermo','GRI-Mech3.0'],
    reactionLibraries = ['Methylformate','Glarborg/highP'],
    seedMechanisms = ['GRI-Mech3.0'],
)

pressureDependence(
    method='modified strong collision',
    minimumGrainSize=(2.0,'kcal/mol'),
    minimumNumberOfGrains=200,
    temperatures=(300,'K',2000,'K',8),
    pressures=(0.01,'bar',100,'bar',5),
    interpolation=('Chebyshev', 6, 4),
)

I have a lot of lines saying:

Calculating densities of states for Network #1...
Calculating phenomenological rate coefficients for Network #1...
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
Exception ZeroDivisionError: 'float division' in 'rmgpy.measure.collision.calculateCollisionFrequency' ignored
...etc...

then it dies with:

Traceback (most recent call last):
  File "rmg.py", line 711, in <module>
    execute(args)
  File "rmg.py", line 275, in execute
    reactionModel.enlarge(spec, database)
  File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 938, in enlarge
    self.updateUnimolecularReactionNetworks(database)
  File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 1399, in updateUnimolecularReactionNetworks
    network.update(self, database)
  File "/Users/rwest/XCodeProjects/RMGpy/RMG-Py/rmgpy/rmg/model.py", line 532, in update
    netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)
  File "reaction.pyx", line 332, in rmgpy.measure.reaction.fitInterpolationModel (build/pyrex/rmgpy/measure/reaction.c:5141)
  File "kinetics.py", line 732, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:12744)
  File "kinetics.py", line 767, in rmgpy.kinetics.Chebyshev.fitToData (build/pyrex/rmgpy/kinetics.c:12358)
ValueError: math domain error

A little up the stack trace is :

 netReaction.kinetics = fitInterpolationModel(netReaction, Tlist, Plist, K[:,:,i,j], model, Tmin, Tmax, Pmin, Pmax)

at which point the matrix K is entirely zero. If it helps, this is the netReaction:

Reaction(index=334, reactants=[Species(index=1, label="methylformate", thermo=MultiNASA(Tmin=100, Tmax=5000, polynomials=[NASA(Tmin=(100,"K"), Tmax=(1133.91,"K"), coeffs=[3.26245,0.0112056,1.91547e-05,-2.51712e-08,8.14166e-12,-44760.2,12.148], comment="""Low temperature range polynomial"""),NASA(Tmin=(1133.91,"K"), Tmax=(5000,"K"), coeffs=[5.13303,0.0182987,-8.34096e-06,1.64356e-09,-1.18236e-13,-46064.6,-0.992883], comment="""High temperature range polynomial""")], comment="""NASA function fitted to Wilhoit function. Weighted RMS error = 0.042_R;(Unweighted) RMS error = 0.046_R;"""), states=StatesModel(modes=[HarmonicOscillator(frequencies=([2782.5,750,1395,475,1775,1000,2750,2800,2850,1350,1500,750,1050,1375,1000,425],"cm^-1")), HinderedRotor(inertia=(0.0877936,"amu_angstrom^2"), barrier=(2.01855,"kJ/mol"), symmetry=1), HinderedRotor(inertia=(0.0654664,"amu_angstrom^2"), barrier=(47.6296,"kJ/mol"), symmetry=1)], spinMultiplicity=1), molecule=[Molecule(SMILES="COC=O")], E0=(-372.379,"kJ/mol"), lennardJones=LennardJones(sigma=(4.687e-10,"m"), epsilon=(7.33678e-21,"J")), molecularWeight=(60.052,"g/mol"))], products=[Species(index=16, label="CH3OH", thermo=MultiNASA(Tmin=100, Tmax=5000, polynomials=[NASA(Tmin=(100,"K"), Tmax=(1143.22,"K"), coeffs=[3.89247,-0.000484239,2.42754e-05,-2.44333e-08,7.4633e-12,-25699.7,5.86289], comment="""Low temperature range polynomial"""),NASA(Tmin=(1143.22,"K"), Tmax=(5000,"K"), coeffs=[5.08762,0.00759196,-2.40448e-06,5.03905e-10,-4.09468e-14,-26774,-3.56642], comment="""High temperature range polynomial""")], comment="""NASA function fitted to Wilhoit function. Weighted RMS error = 0.028_R;(Unweighted) RMS error = 0.030_R;"""), molecule=[Molecule(SMILES="CO")], E0=(-213.698,"kJ/mol"), lennardJones=LennardJones(sigma=(4.443e-10,"m"), epsilon=(1.52838e-21,"J")), molecularWeight=(32.0419,"g/mol")), Species(index=25, label="CO", thermo=MultiNASA(Tmin=100, Tmax=5000, polynomials=[NASA(Tmin=(100,"K"), Tmax=(1552.61,"K"), coeffs=[3.65108,-0.00132484,3.28978e-06,-2.09102e-09,4.32436e-13,-14359.9,3.18576], comment="""Low temperature range polynomial"""),NASA(Tmin=(1552.61,"K"), Tmax=(5000,"K"), coeffs=[3.09568,0.00139239,-5.78124e-07,1.03416e-10,-6.83094e-15,-14342.5,5.61029], comment="""High temperature range polynomial""")], comment="""NASA function fitted to Wilhoit function. Weighted RMS error = 0.006_R;(Unweighted) RMS error = 0.003_R;"""), molecule=[Molecule(SMILES="[C]=O")], E0=(-119.306,"kJ/mol"), lennardJones=LennardJones(sigma=(4.443e-10,"m"), epsilon=(1.52838e-21,"J")), molecularWeight=(28.0101,"g/mol"))], reversible=False)

(edit: title changed after comment below added)

rwest commented 13 years ago

I had no collider bath gas. The problem went away when I added 90 mol% N2, defined as:

species(
    label='N2',
    reactive=False,
    structure=InChI("InChI=1/N2/c1-2"),
)

I think some sanity check like "do you have a bath gas" could help diagnose this.