diff --git a/.github/workflows/tests.yaml b/.github/workflows/tests.yaml index a7d4433ce7..04cf7688e0 100644 --- a/.github/workflows/tests.yaml +++ b/.github/workflows/tests.yaml @@ -344,7 +344,7 @@ jobs: echo "COVERAGE_CORE=sysmon" >> $GITHUB_ENV - name: Run tests with pytest - run: pytest -k 'mlp or newtonnet' --durations=10 --cov=quacc --cov-report=xml + run: pytest -k 'mlp or newtonnet or geodesic' --durations=10 --cov=quacc --cov-report=xml - name: Upload code coverage report to Artifact uses: actions/upload-artifact@v4 diff --git a/docs/user/recipes/recipes_list.md b/docs/user/recipes/recipes_list.md index fe58ce0e73..28b25eb2cd 100644 --- a/docs/user/recipes/recipes_list.md +++ b/docs/user/recipes/recipes_list.md @@ -107,14 +107,16 @@ The list of available quacc recipes is shown below. The "Req'd Extras" column sp
-| Name | Decorator | Documentation | Req'd Extras | -| ------------------- | --------------- | -------------------------------------------- | -------------- | -| NewtonNet Static | `#!Python @job` | [quacc.recipes.newtonnet.core.static_job][] | | -| NewtonNet Relax | `#!Python @job` | [quacc.recipes.newtonnet.core.relax_job][] | | -| NewtonNet Frequency | `#!Python @job` | [quacc.recipes.newtonnet.core.freq_job][] | | -| NewtonNet TS | `#!Python @job` | [quacc.recipes.newtonnet.ts.ts_job][] | `quacc[sella]` | -| NewtonNet IRC | `#!Python @job` | [quacc.recipes.newtonnet.ts.irc_job][] | `quacc[sella]` | -| NewtonNet Quasi IRC | `#!Python @job` | [quacc.recipes.newtonnet.ts.quasi_irc_job][] | `quacc[sella]` | +| Name | Decorator | Documentation | Req'd Extras | +|---------------------| --------------- |-----------------------------------------------| | +| NewtonNet Static | `#!Python @job` | [quacc.recipes.newtonnet.core.static_job][] | | +| NewtonNet Relax | `#!Python @job` | [quacc.recipes.newtonnet.core.relax_job][] | | +| NewtonNet Frequency | `#!Python @job` | [quacc.recipes.newtonnet.core.freq_job][] | | +| NewtonNet TS | `#!Python @job` | [quacc.recipes.newtonnet.ts.ts_job][] | `quacc[sella]` | +| NewtonNet IRC | `#!Python @job` | [quacc.recipes.newtonnet.ts.irc_job][] | `quacc[sella]` | +| NewtonNet Quasi IRC | `#!Python @job` | [quacc.recipes.newtonnet.ts.quasi_irc_job][] | `quacc[sella]` | +| NewtonNet neb | `#!Python @job` | [quacc.recipes.newtonnet.ts.neb_job][] | | +| NewtonNet geodesic | `#!Python @job` | [quacc.recipes.newtonnet.ts.geodesic_job][] | |
diff --git a/src/quacc/atoms/ts.py b/src/quacc/atoms/ts.py new file mode 100644 index 0000000000..cd7329eb99 --- /dev/null +++ b/src/quacc/atoms/ts.py @@ -0,0 +1,115 @@ +from __future__ import annotations + +import logging +from importlib.util import find_spec +from typing import TYPE_CHECKING + +from ase.atoms import Atoms +from monty.dev import requires + +from quacc.atoms.core import copy_atoms + +LOGGER = logging.getLogger(__name__) + +has_geodesic_interpolate = bool(find_spec("geodesic_interpolate")) + +if has_geodesic_interpolate: + from geodesic_interpolate.geodesic import Geodesic + from geodesic_interpolate.interpolation import redistribute + +if TYPE_CHECKING: + from typing import Literal + + +@requires( + has_geodesic_interpolate, + "geodesic-interpolate must be installed. Refer to the quacc documentation.", +) +def geodesic_interpolate_wrapper( + reactant: Atoms, + product: Atoms, + n_images: int = 10, + perform_sweep: bool | Literal["auto"] = "auto", + redistribute_tol: float = 1e-2, + smoother_tol: float = 2e-3, + max_iterations: int = 15, + max_micro_iterations: int = 20, + morse_scaling: float = 1.7, + geometry_friction: float = 1e-2, + distance_cutoff: float = 3.0, + sweep_cutoff_size: int = 35, +) -> list[Atoms]: + """ + Interpolates between two geometries and optimizes the path with the geodesic method. + + Parameters + ---------- + reactant + The ASE Atoms object representing the initial geometry. + product + The ASE Atoms object representing the final geometry. + n_images + Number of images for interpolation. Default is 10. + perform_sweep + Whether to sweep across the path optimizing one image at a time. + Default is to perform sweeping updates if there are more than 35 atoms. + redistribute_tol + the value passed to the tol keyword argument of + geodesic_interpolate.interpolation.redistribute. Default is 1e-2. + smoother_tol + the value passed to the tol keyword argument of geodesic_smoother.smooth + or geodesic_smoother.sweep. Default is 2e-3. + max_iterations + Maximum number of minimization iterations. Default is 15. + max_micro_iterations + Maximum number of micro iterations for the sweeping algorithm. Default is 20. + morse_scaling + Exponential parameter for the Morse potential. Default is 1.7. + geometry_friction + Size of friction term used to prevent very large changes in geometry. Default is 1e-2. + distance_cutoff + Cut-off value for the distance between a pair of atoms to be included in the coordinate system. Default is 3.0. + sweep_cutoff_size + Cut off system size that above which sweep function will be called instead of smooth + in Geodesic. + + Returns + ------- + list[Atoms] + A list of ASE Atoms objects representing the smoothed path between the reactant and product geometries. + """ + reactant = copy_atoms(reactant) + product = copy_atoms(product) + + # Read the initial geometries. + chemical_symbols = reactant.get_chemical_symbols() + + # First redistribute number of images. Perform interpolation if too few and subsampling if too many images are given + raw_interpolated_positions = redistribute( + chemical_symbols, + [reactant.positions, product.positions], + n_images, + tol=redistribute_tol, + ) + + # Perform smoothing by minimizing distance in Cartesian coordinates with redundant internal metric + # to find the appropriate geodesic curve on the hyperspace. + geodesic_smoother = Geodesic( + chemical_symbols, + raw_interpolated_positions, + morse_scaling, + threshold=distance_cutoff, + friction=geometry_friction, + ) + if perform_sweep == "auto": + perform_sweep = len(chemical_symbols) > sweep_cutoff_size + if perform_sweep: + geodesic_smoother.sweep( + tol=smoother_tol, max_iter=max_iterations, micro_iter=max_micro_iterations + ) + else: + geodesic_smoother.smooth(tol=smoother_tol, max_iter=max_iterations) + return [ + Atoms(symbols=chemical_symbols, positions=geom) + for geom in geodesic_smoother.path + ] diff --git a/src/quacc/recipes/newtonnet/core.py b/src/quacc/recipes/newtonnet/core.py index ed63231318..929dcecf15 100644 --- a/src/quacc/recipes/newtonnet/core.py +++ b/src/quacc/recipes/newtonnet/core.py @@ -204,7 +204,7 @@ def freq_job( ) -def _add_stdev_and_hess(summary: dict[str, Any]) -> dict[str, Any]: +def _add_stdev_and_hess(summary: dict[str, Any], **calc_kwargs) -> dict[str, Any]: """ Calculate and add standard deviation values and Hessians to the summary. @@ -218,6 +218,10 @@ def _add_stdev_and_hess(summary: dict[str, Any]) -> dict[str, Any]: ---------- summary A dictionary containing information about the molecular trajectory. + **calc_kwargs + Custom kwargs for the NewtonNet calculator. Set a value to + `quacc.Remove` to remove a pre-existing key entirely. For a list of available + keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator. Returns ------- @@ -226,11 +230,13 @@ def _add_stdev_and_hess(summary: dict[str, Any]) -> dict[str, Any]: Hessian values. """ settings = get_settings() + calc_defaults = { + "model_path": settings.NEWTONNET_MODEL_PATH, + "settings_path": settings.NEWTONNET_CONFIG_PATH, + } + calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs) for i, atoms in enumerate(summary["trajectory"]): - calc = NewtonNet( - model_path=settings.NEWTONNET_MODEL_PATH, - settings_path=settings.NEWTONNET_CONFIG_PATH, - ) + calc = NewtonNet(**calc_flags) results = Runner(atoms, calc).run_calc().calc.results summary["trajectory_results"][i]["hessian"] = results["hessian"] summary["trajectory_results"][i]["energy_std"] = results["energy_disagreement"] diff --git a/src/quacc/recipes/newtonnet/ts.py b/src/quacc/recipes/newtonnet/ts.py index a99a5142a2..b379b4ea21 100644 --- a/src/quacc/recipes/newtonnet/ts.py +++ b/src/quacc/recipes/newtonnet/ts.py @@ -5,6 +5,8 @@ from importlib.util import find_spec from typing import TYPE_CHECKING +import numpy as np +from ase.mep import NEB from monty.dev import requires from quacc import change_settings, get_settings, job, strip_decorator @@ -13,6 +15,7 @@ from quacc.schemas.ase import Summarize from quacc.utils.dicts import recursive_dict_merge +has_geodesic_interpolate = bool(find_spec("geodesic_interpolate")) has_sella = bool(find_spec("sella")) has_newtonnet = bool(find_spec("newtonnet")) @@ -20,7 +23,8 @@ from sella import IRC, Sella if has_newtonnet: from newtonnet.utils.ase_interface import MLAseCalculator as NewtonNet - +if has_geodesic_interpolate: + from quacc.atoms.ts import geodesic_interpolate_wrapper if TYPE_CHECKING: from typing import Any, Literal @@ -29,6 +33,8 @@ from numpy.typing import NDArray from quacc.types import ( + GeodesicSchema, + NebSchema, NewtonNetIRCSchema, NewtonNetQuasiIRCSchema, NewtonNetTSSchema, @@ -85,7 +91,6 @@ def ts_job( calc_defaults = { "model_path": settings.NEWTONNET_MODEL_PATH, "settings_path": settings.NEWTONNET_CONFIG_PATH, - "hess_method": "autograd", } opt_defaults = { "optimizer": Sella, @@ -105,7 +110,7 @@ def ts_job( # Run the TS optimization dyn = Runner(atoms, calc).run_opt(**opt_flags) opt_ts_summary = _add_stdev_and_hess( - Summarize(additional_fields={"name": "NewtonNet TS"}).opt(dyn) + Summarize(additional_fields={"name": "NewtonNet TS"}).opt(dyn), **calc_flags ) # Run a frequency calculation @@ -267,7 +272,193 @@ def quasi_irc_job( return relax_summary -def _get_hessian(atoms: Atoms) -> NDArray: +@job +@requires( + has_newtonnet, "NewtonNet must be installed. Refer to the quacc documentation." +) +def neb_job( + reactant_atoms: Atoms, + product_atoms: Atoms, + interpolation_method: Literal["linear", "idpp", "geodesic"] = "linear", + relax_job_kwargs: dict[str, Any] | None = None, + interpolate_kwargs: dict[str, Any] | None = None, + neb_kwargs: dict[str, Any] | None = None, + **calc_kwargs, +) -> NebSchema: + """ + Perform a nudged elastic band (NEB) calculation to find the minimum energy path (MEP) between the given reactant and product structures. + + Parameters + ---------- + reactant_atoms + The Atoms object representing the reactant structure. + product_atoms + The Atoms object representing the product structure. + interpolation_method + The method to initialize the NEB optimization. There are three choices here, "linear", "idpp" and "geodesic". + Defaults to linear. + relax_job_kwargs + Keyword arguments to use for the [quacc.recipes.newtonnet.core.relax_job][] function. + interpolate_kwargs + Keyword arguments for the interpolate functions (geodesic, linear or idpp). + neb_kwargs + Keyword arguments for the NEB calculation. + **calc_kwargs + Dictionary of custom kwargs for the NewtonNet calculator. Set a value to + `quacc.Remove` to remove a pre-existing key entirely. For a list of available + keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator. + + Returns + ------- + NebSchema + A dictionary containing the NEB results + """ + relax_job_kwargs = relax_job_kwargs or {} + settings = get_settings() + + calc_defaults = { + "model_path": settings.NEWTONNET_MODEL_PATH, + "settings_path": settings.NEWTONNET_CONFIG_PATH, + "hess_method": None, + } + + interpolate_defaults = {"n_images": 20} + + neb_defaults = {"method": "aseneb", "precon": None} + calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs) + interpolate_flags = recursive_dict_merge(interpolate_defaults, interpolate_kwargs) + neb_flags = recursive_dict_merge(neb_defaults, neb_kwargs) + + # Define calculator + calc = NewtonNet(**calc_flags) + + # Run relax job + relax_summary_r = strip_decorator(relax_job)(reactant_atoms, **relax_job_kwargs) + relax_summary_p = strip_decorator(relax_job)(product_atoms, **relax_job_kwargs) + + if interpolation_method == "geodesic": + images = geodesic_interpolate_wrapper( + relax_summary_r["atoms"], relax_summary_p["atoms"], **interpolate_flags + ) + else: + images = [reactant_atoms] + images += [ + reactant_atoms.copy() for _ in range(interpolate_flags["n_images"] - 2) + ] + images += [product_atoms] + neb = NEB(images) + + # Interpolate linearly the positions of the middle images: + neb.interpolate(method=interpolation_method) + images = neb.images + + max_steps = neb_flags.pop("max_steps", None) + dyn = Runner(images, calc).run_neb(max_steps=max_steps, neb_kwargs=neb_flags) + + return { + "relax_reactant": relax_summary_r, + "relax_product": relax_summary_p, + "initial_images": images, + "neb_results": Summarize( + additional_fields={ + "neb_flags": neb_flags, + "calc_flags": calc_flags, + "interpolate_flags": interpolate_flags, + } + ).neb(dyn, len(images)), + } + + +@job +@requires( + has_newtonnet, "NewtonNet must be installed. Refer to the quacc documentation." +) +@requires( + has_geodesic_interpolate, + "geodesic-interpolate must be installed. Refer to the quacc documentation.", +) +def geodesic_job( + reactant_atoms: Atoms, + product_atoms: Atoms, + relax_job_kwargs: dict[str, Any] | None = None, + geodesic_interpolate_kwargs: dict[str, Any] | None = None, + **calc_kwargs, +) -> GeodesicSchema: + """ + Perform a quasi-IRC job using the given reactant and product atoms objects. + + Parameters + ---------- + reactant_atoms + The Atoms object representing the reactant structure. + product_atoms + The Atoms object representing the product structure. + relax_job_kwargs + Keyword arguments to use for [quacc.recipes.newtonnet.core.relax_job][] function. + geodesic_interpolate_kwargs + Keyword arguments for geodesic_interpolate, by default None. + **calc_kwargs + Dictionary of custom kwargs for the NewtonNet calculator. Set a value to + `quacc.Remove` to remove a pre-existing key entirely. For a list of available + keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator. + + Returns + ------- + GeodesicSchema + A dictionary containing the following keys: + - 'relax_reactant': Summary of the relaxed reactant structure. + - 'relax_product': Summary of the relaxed product structure. + - 'initial_images': The interpolated images between reactant and product. + - 'ts_atoms': ASE atoms object for the highest energy structure for the geodesic path + """ + relax_job_kwargs = relax_job_kwargs or {} + settings = get_settings() + + calc_defaults = { + "model_path": settings.NEWTONNET_MODEL_PATH, + "settings_path": settings.NEWTONNET_CONFIG_PATH, + "hess_method": None, + } + + geodesic_defaults = {"n_images": 20} + + calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs) + calc_flags["hess_method"] = None + geodesic_interpolate_flags = recursive_dict_merge( + geodesic_defaults, geodesic_interpolate_kwargs + ) + + # Define calculator + reactant_atoms.calc = NewtonNet(**calc_flags) + product_atoms.calc = NewtonNet(**calc_flags) + + # Run IRC + relax_summary_r = strip_decorator(relax_job)(reactant_atoms, **relax_job_kwargs) + relax_summary_p = strip_decorator(relax_job)(product_atoms, **relax_job_kwargs) + + images = geodesic_interpolate_wrapper( + relax_summary_r["atoms"].copy(), + relax_summary_p["atoms"].copy(), + **geodesic_interpolate_flags, + ) + + potential_energies = [] + for image in images: + image.calc = NewtonNet(**calc_flags) + potential_energies.append(image.get_potential_energy()) + + ts_index = np.argmax(potential_energies) + ts_atoms = images[ts_index] + + return { + "relax_reactant": relax_summary_r, + "relax_product": relax_summary_p, + "initial_images": images, + "ts_atoms": ts_atoms, + } + + +def _get_hessian(atoms: Atoms, **calc_kwargs) -> NDArray: """ Calculate and retrieve the Hessian matrix for the given molecular configuration. @@ -280,6 +471,10 @@ def _get_hessian(atoms: Atoms) -> NDArray: ---------- atoms The ASE Atoms object representing the molecular configuration. + **calc_kwargs + Dictionary of custom kwargs for the NewtonNet calculator. Set a value to + `quacc.Remove` to remove a pre-existing key entirely. For a list of available + keys, refer to the `newtonnet.utils.ase_interface.MLAseCalculator` calculator. Returns ------- @@ -287,10 +482,12 @@ def _get_hessian(atoms: Atoms) -> NDArray: The calculated Hessian matrix, reshaped into a 2D array. """ settings = get_settings() - ml_calculator = NewtonNet( - model_path=settings.NEWTONNET_MODEL_PATH, - settings_path=settings.NEWTONNET_CONFIG_PATH, - ) - ml_calculator.calculate(atoms) - - return ml_calculator.results["hessian"].reshape((-1, 3 * len(atoms))) + calc_defaults = { + "model_path": settings.NEWTONNET_MODEL_PATH, + "settings_path": settings.NEWTONNET_CONFIG_PATH, + "hess_method": "autograd", + } + calc_flags = recursive_dict_merge(calc_defaults, calc_kwargs) + calc = NewtonNet(**calc_flags) + calc.calculate(atoms) + return calc.results["hessian"].reshape((-1, 3 * len(atoms))) diff --git a/src/quacc/runners/ase.py b/src/quacc/runners/ase.py index 4157c87085..74be4a1ced 100644 --- a/src/quacc/runners/ase.py +++ b/src/quacc/runners/ase.py @@ -2,10 +2,13 @@ from __future__ import annotations +from collections.abc import Callable +from copy import deepcopy from importlib.util import find_spec from logging import getLogger +from pathlib import Path from shutil import copy, copytree -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Any import numpy as np from ase.calculators import calculator @@ -17,21 +20,22 @@ Stationary, ZeroRotation, ) -from ase.optimize import BFGS +from ase.mep import NEB +from ase.mep.neb import NEBOptimizer +from ase.optimize import BFGS, BFGSLineSearch from ase.optimize.sciopt import SciPyOptimizer from ase.vibrations import Vibrations from monty.dev import requires from monty.os.path import zpath -from quacc.atoms.core import copy_atoms from quacc.runners._base import BaseRunner -from quacc.runners.prep import terminate +from quacc.runners.prep import calc_cleanup, calc_setup, terminate from quacc.utils.dicts import recursive_dict_merge LOGGER = getLogger(__name__) has_sella = bool(find_spec("sella")) - +has_geodesic_interpolate = bool(find_spec("geodesic_interpolate")) if TYPE_CHECKING: from collections.abc import Callable @@ -40,7 +44,7 @@ from ase.atoms import Atoms from ase.calculators.calculator import Calculator - from ase.optimize.optimize import Dynamics + from ase.optimize.optimize import Dynamics, Optimizer from quacc.types import ( Filenames, @@ -58,7 +62,7 @@ class Runner(BaseRunner): def __init__( self, - atoms: Atoms, + atoms: Atoms | list[Atoms], calculator: Calculator, copy_files: SourceDirectory | dict[SourceDirectory, Filenames] | None = None, ) -> None: @@ -68,7 +72,7 @@ def __init__( Parameters ---------- atoms - The Atoms object to run calculations on. + The Atoms object to run calculations on. A list[Atoms] is used for NEB. calculator The instantiated ASE calculator object to attach to the Atoms object. copy_files @@ -78,10 +82,15 @@ def __init__( ------- None """ - self.atoms = copy_atoms(atoms) - self.atoms.calc = calculator self.copy_files = copy_files - self.setup() + if isinstance(atoms, list): + self.atoms = [image.copy() for image in atoms] + for image in self.atoms: + image.calc = deepcopy(calculator) + else: + self.atoms = atoms.copy() + self.atoms.calc = calculator + self.setup() def run_calc( self, geom_file: str | None = None, properties: list[str] | None = None @@ -347,6 +356,97 @@ def run_md( optimizer_kwargs=dynamics_kwargs, ) + def run_neb( + self, + relax_cell: bool = False, + fmax: float = 0.01, + max_steps: int | None = 1000, + optimizer: Optimizer = NEBOptimizer, + optimizer_kwargs: dict[str, Any] | None = None, + neb_kwargs: dict[str, Any] | None = None, + run_kwargs: dict[str, Any] | None = None, + ) -> Dynamics: + """ + Run an NEB calculation. + + Parameters + ---------- + relax_cell + Whether to relax the unit cell shape and volume. + fmax + Tolerance for the force convergence (in eV/A). + max_steps + Maximum number of steps to take. + optimizer + Optimizer class to use. + optimizer_kwargs + Dictionary of kwargs for the optimizer. Takes all valid kwargs for ASE + Optimizer classes. + neb_kwargs + Dictionary of kwargs for the NEB class. + run_kwargs + Dictionary of kwargs for the `run()` method of the optimizer. + + Returns + ------- + Dynamics + The ASE Dynamics object following an NEB calculation. + """ + images = self.atoms + run_kwargs = run_kwargs or {} + neb_kwargs = neb_kwargs or {} + traj_filename = "opt.traj" + + # Create a parent temporary directory for the NEB run + neb_tmpdir, neb_results_dir = calc_setup(None) + + # Adjust optimizer_kwargs to use the parent directory + optimizer_kwargs = recursive_dict_merge( + { + "logfile": str(neb_tmpdir / "opt.log"), + "restart": str(neb_tmpdir / "opt.json"), + }, + optimizer_kwargs, + ) + + if "trajectory" in optimizer_kwargs: + msg = "Quacc does not support setting the `trajectory` kwarg." + raise ValueError(msg) + + if optimizer == BFGSLineSearch: + raise ValueError("BFGSLineSearch is not allowed as optimizer with NEB.") + + # Copy atoms so we don't modify it in-place + neb = NEB(images, **neb_kwargs) + + # Perform staging operations + for i, image in enumerate(images): + image_tmpdir = neb_tmpdir / f"image_{i}" + image_tmpdir.mkdir() + image.calc.directory = image_tmpdir + + # Define the Trajectory object + traj_file = neb_tmpdir / traj_filename + traj = Trajectory(traj_file, "w", atoms=neb) + + # Set volume relaxation constraints, if relevant + if relax_cell: + for i in range(len(images)): + if images[i].pbc.any(): + images[i] = FrechetCellFilter(images[i]) + + dyn = optimizer(neb, **optimizer_kwargs) + dyn.attach(traj.write) + dyn.run(fmax, max_steps) + traj.close() + dyn.logfile.close() + + calc_cleanup(None, neb_tmpdir, neb_results_dir) + traj.filename = zpath(str(neb_results_dir / traj_filename)) + dyn.trajectory = traj + + return dyn + def _copy_intermediate_files( self, step_number: int, files_to_ignore: list[Path] | None = None ) -> None: diff --git a/src/quacc/schemas/ase.py b/src/quacc/schemas/ase.py index 084fdad8e8..68103ac99d 100644 --- a/src/quacc/schemas/ase.py +++ b/src/quacc/schemas/ase.py @@ -23,6 +23,7 @@ from typing import Any, Literal from ase.atoms import Atoms + from ase.io.trajectory import TrajectoryWriter from ase.md.md import MolecularDynamics from ase.optimize.optimize import Optimizer from ase.vibrations import Vibrations @@ -293,6 +294,89 @@ def md( store=store, ) + def neb( + self, + dyn: Optimizer, + n_images: int, + n_iter_return: int = -1, + trajectory: TrajectoryWriter | list[Atoms] | None = None, + store: Store | None | DefaultSetting = QuaccDefault, + ) -> OptSchema: + """ + Summarize the NEB run results and store them in a database-friendly format. + + Parameters + ---------- + dyn + ASE Optimizer object used for the NEB run. + n_images + Number of images in the NEB run. + n_iter_return + Number of iterations to return. If -1, all iterations are returned. + trajectory + Trajectory of the NEB run, either as a Trajectory object or a list of Atoms objects. + store + Maggma Store object to store the results in. Defaults to `QuaccSettings.STORE`. + + Returns + ------- + OptSchema + A dictionary containing the summarized NEB run results. + """ + store = self._settings.STORE if store == QuaccDefault else store + + # Get trajectory + if trajectory: + atoms_trajectory = trajectory + else: + atoms_trajectory = read(dyn.trajectory.filename, index=":") # type: ignore[union-attr] + + if n_iter_return == -1: + atoms_trajectory = atoms_trajectory[-(n_images):] + else: + atoms_trajectory = _get_nth_iteration( + atoms_trajectory, + int(len(atoms_trajectory) / n_images), + n_images, + n_iter_return, + ) + trajectory_results = [atoms.calc.results for atoms in atoms_trajectory] + ts_index = ( + np.argmax( + [ + result["energy"] + for result in trajectory_results[-(n_images - 1) : -1] + ] + ) + + 1 + ) + ts_atoms = atoms_trajectory[ts_index] + base_task_doc = atoms_to_metadata( + atoms_trajectory[0], charge_and_multiplicity=self.charge_and_multiplicity + ) + + # Clean up the opt parameters + parameters_opt = dyn.todict() + parameters_opt.pop("logfile", None) + parameters_opt.pop("restart", None) + + opt_fields = { + "parameters_opt": parameters_opt, + "trajectory": atoms_trajectory, + "trajectory_results": trajectory_results, + "ts_atoms": ts_atoms, + } + + # Create a dictionary of the inputs/outputs + unsorted_task_doc = base_task_doc | opt_fields | self.additional_fields + + return finalize_dict( + unsorted_task_doc, + directory=None, + gzip_file=self._settings.GZIP_FILES, + store=store, + ) + class VibSummarize: """ @@ -516,3 +600,37 @@ def vib_and_thermo( gzip_file=self._settings.GZIP_FILES, store=store, ) + + +def _get_nth_iteration( + neb_trajectory: list[Atoms], n_iter: int, n_images: int, interval: int +) -> list[Atoms]: + """ + Extract every nth iteration from the NEB trajectory. + + Parameters + ---------- + neb_trajectory + List of configurations (length: n_iter * n_images). + n_iter + Total number of iterations. + n_images + Number of images per iteration. + interval + Interval to get every nth iteration. + + Returns + ------- + list[Atoms] + List of configurations from every nth iteration. + """ + result = [] + start_idx, end_idx = 0, 0 + for i in range(0, n_iter, interval): + start_idx = i * n_images + end_idx = start_idx + n_images + + result.extend(neb_trajectory[start_idx:end_idx]) + if end_idx < len(neb_trajectory) - 1: + result.extend(neb_trajectory[-(n_images):]) + return result diff --git a/src/quacc/types.py b/src/quacc/types.py index f9db8555ef..56de12a0a7 100644 --- a/src/quacc/types.py +++ b/src/quacc/types.py @@ -685,6 +685,31 @@ class NewtonNetQuasiIRCSchema(OptSchema): irc_job: NewtonNetIRCSchema freq_job: VibThermoSchema | None + class NebSchema(TypedDict): + relax_reactant: OptSchema + relax_product: OptSchema + initial_images: list[Atoms] + neb_results: dict + + class NebTsSchema(TypedDict): + relax_reactant: OptSchema + relax_product: OptSchema + initial_images: list[Atoms] + neb_results: dict + ts_results: OptSchema + + class GeodesicSchema(TypedDict): + relax_reactant: OptSchema + relax_product: OptSchema + initial_images: list[Atoms] + ts_atoms: Atoms + + class GeodesicTsSchema(TypedDict): + relax_reactant: OptSchema + relax_product: OptSchema + initial_images: list[Atoms] + ts_results: NewtonNetTSSchema + # ----------- Recipe (Q-Chem) type hints ----------- class QchemQuasiIRCSchema(OptSchema): diff --git a/tests/core/atoms/test_ts.py b/tests/core/atoms/test_ts.py new file mode 100644 index 0000000000..6122188286 --- /dev/null +++ b/tests/core/atoms/test_ts.py @@ -0,0 +1,83 @@ +from __future__ import annotations + +import numpy as np +import pytest + + +@pytest.fixture(scope="module", autouse=True) +def set_seed(): + np.random.seed(42) # noqa: NPY002 + + +from importlib.util import find_spec + +from ase.atoms import Atoms +from ase.build import molecule + +from quacc.atoms.ts import geodesic_interpolate_wrapper + +has_geodesic_interpolate = bool(find_spec("geodesic_interpolate")) + + +@pytest.fixture +def setup_test_environment(tmp_path): + reactant = Atoms( + symbols="CCHHCHH", + positions=[ + [1.4835950817281542, -1.0145410211301968, -0.13209027203235943], + [0.8409564131524673, 0.018549610257914483, -0.07338809662321308], + [-0.6399757891931867, 0.01763740851518944, 0.0581573443268891], + [-1.0005576455546672, 1.0430257532387608, 0.22197240310602892], + [1.402180736662139, 0.944112416574632, -0.12179540364365492], + [-1.1216961389434357, -0.3883639833876232, -0.8769102842015071], + [-0.9645026578514683, -0.6204201840686793, 0.9240543090678239], + ], + ) + + product = Atoms( + symbols="CCHHCHH", + positions=[ + [1.348003553501624, 0.4819311116778978, 0.2752537177143993], + [0.2386618286631742, -0.3433222966734429, 0.37705518940917926], + [-0.9741307940518336, 0.07686022294949588, 0.08710778043683955], + [-1.8314843503320921, -0.5547344604780035, 0.1639037492534953], + [0.3801391040059668, -1.3793340533058087, 0.71035902765307], + [1.9296265384257907, 0.622088341468767, 1.0901733942191298], + [-1.090815880212625, 1.0965111343610956, -0.23791518420660265], + ], + ) + return reactant, product + + +@pytest.mark.skipif( + not has_geodesic_interpolate, + reason="geodesic_interpolate function is not available", +) +def test_geodesic_interpolate_wrapper(setup_test_environment): + reactant, product = setup_test_environment + + # Execute the geodesic_interpolate_wrapper function + smoother_path = geodesic_interpolate_wrapper( + reactant, + product, + n_images=20, + redistribute_tol=5e-4, + smoother_tol=1e-4, + max_iterations=10, + max_micro_iterations=10, + morse_scaling=1.5, + geometry_friction=1e-2, + distance_cutoff=2.5, + ) + assert smoother_path[1].positions[0][0] == pytest.approx(1.378384900, abs=1e-5) + assert smoother_path[5].positions[0][2] == pytest.approx(-0.512075394, abs=1e-5) + + +@pytest.mark.skipif( + not has_geodesic_interpolate, + reason="geodesic_interpolate function is not available", +) +def test_geodesic_interpolate_wrapper_large_system(setup_test_environment): + # Test with large system to trigger sweeping updates + smoother_path = geodesic_interpolate_wrapper(molecule("C60"), molecule("C60")) + assert len(smoother_path) == 10 diff --git a/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py b/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py index ebb0cb6e5e..00bd3a70f3 100644 --- a/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py +++ b/tests/core/recipes/newtonnet_recipes/test_newtonnet_recipes.py @@ -2,16 +2,29 @@ import pytest + +@pytest.fixture(scope="module", autouse=True) +def set_seed(): + np.random.seed(42) # noqa: NPY002 + + pytest.importorskip("sella") pytest.importorskip("newtonnet") from pathlib import Path import numpy as np +from ase import Atoms from ase.build import molecule from quacc import _internally_set_settings from quacc.recipes.newtonnet.core import freq_job, relax_job, static_job -from quacc.recipes.newtonnet.ts import irc_job, quasi_irc_job, ts_job +from quacc.recipes.newtonnet.ts import ( + geodesic_job, + irc_job, + neb_job, + quasi_irc_job, + ts_job, +) current_file_path = Path(__file__).parent @@ -106,7 +119,7 @@ def test_ts_job_with_default_args(tmp_path, monkeypatch): assert "freq_job" in output assert output["results"]["energy"] == pytest.approx(-6.796914263061945) assert output["freq_job"]["results"]["imag_vib_freqs"][0] == pytest.approx( - -2426.7398321816004 + -2426.73983218162, abs=1e-6 ) @@ -282,3 +295,103 @@ def test_quasi_irc_job_with_custom_irc_swaps(tmp_path, monkeypatch): assert output["irc_job"]["results"]["energy"] == pytest.approx(-9.517354965639784) assert output["results"]["energy"] == pytest.approx(-9.517354965639784) assert output["freq_job"]["results"]["energy"] == pytest.approx(-9.517354965639784) + + +@pytest.fixture +def setup_test_environment(): + reactant = Atoms( + symbols="CCHHCHH", + positions=[ + [1.4835950817281542, -1.0145410211301968, -0.13209027203235943], + [0.8409564131524673, 0.018549610257914483, -0.07338809662321308], + [-0.6399757891931867, 0.01763740851518944, 0.0581573443268891], + [-1.0005576455546672, 1.0430257532387608, 0.22197240310602892], + [1.402180736662139, 0.944112416574632, -0.12179540364365492], + [-1.1216961389434357, -0.3883639833876232, -0.8769102842015071], + [-0.9645026578514683, -0.6204201840686793, 0.9240543090678239], + ], + ) + + product = Atoms( + symbols="CCHHCHH", + positions=[ + [1.348003553501624, 0.4819311116778978, 0.2752537177143993], + [0.2386618286631742, -0.3433222966734429, 0.37705518940917926], + [-0.9741307940518336, 0.07686022294949588, 0.08710778043683955], + [-1.8314843503320921, -0.5547344604780035, 0.1639037492534953], + [0.3801391040059668, -1.3793340533058087, 0.71035902765307], + [1.9296265384257907, 0.622088341468767, 1.0901733942191298], + [-1.090815880212625, 1.0965111343610956, -0.23791518420660265], + ], + ) + + return reactant, product + + +def test_neb_job_linear(setup_test_environment): + reactant, product = setup_test_environment + + neb_summary = neb_job( + reactant, product, interpolation_method="linear", neb_kwargs={"max_steps": 10} + ) + + assert len(neb_summary["neb_results"]["trajectory_results"]) == 20 + assert neb_summary["relax_reactant"]["atoms"].positions[0, 0] == pytest.approx( + 0.8815, abs=1e-3 + ) + assert neb_summary["relax_product"]["atoms"].positions[0, 0] == pytest.approx( + 1.117689, abs=1e-3 + ) + + assert neb_summary["neb_results"]["trajectory_results"][1][ + "energy" + ] == pytest.approx(-21.99570976499786, abs=1e-5) + + +def test_neb_job_idpp(setup_test_environment): + reactant, product = setup_test_environment + + neb_summary = neb_job( + reactant, product, interpolation_method="idpp", neb_kwargs={"max_steps": 10} + ) + + assert len(neb_summary["neb_results"]["trajectory_results"]) == 20 + assert neb_summary["relax_reactant"]["atoms"].positions[0, 0] == pytest.approx( + 0.8815, abs=1e-3 + ) + assert neb_summary["relax_product"]["atoms"].positions[0, 0] == pytest.approx( + 1.117689, abs=1e-3 + ) + + assert neb_summary["neb_results"]["trajectory_results"][1][ + "energy" + ] == pytest.approx(-23.58905916097854, abs=1e-5) + + +def test_neb_job_geodesic(setup_test_environment): + reactant, product = setup_test_environment + + neb_summary = neb_job( + reactant, product, interpolation_method="geodesic", neb_kwargs={"max_steps": 10} + ) + + assert len(neb_summary["neb_results"]["trajectory_results"]) == 20 + assert neb_summary["relax_reactant"]["atoms"].positions[0, 0] == pytest.approx( + 0.8815, abs=1e-3 + ) + assert neb_summary["relax_product"]["atoms"].positions[0, 0] == pytest.approx( + 1.117689, abs=1e-3 + ) + + assert neb_summary["neb_results"]["trajectory_results"][1][ + "energy" + ] == pytest.approx(-24.895280838012695, abs=0.05) + + +def test_geodesic_job(setup_test_environment): + reactant, product = setup_test_environment + + geodesic_summary = geodesic_job(reactant, product) + assert geodesic_summary["ts_atoms"].get_potential_energy() == pytest.approx( + -22.597125398584318, abs=1e-6 + ) diff --git a/tests/core/runners/test_ase.py b/tests/core/runners/test_ase.py index 3901cb0058..cc38617ad9 100644 --- a/tests/core/runners/test_ase.py +++ b/tests/core/runners/test_ase.py @@ -2,16 +2,26 @@ import glob import os +from importlib.util import find_spec from logging import WARNING, getLogger from pathlib import Path from shutil import rmtree import numpy as np import pytest + + +@pytest.fixture(scope="module", autouse=True) +def set_seed(): + np.random.seed(42) # noqa: NPY002 + + +from ase import Atoms from ase.build import bulk, molecule from ase.calculators.emt import EMT from ase.calculators.lj import LennardJones from ase.io import read +from ase.mep.neb import NEBOptimizer from ase.optimize import BFGS, BFGSLineSearch from ase.optimize.sciopt import SciPyFminBFGS @@ -19,6 +29,9 @@ from quacc.runners._base import BaseRunner from quacc.runners.ase import Runner +has_geodesic_interpolate = bool(find_spec("geodesic_interpolate")) +test_files_path = Path(__file__).parent / "test_files" + LOGGER = getLogger(__name__) LOGGER.propagate = True @@ -48,6 +61,36 @@ def teardown_function(): os.remove(os.path.join(results_dir, f)) +@pytest.fixture +def setup_test_environment(tmp_path): + reactant = Atoms( + symbols="CCHHCHH", + positions=[ + [1.4835950817281542, -1.0145410211301968, -0.13209027203235943], + [0.8409564131524673, 0.018549610257914483, -0.07338809662321308], + [-0.6399757891931867, 0.01763740851518944, 0.0581573443268891], + [-1.0005576455546672, 1.0430257532387608, 0.22197240310602892], + [1.402180736662139, 0.944112416574632, -0.12179540364365492], + [-1.1216961389434357, -0.3883639833876232, -0.8769102842015071], + [-0.9645026578514683, -0.6204201840686793, 0.9240543090678239], + ], + ) + + product = Atoms( + symbols="CCHHCHH", + positions=[ + [1.348003553501624, 0.4819311116778978, 0.2752537177143993], + [0.2386618286631742, -0.3433222966734429, 0.37705518940917926], + [-0.9741307940518336, 0.07686022294949588, 0.08710778043683955], + [-1.8314843503320921, -0.5547344604780035, 0.1639037492534953], + [0.3801391040059668, -1.3793340533058087, 0.71035902765307], + [1.9296265384257907, 0.622088341468767, 1.0901733942191298], + [-1.090815880212625, 1.0965111343610956, -0.23791518420660265], + ], + ) + return reactant, product + + def test_base_runner(tmp_path, monkeypatch): monkeypatch.chdir(tmp_path) atoms = bulk("Cu") @@ -249,3 +292,40 @@ def fn_hook(dyn): Runner(bulk("Cu"), EMT()).run_opt(fn_hook=fn_hook) with pytest.raises(ValueError, match="Test error"): raise err.value.parent_error + + +def test_run_neb(monkeypatch, tmp_path): + monkeypatch.chdir(tmp_path) + geodesic_path = test_files_path / "geodesic_path.xyz" + images = read(geodesic_path, index=":") + + neb_kwargs = {"method": "aseneb", "precon": None} + dyn = Runner(images, EMT()).run_neb(optimizer=NEBOptimizer, neb_kwargs=neb_kwargs) + traj = read(dyn.trajectory.filename, index=":") + + assert traj[-1].calc.results is not None + assert not os.path.exists(tmp_path / "opt.log") + + +def test_run_neb2(): + geodesic_path = test_files_path / "geodesic_path.xyz" + + images = read(geodesic_path, index=":") + + with pytest.raises( + ValueError, match="BFGSLineSearch is not allowed as optimizer with NEB." + ): + Runner(images, EMT()).run_neb( + optimizer=BFGSLineSearch, neb_kwargs={"method": "aseneb", "precon": None} + ) + + +def test_run_neb_raises_value_error_for_trajectory_kwarg(): + images = [Atoms("H2", positions=[[0, 0, 0], [0, 0, 0.74]])] + + with pytest.raises( + ValueError, match="Quacc does not support setting the `trajectory` kwarg." + ): + Runner(images, EMT()).run_neb( + optimizer=NEBOptimizer, optimizer_kwargs={"trajectory": "some_traj.traj"} + ) diff --git a/tests/core/runners/test_files/geodesic_path.xyz b/tests/core/runners/test_files/geodesic_path.xyz new file mode 100644 index 0000000000..73c288b5b2 --- /dev/null +++ b/tests/core/runners/test_files/geodesic_path.xyz @@ -0,0 +1,90 @@ +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.60210258 -0.85402158 0.29381286 +C -0.28783951 -0.53322413 -0.44372747 +H -0.95727408 0.42884906 0.07350937 +H 0.35481921 0.65326383 0.53310253 +C 0.83497720 -0.13893815 -0.61903514 +H -0.06453974 0.87896748 -0.88297197 +H -0.48224567 -0.43489651 1.04530982 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.59223789 -0.81802021 0.31135853 +C -0.28604018 -0.47152900 -0.42593957 +H -0.92390973 0.45417612 0.20964610 +H 0.38052911 0.69342674 0.54238992 +C 0.85025969 -0.11855947 -0.60382787 +H -0.05043198 0.84149609 -1.07238207 +H -0.56264480 -0.58099027 1.03875497 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.56752580 -0.76523536 0.32900863 +C -0.29444225 -0.36150756 -0.39642003 +H -0.88419505 0.46889917 0.41683865 +H 0.42167940 0.75256538 0.57420807 +C 0.85881363 -0.07344971 -0.57888790 +H -0.01427599 0.74054417 -1.33086389 +H -0.65510555 -0.76181608 0.98611646 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.53153649 -0.70795790 0.32135314 +C -0.31196159 -0.20232825 -0.36082476 +H -0.83211712 0.45957469 0.65330147 +H 0.46707007 0.80230819 0.65159567 +C 0.85547097 0.01215643 -0.55082250 +H 0.02479217 0.53705274 -1.58192589 +H -0.73479098 -0.90080590 0.86732288 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.49710666 -0.64719683 0.28652142 +C -0.32455217 -0.02770666 -0.32588575 +H -0.76903102 0.41846514 0.85063704 +H 0.50991033 0.83413692 0.73959991 +C 0.84936453 0.11862994 -0.53249860 +H 0.03875338 0.26781854 -1.71408411 +H -0.80155171 -0.96414705 0.69571009 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.46929389 -0.58407604 0.23873741 +C -0.32819281 0.13632023 -0.29577230 +H -0.70406684 0.34361913 0.98764228 +H 0.54736483 0.85455598 0.81084680 +C 0.84529882 0.22536273 -0.52705498 +H 0.02725281 -0.01111080 -1.71369497 +H -0.85695070 -0.96467122 0.49929576 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.45166801 -0.52656333 0.20165987 +C -0.32191784 0.26695481 -0.27071285 +H -0.64751808 0.25267214 1.06356077 +H 0.57434755 0.87413911 0.85281977 +C 0.84453688 0.30577859 -0.53229994 +H -0.00510047 -0.24219076 -1.62060929 +H -0.89601604 -0.93079055 0.30558168 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.44176149 -0.48136401 0.18034320 +C -0.31132747 0.36096930 -0.25222570 +H -0.60331694 0.16140671 1.09768002 +H 0.58999854 0.89461629 0.87334640 +C 0.84542407 0.35773560 -0.54324267 +H -0.04447183 -0.41002497 -1.49647966 +H -0.91806787 -0.88333892 0.14057842 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.43559245 -0.44795824 0.16566337 +C -0.30056794 0.42882307 -0.24077310 +H -0.56621759 0.07575027 1.10766820 +H 0.59486995 0.91244138 0.88346800 +C 0.84607809 0.39270844 -0.55602468 +H -0.07914315 -0.53074541 -1.37277954 +H -0.93061181 -0.83101951 0.01277776 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.43009268 -0.42437008 0.15071523 +C -0.29152973 0.47863977 -0.23705145 +H -0.53028723 -0.00274016 1.10374152 +H 0.58906640 0.92480496 0.88885267 +C 0.84668896 0.41932394 -0.56813409 +H -0.10017589 -0.61881108 -1.26384123 +H -0.94385519 -0.77684734 -0.07428264 diff --git a/tests/core/schemas/test_ase.py b/tests/core/schemas/test_ase.py index f32a6cd31c..019b20ae38 100644 --- a/tests/core/schemas/test_ase.py +++ b/tests/core/schemas/test_ase.py @@ -8,6 +8,7 @@ from ase.build import bulk, molecule from ase.calculators.emt import EMT from ase.io import read +from ase.mep import NEB from ase.optimize import BFGS from ase.vibrations import Vibrations from maggma.stores import MemoryStore @@ -429,3 +430,26 @@ def test_errors(tmp_path, monkeypatch): vib.run() with pytest.raises(ValueError, match="Unsupported thermo_method"): VibSummarize(vib, directory=tmp_path).vib_and_thermo("bad") + + +def test_summarize_neb(monkeypatch, tmp_path): + monkeypatch.chdir(tmp_path) + + # Read initial and final states: + images = read(FILE_DIR / "test_files" / "geodesic_path.xyz", index=":") + neb = NEB(images) + for image in images: + image.calc = EMT() + optimizer = BFGS(neb, trajectory="opt.traj") + optimizer.run(fmax=0.05) + + neb_summary = Summarize().neb( + optimizer, len(images), n_iter_return=10, trajectory=read("opt.traj", index=":") + ) + assert neb_summary["trajectory_results"][-2]["energy"] == pytest.approx( + 1.0875330074934446, abs=1e-4 + ) + + ts_atoms = neb_summary["ts_atoms"] + ts_atoms.calc = EMT() + assert ts_atoms.get_potential_energy() == pytest.approx(1.1603536513693768, 1e-4) diff --git a/tests/core/schemas/test_files/geodesic_path.xyz b/tests/core/schemas/test_files/geodesic_path.xyz new file mode 100644 index 0000000000..73c288b5b2 --- /dev/null +++ b/tests/core/schemas/test_files/geodesic_path.xyz @@ -0,0 +1,90 @@ +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.60210258 -0.85402158 0.29381286 +C -0.28783951 -0.53322413 -0.44372747 +H -0.95727408 0.42884906 0.07350937 +H 0.35481921 0.65326383 0.53310253 +C 0.83497720 -0.13893815 -0.61903514 +H -0.06453974 0.87896748 -0.88297197 +H -0.48224567 -0.43489651 1.04530982 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.59223789 -0.81802021 0.31135853 +C -0.28604018 -0.47152900 -0.42593957 +H -0.92390973 0.45417612 0.20964610 +H 0.38052911 0.69342674 0.54238992 +C 0.85025969 -0.11855947 -0.60382787 +H -0.05043198 0.84149609 -1.07238207 +H -0.56264480 -0.58099027 1.03875497 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.56752580 -0.76523536 0.32900863 +C -0.29444225 -0.36150756 -0.39642003 +H -0.88419505 0.46889917 0.41683865 +H 0.42167940 0.75256538 0.57420807 +C 0.85881363 -0.07344971 -0.57888790 +H -0.01427599 0.74054417 -1.33086389 +H -0.65510555 -0.76181608 0.98611646 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.53153649 -0.70795790 0.32135314 +C -0.31196159 -0.20232825 -0.36082476 +H -0.83211712 0.45957469 0.65330147 +H 0.46707007 0.80230819 0.65159567 +C 0.85547097 0.01215643 -0.55082250 +H 0.02479217 0.53705274 -1.58192589 +H -0.73479098 -0.90080590 0.86732288 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.49710666 -0.64719683 0.28652142 +C -0.32455217 -0.02770666 -0.32588575 +H -0.76903102 0.41846514 0.85063704 +H 0.50991033 0.83413692 0.73959991 +C 0.84936453 0.11862994 -0.53249860 +H 0.03875338 0.26781854 -1.71408411 +H -0.80155171 -0.96414705 0.69571009 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.46929389 -0.58407604 0.23873741 +C -0.32819281 0.13632023 -0.29577230 +H -0.70406684 0.34361913 0.98764228 +H 0.54736483 0.85455598 0.81084680 +C 0.84529882 0.22536273 -0.52705498 +H 0.02725281 -0.01111080 -1.71369497 +H -0.85695070 -0.96467122 0.49929576 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.45166801 -0.52656333 0.20165987 +C -0.32191784 0.26695481 -0.27071285 +H -0.64751808 0.25267214 1.06356077 +H 0.57434755 0.87413911 0.85281977 +C 0.84453688 0.30577859 -0.53229994 +H -0.00510047 -0.24219076 -1.62060929 +H -0.89601604 -0.93079055 0.30558168 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.44176149 -0.48136401 0.18034320 +C -0.31132747 0.36096930 -0.25222570 +H -0.60331694 0.16140671 1.09768002 +H 0.58999854 0.89461629 0.87334640 +C 0.84542407 0.35773560 -0.54324267 +H -0.04447183 -0.41002497 -1.49647966 +H -0.91806787 -0.88333892 0.14057842 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.43559245 -0.44795824 0.16566337 +C -0.30056794 0.42882307 -0.24077310 +H -0.56621759 0.07575027 1.10766820 +H 0.59486995 0.91244138 0.88346800 +C 0.84607809 0.39270844 -0.55602468 +H -0.07914315 -0.53074541 -1.37277954 +H -0.93061181 -0.83101951 0.01277776 +7 +Properties=species:S:1:pos:R:3 pbc="F F F" +C 0.43009268 -0.42437008 0.15071523 +C -0.29152973 0.47863977 -0.23705145 +H -0.53028723 -0.00274016 1.10374152 +H 0.58906640 0.92480496 0.88885267 +C 0.84668896 0.41932394 -0.56813409 +H -0.10017589 -0.61881108 -1.26384123 +H -0.94385519 -0.77684734 -0.07428264 diff --git a/tests/requirements-newtonnet.txt b/tests/requirements-newtonnet.txt index 932e4e19dd..5603fb55f6 100644 --- a/tests/requirements-newtonnet.txt +++ b/tests/requirements-newtonnet.txt @@ -1 +1,3 @@ newtonnet==1.1.1 +geodesic-interpolate @ git+https://github.com/virtualzx-nad/geodesic-interpolate.git +sella==2.3.5