Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improved TS energy checks #690

Merged
merged 13 commits into from
Oct 2, 2023
Merged

Improved TS energy checks #690

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
297ec6a
Minor: Simplified Species.populate_ts_checks()
alongd Aug 3, 2023
2f6a6f8
Refactor check_rxn_e_elect() and consolidate rxn energy report
alongd Aug 24, 2023
4f6eb21
Tests: Refactor check_rxn_e_elect()
alongd Aug 3, 2023
802733e
Refactor compute_rxn_e0()
alongd Aug 22, 2023
2b852e3
Tests: Modify the compute_rxn_e0() test
alongd Aug 3, 2023
1f700ac
Modified check_rxn_e0() to change test status w/o returning a value
alongd Aug 3, 2023
c3460c4
Tests: Added check_rxn_e0() test
alongd Aug 20, 2023
3b21f71
Call all TS energy tests from the same place
alongd Aug 5, 2023
d1417c6
Check reaction E0 through check_ts()
alongd Aug 5, 2023
999cd6b
Check species yml_path in sp/geo/freq checks in Scheduler
alongd Aug 4, 2023
397e303
Added species_has_sp_and_freq() into Scheduler
alongd Aug 4, 2023
659bf33
Added copy_e0_values() method to ARCReacrtion
alongd Aug 5, 2023
e8e3fe4
Tests: Scheduler species_has_x() tests for species with `yml_path` attr
alongd Aug 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
177 changes: 101 additions & 76 deletions arc/checks/ts.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,10 +37,16 @@


def check_ts(reaction: 'ARCReaction',
verbose: bool = True,
job: Optional['JobAdapter'] = None,
checks: Optional[List[str]] = None,
rxn_zone_atom_indices: Optional[List[int]] = None,
species_dict: Optional[dict] = None,
project_directory: Optional[str] = None,
kinetics_adapter: Optional[str] = None,
output: Optional[dict] = None,
sp_level: Optional['Level'] = None,
freq_scale_factor: float = 1.0,
verbose: bool = True,
):
"""
Check the TS in terms of energy, normal mode displacement, and IRC.
Expand All @@ -52,19 +58,36 @@

Args:
reaction (ARCReaction): The reaction for which the TS is checked.
verbose (bool, optional): Whether to print logging messages.
job (JobAdapter, optional): The frequency job object instance.
checks (List[str], optional): Specific checks to run. Optional values: 'energy', 'freq', 'IRC', 'rotors'.
rxn_zone_atom_indices (List[int], optional): The 0-indices of atoms identified by the normal mode displacement
as the reaction zone. Automatically determined if not given.
species_dict (dict, optional): The Scheduler species dictionary.
project_directory (str, optional): The path to ARC's project directory.
kinetics_adapter (str, optional): The statmech software to use for kinetic rate coefficient calculations.
output (dict, optional): The Scheduler output dictionary.
sp_level (Level, optional): The single-point energy level of theory.
freq_scale_factor (float, optional): The frequency scaling factor.
verbose (bool, optional): Whether to print logging messages.
"""
checks = checks or list()
for entry in checks:
if entry not in ['energy', 'freq', 'IRC', 'rotors']:
raise ValueError(f"Requested checks could be 'energy', 'freq', 'IRC', or 'rotors', got:\n{checks}")

if 'energy' in checks or not reaction.ts_species.ts_checks['e_elect']:
check_ts_energy(reaction=reaction, verbose=verbose)
if 'energy' in checks:
if not reaction.ts_species.ts_checks['E0']:
rxn_copy = compute_rxn_e0(reaction=reaction,
species_dict=species_dict,
project_directory=project_directory,
kinetics_adapter=kinetics_adapter,
output=output,
sp_level=sp_level,
freq_scale_factor=freq_scale_factor)
reaction.copy_e0_values(rxn_copy)
check_rxn_e0(reaction=reaction, verbose=verbose)
if reaction.ts_species.ts_checks['E0'] is None and not reaction.ts_species.ts_checks['e_elect']:
check_rxn_e_elect(reaction=reaction, verbose=verbose)

Check warning on line 90 in arc/checks/ts.py

View check run for this annotation

Codecov / codecov/patch

arc/checks/ts.py#L90 

Added line #L90 was not covered by tests

if 'freq' in checks or (not reaction.ts_species.ts_checks['normal_mode_displacement'] and job is not None):
check_normal_mode_displacement(reaction, job=job)
Expand Down Expand Up @@ -99,9 +122,9 @@
return True


def check_ts_energy(reaction: 'ARCReaction',
verbose: bool = True,
) -> None:
def check_rxn_e_elect(reaction: 'ARCReaction',
verbose: bool = True,
) -> None:
"""
Check that the TS electronic energy is above both reactant and product wells in a ``reaction``.
Sets the respective energy parameter 'e_elect' in the ``TS.ts_checks`` dictionary.
Expand All @@ -110,60 +133,45 @@
reaction (ARCReaction): The reaction for which the TS is checked.
verbose (bool, optional): Whether to print logging messages.
"""
# Check whether E0 values are already known, e.g. from Arkane species YAML files
check_rxn_e0(reaction=reaction)
check_rxn_e0(reaction=reaction, verbose=verbose)
if reaction.ts_species.ts_checks['E0']:
return

