diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc index e2ccaef75..fc0419ea3 100644 --- a/forte/api/forte_python_module.cc +++ b/forte/api/forte_python_module.cc @@ -181,7 +181,6 @@ PYBIND11_MODULE(_forte, m) { }, "Return the cumulants of the RDMs in a spinorbital basis. Spinorbitals follow the ordering " "abab..."); - m.def("get_gas_occupation", &get_gas_occupation); m.def("get_ci_occupation_patterns", &get_ci_occupation_patterns); @@ -250,6 +249,61 @@ PYBIND11_MODULE(_forte, m) { .def("compute_gradient", &MASTER_DSRG::compute_gradient, "Compute the DSRG gradient") .def("compute_Heff_actv", &MASTER_DSRG::compute_Heff_actv, "Return the DSRG dressed ActiveSpaceIntegrals") + .def("save_Heff_full_ambit", [](MASTER_DSRG& self) { + const auto Heff = self.save_Heff_full(); + return py::make_tuple(Heff.first, Heff.second.at(0), Heff.second.at(1)); + }) + .def("save_Heff_full", [](MASTER_DSRG& self) { + const auto Heff = self.save_Heff_full(); + return py::make_tuple(Heff.first, blockedtensor_to_np(Heff.second.at(0)), blockedtensor_to_np(Heff.second.at(1))); + }) + .def("update_Heff_full", [](MASTER_DSRG& self, double H0, BlockedTensor H1, BlockedTensor H2, std::shared_ptr rdms) { + const auto Heff = self.update_Heff_full(H0, H1, H2, rdms); + return py::make_tuple(Heff.first, blockedtensor_to_np(Heff.second.at(0)), blockedtensor_to_np(Heff.second.at(1))); + }) + .def("compute_mbar", &MASTER_DSRG::compute_mbar, "compute full transformed dipole integrals") + .def("compute_Mbar0_full", &MASTER_DSRG::compute_Mbar0_full, "Return full transformed zero-body dipole integrals") + .def("compute_Mbar1_full", [](MASTER_DSRG& self) { + const auto Mbar1 = self.compute_Mbar1_full(); + return py::make_tuple(blockedtensor_to_np(Mbar1.at(0)), blockedtensor_to_np(Mbar1.at(1)), + blockedtensor_to_np(Mbar1.at(2))); + }) + .def("compute_Mbar2_full", [](MASTER_DSRG& self) { + const auto Mbar2 = self.compute_Mbar2_full(); + return py::make_tuple(blockedtensor_to_np(Mbar2.at(0)), blockedtensor_to_np(Mbar2.at(1)), + blockedtensor_to_np(Mbar2.at(2))); + }) + .def("get_gamma1", [](MASTER_DSRG& self) { + const auto gamma1 = self.get_gamma1(); + return blockedtensor_to_np(gamma1); + }) + .def("get_eta1", [](MASTER_DSRG& self) { + const auto eta1 = self.get_eta1(); + return blockedtensor_to_np(eta1); + }) + .def("get_lambda2", [](MASTER_DSRG& self) { + const auto lambda2 = self.get_lambda2(); + return blockedtensor_to_np(lambda2); + }) + .def("get_lambda3", [](MASTER_DSRG& self) { + const auto lambda3 = self.get_lambda3(); + py::dict pyrdm; + pyrdm[py::str("aaaaaa")] = ambit_to_np(lambda3.at(0)); + pyrdm[py::str("aaAaaA")] = ambit_to_np(lambda3.at(1)); + pyrdm[py::str("aAAaAA")] = ambit_to_np(lambda3.at(2)); + pyrdm[py::str("AAAAAA")] = ambit_to_np(lambda3.at(3)); + return pyrdm; + }) + .def("get_lambda4", [](MASTER_DSRG& self) { + const auto lambda4 = self.get_lambda4(); + py::dict pyrdm; + pyrdm[py::str("aaaaaaaa")] = ambit_to_np(lambda4.at(0)); + pyrdm[py::str("aaaAaaaA")] = ambit_to_np(lambda4.at(1)); + pyrdm[py::str("aaAAaaAA")] = ambit_to_np(lambda4.at(2)); + pyrdm[py::str("aAAAaAAA")] = ambit_to_np(lambda4.at(3)); + pyrdm[py::str("AAAAAAAA")] = ambit_to_np(lambda4.at(4)); + return pyrdm; + }) .def("deGNO_DMbar_actv", &MASTER_DSRG::deGNO_DMbar_actv, "Return the DSRG dressed dipole integrals") .def("nuclear_dipole", &MASTER_DSRG::nuclear_dipole, @@ -268,13 +322,18 @@ PYBIND11_MODULE(_forte, m) { .def("set_ci_vectors", &MASTER_DSRG::set_ci_vectors, "Set the CI eigenvector for DSRG-MRPT2 analytic gradients") .def("set_active_space_solver", &MASTER_DSRG::set_active_space_solver, - "Set the shared pointer for ActiveSpaceSolver"); + "Set the shared pointer for ActiveSpaceSolver") + .def("set_rdms", &MASTER_DSRG::set_rdms, "Set external RDMs (dangerous, use with caution)"); // export SADSRG py::class_(m, "SADSRG") .def("compute_energy", &SADSRG::compute_energy, "Compute the DSRG energy") .def("compute_Heff_actv", &SADSRG::compute_Heff_actv, "Return the DSRG dressed ActiveSpaceIntegrals") + .def("compute_Heff_full", [](SADSRG& self) { + const auto Heff = self.compute_Heff_full(); + return py::make_tuple(blockedtensor_to_np(Heff.at(0)), blockedtensor_to_np(Heff.at(1))); + }) .def("compute_mp_eff_actv", &SADSRG::compute_mp_eff_actv, "Return the DSRG dressed ActiveMultipoleIntegrals") .def("set_Uactv", &SADSRG::set_Uactv, "Ua"_a, @@ -299,7 +358,27 @@ PYBIND11_MODULE(_forte, m) { py::class_(m, "MRDSRG_SO") .def("compute_energy", &MRDSRG_SO::compute_energy, "Compute DSRG energy") .def("compute_Heff_actv", &MRDSRG_SO::compute_Heff_actv, - "Return the DSRG dressed ActiveSpaceIntegrals"); + "Return the DSRG dressed ActiveSpaceIntegrals") + .def("save_Heff_full", [](MRDSRG_SO& self) { + const auto Heff = self.save_Heff_full(); + return py::make_tuple(Heff.first, blockedtensor_to_np(Heff.second.at(0)), blockedtensor_to_np(Heff.second.at(1))); + }) + .def("get_gamma1", [](MRDSRG_SO& self) { + const auto gamma1 = self.get_gamma1(); + return blockedtensor_to_np(gamma1); + }) + .def("get_eta1", [](MRDSRG_SO& self) { + const auto eta1 = self.get_eta1(); + return blockedtensor_to_np(eta1); + }) + .def("get_lambda2", [](MRDSRG_SO& self) { + const auto lambda2 = self.get_lambda2(); + return blockedtensor_to_np(lambda2); + }) + .def("get_lambda3", [](MRDSRG_SO& self) { + const auto lambda3 = self.get_lambda3(); + return blockedtensor_to_np(lambda3); + }); // export DSRG_MRPT spin-adapted code py::class_(m, "DSRG_MRPT") diff --git a/forte/api/integrals_api.cc b/forte/api/integrals_api.cc index 5cf53de43..edb7c765f 100644 --- a/forte/api/integrals_api.cc +++ b/forte/api/integrals_api.cc @@ -98,6 +98,7 @@ void export_ForteIntegrals(py::module& m) { .def( "print_ints", [](const ForteIntegrals& ints) { psi::outfile->Printf("\n%s", ints.repr().c_str()); }, - "Print the integrals"); + "Print the integrals") + .def("mo_dipole_ints", &ForteIntegrals::mo_dipole_ints, "mo_dipole_ints"); } } // namespace forte diff --git a/forte/base_classes/rdms.cc b/forte/base_classes/rdms.cc index 3c48c0a72..30aa7b2ad 100644 --- a/forte/base_classes/rdms.cc +++ b/forte/base_classes/rdms.cc @@ -42,11 +42,13 @@ std::shared_ptr RDMs::build(size_t max_rdm_level, size_t n_orbs, RDMsType std::vector dims1(2, n_orbs); std::vector dims2(4, n_orbs); std::vector dims3(6, n_orbs); + std::vector dims4(8, n_orbs); std::shared_ptr rdms; if (type == RDMsType::spin_dependent) { - ambit::Tensor g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb; + ambit::Tensor g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb, g4aaaa, g4aaab, + g4aabb, g4abbb, g4bbbb; if (max_rdm_level > 0) { g1a = ambit::Tensor::build(ambit::CoreTensor, "g1a", dims1); g1b = ambit::Tensor::build(ambit::CoreTensor, "g1b", dims1); @@ -62,6 +64,13 @@ std::shared_ptr RDMs::build(size_t max_rdm_level, size_t n_orbs, RDMsType g3abb = ambit::Tensor::build(ambit::CoreTensor, "g3abb", dims3); g3bbb = ambit::Tensor::build(ambit::CoreTensor, "g3bbb", dims3); } + if (max_rdm_level > 3) { + g4aaaa = ambit::Tensor::build(ambit::CoreTensor, "g4aaaa", dims4); + g4aaab = ambit::Tensor::build(ambit::CoreTensor, "g4aaab", dims4); + g4aabb = ambit::Tensor::build(ambit::CoreTensor, "g4aabb", dims4); + g4abbb = ambit::Tensor::build(ambit::CoreTensor, "g4abbb", dims4); + g4bbbb = ambit::Tensor::build(ambit::CoreTensor, "g4bbbb", dims4); + } if (max_rdm_level < 1) { rdms = std::make_shared(); @@ -69,9 +78,13 @@ std::shared_ptr RDMs::build(size_t max_rdm_level, size_t n_orbs, RDMsType rdms = std::make_shared(g1a, g1b); } else if (max_rdm_level == 2) { rdms = std::make_shared(g1a, g1b, g2aa, g2ab, g2bb); - } else { + } else if (max_rdm_level == 3) { rdms = std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb); + } else { + rdms = + std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, + g3bbb, g4aaaa, g4aaab, g4aabb, g4abbb, g4bbbb); } } else { ambit::Tensor g1, g2, g3; @@ -91,8 +104,11 @@ std::shared_ptr RDMs::build(size_t max_rdm_level, size_t n_orbs, RDMsType rdms = std::make_shared(g1); } else if (max_rdm_level == 2) { rdms = std::make_shared(g1, g2); - } else { + } else if (max_rdm_level == 3) { rdms = std::make_shared(g1, g2, g3); + } else { + throw std::runtime_error( + "RDMs::build: max_rdm_level > 3 not implemented for spin-free RDMs."); } } return rdms; @@ -105,7 +121,8 @@ std::shared_ptr RDMs::build_from_disk(size_t max_rdm_level, RDMsType type, std::string prefix = (filename_prefix.empty() ? "" : filename_prefix + "."); if (type == RDMsType::spin_dependent) { - ambit::Tensor g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb; + ambit::Tensor g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb, g4aaaa, g4aaab, + g4aabb, g4abbb, g4bbbb; if (max_rdm_level > 0) { g1a = ambit::load_tensor(prefix + "g1a.bin"); g1b = ambit::load_tensor(prefix + "g1b.bin"); @@ -121,6 +138,13 @@ std::shared_ptr RDMs::build_from_disk(size_t max_rdm_level, RDMsType type, g3abb = ambit::load_tensor(prefix + "g3abb.bin"); g3bbb = ambit::load_tensor(prefix + "g3bbb.bin"); } + if (max_rdm_level > 3) { + g4aaaa = ambit::load_tensor(prefix + "g4aaaa.bin"); + g4aaab = ambit::load_tensor(prefix + "g4aaab.bin"); + g4aabb = ambit::load_tensor(prefix + "g4aabb.bin"); + g4abbb = ambit::load_tensor(prefix + "g4abbb.bin"); + g4bbbb = ambit::load_tensor(prefix + "g4bbbb.bin"); + } if (max_rdm_level < 1) { rdms = std::make_shared(); @@ -128,9 +152,13 @@ std::shared_ptr RDMs::build_from_disk(size_t max_rdm_level, RDMsType type, rdms = std::make_shared(g1a, g1b); } else if (max_rdm_level == 2) { rdms = std::make_shared(g1a, g1b, g2aa, g2ab, g2bb); - } else { + } else if (max_rdm_level == 3) { rdms = std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb); + } else { + rdms = + std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, + g3bbb, g4aaaa, g4aaab, g4aabb, g4abbb, g4bbbb); } } else { ambit::Tensor g1, g2, g3; @@ -143,6 +171,10 @@ std::shared_ptr RDMs::build_from_disk(size_t max_rdm_level, RDMsType type, if (max_rdm_level > 2) { g3 = ambit::load_tensor(prefix + "g3.bin"); } + if (max_rdm_level > 3) { + throw std::runtime_error( + "RDMs::build_from_disk: max_rdm_level > 3 not implemented for spin-free RDMs."); + } if (max_rdm_level < 1) { rdms = std::make_shared(); @@ -313,6 +345,353 @@ ambit::Tensor RDMs::make_cumulant_L3abb(const ambit::Tensor& g1a, const ambit::T return L3abb; } +ambit::Tensor RDMs::make_cumulant_L4aaaa(const ambit::Tensor& L1a, const ambit::Tensor& L2aa, + const ambit::Tensor& L3aaa, const ambit::Tensor& g4aaaa) { + timer t("make_cumulant_L4aaaa"); + + auto L4aaaa = g4aaaa.clone(); + /// 1111 contributions (24 terms) + L4aaaa("pqrstuvw") -= L1a("pw") * L1a("qv") * L1a("ru") * L1a("st"); + L4aaaa("pqrstuvw") += L1a("pw") * L1a("qv") * L1a("rt") * L1a("su"); + L4aaaa("pqrstuvw") += L1a("pw") * L1a("qu") * L1a("rv") * L1a("st"); + L4aaaa("pqrstuvw") -= L1a("pw") * L1a("qu") * L1a("rt") * L1a("sv"); + L4aaaa("pqrstuvw") -= L1a("pw") * L1a("qt") * L1a("rv") * L1a("su"); + L4aaaa("pqrstuvw") += L1a("pw") * L1a("qt") * L1a("ru") * L1a("sv"); + L4aaaa("pqrstuvw") += L1a("pv") * L1a("qw") * L1a("ru") * L1a("st"); + L4aaaa("pqrstuvw") -= L1a("pv") * L1a("qw") * L1a("rt") * L1a("su"); + L4aaaa("pqrstuvw") -= L1a("pv") * L1a("qu") * L1a("rw") * L1a("st"); + L4aaaa("pqrstuvw") += L1a("pv") * L1a("qu") * L1a("rt") * L1a("sw"); + L4aaaa("pqrstuvw") += L1a("pv") * L1a("qt") * L1a("rw") * L1a("su"); + L4aaaa("pqrstuvw") -= L1a("pv") * L1a("qt") * L1a("ru") * L1a("sw"); + L4aaaa("pqrstuvw") -= L1a("pu") * L1a("qw") * L1a("rv") * L1a("st"); + L4aaaa("pqrstuvw") += L1a("pu") * L1a("qw") * L1a("rt") * L1a("sv"); + L4aaaa("pqrstuvw") += L1a("pu") * L1a("qv") * L1a("rw") * L1a("st"); + L4aaaa("pqrstuvw") -= L1a("pu") * L1a("qv") * L1a("rt") * L1a("sw"); + L4aaaa("pqrstuvw") -= L1a("pu") * L1a("qt") * L1a("rw") * L1a("sv"); + L4aaaa("pqrstuvw") += L1a("pu") * L1a("qt") * L1a("rv") * L1a("sw"); + L4aaaa("pqrstuvw") += L1a("pt") * L1a("qw") * L1a("rv") * L1a("su"); + L4aaaa("pqrstuvw") -= L1a("pt") * L1a("qw") * L1a("ru") * L1a("sv"); + L4aaaa("pqrstuvw") -= L1a("pt") * L1a("qv") * L1a("rw") * L1a("su"); + L4aaaa("pqrstuvw") += L1a("pt") * L1a("qv") * L1a("ru") * L1a("sw"); + L4aaaa("pqrstuvw") += L1a("pt") * L1a("qu") * L1a("rw") * L1a("sv"); + L4aaaa("pqrstuvw") -= L1a("pt") * L1a("qu") * L1a("rv") * L1a("sw"); + + /// 211 contributions (72 terms) + L4aaaa("pqrstuvw") += L2aa("pqvw") * L1a("ru") * L1a("st"); + L4aaaa("pqrstuvw") -= L2aa("pqvw") * L1a("rt") * L1a("su"); + L4aaaa("pqrstuvw") -= L2aa("pquw") * L1a("rv") * L1a("st"); + L4aaaa("pqrstuvw") += L2aa("pquw") * L1a("rt") * L1a("sv"); + L4aaaa("pqrstuvw") += L2aa("pqtw") * L1a("rv") * L1a("su"); + L4aaaa("pqrstuvw") -= L2aa("pqtw") * L1a("ru") * L1a("sv"); + L4aaaa("pqrstuvw") += L2aa("pquv") * L1a("rw") * L1a("st"); + L4aaaa("pqrstuvw") -= L2aa("pquv") * L1a("rt") * L1a("sw"); + L4aaaa("pqrstuvw") -= L2aa("pqtv") * L1a("rw") * L1a("su"); + L4aaaa("pqrstuvw") += L2aa("pqtv") * L1a("ru") * L1a("sw"); + L4aaaa("pqrstuvw") += L2aa("pqtu") * L1a("rw") * L1a("sv"); + L4aaaa("pqrstuvw") -= L2aa("pqtu") * L1a("rv") * L1a("sw"); + L4aaaa("pqrstuvw") -= L2aa("prvw") * L1a("qu") * L1a("st"); + L4aaaa("pqrstuvw") += L2aa("prvw") * L1a("qt") * L1a("su"); + L4aaaa("pqrstuvw") += L2aa("pruw") * L1a("qv") * L1a("st"); + L4aaaa("pqrstuvw") -= L2aa("pruw") * L1a("qt") * L1a("sv"); + L4aaaa("pqrstuvw") -= L2aa("prtw") * L1a("qv") * L1a("su"); + L4aaaa("pqrstuvw") += L2aa("prtw") * L1a("qu") * L1a("sv"); + L4aaaa("pqrstuvw") -= L2aa("pruv") * L1a("qw") * L1a("st"); + L4aaaa("pqrstuvw") += L2aa("pruv") * L1a("qt") * L1a("sw"); + L4aaaa("pqrstuvw") += L2aa("prtv") * L1a("qw") * L1a("su"); + L4aaaa("pqrstuvw") -= L2aa("prtv") * L1a("qu") * L1a("sw"); + L4aaaa("pqrstuvw") -= L2aa("prtu") * L1a("qw") * L1a("sv"); + L4aaaa("pqrstuvw") += L2aa("prtu") * L1a("qv") * L1a("sw"); + L4aaaa("pqrstuvw") += L2aa("psvw") * L1a("qu") * L1a("rt"); + L4aaaa("pqrstuvw") -= L2aa("psvw") * L1a("qt") * L1a("ru"); + L4aaaa("pqrstuvw") -= L2aa("psuw") * L1a("qv") * L1a("rt"); + L4aaaa("pqrstuvw") += L2aa("psuw") * L1a("qt") * L1a("rv"); + L4aaaa("pqrstuvw") += L2aa("pstw") * L1a("qv") * L1a("ru"); + L4aaaa("pqrstuvw") -= L2aa("pstw") * L1a("qu") * L1a("rv"); + L4aaaa("pqrstuvw") += L2aa("psuv") * L1a("qw") * L1a("rt"); + L4aaaa("pqrstuvw") -= L2aa("psuv") * L1a("qt") * L1a("rw"); + L4aaaa("pqrstuvw") -= L2aa("pstv") * L1a("qw") * L1a("ru"); + L4aaaa("pqrstuvw") += L2aa("pstv") * L1a("qu") * L1a("rw"); + L4aaaa("pqrstuvw") += L2aa("pstu") * L1a("qw") * L1a("rv"); + L4aaaa("pqrstuvw") -= L2aa("pstu") * L1a("qv") * L1a("rw"); + L4aaaa("pqrstuvw") += L2aa("qrvw") * L1a("pu") * L1a("st"); + L4aaaa("pqrstuvw") -= L2aa("qrvw") * L1a("pt") * L1a("su"); + L4aaaa("pqrstuvw") -= L2aa("qruw") * L1a("pv") * L1a("st"); + L4aaaa("pqrstuvw") += L2aa("qruw") * L1a("pt") * L1a("sv"); + L4aaaa("pqrstuvw") += L2aa("qrtw") * L1a("pv") * L1a("su"); + L4aaaa("pqrstuvw") -= L2aa("qrtw") * L1a("pu") * L1a("sv"); + L4aaaa("pqrstuvw") += L2aa("qruv") * L1a("pw") * L1a("st"); + L4aaaa("pqrstuvw") -= L2aa("qruv") * L1a("pt") * L1a("sw"); + L4aaaa("pqrstuvw") -= L2aa("qrtv") * L1a("pw") * L1a("su"); + L4aaaa("pqrstuvw") += L2aa("qrtv") * L1a("pu") * L1a("sw"); + L4aaaa("pqrstuvw") += L2aa("qrtu") * L1a("pw") * L1a("sv"); + L4aaaa("pqrstuvw") -= L2aa("qrtu") * L1a("pv") * L1a("sw"); + L4aaaa("pqrstuvw") -= L2aa("qsvw") * L1a("pu") * L1a("rt"); + L4aaaa("pqrstuvw") += L2aa("qsvw") * L1a("pt") * L1a("ru"); + L4aaaa("pqrstuvw") += L2aa("qsuw") * L1a("pv") * L1a("rt"); + L4aaaa("pqrstuvw") -= L2aa("qsuw") * L1a("pt") * L1a("rv"); + L4aaaa("pqrstuvw") -= L2aa("qstw") * L1a("pv") * L1a("ru"); + L4aaaa("pqrstuvw") += L2aa("qstw") * L1a("pu") * L1a("rv"); + L4aaaa("pqrstuvw") -= L2aa("qsuv") * L1a("pw") * L1a("rt"); + L4aaaa("pqrstuvw") += L2aa("qsuv") * L1a("pt") * L1a("rw"); + L4aaaa("pqrstuvw") += L2aa("qstv") * L1a("pw") * L1a("ru"); + L4aaaa("pqrstuvw") -= L2aa("qstv") * L1a("pu") * L1a("rw"); + L4aaaa("pqrstuvw") -= L2aa("qstu") * L1a("pw") * L1a("rv"); + L4aaaa("pqrstuvw") += L2aa("qstu") * L1a("pv") * L1a("rw"); + L4aaaa("pqrstuvw") += L2aa("rsvw") * L1a("pu") * L1a("qt"); + L4aaaa("pqrstuvw") -= L2aa("rsvw") * L1a("pt") * L1a("qu"); + L4aaaa("pqrstuvw") -= L2aa("rsuw") * L1a("pv") * L1a("qt"); + L4aaaa("pqrstuvw") += L2aa("rsuw") * L1a("pt") * L1a("qv"); + L4aaaa("pqrstuvw") += L2aa("rstw") * L1a("pv") * L1a("qu"); + L4aaaa("pqrstuvw") -= L2aa("rstw") * L1a("pu") * L1a("qv"); + L4aaaa("pqrstuvw") += L2aa("rsuv") * L1a("pw") * L1a("qt"); + L4aaaa("pqrstuvw") -= L2aa("rsuv") * L1a("pt") * L1a("qw"); + L4aaaa("pqrstuvw") -= L2aa("rstv") * L1a("pw") * L1a("qu"); + L4aaaa("pqrstuvw") += L2aa("rstv") * L1a("pu") * L1a("qw"); + L4aaaa("pqrstuvw") += L2aa("rstu") * L1a("pw") * L1a("qv"); + L4aaaa("pqrstuvw") -= L2aa("rstu") * L1a("pv") * L1a("qw"); + + /// 22 contributions (18 terms) + L4aaaa("pqrstuvw") -= L2aa("pqvw") * L2aa("rstu"); + L4aaaa("pqrstuvw") += L2aa("pquw") * L2aa("rstv"); + L4aaaa("pqrstuvw") -= L2aa("pqtw") * L2aa("rsuv"); + L4aaaa("pqrstuvw") -= L2aa("pquv") * L2aa("rstw"); + L4aaaa("pqrstuvw") += L2aa("pqtv") * L2aa("rsuw"); + L4aaaa("pqrstuvw") -= L2aa("pqtu") * L2aa("rsvw"); + L4aaaa("pqrstuvw") += L2aa("prvw") * L2aa("qstu"); + L4aaaa("pqrstuvw") -= L2aa("pruw") * L2aa("qstv"); + L4aaaa("pqrstuvw") += L2aa("prtw") * L2aa("qsuv"); + L4aaaa("pqrstuvw") += L2aa("pruv") * L2aa("qstw"); + L4aaaa("pqrstuvw") -= L2aa("prtv") * L2aa("qsuw"); + L4aaaa("pqrstuvw") += L2aa("prtu") * L2aa("qsvw"); + L4aaaa("pqrstuvw") -= L2aa("psvw") * L2aa("qrtu"); + L4aaaa("pqrstuvw") += L2aa("psuw") * L2aa("qrtv"); + L4aaaa("pqrstuvw") -= L2aa("pstw") * L2aa("qruv"); + L4aaaa("pqrstuvw") -= L2aa("psuv") * L2aa("qrtw"); + L4aaaa("pqrstuvw") += L2aa("pstv") * L2aa("qruw"); + L4aaaa("pqrstuvw") -= L2aa("pstu") * L2aa("qrvw"); + + /// 31 contributions (16 terms) + L4aaaa("pqrstuvw") += L3aaa("pqruvw") * L1a("st"); + L4aaaa("pqrstuvw") -= L3aaa("pqrtvw") * L1a("su"); + L4aaaa("pqrstuvw") += L3aaa("pqrtuw") * L1a("sv"); + L4aaaa("pqrstuvw") -= L3aaa("pqrtuv") * L1a("sw"); + L4aaaa("pqrstuvw") -= L3aaa("pqsuvw") * L1a("rt"); + L4aaaa("pqrstuvw") += L3aaa("pqstvw") * L1a("ru"); + L4aaaa("pqrstuvw") -= L3aaa("pqstuw") * L1a("rv"); + L4aaaa("pqrstuvw") += L3aaa("pqstuv") * L1a("rw"); + L4aaaa("pqrstuvw") += L3aaa("prsuvw") * L1a("qt"); + L4aaaa("pqrstuvw") -= L3aaa("prstvw") * L1a("qu"); + L4aaaa("pqrstuvw") += L3aaa("prstuw") * L1a("qv"); + L4aaaa("pqrstuvw") -= L3aaa("prstuv") * L1a("qw"); + L4aaaa("pqrstuvw") -= L3aaa("qrsuvw") * L1a("pt"); + L4aaaa("pqrstuvw") += L3aaa("qrstvw") * L1a("pu"); + L4aaaa("pqrstuvw") -= L3aaa("qrstuw") * L1a("pv"); + L4aaaa("pqrstuvw") += L3aaa("qrstuv") * L1a("pw"); + + return L4aaaa; +} + +ambit::Tensor RDMs::make_cumulant_L4aaab(const ambit::Tensor& L1a, const ambit::Tensor& L1b, + const ambit::Tensor& L2aa, const ambit::Tensor& L2ab, + const ambit::Tensor& L3aaa, const ambit::Tensor& L3aab, + const ambit::Tensor& g4aaab) { + timer t("make_cumulant_L4aaab"); + + auto L4aaab = g4aaab.clone(); + /// 1111 contributions + L4aaab("pqrStuvW") += L1a("pv") * L1a("qu") * L1a("rt") * L1b("SW"); + L4aaab("pqrStuvW") -= L1a("pv") * L1a("qt") * L1a("ru") * L1b("SW"); + L4aaab("pqrStuvW") -= L1a("pu") * L1a("qv") * L1a("rt") * L1b("SW"); + L4aaab("pqrStuvW") += L1a("pu") * L1a("qt") * L1a("rv") * L1b("SW"); + L4aaab("pqrStuvW") += L1a("pt") * L1a("qv") * L1a("ru") * L1b("SW"); + L4aaab("pqrStuvW") -= L1a("pt") * L1a("qu") * L1a("rv") * L1b("SW"); + + /// 211 contributions + L4aaab("pqrStuvW") -= L2aa("pquv") * L1a("rt") * L1b("SW"); + L4aaab("pqrStuvW") += L2aa("pqtv") * L1a("ru") * L1b("SW"); + L4aaab("pqrStuvW") -= L2aa("pqtu") * L1a("rv") * L1b("SW"); + L4aaab("pqrStuvW") += L2aa("pruv") * L1a("qt") * L1b("SW"); + L4aaab("pqrStuvW") -= L2aa("prtv") * L1a("qu") * L1b("SW"); + L4aaab("pqrStuvW") += L2aa("prtu") * L1a("qv") * L1b("SW"); + L4aaab("pqrStuvW") += L2ab("pSvW") * L1a("qu") * L1a("rt"); + L4aaab("pqrStuvW") -= L2ab("pSvW") * L1a("qt") * L1a("ru"); + L4aaab("pqrStuvW") -= L2ab("pSuW") * L1a("qv") * L1a("rt"); + L4aaab("pqrStuvW") += L2ab("pSuW") * L1a("qt") * L1a("rv"); + L4aaab("pqrStuvW") += L2ab("pStW") * L1a("qv") * L1a("ru"); + L4aaab("pqrStuvW") -= L2ab("pStW") * L1a("qu") * L1a("rv"); + L4aaab("pqrStuvW") -= L2aa("qruv") * L1a("pt") * L1b("SW"); + L4aaab("pqrStuvW") += L2aa("qrtv") * L1a("pu") * L1b("SW"); + L4aaab("pqrStuvW") -= L2aa("qrtu") * L1a("pv") * L1b("SW"); + L4aaab("pqrStuvW") -= L2ab("qSvW") * L1a("pu") * L1a("rt"); + L4aaab("pqrStuvW") += L2ab("qSvW") * L1a("pt") * L1a("ru"); + L4aaab("pqrStuvW") += L2ab("qSuW") * L1a("pv") * L1a("rt"); + L4aaab("pqrStuvW") -= L2ab("qSuW") * L1a("pt") * L1a("rv"); + L4aaab("pqrStuvW") -= L2ab("qStW") * L1a("pv") * L1a("ru"); + L4aaab("pqrStuvW") += L2ab("qStW") * L1a("pu") * L1a("rv"); + L4aaab("pqrStuvW") += L2ab("rSvW") * L1a("pu") * L1a("qt"); + L4aaab("pqrStuvW") -= L2ab("rSvW") * L1a("pt") * L1a("qu"); + L4aaab("pqrStuvW") -= L2ab("rSuW") * L1a("pv") * L1a("qt"); + L4aaab("pqrStuvW") += L2ab("rSuW") * L1a("pt") * L1a("qv"); + L4aaab("pqrStuvW") += L2ab("rStW") * L1a("pv") * L1a("qu"); + L4aaab("pqrStuvW") -= L2ab("rStW") * L1a("pu") * L1a("qv"); + + /// 22 contributions + L4aaab("pqrStuvW") -= L2aa("pquv") * L2ab("rStW"); + L4aaab("pqrStuvW") += L2aa("pqtv") * L2ab("rSuW"); + L4aaab("pqrStuvW") -= L2aa("pqtu") * L2ab("rSvW"); + L4aaab("pqrStuvW") += L2aa("pruv") * L2ab("qStW"); + L4aaab("pqrStuvW") -= L2aa("prtv") * L2ab("qSuW"); + L4aaab("pqrStuvW") += L2aa("prtu") * L2ab("qSvW"); + L4aaab("pqrStuvW") -= L2ab("pSvW") * L2aa("qrtu"); + L4aaab("pqrStuvW") += L2ab("pSuW") * L2aa("qrtv"); + L4aaab("pqrStuvW") -= L2ab("pStW") * L2aa("qruv"); + + /// 31 contributions (16 terms) + L4aaab("pqrStuvW") -= L3aaa("pqrtuv") * L1b("SW"); + L4aaab("pqrStuvW") -= L3aab("pqSuvW") * L1a("rt"); + L4aaab("pqrStuvW") += L3aab("pqStvW") * L1a("ru"); + L4aaab("pqrStuvW") -= L3aab("pqStuW") * L1a("rv"); + L4aaab("pqrStuvW") += L3aab("prSuvW") * L1a("qt"); + L4aaab("pqrStuvW") -= L3aab("prStvW") * L1a("qu"); + L4aaab("pqrStuvW") += L3aab("prStuW") * L1a("qv"); + L4aaab("pqrStuvW") -= L3aab("qrSuvW") * L1a("pt"); + L4aaab("pqrStuvW") += L3aab("qrStvW") * L1a("pu"); + L4aaab("pqrStuvW") -= L3aab("qrStuW") * L1a("pv"); + + return L4aaab; +} + +ambit::Tensor RDMs::make_cumulant_L4aabb(const ambit::Tensor& L1a, const ambit::Tensor& L1b, + const ambit::Tensor& L2aa, const ambit::Tensor& L2ab, + const ambit::Tensor& L2bb, const ambit::Tensor& L3aab, + const ambit::Tensor& L3abb, const ambit::Tensor& g4aabb) { + timer t("make_cumulant_L4aabb"); + + auto L4aabb = g4aabb.clone(); + /// 1111 contributions + L4aabb("pqRStuVW") -= L1a("pu") * L1a("qt") * L1b("RW") * L1b("SV"); + L4aabb("pqRStuVW") += L1a("pu") * L1a("qt") * L1b("RV") * L1b("SW"); + L4aabb("pqRStuVW") += L1a("pt") * L1a("qu") * L1b("RW") * L1b("SV"); + L4aabb("pqRStuVW") -= L1a("pt") * L1a("qu") * L1b("RV") * L1b("SW"); + + /// 211 contributions + L4aabb("pqRStuVW") += L2aa("pqtu") * L1b("RW") * L1b("SV"); + L4aabb("pqRStuVW") -= L2aa("pqtu") * L1b("RV") * L1b("SW"); + L4aabb("pqRStuVW") -= L2ab("pRuW") * L1a("qt") * L1b("SV"); + L4aabb("pqRStuVW") += L2ab("pRtW") * L1a("qu") * L1b("SV"); + L4aabb("pqRStuVW") += L2ab("pRuV") * L1a("qt") * L1b("SW"); + L4aabb("pqRStuVW") -= L2ab("pRtV") * L1a("qu") * L1b("SW"); + L4aabb("pqRStuVW") += L2ab("pSuW") * L1a("qt") * L1b("RV"); + L4aabb("pqRStuVW") -= L2ab("pStW") * L1a("qu") * L1b("RV"); + L4aabb("pqRStuVW") -= L2ab("pSuV") * L1a("qt") * L1b("RW"); + L4aabb("pqRStuVW") += L2ab("pStV") * L1a("qu") * L1b("RW"); + L4aabb("pqRStuVW") += L2ab("qRuW") * L1a("pt") * L1b("SV"); + L4aabb("pqRStuVW") -= L2ab("qRtW") * L1a("pu") * L1b("SV"); + L4aabb("pqRStuVW") -= L2ab("qRuV") * L1a("pt") * L1b("SW"); + L4aabb("pqRStuVW") += L2ab("qRtV") * L1a("pu") * L1b("SW"); + L4aabb("pqRStuVW") -= L2ab("qSuW") * L1a("pt") * L1b("RV"); + L4aabb("pqRStuVW") += L2ab("qStW") * L1a("pu") * L1b("RV"); + L4aabb("pqRStuVW") += L2ab("qSuV") * L1a("pt") * L1b("RW"); + L4aabb("pqRStuVW") -= L2ab("qStV") * L1a("pu") * L1b("RW"); + L4aabb("pqRStuVW") += L2bb("RSVW") * L1a("pu") * L1a("qt"); + L4aabb("pqRStuVW") -= L2bb("RSVW") * L1a("pt") * L1a("qu"); + + /// 22 contributions + L4aabb("pqRStuVW") -= L2aa("pqtu") * L2bb("RSVW"); + L4aabb("pqRStuVW") -= L2ab("pRuW") * L2ab("qStV"); + L4aabb("pqRStuVW") += L2ab("pRtW") * L2ab("qSuV"); + L4aabb("pqRStuVW") += L2ab("pRuV") * L2ab("qStW"); + L4aabb("pqRStuVW") -= L2ab("pRtV") * L2ab("qSuW"); + L4aabb("pqRStuVW") += L2ab("pSuW") * L2ab("qRtV"); + L4aabb("pqRStuVW") -= L2ab("pStW") * L2ab("qRuV"); + L4aabb("pqRStuVW") -= L2ab("pSuV") * L2ab("qRtW"); + L4aabb("pqRStuVW") += L2ab("pStV") * L2ab("qRuW"); + + /// 31 contributions + L4aabb("pqRStuVW") += L3aab("pqRtuW") * L1b("SV"); + L4aabb("pqRStuVW") -= L3aab("pqRtuV") * L1b("SW"); + L4aabb("pqRStuVW") -= L3aab("pqStuW") * L1b("RV"); + L4aabb("pqRStuVW") += L3aab("pqStuV") * L1b("RW"); + L4aabb("pqRStuVW") += L3abb("pRSuVW") * L1a("qt"); + L4aabb("pqRStuVW") -= L3abb("pRStVW") * L1a("qu"); + L4aabb("pqRStuVW") -= L3abb("qRSuVW") * L1a("pt"); + L4aabb("pqRStuVW") += L3abb("qRStVW") * L1a("pu"); + + return L4aabb; +} + +ambit::Tensor RDMs::make_cumulant_L4abbb(const ambit::Tensor& L1a, const ambit::Tensor& L1b, + const ambit::Tensor& L2ab, const ambit::Tensor& L2bb, + const ambit::Tensor& L3abb, const ambit::Tensor& L3bbb, + const ambit::Tensor& g4abbb) { + timer t("make_cumulant_L4abbb"); + + auto L4abbb = g4abbb.clone(); + /// 1111 contributions (24 terms) + L4abbb("pQRStUVW") += L1a("pt") * L1b("QW") * L1b("RV") * L1b("SU"); + L4abbb("pQRStUVW") -= L1a("pt") * L1b("QW") * L1b("RU") * L1b("SV"); + L4abbb("pQRStUVW") -= L1a("pt") * L1b("QV") * L1b("RW") * L1b("SU"); + L4abbb("pQRStUVW") += L1a("pt") * L1b("QV") * L1b("RU") * L1b("SW"); + L4abbb("pQRStUVW") += L1a("pt") * L1b("QU") * L1b("RW") * L1b("SV"); + L4abbb("pQRStUVW") -= L1a("pt") * L1b("QU") * L1b("RV") * L1b("SW"); + + /// 211 contributions (72 terms) + L4abbb("pQRStUVW") += L2ab("pQtW") * L1b("RV") * L1b("SU"); + L4abbb("pQRStUVW") -= L2ab("pQtW") * L1a("RU") * L1b("SV"); + L4abbb("pQRStUVW") -= L2ab("pQtV") * L1b("RW") * L1b("SU"); + L4abbb("pQRStUVW") += L2ab("pQtV") * L1b("RU") * L1b("SW"); + L4abbb("pQRStUVW") += L2ab("pQtU") * L1b("RW") * L1b("SV"); + L4abbb("pQRStUVW") -= L2ab("pQtU") * L1b("RV") * L1b("SW"); + L4abbb("pQRStUVW") -= L2ab("pRtW") * L1b("QV") * L1b("SU"); + L4abbb("pQRStUVW") += L2ab("pRtW") * L1b("QU") * L1b("SV"); + L4abbb("pQRStUVW") += L2ab("pRtV") * L1b("QW") * L1b("SU"); + L4abbb("pQRStUVW") -= L2ab("pRtV") * L1b("QU") * L1b("SW"); + L4abbb("pQRStUVW") -= L2ab("pRtU") * L1b("QW") * L1b("SV"); + L4abbb("pQRStUVW") += L2ab("pRtU") * L1b("QV") * L1b("SW"); + L4abbb("pQRStUVW") += L2ab("pStW") * L1b("QV") * L1b("RU"); + L4abbb("pQRStUVW") -= L2ab("pStW") * L1b("QU") * L1b("RV"); + L4abbb("pQRStUVW") -= L2ab("pStV") * L1b("QW") * L1b("RU"); + L4abbb("pQRStUVW") += L2ab("pStV") * L1b("QU") * L1b("RW"); + L4abbb("pQRStUVW") += L2ab("pStU") * L1b("QW") * L1b("RV"); + L4abbb("pQRStUVW") -= L2ab("pStU") * L1b("QV") * L1b("RW"); + + L4abbb("pQRStUVW") -= L2bb("QRVW") * L1a("pt") * L1b("SU"); + L4abbb("pQRStUVW") += L2bb("QRUW") * L1a("pt") * L1b("SV"); + L4abbb("pQRStUVW") -= L2bb("QRUV") * L1a("pt") * L1b("SW"); + L4abbb("pQRStUVW") += L2bb("QSVW") * L1a("pt") * L1b("RU"); + L4abbb("pQRStUVW") -= L2bb("QSUW") * L1a("pt") * L1b("RV"); + L4abbb("pQRStUVW") += L2bb("QSUV") * L1a("pt") * L1b("RW"); + L4abbb("pQRStUVW") -= L2bb("RSVW") * L1a("pt") * L1b("QU"); + L4abbb("pQRStUVW") += L2bb("RSUW") * L1a("pt") * L1b("QV"); + L4abbb("pQRStUVW") -= L2bb("RSUV") * L1a("pt") * L1b("QW"); + + /// 22 contributions (18 terms) + L4abbb("pQRStUVW") -= L2ab("pQtW") * L2bb("RSUV"); + L4abbb("pQRStUVW") += L2ab("pQtV") * L2bb("RSUW"); + L4abbb("pQRStUVW") -= L2ab("pQtU") * L2bb("RSVW"); + L4abbb("pQRStUVW") += L2ab("pRtW") * L2bb("QSUV"); + L4abbb("pQRStUVW") -= L2ab("pRtV") * L2bb("QSUW"); + L4abbb("pQRStUVW") += L2ab("pRtU") * L2bb("QSVW"); + L4abbb("pQRStUVW") -= L2ab("pStW") * L2bb("QRUV"); + L4abbb("pQRStUVW") += L2ab("pStV") * L2bb("QRUW"); + L4abbb("pQRStUVW") -= L2ab("pStU") * L2bb("QRVW"); + + /// 31 contributions (16 terms) + L4abbb("pQRStUVW") -= L3abb("pQRtVW") * L1b("SU"); + L4abbb("pQRStUVW") += L3abb("pQRtUW") * L1b("SV"); + L4abbb("pQRStUVW") -= L3abb("pQRtUV") * L1b("SW"); + L4abbb("pQRStUVW") += L3abb("pQStVW") * L1b("RU"); + L4abbb("pQRStUVW") -= L3abb("pQStUW") * L1b("RV"); + L4abbb("pQRStUVW") += L3abb("pQStUV") * L1b("RW"); + L4abbb("pQRStUVW") -= L3abb("pRStVW") * L1b("QU"); + L4abbb("pQRStUVW") += L3abb("pRStUW") * L1b("QV"); + L4abbb("pQRStUVW") -= L3abb("pRStUV") * L1b("QW"); + + L4abbb("pQRStUVW") -= L3bbb("QRSUVW") * L1a("pt"); + + return L4abbb; +} + ambit::Tensor RDMs::sf1_to_sd1(const ambit::Tensor& G1) { auto g1 = G1.clone(); g1.scale(0.5); @@ -456,6 +835,34 @@ RDMsSpinDependent::RDMsSpinDependent(ambit::Tensor g1a, ambit::Tensor g1b, ambit _test_rdm_dims(g3bbb, "g3bbb", 6); } +RDMsSpinDependent::RDMsSpinDependent(ambit::Tensor g1a, ambit::Tensor g1b, ambit::Tensor g2aa, + ambit::Tensor g2ab, ambit::Tensor g2bb, ambit::Tensor g3aaa, + ambit::Tensor g3aab, ambit::Tensor g3abb, ambit::Tensor g3bbb, + ambit::Tensor g4aaaa, ambit::Tensor g4aaab, + ambit::Tensor g4aabb, ambit::Tensor g4abbb, + ambit::Tensor g4bbbb) + : g1a_(g1a), g1b_(g1b), g2aa_(g2aa), g2ab_(g2ab), g2bb_(g2bb), g3aaa_(g3aaa), g3aab_(g3aab), + g3abb_(g3abb), g3bbb_(g3bbb), g4aaaa_(g4aaaa), g4aaab_(g4aaab), g4aabb_(g4aabb), + g4abbb_(g4abbb), g4bbbb_(g4bbbb) { + max_rdm_ = 4; + type_ = RDMsType::spin_dependent; + n_orbs_ = g1a.dim(0); + _test_rdm_dims(g1a, "g1a", 2); + _test_rdm_dims(g1b, "g1b", 2); + _test_rdm_dims(g2aa, "g2aa", 4); + _test_rdm_dims(g2ab, "g2ab", 4); + _test_rdm_dims(g2bb, "g2bb", 4); + _test_rdm_dims(g3aaa, "g3aaa", 6); + _test_rdm_dims(g3aab, "g3aab", 6); + _test_rdm_dims(g3abb, "g3abb", 6); + _test_rdm_dims(g3bbb, "g3bbb", 6); + _test_rdm_dims(g4aaaa, "g4aaaa", 8); + _test_rdm_dims(g4aaab, "g4aaab", 8); + _test_rdm_dims(g4aabb, "g4aabb", 8); + _test_rdm_dims(g4abbb, "g4abbb", 8); + _test_rdm_dims(g4bbbb, "g4bbbb", 8); +} + ambit::Tensor RDMsSpinDependent::g1a() const { _test_rdm_level(1, "g1a"); return g1a_; @@ -492,6 +899,26 @@ ambit::Tensor RDMsSpinDependent::g3bbb() const { _test_rdm_level(3, "g3bbb"); return g3bbb_; } +ambit::Tensor RDMsSpinDependent::g4aaaa() const { + _test_rdm_level(4, "g4aaaa"); + return g4aaaa_; +} +ambit::Tensor RDMsSpinDependent::g4aaab() const { + _test_rdm_level(4, "g4aaab"); + return g4aaab_; +} +ambit::Tensor RDMsSpinDependent::g4aabb() const { + _test_rdm_level(4, "g4aabb"); + return g4aabb_; +} +ambit::Tensor RDMsSpinDependent::g4abbb() const { + _test_rdm_level(4, "g4abbb"); + return g4abbb_; +} +ambit::Tensor RDMsSpinDependent::g4bbbb() const { + _test_rdm_level(4, "g4bbbb"); + return g4bbbb_; +} ambit::Tensor RDMsSpinDependent::SF_G1() const { _test_rdm_level(1, "SF_G1"); auto G1 = g1a_.clone(); @@ -577,8 +1004,69 @@ ambit::Tensor RDMsSpinDependent::L3bbb() const { return L3bbb; } +ambit::Tensor RDMsSpinDependent::L4aaaa() const { + _test_rdm_level(4, "L4aaaa"); + auto _L1a = L1a(); + auto _L2aa = L2aa(); + auto _L3aaa = L3aaa(); + auto L4aaaa = make_cumulant_L4aaaa(_L1a, _L2aa, _L3aaa, g4aaaa_); + L4aaaa.set_name("L4aaaa"); + return L4aaaa; +} + +ambit::Tensor RDMsSpinDependent::L4aaab() const { + _test_rdm_level(4, "L4aaab"); + auto _L1a = L1a(); + auto _L1b = L1b(); + auto _L2aa = L2aa(); + auto _L2ab = L2ab(); + auto _L3aaa = L3aaa(); + auto _L3aab = L3aab(); + auto L4aaab = make_cumulant_L4aaab(_L1a, _L1b, _L2aa, _L2ab, _L3aaa, _L3aab, g4aaab_); + L4aaab.set_name("L4aaab"); + return L4aaab; +} + +ambit::Tensor RDMsSpinDependent::L4aabb() const { + _test_rdm_level(4, "L4aabb"); + auto _L1a = L1a(); + auto _L1b = L1b(); + auto _L2aa = L2aa(); + auto _L2ab = L2ab(); + auto _L2bb = L2bb(); + auto _L3aab = L3aab(); + auto _L3abb = L3abb(); + auto L4aabb = make_cumulant_L4aabb(_L1a, _L1b, _L2aa, _L2ab, _L2bb, _L3aab, _L3abb, g4aabb_); + L4aabb.set_name("L4aabb"); + return L4aabb; +} + +ambit::Tensor RDMsSpinDependent::L4abbb() const { + _test_rdm_level(4, "L4abbb"); + auto _L1a = L1a(); + auto _L1b = L1b(); + auto _L2ab = L2ab(); + auto _L2bb = L2bb(); + auto _L3abb = L3abb(); + auto _L3bbb = L3bbb(); + auto L4abbb = make_cumulant_L4abbb(_L1a, _L1b, _L2ab, _L2bb, _L3abb, _L3bbb, g4abbb_); + L4abbb.set_name("L4abbb"); + return L4abbb; +} + +ambit::Tensor RDMsSpinDependent::L4bbbb() const { + _test_rdm_level(4, "L4bbbb"); + auto _L1b = L1b(); + auto _L2bb = L2bb(); + auto _L3bbb = L3bbb(); + auto L4bbbb = make_cumulant_L4aaaa(_L1b, _L2bb, _L3bbb, g4bbbb_); + L4bbbb.set_name("L4bbbb"); + return L4bbbb; +} + std::shared_ptr RDMsSpinDependent::clone() { - ambit::Tensor g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb; + ambit::Tensor g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb, g4aaaa, g4aaab, g4aabb, + g4abbb, g4bbbb; if (max_rdm_ > 0) { g1a = g1a_.clone(); g1b = g1b_.clone(); @@ -594,6 +1082,13 @@ std::shared_ptr RDMsSpinDependent::clone() { g3abb = g3abb_.clone(); g3bbb = g3bbb_.clone(); } + if (max_rdm_ > 3) { + g4aaaa = g4aaaa_.clone(); + g4aaab = g4aaab_.clone(); + g4aabb = g4aabb_.clone(); + g4abbb = g4abbb_.clone(); + g4bbbb = g4bbbb_.clone(); + } std::shared_ptr rdms; @@ -603,9 +1098,12 @@ std::shared_ptr RDMsSpinDependent::clone() { rdms = std::make_shared(g1a, g1b); else if (max_rdm_ == 2) rdms = std::make_shared(g1a, g1b, g2aa, g2ab, g2bb); - else + else if (max_rdm_ == 3) rdms = std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb); + else + rdms = std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, + g3bbb, g4aaaa, g4aaab, g4aabb, g4abbb, g4bbbb); return rdms; } @@ -626,6 +1124,13 @@ void RDMsSpinDependent::scale(double factor) { g3abb_.scale(factor); g3bbb_.scale(factor); } + if (max_rdm_ > 3) { + g4aaaa_.scale(factor); + g4aaab_.scale(factor); + g4aabb_.scale(factor); + g4abbb_.scale(factor); + g4bbbb_.scale(factor); + } } void RDMsSpinDependent::axpy(std::shared_ptr rhs, double a) { @@ -651,6 +1156,13 @@ void RDMsSpinDependent::axpy(std::shared_ptr rhs, double a) { g3abb_("pqrstu") += a * rhs->g3abb()("pqrstu"); g3bbb_("pqrstu") += a * rhs->g3bbb()("pqrstu"); } + if (max_rdm_ > 3) { + g4aaaa_("pqrstuvw") += a * rhs->g4aaaa()("pqrstuvw"); + g4aaab_("pqrstuvw") += a * rhs->g4aaab()("pqrstuvw"); + g4aabb_("pqrstuvw") += a * rhs->g4aabb()("pqrstuvw"); + g4abbb_("pqrstuvw") += a * rhs->g4abbb()("pqrstuvw"); + g4bbbb_("pqrstuvw") += a * rhs->g4bbbb()("pqrstuvw"); + } } void RDMsSpinDependent::rotate(const ambit::Tensor& Ua, const ambit::Tensor& Ub) { @@ -712,6 +1224,34 @@ void RDMsSpinDependent::rotate(const ambit::Tensor& Ua, const ambit::Tensor& Ub) g3bbb_("pqrstu") = g3T("pqrstu"); psi::outfile->Printf("\n Transformed 3 RDMs."); + + if (max_rdm_ == 3) + return; + + // Transform the 4-rdms + auto g4T = ambit::Tensor::build(ambit::CoreTensor, "g4T", g4aaaa_.dims()); + g4T("pqrstuvw") = Ua("ap") * Ua("bq") * Ua("cr") * Ua("ds") * g4aaaa_("abcdijkl") * Ua("it") * + Ua("ju") * Ua("kv") * Ua("lw"); + g4aaaa_("pqrstuvw") = g4T("pqrstuvw"); + + g4T("pqrStuvW") = Ua("ap") * Ua("bq") * Ua("cr") * Ub("DS") * g4aaab_("abcDijkL") * Ua("it") * + Ua("ju") * Ua("kv") * Ub("LW"); + g4aaab_("pqrStuvW") = g4T("pqrStuvW"); + g4T("pqRStuVW") = Ua("ap") * Ua("bq") * Ub("CR") * Ub("DS") * g4aabb_("abCDijKL") * Ua("it") * + Ua("ju") * Ub("KV") * Ub("LW"); + g4aabb_("pqRStuVW") = g4T("pqRStuVW"); + g4T("pQRStUVW") = Ua("ap") * Ub("BQ") * Ub("CR") * Ub("DS") * g4abbb_("aBCDiJKL") * Ua("it") * + Ub("JU") * Ub("KV") * Ub("LW"); + g4abbb_("pQRStUVW") = g4T("pQRStUVW"); + + g4T("PQRSTUVW") = Ub("AP") * Ub("BQ") * Ub("CR") * Ub("DS") * g4bbbb_("ABCDIJKL") * Ub("IT") * + Ub("JU") * Ub("KV") * Ub("LW"); + g4bbbb_("PQRSTUVW") = g4T("PQRSTUVW"); + + psi::outfile->Printf("\n Transformed 4 RDMs."); + + if (max_rdm_ == 4) + return; } void RDMsSpinDependent::dump_to_disk(const std::string& filename_prefix) const { @@ -732,6 +1272,13 @@ void RDMsSpinDependent::dump_to_disk(const std::string& filename_prefix) const { ambit::save(g3abb_, prefix + "g3abb.bin"); ambit::save(g3bbb_, prefix + "g3bbb.bin"); } + if (max_rdm_ > 3) { + ambit::save(g4aaaa_, prefix + "g4aaaa.bin"); + ambit::save(g4aaab_, prefix + "g4aaab.bin"); + ambit::save(g4aabb_, prefix + "g4aabb.bin"); + ambit::save(g4abbb_, prefix + "g4abbb.bin"); + ambit::save(g4bbbb_, prefix + "g4bbbb.bin"); + } } RDMsSpinFree::RDMsSpinFree() { @@ -831,6 +1378,31 @@ ambit::Tensor RDMsSpinFree::g3bbb() const { g3bbb.set_name("g3bbb"); return g3bbb; } +ambit::Tensor RDMsSpinFree::g4aaaa() const { + throw std::runtime_error("RDMsSpinFree::g4aaaa not implemented."); + ambit::Tensor g4aaaa; + return g4aaaa; +} +ambit::Tensor RDMsSpinFree::g4aaab() const { + throw std::runtime_error("RDMsSpinFree::g4aaab not implemented."); + ambit::Tensor g4aaab; + return g4aaab; +} +ambit::Tensor RDMsSpinFree::g4aabb() const { + throw std::runtime_error("RDMsSpinFree::g4aabb not implemented."); + ambit::Tensor g4aabb; + return g4aabb; +} +ambit::Tensor RDMsSpinFree::g4abbb() const { + throw std::runtime_error("RDMsSpinFree::g4abbb not implemented."); + ambit::Tensor g4abbb; + return g4abbb; +} +ambit::Tensor RDMsSpinFree::g4bbbb() const { + throw std::runtime_error("RDMsSpinFree::g4bbbb not implemented."); + ambit::Tensor g4bbbb; + return g4bbbb; +} ambit::Tensor RDMsSpinFree::L1a() const { _test_rdm_level(1, "L1a"); auto L1a = sf1_to_sd1(SF_G1_); @@ -886,6 +1458,36 @@ ambit::Tensor RDMsSpinFree::L3bbb() const { return L3bbb; } +ambit::Tensor RDMsSpinFree::L4aaaa() const { + throw std::runtime_error("RDMsSpinFree::L4aaaa not implemented."); + ambit::Tensor L4aaaa; + return L4aaaa; +} + +ambit::Tensor RDMsSpinFree::L4aaab() const { + throw std::runtime_error("RDMsSpinFree::L4aaab not implemented."); + ambit::Tensor L4aaab; + return L4aaab; +} + +ambit::Tensor RDMsSpinFree::L4aabb() const { + throw std::runtime_error("RDMsSpinFree::L4aabb not implemented."); + ambit::Tensor L4aabb; + return L4aabb; +} + +ambit::Tensor RDMsSpinFree::L4abbb() const { + throw std::runtime_error("RDMsSpinFree::L4abbb not implemented."); + ambit::Tensor L4abbb; + return L4abbb; +} + +ambit::Tensor RDMsSpinFree::L4bbbb() const { + throw std::runtime_error("RDMsSpinFree::L4bbbb not implemented."); + ambit::Tensor L4bbbb; + return L4bbbb; +} + std::shared_ptr RDMsSpinFree::clone() { ambit::Tensor g1, g2, g3; if (max_rdm_ > 0) diff --git a/forte/base_classes/rdms.h b/forte/base_classes/rdms.h index 36cd1c31d..942bb4d00 100644 --- a/forte/base_classes/rdms.h +++ b/forte/base_classes/rdms.h @@ -170,6 +170,28 @@ class RDMs { static ambit::Tensor make_cumulant_L3abb(const ambit::Tensor& g1a, const ambit::Tensor& g1b, const ambit::Tensor& g2ab, const ambit::Tensor& g2bb, const ambit::Tensor& g3abb); + /// Make alpha-alpha-alpha-alpha or beta-beta-beta-beta 4-RDC from 4-RDMs + static ambit::Tensor make_cumulant_L4aaaa(const ambit::Tensor& L1a, const ambit::Tensor& L2aa, + const ambit::Tensor& L3aaa, + const ambit::Tensor& g4aaaa); + /// Make alpha-alpha-alpha-beta 4-RDC from 4-RDMs + static ambit::Tensor make_cumulant_L4aaab(const ambit::Tensor& L1a, const ambit::Tensor& L1b, + const ambit::Tensor& L2aa, const ambit::Tensor& L2ab, + const ambit::Tensor& L3aaa, + const ambit::Tensor& L3aab, + const ambit::Tensor& g4aaab); + /// Make alpha-alpha-beta-beta 4-RDC from 4-RDMs + static ambit::Tensor make_cumulant_L4aabb(const ambit::Tensor& L1a, const ambit::Tensor& L1b, + const ambit::Tensor& L2aa, const ambit::Tensor& L2ab, + const ambit::Tensor& L2bb, const ambit::Tensor& L3aab, + const ambit::Tensor& L3abb, + const ambit::Tensor& g4aabb); + /// Make alpha-beta-beta-beta 4-RDC from 4-RDMs + static ambit::Tensor make_cumulant_L4abbb(const ambit::Tensor& L1a, const ambit::Tensor& L1b, + const ambit::Tensor& L2ab, const ambit::Tensor& L2bb, + const ambit::Tensor& L3abb, + const ambit::Tensor& L3bbb, + const ambit::Tensor& g4abbb); /// Spin-free 1-body to spin-dependent 1-body subject to Ms averaging static ambit::Tensor sf1_to_sd1(const ambit::Tensor& G1); @@ -201,6 +223,16 @@ class RDMs { virtual ambit::Tensor g3abb() const = 0; /// @return the beta-beta-beta 3-RDM virtual ambit::Tensor g3bbb() const = 0; + /// @return the alpha-alpha-alpha-alpha 4-RDM + virtual ambit::Tensor g4aaaa() const = 0; + /// @return the alpha-alpha-alpha-beta 4-RDM + virtual ambit::Tensor g4aaab() const = 0; + /// @return the alpha-alpha-beta-beta 4-RDM + virtual ambit::Tensor g4aabb() const = 0; + /// @return the alpha-beta-beta-beta 4-RDM + virtual ambit::Tensor g4abbb() const = 0; + /// @return the beta-beta-beta-beta 4-RDM + virtual ambit::Tensor g4bbbb() const = 0; // Spin-free RDMs @@ -236,6 +268,16 @@ class RDMs { virtual ambit::Tensor L3abb() const = 0; /// @return the beta-beta-beta 3-RDC virtual ambit::Tensor L3bbb() const = 0; + /// @return the alpha-alpha-alpha-alpha 4-RDC + virtual ambit::Tensor L4aaaa() const = 0; + /// @return the alpha-alpha-alpha-beta 4-RDC + virtual ambit::Tensor L4aaab() const = 0; + /// @return the alpha-alpha-beta-beta 4-RDC + virtual ambit::Tensor L4aabb() const = 0; + /// @return the alpha-beta-beta-beta 4-RDC + virtual ambit::Tensor L4abbb() const = 0; + /// @return the beta-beta-beta-beta 4-RDC + virtual ambit::Tensor L4bbbb() const = 0; // Spin-free density cumulants @@ -281,6 +323,12 @@ class RDMsSpinDependent : public RDMs { RDMsSpinDependent(ambit::Tensor g1a, ambit::Tensor g1b, ambit::Tensor g2aa, ambit::Tensor g2ab, ambit::Tensor g2bb, ambit::Tensor g3aaa, ambit::Tensor g3aab, ambit::Tensor g3abb, ambit::Tensor g3bbb); + /// @brief Construct a RDMsSpinDependent object with the 1-, 2-, 3-, and 4-rdms + RDMsSpinDependent(ambit::Tensor g1a, ambit::Tensor g1b, ambit::Tensor g2aa, ambit::Tensor g2ab, + ambit::Tensor g2bb, ambit::Tensor g3aaa, ambit::Tensor g3aab, + ambit::Tensor g3abb, ambit::Tensor g3bbb, ambit::Tensor g4aaaa, + ambit::Tensor g4aaab, ambit::Tensor g4aabb, ambit::Tensor g4abbb, + ambit::Tensor g4bbbb); /// @return the alpha 1-RDM ambit::Tensor g1a() const override; @@ -300,6 +348,16 @@ class RDMsSpinDependent : public RDMs { ambit::Tensor g3abb() const override; /// @return the beta-beta-beta 3-RDM ambit::Tensor g3bbb() const override; + /// @return the alpha-alpha-alpha-alpha 4-RDM + ambit::Tensor g4aaaa() const override; + /// @return the alpha-alpha-alpha-beta 4-RDM + ambit::Tensor g4aaab() const override; + /// @return the alpha-alpha-beta-beta 4-RDM + ambit::Tensor g4aabb() const override; + /// @return the alpha-beta-beta-beta 4-RDM + ambit::Tensor g4abbb() const override; + /// @return the beta-beta-beta-beta 4-RDM + ambit::Tensor g4bbbb() const override; // Spin-free RDMs @@ -329,6 +387,16 @@ class RDMsSpinDependent : public RDMs { ambit::Tensor L3abb() const override; /// @return the beta-beta-beta 3-RDC ambit::Tensor L3bbb() const override; + /// @return the alpha-alpha-alpha-alpha 4-RDC + ambit::Tensor L4aaaa() const override; + /// @return the alpha-alpha-alpha-beta 4-RDC + ambit::Tensor L4aaab() const override; + /// @return the alpha-alpha-beta-beta 4-RDC + ambit::Tensor L4aabb() const override; + /// @return the alpha-beta-beta-beta 4-RDC + ambit::Tensor L4abbb() const override; + /// @return the beta-beta-beta-beta 4-RDC + ambit::Tensor L4bbbb() const override; // class methods @@ -366,6 +434,16 @@ class RDMsSpinDependent : public RDMs { ambit::Tensor g3abb_; /// The beta-beta-beta 3-RDM ambit::Tensor g3bbb_; + /// The alpha-alpha-alpha-alpha 4-RDM + ambit::Tensor g4aaaa_; + /// The alpha-alpha-alpha-beta 4-RDM + ambit::Tensor g4aaab_; + /// The alpha-alpha-beta-beta 4-RDM + ambit::Tensor g4aabb_; + /// The alpha-beta-beta-beta 4-RDM + ambit::Tensor g4abbb_; + /// The beta-beta-beta-beta 4-RDM + ambit::Tensor g4bbbb_; }; class RDMsSpinFree : public RDMs { @@ -397,6 +475,16 @@ class RDMsSpinFree : public RDMs { ambit::Tensor g3abb() const override; /// @return the beta-beta-beta 3-RDM ambit::Tensor g3bbb() const override; + /// @return the alpha-alpha-alpha-alpha 4-RDM + ambit::Tensor g4aaaa() const override; + /// @return the alpha-alpha-alpha-beta 4-RDM + ambit::Tensor g4aaab() const override; + /// @return the alpha-alpha-beta-beta 4-RDM + ambit::Tensor g4aabb() const override; + /// @return the alpha-beta-beta-beta 4-RDM + ambit::Tensor g4abbb() const override; + /// @return the beta-beta-beta-beta 4-RDM + ambit::Tensor g4bbbb() const override; // Spin-free RDMs @@ -426,6 +514,16 @@ class RDMsSpinFree : public RDMs { ambit::Tensor L3abb() const override; /// @return the beta-beta-beta 3-RDC ambit::Tensor L3bbb() const override; + /// @return the alpha-alpha-alpha-alpha 4-RDC + ambit::Tensor L4aaaa() const override; + /// @return the alpha-alpha-alpha-beta 4-RDC + ambit::Tensor L4aaab() const override; + /// @return the alpha-alpha-beta-beta 4-RDC + ambit::Tensor L4aabb() const override; + /// @return the alpha-beta-beta-beta 4-RDC + ambit::Tensor L4abbb() const override; + /// @return the beta-beta-beta-beta 4-RDC + ambit::Tensor L4bbbb() const override; // class methods diff --git a/forte/fci/string_list_defs.h b/forte/fci/string_list_defs.h index 2f71a60e0..535820849 100644 --- a/forte/fci/string_list_defs.h +++ b/forte/fci/string_list_defs.h @@ -74,6 +74,18 @@ struct H3StringSubstitution { : sign(sign_), p(p_), q(q_), r(r_), J(J_) {} }; +/// 4-hole string substitution +struct H4StringSubstitution { + const int16_t sign; + const int16_t p; + const int16_t q; + const int16_t r; + const int16_t s; + const size_t J; + H4StringSubstitution(int16_t sign_, int16_t p_, int16_t q_, int16_t r_, int16_t s_, size_t J_) + : sign(sign_), p(p_), q(q_), r(r_), s(s_), J(J_) {} +}; + using StringList = std::vector>; /// Maps the integers (p,q,h) to list of strings connected by a^{+}_p a_q, where the string @@ -111,6 +123,10 @@ using H2List = std::map, std::vector, std::vector>; +/// Maps the integers (h_J, add_J, h_I) to list of strings connected by a_p a_q a_r a_s, where the +/// string I belongs to the irrep h_I and J belongs to the irrep h_J and add_J is the address of J +using H4List = std::map, std::vector>; + using Pair = std::pair; using PairList = std::vector>>; diff --git a/forte/genci/genci_solver.cc b/forte/genci/genci_solver.cc index 06619b580..d9ea35d8b 100644 --- a/forte/genci/genci_solver.cc +++ b/forte/genci/genci_solver.cc @@ -377,7 +377,7 @@ void GenCISolver::test_rdms(std::shared_ptr b, std::shared_ptrcompute_rdms(*C_, *C_, max_rdm_level, RDMsType::spin_dependent); C_->test_rdms(*C_, *C_, max_rdm_level, RDMsType::spin_dependent, rdms); } diff --git a/forte/genci/genci_string_lists.cc b/forte/genci/genci_string_lists.cc index f285d940f..2f8995814 100644 --- a/forte/genci/genci_string_lists.cc +++ b/forte/genci/genci_string_lists.cc @@ -124,6 +124,7 @@ void GenCIStringLists::startup(std::shared_ptr mo_space_info) { double h1_list_timer = 0.0; double h2_list_timer = 0.0; double h3_list_timer = 0.0; + double h4_list_timer = 0.0; double vovo_list_timer = 0.0; double vvoo_list_timer = 0.0; @@ -154,6 +155,8 @@ void GenCIStringLists::startup(std::shared_ptr mo_space_info) { gas_beta_2h_occupations_ = generate_1h_occupations(gas_beta_1h_occupations_); gas_alfa_3h_occupations_ = generate_1h_occupations(gas_alfa_2h_occupations_); gas_beta_3h_occupations_ = generate_1h_occupations(gas_beta_2h_occupations_); + gas_alfa_4h_occupations_ = generate_1h_occupations(gas_alfa_3h_occupations_); + gas_beta_4h_occupations_ = generate_1h_occupations(gas_beta_3h_occupations_); if (na_ >= 1) { auto alfa_1h_strings = make_strings_with_occupation( @@ -186,6 +189,16 @@ void GenCIStringLists::startup(std::shared_ptr mo_space_info) { ngas_spaces_, nirrep_, gas_size_, gas_mos_, gas_beta_3h_occupations_, string_class_); beta_address_3h_ = std::make_shared(gas_size_, nb_ - 3, beta_3h_strings); } + if (na_ >= 4) { + auto alfa_4h_strings = make_strings_with_occupation( + ngas_spaces_, nirrep_, gas_size_, gas_mos_, gas_alfa_4h_occupations_, string_class_); + alfa_address_4h_ = std::make_shared(gas_size_, na_ - 4, alfa_4h_strings); + } + if (nb_ >= 4) { + auto beta_4h_strings = make_strings_with_occupation( + ngas_spaces_, nirrep_, gas_size_, gas_mos_, gas_beta_4h_occupations_, string_class_); + beta_address_4h_ = std::make_shared(gas_size_, nb_ - 4, beta_4h_strings); + } nas_ = 0; nbs_ = 0; @@ -250,6 +263,12 @@ void GenCIStringLists::startup(std::shared_ptr mo_space_info) { beta_3h_list = make_3h_list(beta_strings_, beta_address_, beta_address_3h_); h3_list_timer += t.get(); } + { + local_timer t; + alfa_4h_list = make_4h_list(alfa_strings_, alfa_address_, alfa_address_4h_); + beta_4h_list = make_4h_list(beta_strings_, beta_address_, beta_address_4h_); + h4_list_timer += t.get(); + } double total_time = str_list_timer + nn_list_timer + vo_list_timer + oo_list_timer + vvoo_list_timer + vovo_list_timer; @@ -269,6 +288,7 @@ void GenCIStringLists::startup(std::shared_ptr mo_space_info) { {"timing for 1-hole strings", h1_list_timer}, {"timing for 2-hole strings", h2_list_timer}, {"timing for 3-hole strings", h3_list_timer}, + {"timing for 4-hole strings", h4_list_timer}, {"total timing", total_time}}); } std::string table = printer.get_table("String Lists"); @@ -423,4 +443,16 @@ std::vector& GenCIStringLists::get_beta_3h_list(int h_I, s return beta_3h_list[I_tuple]; } +std::vector& GenCIStringLists::get_alfa_4h_list(int h_I, size_t add_I, + int h_J) { + std::tuple I_tuple(h_I, add_I, h_J); + return alfa_4h_list[I_tuple]; +} + +std::vector& GenCIStringLists::get_beta_4h_list(int h_I, size_t add_I, + int h_J) { + std::tuple I_tuple(h_I, add_I, h_J); + return beta_4h_list[I_tuple]; +} + } // namespace forte \ No newline at end of file diff --git a/forte/genci/genci_string_lists.h b/forte/genci/genci_string_lists.h index d6136686e..6d027b8e6 100644 --- a/forte/genci/genci_string_lists.h +++ b/forte/genci/genci_string_lists.h @@ -104,6 +104,10 @@ class GenCIStringLists { auto alfa_address_3h() { return alfa_address_3h_; } /// @return the beta string address object for N - 3 electrons auto beta_address_3h() { return beta_address_3h_; } + /// @return the alpha string address object for N - 4 electrons + auto alfa_address_4h() { return alfa_address_4h_; } + /// @return the beta string address object for N - 4 electrons + auto beta_address_4h() { return beta_address_4h_; } /// @return the address of a determinant in the CI vector size_t determinant_address(const Determinant& d) const; @@ -151,6 +155,9 @@ class GenCIStringLists { std::vector& get_alfa_3h_list(int h_I, size_t add_I, int h_J); std::vector& get_beta_3h_list(int h_I, size_t add_I, int h_J); + std::vector& get_alfa_4h_list(int h_I, size_t add_I, int h_J); + std::vector& get_beta_4h_list(int h_I, size_t add_I, int h_J); + Pair get_pair_list(int h, int n) const { return pair_list_[h][n]; } private: @@ -195,6 +202,8 @@ class GenCIStringLists { std::vector> gas_beta_2h_occupations_; std::vector> gas_alfa_3h_occupations_; std::vector> gas_beta_3h_occupations_; + std::vector> gas_alfa_4h_occupations_; + std::vector> gas_beta_4h_occupations_; std::vector> gas_occupations_; @@ -248,6 +257,9 @@ class GenCIStringLists { /// The 3-hole lists H3List alfa_3h_list; H3List beta_3h_list; + /// The 4-hole lists + H4List alfa_4h_list; + H4List beta_4h_list; /// Addressers /// The alpha string address @@ -266,6 +278,10 @@ class GenCIStringLists { std::shared_ptr alfa_address_3h_; /// The beta string address for N - 3 electrons std::shared_ptr beta_address_3h_; + /// The alpha string address for N - 4 electrons + std::shared_ptr alfa_address_4h_; + /// The beta string address for N - 4 electrons + std::shared_ptr beta_address_4h_; // ==> Class Functions <== diff --git a/forte/genci/genci_vector.h b/forte/genci/genci_vector.h index cfbe68117..6044e53ef 100644 --- a/forte/genci/genci_vector.h +++ b/forte/genci/genci_vector.h @@ -285,6 +285,16 @@ class GenCIVector { return (ncmo * ncmo * ncmo * ncmo * ncmo * p + ncmo * ncmo * ncmo * ncmo * q + ncmo * ncmo * ncmo * r + ncmo * ncmo * s + ncmo * t + u); } + static size_t eight_index(size_t p, size_t q, size_t r, size_t s, size_t t, size_t u, + size_t v, size_t w, size_t ncmo) { + size_t ncmo2 = ncmo * ncmo; + size_t ncmo4 = ncmo2 * ncmo2; + size_t ncmo6 = ncmo4 * ncmo2; + return (ncmo6 * ncmo * p + ncmo6 * q + + ncmo4 * ncmo * r + ncmo4 * s + + ncmo2 * ncmo * t + ncmo2 * u + + ncmo * v + w); + } /// @brief Apply the scalar part of the Hamiltonian to this vector and add it to the result /// @param result The wave function to add the result to @@ -329,7 +339,7 @@ class GenCIVector { // 3-RDM elements are stored in the format // -> rdm[six_index(p,q,r,s,t,u)] - /// Compute the matrix elements of the same spin 3-RDM (with all indices + /// Compute the matrix elements of the same spin 3-RDM (with all indices /// alpha or beta) static ambit::Tensor compute_3rdm_aaa_same_irrep(GenCIVector& C_left, GenCIVector& C_right, bool alfa); @@ -340,6 +350,23 @@ class GenCIVector { /// a_{tb} a_{sa}> static ambit::Tensor compute_3rdm_abb_same_irrep(GenCIVector& C_left, GenCIVector& C_right); + /// 4-RDM elements are stored in the format + /// -> rdm[eight_index(p,q,r,s,t,u,v,w)] + + /// Compute the matrix elements of the same spin 4-RDM + /// (with all indices alpha or beta) + static ambit::Tensor compute_4rdm_aaaa_same_irrep(GenCIVector& C_left, GenCIVector& C_right, + bool alfa); + /// Compute the matrix elements of the alpha-alpha-alpha-beta 4-RDM + static ambit::Tensor compute_4rdm_aaab_same_irrep(GenCIVector& C_left, GenCIVector& C_right); + /// Compute the matrix elements of the alpha-alpha-beta-beta 4-RDM + static ambit::Tensor compute_4rdm_aabb_same_irrep(GenCIVector& C_left, GenCIVector& C_right); + /// Compute the matrix elements of the alpha-beta-beta-beta 4-RDM + static ambit::Tensor compute_4rdm_abbb_same_irrep(GenCIVector& C_left, GenCIVector& C_right); + public: /// @brief Provide a pointer to the a block of the coefficient matrix in such a way that we can /// use its content in several algorithms (sigma vector, RDMs, etc.) diff --git a/forte/genci/genci_vector_rdm.cc b/forte/genci/genci_vector_rdm.cc index 20e355e99..f3defd291 100644 --- a/forte/genci/genci_vector_rdm.cc +++ b/forte/genci/genci_vector_rdm.cc @@ -51,6 +51,13 @@ std::shared_ptr GenCIVector::compute_rdms(GenCIVector& C_left, GenCIVector ambit::Tensor g1a, g1b; ambit::Tensor g2aa, g2ab, g2bb; ambit::Tensor g3aaa, g3aab, g3abb, g3bbb; + ambit::Tensor g4aaaa, g4aaab, g4aabb, g4abbb, g4bbbb; + + if (max_rdm_level >= 5) { + throw std::runtime_error( + "RDMs of order 5 or higher are not implemented in GenCISolver (and " + "more generally in Forte)."); + } if (max_rdm_level >= 1) { local_timer t; @@ -110,6 +117,56 @@ std::shared_ptr GenCIVector::compute_rdms(GenCIVector& C_left, GenCIVector rdm_timing.push_back(t.get()); } + if (max_rdm_level >= 4) { + local_timer t; + if (na >= 4) { + local_timer t_aaaa; + g4aaaa = compute_4rdm_aaaa_same_irrep(C_left, C_right, true); + psi::outfile->Printf("\n Timing for 4-RDM (aaaa): %.3f s", t_aaaa.get()); + } else { + g4aaaa = ambit::Tensor::build(ambit::CoreTensor, "g4aaaa", + {nmo, nmo, nmo, nmo, nmo, nmo, nmo, nmo}); + g4aaaa.zero(); + } + if (nb >= 4) { + local_timer t_bbbb; + g4bbbb = compute_4rdm_aaaa_same_irrep(C_left, C_right, false); + psi::outfile->Printf("\n Timing for 4-RDM (bbbb): %.3f s", t_bbbb.get()); + } else { + g4bbbb = ambit::Tensor::build(ambit::CoreTensor, "g4bbbb", + {nmo, nmo, nmo, nmo, nmo, nmo, nmo, nmo}); + g4bbbb.zero(); + } + if ((na >= 3) and (nb >= 1)) { + local_timer t_aaab; + g4aaab = compute_4rdm_aaab_same_irrep(C_left, C_right); + psi::outfile->Printf("\n Timing for 4-RDM (aaab): %.3f s", t_aaab.get()); + } else { + g4aaab = ambit::Tensor::build(ambit::CoreTensor, "g4aaab", + {nmo, nmo, nmo, nmo, nmo, nmo, nmo, nmo}); + g4aaab.zero(); + } + if ((na >= 2) and (nb >= 2)) { + local_timer t_aabb; + g4aabb = compute_4rdm_aabb_same_irrep(C_left, C_right); + psi::outfile->Printf("\n Timing for 4-RDM (aabb): %.3f s", t_aabb.get()); + } else { + g4aabb = ambit::Tensor::build(ambit::CoreTensor, "g4aabb", + {nmo, nmo, nmo, nmo, nmo, nmo, nmo, nmo}); + g4aabb.zero(); + } + if ((na >= 1) and (nb >= 3)) { + local_timer t_abbb; + g4abbb = compute_4rdm_abbb_same_irrep(C_left, C_right); + psi::outfile->Printf("\n Timing for 4-RDM (abbb): %.3f s", t_abbb.get()); + } else { + g4abbb = ambit::Tensor::build(ambit::CoreTensor, "g4abbb", + {nmo, nmo, nmo, nmo, nmo, nmo, nmo, nmo}); + g4abbb.zero(); + } + rdm_timing.push_back(t.get()); + } + // for (size_t n = 0; n < rdm_timing.size(); ++n) { // psi::outfile->Printf("\n Timing for %d-RDM: %.3f s", n + 1, rdm_timing[n]); // } @@ -125,7 +182,16 @@ std::shared_ptr GenCIVector::compute_rdms(GenCIVector& C_left, GenCIVector return std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, g3abb, g3bbb); } + if (max_rdm_level == 4) { + return std::make_shared(g1a, g1b, g2aa, g2ab, g2bb, g3aaa, g3aab, + g3abb, g3bbb, g4aaaa, g4aaab, g4aabb, g4abbb, + g4bbbb); + } } else { + if (max_rdm_level == 4) { + throw std::runtime_error( + "Spin-free RDMs of order 4 or higher are not implemented in GenCISolver"); + } g1a("pq") += g1b("pq"); if (max_rdm_level > 1) { @@ -144,12 +210,6 @@ std::shared_ptr GenCIVector::compute_rdms(GenCIVector& C_left, GenCIVector if (max_rdm_level == 3) return std::make_shared(g1a, g2aa, g3aaa); } - - if (max_rdm_level >= 4) { - throw std::runtime_error( - "RDMs of order 4 or higher are not implemented in GenCISolver (and " - "more generally in Forte)."); - } return std::make_shared(); } @@ -620,6 +680,475 @@ ambit::Tensor GenCIVector::compute_3rdm_abb_same_irrep(GenCIVector& C_left, GenC return rdm; } +ambit::Tensor GenCIVector::compute_4rdm_aaaa_same_irrep(GenCIVector& C_left, GenCIVector& C_right, + bool alfa) { + size_t ncmo = C_left.ncmo_; + const auto& alfa_address = C_left.alfa_address_; + const auto& beta_address = C_left.beta_address_; + const auto& lists = C_left.lists_; + + auto rdm = ambit::Tensor::build(ambit::CoreTensor, alfa ? "4RDM_AAAA" : "4RDM_BBBB", + {ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo}); + + auto na = alfa_address->nones(); + auto nb = beta_address->nones(); + if ((alfa and (na < 4)) or ((!alfa) and (nb < 4))) + return rdm; + + auto& rdm_data = rdm.data(); + + int num_4h_classes = + alfa ? lists->alfa_address_4h()->nclasses() : lists->beta_address_4h()->nclasses(); + + for (int class_K = 0; class_K < num_4h_classes; ++class_K) { + size_t maxK = alfa ? lists->alfa_address_4h()->strpcls(class_K) + : lists->beta_address_4h()->strpcls(class_K); + + // loop over blocks of matrix C + for (const auto& [nI, class_Ia, class_Ib] : lists->determinant_classes()) { + if (lists->detpblk(nI) == 0) + continue; + + auto Cr = C_right.gather_C_block(CR, alfa, alfa_address, beta_address, class_Ia, + class_Ib, false); + + for (const auto& [nJ, class_Ja, class_Jb] : lists->determinant_classes()) { + // The string class on which we don't act must be the same for I and J + if ((alfa and (class_Ib != class_Jb)) or (not alfa and (class_Ia != class_Ja))) + continue; + if (lists->detpblk(nJ) == 0) + continue; + + // Get a pointer to the correct block of matrix C + auto Cl = C_left.gather_C_block(CL, alfa, alfa_address, beta_address, class_Ja, + class_Jb, false); + + size_t maxL = + alfa ? beta_address->strpcls(class_Ib) : alfa_address->strpcls(class_Ia); + if (maxL > 0) { + for (size_t K = 0; K < maxK; ++K) { + std::vector& Krlist = + alfa ? lists->get_alfa_4h_list(class_K, K, class_Ia) + : lists->get_beta_4h_list(class_K, K, class_Ib); + std::vector& Kllist = + alfa ? lists->get_alfa_4h_list(class_K, K, class_Ja) + : lists->get_beta_4h_list(class_K, K, class_Jb); + for (const auto& [sign_K, p, q, r, s, I] : Krlist) { + for (const auto& [sign_L, t, u, v, w, J] : Kllist) { + rdm_data[eight_index(p, q, r, s, t, u, v, w, ncmo)] += + sign_K * sign_L * psi::C_DDOT(maxL, Cl[J], 1, Cr[I], 1); + } + } + } + } + } + } + } + return rdm; +} + +ambit::Tensor GenCIVector::compute_4rdm_aaab_same_irrep(GenCIVector& C_left, GenCIVector& C_right) { + size_t ncmo = C_left.ncmo_; + const auto& lists = C_left.lists_; + + auto rdm = ambit::Tensor::build(ambit::CoreTensor, "4RDM_AAAB", + {ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo}); + rdm.zero(); + auto& rdm_data = rdm.data(); + + int num_3h_class_Ka = lists->alfa_address_3h()->nclasses(); + int num_1h_class_Kb = lists->beta_address_1h()->nclasses(); + + for (int class_Ka = 0; class_Ka < num_3h_class_Ka; ++class_Ka) { + size_t maxKa = lists->alfa_address_3h()->strpcls(class_Ka); + + for (int class_Kb = 0; class_Kb < num_1h_class_Kb; ++class_Kb) { + size_t maxKb = lists->beta_address_1h()->strpcls(class_Kb); + + // loop over blocks of matrix C + for (const auto& [nI, class_Ia, class_Ib] : lists->determinant_classes()) { + if (lists->detpblk(nI) == 0) + continue; + + const auto Cr = C_right.C_[nI]->pointer(); + + for (const auto& [nJ, class_Ja, class_Jb] : lists->determinant_classes()) { + if (lists->detpblk(nJ) == 0) + continue; + + // Get a pointer to the correct block of matrix C + const auto Cl = C_left.C_[nJ]->pointer(); + + for (size_t Ka = 0; Ka < maxKa; ++Ka) { + auto& Ka_right_list = lists->get_alfa_3h_list(class_Ka, Ka, class_Ia); + auto& Ka_left_list = lists->get_alfa_3h_list(class_Ka, Ka, class_Ja); + for (size_t Kb = 0; Kb < maxKb; ++Kb) { + auto& Kb_right_list = lists->get_beta_1h_list(class_Kb, Kb, class_Ib); + auto& Kb_left_list = lists->get_beta_1h_list(class_Kb, Kb, class_Jb); + for (const auto& [sign_pqr, p, q, r, Ja] : Ka_left_list) { + for (const auto& [sign_s, s, Jb] : Kb_left_list) { + const double ClJ = sign_pqr * sign_s * Cl[Ja][Jb]; + for (const auto& [sign_tuv, t, u, v, Ia] : Ka_right_list) { + const auto CrIa = Cr[Ia]; + for (const auto& [sign_w, w, Ib] : Kb_right_list) { + rdm_data[eight_index(p, q, r, s, t, u, v, w, ncmo)] += + sign_tuv * sign_w * ClJ * CrIa[Ib]; + } + } + } + } + } + } + } + } + } + } + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = 0; q < p; ++q) { + for (size_t r = 0; r < q; ++r) { + for (size_t s = 0; s < ncmo; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = 0; u < t; ++u) { + for (size_t v = 0; v < u; ++v) { + for (size_t w = 0; w < ncmo; ++w) { + const double rdm_element = + rdm_data[eight_index(p, q, r, s, t, u, v, w, ncmo)]; + rdm_data[eight_index(p, q, r, s, t, v, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, r, s, u, t, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, r, s, u, v, t, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, r, s, v, t, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, r, s, v, u, t, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, q, s, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, q, s, t, v, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, q, s, u, t, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, q, s, u, v, t, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, q, s, v, t, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, q, s, v, u, t, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, p, r, s, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, r, s, t, v, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, p, r, s, u, t, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, p, r, s, u, v, t, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, r, s, v, t, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, r, s, v, u, t, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, r, p, s, t, u, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, r, p, s, t, v, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, r, p, s, u, t, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, r, p, s, u, v, t, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, r, p, s, v, t, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, r, p, s, v, u, t, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(r, p, q, s, t, u, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(r, p, q, s, t, v, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(r, p, q, s, u, t, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(r, p, q, s, u, v, t, w, ncmo)] = + rdm_element; + rdm_data[eight_index(r, p, q, s, v, t, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(r, p, q, s, v, u, t, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(r, q, p, s, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(r, q, p, s, t, v, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(r, q, p, s, u, t, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(r, q, p, s, u, v, t, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(r, q, p, s, v, t, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(r, q, p, s, v, u, t, w, ncmo)] = + rdm_element; + } + } + } + } + } + } + } + } + return rdm; +} + +ambit::Tensor GenCIVector::compute_4rdm_aabb_same_irrep(GenCIVector& C_left, GenCIVector& C_right) { + size_t ncmo = C_left.ncmo_; + const auto& lists = C_left.lists_; + + auto rdm = ambit::Tensor::build(ambit::CoreTensor, "4RDM_AABB", + {ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo}); + rdm.zero(); + auto& rdm_data = rdm.data(); + + int num_2h_class_Ka = lists->alfa_address_2h()->nclasses(); + int num_2h_class_Kb = lists->beta_address_2h()->nclasses(); + + for (int class_Ka = 0; class_Ka < num_2h_class_Ka; ++class_Ka) { + size_t maxKa = lists->alfa_address_2h()->strpcls(class_Ka); + + for (int class_Kb = 0; class_Kb < num_2h_class_Kb; ++class_Kb) { + size_t maxKb = lists->beta_address_2h()->strpcls(class_Kb); + + // loop over blocks of matrix C + for (const auto& [nI, class_Ia, class_Ib] : lists->determinant_classes()) { + if (lists->detpblk(nI) == 0) + continue; + + const auto Cr = C_right.C_[nI]->pointer(); + + for (const auto& [nJ, class_Ja, class_Jb] : lists->determinant_classes()) { + if (lists->detpblk(nJ) == 0) + continue; + + // Get a pointer to the correct block of matrix C + const auto Cl = C_left.C_[nJ]->pointer(); + + for (size_t Ka = 0; Ka < maxKa; ++Ka) { + auto& Ka_right_list = lists->get_alfa_2h_list(class_Ka, Ka, class_Ia); + auto& Ka_left_list = lists->get_alfa_2h_list(class_Ka, Ka, class_Ja); + for (size_t Kb = 0; Kb < maxKb; ++Kb) { + auto& Kb_right_list = lists->get_beta_2h_list(class_Kb, Kb, class_Ib); + auto& Kb_left_list = lists->get_beta_2h_list(class_Kb, Kb, class_Jb); + for (const auto& [sign_pq, p, q, Ja] : Ka_left_list) { + for (const auto& [sign_rs, r, s, Jb] : Kb_left_list) { + const double ClJ = sign_pq * sign_rs * Cl[Ja][Jb]; + for (const auto& [sign_tu, t, u, Ia] : Ka_right_list) { + const auto CrIa = Cr[Ia]; + for (const auto& [sign_vw, v, w, Ib] : Kb_right_list) { + rdm_data[eight_index(p, q, r, s, t, u, v, w, ncmo)] += + sign_tu * sign_vw * ClJ * CrIa[Ib]; + } + } + } + } + } + } + } + } + } + } + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = 0; q < p; ++q) { + for (size_t r = 0; r < ncmo; ++r) { + for (size_t s = 0; s < r; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = 0; u < t; ++u) { + for (size_t v = 0; v < ncmo; ++v) { + for (size_t w = 0; w < v; ++w) { + const double rdm_element = + rdm_data[eight_index(p, q, r, s, t, u, v, w, ncmo)]; + rdm_data[eight_index(p, q, r, s, t, u, w, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, r, s, u, t, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, r, s, u, t, w, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, s, r, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, s, r, t, u, w, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, s, r, u, t, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, s, r, u, t, w, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, r, s, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, r, s, t, u, w, v, ncmo)] = + rdm_element; + rdm_data[eight_index(q, p, r, s, u, t, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, p, r, s, u, t, w, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, s, r, t, u, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(q, p, s, r, t, u, w, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, s, r, u, t, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(q, p, s, r, u, t, w, v, ncmo)] = + rdm_element; + } + } + } + } + } + } + } + } + + return rdm; +} + +ambit::Tensor GenCIVector::compute_4rdm_abbb_same_irrep(GenCIVector& C_left, GenCIVector& C_right) { + size_t ncmo = C_left.ncmo_; + const auto& lists = C_left.lists_; + + auto rdm = ambit::Tensor::build(ambit::CoreTensor, "4RDM_ABBB", + {ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo, ncmo}); + rdm.zero(); + auto& rdm_data = rdm.data(); + + int num_1h_class_Ka = lists->alfa_address_1h()->nclasses(); + int num_3h_class_Kb = lists->beta_address_3h()->nclasses(); + + for (int class_Ka = 0; class_Ka < num_1h_class_Ka; ++class_Ka) { + size_t maxKa = lists->alfa_address_1h()->strpcls(class_Ka); + + for (int class_Kb = 0; class_Kb < num_3h_class_Kb; ++class_Kb) { + size_t maxKb = lists->beta_address_3h()->strpcls(class_Kb); + + // loop over blocks of matrix C + for (const auto& [nI, class_Ia, class_Ib] : lists->determinant_classes()) { + if (lists->detpblk(nI) == 0) + continue; + + const auto Cr = C_right.C_[nI]->pointer(); + + for (const auto& [nJ, class_Ja, class_Jb] : lists->determinant_classes()) { + if (lists->detpblk(nJ) == 0) + continue; + + // Get a pointer to the correct block of matrix C + const auto Cl = C_left.C_[nJ]->pointer(); + + for (size_t Ka = 0; Ka < maxKa; ++Ka) { + auto& Ka_right_list = lists->get_alfa_1h_list(class_Ka, Ka, class_Ia); + auto& Ka_left_list = lists->get_alfa_1h_list(class_Ka, Ka, class_Ja); + for (size_t Kb = 0; Kb < maxKb; ++Kb) { + auto& Kb_right_list = lists->get_beta_3h_list(class_Kb, Kb, class_Ib); + auto& Kb_left_list = lists->get_beta_3h_list(class_Kb, Kb, class_Jb); + for (const auto& [sign_p, p, Ja] : Ka_left_list) { + for (const auto& [sign_qrs, q, r, s, Jb] : Kb_left_list) { + const double ClJ = sign_p * sign_qrs * Cl[Ja][Jb]; + for (const auto& [sign_t, t, Ia] : Ka_right_list) { + const auto CrIa = Cr[Ia]; + for (const auto& [sign_uvw, u, v, w, Ib] : Kb_right_list) { + rdm_data[eight_index(p, q, r, s, t, u, v, w, ncmo)] += + sign_t * sign_uvw * ClJ * CrIa[Ib]; + } + } + } + } + } + } + } + } + } + } + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = 0; q < ncmo; ++q) { + for (size_t r = 0; r < q; ++r) { + for (size_t s = 0; s < r; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = 0; u < ncmo; ++u) { + for (size_t v = 0; v < u; ++v) { + for (size_t w = 0; w < v; ++w) { + const double rdm_element = + rdm_data[eight_index(p, q, r, s, t, u, v, w, ncmo)]; + rdm_data[eight_index(p, q, r, s, t, u, w, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, r, s, t, v, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, r, s, t, v, w, u, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, r, s, t, w, u, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, r, s, t, w, v, u, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, s, r, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, s, r, t, u, w, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, s, r, t, v, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, q, s, r, t, v, w, u, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, s, r, t, w, u, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, q, s, r, t, w, v, u, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, q, s, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, q, s, t, u, w, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, q, s, t, v, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, q, s, t, v, w, u, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, q, s, t, w, u, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, q, s, t, w, v, u, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, s, q, t, u, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, s, q, t, u, w, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, s, q, t, v, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, r, s, q, t, v, w, u, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, s, q, t, w, u, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, r, s, q, t, w, v, u, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, s, q, r, t, u, v, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, s, q, r, t, u, w, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, s, q, r, t, v, u, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, s, q, r, t, v, w, u, ncmo)] = + rdm_element; + rdm_data[eight_index(p, s, q, r, t, w, u, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, s, q, r, t, w, v, u, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, s, r, q, t, u, v, w, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, s, r, q, t, u, w, v, ncmo)] = + rdm_element; + rdm_data[eight_index(p, s, r, q, t, v, u, w, ncmo)] = + rdm_element; + rdm_data[eight_index(p, s, r, q, t, v, w, u, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, s, r, q, t, w, u, v, ncmo)] = + -rdm_element; + rdm_data[eight_index(p, s, r, q, t, w, v, u, ncmo)] = + rdm_element; + } + } + } + } + } + } + } + } + return rdm; +} + void GenCIVector::test_rdms(GenCIVector& Cl, GenCIVector& Cr, int max_rdm_level, RDMsType type, std::shared_ptr rdms) { size_t ncmo = Cl.ncmo_; @@ -792,8 +1321,7 @@ void GenCIVector::test_rdms(GenCIVector& Cl, GenCIVector& Cr, int max_rdm_level, if (max_rdm_level >= 3) { auto g3aaa = rdms->g3aaa(); double error_3rdm_aaa = 0.0; - // for (size_t p = 0; p < no_; ++p){ - for (size_t p = 0; p < 1; ++p) { + for (size_t p = 0; p < ncmo; ++p) { for (size_t q = p + 1; q < ncmo; ++q) { for (size_t r = q + 1; r < ncmo; ++r) { for (size_t s = 0; s < ncmo; ++s) { @@ -904,8 +1432,7 @@ void GenCIVector::test_rdms(GenCIVector& Cl, GenCIVector& Cr, int max_rdm_level, auto g3bbb = rdms->g3bbb(); double error_3rdm_bbb = 0.0; - for (size_t p = 0; p < 1; ++p) { - // for (size_t p = 0; p < no_; ++p){ + for (size_t p = 0; p < ncmo; ++p) { for (size_t q = p + 1; q < ncmo; ++q) { for (size_t r = q + 1; r < ncmo; ++r) { for (size_t s = 0; s < ncmo; ++s) { @@ -940,6 +1467,223 @@ void GenCIVector::test_rdms(GenCIVector& Cl, GenCIVector& Cr, int max_rdm_level, psi::Process::environment.globals["BBBBBB 3-RDM ERROR"] = error_3rdm_bbb; psi::outfile->Printf("\n BBBBBB 3-RDM Error : %+e", error_3rdm_bbb); } + + if (max_rdm_level >= 4) { + auto g4aaaa = rdms->g4aaaa(); + double error_4rdm_aaaa = 0.0; + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = p + 1; q < ncmo; ++q) { + for (size_t r = q + 1; r < ncmo; ++r) { + for (size_t s = r + 1; s < ncmo; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = t + 1; u < ncmo; ++u) { + for (size_t v = u + 1; v < ncmo; ++v) { + for (size_t w = v + 1; w < ncmo; ++w) { + double rdm = 0.0; + for (const auto& [I, c_I] : state_vector_r) { + J = I; + double sign = 1.0; + sign *= J.destroy_alfa_bit(t); + sign *= J.destroy_alfa_bit(u); + sign *= J.destroy_alfa_bit(v); + sign *= J.destroy_alfa_bit(w); + sign *= J.create_alfa_bit(s); + sign *= J.create_alfa_bit(r); + sign *= J.create_alfa_bit(q); + sign *= J.create_alfa_bit(p); + if (sign != 0) { + if (state_vector_l.count(J) != 0) { + rdm += sign * to_double(state_vector_l[J] * c_I); + } + } + } + if (std::fabs(rdm) > 1.0e-12) { + double rdm_comp = g4aaaa.at({p, q, r, s, t, u, v, w}); + error_4rdm_aaaa += std::fabs(rdm - rdm_comp); + } + } + } + } + } + } + } + } + } + psi::Process::environment.globals["AAAAAAAA 4-RDM ERROR"] = error_4rdm_aaaa; + psi::outfile->Printf("\n AAAAAAAA 4-RDM Error : %+e", error_4rdm_aaaa); + + auto g4aaab = rdms->g4aaab(); + double error_4rdm_aaab = 0.0; + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = p + 1; q < ncmo; ++q) { + for (size_t r = q + 1; r < ncmo; ++r) { + for (size_t s = 0; s < ncmo; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = t + 1; u < ncmo; ++u) { + for (size_t v = u + 1; v < ncmo; ++v) { + for (size_t w = 0; w < ncmo; ++w) { + double rdm = 0.0; + for (const auto& [I, c_I] : state_vector_r) { + J = I; + double sign = 1.0; + sign *= J.destroy_alfa_bit(t); + sign *= J.destroy_alfa_bit(u); + sign *= J.destroy_alfa_bit(v); + sign *= J.destroy_beta_bit(w); + sign *= J.create_beta_bit(s); + sign *= J.create_alfa_bit(r); + sign *= J.create_alfa_bit(q); + sign *= J.create_alfa_bit(p); + if (sign != 0) { + if (state_vector_l.count(J) != 0) { + rdm += sign * to_double(state_vector_l[J] * c_I); + } + } + } + if (std::fabs(rdm) > 1.0e-12) { + double rdm_comp = g4aaab.at({p, q, r, s, t, u, v, w}); + error_4rdm_aaab += std::fabs(rdm - rdm_comp); + } + } + } + } + } + } + } + } + } + psi::Process::environment.globals["AAABAAAB 4-RDM ERROR"] = error_4rdm_aaab; + psi::outfile->Printf("\n AAABAAAB 4-RDM Error : %+e", error_4rdm_aaab); + + auto g4aabb = rdms->g4aabb(); + double error_4rdm_aabb = 0.0; + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = p + 1; q < ncmo; ++q) { + for (size_t r = 0; r < ncmo; ++r) { + for (size_t s = r + 1; s < ncmo; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = t + 1; u < ncmo; ++u) { + for (size_t v = 0; v < ncmo; ++v) { + for (size_t w = v + 1; w < ncmo; ++w) { + double rdm = 0.0; + for (const auto& [I, c_I] : state_vector_r) { + J = I; + double sign = 1.0; + sign *= J.destroy_alfa_bit(t); + sign *= J.destroy_alfa_bit(u); + sign *= J.destroy_beta_bit(v); + sign *= J.destroy_beta_bit(w); + sign *= J.create_beta_bit(s); + sign *= J.create_beta_bit(r); + sign *= J.create_alfa_bit(q); + sign *= J.create_alfa_bit(p); + if (sign != 0) { + if (state_vector_l.count(J) != 0) { + rdm += sign * to_double(state_vector_l[J] * c_I); + } + } + } + if (std::fabs(rdm) > 1.0e-12) { + double rdm_comp = g4aabb.at({p, q, r, s, t, u, v, w}); + error_4rdm_aabb += std::fabs(rdm - rdm_comp); + } + } + } + } + } + } + } + } + } + psi::Process::environment.globals["AABBAABB 4-RDM ERROR"] = error_4rdm_aabb; + psi::outfile->Printf("\n AABBAABB 4-RDM Error : %+e", error_4rdm_aabb); + + auto g4abbb = rdms->g4abbb(); + double error_4rdm_abbb = 0.0; + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = 0; q < ncmo; ++q) { + for (size_t r = q + 1; r < ncmo; ++r) { + for (size_t s = r + 1; s < ncmo; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = 0; u < ncmo; ++u) { + for (size_t v = u + 1; v < ncmo; ++v) { + for (size_t w = v + 1; w < ncmo; ++w) { + double rdm = 0.0; + for (const auto& [I, c_I] : state_vector_r) { + J = I; + double sign = 1.0; + sign *= J.destroy_alfa_bit(t); + sign *= J.destroy_beta_bit(u); + sign *= J.destroy_beta_bit(v); + sign *= J.destroy_beta_bit(w); + sign *= J.create_beta_bit(s); + sign *= J.create_beta_bit(r); + sign *= J.create_beta_bit(q); + sign *= J.create_alfa_bit(p); + if (sign != 0) { + if (state_vector_l.count(J) != 0) { + rdm += sign * to_double(state_vector_l[J] * c_I); + } + } + } + if (std::fabs(rdm) > 1.0e-12) { + double rdm_comp = g4abbb.at({p, q, r, s, t, u, v, w}); + error_4rdm_abbb += std::fabs(rdm - rdm_comp); + } + } + } + } + } + } + } + } + } + psi::Process::environment.globals["ABBBABBB 4-RDM ERROR"] = error_4rdm_abbb; + psi::outfile->Printf("\n ABBBABBB 4-RDM Error : %+e", error_4rdm_abbb); + + auto g4bbbb = rdms->g4bbbb(); + double error_4rdm_bbbb = 0.0; + for (size_t p = 0; p < ncmo; ++p) { + for (size_t q = p + 1; q < ncmo; ++q) { + for (size_t r = q + 1; r < ncmo; ++r) { + for (size_t s = r + 1; s < ncmo; ++s) { + for (size_t t = 0; t < ncmo; ++t) { + for (size_t u = t + 1; u < ncmo; ++u) { + for (size_t v = u + 1; v < ncmo; ++v) { + for (size_t w = v + 1; w < ncmo; ++w) { + double rdm = 0.0; + for (const auto& [I, c_I] : state_vector_r) { + J = I; + double sign = 1.0; + sign *= J.destroy_beta_bit(t); + sign *= J.destroy_beta_bit(u); + sign *= J.destroy_beta_bit(v); + sign *= J.destroy_beta_bit(w); + sign *= J.create_beta_bit(s); + sign *= J.create_beta_bit(r); + sign *= J.create_beta_bit(q); + sign *= J.create_beta_bit(p); + if (sign != 0) { + if (state_vector_l.count(J) != 0) { + rdm += sign * to_double(state_vector_l[J] * c_I); + } + } + } + if (std::fabs(rdm) > 1.0e-12) { + double rdm_comp = g4bbbb.at({p, q, r, s, t, u, v, w}); + error_4rdm_bbbb += std::fabs(rdm - rdm_comp); + } + } + } + } + } + } + } + } + } + psi::Process::environment.globals["BBBBBBBB 4-RDM ERROR"] = error_4rdm_bbbb; + psi::outfile->Printf("\n BBBBBBBB 4-RDM Error : %+e", error_4rdm_bbbb); + } } } // namespace forte diff --git a/forte/genci/string_lists_makers.cc b/forte/genci/string_lists_makers.cc index 7ec58c412..50af4be48 100644 --- a/forte/genci/string_lists_makers.cc +++ b/forte/genci/string_lists_makers.cc @@ -387,6 +387,97 @@ H3List make_3h_list(const StringList& strings, std::shared_ptr ad return list; } +H4List make_4h_list(const StringList& strings, std::shared_ptr addresser, + std::shared_ptr addresser_4h) { + H4List list; + int n = addresser->nbits(); + int k = addresser->nones(); + size_t nmo = addresser->nbits(); + + if ((k >= 0) and (k <= n)) { + for (const auto& string_class : strings) { + for (const auto& I : string_class) { + const auto& [add_I, class_I] = addresser->address_and_class(I); + for (size_t s = 0; s < nmo; ++s) { + for (size_t r = s + 1; r < nmo; ++r) { + for (size_t q = r + 1; q < nmo; ++q) { + for (size_t p = q + 1; p < nmo; ++p) { + if (I[p] and I[q] and I[r] and I[s]) { + auto J = I; + J[s] = false; + const auto s_sign = J.slater_sign(s); + J[r] = false; + const auto r_sign = J.slater_sign(r); + J[q] = false; + const auto q_sign = J.slater_sign(q); + J[p] = false; + const auto p_sign = J.slater_sign(p); + if (auto it = addresser_4h->find(J); + it != addresser_4h->end()) { + const auto sign = p_sign * q_sign * r_sign * s_sign; + const auto& [add_J, class_J] = it->second; + std::tuple I_tuple(class_J, add_J, + class_I); + list[I_tuple].push_back( + H4StringSubstitution(-sign, p, q, r, s, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, p, q, s, r, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, p, r, q, s, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, p, r, s, q, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, p, s, q, r, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, p, s, r, q, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, q, p, r, s, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, q, p, s, r, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, q, r, p, s, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, q, r, s, p, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, q, s, p, r, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, q, s, r, p, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, r, p, q, s, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, r, p, s, q, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, r, q, p, s, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, r, q, s, p, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, r, s, p, q, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, r, s, q, p, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, s, p, q, r, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, s, p, r, q, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, s, q, p, r, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, s, q, r, p, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(+sign, s, r, p, q, add_I)); + list[I_tuple].push_back( + H4StringSubstitution(-sign, s, r, q, p, add_I)); + } + } + } + } + } + } + } + } + } + return list; +} + std::map, std::vector>> find_string_map(const GenCIStringLists& list_left, const GenCIStringLists& list_right, bool alfa) { std::map, std::vector>> m; diff --git a/forte/genci/string_lists_makers.h b/forte/genci/string_lists_makers.h index e2a3ee6dd..7857dc53f 100644 --- a/forte/genci/string_lists_makers.h +++ b/forte/genci/string_lists_makers.h @@ -79,5 +79,8 @@ H2List make_2h_list(const StringList& strings, std::shared_ptr ad /// Make 3-hole lists (I -> a_p a_q a_r I = sgn J) H3List make_3h_list(const StringList& strings, std::shared_ptr address, std::shared_ptr address_3h); +/// Make 4-hole lists (I -> a_p a_q a_r a_s I = sgn J) +H4List make_4h_list(const StringList& strings, std::shared_ptr address, + std::shared_ptr address_4h); } // namespace forte diff --git a/forte/helpers/helpers.cc b/forte/helpers/helpers.cc index dec29ba58..e093ad2fc 100644 --- a/forte/helpers/helpers.cc +++ b/forte/helpers/helpers.cc @@ -74,6 +74,15 @@ py::array_t vector_to_np(const std::vector& v, const std::vector return py::array_t(dims, &(v.data()[0])); } +py::dict blockedtensor_to_np(ambit::BlockedTensor bt) { + py::dict pybt; + auto labels = bt.block_labels(); + for (const auto& label : labels) { + pybt[py::str(label)] = ambit_to_np(bt.block(label)); + } + return pybt; +} + std::shared_ptr tensor_to_matrix(ambit::Tensor t) { size_t size1 = t.dim(0); size_t size2 = t.dim(1); diff --git a/forte/helpers/helpers.h b/forte/helpers/helpers.h index 3905e4dca..2a7f99a80 100644 --- a/forte/helpers/helpers.h +++ b/forte/helpers/helpers.h @@ -66,6 +66,14 @@ enum class Spin3 { aaa, aab, abb, bbb }; */ py::array_t ambit_to_np(ambit::Tensor t); +/** + * @brief Convert an ambit BlockedTensor to a dictionary of numpy ndarray's. + * The returned tensor is stored according to the C storage convention. + * @param t The input BlockedTensor + * @return A dictionary of numpy arrays + */ +py::dict blockedtensor_to_np(ambit::BlockedTensor bt); + /** * @brief Convert a std::vector to a numpy ndarray. * The returned tensor is stored according to the C storage convention. diff --git a/forte/mrdsrg-so/mrdsrg_so.cc b/forte/mrdsrg-so/mrdsrg_so.cc index 331407a52..184c93d80 100644 --- a/forte/mrdsrg-so/mrdsrg_so.cc +++ b/forte/mrdsrg-so/mrdsrg_so.cc @@ -33,6 +33,9 @@ #include "psi4/libpsi4util/process.h" #include "psi4/libmints/molecule.h" +#include "psi4/libmints/matrix.h" +#include "psi4/libmints/vector.h" + #include "integrals/active_space_integrals.h" #include "base_classes/mo_space_info.h" #include "base_classes/state_info.h" @@ -50,6 +53,7 @@ MRDSRG_SO::MRDSRG_SO(std::shared_ptr rdms, std::shared_ptr scf_in std::shared_ptr mo_space_info) : DynamicCorrelationSolver(rdms, scf_info, options, ints, mo_space_info), BTF_(std::make_shared()), tensor_type_(ambit::CoreTensor) { + BlockedTensor::reset_mo_spaces(); print_method_banner( {"SO-Based Multireference Driven Similarity Renormalization Group", "Chenyang Li"}); startup(); @@ -642,6 +646,10 @@ double MRDSRG_SO::compute_energy() { ++cycle; } while (!converged); + if (foptions_->get_bool("FULL_HBAR")) { + compute_hbar(); + } + outfile->Printf("\n " "----------------------------------------------------------" "----------------------------------------"); @@ -722,7 +730,7 @@ void MRDSRG_SO::compute_hbar() { O2["pqrs"] = C2["pqrs"]; C2["pqrs"] += O2["rspq"]; } - + // Hbar += C Hbar0 += C0; Hbar1["pq"] += C1["pq"]; @@ -1368,4 +1376,13 @@ void MRDSRG_SO::H3_T2_C2(BlockedTensor& H3, BlockedTensor& T2, const double& alp temp["mx"] += 0.5 * T2["myuv"] * Lambda2["uvxy"]; C2["toqr"] -= alpha * temp["mx"] * H3["xtomqr"]; } +std::pair> MRDSRG_SO::save_Heff_full() { + double Edsrg = Eref + Hbar0; + ambit::BlockedTensor Hbar1_copy = BTF_->build(tensor_type_, "Hbar1_copy", {"gg"}); + ambit::BlockedTensor Hbar2_copy = BTF_->build(tensor_type_, "Hbar2_copy", {"gggg"}); + Hbar1_copy["pq"] = Hbar1["pq"]; + Hbar2_copy["pqrs"] = Hbar2["pqrs"]; + std::vector Heff = {Hbar1_copy, Hbar2_copy}; + return std::make_pair(Edsrg, Heff); +} } // namespace forte diff --git a/forte/mrdsrg-so/mrdsrg_so.h b/forte/mrdsrg-so/mrdsrg_so.h index dbdc58595..4885fa11f 100644 --- a/forte/mrdsrg-so/mrdsrg_so.h +++ b/forte/mrdsrg-so/mrdsrg_so.h @@ -281,6 +281,14 @@ class MRDSRG_SO : public DynamicCorrelationSolver { /// Compute the DSRG-MRPT2 energy double compute_energy(); + // double compute_eom(); + + std::pair> save_Heff_full(); + ambit::BlockedTensor get_gamma1() { return Gamma1; } + ambit::BlockedTensor get_eta1() { return Eta1; } + ambit::BlockedTensor get_lambda2() { return Lambda2; } + ambit::BlockedTensor get_lambda3() { return Lambda3; } + /// DSRG transformed Hamiltonian (not implemented) std::shared_ptr compute_Heff_actv(); diff --git a/forte/mrdsrg-spin-adapted/sa_ldsrg2.cc b/forte/mrdsrg-spin-adapted/sa_ldsrg2.cc index 23c9d9004..e3e460ab4 100644 --- a/forte/mrdsrg-spin-adapted/sa_ldsrg2.cc +++ b/forte/mrdsrg-spin-adapted/sa_ldsrg2.cc @@ -192,9 +192,14 @@ double SA_MRDSRG::compute_energy_ldsrg2() { // dump amplitudes to disk dump_amps_to_disk(); + if (full_hbar_) { + int ncomm = compute_hbar(rsc_conv_); + } + final.stop(); Hbar0_ = Ecorr; + return Ecorr; } diff --git a/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc b/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc index 46196f5f4..e67c03fd0 100644 --- a/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc +++ b/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc @@ -95,6 +95,8 @@ void SA_MRDSRG::read_options() { outfile->Printf("\n Changed DSRG_TRANS_TYPE option to UNITARY"); dsrg_trans_type_ = "UNITARY"; } + + full_hbar_ = foptions_->get_bool("FULL_HBAR"); } void SA_MRDSRG::startup() { diff --git a/forte/mrdsrg-spin-adapted/sa_mrdsrg.h b/forte/mrdsrg-spin-adapted/sa_mrdsrg.h index fba0d95b4..690fca8e4 100644 --- a/forte/mrdsrg-spin-adapted/sa_mrdsrg.h +++ b/forte/mrdsrg-spin-adapted/sa_mrdsrg.h @@ -107,6 +107,8 @@ class SA_MRDSRG : public SADSRG { double rsc_conv_adapt_threshold_; double rsc_conv_adapt_delta_e_; + bool full_hbar_; + /// DSRG transformation type std::string dsrg_trans_type_; @@ -165,6 +167,8 @@ class SA_MRDSRG : public SADSRG { int compute_hbar_sequential(double& rsc_conv); /// Compute DSRG-transformed Hamiltonian truncated to 2-nested commutator void compute_hbar_qc(); + // /// Compute EOM-LDSRG2 + // double compute_eom(); /// Add H2's Hermitian conjugate to itself, H2 need to contain gGgG block void add_hermitian_conjugate(BlockedTensor& H2); @@ -196,7 +200,7 @@ class SA_MRDSRG : public SADSRG { /// Clean up for pointers used for DIIS void diis_manager_cleanup(); - /// Compute MR-LDSRG(2) + /// Compute MR-LDSRG(2)s double compute_energy_ldsrg2(); /// Compute DSRG-transformed multipoles diff --git a/forte/mrdsrg-spin-adapted/sadsrg.cc b/forte/mrdsrg-spin-adapted/sadsrg.cc index fdcb3595f..3683e0243 100644 --- a/forte/mrdsrg-spin-adapted/sadsrg.cc +++ b/forte/mrdsrg-spin-adapted/sadsrg.cc @@ -232,7 +232,9 @@ void SADSRG::read_MOSpaceInfo() { core_mos_ = mo_space_info_->corr_absolute_mo("RESTRICTED_DOCC"); actv_mos_ = mo_space_info_->corr_absolute_mo("ACTIVE"); virt_mos_ = mo_space_info_->corr_absolute_mo("RESTRICTED_UOCC"); + gen_mos_ = mo_space_info_->corr_absolute_mo("CORRELATED"); actv_mos_sym_ = mo_space_info_->symmetry("ACTIVE"); + gen_mos_sym_ = mo_space_info_->symmetry("CORRELATED"); if (eri_df_) { aux_mos_ = std::vector(ints_->nthree()); @@ -472,6 +474,22 @@ std::shared_ptr SADSRG::compute_Heff_actv() { return fci_ints; } +std::vector SADSRG::compute_Heff_full() { + std::vector Heff = {Hbar1_, Hbar2_}; + return Heff; +} + +std::pair> SADSRG::compute_Heff_full_degno() { + double Edsrg = Eref_ + Hbar0_; + if (foptions_->get_bool("FORM_HBAR3")) { + throw psi::PSIEXCEPTION("FORM_HBAR3 is not implemented for full Hamiltonian."); + } else { + deGNO_ints_full("Hamiltonian", Edsrg, Hbar1_, Hbar2_); + } + std::vector Heff = {Hbar1_, Hbar2_}; + return std::make_pair(Edsrg, Heff); +} + std::shared_ptr SADSRG::compute_mp_eff_actv() { /** * DSRG transform multipole integrals: Mbar = e^{-A} M e^{A} @@ -707,6 +725,44 @@ void SADSRG::deGNO_ints(const std::string& name, double& H0, BlockedTensor& H1, print_done(t2.get()); } +void SADSRG::deGNO_ints_full(const std::string& name, double& H0, BlockedTensor& H1, BlockedTensor& H2) { + print_h2("De-Normal-Order Full 2-Body DSRG Transformed " + name); + + // compute scalar + local_timer t0; + print_contents("Computing the scalar term"); + + // build a temp["pqrs"] = 2 * H2["pqrs"] - H2["pqsr"] + auto temp = H2.clone(); + temp.scale(2.0); + temp["pqrs"] -= H2["pqsr"]; + + auto L1h = BTF_->build(tensor_type_, "L1h", spin_cases({"hh"})); + L1h.block("aa")("uv") = L1_.block("aa")("uv"); + L1h.block("cc").iterate( + [&](const std::vector& i, double& value) { value = i[0] == i[1] ? 2.0 : 0.0; }); + + // scalar from H1 + double scalar1 = 0.0; + scalar1 -= H1["ji"] * L1h["ij"]; + + // scalar from H2 + double scalar2 = 0.0; + scalar2 += 0.25 * L1h["ij"] * temp["jlik"] * L1h["kl"]; + + scalar2 -= 0.5 * H2["xyuv"] * L2_["uvxy"]; + + H0 += scalar1 + scalar2; + print_done(t0.get()); + + // compute 1-body term + local_timer t1; + print_contents("Computing the 1-body term"); + + H1["pq"] -= 0.5 * temp["piqj"] * L1h["ji"]; + print_done(t1.get()); +} + ambit::BlockedTensor SADSRG::deGNO_Tamp(BlockedTensor& T1, BlockedTensor& T2, BlockedTensor& D1) { BlockedTensor T1eff = BTF_->build(tensor_type_, "T1eff from de-GNO", {"hp"}); diff --git a/forte/mrdsrg-spin-adapted/sadsrg.h b/forte/mrdsrg-spin-adapted/sadsrg.h index c4eaaf8eb..51b389fb7 100644 --- a/forte/mrdsrg-spin-adapted/sadsrg.h +++ b/forte/mrdsrg-spin-adapted/sadsrg.h @@ -62,6 +62,8 @@ class SADSRG : public DynamicCorrelationSolver { /// Compute DSRG transformed Hamiltonian std::shared_ptr compute_Heff_actv() override; + std::vector compute_Heff_full(); + std::pair> compute_Heff_full_degno(); /// Compute DSRG transformed multipole integrals std::shared_ptr compute_mp_eff_actv() override; @@ -181,9 +183,13 @@ class SADSRG : public DynamicCorrelationSolver { std::vector actv_mos_; /// List of virtual MOs std::vector virt_mos_; + /// List of genereal MOs + std::vector gen_mos_; /// List of the symmetry of the active MOs std::vector actv_mos_sym_; - + /// List of the symmetry of the general MOs + std::vector gen_mos_sym_; + /// List of active active occupied MOs (relative to active) std::vector actv_occ_mos_; /// List of active active unoccupied MOs (relative to active) @@ -326,6 +332,9 @@ class SADSRG : public DynamicCorrelationSolver { /// DSRG transformed 3-body Hamiltonian (active only in DSRG-PT, but full in MRDSRG) ambit::BlockedTensor Hbar3_; + // /// EOM-LDSRG2 Hamiltonian: + // std::shared_ptr EOM_Hbar_mat_; + /// Scalar of the DSRG transformed multipoles std::vector Mbar0_; /// DSRG transformed 1-body multipoles @@ -352,6 +361,12 @@ class SADSRG : public DynamicCorrelationSolver { void deGNO_ints(const std::string& name, double& H0, BlockedTensor& H1, BlockedTensor& H2, BlockedTensor& H3); + /** + * De-normal-order a 2-body DSRG transformed integrals (full Hbar!) + * This will change H0 and H1 !!! + */ + void deGNO_ints_full(const std::string& name, double& H0, BlockedTensor& H1, BlockedTensor& H2); + /** * De-normal-order the T1 and T2 amplitudes and return the effective T1 * T1eff = T1 - T2["ivau"] * D1["uv"] diff --git a/forte/mrdsrg-spin-integrated/master_mrdsrg.cc b/forte/mrdsrg-spin-integrated/master_mrdsrg.cc index f3c418a6b..024f3ef2c 100644 --- a/forte/mrdsrg-spin-integrated/master_mrdsrg.cc +++ b/forte/mrdsrg-spin-integrated/master_mrdsrg.cc @@ -154,6 +154,7 @@ void MASTER_DSRG::read_options() { } do_cu3_ = (foptions_->get_str("THREEPDC") != "ZERO"); + do_cu4_ = (foptions_->get_str("FOURPDC") != "ZERO"); outfile->Printf("Done"); } @@ -264,6 +265,15 @@ void MASTER_DSRG::fill_density() { L3abb_ = rdms_->L3abb(); L3bbb_ = rdms_->L3bbb(); } + + // 4-body density cumulants + if (do_cu4_) { + L4aaaa_ = rdms_->L4aaaa(); + L4aaab_ = rdms_->L4aaab(); + L4aabb_ = rdms_->L4aabb(); + L4abbb_ = rdms_->L4abbb(); + L4bbbb_ = rdms_->L4bbbb(); + } } void MASTER_DSRG::init_fock() { @@ -681,6 +691,103 @@ void MASTER_DSRG::deGNO_ints(const std::string& name, double& H0, BlockedTensor& outfile->Printf("Done. Timing %8.3f s", t2.get()); } +void MASTER_DSRG::deGNO_ints_full(const std::string& name, double& H0, BlockedTensor& H1, + BlockedTensor& H2) { + print_h2("De-Normal-Order Full DSRG Transformed " + name); + + // compute scalar + local_timer t0; + outfile->Printf("\n %-40s ... ", "Computing the scalar term"); + + auto L1h = BTF_->build(tensor_type_, "L1h", spin_cases({"hh"})); + L1h["uv"] = Gamma1_["uv"]; + L1h["UV"] = Gamma1_["UV"]; + L1h.block("cc").iterate( + [&](const std::vector& i, double& value) { value = i[0] == i[1] ? 1.0 : 0.0; }); + L1h.block("CC").iterate( + [&](const std::vector& i, double& value) { value = i[0] == i[1] ? 1.0 : 0.0; }); + + // scalar from H1 + double scalar1 = 0.0; + scalar1 -= H1["ji"] * L1h["ij"]; + scalar1 -= H1["JI"] * L1h["IJ"]; + + // scalar from H2 + double scalar2 = 0.0; + scalar2 += 0.5 * L1h["ij"] * H2["jlik"] * L1h["kl"]; + scalar2 += 0.5 * L1h["IJ"] * H2["JLIK"] * L1h["KL"]; + scalar2 += L1h["ij"] * H2["jLiK"] * L1h["KL"]; + + scalar2 -= 0.25 * H2["xyuv"] * Lambda2_["uvxy"]; + scalar2 -= 0.25 * H2["XYUV"] * Lambda2_["UVXY"]; + scalar2 -= H2["xYuV"] * Lambda2_["uVxY"]; + + H0 += scalar1 + scalar2; + outfile->Printf("Done. Timing %8.3f s", t0.get()); + + // compute 1-body term + local_timer t1; + outfile->Printf("\n %-40s ... ", "Computing the 1-body term"); + + H1["pq"] -= H2["piqj"] * L1h["ji"]; + H1["pq"] -= H2["pIqJ"] * L1h["JI"]; + H1["PQ"] -= H2["iPjQ"] * L1h["ji"]; + H1["PQ"] -= H2["PIQJ"] * L1h["JI"]; + outfile->Printf("Done. Timing %8.3f s", t1.get()); +} + +void MASTER_DSRG::GNO_ints_full(const std::string& name, double& H0, BlockedTensor& H1, + BlockedTensor& H2, const std::shared_ptr rdms) { + print_h2("Normal-Order Full DSRG Transformed " + name); + + // compute scalar + local_timer t0; + outfile->Printf("\n %-40s ... ", "Computing the scalar term"); + + auto L1h = BTF_->build(tensor_type_, "L1h", spin_cases({"hh"})); + L1h.block("aa")("pq") = rdms->g1a()("pq"); + L1h.block("AA")("pq") = rdms->g1b()("pq"); + L1h.block("cc").iterate( + [&](const std::vector& i, double& value) { value = i[0] == i[1] ? 1.0 : 0.0; }); + L1h.block("CC").iterate( + [&](const std::vector& i, double& value) { value = i[0] == i[1] ? 1.0 : 0.0; }); + + // compute 1-body term + local_timer t1; + outfile->Printf("\n %-40s ... ", "Computing the 1-body term"); + + H1["pq"] += H2["piqj"] * L1h["ji"]; + H1["pq"] += H2["pIqJ"] * L1h["JI"]; + H1["PQ"] += H2["iPjQ"] * L1h["ji"]; + H1["PQ"] += H2["PIQJ"] * L1h["JI"]; + outfile->Printf("Done. Timing %8.3f s", t1.get()); + + // scalar from H1 + double scalar1 = 0.0; + scalar1 += H1["ji"] * L1h["ij"]; + scalar1 += H1["JI"] * L1h["IJ"]; + + // scalar from H2 + double scalar2 = 0.0; + scalar2 -= 0.5 * L1h["ij"] * H2["jlik"] * L1h["kl"]; + scalar2 -= 0.5 * L1h["IJ"] * H2["JLIK"] * L1h["KL"]; + scalar2 -= L1h["ij"] * H2["jLiK"] * L1h["KL"]; + + auto Lambda2 = BTF_->build(tensor_type_, "Lambda2", spin_cases({"aaaa"})); + + Lambda2.block("aaaa")("pqrs") = rdms->L2aa()("pqrs"); + Lambda2.block("aAaA")("pqrs") = rdms->L2ab()("pqrs"); + Lambda2.block("AAAA")("pqrs") = rdms->L2bb()("pqrs"); + + scalar2 += 0.25 * H2["xyuv"] * Lambda2["uvxy"]; + scalar2 += 0.25 * H2["XYUV"] * Lambda2["UVXY"]; + scalar2 += H2["xYuV"] * Lambda2["uVxY"]; + + H0 += scalar1 + scalar2; + outfile->Printf("Done. Timing %8.3f s", t0.get()); +} + + void MASTER_DSRG::fill_three_index_ints(ambit::BlockedTensor T) { const auto& block_labels = T.block_labels(); for (const std::string& string_block : block_labels) { @@ -2286,4 +2393,39 @@ std::vector MASTER_DSRG::re_two_labels() { return labels; } + +std::pair> MASTER_DSRG::compute_Heff_full_degno() { + double Edsrg = Eref_ + Hbar0_; + if (foptions_->get_bool("FORM_HBAR3")) { + throw psi::PSIEXCEPTION("FORM_HBAR3 is not implemented for full Hamiltonian."); + } else { + deGNO_ints_full("Hamiltonian", Edsrg, Hbar1_, Hbar2_); + } + std::vector Heff = {Hbar1_, Hbar2_}; + return std::make_pair(Edsrg, Heff); +} + +std::pair> MASTER_DSRG::save_Heff_full() { + double Edsrg = Eref_ + Hbar0_; + ambit::BlockedTensor Hbar1_copy = BTF_->build(tensor_type_, "Hbar1_copy", spin_cases({"gg"})); + ambit::BlockedTensor Hbar2_copy = BTF_->build(tensor_type_, "Hbar2_copy", spin_cases({"gggg"})); + Hbar1_copy["pq"] = Hbar1_["pq"]; + Hbar1_copy["PQ"] = Hbar1_["PQ"]; + Hbar2_copy["pqrs"] = Hbar2_["pqrs"]; + Hbar2_copy["pQrS"] = Hbar2_["pQrS"]; + Hbar2_copy["PQRS"] = Hbar2_["PQRS"]; + std::vector Heff = {Hbar1_copy, Hbar2_copy}; + return std::make_pair(Edsrg, Heff); +} + +std::pair> MASTER_DSRG::update_Heff_full(double& H0, BlockedTensor& H1, + BlockedTensor& H2, std::shared_ptr rdms) { + deGNO_ints_full("Hamiltonian", H0, H1, H2); // De-normal ordering based on original RDMs. + outfile->Printf("\n H0 after deGNO: %20.12f", H0); + GNO_ints_full("Hamiltonian", H0, H1, H2, rdms); // Re-normal ordering based on new RDMs. + outfile->Printf("\n H0 after GNO: %20.12f", H0); + std::vector Heff = {H1, H2}; + return std::make_pair(H0, Heff); +} + } // namespace forte diff --git a/forte/mrdsrg-spin-integrated/master_mrdsrg.h b/forte/mrdsrg-spin-integrated/master_mrdsrg.h index f83a217dc..3af9406a1 100644 --- a/forte/mrdsrg-spin-integrated/master_mrdsrg.h +++ b/forte/mrdsrg-spin-integrated/master_mrdsrg.h @@ -47,6 +47,37 @@ class MASTER_DSRG : public DynamicCorrelationSolver { /// Compute DSRG transformed Hamiltonian virtual std::shared_ptr compute_Heff_actv(); + /// Modify the rdms_ member of the parent (base) class. Only used in EOM-DSRG. + virtual void set_rdms(std::shared_ptr rdms) { + throw std::runtime_error("Child class should override this function"); + }; + + /// Save DSRG full transformed Hamiltonian + std::pair> save_Heff_full(); + + /// Update DSRG full transformed Hamiltonian + std::pair> update_Heff_full(double& H0, BlockedTensor& H1, \ + BlockedTensor& H2, std::shared_ptr rdms); + + /// Compute DSRG full transformed Hamiltonian in the de-normal-ordered basis + std::pair> compute_Heff_full_degno(); + + /// Compute DSRG transformed dipole integral + virtual void compute_mbar() { + throw std::runtime_error("Child class should override this function"); + } + std::vector compute_Mbar0_full() { return Mbar0_full_; } + std::vector compute_Mbar1_full() { return Mbar1_full_; } + std::vector compute_Mbar2_full() { return Mbar2_full_; } + + ambit::BlockedTensor get_gamma1() { return Gamma1_; } + ambit::BlockedTensor get_eta1() { return Eta1_; } + ambit::BlockedTensor get_lambda2() { return Lambda2_; } + std::vector get_lambda3() { return {L3aaa_, L3aab_, L3abb_, L3bbb_}; } + std::vector get_lambda4() { + return {L4aaaa_, L4aaab_, L4aabb_, L4abbb_, L4bbbb_}; + } + /// De-normal-order DSRG transformed dipole moment std::vector deGNO_DMbar_actv(); @@ -278,8 +309,15 @@ class MASTER_DSRG : public DynamicCorrelationSolver { ambit::Tensor L3aab_; ambit::Tensor L3abb_; ambit::Tensor L3bbb_; + /// Four-body density cumulants + ambit::Tensor L4aaaa_; + ambit::Tensor L4aaab_; + ambit::Tensor L4aabb_; + ambit::Tensor L4abbb_; + ambit::Tensor L4bbbb_; /// Whether to compute 3-cumulant contributions bool do_cu3_; + bool do_cu4_; // ==> Fock matrix related <== @@ -327,6 +365,18 @@ class MASTER_DSRG : public DynamicCorrelationSolver { void deGNO_ints(const std::string& name, double& H0, BlockedTensor& H1, BlockedTensor& H2, BlockedTensor& H3); + /** + * De-normal-order the full DSRG transformed Hamiltonian + * This will change H0, H1 and H2!!! + */ + void deGNO_ints_full(const std::string& name, double& H0, BlockedTensor& H1, BlockedTensor& H2); + + /** + * Re-Normal-order the full DSRG transformed Hamiltonian + * This will change H0, H1 and H2!!! + */ + void GNO_ints_full(const std::string& name, double& H0, BlockedTensor& H1, BlockedTensor& H2, const std::shared_ptr rdms); + /** * De-normal-order the T1 and T2 amplitudes and return the effective T1 * T1eff = T1 - T2["ivau"] * D1["uv"] @@ -376,6 +426,10 @@ class MASTER_DSRG : public DynamicCorrelationSolver { /// DSRG transformed 3-body dipole integrals (active only) std::array Mbar3_; + std::vector Mbar0_full_; /// This is redundant. + std::vector Mbar1_full_; + std::vector Mbar2_full_; + // ==> commutators <== /** * H1, C1, G1: a rank 2 tensor of all MOs in general diff --git a/forte/mrdsrg-spin-integrated/mrdsrg.h b/forte/mrdsrg-spin-integrated/mrdsrg.h index 8000c675f..96400ec99 100644 --- a/forte/mrdsrg-spin-integrated/mrdsrg.h +++ b/forte/mrdsrg-spin-integrated/mrdsrg.h @@ -260,6 +260,10 @@ class MRDSRG : public MASTER_DSRG { double get_adaptive_rsc_conv(const int& iter, const double& deltaE, const double& rsc_conv, const double& rsc_conv_adapt, const double& rsc_conv_adapt_delta_e, const double& e_conv); + void compute_mbar() override; + void transform_one_body(const std::vector& oetens); + void compute_mbar_ldsrg2(const ambit::BlockedTensor& M, int ind); + /// Temporary one-body Hamiltonian ambit::BlockedTensor O1_; ambit::BlockedTensor C1_; @@ -350,6 +354,9 @@ class MRDSRG : public MASTER_DSRG { /// Compute MR-LDSRG(2) double compute_energy_ldsrg2(); + /// Modify the rdms_ member of the parent (base) class. Only used in EOM-DSRG. + void set_rdms(std::shared_ptr rdms) override; + /// Zeroth-order Hamiltonian ambit::BlockedTensor H0th_; /// DSRG-MRPT2 zeroth-order Hamiltonian diff --git a/forte/mrdsrg-spin-integrated/mrdsrg_nonpt.cc b/forte/mrdsrg-spin-integrated/mrdsrg_nonpt.cc index 10d594fd3..44fd676c8 100644 --- a/forte/mrdsrg-spin-integrated/mrdsrg_nonpt.cc +++ b/forte/mrdsrg-spin-integrated/mrdsrg_nonpt.cc @@ -42,6 +42,7 @@ #include "helpers/timer.h" #include "base_classes/mo_space_info.h" +#include "integrals/one_body_integrals.h" #include "mrdsrg.h" using namespace psi; @@ -830,6 +831,26 @@ double MRDSRG::compute_energy_ldsrg2() { // dump amplitudes to file dump_amps_to_disk(); + + if (foptions_->get_bool("FULL_HBAR")) { + double rsc_conv = rsc_thres; + compute_hbar(rsc_conv); + } + + if (foptions_->get_bool("FULL_MBAR")){ + auto mpints = std::make_shared(ints_, mo_space_info_); + // bare dipoles + std::vector M1; + std::vector dp_dirs{"X", "Y", "Z"}; + for (int z = 0; z < 3; ++z) { + std::string name = "DIPOLE " + dp_dirs[z]; + ambit::BlockedTensor m1 = BTF_->build(tensor_type_, name, spin_cases({"gg"})); + m1.iterate([&](const std::vector& i, const std::vector&, + double& value) { value = mpints->dp_ints_corr(z, i[0], i[1]); }); + M1.push_back(m1); + } + transform_one_body(M1); + } final.stop(); @@ -837,6 +858,177 @@ double MRDSRG::compute_energy_ldsrg2() { return Ecorr; } +void MRDSRG::compute_mbar() { + auto mpints = std::make_shared(ints_, mo_space_info_); + // bare dipoles + std::vector M1; + std::vector dp_dirs{"X", "Y", "Z"}; + for (int z = 0; z < 3; ++z) { + std::string name = "DIPOLE " + dp_dirs[z]; + ambit::BlockedTensor m1 = BTF_->build(tensor_type_, name, spin_cases({"gg"})); + m1.iterate([&](const std::vector& i, const std::vector&, + double& value) { value = mpints->dp_ints_corr(z, i[0], i[1]); }); + M1.push_back(m1); + } + transform_one_body(M1); +} + +void MRDSRG::transform_one_body(const std::vector& oetens) { + int n_tensors = oetens.size(); + Mbar0_full_ = std::vector(n_tensors, 0.0); + Mbar1_full_.resize(n_tensors); + Mbar2_full_.resize(n_tensors); + for (int i = 0; i < n_tensors; ++i) { + Mbar1_full_[i] = BTF_->build(tensor_type_, oetens[i].name() + "1", spin_cases({"gg"})); + Mbar2_full_[i] = BTF_->build(tensor_type_, oetens[i].name() + "2", spin_cases({"gggg"})); + const auto& M = oetens[i]; + compute_mbar_ldsrg2(M, i); + } +} + +void MRDSRG::compute_mbar_ldsrg2(const ambit::BlockedTensor& M, int ind) { + // compute reference multipole + { + Mbar0_full_[ind] = M["uv"] * Gamma1_["vu"]; + Mbar0_full_[ind] += M["UV"] * Gamma1_["VU"]; + auto& M1c = M.block("cc").data(); + for (size_t m = 0, ncore = core_mos_.size(); m < ncore; ++m) { + Mbar0_full_[ind] += M1c[m * ncore + m]; + } + auto& M1C = M.block("CC").data(); + for (size_t m = 0, ncore = core_mos_.size(); m < ncore; ++m) { + Mbar0_full_[ind] += M1C[m * ncore + m]; + } + Mbar1_full_[ind]["pq"] = M["pq"]; + Mbar1_full_[ind]["PQ"] = M["PQ"]; + } + + O1_["pq"] = M["pq"]; + O1_["PQ"] = M["PQ"]; + + bool converged = false; + int maxn = foptions_->get_int("DSRG_RSC_NCOMM"); + double ct_threshold = foptions_->get_double("DSRG_RSC_THRESHOLD"); + std::string dsrg_op = foptions_->get_str("DSRG_TRANS_TYPE"); + + // compute Hbar recursively + for (int n = 1; n <= maxn; ++n) { + // prefactor before n-nested commutator + double factor = 1.0 / n; + + // Compute the commutator C = 1/n [O, T] + double C0 = 0.0; + C1_.zero(); + C2_.zero(); + + // zero-body + H1_T1_C0(O1_, T1_, factor, C0); + H1_T2_C0(O1_, T2_, factor, C0); + if (n != 1) { + H2_T1_C0(O2_, T1_, factor, C0); + H2_T2_C0(O2_, T2_, factor, C0); + } + // one-body + H1_T1_C1(O1_, T1_, factor, C1_); + H1_T2_C1(O1_, T2_, factor, C1_); + if (n != 1) { + H2_T1_C1(O2_, T1_, factor, C1_); + } + if (foptions_->get_str("SRG_COMM") == "STANDARD") { + if (n != 1) { + H2_T2_C1(O2_, T2_, factor, C1_); + } + } else if (foptions_->get_str("SRG_COMM") == "FO") { + BlockedTensor C1p = BTF_->build(tensor_type_, "C1p", spin_cases({"gg"})); + if (n != 1) { + H2_T2_C1(O2_, T2_, factor, C1p); + } + C1p.block("cc").scale(2.0); + C1p.block("aa").scale(2.0); + C1p.block("vv").scale(2.0); + C1p.block("CC").scale(2.0); + C1p.block("AA").scale(2.0); + C1p.block("VV").scale(2.0); + C1_["pq"] += C1p["pq"]; + C1_["PQ"] += C1p["PQ"]; + } + // two-body + if ((foptions_->get_str("SRG_COMM") == "STANDARD") or n < 2) { + H1_T2_C2(O1_, T2_, factor, C2_); + } else if (foptions_->get_str("SRG_COMM") == "FO2") { + O1_.block("cc").scale(2.0); + O1_.block("aa").scale(2.0); + O1_.block("vv").scale(2.0); + O1_.block("CC").scale(2.0); + O1_.block("AA").scale(2.0); + O1_.block("VV").scale(2.0); + H1_T2_C2(O1_, T2_, factor, C2_); + O1_.block("cc").scale(0.5); + O1_.block("aa").scale(0.5); + O1_.block("vv").scale(0.5); + O1_.block("CC").scale(0.5); + O1_.block("AA").scale(0.5); + O1_.block("VV").scale(0.5); + } + if (n != 1) { + H2_T1_C2(O2_, T1_, factor, C2_); + H2_T2_C2(O2_, T2_, factor, C2_); + } + + // printing level + if (print_ > 3) { + std::string dash(38, '-'); + outfile->Printf("\n %s\n", dash.c_str()); + } + + // [H, A] = [H, T] + [H, T]^dagger + if (dsrg_op == "UNITARY") { + C0 *= 2.0; + O1_["pq"] = C1_["pq"]; + O1_["PQ"] = C1_["PQ"]; + C1_["pq"] += O1_["qp"]; + C1_["PQ"] += O1_["QP"]; + O2_["pqrs"] = C2_["pqrs"]; + O2_["pQrS"] = C2_["pQrS"]; + O2_["PQRS"] = C2_["PQRS"]; + C2_["pqrs"] += O2_["rspq"]; + C2_["pQrS"] += O2_["rSpQ"]; + C2_["PQRS"] += O2_["RSPQ"]; + } + + // Hbar += C + Mbar0_full_[ind] += C0; + Mbar1_full_[ind]["pq"] += C1_["pq"]; + Mbar1_full_[ind]["PQ"] += C1_["PQ"]; + Mbar2_full_[ind]["pqrs"] += C2_["pqrs"]; + Mbar2_full_[ind]["pQrS"] += C2_["pQrS"]; + Mbar2_full_[ind]["PQRS"] += C2_["PQRS"]; + + // copy C to O for next level commutator + O1_["pq"] = C1_["pq"]; + O1_["PQ"] = C1_["PQ"]; + O2_["pqrs"] = C2_["pqrs"]; + O2_["pQrS"] = C2_["pQrS"]; + O2_["PQRS"] = C2_["PQRS"]; + + // test convergence of C + double norm_C1 = C1_.norm(); + double norm_C2 = C2_.norm(); + if (print_ > 3) { + outfile->Printf("\n n: %3d, C0: %20.15f, C1 max: %20.15f, C2 max: %20.15f", n, C0, + C1_.norm(0), C2_.norm(0)); + } + if (std::sqrt(norm_C2 * norm_C2 + norm_C1 * norm_C1) < ct_threshold) { + converged = true; + break; + } + } + if (!converged) { + outfile->Printf("\n Warning! Mbar is not converged in %3d-nested commutators!", maxn); + outfile->Printf("\n Please increase DSRG_RSC_NCOMM."); + } +} + void MRDSRG::compute_hbar_qc() { // initialize Hbar with bare H Hbar0_ = 0.0; @@ -1271,4 +1463,12 @@ double MRDSRG::get_adaptive_rsc_conv(const int& iter, const double& deltaE, cons return threshold; } +void MRDSRG::set_rdms(std::shared_ptr rdms) { + // set up RDMS + rdms_ = rdms; + build_density(); + // build_ints(); + // build_fock(H_, V_); +} + } // namespace forte diff --git a/forte/options.yaml b/forte/options.yaml index 1e2d0d916..47cb1c05d 100644 --- a/forte/options.yaml +++ b/forte/options.yaml @@ -1114,6 +1114,14 @@ DSRG: " - Multi-state approach (currently only for MRPT2)\n" " - MS: form 2nd-order Heff_MN = + 0.5 * [ + ]\n" " - XMS: rotate references such that is diagonal before MS procedure" + FULL_HBAR: + type: bool + default: False + help: "Save the full DSRG effective Hamiltonian to disk." + FULL_MBAR: + type: bool + default: False + help: "Save the full DSRG effective dipole integrals to disk." FORM_HBAR3: type: bool default: False @@ -1570,6 +1578,11 @@ Old: default: "MK" choices: ["MK", "MK_DECOMP", "ZERO"] help: "The form of the three-particle density cumulant" + FOURPDC: + type: str + default: "ZERO" + choices: ["MK", "ZERO"] + help: "The form of the four-particle density cumulant" SRG_COMM: type: str default: "STANDARD" diff --git a/forte/orbital-helpers/semi_canonicalize.h b/forte/orbital-helpers/semi_canonicalize.h index 76489f069..90b58c6d9 100644 --- a/forte/orbital-helpers/semi_canonicalize.h +++ b/forte/orbital-helpers/semi_canonicalize.h @@ -158,6 +158,7 @@ class SemiCanonical { std::shared_ptr Ua_; /// Unitary matrix for beta orbital rotation std::shared_ptr Ub_; + /// Unitary matrix for alpha orbital rotation in the active space ambit::Tensor Ua_t_; /// Unitary matrix for beta orbital rotation in the active space diff --git a/forte/proc/dsrg.py b/forte/proc/dsrg.py index 81302ec5e..efa386ad3 100644 --- a/forte/proc/dsrg.py +++ b/forte/proc/dsrg.py @@ -69,7 +69,13 @@ def __init__(self, active_space_solver, state_weights_map, mo_space_info, ints, psi4.core.print_out("\n DSRG relaxation only supports UNITARY transformation. Setting RELAX_REF to NONE.") self.relax_ref = "NONE" - self.max_rdm_level = 3 if options.get_str("THREEPDC") != "ZERO" else 2 + if options.get_str("FOURPDC") != "ZERO": + self.max_rdm_level = 4 + elif options.get_str("THREEPDC") != "ZERO": + self.max_rdm_level = 3 + else: + self.max_rdm_level = 2 + if options.get_str("DSRG_3RDM_ALGORITHM") == "DIRECT": as_type = options.get_str("ACTIVE_SPACE_SOLVER") if as_type == "BLOCK2" and self.solver_type in ["SA-MRDSRG", "SA_MRDSRG"]: @@ -242,6 +248,27 @@ def compute_energy(self): # Spit out energy if reference relaxation not implemented if not self.Heff_implemented: self.relax_maxiter = 0 + + if self.options.get_bool("FULL_HBAR") and self.solver_type in ["MRDSRG","SA-MRDSRG","SA_MRDSRG","MRDSRG_SO", "MRDSRG-SO"] \ + and self.relax_maxiter == 0: + psi4.core.print_out("\n =>** Saving Full Hbar (unrelaxed) **<=\n") + Hbar0, Hbar1, Hbar2 = self.dsrg_solver.save_Heff_full() + psi4.core.print_out(f"\n The Hbar0 term is: {Hbar0}\n") + np.savez("save_Hbar", Hbar0=Hbar0, **Hbar1, **Hbar2) + + if self.options.get_bool("FULL_MBAR") and self.solver_type in ["MRDSRG"] and self.relax_maxiter == 0: + psi4.core.print_out("\n =>** Getting dipole integral (unrelaxed) **<=\n") + self.dsrg_solver.compute_mbar() + Mbar0 = self.dsrg_solver.compute_Mbar0_full() + np.save("Mbar0", Mbar0) + Mbar1 = self.dsrg_solver.compute_Mbar1_full() + Mbar2 = self.dsrg_solver.compute_Mbar2_full() + + for i in range(3): + np.savez(f"Mbar1_{i}", **Mbar1[i]) + np.savez(f"Mbar2_{i}", **Mbar2[i]) + + del Mbar0, Mbar1, Mbar2 # Reference relaxation procedure for n in range(self.relax_maxiter): @@ -252,6 +279,14 @@ def compute_energy(self): # so that the CI vectors are comparable before and after DSRG dressing. # However, the ForteIntegrals object and the dipole integrals always refer to the current semi-canonical basis. # so to compute the dipole moment correctly, we need to make the RDMs and orbital basis consistent + if self.options.get_bool("FULL_HBAR") and self.solver_type in ["SA-MRDSRG","SA_MRDSRG","MRDSRG_SO","MRDSRG-SO"]: + raise NotImplementedError("Relaxed full Hbar not implemented for SO/SA-DSRG") + + if self.options.get_bool("FULL_HBAR") and self.solver_type in ["MRDSRG"] and n == self.relax_maxiter - 1: + psi4.core.print_out("\n =>** Temporarily saving Full Hbar in de-normal-ordered basis (relaxed) **<=\n") + Hbar0, Hbar1, Hbar2 = self.dsrg_solver.save_Heff_full_ambit() + psi4.core.print_out(f"\n The Hbar0 term is: {Hbar0}\n") + ints_dressed = self.dsrg_solver.compute_Heff_actv() if self.fno_pt2_Heff_shift is not None: ints_dressed.add(self.fno_pt2_Heff_shift, 1.0) @@ -374,12 +409,63 @@ def compute_energy(self): # - Compute the DSRG energy self.make_dsrg_solver() self.dsrg_setup() - self.dsrg_solver.set_read_cwd_amps(not self.restart_amps) # don't read from cwd if checkpoint available + # don't read from cwd if checkpoint available + self.dsrg_solver.set_read_cwd_amps(not self.restart_amps) e_dsrg = self.dsrg_solver.compute_energy() + self.fno_pt2_energy_shift if self.fno_pt2_energy_shift != 0.0: psi4.core.print_out(f"\n\n DSRG-MRPT2 FNO energy correction: {self.fno_pt2_energy_shift:20.15f}") psi4.core.print_out(f"\n DSRG-MRPT2 FNO corrected energy: {e_dsrg:20.15f}") + if self.options.get_bool("FULL_HBAR") and self.solver_type in ["MRDSRG","MRDSRG_SO","MRDSRG-SO"]: + if self.relax_maxiter != 0: + self.rdms = self.active_space_solver.compute_average_rdms(self.state_weights_map, self.max_rdm_level, self.rdm_type) + self.rdms.rotate(self.Ua, self.Ub) # To previous semi-canonical basis + + psi4.core.print_out("\n =>** Saving Full Hbar **<=\n") + Hbar0, Hbar1, Hbar2 = self.dsrg_solver.update_Heff_full(Hbar0, Hbar1, Hbar2, self.rdms) + psi4.core.print_out(f"\n The Hbar0 term is: {Hbar0}\n") + np.savez("save_Hbar", Hbar0=Hbar0, **Hbar1, **Hbar2) + + self.dsrg_solver.set_rdms(self.rdms) # Important to update the rdms in the solver + + if self.options.get_bool("FULL_MBAR") and self.solver_type in ["MRDSRG"]: + psi4.core.print_out("\n =>** Getting dipole integral **<=\n") + self.dsrg_solver.compute_mbar() + Mbar0 = self.dsrg_solver.compute_Mbar0_full() + np.save("Mbar0", Mbar0) + Mbar1 = self.dsrg_solver.compute_Mbar1_full() + Mbar2 = self.dsrg_solver.compute_Mbar2_full() + + for i in range(3): + np.savez(f"Mbar1_{i}", **Mbar1[i]) + np.savez(f"Mbar2_{i}", **Mbar2[i]) + + del Mbar0, Mbar1, Mbar2 + + psi4.core.print_out("\n =>** Getting gamma1 **<=\n") + gamma1 = self.dsrg_solver.get_gamma1() + np.savez("save_gamma1", **gamma1) + + psi4.core.print_out("\n =>** Getting eta1 **<=\n") + eta1 = self.dsrg_solver.get_eta1() + np.savez("save_eta1", **eta1) + + psi4.core.print_out("\n =>** Getting lambda2 **<=\n") + lambda2 = self.dsrg_solver.get_lambda2() + np.savez("save_lambda2", **lambda2) + + psi4.core.print_out("\n =>** Getting lambda3 **<=\n") + lambda3 = self.dsrg_solver.get_lambda3() + np.savez("save_lambda3", **lambda3) + + del gamma1, eta1, lambda2, lambda3 + + if self.options.get_str("FOURPDC") != "ZERO" and self.solver_type in ["MRDSRG"]: + psi4.core.print_out("\n =>** Getting lambda4 **<=\n") + lambda4 = self.dsrg_solver.get_lambda4() + np.savez("save_lambda4", **lambda4) + del lambda4 + self.dsrg_cleanup() # dump reference relaxation energies to json file diff --git a/forte/pymodule.py b/forte/pymodule.py index 280b7a3d8..fdd65273b 100644 --- a/forte/pymodule.py +++ b/forte/pymodule.py @@ -108,13 +108,15 @@ def forte_driver(data: ForteData): write_external_rdm_file(data.rdms, data.active_space_solver) if options.get_bool("SPIN_ANALYSIS"): - data = ActiveSpaceRDMs(max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) + data = ActiveSpaceRDMs( + max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) forte.perform_spin_analysis(data.rdms, options, mo_space_info, as_ints) # solver for dynamical correlation from DSRG correlation_solver_type = options.get_str("CORRELATION_SOLVER") if correlation_solver_type != "NONE": - dsrg_proc = ProcedureDSRG(data.active_space_solver, state_weights_map, mo_space_info, ints, options, scf_info) + dsrg_proc = ProcedureDSRG( + data.active_space_solver, state_weights_map, mo_space_info, ints, options, scf_info) return_en = dsrg_proc.compute_energy() dsrg_proc.print_summary() dsrg_proc.push_to_psi4_environment() @@ -125,13 +127,15 @@ def forte_driver(data: ForteData): # DSRG reads consistent CI coefficients before and after SemiCanonical class. # 2. This is OK only when running ground-state calculations state = list(state_map.keys())[0] - psi4.core.print_out(f"\n ==> Coupling Coefficients for {state} <==") + psi4.core.print_out( + f"\n ==> Coupling Coefficients for {state} <==") ci_vectors = data.active_space_solver.eigenvectors(state) dsrg_proc.compute_gradient(ci_vectors) else: psi4.core.print_out("\n Semicanonical orbitals must be used!\n") else: - average_energy = forte.compute_average_state_energy(state_energies_list, state_weights_map) + average_energy = forte.compute_average_state_energy( + state_energies_list, state_weights_map) return_en = average_energy return return_en @@ -187,7 +191,8 @@ def energy_forte(name, **kwargs): # if data.options.get_str("INT_TYPE") == "FCIDUMP": # psi4.core.print_out("\n\n Skipping MCSCF computation. Using integrals from FCIDUMP input\n") if data.options.get_bool("MCSCF_REFERENCE") is False: - psi4.core.print_out("\n\n Skipping MCSCF computation. Using HF or orbitals passed via ref_wfn\n") + psi4.core.print_out( + "\n\n Skipping MCSCF computation. Using HF or orbitals passed via ref_wfn\n") else: active_space_solver_type = data.options.get_str("ACTIVE_SPACE_SOLVER") mcscf_ignore_frozen = data.options.get_bool("MCSCF_IGNORE_FROZEN_ORBS") @@ -236,9 +241,12 @@ def energy_forte(name, **kwargs): psi4.core.set_scalar_variable("CURRENT ENERGY", energy) - psi4.core.print_out(f"\n\n Time to prepare integrals: {start - start_pre_ints:12.3f} seconds") - psi4.core.print_out(f"\n Time to run job : {end - start:12.3f} seconds") - psi4.core.print_out(f"\n Total : {end - start_pre_ints:12.3f} seconds\n") + psi4.core.print_out( + f"\n\n Time to prepare integrals: {start - start_pre_ints:12.3f} seconds") + psi4.core.print_out( + f"\n Time to run job : {end - start:12.3f} seconds") + psi4.core.print_out( + f"\n Total : {end - start_pre_ints:12.3f} seconds\n") if "FCIDUMP" not in data.options.get_str("INT_TYPE"): if data.options.get_bool("DUMP_ORBITALS"): @@ -332,8 +340,10 @@ def gradient_forte(name, **kwargs): ] max_key_size = max(len(k) for k, v in times) for key, value in times: - psi4.core.print_out(f"\n Time to {key:{max_key_size}} : {value:12.3f} seconds") - psi4.core.print_out(f'\n {"Total":{max_key_size + 8}} : {end - time_pre_ints:12.3f} seconds\n') + psi4.core.print_out( + f"\n Time to {key:{max_key_size}} : {value:12.3f} seconds") + psi4.core.print_out( + f'\n {"Total":{max_key_size + 8}} : {end - time_pre_ints:12.3f} seconds\n') # Dump orbitals if needed if data.options.get_bool("DUMP_ORBITALS"): diff --git a/forte/utils/helpers.py b/forte/utils/helpers.py index c88d642f3..e19753b0d 100644 --- a/forte/utils/helpers.py +++ b/forte/utils/helpers.py @@ -184,7 +184,7 @@ def psi4_cubeprop(wfn, path=".", orbs=[], nocc=0, nvir=0, density=False, frontie def prepare_forte_objects( - wfn, mo_spaces=None, active_space="ACTIVE", core_spaces=["RESTRICTED_DOCC"], localize=False, localize_spaces=[] + wfn, mo_spaces, active_space="ACTIVE", core_spaces=["RESTRICTED_DOCC"], localize=False, localize_spaces=[] ): """Take a psi4 wavefunction object and prepare the ForteIntegrals, SCFInfo, and MOSpaceInfo objects @@ -235,15 +235,18 @@ def prepare_forte_objects( point_group = wfn.molecule().point_group().symbol() # create a MOSpaceInfo object - if mo_spaces is None: - mo_space_info = forte.make_mo_space_info(nmopi, point_group, options) - else: - mo_space_info = forte.make_mo_space_info_from_map(nmopi, point_group, mo_spaces) + mo_space_info = forte.make_mo_space_info_from_map(nmopi, point_group, mo_spaces) + + # These variables are needed in make_state_weights_map + nel = wfn.nalpha() + wfn.nbeta() + multiplicity = 1 + options.set_int("NEL", nel) + options.set_int("MULTIPLICITY", multiplicity) state_weights_map = forte.make_state_weights_map(options, mo_space_info) # make a ForteIntegral object - ints = forte.make_ints_from_psi4(wfn, options, mo_space_info) + ints = forte.make_ints_from_psi4(wfn, options, scf_info, mo_space_info) if localize: localizer = forte.Localize(forte.forte_options, ints, mo_space_info) diff --git a/tests/methods/fci-rdms-1/output.ref b/tests/methods/fci-rdms-1/output.ref new file mode 100644 index 000000000..7af108138 --- /dev/null +++ b/tests/methods/fci-rdms-1/output.ref @@ -0,0 +1,928 @@ + + ----------------------------------------------------------------------- + Psi4: An Open-Source Ab Initio Electronic Structure Package + Psi4 1.10a1.dev53 + + Git: Rev {master} 45660d8 dirty + + + D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, + M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, + A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, + R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, + J. S. O'Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, + B. P. Pritchard, B. R. Brooks, H. F. Schaefer III, A. Yu. Sokolov, + K. Patkowski, A. E. DePrince III, U. Bozkaya, R. A. King, + F. A. Evangelista, J. M. Turney, T. D. Crawford, C. D. Sherrill, + J. Chem. Phys. 152(18) 184108 (2020). https://doi.org/10.1063/5.0006002 + + Additional Code Authors + E. T. Seidl, C. L. Janssen, E. F. Valeev, M. L. Leininger, + J. F. Gonthier, R. M. Richard, H. R. McAlexander, M. Saitow, X. Wang, + P. Verma, M. H. Lechner, A. Jiang, S. Behnle, A. G. Heide, + M. F. Herbst, and D. L. Poole + + Previous Authors, Complete List of Code Contributors, + and Citations for Specific Modules + https://github.com/psi4/psi4/blob/master/codemeta.json + https://github.com/psi4/psi4/graphs/contributors + http://psicode.org/psi4manual/master/introduction.html#citing-psifour + + ----------------------------------------------------------------------- + + + Psi4 started on: Tuesday, 17 September 2024 03:25PM + + Process ID: 19158 + Host: Brian-Zs-MBA.local + PSIDATADIR: /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4 + Memory: 500.0 MiB + Threads: 1 + + ==> Input File <== + +-------------------------------------------------------------------------- +#! Generated using commit GITCOMMIT + +import forte + +molecule { +0 1 +Li +H 1 R + +R = 3.0 +units bohr +} + +set { + basis sto-3g + reference rhf + scf_type pk + e_convergence 12 +} + +set forte { + active_space_solver fci + fci_test_rdms true +} + +energy('scf') +energy('forte') + +compare_values(0.0, variable("AA 1-RDM ERROR"),12, "AA 1-RDM") #TEST +compare_values(0.0, variable("BB 1-RDM ERROR"),12, "BB 1-RDM") #TEST +compare_values(0.0, variable("AAAA 2-RDM ERROR"),12, "AAAA 2-RDM") #TEST +compare_values(0.0, variable("BBBB 2-RDM ERROR"),12, "BBBB 2-RDM") #TEST +compare_values(0.0, variable("ABAB 2-RDM ERROR"),12, "ABAB 2-RDM") #TEST +compare_values(0.0, variable("AABAAB 3-RDM ERROR"),12, "AABAAB 3-RDM") #TEST +compare_values(0.0, variable("ABBABB 3-RDM ERROR"),12, "ABBABB 3-RDM") #TEST +compare_values(0.0, variable("AAAAAA 3-RDM ERROR"),12, "AAAAAA 3-RDM") #TEST +compare_values(0.0, variable("BBBBBB 3-RDM ERROR"),12, "BBBBBB 3-RDM") #TEST +-------------------------------------------------------------------------- + +Scratch directory: /tmp/ + => Libint2 <= + + Primary basis highest AM E, G, H: 6, 6, 3 + Auxiliary basis highest AM E, G, H: 7, 7, 4 + Onebody basis highest AM E, G, H: -, -, - + Solid Harmonics ordering: Gaussian + +*** tstart() called on Brian-Zs-MBA.local +*** at Tue Sep 17 15:25:11 2024 + + => Loading Basis Set <= + + Name: STO-3G + Role: ORBITAL + Keyword: BASIS + atoms 1 entry LI line 31 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + atoms 2 entry H line 19 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + + + --------------------------------------------------------- + SCF + by Justin Turney, Rob Parrish, Andy Simmonett + and Daniel G. A. Smith + RHF Reference + 1 Threads, 500 MiB Core + --------------------------------------------------------- + + ==> Geometry <== + + Molecular point group: c2v + Full point group: C_inf_v + + Geometry (in Bohr), charge = 0, multiplicity = 1: + + Center X Y Z Mass + ------------ ----------------- ----------------- ----------------- ----------------- + LI 0.000000000000 0.000000000000 -0.376812030371 7.016003436600 + H 0.000000000000 0.000000000000 2.623187969629 1.007825032230 + + Running in c2v symmetry. + + Rotational constants: A = ************ B = 7.59029 C = 7.59029 [cm^-1] + Rotational constants: A = ************ B = 227551.19787 C = 227551.19787 [MHz] + Nuclear repulsion = 1.000000000000000 + + Charge = 0 + Multiplicity = 1 + Electrons = 4 + Nalpha = 2 + Nbeta = 2 + + ==> Algorithm <== + + SCF Algorithm Type is PK. + DIIS enabled. + MOM disabled. + Fractional occupation disabled. + Guess Type is SAD. + Energy threshold = 1.00e-12 + Density threshold = 1.00e-06 + Integral threshold = 1.00e-12 + + ==> Primary Basis <== + + Basis Set: STO-3G + Blend: STO-3G + Number of shells: 4 + Number of basis functions: 6 + Number of Cartesian functions: 6 + Spherical Harmonics?: true + Max angular momentum: 1 + + ==> Integral Setup <== + + Using in-core PK algorithm. + Calculation information: + Number of atoms: 2 + Number of AO shells: 4 + Number of primitives: 12 + Number of atomic orbitals: 6 + Number of basis functions: 6 + + Integral cutoff 1.00e-12 + Number of threads: 1 + + Performing in-core PK + Using 462 doubles for integral storage. + We computed 55 shell quartets total. + Whereas there are 55 unique shell quartets. + + ==> DiskJK: Disk-Based J/K Matrices <== + + J tasked: Yes + K tasked: Yes + wK tasked: No + Memory [MiB]: 375 + Schwarz Cutoff: 1E-12 + + OpenMP threads: 1 + + Minimum eigenvalue in the overlap matrix is 3.4333995519E-01. + Reciprocal condition number of the overlap matrix is 2.0339047710E-01. + Using symmetric orthogonalization. + + ==> Pre-Iterations <== + + SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF (no occupation information). + + ------------------------- + Irrep Nso Nmo + ------------------------- + A1 4 4 + A2 0 0 + B1 1 1 + B2 1 1 + ------------------------- + Total 6 6 + ------------------------- + + ==> Iterations <== + + Total Energy Delta E RMS |[F,P]| + + @RHF iter SAD: -6.75622365412217 -6.75622e+00 0.00000e+00 + @RHF iter 1: -7.86092861588784 -1.10470e+00 8.57305e-03 DIIS/ADIIS + @RHF iter 2: -7.86223475166447 -1.30614e-03 8.22758e-04 DIIS/ADIIS + @RHF iter 3: -7.86224589388798 -1.11422e-05 7.41119e-05 DIIS + @RHF iter 4: -7.86224624208850 -3.48201e-07 2.61264e-05 DIIS + @RHF iter 5: -7.86224631040868 -6.83202e-08 3.24230e-07 DIIS + @RHF iter 6: -7.86224631041032 -1.64135e-12 1.56779e-08 DIIS + @RHF iter 7: -7.86224631041034 -1.77636e-14 2.79325e-12 DIIS + Energy and wave function converged. + + + ==> Post-Iterations <== + + Orbital Energies [Eh] + --------------------- + + Doubly Occupied: + + 1A1 -2.348477 2A1 -0.286330 + + Virtual: + + 3A1 0.078326 1B1 0.163933 1B2 0.163933 + 4A1 0.551172 + + Final Occupation by Irrep: + A1 A2 B1 B2 + DOCC [ 2, 0, 0, 0 ] + NA [ 2, 0, 0, 0 ] + NB [ 2, 0, 0, 0 ] + + @RHF Final Energy: -7.86224631041034 + + => Energetics <= + + Nuclear Repulsion Energy = 1.0000000000000000 + One-Electron Energy = -12.4548780553254126 + Two-Electron Energy = 3.5926317449150740 + Total Energy = -7.8622463104103382 + +Computation Completed + + +Properties will be evaluated at 0.000000, 0.000000, 0.000000 [a0] + +Properties computed using the SCF density matrix + + + Multipole Moments: + + ------------------------------------------------------------------------------------ + Multipole Electronic (a.u.) Nuclear (a.u.) Total (a.u.) + ------------------------------------------------------------------------------------ + + L = 1. Multiply by 2.5417464519 to convert [e a0] to [Debye] + Dipole X : 0.0000000 0.0000000 0.0000000 + Dipole Y : 0.0000000 0.0000000 0.0000000 + Dipole Z : -3.4031201 1.4927519 -1.9103682 + Magnitude : 1.9103682 + + ------------------------------------------------------------------------------------ + +*** tstop() called on Brian-Zs-MBA.local at Tue Sep 17 15:25:11 2024 +Module time: + user time = 0.10 seconds = 0.00 minutes + system time = 0.01 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes +Total time: + user time = 0.10 seconds = 0.00 minutes + system time = 0.01 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes + +Scratch directory: /tmp/ + + Forte + ---------------------------------------------------------------------------- + A suite of quantum chemistry methods for strongly correlated electrons + + git branch: rdm4 - git commit: 7780d650 + + Developed by: + Francesco A. Evangelista, Chenyang Li, Kevin P. Hannon, + Jeffrey B. Schriber, Tianyuan Zhang, Chenxi Cai, + Nan He, Nicholas Stair, Shuhe Wang, Renke Huang + ---------------------------------------------------------------------------- + + + Preparing forte objects from a Psi4 Wavefunction object + No reference wave function provided for Forte. Computing SCF orbitals using Psi4 ... + => Libint2 <= + + Primary basis highest AM E, G, H: 6, 6, 3 + Auxiliary basis highest AM E, G, H: 7, 7, 4 + Onebody basis highest AM E, G, H: -, -, - + Solid Harmonics ordering: Gaussian + +*** tstart() called on Brian-Zs-MBA.local +*** at Tue Sep 17 15:25:11 2024 + + => Loading Basis Set <= + + Name: STO-3G + Role: ORBITAL + Keyword: BASIS + atoms 1 entry LI line 31 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + atoms 2 entry H line 19 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + + + --------------------------------------------------------- + SCF + by Justin Turney, Rob Parrish, Andy Simmonett + and Daniel G. A. Smith + RHF Reference + 1 Threads, 500 MiB Core + --------------------------------------------------------- + + ==> Geometry <== + + Molecular point group: c2v + Full point group: C_inf_v + + Geometry (in Bohr), charge = 0, multiplicity = 1: + + Center X Y Z Mass + ------------ ----------------- ----------------- ----------------- ----------------- + LI 0.000000000000 0.000000000000 -0.376812030371 7.016003436600 + H 0.000000000000 0.000000000000 2.623187969629 1.007825032230 + + Running in c2v symmetry. + + Rotational constants: A = ************ B = 7.59029 C = 7.59029 [cm^-1] + Rotational constants: A = ************ B = 227551.19787 C = 227551.19787 [MHz] + Nuclear repulsion = 1.000000000000000 + + Charge = 0 + Multiplicity = 1 + Electrons = 4 + Nalpha = 2 + Nbeta = 2 + + ==> Algorithm <== + + SCF Algorithm Type is PK. + DIIS enabled. + MOM disabled. + Fractional occupation disabled. + Guess Type is SAD. + Energy threshold = 1.00e-12 + Density threshold = 1.00e-08 + Integral threshold = 1.00e-12 + + ==> Primary Basis <== + + Basis Set: STO-3G + Blend: STO-3G + Number of shells: 4 + Number of basis functions: 6 + Number of Cartesian functions: 6 + Spherical Harmonics?: true + Max angular momentum: 1 + + ==> Integral Setup <== + + Using in-core PK algorithm. + Calculation information: + Number of atoms: 2 + Number of AO shells: 4 + Number of primitives: 12 + Number of atomic orbitals: 6 + Number of basis functions: 6 + + Integral cutoff 1.00e-12 + Number of threads: 1 + + Performing in-core PK + Using 462 doubles for integral storage. + We computed 55 shell quartets total. + Whereas there are 55 unique shell quartets. + + ==> DiskJK: Disk-Based J/K Matrices <== + + J tasked: Yes + K tasked: Yes + wK tasked: No + Memory [MiB]: 375 + Schwarz Cutoff: 1E-12 + + OpenMP threads: 1 + + Minimum eigenvalue in the overlap matrix is 3.4333995519E-01. + Reciprocal condition number of the overlap matrix is 2.0339047710E-01. + Using symmetric orthogonalization. + + ==> Pre-Iterations <== + + SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF (no occupation information). + + ------------------------- + Irrep Nso Nmo + ------------------------- + A1 4 4 + A2 0 0 + B1 1 1 + B2 1 1 + ------------------------- + Total 6 6 + ------------------------- + + ==> Iterations <== + + Total Energy Delta E RMS |[F,P]| + + @RHF iter SAD: -6.75622365412217 -6.75622e+00 0.00000e+00 + @RHF iter 1: -7.86092861588784 -1.10470e+00 8.57305e-03 DIIS/ADIIS + @RHF iter 2: -7.86223475166447 -1.30614e-03 8.22758e-04 DIIS/ADIIS + @RHF iter 3: -7.86224589388798 -1.11422e-05 7.41119e-05 DIIS + @RHF iter 4: -7.86224624208850 -3.48201e-07 2.61264e-05 DIIS + @RHF iter 5: -7.86224631040868 -6.83202e-08 3.24230e-07 DIIS + @RHF iter 6: -7.86224631041032 -1.64135e-12 1.56779e-08 DIIS + @RHF iter 7: -7.86224631041034 -1.77636e-14 2.79325e-12 DIIS + Energy and wave function converged. + + + ==> Post-Iterations <== + + Orbital Energies [Eh] + --------------------- + + Doubly Occupied: + + 1A1 -2.348477 2A1 -0.286330 + + Virtual: + + 3A1 0.078326 1B1 0.163933 1B2 0.163933 + 4A1 0.551172 + + Final Occupation by Irrep: + A1 A2 B1 B2 + DOCC [ 2, 0, 0, 0 ] + NA [ 2, 0, 0, 0 ] + NB [ 2, 0, 0, 0 ] + + @RHF Final Energy: -7.86224631041034 + + => Energetics <= + + Nuclear Repulsion Energy = 1.0000000000000000 + One-Electron Energy = -12.4548780553254126 + Two-Electron Energy = 3.5926317449150740 + Total Energy = -7.8622463104103382 + +Computation Completed + + +Properties will be evaluated at 0.000000, 0.000000, 0.000000 [a0] + +Properties computed using the SCF density matrix + + + Multipole Moments: + + ------------------------------------------------------------------------------------ + Multipole Electronic (a.u.) Nuclear (a.u.) Total (a.u.) + ------------------------------------------------------------------------------------ + + L = 1. Multiply by 2.5417464519 to convert [e a0] to [Debye] + Dipole X : 0.0000000 0.0000000 0.0000000 + Dipole Y : 0.0000000 0.0000000 0.0000000 + Dipole Z : -3.4031201 1.4927519 -1.9103682 + Magnitude : 1.9103682 + + ------------------------------------------------------------------------------------ + +*** tstop() called on Brian-Zs-MBA.local at Tue Sep 17 15:25:11 2024 +Module time: + user time = 0.08 seconds = 0.00 minutes + system time = 0.00 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes +Total time: + user time = 0.31 seconds = 0.01 minutes + system time = 0.01 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes + + + ==> MO Space Information <== + + ------------------------------------------------- + A1 A2 B1 B2 Sum + ------------------------------------------------- + FROZEN_DOCC 0 0 0 0 0 + RESTRICTED_DOCC 0 0 0 0 0 + GAS1 4 0 1 1 6 + GAS2 0 0 0 0 0 + GAS3 0 0 0 0 0 + GAS4 0 0 0 0 0 + GAS5 0 0 0 0 0 + GAS6 0 0 0 0 0 + RESTRICTED_UOCC 0 0 0 0 0 + FROZEN_UOCC 0 0 0 0 0 + Total 4 0 1 1 6 + ------------------------------------------------- => Loading Basis Set <= + + Name: STO-3G + Role: ORBITAL + Keyword: MINAO_BASIS + atoms 1 entry LI line 31 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + atoms 2 entry H line 19 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + + + State Singlet (Ms = 0) A1 GAS min: 0 0 0 0 0 0 ; GAS max: 12 0 0 0 0 0 ; weights: + 1.000000000000 + Forte will use psi4 integrals + + ==> Primary Basis Set Summary <== + + Basis Set: STO-3G + Blend: STO-3G + Number of shells: 4 + Number of basis functions: 6 + Number of Cartesian functions: 6 + Spherical Harmonics?: true + Max angular momentum: 1 + + + JK created using conventional PK integrals + Using in-core PK algorithm. + Calculation information: + Number of atoms: 2 + Number of AO shells: 4 + Number of primitives: 12 + Number of atomic orbitals: 6 + Number of basis functions: 6 + + Integral cutoff 1.00e-12 + Number of threads: 1 + + Performing in-core PK + Using 462 doubles for integral storage. + We computed 55 shell quartets total. + Whereas there are 55 unique shell quartets. + + ==> DiskJK: Disk-Based J/K Matrices <== + + J tasked: Yes + K tasked: Yes + wK tasked: No + Memory [MiB]: 400 + Schwarz Cutoff: 1E-12 + + OpenMP threads: 1 + + + + ==> Integral Transformation <== + + Number of molecular orbitals: 6 + Number of correlated molecular orbitals: 6 + Number of frozen occupied orbitals: 0 + Number of frozen unoccupied orbitals: 0 + Two-electron integral type: Conventional + + + Computing Conventional Integrals Presorting SO-basis two-electron integrals. + Sorting File: SO Ints (nn|nn) nbuckets = 1 + Constructing frozen core operators + Starting first half-transformation. + Sorting half-transformed integrals. + First half integral transformation complete. + Starting second half-transformation. + Two-electron integral transformation complete. + + Integral transformation done. 0.00423067 s + Reading the two-electron integrals from disk + Size of two-electron integrals: 0.000029 GB + Timing for conventional integral transformation: 0.010 s. + Timing for freezing core and virtual orbitals: 0.000 s. + Timing for computing conventional integrals: 0.010 s. + + ----------------------------------------------------------- + Multi-Configurational Self Consistent Field + Two-Step Approximate Second-Order AO Algorithm + written by Chenyang Li, Kevin P. Hannon, and Shuhe Wang + ----------------------------------------------------------- + + + ==> MCSCF Calculation Information <== + + -------------------------------------------------------- + Print level Default + Integral type CONVENTIONAL + CI solver type FCI + Final orbital type CANONICAL + Derivative type NONE + Optimize orbitals TRUE + Include internal rotations FALSE + Debug printing FALSE + Energy convergence 1.000e-08 + Gradient convergence 1.000e-07 + Max value for rotation 2.000e-01 + Max number of macro iter. 100 + Max number of micro iter. for orbitals 6 + Max number of micro iter. for CI 12 + DIIS start 15 + Min DIIS vectors 3 + Max DIIS vectors 8 + Frequency of DIIS extrapolation 1 + -------------------------------------------------------- + + ==> Independent Orbital Rotations <== + + ORBITAL SPACES A1 A2 B1 B2 + ------------------------------------------------------------- + ACTIVE / RESTRICTED_DOCC 0 0 0 0 + RESTRICTED_UOCC / ACTIVE 0 0 0 0 + RESTRICTED_UOCC / RESTRICTED_DOCC 0 0 0 0 + ------------------------------------------------------------- + + ==> String Lists <== + + -------------------------------------------------------- + number of alpha electrons 2 + number of beta electrons 2 + number of alpha strings 15 + number of beta strings 15 + -------------------------------------------------------- + + ==> FCI Solver <== + + -------------------------------------------------------- + Spin adapt FALSE + Number of determinants 69 + Symmetry 0 + Multiplicity 1 + Number of roots 1 + Target root 0 + -------------------------------------------------------- + + ==> Initial Guess <== + + Initial guess determinants: 50 + + Classification of the initial guess solutions + + Number 2S+1 Selected + ------------------------ + 28 1 * + 21 3 + 1 5 + ------------------------ + + Spin Root Energy Status + ------------------------------------------------------- + singlet 0 -7.882494298809 -0.000000 added + ------------------------------------------------------- + + ==> Davidson-Liu Solver <== + + -------------------------------------------------------- + Print level Default + Energy convergence threshold 1.000e-08 + Residual convergence threshold 1.000e-06 + Schmidt orthogonality threshold 1.000e-12 + Schmidt discard threshold 1.000e-07 + Size of the space 69 + Number of roots 1 + Maximum number of iterations 100 + Collapse subspace size 2 + Maximum subspace size 10 + -------------------------------------------------------- + + Davidson-Liu solver: adding 1 guess vectors + Iteration Average Energy max(∆E) max(Residual) Vectors + --------------------------------------------------------------------------------- + 0 -7.882494298809 7.882494298809 0.008038987350 1 + 1 -7.882504307468 0.000010008659 0.000284247755 2 + 2 -7.882504329166 0.000000021698 0.000021027585 3 + 3 -7.882504329345 0.000000000179 0.000006599970 4 + 4 -7.882504329370 0.000000000025 0.000001340396 5 + 5 -7.882504329372 0.000000000002 0.000000531683 6 + --------------------------------------------------------------------------------- + + ==> Root No. 0 <== + + 2200 0 0 -0.98722460 + 2002 0 0 0.11303965 + 20ba 0 0 0.05889309 + 20ab 0 0 0.05889309 + 2ab0 0 0 -0.03825367 + 2ba0 0 0 -0.03825367 + 2020 0 0 0.03413984 + 2000 2 0 0.02725508 + 2000 0 2 0.02725508 + + Total Energy: -7.882504329372, : 0.000000 + +==> RDMs Test <== + + AA 1-RDM Error : +1.573946e-16 + BB 1-RDM Error : +1.080000e-17 + AAAA 2-RDM Error : +3.705896e-17 + BBBB 2-RDM Error : +6.919200e-19 + ABAB 2-RDM Error : +1.832531e-16 + AAAAAA 3-RDM Error : +0.000000e+00 + AABAAB 3-RDM Error : +0.000000e+00 + ABBABB 3-RDM Error : +0.000000e+00 + BBBBBB 3-RDM Error : +0.000000e+00 + Time for FCI: 0.027383458000 + + ==> Energy Summary <== + + Multi.(2ms) Irrep. No. Energy + -------------------------------------------------------- + 1 ( 0) A1 0 -7.882504329372 0.000000 + -------------------------------------------------------- + +==> RDMs Test <== + + SF 1-RDM Error : +5.214325e-18 + + ==> Natural Orbitals <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B1 0.001534 1B2 0.001534 4A1 0.000059 + + + ==> Dipole Moments [e a0] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State DM_X DM_Y DM_Z |DM| + -------------------------------------------------------------------- + 0A1 0.00000000 0.00000000 -1.81815571 1.81815571 + -------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 1.49275188 1.49275188 + -------------------------------------------------------------------- + +==> RDMs Test <== + + SF 1-RDM Error : +5.214325e-18 + + ==> Natural Orbitals <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B1 0.001534 1B2 0.001534 4A1 0.000059 + + + ==> Quadrupole Moments [e a0^2] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State QM_XX QM_XY QM_XZ QM_YY QM_YZ QM_ZZ + -------------------------------------------------------------------------------------------------- + 0A1 -4.08221296 0.00000000 0.00000000 -4.08221296 0.00000000 -4.94856087 + -------------------------------------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 7.30707704 + -------------------------------------------------------------------------------------------------- + +==> RDMs Test <== + + SF 1-RDM Error : +5.214325e-18 + AAAA 2-RDM Error : +5.248916e-08 + BBBB 2-RDM Error : +5.248916e-08 + ABAB 2-RDM Error : +5.248916e-08 + + ==> Natural Orbitals <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B1 0.001534 1B2 0.001534 4A1 0.000059 + + + ==> String Lists <== + + -------------------------------------------------------- + number of alpha electrons 2 + number of beta electrons 2 + number of alpha strings 15 + number of beta strings 15 + -------------------------------------------------------- + + ==> FCI Solver <== + + -------------------------------------------------------- + Spin adapt FALSE + Number of determinants 69 + Symmetry 0 + Multiplicity 1 + Number of roots 1 + Target root 0 + -------------------------------------------------------- + + ==> Initial Guess <== + + Initial guess determinants: 50 + + Classification of the initial guess solutions + + Number 2S+1 Selected + ------------------------ + 28 1 * + 21 3 + 1 5 + ------------------------ + + Spin Root Energy Status + ------------------------------------------------------- + singlet 0 -7.882494298809 -0.000000 added + ------------------------------------------------------- + + ==> Davidson-Liu Solver <== + + -------------------------------------------------------- + Print level Default + Energy convergence threshold 1.000e-12 + Residual convergence threshold 1.000e-06 + Schmidt orthogonality threshold 1.000e-12 + Schmidt discard threshold 1.000e-07 + Size of the space 69 + Number of roots 1 + Maximum number of iterations 100 + Collapse subspace size 2 + Maximum subspace size 10 + -------------------------------------------------------- + + Davidson-Liu solver: adding 1 guess vectors + Iteration Average Energy max(∆E) max(Residual) Vectors + --------------------------------------------------------------------------------- + 0 -7.882494298809 7.882494298809 0.008038987350 1 + 1 -7.882504307468 0.000010008659 0.000284247755 2 + 2 -7.882504329166 0.000000021698 0.000021027585 3 + 3 -7.882504329345 0.000000000179 0.000006599970 4 + 4 -7.882504329370 0.000000000025 0.000001340396 5 + 5 -7.882504329372 0.000000000002 0.000000531683 6 + 6 -7.882504329372 0.000000000000 0.000000120716 7 + --------------------------------------------------------------------------------- + + ==> Root No. 0 <== + + 2200 0 0 0.98722452 + 2002 0 0 -0.11303969 + 20ba 0 0 -0.05889325 + 20ab 0 0 -0.05889325 + 2ab0 0 0 0.03825444 + 2ba0 0 0 0.03825444 + 2020 0 0 -0.03414006 + 2000 0 2 -0.02725499 + 2000 2 0 -0.02725499 + + Total Energy: -7.882504329372, : 0.000000 + +==> RDMs Test <== + + AA 1-RDM Error : +2.221717e-16 + BB 1-RDM Error : +1.115228e-16 + AAAA 2-RDM Error : +1.406016e-17 + BBBB 2-RDM Error : +5.650255e-17 + ABAB 2-RDM Error : +1.834372e-16 + AAAAAA 3-RDM Error : +0.000000e+00 + AABAAB 3-RDM Error : +0.000000e+00 + ABBABB 3-RDM Error : +0.000000e+00 + BBBBBB 3-RDM Error : +0.000000e+00 + Time for FCI: 0.023231209000 + + ==> Energy Summary <== + + Multi.(2ms) Irrep. No. Energy + -------------------------------------------------------- + 1 ( 0) A1 0 -7.882504329372 0.000000 + -------------------------------------------------------- + +==> RDMs Test <== + + SF 1-RDM Error : +2.221734e-16 + + ==> Natural Orbitals <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B2 0.001534 1B1 0.001534 4A1 0.000059 + + + ==> Dipole Moments [e a0] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State DM_X DM_Y DM_Z |DM| + -------------------------------------------------------------------- + 0A1 0.00000000 0.00000000 -1.81815375 1.81815375 + -------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 1.49275188 1.49275188 + -------------------------------------------------------------------- + +==> RDMs Test <== + + SF 1-RDM Error : +2.221734e-16 + + ==> Natural Orbitals <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B2 0.001534 1B1 0.001534 4A1 0.000059 + + + ==> Quadrupole Moments [e a0^2] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State QM_XX QM_XY QM_XZ QM_YY QM_YZ QM_ZZ + -------------------------------------------------------------------------------------------------- + 0A1 -4.08221616 0.00000000 0.00000000 -4.08221616 0.00000000 -4.94855691 + -------------------------------------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 7.30707704 + -------------------------------------------------------------------------------------------------- + + Time to prepare integrals: 0.122 seconds + Time to run job : 0.061 seconds + Total : 0.184 seconds + AA 1-RDM..............................................................................PASSED + BB 1-RDM..............................................................................PASSED + AAAA 2-RDM............................................................................PASSED + BBBB 2-RDM............................................................................PASSED + ABAB 2-RDM............................................................................PASSED + AABAAB 3-RDM..........................................................................PASSED + ABBABB 3-RDM..........................................................................PASSED + AAAAAA 3-RDM..........................................................................PASSED + BBBBBB 3-RDM..........................................................................PASSED + + Psi4 stopped on: Tuesday, 17 September 2024 03:25PM + Psi4 wall time for execution: 0:00:01.10 + +*** Psi4 exiting successfully. Buy a developer a beer! diff --git a/tests/methods/fci-rdms-3/input.dat b/tests/methods/fci-rdms-3/input.dat new file mode 100644 index 000000000..1da8aa246 --- /dev/null +++ b/tests/methods/fci-rdms-3/input.dat @@ -0,0 +1,42 @@ +#! Generated using commit GITCOMMIT + +import forte + +molecule { +0 1 +Li +H 1 R + +R = 3.0 +units bohr +} + +set { + basis sto-3g + reference rhf + scf_type pk + e_convergence 12 +} + +set forte { + active_space_solver genci + fci_test_rdms true +} + +energy('scf') +energy('forte') + +compare_values(0.0, variable("AA 1-RDM ERROR"),12, "AA 1-RDM") #TEST +compare_values(0.0, variable("BB 1-RDM ERROR"),12, "BB 1-RDM") #TEST +compare_values(0.0, variable("AAAA 2-RDM ERROR"),12, "AAAA 2-RDM") #TEST +compare_values(0.0, variable("BBBB 2-RDM ERROR"),12, "BBBB 2-RDM") #TEST +compare_values(0.0, variable("ABAB 2-RDM ERROR"),12, "ABAB 2-RDM") #TEST +compare_values(0.0, variable("AABAAB 3-RDM ERROR"),12, "AABAAB 3-RDM") #TEST +compare_values(0.0, variable("ABBABB 3-RDM ERROR"),12, "ABBABB 3-RDM") #TEST +compare_values(0.0, variable("AAAAAA 3-RDM ERROR"),12, "AAAAAA 3-RDM") #TEST +compare_values(0.0, variable("BBBBBB 3-RDM ERROR"),12, "BBBBBB 3-RDM") #TEST +compare_values(0.0, variable("AAAAAAAA 4-RDM ERROR"),12, "AAAAAAAA 4-RDM") #TEST +compare_values(0.0, variable("AAABAAAB 4-RDM ERROR"),12, "AAABAAAB 4-RDM") #TEST +compare_values(0.0, variable("AABBAABB 4-RDM ERROR"),12, "AABBAABB 4-RDM") #TEST +compare_values(0.0, variable("ABBBABBB 4-RDM ERROR"),12, "ABBBABBB 4-RDM") #TEST +compare_values(0.0, variable("BBBBBBBB 4-RDM ERROR"),12, "BBBBBBBB 4-RDM") #TEST diff --git a/tests/methods/fci-rdms-3/output.ref b/tests/methods/fci-rdms-3/output.ref new file mode 100644 index 000000000..4a27def6c --- /dev/null +++ b/tests/methods/fci-rdms-3/output.ref @@ -0,0 +1,956 @@ + + ----------------------------------------------------------------------- + Psi4: An Open-Source Ab Initio Electronic Structure Package + Psi4 1.10a1.dev53 + + Git: Rev {master} 45660d8 dirty + + + D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, + M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, + A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, + R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, + J. S. O'Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, + B. P. Pritchard, B. R. Brooks, H. F. Schaefer III, A. Yu. Sokolov, + K. Patkowski, A. E. DePrince III, U. Bozkaya, R. A. King, + F. A. Evangelista, J. M. Turney, T. D. Crawford, C. D. Sherrill, + J. Chem. Phys. 152(18) 184108 (2020). https://doi.org/10.1063/5.0006002 + + Additional Code Authors + E. T. Seidl, C. L. Janssen, E. F. Valeev, M. L. Leininger, + J. F. Gonthier, R. M. Richard, H. R. McAlexander, M. Saitow, X. Wang, + P. Verma, M. H. Lechner, A. Jiang, S. Behnle, A. G. Heide, + M. F. Herbst, and D. L. Poole + + Previous Authors, Complete List of Code Contributors, + and Citations for Specific Modules + https://github.com/psi4/psi4/blob/master/codemeta.json + https://github.com/psi4/psi4/graphs/contributors + http://psicode.org/psi4manual/master/introduction.html#citing-psifour + + ----------------------------------------------------------------------- + + + Psi4 started on: Tuesday, 17 September 2024 03:25PM + + Process ID: 19290 + Host: Brian-Zs-MBA.local + PSIDATADIR: /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4 + Memory: 500.0 MiB + Threads: 1 + + ==> Input File <== + +-------------------------------------------------------------------------- +#! Generated using commit GITCOMMIT + +import forte + +molecule { +0 1 +Li +H 1 R + +R = 3.0 +units bohr +} + +set { + basis sto-3g + reference rhf + scf_type pk + e_convergence 12 +} + +set forte { + active_space_solver genci + fci_test_rdms true +} + +energy('scf') +energy('forte') + +compare_values(0.0, variable("AA 1-RDM ERROR"),12, "AA 1-RDM") #TEST +compare_values(0.0, variable("BB 1-RDM ERROR"),12, "BB 1-RDM") #TEST +compare_values(0.0, variable("AAAA 2-RDM ERROR"),12, "AAAA 2-RDM") #TEST +compare_values(0.0, variable("BBBB 2-RDM ERROR"),12, "BBBB 2-RDM") #TEST +compare_values(0.0, variable("ABAB 2-RDM ERROR"),12, "ABAB 2-RDM") #TEST +compare_values(0.0, variable("AABAAB 3-RDM ERROR"),12, "AABAAB 3-RDM") #TEST +compare_values(0.0, variable("ABBABB 3-RDM ERROR"),12, "ABBABB 3-RDM") #TEST +compare_values(0.0, variable("AAAAAA 3-RDM ERROR"),12, "AAAAAA 3-RDM") #TEST +compare_values(0.0, variable("BBBBBB 3-RDM ERROR"),12, "BBBBBB 3-RDM") #TEST +compare_values(0.0, variable("AAAAAAAA 4-RDM ERROR"),12, "AAAAAAAA 4-RDM") #TEST +compare_values(0.0, variable("AAABAAAB 4-RDM ERROR"),12, "AAABAAAB 4-RDM") #TEST +compare_values(0.0, variable("AABBAABB 4-RDM ERROR"),12, "AABBAABB 4-RDM") #TEST +compare_values(0.0, variable("ABBBABBB 4-RDM ERROR"),12, "ABBBABBB 4-RDM") #TEST +compare_values(0.0, variable("BBBBBBBB 4-RDM ERROR"),12, "BBBBBBBB 4-RDM") #TEST +-------------------------------------------------------------------------- + +Scratch directory: /tmp/ + => Libint2 <= + + Primary basis highest AM E, G, H: 6, 6, 3 + Auxiliary basis highest AM E, G, H: 7, 7, 4 + Onebody basis highest AM E, G, H: -, -, - + Solid Harmonics ordering: Gaussian + +*** tstart() called on Brian-Zs-MBA.local +*** at Tue Sep 17 15:25:33 2024 + + => Loading Basis Set <= + + Name: STO-3G + Role: ORBITAL + Keyword: BASIS + atoms 1 entry LI line 31 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + atoms 2 entry H line 19 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + + + --------------------------------------------------------- + SCF + by Justin Turney, Rob Parrish, Andy Simmonett + and Daniel G. A. Smith + RHF Reference + 1 Threads, 500 MiB Core + --------------------------------------------------------- + + ==> Geometry <== + + Molecular point group: c2v + Full point group: C_inf_v + + Geometry (in Bohr), charge = 0, multiplicity = 1: + + Center X Y Z Mass + ------------ ----------------- ----------------- ----------------- ----------------- + LI 0.000000000000 0.000000000000 -0.376812030371 7.016003436600 + H 0.000000000000 0.000000000000 2.623187969629 1.007825032230 + + Running in c2v symmetry. + + Rotational constants: A = ************ B = 7.59029 C = 7.59029 [cm^-1] + Rotational constants: A = ************ B = 227551.19787 C = 227551.19787 [MHz] + Nuclear repulsion = 1.000000000000000 + + Charge = 0 + Multiplicity = 1 + Electrons = 4 + Nalpha = 2 + Nbeta = 2 + + ==> Algorithm <== + + SCF Algorithm Type is PK. + DIIS enabled. + MOM disabled. + Fractional occupation disabled. + Guess Type is SAD. + Energy threshold = 1.00e-12 + Density threshold = 1.00e-06 + Integral threshold = 1.00e-12 + + ==> Primary Basis <== + + Basis Set: STO-3G + Blend: STO-3G + Number of shells: 4 + Number of basis functions: 6 + Number of Cartesian functions: 6 + Spherical Harmonics?: true + Max angular momentum: 1 + + ==> Integral Setup <== + + Using in-core PK algorithm. + Calculation information: + Number of atoms: 2 + Number of AO shells: 4 + Number of primitives: 12 + Number of atomic orbitals: 6 + Number of basis functions: 6 + + Integral cutoff 1.00e-12 + Number of threads: 1 + + Performing in-core PK + Using 462 doubles for integral storage. + We computed 55 shell quartets total. + Whereas there are 55 unique shell quartets. + + ==> DiskJK: Disk-Based J/K Matrices <== + + J tasked: Yes + K tasked: Yes + wK tasked: No + Memory [MiB]: 375 + Schwarz Cutoff: 1E-12 + + OpenMP threads: 1 + + Minimum eigenvalue in the overlap matrix is 3.4333995519E-01. + Reciprocal condition number of the overlap matrix is 2.0339047710E-01. + Using symmetric orthogonalization. + + ==> Pre-Iterations <== + + SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF (no occupation information). + + ------------------------- + Irrep Nso Nmo + ------------------------- + A1 4 4 + A2 0 0 + B1 1 1 + B2 1 1 + ------------------------- + Total 6 6 + ------------------------- + + ==> Iterations <== + + Total Energy Delta E RMS |[F,P]| + + @RHF iter SAD: -6.75622365412217 -6.75622e+00 0.00000e+00 + @RHF iter 1: -7.86092861588784 -1.10470e+00 8.57305e-03 DIIS/ADIIS + @RHF iter 2: -7.86223475166447 -1.30614e-03 8.22758e-04 DIIS/ADIIS + @RHF iter 3: -7.86224589388798 -1.11422e-05 7.41119e-05 DIIS + @RHF iter 4: -7.86224624208850 -3.48201e-07 2.61264e-05 DIIS + @RHF iter 5: -7.86224631040868 -6.83202e-08 3.24230e-07 DIIS + @RHF iter 6: -7.86224631041032 -1.64135e-12 1.56779e-08 DIIS + @RHF iter 7: -7.86224631041034 -1.77636e-14 2.79325e-12 DIIS + Energy and wave function converged. + + + ==> Post-Iterations <== + + Orbital Energies [Eh] + --------------------- + + Doubly Occupied: + + 1A1 -2.348477 2A1 -0.286330 + + Virtual: + + 3A1 0.078326 1B1 0.163933 1B2 0.163933 + 4A1 0.551172 + + Final Occupation by Irrep: + A1 A2 B1 B2 + DOCC [ 2, 0, 0, 0 ] + NA [ 2, 0, 0, 0 ] + NB [ 2, 0, 0, 0 ] + + @RHF Final Energy: -7.86224631041034 + + => Energetics <= + + Nuclear Repulsion Energy = 1.0000000000000000 + One-Electron Energy = -12.4548780553254126 + Two-Electron Energy = 3.5926317449150740 + Total Energy = -7.8622463104103382 + +Computation Completed + + +Properties will be evaluated at 0.000000, 0.000000, 0.000000 [a0] + +Properties computed using the SCF density matrix + + + Multipole Moments: + + ------------------------------------------------------------------------------------ + Multipole Electronic (a.u.) Nuclear (a.u.) Total (a.u.) + ------------------------------------------------------------------------------------ + + L = 1. Multiply by 2.5417464519 to convert [e a0] to [Debye] + Dipole X : 0.0000000 0.0000000 0.0000000 + Dipole Y : 0.0000000 0.0000000 0.0000000 + Dipole Z : -3.4031201 1.4927519 -1.9103682 + Magnitude : 1.9103682 + + ------------------------------------------------------------------------------------ + +*** tstop() called on Brian-Zs-MBA.local at Tue Sep 17 15:25:33 2024 +Module time: + user time = 0.09 seconds = 0.00 minutes + system time = 0.01 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes +Total time: + user time = 0.09 seconds = 0.00 minutes + system time = 0.01 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes + +Scratch directory: /tmp/ + + Forte + ---------------------------------------------------------------------------- + A suite of quantum chemistry methods for strongly correlated electrons + + git branch: rdm4 - git commit: 7780d650 + + Developed by: + Francesco A. Evangelista, Chenyang Li, Kevin P. Hannon, + Jeffrey B. Schriber, Tianyuan Zhang, Chenxi Cai, + Nan He, Nicholas Stair, Shuhe Wang, Renke Huang + ---------------------------------------------------------------------------- + + + Preparing forte objects from a Psi4 Wavefunction object + No reference wave function provided for Forte. Computing SCF orbitals using Psi4 ... + => Libint2 <= + + Primary basis highest AM E, G, H: 6, 6, 3 + Auxiliary basis highest AM E, G, H: 7, 7, 4 + Onebody basis highest AM E, G, H: -, -, - + Solid Harmonics ordering: Gaussian + +*** tstart() called on Brian-Zs-MBA.local +*** at Tue Sep 17 15:25:33 2024 + + => Loading Basis Set <= + + Name: STO-3G + Role: ORBITAL + Keyword: BASIS + atoms 1 entry LI line 31 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + atoms 2 entry H line 19 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + + + --------------------------------------------------------- + SCF + by Justin Turney, Rob Parrish, Andy Simmonett + and Daniel G. A. Smith + RHF Reference + 1 Threads, 500 MiB Core + --------------------------------------------------------- + + ==> Geometry <== + + Molecular point group: c2v + Full point group: C_inf_v + + Geometry (in Bohr), charge = 0, multiplicity = 1: + + Center X Y Z Mass + ------------ ----------------- ----------------- ----------------- ----------------- + LI 0.000000000000 0.000000000000 -0.376812030371 7.016003436600 + H 0.000000000000 0.000000000000 2.623187969629 1.007825032230 + + Running in c2v symmetry. + + Rotational constants: A = ************ B = 7.59029 C = 7.59029 [cm^-1] + Rotational constants: A = ************ B = 227551.19787 C = 227551.19787 [MHz] + Nuclear repulsion = 1.000000000000000 + + Charge = 0 + Multiplicity = 1 + Electrons = 4 + Nalpha = 2 + Nbeta = 2 + + ==> Algorithm <== + + SCF Algorithm Type is PK. + DIIS enabled. + MOM disabled. + Fractional occupation disabled. + Guess Type is SAD. + Energy threshold = 1.00e-12 + Density threshold = 1.00e-08 + Integral threshold = 1.00e-12 + + ==> Primary Basis <== + + Basis Set: STO-3G + Blend: STO-3G + Number of shells: 4 + Number of basis functions: 6 + Number of Cartesian functions: 6 + Spherical Harmonics?: true + Max angular momentum: 1 + + ==> Integral Setup <== + + Using in-core PK algorithm. + Calculation information: + Number of atoms: 2 + Number of AO shells: 4 + Number of primitives: 12 + Number of atomic orbitals: 6 + Number of basis functions: 6 + + Integral cutoff 1.00e-12 + Number of threads: 1 + + Performing in-core PK + Using 462 doubles for integral storage. + We computed 55 shell quartets total. + Whereas there are 55 unique shell quartets. + + ==> DiskJK: Disk-Based J/K Matrices <== + + J tasked: Yes + K tasked: Yes + wK tasked: No + Memory [MiB]: 375 + Schwarz Cutoff: 1E-12 + + OpenMP threads: 1 + + Minimum eigenvalue in the overlap matrix is 3.4333995519E-01. + Reciprocal condition number of the overlap matrix is 2.0339047710E-01. + Using symmetric orthogonalization. + + ==> Pre-Iterations <== + + SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF (no occupation information). + + ------------------------- + Irrep Nso Nmo + ------------------------- + A1 4 4 + A2 0 0 + B1 1 1 + B2 1 1 + ------------------------- + Total 6 6 + ------------------------- + + ==> Iterations <== + + Total Energy Delta E RMS |[F,P]| + + @RHF iter SAD: -6.75622365412217 -6.75622e+00 0.00000e+00 + @RHF iter 1: -7.86092861588784 -1.10470e+00 8.57305e-03 DIIS/ADIIS + @RHF iter 2: -7.86223475166447 -1.30614e-03 8.22758e-04 DIIS/ADIIS + @RHF iter 3: -7.86224589388798 -1.11422e-05 7.41119e-05 DIIS + @RHF iter 4: -7.86224624208850 -3.48201e-07 2.61264e-05 DIIS + @RHF iter 5: -7.86224631040868 -6.83202e-08 3.24230e-07 DIIS + @RHF iter 6: -7.86224631041032 -1.64135e-12 1.56779e-08 DIIS + @RHF iter 7: -7.86224631041034 -1.77636e-14 2.79325e-12 DIIS + Energy and wave function converged. + + + ==> Post-Iterations <== + + Orbital Energies [Eh] + --------------------- + + Doubly Occupied: + + 1A1 -2.348477 2A1 -0.286330 + + Virtual: + + 3A1 0.078326 1B1 0.163933 1B2 0.163933 + 4A1 0.551172 + + Final Occupation by Irrep: + A1 A2 B1 B2 + DOCC [ 2, 0, 0, 0 ] + NA [ 2, 0, 0, 0 ] + NB [ 2, 0, 0, 0 ] + + @RHF Final Energy: -7.86224631041034 + + => Energetics <= + + Nuclear Repulsion Energy = 1.0000000000000000 + One-Electron Energy = -12.4548780553254126 + Two-Electron Energy = 3.5926317449150740 + Total Energy = -7.8622463104103382 + +Computation Completed + + +Properties will be evaluated at 0.000000, 0.000000, 0.000000 [a0] + +Properties computed using the SCF density matrix + + + Multipole Moments: + + ------------------------------------------------------------------------------------ + Multipole Electronic (a.u.) Nuclear (a.u.) Total (a.u.) + ------------------------------------------------------------------------------------ + + L = 1. Multiply by 2.5417464519 to convert [e a0] to [Debye] + Dipole X : 0.0000000 0.0000000 0.0000000 + Dipole Y : 0.0000000 0.0000000 0.0000000 + Dipole Z : -3.4031201 1.4927519 -1.9103682 + Magnitude : 1.9103682 + + ------------------------------------------------------------------------------------ + +*** tstop() called on Brian-Zs-MBA.local at Tue Sep 17 15:25:33 2024 +Module time: + user time = 0.09 seconds = 0.00 minutes + system time = 0.00 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes +Total time: + user time = 0.31 seconds = 0.01 minutes + system time = 0.01 seconds = 0.00 minutes + total time = 0 seconds = 0.00 minutes + + + ==> MO Space Information <== + + ------------------------------------------------- + A1 A2 B1 B2 Sum + ------------------------------------------------- + FROZEN_DOCC 0 0 0 0 0 + RESTRICTED_DOCC 0 0 0 0 0 + GAS1 4 0 1 1 6 + GAS2 0 0 0 0 0 + GAS3 0 0 0 0 0 + GAS4 0 0 0 0 0 + GAS5 0 0 0 0 0 + GAS6 0 0 0 0 0 + RESTRICTED_UOCC 0 0 0 0 0 + FROZEN_UOCC 0 0 0 0 0 + Total 4 0 1 1 6 + ------------------------------------------------- => Loading Basis Set <= + + Name: STO-3G + Role: ORBITAL + Keyword: MINAO_BASIS + atoms 1 entry LI line 31 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + atoms 2 entry H line 19 file /Users/brianz98/local/psi4/objdir-Release/stage/share/psi4/basis/sto-3g.gbs + + + State Singlet (Ms = 0) A1 GAS min: 0 0 0 0 0 0 ; GAS max: 12 0 0 0 0 0 ; weights: + 1.000000000000 + Forte will use psi4 integrals + + ==> Primary Basis Set Summary <== + + Basis Set: STO-3G + Blend: STO-3G + Number of shells: 4 + Number of basis functions: 6 + Number of Cartesian functions: 6 + Spherical Harmonics?: true + Max angular momentum: 1 + + + JK created using conventional PK integrals + Using in-core PK algorithm. + Calculation information: + Number of atoms: 2 + Number of AO shells: 4 + Number of primitives: 12 + Number of atomic orbitals: 6 + Number of basis functions: 6 + + Integral cutoff 1.00e-12 + Number of threads: 1 + + Performing in-core PK + Using 462 doubles for integral storage. + We computed 55 shell quartets total. + Whereas there are 55 unique shell quartets. + + ==> DiskJK: Disk-Based J/K Matrices <== + + J tasked: Yes + K tasked: Yes + wK tasked: No + Memory [MiB]: 400 + Schwarz Cutoff: 1E-12 + + OpenMP threads: 1 + + + + ==> Integral Transformation <== + + Number of molecular orbitals: 6 + Number of correlated molecular orbitals: 6 + Number of frozen occupied orbitals: 0 + Number of frozen unoccupied orbitals: 0 + Two-electron integral type: Conventional + + + Computing Conventional Integrals Presorting SO-basis two-electron integrals. + Sorting File: SO Ints (nn|nn) nbuckets = 1 + Constructing frozen core operators + Starting first half-transformation. + Sorting half-transformed integrals. + First half integral transformation complete. + Starting second half-transformation. + Two-electron integral transformation complete. + + Integral transformation done. 0.00360425 s + Reading the two-electron integrals from disk + Size of two-electron integrals: 0.000029 GB + Timing for conventional integral transformation: 0.009 s. + Timing for freezing core and virtual orbitals: 0.000 s. + Timing for computing conventional integrals: 0.009 s. + + ----------------------------------------------------------- + Multi-Configurational Self Consistent Field + Two-Step Approximate Second-Order AO Algorithm + written by Chenyang Li, Kevin P. Hannon, and Shuhe Wang + ----------------------------------------------------------- + + + ==> MCSCF Calculation Information <== + + -------------------------------------------------------- + Print level Default + Integral type CONVENTIONAL + CI solver type GENCI + Final orbital type CANONICAL + Derivative type NONE + Optimize orbitals TRUE + Include internal rotations FALSE + Debug printing FALSE + Energy convergence 1.000e-08 + Gradient convergence 1.000e-07 + Max value for rotation 2.000e-01 + Max number of macro iter. 100 + Max number of micro iter. for orbitals 6 + Max number of micro iter. for CI 12 + DIIS start 15 + Min DIIS vectors 3 + Max DIIS vectors 8 + Frequency of DIIS extrapolation 1 + -------------------------------------------------------- + + ==> Independent Orbital Rotations <== + + ORBITAL SPACES A1 A2 B1 B2 + ------------------------------------------------------------- + ACTIVE / RESTRICTED_DOCC 0 0 0 0 + RESTRICTED_UOCC / ACTIVE 0 0 0 0 + RESTRICTED_UOCC / RESTRICTED_DOCC 0 0 0 0 + ------------------------------------------------------------- + + ==> Possible Electron Occupations <== + + Config. Space 1 + α β + ----------------- + 1 2 2 + + ==> String Lists <== + + -------------------------------------------------------- + number of alpha electrons 2 + number of beta electrons 2 + number of alpha strings 15 + number of beta strings 15 + -------------------------------------------------------- + + ==> String-based CI Solver <== + + -------------------------------------------------------- + Print level Default + Spin adapt FALSE + Number of determinants 69 + Symmetry 0 + Multiplicity 1 + Number of roots 1 + Target root 0 + -------------------------------------------------------- + + ==> Initial Guess <== + + Initial guess determinants: 50 + + Classification of the initial guess solutions + + Number 2S+1 Selected + ------------------------ + 28 1 * + 21 3 + 1 5 + ------------------------ + + Spin Root Energy Status + ------------------------------------------------------- + singlet 0 -7.882494298809 -0.000000 added + ------------------------------------------------------- + + ==> Davidson-Liu Solver <== + + -------------------------------------------------------- + Print level Default + Energy convergence threshold 1.000e-08 + Residual convergence threshold 1.000e-06 + Schmidt orthogonality threshold 1.000e-12 + Schmidt discard threshold 1.000e-07 + Size of the space 69 + Number of roots 1 + Maximum number of iterations 100 + Collapse subspace size 2 + Maximum subspace size 10 + -------------------------------------------------------- + + Davidson-Liu solver: adding 1 guess vectors + Iteration Average Energy max(∆E) max(Residual) Vectors + --------------------------------------------------------------------------------- + 0 -7.882494298809 7.882494298809 0.008038987350 1 + 1 -7.882504307468 0.000010008659 0.000284247755 2 + 2 -7.882504329166 0.000000021698 0.000021027585 3 + 3 -7.882504329345 0.000000000179 0.000006599970 4 + 4 -7.882504329370 0.000000000025 0.000001340396 5 + 5 -7.882504329372 0.000000000002 0.000000531683 6 + --------------------------------------------------------------------------------- + Timing for CI: 0.006 s. + + ==> Root No. 0 <== + + 2200 0 0 -0.98722460 + 2002 0 0 0.11303965 + + Total Energy: -7.882504329372, : 0.000000 + Timing for 3-RDM (aab): 0.000 s + Timing for 3-RDM (abb): 0.000 s + Timing for 4-RDM (aabb): 0.005 s + +==> RDMs Test (max level = 4)<== + + AA 1-RDM Error : +2.483883e-16 + BB 1-RDM Error : +1.573397e-16 + AAAA 2-RDM Error : +4.764858e-16 + BBBB 2-RDM Error : +4.767134e-16 + ABAB 2-RDM Error : +2.070125e-16 + AAAAAA 3-RDM Error : +0.000000e+00 + AABAAB 3-RDM Error : +3.585422e-17 + ABBABB 3-RDM Error : +1.834991e-18 + BBBBBB 3-RDM Error : +0.000000e+00 + AAAAAAAA 4-RDM Error : +0.000000e+00 + AAABAAAB 4-RDM Error : +0.000000e+00 + AABBAABB 4-RDM Error : +0.000000e+00 + ABBBABBB 4-RDM Error : +0.000000e+00 + BBBBBBBB 4-RDM Error : +0.000000e+00 + + ==> Energy Summary <== + + Multi.(2ms) Irrep. No. Energy + -------------------------------------------------------- + 1 ( 0) A1 0 -7.882504329372 0.000000 + -------------------------------------------------------- + +==> RDMs Test (max level = 1)<== + + SF 1-RDM Error : +4.975030e-16 + + ==> Natural Orbitals Occupation Numbers <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B1 0.001534 1B2 0.001534 4A1 0.000059 + + + ==> Dipole Moments [e a0] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State DM_X DM_Y DM_Z |DM| + -------------------------------------------------------------------- + 0A1 0.00000000 0.00000000 -1.81815571 1.81815571 + -------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 1.49275188 1.49275188 + -------------------------------------------------------------------- + +==> RDMs Test (max level = 1)<== + + SF 1-RDM Error : +4.975030e-16 + + ==> Natural Orbitals Occupation Numbers <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B1 0.001534 1B2 0.001534 4A1 0.000059 + + + ==> Quadrupole Moments [e a0^2] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State QM_XX QM_XY QM_XZ QM_YY QM_YZ QM_ZZ + -------------------------------------------------------------------------------------------------- + 0A1 -4.08221296 0.00000000 0.00000000 -4.08221296 0.00000000 -4.94856087 + -------------------------------------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 7.30707704 + -------------------------------------------------------------------------------------------------- + +==> RDMs Test (max level = 2)<== + + SF 1-RDM Error : +4.975030e-16 + AAAA 2-RDM Error : +5.248916e-08 + BBBB 2-RDM Error : +5.248916e-08 + ABAB 2-RDM Error : +5.248916e-08 + + ==> Natural Orbitals Occupation Numbers <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B1 0.001534 1B2 0.001534 4A1 0.000059 + + + ==> Possible Electron Occupations <== + + Config. Space 1 + α β + ----------------- + 1 2 2 + + ==> String Lists <== + + -------------------------------------------------------- + number of alpha electrons 2 + number of beta electrons 2 + number of alpha strings 15 + number of beta strings 15 + -------------------------------------------------------- + + ==> String-based CI Solver <== + + -------------------------------------------------------- + Print level Default + Spin adapt FALSE + Number of determinants 69 + Symmetry 0 + Multiplicity 1 + Number of roots 1 + Target root 0 + -------------------------------------------------------- + + ==> Initial Guess <== + + Initial guess determinants: 50 + + Classification of the initial guess solutions + + Number 2S+1 Selected + ------------------------ + 28 1 * + 21 3 + 1 5 + ------------------------ + + Spin Root Energy Status + ------------------------------------------------------- + singlet 0 -7.882494298809 -0.000000 added + ------------------------------------------------------- + + ==> Davidson-Liu Solver <== + + -------------------------------------------------------- + Print level Default + Energy convergence threshold 1.000e-12 + Residual convergence threshold 1.000e-06 + Schmidt orthogonality threshold 1.000e-12 + Schmidt discard threshold 1.000e-07 + Size of the space 69 + Number of roots 1 + Maximum number of iterations 100 + Collapse subspace size 2 + Maximum subspace size 10 + -------------------------------------------------------- + + Davidson-Liu solver: adding 1 guess vectors + Iteration Average Energy max(∆E) max(Residual) Vectors + --------------------------------------------------------------------------------- + 0 -7.882494298809 7.882494298809 0.008038987350 1 + 1 -7.882504307468 0.000010008659 0.000284247755 2 + 2 -7.882504329166 0.000000021698 0.000021027585 3 + 3 -7.882504329345 0.000000000179 0.000006599970 4 + 4 -7.882504329370 0.000000000025 0.000001340396 5 + 5 -7.882504329372 0.000000000002 0.000000531683 6 + 6 -7.882504329372 0.000000000000 0.000000120716 7 + --------------------------------------------------------------------------------- + Timing for CI: 0.001 s. + + ==> Root No. 0 <== + + 2200 0 0 0.98722452 + 2002 0 0 -0.11303969 + + Total Energy: -7.882504329372, : 0.000000 + Timing for 3-RDM (aab): 0.000 s + Timing for 3-RDM (abb): 0.000 s + Timing for 4-RDM (aabb): 0.001 s + +==> RDMs Test (max level = 4)<== + + AA 1-RDM Error : +2.223021e-16 + BB 1-RDM Error : +1.140409e-16 + AAAA 2-RDM Error : +1.112128e-17 + BBBB 2-RDM Error : +6.716114e-17 + ABAB 2-RDM Error : +4.478902e-16 + AAAAAA 3-RDM Error : +0.000000e+00 + AABAAB 3-RDM Error : +7.676617e-18 + ABBABB 3-RDM Error : +7.248220e-19 + BBBBBB 3-RDM Error : +0.000000e+00 + AAAAAAAA 4-RDM Error : +0.000000e+00 + AAABAAAB 4-RDM Error : +0.000000e+00 + AABBAABB 4-RDM Error : +0.000000e+00 + ABBBABBB 4-RDM Error : +0.000000e+00 + BBBBBBBB 4-RDM Error : +0.000000e+00 + + ==> Energy Summary <== + + Multi.(2ms) Irrep. No. Energy + -------------------------------------------------------- + 1 ( 0) A1 0 -7.882504329372 0.000000 + -------------------------------------------------------- + +==> RDMs Test (max level = 1)<== + + SF 1-RDM Error : +2.247255e-16 + + ==> Natural Orbitals Occupation Numbers <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B2 0.001534 1B1 0.001534 4A1 0.000059 + + + ==> Dipole Moments [e a0] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State DM_X DM_Y DM_Z |DM| + -------------------------------------------------------------------- + 0A1 0.00000000 0.00000000 -1.81815375 1.81815375 + -------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 1.49275188 1.49275188 + -------------------------------------------------------------------- + +==> RDMs Test (max level = 1)<== + + SF 1-RDM Error : +2.247255e-16 + + ==> Natural Orbitals Occupation Numbers <== + + 1A1 1.999911 2A1 1.955085 3A1 0.041876 + 1B2 0.001534 1B1 0.001534 4A1 0.000059 + + + ==> Quadrupole Moments [e a0^2] (Nuclear + Electronic) for Singlet (Ms = 0) A1 <== + + State QM_XX QM_XY QM_XZ QM_YY QM_YZ QM_ZZ + -------------------------------------------------------------------------------------------------- + 0A1 -4.08221616 0.00000000 0.00000000 -4.08221616 0.00000000 -4.94855691 + -------------------------------------------------------------------------------------------------- + Nuclear 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 7.30707704 + -------------------------------------------------------------------------------------------------- + + Time to prepare integrals: 0.128 seconds + Time to run job : 0.416 seconds + Total : 0.544 seconds + AA 1-RDM..............................................................................PASSED + BB 1-RDM..............................................................................PASSED + AAAA 2-RDM............................................................................PASSED + BBBB 2-RDM............................................................................PASSED + ABAB 2-RDM............................................................................PASSED + AABAAB 3-RDM..........................................................................PASSED + ABBABB 3-RDM..........................................................................PASSED + AAAAAA 3-RDM..........................................................................PASSED + BBBBBB 3-RDM..........................................................................PASSED + AAAAAAAA 4-RDM........................................................................PASSED + AAABAAAB 4-RDM........................................................................PASSED + AABBAABB 4-RDM........................................................................PASSED + ABBBABBB 4-RDM........................................................................PASSED + BBBBBBBB 4-RDM........................................................................PASSED + + Psi4 stopped on: Tuesday, 17 September 2024 03:25PM + Psi4 wall time for execution: 0:00:01.46 + +*** Psi4 exiting successfully. Buy a developer a beer! diff --git a/tests/methods/tests.yaml b/tests/methods/tests.yaml index 9402f194d..8ea1970d4 100644 --- a/tests/methods/tests.yaml +++ b/tests/methods/tests.yaml @@ -259,6 +259,7 @@ fci: - fci-ecp-1 - fci-ecp-2 - fci-rdms-2 + - fci-rdms-3 - fci-trdms-1 - fci-trdms-2 long: