Add perturbation-based quasi-IRC workflows for ORCA and Q-Chem

docs/about/contributors.md

@@ -26,6 +26,7 @@ Additional contributions were made by the individuals listed [here](https://gith
 - [@ViktoriiaBaib](https://github.com/ViktoriiaBaib): Initial testing of quacc.
 - [@zulissimeta](https://github.com/zulissimeta): Dask support
 - [@yw-fang](https://github.com/yw-fang): VASP Non-SCF recipe
+- [@espottesmith](http://github.com/espottesmith): Some ORCA and Q-Chem workflows, mostly re: IRC
 
 ## Inspiration
 

docs/user/recipes/recipes_list.md

@@ -140,11 +140,14 @@ The list of available quacc recipes is shown below. The "Req'd Extras" column sp
 
 <center>
 
-| Name           | Decorator       | Documentation                             | Req'd Extras |
-| -------------- | --------------- | ----------------------------------------- | ------------ |
-| ORCA Static    | `#!Python @job` | [quacc.recipes.orca.core.static_job][]    |              |
-| ORCA Relax     | `#!Python @job` | [quacc.recipes.orca.core.relax_job][]     |              |
-| ORCA ASE Relax | `#!Python @job` | [quacc.recipes.orca.core.ase_relax_job][] |              |
+| Name                       | Decorator       | Documentation                                         | Req'd Extras |
+| -------------------------- | --------------- | ----------------------------------------------------- | ------------ |
+| ORCA Static                | `#!Python @job` | [quacc.recipes.orca.core.static_job][]                |              |
+| ORCA Relax                 | `#!Python @job` | [quacc.recipes.orca.core.relax_job][]                 |              |
+| ORCA Freq                  | `#!Python @job` | [quacc.recipes.orca.core.freq_job][]                  |              |
+| ORCA ASE Relax             | `#!Python @job` | [quacc.recipes.orca.core.ase_relax_job][]             |              |
+| ORCA ASE Quasi-IRC Perturb | `#!Python @job` | [quacc.recipes.orca.core.ase_quasi_irc_perturb_job][] |              |
+
 
 </center>
 
@@ -170,14 +173,15 @@ The list of available quacc recipes is shown below. The "Req'd Extras" column sp
 
 <center>
 
-| Name             | Decorator       | Documentation                            | Req'd Extras   |
-| ---------------- | --------------- | ---------------------------------------- | -------------- |
-| Q-Chem Static    | `#!Python @job` | [quacc.recipes.qchem.core.static_job][]  |                |
-| Q-Chem Relax     | `#!Python @job` | [quacc.recipes.qchem.core.relax_job][]   |                |
-| Q-Chem Frequency | `#!Python @job` | [quacc.recipes.qchem.core.freq_job][]    |                |
-| Q-Chem TS        | `#!Python @job` | [quacc.recipes.qchem.ts.ts_job][]        | `quacc[sella]` |
-| Q-Chem IRC       | `#!Python @job` | [quacc.recipes.qchem.ts.irc_job][]       | `quacc[sella]` |
-| Q-Chem Quasi IRC | `#!Python @job` | [quacc.recipes.qchem.ts.quasi_irc_job][] | `quacc[sella]` |
+| Name                     | Decorator       | Documentation                                    | Req'd Extras   |
+| ------------------------ | --------------- | ------------------------------------------------ | -------------- |
+| Q-Chem Static            | `#!Python @job` | [quacc.recipes.qchem.core.static_job][]          |                |
+| Q-Chem Relax             | `#!Python @job` | [quacc.recipes.qchem.core.relax_job][]           |                |
+| Q-Chem Frequency         | `#!Python @job` | [quacc.recipes.qchem.core.freq_job][]            |                |
+| Q-Chem TS                | `#!Python @job` | [quacc.recipes.qchem.ts.ts_job][]                | `quacc[sella]` |
+| Q-Chem IRC               | `#!Python @job` | [quacc.recipes.qchem.ts.irc_job][]               | `quacc[sella]` |
+| Q-Chem Quasi IRC         | `#!Python @job` | [quacc.recipes.qchem.ts.quasi_irc_job][]         | `quacc[sella]` |
+| Q-Chem Quasi IRC Perturb | `#!Python @job` | [quacc.recipes.qchem.ts.quasi_irc_perturb_job][] | `quacc[sella]` |
 
 </center>
 

src/quacc/atoms/core.py

@@ -15,6 +15,7 @@
 if TYPE_CHECKING:
     from ase.atoms import Atoms
     from ase.optimize.optimize import Dynamics
+    from numpy.typing import NDArray
 
 logger = logging.getLogger(__name__)
 
@@ -254,3 +255,34 @@ def get_final_atoms_from_dynamics(dynamics: Dynamics) -> Atoms:
     return (
         dynamics.atoms.atoms if isinstance(dynamics.atoms, Filter) else dynamics.atoms
     )
+
+
+def perturb(mol: Atoms, matrix: list[list[float]] | NDArray, scale: float) -> Atoms:
+    """
+    Perturb each atom in a molecule by a (scaled) 1x3 vector, reflecting e.g. a vibrational normal mode.
+
+    Parameters
+    ----------
+    mol
+        ASE Atoms object representing a molecule
+    matrix
+        Nx3 matrix, where N is the numebr of atoms. This means that there is potentially a different translation
+        vector for each atom in the molecule.
+    scale
+        Scaling factor for perturbation
+
+    Returns
+    -------
+    Atoms
+        The input molecule after perturbation
+    """
+
+    mol_copy = copy_atoms(mol)
+    mode = np.asarray(matrix)
+
+    orig_pos = mol_copy.get_positions()
+
+    pos = orig_pos + mode * scale
+    mol_copy.set_positions(pos)
+
+    return mol_copy

src/quacc/recipes/orca/core.py

@@ -7,12 +7,14 @@
 import psutil
 
 from quacc import job
+from quacc.atoms.core import perturb
 from quacc.recipes.orca._base import run_and_summarize, run_and_summarize_opt
 
 if TYPE_CHECKING:
     from typing import Any, Literal
 
     from ase.atoms import Atoms
+    from numpy.typing import NDArray
 
     from quacc.schemas._aliases.cclib import cclibSchema
     from quacc.utils.files import Filenames, SourceDirectory
@@ -81,6 +83,74 @@ def static_job(
     )
 
 
+@job
+def freq_job(
+    atoms: Atoms,
+    charge: int = 0,
+    spin_multiplicity: int = 1,
+    xc: str = "wb97x-d3bj",
+    basis: str = "def2-tzvp",
+    numerical: bool = False,
+    orcasimpleinput: list[str] | None = None,
+    orcablocks: list[str] | None = None,
+    nprocs: int | Literal["max"] = "max",
+    copy_files: SourceDirectory | dict[SourceDirectory, Filenames] | None = None,
+) -> cclibSchema:
+    """
+    Carry out a vibrational frequency analysis calculation.
+
+    Parameters
+    ----------
+    atoms
+        Atoms object
+    charge
+        Charge of the system.
+    spin_multiplicity
+        Multiplicity of the system.
+    xc
+        Exchange-correlation functional
+    basis
+        Basis set
+    numerical
+        If True (default False), a numeric frequency calculation will be requested
+    orcasimpleinput
+        List of `orcasimpleinput` swaps for the calculator. To remove entries
+        from the defaults, put a `#` in front of the name. Refer to the
+        [ase.calculators.orca.ORCA][] calculator for details on `orcasimpleinput`.
+    orcablocks
+        List of `orcablocks` swaps for the calculator. To remove entries
+        from the defaults, put a `#` in front of the name. Refer to the
+        [ase.calculators.orca.ORCA][] calculator for details on `orcablocks`.
+    nprocs
+        Number of processors to use. Defaults to the number of physical cores.
+    copy_files
+        Files to copy (and decompress) from source to the runtime directory.
+
+    Returns
+    -------
+    cclibSchema
+        Dictionary of results from [quacc.schemas.cclib.cclib_summarize_run][].
+        See the type-hint for the data structure.
+    """
+    nprocs = psutil.cpu_count(logical=False) if nprocs == "max" else nprocs
+
+    default_inputs = [xc, basis, "normalprint", "numfreq" if numerical else "freq"]
+
+    default_blocks = [f"%pal nprocs {nprocs} end"]
+
+    return run_and_summarize(
+        atoms,
+        charge,
+        spin_multiplicity,
+        default_inputs=default_inputs,
+        default_blocks=default_blocks,
+        input_swaps=orcasimpleinput,
+        block_swaps=orcablocks,
+        additional_fields={"name": "ORCA vibrational frequency analysis"},
+        copy_files=copy_files,
+    )
+
+
 @job
 def relax_job(
     atoms: Atoms,
@@ -215,3 +285,84 @@ def ase_relax_job(
         additional_fields={"name": "ORCA ASE Relax"},
         copy_files=copy_files,
     )
+
+
+@job
+def ase_quasi_irc_perturb_job(
+    atoms: Atoms,
+    mode: list[list[float]] | NDArray,
+    charge: int = 0,
+    spin_multiplicity: int = 1,
+    perturb_magnitude: float = 0.6,
+    direction: Literal["forward", "reverse"] = "forward",
+    xc: str = "wb97x-d3bj",
+    basis: str = "def2-tzvp",
+    orcasimpleinput: list[str] | None = None,
+    orcablocks: list[str] | None = None,
+    opt_params: dict[str, Any] | None = None,
+    nprocs: int | Literal["max"] = "max",
+    copy_files: SourceDirectory | dict[SourceDirectory, Filenames] | None = None,
+) -> cclibSchema:
+    """
+    Quasi-IRC to optimize a reaction endpoint from a transition-state with known vibrational frequency modes.
+    Perturbs the structure of `atoms` by a finite amount (0.6 * the normalized mode magnitude) along the specified
+    vibrational frequency mode (assumed to be the transition mode), and then performs a `relax_job` on the perturbed
+    structure.
+
+    Parameters
+    ----------
+    atoms
+        Atoms object
+    mode
+        Transition mode. This should be an Nx3 matrix, where N is the number of atoms in `atoms`.
+    charge
+        Charge of the system.
+    spin_multiplicity
+        Multiplicity of the system.
+    perturb_magnitude
+        Factor to multiply the transition mode. Default is 0.6. In some cases, it may be advisable to increase this
+        factor, perhaps to 1.0 or 1.1. Lowering it is not generally found to be helpful.
+    direction
+        Direction of the (Quasi)IRC. Should be "forward" or "reverse".
+    xc
+        Exchange-correlation functional
+    basis
+        Basis set.
+    orcasimpleinput
+        List of `orcasimpleinput` swaps for the calculator. To remove entries
+        from the defaults, put a `#` in front of the name. Refer to the
+        [ase.calculators.orca.ORCA][] calculator for details on `orcasimpleinput`.
+    orcablocks
+        List of `orcablocks` swaps for the calculator. To remove entries
+        from the defaults, put a `#` in front of the name. Refer to the
+        [ase.calculators.orca.ORCA][] calculator for details on `orcablocks`.
+    nprocs
+        Number of processors to use. Defaults to the number of physical cores.
+    copy_files
+        Files to copy (and decompress) from source to the runtime directory.
+
+    Returns
+    -------
+    cclibASEOptSchema
+        Dictionary of results from [quacc.schemas.cclib.cclib_summarize_run][] merged with
+        the results from [quacc.schemas.ase.summarize_opt_run][].
+        See the type-hint for the data structure.
+    """
+    nprocs = psutil.cpu_count(logical=False) if nprocs == "max" else nprocs
+    default_inputs = [xc, basis, "engrad", "normalprint"]
+    default_blocks = [f"%pal nprocs {nprocs} end"]
+
+    scale = perturb_magnitude if direction == "forward" else perturb_magnitude * -1
+
+    return run_and_summarize_opt(
+        perturb(atoms, mode, scale),
+        charge=charge,
+        spin_multiplicity=spin_multiplicity,
+        default_inputs=default_inputs,
+        default_blocks=default_blocks,
+        input_swaps=orcasimpleinput,
+        block_swaps=orcablocks,
+        opt_params=opt_params,
+        additional_fields={"name": "ORCA ASE Quasi-IRC perturbed optimization"},
+        copy_files=copy_files,
+    )