Replace asserts by errors or warnings

iodata/formats/cp2klog.py

@@ -252,7 +252,8 @@ def _read_cp2k_orbital_coeffs(
     next(lit)
     while len(allcoeffs) < len(oe):
         line = next(lit)
-        assert line.startswith("    ORBITAL      L =")
+        if not line.startswith("    ORBITAL      L ="):
+            raise LoadError("Did not find the expected start of a new ORBITAL section.", lit)
         words = line.split()
         angmom = int(words[3])
         state = int(words[6])
@@ -407,7 +408,8 @@ def load_one(lit: LineIterator) -> dict:
     for line in lit:
         if line.startswith("          Core Charge"):
             atcorenum = float(line[70:])
-            assert atcorenum == int(atcorenum)
+            if atcorenum != int(atcorenum):
+                raise LoadError("Inconsistent effective core charge", lit)
             break
         if line.startswith(" Electronic structure"):
             atcorenum = float(atnum)
@@ -447,7 +449,8 @@ def load_one(lit: LineIterator) -> dict:
     # Turn orbital data into a MolecularOrbitals object.
     if restricted:
         norb, nel = _get_norb_nel(oe_alpha)
-        assert nel % 2 == 0
+        if nel % 2 != 0:
+            raise LoadError("Odd number of electrons for a restricted wavefunction.", lit)
         orb_alpha_coeffs = np.zeros([obasis.nbasis, norb])
         orb_alpha_energies = np.zeros(norb)
         orb_alpha_occs = np.zeros(norb)
@@ -466,7 +469,8 @@ def load_one(lit: LineIterator) -> dict:
     else:
         norb_alpha = _get_norb_nel(oe_alpha)[0]
         norb_beta = _get_norb_nel(oe_beta)[0]
-        assert norb_alpha == norb_beta
+        if norb_alpha != norb_beta:
+            raise LoadError("Inconsistent number of alpha and beta orbitals.", lit)
         orb_alpha_coeffs = np.zeros([obasis.nbasis, norb_alpha])
         orb_alpha_energies = np.zeros(norb_alpha)
         orb_alpha_occs = np.zeros(norb_alpha)

iodata/formats/extxyz.py

@@ -34,7 +34,7 @@
 
 from ..docstrings import document_load_many, document_load_one
 from ..periodic import num2sym, sym2num
-from ..utils import LineIterator, amu, angstrom, strtobool
+from ..utils import LineIterator, LoadError, amu, angstrom, strtobool
 from .xyz import load_one as load_one_xyz
 
 __all__ = ()
@@ -76,8 +76,22 @@ def _convert_title_value(value: str):
     return converted_value
 
 
-def _parse_properties(properties: str):
-    """Parse the properties into atom_columns."""
+def _parse_properties(properties: str, lit: LineIterator) -> list[tuple]:
+    """Parse the properties listed in the title into atom_columns.
+
+    Parameters
+    ----------
+    properties
+        The part of the title of an XYZ file containing column definitions.
+    lit
+        The LineIterator reading the file, only used for error messages.
+
+    Returns
+    -------
+    atom_columns
+        A list of tuples specifying the columns in the extended XYZ file.
+        For more details, see ``ATOM_COLUMNS_DOC`` in ``iodata.formats.xyz``.
+    """
     atom_columns = []
     # Maps the dtype to the atom_columns dtype, load_word and dump_word
     dtype_map = {
@@ -124,7 +138,8 @@ def _parse_properties(properties: str):
         (lambda atnum: f"{num2sym[atnum]:2s}"),
     )
     splitted_properties = properties.split(":")
-    assert len(splitted_properties) % 3 == 0
+    if len(splitted_properties) % 3 != 0:
+        raise LoadError(f"Cannot parse property from the title line: {properties}", lit)
     # Each property has 3 values: its name, dtype and shape
     names = splitted_properties[::3]
     dtypes = splitted_properties[1::3]
@@ -145,8 +160,25 @@ def _parse_properties(properties: str):
     return atom_columns
 
 
-def _parse_title(title: str):
-    """Parse the title in an extended xyz file."""
+def _parse_title(title: str, lit: LineIterator) -> tuple[list[tuple], dict[str]]:
+    """Parse the title in an extended xyz file.
+
+    Parameters
+    ----------
+    title
+        The full title line.
+    lit
+        The LineIterator reading the file, only used for error messages.
+
+    Returns
+    -------
+    atom_columns
+        A list of tuples specifying the columns in the extended XYZ file.
+        For more details, see ``ATOM_COLUMNS_DOC`` in ``iodata.formats.xyz``.
+    data
+        Attributes to be included in the ``IOData`` instance,
+        taken from the title line.
+    """
     key_value_pairs = shlex.split(title)
     # A dict of predefined iodata atrributes with their names and dtype convertion functions
 
@@ -163,7 +195,7 @@ def load_cellvecs(word):
         if "=" in key_value_pair:
             key, value = key_value_pair.split("=", 1)
             if key == "Properties":
-                atom_columns = _parse_properties(value)
+                atom_columns = _parse_properties(value, lit)
             elif key in iodata_attrs:
                 data[iodata_attrs[key][0]] = iodata_attrs[key][1](value)
             else:
@@ -184,7 +216,7 @@ def load_one(lit: LineIterator) -> dict:
     atom_line = next(lit)
     title_line = next(lit)
     # parse title
-    atom_columns, title_data = _parse_title(title_line)
+    atom_columns, title_data = _parse_title(title_line, lit)
     lit.back(title_line)
     lit.back(atom_line)
     xyz_data = load_one_xyz(lit, atom_columns)

iodata/formats/fchk.py

@@ -21,6 +21,7 @@
 from collections.abc import Iterator
 from fnmatch import fnmatch
 from typing import Optional, TextIO
+from warnings import warn
 
 import numpy as np
 from numpy.typing import NDArray
@@ -29,7 +30,7 @@
 from ..docstrings import document_dump_one, document_load_many, document_load_one
 from ..iodata import IOData
 from ..orbitals import MolecularOrbitals
-from ..utils import DumpError, LineIterator, LoadError, PrepareDumpError, amu
+from ..utils import DumpError, LineIterator, LoadError, LoadWarning, PrepareDumpError, amu
 
 __all__ = ()
 
@@ -336,7 +337,15 @@ def load_many(lit: LineIterator) -> Iterator[dict]:
                 fchk[f"{prefix} {ipoint + 1:7d} Gradient at each geome"].reshape(-1, natom, 3),
             )
         )
