diff --git a/forte/mrdsrg-so/mrdsrg_so.cc b/forte/mrdsrg-so/mrdsrg_so.cc
index ee7109073..9931c51e1 100644
--- a/forte/mrdsrg-so/mrdsrg_so.cc
+++ b/forte/mrdsrg-so/mrdsrg_so.cc
@@ -108,14 +108,6 @@ void MRDSRG_SO::startup() {
 
     source_ = foptions_->get_str("SOURCE");
 
-    dsrg_trans_type_ = foptions_->get_str("DSRG_TRANS_TYPE");
-    if (dsrg_trans_type_ == "CC" && foptions_->get_str("CORR_LEVEL") == "QDSRG2") {
-        outfile->Printf(
-            "\n  Warning: DSRG_TRANS_TYPE option CC is not supported with CORR_LEVEL QDSRG2.");
-        outfile->Printf("\n  Changed DSRG_TRANS_TYPE option to UNITARY");
-        dsrg_trans_type_ = "UNITARY";
-    }
-
     ntamp_ = foptions_->get_int("NTAMP");
     intruder_tamp_ = foptions_->get_double("INTRUDER_TAMP");
 
@@ -775,8 +767,6 @@ void MRDSRG_SO::compute_hbar() {
     BlockedTensor C1 = ambit::BlockedTensor::build(tensor_type_, "C1", {"gg"});
     BlockedTensor C2 = ambit::BlockedTensor::build(tensor_type_, "C2", {"gggg"});
 
-    bool do_wicked = foptions_->get_bool("DO_WICKED");
-
     // compute Hbar recursively
     for (int n = 1; n <= maxn; ++n) {
         // prefactor before n-nested commutator
@@ -786,6 +776,7 @@ void MRDSRG_SO::compute_hbar() {
         double C0 = 0.0;
         C1.zero();
         C2.zero();
+
         // zero-body
         H1_T1_C0(O1, T1, factor, C0);
         H1_T2_C0(O1, T2, factor, C0);
@@ -808,21 +799,12 @@ void MRDSRG_SO::compute_hbar() {
         //        outfile->Printf("\n  |H2| = %20.12f", C2.norm(1));
         //        outfile->Printf("\n  --------------------------------");
 
-        if (dsrg_trans_type_ == "UNITARY") {
-            // [H, A] = [H, T] + [H, T]^dagger
-            C0 *= 2.0;
-            O1["pq"] = C1["pq"];
-            C1["pq"] += O1["qp"];
-            O2["pqrs"] = C2["pqrs"];
-            C2["pqrs"] += O2["rspq"];
-        }
-
-        if (do_wicked) {
-            O2["pqrs"] = C2["pqrs"];
-            C2["pqrs"] -= O2["qprs"];
-            C2["pqrs"] -= O2["pqsr"];
-            C2["pqrs"] += O2["qpsr"];
-        }
+        // [H, A] = [H, T] + [H, T]^dagger
+        C0 *= 2.0;
+        O1["pq"] = C1["pq"];
+        C1["pq"] += O1["qp"];
+        O2["pqrs"] = C2["pqrs"];
+        C2["pqrs"] += O2["rspq"];
 
         // Hbar += C
         Hbar0 += C0;
diff --git a/forte/mrdsrg-so/mrdsrg_so.h b/forte/mrdsrg-so/mrdsrg_so.h
index da6388a6b..927b126fe 100644
--- a/forte/mrdsrg-so/mrdsrg_so.h
+++ b/forte/mrdsrg-so/mrdsrg_so.h
@@ -120,9 +120,6 @@ class MRDSRG_SO : public DynamicCorrelationSolver {
     /// Order of the Taylor expansion of f(z) = (1-exp(-z^2))/z
     int taylor_order_;
 
-    /// DSRG transformation type
-    std::string dsrg_trans_type_;
-
     std::shared_ptr<BlockedTensorFactory> BTF_;
     TensorType tensor_type_;
 
@@ -213,11 +210,6 @@ class MRDSRG_SO : public DynamicCorrelationSolver {
     void H1_T2_C0(BlockedTensor& H1, BlockedTensor& T2, const double& alpha, double& C0);
     void H2_T2_C0(BlockedTensor& H2, BlockedTensor& T2, const double& alpha, double& C0);
 
-    void H1_T1_C0(BlockedTensor& H1, BlockedTensor& T1, const double& alpha, double& C0);
-    void H2_T1_C0(BlockedTensor& H2, BlockedTensor& T1, const double& alpha, double& C0);
-    void H1_T2_C0(BlockedTensor& H1, BlockedTensor& T2, const double& alpha, double& C0);
-    void H2_T2_C0(BlockedTensor& H2, BlockedTensor& T2, const double& alpha, double& C0);
-
     /// Compute one-body term of commutator [H, T]
     void H1_T1_C1(BlockedTensor& H1, BlockedTensor& T1, const double& alpha, BlockedTensor& C1);
     void H1_T2_C1(BlockedTensor& H1, BlockedTensor& T2, const double& alpha, BlockedTensor& C1);
@@ -226,11 +218,6 @@ class MRDSRG_SO : public DynamicCorrelationSolver {
     void H3_T1_C1(BlockedTensor& H3, BlockedTensor& T1, const double& alpha, BlockedTensor& C1);
     void H3_T2_C1(BlockedTensor& H3, BlockedTensor& T2, const double& alpha, BlockedTensor& C1);
 
-    void H1_T1_C1(BlockedTensor& H1, BlockedTensor& T1, const double& alpha, BlockedTensor& C1);
-    void H1_T2_C1(BlockedTensor& H1, BlockedTensor& T2, const double& alpha, BlockedTensor& C1);
-    void H2_T1_C1(BlockedTensor& H2, BlockedTensor& T1, const double& alpha, BlockedTensor& C1);
-    void H2_T2_C1(BlockedTensor& H2, BlockedTensor& T2, const double& alpha, BlockedTensor& C1);
-
     /// Compute two-body term of commutator [H, T]
     void H2_T1_C2(BlockedTensor& H2, BlockedTensor& T1, const double& alpha, BlockedTensor& C2);
     void H1_T2_C2(BlockedTensor& H1, BlockedTensor& T2, const double& alpha, BlockedTensor& C2);
@@ -238,10 +225,6 @@ class MRDSRG_SO : public DynamicCorrelationSolver {
     void H3_T1_C2(BlockedTensor& H3, BlockedTensor& T1, const double& alpha, BlockedTensor& C2);
     void H3_T2_C2(BlockedTensor& H3, BlockedTensor& T2, const double& alpha, BlockedTensor& C2);
 
-    void H2_T1_C2(BlockedTensor& H2, BlockedTensor& T1, const double& alpha, BlockedTensor& C2);
-    void H1_T2_C2(BlockedTensor& H1, BlockedTensor& T2, const double& alpha, BlockedTensor& C2);
-    void H2_T2_C2(BlockedTensor& H2, BlockedTensor& T2, const double& alpha, BlockedTensor& C2);
-
     /// Compute three-body term of commutator [H, T]
     void H2_T2_C3(BlockedTensor& H2, BlockedTensor& T2, const double& alpha, BlockedTensor& C3);
 
diff --git a/forte/proc/dsrg.py b/forte/proc/dsrg.py
index 3fc96cb7d..b0517aa0a 100644
--- a/forte/proc/dsrg.py
+++ b/forte/proc/dsrg.py
@@ -72,7 +72,13 @@ def __init__(self, active_space_solver, state_weights_map, mo_space_info, ints,
                 "\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"]: