From 05c9e44006a838c72e1adde71c4fac4a2fd0e423 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 16 Jun 2011 15:59:18 +0100 Subject: [PATCH 1/6] Fix RMG PDep bug writing measure input file. The commit e1f226c3d55b3233039f3911ff468c16573496fd changed the method for writing a measure input file, but not the call to it from rmg/model.py --- rmgpy/rmg/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index c328c07b4d..5977e2540d 100755 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -488,8 +488,8 @@ def update(self, reactionModel, database): self.collisionModel = SingleExponentialDownModel(alpha0=4.86 * 4184) # Save input file - rmgpy.measure.output.writeInput(os.path.join(settings.outputDirectory, 'pdep', 'network{0:d}_{1:d}.py'.format(self.index, len(self.isomers))), - self, Tlist, Plist, (grainSize, numGrains), method, model) + rmgpy.measure.output.writeFile(os.path.join(settings.outputDirectory, 'pdep', 'network{0:d}_{1:d}.py'.format(self.index, len(self.isomers))), + self, Tlist, Plist, (grainSize, numGrains), method, model, Tmin, Tmax, Pmin, Pmax) self.printSummary(level=logging.INFO) From ca3b3bc790442430e73fb552d07ee7b243ecd7e6 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 16 Jun 2011 20:45:02 +0100 Subject: [PATCH 2/6] Added (failing) radical detection test to unittest/moleculeTest.py When reading in from the smiles '[CH]' the carbon ends up with a radicalElectrons count of 0 instead of 3. Strangely, [CH2] and [CH3] work fine. This causes problems in data/thermo.py when trying to find the thermo for [CH]. --- rmgpy/data/thermo.py | 4 ++-- unittest/moleculeTest.py | 17 ++++++++++++++++- 2 files changed, 18 insertions(+), 3 deletions(-) diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 4fb9647524..bf13f6fb21 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -582,7 +582,7 @@ def getThermoDataFromGroups(self, molecule): thermoData = None - if sum([atom.radicalElectrons for atom in molecule.atoms]) > 0: + 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) @@ -646,7 +646,7 @@ def getThermoDataFromGroups(self, molecule): # Correct the entropy for the symmetry number - else: + else: # non-radical species # Generate estimate of thermodynamics for atom in molecule.atoms: # Iterate over heavy (non-hydrogen) atoms diff --git a/unittest/moleculeTest.py b/unittest/moleculeTest.py index 8ab51b33c4..936be75848 100644 --- a/unittest/moleculeTest.py +++ b/unittest/moleculeTest.py @@ -856,7 +856,22 @@ def testPickle(self): self.assertEqual(molecule0.getFormula(), molecule.getFormula()) self.assertTrue(molecule0.isIsomorphic(molecule)) self.assertTrue(molecule.isIsomorphic(molecule0)) - + + def testRadicalCH(self): + """ + Test that the species [CH] has three radical electrons and a spin multiplicity of 4. + """ + molecule = Molecule().fromSMILES('[CH]') + self.assertEqual(molecule.atoms[0].radicalElectrons, 3) + self.assertEqual(molecule.atoms[0].spinMultiplicity, 4) + + def testRadicalCH2(self): + """ + Test that the species [CH2] has two radical electrons and a spin multiplicity of 3. + """ + molecule = Molecule().fromSMILES('[CH2]') + self.assertEqual(molecule.atoms[0].radicalElectrons, 2) + self.assertEqual(molecule.atoms[0].spinMultiplicity, 3) ################################################################################ class TestMoleculeSymmetry(unittest.TestCase): From 4a12a62fea572d21706c5bc66edf199caa286ce9 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 16 Jun 2011 21:36:55 +0100 Subject: [PATCH 3/6] Wipe the output.html file from previous run by writing a blank one sooner. As soon as the input file has been read in, we write an output file which has 0 species and 0 reactions. The species pictures etc. from the previous run have been wiped by this point already, but it can be confusing to have the previous output file still there when you're running a new job. --- rmg.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/rmg.py b/rmg.py index fd8429a75b..59a2de34f0 100755 --- a/rmg.py +++ b/rmg.py @@ -251,6 +251,10 @@ def execute(args): # Read input file reactionModel, coreSpecies, reactionSystems, database, seedMechanisms = readInputFile(inputFile) + # Delete previous HTML file + from rmgpy.rmg.output import saveOutputHTML + saveOutputHTML(os.path.join(settings.outputDirectory, 'output.html'), reactionModel) + # Initialize reaction model if args.restart: reactionModel = loadRestartFile(os.path.join(settings.outputDirectory,'restart.pkl.gz')) @@ -335,7 +339,6 @@ def execute(args): # Save the current state of the model core to a pretty HTML file logging.info('Saving latest model core to HTML file...') - from rmgpy.rmg.output import saveOutputHTML saveOutputHTML(os.path.join(settings.outputDirectory, 'output.html'), reactionModel) # Save a Chemkin file containing the current model core From f0d626ff79f37d2a1d77ac9a51e255e87e32db8a Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 16 Jun 2011 21:38:48 +0100 Subject: [PATCH 4/6] Copy the latest chemkin file to chem.inp Whenever the latest chemXXXX.inp file is written, it is copied (or rather hard-linked) to chemkin/chem.inp, so that you can always just look at this file and get the latest. --- rmg.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/rmg.py b/rmg.py index 59a2de34f0..86701c768e 100755 --- a/rmg.py +++ b/rmg.py @@ -343,7 +343,12 @@ def execute(args): # Save a Chemkin file containing the current model core logging.info('Saving latest model core to Chemkin file...') - reactionModel.saveChemkinFile(os.path.join(settings.outputDirectory, 'chemkin', 'chem%04i.inp' % len(reactionModel.core.species))) + this_chemkin_path = os.path.join(settings.outputDirectory, 'chemkin', 'chem%04i.inp' % len(reactionModel.core.species)) + latest_chemkin_path = os.path.join(settings.outputDirectory, 'chemkin','chem.inp') + reactionModel.saveChemkinFile(this_chemkin_path) + if os.path.exists(latest_chemkin_path): + os.unlink(latest_chemkin_path) + os.link(this_chemkin_path,latest_chemkin_path) # Save the restart file if desired if settings.saveRestart: From 9497be35b3e1f8a0713a23cfff71319cbcad2e59 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 16 Jun 2011 22:09:41 +0100 Subject: [PATCH 5/6] RMG-py creates better measure input files. It used to write isomer("HCCOH") reactants("C2H", "OH") instead of isomer("HCCOH(31)") reactants("C2H(3)", "OH(32)") This closes #40 --- rmgpy/measure/output.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/measure/output.py b/rmgpy/measure/output.py index 166e559d62..1b51c5f279 100644 --- a/rmgpy/measure/output.py +++ b/rmgpy/measure/output.py @@ -118,10 +118,10 @@ def writeNetworkConfigurations(f, network): """ # Write isomer configurations for isomer in network.isomers: - f.write('isomer("{0}")\n\n'.format(isomer.label)) + f.write('isomer("{0}")\n\n'.format(isomer)) # Write reactant configurations for reactants in network.reactants: - f.write('reactants("{0}", "{1}")\n\n'.format(reactants[0].label, reactants[1].label)) + f.write('reactants("{0}", "{1}")\n\n'.format(reactants[0], reactants[1])) # No need to write product configurations, as these are assumed ################################################################################ From 060e3c0b33809641b38a75c88aa58acbb2ccbf32 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 16 Jun 2011 22:56:43 +0100 Subject: [PATCH 6/6] Added my methylformate test case to the examples/rmg folder. This is close to something I've been using in RMG-Java. --- examples/rmg/methylformate/input.py | 144 ++++++++++++++++++++++++++++ 1 file changed, 144 insertions(+) create mode 100644 examples/rmg/methylformate/input.py diff --git a/examples/rmg/methylformate/input.py b/examples/rmg/methylformate/input.py new file mode 100644 index 0000000000..5f4fb681ac --- /dev/null +++ b/examples/rmg/methylformate/input.py @@ -0,0 +1,144 @@ +# Data sources +database( + '../RMG-database/input', + thermoLibraries = ['primaryThermoLibrary','DFT_QCI_thermo','GRI-Mech3.0'], + reactionLibraries = ['Methylformate','Glarborg/highP'], + seedMechanisms = ['GRI-Mech3.0'], +) + +# List of species +species( + label='Mfmt', + reactive=True, + structure=SMILES("COC=O"), +) +species( + label='O2', + reactive=True, + structure=SMILES("[O][O]"), +) +species( + label='C2H', + reactive=True, + structure=SMILES("C#[C]"), +) +species( + label='C2H', + reactive=True, + structure=SMILES("C#[C]"), +) +species( + label='CH', + reactive=True, + structure=adjacencyList( + """ + 1 C 3 {2,S} + 2 H 0 {1,S} + """), +) +species( + label='H2O', + reactive=True, + structure=SMILES("O"), +) +species( + label='H2', + reactive=True, + structure=SMILES("[H][H]"), +) +species( + label='CO', + reactive=True, + structure=SMILES("[C]=O"), +) +species( + label='CO2', + reactive=True, + structure=SMILES("C(=O)=O"), +) +species( + label='CH4', + reactive=True, + structure=SMILES("C"), +) +species( + label='CH3', + reactive=True, + structure=SMILES("[CH3]"), +) +species( + label='CH3OH', + reactive=True, + structure=SMILES("CO"), +) +species( + label='C2H4', + reactive=True, + structure=SMILES("C=C"), +) +species( + label='C2H2', + reactive=True, + structure=SMILES("C#C"), +) +species( + label='CH2O', + reactive=True, + structure=SMILES("C=O"), +) +species( + label='CH3CHO', + reactive=True, + structure=SMILES("CC=O"), +) + + +# Bath gas +species( + label='Ar', + reactive=False, + structure=InChI("InChI=1S/Ar"), +) + +# Reaction systems +simpleReactor( + temperature=(1350,'K'), + pressure=(1.0,'bar'), + initialMoleFractions={ + "Mfmt": 0.01, + "O2": 0.02, + "Ar": 0.97, + }, +) + +termination( + time=(0.5,'s'), +) + +simulator( + atol=1e-22, + rtol=1e-8, +) + +model( + toleranceKeepInEdge=0.0, + toleranceMoveToCore=0.005, + toleranceInterruptSimulation=1.0, + maximumEdgeSpecies=100000 +) + +pressureDependence( + method='modified strong collision', # 'reservoir state' + minimumGrainSize=(2.0,'kcal/mol'), + minimumNumberOfGrains=200, + temperatures=(290,'K',3500,'K',8), + pressures=(0.02,'bar',100,'bar',5), + interpolation=('Chebyshev', 6, 4), +) + +options( + units='si', + saveRestart=False, + drawMolecules=False, + generatePlots=True, +)