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

@@ -354,7 +354,11 @@ def build_obasis(
         # 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
 
 

iodata/formats/wfx.py

@@ -128,13 +128,16 @@ def load_data_wfx(lit: LineIterator) -> dict:
     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=" ")