diff --git a/src/kimmdy/analysis.py b/src/kimmdy/analysis.py index aa597e81..4e68277f 100644 --- a/src/kimmdy/analysis.py +++ b/src/kimmdy/analysis.py @@ -396,9 +396,9 @@ def radical_migration( picked_recipes = {} for recipes in run_dir.glob("*decide_recipe/recipes.csv"): task_nr = int(recipes.parents[0].stem.split(sep="_")[0]) - rc, picked_recipe = RecipeCollection.from_csv(recipes) + _, picked_recipe = RecipeCollection.from_csv(recipes) picked_recipes[task_nr] = picked_recipe - sorted_recipes = [val for key, val in sorted(picked_recipes.items())] + sorted_recipes = [v for _, v in sorted(picked_recipes.items())] for sorted_recipe in sorted_recipes: connectivity_difference = {} @@ -598,7 +598,7 @@ def reaction_participation(dir: str, open_plot: bool = False): reaction_count = {"overall": 0} for recipes in run_dir.glob("*decide_recipe/recipes.csv"): # get picked recipe - rc, picked_rp = RecipeCollection.from_csv(recipes) + _, picked_rp = RecipeCollection.from_csv(recipes) assert picked_rp, f"No picked recipe found in {recipes}." # get involved atoms reaction_atom_ids = set() diff --git a/src/kimmdy/config.py b/src/kimmdy/config.py index 8f059494..6d44307d 100644 --- a/src/kimmdy/config.py +++ b/src/kimmdy/config.py @@ -328,9 +328,10 @@ def _validate(self, section: str = "config"): ) name = self.out.name.split("_") out_end = name[-1] - out_start = "_".join(name[:-1]) if out_end.isdigit(): - self.out = self.out.with_name(f"{out_start}_{int(out_end)+1:03}") + self.out = self.out.with_name( + f"{'_'.join(name[:-1])}_{int(out_end)+1:03}" + ) else: self.out = self.out.with_name(self.out.name + "_001") diff --git a/src/kimmdy/coordinates.py b/src/kimmdy/coordinates.py index 7720f0d8..be101a19 100644 --- a/src/kimmdy/coordinates.py +++ b/src/kimmdy/coordinates.py @@ -195,8 +195,12 @@ def merge_dihedrals( parameterizedB = None # construct parameterized Dihedral - if parameterizedA and parameterizedB: + if parameterizedA is not None and parameterizedB is not None: # same + + assert type(parameterizedA) == Dihedral or type(parameterizedA) == DihedralType + assert type(parameterizedB) == Dihedral or type(parameterizedB) == DihedralType + dihedralmerge = Dihedral( *dihedral_key, funct=funct, @@ -207,8 +211,9 @@ def merge_dihedrals( c4=parameterizedB.c1, c5=parameterizedB.periodicity, ) - elif parameterizedA: + elif parameterizedA is not None: # breaking + assert type(parameterizedA) == Dihedral or type(parameterizedA) == DihedralType dihedralmerge = Dihedral( *dihedral_key, funct=funct, @@ -219,8 +224,9 @@ def merge_dihedrals( c4="0.00", c5=parameterizedA.periodicity, ) - elif parameterizedB: + elif parameterizedB is not None: # binding + assert type(parameterizedB) == Dihedral or type(parameterizedB) == DihedralType dihedralmerge = Dihedral( *dihedral_key, funct="9", @@ -242,14 +248,10 @@ def merge_top_moleculetypes_slow_growth( molA: MoleculeType, molB: MoleculeType, ff: FF, - focus_nr: Optional[list[str]] = None, ) -> MoleculeType: """Takes two Topologies and joins them for a smooth free-energy like parameter transition simulation""" hyperparameters = {"morse_well_depth": 300} # [kJ mol-1] - # TODO: - # think about how to bring focus_nr into this - # atoms for nr in molA.atoms.keys(): atomA = molA.atoms[nr] @@ -520,9 +522,7 @@ def merge_top_moleculetypes_slow_growth( return molB -def merge_top_slow_growth( - topA: Topology, topB: Topology, focus_nr: Optional[list[str]] = None -) -> Topology: +def merge_top_slow_growth(topA: Topology, topB: Topology) -> Topology: """Takes two Topologies and joins them for a smooth free-energy like parameter transition simulation. @@ -531,7 +531,7 @@ def merge_top_slow_growth( molA = topA.moleculetypes[REACTIVE_MOLECULEYPE] molB = topB.moleculetypes[REACTIVE_MOLECULEYPE] - molB = merge_top_moleculetypes_slow_growth(molA, molB, topB.ff, focus_nr) + molB = merge_top_moleculetypes_slow_growth(molA, molB, topB.ff) return topB diff --git a/src/kimmdy/kmc.py b/src/kimmdy/kmc.py index c1f7d7b1..2921e04c 100644 --- a/src/kimmdy/kmc.py +++ b/src/kimmdy/kmc.py @@ -235,6 +235,10 @@ def extrande_mod( f"\t\texpected tau: {expected_tau}" ) + b = None + tau = None + l = None + t_max = None for cr_boarders, cr_rates, cr_recipes in zip( pairwise(boarders), rate_windows, recipe_windows ): @@ -262,6 +266,9 @@ def extrande_mod( if chosen_recipe is not None: break + if None in [b, tau, l, t_max]: + logger.error(f"Extrande calculation failed, some variables are None.") + if chosen_recipe is None: logger.info( f"No reaction was chosen\naccepted: 0, rejected: {rejected}, extra: {n_extra}" @@ -343,6 +350,8 @@ def extrande( ) n_extra = 0 + tau = None + l = None while t < t_max: crr_window_idx = np.searchsorted(boarders, t, side="right") - 1 # only 0 rates left -> skip to end @@ -390,6 +399,9 @@ def extrande( f"{n_extra} extra reactions performed during extrande KMC calculation. Try increasing tau_scale" ) + if None in [tau, l]: + logger.error(f"Extrande calculation failed, some variables are None.") + if chosen_recipe is None: logger.info( f"No reaction was chosen\naccepted: 0, rejected: 1, extra: {n_extra}" diff --git a/src/kimmdy/parsing.py b/src/kimmdy/parsing.py index 929f0aec..6607017f 100644 --- a/src/kimmdy/parsing.py +++ b/src/kimmdy/parsing.py @@ -524,15 +524,15 @@ def read_distances_dat(distances_dat: Path) -> dict: class JSONEncoder(json.JSONEncoder): """Encoder that enables writing JSONs with numpy types.""" - def default(self, obj): - if isinstance(obj, np.integer): - return int(obj) - elif isinstance(obj, np.floating): - return float(obj) - elif isinstance(obj, np.ndarray): - return obj.tolist() + def default(self, o): + if isinstance(o, np.integer): + return int(o) + elif isinstance(o, np.floating): + return float(o) + elif isinstance(o, np.ndarray): + return o.tolist() else: - return super(JSONEncoder, self).default(obj) + return super(JSONEncoder, self).default(o) def write_json( diff --git a/src/kimmdy/plugins.py b/src/kimmdy/plugins.py index dcfe959d..66dea168 100644 --- a/src/kimmdy/plugins.py +++ b/src/kimmdy/plugins.py @@ -106,8 +106,8 @@ def parameterize_topology(self, current_topology: Topology) -> Topology: class BasicParameterizer(Parameterizer): """reconstruct base force field state""" - def parameterize_topology(self, current_topology: Topology) -> None: + def parameterize_topology(self, current_topology: Topology) -> Topology: """Do nothing, all necessary actions should already have happened in bind_bond and break_bond of Topology """ - pass + return current_topology diff --git a/src/kimmdy/recipe.py b/src/kimmdy/recipe.py index 0440eaee..95a2341f 100644 --- a/src/kimmdy/recipe.py +++ b/src/kimmdy/recipe.py @@ -7,7 +7,7 @@ import logging from abc import ABC from copy import copy -from dataclasses import dataclass, field +from dataclasses import dataclass from pathlib import Path from typing import TYPE_CHECKING, Callable, Optional @@ -37,10 +37,11 @@ class Relax(RecipeStep): """ -@dataclass class Place(RecipeStep): """Change topology and/or coordinates to place an atom. + Either provide the index (ix_to_place) or the ID (id_to_place) of the atom to place. + Parameters ---------- new_coords : @@ -54,13 +55,43 @@ class Place(RecipeStep): Index of atom to place. 1-based """ - new_coords: tuple[float, float, float] - ix_to_place: Optional[int] = None - id_to_place: Optional[str] = None + def __init__( + self, + new_coords: tuple[float, float, float], + ix_to_place: Optional[int] = None, + id_to_place: Optional[str] = None, + ): + """Create a new Place RecipeStep instance. - _ix_to_place: Optional[int] = field(init=False, repr=False, default=None) + of the first and second atom, either the id (1-based, str) or ix + (0-based, int) can be given. If both are given, the ix is used. - def __post_init__(self): + Parameters + ---------- + new_coords : + New xyz coordinates for atom to place to. Valid for the end point of the recipe timespan. + ix_to_place : + The index of the atom. zero-based, by default None + id_to_place : + The ID of the atom. one-based, by default None + + Raises + ------ + ValueError + If neither an index nor an ID is provided for any of the atoms. + """ + self.new_coords = new_coords + self._ix_to_place: int + if id_to_place is not None: + if type(id_to_place) is not str: + raise ValueError(f"atom_id_1 is {type(id_to_place)}, should be str.") + self._ix_to_place = int(id_to_place) - 1 + if ix_to_place is not None: + if type(ix_to_place) is not int: + raise ValueError(f"atom_ix_1 is {type(ix_to_place)}, should be int.") + self._ix_to_place = ix_to_place + if id_to_place is not None and ix_to_place is not None: + logger.warning(f"Both atom_ix_1 and atom_id_1 are given, using atom_ix_1.") if self._ix_to_place is None: raise ValueError("id_ or ix_ to_place must be provided!") @@ -77,8 +108,6 @@ def ix_to_place(self, value: Optional[int]): @property def id_to_place(self) -> Optional[str]: - if self._ix_to_place is None: - return None return str(self._ix_to_place + 1) @id_to_place.setter @@ -88,12 +117,33 @@ def id_to_place(self, value: Optional[str]): assert isinstance(value, str), f"id_to_place is {type(value)}, should be str." self._ix_to_place = int(value) - 1 + def __eq__(self, other): + """Two Placements are equal if their atom indices are equal and they have the same new coordinates.""" + + if not isinstance(other, Place) or type(self) != type(other): + logger.warning( + f"Comparing RecipeSteps with different types: {type(self)} and {type(other)}. Returning False." + ) + return False + return ( + self._ix_to_place == other._ix_to_place + and self.new_coords == other.new_coords + ) + + def __hash__(self): + return hash((self._ix_to_place, self.new_coords, type(self).__name__)) + + def __repr__(self): + return f"{type(self).__name__}(ix_to_place={self._ix_to_place}, new_coords={self.new_coords})" + + def __str__(self): + return f"{type(self).__name__}({self._ix_to_place}, {self.new_coords})" + -@dataclass class BondOperation(RecipeStep): """Handle a bond operation on the recipe step. - This class takes in either zero-based indices or one-base IDs for two atoms + This class takes in either zero-based indices or one-base IDs for two atoms. Parameters ---------- @@ -117,27 +167,57 @@ class BondOperation(RecipeStep): """ - atom_ix_1: Optional[int] = None - atom_ix_2: Optional[int] = None - atom_id_1: Optional[str] = None - atom_id_2: Optional[str] = None + def __init__( + self, + atom_ix_1: Optional[int] = None, + atom_ix_2: Optional[int] = None, + atom_id_1: Optional[str] = None, + atom_id_2: Optional[str] = None, + ): + """Create a new BondOperation instance. - _atom_ix_1: int = field(init=False, repr=False, default=None) - _atom_ix_2: int = field(init=False, repr=False, default=None) + of the first and second atom, either the id (1-based, str) or ix + (0-based, int) can be given. If both are given, the ix is used. - def __post_init__(self): - e = "" - if not (isinstance(self.atom_ix_1, int) or isinstance(self.atom_id_1, str)): - e += "Exactly on of atom_ix_1 and atom_id_1 must be given!\n" - if not (isinstance(self.atom_ix_2, int) or isinstance(self.atom_id_2, str)): - e += "Exactly on of atom_ix_2 and atom_id_2 must be given!\n" - if len(e) > 0: - raise ValueError(e) + Parameters + ---------- + atom_ix_1 : int, optional + The index of the first atom. zero-based, by default None + atom_ix_2 : int, optional + The index of the second atom. zero-based, by default None + atom_id_1 : str, optional + The ID of the first atom. one-based, by default None + atom_id_2 : str, optional + The ID of the second atom. one-based, by default None + + Raises + ------ + ValueError + If neither an index nor an ID is provided for any of the atoms. + """ + self._atom_ix_1: int + self._atom_ix_2: int + if atom_id_1 is not None: + if type(atom_id_1) is not str: + raise ValueError(f"atom_id_1 is {type(atom_id_1)}, should be str.") + self._atom_ix_1 = int(atom_id_1) - 1 + if atom_id_2 is not None: + if type(atom_id_2) is not str: + raise ValueError(f"atom_id_2 is {type(atom_id_2)}, should be str.") + self._atom_ix_2 = int(atom_id_2) - 1 + if atom_ix_1 is not None: + if type(atom_ix_1) is not int: + raise ValueError(f"atom_ix_1 is {type(atom_ix_1)}, should be int.") + self._atom_ix_1 = atom_ix_1 + if atom_ix_2 is not None: + if type(atom_ix_2) is not int: + raise ValueError(f"atom_ix_2 is {type(atom_ix_2)}, should be int.") + self._atom_ix_2 = atom_ix_2 + if atom_ix_1 is not None and atom_id_1 is not None: + logger.warning(f"Both atom_ix_1 and atom_id_1 are given, using atom_ix_1.") @property def atom_id_1(self) -> str: - if self._atom_ix_1 is None: - return None return str(self._atom_ix_1 + 1) @atom_id_1.setter @@ -160,8 +240,6 @@ def atom_ix_1(self, value: int): @property def atom_id_2(self) -> str: - if self._atom_ix_2 is None: - return None return str(self._atom_ix_2 + 1) @atom_id_2.setter @@ -182,6 +260,27 @@ def atom_ix_2(self, value: int): assert isinstance(value, int), f"atom_ix_2 is {type(value)}, should be int." self._atom_ix_2 = value + def __eq__(self, other): + """Two BondOperations are equal if their atom indices are equal and they are of the same type.""" + + if not isinstance(other, BondOperation) or type(self) != type(other): + logger.warning( + f"Comparing RecipeSteps with different types: {type(self)} and {type(other)}. Returning False." + ) + return False + return ( + self._atom_ix_1 == other._atom_ix_1 and self._atom_ix_2 == other._atom_ix_2 + ) + + def __hash__(self): + return hash((self._atom_ix_1, self._atom_ix_2, type(self).__name__)) + + def __repr__(self): + return f"{type(self).__name__}(atom_ix_1={self._atom_ix_1}, atom_ix_2={self._atom_ix_2})" + + def __str__(self): + return f"{type(self).__name__}({self._atom_ix_1}, {self._atom_ix_2})" + class Break(BondOperation): """Change topology to break a bond @@ -423,10 +522,10 @@ def from_csv(cls, path: Path): with open(path, "r") as f: reader = csv.reader(f) - header = next(reader) # Skip the header row + _ = next(reader) # Skip the header row picked_rp = None for row in reader: - index_s, picked_s, recipe_steps_s, timespans_s, rates_s = row + _, picked_s, recipe_steps_s, timespans_s, rates_s = row picked = eval(picked_s) recipe_steps = eval(recipe_steps_s) diff --git a/tests/test_recipe.py b/tests/test_recipe.py index 8322107b..23a270f8 100644 --- a/tests/test_recipe.py +++ b/tests/test_recipe.py @@ -1,9 +1,12 @@ -from hypothesis import given, strategies as st -from kimmdy import recipe import csv -import pytest from dataclasses import asdict +import pytest +from hypothesis import given +from hypothesis import strategies as st + +from kimmdy import recipe + ## Test RecipeSteps # tests with random inputs @@ -78,10 +81,10 @@ def test_bo_initialization_unequal(): def test_bo_initialization_wrong_type(): # Should raise an error because initialization is with the wrong type - with pytest.raises(AssertionError): - recipe.BondOperation("1", "2") - with pytest.raises(AssertionError): - recipe.BondOperation(atom_id_1=0, atom_id_2=1) + with pytest.raises(ValueError): + recipe.BondOperation("1", "2") # type: ignore + with pytest.raises(ValueError): + recipe.BondOperation(atom_id_1=0, atom_id_2=1) # type: ignore @given( @@ -126,7 +129,7 @@ def test_place_initialization(): with pytest.raises(TypeError): recipe.Place(ix_to_place=1) - with pytest.raises(AssertionError): + with pytest.raises(ValueError): recipe.Place(id_to_place=1, new_coords=(0, 0, 0))