ci & mp: spin != 0

pyblock2/uc/ci.py

@@ -21,12 +21,17 @@
 
 """Arbitrary order Configuration Interaction."""
 
+
 class CI(lib.StreamObject):
     def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None, ci_order=2):
-        from pyscf.scf import hf
+        from pyscf.scf import hf, addons
 
         if isinstance(mf, hf.KohnShamDFT):
             raise RuntimeError("CI Warning: The first argument mf is a DFT object.")
+        if isinstance(mf, scf.rohf.ROHF):
+            lib.logger.warn(mf, 'RCI method does not support ROHF method. ROHF object '
+                            'is converted to UHF object and UCI method is called.')
+            mf = addons.convert_to_uhf(mf)
         if mo_coeff is None:
             mo_coeff = mf.mo_coeff
         if mo_occ is None:
@@ -66,7 +71,9 @@ def kernel(self, ci0=None, max_cycle=50, tol=1e-8):
                 g2e_symm=1,
             )
             nc, nv = (n_elec - spin) // 2, ncas - (n_elec + spin) // 2
-            idx = np.concatenate([np.mgrid[: nc + nv][nc:], np.mgrid[: nc + nv][:nc]])
+            idx = np.concatenate(
+                [np.mgrid[:ncas][ncas - nv :], np.mgrid[:ncas][: ncas - nv]]
+            )
             orb_sym = np.array(orb_sym)[idx]
             h1e = h1e[idx, :][:, idx]
             g2e = g2e[idx, :][:, idx][:, :, idx][:, :, :, idx]
@@ -84,12 +91,18 @@ def kernel(self, ci0=None, max_cycle=50, tol=1e-8):
                 g2e_symm=1,
             )
             nc, nv = (n_elec - abs(spin)) // 2, ncas - (n_elec + abs(spin)) // 2
-            idx = np.concatenate([np.mgrid[: nc + nv][nc:], np.mgrid[: nc + nv][:nc]])
+            idx = np.concatenate(
+                [np.mgrid[:ncas][ncas - nv :], np.mgrid[:ncas][: ncas - nv]]
+            )
             orb_sym = np.array(orb_sym)[idx]
             h1e = tuple(x[idx, :][:, idx] for x in h1e)
             g2e_aa, g2e_ab, g2e_bb = tuple(
                 x[idx, :][:, idx][:, :, idx][:, :, :, idx] for x in g2e
             )
+            h1e = tuple(np.ascontiguousarray(x) for x in h1e)
+            g2e_aa, g2e_ab, g2e_bb = tuple(
+                np.ascontiguousarray(x) for x in (g2e_aa, g2e_ab, g2e_bb)
+            )
             bs = b.sz
             VectorSX = b.VectorSZ
             target = b.SZ(n_elec, spin, 0)
@@ -104,7 +117,13 @@ def kernel(self, ci0=None, max_cycle=50, tol=1e-8):
             VectorSX([target]), False, ncas, b.VectorUInt8(orb_sym), fd, 0
         )
         big.drt = bs.DRT(
-            big.drt.n_sites, big.drt.get_init_qs(), big.drt.orb_sym, nc, nv, n_ex=n_elec
+            big.drt.n_sites,
+            big.drt.get_init_qs(),
+            big.drt.orb_sym,
+            nc,
+            nv,
+            n_ex=n_elec,
+            single_ref=True,
         )
 
         if self.verbose >= 5:
@@ -179,6 +198,14 @@ def diag(self) -> np.ndarray:
             self.ci = kcis if self.nroots > 1 else kcis[0]
             conv = ndav < max_cycle
 
+        if self.ci_order == 0:
+            lib.logger.note(
+                self,
+                "E(%4s) =%s",
+                "CI%d" % 0,
+                "%20.16g" * self.nroots % ((self.e_hf,) * self.nroots),
+            )
+
         return conv, self.e_corr + self.e_hf, self.ci
 
 
@@ -189,4 +216,13 @@ def diag(self) -> np.ndarray:
     mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="ccpvdz", symmetry="d2h", verbose=3)
     mf = scf.RHF(mol).run(conv_tol=1e-14)
     ci.CISD(mf).run()
