Skip to content

Commit

Permalink
Estimating thermo of an adsorbate Species now re-orders Species.molec…
Browse files Browse the repository at this point in the history
…ule.

We do this in the gas phase, so now also adsorbates:
If you have a Species with multiple Molecules (resonance structures)
then the list of Molecules gets re-ordered when you estimate the thermo
so that the most stable (lowest H298) goes first.
This will then be used for drawing pictures, printing SMILES strings, etc.
which might mean it's no the one that you started with.
This may surprise some people in some situations, but I think is the 
way the gas phase works, and it has some benefits, especially for 
machine-generated species.

One thing we don't yet do in adsorbates is prioritize resonance structures coming
from thermo libraries, over those that are merely estimated.
Not sure of the consequences of this.

Other changes in this commit are mostly for clarity.
  • Loading branch information
rwest authored and bjkreitz committed Dec 4, 2024
1 parent 0d4099d commit efe86ca
Showing 1 changed file with 21 additions and 17 deletions.
38 changes: 21 additions & 17 deletions rmgpy/data/thermo.py
Original file line number Diff line number Diff line change
Expand Up @@ -1561,8 +1561,7 @@ def get_thermo_data_for_surface_species(self, species):
"""

# define the comparison function to find the lowest energy

def lowest_energy(species):
def species_enthalpy(species):
if hasattr(species.thermo, 'H298'):
return species.thermo.H298.value_si
else:
Expand All @@ -1577,9 +1576,7 @@ def lowest_energy(species):
resonance_data = []

# iterate over all resonance structures of the species

for i in range(len(species.molecule)):
molecule = species.molecule[i]
for molecule in species.molecule:
# store any labeled atoms to reapply at the end
labeled_atoms = molecule.get_all_labeled_atoms()
molecule.clear_labeled_atoms()
Expand All @@ -1591,7 +1588,10 @@ def lowest_energy(species):
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
# 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:
Expand All @@ -1605,21 +1605,20 @@ def lowest_energy(species):
gas_phase_species_estimates.append(dummy_species)

if gas_phase_species_from_libraries:
gas_phase_species = min(gas_phase_species_from_libraries, key=lowest_energy)
gas_phase_species = min(gas_phase_species_from_libraries, key=species_enthalpy)
else:
gas_phase_species = min(gas_phase_species_estimates, key=lowest_energy)
gas_phase_species = min(gas_phase_species_estimates, key=species_enthalpy)

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

if not isinstance(thermo, ThermoData):
thermo = thermo.to_thermo_data()
find_cp0_and_cpinf(gas_phase_species, thermo)

# 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)"),
Expand All @@ -1639,7 +1638,6 @@ def lowest_energy(species):
# (group_additivity=True means it appends the comments)
add_thermo_data(thermo, adsorption_thermo, group_additivity=True)

resonance_data.append(thermo)
# if the molecule had labels, reapply them
for label, atom in labeled_atoms.items():
if isinstance(atom,list):
Expand All @@ -1648,15 +1646,21 @@ def lowest_energy(species):
else:
atom.label = label

# 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
# a tuple of the H298, the ThermoData object, and the molecule
resonance_data.append((thermo.get_enthalpy(298), thermo, molecule))

# 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 = []
# first element of each tuple is H298, so sort by that
resonance_data = sorted(resonance_data)

for i in range(len(resonance_data)):
enthalpy298.append(resonance_data[i].get_enthalpy(298))
# reorder the resonance structures (molecules) by their H298
species.molecule = [mol for h, thermo, mol in resonance_data]

thermo = resonance_data[np.argmin(enthalpy298)]
thermo = resonance_data[0][1]

if thermo.label:
thermo.label += 'X' * len(surface_sites)
Expand Down

0 comments on commit efe86ca

Please sign in to comment.