1 Subroutine hf3mkr(Axyz,Bxyz,Cxyz,alpha,Gxyz, 2 & RS,GC,ff,R,R0,IJK,Nabc,Lg,Lg3) 3c 4c $Id$ 5c 6 Implicit none 7c::passed 8 integer Nabc, Lg, Lg3 9c--> Cartesian Coordinates for centers a, b, c 10 11 Double Precision Axyz(3),Bxyz(3),Cxyz(3) ! [input] 12 13c--> Exponents (1:3,*) and ES prefactors (4:4,*) 14 15 Double Precision alpha(4,Nabc) ! [input] 16 17c--> Auxiliary Function Integrals & Index 18 19 Double Precision R0(Nabc,Lg3) ! [output] 20 Integer IJK(0:Lg,0:Lg,0:Lg) ! [output] 21 22c--> Scratch Space 23 24 Double Precision Gxyz(3,Nabc), GC(Nabc,3) 25 Double Precision RS(Nabc), ff(2,Nabc), R(Nabc,0:Lg,Lg3) 26c::local 27 double precision PI, PI4 28 Parameter (PI=3.1415926535898D0,PI4=4.D0/PI) 29c 30 double precision a, b, c, abci 31*rak: double precision ab, abi 32 double precision Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz 33*rak: Double Precision Px, Py, Pz 34 double precision GCx, GCy, GCz 35 double precision alpha_t 36*acc_debug: double precision accy 37*acc_debug: integer accy_cnt 38*acc_debug: logical reached 39 integer mp, j, m, n 40 41c 42c Define the auxiliary function integrals necessary to compute the three 43c center nuclear attraction integrals (NAIs). These integrals are scaled 44c by an appropriate factor, RS, defined as 45c 46c / a + b + c \ 1/2 47c RS = | ----------- | 48c \ PI/4 / 49c 50c 51c****************************************************************************** 52* call dfill((Nabc*(Lg+1)*Lg3), 0.0d00 , R, 1) 53* call dfill((Nabc*Lg3),0.0d00, r0, 1) 54 55c Define the center "P" plus C to get "G" center. 56 57 Ax = Axyz(1) 58 Ay = Axyz(2) 59 Az = Axyz(3) 60 61 Bx = Bxyz(1) 62 By = Bxyz(2) 63 Bz = Bxyz(3) 64 65 Cx = Cxyz(1) 66 Cy = Cxyz(2) 67 Cz = Cxyz(3) 68 69 do 00100 mp = 1,Nabc 70 71 a = alpha(1,mp) 72 b = alpha(2,mp) 73 c = alpha(3,mp) 74 75*rak: ab = a + b 76*rak: abi = 1/ab 77*rak: 78*rak: px = abi*(a*Ax + b*Bx) 79*rak: py = abi*(a*Ay + b*By) 80*rak: pz = abi*(a*Az + b*Bz) 81*rak: abci = 1/(ab+c) 82 83 abci = 1/(a+b+c) 84 85*rak: Gxyz(1,mp) = abci*(ab*px + c*Cx) 86*rak: Gxyz(2,mp) = abci*(ab*py + c*Cy) 87*rak: Gxyz(3,mp) = abci*(ab*pz + c*Cz) 88 89 Gxyz(1,mp) = abci*(a*Ax + b*Bx + c*Cx) 90 Gxyz(2,mp) = abci*(a*Ay + b*By + c*Cy) 91 Gxyz(3,mp) = abci*(a*Az + b*Bz + c*Cz) 92 93c Define the scaling factor. 94 95 RS(mp) = sqrt((a+b+c)*PI4) 96 9700100 continue 98 99c Define factors necessary to compute incomplete gamma function and the 100c auxiliary functions. 101 102 do 00200 m = 1,Nabc 103 104 alpha_t = alpha(1,m) + alpha(2,m) + alpha(3,m) 105 106 ff(1,m) = RS(m) 107 ff(2,m) = -2.D0*alpha_t 108 109 GCx = Gxyz(1,m) - Cx 110 GCy = Gxyz(2,m) - Cy 111 GCz = Gxyz(3,m) - Cz 112 113 R(m,0,1) = alpha_t*(GCx*GCx + GCy*GCy + GCz*GCz) 114 115 GC(m,1) = GCx 116 GC(m,2) = GCy 117 GC(m,3) = GCz 118 11900200 continue 120 121c Evaluate the incomplete gamma function. 122 123 call igamma(R,Nabc,Lg) 124 125*acc_debug: accy = 1.0d-30 126*acc_debug: accy_cnt = 0 127*acc_debug: reached = .false. 128*acc_debug:00001 continue 129*acc_debug: call igamma_acc(R,Nabc,Lg,accy,reached) 130*acc_debug: if (.not.reached) then 131*acc_debug: accy_cnt = accy_cnt + 1 132*acc_debug: accy = accy/5.0d00 133*acc_debug: write(6,*)' accy reset to ',accy,' count is ',accy_cnt 134*acc_debug: goto 00001 135*acc_debug: endif 136 137c Define the initial auxiliary functions (i.e., R000j, j=1,Lg). 138 139 do 00300 j = 0,Lg 140 do 00400 m = 1,Nabc 141 R(m,j,1) = ff(1,m)*R(m,j,1) 142 ff(1,m) = ff(1,m)*ff(2,m) 14300400 continue 14400300 continue 145 146c Recursively build the remaining auxiliary functions (i.e., RIJKj, j=0). 147 148 call hfmkr(R,IJK,GC,Nabc,Lg,Lg3) 149 150c Transfer to R0 array. 151 152 do 00500 n = 1,Lg3 153 do 00600 m = 1,Nabc 154 R0(m,n) = R(m,0,n) 15500600 continue 15600500 continue 157 158 end 159