mpo trans symm

pyblock2/driver/core.py

@@ -6798,6 +6798,36 @@ def mps_flip_twos(self, mps):
         zmps = umps.finalize(self.prule)
         return zmps
 
+    def mpo_change_symm(self, mpo, tag="", add_ident=True):
+        """
+        Change symmetry type of MPO.
+        Only works in SAny mode. The resulting MPO should be used in SAny mode.
+
+        Args:
+            mpo : MPO
+                The input MPO.
+            tag : str
+                The tag of the output MPO.
+            add_ident : bool
+                If True, the hidden identity operator will be added into the MPO. Default is True.
+
+        Returns:
+            rmpo : MPO
+                The output MPO.
+        """
+        bw = self.bw
+        assert SymmetryTypes.SAny in bw.symm_type
+        if self.mpi:
+            rmpo = bw.bs.trans_mpo_to_sany(mpo.prim_mpo.prim_mpo, self.ghamil, tag)
+        else:
+            rmpo = bw.bs.trans_mpo_to_sany(mpo.prim_mpo, self.ghamil, tag)
+        rmpo = bw.bs.SimplifiedMPO(rmpo, bw.bs.Rule(), False, False)
+        if add_ident:
+            rmpo = bw.bs.IdentityAddedMPO(rmpo)
+        if self.mpi:
+            rmpo = bw.bs.ParallelMPO(rmpo, self.prule)
+        return rmpo
+
     def mps_change_symm(self, mps, tag, target):
         """
         Change symmetry type of MPS.

src/dmrg/mpo.hpp

@@ -1616,4 +1616,180 @@ template <typename S, typename FL> struct IdentityAddedMPO : MPO<S, FL> {
     }
 };
 
+template <typename S1, typename S2, typename FL> struct TransMPO {
+    typedef typename GMatrix<FL>::FP FP;
+    static shared_ptr<MPO<S2, FL>>
+    forward(const shared_ptr<MPO<S1, FL>> &mpo,
+            const shared_ptr<Hamiltonian<S2, FL>> &hamil,
+            const string &tag = "") {
+        shared_ptr<MPO<S1, FL>> rmpo = make_shared<MPO<S1, FL>>(
+            mpo->n_sites, tag == "" ? "TR-" + mpo->tag : tag);
+        const int n_sites = mpo->n_sites;
+        rmpo->const_e = mpo->const_e;
+        const S2 ref = hamil->vacuum;
+        rmpo->tf = make_shared<TensorFunctions<S2, FL>>(hamil->opf);
+        rmpo->site_op_infos = hamil->site_op_infos;
+        rmpo->basis = hamil->basis;
+        rmpo->hamil = hamil;
+        rmpo->left_operator_names.resize(n_sites, nullptr);
+        rmpo->right_operator_names.resize(n_sites, nullptr);
+        rmpo->tensors.resize(n_sites, nullptr);
+        for (uint16_t m = 0; m < n_sites; m++)
+            rmpo->tensors[m] = make_shared<OperatorTensor<S2, FL>>();
+        rmpo->sparse_form = mpo->sparse_form;
+        rmpo->schemer = nullptr;
+        assert(mpo->schemer == nullptr);
+        shared_ptr<OpExpr<S2>> zero = make_shared<OpExpr<S2>>();
+        auto tr = [&ref](S1 q) -> S2 {
+            return TransStateInfo<S1, S2>::forward(
+                       make_shared<StateInfo<S1>>(q), ref)
+                ->quanta[0];
+        };
+        rmpo->op = make_shared<OpElement<S2, FL>>(
+            mpo->op->name, mpo->op->site_index, tr(mpo->op->q_label),
+            mpo->op->factor);
+        rmpo->left_vacuum = tr(mpo->left_vacuum);
+        for (int ii = 0; ii < n_sites; ii++) {
+            mpo->load_tensor(ii);
+            mpo->load_left_operators(ii);
+            mpo->load_right_operators(ii);
+            shared_ptr<OperatorTensor<S2, FL>> opt = rmpo->tensors[ii];
+            shared_ptr<OperatorTensor<S2, FL>> mopt = mpo->tensors[ii];
+            shared_ptr<Symbolic<S2>> pmat;
+            assert(mopt->lmat == mopt->rmat);
+            if (ii == 0)
+                pmat = make_shared<SymbolicRowVector<S2>>(mopt->lmat->n);
+            else if (ii == n_sites - 1)
+                pmat = make_shared<SymbolicColumnVector<S2>>(mopt->lmat->m);
+            else {
+                pmat = make_shared<SymbolicMatrix<S2>>(mopt->lmat->m,
+                                                       mopt->lmat->n);
+                dynamic_pointer_cast<SymbolicMatrix<S2>>(pmat)->indices =
+                    dynamic_pointer_cast<SymbolicMatrix<S1>>(mopt->lmat)
+                        ->indices;
+                pmat->data.resize(mopt->lmat->data.size());
+            }
+            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<OpElement<S1, FL>>(
+                        mopt->lmat->data[iop]);
+                    pmat->data[iop] = make_shared<OpElement<S2, FL>>(
+                        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<OpSum<S1, FL>>(
+                        mopt->lmat->data[iop]);
+                    vector<shared_ptr<OpExpr<S2>>> strings(p->strings.size());
+                    for (size_t j = 0; j < p->strings.size(); j++) {
+                        const auto pp = dynamic_pointer_cast<OpElement<S1, FL>>(
+                            p->strings[j]);
+                        strings[j] = make_shared<OpElement<S2, FL>>(
+                            pp->name, pp->site_index, tr(pp->q_label),
+                            pp->factor);
+                    }
+                    pmat->data[iop] = sum(strings);
+                } else
+                    assert(false);
+            auto lop = make_shared<SymbolicRowVector<S2>>(
+                mpo->left_operator_names[ii]->n);
+            auto rop = make_shared<SymbolicColumnVector<S2>>(
+                mpo->right_operator_names[ii]->m);
+            for (int iop = 0; iop < lop->data.size(); iop++) {
+                const auto p = dynamic_pointer_cast<OpElement<S1, FL>>(
+                    mpo->left_operator_names[ii]->data[iop]);
+                lop->data[iop] = make_shared<OpElement<S2, FL>>(
+                    p->name, p->site_index, tr(p->q_label), p->factor);
+            }
+            for (int iop = 0; iop < rop->data.size(); iop++) {
+                const auto p = dynamic_pointer_cast<OpElement<S1, FL>>(
+                    mpo->right_operator_names[ii]->data[iop]);
+                rop->data[iop] = make_shared<OpElement<S2, FL>>(
+                    p->name, p->site_index, tr(p->q_label), p->factor);
+            }
+            rmpo->left_operator_names[ii] = lop;
+            rmpo->right_operator_names[ii] = rop;
+            shared_ptr<VectorAllocator<FP>> d_alloc =
+                make_shared<VectorAllocator<FP>>();
+            shared_ptr<StateInfo<S1>> conn =
+                TransStateInfo<S2, S1>::backward_connection(mpo->basis[ii],
+                                                            rmpo->basis[ii]);
+            for (const auto &op : mopt->ops) {
+                const auto p =
+                    dynamic_pointer_cast<OpElement<S1, FL>>(op.first);
+                auto q = make_shared<OpElement<S2, FL>>(
+                    p->name, p->site_index, tr(p->q_label), p->factor);
+                shared_ptr<SparseMatrix<S2, FL>> xmat =
+                    make_shared<SparseMatrix<S2, FL>>(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<FL> r = (*op.second)[k];
+                    shared_ptr<StateInfo<S2>> mls =
+                        TransStateInfo<S1, S2>::forward(
+                            make_shared<StateInfo<S1>>(plu), ref);
+                    shared_ptr<StateInfo<S2>> mrs =
+                        TransStateInfo<S1, S2>::forward(
+                            make_shared<StateInfo<S1>>(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<FL> 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(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);
+                        }
+                }
+                opt->ops[q] = xmat;
+            }
+            mpo->unload_tensor(ii);
+            mpo->unload_left_operators(ii);
+            mpo->unload_right_operators(ii);
+            rmpo->save_tensor(ii);
+            rmpo->unload_tensor(ii);
+            rmpo->save_left_operators(ii);
+            rmpo->unload_left_operators(ii);
+            rmpo->save_right_operators(ii);
+            rmpo->unload_right_operators(ii);
+        }
+        return rmpo;
+    }
+};
+
 } // namespace block2

src/instantiation/block2_dmrg.hpp

@@ -712,6 +712,7 @@ extern template struct block2::MPO<block2::SAny, double>;
 extern template struct block2::DiagonalMPO<block2::SAny, double>;
 extern template struct block2::AncillaMPO<block2::SAny, double>;
 extern template struct block2::IdentityAddedMPO<block2::SAny, double>;
+extern template struct block2::TransMPO<block2::SAny, block2::SAny, double>;
 
 // mpo_fusing.hpp
 extern template struct block2::StackedMPO<block2::SAny, double>;
@@ -1492,6 +1493,7 @@ extern template struct block2::MPO<block2::SAny, complex<double>>;
 extern template struct block2::DiagonalMPO<block2::SAny, complex<double>>;
 extern template struct block2::AncillaMPO<block2::SAny, complex<double>>;
 extern template struct block2::IdentityAddedMPO<block2::SAny, complex<double>>;
+extern template struct block2::TransMPO<block2::SAny, block2::SAny, complex<double>>;
 
 // mpo_fusing.hpp
 extern template struct block2::StackedMPO<block2::SAny, complex<double>>;

src/instantiation/dmrg_a/mpo.cpp

@@ -25,3 +25,4 @@ template struct block2::MPO<block2::SAny, double>;
 template struct block2::DiagonalMPO<block2::SAny, double>;
 template struct block2::AncillaMPO<block2::SAny, double>;
 template struct block2::IdentityAddedMPO<block2::SAny, double>;
+template struct block2::TransMPO<block2::SAny, block2::SAny, double>;

src/instantiation/dmrg_az/mpo.cpp

@@ -24,3 +24,4 @@ template struct block2::MPO<block2::SAny, complex<double>>;
 template struct block2::DiagonalMPO<block2::SAny, complex<double>>;
 template struct block2::AncillaMPO<block2::SAny, complex<double>>;
 template struct block2::IdentityAddedMPO<block2::SAny, complex<double>>;
+template struct block2::TransMPO<block2::SAny, block2::SAny, complex<double>>;

src/pybind.cpp

@@ -306,9 +306,12 @@ PYBIND11_MODULE(block2, m) {
     bind_trans_mps<SAny, SAny>(m_sany, "sany");
     bind_trans_multi_mps<SAny, SAny>(m_sany, "sany");
     bind_fl_trans_mps_spin_specific<SAny, SAny, double>(m_sany, "sany");
+    bind_fl_trans_mpo<SAny, SAny, double>(m_sany, "sany");
 #ifdef _USE_COMPLEX
     bind_dmrg<SAny, complex<double>>(m_sany_cpx, "SAny");
-    bind_fl_trans_mps_spin_specific<SAny, SAny, complex<double>>(m_sany_cpx, "sany");
+    bind_fl_trans_mps_spin_specific<SAny, SAny, complex<double>>(m_sany_cpx,
+                                                                 "sany");
+    bind_fl_trans_mpo<SAny, SAny, complex<double>>(m_sany_cpx, "sany");
 #endif
 #endif
 

src/pybind/dmrg_a/trans_mps.cpp

@@ -27,3 +27,6 @@ template auto
 bind_fl_trans_mps_spin_specific<SAny, SAny, double>(py::module &m,
                                                     const string &aux_name)
     -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));
+template auto bind_fl_trans_mpo<SAny, SAny, double>(py::module &m,
+                                                    const string &aux_name)
+    -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));

src/pybind/dmrg_az/trans_mps.cpp

@@ -23,3 +23,7 @@
 template auto bind_fl_trans_mps_spin_specific<SAny, SAny, complex<double>>(
     py::module &m, const string &aux_name)
     -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));
+template auto
+bind_fl_trans_mpo<SAny, SAny, complex<double>>(py::module &m,
+                                               const string &aux_name)
+    -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));

src/pybind/pybind_dmrg.hpp

@@ -2265,6 +2265,13 @@ void bind_fl_trans_mps_spin_specific(py::module &m, const string &aux_name) {
           &TransUnfusedMPS<S, T, FL>::forward);
 }
 
+template <typename S, typename T, typename FL>
+void bind_fl_trans_mpo(py::module &m, const string &aux_name) {
+
+    m.def(("trans_mpo_to_" + aux_name).c_str(), &TransMPO<S, T, FL>::forward,
+          py::arg("mpo"), py::arg("hamil"), py::arg("tag") = "");
+}
+
 template <typename S = void> void bind_dmrg_types(py::module &m) {
 
     py::enum_<TruncationTypes>(m, "TruncationTypes", py::arithmetic())
@@ -2928,6 +2935,9 @@ extern template auto
 bind_fl_trans_mps_spin_specific<SAny, SAny, double>(py::module &m,
                                                     const string &aux_name)
     -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));
+extern template auto
+bind_fl_trans_mpo<SAny, SAny, double>(py::module &m, const string &aux_name)
+    -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));
 
 #ifdef _USE_COMPLEX
 extern template void bind_fl_general<SAny, complex<double>>(py::module &m);
@@ -2962,6 +2972,10 @@ extern template auto
 bind_fl_trans_mps_spin_specific<SAny, SAny, complex<double>>(
     py::module &m, const string &aux_name)
     -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));
+extern template auto
+bind_fl_trans_mpo<SAny, SAny, complex<double>>(py::module &m,
+                                               const string &aux_name)
+    -> decltype(typename SAny::is_sany_t(typename SAny::is_sany_t()));
 #endif
 
 #endif