non-hermi without left eigen

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

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();

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);

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,

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();

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);
 

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;

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;

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);