PDep rates don't work in Reaction or Kinetics Libraries

[gmagoon]

this essentially on the latest master branch:

Reaction Set Found from Reaction Library [JP10(1) --> MA111(10)]
Reaction Set Found from Reaction Library [C4H3(54) --> C2H2(51) + C2H(41)]
RMG is only utilizing the high-pressure limit parameters for PKL reaction: C4H3(54) --> C4H2(44) + H(38)
RMG is only utilizing the high-pressure limit parameters for PKL reaction: C4H3(54) --> C4H2(44) + H(38)
RMG is only utilizing the high-pressure limit parameters for PKL reaction: C4H3(54) --> C4H2(44) + H(38)
RMG is only utilizing the high-pressure limit parameters for PKL reaction: C4H3(54) --> C4H2(44) + H(38)
Created new intra_H_migration reaction: C4H3(54) --> C4H3J(707)
Created new intra_H_migration reaction: C4H3(54) --> C4H3J(163)
ERROR: java.lang.NullPointerException
    at jing.rxn.PDepNetwork.addReactionToNetworks(PDepNetwork.java:681)
    at jing.rxnSys.ReactionModelGenerator.initializeCoreEdgeModel(ReactionModelGenerator.java:3643)
    at jing.rxnSys.ReactionModelGenerator.initializeCoreEdgeReactionModel(ReactionModelGenerator.java:3867)
    at jing.rxnSys.ReactionModelGenerator.modelGeneration(ReactionModelGenerator.java:1302)
    at RMG.main(RMG.java:96)

Exception in thread "main" java.lang.NullPointerException
        at jing.rxnSys.Logger.log(Logger.java:160)
        at jing.rxnSys.Logger.critical(Logger.java:204)
        at RMG.main(RMG.java:106)

[rwest]

When does this occur?

[gmagoon]

n.b....I have Chebyshevs and PLOGs in pdepreactions.txt of SeedMechanism, ReactionLibrary, and KineticsLibrary

[gmagoon]

This occurs a couple minutes into a PDep RMG run.

[gmagoon]

For anyone interested, I'm working with this condition file:

Database: RMG_database

MaxCarbonNumberPerSpecies: 10 
MaxOxygenNumberPerSpecies: 2
MaxRadicalNumberPerSpecies: 2
//MaxSulfurNumberPerSpecies: 
//MaxSiliconNumberPerSpecies: 
//MaxHeavyAtomPerSpecies:

PrimaryThermoLibrary:
Name: GRIMech3.0
Location: GRI-Mech3.0
Name: RMG-minimal
Location: primaryThermoLibrary
Name: JP10
Location: JP10
END

PrimaryTransportLibrary:
Name: GRIMech3.0
Location: GRI-Mech3.0
END

ForbiddenStructures:
cyclopropeneEndoDB
1 C 0 {2,S} {3,D}
2 C 0 {1,S} {3,S}
3 C 0 {1,D} {2,S}                  

cyclopropeneExoDB
1 C 0 {2,S} {3,S} {4,D}    
2 C 0 {1,S} {3,S}   
3 C 0 {1,S} {2,S}
4 C 0 {1,D}

END

ReadRestart: no
WriteRestart: yes

TemperatureModel: Constant (K) 1250 1500 2000
//note that pressure is slightly different than original mechanism: 1 atm, vs. 1 bar
PressureModel: Constant (atm) 1 10

ThermoMethod: QM both
QMForCyclicsOnly: on
MaxRadNumForQM: 0

InitialStatus:

JP10 (mol/cm3) 1.0 1.0 1.0
1 C 0 {2,S} {9,S} {10,S}
2 C 0 {1,S} {3,S} {6,S}
3 C 0 {2,S} {4,S}
4 C 0 {3,S} {5,S}
5 C 0 {4,S} {6,S}
6 C 0 {5,S} {7,S} {2,S}
7 C 0 {6,S} {8,S} {10,S}
8 C 0 {7,S} {9,S}
9 C 0 {8,S} {1,S}
10 C 0 {1,S} {7,S}

O2 (mol/cm3) 7.0 14.0 28.0
1 O 1 {2,S}
2 O 1 {1,S}

END

InertGas:
N2 (mol/cm3) 52.67 52.67 52.67
//Ar (mol/cm3) 0.0e-6
END


//SpectroscopicDataEstimator: off
//PressureDependence: off
SpectroscopicDataEstimator: FrequencyGroups
PressureDependence: ReservoirState
PDepKineticsModel: Chebyshev 4 4
TRange: (K) 300.0 3000.0 8
PRange: (bar) 0.01 100.0 5

FinishController:
(1) Goal ReactionTime: 0.001 (s)
(2) Error Tolerance: 0.15

DynamicSimulator: DASSL
TimeStep: AUTO
Atol: 1e-18
Rtol: 1e-8

PrimaryKineticLibrary:
Name: JP10kinlib
Location: JP10kinlib
END

ReactionLibrary:
Name: JP10rxnlib
Location: JP10rxnlib
END

SeedMechanism:
Name: Leeds
Location: combustion_core/version5
GenerateReactions: yes
END

ChemkinUnits:
A: moles
Ea: kcal/mol

And I'm working on the JP10_2011 branch of the pharos git repository /home/gmagoon/RMG-java/.git at commit 2d6071e624413fce1f6eabe9703bec156e625b3e

[gmagoon]

Not sure if this is related (maybe it should be a separate issue?) but when I run a NON-Pdep job, it crashes with the following at the first ODE solve, suggesting it has something to do with my PLOG or Chebyshev reactions:

