1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2015
9 * Copyright (c) 2010-2015, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16
17 #include "eri_digest.h"
18 #include "eriworker.h"
19
IntegralDigestor()20 IntegralDigestor::IntegralDigestor() {
21 }
22
~IntegralDigestor()23 IntegralDigestor::~IntegralDigestor() {
24 }
25
JDigestor(const arma::mat & P_)26 JDigestor::JDigestor(const arma::mat & P_) : P(P_) {
27 J.zeros(P.n_rows,P.n_cols);
28 }
29
~JDigestor()30 JDigestor::~JDigestor() {
31 }
32
digest(const std::vector<eripair_t> & shpairs,size_t ip,size_t jp,const std::vector<double> & ints,size_t ioff)33 void JDigestor::digest(const std::vector<eripair_t> & shpairs, size_t ip, size_t jp, const std::vector<double> & ints, size_t ioff) {
34 // Shells in quartet are
35 size_t is=shpairs[ip].is;
36 size_t js=shpairs[ip].js;
37 size_t ks=shpairs[jp].is;
38 size_t ls=shpairs[jp].js;
39
40 // Amount of functions on the first pair is
41 size_t Ni=shpairs[ip].Ni;
42 size_t Nj=shpairs[ip].Nj;
43 // and on second pair is
44 size_t Nk=shpairs[jp].Ni;
45 size_t Nl=shpairs[jp].Nj;
46
47 // First functions on the first pair is
48 size_t i0=shpairs[ip].i0;
49 size_t j0=shpairs[ip].j0;
50 // Second pair is
51 size_t k0=shpairs[jp].i0;
52 size_t l0=shpairs[jp].j0;
53
54 // J_ij = (ij|kl) P_kl
55 {
56 // Work matrix
57 arma::mat Jij(Ni,Nj);
58 Jij.zeros();
59 arma::mat Pkl=P.submat(k0,l0,k0+Nk-1,l0+Nl-1);
60
61 // Degeneracy factor
62 double fac=1.0;
63 if(ks!=ls)
64 fac=2.0;
65
66 // Increment matrix
67 for(size_t ii=0;ii<Ni;ii++)
68 for(size_t jj=0;jj<Nj;jj++) {
69
70 // Matrix element
71 double el=0.0;
72 for(size_t kk=0;kk<Nk;kk++)
73 for(size_t ll=0;ll<Nl;ll++)
74 el+=Pkl(kk,ll)*ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll];
75
76 // Set the element
77 Jij(ii,jj)+=fac*el;
78 }
79
80 // Store the matrix element
81 J.submat(i0,j0,i0+Ni-1,j0+Nj-1)+=Jij;
82 if(is!=js)
83 J.submat(j0,i0,j0+Nj-1,i0+Ni-1)+=arma::trans(Jij);
84 }
85
86 // Permutation: J_kl = (ij|kl) P_ij
87 if(ip!=jp) {
88 // Work matrix
89 arma::mat Jkl(Nk,Nl);
90 Jkl.zeros();
91 arma::mat Pij=P.submat(i0,j0,i0+Ni-1,j0+Nj-1);
92
93 // Degeneracy factor
94 double fac=1.0;
95 if(is!=js)
96 fac=2.0;
97
98 // Increment matrix
99 for(size_t kk=0;kk<Nk;kk++)
100 for(size_t ll=0;ll<Nl;ll++) {
101
102 // Matrix element
103 double el=0.0;
104 for(size_t ii=0;ii<Ni;ii++) {
105 for(size_t jj=0;jj<Nj;jj++) {
106 el+=Pij(ii,jj)*ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll];
107 }
108 }
109
110 // Set the element
111 Jkl(kk,ll)+=fac*el;
112 }
113
114 // Store the matrix element
115 J.submat(k0,l0,k0+Nk-1,l0+Nl-1)+=Jkl;
116 if(ks!=ls)
117 J.submat(l0,k0,l0+Nl-1,k0+Nk-1)+=arma::trans(Jkl);
118 }
119 }
120
get_J() const121 arma::mat JDigestor::get_J() const {
122 return J;
123 }
124
KDigestor(const arma::mat & P_)125 KDigestor::KDigestor(const arma::mat & P_) : P(P_) {
126 K.zeros(P.n_rows,P.n_cols);
127 }
128
~KDigestor()129 KDigestor::~KDigestor() {
130 }
131
digest(const std::vector<eripair_t> & shpairs,size_t ip,size_t jp,const std::vector<double> & ints,size_t ioff)132 void KDigestor::digest(const std::vector<eripair_t> & shpairs, size_t ip, size_t jp, const std::vector<double> & ints, size_t ioff) {
133 size_t is=shpairs[ip].is;
134 size_t js=shpairs[ip].js;
135 size_t ks=shpairs[jp].is;
136 size_t ls=shpairs[jp].js;
137
138 // Amount of functions on the first pair is
139 size_t Ni=shpairs[ip].Ni;
140 size_t Nj=shpairs[ip].Nj;
141 // and on second pair is
142 size_t Nk=shpairs[jp].Ni;
143 size_t Nl=shpairs[jp].Nj;
144
145 // First functions on the first pair is
146 size_t i0=shpairs[ip].i0;
147 size_t j0=shpairs[ip].j0;
148 // Second pair is
149 size_t k0=shpairs[jp].i0;
150 size_t l0=shpairs[jp].j0;
151
152 /*
153 When all indices are different, the
154 following integrals are equivalent:
155 (ij|kl) (ij|lk) (ji|kl) (ji|lk)
156 (kl|ij) (kl|ji) (lk|ij) (lk|ji)
157
158 This translates to
159
160 K(i,k) += (ij|kl) P(j,l) // always
161 K(j,k) += (ij|kl) P(i,l) // if (is!=js)
162 K(i,l) += (ij|kl) P(j,k) // if (ls!=ks)
163 K(j,l) += (ij|kl) P(i,k) // if (is!=js) and (ls!=ks)
164
165 and for ij != kl
166
167 K(k,i) += (ij|kl) P(j,l) // always
168 K(k,j) += (ij|kl) P(i,l) // if (is!=js)
169 K(l,i) += (ij|kl) P(j,k) // if (ks!=ls)
170 K(l,j) += (ij|kl) P(i,k) // if (is!=js) and (ks!=ls)
171
172 However, the latter four permutations just make the
173 exchange matrix symmetric. So the only thing we need to do
174 is do the first four permutations, and at the end we sum up
175 K_ij and K_ji for j>i and set K_ij and K_ji to this
176 value. This makes things a *lot* easier. So:
177 We just need to check if the shells are different, in which
178 case K will get extra increments.
179 */
180
181 // First, do the ik part: K(i,k) += (ij|kl) P(j,l)
182 {
183 arma::mat Kik(Ni,Nk);
184 Kik.zeros();
185 arma::mat Pjl =P.submat(j0,l0,j0+Nj-1,l0+Nl-1);
186
187 // Increment Kik
188 for(size_t ii=0;ii<Ni;ii++)
189 for(size_t kk=0;kk<Nk;kk++)
190 for(size_t ll=0;ll<Nl;ll++)
191 for(size_t jj=0;jj<Nj;jj++)
192 Kik (ii,kk)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pjl (jj,ll);
193
194 // Set elements
195 K.submat(i0,k0,i0+Ni-1,k0+Nk-1)+=Kik;
196 // Symmetrize if necessary
197 if(ip!=jp)
198 K.submat(k0,i0,k0+Nk-1,i0+Ni-1)+=arma::trans(Kik);
199 }
200
201 // Then, the second part: K(j,k) += (ij|kl) P(i,l)
202 if(is!=js) {
203 arma::mat Kjk(Nj,Nk);
204 Kjk.zeros();
205 arma::mat Pil=P.submat(i0,l0,i0+Ni-1,l0+Nl-1);
206
207 // Increment Kjk
208 for(size_t jj=0;jj<Nj;jj++)
209 for(size_t kk=0;kk<Nk;kk++)
210 for(size_t ll=0;ll<Nl;ll++)
211 for(size_t ii=0;ii<Ni;ii++)
212 Kjk(jj,kk)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pil(ii,ll);
213
214 // Set elements
215 K.submat(j0,k0,j0+Nj-1,k0+Nk-1)+=Kjk;
216 // Symmetrize if necessary (take care about possible overlap with next routine)
217 if(ip!=jp) {
218 K.submat(k0,j0,k0+Nk-1,j0+Nj-1)+=arma::trans(Kjk);
219 }
220 }
221
222 // Third part: K(i,l) += (ij|kl) P(j,k)
223 if(ks!=ls) {
224 arma::mat Kil(Ni,Nl);
225 Kil.zeros();
226 arma::mat Pjk=P.submat(j0,k0,j0+Nj-1,k0+Nk-1);
227
228 // Increment Kil
229 for(size_t ii=0;ii<Ni;ii++)
230 for(size_t ll=0;ll<Nl;ll++)
231 for(size_t jj=0;jj<Nj;jj++)
232 for(size_t kk=0;kk<Nk;kk++) {
233 Kil(ii,ll)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pjk(jj,kk);
234 }
235
236 // Set elements
237 K.submat(i0,l0,i0+Ni-1,l0+Nl-1)+=Kil;
238 // Symmetrize if necessary
239 if(ip!=jp)
240 K.submat(l0,i0,l0+Nl-1,i0+Ni-1)+=arma::trans(Kil);
241 }
242
243 // Last permutation: K(j,l) += (ij|kl) P(i,k)
244 if(is!=js && ks!=ls) {
245 arma::mat Kjl(Nj,Nl);
246 Kjl.zeros();
247 arma::mat Pik=P.submat(i0,k0,i0+Ni-1,k0+Nk-1);
248
249 // Increment Kjl
250 for(size_t jj=0;jj<Nj;jj++)
251 for(size_t ll=0;ll<Nl;ll++)
252 for(size_t ii=0;ii<Ni;ii++)
253 for(size_t kk=0;kk<Nk;kk++) {
254 Kjl(jj,ll)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pik(ii,kk);
255 }
256
257 // Set elements
258 K.submat(j0,l0,j0+Nj-1,l0+Nl-1)+=Kjl;
259 // Symmetrize if necessary
260 if (ip!=jp)
261 K.submat(l0,j0,l0+Nl-1,j0+Nj-1)+=arma::trans(Kjl);
262 }
263 }
264
get_K() const265 arma::mat KDigestor::get_K() const {
266 return K;
267 }
268
cxKDigestor(const arma::cx_mat & P_)269 cxKDigestor::cxKDigestor(const arma::cx_mat & P_) : P(P_) {
270 K.zeros(P.n_rows,P.n_cols);
271 }
272
~cxKDigestor()273 cxKDigestor::~cxKDigestor() {
274 }
275
digest(const std::vector<eripair_t> & shpairs,size_t ip,size_t jp,const std::vector<double> & ints,size_t ioff)276 void cxKDigestor::digest(const std::vector<eripair_t> & shpairs, size_t ip, size_t jp, const std::vector<double> & ints, size_t ioff) {
277 size_t is=shpairs[ip].is;
278 size_t js=shpairs[ip].js;
279 size_t ks=shpairs[jp].is;
280 size_t ls=shpairs[jp].js;
281
282 // Amount of functions on the first pair is
283 size_t Ni=shpairs[ip].Ni;
284 size_t Nj=shpairs[ip].Nj;
285 // and on second pair is
286 size_t Nk=shpairs[jp].Ni;
287 size_t Nl=shpairs[jp].Nj;
288
289 // First functions on the first pair is
290 size_t i0=shpairs[ip].i0;
291 size_t j0=shpairs[ip].j0;
292 // Second pair is
293 size_t k0=shpairs[jp].i0;
294 size_t l0=shpairs[jp].j0;
295
296 /*
297 When all indices are different, the
298 following integrals are equivalent:
299 (ij|kl) (ij|lk) (ji|kl) (ji|lk)
300 (kl|ij) (kl|ji) (lk|ij) (lk|ji)
301
302 This translates to
303
304 K(i,k) += (ij|kl) P(j,l) // always
305 K(j,k) += (ij|kl) P(i,l) // if (is!=js)
306 K(i,l) += (ij|kl) P(j,k) // if (ls!=ks)
307 K(j,l) += (ij|kl) P(i,k) // if (is!=js) and (ls!=ks)
308
309 and for ij != kl
310
311 K(k,i) += (ij|kl) P(j,l) // always
312 K(k,j) += (ij|kl) P(i,l) // if (is!=js)
313 K(l,i) += (ij|kl) P(j,k) // if (ks!=ls)
314 K(l,j) += (ij|kl) P(i,k) // if (is!=js) and (ks!=ls)
315
316 However, the latter four permutations just make the
317 exchange matrix symmetric. So the only thing we need to do
318 is do the first four permutations, and at the end we sum up
319 K_ij and K_ji for j>i and set K_ij and K_ji to this
320 value. This makes things a *lot* easier. So:
321 We just need to check if the shells are different, in which
322 case K will get extra increments.
323 */
324
325 // First, do the ik part: K(i,k) += (ij|kl) P(j,l)
326 {
327 arma::cx_mat Kik(Ni,Nk);
328 Kik.zeros();
329 arma::cx_mat Pjl =P.submat(j0,l0,j0+Nj-1,l0+Nl-1);
330
331 // Increment Kik
332 for(size_t ii=0;ii<Ni;ii++)
333 for(size_t kk=0;kk<Nk;kk++)
334 for(size_t ll=0;ll<Nl;ll++)
335 for(size_t jj=0;jj<Nj;jj++)
336 Kik (ii,kk)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pjl (jj,ll);
337
338 // Set elements
339 K.submat(i0,k0,i0+Ni-1,k0+Nk-1)+=Kik;
340 // Symmetrize if necessary
341 if(ip!=jp)
342 K.submat(k0,i0,k0+Nk-1,i0+Ni-1)+=arma::trans(Kik);
343 }
344
345 // Then, the second part: K(j,k) += (ij|kl) P(i,l)
346 if(is!=js) {
347 arma::cx_mat Kjk(Nj,Nk);
348 Kjk.zeros();
349 arma::cx_mat Pil=P.submat(i0,l0,i0+Ni-1,l0+Nl-1);
350
351 // Increment Kjk
352 for(size_t jj=0;jj<Nj;jj++)
353 for(size_t kk=0;kk<Nk;kk++)
354 for(size_t ll=0;ll<Nl;ll++)
355 for(size_t ii=0;ii<Ni;ii++)
356 Kjk(jj,kk)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pil(ii,ll);
357
358 // Set elements
359 K.submat(j0,k0,j0+Nj-1,k0+Nk-1)+=Kjk;
360 // Symmetrize if necessary (take care about possible overlap with next routine)
361 if(ip!=jp) {
362 K.submat(k0,j0,k0+Nk-1,j0+Nj-1)+=arma::trans(Kjk);
363 }
364 }
365
366 // Third part: K(i,l) += (ij|kl) P(j,k)
367 if(ks!=ls) {
368 arma::cx_mat Kil(Ni,Nl);
369 Kil.zeros();
370 arma::cx_mat Pjk=P.submat(j0,k0,j0+Nj-1,k0+Nk-1);
371
372 // Increment Kil
373 for(size_t ii=0;ii<Ni;ii++)
374 for(size_t ll=0;ll<Nl;ll++)
375 for(size_t jj=0;jj<Nj;jj++)
376 for(size_t kk=0;kk<Nk;kk++) {
377 Kil(ii,ll)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pjk(jj,kk);
378 }
379
380 // Set elements
381 K.submat(i0,l0,i0+Ni-1,l0+Nl-1)+=Kil;
382 // Symmetrize if necessary
383 if(ip!=jp)
384 K.submat(l0,i0,l0+Nl-1,i0+Ni-1)+=arma::trans(Kil);
385 }
386
387 // Last permutation: K(j,l) += (ij|kl) P(i,k)
388 if(is!=js && ks!=ls) {
389 arma::cx_mat Kjl(Nj,Nl);
390 Kjl.zeros();
391 arma::cx_mat Pik=P.submat(i0,k0,i0+Ni-1,k0+Nk-1);
392
393 // Increment Kjl
394 for(size_t jj=0;jj<Nj;jj++)
395 for(size_t ll=0;ll<Nl;ll++)
396 for(size_t ii=0;ii<Ni;ii++)
397 for(size_t kk=0;kk<Nk;kk++) {
398 Kjl(jj,ll)+=ints[ioff+((ii*Nj+jj)*Nk+kk)*Nl+ll]*Pik(ii,kk);
399 }
400
401 // Set elements
402 K.submat(j0,l0,j0+Nj-1,l0+Nl-1)+=Kjl;
403 // Symmetrize if necessary
404 if (ip!=jp)
405 K.submat(l0,j0,l0+Nl-1,j0+Nj-1)+=arma::trans(Kjl);
406 }
407 }
408
get_K() const409 arma::cx_mat cxKDigestor::get_K() const {
410 return K;
411 }
412
ForceDigestor()413 ForceDigestor::ForceDigestor() {
414 }
415
~ForceDigestor()416 ForceDigestor::~ForceDigestor() {
417 }
418
JFDigestor(const arma::mat & P_)419 JFDigestor::JFDigestor(const arma::mat & P_) : P(P_) {
420 }
421
~JFDigestor()422 JFDigestor::~JFDigestor() {
423 }
424
digest(const std::vector<eripair_t> & shpairs,size_t ip,size_t jp,dERIWorker & deriw,arma::vec & f)425 void JFDigestor::digest(const std::vector<eripair_t> & shpairs, size_t ip, size_t jp, dERIWorker & deriw, arma::vec & f) {
426 // Shells in question are
427 size_t is=shpairs[ip].is;
428 size_t js=shpairs[ip].js;
429 size_t ks=shpairs[jp].is;
430 size_t ls=shpairs[jp].js;
431
432 // Amount of functions on the first pair is
433 size_t Ni=shpairs[ip].Ni;
434 size_t Nj=shpairs[ip].Nj;
435 // and on second pair is
436 size_t Nk=shpairs[jp].Ni;
437 size_t Nl=shpairs[jp].Nj;
438
439 // First functions on the first pair is
440 size_t i0=shpairs[ip].i0;
441 size_t j0=shpairs[ip].j0;
442 // Second pair is
443 size_t k0=shpairs[jp].i0;
444 size_t l0=shpairs[jp].j0;
445
446 // E_J = P_ij (ij|kl) P_kl. Work matrices
447 arma::mat Pij=P.submat(i0,j0,i0+Ni-1,j0+Nj-1);
448 arma::mat Pkl=P.submat(k0,l0,k0+Nk-1,l0+Nl-1);
449
450 // Degeneracy factor
451 double Jfac=-0.5;
452 if(is!=js)
453 Jfac*=2.0;
454 if(ks!=ls)
455 Jfac*=2.0;
456 if(ip!=jp)
457 Jfac*=2.0;
458
459 // Increment the forces.
460 for(int idx=0;idx<12;idx++) {
461 // Get the integral derivatives
462 const std::vector<double> *erip=deriw.getp(idx);
463
464 // E_J = P_ij (ij|kl) P_kl
465 double el=0.0;
466 for(size_t ii=0;ii<Ni;ii++)
467 for(size_t jj=0;jj<Nj;jj++)
468 for(size_t kk=0;kk<Nk;kk++)
469 for(size_t ll=0;ll<Nl;ll++)
470 el+=Pij(ii,jj)*Pkl(kk,ll)*(*erip)[((ii*Nj+jj)*Nk+kk)*Nl+ll];
471
472 // Increment the element
473 f(idx)+=Jfac*el;
474 }
475 }
476
KFDigestor(const arma::mat & P_,double kfrac_,bool restr)477 KFDigestor::KFDigestor(const arma::mat & P_, double kfrac_, bool restr) : P(P_), kfrac(kfrac_) {
478 fac = restr ? 0.5 : 1.0;
479 }
480
~KFDigestor()481 KFDigestor::~KFDigestor() {
482 }
483
digest(const std::vector<eripair_t> & shpairs,size_t ip,size_t jp,dERIWorker & deriw,arma::vec & f)484 void KFDigestor::digest(const std::vector<eripair_t> & shpairs, size_t ip, size_t jp, dERIWorker & deriw, arma::vec & f) {
485 // Shells on quartet are
486 size_t is=shpairs[ip].is;
487 size_t js=shpairs[ip].js;
488 size_t ks=shpairs[jp].is;
489 size_t ls=shpairs[jp].js;
490
491 // Amount of functions on the first pair is
492 size_t Ni=shpairs[ip].Ni;
493 size_t Nj=shpairs[ip].Nj;
494 // and on second pair is
495 size_t Nk=shpairs[jp].Ni;
496 size_t Nl=shpairs[jp].Nj;
497
498 // First functions on the first pair is
499 size_t i0=shpairs[ip].i0;
500 size_t j0=shpairs[ip].j0;
501 // Second pair is
502 size_t k0=shpairs[jp].i0;
503 size_t l0=shpairs[jp].j0;
504
505 // E_K = P_ik (ij|kl) P_jl
506 arma::mat Pik=P.submat(i0,k0,i0+Ni-1,k0+Nk-1);
507 arma::mat Pjl=P.submat(j0,l0,j0+Nj-1,l0+Nl-1);
508 // + P_jk (ij|kl) P_il
509 arma::mat Pjk=P.submat(j0,k0,j0+Nj-1,k0+Nk-1);
510 arma::mat Pil=P.submat(i0,l0,i0+Ni-1,l0+Nl-1);
511 double K1fac, K2fac;
512 if(is!=js && ks!=ls) {
513 // Get both twice.
514 K1fac=1.0;
515 K2fac=1.0;
516 } else if(is==js && ks==ls) {
517 // Only get the first one, once.
518 K1fac=0.5;
519 K2fac=0.0;
520 } else {
521 // Get both once.
522 K1fac=0.5;
523 K2fac=0.5;
524 }
525 // Switch symmetry
526 if(ip!=jp) {
527 K1fac*=2.0;
528 K2fac*=2.0;
529 }
530 // Restricted calculation?
531 K1fac*=fac*kfrac;
532 K2fac*=fac*kfrac;
533
534 // Increment the forces.
535 for(int idx=0;idx<12;idx++) {
536 // Get the integral derivatives
537 const std::vector<double> * erip=deriw.getp(idx);
538
539 // E_K = P_ik (ij|kl) P_jl
540 double el=0.0;
541 // Increment matrix
542 for(size_t ii=0;ii<Ni;ii++)
543 for(size_t jj=0;jj<Nj;jj++)
544 for(size_t kk=0;kk<Nk;kk++)
545 for(size_t ll=0;ll<Nl;ll++)
546 el+=Pik(ii,kk)*Pjl(jj,ll)*(*erip)[((ii*Nj+jj)*Nk+kk)*Nl+ll];
547
548 // Increment the element
549 f(idx)+=K1fac*el;
550
551 // Second contribution
552 if(K2fac!=0.0) {
553 el=0.0;
554 // Increment matrix
555 for(size_t ii=0;ii<Ni;ii++)
556 for(size_t jj=0;jj<Nj;jj++)
557 for(size_t kk=0;kk<Nk;kk++)
558 for(size_t ll=0;ll<Nl;ll++)
559 el+=Pjk(jj,kk)*Pil(ii,ll)*(*erip)[((ii*Nj+jj)*Nk+kk)*Nl+ll];
560
561 // Increment the element
562 f(idx)+=K2fac*el;
563 }
564 }
565 }
566