diff --git a/pyblock2/uc/ci.py b/pyblock2/uc/ci.py index 0a96d1e3..caf81722 100644 --- a/pyblock2/uc/ci.py +++ b/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() diff --git a/pyblock2/uc/mp.py b/pyblock2/uc/mp.py index a840b839..4af07e85 100644 --- a/pyblock2/uc/mp.py +++ b/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() diff --git a/src/big_site/drt_big_site.hpp b/src/big_site/drt_big_site.hpp index 26be2f62..689a20fa 100644 --- a/src/big_site/drt_big_site.hpp +++ b/src/big_site/drt_big_site.hpp @@ -57,6 +57,7 @@ template ::value> struct DRT { vector> 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 ::value> struct DRT { : DRT(n_sites, vector{q}, orb_sym, n_core, n_virt, n_ex) {} DRT(int n_sites, const vector &init_qs, const vector &orb_sym = vector(), - 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{ @@ -89,7 +90,9 @@ template ::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 ::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 { + 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 { array 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 ::value> struct DRT { } shared_ptr> operator^(int n_ex_new) const { return make_shared>(n_sites, get_init_qs(), orb_sym, n_core, - n_virt, n_ex_new); + n_virt, n_ex_new, single_ref); } vector operator>>(const shared_ptr> &other) const { vector> pbr(2, vector(1, 0)), diff --git a/src/core/integral.hpp b/src/core/integral.hpp index ea582629..d83b2004 100644 --- a/src/core/integral.hpp +++ b/src/core/integral.hpp @@ -554,7 +554,7 @@ template 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 e, const FL *ta, size_t lta, const FL *tb, size_t ltb, const FL *va, size_t lva, @@ -622,10 +622,10 @@ template 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 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 e, + const FL *t, size_t lt, const FL *v, + size_t lv) { params.clear(); ts.clear(); vs.clear(); @@ -666,10 +666,9 @@ template 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 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 e, + const FL *t, size_t lt) { params.clear(); ts.clear(); vs.clear(); @@ -1076,9 +1075,7 @@ template 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")); diff --git a/src/core/integral_compressed.hpp b/src/core/integral_compressed.hpp index f5fa10ce..229b5ef2 100644 --- a/src/core/integral_compressed.hpp +++ b/src/core/integral_compressed.hpp @@ -430,7 +430,7 @@ template struct CompressedFCIDUMP : FCIDUMP { // 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 e, istream &ta, size_t lta, istream &tb, size_t ltb, istream &va, size_t lva, istream &vb, size_t lvb, diff --git a/src/pybind/pybind_big_site.hpp b/src/pybind/pybind_big_site.hpp index bfbd27ad..49549b8e 100644 --- a/src/pybind/pybind_big_site.hpp +++ b/src/pybind/pybind_big_site.hpp @@ -206,6 +206,7 @@ template void bind_drt_big_site(py::module &m) { .def_readwrite("n_core", &DRT::n_core) .def_readwrite("n_virt", &DRT::n_virt) .def_readwrite("n_ex", &DRT::n_ex) + .def_readwrite("single_ref", &DRT::single_ref) .def(py::init<>()) .def(py::init(), py::arg("a"), py::arg("b"), py::arg("c")) @@ -236,6 +237,11 @@ template 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 &, const vector &, + 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::n_rows) .def("initialize", &DRT::initialize) .def("get_init_qs", &DRT::get_init_qs) diff --git a/src/pybind/pybind_core.hpp b/src/pybind/pybind_core.hpp index 9e738d2e..09d9ec50 100644 --- a/src/pybind/pybind_core.hpp +++ b/src/pybind/pybind_core.hpp @@ -2904,7 +2904,7 @@ template void bind_fl_matrix(py::module &m) { }) .def("initialize_sz", [](FCIDUMP *self, uint16_t n_sites, uint16_t n_elec, - uint16_t twos, uint16_t isym, FL e, const py::tuple &t, + int16_t twos, uint16_t isym, FL e, const py::tuple &t, const py::tuple &v) { assert(t.size() == 2 && v.size() == 3); py::array_t ta = t[0].cast>(); @@ -3088,7 +3088,7 @@ template void bind_fl_matrix(py::module &m) { ift.close(); }) .def("initialize_sz", [](CompressedFCIDUMP *self, uint16_t n_sites, - uint16_t n_elec, uint16_t twos, uint16_t isym, + uint16_t n_elec, int16_t twos, uint16_t isym, FL e, const py::tuple &t, const py::tuple &v) { assert(t.size() == 2 && v.size() == 3); if (py::isinstance>(t[0])) {