Ea raised by 12.6 from 9.5 to dHrxn(298K)=22.1 kcal/mol.
Created new Intra_Disproportionation reaction: C3H2(58) --> C3H2(1225)

The model core has 346 reactions and 42 species.
The model edge has 4388 reactions and 1172 species.
The chem.inp file has 42 species (excluding inert gases).
The chem.inp file has 355 reactions (excluding duplicates).
Running time: 2.006 min
Memory used: 54.04 MB / 266.73 MB (20.26%)

Solving reaction system...
ERROR: java.lang.NullPointerException
    at jing.rxnSys.Logger.log(Logger.java:160)
    at jing.rxnSys.Logger.error(Logger.java:213)
    at jing.rxn.PDepReaction.calculateRate(PDepReaction.java:370)
    at jing.rxnSys.JDAS.transferReaction(JDAS.java:652)
    at jing.rxnSys.JDAS.generatePDepODEReactionList(JDAS.java:167)
    at jing.rxnSys.JDASSL.solve(JDASSL.java:161)
    at jing.rxnSys.ReactionSystem.solveReactionSystem(ReactionSystem.java:1499)
    at jing.rxnSys.ReactionModelGenerator.modelGeneration(ReactionModelGenerator.java:1395)
    at RMG.main(RMG.java:96)

Exception in thread "main" java.lang.NullPointerException
        at jing.rxnSys.Logger.log(Logger.java:160)
        at jing.rxnSys.Logger.critical(Logger.java:204)
        at RMG.main(RMG.java:106)

[gmagoon]

Per discussion with rwest, it sounds like this issue may be related to the aforementioned use of PLOG/Chebyshev in the ReactionLibrary and KineticsLibrary, which is likely not supported

[nickvandewiele]

I tried P-dep reactions in KineticsLibrary alone and it gave no error.

I tried the same library folder as ReactionLibrary alone, and it gave the aforementioned error.

Ergo, the bug must be related to the treatment of ReactionLibrary.

[rwest]

The real issue here is that pressure dependent reaction rates from kinetics libraries do not work in ReactionLibrary or KineticsLibrary. (I think it never has worked - it just didn't used to fail in this manner). If you want to use ReactionLibrary or KineticsLibrary, don't have any pdep reactions (all the kinetics are passed to FAME for pressure dependent calculation). For this type of use there is a databases/RMG_database/kinetics_libraries/Glarborg/highP database which has high-pressure rate expressions from Glarborg/C3. Feel free to create others if you like to use Reaction Libraries. If you want to use pressure-dependent rate expressions, you have to use them in a SeedMechanism.

The specific issue of a null pointer exception (PDepNetwork.java:681) is because the reaction in question does not have a reverse, so the line a few lines before the exception sets reaction=null. If you want to investigate further, you can add the line

r.generateReverseReaction();

on line 452 of SeedMechanism.java https://github.com/GreenGroup/RMG-Java/blob/master/source/RMG/jing/rxnSys/SeedMechanism.java#L452 but that seems to just cause problems elsewhere. As any amount of fixing here will not actually make PDep rates work in ReactionLibraries, I personally am not pursuing it further for now.

[gmagoon]

I forgot about this bug and retried again with PLOG in ReactionLibrary...received the following error:

Reaction Set Found from Reaction Library [H2CCCH(53) --> H(41) + C3H2(61), H2CCCH(53) + H2CCCH(53) (included =true) --> C*CC(C#C)*C(87) (included =false), H2CCCH(53) + H2CCCH(53) (included =true) --> FULVENE(21) (included =false), H2CCCH(53) + H2CCCH(53) (included =true) --> BENZENE(69) (included =false), H2CCCH(53) + H2CCCH(53) (included =true) --> PHENYL(70) + H(41) (included =true)]
Created new reverse intra_H_migration reaction: H2CCCH(53) --> C3H3J(515)
ERROR: java.lang.NullPointerException
    at jing.rxn.PDepNetwork.addReactionToNetworks(PDepNetwork.java:681)
    at jing.rxn.PDepNetwork.makeIsomerIncluded(PDepNetwork.java:417)
    at jing.rxn.PDepNetwork.addReactionToNetworks(PDepNetwork.java:783)
    at jing.rxnSys.ReactionModelGenerator.initializeCoreEdgeModel(ReactionModelGenerator.java:3649)
    at jing.rxnSys.ReactionModelGenerator.initializeCoreEdgeReactionModel(ReactionModelGenerator.java:3873)
    at jing.rxnSys.ReactionModelGenerator.modelGeneration(ReactionModelGenerator.java:1308)
    at RMG.main(RMG.java:96)


RMG execution terminated at 2011-09-04 19:48:16

[rwest]

When reading a reaction library which contains PDep reactions we should do something sensible such as:

a) warn the user and discard the reactions b) warn the user and just take the high-P form (in seed mechanisms that I've made this is stored on the same line as the reaction, but some people just put random numbers or 1 0 0 here, as chemkin disregards it) c) stop and tell the user to remove the reactions or put them in a seed mechanism

Or, we should decide how PDep and ReactionLibraries work (at the moment all ReactionLibrary rates are sent as input to FAME, hence they must be high-P limits)

We could also make "high-P" versions of all the default seed mechanisms, such as I had done for Glarborg.

[rwest]

bump! I just ran into this issue again myself. I'm glad I decided to search the issue tracker, because the stack trace is most uninformative.

[rwest]

I have just changed the title of this issue to be more informative. Franklin recently requested this functionality again. I can't remember how hard we decided it would be to fix.