r_e_elect = None if any([spc.e_elect is None for spc in reaction.r_species]) \
else sum(spc.e_elect * reaction.get_species_count(species=spc, well=0) for spc in reaction.r_species)
p_e_elect = None if any([spc.e_elect is None for spc in reaction.p_species]) \
else sum(spc.e_elect * reaction.get_species_count(species=spc, well=1) for spc in reaction.p_species)
ts_e_elect = reaction.ts_species.e_elect

if verbose and all([val is not None for val in [r_e_elect, p_e_elect, ts_e_elect]]):
min_e = extremum_list([r_e_elect, p_e_elect, ts_e_elect], return_min=True)
r_text = f'{r_e_elect - min_e:.2f} kJ/mol' if r_e_elect is not None else 'None'
ts_text = f'{ts_e_elect - min_e:.2f} kJ/mol' if ts_e_elect is not None else 'None'
p_text = f'{p_e_elect - min_e:.2f} kJ/mol' if p_e_elect is not None else 'None'
logger.info(
f'\nReaction {reaction.label} (TS {reaction.ts_label}, TSG {reaction.ts_species.chosen_ts}) '
f'has the following path electronic energy:\n'
f'Reactants: {r_text}\n'
f'TS: {ts_text}\n'
f'Products: {p_text}')

if all([val is not None for val in [r_e_elect, p_e_elect, ts_e_elect]]):
# We have all params, we can make a quantitative decision.
if ts_e_elect > r_e_elect + 1.0 and ts_e_elect > p_e_elect + 1.0:
# TS is above both wells.
reaction.ts_species.ts_checks['e_elect'] = True
return
# TS is not above both wells.
if verbose:
logger.error(f'TS of reaction {reaction.label} has a lower electronic energy value than expected.')
reaction.ts_species.ts_checks['e_elect'] = False
return
# We don't have any params (some are ``None``)
r_ee = sum_list_entries([r.e_elect for r in reaction.r_species],
multipliers=[reaction.get_species_count(species=r, well=0) for r in reaction.r_species])
p_ee = sum_list_entries([p.e_elect for p in reaction.p_species],
multipliers=[reaction.get_species_count(species=p, well=1) for p in reaction.p_species])
ts_ee = reaction.ts_species.e_elect
if verbose:
report_ts_and_wells_energy(r_e=r_ee, p_e=p_ee, ts_e=ts_ee, rxn_label=reaction.label,
ts_label=reaction.ts_label, chosen_ts=reaction.ts_species.chosen_ts,
energy_type='electronic energy')
if all([val is not None for val in [r_ee, p_ee, ts_ee]]):
if ts_ee > r_ee + 1.0 and ts_ee > p_ee + 1.0:
reaction.ts_species.ts_checks['e_elect'] = True
return
if verbose:
logger.error(f'TS of reaction {reaction.label} has a lower electronic energy value than expected.')
reaction.ts_species.ts_checks['e_elect'] = False
return
if verbose:
alongd marked this conversation as resolved.
Show resolved Hide resolved
logger.info('\n')
logger.warning(f"Could not get electronic energy for all species in reaction {reaction.label}.\n")
# We don't really know.
reaction.ts_species.ts_checks['e_elect'] = None
if 'Could not determine TS e_elect relative to the wells; ' not in reaction.ts_species.ts_checks['warnings']:
reaction.ts_species.ts_checks['warnings'] += 'Could not determine TS e_elect relative to the wells; '


def compute_and_check_rxn_e0(reaction: 'ARCReaction',
species_dict: dict,
project_directory: str,
kinetics_adapter: str,
output: dict,
sp_level: 'Level',
freq_scale_factor: float = 1.0,
) -> Optional[bool]:
def compute_rxn_e0(reaction: 'ARCReaction',
species_dict: dict,
project_directory: str,
kinetics_adapter: str,
output: dict,
sp_level: 'Level',
freq_scale_factor: float = 1.0,
) -> Optional['ARCReaction']:
"""
Checking the E0 values between wells and a TS in a ``reaction`` using ZPE from statmech.
This function computed E0 values and populates them in a copy of the given reaction instance.

Args:
reaction (ARCReaction): The reaction for which the TS is checked.
Expand All @@ -175,9 +183,11 @@
freq_scale_factor (float, optional): The frequency scaling factor.

Returns:
Optional[bool]: Whether the test failed and the scheduler should switch to a different TS guess,
``None`` if the test couldn't be performed.
Optional['ARCReaction']: A copy of the reaction object with E0 values populated.
"""
if any(val is None for val in [species_dict, project_directory, kinetics_adapter,
output, sp_level, freq_scale_factor]):
return None

