search

ReactionMechanismGenerator / RMG-Py

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

Problem in Surface_Bidentate_Dissociation causing charged N adsorbate #2039

Closed Tingchenlee closed 3 years ago

Tingchenlee commented 3 years ago

Bug Description

This minimal input file below leads to this crash caused by the reaction family Surface_Bidentate_Dissociation:

For reaction generation 1 process is used.
Error: Could not update atomtypes for this molecule:
multiplicity -187
1 *4 X u0 p0 c0 {2,S}
2 *1 N u0 p1 c-1 {1,S} {3,T}
3 *2 N u0 p1 c-1 {2,T} {4,S}
4 *3 H u0 p0 c0 {3,S}

Error: Problem family: Surface_Bidentate_Dissociation
Error: Problem reactants: (Molecule(smiles="[Pt]"), Molecule(smiles="[Pt]N=N[Pt]"))
Error: <Molecule "[Pt]">
1 *3 H u0 p0 c0 {2,S}
2 *6 X u0 p0 c0 {1,S}

Error: <Molecule "[Pt]N=N[Pt]">
1 *1 N u0 p1 c0 {2,D} {3,S}
2 *2 N u0 p1 c0 {1,D} {4,S}
3 *4 X u0 p0 c0 {1,S}
4 *5 X u0 p0 c0 {2,S}

Traceback (most recent call last):
  File "/Users/lee.ting/Code/RMG-Py/rmg.py", line 118, in <module>
    main()
  File "/Users/lee.ting/Code/RMG-Py/rmg.py", line 112, in main
    rmg.execute(**kwargs)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/main.py", line 922, in execute
    trimolecular_react=self.trimolecular_react)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/model.py", line 599, in enlarge
    procnum=procnum)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/react.py", line 172, in react_all
    return react(spc_fam_tuples, procnum), [fam_tuple[0] for fam_tuple in spc_fam_tuples]
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/react.py", line 66, in react
    reactions = list(map(_react_species_star, spc_fam_tuples))
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/react.py", line 79, in _react_species_star
    return react_species(*args)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/react.py", line 94, in react_species
    reactions = get_db('kinetics').generate_reactions_from_families(species_tuple, only_families=only_families)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/kinetics/database.py", line 534, in generate_reactions_from_families
    prod_resonance=resonance))
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/kinetics/database.py", line 561, in react_molecules
    prod_resonance=prod_resonance))
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1715, in generate_reactions
    self._generate_reactions(reactants, products=products, forward=False, prod_resonance=prod_resonance))
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 2014, in _generate_reactions
    [map_a, map_b], forward)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1588, in _generate_product_structures
    product_structures = self.apply_recipe(reactant_structures, forward=forward)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1517, in apply_recipe
    struct.update(sort_atoms=not self.save_order)
  File "rmgpy/molecule/molecule.py", line 1130, in rmgpy.molecule.molecule.Molecule.update
  File "rmgpy/molecule/molecule.py", line 1355, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/molecule.py", line 1350, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/atomtype.py", line 729, in rmgpy.molecule.atomtype.get_atomtype
  File "rmgpy/molecule/atomtype.py", line 757, in rmgpy.molecule.atomtype.get_atomtype
rmgpy.exceptions.AtomTypeError: Unable to determine atom type for atom N-, which has 1 single bonds, 0 double bonds (0 to O, 0 to S, 0 others), 1 triple bonds, 0 quadruple bonds, 0 benzene bonds, 1 lone pairs, and -1 charge.

How To Reproduce

Try this input file:

database(
    thermoLibraries=['surfaceThermoPt', 'primaryThermoLibrary'],
    reactionLibraries = ['Surface/CPOX_Pt/Deutschmann2006'],
    seedMechanisms = [],
    kineticsDepositories = ['training'],
    kineticsFamilies = ['!Surface_Adsorption_Abstraction_vdW'],
    kineticsEstimator = 'rate rules',

)

