Skip to content

Commit

Permalink
Add perturbation-based quasi-IRC workflows for ORCA and Q-Chem (#2042)
Browse files Browse the repository at this point in the history
## Summary of Changes

Title says most of it. I also added a frequency workflow for ORCA; felt
strange that there wasn't one already.

### Checklist

- [X] I have read the ["Guidelines"
section](https://quantum-accelerators.github.io/quacc/dev/contributing.html#guidelines)
of the contributing guide. Don't lie! 😉
- [X] My PR is on a custom branch and is _not_ named `main`.
- [X] I have added relevant, comprehensive [unit
tests](https://quantum-accelerators.github.io/quacc/dev/contributing.html#unit-tests).

### Notes

- Your PR will likely not be merged without proper and thorough tests.
- If you are an external contributor, you will see a comment from
[@buildbot-princeton](https://github.com/buildbot-princeton). This is
solely for the maintainers.
- When your code is ready for review, ping one of the [active
maintainers](https://quantum-accelerators.github.io/quacc/about/contributors.html#active-maintainers).

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Andrew S. Rosen <asrosen93@gmail.com>
  • Loading branch information
3 people committed Apr 23, 2024
1 parent 23ce801 commit dc41d80
Show file tree
Hide file tree
Showing 11 changed files with 551 additions and 19 deletions.
1 change: 1 addition & 0 deletions docs/about/contributors.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
30 changes: 17 additions & 13 deletions docs/user/recipes/recipes_list.md
Original file line number Diff line number Diff line change
Expand Up @@ -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>

Expand All @@ -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>

Expand Down
32 changes: 32 additions & 0 deletions src/quacc/atoms/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -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__)

Expand Down Expand Up @@ -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
151 changes: 151 additions & 0 deletions src/quacc/recipes/orca/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
)
Loading

0 comments on commit dc41d80

Please sign in to comment.