-    CI(mf, ci_order=4).run()
+    CI(mf, ci_order=2).run()
+
+    mol = gto.M(atom="N 0 0 0; N 0 0 1.1", basis="ccpvdz", symmetry="c1", verbose=3, spin=2)
+    mf = scf.RHF(mol).run(conv_tol=1e-14)
+    ci.CISD(mf).run()
+    CI(mf, ci_order=2).run()
+
+    mf = scf.UHF(mol).run(conv_tol=1e-14)
+    ci.CISD(mf).run()
+    CI(mf, ci_order=2).run()

pyblock2/uc/mp.py

@@ -21,12 +21,20 @@
 
 """Arbitrary order Møller-Plesset Perturbation Theory."""
 
+
 class MP(lib.StreamObject):
     def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None, mp_order=2):
-        from pyscf.scf import hf
+        from pyscf.scf import hf, addons
 
         if isinstance(mf, hf.KohnShamDFT):
             raise RuntimeError("MP Warning: The first argument mf is a DFT object.")
+        if isinstance(mf, scf.rohf.ROHF):
+            lib.logger.warn(
+                mf,
+                "RMP method does not support ROHF method. ROHF object "
+                "is converted to UHF object and UMP method is called.",
+            )
+            mf = addons.convert_to_uhf(mf)
         if mo_coeff is None:
             mo_coeff = mf.mo_coeff
         if mo_occ is None:
@@ -66,7 +74,9 @@ def kernel(self, max_cycle=50, tol=1e-8):
                 g2e_symm=1,
             )
             nc, nv = (n_elec - spin) // 2, ncas - (n_elec + spin) // 2
-            idx = np.concatenate([np.mgrid[: nc + nv][nc:], np.mgrid[: nc + nv][:nc]])
+            idx = np.concatenate(
+                [np.mgrid[:ncas][ncas - nv :], np.mgrid[:ncas][: ncas - nv]]
+            )
             orb_sym = np.array(orb_sym)[idx]
             h1e = h1e[idx, :][:, idx]
             g2e = g2e[idx, :][:, idx][:, :, idx][:, :, :, idx]
@@ -99,7 +109,9 @@ def kernel(self, max_cycle=50, tol=1e-8):
                 g2e_symm=1,
             )
             nc, nv = (n_elec - abs(spin)) // 2, ncas - (n_elec + abs(spin)) // 2