catalystProperties(
    bindingEnergies= {  
                     'H':(-2.75367887E+00, 'eV/molecule'),
                         'C':(-7.02515507E+00, 'eV/molecule'),
                         'N':(-4.63224568E+00, 'eV/molecule'),
                         'O':(-3.81153179E+00, 'eV/molecule'),
                      },
    surfaceSiteDensity=(2.483E-09, 'mol/cm^2'),
)

species(
    label='CH4',
    reactive=True,
    structure=SMILES("[CH4]"),
)

species(
    label='H2',
    reactive=True,
    structure=SMILES("[H][H]"),
)

species(
    label='H2O',
    reactive=True,
    structure=SMILES("O"),
)

species(
    label='N2',
    reactive=True,
    structure=SMILES("N#N"),
)

species(
    label='X',
    reactive=True,
    structure=adjacencyList("1 X u0"),
)

species(
    label='O2',
    reactive=True,
    structure=SMILES("[O][O]"),
)

surfaceReactor(  
    temperature=(723,'K'),
    initialPressure=(50.0, 'bar'),
    initialGasMoleFractions={
        "CH4": 0.0,
        "H2": 0.35,
        "H2O": 0.20,
        "N2": 0.25,
        "O2": 0.20,
    },
    initialSurfaceCoverages={
        "X": 1.0,
    },
    surfaceVolumeRatio=(1.e5, 'm^-1'),
    terminationConversion = { "N2":0.80,},
    terminationTime=(0.1, 's'),
)

simulator(
    atol=1e-18,
    rtol=1e-12,
)

model( 
    toleranceKeepInEdge=0.0,
    toleranceMoveToCore=0.00000001,
    toleranceInterruptSimulation=1,
    maximumEdgeSpecies=100000,
    minCoreSizeForPrune=150,
    toleranceThermoKeepSpeciesInEdge=0.5,
    minSpeciesExistIterationsForPrune=4,
)

options(
    units='si',
    saveRestartPeriod=None,
    generateOutputHTML=True,
    generatePlots=True, 
    saveEdgeSpecies=True,
    saveSimulationProfiles=True,
)

Expected Behavior

It wouldn't crash.

Installation Information

Tingchenlee commented 3 years ago

This issue is probably addressed by: ReactionMechanismGenerator/RMG-database#445

mazeau commented 3 years ago

Made a PR here to the RMG database to correct this issue

rwest commented 3 years ago

Made a PR here to the RMG database to correct this issue

I think that PR addresses #2038 not this one.

Tingchenlee commented 3 years ago

After trying https://github.com/ReactionMechanismGenerator/RMG-database/pull/445, I received the error messages:

Error: Could not update atomtypes for this molecule:
multiplicity 1
1 O u0 p2 c0 {2,S} {4,S}
2 N u1 p0 c0 {1,S} {3,T}
3 N u2 p0 c0 {2,T}
4 H u0 p0 c0 {1,S}
Error: Error attempting to get thermo for species ON(N=[Pt])[Pt](43) with structure 
1 O u0 p2 c0 {2,S} {4,S}
2 N u0 p1 c0 {1,S} {3,S} {5,S}
3 N u0 p1 c0 {2,S} {6,D}
4 H u0 p0 c0 {1,S}
5 X u0 p0 c0 {2,S}
6 X u0 p0 c0 {3,D}
Traceback (most recent call last):
  File "/Users/lee.ting/Code/RMG-Py/rmg.py", line 118, in <module>
    main()
  File "/Users/lee.ting/Code/RMG-Py/rmg.py", line 112, in main
    rmg.execute(**kwargs)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/main.py", line 922, in execute
    trimolecular_react=self.trimolecular_react)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/model.py", line 614, in enlarge
    self.apply_thermo_to_species(procnum)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/model.py", line 824, in apply_thermo_to_species
    self.generate_thermo(spc, rename=True)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/model.py", line 831, in generate_thermo
    submit(spc, self.solvent_name)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/thermo/thermoengine.py", line 175, in submit
    spc.thermo = evaluator(spc, solvent_name=solvent_name)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/thermo/thermoengine.py", line 160, in evaluator
    thermo = generate_thermo_data(spc, solvent_name=solvent_name)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/thermo/thermoengine.py", line 125, in generate_thermo_data
    thermo0 = thermodb.get_thermo_data(spc)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/thermo.py", line 1223, in get_thermo_data
    thermo0 = self.get_thermo_data_for_surface_species(species)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/thermo.py", line 1546, in get_thermo_data_for_surface_species
    dummy_molecule.update()
  File "rmgpy/molecule/molecule.py", line 1130, in rmgpy.molecule.molecule.Molecule.update
  File "rmgpy/molecule/molecule.py", line 1355, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/molecule.py", line 1350, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/atomtype.py", line 729, in rmgpy.molecule.atomtype.get_atomtype
  File "rmgpy/molecule/atomtype.py", line 757, in rmgpy.molecule.atomtype.get_atomtype
