1 Subroutine hf1mkr(Axyz,Bxyz,Cxyz,zan,exinv,ncenters, 2 & alpha,Pxyz,RS,PC,ff,R, 3 & R0,R0C,IJK,NPP,Lp,Lp3,CENTER) 4c $Id$ 5 6 Implicit none 7c::passed 8 Integer ncenters ! number of nuclei 9 Integer NPP ! number of primitive products 10 Integer Lp ! angular momentum sum La+Lb+MXD 11 Integer Lp3 ! number of angular momentum functions 12 Double precision Axyz(3),Bxyz(3) ! Cartesian coordinates of centers 13 Double precision Cxyz(3,ncenters) ! Nuclear cartesian coordinates 14 Double precision zan(ncenters) ! Nuclear charges 15 Double precision exinv(ncenters) ! Nuclear inverse exponents 16 Double precision alpha(2,NPP) ! basis function exponents copy 17 Double precision Pxyz(3,NPP) ! coordinates of product center 18 Double precision RS(NPP) ! scale factor 19 Double precision PC(NPP,3) ! distance from product center to nucleus 20 Double precision ff(2,NPP) ! recursion factors 21 Double precision R(NPP,0:Lp,Lp3) ! R/R000j/RIKJj functions 22 Double precision R0(NPP,Lp3) ! R integrals 23 Double precision R0C(ncenters,NPP,Lp3) ! R integrals per center 24 integer IJK(0:Lp,0:Lp,0:Lp) ! index array 25 Logical CENTER 26c::local 27 integer ic,j,m,n 28 Double precision a,b,f1,f2,alpha_t,beta,PCx,PCy,PCz 29 Double precision PI, PI4 ! pi, 4/pi 30 Parameter (PI=3.1415926535898D0,PI4=4.D0/PI) 31c 32c Define the auxiliary function integrals necessary to compute the nuclear 33c attraction integrals (NAIs). These integrals are scaled by an appropriate 34c factor, RS, defined as 35c 36c / a + b \ 1/2 37c RS = | ------- | 38c \ PI/4 / 39c 40c The scale factor for the Hermite expansion coefficients is assumed to be 41c 42c / PI \ 3/2 / a b __2 \ 43c ES = | ------- | EXP| - ----- AB | 44c \ a + b / \ a + b / 45c 46c Therefore, 47c 48c 2 PI / a b __2 \ 49c ES RS = ------- EXP| - ----- AB | 50c a + b \ a + b / 51c 52c****************************************************************************** 53 54 do 101 m = 1,NPP 55 do 100 j = 1,Lp3 56 R0(m,j) = 0.D0 57 100 continue 58 101 continue 59 60 do 110 m = 1,NPP 61 62c Define the center "P". 63 64 a = alpha(1,m) 65 b = alpha(2,m) 66 67 f1 = a/(a+b) 68 f2 = b/(a+b) 69 70 Pxyz(1,m) = f1*Axyz(1) + f2*Bxyz(1) 71 Pxyz(2,m) = f1*Axyz(2) + f2*Bxyz(2) 72 Pxyz(3,m) = f1*Axyz(3) + f2*Bxyz(3) 73 74c Define the scaling factor. 75 76 RS(m) = sqrt((a+b)*PI4) 77 78 110 continue 79 80c Sum over all centers. 81 82 do 150 ic = 1,ncenters 83 84 beta = exinv(ic) 85 86c Define factors necessary to compute incomplete gamma function and the 87c auxiliary functions. 88 89 do 120 m = 1,NPP 90 91 alpha_t = alpha(1,m) + alpha(2,m) 92 alpha_t = alpha_t/(1.0d0+alpha_t*beta) 93 94 ff(1,m) = RS(m) 95 ff(2,m) = -2.D0*alpha_t 96 97 PCx = Pxyz(1,m) - Cxyz(1,ic) 98 PCy = Pxyz(2,m) - Cxyz(2,ic) 99 PCz = Pxyz(3,m) - Cxyz(3,ic) 100 101 R(m,0,1) = alpha_t*(PCx**2 + PCy**2 + PCz**2) 102 103 PC(m,1) = PCx 104 PC(m,2) = PCy 105 PC(m,3) = PCz 106 107 120 continue 108 109c Evaluate the incomplete gamma function. 110 111 call igamma(R,NPP,Lp) 112 113c Scale for finite nuclear size if necessary 114 115 if (beta .ne. 0.0D0) then 116 do 125 m = 1,NPP 117 alpha_t = sqrt(1.0d0+(alpha(1,m)+alpha(2,m))*beta) 118 do 126 j=0,Lp 119 R(m,j,1) = R(m,j,1)/alpha_t 120 126 continue 121 125 continue 122 end if 123 124c Define the initial auxiliary functions (i.e., R000j, j=1,Lr). 125 126 do 135 j = 0,Lp 127 do 130 m = 1,NPP 128 R(m,j,1) = ff(1,m)*R(m,j,1) 129 ff(1,m) = ff(1,m)*ff(2,m) 130 130 continue 131 135 continue 132 133c Recursively build the remaining auxiliary functions (i.e., RIJKj, j=0). 134 135 call hfmkr(R,IJK,PC,NPP,Lp,Lp3) 136 137c Transfer to R0 array. 138 139 if( CENTER )then 140 do 141 n = 1,Lp3 141 do 140 m = 1,NPP 142 R0C(ic,m,n) = -zan(ic)*R(m,0,n) 143 R0(m,n) = R0(m,n) + R0C(ic,m,n) 144 140 continue 145 141 continue 146 else 147 do 146 n = 1,Lp3 148 do 145 m = 1,NPP 149 R0(m,n) = R0(m,n) - zan(ic)*R(m,0,n) 150 145 continue 151 146 continue 152 end if 153 154 150 continue 155 156 end 157