1
2 /*
3 * This source code is part of
4 *
5 * E R K A L E
6 * -
7 * DFT from Hel
8 *
9 * Written by Susi Lehtola, 2010-2011
10 * Copyright (c) 2010-2011, Susi Lehtola
11 *
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
16 */
17
18 #include "eriscreen.h"
19 #include "eri_digest.h"
20 #include "eriworker.h"
21 #include "mathf.h"
22 #include "timer.h"
23
24 #include <cstdio>
25 // For exceptions
26 #include <sstream>
27 #include <stdexcept>
28
29 #ifdef _OPENMP
30 #include <omp.h>
31 #endif
32
33 // Print out screening table?
34 //#define PRINTOUT
35
36
37 #ifdef CHECKFILL
idx(size_t i,size_t j,size_t k,size_t l)38 size_t idx(size_t i, size_t j, size_t k, size_t l) {
39 if(i<j)
40 std::swap(i,j);
41 if(k<l)
42 std::swap(k,l);
43
44 size_t ij=(i*(i+1))/2+j;
45 size_t kl=(k*(k+1))/2+l;
46
47 if(ij<kl)
48 std::swap(ij,kl);
49
50 return (ij*(ij+1))/2+kl;
51 }
52 #endif
53
54
ERIscreen()55 ERIscreen::ERIscreen() {
56 omega=0.0;
57 alpha=1.0;
58 beta=0.0;
59 basp=NULL;
60 }
61
~ERIscreen()62 ERIscreen::~ERIscreen() {
63 }
64
get_N() const65 size_t ERIscreen::get_N() const {
66 return Nbf;
67 }
68
set_range_separation(double w,double a,double b)69 void ERIscreen::set_range_separation(double w, double a, double b) {
70 omega=w;
71 alpha=a;
72 beta=b;
73 }
74
get_range_separation(double & w,double & a,double & b) const75 void ERIscreen::get_range_separation(double & w, double & a, double & b) const {
76 w=omega;
77 a=alpha;
78 b=beta;
79 }
80
fill(const BasisSet * basisv,double shtol,bool verbose)81 size_t ERIscreen::fill(const BasisSet * basisv, double shtol, bool verbose) {
82 // Form screening table.
83 if(basisv==NULL)
84 return 0;
85
86 basp=basisv;
87 Nbf=basisv->get_Nbf();
88
89 // Form index helper
90 iidx=i_idx(Nbf);
91
92 // Shell-pair list
93 shpairs=basp->get_eripairs(Q,M,shtol,omega,alpha,beta,verbose);
94
95 return shpairs.size();
96 }
97
calculate(std::vector<std::vector<IntegralDigestor * >> & digest,double tol) const98 void ERIscreen::calculate(std::vector< std::vector<IntegralDigestor *> > & digest, double tol) const {
99 // Shells in basis set
100 std::vector<GaussianShell> shells=basp->get_shells();
101 // Get number of shell pairs
102 const size_t Npairs=shpairs.size();
103
104 #ifdef _OPENMP
105 #pragma omp parallel
106 #endif
107 {
108 // ERI worker
109 ERIWorker *eri = (omega==0.0 && alpha==1.0 && beta==0.0) ? new ERIWorker(basp->get_max_am(),basp->get_max_Ncontr()) : new ERIWorker_srlr(basp->get_max_am(),basp->get_max_Ncontr(),omega,alpha,beta);
110
111 // Integral array
112 const std::vector<double> * erip;
113
114 #ifndef _OPENMP
115 int ith=0;
116 #else
117 int ith(omp_get_thread_num());
118 #pragma omp for schedule(dynamic)
119 #endif
120 for(size_t ip=0;ip<Npairs;ip++) {
121 // Loop over second pairs
122 for(size_t jp=0;jp<=ip;jp++) {
123 // Shells on first pair
124 size_t is=shpairs[ip].is;
125 size_t js=shpairs[ip].js;
126 // and those on the second pair
127 size_t ks=shpairs[jp].is;
128 size_t ls=shpairs[jp].js;
129
130 // Schwarz screening estimate
131 double QQ=Q(is,js)*Q(ks,ls);
132 if(QQ<tol) {
133 // Skip due to small value of integral. Because the
134 // integrals have been ordered wrt Q, all the next ones
135 // will be small as well!
136 break;
137 }
138
139 // Distance screening estimate
140 double MM1=M(is,ks)*M(js,ls);
141 double MM2=M(is,ls)*M(js,ks);
142 if(MM1<tol || MM2<tol) {
143 // This pair is negligible
144 continue;
145 }
146
147 // Compute integrals
148 eri->compute(&shells[is],&shells[js],&shells[ks],&shells[ls]);
149 erip=eri->getp();
150
151 // Digest the integrals
152 for(size_t i=0;i<digest[ith].size();i++)
153 digest[ith][i]->digest(shpairs,ip,jp,*erip,0);
154 }
155 }
156
157 delete eri;
158 }
159 }
160
calculate_force(std::vector<std::vector<ForceDigestor * >> & digest,double tol) const161 arma::vec ERIscreen::calculate_force(std::vector< std::vector<ForceDigestor *> > & digest, double tol) const {
162 // Shells
163 std::vector<GaussianShell> shells=basp->get_shells();
164 // Get number of shell pairs
165 const size_t Npairs=shpairs.size();
166
167 // Forces
168 arma::vec F(3*basp->get_Nnuc());
169 F.zeros();
170
171 #ifdef _OPENMP
172 #pragma omp parallel
173 #endif
174 {
175 /// ERI derivative worker
176 dERIWorker *deri = (omega==0.0 && alpha==1.0 && beta==0.0) ? new dERIWorker(basp->get_max_am(),basp->get_max_Ncontr()) : new dERIWorker_srlr(basp->get_max_am(),basp->get_max_Ncontr(),omega,alpha,beta);
177
178 // Shell centers
179 size_t inuc, jnuc, knuc, lnuc;
180
181 #ifndef _OPENMP
182 int ith=0;
183 #else
184 int ith(omp_get_thread_num());
185 arma::vec Fwrk(F);
186 #pragma omp for schedule(dynamic)
187 #endif
188 for(size_t ip=0;ip<Npairs;ip++) {
189 for(size_t jp=0;jp<=ip;jp++) {
190 // Shells on first pair
191 size_t is=shpairs[ip].is;
192 size_t js=shpairs[ip].js;
193 // and those on the second pair
194 size_t ks=shpairs[jp].is;
195 size_t ls=shpairs[jp].js;
196
197 // Shell centers
198 inuc=shells[is].get_center_ind();
199 jnuc=shells[js].get_center_ind();
200 knuc=shells[ks].get_center_ind();
201 lnuc=shells[ls].get_center_ind();
202
203 // Skip when all functions are on the same nucleus - force will vanish
204 if(inuc==jnuc && jnuc==knuc && knuc==lnuc)
205 continue;
206
207 // Schwarz screening estimate
208 double QQ=Q(is,js)*Q(ks,ls);
209 if(QQ<tol) {
210 // Skip due to small value of integral. Because the
211 // integrals have been ordered wrt Q, all the next ones
212 // will be small as well!
213 break;
214 }
215
216 // Distance screening estimate
217 double MM1=M(is,ks)*M(js,ls);
218 double MM2=M(is,ls)*M(js,ks);
219 if(MM1<tol || MM2<tol) {
220 // This pair is negligible
221 continue;
222 }
223
224 // Compute the derivatives.
225 deri->compute(&shells[is],&shells[js],&shells[ks],&shells[ls]);
226
227 // Digest the forces on the nuclei
228 arma::vec f(12);
229 f.zeros();
230
231 // Digest the integrals
232 for(size_t i=0;i<digest[ith].size();i++)
233 digest[ith][i]->digest(shpairs,ip,jp,*deri,f);
234
235 // Increment forces
236 #ifdef _OPENMP
237 Fwrk.subvec(3*inuc,3*inuc+2)+=f.subvec(0,2);
238 Fwrk.subvec(3*jnuc,3*jnuc+2)+=f.subvec(3,5);
239 Fwrk.subvec(3*knuc,3*knuc+2)+=f.subvec(6,8);
240 Fwrk.subvec(3*lnuc,3*lnuc+2)+=f.subvec(9,11);
241 #else
242 F.subvec(3*inuc,3*inuc+2)+=f.subvec(0,2);
243 F.subvec(3*jnuc,3*jnuc+2)+=f.subvec(3,5);
244 F.subvec(3*knuc,3*knuc+2)+=f.subvec(6,8);
245 F.subvec(3*lnuc,3*lnuc+2)+=f.subvec(9,11);
246 #endif
247 }
248 }
249
250 // Collect results
251 #ifdef _OPENMP
252 #pragma omp critical
253 F+=Fwrk;
254 #endif
255
256 delete deri;
257 }
258
259 return F;
260 }
261
calcJ(const arma::mat & P,double tol) const262 arma::mat ERIscreen::calcJ(const arma::mat & P, double tol) const {
263 if(P.n_rows != Nbf || P.n_cols != Nbf) {
264 std::ostringstream oss;
265 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
266 throw std::logic_error(oss.str());
267 }
268
269 #ifdef _OPENMP
270 int nth=omp_get_max_threads();
271 #else
272 int nth=1;
273 #endif
274
275 // Get workers
276 std::vector< std::vector<IntegralDigestor *> > p(nth);
277 #ifdef _OPENMP
278 #pragma omp parallel for
279 #endif
280 for(int i=0;i<nth;i++) {
281 p[i].resize(1);
282 p[i][0]=new JDigestor(P);
283 }
284
285 // Do calculation
286 calculate(p,tol);
287
288 // Collect results
289 arma::mat J(((JDigestor *) p[0][0])->get_J());
290 for(int i=1;i<nth;i++)
291 J+=((JDigestor *) p[i][0])->get_J();
292
293 // Free memory
294 for(size_t i=0;i<p.size();i++)
295 for(size_t j=0;j<p[i].size();j++)
296 delete p[i][j];
297
298 return J;
299 }
300
calcK(const arma::mat & P,double tol) const301 arma::mat ERIscreen::calcK(const arma::mat & P, double tol) const {
302 if(P.n_rows != Nbf || P.n_cols != Nbf) {
303 std::ostringstream oss;
304 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
305 throw std::logic_error(oss.str());
306 }
307
308 #ifdef _OPENMP
309 int nth=omp_get_max_threads();
310 #else
311 int nth=1;
312 #endif
313
314 // Get workers
315 std::vector< std::vector<IntegralDigestor *> > p(nth);
316 #ifdef _OPENMP
317 #pragma omp parallel for
318 #endif
319 for(int i=0;i<nth;i++) {
320 p[i].resize(1);
321 p[i][0]=new KDigestor(P);
322 }
323
324 // Do calculation
325 calculate(p,tol);
326
327 // Collect results
328 arma::mat K(((KDigestor *) p[0][0])->get_K());
329 for(int i=1;i<nth;i++)
330 K+=((KDigestor *) p[i][0])->get_K();
331
332 // Free memory
333 for(size_t i=0;i<p.size();i++)
334 for(size_t j=0;j<p[i].size();j++)
335 delete p[i][j];
336
337 return K;
338 }
339
calcK(const arma::cx_mat & P,double tol) const340 arma::cx_mat ERIscreen::calcK(const arma::cx_mat & P, double tol) const {
341 if(P.n_rows != Nbf || P.n_cols != Nbf) {
342 std::ostringstream oss;
343 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
344 throw std::logic_error(oss.str());
345 }
346
347 #ifdef _OPENMP
348 int nth=omp_get_max_threads();
349 #else
350 int nth=1;
351 #endif
352
353 // Get workers
354 std::vector< std::vector<IntegralDigestor *> > p(nth);
355 #ifdef _OPENMP
356 #pragma omp parallel for
357 #endif
358 for(int i=0;i<nth;i++) {
359 p[i].resize(1);
360 p[i][0]=new cxKDigestor(P);
361 }
362
363 // Do calculation
364 calculate(p,tol);
365
366 // Collect results
367 arma::cx_mat K(((cxKDigestor *) p[0][0])->get_K());
368 for(int i=1;i<nth;i++)
369 K+=((cxKDigestor *) p[i][0])->get_K();
370
371 // Free memory
372 for(size_t i=0;i<p.size();i++)
373 for(size_t j=0;j<p[i].size();j++)
374 delete p[i][j];
375
376 return K;
377 }
378
calcK(const arma::mat & Pa,const arma::mat & Pb,arma::mat & Ka,arma::mat & Kb,double tol) const379 void ERIscreen::calcK(const arma::mat & Pa, const arma::mat & Pb, arma::mat & Ka, arma::mat & Kb, double tol) const {
380 if(Pa.n_rows != Nbf || Pa.n_cols != Nbf) {
381 std::ostringstream oss;
382 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pa.n_rows = " << Pa.n_rows << ", Pa.n_cols = " << Pa.n_cols << "!\n";
383 throw std::logic_error(oss.str());
384 }
385 if(Pb.n_rows != Nbf || Pb.n_cols != Nbf) {
386 std::ostringstream oss;
387 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pb.n_rows = " << Pb.n_rows << ", Pb.n_cols = " << Pb.n_cols << "!\n";
388 throw std::logic_error(oss.str());
389 }
390
391 #ifdef _OPENMP
392 int nth=omp_get_max_threads();
393 #else
394 int nth=1;
395 #endif
396
397 // Get workers
398 std::vector< std::vector<IntegralDigestor *> > p(nth);
399 #ifdef _OPENMP
400 #pragma omp parallel for
401 #endif
402 for(int i=0;i<nth;i++) {
403 p[i].resize(2);
404 p[i][0]=new KDigestor(Pa);
405 p[i][1]=new KDigestor(Pb);
406 }
407
408 // Do calculation
409 calculate(p,tol);
410
411 // Collect results
412 Ka=((KDigestor *) p[0][0])->get_K();
413 Kb=((KDigestor *) p[0][1])->get_K();
414 for(int i=1;i<nth;i++) {
415 Ka+=((KDigestor *) p[i][0])->get_K();
416 Kb+=((KDigestor *) p[i][1])->get_K();
417 }
418
419 // Free memory
420 for(size_t i=0;i<p.size();i++)
421 for(size_t j=0;j<p[i].size();j++)
422 delete p[i][j];
423 }
424
calcK(const arma::cx_mat & Pa,const arma::cx_mat & Pb,arma::cx_mat & Ka,arma::cx_mat & Kb,double tol) const425 void ERIscreen::calcK(const arma::cx_mat & Pa, const arma::cx_mat & Pb, arma::cx_mat & Ka, arma::cx_mat & Kb, double tol) const {
426 if(Pa.n_rows != Nbf || Pa.n_cols != Nbf) {
427 std::ostringstream oss;
428 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pa.n_rows = " << Pa.n_rows << ", Pa.n_cols = " << Pa.n_cols << "!\n";
429 throw std::logic_error(oss.str());
430 }
431 if(Pb.n_rows != Nbf || Pb.n_cols != Nbf) {
432 std::ostringstream oss;
433 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pb.n_rows = " << Pb.n_rows << ", Pb.n_cols = " << Pb.n_cols << "!\n";
434 throw std::logic_error(oss.str());
435 }
436
437 #ifdef _OPENMP
438 int nth=omp_get_max_threads();
439 #else
440 int nth=1;
441 #endif
442
443 // Get workers
444 std::vector< std::vector<IntegralDigestor *> > p(nth);
445 #ifdef _OPENMP
446 #pragma omp parallel for
447 #endif
448 for(int i=0;i<nth;i++) {
449 p[i].resize(2);
450 p[i][0]=new cxKDigestor(Pa);
451 p[i][1]=new cxKDigestor(Pb);
452 }
453
454 // Do calculation
455 calculate(p,tol);
456
457 // Collect results
458 Ka=((cxKDigestor *) p[0][0])->get_K();
459 Kb=((cxKDigestor *) p[0][1])->get_K();
460 for(int i=1;i<nth;i++) {
461 Ka+=((cxKDigestor *) p[i][0])->get_K();
462 Kb+=((cxKDigestor *) p[i][1])->get_K();
463 }
464
465 // Free memory
466 for(size_t i=0;i<p.size();i++)
467 for(size_t j=0;j<p[i].size();j++)
468 delete p[i][j];
469 }
470
calcJK(const arma::mat & P,arma::mat & J,arma::mat & K,double tol) const471 void ERIscreen::calcJK(const arma::mat & P, arma::mat & J, arma::mat & K, double tol) const {
472 if(P.n_rows != Nbf || P.n_cols != Nbf) {
473 std::ostringstream oss;
474 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
475 throw std::logic_error(oss.str());
476 }
477
478 #ifdef _OPENMP
479 int nth=omp_get_max_threads();
480 #else
481 int nth=1;
482 #endif
483
484 // Get workers
485 std::vector< std::vector<IntegralDigestor *> > p(nth);
486 #ifdef _OPENMP
487 #pragma omp parallel for
488 #endif
489 for(int i=0;i<nth;i++) {
490 p[i].resize(2);
491 p[i][0]=new JDigestor(P);
492 p[i][1]=new KDigestor(P);
493 }
494
495 // Do calculation
496 calculate(p,tol);
497
498 // Collect results
499 J=((JDigestor *) p[0][0])->get_J();
500 K=((KDigestor *) p[0][1])->get_K();
501 for(int i=1;i<nth;i++) {
502 J+=((JDigestor *) p[i][0])->get_J();
503 K+=((KDigestor *) p[i][1])->get_K();
504 }
505
506 // Free memory
507 for(size_t i=0;i<p.size();i++)
508 for(size_t j=0;j<p[i].size();j++)
509 delete p[i][j];
510 }
511
calcJK(const arma::cx_mat & P,arma::mat & J,arma::cx_mat & K,double tol) const512 void ERIscreen::calcJK(const arma::cx_mat & P, arma::mat & J, arma::cx_mat & K, double tol) const {
513 if(P.n_rows != Nbf || P.n_cols != Nbf) {
514 std::ostringstream oss;
515 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
516 throw std::logic_error(oss.str());
517 }
518
519 #ifdef _OPENMP
520 int nth=omp_get_max_threads();
521 #else
522 int nth=1;
523 #endif
524
525 // Get workers
526 std::vector< std::vector<IntegralDigestor *> > p(nth);
527 #ifdef _OPENMP
528 #pragma omp parallel for
529 #endif
530 for(int i=0;i<nth;i++) {
531 p[i].resize(2);
532 p[i][0]=new JDigestor(arma::real(P));
533 p[i][1]=new cxKDigestor(P);
534 }
535
536 // Do calculation
537 calculate(p,tol);
538
539 // Collect results
540 J=((JDigestor *) p[0][0])->get_J();
541 K=((cxKDigestor *) p[0][1])->get_K();
542 for(int i=1;i<nth;i++) {
543 J+=((JDigestor *) p[i][0])->get_J();
544 K+=((cxKDigestor *) p[i][1])->get_K();
545 }
546
547 // Free memory
548 for(size_t i=0;i<p.size();i++)
549 for(size_t j=0;j<p[i].size();j++)
550 delete p[i][j];
551 }
552
calcJK(const arma::mat & Pa,const arma::mat & Pb,arma::mat & J,arma::mat & Ka,arma::mat & Kb,double tol) const553 void ERIscreen::calcJK(const arma::mat & Pa, const arma::mat & Pb, arma::mat & J, arma::mat & Ka, arma::mat & Kb, double tol) const {
554 if(Pa.n_rows != Nbf || Pa.n_cols != Nbf) {
555 std::ostringstream oss;
556 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pa.n_rows = " << Pa.n_rows << ", Pa.n_cols = " << Pa.n_cols << "!\n";
557 throw std::logic_error(oss.str());
558 }
559 if(Pb.n_rows != Nbf || Pb.n_cols != Nbf) {
560 std::ostringstream oss;
561 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pb.n_rows = " << Pb.n_rows << ", Pb.n_cols = " << Pb.n_cols << "!\n";
562 throw std::logic_error(oss.str());
563 }
564
565 #ifdef _OPENMP
566 int nth=omp_get_max_threads();
567 #else
568 int nth=1;
569 #endif
570
571 // Get workers
572 std::vector< std::vector<IntegralDigestor *> > p(nth);
573 #ifdef _OPENMP
574 #pragma omp parallel for
575 #endif
576 for(int i=0;i<nth;i++) {
577 p[i].resize(3);
578 p[i][0]=new JDigestor(Pa+Pb);
579 p[i][1]=new KDigestor(Pa);
580 p[i][2]=new KDigestor(Pb);
581 }
582
583 // Do calculation
584 calculate(p,tol);
585
586 // Collect results
587 J=((JDigestor *) p[0][0])->get_J();
588 Ka=((KDigestor *) p[0][1])->get_K();
589 Kb=((KDigestor *) p[0][2])->get_K();
590 for(int i=1;i<nth;i++) {
591 J+=((JDigestor *) p[i][0])->get_J();
592 Ka+=((KDigestor *) p[i][1])->get_K();
593 Kb+=((KDigestor *) p[i][2])->get_K();
594 }
595
596 // Free memory
597 for(size_t i=0;i<p.size();i++)
598 for(size_t j=0;j<p[i].size();j++)
599 delete p[i][j];
600 }
601
calcJK(const arma::cx_mat & Pa,const arma::cx_mat & Pb,arma::mat & J,arma::cx_mat & Ka,arma::cx_mat & Kb,double tol) const602 void ERIscreen::calcJK(const arma::cx_mat & Pa, const arma::cx_mat & Pb, arma::mat & J, arma::cx_mat & Ka, arma::cx_mat & Kb, double tol) const {
603 if(Pa.n_rows != Nbf || Pa.n_cols != Nbf) {
604 std::ostringstream oss;
605 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pa.n_rows = " << Pa.n_rows << ", Pa.n_cols = " << Pa.n_cols << "!\n";
606 throw std::logic_error(oss.str());
607 }
608 if(Pb.n_rows != Nbf || Pb.n_cols != Nbf) {
609 std::ostringstream oss;
610 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pb.n_rows = " << Pb.n_rows << ", Pb.n_cols = " << Pb.n_cols << "!\n";
611 throw std::logic_error(oss.str());
612 }
613
614 #ifdef _OPENMP
615 int nth=omp_get_max_threads();
616 #else
617 int nth=1;
618 #endif
619
620 // Get workers
621 std::vector< std::vector<IntegralDigestor *> > p(nth);
622 #ifdef _OPENMP
623 #pragma omp parallel for
624 #endif
625 for(int i=0;i<nth;i++) {
626 p[i].resize(3);
627 p[i][0]=new JDigestor(arma::real(Pa+Pb));
628 p[i][1]=new cxKDigestor(Pa);
629 p[i][2]=new cxKDigestor(Pb);
630 }
631
632 // Do calculation
633 calculate(p,tol);
634
635 // Collect results
636 J=((JDigestor *) p[0][0])->get_J();
637 Ka=((cxKDigestor *) p[0][1])->get_K();
638 Kb=((cxKDigestor *) p[0][2])->get_K();
639 for(int i=1;i<nth;i++) {
640 J+=((JDigestor *) p[i][0])->get_J();
641 Ka+=((cxKDigestor *) p[i][1])->get_K();
642 Kb+=((cxKDigestor *) p[i][2])->get_K();
643 }
644
645 // Free memory
646 for(size_t i=0;i<p.size();i++)
647 for(size_t j=0;j<p[i].size();j++)
648 delete p[i][j];
649 }
650
calcJK(const std::vector<arma::cx_mat> & P,double jfrac,double kfrac,double tol) const651 std::vector<arma::cx_mat> ERIscreen::calcJK(const std::vector<arma::cx_mat> & P, double jfrac, double kfrac, double tol) const {
652 for(size_t i=0;i<P.size();i++) {
653 if(P[i].n_rows != Nbf || P[i].n_cols != Nbf) {
654 std::ostringstream oss;
655 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P[" << i << "].n_rows = " << P[i].n_rows << ", P[" << i << "].n_cols = " << P[i].n_cols << "!\n";
656 throw std::logic_error(oss.str());
657 }
658 }
659
660 #ifdef _OPENMP
661 int nth=omp_get_max_threads();
662 #else
663 int nth=1;
664 #endif
665
666 bool doj(jfrac!=0.0);
667 bool dok(kfrac!=0.0);
668
669 // Get workers
670 std::vector< std::vector<IntegralDigestor *> > p(nth);
671 #ifdef _OPENMP
672 #pragma omp parallel for
673 #endif
674 for(int i=0;i<nth;i++) {
675 if(doj) {
676 for(size_t j=0;j<P.size();j++)
677 p[i].push_back(new JDigestor(arma::real(P[j])));
678 }
679 if(dok) {
680 for(size_t j=0;j<P.size();j++)
681 p[i].push_back(new cxKDigestor(P[j]));
682 }
683 }
684
685 // Do calculation
686 calculate(p,tol);
687
688 // Collect results
689 std::vector<arma::cx_mat> JK(P.size());
690 for(size_t i=0;i<JK.size();i++)
691 JK[i].zeros(P[i].n_rows,P[i].n_cols);
692
693 size_t joff=0;
694 if(doj) {
695 for(size_t j=0;j<P.size();j++)
696 for(int i=0;i<nth;i++)
697 JK[j]+=jfrac*((JDigestor *) p[i][j+joff])->get_J()*COMPLEX1;
698 joff+=P.size();
699 }
700 // Exchange contribution
701 if(dok) {
702 for(size_t j=0;j<P.size();j++)
703 for(int i=0;i<nth;i++)
704 JK[j]-=kfrac*((cxKDigestor *) p[i][j+joff])->get_K();
705 joff+=P.size();
706 }
707
708 // Free memory
709 for(size_t i=0;i<p.size();i++)
710 for(size_t j=0;j<p[i].size();j++)
711 delete p[i][j];
712
713 return JK;
714 }
715
forceJ(const arma::mat & P,double tol) const716 arma::vec ERIscreen::forceJ(const arma::mat & P, double tol) const {
717 if(P.n_rows != Nbf || P.n_cols != Nbf) {
718 std::ostringstream oss;
719 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
720 throw std::logic_error(oss.str());
721 }
722
723 #ifdef _OPENMP
724 int nth=omp_get_max_threads();
725 #else
726 int nth=1;
727 #endif
728
729 // Get workers
730 std::vector< std::vector<ForceDigestor *> > p(nth);
731 #ifdef _OPENMP
732 #pragma omp parallel for
733 #endif
734 for(int i=0;i<nth;i++) {
735 p[i].resize(1);
736 p[i][0]=new JFDigestor(P);
737 }
738
739 // Do calculation
740 arma::vec f=calculate_force(p,tol);
741
742 // Free memory
743 for(size_t i=0;i<p.size();i++)
744 for(size_t j=0;j<p[i].size();j++)
745 delete p[i][j];
746
747 return f;
748 }
749
forceK(const arma::mat & P,double tol,double kfrac) const750 arma::vec ERIscreen::forceK(const arma::mat & P, double tol, double kfrac) const {
751 if(P.n_rows != Nbf || P.n_cols != Nbf) {
752 std::ostringstream oss;
753 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
754 throw std::logic_error(oss.str());
755 }
756
757 #ifdef _OPENMP
758 int nth=omp_get_max_threads();
759 #else
760 int nth=1;
761 #endif
762
763 // Get workers
764 std::vector< std::vector<ForceDigestor *> > p(nth);
765 #ifdef _OPENMP
766 #pragma omp parallel for
767 #endif
768 for(int i=0;i<nth;i++) {
769 p[i].resize(1);
770 p[i][0]=new KFDigestor(P,kfrac,true);
771 }
772
773 // Do calculation
774 arma::vec f=calculate_force(p,tol);
775
776 // Free memory
777 for(size_t i=0;i<p.size();i++)
778 for(size_t j=0;j<p[i].size();j++)
779 delete p[i][j];
780
781 return f;
782 }
783
forceK(const arma::mat & Pa,const arma::mat & Pb,double tol,double kfrac) const784 arma::vec ERIscreen::forceK(const arma::mat & Pa, const arma::mat & Pb, double tol, double kfrac) const {
785 if(Pa.n_rows != Nbf || Pa.n_cols != Nbf) {
786 std::ostringstream oss;
787 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pa.n_rows = " << Pa.n_rows << ", Pa.n_cols = " << Pa.n_cols << "!\n";
788 throw std::logic_error(oss.str());
789 }
790 if(Pb.n_rows != Nbf || Pb.n_cols != Nbf) {
791 std::ostringstream oss;
792 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pb.n_rows = " << Pb.n_rows << ", Pb.n_cols = " << Pb.n_cols << "!\n";
793 throw std::logic_error(oss.str());
794 }
795
796 #ifdef _OPENMP
797 int nth=omp_get_max_threads();
798 #else
799 int nth=1;
800 #endif
801
802 // Get workers
803 std::vector< std::vector<ForceDigestor *> > p(nth);
804 #ifdef _OPENMP
805 #pragma omp parallel for
806 #endif
807 for(int i=0;i<nth;i++) {
808 p[i].resize(3);
809 p[i][0]=new JFDigestor(Pa+Pb);
810 p[i][1]=new KFDigestor(Pa,kfrac,false);
811 p[i][2]=new KFDigestor(Pb,kfrac,false);
812 }
813
814 // Do calculation
815 arma::vec f=calculate_force(p,tol);
816
817 // Free memory
818 for(size_t i=0;i<p.size();i++)
819 for(size_t j=0;j<p[i].size();j++)
820 delete p[i][j];
821
822 return f;
823 }
824
forceJK(const arma::mat & P,double tol,double kfrac) const825 arma::vec ERIscreen::forceJK(const arma::mat & P, double tol, double kfrac) const {
826 if(P.n_rows != Nbf || P.n_cols != Nbf) {
827 std::ostringstream oss;
828 oss << "Error in ERIscreen: Nbf = " << Nbf << ", P.n_rows = " << P.n_rows << ", P.n_cols = " << P.n_cols << "!\n";
829 throw std::logic_error(oss.str());
830 }
831
832 #ifdef _OPENMP
833 int nth=omp_get_max_threads();
834 #else
835 int nth=1;
836 #endif
837
838 // Get workers
839 std::vector< std::vector<ForceDigestor *> > p(nth);
840 #ifdef _OPENMP
841 #pragma omp parallel for
842 #endif
843 for(int i=0;i<nth;i++) {
844 p[i].resize(2);
845 p[i][0]=new JFDigestor(P);
846 p[i][1]=new KFDigestor(P,kfrac,true);
847 }
848
849 // Do calculation
850 arma::vec f=calculate_force(p,tol);
851
852 // Free memory
853 for(size_t i=0;i<p.size();i++)
854 for(size_t j=0;j<p[i].size();j++)
855 delete p[i][j];
856
857 return f;
858 }
859
forceJK(const arma::mat & Pa,const arma::mat & Pb,double tol,double kfrac) const860 arma::vec ERIscreen::forceJK(const arma::mat & Pa, const arma::mat & Pb, double tol, double kfrac) const {
861 if(Pa.n_rows != Nbf || Pa.n_cols != Nbf) {
862 std::ostringstream oss;
863 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pa.n_rows = " << Pa.n_rows << ", Pa.n_cols = " << Pa.n_cols << "!\n";
864 throw std::logic_error(oss.str());
865 }
866 if(Pb.n_rows != Nbf || Pb.n_cols != Nbf) {
867 std::ostringstream oss;
868 oss << "Error in ERIscreen: Nbf = " << Nbf << ", Pb.n_rows = " << Pb.n_rows << ", Pb.n_cols = " << Pb.n_cols << "!\n";
869 throw std::logic_error(oss.str());
870 }
871
872 #ifdef _OPENMP
873 int nth=omp_get_max_threads();
874 #else
875 int nth=1;
876 #endif
877
878 // Get workers
879 std::vector< std::vector<ForceDigestor *> > p(nth);
880 #ifdef _OPENMP
881 #pragma omp parallel for
882 #endif
883 for(int i=0;i<nth;i++) {
884 p[i].resize(3);
885 p[i][0]=new JFDigestor(Pa+Pb);
886 p[i][1]=new KFDigestor(Pa,kfrac,false);
887 p[i][2]=new KFDigestor(Pb,kfrac,false);
888 }
889
890 // Do calculation
891 arma::vec f=calculate_force(p,tol);
892
893 // Free memory
894 for(size_t i=0;i<p.size();i++)
895 for(size_t j=0;j<p[i].size();j++)
896 delete p[i][j];
897
898 return f;
899 }
900