1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: asd_dmrg/block_ops_1.cc
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 #include <src/asd/dmrg/block_operators.h>
26 #include <src/asd/dmrg/dmrg_block.h>
27
28 using namespace std;
29 using namespace bagel;
30
BlockOperators1(shared_ptr<const DMRG_Block1> left,shared_ptr<DimerJop> jop,const double thresh)31 BlockOperators1::BlockOperators1(shared_ptr<const DMRG_Block1> left, shared_ptr<DimerJop> jop, const double thresh) : BlockOperators(thresh), left_(left) {
32 const int rnorb = jop->monomer_jop<0>()->nocc();
33 const int lnorb = jop->monomer_jop<1>()->nocc();
34
35 for (auto& binfo : left->blocks()) {
36 const BlockKey bk = binfo.key();
37 { // build pure parts
38 shared_ptr<const btas::Tensor3<double>> gamma_aa = left->coupling({GammaSQ::CreateAlpha,GammaSQ::AnnihilateAlpha}).at({bk,bk}).data;
39 shared_ptr<const btas::Tensor3<double>> gamma_bb = left->coupling({GammaSQ::CreateBeta, GammaSQ::AnnihilateBeta }).at({bk,bk}).data;
40
41 const int gsize = gamma_aa->extent(0) * gamma_aa->extent(1);
42 assert(gsize == binfo.nstates*binfo.nstates);
43
44 // 2e part
45 shared_ptr<Matrix> ham = left->h2e(bk)->copy();
46
47 // 1e part
48 const btas::Tensor3<double> gamma_aa_plus_bb(*gamma_aa + *gamma_bb);
49 assert(ham->size()==gsize);
50 shared_ptr<const Matrix> mo1e = jop->monomer_jop<1>()->mo1e()->matrix();
51 dgemv_("N", gsize, mo1e->size(), 1.0, gamma_aa_plus_bb.data(), gsize, mo1e->data(), 1, 1.0, ham->data(), 1);
52 ham_.emplace(bk, ham);
53
54 // now for Q parts
55 const Matrix& coulomb = *jop->coulomb_matrix<1,0,1,0>();
56 const Matrix& exchange = *jop->coulomb_matrix<0,1,1,0>();
57
58 auto Qaa = make_shared<btas::Tensor4<double>>(gamma_aa->extent(0), gamma_aa->extent(1), rnorb, rnorb);
59 dgemm_("N", "T", gsize, rnorb*rnorb, lnorb*lnorb, 1.0, gamma_aa_plus_bb.data(), gsize, coulomb.data(), coulomb.ndim(),
60 0.0, Qaa->data(), gsize);
61 dgemm_("N", "T", gsize, rnorb*rnorb, lnorb*lnorb, -1.0, gamma_aa->data(), gsize, exchange.data(), exchange.ndim(),
62 1.0, Qaa->data(), gsize);
63 Q_aa_.emplace(bk, Qaa);
64
65 auto Qbb = make_shared<btas::Tensor4<double>>(gamma_bb->extent(0), gamma_bb->extent(1), rnorb, rnorb);
66 dgemm_("N", "T", gsize, rnorb*rnorb, lnorb*lnorb, 1.0, gamma_aa_plus_bb.data(), gsize, coulomb.data(), coulomb.ndim(),
67 0.0, Qbb->data(), gsize);
68 dgemm_("N", "T", gsize, rnorb*rnorb, lnorb*lnorb, -1.0, gamma_bb->data(), gsize, exchange.data(), exchange.ndim(),
69 1.0, Qbb->data(), gsize);
70 Q_bb_.emplace(bk, Qbb);
71
72 const BlockKey abkey(bk.nelea-1, bk.neleb+1);
73 if (left->contains(abkey)) {
74 shared_ptr<btas::Tensor3<double>> gamma_ab = left->coupling({GammaSQ::CreateBeta, GammaSQ::AnnihilateAlpha}).at({abkey,bk}).data;
75 auto Qab = make_shared<btas::Tensor4<double>>(gamma_ab->extent(0), gamma_ab->extent(1), rnorb, rnorb);
76 const int gabsize = gamma_ab->extent(0)*gamma_ab->extent(1);
77 dgemm_("N", "T", gabsize, rnorb*rnorb, lnorb*lnorb, -1.0, gamma_ab->data(), gabsize, exchange.data(), exchange.ndim(),
78 0.0, Qab->data(), gabsize);
79 Q_ab_.emplace(bk, Qab);
80 }
81 }
82
83 { // S_a and S_b
84 const Matrix& J_1101 = *jop->coulomb_matrix<1,1,0,1>();
85 const Matrix& h01 = *jop->cross_mo1e();
86
87 const BlockKey akey(bk.nelea+1, bk.neleb);
88 if (left->contains(akey)) {
89 shared_ptr<const btas::Tensor3<double>> gamma_a = left->coupling({GammaSQ::CreateAlpha}).at({akey,bk}).data;
90 shared_ptr<const btas::Tensor3<double>> gamma_aaa = left->coupling({GammaSQ::CreateAlpha, GammaSQ::CreateAlpha, GammaSQ::AnnihilateAlpha}).at({akey,bk}).data;
91 shared_ptr<const btas::Tensor3<double>> gamma_abb = left->coupling({GammaSQ::CreateAlpha, GammaSQ::CreateBeta, GammaSQ::AnnihilateBeta}).at({akey,bk}).data;
92
93 const int asize = gamma_aaa->extent(0) * gamma_aaa->extent(1);
94
95 // S_a
96 btas::Tensor3<double> gamma_akk(*gamma_aaa + *gamma_abb);
97 auto Sa = make_shared<btas::Tensor3<double>>(gamma_akk.extent(0), gamma_akk.extent(1), rnorb);
98 dgemm_("N", "T", asize, rnorb, lnorb, 1.0, gamma_a->data(), asize, h01.data(), h01.ndim(), 0.0, Sa->data(), asize);
99 dgemm_("N", "T", asize, rnorb, lnorb*lnorb*lnorb, 1.0, gamma_akk.data(), asize, J_1101.data(), J_1101.ndim(), 1.0, Sa->data(), asize);
100 S_a_.emplace(bk, Sa);
101 }
102
103 const BlockKey bkey(bk.nelea, bk.neleb+1);
104 if (left->contains(bkey)) {
105 shared_ptr<const btas::Tensor3<double>> gamma_b = left->coupling({GammaSQ::CreateBeta}).at({bkey,bk}).data;
106 shared_ptr<const btas::Tensor3<double>> gamma_baa = left->coupling({GammaSQ::CreateBeta, GammaSQ::CreateAlpha, GammaSQ::AnnihilateAlpha}).at({bkey, bk}).data;
107 shared_ptr<const btas::Tensor3<double>> gamma_bbb = left->coupling({GammaSQ::CreateBeta, GammaSQ::CreateBeta, GammaSQ::AnnihilateBeta}).at({bkey, bk}).data;
108
109 const int bsize = gamma_bbb->extent(0) * gamma_bbb->extent(1);
110
111 // S_b
112 btas::Tensor3<double> gamma_bkk(*gamma_baa + *gamma_bbb);
113 auto Sb = make_shared<btas::Tensor3<double>>(gamma_bkk.extent(0), gamma_bkk.extent(1), rnorb);
114 dgemm_("N", "T", bsize, rnorb, lnorb, 1.0, gamma_b->data(), bsize, h01.data(), h01.ndim(), 0.0, Sb->data(), bsize);
115 dgemm_("N", "T", bsize, rnorb, lnorb*lnorb*lnorb, 1.0, gamma_bkk.data(), bsize, J_1101.data(), J_1101.ndim(), 1.0, Sb->data(), bsize);
116 S_b_.emplace(bk, Sb);
117 }
118 }
119
120 { // P_aa, P_bb, P_ab
121 const Matrix& J_0110 = *jop->coulomb_matrix<0,1,1,0>();
122 auto compute_Pxx = [&J_0110, &lnorb, &rnorb, &left, &bk] (const BlockKey new_key, list<GammaSQ> ops, const double fac) {
123 shared_ptr<const btas::Tensor4<double>> gamma = left->coupling(ops).at({new_key, bk}).data;
124 auto Pxx = make_shared<btas::Tensor4<double>>(gamma->extent(0), gamma->extent(1), rnorb, rnorb);
125 const int gsize = gamma->extent(0)*gamma->extent(1);
126 dgemm_("N", "T", gsize, rnorb*rnorb, lnorb*lnorb, fac, gamma->data(), gsize, J_0110.data(), J_0110.ndim(), 0.0, Pxx->data(), gsize);
127 return Pxx;
128 };
129
130 const BlockKey aakey(bk.nelea-2,bk.neleb);
131 if (left->contains(aakey))
132 P_aa_.emplace(bk, compute_Pxx(aakey, {GammaSQ::AnnihilateAlpha, GammaSQ::AnnihilateAlpha}, 0.5));
133
134 const BlockKey bbkey(bk.nelea,bk.neleb-2);
135 if (left->contains(bbkey))
136 P_bb_.emplace(bk, compute_Pxx(bbkey, {GammaSQ::AnnihilateBeta, GammaSQ::AnnihilateBeta}, 0.5));
137
138 const BlockKey abkey(bk.nelea-1, bk.neleb-1);
139 if (left->contains(abkey))
140 P_ab_.emplace(bk, compute_Pxx(abkey, {GammaSQ::AnnihilateBeta, GammaSQ::AnnihilateAlpha}, 1.0));
141 }
142
143 { // D_a, D_b
144 const Matrix& J_0100 = *jop->coulomb_matrix<0,1,0,0>();
145
146 const BlockKey akey(bk.nelea+1, bk.neleb);
147 if (left->contains(akey)) {
148 shared_ptr<btas::Tensor3<double>> gamma_a = left->coupling({GammaSQ::CreateAlpha}).at({akey,bk}).data;
149 auto Da = make_shared<btas::TensorN<double,5>>(gamma_a->extent(0), gamma_a->extent(1), rnorb, rnorb, rnorb);
150 const int gsize = Da->extent(0)*Da->extent(1);
151 dgemm_("N", "T", gsize, rnorb*rnorb*rnorb, lnorb, 1.0, gamma_a->data(), gsize, J_0100.data(), J_0100.ndim(), 0.0, Da->data(), gsize);
152 D_a_.emplace(bk, Da);
153 }
154
155 const BlockKey bkey(bk.nelea, bk.neleb+1);
156 if (left->contains(bkey)) {
157 shared_ptr<btas::Tensor3<double>> gamma_b = left->coupling({GammaSQ::CreateBeta}).at({bkey,bk}).data;
158 auto Db = make_shared<btas::TensorN<double,5>>(gamma_b->extent(0), gamma_b->extent(1), rnorb, rnorb, rnorb);
159 const int gsize = Db->extent(0)*Db->extent(1);
160 dgemm_("N", "T", gsize, rnorb*rnorb*rnorb, lnorb, 1.0, gamma_b->data(), gsize, J_0100.data(), J_0100.ndim(), 0.0, Db->data(), gsize);
161 D_b_.emplace(bk, Db);
162 }
163 }
164 }
165 }
166
gamma_a(const BlockKey bk,int i) const167 shared_ptr<BlockSparseMatrix> BlockOperators1::gamma_a(const BlockKey bk, int i) const {
168 return get_sparse_mat_block(left_->coupling({GammaSQ::CreateAlpha}).at({BlockKey(bk.nelea+1,bk.neleb), bk}).data, i);
169 }
170
gamma_b(const BlockKey bk,int i) const171 shared_ptr<BlockSparseMatrix> BlockOperators1::gamma_b(const BlockKey bk, int i) const {
172 return get_sparse_mat_block(left_->coupling({GammaSQ::CreateBeta}).at({BlockKey(bk.nelea,bk.neleb+1), bk}).data, i);
173 }
174
gamma_a_as_matrix(const BlockKey bk,int i) const175 shared_ptr<Matrix> BlockOperators1::gamma_a_as_matrix(const BlockKey bk, int i) const {
176 return get_mat_block(left_->coupling({GammaSQ::CreateAlpha}).at({BlockKey(bk.nelea+1,bk.neleb), bk}).data, i);
177 }
178
gamma_b_as_matrix(const BlockKey bk,int i) const179 shared_ptr<Matrix> BlockOperators1::gamma_b_as_matrix(const BlockKey bk, int i) const {
180 return get_mat_block(left_->coupling({GammaSQ::CreateBeta}).at({BlockKey(bk.nelea,bk.neleb+1), bk}).data, i);
181 }
182