-            idx = np.concatenate([np.mgrid[: nc + nv][nc:], np.mgrid[: nc + nv][:nc]])
+            idx = np.concatenate(
+                [np.mgrid[:ncas][ncas - nv :], np.mgrid[:ncas][: ncas - nv]]
+            )
             orb_sym = np.array(orb_sym)[idx]
             h1e = tuple(x[idx, :][:, idx] for x in h1e)
             g2e_aa, g2e_ab, g2e_bb = tuple(
@@ -110,35 +122,51 @@ def kernel(self, max_cycle=50, tol=1e-8):
             target = b.SZ(n_elec, spin, 0)
             vacuum = b.SZ(0)
             g2e = g2e_aa, g2e_bb, g2e_ab
+            h1e = tuple(np.ascontiguousarray(x) for x in h1e)
+            g2e = tuple(
+                np.ascontiguousarray(x) for x in g2e
+            )
             fd = b.FCIDUMP()
             fd.initialize_sz(ncas, n_elec, spin, 0, ecore, h1e, g2e)
 
+            if spin == 0:
+                nva = slice(nv, None)
+                nvb = slice(nv, None)
+            elif spin > 0:
+                nva = slice(nv, None)
+                nvb = slice(nv, nv + nc)
+            else:
+                nva = slice(nv, nv + nc)
+                nvb = slice(nv, None)
+
             f1e = (
                 h1e[0]
-                + 0.5 * np.einsum("mnjj->mn", g2e[0][:, :, nv:, nv:])
-                + 0.5 * np.einsum("mnjj->mn", g2e[2][:, :, nv:, nv:])
-                - 0.5 * np.einsum("mjjn->mn", g2e[0][:, nv:, nv:, :]),
+                + np.einsum("mnjj->mn", g2e[0][:, :, nva, nva])
+                + np.einsum("mnjj->mn", g2e[2][:, :, nvb, nvb])
+                - np.einsum("mjjn->mn", g2e[0][:, nva, nva, :]),
                 h1e[1]
-                + 0.5 * np.einsum("mnjj->mn", g2e[1][:, :, nv:, nv:])
-                + 0.5 * np.einsum("jjmn->mn", g2e[2][nv:, nv:, :, :])
-                - 0.5 * np.einsum("mjjn->mn", g2e[1][:, nv:, nv:, :]),
+                + np.einsum("mnjj->mn", g2e[1][:, :, nvb, nvb])
+                + np.einsum("jjmn->mn", g2e[2][nva, nva, :, :])
+                - np.einsum("mjjn->mn", g2e[1][:, nvb, nvb, :]),
             )
             e0 = (
                 ecore
-                + np.einsum("jj->", f1e[0][nv:, nv:])
-                + np.einsum("jj->", f1e[1][nv:, nv:])
+                + np.einsum("jj->", f1e[0][nva, nva])
+                + np.einsum("jj->", f1e[1][nvb, nvb])
             )
             e1 = (
-                -0.5 * np.einsum("iijj->", g2e[0][nv:, nv:, nv:, nv:])
-                - 0.5 * np.einsum("iijj->", g2e[1][nv:, nv:, nv:, nv:])
-                - 1.0 * np.einsum("iijj->", g2e[2][nv:, nv:, nv:, nv:])
-                + 0.5 * np.einsum("ijji->", g2e[0][nv:, nv:, nv:, nv:])
-                + 0.5 * np.einsum("ijji->", g2e[1][nv:, nv:, nv:, nv:])
+                -0.5 * np.einsum("iijj->", g2e[0][nva, nva, nva, nva])
+                - 0.5 * np.einsum("iijj->", g2e[1][nvb, nvb, nvb, nvb])
+                - 1.0 * np.einsum("iijj->", g2e[2][nva, nva, nvb, nvb])
+                + 0.5 * np.einsum("ijji->", g2e[0][nva, nva, nva, nva])
+                + 0.5 * np.einsum("ijji->", g2e[1][nvb, nvb, nvb, nvb])
             )
 
+            f1e = tuple(np.ascontiguousarray(x) for x in f1e)
+
             fd0 = b.FCIDUMP()
             fd0.initialize_sz(
-                ncas, n_elec, spin, 0, ecore - e0, f1e, (0.0 * x for x in g2e)
+                ncas, n_elec, spin, 0, ecore - e0, f1e, tuple(0.0 * x for x in g2e)
             )
             fd1 = b.FCIDUMP()
             fd1.initialize_sz(
@@ -152,7 +180,13 @@ def kernel(self, max_cycle=50, tol=1e-8):
             VectorSX([target]), False, ncas, b.VectorUInt8(orb_sym), fd, 0
         )
         big.drt = bs.DRT(
-            big.drt.n_sites, big.drt.get_init_qs(), big.drt.orb_sym, nc, nv, n_ex=n_elec
+            big.drt.n_sites,
+            big.drt.get_init_qs(),
+            big.drt.orb_sym,
+            nc,
+            nv,
+            n_ex=n_elec,
+            single_ref=True,
         )
 
         big0 = bs.DRTBigSite(
@@ -258,6 +292,14 @@ def diag(self) -> np.ndarray:
             self.e_corr = sum(ecorrs[: self.mp_order + 1]) - self.e_hf
             conv = niter < max_cycle
             self.ci = wfns
+
+        if self.mp_order == 1:
+            lib.logger.note(
+                self,
+                "E(%4s) =%s",
+                "MP%d" % 1,
+                "%20.16g" % sum(ecorrs[:2]),
+            )
 
         return conv, self.e_corr + self.e_hf, self.ci
 
@@ -270,3 +312,17 @@ def diag(self) -> np.ndarray:
     mf = scf.RHF(mol).run(conv_tol=1e-14)
     mp.MP2(mf).run()
     MP(mf, mp_order=2).run()
+
+    mol = gto.M(
+        atom="N 0 0 0; N 0 0 1.1", basis="ccpvdz", symmetry="c1", verbose=3, spin=2
+    )
+    mf = scf.RHF(mol).run(conv_tol=1e-14)
+    mp.MP2(mf).run()
+    MP(mf, mp_order=2).run() # not okay
+
+    mol = gto.M(
+        atom="N 0 0 0; N 0 0 1.1", basis="ccpvdz", symmetry="c1", verbose=3, spin=2
+    )
+    mf = scf.UHF(mol).run(conv_tol=1e-14)
+    mp.MP2(mf).run()
+    MP(mf, mp_order=2).run()

src/big_site/drt_big_site.hpp

@@ -57,6 +57,7 @@ template <typename S, ElemOpTypes T = ElemT<S>::value> struct DRT {
     vector<array<LL, 5>> xs;
     int n_sites, n_init_qs;
     int n_core, n_virt, n_ex;
+    bool single_ref;
     DRT() : n_sites(0), n_init_qs(0), n_core(0), n_virt(0), n_ex(0) {}
     DRT(int16_t a, int16_t b, int16_t c,
         typename S::pg_t ipg = (typename S::pg_t)0,
@@ -70,9 +71,9 @@ template <typename S, ElemOpTypes T = ElemT<S>::value> struct DRT {
         : DRT(n_sites, vector<S>{q}, orb_sym, n_core, n_virt, n_ex) {}
     DRT(int n_sites, const vector<S> &init_qs,
         const vector<typename S::pg_t> &orb_sym = vector<typename S::pg_t>(),
-        int n_core = 0, int n_virt = 0, int n_ex = 0)
+        int n_core = 0, int n_virt = 0, int n_ex = 0, bool single_ref = false)
         : n_sites(n_sites), orb_sym(orb_sym), n_init_qs((int)init_qs.size()),
-          n_core(n_core), n_virt(n_virt), n_ex(n_ex) {
+          n_core(n_core), n_virt(n_virt), n_ex(n_ex), single_ref(single_ref) {
         if (T == ElemOpTypes::SU2 || T == ElemOpTypes::SZ) {
             for (auto &q : init_qs) {
                 abc.push_back(array<int16_t, 4>{
@@ -89,7 +90,9 @@ template <typename S, ElemOpTypes T = ElemT<S>::value> struct DRT {
     virtual ~DRT() = default;
     int n_rows() const { return (int)abc.size(); }
     void initialize() {
-        int nc = this->n_core, nv = this->n_virt, nx = this->n_ex;
+        const int nc = this->n_core, nv = this->n_virt, nx = this->n_ex;
+        const bool sr = this->single_ref;
+        const int sp = (int)(n_init_qs > 0 && abc[0][1] > 0);
         abc.resize(n_init_qs);
         pgs.resize(n_init_qs);
         auto make_abc = [](int16_t a, int16_t b, int16_t c, int16_t t,
@@ -126,16 +129,19 @@ template <typename S, ElemOpTypes T = ElemT<S>::value> struct DRT {
                 return false;
             }
         };
-        auto make_abct = [&make_abc, &nc, &nv](int k, int16_t a, int16_t b,
-                                               int16_t c, int16_t t,
-                                               int16_t d) -> array<int16_t, 4> {
+        auto make_abct = [&make_abc, &nc, &nv, &sr, &sp](int k, int16_t a, int16_t b, int16_t c,
+                               int16_t t, int16_t d) -> array<int16_t, 4> {
             array<int16_t, 4> r = make_abc(a, b, c, t, d);
-            r[3] = (int16_t)(k < nv || k > nc + nv
-                                 ? 0
-                                 : (k < nc + nv ? (int)t
-                                                : max(0, nc + nc -
-                                                             (a + a + abs(b) -
-                                                              (d + 1) / 2))));
+            if (sr && k >= nc + nv)
+                    r[3] = (int16_t)(t + (sp ? d >= 2 : d == 1 || d == 3));
+            else
+                r[3] =
+                    (int16_t)(k < nv || k > nc + nv
+                                  ? 0
+                                  : (k < nc + nv ? (int)t
+                                                 : max(0, nc + nc -
+                                                              (a + a + abs(b) -
+                                                               (d + 1) / 2))));
             return r;
         };
         auto allow_abct = [&allow_abc, &nc, &nv,
@@ -326,7 +332,7 @@ template <typename S, ElemOpTypes T = ElemT<S>::value> struct DRT {
     }
     shared_ptr<DRT<S>> operator^(int n_ex_new) const {
         return make_shared<DRT<S>>(n_sites, get_init_qs(), orb_sym, n_core,
-                                   n_virt, n_ex_new);
+                                   n_virt, n_ex_new, single_ref);
     }
     vector<LL> operator>>(const shared_ptr<DRT<S>> &other) const {
         vector<vector<int>> pbr(2, vector<int>(1, 0)),

src/core/integral.hpp

@@ -554,7 +554,7 @@ template <typename FL> struct FCIDUMP {
     // Two-electron integrals can be three general rank-4 arrays
     // or 8-fold, 8-fold, 4-fold rank-1 arrays
     virtual ~FCIDUMP() = default;
-    virtual void initialize_sz(uint16_t n_sites, uint16_t n_elec, uint16_t twos,
+    virtual void initialize_sz(uint16_t n_sites, uint16_t n_elec, int16_t twos,
                                uint16_t isym, typename const_fl_type<FL>::FL e,
                                const FL *ta, size_t lta, const FL *tb,
                                size_t ltb, const FL *va, size_t lva,
@@ -622,10 +622,10 @@ template <typename FL> struct FCIDUMP {
     }
     // Initialize integrals: SU(2) case
     // Two-electron integrals can be general rank-4 array or 8-fold rank-1 array
-    virtual void initialize_su2(uint16_t n_sites, uint16_t n_elec,
-                                uint16_t twos, uint16_t isym,
-                                typename const_fl_type<FL>::FL e, const FL *t,
-                                size_t lt, const FL *v, size_t lv) {
+    virtual void initialize_su2(uint16_t n_sites, uint16_t n_elec, int16_t twos,
+                                uint16_t isym, typename const_fl_type<FL>::FL e,
+                                const FL *t, size_t lt, const FL *v,
+                                size_t lv) {
         params.clear();
         ts.clear();
         vs.clear();
@@ -666,10 +666,9 @@ template <typename FL> struct FCIDUMP {
         uhf = false;
     }
     // Initialize with only h1e integral
-    virtual void initialize_h1e(uint16_t n_sites, uint16_t n_elec,
-                                uint16_t twos, uint16_t isym,
-                                typename const_fl_type<FL>::FL e, const FL *t,
-                                size_t lt) {
+    virtual void initialize_h1e(uint16_t n_sites, uint16_t n_elec, int16_t twos,
+                                uint16_t isym, typename const_fl_type<FL>::FL e,
+                                const FL *t, size_t lt) {
         params.clear();
         ts.clear();
         vs.clear();
@@ -1076,9 +1075,7 @@ template <typename FL> struct FCIDUMP {
         return error;
     }
     // Target 2S or 2Sz
-    uint16_t twos() const {
-        return (uint16_t)Parsing::to_int(params.at("ms2"));
-    }
+    int16_t twos() const { return (int16_t)Parsing::to_int(params.at("ms2")); }
     // Number of sites
     uint16_t n_sites() const {
         return (uint16_t)Parsing::to_int(params.at("norb"));

src/core/integral_compressed.hpp

@@ -430,7 +430,7 @@ template <typename FL> struct CompressedFCIDUMP : FCIDUMP<FL> {
     // Initialize integrals: U(1) case
     // Two-electron integrals can be three general rank-4 arrays
     // or 8-fold, 8-fold, 4-fold rank-1 arrays
-    void initialize_sz(uint16_t n_sites, uint16_t n_elec, uint16_t twos,
+    void initialize_sz(uint16_t n_sites, uint16_t n_elec, int16_t twos,
                        uint16_t isym, typename const_fl_type<FL>::FL e,
                        istream &ta, size_t lta, istream &tb, size_t ltb,
                        istream &va, size_t lva, istream &vb, size_t lvb,

src/pybind/pybind_big_site.hpp

@@ -206,6 +206,7 @@ template <typename S> void bind_drt_big_site(py::module &m) {
         .def_readwrite("n_core", &DRT<S>::n_core)
         .def_readwrite("n_virt", &DRT<S>::n_virt)
         .def_readwrite("n_ex", &DRT<S>::n_ex)
+        .def_readwrite("single_ref", &DRT<S>::single_ref)
         .def(py::init<>())
         .def(py::init<int16_t, int16_t, int16_t>(), py::arg("a"), py::arg("b"),
              py::arg("c"))
@@ -236,6 +237,11 @@ template <typename S> void bind_drt_big_site(py::module &m) {
                       int, int, int>(),
              py::arg("n_sites"), py::arg("init_qs"), py::arg("orb_sym"),
              py::arg("n_core"), py::arg("n_virt"), py::arg("n_ex"))
+        .def(py::init<int, const vector<S> &, const vector<typename S::pg_t> &,
+                      int, int, int, bool>(),
+             py::arg("n_sites"), py::arg("init_qs"), py::arg("orb_sym"),
+             py::arg("n_core"), py::arg("n_virt"), py::arg("n_ex"),
+             py::arg("single_ref"))
         .def_property_readonly("n_rows", &DRT<S>::n_rows)
         .def("initialize", &DRT<S>::initialize)
         .def("get_init_qs", &DRT<S>::get_init_qs)