diff --git a/src/dmrg/mpo.hpp b/src/dmrg/mpo.hpp index 4ab3c53c..6edb0355 100644 --- a/src/dmrg/mpo.hpp +++ b/src/dmrg/mpo.hpp @@ -301,6 +301,15 @@ template struct MPO { } return make_tuple(max(lnnz, rnnz), max(lsz, rsz), bdim); } + vector get_bond_dims() { + vector r; + for (int ii = 0; ii < n_sites; ii++) { + this->load_left_operators(ii); + r.push_back(left_operator_names[ii]->n); + this->unload_left_operators(ii); + } + return r; + } string get_filename(int i, int ixtag, const string &dir = "") const { const static string xtag[] = {"TENSOR", "LEFT.OP", "RIGHT.OP", "MIDDLE.OP"}; @@ -1621,8 +1630,8 @@ template struct TransMPO { static shared_ptr> forward(const shared_ptr> &mpo, const shared_ptr> &hamil, - const string &tag = "") { - shared_ptr> rmpo = make_shared>( + const string &tag = "", S2 left_vacuum = S2(S2::invalid)) { + shared_ptr> rmpo = make_shared>( mpo->n_sites, tag == "" ? "TR-" + mpo->tag : tag); const int n_sites = mpo->n_sites; rmpo->const_e = mpo->const_e; @@ -1639,145 +1648,416 @@ template struct TransMPO { rmpo->sparse_form = mpo->sparse_form; rmpo->schemer = nullptr; assert(mpo->schemer == nullptr); + shared_ptr> cg = mpo->hamil->opf->cg; shared_ptr> zero = make_shared>(); - auto tr = [&ref](S1 q) -> S2 { + auto trx = [&ref](S1 q) -> shared_ptr> { return TransStateInfo::forward( - make_shared>(q), ref) - ->quanta[0]; + make_shared>(q), ref); }; - rmpo->op = make_shared>( - mpo->op->name, mpo->op->site_index, tr(mpo->op->q_label), - mpo->op->factor); - rmpo->left_vacuum = tr(mpo->left_vacuum); + auto tr_opq = trx(mpo->op->q_label), tr_vac = trx(mpo->left_vacuum); + assert(tr_opq->n == 1); + rmpo->op = + make_shared>(mpo->op->name, mpo->op->site_index, + tr_opq->quanta[0], mpo->op->factor); + if (left_vacuum == S2(S2::invalid)) { + assert(tr_vac->n == 1); + rmpo->left_vacuum = tr_vac->quanta[0]; + } else { + assert(tr_vac->find_state(left_vacuum) != -1); + rmpo->left_vacuum = left_vacuum; + } + vector>> llopd; + llopd.push_back(make_shared>(OpNames::I, SiteIndex(), + rmpo->left_vacuum)); for (int ii = 0; ii < n_sites; ii++) { mpo->load_tensor(ii); mpo->load_left_operators(ii); mpo->load_right_operators(ii); + if (ii != 0) + mpo->load_left_operators(ii - 1); + vector l_multi(mpo->left_operator_names[ii]->n); + vector>> lopd; + lopd.reserve(mpo->left_operator_names[ii]->n); + for (int iop = 0, jop = 0; iop < mpo->left_operator_names[ii]->n; + iop++) { + const auto p = dynamic_pointer_cast>( + mpo->left_operator_names[ii]->data[iop]); + auto xqs = trx(p->q_label); + for (int kop = 0; kop < xqs->n; kop++) { + FL factor = + ii == 0 ? (FL)cg->cg(mpo->left_vacuum, p->q_label, + p->q_label, rmpo->left_vacuum, + xqs->quanta[kop], xqs->quanta[kop]) + : (FL)1.0; + lopd.push_back(make_shared>( + OpNames::XL, + SiteIndex({(uint16_t)(iop / 1000 / 1000), + (uint16_t)(iop / 1000 % 1000), + (uint16_t)(iop % 1000), (uint16_t)kop}, + {}), + xqs->quanta[kop], p->factor * factor)); + } + jop += (l_multi[iop] = xqs->n); + } + vector r_multi(mpo->right_operator_names[ii]->m); + vector>> ropd; + ropd.reserve(mpo->right_operator_names[ii]->m); + for (int iop = 0, jop = 0; iop < mpo->right_operator_names[ii]->m; + iop++) { + const auto p = dynamic_pointer_cast>( + mpo->right_operator_names[ii]->data[iop]); + auto xqs = trx(p->q_label); + for (int kop = xqs->n - 1; kop >= 0; kop--) { + FL factor = + ii == n_sites - 1 + ? (FL)cg->cg( + dynamic_pointer_cast>( + mpo->left_operator_names[ii - 1] + ->data[iop]) + ->q_label, + p->q_label, mpo->op->q_label, + dynamic_pointer_cast>( + llopd[jop + xqs->n - 1 - kop]) + ->q_label, + xqs->quanta[kop], tr_opq->quanta[0]) + : (FL)1.0; + ropd.push_back(make_shared>( + OpNames::XR, + SiteIndex({(uint16_t)(iop / 1000 / 1000), + (uint16_t)(iop / 1000 % 1000), + (uint16_t)(iop % 1000), (uint16_t)kop}, + {}), + xqs->quanta[kop], p->factor * factor)); + } + jop += (r_multi[iop] = xqs->n); + } + auto lop = make_shared>((int)lopd.size()); + auto rop = make_shared>((int)ropd.size()); + lop->data = lopd; + rop->data = ropd; + rmpo->left_operator_names[ii] = lop; + rmpo->right_operator_names[ii] = rop; + vector l_accu(l_multi.size(), 0); + vector r_accu(r_multi.size(), 0); + for (size_t iop = 0; iop < l_multi.size() - 1; iop++) + l_accu[iop + 1] = l_accu[iop] + l_multi[iop]; + for (size_t iop = 0; iop < r_multi.size() - 1; iop++) + r_accu[iop + 1] = r_accu[iop] + r_multi[iop]; shared_ptr> opt = rmpo->tensors[ii]; shared_ptr> mopt = mpo->tensors[ii]; - shared_ptr> pmat; assert(mopt->lmat == mopt->rmat); + shared_ptr> pmat; if (ii == 0) - pmat = make_shared>(mopt->lmat->n); + pmat = make_shared>(lop->n); else if (ii == n_sites - 1) - pmat = make_shared>(mopt->lmat->m); - else { - pmat = make_shared>(mopt->lmat->m, - mopt->lmat->n); - dynamic_pointer_cast>(pmat)->indices = - dynamic_pointer_cast>(mopt->lmat) - ->indices; - pmat->data.resize(mopt->lmat->data.size()); - } + pmat = make_shared>(rop->m); + else + pmat = make_shared>(rop->m, lop->n); + unordered_map>, + vector>>> + op_map; opt->lmat = opt->rmat = pmat; - for (int iop = 0; iop < mopt->lmat->data.size(); iop++) - if (mopt->lmat->data[iop]->get_type() == OpTypes::Zero) - pmat->data[iop] = zero; - else if (mopt->lmat->data[iop]->get_type() == OpTypes::Elem) { - const auto p = dynamic_pointer_cast>( - mopt->lmat->data[iop]); - pmat->data[iop] = make_shared>( - p->name, p->site_index, tr(p->q_label), p->factor); - } else if (mopt->lmat->data[iop]->get_type() == OpTypes::Sum) { - const auto p = dynamic_pointer_cast>( - mopt->lmat->data[iop]); - vector>> strings(p->strings.size()); - for (size_t j = 0; j < p->strings.size(); j++) { - const auto pp = dynamic_pointer_cast>( - p->strings[j]); - strings[j] = make_shared>( - pp->name, pp->site_index, tr(pp->q_label), - pp->factor); + size_t imxx = 0; + if (ii == 0 || ii == n_sites - 1) { + const vector &xmulti = ii == 0 ? l_multi : r_multi; + unordered_map>, int> xmap; + xmap.reserve(pmat->data.size()); + for (int iop = 0, jop = 0; iop < mopt->lmat->data.size(); + jop += xmulti[iop++]) + if (mopt->lmat->data[iop]->get_type() == OpTypes::Zero) { + for (int kop = 0; kop < xmulti[iop]; kop++) + pmat->data[jop + kop] = zero; + } else { + vector>> pps; + if (mopt->lmat->data[iop]->get_type() == OpTypes::Elem) + pps.push_back( + dynamic_pointer_cast>( + mopt->lmat->data[iop])); + else if (mopt->lmat->data[iop]->get_type() == + OpTypes::Sum) { + const auto p = dynamic_pointer_cast>( + mopt->lmat->data[iop]); + for (size_t j = 0; j < p->strings.size(); j++) { + const auto pp = + dynamic_pointer_cast>( + p->strings[j]); + } + } else + assert(false); + const S1 gxll = + ii == 0 ? S1(S1::invalid) + : dynamic_pointer_cast>( + mpo->left_operator_names[ii - 1] + ->data[iop]) + ->q_label; + for (int kop = 0; kop < xmulti[iop]; kop++) { + const size_t imx = jop + kop; + const S1 gx = + dynamic_pointer_cast>( + ii == 0 ? mpo->left_operator_names[ii] + ->data[iop] + : mpo->right_operator_names[ii] + ->data[iop]) + ->q_label; + const S2 gq = + dynamic_pointer_cast>( + ii == 0 ? lopd[jop + kop] : ropd[jop + kop]) + ->q_label; + const S2 gqll = + ii == 0 + ? S2(S2::invalid) + : dynamic_pointer_cast>( + llopd[jop + kop]) + ->q_label; + vector>> exprs; + for (const auto &p : pps) { + auto ap = abs_value((shared_ptr>)p); + auto xqs = trx(p->q_label); + if (!op_map.count(ap)) { + for (int mop = 0; mop < xqs->n; mop++) + op_map[ap].push_back(make_shared< + OpElement>( + ii == 0 ? OpNames::XL : OpNames::XR, + SiteIndex( + {(uint16_t)(iop / 1000 / 1000), + (uint16_t)(iop / 1000 % 1000), + (uint16_t)(iop % 1000), + (uint16_t)mop}, + {}), + xqs->quanta[mop])); + } + for (int mop = 0; mop < xqs->n; mop++) { + if (ii == 0 && xqs->quanta[mop].combine( + gq, rmpo->left_vacuum) == + S2(S2::invalid)) + continue; + if (ii == n_sites - 1 && + xqs->quanta[mop] != gq) + continue; + assert(ii == n_sites - 1 || mop == kop); + FL factor = + ii == 0 + ? (FL)cg->cg(mpo->left_vacuum, + p->q_label, gx, + rmpo->left_vacuum, + xqs->quanta[mop], gq) + : (FL)cg->cg(gxll, p->q_label, + mpo->op->q_label, gqll, + xqs->quanta[mop], + tr_opq->quanta[0]); + if (abs(factor) < TINY) + continue; + exprs.push_back( + op_map.at(ap).at(mop)->scalar_multiply( + factor * p->factor)); + } + } + if (exprs.size() == 0) + pmat->data[jop + kop] = zero; + else if (exprs.size() == 1) + pmat->data[jop + kop] = exprs[0]; + else + pmat->data[jop + kop] = sum(exprs); + } } - pmat->data[iop] = sum(strings); - } else - assert(false); - auto lop = make_shared>( - mpo->left_operator_names[ii]->n); - auto rop = make_shared>( - mpo->right_operator_names[ii]->m); - for (int iop = 0; iop < lop->data.size(); iop++) { - const auto p = dynamic_pointer_cast>( - mpo->left_operator_names[ii]->data[iop]); - lop->data[iop] = make_shared>( - p->name, p->site_index, tr(p->q_label), p->factor); + } else { + shared_ptr> ppmat = + dynamic_pointer_cast>(pmat); + shared_ptr> mpmat = + dynamic_pointer_cast>(mopt->lmat); + ppmat->indices.reserve(mpmat->indices.size()); + ppmat->data.reserve(mpmat->data.size()); + for (int iiop = 0; iiop < mpmat->data.size(); iiop++) { + const int irop = mpmat->indices[iiop].first; + const int ilop = mpmat->indices[iiop].second; + if (mpmat->data[iiop]->get_type() == OpTypes::Zero) + continue; + vector>> pps; + if (mpmat->data[iiop]->get_type() == OpTypes::Elem) + pps.push_back(dynamic_pointer_cast>( + mpmat->data[iiop])); + else if (mpmat->data[iiop]->get_type() == OpTypes::Sum) { + const auto p = dynamic_pointer_cast>( + mpmat->data[iiop]); + for (size_t j = 0; j < p->strings.size(); j++) + pps.push_back( + dynamic_pointer_cast>( + p->strings[j]->get_op()->scalar_multiply( + p->strings[j]->factor))); + } else + assert(false); + for (int krop = 0; krop < r_multi[irop]; krop++) { + const size_t irmx = r_accu[irop] + krop; + const S1 gxll = + dynamic_pointer_cast>( + mpo->left_operator_names[ii - 1]->data[irop]) + ->q_label; + const S2 gqll = + dynamic_pointer_cast>(llopd[irmx]) + ->q_label; + for (int klop = 0; klop < l_multi[ilop]; klop++) { + const size_t ilmx = l_accu[ilop] + klop; + const S1 gxl = + dynamic_pointer_cast>( + mpo->left_operator_names[ii]->data[ilop]) + ->q_label; + const S2 gql = + dynamic_pointer_cast>( + lopd[ilmx]) + ->q_label; + vector>> exprs; + for (const auto &p : pps) { + auto ap = abs_value((shared_ptr>)p); + auto xqs = trx(p->q_label); + if (!op_map.count(ap)) { + if (dynamic_pointer_cast>( + ap) + ->name == OpNames::I) { + assert(xqs->n == 1); + op_map[ap].push_back( + make_shared>( + p->name, p->site_index, + xqs->quanta[0])); + } else { + for (int mop = 0; mop < xqs->n; mop++) + op_map[ap].push_back( + make_shared>( + OpNames::X, + SiteIndex( + {(uint16_t)(imxx / 4000 / + 4000), + (uint16_t)(imxx / 4000 % + 4000), + (uint16_t)(imxx % 4000), + (uint16_t)mop}, + {}), + xqs->quanta[mop])); + imxx++; + } + } + for (int mop = 0; mop < xqs->n; mop++) { + if (xqs->quanta[mop].combine(gql, gqll) == + S2(S2::invalid)) + continue; + FL factor = + (FL)cg->cg(gxll, p->q_label, gxl, gqll, + xqs->quanta[mop], gql); + if (abs(factor) < TINY) + continue; + exprs.push_back( + op_map.at(ap).at(mop)->scalar_multiply( + factor * p->factor)); + } + } + if (exprs.size() != 0) + ppmat->indices.push_back( + make_pair((int)irmx, (int)ilmx)); + if (exprs.size() == 1) + ppmat->data.push_back(exprs[0]); + else if (exprs.size() != 0) + ppmat->data.push_back(sum(exprs)); + } + } + } } - for (int iop = 0; iop < rop->data.size(); iop++) { - const auto p = dynamic_pointer_cast>( - mpo->right_operator_names[ii]->data[iop]); - rop->data[iop] = make_shared>( - p->name, p->site_index, tr(p->q_label), p->factor); + for (const auto &op : mopt->ops) { + if (!op_map.count(op.first)) { + const auto p = + dynamic_pointer_cast>(op.first); + auto xqs = trx(op.second->info->delta_quantum); + if (xqs->n == 1 && p->name == OpNames::I) + for (int mop = 0; mop < xqs->n; mop++) + op_map[op.first].push_back( + make_shared>( + p->name, p->site_index, xqs->quanta[0], + p->factor)); + } } - rmpo->left_operator_names[ii] = lop; - rmpo->right_operator_names[ii] = rop; shared_ptr> d_alloc = make_shared>(); shared_ptr> conn = TransStateInfo::backward_connection(mpo->basis[ii], rmpo->basis[ii]); - for (const auto &op : mopt->ops) { + for (const auto &opx : op_map) { const auto p = - dynamic_pointer_cast>(op.first); - auto q = make_shared>( - p->name, p->site_index, tr(p->q_label), p->factor); - shared_ptr> xmat = - make_shared>(d_alloc); - xmat->allocate(hamil->find_site_op_info(ii, q->q_label)); - xmat->factor = op.second->factor; - for (int k = 0; k < op.second->info->n; k++) { - S1 plu = op.second->info->quanta[k].get_bra( - op.second->info->delta_quantum); - S1 pru = op.second->info->quanta[k].get_ket(); - GMatrix r = (*op.second)[k]; - shared_ptr> mls = - TransStateInfo::forward( - make_shared>(plu), ref); - shared_ptr> mrs = - TransStateInfo::forward( - make_shared>(pru), ref); - for (int iln = 0; iln < mls->n; iln++) - for (int irn = 0; irn < mrs->n; irn++) { - S2 lqn = mls->quanta[iln], rqn = mrs->quanta[irn]; - GMatrix xr = - (*xmat)[q->q_label.combine(lqn, rqn)]; - int il = rmpo->basis[ii]->find_state(lqn); - int ir = rmpo->basis[ii]->find_state(rqn); - MKL_INT zl = rmpo->basis[ii]->n_states[il], - zr = rmpo->basis[ii]->n_states[ir]; - int klst = conn->n_states[il]; - int krst = conn->n_states[ir]; - int kled = il == rmpo->basis[ii]->n - 1 - ? conn->n - : conn->n_states[il + 1]; - int kred = ir == rmpo->basis[ii]->n - 1 - ? conn->n - : conn->n_states[ir + 1]; - size_t lsh = 0, rsh = 0; - for (int ilp = klst; - ilp < kled && conn->quanta[ilp] != plu; ilp++) - lsh += - mpo->basis[ii] + dynamic_pointer_cast>(opx.first); + const auto pmat = mopt->ops.at(opx.first); + for (size_t mop = 0; mop < opx.second.size(); mop++) { + auto q = opx.second[mop]; + shared_ptr> xmat = + make_shared>(d_alloc); + xmat->allocate(hamil->find_site_op_info(ii, q->q_label)); + xmat->factor = pmat->factor; + for (int k = 0; k < pmat->info->n; k++) { + S1 plu = pmat->info->quanta[k].get_bra( + pmat->info->delta_quantum); + S1 pru = pmat->info->quanta[k].get_ket(); + GMatrix r = (*pmat)[k]; + shared_ptr> mls = + TransStateInfo::forward( + make_shared>(plu), ref); + shared_ptr> mrs = + TransStateInfo::forward( + make_shared>(pru), ref); + for (int iln = 0; iln < mls->n; iln++) + for (int irn = 0; irn < mrs->n; irn++) { + S2 lqn = mls->quanta[iln], + rqn = mrs->quanta[irn]; + if (q->q_label.combine(lqn, rqn) == + S2(S2::invalid)) + continue; + FL factor = (FL)cg->cg(pru, p->q_label, plu, + rqn, q->q_label, lqn); + if (abs(factor) < TINY) + continue; + GMatrix xr = + (*xmat)[q->q_label.combine(lqn, rqn)]; + int il = rmpo->basis[ii]->find_state(lqn); + int ir = rmpo->basis[ii]->find_state(rqn); + MKL_INT zl = rmpo->basis[ii]->n_states[il], + zr = rmpo->basis[ii]->n_states[ir]; + int klst = conn->n_states[il]; + int krst = conn->n_states[ir]; + int kled = il == rmpo->basis[ii]->n - 1 + ? conn->n + : conn->n_states[il + 1]; + int kred = ir == rmpo->basis[ii]->n - 1 + ? conn->n + : conn->n_states[ir + 1]; + size_t lsh = 0, rsh = 0; + for (int ilp = klst; + ilp < kled && conn->quanta[ilp] != plu; + ilp++) + lsh += mpo->basis[ii]->n_states + [mpo->basis[ii]->find_state( + conn->quanta[ilp])]; + for (int irp = krst; + irp < kred && conn->quanta[irp] != pru; + irp++) + rsh += mpo->basis[ii]->n_states + [mpo->basis[ii]->find_state( + conn->quanta[irp])]; + MKL_INT kl = + (MKL_INT)mpo->basis[ii] ->n_states[mpo->basis[ii]->find_state( - conn->quanta[ilp])]; - for (int irp = krst; - irp < kred && conn->quanta[irp] != pru; irp++) - rsh += - mpo->basis[ii] + plu)]; + MKL_INT kr = + (MKL_INT)mpo->basis[ii] ->n_states[mpo->basis[ii]->find_state( - conn->quanta[irp])]; - MKL_INT kl = - (MKL_INT)mpo->basis[ii] - ->n_states[mpo->basis[ii]->find_state(plu)]; - MKL_INT kr = - (MKL_INT)mpo->basis[ii] - ->n_states[mpo->basis[ii]->find_state(pru)]; - for (MKL_INT ikl = 0; ikl < kl; ikl++) - for (MKL_INT ikr = 0; ikr < kr; ikr++) - xr(ikl + lsh, ikr + rsh) = r(ikl, ikr); - } + pru)]; + for (MKL_INT ikl = 0; ikl < kl; ikl++) + for (MKL_INT ikr = 0; ikr < kr; ikr++) + xr(ikl + lsh, ikr + rsh) = + (FL)factor * r(ikl, ikr); + } + } + assert(q->q_label == xmat->info->delta_quantum); + opt->ops[q] = xmat; } - opt->ops[q] = xmat; } + llopd = lopd; + if (ii != 0) + mpo->unload_left_operators(ii - 1); mpo->unload_tensor(ii); mpo->unload_left_operators(ii); mpo->unload_right_operators(ii); diff --git a/src/dmrg/mpo_fusing.hpp b/src/dmrg/mpo_fusing.hpp index 9e5bdee1..0a9cf9e8 100644 --- a/src/dmrg/mpo_fusing.hpp +++ b/src/dmrg/mpo_fusing.hpp @@ -147,6 +147,10 @@ template struct StackedMPO : MPO { for (int imb = 0; imb < mpob->left_operator_names[m]->n; imb++) { const int imx = ima * mpob->left_operator_names[m]->n + imb; + const auto opxa = dynamic_pointer_cast>( + mpoa->left_operator_names[m]->data[ima]); + const auto opxb = dynamic_pointer_cast>( + mpob->left_operator_names[m]->data[imb]); left_operator_names[m]->data[imx] = make_shared>( OpNames::XL, @@ -154,18 +158,18 @@ template struct StackedMPO : MPO { (uint16_t)(imx / 1000 % 1000), (uint16_t)(imx % 1000)}, {}), - dynamic_pointer_cast>( - mpoa->left_operator_names[m]->data[ima]) - ->q_label + - dynamic_pointer_cast>( - mpob->left_operator_names[m]->data[imb]) - ->q_label); + opxa->q_label + opxb->q_label, + opxa->factor * opxb->factor); } for (int ima = 0; ima < mpoa->right_operator_names[m]->m; ima++) for (int imb = 0; imb < mpob->right_operator_names[m]->m; imb++) { const int imx = ima * mpob->right_operator_names[m]->m + imb; + const auto opxa = dynamic_pointer_cast>( + mpoa->right_operator_names[m]->data[ima]); + const auto opxb = dynamic_pointer_cast>( + mpob->right_operator_names[m]->data[imb]); right_operator_names[m]->data[imx] = make_shared>( OpNames::XR, @@ -173,12 +177,8 @@ template struct StackedMPO : MPO { (uint16_t)(imx / 1000 % 1000), (uint16_t)(imx % 1000)}, {}), - dynamic_pointer_cast>( - mpoa->right_operator_names[m]->data[ima]) - ->q_label + - dynamic_pointer_cast>( - mpob->right_operator_names[m]->data[imb]) - ->q_label); + opxa->q_label + opxb->q_label, + opxa->factor * opxb->factor); } if (m == 0) { shared_ptr> pmata = @@ -341,6 +341,40 @@ template struct StackedMPO : MPO { pmat->data[imx] = zero; continue; } + vector>> ppas; + vector>> ppbs; + if (pmata->data[ima]->get_type() == OpTypes::Elem) + ppas.push_back( + dynamic_pointer_cast>( + pmata->data[ima])); + else if (pmata->data[ima]->get_type() == OpTypes::Sum) { + const auto p = dynamic_pointer_cast>( + pmata->data[ima]); + for (size_t j = 0; j < p->strings.size(); j++) + ppas.push_back( + dynamic_pointer_cast>( + p->strings[j] + ->get_op() + ->scalar_multiply( + p->strings[j]->factor))); + } else + assert(false); + if (pmatb->data[imb]->get_type() == OpTypes::Elem) + ppbs.push_back( + dynamic_pointer_cast>( + pmatb->data[imb])); + else if (pmatb->data[imb]->get_type() == OpTypes::Sum) { + const auto p = dynamic_pointer_cast>( + pmatb->data[imb]); + for (size_t j = 0; j < p->strings.size(); j++) + ppbs.push_back( + dynamic_pointer_cast>( + p->strings[j] + ->get_op() + ->scalar_multiply( + p->strings[j]->factor))); + } else + assert(false); shared_ptr> opel = make_shared>( OpNames::X, @@ -349,40 +383,34 @@ template struct StackedMPO : MPO { (uint16_t)(imx / 1000 % 1000), (uint16_t)(imx % 1000)}, {}), - dynamic_pointer_cast>( - pmata->data[ima]) - ->q_label + - dynamic_pointer_cast>( - pmatb->data[imb]) - ->q_label); - shared_ptr> mata = - mpoa->tensors[m]->ops.at( - abs_value(pmata->data[ima])); - shared_ptr> matb = - mpob->tensors[m]->ops.at( - abs_value(pmatb->data[imb])); + ppas[0]->q_label + ppbs[0]->q_label); shared_ptr> xmat = make_shared>(xd_alloc); - const S q = (mata->info->delta_quantum + - matb->info->delta_quantum)[0]; - const FL phase = - mata->info->delta_quantum.is_fermion() && - dynamic_pointer_cast>( - mpob->right_operator_names[m] - ->data[pmatb->indices[imb].first]) - ->q_label.is_fermion() - ? -1 - : 1; + const S q = (ppas[0]->q_label + ppbs[0]->q_label)[0]; xmat->allocate(site_op_infos_mp[m].at(q)); - MPO::tf->opf->product(0, mata, matb, xmat, - phase); - pmat->data[imx] = opel->scalar_multiply( - dynamic_pointer_cast>( - pmata->data[ima]) - ->factor * - dynamic_pointer_cast>( - pmatb->data[imb]) - ->factor); + for (const auto xpa : ppas) + for (const auto xpb : ppbs) { + shared_ptr> mata = + mpoa->tensors[m]->ops.at( + abs_value((shared_ptr>)xpa)); + shared_ptr> matb = + mpob->tensors[m]->ops.at( + abs_value((shared_ptr>)xpb)); + const FL phase = + mata->info->delta_quantum.is_fermion() && + dynamic_pointer_cast< + OpElement>( + mpob->right_operator_names[m] + ->data[pmatb->indices[imb] + .first]) + ->q_label.is_fermion() + ? -1 + : 1; + MPO::tf->opf->product( + 0, mata, matb, xmat, + phase * xpa->factor * xpb->factor); + } + pmat->data[imx] = opel; mats[imx] = xmat; } } diff --git a/src/pybind/pybind_core.hpp b/src/pybind/pybind_core.hpp index 101bf636..02d295c8 100644 --- a/src/pybind/pybind_core.hpp +++ b/src/pybind/pybind_core.hpp @@ -618,6 +618,7 @@ template void bind_fl_expr(py::module &m) { .def_readwrite("conj", &OpProduct::conj) .def_readwrite("a", &OpProduct::a) .def_readwrite("b", &OpProduct::b) + .def("get_op", &OpProduct::get_op) .def("__hash__", &OpProduct::hash); py::class_, shared_ptr>, diff --git a/src/pybind/pybind_dmrg.hpp b/src/pybind/pybind_dmrg.hpp index 1c4e0d23..b0b69ccb 100644 --- a/src/pybind/pybind_dmrg.hpp +++ b/src/pybind/pybind_dmrg.hpp @@ -1814,6 +1814,7 @@ template void bind_fl_mpo(py::module &m) { .def_readwrite("tread", &MPO::tread) .def_readwrite("twrite", &MPO::twrite) .def("get_summary", &MPO::get_summary) + .def("get_bond_dims", &MPO::get_bond_dims) .def("get_filename", &MPO::get_filename) .def("load_left_operators", &MPO::load_left_operators) .def("save_left_operators", &MPO::save_left_operators) @@ -2269,7 +2270,8 @@ template void bind_fl_trans_mpo(py::module &m, const string &aux_name) { m.def(("trans_mpo_to_" + aux_name).c_str(), &TransMPO::forward, - py::arg("mpo"), py::arg("hamil"), py::arg("tag") = ""); + py::arg("mpo"), py::arg("hamil"), py::arg("tag") = "", + py::arg("left_vacuum") = T(T::invalid)); } template void bind_dmrg_types(py::module &m) {