diff --git a/forte/CMakeLists.txt b/forte/CMakeLists.txt
index 6568069a6..7a74cea42 100644
--- a/forte/CMakeLists.txt
+++ b/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
diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc
index d46333ad0..039e3635f 100644
--- a/forte/api/forte_python_module.cc
+++ b/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")
diff --git a/forte/api/integrals_api.cc b/forte/api/integrals_api.cc
index 561696390..a6888be78 100644
--- a/forte/api/integrals_api.cc
+++ b/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")
diff --git a/forte/mrdsrg-spin-adapted/sa_dsrgpt.cc b/forte/mrdsrg-spin-adapted/sa_dsrgpt.cc
index f42019656..8116bf465 100644
--- a/forte/mrdsrg-spin-adapted/sa_dsrgpt.cc
+++ b/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
diff --git a/forte/mrdsrg-spin-adapted/sa_dsrgpt.h b/forte/mrdsrg-spin-adapted/sa_dsrgpt.h
index fc757741b..f427d0fa0 100644
--- a/forte/mrdsrg-spin-adapted/sa_dsrgpt.h
+++ b/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_;
 
diff --git a/forte/mrdsrg-spin-adapted/sa_ldsrg2.cc b/forte/mrdsrg-spin-adapted/sa_ldsrg2.cc
index 1479a93dd..c604b250a 100644
--- a/forte/mrdsrg-spin-adapted/sa_ldsrg2.cc
+++ b/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);
 
diff --git a/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc b/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc
index 68e7f7343..6c0d08a72 100644
--- a/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc
+++ b/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;
 }
 
