1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: nevpt2_mat.cc
4 // Copyright (C) 2014 Toru Shiozaki
5 //
6 // Author: Toru Shiozaki <shiozaki@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 #ifdef NEVPT2IMPL
26
27 #ifndef NDEBUG
is_hermitian_conjugate(shared_ptr<const T> t,shared_ptr<const T> u)28 template<typename T> inline bool is_hermitian_conjugate(shared_ptr<const T> t, shared_ptr<const T> u) { assert(false); return false; }
is_hermitian_conjugate(shared_ptr<const Matrix> t,shared_ptr<const Matrix> u)29 template<> inline bool is_hermitian_conjugate<Matrix>(shared_ptr<const Matrix> t, shared_ptr<const Matrix> u) {
30 return abs((*t->transpose()-*u).rms()) < 1.0e-8;
31 }
is_hermitian_conjugate(shared_ptr<const ZMatrix> t,shared_ptr<const ZMatrix> u)32 template<> inline bool is_hermitian_conjugate<ZMatrix>(shared_ptr<const ZMatrix> t, shared_ptr<const ZMatrix> u) {
33 return abs((*t->transpose_conjg()-*u).rms()) < 1.0e-8;
34 }
35 #endif
36
37 template<typename DataType>
compute_kmat()38 void NEVPT2<DataType>::compute_kmat() {
39 const double fac2 = is_same<DataType,double>::value ? 2.0 : 1.0;
40 {
41 // Eq. (27)
42 auto kmat = make_shared<MatType>(nact_, nact_, true);
43 btas::contract(1.0, *fockact_c_, {0,1}, *rdm1_, {2,1}, 0.0, *kmat, {0,2});
44 *kmat += *qvec_;
45 // transpose (assuming that kmat is hermitian...)
46 kmat = kmat->get_conjg();
47 kmat->localize();
48
49 // Eq. (A4)
50 auto kmatp = make_shared<MatType>(*kmat * (-1.0));
51 *kmatp += *fockact_->get_conjg() * fac2;
52 kmatp->localize();
53
54 kmat_ = kmat;
55 kmatp_ = kmatp;
56 }
57 {
58 auto compute_kmat = [this,&fac2](shared_ptr<const MatType> rdm2, shared_ptr<const MatType> rdm3, shared_ptr<const MatType> fock, const double sign) {
59 auto out = rdm2->clone();
60 // Eq. (A7) and (A9)
61 for (int b = 0; b != nact_; ++b)
62 for (int a = 0; a != nact_; ++a)
63 for (int bp = 0; bp != nact_; ++bp)
64 for (int ap = 0; ap != nact_; ++ap)
65 for (int c = 0; c != nact_; ++c) {
66 out->element(ap+nact_*bp, a+nact_*b) += rdm2->element(ap+nact_*bp, c+nact_*b) * fock->element(a, c)
67 + rdm2->element(ap+nact_*bp, a+nact_*c) * fock->element(b, c);
68 for (int d = 0; d != nact_; ++d)
69 for (int e = 0; e != nact_; ++e) {
70 out->element(ap+nact_*bp, a+nact_*b) += 0.5 * ints2_->element(e+nact_*a, c+nact_*d)
71 * (sign*2.0 * rdm3->element(ap+nact_*(bp+nact_*e), d+nact_*(b+nact_*c))
72 +sign*(b == e ? rdm2->element(ap+nact_*bp, d+nact_*c) : 0.0))
73 + 0.5 * ints2_->element(e+nact_*b, c+nact_*d)
74 * (sign*2.0 * rdm3->element(ap+nact_*(bp+nact_*e), a+nact_*(d+nact_*c))
75 +sign*(a == e ? rdm2->element(ap+nact_*bp, c+nact_*d) : 0.0));
76 }
77 }
78 return out;
79 };
80
81 kmat2_ = compute_kmat( rdm2_, rdm3_, fockact_c_, 1.0);
82 kmatp2_ = compute_kmat(hrdm2_, hrdm3_, fockact_h_, -1.0);
83 // for CAS references, these are Hermitian; if not, suspect the factor of 2.0 above
84 assert(kmat_->is_hermitian());
85 assert(kmatp_->is_hermitian());
86 assert(kmat2_->is_hermitian());
87 assert(kmatp2_->is_hermitian());
88 }
89 }
90
91
92 template<typename DataType>
compute_abcd()93 void NEVPT2<DataType>::compute_abcd() {
94 auto id2 = [this]( const int k, const int l) { return ( (k+nact_*l)); };
95 auto id3 = [this]( const int j, const int k, const int l) { return (j+nact_*(k+nact_*l)); };
96 auto id4 = [this](const int i, const int j, const int k, const int l) { return i+nact_*(j+nact_*(k+nact_*l)); };
97
98 const double fac2 = is_same<DataType,double>::value ? 2.0 : 1.0;
99 // A matrices
100 {
101 shared_ptr<MatType> amat2 = rdm2_->clone();
102 {
103 for (int b = 0; b != nact_; ++b)
104 for (int a = 0; a != nact_; ++a)
105 for (int bp = 0; bp != nact_; ++bp)
106 for (int ap = 0; ap != nact_; ++ap)
107 for (int c = 0; c != nact_; ++c) {
108 amat2->element(ap+nact_*bp,a+nact_*b) += fockact_p_->element(c,a) * ardm2_->element(bp+nact_*ap,c+nact_*b)
109 - fockact_p_->element(b,c) * ardm2_->element(bp+nact_*ap,a+nact_*c);
110 for (int d = 0; d != nact_; ++d)
111 for (int e = 0; e != nact_; ++e)
112 amat2->element(ap+nact_*bp,a+nact_*b) += 0.5 * ints2_->element(c+nact_*d,e+nact_*a) * (ardm3_->element(id3(bp,ap,c),id3(e,d,b))
113 + ardm3_->element(id3(bp,ap,d),id3(b,c,e)))
114 - 0.5 * ints2_->element(b+nact_*c,e+nact_*d) * (ardm3_->element(id3(bp,ap,a),id3(e,c,d))
115 + ardm3_->element(id3(bp,ap,c),id3(d,a,e)));
116 }
117 shared_ptr<MatType> tmp = amat2->copy();
118 sort_indices<1,0,3,2,0,1,1,1>(tmp->data(), amat2->data(), nact_, nact_, nact_, nact_);
119 }
120 shared_ptr<MatType> amat3 = rdm3_->clone();
121 shared_ptr<MatType> amat3t = rdm3_->clone();
122 {
123 const MatType fock_pp = *fockact_p_ * 2.0 - *fockact_c_;
124 for (int c = 0; c != nact_; ++c)
125 for (int b = 0; b != nact_; ++b)
126 for (int a = 0; a != nact_; ++a)
127 for (int cp = 0; cp != nact_; ++cp)
128 for (int bp = 0; bp != nact_; ++bp)
129 for (int ap = 0; ap != nact_; ++ap)
130 for (int d = 0; d != nact_; ++d) {
131 amat3->element(id3(ap,bp,cp),id3(a,b,c)) += fock_pp.element(d,a)*ardm3_->element(id3(cp,ap,bp),id3(b,d,c))
132 - fockact_c_->element(c,d)*ardm3_->element(id3(cp,ap,bp),id3(b,a,d))
133 - fock_pp.element(b,d)*ardm3_->element(id3(cp,ap,bp),id3(d,a,c));
134 amat3t->element(id3(ap,bp,cp),id3(a,b,c))+= fock_pp.element(d,a)*srdm3_->element(id3(cp,ap,bp),id3(b,d,c))
135 - fockact_c_->element(c,d)*srdm3_->element(id3(cp,ap,bp),id3(b,a,d))
136 + fockact_c_->element(d,b)*srdm3_->element(id3(cp,ap,bp),id3(d,a,c));
137 for (int e = 0; e != nact_; ++e) {
138 amat3->element(id3(ap,bp,cp),id3(a,b,c)) += ints2_->element(id2(c,d),id2(e,a))*ardm3_->element(id3(cp,ap,bp),id3(b,d,e));
139 amat3t->element(id3(ap,bp,cp),id3(a,b,c))+= ints2_->element(id2(c,d),id2(e,a))*srdm3_->element(id3(cp,ap,bp),id3(b,d,e));
140 for (int f = 0; f != nact_; ++f) {
141 amat3->element(id3(ap,bp,cp),id3(a,b,c)) += ints2_->element(id2(d,e),id2(f,a))*ardm4_->element(id4(cp,ap,bp,b),id4(d,f,e,c))
142 - ints2_->element(id2(d,c),id2(f,e))*ardm4_->element(id4(cp,ap,bp,b),id4(d,f,a,e))
143 - ints2_->element(id2(d,b),id2(f,e))*ardm4_->element(id4(cp,ap,bp,e),id4(d,f,a,c));
144 amat3t->element(id3(ap,bp,cp),id3(a,b,c))+= ints2_->element(id2(d,e),id2(f,a))
145 *((b == bp ? fac2 : 0.0)*ardm3_->element(id3(cp,ap,d),id3(f,e,c)) - ardm4_->element(id4(cp,ap,b,bp),id4(d,f,e,c)))
146 - ints2_->element(id2(d,c),id2(f,e))
147 *((b == bp ? fac2 : 0.0)*ardm3_->element(id3(cp,ap,d),id3(f,a,e)) - ardm4_->element(id4(cp,ap,b,bp),id4(d,f,a,e)))
148 + ints2_->element(id2(d,e),id2(f,b))
149 *((bp == e ? fac2 : 0.0)*ardm3_->element(id3(cp,ap,d),id3(f,a,c)) - ardm4_->element(id4(cp,ap,e,bp),id4(d,f,a,c)));
150 }
151 }
152 }
153 }
154 amat2_ = amat2;
155 assert(amat2_->is_hermitian());
156 amat3_ = amat3;
157 assert(amat3_->is_hermitian());
158 amat3t_ = amat3t;
159 assert(amat3t_->is_hermitian());
160 }
161 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
162 // B matrices
163 {
164 auto bmat2 = make_shared<MatType>(nact_*nact_*nact_, nact_, true);
165 auto bmat2t = make_shared<MatType>(nact_*nact_*nact_, nact_, true);
166 for (int a = 0; a != nact_; ++a)
167 for (int cp = 0; cp != nact_; ++cp)
168 for (int bp = 0; bp != nact_; ++bp)
169 for (int ap = 0; ap != nact_; ++ap)
170 for (int c = 0; c != nact_; ++c) {
171 bmat2->element(id3(ap,bp,cp),a) -= (fockact_p_->element(a,c)*2.0-fockact_c_->element(a,c)) * ardm2_->element(id2(cp,ap),id2(bp,c));
172 bmat2t->element(id3(ap,bp,cp),a) += fockact_c_->element(c,a) * srdm2_->element(id2(cp,ap),id2(bp,c));
173 for (int e = 0; e != nact_; ++e)
174 for (int f = 0; f != nact_; ++f) {
175 bmat2->element(id3(ap,bp,cp),a) -= ints2_->element(id2(a,c),id2(e,f)) * ardm3_->element(id3(cp,ap,bp),id3(e,c,f));
176 bmat2t->element(id3(ap,bp,cp),a) += ints2_->element(id2(e,c),id2(a,f)) * srdm3_->element(id3(cp,ap,bp),id3(e,c,f));
177 }
178 }
179 bmat2_ = bmat2;
180 bmat2t_ = bmat2t;
181 }
182 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
183 // C matrices
184 {
185 auto cmat2 = make_shared<MatType>(nact_, nact_*nact_*nact_, true);
186 auto cmat2t = make_shared<MatType>(nact_, nact_*nact_*nact_, true);
187 // <d c+ b+ a>
188 shared_ptr<MatType> s2rdm2 = ardm2_->clone();
189 for (int a = 0; a != nact_; ++a)
190 for (int b = 0; b != nact_; ++b)
191 for (int c = 0; c != nact_; ++c)
192 for (int d = 0; d != nact_; ++d)
193 s2rdm2->element(id2(d,c),id2(b,a)) += (d == c ? fac2 : 0.0) * rdm1_->element(b,a) - ardm2_->element(id2(c,d), id2(b,a));
194
195 const MatType fock_pp = *fockact_p_ * 2.0 - *fockact_c_;
196 for (int c = 0; c != nact_; ++c)
197 for (int a = 0; a != nact_; ++a)
198 for (int b = 0; b != nact_; ++b)
199 for (int ap = 0; ap != nact_; ++ap)
200 for (int d = 0; d != nact_; ++d) {
201 cmat2->element(ap, id3(a,b,c)) += fock_pp.element(d,a)*ardm2_->element(id2(ap,b),id2(d,c)) - fockact_c_->element(c,d)*ardm2_->element(id2(ap,b),id2(a,d))
202 - fock_pp.element(b,d)*ardm2_->element(id2(ap,d),id2(a,c));
203 cmat2t->element(ap, id3(a,b,c)) += fock_pp.element(d,a)*s2rdm2->element(id2(ap,b),id2(d,c)) - fockact_c_->element(c,d)*s2rdm2->element(id2(ap,b),id2(a,d))
204 + fockact_c_->element(d,b)*s2rdm2->element(id2(ap,d),id2(a,c));
205 for (int e = 0; e != nact_; ++e) {
206 for (int f = 0; f != nact_; ++f) {
207 cmat2->element(ap, id3(a,b,c)) += ints2_->element(id2(d,e),id2(f,a))*ardm3_->element(id3(ap,b,d),id3(f,e,c))
208 - ints2_->element(id2(d,c),id2(f,e))*ardm3_->element(id3(ap,b,d),id3(f,a,e))
209 - ints2_->element(id2(d,b),id2(f,e))*ardm3_->element(id3(ap,e,d),id3(f,a,c));
210 cmat2t->element(ap, id3(a,b,c))+= ints2_->element(id2(d,e),id2(f,a))*((ap == b ? fac2 : 0.0)*ardm2_->element(id2(d,f),id2(e,c)) - ardm3_->element(id3(b,ap,d),id3(f,e,c)))
211 - ints2_->element(id2(d,c),id2(f,e))*((ap == b ? fac2 : 0.0)*ardm2_->element(id2(d,f),id2(a,e)) - ardm3_->element(id3(b,ap,d),id3(f,a,e)))
212 + ints2_->element(id2(d,e),id2(f,b))*((ap == e ? fac2 : 0.0)*ardm2_->element(id2(d,f),id2(a,c)) - ardm3_->element(id3(e,ap,d),id3(f,a,c)));
213 }
214 cmat2->element(ap, id3(a,b,c)) += ints2_->element(id2(c,d),id2(e,a))*ardm2_->element(id2(ap,b),id2(d,e));
215 cmat2t->element(ap, id3(a,b,c))+= ints2_->element(id2(c,d),id2(e,a))*s2rdm2->element(id2(ap,b),id2(d,e));
216 }
217 }
218 cmat2_ = cmat2;
219 cmat2t_ = cmat2t;
220 }
221
222 // checking
223 assert(is_hermitian_conjugate(bmat2_, cmat2_));
224 assert(is_hermitian_conjugate(bmat2t_, cmat2t_));
225
226 /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
227 // D matrices
228 {
229 shared_ptr<MatType> s2rdm2 = ardm2_->clone();
230 shared_ptr<MatType> s2rdm3 = ardm3_->clone();
231 for (int a = 0; a != nact_; ++a)
232 for (int b = 0; b != nact_; ++b)
233 for (int c = 0; c != nact_; ++c)
234 for (int d = 0; d != nact_; ++d)
235 s2rdm2->element(id2(d,c),id2(b,a)) += (d == c ? fac2 : 0.0) * rdm1_->element(b,a) - ardm2_->element(id2(c,d), id2(b,a));
236 for (int a = 0; a != nact_; ++a)
237 for (int b = 0; b != nact_; ++b)
238 for (int c = 0; c != nact_; ++c)
239 for (int d = 0; d != nact_; ++d)
240 for (int e = 0; e != nact_; ++e)
241 for (int f = 0; f != nact_; ++f)
242 s2rdm3->element(id3(f,e,d),id3(c,b,a)) += (f == e ? fac2 : 0.0) * ardm2_->element(id2(d,c),id2(b,a)) - ardm3_->element(id3(e,f,d), id3(c,b,a));
243
244 shared_ptr<MatType> dmat2 = rdm2_->clone();
245 for (int b = 0; b != nact_; ++b)
246 for (int a = 0; a != nact_; ++a)
247 for (int bp = 0; bp != nact_; ++bp)
248 for (int ap = 0; ap != nact_; ++ap)
249 for (int c = 0; c != nact_; ++c) {
250 dmat2->element(ap+nact_*bp,a+nact_*b) -= fockact_p_->element(b,c) * ((a == ap ? fac2: 0.0) * rdm1_->element(bp,c) - rdm2_->element(a+nact_*bp,ap+nact_*c));
251 dmat2->element(ap+nact_*bp,a+nact_*b) += fockact_p_->element(c,a) * ((c == ap ? fac2: 0.0) * rdm1_->element(bp,b) - rdm2_->element(c+nact_*bp,ap+nact_*b));
252 for (int d = 0; d != nact_; ++d)
253 for (int e = 0; e != nact_; ++e) {
254 dmat2->element(ap+nact_*bp,a+nact_*b) -= 0.5 * ints2_->element(c+nact_*b,e+nact_*d)
255 * (srdm3_->element(id3(c,e,ap),id3(a,bp,d)) + s2rdm3->element(id3(ap,a,bp),id3(d,c,e))
256 + (ap == bp ? 1.0 : 0.0) *(ardm2_->element(c+nact_*e,a+nact_*d) + ardm2_->element(a+nact_*d,c+nact_*e))
257 + (c == ap ? 1.0 : 0.0) *s2rdm2->element(e+nact_*a,bp+nact_*d)
258 - (bp == e ? 1.0 : 0.0) *s2rdm2->element(ap+nact_*a,c+nact_*d));
259 dmat2->element(ap+nact_*bp,a+nact_*b) += 0.5 * ints2_->element(c+nact_*d,e+nact_*a)
260 * (srdm3_->element(id3(c,e,ap),id3(d,bp,b)) + s2rdm3->element(id3(ap,d,bp),id3(b,c,e))
261 + (ap == bp ? 1.0 : 0.0) *(ardm2_->element(c+nact_*e,d+nact_*b) + ardm2_->element(d+nact_*b,c+nact_*e))
262 + (c == ap ? 1.0 : 0.0) * s2rdm2->element(id2(e,d), id2(bp,b))
263 - (bp == e ? 1.0 : 0.0) * s2rdm2->element(id2(ap,d), id2(c,b)));
264
265 }
266 }
267 shared_ptr<MatType> tmp = dmat2->copy();
268 sort_indices<1,0,3,2,0,1,1,1>(tmp->data(), dmat2->data(), nact_, nact_, nact_, nact_);
269 dmat2_ = dmat2;
270 assert(dmat2_->is_hermitian());
271 }
272 {
273 shared_ptr<MatType> dmat1 = rdm1_->clone();
274 shared_ptr<MatType> dmat1t = rdm1_->clone();
275 for (int a = 0; a != nact_; ++a)
276 for (int ap = 0; ap != nact_; ++ap)
277 for (int c = 0; c != nact_; ++c) {
278 dmat1->element(ap,a) += -(fockact_p_->element(a,c)*2.0-fockact_c_->element(a,c)) * rdm1_->element(ap,c);
279 dmat1t->element(ap,a) += fockact_c_->element(c,a) * hrdm1_->element(c,ap);
280 for (int e = 0; e != nact_; ++e)
281 for (int f = 0; f != nact_; ++f) {
282 dmat1->element(ap,a) += - ints2_->element(id2(a,c),id2(e,f)) * ardm2_->element(id2(ap,e),id2(c,f));
283 dmat1t->element(ap,a) += ints2_->element(id2(e,c),id2(a,f)) * ((ap == e ? fac2 : 0.0)*rdm1_->element(c,f) - ardm2_->element(id2(e,ap),id2(c,f)));
284 }
285 }
286 dmat1_ = dmat1;
287 assert(dmat1_->is_hermitian());
288 dmat1t_ = dmat1t;
289 assert(dmat1t_->is_hermitian());
290 }
291 }
292
293 #endif
294