rmgpy.exceptions.AtomTypeError: Unable to determine atom type for atom N., which has 1 single bonds, 0 double bonds (0 to O, 0 to S, 0 others), 1 triple bonds, 0 quadruple bonds, 0 benzene bonds, 0 lone pairs, and +0 charge.
rwest commented 3 years ago

I think it's trying to get the thermo for the bidentate

Error: Error attempting to get thermo for species ON(N=[Pt])[Pt](43) with structure 
1 O u0 p2 c0 {2,S} {4,S}
2 N u0 p1 c0 {1,S} {3,S} {5,S}
3 N u0 p1 c0 {2,S} {6,D}
4 H u0 p0 c0 {1,S}
5 X u0 p0 c0 {2,S}
6 X u0 p0 c0 {3,D}

but to do so creates the gas phase counterpart, and it comes up with 

multiplicity 1
1 O u0 p2 c0 {2,S} {4,S}
2 N u1 p0 c0 {1,S} {3,T}
3 N u2 p0 c0 {2,T}
4 H u0 p0 c0 {1,S}

which as @mazeau pointed out, has two unpaired electrons on atom 3. Maybe the thermo estimation of bidentate adsorption things is the problem? (But the title of this issue is Surface_Bidentate_Dissociation )

Tingchenlee commented 3 years ago

Should we discuss this in other issues? I encountered another problem species when I was testing my new library(metal=Ru):

Error: Could not update atomtypes for this molecule:
multiplicity 1
1 N u1 p0 c0 {2,S} {3,T}
2 N u0 p1 c0 {1,S} {4,D}
3 N u2 p0 c0 {1,T}
4 N u0 p1 c0 {2,D} {5,S}
5 H u0 p0 c0 {4,S}

Error: Error attempting to get thermo for species N=NN(N=[Pt])[Pt](43) with structure 
1 N u0 p1 c0 {2,S} {3,S} {6,S}
2 N u0 p1 c0 {1,S} {4,D}
3 N u0 p1 c0 {1,S} {7,D}
4 N u0 p1 c0 {2,D} {5,S}
5 H u0 p0 c0 {4,S}
6 X u0 p0 c0 {1,S}
7 X u0 p0 c0 {3,D}

Traceback (most recent call last):
  File "/Users/lee.ting/Code/RMG-Py/rmg.py", line 118, in <module>
    main()
  File "/Users/lee.ting/Code/RMG-Py/rmg.py", line 112, in main
    rmg.execute(**kwargs)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/main.py", line 922, in execute
    trimolecular_react=self.trimolecular_react)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/model.py", line 614, in enlarge
    self.apply_thermo_to_species(procnum)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/model.py", line 824, in apply_thermo_to_species
    self.generate_thermo(spc, rename=True)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/rmg/model.py", line 831, in generate_thermo
    submit(spc, self.solvent_name)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/thermo/thermoengine.py", line 175, in submit
    spc.thermo = evaluator(spc, solvent_name=solvent_name)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/thermo/thermoengine.py", line 160, in evaluator
    thermo = generate_thermo_data(spc, solvent_name=solvent_name)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/thermo/thermoengine.py", line 125, in generate_thermo_data
    thermo0 = thermodb.get_thermo_data(spc)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/thermo.py", line 1223, in get_thermo_data
    thermo0 = self.get_thermo_data_for_surface_species(species)
  File "/Users/lee.ting/Code/RMG-Py/rmgpy/data/thermo.py", line 1546, in get_thermo_data_for_surface_species
    dummy_molecule.update()
  File "rmgpy/molecule/molecule.py", line 1130, in rmgpy.molecule.molecule.Molecule.update
  File "rmgpy/molecule/molecule.py", line 1355, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/molecule.py", line 1350, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/atomtype.py", line 729, in rmgpy.molecule.atomtype.get_atomtype
  File "rmgpy/molecule/atomtype.py", line 757, in rmgpy.molecule.atomtype.get_atomtype
