FCIDump Driver

.pylintdict

@@ -140,6 +140,7 @@ eigenstate
 eigrange
 eigs
 eigvecs
+einact
 einsum
 els
 endif
@@ -159,6 +160,7 @@ excitations
 expr
 factorizers
 factr
+fcidump
 fcompiler
 fermionic
 FermionicOperator
@@ -224,6 +226,7 @@ indvar
 init
 initializer
 initio
+inline
 innerproduct
 inreg
 instantiations
@@ -239,6 +242,7 @@ ising
 isinstance
 iso
 isub
+isym
 iteratively
 izaac
 jac
@@ -249,6 +253,7 @@ jt
 jth
 jw
 kaicher
+ket
 killoran
 kingma
 kitaev
@@ -317,6 +322,7 @@ msq
 multiclass
 multinomial
 multiprocess
+namelist
 nan
 narray
 nasdaq
@@ -326,6 +332,7 @@ ndarray
 ndarray's
 negoro
 nelder
+nelec
 neq
 nevals
 newtons's
@@ -338,6 +345,7 @@ nn
 noancilla
 noint
 noqa
+norb
 norbs
 nosignatures
 np
@@ -350,14 +358,16 @@ objval
 occ
 oe
 ok
-oneee
+onee
 online
 onwards
 oplus
 optim
 optimizer's
 optimizers
+orbsym
 org
+outpath
 overfit
 params
 parentname
@@ -581,4 +591,4 @@ zzz
 ucc
 uccd
 UCCS
-vir
+vir

qiskit/chemistry/drivers/__init__.py

@@ -2,7 +2,7 @@
 
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2018, 2019.
+# (C) Copyright IBM 2018, 2020.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -78,6 +78,7 @@
 
 """
 from ._basedriver import BaseDriver, UnitsType, HFMethodType
+from .fcidumpd import FCIDumpDriver
 from .gaussiand import GaussianDriver
 from .hdf5d import HDF5Driver
 from .psi4d import PSI4Driver
@@ -87,6 +88,7 @@
 __all__ = ['BaseDriver',
            'UnitsType',
            'HFMethodType',
+           'FCIDumpDriver',
            'GaussianDriver',
            'HDF5Driver',
            'PSI4Driver',

qiskit/chemistry/drivers/fcidumpd/__init__.py

@@ -0,0 +1,22 @@
+# -*- coding: utf-8 -*-
+
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""FCIDump driver package.
+
+Contains tools to parse and dump FCIDump files.
+"""
+
+from .fcidumpdriver import FCIDumpDriver
+
+__all__ = ['FCIDumpDriver']

qiskit/chemistry/drivers/fcidumpd/dumper.py

@@ -0,0 +1,122 @@
+# -*- coding: utf-8 -*-
+
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""FCIDump dumper."""
+
+from typing import List, Optional
+from io import TextIOWrapper
+import itertools
+import numpy as np
+
+
+def dump(outpath: str, norb: int, nelec: int, hijs: List[float], hijkls: List[float], einact: float,
+         ms2: int = 0, orbsym: Optional[List[int]] = None, isym: int = 1
+         ) -> None:
+    # pylint: disable=wrong-spelling-in-docstring
+    """Generates a FCIDump output.
+
+    Args:
+        outpath: Path to the output file.
+        norb: The number of orbitals.
+        nelec: The number of electrons.
+        hijs: The pair of alpha and beta 1-electron integrals. The latter may be None.
+        hijkls: The triplet of alpha/alpha, beta/alpha and beta/beta 2-electron integrals. The
+            latter two may be None.
+        einact: The inactive energy.
+        ms2: 2*S, where S is the spin quantum number.
+        orbsym: A list of spatial symmetries of the orbitals.
+        isym: The spatial symmetry of the wave function.
+    """
+    hij, hij_b = hijs
+    hijkl, hijkl_ba, hijkl_bb = hijkls
+    # assert that either all beta variables are None or all of them are not
+    assert all([h is None for h in [hij_b, hijkl_ba, hijkl_bb]]) \
+        or all([h is not None for h in [hij_b, hijkl_ba, hijkl_bb]])
+    assert norb == hij.shape[0] == hijkl.shape[0]
+    mos = range(norb)
+    with open(outpath, 'w') as outfile:
+        # print header
+        outfile.write('&FCI NORB={:4d},NELEC={:4d},MS2={:4d}\n'.format(norb, nelec, ms2))
+        if orbsym is None:
+            outfile.write(' ORBSYM=' + '1,'*norb + '\n')
+        else:
+            assert len(orbsym) == norb
+            outfile.write(' ORBSYM=' + ','.join(orbsym) + '\n')
+        outfile.write(' ISYM={:d},\n/&END\n'.format(isym))
+        # append 2e integrals
+        _dump_2e_ints(hijkl, mos, outfile)
+        if hijkl_ba is not None:
+            _dump_2e_ints(hijkl_ba.transpose(), mos, outfile, beta=1)
+        if hijkl_bb is not None:
+            _dump_2e_ints(hijkl_bb, mos, outfile, beta=2)
+        # append 1e integrals
+        _dump_1e_ints(hij, mos, outfile)
+        if hij_b is not None:
+            _dump_1e_ints(hij_b, mos, outfile, beta=True)
+        # TODO append MO energies (last three indices are 0)
+        # append inactive energy
+        _write_to_outfile(outfile, einact, (0, 0, 0, 0))
+
+
+def _dump_1e_ints(hij: List[float], mos: List[int], outfile: TextIOWrapper,
+                  beta: bool = False) -> None:
+    idx_offset = 1 if not beta else 1+len(mos)
+    hij_elements = set()
+    for i, j in itertools.product(mos, repeat=2):
+        if i == j:
+            _write_to_outfile(outfile, hij[i][j], (i+idx_offset, j+idx_offset, 0, 0))
+            continue
+        if (j, i) in hij_elements and np.isclose(hij[i][j], hij[j][i]):
+            continue
+        _write_to_outfile(outfile, hij[i][j], (i+idx_offset, j+idx_offset, 0, 0))
+        hij_elements.add((i, j))
+
+
+def _dump_2e_ints(hijkl: List[float], mos: List[int], outfile: TextIOWrapper,
+                  beta: int = 0) -> None:
+    idx_offsets = [1, 1]
+    for b in range(beta):
+        idx_offsets[1-b] += len(mos)
+    hijkl_elements = set()
+    # pylint: disable=invalid-name
+    for elem in itertools.product(mos, repeat=4):
+        if np.isclose(hijkl[elem], 0.0, atol=1e-14):
+            continue
+        if len(set(elem)) == 1:
+            _write_to_outfile(outfile, hijkl[elem], (*[e+idx_offsets[0] for e in elem[:2]],
+                                                     *[e+idx_offsets[1] for e in elem[2:]]))
+            continue
+        if beta != 1 and elem[::-1] in hijkl_elements and \
+                np.isclose(hijkl[elem], hijkl[elem[::-1]]):
+            continue
+        bra_perms = set(itertools.permutations(elem[:2]))
+        ket_perms = set(itertools.permutations(elem[2:]))
+        if beta == 1:
+            permutations = itertools.product(bra_perms, ket_perms)
+        else:
+            permutations = itertools.chain(
+                itertools.product(bra_perms, ket_perms),
+                itertools.product(ket_perms, bra_perms)
+            )
+        for perm in {e1 + e2 for e1, e2 in permutations}:
+            if perm in hijkl_elements and np.isclose(hijkl[elem], hijkl[perm]):
+                break
+        else:
+            _write_to_outfile(outfile, hijkl[elem], (*[e+idx_offsets[0] for e in elem[:2]],
+                                                     *[e+idx_offsets[1] for e in elem[2:]]))
+            hijkl_elements.add(elem)
+
+
+def _write_to_outfile(outfile: str, value: float, indices: List[int]):
+    outfile.write('{:23.16E}{:4d}{:4d}{:4d}{:4d}\n'.format(value, *indices))

qiskit/chemistry/drivers/fcidumpd/fcidumpdriver.py

@@ -0,0 +1,106 @@
+# -*- coding: utf-8 -*-
+
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2020.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""FCIDump Driver."""
+
+from typing import List, Optional
+from qiskit.chemistry.drivers import BaseDriver
+from qiskit.chemistry import QiskitChemistryError, QMolecule
+from .dumper import dump
+from .parser import parse
+
+
+class FCIDumpDriver(BaseDriver):
+    """Python implementation of an FCIDump driver.
+
+    The FCIDump format is partially defined in Knowles1989.
+
+    References:
+        Knowles1989: Peter J. Knowles, Nicholas C. Handy,
+            A determinant based full configuration interaction program,
+            Computer Physics Communications, Volume 54, Issue 1, 1989, Pages 75-83,
+            ISSN 0010-4655, https://doi.org/10.1016/0010-4655(89)90033-7.
+    """
+
+    def __init__(self, fcidump_input: str, atoms: Optional[List[str]] = None) -> None:
+        """
+        Args:
+            fcidump_input: Path to the FCIDump file.
+            atoms: Allows to specify the atom list of the molecule. If it is provided, the created
+                QMolecule instance will permit frozen core Hamiltonians. This list must consist of
+                valid atom symbols.
+
+        Raises:
+            QiskitChemistryError: If ``fcidump_input`` is not a string or if ``atoms`` is not a list
+                of valid atomic symbols as specified in ``QMolecule``.
+        """
+        super().__init__()
+
+        if not isinstance(fcidump_input, str):
+            raise QiskitChemistryError(
+                "The fcidump_input must be str, not '{}'".format(fcidump_input))
+        self._fcidump_input = fcidump_input
+
+        if atoms and not isinstance(atoms, list) \
+                and not all([sym in QMolecule.symbols for sym in atoms]):
+            raise QiskitChemistryError(
+                "The atoms must be a list of valid atomic symbols, not '{}'".format(atoms))
+        self.atoms = atoms
+
+    def run(self) -> QMolecule:
+        """Constructs a QMolecule instance out of a FCIDump file.
+
+        Returns:
+            A QMolecule instance populated with a minimal set of required data.
+        """
+        fcidump_data = parse(self._fcidump_input)
+
+        q_mol = QMolecule()
+
+        q_mol.nuclear_repulsion_energy = fcidump_data.get('ecore', None)
+        q_mol.num_orbitals = fcidump_data.get('NORB')
+        q_mol.multiplicity = fcidump_data.get('MS2', 0) + 1
+        q_mol.charge = 0  # ensures QMolecule.log() works
+        q_mol.num_beta = (fcidump_data.get('NELEC') - (q_mol.multiplicity - 1)) // 2
+        q_mol.num_alpha = fcidump_data.get('NELEC') - q_mol.num_beta
+        if self.atoms is not None:
+            q_mol.num_atoms = len(self.atoms)
+            q_mol.atom_symbol = self.atoms
+            q_mol.atom_xyz = [[float('NaN')] * 3] * len(self.atoms)  # ensures QMolecule.log() works
+
+        q_mol.mo_onee_ints = fcidump_data.get('hij', None)
+        q_mol.mo_onee_ints_b = fcidump_data.get('hij_b', None)
+        q_mol.mo_eri_ints = fcidump_data.get('hijkl', None)
+        q_mol.mo_eri_ints_bb = fcidump_data.get('hijkl_bb', None)
+        q_mol.mo_eri_ints_ba = fcidump_data.get('hijkl_ba', None)
+
+        return q_mol
+
+    @staticmethod
+    def dump(q_mol: QMolecule, outpath: str, orbsym: Optional[List[int]] = None,
+             isym: int = 1) -> None:
+        """Convenience method to produce an FCIDump output file.
+
+        Args:
+            outpath: Path to the output file.
+            q_mol: QMolecule data to be dumped. It is assumed that the nuclear_repulsion_energy in
+                this QMolecule instance contains the inactive core energy.
+            orbsym: A list of spatial symmetries of the orbitals.
+            isym: The spatial symmetry of the wave function.
+        """
+        dump(outpath,
+             q_mol.num_orbitals, q_mol.num_alpha + q_mol.num_beta,
+             (q_mol.mo_onee_ints, q_mol.mo_onee_ints_b),
+             (q_mol.mo_eri_ints, q_mol.mo_eri_ints_ba, q_mol.mo_eri_ints_bb),
+             q_mol.nuclear_repulsion_energy, ms2=q_mol.multiplicity - 1, orbsym=orbsym, isym=isym)