1 // 2 // BAGEL - Brilliantly Advanced General Electronic Structure Library 3 // Filename: asd_dmrg/block_operators.h 4 // Copyright (C) 2014 Toru Shiozaki 5 // 6 // Author: Shane Parker <shane.parker@u.northwestern.edu> 7 // Maintainer: Shiozaki group 8 // 9 // This file is part of the BAGEL package. 10 // 11 // This program is free software: you can redistribute it and/or modify 12 // it under the terms of the GNU General Public License as published by 13 // the Free Software Foundation, either version 3 of the License, or 14 // (at your option) any later version. 15 // 16 // This program is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU General Public License for more details. 20 // 21 // You should have received a copy of the GNU General Public License 22 // along with this program. If not, see <http://www.gnu.org/licenses/>. 23 // 24 25 #ifndef __BAGEL_ASD_DMRG_BLOCK_OPERATORS_H 26 #define __BAGEL_ASD_DMRG_BLOCK_OPERATORS_H 27 28 #include <src/asd/dimer/dimer_jop.h> 29 #include <src/asd/dmrg/block_key.h> 30 #include <src/util/math/blocksparsematrix.h> 31 32 namespace bagel { 33 34 class DMRG_Block; 35 class DMRG_Block1; 36 class DMRG_Block2; 37 38 /// Base class for access to contracted operators. All operators are keyed by the ket of the operation. The contracted operators that get built are the following: 39 /** \f[ \hat H^{L'L} = \sum_{p,q} h_{pq} \hat p^\dagger \hat q + \frac{1}{2} \sum_{pqrs} (pq|rs)p^\dagger r^\dagger s q \f] 40 \f[ \hat Q_{r_\alpha^\dagger s_\alpha}^{L'L} = \sum_{p,q} \left\{ \Gamma^{L'L}_{p^\dagger_\alpha q_\alpha} \left[(pq|rs) - (ps|qr)\right] + \Gamma^{L'L}_{p^\dagger_\beta q_\beta}(pq|rs)\right\} \f] 41 \f[ \hat Q_{r_\beta^\dagger s_\beta}^{L'L} = \sum_{p,q} \left\{ \Gamma^{L'L}_{p^\dagger_\beta q_\beta} \left[(pq|rs) - (ps|qr)\right] + \Gamma^{L'L}_{p^\dagger_\alpha q_\alpha}(pq|rs)\right\} \f] 42 \f[ \hat Q_{r_\alpha^\dagger s_\beta}^{L'L} = - \sum_{p,q} \Gamma^{L'L}_{p^\dagger_\beta q_\alpha} (ps|qr) \f] 43 \f[ \hat S_{p_\alpha}^{L'L} = \sum_q h_{pq} \Gamma^{L'L}_{q^\dagger_\alpha} + \sum_{\tau}\sum_{q,r,s} \Gamma^{L'L}_{s^\dagger_\alpha r^\dagger_\tau q_\tau} (ps|rq) \f] 44 \f[ \hat S_{p_\beta}^{L'L} = \sum_q h_{pq} \Gamma^{L'L}_{q^\dagger_\beta} + \sum_{\tau}\sum_{q,r,s} \Gamma^{L'L}_{s^\dagger_\beta r^\dagger_\tau q_\tau} (ps|rq) \f] 45 \f[ \hat P_{r_\alpha^\dagger s_\alpha^\dagger}^{L'L} = \sum_{p,q} \Gamma^{L'L}_{p_\alpha q_\alpha} (pr|qs) \f] 46 \f[ \hat P_{r_\beta^\dagger s_\beta^\dagger}^{L'L} = \sum_{p,q} \Gamma^{L'L}_{p_\beta q_\beta} (pr|qs) \f] 47 \f[ \hat P_{r_\alpha^\dagger s_\beta^\dagger}^{L'L} = \sum_{p,q} \Gamma^{L'L}_{p_\beta q_\alpha} (pr|qs) \f] 48 \f[ \hat D_{q_\tau^\dagger r_\tau s_\alpha}^{L'L} = \sum_s \Gamma^{L'L}_{p_\alpha^\dagger} (ps|rq) \f] 49 \f[ \hat D_{q_\tau^\dagger r_\tau s_\beta}^{L'L} = \sum_s \Gamma^{L'L}_{p_\beta^\dagger} (ps|rq) \f] */ 50 class BlockOperators { 51 protected: 52 double thresh_; ///< threshold for building BlockSparseMatrix objects 53 54 public: BlockOperators(const double thresh)55 BlockOperators(const double thresh) : thresh_(thresh) {} 56 57 virtual std::shared_ptr<Matrix> ham(const BlockKey bk) const = 0; 58 virtual std::shared_ptr<BlockSparseMatrix> Q_aa(const BlockKey bk, int i, int j) const = 0; 59 virtual std::shared_ptr<BlockSparseMatrix> Q_bb(const BlockKey bk, int i, int j) const = 0; 60 virtual std::shared_ptr<BlockSparseMatrix> Q_ab(const BlockKey bk, int i, int j) const = 0; 61 virtual std::shared_ptr<BlockSparseMatrix> S_a(const BlockKey bk, int i) const = 0; 62 virtual std::shared_ptr<BlockSparseMatrix> S_b(const BlockKey bk, int i) const = 0; 63 virtual std::shared_ptr<BlockSparseMatrix> P_aa(const BlockKey bk, int i, int j) const = 0; 64 virtual std::shared_ptr<BlockSparseMatrix> P_bb(const BlockKey bk, int i, int j) const = 0; 65 virtual std::shared_ptr<BlockSparseMatrix> P_ab(const BlockKey bk, int i, int j) const = 0; 66 virtual std::shared_ptr<BlockSparseMatrix> gamma_a(const BlockKey bk, int i) const = 0; 67 virtual std::shared_ptr<BlockSparseMatrix> gamma_b(const BlockKey bk, int i) const = 0; 68 69 virtual double ham(const BlockKey bk, int brastate, int ketstate) const = 0; 70 virtual double Q_aa(const BlockKey bk, int brastate, int ketstate, int i, int j) const = 0; 71 virtual double Q_bb(const BlockKey bk, int brastate, int ketstate, int i, int j) const = 0; 72 virtual double Q_ab(const BlockKey bk, int brastate, int ketstate, int i, int j) const = 0; 73 virtual double S_a(const BlockKey bk, int brastate, int ketstate, int i) const = 0; 74 virtual double S_b(const BlockKey bk, int brastate, int ketstate, int i) const = 0; 75 virtual double P_aa(const BlockKey bk, int brastate, int ketstate, int i, int j) const = 0; 76 virtual double P_bb(const BlockKey bk, int brastate, int ketstate, int i, int j) const = 0; 77 virtual double P_ab(const BlockKey bk, int brastate, int ketstate, int i, int j) const = 0; 78 virtual double D_a(const BlockKey bk, int brastate, int ketstate, int i, int j, int k) const = 0; 79 virtual double D_b(const BlockKey bk, int brastate, int ketstate, int i, int j, int k) const = 0; 80 }; 81 82 /// Operators for a single block 83 class BlockOperators1 : public BlockOperators { 84 protected: 85 std::shared_ptr<const DMRG_Block1> left_; 86 87 std::unordered_map<BlockKey, std::shared_ptr<const Matrix>> ham_; 88 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor4<double>>> Q_aa_; 89 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor4<double>>> Q_bb_; 90 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor4<double>>> Q_ab_; 91 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor3<double>>> S_a_; 92 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor3<double>>> S_b_; 93 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor4<double>>> P_aa_; 94 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor4<double>>> P_bb_; 95 std::unordered_map<BlockKey, std::shared_ptr<const btas::Tensor4<double>>> P_ab_; 96 std::unordered_map<BlockKey, std::shared_ptr<const btas::TensorN<double,5>>> D_a_; 97 std::unordered_map<BlockKey, std::shared_ptr<const btas::TensorN<double,5>>> D_b_; 98 99 public: 100 BlockOperators1(std::shared_ptr<const DMRG_Block1> left, std::shared_ptr<DimerJop> jop, const double thresh = 1.0e-12); 101 ham(const BlockKey bk)102 std::shared_ptr<Matrix> ham(const BlockKey bk) const override { return ham_.at(bk)->copy(); } Q_aa(const BlockKey bk,int i,int j)103 std::shared_ptr<BlockSparseMatrix> Q_aa(const BlockKey bk, int i, int j) const override { return get_sparse_mat_block(Q_aa_.at(bk), i, j); } Q_bb(const BlockKey bk,int i,int j)104 std::shared_ptr<BlockSparseMatrix> Q_bb(const BlockKey bk, int i, int j) const override { return get_sparse_mat_block(Q_bb_.at(bk), i, j); } Q_ab(const BlockKey bk,int i,int j)105 std::shared_ptr<BlockSparseMatrix> Q_ab(const BlockKey bk, int i, int j) const override { return get_sparse_mat_block(Q_ab_.at(bk), i, j); } S_a(const BlockKey bk,int i)106 std::shared_ptr<BlockSparseMatrix> S_a(const BlockKey bk, int i) const override { return get_sparse_mat_block(S_a_.at(bk), i); } S_b(const BlockKey bk,int i)107 std::shared_ptr<BlockSparseMatrix> S_b(const BlockKey bk, int i) const override { return get_sparse_mat_block(S_b_.at(bk), i); } P_aa(const BlockKey bk,int i,int j)108 std::shared_ptr<BlockSparseMatrix> P_aa(const BlockKey bk, int i, int j) const override { return get_sparse_mat_block(P_aa_.at(bk), i, j); } P_bb(const BlockKey bk,int i,int j)109 std::shared_ptr<BlockSparseMatrix> P_bb(const BlockKey bk, int i, int j) const override { return get_sparse_mat_block(P_bb_.at(bk), i, j); } P_ab(const BlockKey bk,int i,int j)110 std::shared_ptr<BlockSparseMatrix> P_ab(const BlockKey bk, int i, int j) const override { return get_sparse_mat_block(P_ab_.at(bk), i, j); } 111 std::shared_ptr<BlockSparseMatrix> gamma_a(const BlockKey bk, int i) const override; 112 std::shared_ptr<BlockSparseMatrix> gamma_b(const BlockKey bk, int i) const override; 113 114 std::shared_ptr<Matrix> gamma_a_as_matrix(const BlockKey bk, int i) const; 115 std::shared_ptr<Matrix> gamma_b_as_matrix(const BlockKey bk, int i) const; 116 S_a_copy_to(double * target,const BlockKey bk,int i)117 void S_a_copy_to(double* target, const BlockKey bk, int i) const { copy_to(target, S_a_.at(bk), i); } S_b_copy_to(double * target,const BlockKey bk,int i)118 void S_b_copy_to(double* target, const BlockKey bk, int i) const { copy_to(target, S_b_.at(bk), i); } Q_aa_copy_to(double * target,const BlockKey bk,int i,int j)119 void Q_aa_copy_to(double* target, const BlockKey bk, int i, int j) const { copy_to(target, Q_aa_.at(bk), i, j); } Q_bb_copy_to(double * target,const BlockKey bk,int i,int j)120 void Q_bb_copy_to(double* target, const BlockKey bk, int i, int j) const { copy_to(target, Q_bb_.at(bk), i, j); } Q_ab_copy_to(double * target,const BlockKey bk,int i,int j)121 void Q_ab_copy_to(double* target, const BlockKey bk, int i, int j) const { copy_to(target, Q_ab_.at(bk), i, j); } P_aa_copy_to(double * target,const BlockKey bk,int i,int j)122 void P_aa_copy_to(double* target, const BlockKey bk, int i, int j) const { copy_to(target, P_aa_.at(bk), i, j); } P_bb_copy_to(double * target,const BlockKey bk,int i,int j)123 void P_bb_copy_to(double* target, const BlockKey bk, int i, int j) const { copy_to(target, P_bb_.at(bk), i, j); } P_ab_copy_to(double * target,const BlockKey bk,int i,int j)124 void P_ab_copy_to(double* target, const BlockKey bk, int i, int j) const { copy_to(target, P_ab_.at(bk), i, j); } 125 ham(const BlockKey bk,int brastate,int ketstate)126 double ham(const BlockKey bk, int brastate, int ketstate) const override { return ham_.at(bk)->element(brastate,ketstate); } Q_aa(const BlockKey bk,int brastate,int ketstate,int i,int j)127 double Q_aa(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return (*Q_aa_.at(bk))(brastate, ketstate, i, j); } Q_bb(const BlockKey bk,int brastate,int ketstate,int i,int j)128 double Q_bb(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return (*Q_bb_.at(bk))(brastate, ketstate, i, j); } Q_ab(const BlockKey bk,int brastate,int ketstate,int i,int j)129 double Q_ab(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return (*Q_ab_.at(bk))(brastate, ketstate, i, j); } S_a(const BlockKey bk,int brastate,int ketstate,int i)130 double S_a(const BlockKey bk, int brastate, int ketstate, int i) const override { return (*S_a_.at(bk))(brastate, ketstate, i); } S_b(const BlockKey bk,int brastate,int ketstate,int i)131 double S_b(const BlockKey bk, int brastate, int ketstate, int i) const override { return (*S_b_.at(bk))(brastate, ketstate, i); } P_aa(const BlockKey bk,int brastate,int ketstate,int i,int j)132 double P_aa(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return (*P_aa_.at(bk))(brastate, ketstate, i, j); } P_bb(const BlockKey bk,int brastate,int ketstate,int i,int j)133 double P_bb(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return (*P_bb_.at(bk))(brastate, ketstate, i, j); } P_ab(const BlockKey bk,int brastate,int ketstate,int i,int j)134 double P_ab(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return (*P_ab_.at(bk))(brastate, ketstate, i, j); } D_a(const BlockKey bk,int brastate,int ketstate,int i,int j,int k)135 double D_a(const BlockKey bk, int brastate, int ketstate, int i, int j, int k) const override { return (*D_a_.at(bk))(brastate, ketstate, i, j, k); } D_b(const BlockKey bk,int brastate,int ketstate,int i,int j,int k)136 double D_b(const BlockKey bk, int brastate, int ketstate, int i, int j, int k) const override { return (*D_b_.at(bk))(brastate, ketstate, i, j, k); } 137 Q_aa_as_matview(const BlockKey bk)138 const MatView Q_aa_as_matview(const BlockKey bk) const { return get_as_matview(Q_aa_.at(bk)); } Q_bb_as_matview(const BlockKey bk)139 const MatView Q_bb_as_matview(const BlockKey bk) const { return get_as_matview(Q_bb_.at(bk)); } Q_ab_as_matview(const BlockKey bk)140 const MatView Q_ab_as_matview(const BlockKey bk) const { return get_as_matview(Q_ab_.at(bk)); } S_a_as_matview(const BlockKey bk)141 const MatView S_a_as_matview(const BlockKey bk) const { return get_as_matview(S_a_.at(bk)); } S_b_as_matview(const BlockKey bk)142 const MatView S_b_as_matview(const BlockKey bk) const { return get_as_matview(S_b_.at(bk)); } P_aa_as_matview(const BlockKey bk)143 const MatView P_aa_as_matview(const BlockKey bk) const { return get_as_matview(P_aa_.at(bk)); } P_bb_as_matview(const BlockKey bk)144 const MatView P_bb_as_matview(const BlockKey bk) const { return get_as_matview(P_bb_.at(bk)); } P_ab_as_matview(const BlockKey bk)145 const MatView P_ab_as_matview(const BlockKey bk) const { return get_as_matview(P_ab_.at(bk)); } D_a_as_matview(const BlockKey bk)146 const MatView D_a_as_matview(const BlockKey bk) const { return get_as_matview(D_a_.at(bk)); } D_b_as_matview(const BlockKey bk)147 const MatView D_b_as_matview(const BlockKey bk) const { return get_as_matview(D_b_.at(bk)); } 148 149 private: 150 template <typename TensorType> get_as_matview(std::shared_ptr<const TensorType> t)151 const MatView get_as_matview(std::shared_ptr<const TensorType> t) const { 152 return MatView(btas::make_view(btas::CRange<2>(t->extent(0)*t->extent(1), t->size()/(t->extent(0)*t->extent(1))), t->storage()), true); 153 } 154 155 template <typename... Args> get_sparse_mat_block(Args...args)156 std::shared_ptr<BlockSparseMatrix> get_sparse_mat_block(Args... args) const { 157 return std::make_shared<BlockSparseMatrix>(get_mat_block(args...)); 158 } 159 get_mat_block(std::shared_ptr<const btas::Tensor3<double>> g,const int i)160 std::shared_ptr<Matrix> get_mat_block(std::shared_ptr<const btas::Tensor3<double>> g, const int i) const { 161 auto out = std::make_shared<Matrix>(g->extent(0), g->extent(1), true); 162 std::copy_n(&(*g)(0, 0, i), out->size(), out->data()); 163 return out; 164 } 165 get_mat_block(std::shared_ptr<const btas::Tensor4<double>> g,const int i,const int j)166 std::shared_ptr<Matrix> get_mat_block(std::shared_ptr<const btas::Tensor4<double>> g, const int i, const int j) const { 167 auto out = std::make_shared<Matrix>(g->extent(0), g->extent(1), true); 168 std::copy_n(&(*g)(0, 0, i, j), out->size(), out->data()); 169 return out; 170 } 171 get_mat_block(std::shared_ptr<const btas::TensorN<double,5>> g,const int i,const int j,const int k)172 std::shared_ptr<Matrix> get_mat_block(std::shared_ptr<const btas::TensorN<double,5>> g, const int i, const int j, const int k) const { 173 auto out = std::make_shared<Matrix>(g->extent(0), g->extent(1), true); 174 std::copy_n(&(*g)(0, 0, i, j, k), out->size(), out->data()); 175 return out; 176 } 177 copy_to(double * target,std::shared_ptr<const btas::Tensor3<double>> g,const int i)178 void copy_to(double* target, std::shared_ptr<const btas::Tensor3<double>> g, const int i) const { 179 std::copy_n(&(*g)(0, 0, i), g->extent(0)*g->extent(1), target); 180 } 181 copy_to(double * target,std::shared_ptr<const btas::Tensor4<double>> g,const int i,const int j)182 void copy_to(double* target, std::shared_ptr<const btas::Tensor4<double>> g, const int i, const int j) const { 183 std::copy_n(&(*g)(0, 0, i, j), g->extent(0)*g->extent(1), target); 184 } 185 186 }; 187 188 /// Operators for two blocks 189 class BlockOperators2 : public BlockOperators { 190 protected: 191 std::shared_ptr<const DMRG_Block2> blocks_; 192 std::shared_ptr<const BlockOperators1> left_ops_; 193 std::shared_ptr<const BlockOperators1> right_ops_; 194 std::shared_ptr<const BlockOperators1> intra_ops_; 195 std::shared_ptr<DimerJop> jop_; 196 197 public: 198 BlockOperators2(std::shared_ptr<const DMRG_Block2> blocks, std::shared_ptr<DimerJop> jop, const double thresh = 1.0e-13); 199 200 std::shared_ptr<Matrix> ham(const BlockKey bk) const override; 201 std::shared_ptr<BlockSparseMatrix> Q_aa(const BlockKey bk, int i, int j) const override; 202 std::shared_ptr<BlockSparseMatrix> Q_bb(const BlockKey bk, int i, int j) const override; 203 std::shared_ptr<BlockSparseMatrix> Q_ab(const BlockKey bk, int i, int j) const override; 204 std::shared_ptr<BlockSparseMatrix> S_a(const BlockKey bk, int i) const override; 205 std::shared_ptr<BlockSparseMatrix> S_b(const BlockKey bk, int i) const override; 206 std::shared_ptr<BlockSparseMatrix> P_aa(const BlockKey bk, int i, int j) const override; 207 std::shared_ptr<BlockSparseMatrix> P_bb(const BlockKey bk, int i, int j) const override; 208 std::shared_ptr<BlockSparseMatrix> P_ab(const BlockKey bk, int i, int j) const override; 209 std::shared_ptr<BlockSparseMatrix> gamma_a(const BlockKey bk, int i) const override; 210 std::shared_ptr<BlockSparseMatrix> gamma_b(const BlockKey bk, int i) const override; 211 212 // these should all be replaced with more direct functions ham(const BlockKey bk,int brastate,int ketstate)213 double ham(const BlockKey bk, int brastate, int ketstate) const override { return ham(bk)->element(brastate,ketstate); } Q_aa(const BlockKey bk,int brastate,int ketstate,int i,int j)214 double Q_aa(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return Q_aa(bk, i, j)->element(brastate, ketstate); } Q_bb(const BlockKey bk,int brastate,int ketstate,int i,int j)215 double Q_bb(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return Q_bb(bk, i, j)->element(brastate, ketstate); } Q_ab(const BlockKey bk,int brastate,int ketstate,int i,int j)216 double Q_ab(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return Q_ab(bk, i, j)->element(brastate, ketstate); } S_a(const BlockKey bk,int brastate,int ketstate,int i)217 double S_a(const BlockKey bk, int brastate, int ketstate, int i) const override { return S_a(bk, i)->element(brastate, ketstate); } S_b(const BlockKey bk,int brastate,int ketstate,int i)218 double S_b(const BlockKey bk, int brastate, int ketstate, int i) const override { return S_b(bk, i)->element(brastate, ketstate); } P_aa(const BlockKey bk,int brastate,int ketstate,int i,int j)219 double P_aa(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return P_aa(bk, i, j)->element(brastate, ketstate); } P_bb(const BlockKey bk,int brastate,int ketstate,int i,int j)220 double P_bb(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return P_bb(bk, i, j)->element(brastate, ketstate); } P_ab(const BlockKey bk,int brastate,int ketstate,int i,int j)221 double P_ab(const BlockKey bk, int brastate, int ketstate, int i, int j) const override { return P_ab(bk, i, j)->element(brastate, ketstate); } 222 double D_a(const BlockKey bk, int brastate, int ketstate, int i, int j, int k) const override; 223 double D_b(const BlockKey bk, int brastate, int ketstate, int i, int j, int k) const override; 224 }; 225 226 } 227 228 #endif 229 230