rmgpy.exceptions.AtomTypeError: Unable to determine atom type for atom N., which has 1 single bonds, 0 double bonds (0 to O, 0 to S, 0 others), 1 triple bonds, 0 quadruple bonds, 0 benzene bonds, 0 lone pairs, and +0 charge.
Tingchenlee commented 3 years ago

Katrin referred that it looks like RMG isn't updating the bonds and lone pairs correctly when removing it from the surface, like the #2034 Other testings of Rh and Pd surface have the same problem of species:

Error: Error attempting to get thermo for species ON(N=[Pt])[Pt](48) with structure 
1 O u0 p2 c0 {2,S} {4,S}
2 N u0 p1 c0 {1,S} {3,S} {5,S}
3 N u0 p1 c0 {2,S} {6,D}
4 H u0 p0 c0 {1,S}
5 X u0 p0 c0 {2,S}
6 X u0 p0 c0 {3,D}
josefedupe commented 3 years ago

Hello there, I was a facing a similar scenario with the thermo of the bidentate adsorbate for adsorption of [O][O] over Ni:

For reaction generation 1 process is used.
Generating thermo for new species...
Error: Could not update atomtypes for this molecule:
1 O u0 p1 c+1 {2,T}
2 O u0 p2 c-1 {1,T}

Error: Error attempting to get thermo for species [Pt]OO[Pt](18) with structure
1 O u0 p2 c0 {2,S} {3,S}
2 O u0 p2 c0 {1,S} {4,S}
3 X u0 p0 c0 {1,S}
4 X u0 p0 c0 {2,S}

Traceback (most recent call last):
  File "/mnt/c/users/josefedupe/RMG-Py/rmg.py", line 118, in <module>
    main()
  File "/mnt/c/users/josefedupe/RMG-Py/rmg.py", line 112, in main
    rmg.execute(**kwargs)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/rmg/main.py", line 922, in execute
    trimolecular_react=self.trimolecular_react)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/rmg/model.py", line 610, in enlarge
    self.apply_thermo_to_species(procnum)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/rmg/model.py", line 820, in apply_thermo_to_species
    self.generate_thermo(spc, rename=True)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/rmg/model.py", line 827, in generate_thermo
    submit(spc, self.solvent_name)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/thermo/thermoengine.py", line 175, in submit
    spc.thermo = evaluator(spc, solvent_name=solvent_name)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/thermo/thermoengine.py", line 160, in evaluator
    thermo = generate_thermo_data(spc, solvent_name=solvent_name)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/thermo/thermoengine.py", line 125, in generate_thermo_data
    thermo0 = thermodb.get_thermo_data(spc)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/data/thermo.py", line 1229, in get_thermo_data
    thermo0 = self.get_thermo_data_for_surface_species(species)
  File "/mnt/c/users/josefedupe/RMG-Py/rmgpy/data/thermo.py", line 1562, in get_thermo_data_for_surface_species
    dummy_molecule.update()
  File "rmgpy/molecule/molecule.py", line 1130, in rmgpy.molecule.molecule.Molecule.update
  File "rmgpy/molecule/molecule.py", line 1355, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/molecule.py", line 1350, in rmgpy.molecule.molecule.Molecule.update_atomtypes
  File "rmgpy/molecule/atomtype.py", line 729, in rmgpy.molecule.atomtype.get_atomtype
  File "rmgpy/molecule/atomtype.py", line 757, in rmgpy.molecule.atomtype.get_atomtype
rmgpy.exceptions.AtomTypeError: Unable to determine atom type for atom O-, which has 0 single bonds, 0 double bonds (0 to O, 0 to S, 0 others), 1 triple bonds, 0 quadruple bonds, 0 benzene bonds, 2 lone pairs, and -1 charge.

I think, it was result of the treatment of thermo.py for bidentate adsorption where the bidentate adsorption of CO, X-C=O-X, requires to generate C[-1]#O[+1]. With this in mind, I add the following “try: - except:” to treat [Pt]OO[Pt] and generate [O][O]:

if len(adsorbed_atoms) == 2:
            # Bidentate adsorption.
            # Conserve the old molecular structure after sites removal in case
            # of generate an invalid molecule, i.e.  O[-1]#O[+1]
            dummy_molecule_old = dummy_molecule.copy(deep=True)
            # Try to turn adjacent biradical into a bond.
            try:
                bond = adsorbed_atoms[0].bonds[adsorbed_atoms[1]]
            except KeyError:
                pass # the two adsorbed atoms are not bonded to each other
            else:
                if bond.order < 3:
                    bond.increment_order()
                    adsorbed_atoms[0].decrement_radical()
                    adsorbed_atoms[1].decrement_radical()
                    if (adsorbed_atoms[0].radical_electrons and
                            adsorbed_atoms[1].radical_electrons and
                            bond.order < 3):
                        # There are still spare adjacenct radicals, so do it again
                        bond.increment_order()
                        adsorbed_atoms[0].decrement_radical()
                        adsorbed_atoms[1].decrement_radical()
                    if (adsorbed_atoms[0].lone_pairs and
                            adsorbed_atoms[1].lone_pairs and 
                            bond.order < 3):
                        # X#C-C#X will end up with .:C-C:. in gas phase
                        # and we want to get to .C#C. but not :C=C:
                        bond.increment_order()
                        adsorbed_atoms[0].decrement_lone_pairs()
                        adsorbed_atoms[0].increment_radical()
                        adsorbed_atoms[1].decrement_lone_pairs()
                        adsorbed_atoms[1].increment_radical()
                #For bidentate CO because we want C[-1]#O[+1] but not .C#O.
                #and for bidentate OO we want [O][O] but not O[-1]#O[+1], see
                #line 1556
                if (bond.order == 3 and adsorbed_atoms[0].radical_electrons and 
                    adsorbed_atoms[1].radical_electrons and 
                    (adsorbed_atoms[0].lone_pairs or adsorbed_atoms[1].lone_pairs)):
                    adsorbed_atoms[0].decrement_radical()
                    adsorbed_atoms[1].decrement_radical()
                    if adsorbed_atoms[0].lone_pairs:
                        adsorbed_atoms[1].increment_lone_pairs()
                    else:
                        adsorbed_atoms[0].increment_lone_pairs()
            try:
                dummy_molecule.update()
            except rmgpy.exceptions.AtomTypeError:
                dummy_molecule = dummy_molecule_old
                logging.exception("The above molecule was reestructured to: \n{0}".format(dummy_molecule.to_adjacency_list()))
                pass

        dummy_molecule.update_connectivity_values()
        dummy_molecule.update()
rwest commented 3 years ago

Hmm. This is getting complicated. Seems we want C*O* to become C[-1]#O[+1] O*O* to become [O][O] C*C* to become [C]#[C] N*N*O to become N#N=O (the middle Nitrogen has valence 5) etc...

Do we treat CO and OO as special known cases, and then everything else, a more generic arroach like: generate everything feasible, try them all, and use the lowest energy that doesn't fail?

rwest commented 3 years ago

Starting with 

1 O u0 p2 c0 {2,S} {4,S}
2 N u0 p1 c0 {1,S} {3,S} {5,S}
3 N u0 p1 c0 {2,S} {6,D}
4 H u0 p0 c0 {1,S}
5 X u0 p0 c0 {2,S}
6 X u0 p0 c0 {3,D}

