B-MRPT2

forte/CMakeLists.txt

@@ -154,6 +154,7 @@ mrdsrg-spin-adapted/sadsrg.cc
 mrdsrg-spin-adapted/sadsrg_amps_analysis.cc
 mrdsrg-spin-adapted/sadsrg_block_labels.cc
 mrdsrg-spin-adapted/sadsrg_comm.cc
+mrdsrg-spin-adapted/sadsrg_orbital_rotation.cc
 mrdsrg-spin-adapted/sa_dsrgpt.cc
 mrdsrg-spin-adapted/sa_ldsrg2.cc
 mrdsrg-spin-adapted/sa_mrdsrg.cc

forte/api/forte_python_module.cc

@@ -366,7 +366,9 @@ PYBIND11_MODULE(_forte, m) {
              "Set the map from state to the weights of all computed roots")
         .def("set_read_cwd_amps", &SADSRG::set_read_amps_cwd,
              "Set if reading amplitudes in the current directory or not")
-        .def("clean_checkpoints", &SADSRG::clean_checkpoints, "Delete amplitudes checkpoint files");
+        .def("clean_checkpoints", &SADSRG::clean_checkpoints, "Delete amplitudes checkpoint files")
+        .def("is_brueckner_converged", &SADSRG::is_brueckner_converged,
+             "Return true if DSRG Brueckner orbital rotation is converged");
 
     // export MRDSRG_SO
     py::class_<MRDSRG_SO>(m, "MRDSRG_SO")

forte/api/integrals_api.cc