Check warning on line 190 in arc/checks/ts.py

View check run for this annotation

Codecov / codecov/patch

arc/checks/ts.py#L190 

Added line #L190 was not covered by tests
for spc in reaction.r_species + reaction.p_species + [reaction.ts_species]:
folder = 'rxns' if species_dict[spc.label].is_ts else 'Species'
freq_path = os.path.join(project_directory, 'output', folder, spc.label, 'geometry', 'freq.out')
Expand All @@ -201,53 +211,68 @@
e0_only=True,
skip_rotors=True,
)
e0_pass = check_rxn_e0(reaction=rxn_copy)
if not e0_pass:
if rxn_copy.ts_species.ts_guesses_exhausted:
return False
return True # Switch TS.
return False # Don't switch TS.
return rxn_copy
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is not clear from the name and use of this function that it should return a reaction object. The name contains compute, so it should either return a computed result, or populate an object property. Could you elaborate and consider renaming?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a comment in the docstring. It's indeed unclear from the name, but the use does show it, also the type hint and now the docstring.



def check_rxn_e0(reaction: 'ARCReaction',
verbose: bool = True,
) -> Optional[bool]:
):
"""
Check the E0 values between wells and a TS in a ``reaction``, assuming that E0 values are available.

Args:
reaction (ARCReaction): The reaction to consider.
verbose (bool, optional): Whether to print logging messages.

Returns:
Optional[bool]: Whether the Ts E0 is above both R and P E0 wells.
"""
if reaction.ts_species.ts_checks['E0']:
return True
return

Check warning on line 228 in arc/checks/ts.py

View check run for this annotation

Codecov / codecov/patch

arc/checks/ts.py#L228 

Added line #L228 was not covered by tests
r_e0 = sum_list_entries([r.e0 for r in reaction.r_species],
multipliers=[reaction.get_species_count(species=r, well=0) for r in reaction.r_species])
p_e0 = sum_list_entries([p.e0 for p in reaction.p_species],
multipliers=[reaction.get_species_count(species=p, well=1) for p in reaction.p_species])
ts_e0 = reaction.ts_species.e0
if any(e0 is None for e0 in [r_e0, p_e0, ts_e0]):
reaction.ts_species.ts_checks['E0'] = None
else:
if verbose:
report_ts_and_wells_energy(r_e=r_e0, p_e=p_e0, ts_e=ts_e0, rxn_label=reaction.label, ts_label=reaction.ts_label,
chosen_ts=reaction.ts_species.chosen_ts, energy_type='E0')
if r_e0 >= ts_e0 or p_e0 >= ts_e0:
reaction.ts_species.ts_checks['E0'] = False
else:
reaction.ts_species.ts_checks['E0'] = True


def report_ts_and_wells_energy(r_e: float,
p_e: float,
ts_e: float,
rxn_label: str,
ts_label: str,
chosen_ts: int,
energy_type: str = 'electronic energy',
):
"""
Report the relative R/TS/P energies.

if verbose and all([val is not None for val in [r_e0, p_e0, ts_e0]]):
min_e0 = extremum_list([r_e0, p_e0, ts_e0], return_min=True)
r_text = f'{r_e0 - min_e0:.2f} kJ/mol' if r_e0 is not None else 'None'
ts_text = f'{ts_e0 - min_e0:.2f} kJ/mol' if ts_e0 is not None else 'None'
p_text = f'{p_e0 - min_e0:.2f} kJ/mol' if p_e0 is not None else 'None'
Args:
r_e (float): The reactant energy.
p_e (float): The product energy.
ts_e (float): The TS energy.
rxn_label (str): The reaction label.
ts_label (str): The TS label.
chosen_ts (int): The TSG number.
energy_type (str): The energy type: 'electronic energy' or 'E0'.
"""
if all([val is not None for val in [r_e, p_e, ts_e]]):
min_e = extremum_list(lst=[r_e, p_e, ts_e], return_min=True)
r_text = f'{r_e - min_e:.2f} kJ/mol'
ts_text = f'{ts_e - min_e:.2f} kJ/mol'
p_text = f'{p_e - min_e:.2f} kJ/mol'
logger.info(
f'\nReaction {reaction.label} (TS {reaction.ts_label}, TSG {reaction.ts_species.chosen_ts}) '
f'has the following path E0 values:\n'
f'\nReaction {rxn_label} (TS {ts_label}, TSG {chosen_ts}) has the following {energy_type} values:\n'
f'Reactants: {r_text}\n'
f'TS: {ts_text}\n'
f'Products: {p_text}')
if any(e0 is None for e0 in [r_e0, p_e0, ts_e0]):
return None
if r_e0 >= ts_e0 or p_e0 >= ts_e0:
reaction.ts_species.ts_checks['E0'] = False
return False
reaction.ts_species.ts_checks['E0'] = True
return True


def check_normal_mode_displacement(reaction: 'ARCReaction',
Expand Down
Loading
Loading