it first creates the gas phase counterpart

multiplicity 1
1 O u0 p2 c0 {2,S} {4,S}
2 N u1 p1 c0 {1,S} {3,S}
3 N u2 p1 c0 {2,S}
4 H u0 p0 c0 {1,S}

note the three unpaired electrons total then it turns two of them into a double bond

                if bond.order < 3:
                    bond.increment_order()
                    adsorbed_atoms[0].decrement_radical()
                    adsorbed_atoms[1].decrement_radical()

to get 

multiplicity 1
1 O u0 p2 c0 {2,S} {4,S}
2 N u0 p1 c0 {1,S} {3,D}
3 N u1 p1 c0 {2,D}
4 H u0 p0 c0 {1,S}

then it tries to turn the lone pair into another bond

                    if (adsorbed_atoms[0].lone_pairs and
                            adsorbed_atoms[1].lone_pairs and 
                            bond.order < 3):
                        # X#C-C#X will end up with .:C-C:. in gas phase
                        # and we want to get to .C#C. but not :C=C:
                        bond.increment_order()
                        adsorbed_atoms[0].decrement_lone_pairs()
                        adsorbed_atoms[0].increment_radical()
                        adsorbed_atoms[1].decrement_lone_pairs()
                        adsorbed_atoms[1].increment_radical()

and gets

multiplicity 1
1 O u0 p2 c0 {2,S} {4,S}
2 N u1 p0 c+1 {1,S} {3,T}
3 N u2 p0 c+1 {2,T}
4 H u0 p0 c0 {1,S}
josefedupe commented 3 years ago

Yes, I think the problem would be to decide which molecular structure for the gas phase counterpart should be used. A practical solution could be to test after each “if” statement the molecular structure, if it does not have “AtomTypeError”. But if the last is applied the final gas phase counter part of [Pt]OO[Pt] will be O=O, which, by default is set to not be allowed in the field “generatedSpeciesConstrains” of the “input.py” file. Then, it could be tested the H0 of each molecular structure, and select the one with lower H0 and without “AtomTypeErrors”. With that in mind, may be the lines to be added could be as follows,