@@ -39,8 +39,8 @@ namespace forte {
 /// export ForteIntegrals
 void export_ForteIntegrals(py::module& m) {
     py::class_<ForteIntegrals, std::shared_ptr<ForteIntegrals>>(m, "ForteIntegrals")
-        .def("rotate_orbitals", &ForteIntegrals::rotate_orbitals, "Rotate MOs during contructor")
-        .def("nmo", &ForteIntegrals::nmo, "Return the total number of moleuclar orbitals")
+        .def("rotate_orbitals", &ForteIntegrals::rotate_orbitals, "Rotate MOs during constructor")
+        .def("nmo", &ForteIntegrals::nmo, "Return the total number of molecular orbitals")
         .def("ncmo", &ForteIntegrals::ncmo, "Return the number of correlated orbitals")
         .def("frozen_core_energy", &ForteIntegrals::frozen_core_energy,
              "Return the frozen core energy")

forte/mrdsrg-spin-adapted/sa_dsrgpt.cc

@@ -27,8 +27,6 @@
  * @END LICENSE
  */
 
-#include "psi4/libpsi4util/PsiOutStream.h"
-
 #include "helpers/disk_io.h"
 #include "helpers/printing.h"
 #include "helpers/timer.h"
@@ -59,7 +57,7 @@ void SA_DSRGPT::print_options() {
 
     if (ints_->integral_type() == Cholesky) {
         auto cholesky_threshold = foptions_->get_double("CHOLESKY_TOLERANCE");
-        calculation_info_double.push_back({"Cholesky tolerance", cholesky_threshold});
+        calculation_info_double.emplace_back("Cholesky tolerance", cholesky_threshold);
     }
 
     std::vector<std::pair<std::string, std::string>> calculation_info_string{
@@ -71,15 +69,15 @@ void SA_DSRGPT::print_options() {
         {"Internal amplitudes", internal_amp_}};
 
     if (multi_state_) {
-        calculation_info_string.push_back({"State type", "MULTIPLE STATES"});
-        calculation_info_string.push_back({"Multi-state type", multi_state_algorithm_});
+        calculation_info_string.emplace_back("State type", "MULTIPLE STATES");
+        calculation_info_string.emplace_back("Multi-state type", multi_state_algorithm_);
     } else {
-        calculation_info_string.push_back({"State type", "SINGLE STATE"});
+        calculation_info_string.emplace_back("State type", "SINGLE STATE");
     }
 
     if (internal_amp_ != "NONE") {
-        calculation_info_string.push_back({"Internal amplitudes levels", internal_amp_});
-        calculation_info_string.push_back({"Internal amplitudes selection", internal_amp_select_});
+        calculation_info_string.emplace_back("Internal amplitudes levels", internal_amp_);
+        calculation_info_string.emplace_back("Internal amplitudes selection", internal_amp_select_);
     }
 
     // Print some information

forte/mrdsrg-spin-adapted/sa_dsrgpt.h

@@ -66,10 +66,10 @@ class SA_DSRGPT : public SADSRG {
     ambit::BlockedTensor F0th_;
     /// Generalized Fock matrix (bare off-diagonal blocks)
     ambit::BlockedTensor F1st_;
-    /// Single excitation amplitude
-    ambit::BlockedTensor T1_;
-    /// Double excitation amplitude
-    ambit::BlockedTensor T2_;
+//    /// Single excitation amplitude
+//    ambit::BlockedTensor T1_;
+//    /// Double excitation amplitude
+//    ambit::BlockedTensor T2_;
     /// Double excitation amplitude (2 * J - K)
     ambit::BlockedTensor S2_;
 

forte/mrdsrg-spin-adapted/sa_ldsrg2.cc

@@ -93,13 +93,15 @@ double SA_MRDSRG::compute_energy_ldsrg2() {
         local_timer t_hbar;
         timer hbar("Compute Hbar");
         if (corrlv_string_ == "LDSRG2_QC") {
-            compute_hbar_qc();
+            if (sequential_Hbar_)
+                compute_hbar_qc_sequential();
+            else
+                compute_hbar_qc();
         } else {
-            if (sequential_Hbar_) {
+            if (sequential_Hbar_)
                 compute_hbar_sequential();
-            } else {
+            else
                 compute_hbar();
-            }
         }
         hbar.stop();
         double Edelta = Hbar0_ - Ecorr;
@@ -162,9 +164,9 @@ double SA_MRDSRG::compute_energy_ldsrg2() {
     outfile->Printf("\n    %s", dash.c_str());
     outfile->Printf("\n\n  ==> MR-LDSRG(2)%s Energy Summary <==\n", level.c_str());
     std::vector<std::pair<std::string, double>> energy;
-    energy.push_back({"E0 (reference)", Eref_});
-    energy.push_back({"MR-LDSRG(2) correlation energy", Ecorr});
-    energy.push_back({"MR-LDSRG(2) total energy", Eref_ + Ecorr});
+    energy.emplace_back("E0 (reference)", Eref_);
+    energy.emplace_back("MR-LDSRG(2) correlation energy", Ecorr);
+    energy.emplace_back("MR-LDSRG(2) total energy", Eref_ + Ecorr);
     for (auto& str_dim : energy) {
         outfile->Printf("\n    %-30s = %23.15f", str_dim.first.c_str(), str_dim.second);
     }
@@ -304,20 +306,7 @@ void SA_MRDSRG::compute_hbar_sequential() {
 
     timer rotation("Hbar T1 rotation");
 
-    ambit::BlockedTensor A1;
-    A1 = BTF_->build(tensor_type_, "A1 Amplitudes", {"gg"}, true);
-    A1["ia"] = T1_["ia"];
-    A1["ai"] -= T1_["ia"];
-
-    size_t ncmo = core_mos_.size() + actv_mos_.size() + virt_mos_.size();
-
-    psi::SharedMatrix A1_m(new psi::Matrix("A1 alpha", ncmo, ncmo));
-    A1.iterate([&](const std::vector<size_t>& i, const std::vector<SpinType>&, double& value) {
-        A1_m->set(i[0], i[1], value);
-    });
-
-    // >=3 is required for high energy convergence
-    A1_m->expm(3);
+    auto A1_m = expA1(T1_, false);
 
     ambit::BlockedTensor U1;
     U1 = BTF_->build(tensor_type_, "Transformer", {"gg"}, true);
@@ -570,6 +559,129 @@ void SA_MRDSRG::compute_hbar_qc() {
     Hbar1_["ia"] += temp["ai"];
 }
 
+void SA_MRDSRG::compute_hbar_qc_sequential() {
+    if (!eri_df_)
+        throw std::runtime_error("Not available for conventional integrals!");
+
+    timer rotation("Hbar T1 rotation");
+
+    auto A1_m = expA1(T1_, false);
+
+    auto U1 = BTF_->build(tensor_type_, "Transformer", {"gg"}, true);
+    U1.iterate([&](const std::vector<size_t>& i, const std::vector<SpinType>&, double& value) {
+        value = A1_m->get(i[0], i[1]);
+    });
+
+    /// Recompute Hbar0 (ref. energy + T1 correlation) and F (Fock)
+    /// E = 0.5 * ( H["ji"] + F["ji] ) * D1["ij"] + 0.25 * V["xyuv"] * L2["uvxy"]
+
+    // F is now "bare H1"
+    auto F = BTF_->build(tensor_type_, "Fock (A1 dressed)", {"gg"});
+    F["rs"] = U1["rp"] * H_["pq"] * U1["sq"];
+
+    Hbar0_ = 0.0;
+    F.block("cc").iterate([&](const std::vector<size_t>& i, double& value) {
+        if (i[0] == i[1])
+            Hbar0_ += value;
+    });
+    Hbar0_ += 0.5 * F["uv"] * L1_["vu"];
+
+    // for simplicity, create a core-core density matrix
+    auto D1c = BTF_->build(tensor_type_, "L1 core", {"cc"});
+    for (size_t m = 0, nc = core_mos_.size(); m < nc; ++m) {
+        D1c.block("cc").data()[m * nc + m] = 2.0;
+    }
+
+    // F becomes "Fock"
+    auto B = BTF_->build(tensor_type_, "B 3-idx", {"Lgg"}, true);
+    B["grs"] = U1["rp"] * B_["gpq"] * U1["sq"];
+
+    auto temp = BTF_->build(tensor_type_, "B temp", {"L"}, true);
+    temp["g"] = B["gmn"] * D1c["mn"];
+    temp["g"] += B["guv"] * L1_["uv"];
+    F["pq"] += temp["g"] * B["gpq"];
+
+    F["pq"] -= 0.5 * B["gpm"] * B["gnq"] * D1c["mn"];
+    F["pq"] -= 0.5 * B["gpu"] * B["gvq"] * L1_["uv"];
+
+    // compute fully contracted term from T1
+    F.block("cc").iterate([&](const std::vector<size_t>& i, double& value) {
+        if (i[0] == i[1])
+            Hbar0_ += value;
+    });
+    Hbar0_ += 0.5 * F["uv"] * L1_["vu"];
+
+    Hbar0_ += 0.5 * B["gux"] * B["gvy"] * L2_["xyuv"];
+
+    Hbar0_ += Efrzc_ + Enuc_ - Eref_;
+
+    rotation.stop();
+
+    ////////////////////////////////////////////////////////////////////////////////////
+
+    // iteration variables
+    timer comm("Hbar T2 commutator");
+
+    Hbar1_["ia"] = F["ia"];
+    Hbar2_["ijab"] = B["gia"] * B["gjb"];
+
+    // diagonal blocks of fock
+    auto fock = ambit::BlockedTensor::build(tensor_type_, "fock diag", {"cc", "aa", "vv"});
+    fock["pq"] = F_["pq"];
+
+    // final second-order Hbar in the active space
+    auto hbar2 = ambit::BlockedTensor::build(tensor_type_, "hbar2", {"aaaa"});
+    hbar2["pqrs"] = Hbar2_["pqrs"];
+
+    // S1 = [H0th, A2]_1
+    auto S1 = BTF_->build(tensor_type_, "S1", {"gg"}, true);
+    auto tmp1 = BTF_->build(tensor_type_, "tmp1", {"gg"}, true);
+    H1_T2_C1(fock, T2_, 1.0, tmp1);
+    S1["pq"] += tmp1["pq"];
+    S1["pq"] += tmp1["qp"];
+
+    S1["pq"] += 2.0 * F["pq"];
+
+    // S2 = [H0th, A2]_2
+    auto S2 = BTF_->build(tensor_type_, "S2", od_two_labels(), true);
+    auto tmp2 = BTF_->build(tensor_type_, "tmp2", {"hhpp"}, true);
+    H1_T2_C2(fock, T2_, 1.0, tmp2);
+    S2["pqrs"] += tmp2["pqrs"];
+    S2["pqrs"] += tmp2["rspq"];
+    Hbar2_["pqrs"] += S2["pqrs"];
+
+    S2["pqrs"] += 2.0 * B["gpr"] * B["gqs"];
+
+    tmp1.zero();
+    tmp2 = BTF_->build(tensor_type_, "tmp2", {"aaaa"}, true);
+
+    // 0.5 * [S1 + H1, A2]_0,1,2
+    H1_T2_C0(S1, T2_, 1.0, Hbar0_);
+    H1_T2_C1(S1, T2_, 0.5, tmp1);
+    H1_T2_C2(S1, T2_, 0.5, tmp2);
+
+    // 0.5 * [S2 + H2, A2]_0,1,2
+    H2_T2_C0(S2, T2_, DT2_, 1.0, Hbar0_);
+    H2_T2_C2(S2, T2_, DT2_, 0.5, tmp2);
+    H2_T2_C1(S2, T2_, DT2_, 0.5, tmp1);
+
+    S2 = BTF_->build(tensor_type_, "S2", diag_two_labels(), true);
+    S2["pqrs"] += B["gpr"] * B["gqs"];
+    H2_T2_C1(S2, T2_, DT2_, 1.0, tmp1);
+
+    // [H2, A2]_0,1,2
+//    V_T2_C0_DF(B, T2_, DT2_, 2.0, Hbar0_);
+//    V_T2_C1_DF(B, T2_, DT2_, 1.0, tmp1);
+//    V_T2_C2_DF(B, T2_, DT2_, 1.0, tmp2);
+
+    Hbar1_["pq"] += tmp1["pq"];
+    Hbar1_["pq"] += tmp1["qp"];
+    hbar2["uvxy"] += tmp2["uvxy"];
+    hbar2["uvxy"] += tmp2["xyuv"];
+
+    Hbar2_["uvxy"] = hbar2["uvxy"];
+}
+
 void SA_MRDSRG::setup_ldsrg2_tensors() {
     BlockedTensor::set_expert_mode(true);
 

forte/mrdsrg-spin-adapted/sa_mrdsrg.cc

@@ -34,6 +34,7 @@
 #include "psi4/libpsi4util/PsiOutStream.h"
 #include "psi4/libpsio/psio.hpp"
 
+#include "helpers/helpers.h"
 #include "helpers/printing.h"
 #include "helpers/timer.h"
 #include "sa_mrdsrg.h"
@@ -54,15 +55,15 @@ SA_MRDSRG::SA_MRDSRG(std::shared_ptr<RDMs> rdms, std::shared_ptr<SCFInfo> scf_in
 
 void SA_MRDSRG::read_options() {
     corrlv_string_ = foptions_->get_str("CORR_LEVEL");
-    std::vector<std::string> available{"LDSRG2", "LDSRG2_QC"};
+    std::vector<std::string> available{"LDSRG2", "LDSRG2_QC", "CC2"};
     if (std::find(available.begin(), available.end(), corrlv_string_) == available.end()) {
         outfile->Printf("\n  Warning: CORR_LEVEL option %s is not implemented.",
                         corrlv_string_.c_str());
         outfile->Printf("\n  Changed CORR_LEVEL option to LDSRG2_QC");
 
         corrlv_string_ = "LDSRG2_QC";
-        warnings_.push_back(std::make_tuple("Unsupported CORR_LEVEL", "Change to LDSRG2_QC",
-                                            "Change options in input.dat"));
+        warnings_.emplace_back("Unsupported CORR_LEVEL", "Change to LDSRG2_QC",
+                               "Change options in input.dat");
     }
 
     sequential_Hbar_ = foptions_->get_bool("DSRG_HBAR_SEQ");
@@ -71,7 +72,7 @@ void SA_MRDSRG::read_options() {
     rsc_ncomm_ = foptions_->get_int("DSRG_RSC_NCOMM");
     rsc_conv_ = foptions_->get_double("DSRG_RSC_THRESHOLD");
 
-    maxiter_ = foptions_->get_int("MAXITER");
+    maxiter_ = foptions_->get_int("DSRG_MAXITER");
     e_conv_ = foptions_->get_double("E_CONVERGENCE");
     r_conv_ = foptions_->get_double("R_CONVERGENCE");
 
@@ -106,6 +107,17 @@ void SA_MRDSRG::startup() {
 
     t1_file_cwd_ = "forte.mrdsrg.adapted.t1.bin";
     t2_file_cwd_ = "forte.mrdsrg.adapted.t2.bin";
+
+    // build transformation matrix to orthogonalize T1 excited basis
+    if (t1_type_ == "PROJECT") {
+        // allocate residuals in orthogonal basis
+        Oca_ = ambit::Tensor::build(tensor_type_, "Omega ca", {core_mos_.size(), Xca_.dim(1)});
+        Oav_ = ambit::Tensor::build(tensor_type_, "Omega av", {Xav_.dim(1), virt_mos_.size()});
+
+        // compute denominators
+        compute_proj_denom_ca();
+        compute_proj_denom_av();
+    }
 }
 
 void SA_MRDSRG::print_options() {
@@ -134,7 +146,8 @@ void SA_MRDSRG::print_options() {
         {"Reference relaxation", relax_ref_},
         {"3RDM algorithm", L3_algorithm_},
         {"Core-Virtual source type", ccvv_source_},
-        {"T1 amplitudes initial guess", t1_guess_}};
+        {"T1 amplitudes initial guess", t1_guess_},
+        {"T1 amplitudes type", t1_type_}};
 
     if (internal_amp_ != "NONE") {
         calculation_info_string.emplace_back("Internal amplitudes levels", internal_amp_);
@@ -148,9 +161,19 @@ void SA_MRDSRG::print_options() {
         {"Read amplitudes from current dir", read_amps_cwd_},
         {"Write amplitudes to current dir", dump_amps_cwd_}};
 
+    if (brueckner_) {
+        calculation_info_bool.emplace_back("DSRG Brueckner orbitals", brueckner_);
+        calculation_info_double.emplace_back("Brueckner convergence", brueckner_conv_);
+    }
+
     // print information
     print_selected_options("Computation Information", calculation_info_string,
                            calculation_info_bool, calculation_info_double, calculation_info_int);
+
+    // stop if there are some conflicts
+    if (corrlv_string_ == "LDSRG2_QC" and !eri_df_ and sequential_Hbar_) {
+        throw std::runtime_error("Sequential LDSRG2_QC only available with DF/CD integrals!");
+    }
 }
 
 void SA_MRDSRG::check_memory() {
@@ -221,7 +244,7 @@ void SA_MRDSRG::build_ints() {
 
 double SA_MRDSRG::compute_energy() {
     // build initial amplitudes
-    T1_ = BTF_->build(tensor_type_, "T1 Amplitudes", {"hp"});
+    //    T1_ = BTF_->build(tensor_type_, "T1 Amplitudes", {"hp"});
     T2_ = BTF_->build(tensor_type_, "T2 Amplitudes", {"hhpp"});
     guess_t(V_, T2_, F_, T1_, B_);
 
@@ -238,6 +261,10 @@ double SA_MRDSRG::compute_energy() {
     //    default: { Etotal += compute_energy_ldsrg2_qc(); }
     //    }
 
+    if (brueckner_) {
+        brueckner_orbital_rotation(T1_);
+    }
+
     return Etotal;
 }