NonSingletCI Class, AP1roGSDGeneralizedSen-o Wfn

pyci/__init__.py

@@ -18,7 +18,7 @@
 from .pyci import __version__, c_long, c_ulong, c_double, sparse_op
 
 from .pyci import secondquant_op, wavefunction, one_spin_wfn, two_spin_wfn
-from .pyci import doci_wfn, fullci_wfn, genci_wfn, sparse_op
+from .pyci import doci_wfn, fullci_wfn, genci_wfn, nonsingletci_wfn, sparse_op
 from .pyci import get_num_threads, set_num_threads, popcnt, ctz
 from .pyci import compute_overlap, compute_rdms, compute_transition_rdms
 from .pyci import add_hci, compute_enpt2
@@ -45,6 +45,7 @@
     "doci_wfn",
     "fullci_wfn",
     "genci_wfn",
+    "nonsingletci_wfn",
     "sparse_op",
     "get_num_threads",
     "set_num_threads",

pyci/fanci/__init__.py

@@ -10,6 +10,7 @@
     "AP1roG",
     "DetRatio",
     "pCCDS",
+    "AP1roGeneralizedSeno",
 ]
 
 
@@ -18,3 +19,4 @@
 from .ap1rog import AP1roG
 from .detratio import DetRatio
 from .pccds import pCCDS
+from .ap1rogen import AP1roGeneralizedSeno

pyci/fanci/ap1rogen.py

