1 Subroutine hferi(Ep,Eq,R0,IJK,ERI,E3,sum,MPP,MPQ, 2 & NPP,NPQ,Nint,La,Lb,Lc,Ld,Lr,MXD, 3 & canAB,canCD,canPQ) 4c $Id$ 5 6 Implicit none 7#include "sh_order.fh" 8 integer mpp,mpq,npp,npq,nint 9 integer la,lb,lc,ld,lr,mxd 10 integer ld2,lqmax,lqmax3 11 integer iq,jq,nn,ia,ja,ka,mb_limit 12 integer ij,mb,ib,jb,kb,ip,jp,kp,mp,ir,jr,lq,nq 13 integer kr,nr,mq,mc,ic,jc,kc,md_limit,kl,md,id,jd,kd 14 integer lb2,lc2,kq,ma,nc,la6,lb6,lc6,ld6 15 integer nq1,iq1,iqfin,iqlast 16 Logical canAB,canCD,canPQ 17 18c--> Hermite Linear Expansion Coefficients 19 20 Double precision Ep(3,NPP,0:MXD,0:(La+Lb),0:La,0:Lb) 21 Double precision Eq(3,NPQ,0:MXD,0:(Lc+Ld),0:Lc,0:Ld) 22 23c--> Auxiliary Function Integrals & Index 24 25 Double precision R0(MPQ,MPP,*) 26 Integer IJK(0:Lr,0:Lr,0:Lr) 27 28c--> ERI 29 30 Double precision ERI(Nint) 31 32c--> Scratch Space 33 34 Double precision E3(*),sum(MPQ,*),erinn,zot, 35 , zot1,zot2 36 Integer Nxyz(3),nqmx 37c 38c Compute electron repulsion integrals (ERI). 39c 40c Formula: 41c 42c __ 43c \ Ic,Id;n10 Jc,Jd;n11 Kc,Kd;n12 44c ERI = / Ex Ey Ez SUM 45c -- Iq Jq Kq Iq,Jq,Kq 46c Iq,Jq,Kq 47c 48c __ 49c \ Lq Ia,Ib;n7 Ja,Jb;n8 Ka,Kb;n9 50c SUM = / (-1) Ex Ey Ez R 51c Iq,Jq,Kq -- Ip Jp Kp Ir,Jr,Kr 52c Ip,Jp,Kp 53c 54c Ir = (Ip+n1) + (Iq+n4) 55c where Jr = (Jp+n2) + (Jq+n5) 56c Kr = (Kp+n3) + (Kq+n6) 57c 58c and Lq = (Iq+n4) + (Jq+n5) + (Kq+n6) 59c 60c N.B. For simple ERI (i.e., no derivative integrals) n[1-12] = 0! 61c 62c****************************************************************************** 63*rak: integer num_pg, num_tot, num_qg 64*rak: save num_pg, num_tot, num_qg 65*rak: data num_pg /0/ 66*rak: data num_qg /0/ 67*rak: data num_tot /0/ 68*rak: 69*rak: if(MPP.ge.MPQ) then 70*rak: num_pg = num_pg + 1 71*rak: else 72*rak: num_qg = num_qg + 1 73*rak: endif 74*rak: num_tot = num_tot + 1 75*rak: if (num_tot.lt.3200.or.mod(num_tot,1000).eq.0) then 76*rak: write(6,*)' num_pg ',num_pg 77*rak: write(6,*)' num_qg ',num_qg 78*rak: write(6,*)' num_tot ',num_tot 79*rak: endif 80c 81c General case: [ab|cd] 82 83c Define the number of shell components on each center. 84 85 Lb2 = ((Lb+1)*(Lb+2))/2 86 Lc2 = ((Lc+1)*(Lc+2))/2 87 Ld2 = ((Ld+1)*(Ld+2))/2 88 la6=(la*(la+1)*(la+2))/6 89 lb6=(lb*(lb+1)*(lb+2))/6 90 lc6=(lc*(lc+1)*(lc+2))/6 91 ld6=(ld*(ld+1)*(ld+2))/6 92 md_limit = Ld2 93 94c Initialize the block of ERIs. 95 96 Lqmax = Lc + Ld 97 Lqmax3 = ((Lqmax+1)*(Lqmax+2)*(Lqmax+3))/6 98 99c Loop over the components of the "A" and "B" shells. 100 nqmx=0 101 do Iq = 0,Lqmax 102 do Jq = 0,Lqmax-Iq 103!DEC$ LOOP COUNT MAX=10, MIN=1 104 do Kq = 0,Lqmax-Iq-Jq 105 nqmx = max(IJK(Iq,Jq,Kq),nqmx) 106 enddo 107 enddo 108 enddo 109 nn = 0 110 do ma = 1,((La+1)*(La+2))/2 111 112 nc = la6 + ma 113 ia = Ixyz(1,nc) 114 ja = Ixyz(2,nc) 115 ka = Ixyz(3,nc) 116 117 118 if( canAB )then 119 mb_limit = ma 120 ij = (ma*(ma-1))/2 121 else 122 mb_limit = Lb2 123 ij = (ma-1)*Lb2 124 end if 125 126 do mb = 1,mb_limit 127 ij=ij+1 128 129 130 nc = lb6 + mb 131 ib = Ixyz(1,nc) 132 jb = Ixyz(2,nc) 133 kb = Ixyz(3,nc) 134 135c Sum across (Ip,Jp,Kp) for each value of (Iq,Jq,Kq). 136 call dcopy(mpq*nqmx,0d0,0,sum,1) 137 138 do Ip = 0,Ia+Ib 139 do Jp = 0,Ja+Jb 140 do Kp = 0,Ka+Kb 141 if(MPP.eq.1) then 142 E3(1) = Ep(1,1,0,Ip,Ia,Ib)* 143 & Ep(2,1,0,Jp,Ja,Jb)* 144 & Ep(3,1,0,Kp,Ka,Kb) 145 else 146 147c Define the product of the Hermite expansions coefficients for 148c overlap distribution "P". 149!DEC$ LOOP COUNT MAX=30, MIN=1 150 do mp = 1,MPP 151 E3(mp) = Ep(1,mp,0,Ip,Ia,Ib)* 152 & Ep(2,mp,0,Jp,Ja,Jb)* 153 & Ep(3,mp,0,Kp,Ka,Kb) 154 end do 155 endif 156 157 do Iq = 0,Lqmax 158 Ir = Ip + Iq 159 do Jq = 0,Lqmax-Iq 160 Jr = Jp + Jq 161 Lq = Iq + Jq -1 162 do Kq = 0,Lqmax-Iq-Jq 163 164 nq = IJK(Iq,Jq,Kq) 165 166 Kr = Kp + Kq 167 168 nr = IJK(Ir,Jr,Kr) 169 170c Include the factor of (-1)**(Iq+Jq+Kq). 171 Lq=Lq+1 172#if defined(GCC4) || defined(PGLINUX) 173 if(IAND(Lq,1).eq.1)then 174#else 175 if(AND(Lq,1).eq.1)then 176#endif 177 if(MPQ.eq.1) then 178!DEC$ LOOP COUNT MAX=30, MIN=1 179 do mp = 1,MPP 180 sum(1,nq) = sum(1,nq)-E3(mp)*R0(1,mp,nr) 181 end do 182 else 183 do mp = 1,MPP 184 zot=-E3(mp) 185!DEC$ LOOP COUNT MAX=30, MIN=1 186 do mq = 1,MPQ 187 sum(mq,nq) = sum(mq,nq)+zot*R0(mq,mp,nr) 188 end do 189 end do 190 endif 191 else 192 if(MPQ.eq.1) then 193!DEC$ LOOP COUNT MAX=30, MIN=1 194 do mp = 1,MPP 195 sum(1,nq) = sum(1,nq)+E3(mp)*R0(1,mp,nr) 196 end do 197 else 198 do mp = 1,MPP 199 zot=E3(mp) 200!DEC$ LOOP COUNT MAX=30, MIN=1 201 do mq = 1,MPQ 202 sum(mq,nq) = sum(mq,nq)+zot*R0(mq,mp,nr) 203 end do 204 end do 205 endif 206 end if 207 208 end do 209 end do 210 end do 211 212 end do 213 end do 214 end do 215 216c Loop over the components of the "C" and "D" shells. 217 218 do mc = 1,Lc2 219 220 nc = lc6 + mc 221 ic = Ixyz(1,nc) 222 jc = Ixyz(2,nc) 223 kc = Ixyz(3,nc) 224 225 if( canCD ) md_limit = mc 226 227 228 if( canAB )then 229 kl = (mc*(mc-1))/2 230 else 231 kl = (mc-1)*Ld2 232 end if 233 do md = 1,md_limit 234 235 if( canPQ )then 236 kl=kl+1 237 if( kl.gt.ij ) go to 480 238 end if 239 240 nc = ld6 + md 241 id = Ixyz(1,nc) 242 jd = Ixyz(2,nc) 243 kd = Ixyz(3,nc) 244 245 nn = nn + 1 246 247c Sum across (Iq,Jq,Kq). 248c Define the product of the Hermite expansion coefficients for 249c overlap distribution "Q" and calculate eri 250 erinn=0d0 251 252 if(MPQ.eq.1) then 253 do Iq = 0,Ic+Id 254 do Jq = 0,Jc+Jd 255!DEC$ LOOP COUNT MAX=30, MIN=1 256 do Kq = 0,Kc+Kd 257 nq = IJK(Iq,Jq,Kq) 258 ERInn = ERInn + 259 & Eq(1,1,0,Iq,Ic,Id)* 260 & Eq(2,1,0,Jq,Jc,Jd)* 261 & Eq(3,1,0,Kq,Kc,Kd)*sum(1,nq) 262 enddo 263 enddo 264 enddo 265 else 266 iqfin=0 267 iqlast=ic+id 268 if(iqlast.gt.0) then 269 iqfin=-1 270#if defined(GCC4) || defined(PGLINUX) 271 if(iand(iqlast,1).eq.0) then 272#else 273 if(and(iqlast,1).eq.0) then 274#endif 275 iqfin=ic+id 276 iqlast=iqlast-2 277 endif 278 do Iq = 0,iqlast,2 279 iq1=iq+1 280 do Jq = 0,Jc+Jd 281 do Kq = 0,Kc+Kd 282 nq = IJK(Iq,Jq,Kq) 283 nq1 = IJK(Iq1,Jq,Kq) 284!DEC$ LOOP COUNT MAX=30, MIN=1 285 do mq = 1,MPQ 286 ERInn = ERInn + 287 & (Eq(1,mq,0,Iq,Ic,Id)*sum(mq,nq)+ 288 & Eq(1,mq,0,Iq1,Ic,Id)*sum(mq,nq1))* 289 & Eq(2,mq,0,Jq,Jc,Jd)* 290 & Eq(3,mq,0,Kq,Kc,Kd) 291 end do 292 end do 293 enddo 294 enddo 295 endif 296 if(iqfin.ne.-1) then 297 do Jq = 0,Jc+Jd 298 do Kq = 0,Kc+Kd 299 nq = IJK(Iqfin,Jq,Kq) 300!DEC$ LOOP COUNT MAX=30, MIN=1 301 do mq = 1,MPQ 302 ERInn = ERInn + 303 & Eq(1,mq,0,Iqfin,Ic,Id)* 304 & Eq(2,mq,0,Jq,Jc,Jd)* 305 & Eq(3,mq,0,Kq,Kc,Kd)*sum(mq,nq) 306 end do 307 end do 308 enddo 309 endif 310 endif 311 eri(nn)=erinn 312 313 end do ! md 314 480 continue 315 end do ! mc 316 end do ! mb 317 end do ! ma 318 319 end 320