Source code for qm_atlas.workflows.score_tautomers
import functools
import logging
from pathlib import Path
import numpy as np
from ppqm.utils import functools as ppqm_functools
from rdkit import Chem
from qm_atlas.constants import DEFAULT_SCR
from qm_atlas.tasks import calculate_properties, conformer_generation, optimize
from qm_atlas.workflows import fast_conformers
from qm_atlas.wrappers import unicon
from qm_atlas.wrappers.cosmotherm import cosmo_tasks, cosmo_utils
from qm_atlas.wrappers.turbomole import utils as tm_utils
_logger = logging.getLogger(__name__)
TAUTOMER_TM_SP_CONFIG = calculate_properties.TurbomoleSinglePointOptions(
basis="def2-mTZVP",
functional=["b97-3c"],
control_template=None,
)
[docs]
def generate_tautomers(
mol: Chem.Mol,
energy_threshold: float = 7,
solvent_name: str = "h2o",
single_point_config: calculate_properties.TurbomoleSinglePointOptions = TAUTOMER_TM_SP_CONFIG,
cosmo_options: cosmo_utils.CosmoSettings = cosmo_tasks.DESCRIPTORS_SETTINGS_SIMPLIFIED,
solvents_dirs: list[Path] | None = None,
use_config_cosmo_dir: bool = True,
conformer_generation_options: conformer_generation.ConformerGenerationOptions = fast_conformers.FAST_DEFAULT_OPTIONS,
optimization_config: list[
optimize.XtbOptions | optimize.JobexOptions | optimize.XtbTurbomoleOptions
] = fast_conformers.DEFAULT_OPTIMIZATION_CONFIG,
final_optimization_config: list[
optimize.XtbOptions | optimize.JobexOptions | optimize.XtbTurbomoleOptions
] = fast_conformers.DEFAULT_FINAL_OPTIMIZATION_CONFIG,
n_cores: int = 1,
scr: Path = DEFAULT_SCR,
) -> list[Chem.Mol]:
"""Enumerate the tautomers of a molecule and keep the energetically relevant ones.
Tautomers are generated with UniCon. Each tautomer is expanded into a conformer
ensemble and geometry-optimized with
:func:`~qm_atlas.workflows.fast_conformers.generate_fast_conformers`, and every
conformer receives a DFT single-point calculation from which COSMOtherm free
energies are extracted. Each tautomer is scored by the lowest ``E_COSMO+dE+Mu``
free energy across its conformers in ``solvent_name``. Tautomers whose free
energy lies within ``energy_threshold`` of the most stable tautomer are returned.
Args:
mol (Chem.Mol):
The molecule whose tautomers should be enumerated and scored.
energy_threshold (float):
Free-energy window in kcal/mol relative to the most stable tautomer.
Tautomers with a relative energy below this value are kept. Defaults to 7.
solvent_name (str):
Name of the solvent used for the COSMOtherm free-energy evaluation.
Defaults to "h2o".
single_point_config (TurbomoleSinglePointOptions):
Configuration for the DFT single-point calculations used to obtain the
COSMOtherm properties. Defaults to
:data:`~qm_atlas.workflows.score_tautomers.TAUTOMER_TM_SP_CONFIG`.
cosmo_options (cosmo_utils.CosmoSettings):
Settings for the COSMOtherm descriptor calculation. Defaults to
``cosmo_tasks.DESCRIPTORS_SETTINGS_SIMPLIFIED``.
solvents_dirs (list[Path] | None):
Additional directories in which to look for the solvent ``.cosmo`` files.
Defaults to None (use the COSMOtherm database and configured directories).
use_config_cosmo_dir (bool):
Whether to prepend the solvent directories from the software configuration
to ``solvents_dirs``. Defaults to True.
conformer_generation_options (ConformerGenerationOptions):
Force-field conformer generator settings used during conformer expansion.
Defaults to
:data:`~qm_atlas.workflows.fast_conformers.FAST_DEFAULT_OPTIONS`.
optimization_config (list[XtbOptions | JobexOptions | XtbTurbomoleOptions]):
Sequence of optimization steps applied during conformer expansion.
Defaults to
:data:`~qm_atlas.workflows.fast_conformers.DEFAULT_OPTIMIZATION_CONFIG`.
final_optimization_config (list[XtbOptions | JobexOptions | XtbTurbomoleOptions]):
Optional final optimization steps applied to the expanded conformers.
Defaults to
:data:`~qm_atlas.workflows.fast_conformers.DEFAULT_FINAL_OPTIMIZATION_CONFIG`.
n_cores (int):
Number of CPU cores to use. Conformer expansion is parallelized over the
tautomers. Defaults to 1.
scr (Path):
Scratch directory for the calculations. Defaults to DEFAULT_SCR.
Returns:
list[Chem.Mol]:
The tautomers (as generated by UniCon) whose lowest COSMOtherm free
energy lies within ``energy_threshold`` of the most stable tautomer.
"""
# generate tautomers using unicon
tautomers = unicon.get_tautomers(mol)
# TODO: explore using other tautomer generators here, not just unicon.
_logger.info(f"Unicon generated {len(tautomers)} tautomers")
conformer_generator = functools.partial(
fast_conformers.generate_fast_conformers,
conformer_generation_options=conformer_generation_options,
optimization_config=optimization_config,
final_optimization_config=final_optimization_config,
scr=scr,
)
# conformer expansion for all tautomers
n_procs = min(n_cores, len(tautomers))
if n_procs == 1:
tautomers_3d = [conformer_generator(tautomer, n_cores=n_cores) for tautomer in tautomers]
elif n_procs > 1:
tautomers_3d = ppqm_functools.func_parallel(
conformer_generator, tautomers, n_cores=n_procs, title="Conf_Exp_Tautomers"
)
else:
raise ValueError(f"Received {n_cores} for number of cores")
# DFT energies for all conformers
tautomer_energies = []
for tautomer_3d in tautomers_3d:
single_point_properties = calculate_properties.calculate_property_mol(
tautomer_3d,
single_point_config,
n_cores=n_cores,
scr=scr,
)
cosmo_outputs = [
calculate_properties.extract_cosmo_output(sp_prop, remove_from_dict=True)
for sp_prop in single_point_properties
]
_, cosmo_outputs, _ = tm_utils.extract_converged(tautomer_3d, cosmo_outputs)
# catch unconverged calculations
if len(cosmo_outputs) == 0:
tautomer_energies.append(np.inf)
smi = Chem.MolToSmiles(tautomer_3d)
_logger.error(f"Single Point Calculations failed for Tauromer {smi}")
continue
cosmo_df = cosmo_tasks.calculate_descriptors_df(
cosmo_outputs,
solvent_name=solvent_name,
level="bp-tzvp",
solvents_dirs=solvents_dirs,
options=cosmo_options,
use_config_cosmo_dir=use_config_cosmo_dir,
)
tautomer_energies.append(cosmo_df["E_COSMO+dE+Mu"].min())
energy_differences = np.array(tautomer_energies) - np.min(tautomer_energies)
stable_indices = energy_differences < energy_threshold
selected_tautomers = [tautomers[idx] for idx, stable in enumerate(stable_indices) if stable]
n = len(selected_tautomers)
_logger.info(f"Keeping {n} after applying energy cutoff: {energy_threshold} kcal/mol")
return selected_tautomers