diff --git a/.gitignore b/.gitignore index 6ccf65a08..3e370f368 100755 --- a/.gitignore +++ b/.gitignore @@ -11,7 +11,7 @@ autode\.egg-info/ doc/_build/ **.DS_Store **htmlcov/ -.coverage +.coverage* coverage.xml tests/*.xyz **/tmp.xyz @@ -20,3 +20,4 @@ tests/*.xyz *.inp *.pyc .autode_calculations +tests/initial_path/* diff --git a/autode/__init__.py b/autode/__init__.py index a71a3171d..15865ec2d 100644 --- a/autode/__init__.py +++ b/autode/__init__.py @@ -5,6 +5,7 @@ from autode import neb from autode import mol_graphs from autode import hessians +from autode.neb import NEB, CINEB from autode.reactions.reaction import Reaction from autode.reactions.multistep import MultiStepReaction from autode.transition_states.transition_state import TransitionState @@ -61,6 +62,8 @@ "NCIComplex", "Config", "Calculation", + "NEB", + "CINEB", "pes", "neb", "geom", diff --git a/autode/calculations/calculation.py b/autode/calculations/calculation.py index 9318f0e25..854e9d4b7 100644 --- a/autode/calculations/calculation.py +++ b/autode/calculations/calculation.py @@ -3,12 +3,8 @@ from copy import deepcopy from typing import Optional, List -from autode.utils import deprecated -from autode.atoms import Atoms from autode.point_charges import PointCharge from autode.log import logger -from autode.hessians import Hessian -from autode.values import PotentialEnergy, Gradient from autode.calculations.types import CalculationType from autode.calculations.executors import ( CalculationExecutor, @@ -305,47 +301,6 @@ def _check_properties_exist(self) -> None: return None - # ----------------- Deprecated functions ---------------------- - - @deprecated - def get_energy(self) -> PotentialEnergy: - self._executor.set_properties() - return self.molecule.energy - - @deprecated - def optimisation_converged(self) -> bool: - self._executor.set_properties() - return self.optimiser.converged - - @deprecated - def optimisation_nearly_converged(self) -> bool: - self._executor.set_properties() - tol = PotentialEnergy(0.1, units="kcal mol-1") - return ( - not self.optimiser.converged - and self.optimiser.last_energy_change < tol - ) - - @deprecated - def get_final_atoms(self) -> Atoms: - self._executor.set_properties() - return self.molecule.atoms - - @deprecated - def get_atomic_charges(self) -> List[float]: - self._executor.set_properties() - return self.molecule.partial_charges - - @deprecated - def get_gradients(self) -> Gradient: - self._executor.set_properties() - return self.molecule.gradient - - @deprecated - def get_hessian(self) -> Hessian: - self._executor.set_properties() - return self.molecule.hessian - def _are_opt(keywords) -> bool: return isinstance(keywords, kws.OptKeywords) diff --git a/autode/constraints.py b/autode/constraints.py index 6d9290603..77cda20ef 100644 --- a/autode/constraints.py +++ b/autode/constraints.py @@ -1,5 +1,7 @@ from collections.abc import MutableMapping from typing import Optional, Dict, List +from copy import deepcopy + from autode.values import Distance from autode.log import logger @@ -122,6 +124,9 @@ def update( return None + def copy(self) -> "Constraints": + return deepcopy(self) + class DistanceConstraints(MutableMapping): def __init__(self, *args, **kwargs): @@ -178,3 +183,6 @@ def __setitem__(self, key, value): ) self._store[self._key_transform(key)] = Distance(value) + + def copy(self) -> "DistanceConstraints": + return deepcopy(self) diff --git a/autode/neb/ci.py b/autode/neb/ci.py index 692a1868f..d8c859f17 100644 --- a/autode/neb/ci.py +++ b/autode/neb/ci.py @@ -3,13 +3,14 @@ https://doi.org/10.1063/1.1329672 """ import numpy as np +from scipy.optimize import OptimizeResult from typing import Optional from autode.neb.original import NEB, Images, Image from autode.log import logger class CImage(Image): - def __init__(self, image): + def __init__(self, image: Image): """ Construct a climbing image from a non-climbing one @@ -17,11 +18,11 @@ def __init__(self, image): Arguments: image (autode.neb.Image): """ - super().__init__(image.name, image.k) + super().__init__(species=image, name=image.name, k=image.k) # Set all the current attributes from the regular image self.__dict__.update(image.__dict__) - def get_force(self, im_l, im_r) -> np.ndarray: + def get_force(self, im_l: Image, im_r: Image) -> np.ndarray: """ Compute F_m @@ -34,14 +35,13 @@ def get_force(self, im_l, im_r) -> np.ndarray: hat_tau, x_l, x, x_r = self._tau_xl_x_xr(im_l, im_r) # F_m = ∇V(x_m) + (2∇V(x_m).τ)τ - return -self.grad + 2.0 * np.dot(self.grad, hat_tau) * hat_tau + return -self.gradient + 2.0 * np.dot(self.gradient, hat_tau) * hat_tau class CImages(Images): def __init__( self, images: Images, - num: int = None, wait_iterations: int = 4, init_k: Optional[float] = None, ): @@ -55,11 +55,10 @@ def __init__( wait_iterations (int): Number of iterations to wait before turning on the climbing image - init_k: Initial force constant + init_k: Initial force constant. Must be defined if the images are empty """ super().__init__( - num=num if num is not None else len(images), - init_k=init_k if init_k is not None else images[0].k, + init_k=init_k if init_k is not None else images.init_k, ) self.wait_iteration = wait_iterations @@ -102,9 +101,7 @@ def __init__(self, *args, **kwargs): self.images = CImages(self.images) - def _minimise( - self, method, n_cores, etol, max_n=30 - ) -> "scipy.optimize.OptimizeResult": + def _minimise(self, method, n_cores, etol, max_n=30) -> OptimizeResult: """Minimise th energy of every image in the NEB""" logger.info(f"Minimising to ∆E < {etol:.4f} Ha on all NEB coordinates") result = super()._minimise(method, n_cores, etol, max_n) diff --git a/autode/neb/idpp.py b/autode/neb/idpp.py index 9b062abea..1b887247d 100644 --- a/autode/neb/idpp.py +++ b/autode/neb/idpp.py @@ -1,5 +1,10 @@ import numpy as np + from scipy.spatial import distance_matrix +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + from autode.neb.original import Images, Image class IDPP: @@ -17,7 +22,7 @@ class IDPP: as suggested in the paper. """ - def __init__(self, images: "autode.neb.original.Images"): + def __init__(self, images: "Images"): """Initialise a IDPP potential from a set of NEB images""" if len(images) < 2: @@ -29,7 +34,7 @@ def __init__(self, images: "autode.neb.original.Images"): self._set_distance_matrices(images) - def __call__(self, image: "autode.neb.original.Image") -> float: + def __call__(self, image: "Image") -> float: r""" Value of the IDPP objective function for a single image defined by, @@ -52,7 +57,7 @@ def __call__(self, image: "autode.neb.original.Image") -> float: return 0.5 * np.sum(w * (r_k - r) ** 2) - def grad(self, image: "autode.neb.original.Image") -> np.ndarray: + def grad(self, image: "Image") -> np.ndarray: r""" Gradient of the potential with respect to displacement of the Cartesian components: :math:`\nabla S = (dS/dx_0, dS/dy_0, dS/dz_0, @@ -75,7 +80,7 @@ def grad(self, image: "autode.neb.original.Image") -> np.ndarray: (np.ndarray): :math:`\nabla S` """ - x = np.array(image.species.coordinates).flatten() + x = np.array(image.coordinates).flatten() grad = np.zeros_like(x) r = self._distance_matrix(image, unity_diagonal=True) @@ -105,9 +110,9 @@ def grad(self, image: "autode.neb.original.Image") -> np.ndarray: grad[1::3] = np.sum(a * delta[1::3, 1::3], axis=1) # y grad[2::3] = np.sum(a * delta[2::3, 2::3], axis=1) # z - return grad + return grad.reshape((-1, 3)) - def _set_distance_matrices(self, images) -> None: + def _set_distance_matrices(self, images: "Images") -> None: """ For each image determine the optimum distance matrix using @@ -130,14 +135,16 @@ def _set_distance_matrices(self, images) -> None: self._diagonal_distance_matrix_idxs = np.diag_indices_from(delta) return None - def _req_distance_matrix(self, image): + def _req_distance_matrix(self, image: "Image"): """Required distance matrix for an image, with elements r_{ij}^k""" return self._dists[image.name] - def _distance_matrix(self, image, unity_diagonal=False) -> np.ndarray: + def _distance_matrix( + self, image: "Image", unity_diagonal: bool = False + ) -> np.ndarray: """Distance matrix for an image""" - x = image.species.coordinates + x = image.coordinates r = distance_matrix(x, x) if unity_diagonal: @@ -145,7 +152,7 @@ def _distance_matrix(self, image, unity_diagonal=False) -> np.ndarray: return r - def _weight_matrix(self, image) -> np.ndarray: + def _weight_matrix(self, image: "Image") -> np.ndarray: r""" Weight matrix with elements diff --git a/autode/neb/neb.py b/autode/neb/neb.py index 754f4f24a..3f5fe2228 100644 --- a/autode/neb/neb.py +++ b/autode/neb/neb.py @@ -1,18 +1,25 @@ import autode.exceptions as ex + from autode.config import Config from autode.log import logger from autode.neb.ci import CINEB from autode.transition_states.ts_guess import TSguess from autode.utils import work_in +from typing import TYPE_CHECKING, Optional + +if TYPE_CHECKING: + from autode.species.species import Species + from autode.wrappers.methods import Method + def get_ts_guess_neb( - reactant: "autode.species.Species", - product: "autode.species.Species", - method: "autode.wrappers.methods.Method", + reactant: "Species", + product: "Species", + method: "Method", name: str = "neb", n: int = 10, -): +) -> Optional[TSguess]: """ Get a transition state guess using a nudged elastic band calculation. The geometry of the reactant is used as the fixed initial point and the final @@ -37,17 +44,13 @@ def get_ts_guess_neb( assert n is not None logger.info("Generating a TS guess using a nudged elastic band") - neb = CINEB( - initial_species=reactant.copy(), final_species=product.copy(), num=n - ) + neb = CINEB.from_end_points(reactant, product, num=n) - # Calculate and generate the TS guess, in a unique directory - try: - - @work_in(name) - def calculate(): - return neb.calculate(method=method, n_cores=Config.n_cores) + @work_in(name) + def calculate(): + return neb.calculate(method=method, n_cores=Config.n_cores) + try: calculate() except ex.CalculationException: @@ -55,7 +58,7 @@ def calculate(): return None return TSguess( - atoms=neb.get_species_saddle_point().atoms, + atoms=neb.peak_species.atoms, reactant=reactant, product=product, name=name, diff --git a/autode/neb/original.py b/autode/neb/original.py index 3d161f02c..88d15f59c 100644 --- a/autode/neb/original.py +++ b/autode/neb/original.py @@ -2,22 +2,26 @@ The theory behind this original NEB implementation is taken from Henkelman and H. J ́onsson, J. Chem. Phys. 113, 9978 (2000) """ -from typing import Optional, Sequence +import numpy as np +import matplotlib.pyplot as plt + +from typing import Optional, Sequence, List, Any, TYPE_CHECKING +from copy import deepcopy + from autode.log import logger from autode.calculations import Calculation from autode.wrappers.methods import Method +from autode.species.species import Species from autode.input_output import xyz_file_to_molecules from autode.path import Path from autode.utils import work_in, ProcessPool from autode.config import Config from autode.neb.idpp import IDPP from scipy.optimize import minimize -from autode.values import Distance, PotentialEnergy -from copy import deepcopy -import numpy as np -import matplotlib.pyplot as plt +from autode.values import Distance, PotentialEnergy, ForceConstant, Gradient -blues = plt.get_cmap("Blues") +if TYPE_CHECKING: + from autode.wrappers.methods import Method def energy_gradient(image, method, n_cores): @@ -39,7 +43,7 @@ def _est_energy_gradient(image, est_method, n_cores): """Electronic structure energy and gradint""" calc = Calculation( name=f"{image.name}_{image.iteration}", - molecule=image.species, + molecule=image, method=est_method, keywords=est_method.keywords.grad, n_cores=n_cores, @@ -48,19 +52,16 @@ def _est_energy_gradient(image, est_method, n_cores): @work_in(image.name) def run(): calc.run() - image.grad = image.species.gradient.flatten() - image.energy = image.species.energy - return None run() return image def _idpp_energy_gradient( - image: "autode.neb.original.Image", - idpp: "autode.neb.idpp.IDPP", + image: "Image", + idpp: IDPP, n_cores: int, -) -> "autode.neb.original.Image": +) -> "Image": """ Evaluate the energy and gradient of an image using an image dependent pair potential IDDP instance and set the energy and gradient on the image @@ -77,7 +78,7 @@ def _idpp_energy_gradient( (autode.neb.original.Image): Image """ image.energy = idpp(image) - image.grad = idpp.grad(image) + image.gradient = idpp.grad(image) return image @@ -127,7 +128,7 @@ def derivative(flat_coords, images, method, n_cores, plot_energies): """ # Forces for the first image are fixed at zero - forces = np.array(images[0].grad) + forces = np.zeros(shape=images[0].gradient.shape) # No need to calculate gradient as should already be there from energy eval for i in range(1, len(images) - 1): @@ -135,37 +136,51 @@ def derivative(flat_coords, images, method, n_cores, plot_energies): forces = np.append(forces, force) # Final zero set of forces - forces = np.append(forces, images[-1].grad) + forces = np.append(forces, np.zeros(shape=images[-1].gradient.shape)) # dV/dx is negative of the force logger.info(f"|F| = {np.linalg.norm(forces):.4f} Ha Å-1") return -forces -class Image: - def __init__(self, name: str, k: float): +class Image(Species): + def __init__( + self, + species: Species, + name: str, + k: ForceConstant, + ): """ Image in a NEB + -------------------------------------------------------------------------------- Arguments: - name (str): - """ - self.name = name + species (Species): Molecule for which this image represents + + name (str): Name of this image - # Current optimisation iteration of this image - self.iteration = 0 + k (ForceConstant): Force constant of the harmonic potential joining this + image to its neighbour(s) + """ + super().__init__( + name=name, + charge=species.charge, + mult=species.mult, + atoms=species.atoms.copy(), + ) + self.solvent = deepcopy(species.solvent) + self.energy = deepcopy(species.energy) - # Force constant in Eh/Å^2 + self.iteration = 0 #: Current optimisation iteration of this image self.k = k - self.species: Optional["autode.species.Species"] = None - self.energy: Optional["autode.values.Energy"] = None - self.grad: Optional[np.ndarray] = None + def _generate_conformers(self, *args, **kwargs): + raise RuntimeError("Cannot create conformers of an image") def _tau_xl_x_xr( self, - im_l: "autode.neb.original.Image", - im_r: "autode.neb.original.Image", + im_l: "Image", + im_r: "Image", ) -> tuple: """ Calculate the normalised τ vector, along with the coordinates of the @@ -194,7 +209,7 @@ def _tau_xl_x_xr( # x_i-1, x_i, x_i+1 x_l, x, x_r = [ - image.species.coordinates.flatten() for image in (im_l, self, im_r) + image.coordinates.flatten() for image in (im_l, self, im_r) ] # τ_i+ tau_plus = x_r - x @@ -221,8 +236,8 @@ def _tau_xl_x_xr( def get_force( self, - im_l: "autode.neb.original.Image", - im_r: "autode.neb.original.Image", + im_l: "Image", + im_r: "Image", ) -> np.ndarray: """ Compute F_i. Notation from: @@ -244,32 +259,47 @@ def get_force( ) * hat_tau # ∇V(x)_i|_|_ = ∇V(x)_i - (∇V(x)_i•τ) τ - grad_perp = self.grad - np.dot(self.grad, hat_tau) * hat_tau + grad_perp = self.gradient - np.dot(self.gradient, hat_tau) * hat_tau # F_i = F_i^s|| - ∇V(x)_i|_|_ return f_parallel - grad_perp + @property + def gradient(self) -> Optional[np.ndarray]: + return None if self._grad is None else self._grad.flatten() + + @gradient.setter + def gradient(self, value: Optional[np.ndarray]): + self._grad = None if value is None else value.flatten() + class Images(Path): - def __init__(self, num, init_k, min_k=None, max_k=None): + def __init__( + self, + init_k: ForceConstant, + min_k: Optional[ForceConstant] = None, + max_k: Optional[ForceConstant] = None, + ): """ Set of images joined by harmonic springs with force constant k ----------------------------------------------------------------------- Arguments: - num (int): Number of images + init_k (ForceConstant): Initial force constant - init_k (float): Initial force constant (Ha Å^-2) + min_k (ForceConstant | None): Minimum value of k - min_k (None | float): Minimum value of k - - max_k (None | float): Maximum value of k + max_k (ForceConstant | None): Maximum value of k """ - super().__init__(*(Image(name=str(i), k=init_k) for i in range(num))) + super().__init__() - self.min_k = init_k / 10 if min_k is None else float(min_k) - self.max_k = 2 * init_k if max_k is None else float(max_k) + self.init_k = init_k + self.min_k = init_k / 10 if min_k is None else min_k + self.max_k = 2 * init_k if max_k is None else max_k + assert ( + self.max_k > self.min_k + ), "Can't set the min force constant above the max" def __eq__(self, other): """Equality od two climbing image NEB paths""" @@ -310,7 +340,7 @@ def increment(self): image.k = self.min_k else: - image.k = self.max_k - delta_k * ( + image.k = self.max_k - delta_k * float( (e_max - image.energy) / (e_max - e_ref) ) return None @@ -319,6 +349,8 @@ def plot_energies( self, save=False, name="None", color=None, xlabel="NEB coordinate" ): """Plot the NEB surface""" + blues = plt.get_cmap("Blues") + color = ( blues((self[0].iteration + 1) / 20) if color is None @@ -331,7 +363,7 @@ def coords(self): coords = np.array([]) for image in self: - coords = np.append(coords, image.species.coordinates.flatten()) + coords = np.append(coords, image.coordinates.flatten()) return coords def set_coords(self, coords): @@ -343,14 +375,20 @@ def set_coords(self, coords): coords (np.ndarray): shape (num x n x 3,) """ - n_atoms = self[0].species.n_atoms + n_atoms = self[0].n_atoms coords = coords.reshape((len(self), n_atoms, 3)) for i, image in enumerate(self): - image.species.coordinates = coords[i] + image.coordinates = coords[i] return None + def append_species(self, species: Species) -> None: + """Add a species to the list of images""" + super().append( + Image(species=species, name=f"{len(self)}", k=self.init_k) + ) + def copy(self) -> "Images": return deepcopy(self) @@ -361,51 +399,33 @@ class NEB: def __init__( self, - initial_species: Optional["autode.Species"] = None, - final_species: Optional["autode.Species"] = None, - num: int = 8, - species_list: Optional[list] = None, - k: float = 0.1, + init_k: ForceConstant = ForceConstant(0.1, units="Ha / Å^2"), + **kwargs, ): """ - Nudged elastic band + Nudged elastic band (NEB) + + Warning: The initial/final species or those in a species list must have + the same atom ordering. ----------------------------------------------------------------------- Arguments: - initial_species (autode.species.Species): - - final_species (autode.species.Species): - - num (int): Number of images in the NEB - - species_list (list(autode.species.Species)): Intermediate images - along the NEB + init_k: Initial force constant between each image """ - self._init_k = k - - if species_list is not None: - self.images = Images(num=len(species_list), init_k=k) - - for image, species in zip(self.images, species_list): - image.species = species - image.energy = species.energy - - else: - self.images = Images(num=num, init_k=k) - self._init_from_end_points(initial_species, final_species) - - logger.info(f"Initialised a NEB with {len(self.images)} images") + self._raise_exception_if_any(kwargs) + self._init_k = init_k + self.images = Images(init_k=init_k) @property - def init_k(self) -> float: - """Initial force constant used to in this NEB. Units of Ha Å^-1""" + def init_k(self) -> ForceConstant: + """Initial force constant used to in this NEB""" return self._init_k @classmethod def from_file( cls, filename: str, - k: Optional[float] = None, + init_k: Optional[float] = None, ) -> "NEB": """ Create a nudged elastic band from a .xyz file containing multiple @@ -413,7 +433,7 @@ def from_file( """ molecules = xyz_file_to_molecules(filename) - if k is None and all(m.energy is not None for m in molecules): + if init_k is None and all(m.energy is not None for m in molecules): logger.info( "Have a set of energies from file. Can adaptively " "choose a sensible force constant (k)" @@ -426,18 +446,86 @@ def from_file( # TODO: test reasonableness of this function... # use a shifted tanh to interpolate in [0.005, 0.2005] - k = ( + init_k = ForceConstant( 0.1 * (np.tanh((max_de.to("kcal mol-1") - 40) / 20) + 1) - + 0.005 + + 0.005, + units="Ha / Å^2", + ) + + if init_k is None: # choose a sensible default + init_k = 0.1 + + logger.info( + f"Using k = {init_k:.6f} Ha Å^-1 as the NEB force constant" + ) + return cls.from_list(molecules, init_k=init_k) + + @classmethod + def from_list( + cls, + species_list: List[Species], + init_k: ForceConstant = ForceConstant(0.1, units="Ha / Å^2"), + ) -> "NEB": + """ + Nudged elastic band constructed from list of species + + ----------------------------------------------------------------------- + Arguments: + species_list: Full set of initial images that will form the while NEB + + init_k: Force constant + + Returns: + (NEB): + """ + neb = cls(init_k=init_k) + + for species in species_list: + neb.images.append_species(species) + + logger.info(f"Initialised a NEB with {len(neb.images)} images") + return neb + + @classmethod + def from_end_points( + cls, + initial: Species, + final: Species, + num: int, + init_k: ForceConstant = ForceConstant(0.1, units="Ha / Å^2"), + ) -> "NEB": + """ + Construct a nudged elastic band from only the endpoints. The atomic + ordering must be identical in the initial and final species + + ----------------------------------------------------------------------- + Arguments: + initial: Initial/left-most species in the NEB + + final: Final/right-most species in the NEB + + num: Number of images to create + + init_k: Initial force constant + + Returns: + (NEB): + """ + + if initial.sorted_atomic_symbols != final.sorted_atomic_symbols: + raise ValueError( + "Cannot construct a NEB from species with different atoms" ) - if k is None: # choose a sensible default - k = 0.1 + neb = cls.from_list( + species_list=cls._interpolated_species(initial, final, n=num), + init_k=init_k, + ) + neb.idpp_relax() - logger.info(f"Using k = {k:.6f} Ha Å^-1 as the NEB force constant") - return cls(species_list=molecules, k=k) + return neb - def _minimise(self, method, n_cores, etol, max_n=30) -> None: + def _minimise(self, method, n_cores, etol, max_n=30) -> Any: """Minimise the energy of every image in the NEB""" logger.info(f"Minimising to ∆E < {etol:.4f} Ha on all NEB coordinates") @@ -454,10 +542,6 @@ def _minimise(self, method, n_cores, etol, max_n=30) -> None: logger.info(f"NEB path energy = {result.fun:.5f} Ha, {result.message}") return result - def contains_peak(self) -> bool: - """Does this nudged elastic band calculation contain an energy peak?""" - return self.images.contains_peak - def partition( self, max_delta: Distance, @@ -483,84 +567,77 @@ def partition( logger.info("Interpolating") assert len(self.images) > 1 - init_k = self.images[0].k _list = [] for i, left_image in enumerate(self.images[:-1]): right_image = self.images[i + 1] - left_species, right_species = ( - left_image.species, - right_image.species, - ) - sub_neb = NEB(species_list=[left_species, right_species]) - n = 2 + sub_neb = NEB.from_end_points(left_image, right_image, num=n) while ( sub_neb._max_atom_distance_between_images(distance_idxs) > max_delta ): - sub_neb.images = self._images_type(num=n, init_k=init_k) - sub_neb.images[0].species = left_species - sub_neb.images[-1].species = right_species - sub_neb.interpolate_geometries() - try: - sub_neb.idpp_relax() + sub_neb = NEB.from_end_points( + left_image, right_image, num=n + ) except RuntimeError: logger.warning("Failed to IDPP relax the interpolated NEB") n += 1 for image in sub_neb.images[:-1]: # add all apart from the last - _list.append(image.species) + _list.append(image) - _list.append(self._species_at(-1)) # end with the last - self.images = self._images_type(num=len(_list), init_k=init_k) + _list.append(self.images[-1]) # end with the last + self.images.clear() - for i, image in enumerate(self.images): - image.species = _list[i] + for image in _list: + self.images.append_species(image) logger.info( f"Partition successful – now have {len(self.images)} " f"images" ) return None - def _species_at(self, idx: int) -> "autode.species.Species": - """Return the species associated with a particular image index""" - return self.images[idx].species - def print_geometries(self, name="neb") -> None: return self.images.print_geometries(name) - def interpolate_geometries(self) -> None: + @staticmethod + def _interpolated_species( + initial: Species, final: Species, n: int + ) -> List[Species]: """Generate simple interpolated coordinates for these set of images in Cartesian coordinates""" - n = len(self.images) + + intermediate_species = [] # Interpolate images between the starting point i=0 and end point i=n-1 for i in range(1, n - 1): # Use a copy of the starting point for atoms, charge etc. - self.images[i].species = self.images[0].species.copy() + species = initial.copy() # For all the atoms in the species translate an amount so the # spacing is even between the initial and final points - for j, atom in enumerate(self.images[i].species.atoms): + for j, atom in enumerate(species.atoms): # Shift vector is final minus current - shift = self.images[-1].species.atoms[j].coord - atom.coord + shift = final.atoms[j].coord - atom.coord # then an equal spacing is the i-th point in the grid atom.translate(vec=shift * (i / n)) - return None + intermediate_species.append(species) + + return [initial.copy()] + intermediate_species + [final.copy()] @work_in("neb") def calculate( self, - method: "autode.wrappers.methods.Method", + method: "Method", n_cores: int, name_prefix: str = "", etol_per_image: PotentialEnergy = PotentialEnergy( @@ -585,11 +662,9 @@ def calculate( """ self.print_geometries(name=f"{name_prefix}neb_init") - # Calculate energy on the first and final points + # Calculate energy on the first and final points as these will not be recalc-ed for idx in [0, -1]: energy_gradient(self.images[idx], method=method, n_cores=n_cores) - # Zero the forces so the end points don't move - self.images[idx].grad = np.zeros(shape=self.images[idx].grad.shape) if isinstance(etol_per_image, PotentialEnergy): etol_per_image = float( @@ -609,23 +684,15 @@ def calculate( plt.close() return None - def get_species_saddle_point(self) -> "autode.species.Species": - """Find a TS guess species for this NEB: highest energy saddle point""" - if self.images.peak_idx is None: + @property + def peak_species(self) -> Optional[Species]: + """TS guess species for this NEB: highest energy saddle point""" + if not self.images.contains_peak: logger.warning("Found no peaks in the NEB") return None - return self.images[self.images.peak_idx].species - - def _init_from_end_points(self, initial, final) -> None: - """Initialise from the start and finish points of the NEB""" - - self.images[0].species = initial - self.images[-1].species = final - self.interpolate_geometries() - self.idpp_relax() - - return None + image = self.images[self.images.peak_idx] + return image.new_species() def idpp_relax(self) -> None: """ @@ -636,18 +703,18 @@ def idpp_relax(self) -> None: :py:meth:`IDPP ` """ logger.info(f"Minimising NEB with IDPP potential") - images = self.images.copy() - idpp = IDPP(images) - images.min_k = images.max_k = 0.1 + images = self.images.copy() + images.min_k = images.max_k = ForceConstant(0.1, units="Ha / Å^2") + idpp = IDPP(images=images) for i, image in enumerate(images): image.energy = idpp(image) - image.grad = idpp.grad(image) + image.gradient = idpp.grad(image) # Initial and final images are fixed, with zero gradient if i == 0 or i == len(images) - 1: - image.grad[:] = 0.0 + image.gradient[:] = 0.0 result = minimize( total_energy, @@ -670,7 +737,7 @@ def _max_atom_distance_between_images( Calculate the maximum atom-atom distance between two consecutive images """ if idxs is None: # Use all pairwise distances - idxs = np.arange(self.images[0].species.n_atoms) + idxs = np.arange(self.images[0].n_atoms) else: idxs = np.array(idxs) @@ -678,8 +745,8 @@ def _max_atom_distance_between_images( for i in range(len(self.images) // 2): k = 2 * i - x_i = self._species_at(k).coordinates - x_j = self._species_at(k + 1).coordinates + x_i = self.images[k].coordinates + x_j = self.images[k + 1].coordinates max_distance = np.max(np.linalg.norm(x_i - x_j, axis=1)[idxs]) if max_distance > overall_max_distance: @@ -690,3 +757,22 @@ def _max_atom_distance_between_images( @property def max_atom_distance_between_images(self) -> Distance: return self._max_atom_distance_between_images(idxs=None) + + @staticmethod + def _raise_exception_if_any(kwargs: dict) -> None: + + if len(kwargs) == 0: + return + elif any( + arg in kwargs + for arg in ("initial_species", "final_species", "num") + ): + raise ValueError( + "Cannot construct a NEB. Please use NEB.from_endpoints()" + ) + elif "species_list" in kwargs: + raise ValueError( + "Cannot construct a NEB from a species list. Please use NEB.from_list()" + ) + else: + raise ValueError("Unrecognised keyword argument") diff --git a/autode/opt/optimisers/dimer.py b/autode/opt/optimisers/dimer.py index 92e592b17..cefd8f15c 100644 --- a/autode/opt/optimisers/dimer.py +++ b/autode/opt/optimisers/dimer.py @@ -273,11 +273,10 @@ def _update_gradient_at(self, point: DimerPoint) -> None: ) calc.run() - self._coords.e = self._species.energy = calc.get_energy() - self._species.gradient = calc.get_gradients() + self._coords.e = self._species.energy self._coords.set_g_at( - point, calc.get_gradients().flatten(), mass_weighted=False + point, self._species.gradient.flatten(), mass_weighted=False ) calc.clean_up(force=True, everything=True) diff --git a/autode/path/__init__.py b/autode/path/__init__.py index 2136e25d2..d90d12446 100644 --- a/autode/path/__init__.py +++ b/autode/path/__init__.py @@ -1,6 +1,5 @@ from autode.path.path import Path from autode.path.adaptive import AdaptivePath -from autode.path.adaptive import PathPoint -__all__ = ["Path", "AdaptivePath", "PathPoint"] +__all__ = ["Path", "AdaptivePath"] diff --git a/autode/path/adaptive.py b/autode/path/adaptive.py index c79613584..78fde9f04 100644 --- a/autode/path/adaptive.py +++ b/autode/path/adaptive.py @@ -1,20 +1,28 @@ import autode as ade import numpy as np -from copy import deepcopy + from autode.log import logger -from autode.mol_graphs import make_graph from autode.path.path import Path from autode.transition_states.ts_guess import TSguess from autode.utils import work_in +from autode.bonds import ScannedBond + +from typing import TYPE_CHECKING, List, Optional + + +if TYPE_CHECKING: + from autode.species import ReactantComplex, ProductComplex, Species + from autode.wrappers.methods import Method + from autode.bond_rearrangement import BondRearrangement def get_ts_adaptive_path( - reactant: "autode.species.ReactantComplex", - product: "autode.species.ProductComplex", - method: "autode.wrappers.methods.Method", - bond_rearr: "autode.bond_rearrangement.BondRearrangement", + reactant: "ReactantComplex", + product: "ProductComplex", + method: "Method", + bond_rearr: "BondRearrangement", name: str = "adaptive", -) -> "autode.transition_states.TSguess ": +) -> Optional[TSguess]: """ Generate a TS guess geometry based on an adaptive path along multiple breaking and/or forming bonds @@ -49,7 +57,7 @@ def get_ts_adaptive_path( return None ts_guess = TSguess( - atoms=ts_path[ts_path.peak_idx].species.atoms, + atoms=ts_path[ts_path.peak_idx].atoms, reactant=reactant, product=product, bond_rearr=bond_rearr, @@ -58,7 +66,9 @@ def get_ts_adaptive_path( return ts_guess -def pruned_active_bonds(reactant, fbonds, bbonds): +def pruned_active_bonds( + reactant: "ReactantComplex", fbonds: list, bbonds: list +) -> List[ScannedBond]: """ Prune the set of forming and breaking bonds for special cases @@ -126,38 +136,14 @@ def pruned_active_bonds(reactant, fbonds, bbonds): return fbonds + bbonds -class PathPoint: - def __init__(self, species, constraints): - """ - Point on a PES path - - ----------------------------------------------------------------------- - Arguments: - species (autode.species.Species): - - constraints (dict): Distance constraints keyed with atom indexes as - a tuple with distances (floats) as values - """ - self.species = species - self.species.constraints.distance = constraints - - self.energy = None # Ha - self.grad = None # Ha Å^-1 - - @property - def constraints(self) -> dict: - return self.species.constraints.distance - - def copy(self): - """Return a copy of this point""" - return PathPoint( - species=self.species.new_species(), - constraints=deepcopy(self.species.constraints.distance), - ) - - class AdaptivePath(Path): - def __init__(self, init_species, bonds, method, final_species=None): + def __init__( + self, + bonds: List[ScannedBond], + method: "Method", + init_species: Optional["Species"] = None, + final_species: Optional["Species"] = None, + ): """ PES Path @@ -177,18 +163,15 @@ def __init__(self, init_species, bonds, method, final_species=None): self.bonds = bonds self.final_species = final_species - # Bonds need to have the initial and final dists to drive along them - for bond in bonds: - assert bond.curr_dist is not None and bond.final_dist is not None - # Add the first point - will run a constrained minimisation if possible - init_point = PathPoint( - species=init_species, - constraints={bond.atom_indexes: bond.curr_dist for bond in bonds}, - ) + if init_species is not None: + point = init_species.new_species() + point.constraints.distance = { + b.atom_indexes: b.curr_dist for b in bonds + } + self.append(point) - if init_point.species.n_atoms > 0: - self.append(init_point) + self._check_bonds_have_initial_and_final_distances() def __eq__(self, other): """Equality of two adaptive paths""" @@ -197,19 +180,22 @@ def __eq__(self, other): return super().__eq__(other) + def _check_bonds_have_initial_and_final_distances(self) -> None: + for bond in self.bonds: + assert bond.curr_dist is not None and bond.final_dist is not None + @work_in("initial_path") def append(self, point) -> None: """ - Append a point to the path and optimise it + Append a point to the path and optimise it ----------------------------------------------------------------------- Arguments: - point (PathPoint): Point on a path + point (Species): Point on a path Raises: (autode.exceptions.CalculationException): """ - super().append(point) idx = len(self) - 1 keywords = self.method.keywords.low_opt.copy() @@ -217,38 +203,33 @@ def append(self, point) -> None: calc = ade.Calculation( name=f"path_opt{idx}", - molecule=self[idx].species, + molecule=point, method=self.method, keywords=keywords, n_cores=ade.Config.n_cores, ) calc.run() - - # Set the required properties from the calculation - self[idx].species.atoms = calc.get_final_atoms() - make_graph(self[idx].species) - self[idx].energy = calc.get_energy() + point.reset_graph() if self.method.name == "xtb" or self.method.name == "mopac": # XTB prints gradients including the constraints, which are ~0 # the gradient here is just the derivative of the electronic energy # so rerun a gradient calculation, which should be very fast # while MOPAC doesn't print gradients for a constrained opt + tmp_point_for_grad = point.new_species() calc = ade.Calculation( name=f"path_grad{idx}", - molecule=self[idx].species.new_species(), + molecule=tmp_point_for_grad, method=self.method, keywords=self.method.keywords.grad, n_cores=ade.Config.n_cores, ) calc.run() - self[idx].grad = calc.get_gradients() calc.clean_up(force=True, everything=True) + assert tmp_point_for_grad.gradient is not None + point.gradient = tmp_point_for_grad.gradient - else: - self[idx].grad = calc.get_gradients() - - return None + return super().append(point) def plot_energies( self, save=True, name="init_path", color="k", xlabel="ζ" @@ -265,10 +246,9 @@ def contains_suitable_peak(self) -> bool: logger.info("Products made and have a peak. Assuming suitable!") return True - # Products aren't made by isomorphism but we may still have a suitable - # peak.. + # Products aren't made by isomorphism, but we may still have a suitable peak if any( - self[-1].constraints[b.atom_indexes] == b.final_dist + self[-1].constraints.distance[b.atom_indexes] == b.final_dist for b in self.bonds ): logger.warning( @@ -296,7 +276,7 @@ def _adjust_constraints(self, point): max_step, min_step = ade.Config.max_step_size, ade.Config.min_step_size for bond in self.bonds: - (i, j), coords = bond.atom_indexes, self[-1].species.coordinates + (i, j), coords = bond.atom_indexes, self[-1].coordinates # Normalised r_ij vector vec = coords[j] - coords[i] @@ -305,8 +285,8 @@ def _adjust_constraints(self, point): # Calculate |∇E_i·r| i.e. the gradient along the bond. Positive # values are downhill in energy to form the bond and negative # downhill to break it - gradi = np.dot(self[-1].grad[i], vec) # |∇E_i·r| bond midpoint - gradj = np.dot(self[-1].grad[j], -vec) + gradi = np.dot(self[-1].gradient[i], vec) # |∇E_i·r| bond midpoint + gradj = np.dot(self[-1].gradient[j], -vec) # Exclude gradients from atoms that are being substituted if atom_idxs.count(i) > 1: @@ -330,7 +310,7 @@ def _adjust_constraints(self, point): ) + min_step dr *= np.sign(bond.dr) - new_dist = point.species.distance(*bond.atom_indexes) + dr + new_dist = point.distance(*bond.atom_indexes) + dr # No need to go exceed final distances on forming/breaking bonds if bond.forming and new_dist < bond.final_dist: @@ -342,7 +322,7 @@ def _adjust_constraints(self, point): else: logger.info(f"Using step {dr:.3f} Å on bond: {bond}") - point.constraints[bond.atom_indexes] = new_dist + point.constraints.distance[bond.atom_indexes] = new_dist return None @@ -362,13 +342,13 @@ def generate(self, init_step_size=0.2, name="initial") -> None: # Always perform an initial step linear in all bonds logger.info("Performing a linear step and calculating gradients") - point = self[0].copy() + point = self[0].new_species(with_constraints=True) for bond in self.bonds: # Shift will be -min_step_size if ∆r is negative and larger than # the minimum step size dr = np.sign(bond.dr) * min(init_step_size, np.abs(bond.dr)) - point.constraints[bond.atom_indexes] += dr + point.constraints.distance[bond.atom_indexes] += dr self.append(point) logger.info("First point found") @@ -376,14 +356,14 @@ def generate(self, init_step_size=0.2, name="initial") -> None: def reached_final_point(): """Are there any more points to add?""" return all( - point.constraints[b.atom_indexes] == b.final_dist + point.constraints.distance[b.atom_indexes] == b.final_dist for b in self.bonds ) logger.info("Adaptively adding points to the path") while not (reached_final_point() or self.contains_suitable_peak()): - point = self[-1].copy() + point = self[-1].new_species(with_constraints=True) self._adjust_constraints(point=point) self.append(point) diff --git a/autode/path/path.py b/autode/path/path.py index 6378b08c3..b1a2ebcda 100644 --- a/autode/path/path.py +++ b/autode/path/path.py @@ -1,21 +1,23 @@ import numpy as np import matplotlib.pyplot as plt -from typing import Optional -from autode import mol_graphs + +from autode.species import Species from autode.input_output import atoms_to_xyz_file from autode.log import logger from autode.units import KcalMol +from typing import Optional + class Path(list): - def __init__(self, *args, units=KcalMol): + def __init__(self, *args: Species, units=KcalMol): """ Base path class that may be populated with species or nudged elastic band images, *must* have .energy attributes ----------------------------------------------------------------------- Arguments: - args (list(autode.path.PathPoint | autode.neb.Image)): + args (autode.species.species.Species): Keyword Arguments: units (autode.units.Unit): @@ -23,12 +25,7 @@ def __init__(self, *args, units=KcalMol): super().__init__() for arg in args: - if not (hasattr(arg, "energy") or hasattr(arg, "species")): - raise ValueError( - "A Path must be initialised from a class " - "with both energy and species attributes" - ) - + assert isinstance(arg, Species) self.append(arg) self.units = units @@ -95,7 +92,7 @@ def peak_idx(self) -> Optional[int]: def contains_peak(self) -> bool: return self.peak_idx is not None - def product_idx(self, product) -> Optional[int]: + def product_idx(self, product: "Species") -> Optional[int]: """ Get the index of the point in the path at which products are made. If they are not made or they cannot be checked then return None @@ -112,9 +109,7 @@ def product_idx(self, product) -> Optional[int]: return None for i, point in enumerate(self): - if mol_graphs.is_isomorphic( - graph1=point.species.graph, graph2=product.graph - ): + if product.graph.is_isomorphic_to(point.graph): logger.info(f"Products made at point {i}") return i @@ -175,20 +170,19 @@ def print_geometries(self, name) -> None: open(f"{name}.xyz", "w").close() # Empty the file for i, image in enumerate(self): - assert image.species is not None energy = image.energy if image.energy is not None else "none" title_line = ( f"autodE path point {i}. E = {energy} " - f"charge = {image.species.charge} " - f"mult = {image.species.mult} " + f"charge = {image.charge} " + f"mult = {image.mult} " ) - if image.species.solvent is not None: - title_line += f"solvent = {image.species.solvent.name} " + if image.solvent is not None: + title_line += f"solvent = {image.solvent.name} " atoms_to_xyz_file( - image.species.atoms, + image.atoms, f"{name}.xyz", title_line=title_line, append=True, diff --git a/autode/species/species.py b/autode/species/species.py index d63060573..77dfa5a6b 100644 --- a/autode/species/species.py +++ b/autode/species/species.py @@ -126,22 +126,30 @@ def copy(self) -> "Species": """Copy this whole species""" return deepcopy(self) - def new_species(self, name="species") -> "Species": + def new_species( + self, name="species", with_constraints: bool = False + ) -> "Species": """ A new version of this species, identical properties without any energies, gradients, hessian, conformers or constraints. ----------------------------------------------------------------------- - Keyword Arguments: + Arguments: name (str): Name of the new species + with_constraints (bool): Should the constraints from this species be copied + into the new one + Returns: (autode.species.Species): """ - species = Species(name, deepcopy(self.atoms), self.charge, self.mult) + species = Species(name, self.atoms.copy(), self.charge, self.mult) species.graph = None if self.graph is None else self.graph.copy() species.solvent = None if self.solvent is None else self.solvent.copy() + if with_constraints: + species.constraints = self.constraints.copy() + return species @property diff --git a/autode/transition_states/locate_tss.py b/autode/transition_states/locate_tss.py index ec8178d18..e1163f834 100644 --- a/autode/transition_states/locate_tss.py +++ b/autode/transition_states/locate_tss.py @@ -340,7 +340,7 @@ def _get_ts_neb_from_adaptive_path( neb = NEB.from_file(f"{ad_name}_path.xyz") - if not neb.contains_peak(): + if not neb.images.contains_peak: logger.info("Adaptive path had no peak – not running a NEB") return None @@ -356,9 +356,8 @@ def _get_ts_neb_from_adaptive_path( ) if neb.images.contains_peak: - peak_idx = neb.images.peak_idx ts_guess = TSguess( - atoms=neb.images[peak_idx].species.atoms, + atoms=neb.peak_species.atoms, reactant=reactant, product=product, bond_rearr=bond_rearr, diff --git a/autode/transition_states/transition_state.py b/autode/transition_states/transition_state.py index 792a4172d..e02a0c024 100644 --- a/autode/transition_states/transition_state.py +++ b/autode/transition_states/transition_state.py @@ -92,7 +92,7 @@ def _run_opt_ts_calc(self, method, name_ext): try: optts_calc.run() - if not optts_calc.optimisation_converged(): + if not optts_calc.optimiser.converged: self._reoptimise(optts_calc, name_ext, method) except CalculationException: @@ -100,32 +100,34 @@ def _run_opt_ts_calc(self, method, name_ext): return - def _reoptimise(self, calc, name_ext, method) -> Calculation: + def _reoptimise( + self, calc: Calculation, name_ext: str, method: "Method" + ) -> Calculation: """Rerun a calculation for more steps""" - if not calc.optimisation_nearly_converged(): + if calc.optimiser.last_energy_change.to("kcal mol-1") > 0.1: self.warnings += f"TS for {self.name} was not fully converged." logger.info("Optimisation did not converge") return calc logger.info("Optimisation nearly converged") - if self.could_have_correct_imag_mode: - logger.info( - "Still have correct imaginary mode, trying " - "more optimisation steps" - ) + if not self.could_have_correct_imag_mode: + logger.warning("Lost imaginary mode") + return calc - self.atoms = calc.get_final_atoms() - calc = Calculation( - name=f"{self.name}_{name_ext}_reopt", - molecule=self, - method=method, - n_cores=Config.n_cores, - keywords=method.keywords.opt_ts, - ) - calc.run() - else: - logger.info("Lost imaginary mode") + logger.info( + "Still have correct imaginary mode, trying " + "more optimisation steps" + ) + + calc = Calculation( + name=f"{self.name}_{name_ext}_reopt", + molecule=self, + method=method, + n_cores=Config.n_cores, + keywords=method.keywords.opt_ts, + ) + calc.run() return calc diff --git a/autode/values.py b/autode/values.py index 77fc41ca1..8de94822a 100644 --- a/autode/values.py +++ b/autode/values.py @@ -1,40 +1,24 @@ import numpy as np from abc import ABC, abstractmethod -from typing import Any, Union, Type, Optional, Sequence, Tuple +from typing import Any, Union, Type, Optional, Sequence, TYPE_CHECKING from copy import deepcopy from collections.abc import Iterable from autode.log import logger + +# fmt: off from autode.units import ( - Unit, - ha, - kjmol, - kcalmol, - ev, - J, - ang, - a0, - nm, - pm, - m, - ang_amu_half, - rad, - deg, - wavenumber, - hz, - amu, - kg, - m_e, - amu_ang_sq, - kg_m_sq, - ha_per_ang, - ha_per_a0, - ev_per_ang, - kcalmol_per_ang, - byte, - MB, - GB, - TB, + Unit, ha, m, ang_amu_half, ha_per_a0, ev_per_ang, + kjmol, kcalmol, rad, deg, kcalmol_per_ang, byte, + ev, J, wavenumber, hz, MB, ha_per_ang, + ang, a0, amu, kg, GB, kg_m_sq, + nm, pm, m_e, amu_ang_sq, TB, ha_per_a0_sq, + ha_per_ang_sq, J_per_m_sq, J_per_ang_sq, J_per_ang_sq_kg, ) +# fmt: on + +if TYPE_CHECKING: + from autode.wrappers.methods import Method + from autode.wrappers.keywords.keywords import Keywords def _to( @@ -285,10 +269,10 @@ class Energy(Value): def __init__( self, - value, + value: Any, units: Union[Unit, str] = ha, - method: Optional["autode.wrappers.methods.Method"] = None, - keywords: Optional["autode.wrappers.keywords.Keywords"] = None, + method: Optional["Method"] = None, + keywords: Optional["Keywords"] = None, estimated: bool = False, ): """ @@ -340,8 +324,8 @@ def __eq__(self, other): def set_method_str( self, - method: Optional["autode.wrappers.methods.Method"], - keywords: Optional["autode.wrappers.keywords.Keywords"], + method: Optional["Method"], + keywords: Optional["Keywords"], ) -> None: self.method_str = method_string(method, keywords) @@ -575,6 +559,23 @@ def __init__(self, value, units=amu): super().__init__(value, units=units) +class ForceConstant(Value): + + implemented_units = [ + ha_per_ang_sq, + ha_per_a0_sq, + J_per_m_sq, + J_per_ang_sq, + J_per_ang_sq_kg, + ] + + def __repr__(self): + return f"Force constant({round(self, 5)} {self.units.name})" + + def __init__(self, value, units=ha_per_ang_sq): + super().__init__(value, units=units) + + class ValueArray(ABC, np.ndarray): """ Abstract base class for an array of values, e.g. gradients or a Hessian @@ -782,8 +783,8 @@ def __repr__(self): def method_string( - method: Optional["autode.wrappers.methods.Method"], - keywords: Optional["autode.wrappers.keywords.Keywords"], + method: Optional["Method"], + keywords: Optional["Keywords"], ) -> str: """ Create a method string for a method and the keywords diff --git a/doc/changelog.rst b/doc/changelog.rst index eebce275e..c5fab9841 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -8,13 +8,19 @@ Changelog Functionality improvements ************************** - Adds :code:`temporary_config()` context manager for temporary configuration changes +- Adds a :code:`ForceConstant` value Usability improvements/Changes ****************************** - Adds full usability of autodE on Windows, including parallelisation with :code:`loky` - Optional timeout for graph isomorphism test in Windows, turned on by :code:`Config.use_experimental_timeout=True` (default behaviour kept for Linux/macOS) - +- A NEB :code:`Image` now derives from a :code:`Species` superclass +- Modifies NEB :code:`Image` constructor to be formed from an image +- Defines named constructors (:code:`from_endpoints(...)`, :code:`from_list(...)`) for :code:`NEB` +- Removes :code:`NEB().contains_peak()` in favour of :code:`NEB().images.contains_peak` +- Modifies the :code:`CImages` constructor to ensure it's constructed from an :code:`Images` instance +- Removes :code:`NEB.get_species_saddle_point()` in favour of :code:`NEB.peak_species` Bug Fixes ********* diff --git a/doc/common/claisen_cineb.py b/doc/common/claisen_cineb.py index 76d4b8b3f..443c4d25a 100644 --- a/doc/common/claisen_cineb.py +++ b/doc/common/claisen_cineb.py @@ -5,10 +5,10 @@ # Create an 8 image nudged elastic band with intermediate images interpolated # from the final points, thus they must be structurally similar -neb = ade.neb.CINEB(initial_species=reac, final_species=prod, num=8) +neb = ade.CINEB.from_end_points(reac, prod, num=8) # minimise with XTB neb.calculate(method=ade.methods.XTB(), n_cores=4) # print the geometry of the peak species -peak_species = neb.get_species_saddle_point() -peak_species.print_xyz_file(filename="peak.xyz") +print("Found a peak: ", neb.images.contains_peak) +neb.peak_species.print_xyz_file(filename="peak.xyz") diff --git a/doc/common/nci_FF_example.py b/doc/common/nci_FF_example.py index e966b30ba..0d865d652 100644 --- a/doc/common/nci_FF_example.py +++ b/doc/common/nci_FF_example.py @@ -1,8 +1,7 @@ import autode as ade -from autode.methods import XTB -from autode.calculations import Calculation -from scipy.optimize import minimize import numpy as np + +from scipy.optimize import minimize from scipy.spatial import distance_matrix ade.Config.max_num_complex_conformers = 10 @@ -85,17 +84,16 @@ def rotation_translation( def set_charges_vdw(species): """Calculate the partial atomic charges to atoms with XTB""" - calc = Calculation( + calc = ade.Calculation( name="tmp", molecule=species, - method=XTB(), - keywords=ade.SinglePointKeywords([]), + method=ade.methods.XTB(), + keywords=ade.SinglePointKeywords(), ) calc.run() - charges = calc.get_atomic_charges() for i, atom in enumerate(species.atoms): - atom.charge = charges[i] + atom.charge = atom.partial_charge atom.vdw = float(atom.vdw_radius) return None diff --git a/examples/tutorials/j_NEB.py b/examples/tutorials/j_NEB.py index 26d3e15a6..0fdca4105 100644 --- a/examples/tutorials/j_NEB.py +++ b/examples/tutorials/j_NEB.py @@ -13,14 +13,13 @@ # that support gradient evaluations (all of them!). For example, to # set up a set of images and relax to the ~minimum energy path for a # Diels Alder reaction between benzoquinone and cyclopentadiene -neb = ade.neb.NEB( - initial_species=ade.Molecule("_data/DielsAlder/reactant.xyz"), - final_species=ade.Molecule("_data/DielsAlder/product.xyz"), +neb = ade.NEB.from_end_points( + ade.Molecule("_data/DielsAlder/reactant.xyz"), + ade.Molecule("_data/DielsAlder/product.xyz"), num=5, ) - neb.calculate(method=orca, n_cores=10) # will have generated a plot of the relaxation, along with a .xyz # trajectory of the initial and final NEB path -# To use a climbing image NEB simply replace ade.neb.NEB with ade.neb.CINEB +# To use a climbing image NEB replace ade.NEB with ade.CINEB diff --git a/pyproject.toml b/pyproject.toml index 83979720a..13f1eb81f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,3 +10,9 @@ target-version = ['py39'] omit = [ "setup.py" ] + +[tool.coverage.report] +exclude_lines = [ + "pragma: no cover", + "if\\s+(typing\\.)?TYPE_CHECKING:" +] diff --git a/tests/requirements.txt b/tests/requirements.txt index 8d8087843..6e7a91e34 100644 --- a/tests/requirements.txt +++ b/tests/requirements.txt @@ -1,4 +1,4 @@ -coverage=6.2 +coverage pytest pytest-cov codecov \ No newline at end of file diff --git a/tests/test_calculation.py b/tests/test_calculation.py index bce7def83..b401013da 100644 --- a/tests/test_calculation.py +++ b/tests/test_calculation.py @@ -29,17 +29,19 @@ ) -test_mol = Molecule(smiles="O", name="test_mol") - - def h_atom() -> Molecule: return Molecule(atoms=[Atom("H")], mult=2) +def h2o() -> Molecule: + return Molecule(smiles="O", name="test_mol") + + @work_in_tmp_dir() def test_calc_class(): xtb = XTB() + test_mol = h2o() calc = Calculation( name="-tmp", molecule=test_mol, method=xtb, keywords=xtb.keywords.sp @@ -51,17 +53,10 @@ def test_calc_class(): assert calc.method.name == "xtb" assert len(calc.input.filenames) == 0 - with pytest.raises(ex.CouldNotGetProperty): - _ = calc.get_energy() - - with pytest.raises(ex.CouldNotGetProperty): - _ = calc.get_final_atoms() - - with pytest.raises(ex.CouldNotGetProperty): - _ = calc.get_gradients() - - with pytest.raises(ex.CouldNotGetProperty): - _ = calc.get_atomic_charges() + # Without + assert test_mol.energy is None + assert test_mol.gradient is None + assert test_mol.hessian is None # With a filename that doesn't exist a NoOutput exception should be raised calc.output.filename = "/a/path/that/does/not/exist/tmp" @@ -69,12 +64,12 @@ def test_calc_class(): _ = calc.output.file_lines # With no output should not be able to get properties - calc.output.filename = "tmp.out" - with open(calc.output.filename, "w") as output_file: + with open("tmp.out", "w") as output_file: print("some\ntest\noutput", file=output_file) - with pytest.raises(ex.CouldNotGetProperty): - _ = calc.get_atomic_charges() + with pytest.raises(ex.CalculationException): + # Cannot set even the energy from an invalid output file + calc.set_output_filename("tmp.out") # Should default to a single core assert calc.n_cores == 1 @@ -105,6 +100,8 @@ def test_calc_class(): def test_calc_copy(): orca = ORCA() + test_mol = h2o() + calc = Calculation( name="tmp", molecule=test_mol, method=orca, keywords=orca.keywords.sp ) @@ -160,6 +157,7 @@ def test_distance_const_check(): def test_calc_string(): xtb = XTB() + test_mol = h2o() a = test_mol.copy() no_const = Calculation( @@ -197,6 +195,7 @@ def test_fix_unique(): autodE checks the input of each previously run calc with the name name""" orca = ORCA() + test_mol = h2o() calc = CalculationExecutor( name="tmp", molecule=test_mol, method=orca, keywords=orca.keywords.sp @@ -225,7 +224,7 @@ def test_fix_unique(): def test_solvent_get(): xtb = XTB() - _test_mol = Molecule(smiles="O", name="test_mol") + _test_mol = h2o() # Can't get the name of a solvent if molecule.solvent is not a string with pytest.raises(ex.SolventUnavailable): @@ -240,7 +239,7 @@ def test_solvent_get(): # Currently iodoethane is not in XTB - might be in the future _test_mol.solvent = "iodoethane" - assert not hasattr(test_mol.solvent, "xtb") + assert not hasattr(_test_mol.solvent, "xtb") assert _test_mol.solvent.is_implicit with pytest.raises(ex.SolventUnavailable): @@ -253,6 +252,7 @@ def test_solvent_get(): def test_input_gen(): xtb = XTB() + test_mol = h2o() calc = Calculation( name="tmp", molecule=test_mol, method=xtb, keywords=xtb.keywords.sp @@ -286,6 +286,8 @@ def test_input_gen(): def test_exec_not_avail_method(): orca = ORCA() + test_mol = h2o() + orca.path = "/a/non/existent/path" assert not orca.is_available diff --git a/tests/test_constraints.py b/tests/test_constraints.py index a1a0f95fd..cd7f6b0f7 100644 --- a/tests/test_constraints.py +++ b/tests/test_constraints.py @@ -1,5 +1,5 @@ import pytest -from autode.constraints import Constraints +from autode.constraints import Constraints, DistanceConstraints def test_base_properties(): @@ -110,3 +110,12 @@ def test_clear(): consts.update(cartesian=[0]) assert len(consts.cartesian) == 1 + + +def test_copy(): + + constraints = DistanceConstraints({(0, 1): 1.0}) + copied_constraints = constraints.copy() + constraints[(0, 1)] = 1.1 + + assert abs(copied_constraints[(0, 1)] - 1.0) < 1e-10 diff --git a/tests/test_neb.py b/tests/test_neb.py index 4019664e8..e674091b2 100644 --- a/tests/test_neb.py +++ b/tests/test_neb.py @@ -4,7 +4,7 @@ import pytest from autode.path import Path from autode.neb import NEB, CINEB -from autode.values import Distance +from autode.values import Distance, ForceConstant from autode.neb.ci import Images, CImages, Image from autode.neb.idpp import IDPP from autode.species.molecule import Species, Molecule @@ -40,30 +40,35 @@ def test_neb_properties(): atoms=[Atom("H"), Atom("H", x=1.3), Atom("H", x=2.0)], ) - neb = NEB(initial_species=reac, final_species=prod, num=3) + neb = NEB.from_end_points(reac, prod, num=3) assert len(neb.images) == 3 - assert neb.get_species_saddle_point() is None - assert not neb.contains_peak() + assert neb.peak_species is None + assert not neb.images.contains_peak # Should move monotonically from 0.7 -> 1.3 Angstroms for i in range(1, len(neb.images)): - prev_bb_dist = neb.images[i - 1].species.distance(0, 1) - curr_bb_dist = neb.images[i].species.distance(0, 1) + prev_bb_dist = neb.images[i - 1].distance(0, 1) + curr_bb_dist = neb.images[i].distance(0, 1) assert curr_bb_dist > prev_bb_dist def test_image_properties(): - images = CImages(images=[Image(name="tmp", k=0.1)]) + k = ForceConstant(0.1) + images = CImages(images=Images(init_k=k)) assert images != 0 assert images == images - images = Images(init_k=0.1, num=1) + images = Images(init_k=k) assert images != 0 assert images == images + image = Image(Molecule(smiles="CC"), k=ForceConstant(1.0), name="tmp") + with pytest.raises(Exception): + image._generate_conformers() + def test_contains_peak(): @@ -90,15 +95,15 @@ def test_contains_peak(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "neb.zip")) def test_full_calc_with_xtb(): - sn2_neb = NEB( - initial_species=Species( + sn2_neb = NEB.from_end_points( + initial=Species( name="inital", charge=-1, mult=1, atoms=xyz_file_to_atoms("sn2_init.xyz"), solvent_name="water", ), - final_species=Species( + final=Species( name="final", charge=-1, mult=1, @@ -108,16 +113,10 @@ def test_full_calc_with_xtb(): num=14, ) - sn2_neb.interpolate_geometries() - - xtb = XTB() - - xtb.path = shutil.which("xtb") - sn2_neb.calculate(method=xtb, n_cores=2) + sn2_neb.calculate(method=XTB(), n_cores=2) # There should be a peak in this surface - peak_species = sn2_neb.get_species_saddle_point() - assert peak_species is not None + assert sn2_neb.peak_species is not None assert all(image.energy is not None for image in sn2_neb.images) @@ -174,7 +173,10 @@ def test_get_ts_guess_neb(): def test_climbing_image(): - images = CImages(images=[Image(name="tmp", k=0.1)]) + k = ForceConstant(0.1) + images = CImages(images=Images(init_k=k)) + images.append_species(Molecule(atoms=[Atom("H")], mult=2)) + assert images.peak_idx is None assert images[0].iteration == 0 images[0].iteration = 10 @@ -183,31 +185,35 @@ def test_climbing_image(): def _simple_h2_images(num, shift, increment): """Simple set of images for a n-image NEB for H2""" - images = Images(num=num, init_k=1.0) + images = Images(init_k=ForceConstant(1.0)) for i in range(num): - images[i].species = Molecule( - atoms=[Atom("H"), Atom("H", x=shift + i * increment)] - ) + mol = Molecule(atoms=[Atom("H"), Atom("H", x=shift + i * increment)]) + images.append_species(mol) return images def test_energy_gradient_type(): + k = ForceConstant(1.0) + image = Image(species=Molecule(atoms=[Atom("H")], mult=2), name="tmp", k=k) + # Energy and gradient must have a method (EST or IDPP) with pytest.raises(ValueError): - _ = energy_gradient(image=Image("tmp", k=1.0), method=None, n_cores=1) + _ = energy_gradient(image=image, method=None, n_cores=1) def test_iddp_init(): """IDPP requires at least 2 images""" + k = ForceConstant(0.1) + with pytest.raises(ValueError): - _ = IDPP(Images(num=0, init_k=0.1)) + _ = IDPP(Images(init_k=k)) with pytest.raises(ValueError): - _ = IDPP(Images(num=1, init_k=0.1)) + _ = IDPP(Images(init_k=k)) def test_iddp_energy(): @@ -236,7 +242,7 @@ def test_iddp_gradient(): value = idpp(image) # and the gradient calculable - grad = idpp.grad(image) + grad = idpp.grad(image).flatten() assert grad is not None # And the gradient be close to the numerical analogue @@ -246,9 +252,9 @@ def num_grad(n, h=1e-8): shift_vec = np.zeros(3) shift_vec[k] = h - image.species.atoms[i].translate(shift_vec) + image.atoms[i].translate(shift_vec) new_value = idpp(image) - image.species.atoms[i].translate(-shift_vec) + image.atoms[i].translate(-shift_vec) return (new_value - value) / h @@ -260,7 +266,7 @@ def num_grad(n, h=1e-8): assert np.isclose(analytic_value, num_grad(i), atol=1e-5) -@work_in_tmp_dir(filenames_to_copy=[], kept_file_exts=[]) +@work_in_tmp_dir() def test_neb_interpolate_and_idpp_relax(): mol = Molecule( @@ -277,12 +283,10 @@ def test_neb_interpolate_and_idpp_relax(): rot_mol = mol.copy() rot_mol.rotate(axis=[1.0, 0.0, 0.0], theta=1.5) - neb = NEB(initial_species=mol, final_species=rot_mol, num=10) - neb.interpolate_geometries() - neb.idpp_relax() + neb = NEB.from_end_points(initial=mol, final=rot_mol, num=10) for image in neb.images: - assert are_coords_reasonable(image.species.coordinates) + assert are_coords_reasonable(image.coordinates) def test_max_delta_between_images(): @@ -293,12 +297,12 @@ def test_max_delta_between_images(): ] assert np.isclose( - NEB(species_list=_list).max_atom_distance_between_images, 1.0 + NEB.from_list(_list).max_atom_distance_between_images, 1.0 ) _list[0].atoms[1].coord[0] = 1.7 # x coordinate of the second atom assert np.isclose( - NEB(species_list=_list).max_atom_distance_between_images, 0.0 + NEB.from_list(_list).max_atom_distance_between_images, 0.0 ) @@ -309,7 +313,7 @@ def test_max_delta_between_images_h3(): Molecule(atoms=[Atom("H"), Atom("H", x=0.70657), Atom("H", x=2.7)]), ] - neb = NEB(species_list=_list) + neb = NEB.from_list(_list) assert np.isclose(neb.max_atom_distance_between_images, 0.00657) assert np.isclose( @@ -327,7 +331,7 @@ def test_partition_max_delta(): Molecule(atoms=[Atom("H"), Atom("H", x=2.0), Atom("H", x=2.7)]), ] - h2_h = NEB(species_list=_list) + h2_h = NEB.from_list(_list) max_delta = Distance(0.1, units="Å") assert ( @@ -343,8 +347,7 @@ def test_partition_max_delta(): assert ( np.max( np.linalg.norm( - h2_h.images[i].species.coordinates - - h2_h.images[j].species.coordinates, + h2_h.images[i].coordinates - h2_h.images[j].coordinates, axis=1, ) ) @@ -373,7 +376,7 @@ def test_init_from_file_sets_force_constant(): ) # Should be able to set the initial force constant - neb = NEB.from_file("tmp.xyz", k=0.234) + neb = NEB.from_file("tmp.xyz", init_k=0.234) assert len(neb.images) == 3 assert np.isclose(neb.init_k, 0.234) @@ -404,3 +407,17 @@ def test_init_from_file_sets_force_constant_no_energies(): neb = NEB.from_file("tmp.xyz") # Estimated value should be reasonable even without energies assert 0.001 < neb.init_k < 0.2 + + +def test_neb_constructor_with_kwargs_raises(): + + with pytest.raises(Exception): + _ = NEB(init_k=ForceConstant(0.1), another_arg="a string") + + +def test_constructing_neb_from_endpoints_with_different_atoms_raises(): + + with pytest.raises(Exception): + _ = NEB.from_end_points( + Molecule(smiles="O"), Molecule(smiles="C"), num=4 + ) diff --git a/tests/test_path.py b/tests/test_path.py index d6fa2d64e..237beffa6 100644 --- a/tests/test_path.py +++ b/tests/test_path.py @@ -4,11 +4,12 @@ from autode.atoms import Atom from autode.methods import XTB -from autode.path import Path, AdaptivePath, PathPoint +from autode.path import Path, AdaptivePath from autode.path.adaptive import pruned_active_bonds from autode.bonds import FormingBond, BreakingBond from autode.species import Species, Molecule from autode.units import Unit, KcalMol +from . import testutils test_species = Species(name="tmp", charge=0, mult=1, atoms=[Atom("He")]) test_mol = Molecule(smiles="O") @@ -24,7 +25,7 @@ def test_path_properties_empty(): assert path == Path() # should be able to compare paths assert path != 0 - with pytest.raises(ValueError): + with pytest.raises(Exception): _ = Path("does not have correct attributes") # With no species there should be no peak/saddle/energies @@ -41,16 +42,16 @@ def test_path_properties_empty(): def test_path_properties(): - p1 = PathPoint(test_species.copy(), constraints={}) + p1 = test_species.copy() p1.energy = -3 - p2 = PathPoint(test_species.copy(), constraints={}) + p2 = test_species.copy() p2.energy = -2 path = Path(p1, p2, units=KcalMol) assert all(np.isclose(path.energies, np.array([-3, -2]))) assert all(np.isclose(path.rel_energies, 627.509 * np.array([0, 1]))) - p3 = PathPoint(test_species.copy(), constraints={}) + p3 = test_species.copy() path = Path(p1, p2, p3) # There is an energy not set, should not be able to find a peak @@ -77,20 +78,12 @@ def test_path_properties(): def test_point_properties(): - point = PathPoint(species=test_species.copy(), constraints={}) + point = test_species.copy() assert point.energy is None - assert point.grad is None - assert not point.species.constraints.any - assert point.species.name == "tmp" - - point.energy = 1 - point.grad = np.zeros(shape=(1, 3)) - - # Copies should not copy energy or gradients - new_point = point.copy() - assert new_point.energy is None - assert new_point.grad is None + assert point.gradient is None + assert not point.constraints.any + assert point.name == "tmp" def test_pruning_bonds(): @@ -172,7 +165,7 @@ def test_pruning_bonds2(): def test_products_made(): - path = Path(PathPoint(species=test_mol, constraints={})) + path = Path(test_mol) assert not path.products_made(product=None) # Species have no graphs @@ -187,11 +180,20 @@ def test_products_made(): assert not path.products_made(product=diff_mol) +@testutils.requires_with_working_xtb_install def test_adaptive_path(): species_no_atoms = Species(name="tmp", charge=0, mult=1, atoms=[]) - path1 = AdaptivePath(init_species=species_no_atoms, bonds=[], method=XTB()) - assert len(path1) == 0 + + with pytest.raises(Exception): + # cannot create a path with a molecule with no atoms + _ = AdaptivePath(init_species=species_no_atoms, bonds=[], method=XTB()) + + path1 = AdaptivePath( + init_species=Molecule(smiles="O"), bonds=[], method=XTB() + ) + + assert len(path1) == 1 assert path1.method.name == "xtb" assert len(path1.bonds) == 0 diff --git a/tests/test_pes/test_calculate.py b/tests/test_pes/test_calculate.py index 94a3a8d0f..3f2f3ee4d 100644 --- a/tests/test_pes/test_calculate.py +++ b/tests/test_pes/test_calculate.py @@ -24,7 +24,7 @@ def test_calculate_no_species(): @testutils.requires_with_working_xtb_install @work_in_tmp_dir(filenames_to_copy=[], kept_file_exts=[]) def test_calculate_1d(): - pes = RelaxedPESnD(species=h2(), rs={(0, 1): (1.7, 10)}) + pes = RelaxedPESnD(species=h2(), rs={(0, 1): (1.5, 10)}) pes.calculate(method=XTB()) diff --git a/tests/test_qrc.py b/tests/test_qrc.py index 65ec0f792..ae1787b13 100644 --- a/tests/test_qrc.py +++ b/tests/test_qrc.py @@ -14,15 +14,14 @@ def test_hshift_displacement(): orca = ORCA() - reac = Reactant(name="reactant", smiles="CC[C]([H])[H]") + ts = Reactant(name="reactant", smiles="CC[C]([H])[H]") calc = Calculation( - name="tmp", molecule=reac, method=orca, keywords=orca.keywords.opt_ts + name="tmp", molecule=ts, method=orca, keywords=orca.keywords.opt_ts ) calc.set_output_filename("TS_hshift.out") assert calc.terminated_normally - ts = Molecule(atoms=calc.get_final_atoms()) - ts.hessian = calc.get_hessian() + assert ts.hessian is not None f_disp_mol = displaced_species_along_mode( ts, mode_number=6, disp_factor=1.0 # TS mode diff --git a/tests/test_ts/test_mode_checking.py b/tests/test_ts/test_mode_checking.py index 30899643c..199386287 100644 --- a/tests/test_ts/test_mode_checking.py +++ b/tests/test_ts/test_mode_checking.py @@ -45,28 +45,24 @@ def test_imag_modes(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "mode_checking.zip")) def test_graph_no_other_bonds(): - reac = Reactant( - name="r", + ts = TSbase( atoms=xyz_file_to_atoms("h_shift_correct_ts_mode.xyz"), mult=2, + bond_rearr=BondRearrangement( + breaking_bonds=[(1, 10)], forming_bonds=[(5, 10)] + ), ) calc = Calculation( name="h_shift", - molecule=reac, + molecule=ts, method=orca, keywords=orca.keywords.opt_ts, n_cores=1, ) calc.set_output_filename("h_shift_correct_ts_mode.out") - ts = TSbase( - atoms=calc.get_final_atoms(), - bond_rearr=BondRearrangement( - breaking_bonds=[(1, 10)], forming_bonds=[(5, 10)] - ), - ) - ts.hessian = calc.get_hessian() + assert ts.hessian is not None f_ts = displaced_species_along_mode(ts, mode_number=6, disp_factor=1.0) b_ts = displaced_species_along_mode(ts, mode_number=6, disp_factor=-1.0) @@ -85,7 +81,8 @@ def has_correct_mode(name, fbonds, bbonds): keywords=orca.keywords.opt_ts, n_cores=1, ) - calc.molecule = reac = Reactant( + # need to bypass the pre-calculation checks on the molecule. e.g. valid spin state + calc.molecule = reactant = Reactant( name="r", atoms=xyz_file_to_atoms(f"{name}.xyz") ) @@ -93,11 +90,11 @@ def has_correct_mode(name, fbonds, bbonds): # Don't require all bonds to be breaking/making in a 'could be ts' function ts = TSbase( - atoms=reac.atoms, + atoms=reactant.atoms, bond_rearr=BondRearrangement( breaking_bonds=bbonds, forming_bonds=fbonds ), ) - ts.hessian = calc.get_hessian() + ts.hessian = reactant.hessian return ts.imag_mode_has_correct_displacement(req_all=False) diff --git a/tests/test_ts/test_ts_adapt_neb.py b/tests/test_ts/test_ts_adapt_neb.py index bff7319ab..448d92edb 100644 --- a/tests/test_ts/test_ts_adapt_neb.py +++ b/tests/test_ts/test_ts_adapt_neb.py @@ -43,8 +43,9 @@ def test_ts_from_neb_optimised_after_adapt(): rxn = _sn2_reaction() neb = NEB.from_file("dOr2Us_ll_ad_0-5_0-1_path.xyz") + assert len(neb.images) > 0 - assert neb.images[0].species.solvent is not None + assert neb.images[0].solvent is not None def get_ts_guess(): return _get_ts_neb_from_adaptive_path( diff --git a/tests/test_values.py b/tests/test_values.py index dd8059111..92e47bc9b 100644 --- a/tests/test_values.py +++ b/tests/test_values.py @@ -8,6 +8,7 @@ Coordinate, Coordinates, MomentOfInertia, + ForceConstant, _to, ) @@ -125,3 +126,12 @@ def test_copy_conversion(): assert not np.allclose(x, y) assert np.allclose(x, np.ones(shape=(1, 3))) + + +def test_force_constant(): + + fc = ForceConstant(0.1) + assert "force" in repr(fc).lower() + + # should be able to convert to Ha/a0^2 without any problems + _ = fc.to("Ha/a0^2") diff --git a/tests/test_wrappers/test_gaussian.py b/tests/test_wrappers/test_gaussian.py index 1cb184a05..3b96756de 100644 --- a/tests/test_wrappers/test_gaussian.py +++ b/tests/test_wrappers/test_gaussian.py @@ -23,7 +23,6 @@ from .. import testutils here = os.path.dirname(os.path.abspath(__file__)) -test_mol = Molecule(name="methane", smiles="C") method = G09() opt_keywords = OptKeywords(["PBE1PBE/Def2SVP", "Opt"]) @@ -39,6 +38,10 @@ sp_keywords = SinglePointKeywords(["PBE1PBE/Def2SVP"]) +def methane(): + return Molecule(name="methane", smiles="C") + + def test_printing_ecp(): tmp_file = open("tmp.com", "w") @@ -90,7 +93,7 @@ def test_input_print_max_opt(): keywds = opt_keywords.copy() keywds.max_opt_cycles = 10 - str_keywords = _get_keywords(CalculationInput(keywds), molecule=test_mol) + str_keywords = _get_keywords(CalculationInput(keywds), molecule=methane()) # Should be only a single instance of the maxcycles declaration assert sum("maxcycles=10" in kw.lower() for kw in str_keywords) == 1 @@ -145,8 +148,8 @@ def test_gauss_opt_calc(): assert os.path.exists("opt_g09.com") assert os.path.exists("opt_g09.log") - assert len(calc.get_final_atoms()) == 5 - assert calc.get_energy() == -499.729222331 + assert len(methylchloride.atoms) == 5 + assert methylchloride.energy.to("Ha") == -499.729222331 assert calc.output.exists assert calc.output.file_lines is not None assert methylchloride.imaginary_frequencies is None @@ -154,16 +157,15 @@ def test_gauss_opt_calc(): assert calc.input.filename == "opt_g09.com" assert calc.output.filename == "opt_g09.log" assert calc.terminated_normally - assert calc.optimisation_converged() - assert calc.optimisation_nearly_converged() is False + assert calc.optimiser.converged - charges = calc.get_atomic_charges() + charges = methylchloride.partial_charges assert len(charges) == methylchloride.n_atoms # Should be no very large atomic charges in this molecule assert all(-1.0 < c < 1.0 for c in charges) - gradients = calc.get_gradients() + gradients = methylchloride.gradient assert len(gradients) == methylchloride.n_atoms assert len(gradients[0]) == 3 @@ -202,23 +204,21 @@ def test_gauss_optts_calc(): assert bond_added - mol = Molecule(atoms=calc.get_final_atoms()) - mol.calc_thermo(calc=calc, ss="1atm", lfm_method="igm") + test_mol.calc_thermo(calc=calc, ss="1atm", lfm_method="igm") assert calc.terminated_normally - assert calc.optimisation_converged() - assert calc.optimisation_nearly_converged() is False - assert len(mol.imaginary_frequencies) == 1 + assert calc.optimiser.converged + assert len(test_mol.imaginary_frequencies) == 1 - assert -40.324 < mol.free_energy < -40.322 - assert -40.301 < mol.enthalpy < -40.298 + assert -40.324 < test_mol.free_energy < -40.322 + assert -40.301 < test_mol.enthalpy < -40.298 def test_bad_gauss_output(): calc = Calculation( name="no_output", - molecule=test_mol, + molecule=methane(), method=method, keywords=opt_keywords, ) @@ -226,10 +226,7 @@ def test_bad_gauss_output(): calc.rev_output_file_lines = [] with pytest.raises(CalculationException): - _ = calc.get_energy() - - with pytest.raises(CalculationException): - calc.get_final_atoms() + calc.set_output_filename("no_output") @testutils.work_in_zipped_dir(os.path.join(here, "data", "g09.zip")) @@ -254,26 +251,26 @@ def test_fix_angle_error(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "g09.zip")) def test_constraints(): - a = test_mol.copy() + a = methane() a.constraints.distance = {(0, 1): 1.2} calc = Calculation( name="const_dist_opt", molecule=a, method=method, keywords=opt_keywords ) calc.run() - opt_atoms = calc.get_final_atoms() + opt_atoms = a.atoms assert ( 1.199 < np.linalg.norm(opt_atoms[0].coord - opt_atoms[1].coord) < 1.201 ) - b = test_mol.copy() + b = methane() b.constraints.cartesian = [0] calc = Calculation( name="const_cart_opt", molecule=b, method=method, keywords=opt_keywords ) calc.run() - opt_atoms = calc.get_final_atoms() - assert np.linalg.norm(test_mol.atoms[0].coord - opt_atoms[0].coord) < 1e-3 + opt_atoms = b.atoms + assert np.linalg.norm(methane().atoms[0].coord - opt_atoms[0].coord) < 1e-3 @testutils.work_in_zipped_dir(os.path.join(here, "data", "g09.zip")) @@ -305,9 +302,10 @@ def test_point_charge_calc(): # Methane single point using a point charge with a unit positive charge # located at (10, 10, 10) + mol = methane() calc = Calculation( name="methane_point_charge", - molecule=test_mol, + molecule=mol, method=method, keywords=sp_keywords, point_charges=[PointCharge(charge=1.0, x=10.0, y=10.0, z=10.0)], @@ -329,14 +327,14 @@ def test_point_charge_calc(): assert float(z) == 10.0 assert float(charge) == 1.0 - assert -40.428 < calc.get_energy() < -40.427 + assert -40.428 < mol.energy < -40.427 # Gaussian needs x-matrix and nosymm in the input line to run optimisations # with point charges.. for opt_keyword in ["Opt", "Opt=Tight", "Opt=(Tight)"]: calc = Calculation( name="methane_point_charge_o", - molecule=test_mol, + molecule=methane(), method=method, keywords=OptKeywords(["PBE1PBE/Def2SVP", opt_keyword]), point_charges=[PointCharge(charge=1.0, x=3.0, y=3.0, z=3.0)], @@ -423,5 +421,5 @@ def test_xtb_optts(): calc.run() # Even though a Hessian is not requested it should be added - assert calc.get_hessian() is not None - assert np.isclose(calc.get_energy().to("Ha"), -13.1297380, atol=1e-5) + assert orca_ts.hessian is not None + assert np.isclose(orca_ts.energy.to("Ha"), -13.1297380, atol=1e-5) diff --git a/tests/test_wrappers/test_mopac.py b/tests/test_wrappers/test_mopac.py index 8f126522f..e4cc07d68 100644 --- a/tests/test_wrappers/test_mopac.py +++ b/tests/test_wrappers/test_mopac.py @@ -42,11 +42,11 @@ def test_mopac_opt_calculation(): assert os.path.exists("opt_mopac.mop") is True assert os.path.exists("opt_mopac.out") is True - assert len(calc.get_final_atoms()) == 5 + assert mol.n_atoms == 5 # Actual energy in Hartrees energy = Constants.eV_to_ha * -430.43191 - assert energy - 0.0001 < calc.get_energy() < energy + 0.0001 + assert energy - 0.0001 < mol.energy < energy + 0.0001 assert calc.output.exists assert calc.output.file_lines is not None @@ -60,9 +60,10 @@ def test_mopac_opt_calculation(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "mopac.zip")) def test_mopac_with_pc(): + mol = mecl() calc = Calculation( name="opt_pc", - molecule=mecl(), + molecule=mol, method=method, keywords=Config.MOPAC.keywords.opt, point_charges=[PointCharge(1, x=4, y=4, z=4)], @@ -71,11 +72,11 @@ def test_mopac_with_pc(): assert os.path.exists("opt_pc_mopac.mop") is True assert os.path.exists("opt_pc_mopac.out") is True - assert len(calc.get_final_atoms()) == 5 + assert len(mol.atoms) == 5 # Actual energy in Hartrees without any point charges energy = Constants.eV_to_ha * -430.43191 - assert np.abs(calc.get_energy() - energy) > 0.0001 + assert np.abs(mol.energy - energy) > 0.0001 @testutils.work_in_zipped_dir(os.path.join(here, "data", "mopac.zip")) @@ -91,7 +92,6 @@ def test_other_spin_states(): keywords=Config.MOPAC.keywords.sp, ) calc.run() - singlet_energy = calc.get_energy() o_triplet = Molecule(atoms=[Atom("O")], mult=3) o_triplet.name = "molecule" @@ -103,9 +103,8 @@ def test_other_spin_states(): keywords=Config.MOPAC.keywords.sp, ) calc.run() - triplet_energy = calc.get_energy() - assert triplet_energy < singlet_energy + assert o_triplet.energy < o_singlet.energy h_doublet = Molecule(atoms=[Atom("H")], mult=2) h_doublet.name = "molecule" @@ -119,7 +118,7 @@ def test_other_spin_states(): calc.run() # Open shell doublet should work - assert calc.get_energy() is not None + assert h_doublet.energy is not None h_quin = Molecule(atoms=[Atom("H")], mult=5) h_quin.name = "molecule" @@ -147,11 +146,9 @@ def test_bad_geometry(): keywords=Config.MOPAC.keywords.opt, ) - calc.output.filename = "h2_overlap_opt_mopac.out" - assert not calc.terminated_normally - with pytest.raises(Exception): - _ = calc.get_energy() + # cannot even get the energy from the output file + calc.set_output_filename("h2_overlap_opt_mopac.out") assert not method.optimiser_from(calc).converged @@ -168,21 +165,23 @@ def test_constrained_opt(): keywords=Config.MOPAC.keywords.opt, ) calc.run() - opt_energy = calc.get_energy() + opt_energy = methane.energy # Constrained optimisation with a C–H distance of 1.2 Å # (carbon is the first atom in the file) + constrained_methane = methane.copy() methane.constraints.distance = {(0, 1): 1.2} + const = Calculation( name="methane_const", - molecule=methane, + molecule=constrained_methane, method=method, keywords=Config.MOPAC.keywords.opt, ) const.run() - assert opt_energy < const.get_energy() - assert calc.get_hessian() is None + assert methane.energy < constrained_methane.energy + assert methane.hessian is None @testutils.work_in_zipped_dir(os.path.join(here, "data", "mopac.zip")) @@ -197,10 +196,9 @@ def test_grad(): keywords=Config.MOPAC.keywords.grad, ) grad_calc.run() - energy = grad_calc.get_energy() - assert energy is not None + assert h2.energy is not None - gradients = grad_calc.get_gradients() + gradients = h2.gradient assert gradients.shape == (2, 3) delta_r = 1e-5 @@ -209,7 +207,7 @@ def test_grad(): ) h2_disp.single_point(method) - delta_energy = h2_disp.energy - energy # Ha] + delta_energy = h2_disp.energy - h2.energy grad = delta_energy / delta_r # Ha A^-1 # Difference between the absolute and finite difference approximation @@ -226,10 +224,9 @@ def test_broken_grad(): method=method, keywords=Config.MOPAC.keywords.grad, ) - grad_calc_broken.output.filename = "h2_grad_broken.out" with pytest.raises(CouldNotGetProperty): - _ = grad_calc_broken.get_gradients() + grad_calc_broken.set_output_filename("h2_grad_broken.out") @testutils.work_in_zipped_dir(os.path.join(here, "data", "mopac.zip")) diff --git a/tests/test_wrappers/test_nwchem.py b/tests/test_wrappers/test_nwchem.py index 0e1bfb2bf..a56283c1b 100644 --- a/tests/test_wrappers/test_nwchem.py +++ b/tests/test_wrappers/test_nwchem.py @@ -1,4 +1,7 @@ -from autode import utils +import numpy as np +import pytest +import os + from autode.wrappers.NWChem import NWChem, ecp_block from autode.point_charges import PointCharge from autode.calculations import Calculation @@ -11,12 +14,9 @@ from autode.config import Config from autode.atoms import Atom from .. import testutils -import numpy as np -import pytest -import os + here = os.path.dirname(os.path.abspath(__file__)) -test_mol = Molecule(name="methane", smiles="C") method = NWChem() method.path = here # spoof install @@ -32,24 +32,24 @@ @testutils.work_in_zipped_dir(os.path.join(here, "data", "nwchem.zip")) def test_opt_calc(): + test_mol = Molecule(name="methane", smiles="C") + calc = Calculation( name="opt", molecule=test_mol, method=NWChem(), keywords=opt_keywords ) calc.run() - final_atoms = calc.get_final_atoms() - assert len(final_atoms) == 5 - assert type(final_atoms[0]) is Atom - assert -40.4165 < calc.get_energy() < -40.4164 + assert len(test_mol.atoms) == 5 + assert type(test_mol.atoms[0]) is Atom + assert -40.4165 < test_mol.energy < -40.4164 assert calc.terminated_normally - assert calc.optimisation_converged() - assert calc.optimisation_nearly_converged() is False + assert calc.optimiser.converged # No Hessian is computed for an optimisation calculation - assert calc.get_hessian() is None + assert test_mol.hessian is None # Optimisation should result in small gradients - gradients = calc.get_gradients() + gradients = test_mol.gradient assert len(gradients) == 5 assert all(-0.1 < np.linalg.norm(g) < 0.1 for g in gradients) @@ -154,7 +154,7 @@ def test_hessian_extract_butane(): ) calc.set_output_filename("butane_hess_nwchem.out") - hess = calc.get_hessian() + hess = butane.hessian assert hess is not None # bottom right corner element should be positive @@ -218,11 +218,11 @@ def test_point_charge_calculation(): ) calc.run() - assert calc.get_energy() is not None + assert h.energy is not None # H atom energy with a point charge should be different from the # isolated atoms HF energy - assert not np.isclose(calc.get_energy(), -0.5, atol=0.001) + assert not np.isclose(h.energy.to("Ha"), -0.5, atol=0.001) @testutils.work_in_zipped_dir(os.path.join(here, "data", "nwchem.zip")) diff --git a/tests/test_wrappers/test_orca.py b/tests/test_wrappers/test_orca.py index 3c141ada9..2b5734ccb 100644 --- a/tests/test_wrappers/test_orca.py +++ b/tests/test_wrappers/test_orca.py @@ -53,20 +53,18 @@ def test_orca_opt_calculation(): assert os.path.exists("opt_orca.inp") is True assert os.path.exists("opt_orca.out") is True - assert len(calc.get_final_atoms()) == 5 - assert -499.735 < calc.get_energy() < -499.730 + assert len(methylchloride.atoms) == 5 + assert -499.735 < methylchloride.energy < -499.730 assert calc.output.exists assert calc.output.file_lines is not None assert calc.input.filename == "opt_orca.inp" assert calc.output.filename == "opt_orca.out" assert calc.terminated_normally - assert calc.optimisation_converged() - - assert calc.optimisation_nearly_converged() is False + assert calc.optimiser.converged # Should have a partial atomic charge for every atom - charges = calc.get_atomic_charges() + charges = methylchloride.partial_charges assert charges == [-0.006954, -0.147352, 0.052983, 0.052943, 0.053457] calc = Calculation( @@ -113,12 +111,11 @@ def test_orca_optts_calculation(): assert ts.normal_mode(mode_number=6) is not None assert calc.terminated_normally - assert calc.optimisation_converged() - assert calc.optimisation_nearly_converged() is False + assert calc.optimiser.converged assert len(ts.imaginary_frequencies) == 1 # Gradients should be an n_atom x 3 array - gradients = calc.get_gradients() + gradients = ts.gradient assert gradients.shape == (ts.n_atoms, 3) assert -599.437 < ts.enthalpy < -599.436 @@ -135,10 +132,7 @@ def test_bad_orca_output(): ) with pytest.raises(ex.CouldNotGetProperty): - _ = calc.get_energy() - - with pytest.raises(ex.CouldNotGetProperty): - _ = calc.get_final_atoms() + calc.set_output_filename("no_output") calc.output_file_lines = None assert calc.terminated_normally is False @@ -217,7 +211,6 @@ def test_gradients(): keywords=method.keywords.grad, ) calc.run() - h2.energy = calc.get_energy() delta_r = 1e-8 @@ -232,7 +225,6 @@ def test_gradients(): keywords=method.keywords.grad, ) calc.run() - h2_disp.energy = calc.get_energy() delta_energy = h2_disp.energy - h2.energy # Ha grad = delta_energy / delta_r # Ha A^-1 @@ -246,7 +238,7 @@ def test_gradients(): calc.run() - diff = calc.get_gradients()[1, 0] - grad # Ha A^-1 + diff = h2.gradient[1, 0] - grad # Ha A^-1 # Difference between the absolute and finite difference approximation assert np.abs(diff) < 1e-3 @@ -255,15 +247,16 @@ def test_gradients(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "orca.zip")) def test_mp2_numerical_gradients(): + mol = Molecule("tmp_orca.xyz", charge=-1) calc = Calculation( name="tmp", - molecule=Molecule(atoms=xyz_file_to_atoms("tmp_orca.xyz"), charge=-1), + molecule=mol, method=method, keywords=method.keywords.grad, ) calc.set_output_filename(filename="tmp_orca.out") - gradients = calc.get_gradients() + gradients = mol.gradient assert len(gradients) == 6 expected = ( np.array([-0.00971201, -0.00773534, -0.02473580]) / Constants.a0_to_ang @@ -274,7 +267,7 @@ def test_mp2_numerical_gradients(): calc.set_output_filename(filename="numerical_orca.out") assert calc.output.filename == "numerical_orca.out" - gradients = calc.get_gradients() + gradients = mol.gradient assert len(gradients) == 6 expected = ( np.array([0.012397372, 0.071726232, -0.070942743]) @@ -331,15 +324,14 @@ def test_keyword_setting(): @testutils.work_in_zipped_dir(os.path.join(here, "data", "orca.zip")) def test_hessian_extraction(): + h2o = Molecule(smiles="O") calc = Calculation( name="tmp", - molecule=Molecule(smiles="O"), + molecule=h2o, method=method, keywords=method.keywords.hess, ) - calc.output.filename = "H2O_hess_orca.out" - with open("H2O_hess_orca.xyz", "w") as xyz_file: print( "3\n", @@ -350,7 +342,9 @@ def test_hessian_extraction(): file=xyz_file, ) - hessian = calc.get_hessian() + calc.set_output_filename("H2O_hess_orca.out") + + hessian = h2o.hessian assert hessian.shape == (9, 9) # should not have any very large values assert np.sum(np.abs(hessian)) < 100 @@ -372,12 +366,10 @@ def test_charges_from_v5_output_file(): method=method, keywords=method.keywords.sp, ) - calc.output.filename = "h2o_orca_v5_charges.out" - + calc.set_output_filename("h2o_orca_v5_charges.out") assert calc.output.exists - - # q_O q_H q_H - assert calc.get_atomic_charges() == [-0.313189, 0.156594, 0.156594] + # q_O q_H q_H + assert water.partial_charges == [-0.313189, 0.156594, 0.156594] def test_unsupported_freq_scaling(): diff --git a/tests/test_wrappers/test_qchem.py b/tests/test_wrappers/test_qchem.py index ef15abbe0..176c7b1bf 100644 --- a/tests/test_wrappers/test_qchem.py +++ b/tests/test_wrappers/test_qchem.py @@ -35,7 +35,7 @@ def _blank_calc(name="test"): def _completed_thf_calc(): calc = _blank_calc() - calc.output.filename = "smd_thf.out" + calc.set_output_filename("smd_thf.out") assert calc.output.exists assert len(calc.output.file_lines) > 0 @@ -204,9 +204,8 @@ def test_simple_input_generation(): def test_energy_extraction(): calc = _completed_thf_calc() - energy = calc.get_energy() - assert np.isclose(energy.to("Ha"), -232.45463628, atol=1e-8) + assert np.isclose(calc.molecule.energy.to("Ha"), -232.45463628, atol=1e-8) for calc in (_blank_calc(), _broken_output_calc(), _broken_output_calc2()): with pytest.raises(CalculationException): @@ -294,9 +293,11 @@ def test_h2o_opt(): @work_in_zipped_dir(qchem_data_zip_path) def test_gradient_extraction_h2o(): + + h2o = Molecule(smiles="O") calc = Calculation( name="test", - molecule=Molecule(smiles="O"), + molecule=h2o, method=method, keywords=OptKeywords(), ) @@ -305,15 +306,13 @@ def test_gradient_extraction_h2o(): assert calc.output.exists - grad = calc.get_gradients() - assert grad.shape == (3, 3) + assert h2o.gradient.shape == (3, 3) # The minimum should have a gradient close to zero - assert np.allclose(grad, np.zeros(shape=(3, 3)), atol=1e-4) + assert np.allclose(h2o.gradient, np.zeros(shape=(3, 3)), atol=1e-4) # also for this calculation the optimisation has converged assert calc.optimiser.converged - assert not calc.optimisation_nearly_converged() @work_in_zipped_dir(qchem_data_zip_path) @@ -321,10 +320,9 @@ def test_gradient_extraction_h2(): calc = _blank_calc() calc.molecule = Molecule(atoms=[Atom("H"), Atom("H", x=0.77)]) - calc.output.filename = "H2_qchem.out" + calc.set_output_filename("H2_qchem.out") - grad = calc.get_gradients() - assert grad.shape == (2, 3) + assert calc.molecule.gradient.shape == (2, 3) @work_in_zipped_dir(qchem_data_zip_path) @@ -343,10 +341,11 @@ def test_butane_gradient_extraction(): @work_in_zipped_dir(qchem_data_zip_path) def test_h2o_hessian_extraction(): + h2o = Molecule(smiles="O") calc = _blank_calc() calc.input.keywords = method.keywords.hess - calc.output.filename = "H2O_hess_qchem.out" - calc.molecule = Molecule(smiles="O") + calc.molecule = h2o + calc.set_output_filename("H2O_hess_qchem.out") hess = method.hessian_from(calc) assert hess.shape == (9, 9) @@ -358,7 +357,7 @@ def test_h2o_hessian_extraction(): # Final atoms are available, as the same ones input assert np.allclose( calc.molecule.coordinates, - calc.get_final_atoms().coordinates, + h2o.coordinates, atol=1e-8, ) @@ -446,9 +445,7 @@ def test_ts_opt(): # Should skip calculation for already completed and saved calculation calc.run() - assert np.isclose(calc.get_energy().to("Ha"), -599.4788133790, atol=1e-8) - - ts_mol.hessian = calc.get_hessian() + assert np.isclose(ts_mol.energy.to("Ha"), -599.4788133790, atol=1e-8) assert ts_mol.hessian is not None assert sum(freq.is_imaginary for freq in ts_mol.vib_frequencies) == 1