-        assert len(trajectory) == nstep
+        if len(trajectory) != nstep:
+            warn(
+                LoadWarning(
+                    "The size of the optimization trajectory is inconsistent with the values in the"
+                    "field 'Optimization Number of geometries'. The latter is ignored.",
+                    lit,
+                ),
+                stacklevel=2,
+            )
         for istep, (energy, recor, atcoords, gradients) in enumerate(trajectory):
             data = {
                 "title": fchk["title"],
@@ -349,7 +358,7 @@ def load_many(lit: LineIterator) -> Iterator[dict]:
                     "ipoint": ipoint,
                     "npoint": len(nsteps),
                     "istep": istep,
-                    "nstep": nstep,
+                    "nstep": len(trajectory),
                 },
             }
             if prefix == "IRC point":

iodata/formats/fcidump.py

@@ -29,12 +29,13 @@
 """
 
 from typing import TextIO
+from warnings import warn
 
 import numpy as np
 
 from ..docstrings import document_dump_one, document_load_one
 from ..iodata import IOData
-from ..utils import LineIterator, LoadError, set_four_index_element
+from ..utils import LineIterator, LoadError, LoadWarning, set_four_index_element
 
 __all__ = ()
 
@@ -86,9 +87,11 @@ def load_one(lit: LineIterator) -> dict:
             ij = int(words[2]) - 1
             ik = int(words[3]) - 1
             il = int(words[4]) - 1
-            # Uncomment the following line if you want to assert that the
-            # FCIDUMP file does not contain duplicate 4-index entries.
-            # assert two_mo.get_element(ii,ik,ij,il) == 0.0
+            if two_mo[ii, ik, ij, il] != 0.0:
+                warn(
+                    LoadWarning("Duplicate entries in the FCIDUMP file are ignored", lit),
+                    stacklevel=2,
+                )
             set_four_index_element(two_mo, ii, ik, ij, il, value)
         elif words[1] != "0":
             ii = int(words[1]) - 1

iodata/formats/molden.py

@@ -579,8 +579,6 @@ def _fix_mo_coeffs_psi4(obasis: MolecularBasis) -> Union[MolecularBasis, None]:
     correction = []
     corrected = False
     for shell in obasis.shells:
-        # We can safely assume segmented shells.
-        assert shell.ncon == 1
         angmom = shell.angmoms[0]
         kind = shell.kinds[0]
         factors = None
@@ -594,7 +592,8 @@ def _fix_mo_coeffs_psi4(obasis: MolecularBasis) -> Union[MolecularBasis, None]:
         if factors is None:
             factors = np.ones(shell.nbasis)
         else:
-            assert len(factors) == shell.nbasis
+            if len(factors) != shell.nbasis:
+                raise AssertionError("Internal error when correcting for Psi4 errors.")
             corrected = True
         correction.append(factors)
     if corrected:
@@ -612,8 +611,6 @@ def _fix_mo_coeffs_cfour(obasis: MolecularBasis) -> Union[MolecularBasis, None]:
     correction = []
     corrected = False
     for shell in obasis.shells:
-        # We can safely assume segmented shells.
-        assert shell.ncon == 1
         angmom = shell.angmoms[0]
         kind = shell.kinds[0]
         factors = None
@@ -632,7 +629,8 @@ def _fix_mo_coeffs_cfour(obasis: MolecularBasis) -> Union[MolecularBasis, None]:
         if factors is None:
             factors = np.ones(shell.nbasis)
         else:
-            assert len(factors) == shell.nbasis
+            if len(factors) != shell.nbasis:
+                raise AssertionError("Internal error when correcting for CFOUR errors.")
             corrected = True
         correction.append(factors)
     if corrected:
@@ -670,6 +668,8 @@ def _fix_molden_from_buggy_codes(result: dict, lit: LineIterator, norm_threshold
         coeffsb = result["mo"].coeffsb
     else:
         raise LoadError(f"Molecular orbital kind={result['mo'].kind} not recognized.", lit)
+    if any(shell.ncon != 1 for shell in obasis.shells):
+        raise LoadError("Generalized contractions are not supported", lit)
 
     if _is_normalized_properly(obasis, atcoords, coeffsa, coeffsb, norm_threshold):
         # The file is good. No need to change obasis.

iodata/formats/molekel.py

@@ -125,21 +125,26 @@ def _load_helper_coeffs(
             # read a1g line
             words = line.split()
             ncol = len(words)
-            assert ncol > 0
+            if ncol == 0:
+                raise LoadError("Expect irrep, got empty line", line)
             irreps.extend(words)
             cols = [np.zeros((nbasis, 1), float) for _ in range(ncol)]
             in_orb = 1
         elif in_orb == 1:
             # read energies
             words = line.split()
-            assert len(words) == ncol
+            if len(words) != ncol:
+                raise LoadError(f"Wrong number of energies: expected {ncol}, got {len(words)}", lit)
             energies.extend(float(word) for word in words)
             in_orb = 2
             ibasis = 0
         elif in_orb == 2:
             # read expansion coefficients
             words = line.split()
-            assert len(words) == ncol
+            if len(words) != ncol:
+                raise LoadError(
+                    f"Wrong number of coefficients: expected {ncol}, got {len(words)}", lit
+                )
             for icol in range(ncol):
                 cols[icol][ibasis] = float(words[icol])
             ibasis += 1
@@ -224,8 +229,10 @@ def load_one(lit: LineIterator, norm_threshold: float = 1e-4) -> dict:
     nelec = atnums.sum() - charge
     if coeffsb is None:
         # restricted closed-shell
-        assert nelec % 2 == 0
-        assert abs(occsa.sum() - nelec) < 1e-7
+        if nelec % 2 != 0:
+            raise LoadError("Odd number of electrons found in restricted case.", lit)
+        if abs(occsa.sum() - nelec) > 1e-7:
+            raise LoadError("Occupation numbers are inconsistent with number of electrons", lit)
         mo = MolecularOrbitals(
             "restricted", coeffsa.shape[1], coeffsa.shape[1], occsa, coeffsa, energiesa, irrepsa
         )
@@ -235,22 +242,20 @@ def load_one(lit: LineIterator, norm_threshold: float = 1e-4) -> dict:
                 "Beta occupation numbers not found in mkl file while beta orbitals were present.",
                 lit,
             )
-        nalpha = int(np.round(occsa.sum()))
-        nbeta = int(np.round(occsb.sum()))
+        nalpha = occsa.sum()
+        nbeta = occsb.sum()
         if abs(spinpol - abs(nalpha - nbeta)) > 1e-7:
             warn(
                 LoadWarning(
-                    f"The spin polarization ({spinpol}) is inconsistent with the"
+                    f"The spin polarization ({spinpol}) is inconsistent with the "
                     f"difference between alpha and beta occupation numbers ({nalpha} - {nbeta}). "
                     "The spin polarization will be rederived from the occupation numbers.",
                     lit,
                 ),
                 stacklevel=2,
             )
-        assert nelec == nalpha + nbeta
-        assert coeffsa.shape == coeffsb.shape
-        assert energiesa.shape == energiesb.shape
-        assert occsa.shape == occsb.shape
+        if abs(nelec - (nalpha + nbeta)) > 1e-7:
+            raise LoadError("Occupation numbers are inconsistent with number of electrons", lit)
         mo = MolecularOrbitals(
             "unrestricted",
             coeffsa.shape[1],

iodata/formats/mwfn.py

@@ -208,7 +208,8 @@ def _load_helper_mo(lit: LineIterator, n_basis: int, n_mo: int) -> dict:
         line = next(lit)
         while "Index" not in line:
             line = next(lit)
-        assert line.startswith("Index")
+        if not line.startswith("Index"):
+            raise LoadError(f"Expecting line starting with 'Index', got '{line}'", lit)
         data["mo_numbers"][index] = line.split()[1]
         data["mo_type"][index] = next(lit).split()[1]
         data["mo_energies"][index] = next(lit).split()[1]
@@ -324,7 +325,8 @@ def load_one(lit: LineIterator) -> dict:
         # 2: beta
         norba = (inp["mo_type"] == 1).sum()
         norbb = (inp["mo_type"] == 2).sum()
-        assert (inp["mo_type"] == 0).sum() == 0
+        if (inp["mo_type"] == 0).sum() != 0:
+            raise LoadError("Restricted orbtials found in unrestricted wavefunction.", lit)
     # Build MolecularOrbitals instance
     mo = MolecularOrbitals(
         inp["mo_kind"], norba, norbb, inp["mo_occs"], inp["mo_coeffs"], inp["mo_energies"], None

iodata/formats/qchemlog.py

@@ -27,7 +27,7 @@
 from ..docstrings import document_load_one
 from ..orbitals import MolecularOrbitals
 from ..periodic import sym2num
-from ..utils import LineIterator, angstrom, calmol, kcalmol, kjmol, strtobool
+from ..utils import LineIterator, LoadError, angstrom, calmol, kcalmol, kjmol, strtobool
 
 __all__ = ()
 
@@ -109,8 +109,12 @@ def load_one(lit: LineIterator) -> dict:
         # Convert to alphabetical ordering: xx, xy, xz, yy, yz, zz
         moments[(2, "c")] = data["quadrupole"][[0, 1, 3, 2, 4, 5]]
     # check total dipole parsed
-    if "dipole_tol" in data and "dipole" in data:
-        assert abs(np.linalg.norm(data["dipole"]) - data["dipole_tol"]) < 1.0e-4
+    if (
+        "dipole_tol" in data
+        and "dipole" in data
+        and abs(np.linalg.norm(data["dipole"]) - data["dipole_tol"]) > 1.0e-4
+    ):
+        raise LoadError("Inconsistent dipole and dipole_tol", lit.filename)
     if moments:
         result["moments"] = moments
 

iodata/formats/wfn.py

@@ -259,103 +259,107 @@
     )
 
 
 def build_obasis(
     icenters: NDArray[int],
     type_assignments: NDArray[int],
     exponents: NDArray[float],
     lit: LineIterator,
 ) -> tuple[MolecularBasis, NDArray[int]]:
     """Construct a basis set using the arrays read from a WFN or WFX file.
 
     Parameters
     ----------
     icenters
         The center indices for all basis functions. shape=(nbasis,). Lowest
         index is zero.
     type_assignments
         Integer codes for basis function names. shape=(nbasis,). Lowest index
         is zero.
     exponents
         The Gaussian exponents of all basis functions. shape=(nbasis,)
 
     """
     # Build the basis set, keeping track of permutations in case there are
     # deviations from the default ordering of primitives in a WFN file.
     shells = []
     ibasis = 0
     nbasis = len(icenters)
     permutation = np.zeros(nbasis, dtype=int)
     # Loop over all (batches of primitive) basis functions and extract shells.
     while ibasis < nbasis:
         # Determine the angular moment of the shell
         type_assignment = type_assignments[ibasis]
         # multiple different type assignments (codes for individual basis
         # functions) can match one angular momentum.
         angmom = 0 if type_assignment == 0 else len(PRIMITIVE_NAMES[type_assignments[ibasis]])
         # The number of cartesian functions for the current angular momentum
         ncart = len(CONVENTIONS[(angmom, "c")])
         # Determine how many shells are to be read in one batch. E.g. for a
         # contracted p shell, the WFN format contains first all px basis
         # functions, the all py, finally all pz. These need to be regrouped into
         # shells.
         # This pattern can almost be used to reverse-engineer contractions.
         # One should also check (i) if the corresponding mo-coefficients are the
         # same (after fixing them for normalization) and (ii) if the functions
         # are centered on the same atom.
         # For now, this implementation makes no attempt to reverse-engineer
         # contractions, but it can be done.
         ncon = 1  # the contraction length
         if angmom > 0:
             # batches for s-type functions are not necessary and may result in
             # multiple centers being pulled into one batch.
             while (
                 ibasis + ncon < len(type_assignments)
                 and type_assignments[ibasis + ncon] == type_assignment
             ):
                 ncon += 1
         # Check if the type assignment is consistent for remaining basis
         # functions in this batch.
         for ifn in range(ncart):
             if not (
                 type_assignments[ibasis + ncon * ifn : ibasis + ncon * (ifn + 1)]
                 == type_assignments[ibasis + ncon * ifn]
             ).all():
                 raise LoadError("Inconcsistent type assignments in current batch of shells.", lit)
         # Check if all basis functions in the current batch sit on
         # the same center. If not, IOData cannot read this file.
         icenter = icenters[ibasis]
         if not (icenters[ibasis : ibasis + ncon * ncart] == icenter).all():
             raise LoadError("Incomplete shells in WFN file not supported by IOData.", lit)
         # Check if the same exponent is used for corresponding basis functions.
         batch_exponents = exponents[ibasis : ibasis + ncon]
         for ifn in range(ncart):
             if not (
                 exponents[ibasis + ncon * ifn : ibasis + ncon * (ifn + 1)] == batch_exponents
             ).all():
                 raise LoadError(
                     "Exponents must be the same for corresponding basis functions.", lit
                 )
         # A permutation is needed because we need to regroup basis functions
         # into shells.
         batch_primitive_names = [
             PRIMITIVE_NAMES[type_assignments[ibasis + ifn * ncon]] for ifn in range(ncart)
         ]
         for irep in range(ncon):
             for i, primitive_name in enumerate(batch_primitive_names):
                 ifn = CONVENTIONS[(angmom, "c")].index(primitive_name)
                 permutation[ibasis + irep * ncart + ifn] = ibasis + irep + i * ncon
         # WFN uses non-normalized primitives, which will be corrected for
         # when processing the MO coefficients. Normalized primitives will
         # be used here. No attempt is made here to reconstruct the contraction.
         shells.extend(
             Shell(icenter, [angmom], ["c"], np.array([exponent]), np.array([[1.0]]))
             for exponent in batch_exponents
         )
         # Move on to the next contraction
         ibasis += ncart * ncon
     obasis = MolecularBasis(shells, CONVENTIONS, "L2")
-    assert obasis.nbasis == nbasis
+    if obasis.nbasis != nbasis:
+        raise LoadError(
+            "The basis set size derived from icenters is inconsistent with "
+            "the number of primitives."
+        )
     return obasis, permutation
 
 
 def get_mocoeff_scales(obasis: MolecularBasis) -> NDArray[float]:

iodata/formats/wfx.py

@@ -117,55 +117,58 @@
     )
 
 
 def load_data_wfx(lit: LineIterator) -> dict:
     """Process loaded WFX data."""
     # get all section labels and required labels for WFX files
     lbs_str, lbs_int, lbs_float, lbs_aint, lbs_afloat, lbs_other, required_tags = _wfx_labels()
     # load sections in WFX and check required tags exists
     data = parse_wfx(lit, required_tags)
 
     # process raw data to convert them to proper type based on their label
     result = {}
     for key, value in data.items():
         if key in lbs_str:
-            assert len(value) == 1
+            if len(value) != 1:
+                raise LoadError(f"Expecting one value, got {len(value)}", lit)
             result[lbs_str[key]] = value[0]
         elif key in lbs_int:
-            assert len(value) == 1
+            if len(value) != 1:
+                raise LoadError(f"Expecting one value, got {len(value)}", lit)
             result[lbs_int[key]] = int(value[0])
         elif key in lbs_float:
-            assert len(value) == 1
+            if len(value) != 1:
+                raise LoadError(f"Expecting one value, got {len(value)}", lit)
             result[lbs_float[key]] = float(value[0])
         elif key in lbs_afloat:
             result[lbs_afloat[key]] = np.fromstring(" ".join(value), dtype=float, sep=" ")
         elif key in lbs_aint:
             result[lbs_aint[key]] = np.fromstring(" ".join(value), dtype=int, sep=" ")
         elif key in lbs_other:
             result[lbs_other[key]] = value
         else:
             warn(LoadWarning(f"Not recognized section label, skip {key}", lit), stacklevel=2)
 
     # reshape some arrays
     result["atcoords"] = result["atcoords"].reshape(-1, 3)
     result["mo_coeffs"] = result["mo_coeffs"].reshape(result["num_primitives"], -1, order="F")
     # process nuclear gradient, if present
     if "nuclear_gradient" in result:
         gradient_mix = np.array([i.split() for i in result.pop("nuclear_gradient")]).reshape(-1, 4)
         gradient_atoms = gradient_mix[:, 0].astype(np.str_)
         index = [result["nuclear_names"].index(atom) for atom in gradient_atoms]
         result["atgradient"] = np.full((len(result["nuclear_names"]), 3), np.nan)
         result["atgradient"][index] = gradient_mix[:, 1:].astype(float)
     # check keywords & number of perturbations
     perturbation_check = {"GTO": 0, "GIAO": 3, "CGST": 6}
     key = result["keywords"]
     num = result["num_perturbations"]
     if key not in perturbation_check:
         raise LoadError(f"The keywords is {key}, but it should be either GTO, GIAO or CGST.", lit)
     if num != perturbation_check[key]:
         raise LoadError(
             f"Number of perturbations of {key} is {num}, expected {perturbation_check[key]}.", lit
         )
     return result
 
 
 def parse_wfx(lit: LineIterator, required_tags: Optional[list] = None) -> dict: