diff --git a/ml_peg/calcs/bulk_crystal/elasticity/calc_elasticity.py b/ml_peg/calcs/bulk_crystal/elasticity/calc_elasticity.py index 90628dad1..2886f237d 100644 --- a/ml_peg/calcs/bulk_crystal/elasticity/calc_elasticity.py +++ b/ml_peg/calcs/bulk_crystal/elasticity/calc_elasticity.py @@ -141,7 +141,6 @@ def process_result(self, result: dict | None, model_name: str) -> dict: # noqa: dict Dictionary of processed elastic properties. """ - key = get_crystal_system(result["final_structure"]) return { f"bulk_modulus_vrh_{model_name}": ( result["bulk_modulus_vrh"] * eVA3ToGPa @@ -158,7 +157,11 @@ def process_result(self, result: dict | None, model_name: str) -> dict: # noqa: if result is not None else float("nan") ), - f"crystal_system_{model_name}": key, + f"crystal_system_{model_name}": get_crystal_system( + result["final_structure"] + ) + if result is not None + else None, } diff --git a/ml_peg/calcs/bulk_crystal/high_pressure_relaxation/calc_high_pressure_relaxation.py b/ml_peg/calcs/bulk_crystal/high_pressure_relaxation/calc_high_pressure_relaxation.py index 7c274bf89..22e46e4ee 100644 --- a/ml_peg/calcs/bulk_crystal/high_pressure_relaxation/calc_high_pressure_relaxation.py +++ b/ml_peg/calcs/bulk_crystal/high_pressure_relaxation/calc_high_pressure_relaxation.py @@ -9,12 +9,14 @@ import random from typing import Any import urllib.request +from warnings import warn from ase import Atoms from ase.constraints import FixSymmetry from ase.io import write as ase_write from ase.units import GPa from janus_core.calculations.geom_opt import GeomOpt +import numpy as np import pandas as pd from pymatgen.entries.computed_entries import ComputedStructureEntry from pymatgen.io.ase import AseAtomsAdaptor @@ -200,13 +202,17 @@ def relax_with_pressure( converged = True counter += 1 atoms = relaxed - except Exception as e: - print(f"Relaxation failed: {e}") - return None, False, None + except Exception as exc: + warn(f"Relaxation failed: {exc}", stacklevel=2) + return None, False, np.nan if not converged: - return None, False, None + return None, False, np.nan # Calculate enthalpy: H = E + PV - energy = relaxed.get_potential_energy() + try: + energy = relaxed.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + energy = np.nan volume = relaxed.get_volume() enthalpy = energy + pressure_gpa * GPa * volume n_atoms = len(relaxed) diff --git a/ml_peg/calcs/bulk_crystal/iron_properties/calc_iron_properties.py b/ml_peg/calcs/bulk_crystal/iron_properties/calc_iron_properties.py index 86b7f0321..e80e1dc28 100644 --- a/ml_peg/calcs/bulk_crystal/iron_properties/calc_iron_properties.py +++ b/ml_peg/calcs/bulk_crystal/iron_properties/calc_iron_properties.py @@ -16,9 +16,9 @@ from __future__ import annotations import json -import logging from pathlib import Path from typing import TYPE_CHECKING, Any +from warnings import warn from ase.build import bulk from ase.constraints import FixedLine @@ -52,8 +52,6 @@ from ase import Atoms from ase.calculators.calculator import Calculator -logger = logging.getLogger(__name__) - MODELS = load_models(current_models) # Local directory for calculator outputs @@ -139,10 +137,17 @@ def run_eos_calculation(calc: Calculator) -> dict[str, Any]: # Relax atomic positions at fixed cell volume # (matches LAMMPS minimize behavior) - opt = BFGS(atoms, logfile=None) - opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + try: + opt = BFGS(atoms, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception as exc: + warn(f"EoS relaxation failed: {exc}", stacklevel=2) - energy = atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn(f"Energy calculation failed: {exc}", stacklevel=2) + energy = np.nan volume = atoms.get_volume() n_atoms = len(atoms) @@ -196,9 +201,12 @@ def run_elastic_calculation( atoms_ref.calc = calc # Box relaxation - ecf = ExpCellFilter(atoms_ref) - opt = BFGS(ecf, logfile=None) - opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + try: + ecf = ExpCellFilter(atoms_ref) + opt = BFGS(ecf, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception as exc: + warn(f"Elasticity relaxation failed: {exc}", stacklevel=2) # Apply random jiggle to atoms to prevent staying on saddle points rng = np.random.default_rng(seed=87287) @@ -219,8 +227,16 @@ def run_elastic_calculation( atoms_pos.calc = calc opt_pos = BFGS(atoms_pos, logfile=None) - opt_pos.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) - stress_pos = atoms_pos.get_stress(voigt=True) + try: + opt_pos.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception as exc: + warn(f"Positive strain relaxation failed: {exc}", stacklevel=2) + + try: + stress_pos = atoms_pos.get_stress(voigt=True) + except Exception as exc: + warn(f"Positive strain stress calculation failed: {exc}", stacklevel=2) + stress_pos = np.nan * np.ones(6) # Negative strain with off-diagonal cell adjustment atoms_neg = apply_voigt_strain(atoms_ref.copy(), direction, -ELASTIC_STRAIN) @@ -228,8 +244,15 @@ def run_elastic_calculation( atoms_neg.calc = calc opt_neg = BFGS(atoms_neg, logfile=None) - opt_neg.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) - stress_neg = atoms_neg.get_stress(voigt=True) + try: + opt_neg.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception as exc: + warn(f"Negative strain relaxation failed: {exc}", stacklevel=2) + try: + stress_neg = atoms_neg.get_stress(voigt=True) + except Exception as exc: + warn(f"Negative strain stress calculation failed: {exc}", stacklevel=2) + stress_neg = np.nan * np.ones(6) # Compute elastic constants using stress differences # C_ij = dσ_i / dε_j = (σ_pos - σ_neg) / (2 * ε) @@ -308,15 +331,21 @@ def run_bain_path_calculation( atoms_relaxed = relax_volume_isotropic(atoms, calc) # Step 2: Atomic position relaxation at fixed cell - opt = BFGS(atoms_relaxed, logfile=None) try: + opt = BFGS(atoms_relaxed, logfile=None) opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) - except Exception: - logger.warning( - "BFGS relaxation failed for Bain path c/a=%.3f", ratio, exc_info=True + except Exception as exc: + warn(f"Relaxation failed for Bain path c/a={ratio}: {exc}", stacklevel=2) + + try: + energy = atoms_relaxed.get_potential_energy() + except Exception as exc: + warn( + f"Energy calculation failed for Bain path c/a={ratio}: {exc}", + stacklevel=2, ) + energy = np.nan - energy = atoms_relaxed.get_potential_energy() cell = atoms_relaxed.get_cell() ca_actual = cell[2, 2] / cell[1, 1] @@ -376,7 +405,12 @@ def run_vacancy_calculation( atoms_perfect.calc = calc n_atoms = len(atoms_perfect) - E_perfect = atoms_perfect.get_potential_energy() # noqa: N806 + try: + E_perfect = atoms_perfect.get_potential_energy() # noqa: N806 + except Exception as exc: + warn(f"Energy calculation failed for perfect BCC: {exc}", stacklevel=2) + E_perfect = np.nan # noqa: N806 + E_coh = E_perfect / n_atoms # noqa: N806 atoms_defect = atoms_perfect.copy() @@ -384,9 +418,18 @@ def run_vacancy_calculation( _set_iron_info(atoms_defect) atoms_defect.calc = calc - opt = BFGS(atoms_defect, logfile=None) - opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) - E_defect = atoms_defect.get_potential_energy() # noqa: N806 + try: + opt = BFGS(atoms_defect, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception as exc: + warn(f"Relaxation failed for vacancy defect: {exc}", stacklevel=2) + + try: + E_defect = atoms_defect.get_potential_energy() # noqa: N806 + except Exception as exc: + warn(f"Energy calculation failed for vacancy defect: {exc}", stacklevel=2) + E_defect = np.nan # noqa: N806 + E_vac = (E_defect - E_perfect) + E_coh # noqa: N806 return { @@ -452,6 +495,7 @@ def run_surface_calculations( for name, cfg in SURFACE_CONFIG.items(): print(f"Running surface calculation for {name}...") + create_fn = cfg["create_fn"] area_axes = cfg["area_axes"] vacuum = cfg["vacuum"] @@ -465,10 +509,23 @@ def run_surface_calculations( slab_kwargs = {"layers": cfg["layers"], "vacuum": vacuum} # Bulk reference - atoms_bulk = create_fn(lattice_parameter, **bulk_kwargs) + try: + atoms_bulk = create_fn(lattice_parameter, **bulk_kwargs) + except ValueError as err: + warn(f"Failed to create bulk structure for {name}: {err}", stacklevel=2) + surfaces[name] = np.nan + continue + _set_iron_info(atoms_bulk) atoms_bulk.calc = calc - e_bulk = atoms_bulk.get_potential_energy() + try: + e_bulk = atoms_bulk.get_potential_energy() + except Exception as exc: + warn( + f"Energy calculation failed for bulk reference {name}: {exc}", + stacklevel=2, + ) + e_bulk = np.nan cell = atoms_bulk.get_cell() area = np.linalg.norm(np.cross(cell[area_axes[0]], cell[area_axes[1]])) @@ -476,9 +533,17 @@ def run_surface_calculations( atoms_slab = create_fn(lattice_parameter, **slab_kwargs) _set_iron_info(atoms_slab) atoms_slab.calc = calc - opt = BFGS(atoms_slab, logfile=None) - opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) - e_slab = atoms_slab.get_potential_energy() + try: + opt = BFGS(atoms_slab, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception as exc: + warn(f"Relaxation failed for {name}: {exc}", stacklevel=2) + + try: + e_slab = atoms_slab.get_potential_energy() + except Exception as exc: + warn(f"Energy calculation failed for {name}: {exc}", stacklevel=2) + e_slab = np.nan surfaces[name] = calculate_surface_energy(e_slab, e_bulk, area) @@ -516,23 +581,41 @@ def run_sfe_calculation( dict[str, Any] Dictionary with displacements, sfe_J_per_m2, and max_sfe. """ - config = SFE_CONFIG[sfe_type] - atoms = config["create_fn"](lattice_parameter) - _set_iron_info(atoms) - atoms.calc = calc - # Calculate Burgers vector magnitude: b = a * sqrt(3) / 2 burgers_vector = lattice_parameter * np.sqrt(3) / 2 step_size = burgers_vector / SFE_STEPS + config = SFE_CONFIG[sfe_type] + + try: + atoms = config["create_fn"](lattice_parameter) + except Exception as err: + warn(f"Failed to create SFE structure for {sfe_type}: {err}", stacklevel=2) + return { + "displacements": [step * step_size for step in range(SFE_STEPS + 1)], + "sfe_J_per_m2": [np.nan for _ in range(SFE_STEPS + 1)], + "max_sfe": np.nan, + } + + _set_iron_info(atoms) + atoms.calc = calc + cell = atoms.get_cell() ly = cell[1, 1] lz = cell[2, 2] area = ly * lz - opt = BFGS(atoms, logfile=None) - opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) - e0 = atoms.get_potential_energy() + try: + opt = BFGS(atoms, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception as exc: + warn(f"Relaxation failed for SFE {sfe_type}: {exc}", stacklevel=2) + + try: + e0 = atoms.get_potential_energy() + except Exception as exc: + warn(f"Energy calculation failed for {sfe_type}: {exc}", stacklevel=2) + e0 = np.nan positions = atoms.get_positions() x_mid = (positions[:, 0].min() + positions[:, 0].max()) / 2 + 0.1 @@ -552,19 +635,23 @@ def run_sfe_calculation( atoms.set_constraint(constraints) - opt = BFGS(atoms, logfile=None) try: + opt = BFGS(atoms, logfile=None) opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) - except Exception: - logger.warning( - "BFGS relaxation failed for SFE %s step %d", - sfe_type, - step, - exc_info=True, - ) + except Exception as exc: + warn(f"Relaxation failed for SFE {sfe_type}: {exc}", stacklevel=2) atoms.set_constraint() - energy = atoms.get_potential_energy() + + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn( + f"Energy calculation failed for SFE {sfe_type} step {step}: {exc}", + stacklevel=2, + ) + energy = np.nan + sfe = (energy - e0) / (2 * area) * EV_PER_A2_TO_J_PER_M2 displacements.append(step * step_size) @@ -650,13 +737,22 @@ def run_ts_calculation( atoms.set_positions(positions) # Calculate energy (no relaxation!) - energy = atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn(f"Energy calculation failed for TS step {dd}: {exc}", stacklevel=2) + energy = np.nan # Calculate forces for stress - forces = atoms.get_forces() + try: + forces = atoms.get_forces() - # Sum of z-forces on upper region - fz_upper = np.sum(forces[upper_indices, 2]) + # Sum of z-forces on upper region + fz_upper = np.sum(forces[upper_indices, 2]) + except Exception as exc: + warn(f"Force calculation failed for TS step {dd}: {exc}", stacklevel=2) + forces = np.nan + fz_upper = np.nan # Convert to stress (GPa): σ = F / A # Negate because forces on upper atoms point downward (negative z) @@ -836,14 +932,18 @@ def run_iron_properties(model_name: str, model: Any) -> None: write(structures_dir / "vacancy.extxyz", vac_atoms) for name, cfg in SURFACE_CONFIG.items(): - create_fn = cfg["create_fn"] - if "size" in cfg: - slab = create_fn(a0, size=cfg["size"], vacuum=cfg["vacuum"]) - else: - slab = create_fn(a0, layers=cfg["layers"], vacuum=cfg["vacuum"]) - _set_iron_info(slab) - slab.info["description"] = f"BCC Fe ({name}) surface slab" - write(structures_dir / f"surface_{name}.extxyz", slab) + try: + create_fn = cfg["create_fn"] + if "size" in cfg: + slab = create_fn(a0, size=cfg["size"], vacuum=cfg["vacuum"]) + else: + slab = create_fn(a0, layers=cfg["layers"], vacuum=cfg["vacuum"]) + _set_iron_info(slab) + slab.info["description"] = f"BCC Fe ({name}) surface slab" + write(structures_dir / f"surface_{name}.extxyz", slab) + except Exception as exc: + warn(f"Failed to create structure for surface {name}: {exc}", stacklevel=2) + continue # Save all results as JSON (write_dir / "results.json").write_text(json.dumps(results, indent=2, default=str)) diff --git a/ml_peg/calcs/bulk_crystal/iron_properties/iron_utils.py b/ml_peg/calcs/bulk_crystal/iron_properties/iron_utils.py index 6b015d70e..c2753dc6b 100644 --- a/ml_peg/calcs/bulk_crystal/iron_properties/iron_utils.py +++ b/ml_peg/calcs/bulk_crystal/iron_properties/iron_utils.py @@ -7,6 +7,7 @@ from __future__ import annotations from typing import TYPE_CHECKING +from warnings import warn from ase import Atoms from ase.build import bulk @@ -75,12 +76,15 @@ def fit_eos( - V0: Equilibrium volume per atom (Angstrom^3) - a0: Equilibrium lattice parameter (Angstrom) - for BCC """ - eos = EquationOfState(vol, ene, eos="birchmurnaghan") - V0, E0, B0_raw = eos.fit() # noqa: N806 - - B0_GPa = B0_raw * EV_PER_A3_TO_GPA # noqa: N806 - Bp = eos.eos_parameters[2] # noqa: N806 - a0 = (V0 * 2) ** (1.0 / 3.0) + try: + eos = EquationOfState(vol, ene, eos="birchmurnaghan") + V0, E0, B0_raw = eos.fit() # noqa: N806 + B0_GPa = B0_raw * EV_PER_A3_TO_GPA # noqa: N806 + Bp = eos.eos_parameters[2] # noqa: N806 + a0 = (V0 * 2) ** (1.0 / 3.0) + except Exception as exc: + warn(f"EOS fitting failed: {exc}", stacklevel=2) + return dict.fromkeys(["E0", "B0", "Bp", "V0", "a0"], np.nan) return {"E0": E0, "B0": B0_GPa, "Bp": Bp, "V0": V0, "a0": a0} @@ -150,7 +154,14 @@ def energy_at_scale(scale: float) -> float: test_atoms = atoms.copy() test_atoms.set_cell(original_cell * scale, scale_atoms=True) test_atoms.calc = calc - return test_atoms.get_potential_energy() + try: + return test_atoms.get_potential_energy() + except Exception as exc: + warn( + f"Energy calculation failed for scale factor {scale}: {exc}", + stacklevel=2, + ) + return np.nan # Find optimal scale factor that minimizes energy result = minimize_scalar( diff --git a/ml_peg/calcs/bulk_crystal/lattice_constants/calc_lattice_constants.py b/ml_peg/calcs/bulk_crystal/lattice_constants/calc_lattice_constants.py index deb3b8c88..2c0a87b93 100644 --- a/ml_peg/calcs/bulk_crystal/lattice_constants/calc_lattice_constants.py +++ b/ml_peg/calcs/bulk_crystal/lattice_constants/calc_lattice_constants.py @@ -6,6 +6,7 @@ import json from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms from ase.build import bulk @@ -136,9 +137,12 @@ def test_lattice_consts(mlip: tuple[str, Any]) -> None: for name, crystal in crystals.items(): crystal.calc = copy(calc) - GeomOpt( - struct=crystal, - fmax=0.03, - write_traj=True, - file_prefix=OUT_PATH / model_name / name, - ).run() + try: + GeomOpt( + struct=crystal, + fmax=0.03, + write_traj=True, + file_prefix=OUT_PATH / model_name / name, + ).run() + except Exception as exc: + warn(f"Error during optimisation for {name}: {exc}", stacklevel=2) diff --git a/ml_peg/calcs/bulk_crystal/low_dimensional_relaxation/calc_low_dimensional_relaxation.py b/ml_peg/calcs/bulk_crystal/low_dimensional_relaxation/calc_low_dimensional_relaxation.py index 6f10dc914..9fb30d558 100644 --- a/ml_peg/calcs/bulk_crystal/low_dimensional_relaxation/calc_low_dimensional_relaxation.py +++ b/ml_peg/calcs/bulk_crystal/low_dimensional_relaxation/calc_low_dimensional_relaxation.py @@ -9,6 +9,7 @@ import random from typing import Any import urllib.request +from warnings import warn from ase import Atoms from ase.io import write as ase_write @@ -335,14 +336,21 @@ def relax_low_dimensional( converged = True counter += 1 atoms = relaxed - except Exception as e: - print(f"Relaxation failed: {e}") - return None, False, None, None + except Exception as exc: + warn(f"Relaxation failed: {exc}", stacklevel=2) + return None, False, np.nan, None if relaxed is None: - return None, False, None, max_force + return None, False, np.nan, max_force - energy_per_atom = relaxed.get_potential_energy() / len(relaxed) + try: + energy_per_atom = relaxed.get_potential_energy() / len(relaxed) + except Exception as exc: + warn( + f"Error calculating energy: {exc}", + stacklevel=2, + ) + energy_per_atom = np.nan return relaxed, converged, energy_per_atom, max_force diff --git a/ml_peg/calcs/conformers/37Conf8/calc_37Conf8.py b/ml_peg/calcs/conformers/37Conf8/calc_37Conf8.py index d9a4dce11..ac7d2e7e6 100644 --- a/ml_peg/calcs/conformers/37Conf8/calc_37Conf8.py +++ b/ml_peg/calcs/conformers/37Conf8/calc_37Conf8.py @@ -8,9 +8,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pandas as pd import pytest from tqdm import tqdm @@ -65,14 +67,26 @@ def test_37conf8_conformer_energies(mlip: tuple[str, Any]) -> None: zero_conf.info["charge"] = 0 zero_conf.info["spin"] = 1 zero_conf.calc = calc - e_model_zero_conf = zero_conf.get_potential_energy() + try: + e_model_zero_conf = zero_conf.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + e_model_zero_conf = np.nan else: atoms = read(data_path / "PBEPBE-D3" / f"{label}_PBEPBE-D3.xyz") atoms.info["charge"] = 0 atoms.info["spin"] = 1 atoms.calc = calc - atoms.info["model_rel_energy"] = ( - atoms.get_potential_energy() - e_model_zero_conf - ) + try: + atoms.info["model_rel_energy"] = ( + atoms.get_potential_energy() - e_model_zero_conf + ) + except Exception as exc: + warn( + f"Error calculating relative energy for {label}: {exc}", + stacklevel=2, + ) + atoms.info["model_rel_energy"] = np.nan + atoms.info["ref_energy"] = float(df.iloc[i][2]) * KCAL_TO_EV write(write_dir / f"{label}.xyz", atoms) diff --git a/ml_peg/calcs/conformers/ACONFL/calc_ACONFL.py b/ml_peg/calcs/conformers/ACONFL/calc_ACONFL.py index 9324a0485..7233f8e3c 100644 --- a/ml_peg/calcs/conformers/ACONFL/calc_ACONFL.py +++ b/ml_peg/calcs/conformers/ACONFL/calc_ACONFL.py @@ -11,9 +11,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -65,9 +67,17 @@ def test_aconfl_conformer_energies(mlip: tuple[str, Any]) -> None: zero_atoms = read(data_path / zero_atoms_label / "struc.xyz") zero_atoms.calc = calc zero_atoms.info.update({"charge": 0, "spin": 1}) - atoms.info["model_rel_energy"] = ( - atoms.get_potential_energy() - zero_atoms.get_potential_energy() - ) + try: + atoms.info["model_rel_energy"] = ( + atoms.get_potential_energy() - zero_atoms.get_potential_energy() + ) + except Exception as exc: + warn( + f"Error calculating energy for {atoms_label} or " + f"{zero_atoms_label}: {exc}", + stacklevel=2, + ) + atoms.info["model_rel_energy"] = np.nan atoms.info["ref_rel_energy"] = ref_rel_energy write_dir = OUT_PATH / model_name diff --git a/ml_peg/calcs/conformers/DipCONFS/calc_DipCONFS.py b/ml_peg/calcs/conformers/DipCONFS/calc_DipCONFS.py index 313733ac6..a90186d84 100644 --- a/ml_peg/calcs/conformers/DipCONFS/calc_DipCONFS.py +++ b/ml_peg/calcs/conformers/DipCONFS/calc_DipCONFS.py @@ -13,9 +13,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pandas as pd import pytest from tqdm import tqdm @@ -97,14 +99,22 @@ def test_dipconfs(mlip: tuple[str, Any]) -> None: # Get zero ref conformer model energy zero_conf = get_atoms(data_path / zero_conf_label / "struc.xyz") zero_conf.calc = calc - e_model_zero_conf = zero_conf.get_potential_energy() + try: + e_model_zero_conf = zero_conf.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {zero_conf_label}: {exc}", stacklevel=2) + e_model_zero_conf = np.nan # Get current conformer model energy atoms = get_atoms(data_path / label / "struc.xyz") atoms.calc = calc - atoms.info["model_rel_energy"] = ( - atoms.get_potential_energy() - e_model_zero_conf - ) + try: + atoms.info["model_rel_energy"] = ( + atoms.get_potential_energy() - e_model_zero_conf + ) + except Exception as exc: + warn(f"Error calculating relative energy for {label}: {exc}", stacklevel=2) + atoms.info["model_rel_energy"] = np.nan atoms.info["ref_energy"] = e_rel_ref write_dir = OUT_PATH / model_name diff --git a/ml_peg/calcs/conformers/Glucose205/calc_Glucose205.py b/ml_peg/calcs/conformers/Glucose205/calc_Glucose205.py index 2c08c03d9..7e0648b47 100644 --- a/ml_peg/calcs/conformers/Glucose205/calc_Glucose205.py +++ b/ml_peg/calcs/conformers/Glucose205/calc_Glucose205.py @@ -10,9 +10,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write +import numpy as np import pandas as pd import pytest from tqdm import tqdm @@ -122,7 +124,11 @@ def test_glucose205(mlip: tuple[str, Any]) -> None: data_path / "Glucose_structures" / f"{lowest_conf_label}.xyz" ) conf_lowest.calc = calc - e_conf_lowest_model = conf_lowest.get_potential_energy() + try: + e_conf_lowest_model = conf_lowest.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {lowest_conf_label}: {exc}", stacklevel=2) + e_conf_lowest_model = np.nan for label, e_ref in tqdm(ref_energies.items()): # Skip the reference conformer for which the error is automatically zero @@ -131,9 +137,13 @@ def test_glucose205(mlip: tuple[str, Any]) -> None: atoms = get_atoms(data_path / "Glucose_structures" / f"{label}.xyz") atoms.calc = calc - atoms.info["model_rel_energy"] = ( - atoms.get_potential_energy() - e_conf_lowest_model - ) + try: + atoms.info["model_rel_energy"] = ( + atoms.get_potential_energy() - e_conf_lowest_model + ) + except Exception as exc: + warn(f"Error calculating relative energy for {label}: {exc}", stacklevel=2) + atoms.info["model_rel_energy"] = np.nan atoms.info["ref_energy"] = e_ref write_dir = OUT_PATH / model_name diff --git a/ml_peg/calcs/conformers/MPCONF196/calc_MPCONF196.py b/ml_peg/calcs/conformers/MPCONF196/calc_MPCONF196.py index 145bd924a..a033fabf2 100644 --- a/ml_peg/calcs/conformers/MPCONF196/calc_MPCONF196.py +++ b/ml_peg/calcs/conformers/MPCONF196/calc_MPCONF196.py @@ -9,6 +9,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write @@ -135,7 +136,11 @@ def test_mpconf196(mlip: tuple[str, Any]) -> None: atoms = get_atoms(data_path / xyz_fname) atoms.translate(-atoms.get_center_of_mass()) atoms.calc = calc - model_abs_energies.append(atoms.get_potential_energy()) + try: + model_abs_energies.append(atoms.get_potential_energy()) + except Exception as exc: + warn(f"Error calculating energy for {xyz_fname}: {exc}", stacklevel=2) + model_abs_energies.append(np.nan) ref_abs_energies.append(e_ref) current_molecule_labels.append(label) diff --git a/ml_peg/calcs/conformers/Maltose222/calc_Maltose222.py b/ml_peg/calcs/conformers/Maltose222/calc_Maltose222.py index 0b983b612..f3c503f64 100644 --- a/ml_peg/calcs/conformers/Maltose222/calc_Maltose222.py +++ b/ml_peg/calcs/conformers/Maltose222/calc_Maltose222.py @@ -10,9 +10,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pandas as pd import pytest from tqdm import tqdm @@ -122,7 +124,11 @@ def test_maltose222(mlip: tuple[str, Any]) -> None: data_path / "Maltose_structures" / f"{lowest_conf_label}.xyz" ) conf_lowest.calc = calc - e_conf_lowest_model = conf_lowest.get_potential_energy() + try: + e_conf_lowest_model = conf_lowest.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {lowest_conf_label}: {exc}", stacklevel=2) + e_conf_lowest_model = np.nan for label, e_ref in tqdm(ref_energies.items()): # Skip the reference conformer for which the error is automatically zero @@ -131,9 +137,13 @@ def test_maltose222(mlip: tuple[str, Any]) -> None: atoms = get_atoms(data_path / "Maltose_structures" / f"{label}.xyz") atoms.calc = calc - atoms.info["model_rel_energy"] = ( - atoms.get_potential_energy() - e_conf_lowest_model - ) + try: + atoms.info["model_rel_energy"] = ( + atoms.get_potential_energy() - e_conf_lowest_model + ) + except Exception as exc: + warn(f"Error calculating relative energy for {label}: {exc}", stacklevel=2) + atoms.info["model_rel_energy"] = np.nan atoms.info["ref_energy"] = e_ref write_dir = OUT_PATH / model_name diff --git a/ml_peg/calcs/conformers/OpenFF_Tors/calc_OpenFF_Tors.py b/ml_peg/calcs/conformers/OpenFF_Tors/calc_OpenFF_Tors.py index f8f164623..fb4a42180 100644 --- a/ml_peg/calcs/conformers/OpenFF_Tors/calc_OpenFF_Tors.py +++ b/ml_peg/calcs/conformers/OpenFF_Tors/calc_OpenFF_Tors.py @@ -10,6 +10,7 @@ import json from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import write @@ -81,14 +82,25 @@ def test_openff_tors(mlip: tuple[str, Any]) -> None: if i == 0: e_ref_zero_conf = ref_energy * units.Hartree - e_model_zero_conf = atoms.get_potential_energy() + try: + e_model_zero_conf = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + e_model_zero_conf = np.nan else: atoms.info["ref_rel_energy"] = ( ref_energy * units.Hartree - e_ref_zero_conf ) - atoms.info["model_rel_energy"] = ( - atoms.get_potential_energy() - e_model_zero_conf - ) + try: + atoms.info["model_rel_energy"] = ( + atoms.get_potential_energy() - e_model_zero_conf + ) + except Exception as exc: + warn( + f"Error calculating relative energy for {label}: {exc}", + stacklevel=2, + ) + atoms.info["model_rel_energy"] = np.nan write_dir = OUT_PATH / model_name write_dir.mkdir(parents=True, exist_ok=True) write(write_dir / f"{label}.xyz", atoms) diff --git a/ml_peg/calcs/conformers/UpU46/calc_UpU46.py b/ml_peg/calcs/conformers/UpU46/calc_UpU46.py index f7d8556cf..0096ec6fe 100644 --- a/ml_peg/calcs/conformers/UpU46/calc_UpU46.py +++ b/ml_peg/calcs/conformers/UpU46/calc_UpU46.py @@ -10,9 +10,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -104,7 +106,11 @@ def test_upu46(mlip: tuple[str, Any]) -> None: conf_lowest = get_atoms(data_path / f"{zero_conf_label}.xyz") conf_lowest.calc = calc - e_conf_lowest_model = conf_lowest.get_potential_energy() + try: + e_conf_lowest_model = conf_lowest.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {zero_conf_label}: {exc}", stacklevel=2) + e_conf_lowest_model = np.nan for label, e_ref in tqdm(ref_energies.items()): # Skip the reference conformer for @@ -114,9 +120,13 @@ def test_upu46(mlip: tuple[str, Any]) -> None: atoms = get_atoms(data_path / f"{label}.xyz") atoms.calc = calc - atoms.info["model_rel_energy"] = ( - atoms.get_potential_energy() - e_conf_lowest_model - ) + try: + atoms.info["model_rel_energy"] = ( + atoms.get_potential_energy() - e_conf_lowest_model + ) + except Exception as exc: + warn(f"Error calculating relative energy for {label}: {exc}", stacklevel=2) + atoms.info["model_rel_energy"] = np.nan atoms.info["ref_energy"] = e_ref atoms.calc = None diff --git a/ml_peg/calcs/conformers/solvMPCONF196/calc_solvMPCONF196.py b/ml_peg/calcs/conformers/solvMPCONF196/calc_solvMPCONF196.py index c6cc078dd..be51974af 100644 --- a/ml_peg/calcs/conformers/solvMPCONF196/calc_solvMPCONF196.py +++ b/ml_peg/calcs/conformers/solvMPCONF196/calc_solvMPCONF196.py @@ -9,6 +9,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write @@ -139,7 +140,11 @@ def test_solvmpconf196(mlip: tuple[str, Any]) -> None: ) atoms.translate(-atoms.get_center_of_mass()) atoms.calc = calc - model_abs_energies.append(atoms.get_potential_energy()) + try: + model_abs_energies.append(atoms.get_potential_energy()) + except Exception as exc: + warn(f"Error calculating energy for {xyz_label}: {exc}", stacklevel=2) + model_abs_energies.append(np.nan) ref_abs_energies.append(e_ref) current_molecule_labels.append(label) diff --git a/ml_peg/calcs/defect/Defectstab/calc_Defectstab.py b/ml_peg/calcs/defect/Defectstab/calc_Defectstab.py index 7f58a9cee..21477f736 100644 --- a/ml_peg/calcs/defect/Defectstab/calc_Defectstab.py +++ b/ml_peg/calcs/defect/Defectstab/calc_Defectstab.py @@ -4,8 +4,10 @@ from pathlib import Path from typing import Any +from warnings import warn from ase.io import read, write +import numpy as np import pytest from ml_peg.calcs.utils.utils import download_s3_data @@ -105,7 +107,14 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: # Assuming existing logic for fe_sia where ref.poscar header has E_bulk # And we need to calculate E_bulk_pred atoms_ref.calc = calc - e_ref_bulk_pred = atoms_ref.get_potential_energy() + try: + e_ref_bulk_pred = atoms_ref.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating bulk reference energy for {ref_file}: {exc}", + stacklevel=2, + ) + e_ref_bulk_pred = np.nan n_bulk = len(atoms_ref) # Save ref calculation @@ -126,7 +135,11 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: e_dft = get_ref_energy(poscar_path) atoms.calc = calc - e_pred = atoms.get_potential_energy() + try: + e_pred = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {poscar_path}: {exc}", stacklevel=2) + e_pred = np.nan atoms.info["energy_dft"] = e_dft atoms.info["energy_pred"] = e_pred @@ -175,7 +188,11 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: atoms.info.setdefault("spin", 1) e_dft = get_ref_energy(p_path) atoms.calc = calc - e_pred = atoms.get_potential_energy() # Just to verify/store if needed + try: + e_pred = atoms.get_potential_energy() # Just to verify/store if needed + except Exception as exc: + warn(f"Error calculating energy for {p_path}: {exc}", stacklevel=2) + e_pred = np.nan n_atoms = len(atoms) ref_energies[f"E_{el}_dft"] = e_dft / n_atoms @@ -200,7 +217,11 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: atoms_nd.info.setdefault("charge", 0) atoms_nd.info.setdefault("spin", 1) atoms_nd.calc = calc - e_no_defect_pred = atoms_nd.get_potential_energy() + try: + e_no_defect_pred = atoms_nd.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {f}: {exc}", stacklevel=2) + e_no_defect_pred = np.nan break # Process Boron Carbide Subsets @@ -217,7 +238,11 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: e_dft = get_ref_energy(poscar_path) atoms.calc = calc - e_pred = atoms.get_potential_energy() + try: + e_pred = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {poscar_path}: {exc}", stacklevel=2) + e_pred = np.nan atoms.info["energy_dft"] = e_dft atoms.info["energy_pred"] = e_pred @@ -276,7 +301,11 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: atoms_mai.info.setdefault("charge", 0) atoms_mai.info.setdefault("spin", 1) atoms_mai.calc = calc - e_mai_pred = atoms_mai.get_potential_energy() + try: + e_mai_pred = atoms_mai.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {p}: {exc}", stacklevel=2) + e_mai_pred = np.nan elif "pristine" in p.name: e_pris_dft = get_ref_energy(p) atoms_pris = read(p) @@ -284,7 +313,11 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: atoms_pris.info.setdefault("charge", 0) atoms_pris.info.setdefault("spin", 1) atoms_pris.calc = calc - e_pris_pred = atoms_pris.get_potential_energy() + try: + e_pris_pred = atoms_pris.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {p}: {exc}", stacklevel=2) + e_pris_pred = np.nan for poscar_path in mapi_files: atoms = read(poscar_path) @@ -294,7 +327,11 @@ def test_defectstab(mlip: tuple[str, Any]) -> None: e_dft = get_ref_energy(poscar_path) atoms.calc = calc - e_pred = atoms.get_potential_energy() + try: + e_pred = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {poscar_path}: {exc}", stacklevel=2) + e_pred = np.nan atoms.info["energy_dft"] = e_dft atoms.info["energy_pred"] = e_pred diff --git a/ml_peg/calcs/defect/Relastab/calc_Relastab.py b/ml_peg/calcs/defect/Relastab/calc_Relastab.py index 1d6df0a48..ba7fbcfaf 100644 --- a/ml_peg/calcs/defect/Relastab/calc_Relastab.py +++ b/ml_peg/calcs/defect/Relastab/calc_Relastab.py @@ -4,8 +4,10 @@ from pathlib import Path from typing import Any +from warnings import warn from ase.io import read, write +import numpy as np import pytest from ml_peg.calcs.utils.utils import download_s3_data @@ -71,8 +73,12 @@ def test_relastab(mlip: tuple[str, Any]) -> None: # Calculate atoms.calc = calc - - atoms.get_potential_energy() + try: + e_pred = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {poscar}: {exc}", stacklevel=2) + e_pred = np.nan + atoms.info["energy_pred"] = e_pred atoms.info["system"] = poscar.stem atoms.info["subset"] = subset_name diff --git a/ml_peg/calcs/lanthanides/isomer_complexes/calc_isomer_complexes.py b/ml_peg/calcs/lanthanides/isomer_complexes/calc_isomer_complexes.py index f662009e0..afc08b958 100644 --- a/ml_peg/calcs/lanthanides/isomer_complexes/calc_isomer_complexes.py +++ b/ml_peg/calcs/lanthanides/isomer_complexes/calc_isomer_complexes.py @@ -5,9 +5,11 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -125,7 +127,11 @@ def test_isomer_complexes(mlip: tuple[str, Any]) -> None: atoms.info["spin_multiplicity"] = entry["multiplicity"] atoms.info["spin"] = entry["multiplicity"] atoms.calc = copy(calc) - atoms.info["model_energy"] = atoms.get_potential_energy() + try: + atoms.info["model_energy"] = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {entry['xyz']}: {exc}", stacklevel=2) + atoms.info["model_energy"] = np.nan atoms.info["model"] = model_name atoms.info["system"] = entry["system"] diff --git a/ml_peg/calcs/molecular/BMIMCl_RDF/calc_BMIMCl_RDF.py b/ml_peg/calcs/molecular/BMIMCl_RDF/calc_BMIMCl_RDF.py index 850f401d9..b05a95f25 100644 --- a/ml_peg/calcs/molecular/BMIMCl_RDF/calc_BMIMCl_RDF.py +++ b/ml_peg/calcs/molecular/BMIMCl_RDF/calc_BMIMCl_RDF.py @@ -10,6 +10,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import write @@ -17,6 +18,7 @@ from ase.md.velocitydistribution import MaxwellBoltzmannDistribution from ase.optimize import LBFGS import molify +import numpy as np import pytest from tqdm import tqdm @@ -59,7 +61,10 @@ def test_bmimcl_md(mlip: tuple[str, Any]) -> None: box.calc = calc opt = LBFGS(box) - opt.run(fmax=0.1, steps=1000) + try: + opt.run(fmax=0.1, steps=1000) + except Exception as exc: + warn(f"Error during optimisation: {exc}", stacklevel=2) MaxwellBoltzmannDistribution(box, temperature_K=TEMPERATURE) @@ -78,5 +83,10 @@ def test_bmimcl_md(mlip: tuple[str, Any]) -> None: traj_file.unlink() for _ in tqdm(range(STEPS), desc=f"{model_name} MD"): - dyn.run(1) + try: + dyn.run(1) + except Exception as exc: + warn(f"Error during MD step for {model_name}: {exc}", stacklevel=2) + box.info["energy"] = np.nan + break write(traj_file, box, append=True) diff --git a/ml_peg/calcs/molecular/GMTKN55/calc_GMTKN55.py b/ml_peg/calcs/molecular/GMTKN55/calc_GMTKN55.py index 42de08e75..e75a093ad 100644 --- a/ml_peg/calcs/molecular/GMTKN55/calc_GMTKN55.py +++ b/ml_peg/calcs/molecular/GMTKN55/calc_GMTKN55.py @@ -5,6 +5,7 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms from ase.io import write @@ -97,7 +98,17 @@ def test_gmtkn55(mlip: tuple[str, Any]) -> None: atoms.pbc = False atoms.calc = copy(calc) - atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {species_name}: {exc}", + stacklevel=2, + ) + energy = np.nan + + atoms.info["energy"] = energy + atoms.calc = None system_structs.append(atoms) diff --git a/ml_peg/calcs/molecular/Wiggle150/calc_Wiggle150.py b/ml_peg/calcs/molecular/Wiggle150/calc_Wiggle150.py index 764005ab1..b9297a6ac 100644 --- a/ml_peg/calcs/molecular/Wiggle150/calc_Wiggle150.py +++ b/ml_peg/calcs/molecular/Wiggle150/calc_Wiggle150.py @@ -5,10 +5,12 @@ from collections.abc import Iterable from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.calculators.calculator import Calculator from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -111,7 +113,11 @@ def get_energy(atoms: Atoms, calc: Calculator) -> float: """ atoms_copy = atoms.copy() atoms_copy.calc = calc - energy = atoms_copy.get_potential_energy() + try: + energy = atoms_copy.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + return np.nan return float(energy) diff --git a/ml_peg/calcs/molecular_crystal/CPOSS209/calc_CPOSS209.py b/ml_peg/calcs/molecular_crystal/CPOSS209/calc_CPOSS209.py index e1a3a5d06..df0df28aa 100644 --- a/ml_peg/calcs/molecular_crystal/CPOSS209/calc_CPOSS209.py +++ b/ml_peg/calcs/molecular_crystal/CPOSS209/calc_CPOSS209.py @@ -5,6 +5,7 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase.io import read, write import numpy as np @@ -74,7 +75,13 @@ def test_lattice_energy(mlip: tuple[str, Any]) -> None: solid.info.setdefault("charge", 0) solid.info.setdefault("spin", 1) solid.calc = copy(calc) - solid.get_potential_energy() + try: + solid.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {crystal_file}: {exc}", stacklevel=2 + ) + solid.info["energy"] = np.nan solid.info["ref"] = ref_crystal solid.info["num_molecules"] = num_mol @@ -132,7 +139,13 @@ def test_lattice_energy(mlip: tuple[str, Any]) -> None: molecule.info.setdefault("charge", 0) molecule.info.setdefault("spin", 1) molecule.calc = copy(calc) - molecule.get_potential_energy() + try: + molecule.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {molecule_file}: {exc}", stacklevel=2 + ) + molecule.info["energy"] = np.nan # One molecule in each gas phase file molecule.info["num_molecules"] = 1 diff --git a/ml_peg/calcs/molecular_crystal/DMC_ICE13/calc_DMC_ICE13.py b/ml_peg/calcs/molecular_crystal/DMC_ICE13/calc_DMC_ICE13.py index 4d6ae4552..0c83e9b91 100644 --- a/ml_peg/calcs/molecular_crystal/DMC_ICE13/calc_DMC_ICE13.py +++ b/ml_peg/calcs/molecular_crystal/DMC_ICE13/calc_DMC_ICE13.py @@ -6,8 +6,10 @@ import json from pathlib import Path from typing import Any +from warnings import warn from ase.io import read, write +import numpy as np import pytest from ml_peg.calcs.utils.utils import download_s3_data @@ -58,7 +60,11 @@ def test_lattice_energy(mlip: tuple[str, Any]) -> None: # Set default charge and spin water.info.setdefault("charge", 0) water.info.setdefault("spin", 1) - water.get_potential_energy() + try: + water.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + water.info["energy"] = np.nan for polymorph in polymorphs: polymorph_path = data_dir / polymorph / "POSCAR" @@ -69,7 +75,11 @@ def test_lattice_energy(mlip: tuple[str, Any]) -> None: # Set default charge and spin struct.info.setdefault("charge", 0) struct.info.setdefault("spin", 1) - struct.get_potential_energy() + try: + struct.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {polymorph}: {exc}", stacklevel=2) + struct.info["energy"] = np.nan struct.info["ref"] = ref struct.info["polymorph"] = polymorph diff --git a/ml_peg/calcs/molecular_crystal/X23/calc_X23.py b/ml_peg/calcs/molecular_crystal/X23/calc_X23.py index df5aa92b5..d9680e881 100644 --- a/ml_peg/calcs/molecular_crystal/X23/calc_X23.py +++ b/ml_peg/calcs/molecular_crystal/X23/calc_X23.py @@ -5,6 +5,7 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write @@ -63,14 +64,25 @@ def test_lattice_energy(mlip: tuple[str, Any]) -> None: # Set default charge and spin molecule.info.setdefault("charge", 0) molecule.info.setdefault("spin", 1) - molecule.get_potential_energy() + try: + molecule.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {system}: {exc}", + stacklevel=2, + ) + molecule.info["energy"] = np.nan solid = read(solid_path, index=0, format="vasp") solid.calc = copy(calc) # Set default charge and spin solid.info.setdefault("charge", 0) solid.info.setdefault("spin", 1) - solid.get_potential_energy() + try: + solid.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {system}: {exc}", stacklevel=2) + solid.info["energy"] = np.nan ref = np.loadtxt(ref_path)[0] num_molecules = np.loadtxt(num_molecules_path) diff --git a/ml_peg/calcs/molecular_dynamics/liquid_densities/calc_liquid_densities.py b/ml_peg/calcs/molecular_dynamics/liquid_densities/calc_liquid_densities.py index e548e46aa..7d6215ce0 100644 --- a/ml_peg/calcs/molecular_dynamics/liquid_densities/calc_liquid_densities.py +++ b/ml_peg/calcs/molecular_dynamics/liquid_densities/calc_liquid_densities.py @@ -11,10 +11,12 @@ from pathlib import Path import time from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import Trajectory, read from ase.md.nose_hoover_chain import IsotropicMTKNPT +import numpy as np import pytest from ml_peg.calcs.utils.utils import download_s3_data @@ -121,7 +123,11 @@ def run_npt(atoms, calc, output_fname): ) dyn.nsteps = nsteps dyn.attach(log_md, interval=LOG_INTERVAL, dyn=dyn, start_time=time.time()) - dyn.run(steps=NUM_MD_STEPS - nsteps) + try: + dyn.run(steps=NUM_MD_STEPS - nsteps) + except Exception as exc: + warn(f"Error running MD: {exc}", stacklevel=2) + dyn.atoms.info["energy"] = np.nan @pytest.mark.very_slow diff --git a/ml_peg/calcs/molecular_dynamics/water_density/calc_water_density.py b/ml_peg/calcs/molecular_dynamics/water_density/calc_water_density.py index 6f92b5031..daed81224 100644 --- a/ml_peg/calcs/molecular_dynamics/water_density/calc_water_density.py +++ b/ml_peg/calcs/molecular_dynamics/water_density/calc_water_density.py @@ -7,10 +7,12 @@ from pathlib import Path import time from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import Trajectory, read from ase.md.nose_hoover_chain import IsotropicMTKNPT +import numpy as np import pytest from ml_peg.calcs.utils.utils import download_s3_data @@ -117,7 +119,11 @@ def run_npt(atoms, calc, output_fname, temperature): ) dyn.nsteps = nsteps dyn.attach(log_md, interval=LOG_INTERVAL, dyn=dyn, start_time=time.time()) - dyn.run(steps=NUM_MD_STEPS - nsteps) + try: + dyn.run(steps=NUM_MD_STEPS - nsteps) + except Exception as exc: + warn(f"Error running MD: {exc}", stacklevel=2) + dyn.atoms.info["energy"] = np.nan @pytest.mark.very_slow diff --git a/ml_peg/calcs/molecular_reactions/BH2O_36/calc_BH2O_36.py b/ml_peg/calcs/molecular_reactions/BH2O_36/calc_BH2O_36.py index 7f02fbd64..acf107a69 100644 --- a/ml_peg/calcs/molecular_reactions/BH2O_36/calc_BH2O_36.py +++ b/ml_peg/calcs/molecular_reactions/BH2O_36/calc_BH2O_36.py @@ -10,8 +10,10 @@ import json from pathlib import Path from typing import Any +from warnings import warn from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -106,21 +108,42 @@ def test_bh2o_36(mlip: tuple[str, Any]) -> None: atoms_rct.info["spin"] = 1 atoms_rct.info["ref_energy"] = system["rct"]["energy"] atoms_rct.calc = calc - atoms_rct.info["pred_energy"] = atoms_rct.get_potential_energy() + try: + atoms_rct.info["pred_energy"] = atoms_rct.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {identifier}: {exc}", + stacklevel=2, + ) + atoms_rct.info["pred_energy"] = np.nan atoms_pro = read(system["pro"]["xyz_path"]) atoms_pro.info["charge"] = int(system["pro"]["charge"]) atoms_pro.info["spin"] = 1 atoms_pro.info["ref_energy"] = system["pro"]["energy"] atoms_pro.calc = calc - atoms_pro.info["pred_energy"] = atoms_pro.get_potential_energy() + try: + atoms_pro.info["pred_energy"] = atoms_pro.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {identifier}: {exc}", + stacklevel=2, + ) + atoms_pro.info["pred_energy"] = np.nan atoms_ts = read(system["ts"]["xyz_path"]) atoms_ts.info["charge"] = int(system["ts"]["charge"]) atoms_ts.info["spin"] = 1 atoms_ts.info["ref_energy"] = system["ts"]["energy"] atoms_ts.calc = calc - atoms_ts.info["pred_energy"] = atoms_ts.get_potential_energy() + try: + atoms_ts.info["pred_energy"] = atoms_ts.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {identifier}: {exc}", + stacklevel=2, + ) + atoms_ts.info["pred_energy"] = np.nan write_dir = OUT_PATH / model_name write_dir.mkdir(parents=True, exist_ok=True) diff --git a/ml_peg/calcs/molecular_reactions/BH9/calc_BH9.py b/ml_peg/calcs/molecular_reactions/BH9/calc_BH9.py index d49a7e300..5a68e23b6 100644 --- a/ml_peg/calcs/molecular_reactions/BH9/calc_BH9.py +++ b/ml_peg/calcs/molecular_reactions/BH9/calc_BH9.py @@ -10,9 +10,11 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -163,7 +165,11 @@ def test_bh9(mlip: tuple[str, Any]) -> None: # Create list for TS and reactant atoms structs = [process_atoms(xyz_path / f"{label}TS.xyz")] structs[0].calc = calc - structs[0].info["model_energy"] = structs[0].get_potential_energy() + try: + structs[0].info["model_energy"] = structs[0].get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + structs[0].info["model_energy"] = np.nan structs[0].info["label"] = label # Write both forward and reverse barriers, only forward used in analysis here @@ -173,7 +179,16 @@ def test_bh9(mlip: tuple[str, Any]) -> None: for file in xyz_path.glob(f"{label}R*.xyz"): reactant_atoms = process_atoms(file) reactant_atoms.calc = copy(calc) - reactant_atoms.info["model_energy"] = reactant_atoms.get_potential_energy() + try: + reactant_atoms.info["model_energy"] = ( + reactant_atoms.get_potential_energy() + ) + except Exception as exc: + warn( + f"Error calculating energy for {file}: {exc}", + stacklevel=2, + ) + reactant_atoms.info["model_energy"] = np.nan reactant_atoms.info["label"] = label structs.append(reactant_atoms) diff --git a/ml_peg/calcs/molecular_reactions/CYCLO70/calc_CYCLO70.py b/ml_peg/calcs/molecular_reactions/CYCLO70/calc_CYCLO70.py index 2423bdac2..fde267d8a 100644 --- a/ml_peg/calcs/molecular_reactions/CYCLO70/calc_CYCLO70.py +++ b/ml_peg/calcs/molecular_reactions/CYCLO70/calc_CYCLO70.py @@ -12,9 +12,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -94,7 +96,15 @@ def test_cyclo70(mlip: tuple[str, Any]) -> None: atoms.info["charge"] = int(atoms.info["charge"]) else: atoms.info["charge"] = 0 - bh_forward_model -= atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {atoms_label} in {rxn}: {exc}", + stacklevel=2, + ) + energy = np.nan + bh_forward_model -= energy atoms.info["label"] = atoms_label atoms.calc = None structs_forward.append(atoms) @@ -110,7 +120,15 @@ def test_cyclo70(mlip: tuple[str, Any]) -> None: atoms.info["charge"] = int(atoms.info["charge"]) else: atoms.info["charge"] = 0 - bh_reverse_model -= atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {atoms_label} in {rxn}: {exc}", + stacklevel=2, + ) + energy = np.nan + bh_reverse_model -= energy atoms.info["label"] = atoms_label atoms.calc = None structs_reverse.append(atoms) @@ -126,8 +144,16 @@ def test_cyclo70(mlip: tuple[str, Any]) -> None: atoms.info["charge"] = int(atoms.info["charge"]) else: atoms.info["charge"] = 0 - bh_forward_model += atoms.get_potential_energy() - bh_reverse_model += atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {atoms_label} in {rxn}: {exc}", + stacklevel=2, + ) + energy = np.nan + bh_forward_model += energy + bh_reverse_model += energy atoms.info["ref_forward_bh"] = bh_forward_ref atoms.info["ref_reverse_bh"] = bh_reverse_ref diff --git a/ml_peg/calcs/molecular_reactions/Criegee22/calc_Criegee22.py b/ml_peg/calcs/molecular_reactions/Criegee22/calc_Criegee22.py index a7fbc1ca8..93f3b59d4 100644 --- a/ml_peg/calcs/molecular_reactions/Criegee22/calc_Criegee22.py +++ b/ml_peg/calcs/molecular_reactions/Criegee22/calc_Criegee22.py @@ -9,9 +9,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -61,14 +63,28 @@ def test_criegee22(mlip: tuple[str, Any]) -> None: atoms_reac.calc = calc atoms_reac.info["charge"] = 0 atoms_reac.info["spin"] = 1 - atoms_reac.info["model_energy"] = atoms_reac.get_potential_energy() + try: + atoms_reac.info["model_energy"] = atoms_reac.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating model energy for {label}: {exc}", + stacklevel=2, + ) + atoms_reac.info["model_energy"] = np.nan atoms_reac.info["ref_energy"] = 0 atoms_ts = read(data_path / "structures" / f"{label}-TS.xyz") atoms_ts.calc = calc atoms_ts.info["charge"] = 0 atoms_ts.info["spin"] = 1 - atoms_ts.info["model_energy"] = atoms_ts.get_potential_energy() + try: + atoms_ts.info["model_energy"] = atoms_ts.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating model energy for {label}: {exc}", + stacklevel=2, + ) + atoms_ts.info["model_energy"] = np.nan atoms_ts.info["ref_energy"] = bh_ref write_dir = OUT_PATH / model_name diff --git a/ml_peg/calcs/molecular_reactions/RDB7/calc_RDB7.py b/ml_peg/calcs/molecular_reactions/RDB7/calc_RDB7.py index 06c2c6371..9c2e1f476 100644 --- a/ml_peg/calcs/molecular_reactions/RDB7/calc_RDB7.py +++ b/ml_peg/calcs/molecular_reactions/RDB7/calc_RDB7.py @@ -12,6 +12,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atom, Atoms, units from ase.io import write @@ -118,12 +119,25 @@ def test_rdb87(mlip: tuple[str, Any]) -> None: bh_forward_ref -= get_cc_energy(qm_path) atoms = get_atoms_from_molpro(qm_path) atoms.calc = calc - bh_forward_model -= atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for {qm_path}: {exc}", + stacklevel=2, + ) + energy = np.nan + bh_forward_model -= energy for qm_path in (data_path / "qm_logs" / f"rxn{label}").glob("ts*"): bh_forward_ref += get_cc_energy(qm_path) atoms = get_atoms_from_molpro(qm_path) atoms.calc = calc - bh_forward_model += atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {qm_path}: {exc}", stacklevel=2) + energy = np.nan + bh_forward_model += energy atoms.info["model_forward_barrier"] = bh_forward_model atoms.info["ref_forward_barrier"] = bh_forward_ref diff --git a/ml_peg/calcs/nebs/li_diffusion/calc_li_diffusion.py b/ml_peg/calcs/nebs/li_diffusion/calc_li_diffusion.py index 5fad0ebc3..bded9ec8b 100644 --- a/ml_peg/calcs/nebs/li_diffusion/calc_li_diffusion.py +++ b/ml_peg/calcs/nebs/li_diffusion/calc_li_diffusion.py @@ -3,11 +3,13 @@ from __future__ import annotations from pathlib import Path +from warnings import warn from ase import Atoms -from ase.io import read +from ase.io import read, write from janus_core.calculations.geom_opt import GeomOpt from janus_core.calculations.neb import NEB +import numpy as np import pytest from ml_peg.models import current_models @@ -43,11 +45,19 @@ def relaxed_structs() -> dict[str, Atoms]: geomopt = GeomOpt( struct=struct, write_results=True, - file_prefix=OUT_PATH / f"{struct_name}-{model_name}", + file_prefix=OUT_PATH / model_name / struct_name, filter_class=None, + steps=25, ) - geomopt.run() - relaxed_structs[f"{struct_name}-{model_name}"] = geomopt.struct + try: + geomopt.run() + except Exception as exc: + warn( + f"Error for {struct_name} with {model_name}: {exc}", + stacklevel=2, + ) + struct.info["energy"] = np.nan + relaxed_structs[f"{model_name}-{struct_name}"] = geomopt.struct return relaxed_structs @@ -64,23 +74,38 @@ def test_li_diffusion_b(relaxed_structs: dict[str, Atoms], model_name: str) -> N model_name Name of model to use. """ + init_struct = relaxed_structs[f"{model_name}-LiFePO4_start_bc.cif"] + final_struct = relaxed_structs[f"{model_name}-LiFePO4_end_b.cif"] neb = NEB( - init_struct=relaxed_structs[f"LiFePO4_start_bc.cif-{model_name}"], - final_struct=relaxed_structs[f"LiFePO4_end_b.cif-{model_name}"], + init_struct=init_struct, + final_struct=final_struct, n_images=11, interpolator="pymatgen", minimize=True, plot_band=True, write_band=True, - file_prefix=OUT_PATH / f"li_diffusion_b-{model_name}", + file_prefix=OUT_PATH / model_name / "li_diffusion_b", + steps=100, ) # Set default charge and spin for all images - neb.interpolate() - neb.interpolator = None - for image in neb.images: - image.info.setdefault("charge", 0) - image.info.setdefault("spin", 1) - neb.run() + try: + neb.interpolate() + neb.interpolator = None + for image in neb.images: + image.info.setdefault("charge", 0) + image.info.setdefault("spin", 1) + neb.run() + except Exception as e: + print(f"Error running NEB for {model_name} path B: {e}") + # Write out equivalent data + out_dir = OUT_PATH / model_name + out_dir.mkdir(exist_ok=True, parents=True) + write(out_dir / "li_diffusion_b-neb-band.extxyz", [init_struct, final_struct]) + with open( + out_dir / "li_diffusion_b-neb-results.dat", "w", encoding="utf8" + ) as out: + print("#Barrier [eV] | delta E [eV] | Max force [eV/Å] ", file=out) + print(np.nan, np.nan, np.nan, file=out) @pytest.mark.slow @@ -96,20 +121,35 @@ def test_li_diffusion_c(relaxed_structs: dict[str, Atoms], model_name: str) -> N model_name Name of model to use. """ + init_struct = relaxed_structs[f"{model_name}-LiFePO4_start_bc.cif"] + final_struct = relaxed_structs[f"{model_name}-LiFePO4_end_c.cif"] neb = NEB( - init_struct=relaxed_structs[f"LiFePO4_start_bc.cif-{model_name}"], - final_struct=relaxed_structs[f"LiFePO4_end_c.cif-{model_name}"], + init_struct=init_struct, + final_struct=final_struct, n_images=11, interpolator="pymatgen", minimize=True, plot_band=True, write_band=True, - file_prefix=OUT_PATH / f"li_diffusion_c-{model_name}", + file_prefix=OUT_PATH / model_name / "li_diffusion_c", + steps=500, ) # Set default charge and spin for all images - neb.interpolate() - neb.interpolator = None - for image in neb.images: - image.info.setdefault("charge", 0) - image.info.setdefault("spin", 1) - neb.run() + try: + neb.interpolate() + neb.interpolator = None + for image in neb.images: + image.info.setdefault("charge", 0) + image.info.setdefault("spin", 1) + neb.run() + except Exception as e: + print(f"Error running NEB for {model_name} path C: {e}") + # Write out equivalent data + out_dir = OUT_PATH / model_name + out_dir.mkdir(exist_ok=True, parents=True) + write(out_dir / "li_diffusion_c-neb-band.extxyz", [init_struct, final_struct]) + with open( + out_dir / "li_diffusion_c-neb-results.dat", "w", encoding="utf8" + ) as out: + print("#Barrier [eV] | delta E [eV] | Max force [eV/Å] ", file=out) + print(np.nan, np.nan, np.nan, file=out) diff --git a/ml_peg/calcs/nebs/si_defects/calc_si_defects.py b/ml_peg/calcs/nebs/si_defects/calc_si_defects.py index e77ceff34..4d8ea2f0e 100644 --- a/ml_peg/calcs/nebs/si_defects/calc_si_defects.py +++ b/ml_peg/calcs/nebs/si_defects/calc_si_defects.py @@ -6,6 +6,7 @@ from pathlib import Path import sys from typing import Any +from warnings import warn from ase.atoms import Atoms from ase.io import read, write @@ -151,10 +152,18 @@ def test_si_defects(mlip: tuple[str, Any]) -> None: out_atoms = atoms.copy() out_atoms.info["ref_energy_ev"] = ref_energy_ev out_atoms.arrays["ref_forces"] = ref_forces - out_atoms.info["pred_energy_ev"] = float(atoms_pred.get_potential_energy()) - out_atoms.arrays["pred_forces"] = np.asarray( - atoms_pred.get_forces(), dtype=float - ) + try: + pred_energy = float(atoms_pred.get_potential_energy()) + pred_forces = np.asarray(atoms_pred.get_forces(), dtype=float) + except Exception as exc: + warn( + f"Error calculating energy for {case}: {exc}", + stacklevel=2, + ) + pred_energy = np.nan + pred_forces = np.full((len(atoms_pred), 3), np.nan) + out_atoms.info["pred_energy_ev"] = pred_energy + out_atoms.arrays["pred_forces"] = pred_forces results.append(out_atoms) write(out_dir / "si_defects.extxyz", results) diff --git a/ml_peg/calcs/non_covalent_interactions/IONPI19/calc_IONPI19.py b/ml_peg/calcs/non_covalent_interactions/IONPI19/calc_IONPI19.py index f23c27240..34a25af7e 100644 --- a/ml_peg/calcs/non_covalent_interactions/IONPI19/calc_IONPI19.py +++ b/ml_peg/calcs/non_covalent_interactions/IONPI19/calc_IONPI19.py @@ -9,9 +9,11 @@ import os from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -135,7 +137,11 @@ def test_ionpi19(mlip: tuple[str, Any]) -> None: atoms = get_atoms(data_path / label) atoms.info["ref_int_energy"] = ref_energies[system_id] atoms.calc = calc - atoms.info["model_energy"] = atoms.get_potential_energy() + try: + atoms.info["model_energy"] = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_energy"] = np.nan write_dir = OUT_PATH / model_name write_dir.mkdir(parents=True, exist_ok=True) diff --git a/ml_peg/calcs/non_covalent_interactions/NCIA_D1200/calc_NCIA_D1200.py b/ml_peg/calcs/non_covalent_interactions/NCIA_D1200/calc_NCIA_D1200.py index bdf56367f..e3f3b1ae0 100644 --- a/ml_peg/calcs/non_covalent_interactions/NCIA_D1200/calc_NCIA_D1200.py +++ b/ml_peg/calcs/non_covalent_interactions/NCIA_D1200/calc_NCIA_D1200.py @@ -9,6 +9,7 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write @@ -135,11 +136,15 @@ def test_ncia_d1200(mlip: tuple[str, Any]) -> None: atoms_a.calc = copy(calc) atoms_b.calc = copy(calc) - atoms.info["model_int_energy"] = ( - atoms.get_potential_energy() - - atoms_a.get_potential_energy() - - atoms_b.get_potential_energy() - ) + try: + atoms.info["model_int_energy"] = ( + atoms.get_potential_energy() + - atoms_a.get_potential_energy() + - atoms_b.get_potential_energy() + ) + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_int_energy"] = np.nan atoms.info["ref_int_energy"] = ref_energy atoms.calc = None diff --git a/ml_peg/calcs/non_covalent_interactions/NCIA_D442x10/calc_NCIA_D442x10.py b/ml_peg/calcs/non_covalent_interactions/NCIA_D442x10/calc_NCIA_D442x10.py index 8f1529b48..b46178fae 100644 --- a/ml_peg/calcs/non_covalent_interactions/NCIA_D442x10/calc_NCIA_D442x10.py +++ b/ml_peg/calcs/non_covalent_interactions/NCIA_D442x10/calc_NCIA_D442x10.py @@ -9,6 +9,7 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write @@ -133,11 +134,15 @@ def test_ncia_d442x10(mlip: tuple[str, Any]) -> None: atoms_a.calc = copy(calc) atoms_b.calc = copy(calc) - atoms.info["model_int_energy"] = ( - atoms.get_potential_energy() - - atoms_a.get_potential_energy() - - atoms_b.get_potential_energy() - ) + try: + atoms.info["model_int_energy"] = ( + atoms.get_potential_energy() + - atoms_a.get_potential_energy() + - atoms_b.get_potential_energy() + ) + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_int_energy"] = np.nan atoms.info["ref_int_energy"] = ref_energy atoms.calc = None diff --git a/ml_peg/calcs/non_covalent_interactions/NCIA_HB300SPXx10/calc_NCIA_HB300SPXx10.py b/ml_peg/calcs/non_covalent_interactions/NCIA_HB300SPXx10/calc_NCIA_HB300SPXx10.py index 7ced2ae54..c99bb4f3b 100644 --- a/ml_peg/calcs/non_covalent_interactions/NCIA_HB300SPXx10/calc_NCIA_HB300SPXx10.py +++ b/ml_peg/calcs/non_covalent_interactions/NCIA_HB300SPXx10/calc_NCIA_HB300SPXx10.py @@ -8,9 +8,11 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -126,11 +128,15 @@ def test_ncia_hb300spxx10(mlip: tuple[str, Any]) -> None: atoms_a.calc = calc atoms_b.calc = calc - atoms.info["model_int_energy"] = ( - atoms.get_potential_energy() - - atoms_a.get_potential_energy() - - atoms_b.get_potential_energy() - ) + try: + atoms.info["model_int_energy"] = ( + atoms.get_potential_energy() + - atoms_a.get_potential_energy() + - atoms_b.get_potential_energy() + ) + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_int_energy"] = np.nan atoms.info["ref_int_energy"] = ref_energy atoms.calc = None diff --git a/ml_peg/calcs/non_covalent_interactions/NCIA_HB375x10/calc_NCIA_HB375x10.py b/ml_peg/calcs/non_covalent_interactions/NCIA_HB375x10/calc_NCIA_HB375x10.py index 459b1c967..505f3ce2c 100644 --- a/ml_peg/calcs/non_covalent_interactions/NCIA_HB375x10/calc_NCIA_HB375x10.py +++ b/ml_peg/calcs/non_covalent_interactions/NCIA_HB375x10/calc_NCIA_HB375x10.py @@ -9,9 +9,11 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -132,11 +134,15 @@ def test_ncia_hb375x10(mlip: tuple[str, Any]) -> None: atoms_a.calc = copy(calc) atoms_b.calc = copy(calc) - atoms.info["model_int_energy"] = ( - atoms.get_potential_energy() - - atoms_a.get_potential_energy() - - atoms_b.get_potential_energy() - ) + try: + atoms.info["model_int_energy"] = ( + atoms.get_potential_energy() + - atoms_a.get_potential_energy() + - atoms_b.get_potential_energy() + ) + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_int_energy"] = np.nan atoms.info["ref_int_energy"] = ref_energy atoms.calc = None diff --git a/ml_peg/calcs/non_covalent_interactions/NCIA_IHB100x10/calc_NCIA_IHB100x10.py b/ml_peg/calcs/non_covalent_interactions/NCIA_IHB100x10/calc_NCIA_IHB100x10.py index c067fbe1e..eb44b1653 100644 --- a/ml_peg/calcs/non_covalent_interactions/NCIA_IHB100x10/calc_NCIA_IHB100x10.py +++ b/ml_peg/calcs/non_covalent_interactions/NCIA_IHB100x10/calc_NCIA_IHB100x10.py @@ -9,9 +9,11 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -129,11 +131,15 @@ def test_ncia_ihb100x10(mlip: tuple[str, Any]) -> None: atoms_a.calc = copy(calc) atoms_b.calc = copy(calc) - atoms.info["model_int_energy"] = ( - atoms.get_potential_energy() - - atoms_a.get_potential_energy() - - atoms_b.get_potential_energy() - ) + try: + atoms.info["model_int_energy"] = ( + atoms.get_potential_energy() + - atoms_a.get_potential_energy() + - atoms_b.get_potential_energy() + ) + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_int_energy"] = np.nan atoms.info["ref_int_energy"] = ref_energy atoms.calc = None diff --git a/ml_peg/calcs/non_covalent_interactions/NCIA_R739x5/calc_NCIA_R739x5.py b/ml_peg/calcs/non_covalent_interactions/NCIA_R739x5/calc_NCIA_R739x5.py index 9e23dd278..372511433 100644 --- a/ml_peg/calcs/non_covalent_interactions/NCIA_R739x5/calc_NCIA_R739x5.py +++ b/ml_peg/calcs/non_covalent_interactions/NCIA_R739x5/calc_NCIA_R739x5.py @@ -8,6 +8,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write @@ -133,11 +134,15 @@ def test_ncia_r739x5(mlip: tuple[str, Any]) -> None: atoms_a.calc = calc atoms_b.calc = calc - atoms.info["model_int_energy"] = ( - atoms.get_potential_energy() - - atoms_a.get_potential_energy() - - atoms_b.get_potential_energy() - ) + try: + atoms.info["model_int_energy"] = ( + atoms.get_potential_energy() + - atoms_a.get_potential_energy() + - atoms_b.get_potential_energy() + ) + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_int_energy"] = np.nan atoms.info["ref_int_energy"] = ref_energy atoms.calc = None diff --git a/ml_peg/calcs/non_covalent_interactions/NCIA_SH250x10/calc_NCIA_SH250x10.py b/ml_peg/calcs/non_covalent_interactions/NCIA_SH250x10/calc_NCIA_SH250x10.py index 590f7ceee..7e997a91e 100644 --- a/ml_peg/calcs/non_covalent_interactions/NCIA_SH250x10/calc_NCIA_SH250x10.py +++ b/ml_peg/calcs/non_covalent_interactions/NCIA_SH250x10/calc_NCIA_SH250x10.py @@ -9,9 +9,11 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -126,11 +128,15 @@ def test_lattice_energy(mlip: tuple[str, Any]) -> None: atoms_a.calc = copy(calc) atoms_b.calc = copy(calc) - atoms.info["model_int_energy"] = ( - atoms.get_potential_energy() - - atoms_a.get_potential_energy() - - atoms_b.get_potential_energy() - ) + try: + atoms.info["model_int_energy"] = ( + atoms.get_potential_energy() + - atoms_a.get_potential_energy() + - atoms_b.get_potential_energy() + ) + except Exception as exc: + warn(f"Error calculating energy for {label}: {exc}", stacklevel=2) + atoms.info["model_int_energy"] = np.nan atoms.info["ref_int_energy"] = ref_energy atoms.calc = None diff --git a/ml_peg/calcs/non_covalent_interactions/QUID/calc_QUID.py b/ml_peg/calcs/non_covalent_interactions/QUID/calc_QUID.py index 9080cc61e..af98e38ea 100644 --- a/ml_peg/calcs/non_covalent_interactions/QUID/calc_QUID.py +++ b/ml_peg/calcs/non_covalent_interactions/QUID/calc_QUID.py @@ -11,6 +11,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms from ase.io import write @@ -104,6 +105,7 @@ def compute_interaction_energy(dataset, label, calc): # List to store dimer and monomers. atoms_list = [] model_int_energy = 0 + for atoms_name, stoich in zip( ["dimer", "small_monomer", "big_monomer"], [1, -1, -1], strict=False ): @@ -112,9 +114,16 @@ def compute_interaction_energy(dataset, label, calc): atoms = Atoms(numbers=atomic_numbers, positions=positions) atoms.info.update({"charge": 0, "spin": 1}) atoms.calc = calc - model_int_energy += atoms.get_potential_energy() * stoich - atoms.calc = None atoms_list.append(atoms) + try: + model_int_energy += atoms.get_potential_energy() * stoich + except Exception as exc: + warn(f"Error calculating energy for {atoms_name}: {exc}", stacklevel=2) + model_int_energy = np.nan + break + + for atoms in atoms_list: + atoms.calc = None atoms_list[0].info.update( { "ref_int_energy": ref_int_energy, diff --git a/ml_peg/calcs/physicality/diatomics/calc_diatomics.py b/ml_peg/calcs/physicality/diatomics/calc_diatomics.py index 6033b7976..be8f6c9e0 100644 --- a/ml_peg/calcs/physicality/diatomics/calc_diatomics.py +++ b/ml_peg/calcs/physicality/diatomics/calc_diatomics.py @@ -6,6 +6,7 @@ import itertools import json from pathlib import Path +from warnings import warn from ase import Atoms from ase.data import chemical_symbols @@ -176,8 +177,16 @@ def run_diatomics(model_name: str, model) -> None: atoms.info.setdefault("spin", 1) atoms.calc = calc - energy = float(atoms.get_potential_energy()) - forces = atoms.get_forces() + try: + energy = float(atoms.get_potential_energy()) + forces = atoms.get_forces() + except Exception as exc: + warn( + f"Error calculating energy for {pair_label}: {exc}", + stacklevel=2, + ) + energy = np.nan + forces = np.full((len(atoms), 3), np.nan) bond_vector = atoms.positions[1] - atoms.positions[0] force_parallel = _project_force(forces, bond_vector) diff --git a/ml_peg/calcs/physicality/extensivity/calc_extensivity.py b/ml_peg/calcs/physicality/extensivity/calc_extensivity.py index 44792a7e9..f94114a8e 100644 --- a/ml_peg/calcs/physicality/extensivity/calc_extensivity.py +++ b/ml_peg/calcs/physicality/extensivity/calc_extensivity.py @@ -5,10 +5,12 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms from ase.build import bulk, surface from ase.io import write +import numpy as np import pytest from ml_peg.models import current_models @@ -87,14 +89,15 @@ def test_extensivity(mlip: tuple[str, Any]) -> None: slab2_big = slab2.copy() slab2_big.set_cell(tall_cell, scale_atoms=False) - slab1_big.calc = calc - slab1_big.get_potential_energy() - - slab2_big.calc = copy(calc) - slab2_big.get_potential_energy() - - combined.calc = copy(calc) - combined.get_potential_energy() + for struct in (slab1_big, slab2_big, combined): + struct.calc = copy(calc) + try: + energy = struct.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + energy = np.nan + struct.info["energy"] = energy + struct.calc = None # Write output structures write_dir = OUT_PATH / model_name diff --git a/ml_peg/calcs/physicality/locality/calc_locality.py b/ml_peg/calcs/physicality/locality/calc_locality.py index 419cf3169..a89159bda 100644 --- a/ml_peg/calcs/physicality/locality/calc_locality.py +++ b/ml_peg/calcs/physicality/locality/calc_locality.py @@ -4,6 +4,7 @@ from copy import copy from pathlib import Path +from warnings import warn from ase import Atoms from ase.io import write @@ -152,12 +153,12 @@ def prepared_solute() -> dict[str, Atoms]: solutes = {} for model_name, calc in MODELS.items(): solute = solute.copy() + solute.calc = calc.get_calculator(precision="high") + solutes[model_name] = solute try: - solute.calc = calc.get_calculator(precision="high") solute.get_forces() - solutes[model_name] = solute - # If a model fails, don't block other model tests - except (ModuleNotFoundError, RuntimeError, TypeError): + except Exception as exc: + warn(f"Error preparing solute: {exc} for {model_name}", stacklevel=2) continue return solutes @@ -182,7 +183,11 @@ def test_ghost_atoms(prepared_solute: dict[str, Atoms], model_name: str) -> None system_ghost = prepare_ghost_system(solute, ghost_num, ghost_dist) system_ghost.calc = copy(solute.calc) - system_ghost.get_forces() + try: + system_ghost.get_forces() + except Exception as exc: + warn(f"Error calculating forces: {exc}", stacklevel=2) + system_ghost.arrays["forces"] = np.full((len(system_ghost), 3), np.nan) # Write output structures write_dir = OUT_PATH / model_name @@ -214,7 +219,11 @@ def test_rand_h(prepared_solute: dict[str, Atoms], model_name: str) -> None: for _ in range(rand_trials): system_rand_h = add_random_h(solute, rand_min_dist, rand_max_dist, rng) system_rand_h.calc = copy(solute.calc) - system_rand_h.get_forces() + try: + system_rand_h.get_forces() + except Exception as exc: + warn(f"Error calculating forces: {exc}", stacklevel=2) + system_rand_h.arrays["forces"] = np.full((len(system_rand_h), 3), np.nan) rand_h_structures.append(system_rand_h) # Write output structures diff --git a/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py b/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py index 9bf9101c9..4a825ab18 100644 --- a/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py +++ b/ml_peg/calcs/physicality/oxidation_states/calc_oxidation_states.py @@ -4,6 +4,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase.geometry.rdf import get_rdf from ase.io import read @@ -67,7 +68,10 @@ def test_iron_oxidation_state_md(mlip: tuple[str, Any]) -> None: file_prefix=out_dir / f"{salt}_{model_name}", restart_every=1000, ) - npt.run() + try: + npt.run() + except Exception as exc: + warn(f"Error during MD for {salt}: {exc}", stacklevel=2) @pytest.mark.parametrize("model_name", MODELS) diff --git a/ml_peg/calcs/physicality/water_slab_dipoles/calc_water_slab_dipoles.py b/ml_peg/calcs/physicality/water_slab_dipoles/calc_water_slab_dipoles.py index f81e78c59..3f6722aef 100644 --- a/ml_peg/calcs/physicality/water_slab_dipoles/calc_water_slab_dipoles.py +++ b/ml_peg/calcs/physicality/water_slab_dipoles/calc_water_slab_dipoles.py @@ -4,6 +4,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read @@ -84,6 +85,9 @@ def test_water_dipole(mlip: tuple[str, Any]) -> None: print("Starting MD") - md.run() + try: + md.run() + except Exception as exc: + warn(f"Error during MD: {exc}", stacklevel=2) print("Finished MD") diff --git a/ml_peg/calcs/supramolecular/LNCI16/calc_LNCI16.py b/ml_peg/calcs/supramolecular/LNCI16/calc_LNCI16.py index 2a69deafa..1f39c7671 100644 --- a/ml_peg/calcs/supramolecular/LNCI16/calc_LNCI16.py +++ b/ml_peg/calcs/supramolecular/LNCI16/calc_LNCI16.py @@ -3,12 +3,14 @@ from __future__ import annotations from pathlib import Path +from warnings import warn from ase import Atoms, units from ase.calculators.calculator import Calculator from ase.io import read, write import mlipx from mlipx.abc import NodeWithCalculator +import numpy as np from tqdm import tqdm import zntrack @@ -177,12 +179,16 @@ def interaction_energy(frags: dict[str, Atoms], calc: Calculator) -> float: guest_copy = frags["guest"].copy() complex_copy.calc = calc - e_complex = complex_copy.get_potential_energy() host_copy.calc = calc - e_host = host_copy.get_potential_energy() guest_copy.calc = calc - e_guest = guest_copy.get_potential_energy() - return e_complex - e_host - e_guest + try: + e_complex = complex_copy.get_potential_energy() + e_host = host_copy.get_potential_energy() + e_guest = guest_copy.get_potential_energy() + return e_complex - e_host - e_guest + except Exception as exc: + warn(f"Error calculating energies: {exc}", stacklevel=2) + return np.nan @staticmethod def benchmark_lnci16( diff --git a/ml_peg/calcs/supramolecular/S30L/calc_S30L.py b/ml_peg/calcs/supramolecular/S30L/calc_S30L.py index ba899644d..e51ac88ec 100644 --- a/ml_peg/calcs/supramolecular/S30L/calc_S30L.py +++ b/ml_peg/calcs/supramolecular/S30L/calc_S30L.py @@ -3,13 +3,14 @@ from __future__ import annotations from pathlib import Path -import warnings +from warnings import warn from ase import Atoms, units from ase.calculators.calculator import Calculator from ase.io import read, write import mlipx from mlipx.abc import NodeWithCalculator +import numpy as np from tqdm import tqdm import zntrack @@ -59,9 +60,7 @@ def read_charge(folder: Path) -> int: try: return int(f.read_text().strip()) except ValueError: - warnings.warn( - f"Invalid charge in {f} - assuming neutral.", stacklevel=2 - ) + warn(f"Invalid charge in {f} - assuming neutral.", stacklevel=2) return 0 def read_atoms(self, folder: Path, ident: str) -> Atoms: @@ -134,12 +133,16 @@ def interaction_energy(frags: dict[str, Atoms], calc: Calculator) -> float: Interaction energy in eV. """ frags["complex"].calc = calc - e_complex = frags["complex"].get_potential_energy() frags["host"].calc = calc - e_host = frags["host"].get_potential_energy() frags["guest"].calc = calc - e_guest = frags["guest"].get_potential_energy() - return e_complex - e_host - e_guest + try: + e_complex = frags["complex"].get_potential_energy() + e_host = frags["host"].get_potential_energy() + e_guest = frags["guest"].get_potential_energy() + return e_complex - e_host - e_guest + except Exception as exc: + warn(f"Error calculating energies: {exc}", stacklevel=2) + return np.nan @staticmethod def parse_references(path: Path) -> dict[int, float]: @@ -218,8 +221,8 @@ def benchmark_s30l( complex_atoms_list.append(complex_atoms) - except Exception as e: - print(f"Error processing system {idx}: {e}") + except Exception as exc: + warn(f"Error processing system {idx}: {exc}", stacklevel=2) continue return complex_atoms_list diff --git a/ml_peg/calcs/supramolecular/utils/plf547_pla15_utils.py b/ml_peg/calcs/supramolecular/utils/plf547_pla15_utils.py index afecdc7d8..14a6f0d2a 100644 --- a/ml_peg/calcs/supramolecular/utils/plf547_pla15_utils.py +++ b/ml_peg/calcs/supramolecular/utils/plf547_pla15_utils.py @@ -4,6 +4,7 @@ import logging from pathlib import Path +from warnings import warn from ase import Atoms, units from ase.calculators.calculator import Calculator @@ -122,18 +123,22 @@ def get_interaction_energy(fragments: dict[str, Atoms], calc: Calculator) -> flo Interaction energy in eV. """ fragments["complex"].calc = calc - e_complex = fragments["complex"].get_potential_energy() - fragments["complex"].calc = None - fragments["protein"].calc = calc - e_protein = fragments["protein"].get_potential_energy() - fragments["protein"].calc = None - fragments["ligand"].calc = calc - e_ligand = fragments["ligand"].get_potential_energy() - fragments["ligand"].calc = None - - return e_complex - e_protein - e_ligand + try: + e_complex = fragments["complex"].get_potential_energy() + e_protein = fragments["protein"].get_potential_energy() + e_ligand = fragments["ligand"].get_potential_energy() + result = e_complex - e_protein - e_ligand + except Exception as exc: + warn(f"Error calculating energies: {exc}", stacklevel=2) + result = np.nan + finally: + fragments["complex"].calc = None + fragments["protein"].calc = None + fragments["ligand"].calc = None + + return result def mda_atoms_to_ase(atom_list, charge: float, identifier: str) -> Atoms: diff --git a/ml_peg/calcs/surfaces/CMRAds200/calc_CMRAds200.py b/ml_peg/calcs/surfaces/CMRAds200/calc_CMRAds200.py index b1a0eb1a6..b65d84283 100644 --- a/ml_peg/calcs/surfaces/CMRAds200/calc_CMRAds200.py +++ b/ml_peg/calcs/surfaces/CMRAds200/calc_CMRAds200.py @@ -5,6 +5,7 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase.io import read, write import numpy as np @@ -84,10 +85,18 @@ def test_cmrads200(mlip: tuple[str, Any]) -> None: stoichiometry_list = np.array( cmrads_reactions_dict[reaction_id]["Stoichiometry"].split(",") ).astype(float) - systems_energies = np.array( - [structs_energy_dict[sys].get_potential_energy() for sys in systems_list] - ) - reaction_energy = np.sum(systems_energies * stoichiometry_list) + try: + systems_energies = np.array( + [ + structs_energy_dict[sys].get_potential_energy() + for sys in systems_list + ] + ) + reaction_energy = np.sum(systems_energies * stoichiometry_list) + except Exception as exc: + warn(f"Energy calculation failed for {reaction_id}: {exc}", stacklevel=2) + reaction_energy = np.nan + mol_surface = structs_energy_dict[systems_list[0]].copy() mol_surface.info["pred_adsorption_energy"] = reaction_energy mol_surface.info["PBE_adsorption_energy"] = cmrads_reactions_dict[reaction_id][ diff --git a/ml_peg/calcs/surfaces/OC157/calc_OC157.py b/ml_peg/calcs/surfaces/OC157/calc_OC157.py index 9110c5e42..fb74d6e84 100644 --- a/ml_peg/calcs/surfaces/OC157/calc_OC157.py +++ b/ml_peg/calcs/surfaces/OC157/calc_OC157.py @@ -4,12 +4,14 @@ from copy import copy from pathlib import Path +from warnings import warn from ase import Atoms from ase.calculators.calculator import Calculator from ase.io import read, write import mlipx from mlipx.abc import NodeWithCalculator +import numpy as np from tqdm import tqdm import zntrack @@ -123,7 +125,11 @@ def evaluate_energies(triplet: list[Atoms], calc: Calculator) -> None: atoms.info.setdefault("charge", 0) atoms.info.setdefault("spin", 1) - atoms.get_potential_energy() + try: + atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + atoms.info["energy"] = np.nan def run(self): """Run OC157 energy calculations.""" diff --git a/ml_peg/calcs/surfaces/S24/calc_S24.py b/ml_peg/calcs/surfaces/S24/calc_S24.py index 5709ad060..0f948d5a1 100644 --- a/ml_peg/calcs/surfaces/S24/calc_S24.py +++ b/ml_peg/calcs/surfaces/S24/calc_S24.py @@ -4,12 +4,14 @@ from copy import copy from pathlib import Path +from warnings import warn from ase import Atoms from ase.calculators.calculator import Calculator from ase.io import read, write import mlipx from mlipx.abc import NodeWithCalculator +import numpy as np from tqdm import tqdm import zntrack @@ -114,7 +116,11 @@ def evaluate_energies(atoms_list: list[Atoms], calc: Calculator) -> None: atoms.info.setdefault("charge", 0) atoms.info.setdefault("spin", 1) - atoms.get_potential_energy() + try: + atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + atoms.info["energy"] = np.nan def run(self): """Run S24 energy calculations.""" @@ -140,9 +146,21 @@ def run(self): sys_id = f"{(i // 3) + 1:03d}" # Store reference energies - surface.info["ref_energy"] = surface.get_potential_energy() - mol_surface.info["ref_energy"] = mol_surface.get_potential_energy() - molecule.info["ref_energy"] = molecule.get_potential_energy() + try: + surface.info["ref_energy"] = surface.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {sys_id}: {exc}", stacklevel=2) + surface.info["ref_energy"] = np.nan + try: + mol_surface.info["ref_energy"] = mol_surface.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {sys_id}: {exc}", stacklevel=2) + mol_surface.info["ref_energy"] = np.nan + try: + molecule.info["ref_energy"] = molecule.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {sys_id}: {exc}", stacklevel=2) + molecule.info["ref_energy"] = np.nan # Store system information surface_formula = surface.get_chemical_formula() @@ -166,12 +184,16 @@ def run(self): self.evaluate_energies(triplet, calc) # Calculate and store adsorption energies - surface_energy = surface.get_potential_energy() - mol_surf_energy = mol_surface.get_potential_energy() - molecule_energy = molecule.get_potential_energy() - pred_ads_energy = self.compute_adsorption_energy( - surface_energy, mol_surf_energy, molecule_energy, adsorbate_count - ) + try: + surface_energy = surface.get_potential_energy() + mol_surf_energy = mol_surface.get_potential_energy() + molecule_energy = molecule.get_potential_energy() + pred_ads_energy = self.compute_adsorption_energy( + surface_energy, mol_surf_energy, molecule_energy, adsorbate_count + ) + except Exception as exc: + warn(f"Error calculating energy for {sys_id}: {exc}", stacklevel=2) + pred_ads_energy = np.nan ref_surface_energy = surface.info["ref_energy"] ref_mol_surf_energy = mol_surface.info["ref_energy"] diff --git a/ml_peg/calcs/surfaces/SBH17/calc_SBH17.py b/ml_peg/calcs/surfaces/SBH17/calc_SBH17.py index eb9b44612..00d0315fc 100644 --- a/ml_peg/calcs/surfaces/SBH17/calc_SBH17.py +++ b/ml_peg/calcs/surfaces/SBH17/calc_SBH17.py @@ -5,6 +5,7 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase.io import read, write import numpy as np @@ -61,14 +62,22 @@ def test_surface_barrier(mlip: tuple[str, Any]) -> None: gp.info.setdefault("spin", 1) gp.info["spin"] = int(round(gp.info["spin"])) gp.calc = calc - gp.get_potential_energy() + try: + gp.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {system}: {exc}", stacklevel=2) + gp.info["energy"] = np.nan ts = read(ts_path, index=0) ts.info.setdefault("charge", 0) ts.info.setdefault("spin", 1) ts.info["spin"] = int(round(ts.info["spin"])) ts.calc = copy(calc) - ts.get_potential_energy() + try: + ts.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {system}: {exc}", stacklevel=2) + ts.info["energy"] = np.nan ref = np.loadtxt(ref_path) diff --git a/ml_peg/calcs/surfaces/elemental_slab_oxygen_adsorption/calc_elemental_slab_oxygen_adsorption.py b/ml_peg/calcs/surfaces/elemental_slab_oxygen_adsorption/calc_elemental_slab_oxygen_adsorption.py index c6a379220..a4e610ea2 100644 --- a/ml_peg/calcs/surfaces/elemental_slab_oxygen_adsorption/calc_elemental_slab_oxygen_adsorption.py +++ b/ml_peg/calcs/surfaces/elemental_slab_oxygen_adsorption/calc_elemental_slab_oxygen_adsorption.py @@ -4,12 +4,14 @@ from copy import copy from pathlib import Path +from warnings import warn from ase import Atoms from ase.calculators.calculator import Calculator from ase.io import read, write import mlipx from mlipx.abc import NodeWithCalculator +import numpy as np from tqdm import tqdm import zntrack @@ -105,7 +107,11 @@ def evaluate_energies(atoms_list: list[Atoms], calc: Calculator) -> None: atoms.info.setdefault("charge", 0) atoms.info.setdefault("spin", 1) - atoms.get_potential_energy() + try: + atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + atoms.info["energy"] = np.nan def run(self): """Run oxygen adsorption energy calculations.""" diff --git a/ml_peg/calcs/surfaces/graphene_wetting_under_strain/calc_graphene_wetting_under_strain.py b/ml_peg/calcs/surfaces/graphene_wetting_under_strain/calc_graphene_wetting_under_strain.py index 88eb8439d..84bcb0c55 100644 --- a/ml_peg/calcs/surfaces/graphene_wetting_under_strain/calc_graphene_wetting_under_strain.py +++ b/ml_peg/calcs/surfaces/graphene_wetting_under_strain/calc_graphene_wetting_under_strain.py @@ -4,8 +4,10 @@ from pathlib import Path from typing import Any +from warnings import warn import ase.io +import numpy as np import pytest from tqdm import tqdm import yaml @@ -66,7 +68,11 @@ def test_graphene_wetting_energy(mlip: tuple[str, Any]) -> None: atoms.calc = calc atoms.info.setdefault("charge", 0) atoms.info.setdefault("spin", 1) - water_energy = atoms.get_potential_energy() + try: + water_energy = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + water_energy = np.nan # Iterate through strain conditions for strain in strains: @@ -74,7 +80,14 @@ def test_graphene_wetting_energy(mlip: tuple[str, Any]) -> None: atoms.calc = calc atoms.info.setdefault("charge", 0) atoms.info.setdefault("spin", 1) - graphene_energy = atoms.get_potential_energy() + try: + graphene_energy = atoms.get_potential_energy() + except Exception as exc: + warn( + f"Error calculating energy for strain {strain}: {exc}", + stacklevel=2, + ) + graphene_energy = np.nan # Iterate through orientations for orientation in orientations: @@ -89,7 +102,11 @@ def test_graphene_wetting_energy(mlip: tuple[str, Any]) -> None: atoms.calc = calc atoms.info.setdefault("charge", 0) atoms.info.setdefault("spin", 1) - mlip_potential_energy = atoms.get_potential_energy() + try: + mlip_potential_energy = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy: {exc}", stacklevel=2) + mlip_potential_energy = np.nan mlip_adsorption_energy = ( mlip_potential_energy - graphene_energy - water_energy ) diff --git a/ml_peg/calcs/tm_complexes/3dTMV/calc_3dTMV.py b/ml_peg/calcs/tm_complexes/3dTMV/calc_3dTMV.py index 3182fe8af..37fb178a2 100644 --- a/ml_peg/calcs/tm_complexes/3dTMV/calc_3dTMV.py +++ b/ml_peg/calcs/tm_complexes/3dTMV/calc_3dTMV.py @@ -9,9 +9,11 @@ from copy import copy from pathlib import Path from typing import Any +from warnings import warn from ase import units from ase.io import read, write +import numpy as np import pytest from tqdm import tqdm @@ -140,13 +142,21 @@ def test_3dtmv(mlip: tuple[str, Any]) -> None: atoms_ox.info["charge"] = MOLECULAR_DATA[complex_id]["charge_ox"] atoms_ox.info["spin"] = MOLECULAR_DATA[complex_id]["mult_ox"] atoms_ox.calc = copy(calc) - oxidized_energy = atoms_ox.get_potential_energy() + try: + oxidized_energy = atoms_ox.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {complex_id}: {exc}", stacklevel=2) + oxidized_energy = np.nan atoms_in = atoms.copy() atoms_in.info["charge"] = MOLECULAR_DATA[complex_id]["charge_in"] atoms_in.info["spin"] = MOLECULAR_DATA[complex_id]["mult_in"] atoms_in.calc = copy(calc) - initial_energy = atoms_in.get_potential_energy() + try: + initial_energy = atoms_in.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {complex_id}: {exc}", stacklevel=2) + initial_energy = np.nan model_ion_energy = oxidized_energy - initial_energy diff --git a/ml_peg/calcs/utils/gscdb138.py b/ml_peg/calcs/utils/gscdb138.py index 4d0a019ca..0fc26c1e0 100644 --- a/ml_peg/calcs/utils/gscdb138.py +++ b/ml_peg/calcs/utils/gscdb138.py @@ -4,6 +4,7 @@ from pathlib import Path from typing import Any +from warnings import warn from ase import Atoms, units from ase.io import read, write @@ -122,7 +123,11 @@ def run_gscdb138( print(f"Skipping {specy}") continue atoms.calc = calc - energy = atoms.get_potential_energy() + try: + energy = atoms.get_potential_energy() + except Exception as exc: + warn(f"Error calculating energy for {specy}: {exc}", stacklevel=2) + energy = np.nan e_rel_model += stoi * energy atoms.info["model_energy"] = energy atoms.calc = None