if len(adsorbed_atoms) == 2:
            # Bidentate adsorption.
            # Conserve the old molecular structure after sites removal in case
            # of generate an invalid molecule, i.e.  O[-1]#O[+1]
            #print('here')
            dummy_molecule_old = dummy_molecule.copy(deep=True)
            dummy_molecule_old.update_connectivity_values()
            dummy_molecule_old.update()
            dummy_species_old = Species()
            dummy_species_old.molecule.append(dummy_molecule_old)
            dummy_species_old.generate_resonance_structures()
            H298_old = []
            H298_old.append(self.get_thermo_data(dummy_species_old).get_enthalpy(298.))
            molecule_test = []
            molecule_test.append(dummy_molecule_old.to_smiles())
            # Try to turn adjacent biradical into a bond.
            try:
                bond = adsorbed_atoms[0].bonds[adsorbed_atoms[1]]
            except KeyError:
                pass # the two adsorbed atoms are not bonded to each other
            else:
                if bond.order < 3:
                    #print('in here 1533')
                    bond.increment_order()
                    adsorbed_atoms[0].decrement_radical()
                    adsorbed_atoms[1].decrement_radical()
                    try:
                        dummy_molecule.update_connectivity_values()
                        dummy_molecule.update()
                        dummy_species_old = Species()
                        dummy_species_old.molecule.append(dummy_molecule)
                        dummy_species_old.generate_resonance_structures()
                        molecule_test.append(dummy_molecule.to_smiles())
                        H298_old.append(self.get_thermo_data(dummy_species_old).get_enthalpy(298.))
                        # print(dummy_species_old)
                    except rmgpy.exceptions.AtomTypeError:
                        #dummy_molecule = dummy_molecule_old
                        #logging.exception("The above molecule was reestructured to: \n{0}".format(dummy_molecule_old.to_adjacency_list()))
                        pass
                    if (adsorbed_atoms[0].radical_electrons and
                            adsorbed_atoms[1].radical_electrons and
                            bond.order < 3):
                        # There are still spare adjacenct radicals, so do it again
                        print('in here 1551')
                        bond.increment_order()
                        adsorbed_atoms[0].decrement_radical()
                        adsorbed_atoms[1].decrement_radical()
                        try:
                            dummy_molecule.update_connectivity_values()
                            dummy_molecule.update()
                            dummy_species_old = Species()
                            dummy_species_old.molecule.append(dummy_molecule)
                            dummy_species_old.generate_resonance_structures()
                            molecule_test.append(dummy_molecule.to_smiles())
                            H298_old.append(self.get_thermo_data(dummy_species_old).get_enthalpy(298.))
                        except rmgpy.exceptions.AtomTypeError:
                            #print('exception 1563')
                            pass
                    if (adsorbed_atoms[0].lone_pairs and
                            adsorbed_atoms[1].lone_pairs and 
                            bond.order < 3):
                        # X#C-C#X will end up with .:C-C:. in gas phase
                        # and we want to get to .C#C. but not :C=C:
                        #print('in here 1571')
                        bond.increment_order()
                        adsorbed_atoms[0].decrement_lone_pairs()
                        adsorbed_atoms[0].increment_radical()
                        adsorbed_atoms[1].decrement_lone_pairs()
                        adsorbed_atoms[1].increment_radical()
                        try:
                            dummy_molecule.update_connectivity_values()
                            dummy_molecule.update()
                            dummy_species_old = Species()
                            dummy_species_old.molecule.append(dummy_molecule)
                            dummy_species_old.generate_resonance_structures()
                            molecule_test.append(dummy_molecule.to_smiles())
                            H298_old.append(self.get_thermo_data(dummy_species_old).get_enthalpy(298.))
                        except rmgpy.exceptions.AtomTypeError or rmgpy.exceptions.ResonanceError:
                            #print('exception 1585')
                            pass
                #For bidentate CO because we want C[-1]#O[+1] but not .C#O.
                #and for bidentate OO we want [O][O] but not O[-1]#O[+1], see
                #line 1556
                if (bond.order == 3 and adsorbed_atoms[0].radical_electrons and 
                    adsorbed_atoms[1].radical_electrons and 
                    (adsorbed_atoms[0].lone_pairs or adsorbed_atoms[1].lone_pairs)):
                    #print('in here 1594')
                    adsorbed_atoms[0].decrement_radical()
                    adsorbed_atoms[1].decrement_radical()
                    if adsorbed_atoms[0].lone_pairs:
                        adsorbed_atoms[1].increment_lone_pairs()
                    else:
                        adsorbed_atoms[0].increment_lone_pairs()
                try:
                    dummy_molecule.update_connectivity_values()
                    dummy_molecule.update()
                    dummy_species_old = Species()
                    dummy_species_old.molecule.append(dummy_molecule)
                    dummy_species_old.generate_resonance_structures()
                    molecule_test.append(dummy_molecule.to_smiles())
                    H298_old.append(self.get_thermo_data(dummy_species_old).get_enthalpy(298.))
                except rmgpy.exceptions.AtomTypeError or rmgpy.exceptions.ResonanceError:
                    #print('in here 1619')
                    index_min_H298 = min(range(len(H298_old)), key=H298_old.__getitem__)
                    dummy_molecule.from_smiles(molecule_test[index_min_H298])
                    logging.exception("The above molecule was reestructured to: \n{0}".format(dummy_molecule_old.to_adjacency_list()))
                    pass

I think it is just the idea, I am not sure if it could work generally. And sorry about my awkward way to program, I am new in python.

rwest commented 3 years ago

Thanks for the help, @josefedupe! Maybe you could help us review pull request #2045, see if it works for your scenarios.

josefedupe commented 3 years ago

Thanks to you @rwest, RMG-CAT is an amazing code! The scenario was presented in #2045. And sorry for the late response!