ReactionMechanismGenerator / RMG-Py

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

Can rmg.py generate mechanism for reactive Nitrogen-containing molecules ? #645

Closed OUchemical closed 8 years ago

OUchemical commented 8 years ago

Hi guys, I'm trying to generate n-phenyl-1-naphthylamines liquid phase oxidation mechanism in hexane solvent, and rmg reported that "PAN" is global forbidden species, so I went back to "ForbiddenStructures.py" and commented out species "N_birad_singlet_2singleBonds", systems started reporting error message: "Invalid valency for atom C with 0 unpaired electrons, 0 pairs of electrons, and 0 charge. "

I'm wondering if rmg.py can generate mechanism for reactive N-containing molecules, or seed mechanisms files or reaction libraries need to be pre-defined by the user in order to proceed?

Error messages are pasted below:

Traceback (most recent call last): File "C:\Anaconda2\RMGPycode\RMG-Py\rmg.py", line 165, in <module> rmg.execute(inputFile, output_dir, **kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\main.py", line 522, in execute self.initialize(inputFile, output_directory, **kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\main.py", line 434, in initialize raise ForbiddenStructureException("Input species {0} is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.".format(spec.label)) rmgpy.data.base.ForbiddenStructureException: Input species PAN is globally forbidden. You may explicitly allow it, but it will remain inert unless found in a seed mechanism or reaction library.

Generating thermodynamics for new species... 1 radicals on [C]1=CC=CC2=CC=CC=C12 exceeds limit of 0. Using HBI method. Reading existing thermo file QMfiles\pm3\UFWIBTONFRDIAS-UHFFFAOYSA.thermo Traceback (most recent call last): File "C:\Anaconda2\RMGPycode\RMG-Py\rmg.py", line 165, in <module> rmg.execute(inputFile, output_dir, **kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\main.py", line 549, in execute bimolecularReact=self.bimolecularReact) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 752, in enlarge spec.generateThermoData(database, quantumMechanics=self.quantumMechanics) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 98, in generateThermoData thermo0 = database.thermo.getThermoData(self, trainingSet=None, quantumMechanics=quantumMechanics) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\thermo.py", line 755, in getThermoData tdata = self.estimateRadicalThermoViaHBI(molecule, quantumMechanics.getThermoData) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\thermo.py", line 1056, in estimateRadicalThermoViaHBI thermoData_sat = stableThermoEstimator(saturatedStruct) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\qm\main.py", line 194, in getThermoData thermo0 = qm_molecule_calculator.generateThermoData() File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\qm\molecule.py", line 324, in generateThermoData if self.loadThermoData(): File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\qm\molecule.py", line 374, in loadThermoData loadedMolecule = rmgpy.molecule.Molecule().fromAdjacencyList(local_context['adjacencyList']) File "rmgpy\molecule\molecule.py", line 1183, in rmgpy.molecule.molecule.Molecule.fromAdjacencyList (build\Release\pyrex\rmgpy\molecule\mo lecule.c:22764) File "rmgpy\molecule\molecule.py", line 1191, in rmgpy.molecule.molecule.Molecule.fromAdjacencyList (build\Release\pyrex\rmgpy\molecule\mo lecule.c:22339) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\molecule\adjlist.py", line 710, in fromAdjacencyList for atom in atoms: ConsistencyChecker.check_partial_charge(atom) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\molecule\adjlist.py", line 113, in check_partial_charge .format(symbol=atom.symbol, radicals=atom.radicalElectrons, lonePairs=atom.lonePairs, charge=atom.charge)) rmgpy.molecule.adjlist.InvalidAdjacencyListError: Invalid valency for atom C with 0 unpaired electrons, 0 pairs of electrons, and 0 charge.

alongd commented 8 years ago

Hi OUchemical,

Yes, RMG can generate mechanisms for reacting N species. Can you attach your input file so I can regenerate your error?

OUchemical commented 8 years ago

Thanks a lot alongd! Here's my input file: input.txt

alongd commented 8 years ago

Try using the following forbiddenStructures.py file: https://github.com/alongd/RMG-database/blob/nitrogen/input/forbiddenStructures.py

(see issue #514 for details on commented structures)

OUchemical commented 8 years ago

I just tried using that forbiddenStructures file, and now I'm getting this error:

Error: <Molecule "C12(C3C4(C2=CC1C=C4)C=C3)Nc1ccccc1"> Error: <Molecule "c12c(Nc3ccccc3)cccc1cccc2"> Traceback (most recent call last): File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\scoop_framework\util.py", line 108, in call return self.myfn(_args, *_kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\react.py", line 83, in reactMolecules rxns = family.generateReactions(molecules) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\kinetics\family.py", line 1347, in generateReactions reactionList.extend(self.generateReactions(reactants, forward=False)) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\kinetics\family.py", line 1566, in generateReactions reaction.degeneracy = self.calculateDegeneracy(reaction) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\kinetics\family.py", line 1364, in calculateDegeneracy 'but generated {2}').format(reaction, self.label, len(reactions))) KineticsError: Unable to calculate degeneracy for reaction <Molecule "C12(C3C4(C2=CC1C=C4)C=C3)Nc1ccccc1"> <=> <Molecule "c12c(Nc3ccccc3)ccc c1cccc2"> in reaction family Intra_Diels_alder. Expected 1 reaction but generated 2

Traceback (most recent call last): File "C:\Anaconda2\RMGPycode\RMG-Py\rmg.py", line 165, in rmg.execute(inputFile, output_dir, _kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\main.py", line 549, in execute bimolecularReact=self.bimolecularReact) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 731, in enlarge reactions = list(react(self.core.species[i].copy(deep=True))) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\react.py", line 62, in react combos File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\scoopframework\util.py", line 150, in map return map(_args, _kwargs) File "C:\Anaconda2\envs\rmg_env\lib\site-packages\scoop-0.7.2.0-py2.7.egg\scoop\fallbacks.py", line 49, in wrapper File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\scoop_framework\util.py", line 108, in call return self.myfn(_args, **kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\react.py", line 83, in reactMolecules rxns = family.generateReactions(molecules) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\kinetics\family.py", line 1347, in generateReactions reactionList.extend(self.generateReactions(reactants, forward=False)) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\kinetics\family.py", line 1566, in generateReactions reaction.degeneracy = self.calculateDegeneracy(reaction) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\kinetics\family.py", line 1364, in calculateDegeneracy 'but generated {2}').format(reaction, self.label, len(reactions))) rmgpy.data.kinetics.common.KineticsError: Unable to calculate degeneracy for reaction <Molecule "C12(C3C4(C2=CC1C=C4)C=C3)Nc1ccccc1"> <=> <M olecule "c12c(Nc3ccccc3)cccc1cccc2"> in reaction family Intra_Diels_alder. Expected 1 reaction but generated 2

alongd commented 8 years ago

See issue #617 where a similar error was reported. Try replacing RMG-Py/rmgpy/data/kinetics/family.py with the following file: https://github.com/alongd/RMG-Py/blob/nitrogen/rmgpy/data/kinetics/family.py

OUchemical commented 8 years ago

Now getting this error:

Traceback (most recent call last): File "C:\Anaconda2\RMGPycode\RMG-Py\rmg.py", line 165, in rmg.execute(inputFile, output_dir, **kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\main.py", line 549, in execute bimolecularReact=self.bimolecularReact) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 733, in enlarge self.processNewReactions(reactions, self.core.species[i], None) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 830, in processNewReactions rxn, isNew = self.makeNewReaction(rxn) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 568, in makeNewReaction reactants = [self.makeNewSpecies(reactant)[0] for reactant in forward.reactants] File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 424, in makeNewSpecies molecule.clearLabeledAtoms() AttributeError: 'rmgpy.species.Species' object has no attribute 'clearLabeledAtoms'

was there something I needed to change alongside family.py?

alongd commented 8 years ago

What you're describing was solved back in 2011 in issue #35 . I don't suppose your RMG version was installed prior to 2011 and wasn't updated since, right? :)

I don't think I can help with this, maybe a more experienced user. @KEHANG , do you have any insights on the error above?

BTW, @OUchemical , I had no problem running your input file after implementing the fixes we discussed above, so worse case you can download my nitrogen branch using git.

OUchemical commented 8 years ago

Hi alongd,

I've copied files from you nitrogen branch, and still getting the following error. You can run my input file without problem? Can you send me the output file you generated from your end?

Thank you so much for your help!

Error: <Molecule "C123C(C4C(C=C14)C=C3)C=C2Nc1ccccc1"> Error: <Molecule "c12c(Nc3ccccc3)cccc1cccc2"> SKIP raise KineticsError for degeneracy flag Generating thermodynamics for new species... 1 radicals on [C]1=CC=CC2=CC=CC=C12 exceeds limit of 0. Using HBI method. Reading existing thermo file QMfiles\pm3\UFWIBTONFRDIAS-UHFFFAOYSA.thermo Traceback (most recent call last): File "C:\Anaconda2\RMGPycode\RMG-Py\rmg.py", line 165, in rmg.execute(inputFile, output_dir, **kwargs) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\main.py", line 539, in execute bimolecularReact=self.bimolecularReact) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 753, in enlarge spec.generateThermoData(database, quantumMechanics=self.quantumMechanics) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 99, in generateThermoData thermo0 = database.thermo.getThermoData(self, trainingSet=None, quantumMechanics=quantumMechanics) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\thermo.py", line 754, in getThermoData tdata = self.estimateRadicalThermoViaHBI(molecule, quantumMechanics.getThermoData) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\data\thermo.py", line 1055, in estimateRadicalThermoViaHBI thermoData_sat = stableThermoEstimator(saturatedStruct) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\qm\main.py", line 194, in getThermoData thermo0 = qm_molecule_calculator.generateThermoData() File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\qm\molecule.py", line 324, in generateThermoData if self.loadThermoData(): File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\qm\molecule.py", line 374, in loadThermoData loadedMolecule = rmgpy.molecule.Molecule().fromAdjacencyList(local_context['adjacencyList']) File "rmgpy\molecule\molecule.py", line 1183, in rmgpy.molecule.molecule.Molecule.fromAdjacencyList (build\Release\pyrex\rmgpy\molecule\mo lecule.c:22764) File "rmgpy\molecule\molecule.py", line 1191, in rmgpy.molecule.molecule.Molecule.fromAdjacencyList (build\Release\pyrex\rmgpy\molecule\mo lecule.c:22339) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\molecule\adjlist.py", line 710, in fromAdjacencyList for atom in atoms: ConsistencyChecker.check_partial_charge(atom) File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\molecule\adjlist.py", line 113, in check_partial_charge .format(symbol=atom.symbol, radicals=atom.radicalElectrons, lonePairs=atom.lonePairs, charge=atom.charge)) rmgpy.molecule.adjlist.InvalidAdjacencyListError: Invalid valency for atom C with 0 unpaired electrons, 0 pairs of electrons, and 0 charge.

alongd commented 8 years ago

Hi, My run of your input file went till the part when RMG requires MOPAC. I could reproduce most of the errors you described, and solved them by running on my nitrogen branch. Anyway, I do not have the output for your model (haven't installed the MOPAC password yet).

I suggest to focus on solving your error:

Traceback (most recent call last):
File "C:\Anaconda2\RMGPycode\RMG-Py\rmg.py", line 165, in
rmg.execute(inputFile, output_dir, **kwargs)
File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\main.py", line 549, in execute
bimolecularReact=self.bimolecularReact)
File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 733, in enlarge
self.processNewReactions(reactions, self.core.species[i], None)
File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 830, in processNewReactions
rxn, isNew = self.makeNewReaction(rxn)
File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 568, in makeNewReaction
reactants = [self.makeNewSpecies(reactant)[0] for reactant in forward.reactants]
File "C:\Anaconda2\RMGPycode\RMG-Py\rmgpy\rmg\model.py", line 424, in makeNewSpecies
molecule.clearLabeledAtoms()
AttributeError: 'rmgpy.species.Species' object has no attribute 'clearLabeledAtoms'

@nyee , can you perhaps comment on that?

@OUchemical , you can also try working on my nitrogen branch (use GIT to pull it), but I can't say whether or not more errors will pop later in the run.

OUchemical commented 8 years ago

Thanks @alongd , what I did is that I commented out adjlist.py line 111-113, where raise InvalidAdjacencyListError, and my program run for more than 24 hours now, and still working on thermochemistry calculation of different potential products by Mopac. Lots of complex structures have been generated. I'm not sure if this will eventually generate reactions model at all, do you have some comments on that? Many thanks!

alongd commented 8 years ago

Good luck, let us know how it goes!

OUchemical commented 8 years ago

My program run for 4.5 days, and stopped with the following error. I have also tried liquid phase phenylamine oxidation as well, same error appeared and stopped the calculation. Do you have comments on that?

Thanks,

File "C:\code\RMG-Py\rmgpy\data\thermo.py", line 781, in getThermoData shortDesc = thermo0.comment File "C:\code\RMG-Py\rmgpy\data\thermo.py", line 321, in loadEntry raise DatabaseError('Found a duplicate molecule with label {0} in the thermo library. Please correct your library.'.format(label)) rmgpy.data.base.DatabaseError: Found a duplicate molecule with label C12(C3C(C(c 4c1cccc4)C=[C]2)c1c(C(=C3)Nc2ccccc2)cccc1)Nc1ccccc1_(D) in the thermo library. Please correct your library.

rwest commented 8 years ago

Thanks for reporting back. Pllease could you provide the full stack trace? My guess is that this is to do with QMThermo.

Out of interest, how does the incomplete model perform? After 4.5 days it must be relatively large and hopefully has many important pathways. (You can get the incomplete model at any time from the "CHEMKIN" folder - you don't need to wait for RMG to finish)

OUchemical commented 8 years ago

Thank @rwest for replying, the log file is over 8 MB, most are QM Thermo calculation, since reactants are n-phenyl-1-naphthylamine (antioxidant), and oxygen in liquid reactor, RMG generated 203 edge species, and 210 edge reactions, 0 core species and 0 core reaction, therefore no kinetic data reported.

Traceback (most recent call last): File "C:\code\RMG-Py\rmg.py", line 165, in rmg.execute(inputFile, outputdir, **kwargs) File "C:\code\RMG-Py\rmgpy\rmg\main.py", line 640, in execute bimolecularReact=self.bimolecularReact) File "C:\code\RMG-Py\rmgpy\rmg\model.py", line 753, in enlarge spec.generateThermoData(database, quantumMechanics=self.quantumMechanics) File "C:\code\RMG-Py\rmgpy\rmg\model.py", line 99, in generateThermoData thermo0 = database.thermo.getThermoData(self, trainingSet=None, quantumMecha nics=quantumMechanics) File "C:\code\RMG-Py\rmgpy\data\thermo.py", line 781, in getThermoData shortDesc = thermo0.comment File "C:\code\RMG-Py\rmgpy\data\thermo.py", line 321, in loadEntry raise DatabaseError('Found a duplicate molecule with label {0} in the thermo library. Please correct your library.'.format(label)) rmgpy.data.base.DatabaseError: Found a duplicate molecule with label C12(C3C(C(c 4c1cccc4)C=[C]2)c1c(C(=C3)Nc2ccccc2)cccc1)Nc1ccccc1(D) in the thermo library. Please correct your library.

OUchemical commented 8 years ago

My program stuck at this point overnight, it calculated 1 new core species, and tried to add 1 new core reaction, but at the end, shows on the screen: "added 0 new core reactions"

Has anyone come across similar problem before?

Thanks,

Summary of Model Enlargement

Added 1 new core species C1(C2[C]3C(=CC=C2)C(=CC=C3)Nc2ccccc2)[C]2C(=CC=CC2=CC=C1)Nc1ccccc1(170) Created 0 new edge species Moved 1 reactions from edge to core C1(C2[C]3C(=CC=C2)C(=CC=C3)Nc2ccccc2)[C]2C(=CC=CC2=CC=C1)Nc1ccccc1(170) <=> PAN(3) + PAN(3) Added 0 new core reactions Created 0 new edge reactions

After model enlargement: The model core has 8 species and 1 reactions The model edge has 272 species and 394 reactions

alongd commented 8 years ago

did it give any error Traceback?

OUchemical commented 8 years ago

There's no error reported. Program is still running, but seems stuck, there's no more print over past 10 hours.

alongd commented 8 years ago

OK, could be a memory issue. Look for the latest appearance of the memory report, should look like:

Chemkin file contains 2201 reactions.
Updating RMG execution statistics...
    Execution time (DD:HH:MM:SS): 00:02:59:46
    Memory used: 2903.03 MB

Plus, how much memory do you have on the computer running the simulation?

OUchemical commented 8 years ago

Thanks for you reply, my computer has 32 GB of RAM, and is only used less than 9 GB. So probably not memory limiting issue.

nickvandewiele commented 8 years ago

I noticed that the simulation involves a molecule that is, to RMG's standards, quite large and complex. I believe that although it seems that RMG has crashed, it is still simulating.

RMG is in the process of generating all possible reactions between the molecules in the core, against the given list of reaction families. That process may take a very long time, especially with molecules with many reacting sites.

In addition, the presence of aromatic rings will lead to structures with many resonance isomers, which also adds up.

I don't know the specifics of your RMG simulation, but this simulation may take a very long time (weeks), and memory (10s of GB) to complete. Be aware that many options in the input file may drastically increase the computational requirements of a simulation. Some guidelines are provided here.

alongd commented 8 years ago

@OUchemical, let us know how it's going. Perhaps we should close this issue?

OUchemical commented 8 years ago

Sorry for the late reply, and yes we could close the issue I think. I've done through various MOPAC calculation in terms of obtaining the activation energy for pre-picked reactions, reaction rate will be verified by experiments, and I'm planning on building a specific PAN oxidation library for RMG, I'll let you know with updates

On Tue, Aug 9, 2016 at 4:13 PM, Alon Grinberg Dana <notifications@github.com

wrote:

@OUchemical https://github.com/OUchemical, let us know how it's going. Perhaps we should close this issue?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ReactionMechanismGenerator/RMG-Py/issues/645#issuecomment-238676211, or mute the thread https://github.com/notifications/unsubscribe-auth/ASZX5x-JTBYeeRsS6bDfEtXQk_3uBI4rks5qeN9mgaJpZM4IbSK2 .

Yifan Yang Postdoctoral researcher

Office: Canal building 3105 Department of Mechanical and Aerospace Engineering Carleton University, 1125 Colonel By Drive, Ottawa, Ontario, Canada K1S 5B6 E-mail: yifanyang@cunet.carleton.ca yifanyang@cunet.carleton.ca;

alongd commented 8 years ago

Great, please do update us, and post issues you encounter.