1c $Id$ 2 Subroutine hferi_gen(Ep,Eq,PairP,PairQ,R0,IJK,ERI,E3,t1,t2,t3,t4, 3 & MPP,MPQ,NPP,NPQ,La,Lb,Lc,Ld,La2,Lb2,Lc2,Ld2,Lqmax,Lqmax3,Lr, 4 & Acoefs,Bcoefs,Ccoefs,Dcoefs,NPA,NPB,NPC,NPD,NCA,NCB,NCC,NCD, 5 & MXD,canAB,canCD,canPQ) 6 7 Implicit none 8#include "errquit.fh" 9 10 Logical canAB,canCD,canPQ 11 12c--> Number of pairs in pass, total number 13 14 Integer MPP,MPQ, NPP,NPQ 15 16c--> Angular momenta and number of cartesian components 17 18 Integer La,Lb,Lc,Ld, La2,Lb2,Lc2,Ld2, Lqmax,Lqmax3,Lr 19 20c--> Derivative order 21 22 Integer MXD 23 24c--> Number of primitives, number of contracted functions 25 26 Integer NPA,NPB,NPC,NPD, NCA,NCB,NCC,NCD 27 28c--> Hermite Linear Expansion Coefficients 29 30 Double precision Ep(3,NPP,0:MXD,0:(La+Lb),0:La,0:Lb) 31 Double precision Eq(3,NPQ,0:MXD,0:(Lc+Ld),0:Lc,0:Ld) 32 33c--> Pair index arrays 34 35 Integer PairP(2,NPP) 36 Integer PairQ(2,NPQ) 37 38c--> Auxiliary Function Integrals & Index 39 40 Double precision R0(MPQ,MPP,*) 41 Integer IJK(0:Lr,0:Lr,0:Lr) 42 43c--> Contraction coefficients 44 45 Double precision Acoefs(NPA,NCA) 46 Double precision Bcoefs(NPB,NCB) 47 Double precision Ccoefs(NPC,NCC) 48 Double precision Dcoefs(NPD,NCD) 49 50c--> ERI 51 52 Double precision ERI(Ld2,NCD,Lc2,NCC,Lb2*NCB,La2*NCA) 53 54c--> Scratch Space 55 56 Double precision E3(*) 57 Double precision t1(MPQ,Lqmax3,NPB,NCA) ! 1st 1/4 trans 58 Double precision t2(MPQ,Lqmax3,MPP) ! 2nd 1/4 trans, primitives 59 Double precision t3(MPQ) ! 3rd 1/4 trans 60 Double precision t4(NPD) ! 4th 1/4 trans 61 Integer Nxyz(3) 62 63c--> Local variables 64 65 Integer Lq,ninti 66 Integer ma,mb,mc,md, mp,mq, nq,nr 67 Integer Ia,Ja,Ka, Ib,Jb,Kb, Ic,Jc,Kc, Id,Jd,Kd 68 Integer Ip,Jp,Kp, Iq,Jq,Kq, Ir,Jr,Kr 69 Integer ipa,ipb,ipc,ipd, ica,icb,icc,icd 70 Integer ie,icab,NCAB 71 Integer ind_ca,ind_cb 72 Double precision ca,cb 73 74c 75c Compute electron repulsion integrals (ERI). 76c 77c Formula: 78c 79c __ 80c \ Ic,Id;n10 Jc,Jd;n11 Kc,Kd;n12 81c ERI = / Ex Ey Ez SUM 82c -- Iq Jq Kq Iq,Jq,Kq 83c Iq,Jq,Kq 84c 85c __ 86c \ Lq Ia,Ib;n7 Ja,Jb;n8 Ka,Kb;n9 87c SUM = / (-1) Ex Ey Ez R 88c Iq,Jq,Kq -- Ip Jp Kp Ir,Jr,Kr 89c Ip,Jp,Kp 90c 91c Ir = (Ip+n1) + (Iq+n4) 92c where Jr = (Jp+n2) + (Jq+n5) 93c Kr = (Kp+n3) + (Kq+n6) 94c 95c and Lq = (Iq+n4) + (Jq+n5) + (Kq+n6) 96c 97c N.B. For simple ERI (i.e., no derivative integrals) n[1-12] = 0! 98c 99c****************************************************************************** 100c General case: [ab|cd] 101 102 if (canAB .or. canCD .or. canPQ) 103 & call errquit ('hferi_gen:can''t do canonical integrals',911, 104 & INT_ERR) 105 106c Initialize the block of ERIs. 107 108 NCAB = NCA*NCB 109 ninti = La2*Lb2*Lc2*Ld2*NCA*NCB*NCC*NCD 110 call dfill (ninti,0.0d0,eri,1) 111 112c Loop over the components of the "A" and "B" shells. 113 114 do ma = 1,La2 115 116 call getNxyz(La,ma,Nxyz) 117 Ia = Nxyz(1) 118 Ja = Nxyz(2) 119 Ka = Nxyz(3) 120 121 do mb = 1,Lb2 122 123 call getNxyz(Lb,mb,Nxyz) 124 Ib = Nxyz(1) 125 Jb = Nxyz(2) 126 Kb = Nxyz(3) 127 128c Sum across (Ip,Jp,Kp) for each value of (Iq,Jq,Kq). 129 130 call dfill(MPP*MPQ*Lqmax3,0.0d00,t2,1) 131 132 do Ip = 0,Ia+Ib 133 do Jp = 0,Ja+Jb 134 do Kp = 0,Ka+Kb 135 136c Define the product of the Hermite expansion coefficients for 137c overlap distribution "P". 138 139 do mp = 1,MPP 140 E3(mp) = Ep(1,mp,0,Ip,Ia,Ib) 141 & *Ep(2,mp,0,Jp,Ja,Jb) 142 & *Ep(3,mp,0,Kp,Ka,Kb) 143 end do 144 145 do Iq = 0,Lqmax 146 do Jq = 0,Lqmax-Iq 147 do Kq = 0,Lqmax-Iq-Jq 148 149 nq = IJK(Iq,Jq,Kq) 150 Ir = Ip + Iq 151 Jr = Jp + Jq 152 Kr = Kp + Kq 153 nr = IJK(Ir,Jr,Kr) 154 155c Include the factor of (-1)**(Iq+Jq+Kq). 156c Sum over Hermite functions 157 158 Lq = Iq + Jq + Kq 159 if( mod(Lq,2).eq.0 )then 160 do mp = 1,MPP 161 do mq = 1,MPQ 162 t2(mq,nq,mp) = t2(mq,nq,mp) 163 & + E3(mp)*R0(mq,mp,nr) 164 end do 165 end do 166 else 167 do mp = 1,MPP 168 do mq = 1,MPQ 169 t2(mq,nq,mp) = t2(mq,nq,mp) 170 & - E3(mp)*R0(mq,mp,nr) 171 end do 172 end do 173 end if 174 175 end do 176 end do 177 end do 178 179 end do 180 end do 181 end do 182 183c Contract over shell a 184 185 call dfill(MPQ*Lqmax3*NPB*NCA,0.0d0,t1,1) 186 do ica = 1,NCA 187 do mp = 1,MPP 188 ipa = PairP(1,mp) 189 ipb = PairP(2,mp) 190 ca = Acoefs(ipa,ica) 191 do nq = 1,Lqmax3 192 do mq = 1,MPQ 193 t1(mq,nq,ipb,ica) = t1(mq,nq,ipb,ica) + 194 & t2(mq,nq,mp)*ca 195 end do 196 end do 197 end do 198 end do 199 200c Contract over shell b 201 202 call dfill(MPQ*Lqmax3*NCAB,0.0d0,t2,1) 203 icab = 0 204 do ica = 1,NCA 205 do icb = 1,NCB 206 icab = icab+1 207 do ipb = 1,NPB 208 cb = Bcoefs(ipb,icb) 209 do nq = 1,Lqmax3 210 do mq = 1,MPQ 211 t2(mq,nq,icab) = t2(mq,nq,icab) 212 & +t1(mq,nq,ipb,ica)*cb 213 end do 214 end do 215 end do 216 end do 217 end do 218 219 220c Loop over the components of the "C" and "D" shells. 221 222 do mc = 1,Lc2 223 224 call getNxyz(Lc,mc,Nxyz) 225 Ic = Nxyz(1) 226 Jc = Nxyz(2) 227 Kc = Nxyz(3) 228 229 do md = 1,Ld2 230 231 call getNxyz(Ld,md,Nxyz) 232 Id = Nxyz(1) 233 Jd = Nxyz(2) 234 Kd = Nxyz(3) 235 236c Define the product of the Hermite expansion coefficients for 237c overlap distribution "Q". 238 239 ie = 0 240 do Iq = 0,Ic+Id 241 do Jq = 0,Jc+Jd 242 do Kq = 0,Kc+Kd 243 do mq = 1,MPQ 244 ie = ie+1 245 E3(ie) = Eq(1,mq,0,Iq,Ic,Id) 246 & *Eq(2,mq,0,Jq,Jc,Jd) 247 & *Eq(3,mq,0,Kq,Kc,Kd) 248 end do 249 end do 250 end do 251 end do 252 253c Sum across (Iq,Jq,Kq). 254 255 icab = 0 256 ind_ca = ma 257 do ica = 1,NCA 258 ind_cb = mb 259 do icb = 1,NCB 260 call dfill (MPQ,0.0d0,t3,1) 261 icab = icab+1 262 ie = 0 263 do Iq = 0,Ic+Id 264 do Jq = 0,Jc+Jd 265 do Kq = 0,Kc+Kd 266 267 nq = IJK(Iq,Jq,Kq) 268 269c Sum over Hermite functions to give integrals 270 271 do mq = 1,MPQ 272 ie = ie+1 273 t3(mq) = t3(mq) + E3(ie)*t2(mq,nq,icab) 274 end do 275 276 end do 277 end do 278 end do 279 280c Contract over C shell 281 282 do icc = 1,NCC 283 call dfill (NPD,0.0d0,t4,1) 284 do mq = 1,MPQ 285 ipc = PairQ(1,mq) 286 ipd = PairQ(2,mq) 287 t4(ipd) = t4(ipd) + t3(mq)*Ccoefs(ipc,icc) 288 end do 289 290c Contract over D shell 291 292 do icd = 1,NCD 293 do ipd = 1,NPD 294 ERI(md,icd,mc,icc,ind_cb,ind_ca) = 295 & ERI(md,icd,mc,icc,ind_cb,ind_ca) + 296 & t4(ipd)*Dcoefs(ipd,icd) 297 end do 298 end do 299 end do 300 301 ind_cb = ind_cb+Lb2 302 end do ! icb 303 ind_ca = ind_ca+La2 304 end do ! ica 305 306 end do 307 end do 308 309 end do 310 end do 311 312 end 313