@@ -0,0 +1,195 @@
+r"""
+FanCI AP1roGeneralizedseno module for AP1roGSDGeneralized_sen-o wavefunction.
+
+"""
+
+from typing import Any, Union
+
+import numpy as np
+
+import pyci
+
+from ..pyci import AP1roGeneralizedSenoObjective
+from .fanci import FanCI
+import pdb
+
+
+__all___ = [
+    "AP1roGeneralizedSeno",
+]
+
+
+class AP1roGeneralizedSeno(FanCI):
+    r"""
+    DOC
+    """
+
+    def __init__(
+        self,
+        ham: pyci.hamiltonian,
+        nocc: int,
+        nproj: int = None,
+        wfn: pyci.nonsingletci_wfn = None,
+        fill: str = 'excitation',
+        **kwargs: Any,
+    ) -> None:
+        r"""
+        Initialize the FanCI problem.
+
+        Parameters
+        ----------
+        ham : pyci.hamiltonian
+            PyCI Hamiltonian.
+        nocc : int
+            Number of occupied orbitals.
+        nproj : int, optional
+            Number of determinants in projection ("P") space.
+        wfn : pyci.nonsingletci_wfn
+            If specified, this PyCI wave function defines the projection ("P") space.
+        fill : 
+        kwargs : Any, optional
+            Additional keyword arguments for base FanCI class.
+
+        """
+        # SEE OTHER FANCI WFNS
+
+        if not isinstance(ham, pyci.hamiltonian):
+            raise TypeError(f"Invalid `ham` type `{type(ham)}`; must be `pyci.hamiltonian`")
+
+        nparam = nocc * (ham.nbasis - nocc) + (2 * nocc) * (2 * (ham.nbasis - nocc)) + 1 #less params considering we added singles as well
+        nproj = nparam if nproj is None else nproj
+
+        # if nproj > nparam:
+        #     raise ValueError("nproj cannot be greater than the size of the space")
+
+        if wfn is None:
+            # wfn = pyci.doci_wfn(ham.nbasis, nocc, nocc)
+            # wfn.add_excited_dets(1) # add pair excited determinants
+            wfn = pyci.fullci_wfn(ham.nbasis, nocc, nocc)
+            # exc=0 ensures addint HF determinannt first
+            pyci.add_excitations(wfn,1,2,3,4)
+            print("Printing FCI wfn dets: ")
+            for i, sd in enumerate(wfn.to_occ_array()):
+                sd = np.array(sd)
+                print(wfn.to_det_array()[i],np.concatenate((sd[0],sd[1]+ham.nbasis)))
+
+            wfn = pyci.nonsingletci_wfn(wfn)
+
+        elif not isinstance(wfn, pyci.nonsingletci_wfn):
+            raise TypeError(f"Invalid `wfn` type `{type(wfn)}`; must be `pyci.nonsingletci_wfn`")
+        elif wfn.nocc != nocc:
+            raise ValueError(f"wfn.nocc does not match `nocc={nocc}` parameter")
+
+
+        # Initialize base class
+        FanCI.__init__(self, ham, wfn, nproj, nparam, **kwargs)
+        print("Nonsingletci nbasis: ",wfn.nbasis)
+        print("\n\nPrinting nonsingletci wfn dets: ")
+        for sd in (self._sspace):
+            print(sd)
+            # sd = np.array(sd)
+            # print(np.concatenate((sd[0],sd[1]+ham.nbasis)))
+        print("Done printing\n\n")
+
+        # Assign reference occupations
+        #ref_occs_up = np.arange(nocc_up, dtype=pyci.c_long)
+        #ref_occs_dn = np.arange(nocc_dn, dtype=pyci.c_long)
+
+        # Save sub-class-specific attributes
+        #self._ref_occs = [ref_occs_up, ref_occs_dn]
+
+        # Initiazlize C++ extension
+        try:
+            norm_det = kwargs["norm_det"]
+            idx_det_cons = np.asarray([elem[0] for elem in norm_det], dtype=int)
+            det_cons = np.asarray([elem[1] for elem in norm_det], dtype=int)
+        except KeyError:
+            idx_det_cons = None
+            det_cons = None
+
+        try:
+            norm_param = kwargs["norm_param"]
+            idx_param_cons = np.asarray([elem[0] for elem in norm_param], dtype=int)
+            param_cons = np.asarray([elem[1] for elem in norm_param], dtype=int)
+        except KeyError:
+            idx_param_cons = None
+            param_cons = None
+
+        self._cext = AP1roGeneralizedSenoObjective(
+            self._ci_op, self._wfn,
+            idx_det_cons=idx_det_cons, det_cons=det_cons,
+            idx_param_cons=idx_param_cons, param_cons=param_cons,
+        )
+
+    def compute_overlap(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Compute the FanCI overlap vector.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Parameter array, [p_0, p_1, ..., p_n].
+
+        Returns
+        -------
+        ovlp : np.ndarray
+            Overlap array.
+
+        """
+        return self._cext.overlap(x)
+
+    def compute_overlap_deriv(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Compute the FanCI overlap derivative matrix.
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Parameter array, [p_0, p_1, ..., p_n].
+
+        Returns
+        -------
+        ovlp : np.ndarray
+            Overlap derivative array.
+
+        """
+        return self._cext.d_overlap(x)
+
+    def compute_objective(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Compute the FanCI objective function.
+
+            f : x[k] -> y[n]
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Parameter array, [p_0, p_1, ..., p_n, E].
+
+        Returns
+        -------
+        obj : np.ndarray
+            Objective vector.
+
+        """
+        return self._cext.objective(self._ci_op, x)
+
+    def compute_jacobian(self, x: np.ndarray) -> np.ndarray:
+        r"""
+        Compute the Jacobian of the FanCI objective function.
+
+            j : x[k] -> y[n, k]
+
+        Parameters
+        ----------
+        x : np.ndarray
+            Parameter array, [p_0, p_1, ..., p_n, E].
+
+        Returns
+        -------
+        jac : np.ndarray
+            Jacobian matrix.
+
+        """
+        return self._cext.jacobian(self._ci_op, x)
+

pyci/fanci/fanci.py

@@ -13,6 +13,8 @@
 
 from scipy.optimize import OptimizeResult, least_squares, root
 
+from scipy.optimize import minimize
+
 import pyci
 
 
@@ -200,11 +202,28 @@ def __init__(
         wfn = fill_wavefunction(wfn, nproj, fill)
 
         # Compute CI matrix operator with nproj rows and len(wfn) columns
-        ci_op = pyci.sparse_op(ham, wfn, nrow=nproj, ncol=len(wfn), symmetric=False)
+        if type(wfn).__name__ == "nonsingletci_wfn":
+            ci_op = pyci.sparse_op(ham, wfn, nrow=nproj, ncol=len(wfn), symmetric=False, wfntype="nonsingletci")
+        else:
+            ci_op = pyci.sparse_op(ham, wfn, nrow=nproj, ncol=len(wfn), symmetric=False)
 
         # Compute arrays of occupations
         sspace = wfn.to_occ_array()
+
+        # if type(wfn).__name__ == "nonsingletci_wfn":
+        #     pspace = sspace[np.r_[0:13, 21:29, 37:45, 53:60, 66:72, 75:77, 81, 99:102, 104:108, 109:117, 126:129, 131:135, 
+        #               136:144, 171:183, 189:195, 201:206, 209, 221:224, 225:231, 2316:239, 240:246, 261:267, 
+        #               273:278, 282:286, 288, 300:302, 303:309, 314:316, 317:323, 337:344, 348:352, 356:359, 
+        #               365:369, 371:375, 381:384, 388:391, 393:395, 400:403, 405:408, 413:416, 418:420, 422, 
+        #               424:426, 427, 430]]
+        #     # nproj = len(pspace)
+        # else:
         pspace = sspace[:nproj]
+        print("len(sspace), len(pspace): ", len(sspace), len(pspace))
+        # dec_sspace = wfn.to_det_array()
+        # print("Printing sspace")
+        # for dec, bin in zip(dec_sspace, sspace):
+        #     print("\n", dec, bin)
 
         # Assign attributes to instance
         self._nequation = nequation
@@ -271,6 +290,10 @@ def optimize(
         elif mode == "root":
             opt_args = f, x0
             optimizer = root
+        # FIXME: For BFGS to work, objective function must return a scalar value
+        elif mode == "bfgs":
+            opt_args = f, x0
+            optimizer = minimize
         else:
             raise ValueError("invalid mode parameter")