1 Subroutine hf2mkr(Axyz,Bxyz,Cxyz,Dxyz,alpha_p,alpha_q,R0,IJK, 2 & Pxyz,Qxyz,PQ,ff,R,ffscr,rscr,NPP,NPQ,Lr,Lr3) 3c 4c $Id$ 5 6 implicit none 7#include "case.fh" 8c 9 double precision PI,P1 10 Parameter (PI=3.1415926535898D0,P1=4.D0/PI) 11 12c--> ints 13 integer npp,npq 14 integer Lr, Lr3 15 integer IJK(0:Lr,0:Lr,0:Lr) 16 17c--> Cartesian Coordinates 18 double precision Axyz(3),Bxyz(3),Cxyz(3),Dxyz(3) 19 20c--> Exponents 21 22 double precision alpha_p(2,NPP),alpha_q(2,NPQ) 23 24c--> Auxiliary Function Integrals & Index 25 26 double precision R0((NPP*NPQ),Lr3) 27 28c--> Scratch Space 29 30 double precision Pxyz(3,NPP),Qxyz(3,NPQ),PQ((NPP*NPQ),3), 31 & ff(2,(NPP*NPQ)), R((NPP*NPQ),0:Lr,Lr3) 32c 33 double precision ffscr(2,(NPP*NPQ)), Rscr((NPP*NPQ),0:Lr,Lr3), 34 & rhoscaled,radd 35 36 integer mp, mq, mr, j, n 37 double precision a,b,c,d,f1,f2,p,q,rho 38 double precision PQx, PQy, PQz 39c 40c Evaluate the auxiliary function integrals. These integrals are scaled by a 41c factor appropriate for ERIs, defined as 42c 43c 1/2 44c ES = ((4/PI)*(pq/(p+q))) where p = a + b and q = c + d. 45c 46c****************************************************************************** 47c 48c Define the center "P". 49 50 do 100 mp = 1,NPP 51 52 a = alpha_p(1,mp) 53 b = alpha_p(2,mp) 54 55 f1 = a/(a+b) 56 f2 = b/(a+b) 57 58 Pxyz(1,mp) = f1*Axyz(1) + f2*Bxyz(1) 59 Pxyz(2,mp) = f1*Axyz(2) + f2*Bxyz(2) 60 Pxyz(3,mp) = f1*Axyz(3) + f2*Bxyz(3) 61 62 100 continue 63 64c Define the center "Q". 65 66 do 110 mq = 1,NPQ 67 68 c = alpha_q(1,mq) 69 d = alpha_q(2,mq) 70 71 f1 = c/(c+d) 72 f2 = d/(c+d) 73 74 Qxyz(1,mq) = f1*Cxyz(1) + f2*Dxyz(1) 75 Qxyz(2,mq) = f1*Cxyz(2) + f2*Dxyz(2) 76 Qxyz(3,mq) = f1*Cxyz(3) + f2*Dxyz(3) 77 78 110 continue 79 80c Define factors necessary to compute incomplete gamma function and the 81c auxiliary functions. 82 83 mr = 0 84 do 125 mp = 1,NPP 85 do 120 mq = 1,NPQ 86 mr = mr + 1 87 88 p = alpha_p(1,mp) + alpha_p(2,mp) 89 q = alpha_q(1,mq) + alpha_q(2,mq) 90 91 rho = p*q/(p+q) 92 93 PQx = Pxyz(1,mp) - Qxyz(1,mq) 94 PQy = Pxyz(2,mp) - Qxyz(2,mq) 95 PQz = Pxyz(3,mp) - Qxyz(3,mq) 96 97 PQ(mr,1) = PQx 98 PQ(mr,2) = PQy 99 PQ(mr,3) = PQz 100c 101 R(mr,0,1) = rho*(PQx**2 + PQy**2 + PQz**2) 102 ff(1,mr) = 2.d0*dsqrt(rho/pi) 103 ff(2,mr) = -2.D0*rho 104c 105 if (doscreen) then 106 rhoscaled = rho 107 call case_md(rhoscaled) 108 Rscr(mr,0,1) = rhoscaled*(PQx**2 + PQy**2 + PQz**2) 109 ffscr(1,mr) = 2.d0*dsqrt(rhoscaled/pi) 110 ffscr(2,mr) = -2.D0*rhoscaled 111 end if 112c 113 120 continue 114 125 continue 115 116c Evaluate the incomplete gamma function. 117 118 call igamma(R,(NPP*NPQ),Lr) 119 if (doscreen) call igamma(Rscr,(NPP*NPQ),Lr) 120 121c Define the initial auxiliary functions (i.e., R000j, j=1,Lr). 122 123 do 135 j = 0,Lr 124 do 130 mr = 1,(NPP*NPQ) 125 R(mr,j,1) = ff(1,mr)*R(mr,j,1) 126 ff(1,mr) = ff(1,mr)*ff(2,mr) 127 130 continue 128 if (doscreen) then 129 do mr = 1,(NPP*NPQ) 130 Rscr(mr,j,1) = ffscr(1,mr)*Rscr(mr,j,1) 131 ffscr(1,mr) = ffscr(1,mr)*ffscr(2,mr) 132 enddo 133 endif 134 135 continue 135 136c Recursively build the remaining auxiliary functions (i.e., RIJKj, j=0). 137 138 call hfmkr(R,IJK,PQ,(NPP*NPQ),Lr,Lr3) 139 if (doscreen) call hfmkr(Rscr,IJK,PQ,(NPP*NPQ),Lr,Lr3) 140 141c Transfer to R0 array. 142 143c write(6,*) "cam_alpha:",cam_alpha 144c write(6,*) "cam_beta:",cam_beta 145c write(6,*) "doscreen:",doscreen 146c write(6,*) "cam_srhf:",cam_srhf 147c 148 if (doscreen) then 149 do n = 1,Lr3 150 if (cam_srhf) then ! for pure short-range HF (1-erf(r)/r) 151 do mr = 1,(NPP*NPQ) 152 R0(mr,n) = R(mr,0,n)-Rscr(mr,0,n) 153 enddo 154 else 155 do mr = 1,(NPP*NPQ) 156 R0(mr,n) = cam_alpha*R(mr,0,n)+cam_beta*Rscr(mr,0,n) 157 enddo 158 endif 159 enddo 160 else 161 do n = 1,Lr3 162 do mr = 1,(NPP*NPQ) 163 R0(mr,n) = R(mr,0,n) 164 enddo 165 enddo 166 endif 167 168 end 169