1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * HF/DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2013
9 * Copyright (c) 2010-2013, 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 "../eriworker.h"
18 #include "../gaunt.h"
19 #include "../linalg.h"
20 #include "../mathf.h"
21 #include "integrals.h"
22 #include "../basis.h"
23 #include "../slaterfit/form_exponents.h"
24
25 /// A. Kumar and P. C. Mishra, Pramana 29 (1987), pp. 385-390.
Ul(int l,int n12,int n34,double z12,double z34)26 double Ul(int l, int n12, int n34, double z12, double z34) {
27 // Prefactor
28 double prefac=fact(n34+l)/pow(z34,n34+l+1);
29
30 // First term
31 double term1=fact(n12-l-1)/pow(z12,n12-l);
32
33 // Second term
34 double term2=0.0;
35 for(int lp=1;lp<=n34+l+1;lp++)
36 term2-=pow(z34,n34+l-lp+1)/fact(n34+l-lp+1) * fact(n12+n34-lp)/pow(z12+z34,n12+n34-lp+1);
37
38 // Third term
39 double term3=0.0;
40 for(int lp=1;lp<=n34-l;lp++)
41 term3+=pow(z34,n34+l-lp+1)*fact(n12+n34-lp) / (fact(n34-l-lp)*pow(z12+z34,n12+n34-lp+1));
42 term3*=fact(n34-l-1)/fact(n34+l);
43
44 return prefac*(term1+term2+term3);
45 }
46
47 /// Gaunt factor. E. U. Condon and G. H. Shortley, The theory of atomic spectra, Cambridge University Press 1963, page 176.
gaunt(int l,int m,int lp,int mp,int L)48 double gaunt(int l, int m, int lp, int mp, int L) {
49 // Cannot couple if
50 if((l+lp+L)%2==1)
51 return 0.0;
52
53 int g=(l+lp+L)/2;
54 if(g<l || g<lp || g<L)
55 return 0.0;
56
57 // printf("1fact with arguments %i, %i, %i, %i, %i, %i\n",2*g-2*lp,g,g-l,g-lp,g-L,2*g+1);
58 double fac1=pow(-1.0,g-l-mp)*fact(2*g-2*lp)*fact(g)/(fact(g-l)*fact(g-lp)*fact(g-L)*fact(2*g+1));
59 // printf("fac1 computed\n"); fflush(stdout);
60
61 // printf("2fact with arguments %i, %i, %i, %i, %i, %i\n",L-m-mp,l+m,lp+mp,lp-mp,L+m+mp,l-m);
62 double fac2=sqrt((2*l+1)*(2*lp+1)*(2*L+1)*fact(L-m-mp)*fact(l+m)*fact(lp+mp)*fact(lp-mp)/(2*(fact(L+m+mp)*fact(l-m))));
63 // printf("fac2 computed\n"); fflush(stdout);
64
65 // Third term,
66 double fac3=0.0;
67
68 // Minimum of t
69 int tmin=0;
70 int tmin1=-(L+m+mp);
71 int tmin2=-(l-lp+m+mp);
72 if(tmin<tmin1)
73 tmin=tmin1;
74 if(tmin<tmin2)
75 tmin=tmin2;
76
77 int tmax=lp-mp;
78 int tmax2=l+lp-m-mp;
79 int tmax3=L-m-mp;
80 if(tmax>tmax2)
81 tmax=tmax2;
82 if(tmax>tmax3)
83 tmax=tmax3;
84
85 for(int t=tmin;t<=tmax;t++) {
86 // printf("3fact with arguments %i, %i, %i, %i, %i, %i\n",L+m+mp+t,l+lp-m-mp-t,L-m-mp-t,l-lp+m+mp+t,lp-mp-t,t);
87 // fflush(stdout);
88 fac3+=pow(-1.0,t)*fact(L+m+mp+t)*fact(l+lp-m-mp-t)/(fact(L-m-mp-t)*fact(l-lp+m+mp+t)*fact(lp-mp-t)*fact(t));
89 }
90 return fac1*fac2*fac3;
91 }
92
93 // Angular factor, Condon-Shortley page 175
ck(int k,int l,int m,int lp,int mp)94 double ck(int k, int l, int m, int lp, int mp) {
95 if(k<abs(m-mp))
96 return 0.0;
97
98 return sqrt(2.0/(2*k+1))*gaunt(k,m-mp,lp,mp,l);
99 }
100
101 // Normalization coefficient
normalization(int n,double z)102 double normalization(int n, double z) {
103 return sqrt(pow(2.0*z,2*n+1)/fact(2*n));
104 }
105
106 // Primitive overlap integral
overlap_primitive(int nab,double zab)107 double overlap_primitive(int nab, double zab) {
108 return fact(nab)/pow(zab,nab+1);
109 }
110
111 // Repulsion integral (ab|cd) = \int \phi_a(r_1) \phi_b(r_1) r^{-1}_{12} \phi_c(r_2) \phi_d(r_2) d^3 r_1 d^3 r_2
ERI_unnormalized(int na,int nb,int nc,int nd,double za,double zb,double zc,double zd,int la,int ma,int lb,int mb,int lc,int mc,int ld,int md)112 double ERI_unnormalized(int na, int nb, int nc, int nd, double za, double zb, double zc, double zd, int la, int ma, int lb, int mb, int lc, int mc, int ld, int md) {
113 // Remember complex conjugation in ERI
114 if(mb-ma!=md-mc) {
115 // printf("m different - ma=%i, mb=%i => %i, mc=%i, md=%i => %i\n",ma,mb,mb-ma,mc,md,md-mc);
116 return 0.0;
117 }
118
119 // What is maximum angular momentum we can couple to?
120 // int lmax=la+lc;
121 /*
122 if(lb+ld<lmax)
123 lmax=lb+ld;
124 */
125
126 // The functions la and lb can couple to |la-lb| <= l <= la+lb
127 // but we also must require that |lc-ld| <= l <= lc+ld
128 // and |m|<=l
129 int lmin1=std::abs(la-lb);
130 int lmin2=std::abs(lc-ld);
131 int lmin3=abs(mb-ma);
132 int lmin=std::max(lmin1,std::max(lmin2,lmin3));
133
134 int lmax=std::min(la+lb,lc+ld);
135
136 double eri=0.0;
137
138 for(int l=lmin;l<=lmax;l++)
139 eri+=Ul(l,na+nb,nc+nd,za+zb,zc+zd)*ck(l,la,ma,lb,mb)*ck(l,lc,mc,ld,md);
140
141 return eri;
142 }
143
ERI(int na,int nb,int nc,int nd,double za,double zb,double zc,double zd,int la,int ma,int lb,int mb,int lc,int mc,int ld,int md)144 double ERI(int na, int nb, int nc, int nd, double za, double zb, double zc, double zd, int la, int ma, int lb, int mb, int lc, int mc, int ld, int md) {
145 double eri=ERI_unnormalized(na,nb,nc,nd,za,zb,zc,zd,la,ma,lb,mb,lc,mc,ld,md);
146 // Plug in normalization
147 eri*=normalization(na,za)*normalization(nb,zb)*normalization(nc,zc)*normalization(nd,zd);
148 return eri;
149 }
150
coulomb_overlap(const bf_t & a,const bf_t & b)151 double coulomb_overlap(const bf_t & a, const bf_t & b) {
152 bf_t dummy;
153 dummy.zeta=0.0;
154 dummy.n=1;
155 dummy.l=0;
156 dummy.m=0;
157
158 // Overlap wrt overlap-normalized functions
159 // return ERI_unnormalized(a,dummy,b,dummy)*normalization(a.n,a.zeta)*normalization(b.n,b.zeta);
160
161 // Overlap wrt ERI-normalized functions
162 return ERI_unnormalized(a,dummy,b,dummy)/sqrt(ERI_unnormalized(a,dummy,a,dummy)*ERI_unnormalized(b,dummy,b,dummy));
163 }
164
three_overlap(int na,int nc,int nd,int la,int ma,int lc,int mc,int ld,int md,double za,double zc,double zd)165 double three_overlap(int na, int nc, int nd, int la, int ma, int lc, int mc, int ld, int md, double za, double zc, double zd) {
166 if(ma!=md-mc) {
167 return 0.0;
168 }
169 // The functions lc and ld can couple to |lc-ld| <= l <= lc+ld
170 // which must thus be the same as la
171 if(la<abs(lc-ld) || la>lc+ld)
172 return 0.0;
173
174 //return overlap_primitive(na+nc+nd,za+zc+zd)*gaunt_coefficient(la,ma,lc,mc,ld,md)*normalization(na,za)*normalization(nc,zc)*normalization(nd,zd);
175 return overlap_primitive(na+nc+nd,za+zc+zd)*ck(la,lc,mc,ld,md)*normalization(na,za)*normalization(nc,zc)*normalization(nd,zd);
176 }
177
178
179 // Do Gaussian expansion of ERI
gaussian_ERI(int la,int ma,int lb,int mb,int lc,int mc,int ld,int md,double za,double zb,double zc,double zd,int nfit)180 double gaussian_ERI(int la, int ma, int lb, int mb, int lc, int mc, int ld, int md, double za, double zb, double zc, double zd, int nfit) {
181 ERIWorker eri(std::max(la,lb),nfit);
182 std::vector<double> eris;
183
184 std::vector<contr_t> ca=slater_fit(za,la,nfit,false);
185 std::vector<contr_t> cb=slater_fit(zb,lb,nfit,false);
186 std::vector<contr_t> cc=slater_fit(zc,lc,nfit,false);
187 std::vector<contr_t> cd=slater_fit(zd,ld,nfit,false);
188
189 GaussianShell ash(la,true,ca);
190 GaussianShell bsh(lb,true,cb);
191 GaussianShell csh(lc,true,cc);
192 GaussianShell dsh(ld,true,cd);
193
194 coords_t cen={0.0, 0.0, 0.0};
195
196 ash.set_first_ind(0);
197 ash.set_center(cen,0);
198
199 bsh.set_first_ind(ash.get_last_ind()+1);
200 bsh.set_center(cen,0);
201
202 csh.set_first_ind(bsh.get_last_ind()+1);
203 csh.set_center(cen,0);
204
205 dsh.set_first_ind(csh.get_last_ind()+1);
206 dsh.set_center(cen,0);
207
208
209 ash.convert_contraction();
210 ash.normalize();
211 bsh.convert_contraction();
212 bsh.normalize();
213 csh.convert_contraction();
214 csh.normalize();
215 dsh.convert_contraction();
216 dsh.normalize();
217
218 eri.compute(&ash,&bsh,&csh,&dsh);
219 eris=eri.get();
220
221 // Lengths
222 int bn=2*lb+1;
223 int cn=2*lc+1;
224 int dn=2*ld+1;
225
226 // Indices
227 int ai=la+ma;
228 int bi=lb+mb;
229 int ci=lc+mc;
230 int di=ld+md;
231
232 return eris[((ai*bn+bi)*cn+ci)*dn+di];
233 }
234
235 // Repulsion integral (ab|cd) = \int \phi_a(r_1) \phi_b(r_1) r^{-1}_{12} \phi_c(r_2) \phi_d(r_2) d^3 r_1 d^3 r_2
overlap(int na,int nb,double za,double zb,int la,int ma,int lb,int mb)236 double overlap(int na, int nb, double za, double zb, int la, int ma, int lb, int mb) {
237 if(la!=lb || ma!=mb)
238 return 0.0;
239
240 return normalization(na,za)*normalization(nb,zb)*overlap_primitive(na+nb,za+zb);
241 }
242
243 // Nuclear attraction integral
nuclear(int na,int nb,double za,double zb,int la,int ma,int lb,int mb)244 double nuclear(int na, int nb, double za, double zb, int la, int ma, int lb, int mb) {
245 if(la!=lb || ma!=mb)
246 return 0.0;
247
248 return -normalization(na,za)*normalization(nb,zb)*overlap_primitive(na+nb-1,za+zb);
249 }
250
251 // Kinetic energy integral
kinetic(int na,int nb,double za,double zb,int la,int ma,int lb,int mb)252 double kinetic(int na, int nb, double za, double zb, int la, int ma, int lb, int mb) {
253 if(la!=lb || ma!=mb)
254 return 0.0;
255
256 // First term is
257 double term1=(lb*(lb+1)-nb*(nb-1))*overlap_primitive(na+nb-2,za+zb);
258 // Second term is
259 double term2=2*zb*nb*overlap_primitive(na+nb-1,za+zb);
260 // Third term is
261 double term3=-zb*zb*overlap_primitive(na+nb,za+zb);
262
263 return 0.5*normalization(na,za)*normalization(nb,zb)*(term1+term2+term3);
264 }
265
266 /// Form overlap matrix
overlap(const std::vector<bf_t> & basis)267 arma::mat overlap(const std::vector<bf_t> & basis) {
268 arma::mat S(basis.size(),basis.size());
269 S.zeros();
270 #ifdef _OPENMP
271 #pragma omp parallel for schedule(dynamic)
272 #endif
273 for(size_t i=0;i<basis.size();i++)
274 for(size_t j=0;j<=i;j++) {
275 double el=overlap(basis[i].n,basis[j].n,basis[i].zeta,basis[j].zeta,basis[i].l,basis[i].m,basis[j].l,basis[j].m);
276 S(i,j)=el;
277 S(j,i)=el;
278 }
279
280 return S;
281 }
282
kinetic(const std::vector<bf_t> & basis)283 arma::mat kinetic(const std::vector<bf_t> & basis) {
284 arma::mat T(basis.size(),basis.size());
285 T.zeros();
286 #ifdef _OPENMP
287 #pragma omp parallel for schedule(dynamic)
288 #endif
289 for(size_t i=0;i<basis.size();i++)
290 for(size_t j=0;j<=i;j++) {
291 double el=kinetic(basis[i].n,basis[j].n,basis[i].zeta,basis[j].zeta,basis[i].l,basis[i].m,basis[j].l,basis[j].m);
292 T(i,j)=el;
293 T(j,i)=el;
294 }
295
296 return T;
297 }
298
nuclear(const std::vector<bf_t> & basis,int Z)299 arma::mat nuclear(const std::vector<bf_t> & basis, int Z) {
300 arma::mat V(basis.size(),basis.size());
301 V.zeros();
302 #ifdef _OPENMP
303 #pragma omp parallel for schedule(dynamic)
304 #endif
305 for(size_t i=0;i<basis.size();i++)
306 for(size_t j=0;j<=i;j++) {
307 double el=nuclear(basis[i].n,basis[j].n,basis[i].zeta,basis[j].zeta,basis[i].l,basis[i].m,basis[j].l,basis[j].m);
308 V(i,j)=el;
309 V(j,i)=el;
310 }
311
312 return Z*V;
313 }
314
coulomb(const std::vector<bf_t> & basis,const arma::mat & P)315 arma::mat coulomb(const std::vector<bf_t> & basis, const arma::mat & P) {
316 arma::mat J(basis.size(),basis.size());
317 J.zeros();
318 #ifdef _OPENMP
319 #pragma omp parallel for schedule(dynamic)
320 #endif
321 for(size_t i=0;i<basis.size();i++)
322 for(size_t j=0;j<=i;j++) {
323 double el=0.0;
324
325 for(size_t k=0;k<basis.size();k++)
326 for(size_t l=0;l<basis.size();l++)
327 el+=P(k,l)*ERI(basis[i],basis[j],basis[k],basis[l]);
328
329 J(i,j)=el;
330 J(j,i)=el;
331 }
332
333 return J;
334 }
335
integral(const std::vector<bf_t> & fitbas,bool coulomb=false)336 arma::vec integral(const std::vector<bf_t> & fitbas, bool coulomb=false) {
337 arma::vec r(fitbas.size());
338 r.zeros();
339 for(size_t i=0;i<fitbas.size();i++)
340 if(fitbas[i].l==0) {
341 // Integral is
342 r(i)=sqrt(4.0*M_PI)*pow(fitbas[i].zeta,-fitbas[i].n-2)*fact(fitbas[i].n+1);
343
344 if(!coulomb) {
345 r(i)*=normalization(fitbas[i].n,fitbas[i].zeta);
346 } else {
347 bf_t dummy;
348 dummy.zeta=0.0;
349 dummy.n=1;
350 dummy.l=0;
351 dummy.m=0;
352
353 r(i)/=sqrt(ERI_unnormalized(fitbas[i],dummy,fitbas[i],dummy));
354 }
355 }
356
357 return r;
358 }
359
coulomb_coul(const std::vector<bf_t> & basis,const std::vector<bf_t> & fitbas,const arma::mat & P)360 arma::mat coulomb_coul(const std::vector<bf_t> & basis, const std::vector<bf_t> & fitbas, const arma::mat & P) {
361 bf_t dummy;
362 dummy.zeta=0.0;
363 dummy.n=1;
364 dummy.l=0;
365 dummy.m=0;
366
367 // Coulomb overlap matrix
368 arma::mat V(fitbas.size(),fitbas.size());
369 #ifdef _OPENMP
370 #pragma omp parallel for schedule(dynamic)
371 #endif
372 for(size_t i=0;i<fitbas.size();i++)
373 for(size_t j=0;j<=i;j++) {
374 double el=ERI_unnormalized(fitbas[i],dummy,fitbas[j],dummy)/sqrt(ERI_unnormalized(fitbas[i],dummy,fitbas[i],dummy)*ERI_unnormalized(fitbas[j],dummy,fitbas[j],dummy));
375 V(i,j)=el;
376 V(j,i)=el;
377 }
378 arma::mat Vinvh=BasOrth(V,false);
379 arma::mat Vinv=Vinvh*Vinvh;
380
381 // Three-center integrals
382 arma::cube thr(fitbas.size(),basis.size(),basis.size());
383 #ifdef _OPENMP
384 #pragma omp parallel for schedule(dynamic)
385 #endif
386 for(size_t i=0;i<fitbas.size();i++)
387 for(size_t k=0;k<basis.size();k++)
388 for(size_t l=0;l<basis.size();l++) {
389 thr(i,k,l)=ERI_unnormalized(fitbas[i],dummy,basis[k],basis[l])/sqrt(ERI_unnormalized(fitbas[i],dummy,fitbas[i],dummy))*normalization(basis[k].n,basis[k].zeta)*normalization(basis[l].n,basis[l].zeta);
390 }
391
392 // Compute fitting coefficients
393 arma::vec gamma(fitbas.size());
394 gamma.zeros();
395 #ifdef _OPENMP
396 #pragma omp parallel for schedule(dynamic)
397 #endif
398 for(size_t i=0;i<fitbas.size();i++)
399 for(size_t k=0;k<basis.size();k++)
400 for(size_t l=0;l<basis.size();l++) {
401 gamma(i)+=P(k,l)*thr(i,k,l);
402 }
403
404 // printf("Norm of expanded density is %.12f (coul).\n",arma::dot(gamma,integral(fitbas,true)));
405 //printf("Norm of expanded density is %.12f (ovl).\n",arma::dot(gamma,integral(fitbas,false)));
406
407 // Translate coefficients
408 gamma=Vinv*gamma;
409 // printf("Norm of expanded density is %.12f (coul).\n",arma::dot(gamma,integral(fitbas,true)));
410 // printf("Norm of expanded density is %.12f (ovl).\n",arma::dot(gamma,integral(fitbas,false)));
411
412 // Collect Coulomb matrix
413 arma::mat J(basis.size(),basis.size());
414 J.zeros();
415 #ifdef _OPENMP
416 #pragma omp parallel for schedule(dynamic)
417 #endif
418 for(size_t i=0;i<basis.size();i++)
419 for(size_t j=0;j<=i;j++) {
420 double el=0.0;
421
422 for(size_t k=0;k<fitbas.size();k++)
423 el+=gamma(k)*thr(k,i,j);
424
425 J(i,j)=el;
426 J(j,i)=el;
427 }
428
429 return J;
430 }
431
coulomb_ovl(const std::vector<bf_t> & basis,const std::vector<bf_t> & fitbas,const arma::mat & P)432 arma::mat coulomb_ovl(const std::vector<bf_t> & basis, const std::vector<bf_t> & fitbas, const arma::mat & P) {
433 // Form overlap matrix
434 arma::mat S(fitbas.size(),fitbas.size());
435 #ifdef _OPENMP
436 #pragma omp parallel for schedule(dynamic)
437 #endif
438 for(size_t i=0;i<fitbas.size();i++)
439 for(size_t j=0;j<=i;j++) {
440 double el=overlap(fitbas[i],fitbas[j]);
441 S(i,j)=el;
442 S(j,i)=el;
443 }
444
445 // Inverse overlap
446 arma::mat Sinvh=BasOrth(S,false);
447 arma::mat Sinv=Sinvh*Sinvh;
448
449 // Coulomb overlap matrix
450 arma::mat V(fitbas.size(),fitbas.size());
451 #ifdef _OPENMP
452 #pragma omp parallel for schedule(dynamic)
453 #endif
454 for(size_t i=0;i<fitbas.size();i++)
455 for(size_t j=0;j<=i;j++) {
456 double el=coulomb_overlap(fitbas[i],fitbas[j]);
457 V(i,j)=el;
458 V(j,i)=el;
459 }
460
461 // Three-center integrals
462 arma::cube thr(fitbas.size(),basis.size(),basis.size());
463 #ifdef _OPENMP
464 #pragma omp parallel for schedule(dynamic)
465 #endif
466 for(size_t i=0;i<fitbas.size();i++)
467 for(size_t k=0;k<basis.size();k++)
468 for(size_t l=0;l<basis.size();l++) {
469 thr(i,k,l)=three_overlap(fitbas[i],basis[k],basis[l]);
470 }
471
472 // Compute fitting coefficients
473 arma::vec gamma(fitbas.size());
474 gamma.zeros();
475 #ifdef _OPENMP
476 #pragma omp parallel for schedule(dynamic)
477 #endif
478 for(size_t i=0;i<fitbas.size();i++)
479 for(size_t k=0;k<basis.size();k++)
480 for(size_t l=0;l<basis.size();l++) {
481 gamma(i)+=P(k,l)*thr(i,k,l);
482 }
483 // printf("Norm of expanded density was %.12f.\n",arma::dot(gamma,integral(fitbas)));
484 gamma=Sinv*gamma;
485
486 // printf("Norm of expanded density was %.12f.\n",arma::dot(gamma,integral(fitbas)));
487
488 /*
489 double Nelfit=arma::dot(gamma,integral(fitbas));
490 double Nel=arma::trace(P*overlap(basis));
491 // Renormalize
492 gamma*=Nel/Nelfit;
493 */
494
495 // Translate coefficients
496 gamma=Sinv*V*gamma;
497
498 // Collect Coulomb matrix
499 arma::mat J(basis.size(),basis.size());
500 J.zeros();
501 #ifdef _OPENMP
502 #pragma omp parallel for schedule(dynamic)
503 #endif
504 for(size_t i=0;i<basis.size();i++)
505 for(size_t j=0;j<=i;j++) {
506 double el=0.0;
507
508 for(size_t k=0;k<fitbas.size();k++)
509 el+=gamma(k)*thr(k,i,j);
510
511 J(i,j)=el;
512 J(j,i)=el;
513 }
514
515 return J;
516 }
517
coulomb(const std::vector<bf_t> & basis,const std::vector<bf_t> & fitbas,const arma::mat & P)518 arma::mat coulomb(const std::vector<bf_t> & basis, const std::vector<bf_t> & fitbas, const arma::mat & P) {
519 return coulomb_coul(basis,fitbas,P);
520 // return coulomb_ovl(basis,fitbas,P);
521 }
522
exchange(const std::vector<bf_t> & basis,const arma::mat & P)523 arma::mat exchange(const std::vector<bf_t> & basis, const arma::mat & P) {
524 arma::mat K(basis.size(),basis.size());
525 K.zeros();
526
527 #ifdef _OPENMP
528 #pragma omp parallel for schedule(dynamic)
529 #endif
530 for(size_t i=0;i<basis.size();i++)
531 for(size_t j=0;j<=i;j++) {
532 double el=0.0;
533
534 for(size_t k=0;k<basis.size();k++)
535 for(size_t l=0;l<basis.size();l++)
536 el+=P(k,l)*ERI(basis[i],basis[k],basis[j],basis[l]);
537
538 K(i,j)=el;
539 K(j,i)=el;
540 }
541
542 return K;
543 }
544