diff --git a/CHANGELOG.md b/CHANGELOG.md index 3ef7c55c..3b7431e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ - Added `calc_rotational_constants` to calculate molecular rotational constants to `Structure` class (#230). - Added `calc_rotor_type` to classify a molecule's rotor type to `Structure` class (#230). - Add type checks for ncores and memory setters in `Input` and checks for invalid negative values (#261). +- Added `Output.get_fcidump()` and the `Fcidump` class for parsing FCIDUMP files generated by ORCA (#253). ### Changed - Refactored methods from Runner into BaseRunner (#193) diff --git a/examples/exmp054_casscf_fcidump/job.py b/examples/exmp054_casscf_fcidump/job.py index fcbad4e5..d4b18d90 100755 --- a/examples/exmp054_casscf_fcidump/job.py +++ b/examples/exmp054_casscf_fcidump/job.py @@ -23,8 +23,8 @@ def run_exmp054( calc = Calculator(basename="job", working_dir=working_dir) calc.structure = structure - - calc.input.add_blocks(BlockCasscf(nel=2, norb=2), BlockOutput(dumpactints=True)) + casscf_block = BlockCasscf(nel=2, norb=2) + calc.input.add_blocks(casscf_block, BlockOutput(dumpactints=True)) calc.input.ncores = 4 calc.write_input() @@ -48,13 +48,18 @@ def run_exmp054( print("CASSCF energy") print(output.results_properties.geometries[0].energy[0].totalenergy[0][0]) - # > retrieve the path to the fcidump file - fcidump_file = output.get_outfile().with_suffix(".fcidump") + # > parse the fcidump file + fcidump = output.get_fcidump() + + if fcidump is None: + print(f"FCIDUMP generation failed, see output file: {output.get_outfile()}") + print(output.error_message()) + sys.exit(1) - # > Read and print the file - print("FCIDUMP file:") - with open(fcidump_file, "r") as f: - print(f.read()) + print("hcore matrix") + print(fcidump.hcore_matrix) + print("ERI tensor") + print(fcidump.eri_tensor) return output diff --git a/src/opi/output/core.py b/src/opi/output/core.py index faab0b02..4daae395 100644 --- a/src/opi/output/core.py +++ b/src/opi/output/core.py @@ -16,6 +16,7 @@ from opi.execution.core import Runner from opi.input.structures import Atom, Coordinates, Structure from opi.output.cube import CubeOutput +from opi.output.fcidump import Fcidump from opi.output.gbw_suffix import GbwSuffix from opi.output.grepper.recipes import ( get_error_message, @@ -2742,3 +2743,33 @@ def get_free_solvation_energy_opencosmors(self) -> float | None: except (KeyError, IndexError): return None return float(free_solv_energy) + + def get_fcidump(self) -> Fcidump | None: + """ + Parse the FCIDUMP file generated by ORCA and return its data in the `Fcidump` data class. + The FCIDUMP file has to be generated by the ORCA job and cannot be generated on-the-fly after the calculation. + To generate FCIDUMP files, set the following options: + ``` + %output + dumpactints true + end + ``` + See `BlockOutput` for the related block option. + + Returns + ------- + fcidump_data: Fcidump | None + The parsed FCIDUMP data or None if the file is not present or could not be parsed. + """ + + # > Get path to the FCIDUMP file + fci_path = self.get_file(".fcidump") + + # > If there is no file we return None + if not fci_path.is_file(): + return None + + try: + return Fcidump.from_file(fci_path) + except ValueError: + return None diff --git a/src/opi/output/fcidump.py b/src/opi/output/fcidump.py new file mode 100644 index 00000000..40b7b794 --- /dev/null +++ b/src/opi/output/fcidump.py @@ -0,0 +1,184 @@ +"""Parse a potential FCIDUMP file""" + +import re +from dataclasses import dataclass, field +from functools import cached_property +from pathlib import Path +from typing import Self + +import numpy as np + + +@dataclass +class Fcidump: + """ + Reads and stores data from a FCIDUMP file. One and two-electrons integrals are stored as dicts and can be + accessed as numpy arrays via `hcore_matrix` and `eri_tensor`. + + Attributes + -------- + norb: int + Number of active orbitals. + nelec: int + Number of active electrons. + ms2: int + Twice the total spin projection, i.e. the difference of alpha and beta electrons. + orbsym: list[int] + Symmetry labels of the orbitals. + isym: int + Overall symmetry of the electronic structure. + one_electron: dict[tuple[int, int], float] + Dictionary that contains the one-electron integrals. + two_electron: dict[tuple[int, int, int, int], float] + Dictionary that contains the two-electron integrals. + orbital_energies: dict[int, float] + Dictionary that contains the orbital energies, if present in the FCIDUMP file. + e_nuc: float + Core energy contribution. Contains the contracted energy of the inactive space. + path: Path + Path to the FCIDUMP file. + """ + + norb: int + nelec: int + ms2: int + orbsym: list[int] + isym: int + one_electron: dict[tuple[int, int], float] = field(default_factory=dict) + two_electron: dict[tuple[int, int, int, int], float] = field(default_factory=dict) + orbital_energies: dict[int, float] = field(default_factory=dict) + e_nuc: float = 0.0 + path: Path = field(default_factory=Path) + + @cached_property + def hcore_matrix(self) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]: + """Return the one-electron integrals as a symmetric (norb, norb) numpy array.""" + mat = np.zeros((self.norb, self.norb)) + for (i, j), val in self.one_electron.items(): + mat[i - 1, j - 1] = val + mat[j - 1, i - 1] = val + return mat + + @cached_property + def eri_tensor(self) -> np.ndarray[tuple[int, int, int, int], np.dtype[np.float64]]: + """Return the two-electron integrals as a (norb, norb, norb, norb) numpy array. + + Uses chemist's notation (ij|kl) with 8-fold permutation symmetry applied. + """ + tensor = np.zeros((self.norb,) * 4) + if not self.two_electron: + return tensor + + # > Pull all stored indices/values into arrays once + idx = np.array(list(self.two_electron.keys()), dtype=np.int64) - 1 # (n_integrals, 4) + vals = np.array(list(self.two_electron.values()), dtype=np.float64) # (n_integrals,) + a, b, c, d = idx[:, 0], idx[:, 1], idx[:, 2], idx[:, 3] + + # > Vectorized assignment for each of the 8 symmetry-equivalent index permutations + for p, q, r, s in ( + (a, b, c, d), + (b, a, c, d), + (a, b, d, c), + (b, a, d, c), + (c, d, a, b), + (d, c, a, b), + (c, d, b, a), + (d, c, b, a), + ): + tensor[p, q, r, s] = vals + return tensor + + @classmethod + def from_file(cls, path: Path | str) -> Self: + """ + Parse a FCIDUMP file and return the populated `Fcidump` object. + The FCIDUMP file is documented in the paper: + Knowles, P. J.; Handy, N. C. A Determinant Based Full Configuration Interaction Program. + Computer Physics Communications 1989, 54, 75–83. https://doi.org/10.1016/0010-4655(89)90033-7 + + Raises + ------- + ValueError + If the FCIDUMP file cannot be parsed. + FileNotFoundError + If the FCIDUMP file cannot be found at the given path. + """ + path = Path(path) + + if not path.is_file(): + raise FileNotFoundError(f"{cls.__name__}: FCIDUMP file not found at {path}") + + text = path.read_text() + + # > Split header and body + end_match = re.search(r"&END|/", text, re.IGNORECASE) + if end_match is None: + raise ValueError( + f"{cls.__name__}: Could not find header terminator (&END or /) in {path}" + ) + header = text[: end_match.end()] + body = text[end_match.end() :] + + # > Parse the header + dump = cls( + norb=cls._get_int("NORB", header), + nelec=cls._get_int("NELEC", header), + ms2=cls._get_int("MS2", header), + orbsym=cls._get_int_list("ORBSYM", header), + isym=cls._get_int("ISYM", header), + path=Path(path), + ) + + # > Parse the integrals from the body + for line in body.splitlines(): + parts = line.split() + if not parts: + continue + if len(parts) != 5: + raise ValueError(f"{cls.__name__}: Could not parse {line} in {path}") + try: + val, i, j, k, ll = ( + float(parts[0]), + int(parts[1]), + int(parts[2]), + int(parts[3]), + int(parts[4]), + ) + except ValueError: + raise ValueError(f"{cls.__name__}: Could not parse {line} in {path}") + # > Inactive contribution + if i == 0 and j == 0 and k == 0 and ll == 0: + dump.e_nuc = val + # > Orbital energies (not written by ORCA, but by some other programs) + elif j == 0 and k == 0 and ll == 0: + dump.orbital_energies[i] = val + # > One-electron matrix + elif k == 0 and ll == 0: + dump.one_electron[(i, j)] = val + # > Two-electron tensor + else: + dump.two_electron[(i, j, k, ll)] = val + + return dump + + @classmethod + def _get_int(cls, key: str, header: str) -> int: + """Return the positive integer value of the given key.""" + m = re.search(rf"{key}\s*=\s*(\d+)", header, re.IGNORECASE) + if m is None: + raise ValueError(f"{cls.__name__}: Could not parse {key}") + return int(m.group(1)) + + @classmethod + def _get_int_list(cls, key: str, header: str) -> list[int]: + """Return a list of non-negative integers corresponding to the given key.""" + m = re.search(rf"{key}\s*=\s*([-\d\s,]+)", header, re.IGNORECASE) + if m is None: + raise ValueError(f"{cls.__name__}: Could not parse {key}") + try: + values = [int(x) for x in re.split(r"[,\s]+", m.group(1).strip()) if x] + except ValueError: + raise ValueError(f"{cls.__name__}: Could not parse {key}") + if not values or any(v < 0 for v in values): + raise ValueError(f"{cls.__name__}: Could not parse {key}") + return values diff --git a/tests/examples/test_exmp054_casscf_fcidump.py b/tests/examples/test_exmp054_casscf_fcidump.py index 40a081fb..584eaa14 100644 --- a/tests/examples/test_exmp054_casscf_fcidump.py +++ b/tests/examples/test_exmp054_casscf_fcidump.py @@ -1,9 +1,29 @@ +import numpy as np import pytest from examples.exmp054_casscf_fcidump.job import run_exmp054 from opi.input.structures import Structure +def _naive_eri(two_electron: dict, norb: int) -> np.ndarray: + """Independent, unvectorized reference implementation for differential testing.""" + tensor = np.zeros((norb,) * 4) + for (i, j, k, ll), val in two_electron.items(): + a, b, c, d = i - 1, j - 1, k - 1, ll - 1 + for p, q, r, s in [ + (a, b, c, d), + (b, a, c, d), + (a, b, d, c), + (b, a, d, c), + (c, d, a, b), + (d, c, a, b), + (c, d, b, a), + (d, c, b, a), + ]: + tensor[p, q, r, s] = val + return tensor + + @pytest.mark.examples @pytest.mark.orca def test_exmp054_casscf_fcidump(example_input_file, tmp_path) -> None: @@ -15,6 +35,20 @@ def test_exmp054_casscf_fcidump(example_input_file, tmp_path) -> None: # Run the example in tmp_path output = run_exmp054(structure=structure, working_dir=tmp_path) - fcidump_file = output.get_outfile().with_suffix(".fcidump") + fcidump = output.get_fcidump() + assert fcidump is not None + fcidump_file = fcidump.path assert fcidump_file.exists() and fcidump_file.is_file() + + norb = fcidump.norb + assert fcidump.hcore_matrix.shape == (norb, norb) + assert fcidump.eri_tensor.shape == (norb, norb, norb, norb) + + # hcore must be symmetric + assert np.allclose(fcidump.hcore_matrix, fcidump.hcore_matrix.T) + + # differential test: vectorized eri_tensor must match an independently + # written, unvectorized reference implementation on this same real data; + # catches an all-zero or misplaced tensor at every symmetry-equivalent position + assert np.allclose(fcidump.eri_tensor, _naive_eri(fcidump.two_electron, norb)) diff --git a/tests/unit/test_output_fcidump_parser.py b/tests/unit/test_output_fcidump_parser.py new file mode 100644 index 00000000..4c3d6241 --- /dev/null +++ b/tests/unit/test_output_fcidump_parser.py @@ -0,0 +1,146 @@ +from pathlib import Path + +import pytest + +from opi.output.fcidump import Fcidump + + +@pytest.mark.unit +@pytest.mark.output +def test_parse_fcidump_header(tmp_path: Path) -> None: + fcidump_text = """ + &FCI + NORB= 2, NELEC= 2, MS2= 0, + ORBSYM=1,1, + ISYM=0, + / + 0.5000000000 1 1 1 1 + 0.1000000000 2 1 1 1 + 0.2000000000 2 2 1 1 + 0.3000000000 2 2 2 2 + -0.1234567890 1 1 0 0 + 0.0987654321 2 1 0 0 + 0.0500000000 2 2 0 0 + -1.2345678901 0 0 0 0 + """ + fci_file = tmp_path / "test.fcidump" + fci_file.write_text(fcidump_text) + + dump = Fcidump.from_file(fci_file) + + assert dump.norb == 2 + assert dump.nelec == 2 + assert dump.ms2 == 0 + assert dump.orbsym == [1, 1] + assert dump.isym == 0 + assert pytest.approx(dump.e_nuc) == -1.2345678901 + + +@pytest.mark.unit +@pytest.mark.output +def test_get_int() -> None: + """`_get_int()` extracts integer header values, tolerating whitespace and mixed-case keys.""" + header = "&FCI\n NORB= 12, NELEC=8,\n ms2 = 0,\n/" + + assert Fcidump._get_int("NORB", header) == 12 + assert Fcidump._get_int("NELEC", header) == 8 + # key matching is case-insensitive + assert Fcidump._get_int("MS2", header) == 0 + + +@pytest.mark.unit +@pytest.mark.output +def test_get_int_missing_key_raises() -> None: + """`_get_int()` raises a ValueError naming the key when it is absent from the header.""" + with pytest.raises(ValueError, match="ISYM"): + Fcidump._get_int("ISYM", "&FCI\n NORB= 2,\n/") + + +@pytest.mark.unit +@pytest.mark.output +def test_get_int_negative_value_raises() -> None: + """`_get_int()` only accepts non-negative integers, so a negative value raises a ValueError.""" + with pytest.raises(ValueError, match="MS2"): + Fcidump._get_int("MS2", "&FCI\n MS2= -2,\n/") + + +@pytest.mark.unit +@pytest.mark.output +def test_get_int_list() -> None: + """`_get_int_list()` parses comma-separated integers, tolerating whitespace and mixed case.""" + assert Fcidump._get_int_list("ORBSYM", "ORBSYM=1,1,2,1,") == [1, 1, 2, 1] + # tolerates whitespace around the separators and is case-insensitive + assert Fcidump._get_int_list("ORBSYM", "orbsym = 1 , 2 , 3,") == [1, 2, 3] + # a single-entry list (norb=1) must parse as well + assert Fcidump._get_int_list("ORBSYM", "ORBSYM=1,") == [1] + + +@pytest.mark.unit +@pytest.mark.output +def test_get_int_list_missing_key_raises() -> None: + """`_get_int_list()` raises a ValueError naming the key when it is absent from the header.""" + with pytest.raises(ValueError, match="ORBSYM"): + Fcidump._get_int_list("ORBSYM", "&FCI\n NORB= 2,\n/") + + +@pytest.mark.unit +@pytest.mark.output +def test_get_int_list_negative_value_raises() -> None: + """`_get_int_list()` only accepts non-negative integers, so a negative entry raises a ValueError.""" + with pytest.raises(ValueError, match="ORBSYM"): + Fcidump._get_int_list("ORBSYM", "ORBSYM=-1,1,") + + +@pytest.mark.unit +@pytest.mark.output +def test_hcore_matrix_shape_and_symmetry() -> None: + dump = Fcidump( + norb=2, + nelec=2, + ms2=0, + orbsym=[1, 1], + isym=0, + one_electron={(1, 1): -0.5, (2, 1): 0.1, (2, 2): -0.3}, + ) + mat = dump.hcore_matrix + + assert mat.shape == (2, 2) + assert pytest.approx(mat[0, 0]) == -0.5 + assert pytest.approx(mat[1, 1]) == -0.3 + # off-diagonal must be symmetric + assert pytest.approx(mat[1, 0]) == 0.1 + assert pytest.approx(mat[0, 1]) == 0.1 + + +@pytest.mark.unit +@pytest.mark.output +def test_eri_tensor_empty_two_electron() -> None: + """Without two-electron integrals, eri_tensor must be all-zero instead of raising.""" + dump = Fcidump(norb=2, nelec=2, ms2=0, orbsym=[1, 1], isym=0) + + assert dump.eri_tensor.shape == (2, 2, 2, 2) + assert not dump.eri_tensor.any() + + +@pytest.mark.unit +@pytest.mark.output +def test_eri_tensor_shape_and_symmetry() -> None: + dump = Fcidump( + norb=2, + nelec=2, + ms2=0, + orbsym=[1, 1], + isym=0, + two_electron={(1, 1, 1, 1): 0.5, (2, 1, 1, 1): 0.1, (2, 2, 1, 1): 0.2, (2, 2, 2, 2): 0.3}, + ) + tensor = dump.eri_tensor + + assert tensor.shape == (2, 2, 2, 2) + assert pytest.approx(tensor[0, 0, 0, 0]) == 0.5 + assert pytest.approx(tensor[1, 1, 1, 1]) == 0.3 + # check 8-fold symmetry for (2,1,1,1) -> index (1,0,0,0) + val = 0.1 + assert pytest.approx(tensor[1, 0, 0, 0]) == val # (ij|kl) + assert pytest.approx(tensor[0, 1, 0, 0]) == val # (ji|kl) + assert pytest.approx(tensor[0, 0, 1, 0]) == val # (kl|ij) + assert pytest.approx(tensor[0, 0, 0, 1]) == val # (lk|ij)