Merge pull request #292 from SCM-NV/orb

BUG: Remove MO padding when the requested MO range is larger than the number of available MOs

src/qmflows/parsers/cp2KParser.py

@@ -25,7 +25,7 @@
 from ..type_hints import Literal as Literal_
 from ..type_hints import PathLike, T, WarnDict, WarnMap
 from ..utils import file_to_context
-from ..warnings_qmflows import QMFlows_Warning
+from ..warnings_qmflows import QMFlows_Warning, Orbital_Warning
 from .parser import minusOrplus, natural, point, try_search_pattern
 from .xyzParser import manyXYZ, tuplesXYZ_to_plams
 from ._cp2k_basis_parser import readCp2KBasis
@@ -213,23 +213,37 @@ def read_coefficients(
     energies = np.empty(norbitals)
     coefficients = np.empty((norbital_functions, norbitals))
 
+    n_orb_actual = 0
     for i, lines in enumerate(chunks):
         j = 2 * i
         es = lines[1]
         # There could be a single or double column
         css = [k[4:] for k in lines[3:]]
         # There is an odd number of MO and this is the last one
         if len(es) == 1:
-            energies[-1] = float(es[0])
-            coefficients[:, -1] = np.concatenate(css)
+            n_orb_actual = j + 1
+            energies[j] = float(es[0])
+            coefficients[:, j] = np.concatenate(css)
         else:
             # rearrange the coefficients
+            n_orb_actual = j + 2
             css2 = np.transpose(css)
             energies[j: j + 2] = es
             coefficients[:, j] = css2[0]
             coefficients[:, j + 1] = css2[1]
 
-    return InfoMO(energies, coefficients)
+    # It's possible that the range of requested MOs is larger than the
+    # number of available MOs for a given basis set; the arrays will require
+    # some trimming in that case
+    if norbitals != n_orb_actual:
+        warnings.warn(
+            f"Trimming MOs, the requested number of MOs ({norbitals}) is larger "
+            f"than the available amount ({n_orb_actual})",
+            Orbital_Warning,
+        )
+        return InfoMO(energies[:n_orb_actual], coefficients[:, :n_orb_actual])
+    else:
+        return InfoMO(energies, coefficients)
 
 
 def _get_mos(path: PathLike) -> List[List[str]]:

src/qmflows/warnings_qmflows.py

@@ -11,7 +11,7 @@
 __all__ = [
     'QMFlows_Warning', 'Key_Warning', 'Assertion_Warning',
     'SCF_Convergence_Warning', 'Geometry_Convergence_Warning',
-    'Parameter_Warning', 'Charge_Warning',
+    'Parameter_Warning', 'Charge_Warning', 'Orbital_Warning',
     'cp2k_warnings'
 ]
 
@@ -44,6 +44,10 @@ class Charge_Warning(Parameter_Warning):
     """Warning class for charges in classical forcefields."""
 
 
+class Orbital_Warning(QMFlows_Warning):
+    """Warning class for orbital-related issues."""
+
+
 def _eval_charge(msg: str, tolerance: float = 0.1) -> Optional[str]:
     """Check of the total molecular charge is integer within a given *tolerance*."""
     charge = float(msg.rsplit(maxsplit=1)[1])

test/test_using_plams/test_cp2k.py

@@ -1,14 +1,16 @@
 """Run an small cp2k calculation."""
 
+import warnings
 from pathlib import Path
 
+import pytest
 import numpy as np
 import h5py
-import pytest
 from assertionlib import assertion
 from scm.plams import Molecule
 
 from qmflows import Settings, cp2k, run, templates
+from qmflows.warnings_qmflows import Orbital_Warning
 from qmflows.test_utils import (
     PATH,
     PATH_MOLECULES,
@@ -40,7 +42,7 @@ def test_cp2k_opt(tmp_path: Path) -> None:
 
 @requires_cp2k
 @pytest.mark.slow
-@pytest.mark.parametrize("mo_index_range", ["1_4", "1_2", "3_4"])
+@pytest.mark.parametrize("mo_index_range", ["1_4", "1_2", "3_4", "0_99"])
 def test_cp2k_singlepoint(tmp_path: Path, mo_index_range: str) -> None:
     """Run a simple single point."""
     mol = Molecule(PATH_MOLECULES / "h2o.xyz", 'xyz', charge=0, multiplicity=1)
@@ -58,13 +60,20 @@ def test_cp2k_singlepoint(tmp_path: Path, mo_index_range: str) -> None:
     )
     job = cp2k(s, mol)
     result = run(job, path=tmp_path, folder="test_cp2k_singlepoint")
+    validate_status(result)
 
     with h5py.File(PATH / "test_output.hdf5", "r") as f:
         key = f"test_using_plams/test_cp2k/test_cp2k_singlepoint/{mo_index_range}"
         ref = f[key][...].view(np.recarray)
 
-    validate_status(result)
-    orbitals = result.orbitals
+    if mo_index_range == "0_99":
+        with pytest.warns(Orbital_Warning):
+            orbitals = result.orbitals
+    else:
+        with warnings.catch_warnings():
+            warnings.simplefilter("error", Orbital_Warning)
+            orbitals = result.orbitals
+
     assertion.is_not(orbitals, None)
     np.testing.assert_allclose(orbitals.eigenvalues, ref.eigenvalues)
     np.testing.assert_allclose(np.abs(orbitals.eigenvectors).T, ref.eigenvectors)