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