Skip to content

Commit

Permalink
add resonance to get_thermo_for_surface_species
Browse files Browse the repository at this point in the history
  • Loading branch information
@bjkreitz
bjkreitz committed Sep 6, 2023
1 parent 8e2bee1 commit da17918
Showing 1 changed file with 77 additions and 58 deletions.
135 changes: 77 additions & 58 deletions rmgpy/data/thermo.py
View file
Expand Up @@ -1542,77 +1542,96 @@ def get_thermo_data_for_surface_species(self, species):
Returns a :class:`ThermoData` object, with no Cp0 or CpInf
"""

# define the comparison function to find the lowest energy

def lowest_energy(species):
if hasattr(species.thermo, 'H298'):
return species.thermo.H298.value_si
else:
return species.thermo.get_enthalpy(298.0)

if species.is_surface_site():
raise DatabaseError("Can't estimate thermo of vacant site. Should be in library (and should be 0).")

logging.debug("Trying to generate thermo for surface species using first of %d resonance isomer(s):",
len(species.molecule))
molecule = species.molecule[0]
# store any labeled atoms to reapply at the end
labeled_atoms = molecule.get_all_labeled_atoms()
molecule.clear_labeled_atoms()
logging.debug("Before removing from surface:\n" + molecule.to_adjacency_list())
# only want/need to do one resonance structure,
# because will need to regenerate others in gas phase
dummy_molecules = molecule.get_desorbed_molecules()
for mol in dummy_molecules:
mol.clear_labeled_atoms()
if len(dummy_molecules) == 0:
raise RuntimeError(f"Cannot get thermo for gas-phase molecule. No valid dummy molecules from original molecule:\n{molecule.to_adjacency_list()}")

resonance_data = []

# iterate over all resonance structures of the species

for i in range(len(species.molecule)):
molecule = species.molecule[i]
# store any labeled atoms to reapply at the end
labeled_atoms = molecule.get_all_labeled_atoms()
molecule.clear_labeled_atoms()
logging.debug("Before removing from surface:\n" + molecule.to_adjacency_list())
# get all desorbed molecules for a given adsorbate
dummy_molecules = molecule.get_desorbed_molecules()
for mol in dummy_molecules:
mol.clear_labeled_atoms()
if len(dummy_molecules) == 0:
raise RuntimeError(f"Cannot get thermo for gas-phase molecule. No valid dummy molecules from original molecule:\n{molecule.to_adjacency_list()}")

# if len(molecule) > 1, it will assume all resonance structures have already been generated when it tries to generate them, so evaluate each configuration separately and pick the lowest energy one by H298 value
gas_phase_species_from_libraries = []
gas_phase_species_estimates = []
for dummy_molecule in dummy_molecules:
dummy_species = Species()
dummy_species.molecule = [dummy_molecule]
dummy_species.generate_resonance_structures()
dummy_species.thermo = self.get_thermo_data(dummy_species)
if dummy_species.thermo.label:
gas_phase_species_from_libraries.append(dummy_species)
else:
gas_phase_species_estimates.append(dummy_species)
# if len(molecule) > 1, it will assume all resonance structures have already been generated when it tries to generate them, so evaluate each configuration separately and pick the lowest energy one by H298 value
gas_phase_species_from_libraries = []
gas_phase_species_estimates = []
for dummy_molecule in dummy_molecules:
dummy_species = Species()
dummy_species.molecule = [dummy_molecule]
dummy_species.generate_resonance_structures()
dummy_species.thermo = self.get_thermo_data(dummy_species)
if dummy_species.thermo.label:
gas_phase_species_from_libraries.append(dummy_species)
else:
gas_phase_species_estimates.append(dummy_species)

# define the comparison function to find the lowest energy
def lowest_energy(species):
if hasattr(species.thermo, 'H298'):
return species.thermo.H298.value_si
if gas_phase_species_from_libraries:
gas_phase_species = min(gas_phase_species_from_libraries, key=lowest_energy)
else:
return species.thermo.get_enthalpy(298.0)
gas_phase_species = min(gas_phase_species_estimates, key=lowest_energy)

if gas_phase_species_from_libraries:
species = min(gas_phase_species_from_libraries, key=lowest_energy)
else:
species = min(gas_phase_species_estimates, key=lowest_energy)
thermo = gas_phase_species.thermo
thermo.comment = f"Gas phase thermo for {thermo.label or gas_phase_species.molecule[i].to_smiles()} from {thermo.comment}. Adsorption correction:"
logging.debug("Using thermo from gas phase for species {}\n".format(gas_phase_species.label) + repr(thermo))

thermo = species.thermo
thermo.comment = f"Gas phase thermo for {thermo.label or species.molecule[0].to_smiles()} from {thermo.comment}. Adsorption correction:"
logging.debug("Using thermo from gas phase for species {}\n".format(species.label) + repr(thermo))
if not isinstance(thermo, ThermoData):
thermo = thermo.to_thermo_data()
find_cp0_and_cpinf(gas_phase_species, thermo)

if not isinstance(thermo, ThermoData):
thermo = thermo.to_thermo_data()
find_cp0_and_cpinf(species, thermo)
# Get the adsorption energy
# Create the ThermoData object

# Get the adsorption energy
# Create the ThermoData object
adsorption_thermo = ThermoData(
Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
Cpdata=([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "J/(mol*K)"),
H298=(0.0, "kJ/mol"),
S298=(0.0, "J/(mol*K)"),
)
adsorption_thermo = ThermoData(
Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
Cpdata=([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "J/(mol*K)"),
H298=(0.0, "kJ/mol"),
S298=(0.0, "J/(mol*K)"),
)

surface_sites = molecule.get_surface_sites()
try:
self._add_adsorption_correction(adsorption_thermo, self.groups['adsorptionPt111'], molecule, surface_sites)
except (KeyError, DatabaseError):
logging.error("Couldn't find in adsorption thermo database:")
logging.error(molecule)
logging.error(molecule.to_adjacency_list())
raise
surface_sites = molecule.get_surface_sites()
try:
self._add_adsorption_correction(adsorption_thermo, self.groups['adsorptionPt111'], molecule, surface_sites)
except (KeyError, DatabaseError):
logging.error("Couldn't find in adsorption thermo database:")
logging.error(molecule)
logging.error(molecule.to_adjacency_list())
raise

# (group_additivity=True means it appends the comments)
add_thermo_data(thermo, adsorption_thermo, group_additivity=True)

resonance_data.append(thermo)

# (group_additivity=True means it appends the comments)
add_thermo_data(thermo, adsorption_thermo, group_additivity=True)
# Get the enthalpy of formation of every adsorbate at 298 K and determine the resonance structure with the lowest enthalpy of formation
# We assume that the lowest enthalpy of formation is the correct thermodynamic property for the adsorbate

enthalpy298 = []

for i in range(len(resonance_data)):
enthalpy298.append(resonance_data[i].get_enthalpy(298))

thermo = resonance_data[np.argmin(enthalpy298)]

if thermo.label:
thermo.label += 'X' * len(surface_sites)
Expand Down Expand Up @@ -2849,7 +2868,6 @@ def register_in_central_thermo_db(self, species):
except ValueError:
logging.info('Fail to generate inchi/smiles for species below:\n{0}'.format(species.to_adjacency_list()))


def find_cp0_and_cpinf(species, heat_capacity):
"""
Calculate the Cp0 and CpInf values, and add them to the HeatCapacityModel object.
Expand All @@ -2860,3 +2878,4 @@ def find_cp0_and_cpinf(species, heat_capacity):
if heat_capacity.CpInf is None:
cp_inf = species.calculate_cpinf()
heat_capacity.CpInf = (cp_inf, "J/(mol*K)")

0 comments on commit da17918

Please sign in to comment.