Skip to content

Commit

Permalink
Adding metal attributes to thermo
Browse files Browse the repository at this point in the history 
The Pt111 values in the metal database are not the same as the binding energies used to calculate the Pt111 surface thermo library, so they are hard coded in here
  • Loading branch information
@mazeau
mazeau committed Oct 5, 2020
1 parent 768749f commit 0283e3d
Showing 1 changed file with 58 additions and 34 deletions.
92 changes: 58 additions & 34 deletions rmgpy/data/thermo.py
View file
Expand Up @@ -48,6 +48,8 @@
from rmgpy.molecule import Molecule, Bond, Group
from rmgpy.species import Species
from rmgpy.thermo import NASAPolynomial, NASA, ThermoData, Wilhoit
from rmgpy.data.surface import MetalDatabase
from rmgpy import settings

#: This dictionary is used to add multiplicity to species label
_multiplicity_labels = {1: 'S', 2: 'D', 3: 'T', 4: 'Q', 5: 'V'}
Expand Down Expand Up @@ -845,7 +847,7 @@ def __init__(self):
self.global_context = {}

# Catalyst properties
self.set_delta_atomic_adsorption_energies()
self.set_binding_energies()

def __reduce__(self):
"""
Expand Down Expand Up @@ -1192,7 +1194,7 @@ def record_ring_generic_nodes(self):
continue
self.groups['ring'].generic_nodes.append(label)

def get_thermo_data(self, species, training_set=None):
def get_thermo_data(self, species, metal_to_scale_from=None, metal_to_scale_to=None, training_set=None):
"""
Return the thermodynamic parameters for a given :class:`Species`
object `species`. This function first searches the loaded libraries
Expand All @@ -1202,25 +1204,31 @@ def get_thermo_data(self, species, training_set=None):
The method corrects for symmetry when the molecule uses machine
learning or group additivity. Libraries and direct QM calculations
are already corrected.
If either metal to scale to or from is not specified, assume the binding energies given in the input file
Returns: ThermoData
"""
from rmgpy.rmg.input import get_input

thermo0 = self.get_thermo_data_from_libraries(species)

if thermo0 is not None:
if thermo0 is not None: # was able to find thermodata in the loaded libraries
if len(thermo0) != 3:
raise RuntimeError("thermo0 should be a tuple (thermo_data, library, entry), not {0}".format(thermo0))
thermo0 = thermo0[0]

if species.contains_surface_site():
thermo0 = self.correct_binding_energy(thermo0, species)
if 'surface' or 'adsorption' not in thermo0.comment: # the first time it is being loaded in?
thermo0 = self.correct_binding_energy(thermo0, species, metal_to_scale_from=thermo0.metal, metal_to_scale_to=metal_to_scale_to)
else: # if the thermo has already been corrected, scale from the input binding energies
# if the thermo has already been corrected, scale from that point
thermo0 = self.correct_binding_energy(thermo0, species, metal_to_scale_from=None, metal_to_scale_to=metal_to_scale_to)
return thermo0

if species.contains_surface_site():
thermo0 = self.get_thermo_data_for_surface_species(species)
thermo0 = self.correct_binding_energy(thermo0, species)
thermo0 = self.correct_binding_energy(thermo0, species, metal_to_scale_to=metal_to_scale_to)
return thermo0

try:
Expand Down Expand Up @@ -1351,54 +1359,70 @@ def get_thermo_data(self, species, training_set=None):
# Return the resulting thermo parameters
return thermo0

def set_delta_atomic_adsorption_energies(self, binding_energies=None):
def set_binding_energies(self, binding_energies=None):
"""
Sets and stores the change in atomic binding energy between
the desired and the Pt(111) default.
This depends on the two metal surfaces: the reference one used in
the database of adsorption energies, and the desired surface.
If binding_energies are not provided, resets the values to those
of the Pt(111) default.
Sets and stores the atomic binding energies specified in the input file.
Args:
binding_energies (dict, optional): the desired binding energies with
elements as keys and binding energy/unit tuples as values
Returns:
None, stores result in self.delta_atomic_adsorption_energy
None, stores result in self.binding_energies
"""
reference_binding_energies = {
'C': rmgpy.quantity.Energy(-7.025, 'eV/molecule'),
'H': rmgpy.quantity.Energy(-2.754, 'eV/molecule'),
'O': rmgpy.quantity.Energy(-3.811, 'eV/molecule'),
'N': rmgpy.quantity.Energy(-4.632, 'eV/molecule'),
}
metal_db = MetalDatabase()
metal_db.load(os.path.join(settings['database.directory'], 'surface'))

# Use Pt(111) reference if no binding energies are provided
if binding_energies is None:
binding_energies = reference_binding_energies
# Use Pt(111) reference that surfaceThermoPt111 was calculated on if no binding energies are provided
binding_energies = {
'C': rmgpy.quantity.Energy(-7.025, 'eV/molecule'),
'H': rmgpy.quantity.Energy(-2.754, 'eV/molecule'),
'O': rmgpy.quantity.Energy(-3.811, 'eV/molecule'),
'N': rmgpy.quantity.Energy(-4.632, 'eV/molecule'),
}
elif binding_energies is str:
# Want to load the binding energies from the database
metal_db.find_binding_energies(binding_energies)

self.delta_atomic_adsorption_energy = {
'C': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
'H': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
'O': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
'N': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
}
for element, energy in binding_energies.items():
binding_energies[element] = rmgpy.quantity.Energy(binding_energies[element])

for element, deltaEnergy in self.delta_atomic_adsorption_energy.items():
deltaEnergy.value_si = binding_energies[element].value_si - reference_binding_energies[element].value_si
self.binding_energies = binding_energies

def correct_binding_energy(self, thermo, species):
def correct_binding_energy(self, thermo, species, metal_to_scale_to=None, metal_to_scale_from=None):
"""
Changes the provided thermo, by applying a linear scaling relation
to correct the adsorption energy.
:param thermo: starting thermo data
:param species: the species (which is an adsorbate)
:param metal: a tuple of the metal you're scaling from and the metal you want to scale to
:return: corrected thermo
"""
metal_db = MetalDatabase()
metal_db.load(os.path.join(settings['database.directory'], 'surface'))

if metal_to_scale_to is None:
metal_to_scale_to_binding_energies = self.binding_energies
else:
metal_to_scale_to_binding_energies = metal_db.find_binding_energies(metal_to_scale_to)

if metal_to_scale_from is None:
metal_to_scale_to_binding_energies = self.binding_energies
else:
metal_to_scale_from_binding_energies = metal_db.find_binding_energies(metal_to_scale_from)

delta_atomic_adsorption_energy = {
'C': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
'H': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
'O': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
'N': rmgpy.quantity.Energy(0.0, 'eV/molecule'),
}

for element, deltaEnergy in delta_atomic_adsorption_energy.items():
deltaEnergy.value_si = metal_to_scale_to_binding_energies[element].value_si - metal_to_scale_from_binding_energies[element].value_si

molecule = species.molecule[0]
# only want/need to do one resonance structure
surface_sites = []
Expand Down Expand Up @@ -1438,7 +1462,7 @@ def correct_binding_energy(self, thermo, species):
comments = []
for element in 'CHON':
if normalized_bonds[element]:
change_in_binding_energy = self.delta_atomic_adsorption_energy[element].value_si * normalized_bonds[element]
change_in_binding_energy = delta_atomic_adsorption_energy[element].value_si * normalized_bonds[element]
thermo.H298.value_si += change_in_binding_energy
comments.append(f'{normalized_bonds[element]:.2f}{element}')
thermo.comment += " Binding energy corrected by LSR ({})".format('+'.join(comments))
Expand Down Expand Up @@ -1631,7 +1655,7 @@ def get_thermo_data_from_libraries(self, species, training_set=None):

return None

def get_all_thermo_data(self, species):
def get_all_thermo_data(self, species, metal=None):
"""
Return all possible sets of thermodynamic parameters for a given
:class:`Species` object `species`. The hits from the depository come
Expand Down

0 comments on commit 0283e3d

Please sign in to comment.