diff --git a/forte/mrdsrg-spin-adapted/sa_mrdsrg.h b/forte/mrdsrg-spin-adapted/sa_mrdsrg.h
index 3b4869e83..236bc0af0 100644
--- a/forte/mrdsrg-spin-adapted/sa_mrdsrg.h
+++ b/forte/mrdsrg-spin-adapted/sa_mrdsrg.h
@@ -113,10 +113,6 @@ class SA_MRDSRG : public SADSRG {
     ambit::BlockedTensor B_;
     /// Generalized Fock matrix
     ambit::BlockedTensor F_;
-    /// Single excitation amplitude
-    ambit::BlockedTensor T1_;
-    /// Double excitation amplitude
-    ambit::BlockedTensor T2_;
     /// Difference of consecutive singles
     ambit::BlockedTensor DT1_;
     /// Difference of consecutive doubles
@@ -153,13 +149,31 @@ class SA_MRDSRG : public SADSRG {
     void guess_t1(BlockedTensor& F, BlockedTensor& T2, BlockedTensor& T1);
     /// Update T1 in every iteration
     void update_t1();
+    /// Update T1 using projective conditions
+    void update_t1_proj();
+
+    /// Denominator of T1 core-actv block in orthogonalized basis
+    ambit::Tensor Aca_;
+    /// Residual of T1 core-actv block update in orthogonalized basis
+    ambit::Tensor Oca_;
+    /// Denominator of T1 actv-virt block in orthogonalized basis
+    ambit::Tensor Aav_;
+    /// Residual of T1 actv-virt block update in orthogonalized basis
+    ambit::Tensor Oav_;
+
+    /// Compute core-actv block denominator for projective T1 update
+    void compute_proj_denom_ca();
+    /// Compute actv-virt block denominator for projective T1 update
+    void compute_proj_denom_av();
 
     /// Compute DSRG-transformed Hamiltonian
     void compute_hbar();
-    /// Compute DSRG-transformed Hamiltonian Hbar sequentially
+    /// Compute DSRG-transformed Hamiltonian sequentially
     void compute_hbar_sequential();
     /// Compute DSRG-transformed Hamiltonian truncated to 2-nested commutator
     void compute_hbar_qc();
+    /// Compute DSRG-transformed Hamiltonian truncated to 2-nested commutator sequentially
+    void compute_hbar_qc_sequential();
 
     /// Add H2's Hermitian conjugate to itself, H2 need to contain gGgG block
     void add_hermitian_conjugate(BlockedTensor& H2);
diff --git a/forte/mrdsrg-spin-adapted/sa_mrdsrg_amps.cc b/forte/mrdsrg-spin-adapted/sa_mrdsrg_amps.cc
index 086da47e4..8c08de94b 100644
--- a/forte/mrdsrg-spin-adapted/sa_mrdsrg_amps.cc
+++ b/forte/mrdsrg-spin-adapted/sa_mrdsrg_amps.cc
@@ -51,8 +51,13 @@ void SA_MRDSRG::guess_t(BlockedTensor& V, BlockedTensor& T2, BlockedTensor& F, B
 }
 
 void SA_MRDSRG::update_t() {
-    update_t2();
-    update_t1();
+    if (t1_type_ == "PROJECT") {
+        update_t1_proj();
+        update_t2();
+    } else {
+        update_t2();
+        update_t1();
+    }
 }
 
 void SA_MRDSRG::guess_t2(BlockedTensor& V, BlockedTensor& T2, BlockedTensor& B) {
@@ -418,6 +423,168 @@ void SA_MRDSRG::update_t1() {
     Hbar1_["uv"] = Hbar1copy["uv"];
 }
 
+void SA_MRDSRG::update_t1_proj() {
+    /// update T1 according to <{a^i_a} Hbar^a_i> = 0
+    DT1_["ia"] = T1_["ia"];
+
+    // core-virt block
+    const auto& Hbar1cv_data = Hbar1_.block("cv").data();
+    auto dim_v = Hbar1_.block("cv").dim(1);
+    T1_.block("cv").iterate([&](const std::vector<size_t>& i, double& value) {
+        auto i0 = label_to_spacemo_['c'][i[0]];
+        auto i1 = label_to_spacemo_['v'][i[1]];
+        value += Hbar1cv_data[i[0] * dim_v + i[1]] / (Fdiag_[i0] - Fdiag_[i1]);
+    });
+
+    // core-actv block
+    // 1. residual in original basis
+    auto Nca = ambit::BlockedTensor::build(tensor_type_, "Nca", {"ca"});
+    Nca["mz"] += 0.5 * Hbar1_["mu"] * Eta1_["uz"];
+    Nca["mz"] -= 0.5 * Hbar2_["xmuv"] * L2_["uvxz"];
+
+    // 2. transform residual to orthogonal basis
+    Oca_("mz") = Nca.block("ca")("mw") * Xca_("wz");
+
+    // 3. update in orthogonal basis
+    auto Oca = Oca_.clone();
+    Oca_("mz") = Oca("mz") * Aca_("mz");
+    // auto& Oca_data = Oca_.data();
+    // const auto& Aca_data = Aca_.data();
+    // std::transform(Oca_data.begin(), Oca_data.end(), Aca_data.begin(), Oca_data.begin(),
+    //                [](double n, double d) { return n * d; });
+    // Oca_.print();
+
+    // 4. transform update to original basis
+    T1_.block("ca")("mu") -= Oca_("mz") * Xca_("uz");
+
+    // actv-virt block
+    // 1. residual in original basis
+    auto Nav = ambit::BlockedTensor::build(tensor_type_, "Nav", {"av"});
+    Nav["ze"] += 0.5 * Hbar1_["ue"] * L1_["zu"];
+    Nav["ze"] += 0.5 * Hbar2_["uvex"] * L2_["zxuv"];
+
+    // 2. transform residual to orthogonal basis
+    Oav_("ze") = Nav.block("av")("we") * Xav_("wz");
+
+    // 3. update in orthogonal basis
+    auto Oav = Oav_.clone();
+    Oav_("ze") = Oav("ze") * Aav_("ze");
+    // auto& Oav_data = Oav_.data();
+    // const auto& Aav_data = Aav_.data();
+    // std::transform(Oav_data.begin(), Oav_data.end(), Aav_data.begin(), Oav_data.begin(),
+    //                [](double n, double d) { return n * d; });
+    // Oav_.print();
+
+    // 4. transform update to original basis
+    T1_.block("av")("ue") -= Oav_("ze") * Xav_("uz");
+
+    // compute difference
+    DT1_["ia"] -= T1_["ia"];
+    DT1_.scale(-1.0);
+    T1rms_ = DT1_.norm();
+    // DT1_.print();
+
+    // compute norm and find maximum
+    T1max_ = T1_.norm(0);
+    T1norm_ = T1_.norm(2);
+}
+
+void SA_MRDSRG::compute_proj_denom_ca() {
+    // compute diagonal elements of Jacobian in orthogonal basis
+    // <{m^+ v} [H, {u^+ n} - {n^+ u}]> X_uu' X_vu' δ_mn
+
+    auto I = ambit::BlockedTensor::build(tensor_type_, "I", {"cc"});
+    I.iterate([&](const std::vector<size_t>& i, const std::vector<SpinType>&, double& value) {
+        if (i[0] == i[1]) {
+            value = 1.0;
+        }
+    });
+
+    // TODO: need to optimize
+    // 1. directly form caa
+    // 2. density fitting integrals
+    auto Dccaa = ambit::BlockedTensor::build(tensor_type_, "Dccaa", {"ccaa"});
+
+    Dccaa["mnuv"] += 0.5 * F_["vx"] * Eta1_["xu"] * I["nm"];
+    Dccaa["mnuv"] -= 0.5 * F_["nm"] * Eta1_["vu"];
+    Dccaa["mnuv"] += 0.25 * V_["xnmy"] * Eta1_["vx"] * Eta1_["yu"];
+    Dccaa["mnuv"] -= 0.25 * V_["nxmy"] * Eta1_["vx"] * Eta1_["yu"];
+    Dccaa["mnuv"] -= 0.5 * V_["vwxy"] * L2_["xyuw"] * I["mn"];
+    Dccaa["mnuv"] += 0.5 * V_["nxmy"] * L2_["vyux"];
+    Dccaa["mnuv"] -= (1.0 / 6.0) * V_["xnmy"] * L2_["vyux"];
+    Dccaa["mnuv"] += (1.0 / 6.0) * V_["xnmy"] * L2_["yvux"];
+
+    // only keep for m = n
+    auto ncore = core_mos_.size();
+    auto nactv = actv_mos_.size();
+    auto Dcaa = ambit::Tensor::build(tensor_type_, "Dcaa", {ncore, nactv, nactv});
+
+    const auto& Dccaa_data = Dccaa.block("ccaa").data();
+    auto ccaa_dim3 = ncore * nactv * nactv;
+    auto ccaa_dim2 = nactv * nactv;
+
+    Dcaa.iterate([&](const std::vector<size_t>& i, double& value) {
+        value = Dccaa_data[i[0] * ccaa_dim3 + i[0] * ccaa_dim2 + i[1] * nactv + i[2]];
+    });
+
+    // transform to orthogonal T1 basis
+    auto nactv_h = Xca_.dim(1);
+    Aca_ = ambit::Tensor::build(tensor_type_, "Inverse of Diag. Jacobian ca", {ncore, nactv_h});
+    Aca_("mx") = Dcaa("muv") * Xca_("ux") * Xca_("vx");
+
+    // TODO: decide we store inverse or the actual denominator
+    Aca_.iterate([&](const std::vector<size_t>& i, double& value) { value = 1.0 / value; });
+    // Aca_.print();
+}
+
+void SA_MRDSRG::compute_proj_denom_av() {
+    // compute diagonal elements of Jacobian in orthogonal basis
+    // <{v^+ f} [H, {e^+ u} - {u^+ e}]> X_uu' X_vu' δ_ef
+
+    auto I = ambit::BlockedTensor::build(tensor_type_, "I", {"vv"});
+    I.iterate([&](const std::vector<size_t>& i, const std::vector<SpinType>&, double& value) {
+        if (i[0] == i[1]) {
+            value = 1.0;
+        }
+    });
+
+    // TODO: need to optimize
+    // 1. directly form vaa
+    // 2. density fitting integrals
+    auto Dvvaa = ambit::BlockedTensor::build(tensor_type_, "Dvvaa", {"vvaa"});
+
+    Dvvaa["efuv"] += 0.5 * F_["ef"] * L1_["vu"];
+    Dvvaa["efuv"] -= 0.5 * F_["xu"] * L1_["vx"] * I["ef"];
+    Dvvaa["efuv"] -= 0.25 * V_["eyfx"] * L1_["xu"] * L1_["vy"];
+    Dvvaa["efuv"] += 0.25 * V_["yefx"] * L1_["xu"] * L1_["vy"];
+    Dvvaa["efuv"] -= 0.5 * V_["xyuw"] * L2_["vwxy"] * I["ef"];
+    Dvvaa["efuv"] += 0.5 * V_["exfy"] * L2_["vyux"];
+    Dvvaa["efuv"] -= (1.0 / 6.0) * V_["xefy"] * L2_["vyux"];
+    Dvvaa["efuv"] += (1.0 / 6.0) * V_["xefy"] * L2_["yvux"];
+
+    // only keep for e = f
+    auto nvirt = virt_mos_.size();
+    auto nactv = actv_mos_.size();
+    auto Dvaa = ambit::Tensor::build(tensor_type_, "Dvaa", {nvirt, nactv, nactv});
+
+    const auto& Dvvaa_data = Dvvaa.block("vvaa").data();
+    auto vvaa_dim3 = nvirt * nactv * nactv;
+    auto vvaa_dim2 = nactv * nactv;
+
+    Dvaa.iterate([&](const std::vector<size_t>& i, double& value) {
+        value = Dvvaa_data[i[0] * vvaa_dim3 + i[0] * vvaa_dim2 + i[1] * nactv + i[2]];
+    });
+
+    // transform to orthogonal T1 basis
+    auto nactv_p = Xav_.dim(1);
+    Aav_ = ambit::Tensor::build(tensor_type_, "Inverse of Diag. Jacobian av", {nactv_p, nvirt});
+    Aav_("xe") = Dvaa("euv") * Xav_("ux") * Xav_("vx");
+
+    // TODO: decide we store inverse or the actual denominator
+    Aav_.iterate([&](const std::vector<size_t>& i, double& value) { value = 1.0 / value; });
+    // Aav_.print();
+}
+
 void SA_MRDSRG::dump_amps_to_disk() {
     // dump to psi4 scratch directory for reference relaxation
     if (restart_amps_ and (relax_ref_ != "NONE")) {
diff --git a/forte/mrdsrg-spin-adapted/sa_mrpt2.cc b/forte/mrdsrg-spin-adapted/sa_mrpt2.cc
index 85d4843e0..15fd4ba5b 100644
--- a/forte/mrdsrg-spin-adapted/sa_mrpt2.cc
+++ b/forte/mrdsrg-spin-adapted/sa_mrpt2.cc
@@ -207,7 +207,7 @@ void SA_MRPT2::init_amps() {
     timer t("Initialize T1 and T2");
     print_contents("Allocating amplitudes");
 
-    T1_ = BTF_->build(tensor_type_, "T1 Amplitudes", {"hp"});
+//    T1_ = BTF_->build(tensor_type_, "T1 Amplitudes", {"hp"});
 
     if (eri_df_) {
         std::vector<std::string> blocks{"aavv", "ccaa", "caav", "acav", "aava", "caaa", "aaaa"};
diff --git a/forte/mrdsrg-spin-adapted/sa_mrpt3.cc b/forte/mrdsrg-spin-adapted/sa_mrpt3.cc
index 7351995e2..41aa4a45d 100644
--- a/forte/mrdsrg-spin-adapted/sa_mrpt3.cc
+++ b/forte/mrdsrg-spin-adapted/sa_mrpt3.cc
@@ -100,7 +100,7 @@ void SA_MRPT3::build_ints() {
 void SA_MRPT3::init_amps() {
     timer t("Initialize T1 and T2");
     print_contents("Allocating 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"});
     S2_ = BTF_->build(tensor_type_, "S2 Amplitudes", {"hhpp"});
     print_done(t.stop());
diff --git a/forte/mrdsrg-spin-adapted/sadsrg.cc b/forte/mrdsrg-spin-adapted/sadsrg.cc
index eb107eb4e..994f6b4e4 100644
--- a/forte/mrdsrg-spin-adapted/sadsrg.cc
+++ b/forte/mrdsrg-spin-adapted/sadsrg.cc
@@ -138,6 +138,12 @@ void SADSRG::startup() {
     chk_filename_prefix_ += "forte." + std::to_string(getpid());
     chk_filename_prefix_ += "." + psi::Process::environment.molecule()->name();
     Bcan_files_.clear();
+
+    // all methods need T1 amplitudes
+    T1_ = BTF_->build(tensor_type_, "T1 Amplitudes", {"hp"});
+    if (t1_type_ == "PROJECT") {
+        build_transformations_orthogonal_t1();
+    }
 }
 
 void SADSRG::build_fock_from_ints() {
@@ -199,6 +205,12 @@ void SADSRG::read_options() {
 
     relax_ref_ = foptions_->get_str("RELAX_REF");
 
+    brueckner_ = foptions_->get_bool("DSRG_BRUECKNER");
+    brueckner_conv_ = foptions_->get_double("BRUECKNER_CONVERGENCE");
+
+    t1_type_ = foptions_->get_str("DSRG_T1_AMPS_TYPE");
+    t1_proj_cutoff_ = foptions_->get_double("DSRG_T1_AMPS_PROJ_CUTOFF");
+
     multi_state_ = foptions_->get_gen_list("AVG_STATE").size() != 0;
     multi_state_algorithm_ = foptions_->get_str("DSRG_MULTI_STATE");
 
@@ -851,6 +863,90 @@ void SADSRG::print_cumulant_summary() {
     outfile->Printf("\n    %s", dash.c_str());
 }
 
+void SADSRG::build_transformations_orthogonal_t1() {
+    auto nirrep = mo_space_info_->nirrep();
+    auto nactv = mo_space_info_->size("ACTIVE");
+    auto dim_actv = mo_space_info_->dimension("ACTIVE");
+
+    // put 1RDM in psi4 matrix form
+    auto& L1data = L1_.block("aa").data();
+    auto L1mat = std::make_shared<psi::Matrix>("1RDM", dim_actv, dim_actv);
+    for (size_t h = 0, offset = 0; h < nirrep; ++h) {
+        for (size_t p = 0; p < dim_actv[h]; ++p) {
+            for (size_t q = 0; q < dim_actv[h]; ++q) {
+                L1mat->set(h, p, q, 0.5 * L1data[(p + offset) * nactv + q + offset]);
+            }
+        }
+        offset += dim_actv[h];
+    }
+
+    // diagonalize 1RDM
+    auto U = std::make_shared<psi::Matrix>("U NO", dim_actv, dim_actv);
+    auto s = std::make_shared<psi::Vector>("s", dim_actv);
+    L1mat->diagonalize(U, s, descending);
+
+    // apply linear dependency cutoff
+    psi::Dimension rank_part(nirrep), rank_hole(nirrep);
+    for (int h = 0; h < nirrep; ++h) {
+        if (!dim_actv[h])
+            continue;
+        auto keep_part = 0, keep_hole = 0;
+        for (int i = 0; i < dim_actv[h]; ++i) {
+            keep_part += (s->get(h, i) > t1_proj_cutoff_ ? 1 : 0);
+            keep_hole += ((1.0 - s->get(h, i)) > t1_proj_cutoff_ ? 1 : 0);
+        }
+        rank_part[h] = keep_part;
+        rank_hole[h] = keep_hole;
+    }
+
+    // compute transformation matrix X = U s^-1/2
+    auto Xp = std::make_shared<psi::Matrix>("Xp = U n^-1/2", dim_actv, rank_part);
+    for (int h = 0; h < nirrep; ++h) {
+        if (!dim_actv[h] or !rank_part[h])
+            continue;
+        for (int i = 0; i < rank_part[h]; ++i) {
+            auto col = U->get_column(h, i);
+            col->scale(pow(s->get(h, i), -0.5));
+            Xp->set_column(h, i, col);
+        }
+    }
+
+    auto Xh = std::make_shared<psi::Matrix>("Xh = U (1 - n)^-1/2", dim_actv, rank_hole);
+    for (int h = 0; h < nirrep; ++h) {
+        if (!dim_actv[h] or !rank_hole[h])
+            continue;
+        for (int i = 0; i < rank_hole[h]; ++i) {
+            auto col = U->get_column(h, dim_actv[h] - i - 1);
+            col->scale(pow(1.0 - s->get(h, dim_actv[h] - i - 1), -0.5));
+            Xh->set_column(h, i, col);
+        }
+    }
+
+    // put transformation matrix in ambit form
+    size_t psize = rank_part.sum();
+    size_t hsize = rank_hole.sum();
+
+    Xav_ = ambit::Tensor::build(tensor_type_, "Xav", {nactv, psize});
+    Xca_ = ambit::Tensor::build(tensor_type_, "Xca", {nactv, hsize});
+
+    auto& Xav_data = Xav_.data();
+    auto& Xca_data = Xca_.data();
+
+    for (size_t h = 0, offset_p = 0, offset_part = 0, offset_hole = 0; h < nirrep; ++h) {
+        for (size_t p = 0; p < dim_actv[h]; ++p) {
+            for (size_t q = 0; q < rank_part[h]; ++q) {
+                Xav_data[(p + offset_p) * psize + q + offset_part] = Xp->get(h, p, q);
+            }
+            for (size_t q = 0; q < rank_hole[h]; ++q) {
+                Xca_data[(p + offset_p) * hsize + q + offset_part] = Xh->get(h, p, q);
+            }
+        }
+        offset_p += dim_actv[h];
+        offset_part += rank_part[h];
+        offset_hole += rank_hole[h];
+    }
+}
+
 std::vector<double> SADSRG::diagonalize_Fock_diagblocks(BlockedTensor& U) {
     // set U to identity and output diagonal Fock
     U.iterate([&](const std::vector<size_t>& i, const std::vector<SpinType>&, double& value) {
@@ -903,6 +999,8 @@ std::vector<double> SADSRG::diagonalize_Fock_diagblocks(BlockedTensor& U) {
         auto Usub = std::make_shared<psi::Matrix>("U " + space, dim, dim);
         auto evals = std::make_shared<psi::Vector>("evals " + space, dim);
         Fd->diagonalize(Usub, evals);
+        outfile->Printf("\n  space: %s\n", space.c_str());
+        evals->print();
 
         auto& Udata = U.block(block).data();
         auto corr_abs_indices = mo_space_info_->corr_absolute_mo(space);
diff --git a/forte/mrdsrg-spin-adapted/sadsrg.h b/forte/mrdsrg-spin-adapted/sadsrg.h
index 6bb07bf55..d6a6d6a07 100644
--- a/forte/mrdsrg-spin-adapted/sadsrg.h
+++ b/forte/mrdsrg-spin-adapted/sadsrg.h
@@ -70,6 +70,9 @@ class SADSRG : public DynamicCorrelationSolver {
     /// Set unitary matrix (in active space) from original to semicanonical
     void set_Uactv(ambit::Tensor& U);
 
+    /// Return if DSRG Brueckner orbitals are converged
+    bool is_brueckner_converged() { return brueckner_absmax_ < brueckner_conv_; }
+
   protected:
     /// Startup function called in constructor
     void startup();
@@ -179,11 +182,6 @@ class SADSRG : public DynamicCorrelationSolver {
     /// List of the symmetry of the active MOs
     std::vector<int> actv_mos_sym_;
 
-    /// List of active active occupied MOs (relative to active)
-    std::vector<size_t> actv_occ_mos_;
-    /// List of active active unoccupied MOs (relative to active)
-    std::vector<size_t> actv_uocc_mos_;
-
     /// List of auxiliary MOs when DF/CD
     std::vector<size_t> aux_mos_;
 
@@ -257,9 +255,9 @@ class SADSRG : public DynamicCorrelationSolver {
     /// Fill in diagonal elements of Fock matrix to Fdiag
     void fill_Fdiag(BlockedTensor& F, std::vector<double>& Fdiag);
 
-    /// Check orbitals if semicanonical
+    /// Check orbitals if semi-canonical
     bool check_semi_orbs();
-    /// Are orbitals semi-canonicalized?
+    /// Are orbitals semi-canonical?
     bool semi_canonical_;
     /// Checked results of each block of Fock matrix
     std::map<std::string, bool> semi_checked_results_;
@@ -400,7 +398,7 @@ class SADSRG : public DynamicCorrelationSolver {
     /// Compute zero-body term of commutator [V, T1], V is constructed from B (DF/CD)
     void V_T1_C0_DF(BlockedTensor& B, BlockedTensor& T1, const double& alpha, double& C0);
     /// Compute zero-body term of commutator [V, T2], V is constructed from B (DF/CD)
-    std::vector<double> V_T2_C0_DF(BlockedTensor& B, BlockedTensor& T1, BlockedTensor& S2,
+    std::vector<double> V_T2_C0_DF(BlockedTensor& B, BlockedTensor& T2, BlockedTensor& S2,
                                    const double& alpha, double& C0);
 
     /// Compute one-body term of commutator [V, T1], V is constructed from B (DF/CD)
@@ -471,8 +469,41 @@ class SADSRG : public DynamicCorrelationSolver {
     /// Print done and timing
     void print_done(double t, const std::string& done="Done");
 
+    // ==> orbital rotations <==
+
+    /// Perform orbital rotations using an unitary matrix
+    void brueckner_orbital_rotation(ambit::BlockedTensor T1);
+
+    /// Whether perform orbital rotation to make T1 vanishing
+    bool brueckner_;
+    /// Convergence threshold for Brueckner orbitals
+    double brueckner_conv_;
+    /// Max element of T1 to compare against Brueckner convergence threshold
+    double brueckner_absmax_;
+
+    /// Compute the exponential of exp(T1 - T1^+) and return a Psi4 SharedMatrix
+    /// @param T1 the T1 matrix
+    /// @param with_symmetry if in blocked form for the returned SharedMatrix
+    psi::SharedMatrix expA1(ambit::BlockedTensor T1, bool with_symmetry);
+
     // ==> common amplitudes analysis and printing <==
 
+    /// Condition for T1 amplitudes
+    std::string t1_type_;
+    /// Cutoff to remove internally contracted singles
+    double t1_proj_cutoff_;
+    /// Transformation matrix to orthogonalize T1 (core-actv)
+    ambit::Tensor Xca_;
+    /// Transformation matrix to orthogonalize T1 (actv-virt)
+    ambit::Tensor Xav_;
+    /// Build transformation matrix to orthogonalize T1
+    void build_transformations_orthogonal_t1();
+
+    /// Single excitation amplitudes
+    ambit::BlockedTensor T1_;
+    /// Double excitation amplitudes
+    ambit::BlockedTensor T2_;
+
     /// Prune internal amplitudes for T1
     void internal_amps_T1(BlockedTensor& T1);
     /// Prune internal amplitudes for T2
diff --git a/forte/proc/dsrg.py b/forte/proc/dsrg.py
index 3b8a72e5d..0fe7e9956 100644
--- a/forte/proc/dsrg.py
+++ b/forte/proc/dsrg.py
@@ -89,6 +89,10 @@ def __init__(self, active_space_solver, state_weights_map, mo_space_info, ints,
 
         self.save_relax_energies = options.get_bool("DSRG_DUMP_RELAXED_ENERGIES")
 
+        # Orbital rotations
+        self.do_brueckner = options.get_bool("DSRG_BRUECKNER")
+        self.brueckner_maxiter = options.get_int("BRUECKNER_MAXITER")
+
         # Filter out some ms-dsrg algorithms
         ms_dsrg_algorithm = options.get_str("DSRG_MULTI_STATE")
         if self.do_multi_state and ("SA" not in ms_dsrg_algorithm):
@@ -198,8 +202,11 @@ def compute_energy(self):
         self.make_dsrg_solver()
         self.dsrg_setup()
         e_dsrg = self.dsrg_solver.compute_energy()
-        psi4.core.set_scalar_variable("UNRELAXED ENERGY", e_dsrg)
 
+        if self.do_brueckner:
+            e_dsrg = self.compute_brueckner_iterations()
+
+        psi4.core.set_scalar_variable("UNRELAXED ENERGY", e_dsrg)
         self.energies_environment[0] = {k: v for k, v in psi4.core.variables().items() if 'ROOT' in k}
 
         # Spit out energy if reference relaxation not implemented
@@ -376,7 +383,7 @@ def reorder_weights(self, state_ci_wfn_map):
             overlap_np = np.abs(overlap.to_array())
             max_values = np.max(overlap_np, axis=1)
             permutation = np.argmax(overlap_np, axis=1)
-            check_pass = len(permutation) == len(set(permutation)) and np.all(max_values > 0.5)
+            check_pass = len(permutation) == len(set(permutation)) and np.all(max_values > 0.8)
 
             if not check_pass:
                 msg = "Relaxed states are likely wrong. Please increase the number of roots."
@@ -479,3 +486,38 @@ def push_to_psi4_environment(self):
 
                     if self.do_dipole and (not self.do_multi_state):
                         psi4.core.set_scalar_variable('FULLY RELAXED DIPOLE', self.dipoles[-1][1][-1])
+
+    def compute_brueckner_iterations(self):
+        """ Iterations to obtain DSRG Brueckner orbitals. """
+        e_dsrg = 0.0
+        converged = False
+
+        for i in range(self.brueckner_maxiter):
+            # form new active-space integrals using rotated ForteIntegrals
+            as_ints = forte.make_active_space_ints(self.mo_space_info, self.ints, "ACTIVE", ["RESTRICTED_DOCC"])
+
+            # solve for active space
+            self.active_space_solver.set_active_space_integrals(as_ints)
+            state_energies_list = self.active_space_solver.compute_energy()
+            forte.compute_average_state_energy(state_energies_list, self.state_weights_map)
+            self.rdms = self.active_space_solver.compute_average_rdms(self.state_weights_map, self.max_rdm_level,
+                                                                      self.rdm_type)
+
+            # semi-canonicalize orbitals
+            if self.do_semicanonical:
+                self.semi.semicanonicalize(self.rdms)
+
+            # solve for DSRG
+            self.make_dsrg_solver()
+            self.dsrg_setup()
+            e_dsrg = self.dsrg_solver.compute_energy()
+
+            if self.dsrg_solver.is_brueckner_converged():
+                psi4.core.print_out("\n  DSRG orbital update converged.")
+                converged = True
+                break
+
+        if not converged:
+            raise psi4.p4util.PsiException(f"DSRG orbital update did not converge in {self.brueckner_maxiter} cycles!")
+
+        return e_dsrg
diff --git a/forte/register_forte_options.py b/forte/register_forte_options.py
index ee4159b63..016adcc75 100644
--- a/forte/register_forte_options.py
+++ b/forte/register_forte_options.py
@@ -591,7 +591,7 @@ def register_dsrg_options(options):
 
     options.add_str(
         "CORR_LEVEL", "PT2",
-        ["PT2", "PT3", "LDSRG2", "LDSRG2_QC", "LSRG2", "SRG_PT2", "QDSRG2", "LDSRG2_P3", "QDSRG2_P3"],
+        ["PT2", "PT3", "CC2", "LDSRG2", "LDSRG2_QC", "LSRG2", "SRG_PT2", "QDSRG2", "LDSRG2_P3", "QDSRG2_P3"],
         "Correlation level of MR-DSRG (used in mrdsrg code, "
         "LDSRG2_P3 and QDSRG2_P3 not implemented)"
     )
@@ -770,6 +770,21 @@ def register_dsrg_options(options):
 
     options.add_bool("DSRG_RDM_MS_AVG", False, "Form Ms-averaged density if true")
 
+    options.add_str("DSRG_T1_AMPS_TYPE", "MANY_BODY", ["MANY_BODY", "PROJECT"],
+                    "T1 conditions for MRDSRG")
+
+    options.add_double("DSRG_T1_AMPS_PROJ_CUTOFF", 0.05,
+                       "Cutoff for removal linear dependencies of internally contracted singles")
+
+    options.add_bool("DSRG_BRUECKNER", False,
+                     "Rotate orbitals such that MRDSRG T1 amplitudes become zero")
+
+    options.add_double("BRUECKNER_CONVERGENCE", 1.0e-6,
+                       "Threshold to consider Brueckner orbitals converged")
+
+    options.add_int("BRUECKNER_MAXITER", 40,
+                    "The max number of iterations for Brueckner orbital update")
+
 
 def register_dwms_options(options):
     options.set_group("DWMS")