diff --git a/pyblock2/driver/core.py b/pyblock2/driver/core.py
index b450db32..1be55f1d 100644
--- a/pyblock2/driver/core.py
+++ b/pyblock2/driver/core.py
@@ -154,6 +154,15 @@ def get_mpo(self, expr, iprint=0):
         mpo = bw.bs.SimplifiedMPO(mpo, bw.bs.Rule(), False, False)
         return mpo
 
+    def orbital_reordering(self, h1e, g2e):
+        bw = self.bw
+        import numpy as np
+        xmat = np.abs(np.einsum('ijji->ij', g2e, optimize=True))
+        kmat = np.abs(h1e) * 1E-7 + xmat
+        kmat = bw.b.VectorDouble(kmat.flatten())
+        idx = bw.b.OrbitalOrdering.fiedler(len(h1e), kmat)
+        return np.array(idx)
+
     def dmrg(
         self,
         mpo,
@@ -176,7 +185,7 @@ def dmrg(
                 thrds = [1e-6] * 4 + [1e-7] * 1
             else:
                 thrds = [1e-5] * 4 + [5e-6] * 1
-        if dav_type is not None and "NonHermitian" in dav_type:
+        if dav_type is not None and "LeftEigen" in dav_type:
             bra = ket.deep_copy("BRA")
         else:
             bra = ket
diff --git a/src/core/iterative_matrix_functions.hpp b/src/core/iterative_matrix_functions.hpp
index e586927d..710fbb6b 100644
--- a/src/core/iterative_matrix_functions.hpp
+++ b/src/core/iterative_matrix_functions.hpp
@@ -108,7 +108,7 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
                                             FP imag_cutoff = (FP)1E-3) {
         int k = (int)vs.size();
         MKL_INT vm = vs[0].m, vn = vs[0].n, n = vm * vn;
-        if (davidson_type & DavidsonTypes::NonHermitian) {
+        if (davidson_type & DavidsonTypes::LeftEigen) {
             assert(k % 2 == 0);
             k = k / 2;
         }
@@ -137,7 +137,8 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
             threading->activate_global_mkl();
             if (davidson_type & DavidsonTypes::NonHermitian) {
                 tl.resize((size_t)n * n);
-                left.data = tl.data();
+                if (davidson_type & DavidsonTypes::LeftEigen)
+                    left.data = tl.data();
                 eig(h, ld, ld_imag, left);
             } else
                 eigs(h, ld);
@@ -156,7 +157,7 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
                 eigvals[i] = ld.data[eigval_idxs[i]];
                 copy(vs[i], GMatrix<FL>(h.data + eigval_idxs[i] * n, vm, vn));
             }
-            if (davidson_type & DavidsonTypes::NonHermitian)
+            if (davidson_type & DavidsonTypes::LeftEigen)
                 for (int i = 0; i < k; i++)
                     copy(vs[i + k],
                          GMatrix<FL>(left.data + eigval_idxs[i] * n, vm, vn));
@@ -165,7 +166,7 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
             pcomm->broadcast(eigvals.data(), eigvals.size(), pcomm->root);
             for (int j = 0; j < k; j++)
                 pcomm->broadcast(vs[j].data, vs[j].size(), pcomm->root);
-            if (davidson_type & DavidsonTypes::NonHermitian)
+            if (davidson_type & DavidsonTypes::LeftEigen)
                 for (int j = 0; j < k; j++)
                     pcomm->broadcast(vs[j + vs.size() / 2].data,
                                      vs[j + vs.size() / 2].size(), pcomm->root);
@@ -186,8 +187,10 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
         shared_ptr<VectorAllocator<FL>> d_alloc =
             make_shared<VectorAllocator<FL>>();
         int k = (int)vs.size();
-        assert(k % 2 == 0);
-        k = k / 2;
+        if (davidson_type & DavidsonTypes::LeftEigen) {
+            assert(k % 2 == 0);
+            k = k / 2;
+        }
         if (deflation_min_size < k)
             deflation_min_size = k;
         if (deflation_max_size < k + k / 2)
@@ -244,7 +247,9 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
                 sigmas[i].clear();
                 op(bs[i], sigmas[i]);
             }
-            GMatrix<FL> left(leftd.data(), m, m);
+            GMatrix<FL> left(nullptr, m, m);
+            if (davidson_type & DavidsonTypes::LeftEigen)
+                left.data = leftd.data();
             if (pcomm == nullptr || pcomm->root == pcomm->rank) {
                 GDiagonalMatrix<FP> ld(nullptr, m), ld_imag(nullptr, m);
                 GMatrix<FL> alpha(nullptr, m, m);
@@ -325,7 +330,8 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
                 }
                 for (int i = ck; i < k; i++) {
                     copy(qs[i - ck], sigmaxs[i]);
-                    iadd(qs[i - ck], bxs[i], -ld(eigval_idxs[i], eigval_idxs[i]));
+                    iadd(qs[i - ck], bxs[i],
+                         -ld(eigval_idxs[i], eigval_idxs[i]));
                 }
                 qq = complex_dot(qs[0], qs[0]);
                 if (iprint) {
@@ -339,8 +345,8 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
                 }
                 for (int i = ck; i < k; i++) {
                     if (davidson_type & DavidsonTypes::DavidsonPrecond)
-                        davidson_precondition(qs[i - ck], ld.data[eigval_idxs[ck]],
-                                              aa);
+                        davidson_precondition(qs[i - ck],
+                                              ld.data[eigval_idxs[ck]], aa);
                     else if (!(davidson_type & DavidsonTypes::NoPrecond))
                         olsen_precondition(qs[i - ck], bxs[i],
                                            ld.data[eigval_idxs[ck]], aa);
@@ -359,7 +365,8 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
             if (abs(qq) < conv_thrd) {
                 ck++;
                 if (ck == k) {
-                    if (pcomm == nullptr || pcomm->root == pcomm->rank) {
+                    if ((davidson_type & DavidsonTypes::LeftEigen) &&
+                        (pcomm == nullptr || pcomm->root == pcomm->rank)) {
                         int ntg = threading->activate_global();
 #pragma omp parallel num_threads(ntg)
                         {
@@ -406,7 +413,8 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
                     pcomm->broadcast(&m, 1, pcomm->root);
             }
             if (xiter == soft_max_iter) {
-                if (pcomm == nullptr || pcomm->root == pcomm->rank) {
+                if ((davidson_type & DavidsonTypes::LeftEigen) &&
+                    (pcomm == nullptr || pcomm->root == pcomm->rank)) {
                     int ntg = threading->activate_global();
 #pragma omp parallel num_threads(ntg)
                     {
@@ -438,16 +446,18 @@ template <typename FL> struct IterativeMatrixFunctions : GMatrixFunctions<FL> {
         if (pcomm == nullptr || pcomm->root == pcomm->rank) {
             for (int i = 0; i < k; i++)
                 copy(vs[i], bxs[i]);
-            for (int i = 0; i < k; i++)
-                copy(vs[i + vs.size() / 2], sigmaxs[i]);
+            if (davidson_type & DavidsonTypes::LeftEigen)
+                for (int i = 0; i < k; i++)
+                    copy(vs[i + vs.size() / 2], sigmaxs[i]);
         }
         if (pcomm != nullptr) {
             pcomm->broadcast(eigvals.data(), eigvals.size(), pcomm->root);
             for (int j = 0; j < k; j++)
                 pcomm->broadcast(vs[j].data, vs[j].size(), pcomm->root);
-            for (int j = 0; j < k; j++)
-                pcomm->broadcast(vs[j + vs.size() / 2].data,
-                                 vs[j + vs.size() / 2].size(), pcomm->root);
+            if (davidson_type & DavidsonTypes::LeftEigen)
+                for (int j = 0; j < k; j++)
+                    pcomm->broadcast(vs[j + vs.size() / 2].data,
+                                     vs[j + vs.size() / 2].size(), pcomm->root);
         }
         if (pcomm == nullptr || pcomm->root == pcomm->rank)
             q.deallocate();
diff --git a/src/core/matrix_functions.hpp b/src/core/matrix_functions.hpp
index add9869b..b7d3aabb 100644
--- a/src/core/matrix_functions.hpp
+++ b/src/core/matrix_functions.hpp
@@ -248,8 +248,12 @@ enum struct DavidsonTypes : uint16_t {
     NoPrecond = 64,
     NonHermitian = 128,
     Exact = 256,
+    LeftEigen = 512,
     ExactNonHermitian = 128 | 256,
-    NonHermitianDavidsonPrecond = 128 | 32
+    ExactNonHermitianLeftEigen = 128 | 256 | 512,
+    NonHermitianDavidsonPrecond = 128 | 32,
+    NonHermitianDavidsonPrecondLeftEigen = 128 | 32 | 512,
+    NonHermitianLeftEigen = 128 | 512
 };
 
 inline bool operator&(DavidsonTypes a, DavidsonTypes b) {
@@ -1352,8 +1356,9 @@ struct GMatrixFunctions<
         MKL_INT lwork = 34 * a.n, info;
         FL *work = d_alloc->allocate(lwork);
         FL *vr = d_alloc->allocate(a.m * a.n);
-        xgeev<FL>("V", lv.data == 0 ? "N" : "V", &a.n, a.data, &a.n, wr.data,
-                  vr, &a.n, lv.data, &a.n, work, &lwork, wi.data, &info);
+        xgeev<FL>("V", lv.data == nullptr ? "N" : "V", &a.n, a.data, &a.n,
+                  wr.data, vr, &a.n, lv.data, &a.n, work, &lwork, wi.data,
+                  &info);
         assert(info == 0);
         uint8_t tag = 0;
         copy(a, GMatrix<FL>(vr, a.m, a.n));
@@ -1362,8 +1367,9 @@ struct GMatrixFunctions<
                 k++;
                 for (MKL_INT j = 0; j < a.n; j++)
                     a(k, j) = -a(k, j);
-                for (MKL_INT j = 0; j < a.n; j++)
-                    lv(k, j) = -lv(k, j);
+                if (lv.data != nullptr)
+                    for (MKL_INT j = 0; j < a.n; j++)
+                        lv(k, j) = -lv(k, j);
             }
         d_alloc->deallocate(vr, a.m * a.n);
         d_alloc->deallocate(work, lwork);
diff --git a/src/dmrg/effective_hamiltonian.hpp b/src/dmrg/effective_hamiltonian.hpp
index e25f082d..afaab488 100644
--- a/src/dmrg/effective_hamiltonian.hpp
+++ b/src/dmrg/effective_hamiltonian.hpp
@@ -374,7 +374,7 @@ struct EffectiveHamiltonian<S, FL, MPS<S, FL>> {
         GDiagonalMatrix<FL> aa(diag->data, (MKL_INT)diag->total_memory);
         vector<GMatrix<FL>> bs = vector<GMatrix<FL>>{
             GMatrix<FL>(ket->data, (MKL_INT)ket->total_memory, 1)};
-        if (davidson_type & DavidsonTypes::NonHermitian)
+        if (davidson_type & DavidsonTypes::LeftEigen)
             bs.push_back(GMatrix<FL>(bra->data, (MKL_INT)bra->total_memory, 1));
         vector<GMatrix<FL>> ors =
             vector<GMatrix<FL>>(ortho_bra.size(), GMatrix<FL>(nullptr, 0, 0));
@@ -1357,7 +1357,7 @@ struct EffectiveHamiltonian<S, FL, MultiMPS<S, FL>> {
         for (int i = 0; i < (int)min((MKL_INT)ket.size(), (MKL_INT)aa.n); i++)
             bs.push_back(
                 GMatrix<FL>(ket[i]->data, (MKL_INT)ket[i]->total_memory, 1));
-        if (davidson_type & DavidsonTypes::NonHermitian)
+        if (davidson_type & DavidsonTypes::LeftEigen)
             for (int i = 0; i < (int)min((MKL_INT)bra.size(), (MKL_INT)aa.n);
                  i++)
                 bs.push_back(GMatrix<FL>(bra[i]->data,
diff --git a/src/dmrg/sweep_algorithm.hpp b/src/dmrg/sweep_algorithm.hpp
index d804ea72..03dc6549 100644
--- a/src/dmrg/sweep_algorithm.hpp
+++ b/src/dmrg/sweep_algorithm.hpp
@@ -2241,11 +2241,14 @@ template <typename S, typename FL, typename FLS> struct DMRG {
             throw runtime_error(
                 "Parallel MPS and different bra and ket is not yet supported!");
         if ((me->bra != me->ket &&
-             !(davidson_type & DavidsonTypes::NonHermitian)) ||
+             !((davidson_type & DavidsonTypes::NonHermitian) &&
+               (davidson_type & DavidsonTypes::LeftEigen))) ||
             (me->bra == me->ket &&
-             (davidson_type & DavidsonTypes::NonHermitian)))
-            throw runtime_error("Different BRA and KET must be used together "
-                                "with non-hermitian Hamiltonian!");
+             ((davidson_type & DavidsonTypes::NonHermitian) &&
+              (davidson_type & DavidsonTypes::LeftEigen))))
+            throw runtime_error(
+                "Different BRA and KET must be used together "
+                "with non-hermitian Hamiltonian and left eigen vector!");
         Timer start, current;
         start.get_time();
         current.get_time();
diff --git a/src/pybind/pybind_core.hpp b/src/pybind/pybind_core.hpp
index 1576e1b2..1400ac95 100644
--- a/src/pybind/pybind_core.hpp
+++ b/src/pybind/pybind_core.hpp
@@ -1523,6 +1523,12 @@ template <typename S = void> void bind_types(py::module &m) {
         .value("ExactNonHermitian", DavidsonTypes::ExactNonHermitian)
         .value("NonHermitianDavidsonPrecond",
                DavidsonTypes::NonHermitianDavidsonPrecond)
+        .value("LeftEigen", DavidsonTypes::LeftEigen)
+        .value("ExactNonHermitianLeftEigen",
+               DavidsonTypes::ExactNonHermitianLeftEigen)
+        .value("NonHermitianDavidsonPrecondLeftEigen",
+               DavidsonTypes::NonHermitianDavidsonPrecondLeftEigen)
+        .value("NonHermitianLeftEigen", DavidsonTypes::NonHermitianLeftEigen)
         .def(py::self & py::self)
         .def(py::self | py::self);
 
diff --git a/unit_test/debug_test_dmrg_n2_sto3g_nh.cpp b/unit_test/debug_test_dmrg_n2_sto3g_nh.cpp
index a8532fcb..548aed4f 100644
--- a/unit_test/debug_test_dmrg_n2_sto3g_nh.cpp
+++ b/unit_test/debug_test_dmrg_n2_sto3g_nh.cpp
@@ -98,9 +98,9 @@ void TestNonHermDMRGN2STO3G<FL>::test_dmrg(
             mps->save_data();
 
             // ME
-            shared_ptr<MPS<S, FL>> bra = mps->deep_copy("BRA");
+            // shared_ptr<MPS<S, FL>> bra = mps->deep_copy("BRA");
             shared_ptr<MovingEnvironment<S, FL, FL>> me =
-                make_shared<MovingEnvironment<S, FL, FL>>(mpo, bra, mps,
+                make_shared<MovingEnvironment<S, FL, FL>>(mpo, mps, mps,
                                                           "DMRG");
             me->init_environments(false);
             if (!condense)
@@ -110,7 +110,7 @@ void TestNonHermDMRGN2STO3G<FL>::test_dmrg(
             // DMRG
             shared_ptr<DMRG<S, FL, FL>> dmrg =
                 make_shared<DMRG<S, FL, FL>>(me, bdims, noises);
-            dmrg->davidson_type = DavidsonTypes::ExactNonHermitian;
+            dmrg->davidson_type = DavidsonTypes::NonHermitian;
             dmrg->iprint = 2;
             dmrg->decomp_type = dt;
             dmrg->noise_type = nt;
@@ -162,6 +162,7 @@ TYPED_TEST_CASE(TestNonHermDMRGN2STO3G, TestFL);
 
 TYPED_TEST(TestNonHermDMRGN2STO3G, TestSU2) {
     using FL = TypeParam;
+    using FLL = typename GMatrix<FL>::FL;
 
     shared_ptr<FCIDUMP<FL>> fcidump = make_shared<FCIDUMP<FL>>();
     PGTypes pg = PGTypes::D2H;
@@ -234,6 +235,7 @@ TYPED_TEST(TestNonHermDMRGN2STO3G, TestSU2) {
 
 TYPED_TEST(TestNonHermDMRGN2STO3G, TestSZ) {
     using FL = TypeParam;
+    using FLL = typename GMatrix<FL>::FL;
 
     shared_ptr<FCIDUMP<FL>> fcidump = make_shared<FCIDUMP<FL>>();
     PGTypes pg = PGTypes::D2H;
@@ -319,6 +321,7 @@ TYPED_TEST(TestNonHermDMRGN2STO3G, TestSZ) {
 
 TYPED_TEST(TestNonHermDMRGN2STO3G, TestSGF) {
     using FL = TypeParam;
+    using FLL = typename GMatrix<FL>::FL;
 
     shared_ptr<FCIDUMP<FL>> fcidump = make_shared<FCIDUMP<FL>>();
     PGTypes pg = PGTypes::D2H;
diff --git a/unit_test/debug_test_dmrg_sa_n2_sto3g_nh.cpp b/unit_test/debug_test_dmrg_sa_n2_sto3g_nh.cpp
index 90861d40..da31948b 100644
--- a/unit_test/debug_test_dmrg_sa_n2_sto3g_nh.cpp
+++ b/unit_test/debug_test_dmrg_sa_n2_sto3g_nh.cpp
@@ -51,8 +51,8 @@ void TestDMRGN2STO3GSA<FL>::test_dmrg(
     // MPO simplification
     cout << "MPO simplification start" << endl;
     mpo = make_shared<SimplifiedMPO<S, FL>>(
-        mpo, make_shared<NoTransposeRule<S, FL>>(make_shared<RuleQC<S, FL>>()), true, true,
-        OpNamesSet({OpNames::R, OpNames::RD}));
+        mpo, make_shared<NoTransposeRule<S, FL>>(make_shared<RuleQC<S, FL>>()),
+        true, true, OpNamesSet({OpNames::R, OpNames::RD}));
     cout << "MPO simplification end .. T = " << t.get_time() << endl;
 
     bond_dim = 120;
@@ -82,9 +82,9 @@ void TestDMRGN2STO3GSA<FL>::test_dmrg(
     mps->save_data();
 
     // ME
-    shared_ptr<MPS<S, FL>> bra = mps->deep_copy("BRA");
+    // shared_ptr<MPS<S, FL>> bra = mps->deep_copy("BRA");
     shared_ptr<MovingEnvironment<S, FL, FL>> me =
-        make_shared<MovingEnvironment<S, FL, FL>>(mpo, bra, mps, "DMRG");
+        make_shared<MovingEnvironment<S, FL, FL>>(mpo, mps, mps, "DMRG");
     me->init_environments(true);
     me->delayed_contraction = OpNamesSet::normal_ops();
     me->cached_contraction = true;
diff --git a/unit_test/test_matrix.cpp b/unit_test/test_matrix.cpp
index d782ae51..2ca251c6 100644
--- a/unit_test/test_matrix.cpp
+++ b/unit_test/test_matrix.cpp
@@ -1123,7 +1123,9 @@ TYPED_TEST(TestMatrix, TestDavidsonRealNonSymmetricExact) {
         }
         MatMul mop(a);
         vector<FL> vw = IterativeMatrixFunctions<FL>::davidson(
-            mop, aa, bs, 0, DavidsonTypes::NonHermitian | DavidsonTypes::Exact,
+            mop, aa, bs, 0,
+            DavidsonTypes::NonHermitian | DavidsonTypes::Exact |
+                DavidsonTypes::LeftEigen,
             ndav, true, (shared_ptr<ParallelCommunicator<SZ>>)nullptr, conv,
             n * k * 5, n * k * 4, k * 2, max((MKL_INT)5, k + 10));
         ASSERT_EQ((int)vw.size(), k);
@@ -1209,8 +1211,9 @@ TYPED_TEST(TestMatrix, TestDavidsonRealNonSymmetric) {
         MatMul mop(a);
         vector<FL> vw = IterativeMatrixFunctions<FL>::davidson(
             mop, aa, bs, 0,
-            DavidsonTypes::NonHermitian | DavidsonTypes::DavidsonPrecond, ndav,
-            false, (shared_ptr<ParallelCommunicator<SZ>>)nullptr, conv,
+            DavidsonTypes::NonHermitian | DavidsonTypes::DavidsonPrecond |
+                DavidsonTypes::LeftEigen,
+            ndav, false, (shared_ptr<ParallelCommunicator<SZ>>)nullptr, conv,
             n * k * 50, n * k * 40, k * 2, min(max((MKL_INT)5, k * 5 + 10), n));
         ASSERT_EQ((int)vw.size(), k);
         GDiagonalMatrix<FL> w(&vw[0], k);