1c $Id$ 2 Subroutine hf2(Axyz,Aprims,Acoefs,NPA,NCA,La, 3 & Bxyz,Bprims,Bcoefs,NPB,NCB,Lb, 4 & Cxyz,Cprims,Ccoefs,NPC,NCC,Lc, 5 & Dxyz,Dprims,Dcoefs,NPD,NCD,Ld, 6 & bERI,Nints,canAB,canCD,canPQ, 7 & DryRun,W0,maxW0) 8 9 Implicit none 10#include "stdio.fh" 11#include "mafdecls.fh" 12#include "errquit.fh" 13#include "case.fh" 14 15 Integer NPA,NCA,La 16 Integer NPB,NCB,Lb 17 Integer NPC,NCC,Lc 18 Integer NPD,NCD,Ld 19 Integer Nints,maxW0 20 21 Logical canAB,canCD,canPQ 22 23 Logical GenCon,DryRun 24 25c--> Cartesian Coordinates, Primitives & Contraction Coefficients 26 27 Double Precision Axyz(3),Aprims(NPA),Acoefs(NPA,NCA) 28 Double Precision Bxyz(3),Bprims(NPB),Bcoefs(NPB,NCB) 29 Double Precision Cxyz(3),Cprims(NPC),Ccoefs(NPC,NCC) 30 Double Precision Dxyz(3),Dprims(NPD),Dcoefs(NPD,NCD) 31c 32c--> Block of Electron Repulsion Integrals 33 34 Double Precision bERI(Nints) 35 36c--> Scratch Space 37 38 Double Precision W0(maxW0) 39 40c--> Local variables 41 42 Integer NCP,NCQ,MXD,nd 43 Integer La2,Lb2,Lc2,Ld2 44 Integer lp,lq,lr,lp3,lq3,lr3,lpq3 45 Integer i_ALPHAp,i_IPAIRp,i_ESp,i_Ep 46 Integer i_ALPHAq,i_IPAIRq,i_ESq,i_Eq 47 Integer i_left,i_right,i_pf 48 Integer MaxMem,MaxAvail,MemPerPair,MaxPairs,NPass 49 Integer maxpp,maxpq,ipp,ipq 50 Integer MPP,MPQ,MPR,NPP,NPQ,NPR 51 Integer lszp,lszq,ninti 52 Integer i_R0,i_IJK,i_P,i_Q,i_PQ,i_ff,i_Rj 53 Integer i_E3,i_t1,i_t2,i_t3,i_t4,i_eri,i_top 54 Integer i,j,istep 55 integer i_ffscr,l_ffscr,i_rscr,l_rscr,nvals 56 57 Logical doit_PQ 58 59c 60c Compute 4-ctr electron repulsion integrals (ERI) for four shells of 61c contracted Gaussian functions. 62c 63c****************************************************************************** 64 65 MXD = 0 66 istep = MA_sizeof(MT_DBL,1,MT_INT) 67 if (istep .gt. 2) call errquit( 68 & 'hf2:Too many integers per real for Ipair addressing',911, 69 & INT_ERR) 70 istep = 2/istep 71 72c Determine whether general or segmented contraction is used. 73 74 NCP = NCA*NCB 75 NCQ = NCC*NCD 76 77 GenCon = (NCP.ne.1) .or. (NCQ.ne.1) 78 79c Define the angular momentum of the overlap distributions. 80 81 Lp = La + Lb 82 Lq = Lc + Ld 83 Lr = Lp + Lq 84 85c Increment "Lr" to account for the order of differentiation. 86 87 Lr = Lr + MXD 88 89c Define the accumulated number of angular momentum functions <= Lr. 90 91 Lp3 = ((Lp+1)*(Lp+2)*(Lp+3))/6 92 Lq3 = ((Lq+1)*(Lq+2)*(Lq+3))/6 93 Lpq3 = max(Lp3,Lq3) 94 Lr3 = ((Lr+1)*(Lr+2)*(Lr+3))/6 95 96c Define the prefactor of the overlap distribution "P". 97 98c Assign pointers to scratch space. 99 100 i_ALPHAp = 1 101 i_IPAIRp = i_ALPHAp + 2*(NPA*NPB) 102 i_left = i_IPAIRp + 2*(NPA*NPB) - 1 103 i_ESp = (maxW0+1) - 3*(NPA*NPB) 104 i_right = i_ESp 105 106 MaxMem = 1 ! take care of compiler warnings 107 if (DryRun) then 108 MaxMem = i_left + (maxW0 - (i_right-1)) 109 NPP = NPA*NPB 110 else 111 if (i_left.ge.i_right) then 112 write(LuOut,*) 'HF2: Insufficient scratch space.' 113 write(LuOut,*) ' needed ',i_left + (maxW0-(i_right-1)) 114 write(LuOut,*) ' allocated ',maxW0 115 write(LuOut,*) 'From the left ' 116 write(LuOut,*) 'ALPHAp: ',i_ALPHAp 117 write(LuOut,*) 'IPAIRp: ',i_IPAIRp 118 write(LuOut,*) 'From the right ' 119 write(LuOut,*) 'ESp : ',i_ESp 120 call errquit('hf2: Not enough scratch space for hfset AB',911, 121 & MEM_ERR) 122 end if 123 call hfset(Axyz,Aprims,Acoefs,NPA,NCA, 124 & Bxyz,Bprims,Bcoefs,NPB,NCB, 125 & GenCon,W0(i_ALPHAp),W0(i_IPAIRp),W0(i_ESp),NPP) 126 end if 127 128c Define the prefactor of the overlap distribution "Q". 129 130c Assign pointers to scratch space. 131 132 i_ALPHAq = i_IPAIRp + 2*(NPA*NPB) 133 i_IPAIRq = i_ALPHAq + 2*(NPC*NPD) 134 i_left = i_IPAIRq + 2*(NPC*NPD) - 1 135 i_ESq = i_right - 3*(NPC*NPD) 136 i_right = i_ESq 137 138 if( DryRun )then 139 MaxMem = max( MaxMem, i_left + (maxW0 - (i_right-1)) ) 140 NPQ = NPC*NPD 141 else 142 if (i_left.ge.i_right) then 143 write(LuOut,*) 'HF2: Insufficient scratch space.' 144 write(LuOut,*) ' needed ',i_left + (maxW0-(i_right-1)) 145 write(LuOut,*) ' allocated ',maxW0 146 write(LuOut,*) 'From the left ' 147 write(LuOut,*) 'ALPHAp: ',i_ALPHAp 148 write(LuOut,*) 'IPAIRp: ',i_IPAIRp 149 write(LuOut,*) 'ALPHAq: ',i_ALPHAq 150 write(LuOut,*) 'IPAIRq: ',i_IPAIRq 151 write(LuOut,*) 'From the right ' 152 write(LuOut,*) 'ESp : ',i_ESp 153 write(LuOut,*) 'ESq : ',i_ESq 154 call errquit('hf2: Not enough scratch space for hfset CD',911, 155 & MEM_ERR) 156 end if 157 call hfset(Cxyz,Cprims,Ccoefs,NPC,NCC, 158 & Dxyz,Dprims,Dcoefs,NPD,NCD, 159 & GenCon,W0(i_ALPHAq),W0(i_IPAIRq),W0(i_ESq),NPQ) 160 end if 161 162c Zero out integral block. Return if NPR = 0, i.e. no integrals 163 164 NPR = NPP*NPQ 165 if (NPR .eq. 0) then 166! if (.not.DryRun) call dfill(Nints,0.0d00,beri,1) 167 if (.not.DryRun) call dcopy(Nints,0d0,0,beri,1) 168 go to 99 169 endif 170 171c Define cartesian components for shells, bra and ket length and 172c number of integrals. 173 174 La2 = (la+1)*(la+2)/2 175 Lb2 = (lb+1)*(lb+2)/2 176 lszp = La2*Lb2*NCP 177 178 Lc2 = (lc+1)*(lc+2)/2 179 Ld2 = (ld+1)*(ld+2)/2 180 lszq = Lc2*Ld2*NCQ 181 182 ninti = lszp*lszq 183 if (.not.DryRun) call dfill(Ninti,0.0d00,beri,1) 184 185c Define the Hermite linear expansion coefficients. 186 187c Assign pointers to scratch space. 188 189 i_Ep = i_IPAIRq + 2*(NPC*NPD) 190 i_Eq = i_Ep + 3*NPP*(MXD+1)*((Lp+1)*(La+1)*(Lb+1)) 191 i_pf = i_Eq + 3*NPQ*(MXD+1)*((Lq+1)*(Lc+1)*(Ld+1)) 192 i_left = i_pf + 2*max(NPP,NPQ) - 1 193 194 if( DryRun )then 195 MaxMem = max( MaxMem, i_left + (maxW0 - (i_right-1)) ) 196 else 197 if (i_left.ge.i_right) then 198 write(LuOut,*) 'HF2: Insufficient scratch space.' 199 write(LuOut,*) ' needed ',i_left + (maxW0-(i_right-1)) 200 write(LuOut,*) ' allocated ',maxW0 201 write(LuOut,*) 'From the left ' 202 write(LuOut,*) 'ALPHAp: ',i_ALPHAp 203 write(LuOut,*) 'IPAIRp: ',i_IPAIRp 204 write(LuOut,*) 'ALPHAq: ',i_ALPHAq 205 write(LuOut,*) 'IPAIRq: ',i_IPAIRq 206 write(LuOut,*) 'Ep : ',i_Ep 207 write(LuOut,*) 'Eq : ',i_Eq 208 write(LuOut,*) 'pf : ',i_pf 209 write(LuOut,*) 'From the right ' 210 write(LuOut,*) 'ESp : ',i_ESp 211 write(LuOut,*) 'ESq : ',i_ESq 212 call errquit('hf2: Not enough scratch space for hfmke',911, 213 & MEM_ERR) 214 end if 215 216 do nd = 0,MXD 217 call hfmke(Axyz,Bxyz,W0(i_ALPHAp),W0(i_ESp),W0(i_Ep),W0(i_pf), 218 & nd,NPP,MXD,La,Lb) 219 call hfmke(Cxyz,Dxyz,W0(i_ALPHAq),W0(i_ESq),W0(i_Eq),W0(i_pf), 220 & nd,NPQ,MXD,Lc,Ld) 221 end do 222 223 end if 224 225c Evaluate the auxiliary function integrals. 226 227c Determine memory requirements and assign pointers to scratch space. 228c If insufficient memory, do multipassing. 229 230 i_IJK = i_Eq + 3*NPQ*(MXD+1)*((Lq+1)*(Lc+1)*(Ld+1)) 231 i_R0 = i_IJK + (Lr+1)**3 232 233 MaxAvail = maxW0-i_R0-3*(NPP+NPQ)+1 234 MemPerPair = (Lr+2)*Lr3 + 5 235 if (MaxAvail-MemPerPair .le. 0) then 236 write(LuOut,*) 'HF2: Insufficient scratch space.' 237 write(LuOut,*) ' needed ',i_R0+3*(NPP+NPQ) 238 & +MemPerPair*NPR 239 write(LuOut,*) ' allocated ',maxW0 240 write(LuOut,*) 'ALPHAp: ',i_ALPHAp 241 write(LuOut,*) 'IPAIRp: ',i_IPAIRp 242 write(LuOut,*) 'ALPHAq: ',i_ALPHAq 243 write(LuOut,*) 'IPAIRq: ',i_IPAIRq 244 write(LuOut,*) 'Ep : ',i_Ep 245 write(LuOut,*) 'Eq : ',i_Eq 246 write(LuOut,*) 'R0 : ',i_R0 247 write(LuOut,*) 'IJK : ',i_IJK 248 call errquit('hf2: Not enough scratch space for hf2mkr',911, 249 & MEM_ERR) 250 end if 251 252 MaxPairs = Maxavail/MemPerPair 253 254 if (MaxPairs .ge. NPR) then 255 maxpp = NPP 256 maxpq = NPQ 257 else if (NPP .ge. NPQ) then 258 if (MaxPairs .ge. NPP) then 259 maxpp = NPP 260 maxpq = MaxPairs/NPP 261 else if (MaxPairs .ge. NPQ) then 262 maxpq = NPQ 263 maxpp = MaxPairs/NPQ 264 else 265 NPass = (NPP-1)/MaxPairs+1 266 maxpp = min((NPP-1)/NPass+1,MaxPairs) 267 maxpq = 1 268 end if 269 else 270 if (MaxPairs .ge. NPQ) then 271 maxpq = NPQ 272 maxpp = MaxPairs/NPQ 273 else if (MaxPairs .ge. NPP) then 274 maxpp = NPP 275 maxpq = MaxPairs/NPP 276 else 277 NPass = (NPP-1)/MaxPairs+1 278 maxpp = min((NPP-1)/NPass+1,MaxPairs) 279 maxpq = 1 280 end if 281 end if 282 283 do ipp = 0,NPP-1,maxpp 284 MPP = min(maxpp,NPP-ipp) 285 do ipq = 0,NPQ-1,maxpq 286 MPQ = min(maxpq,NPQ-ipq) 287 MPR = MPP*MPQ 288 i_P = i_R0 + MPR*Lr3 289 i_Q = i_P + 3*MPP 290 i_PQ = i_Q + 3*MPQ 291 i_ff = i_PQ + MPR*3 292 i_Rj = i_ff + 2*MPR 293 i_top = i_Rj + MPR*(Lr+1)*Lr3 - 1 294 295 if (Lp.eq.Lq) then 296 doit_PQ = MPP.ge.MPQ 297 else if (abs((MPP-MPQ)).ge.4) then 298 doit_PQ = MPP.ge.MPQ 299 else if (Lp.lt.Lq) then 300 doit_PQ = .true. 301 else 302 doit_PQ = .false. 303 endif 304 305 if( DryRun )then 306 MaxMem = max( MaxMem, i_top ) 307 else 308 if(doscreen) then 309 nvals=2*MPP*MPQ 310 if (.not.MA_Push_get(mt_dbl,nvals,'ffscr', 311 C l_ffscr,i_ffscr)) call errquit( 312 E 'hf2: ran out of MA memory',0,MEM_ERR) 313 nvals=2*MPP*MPQ*lr3*(lr+1) 314 if (.not.MA_Push_get(mt_dbl,nvals,'rscr', 315 C l_rscr,i_rscr)) call errquit( 316 E 'hf2: ran out of MA memory',1,MEM_ERR) 317 else 318 nvals=1 319 i_rscr=0 320 i_ffscr=0 321 endif 322 if (doit_PQ) then 323 call hf2mkr(Axyz,Bxyz,Cxyz,Dxyz, 324 & W0(i_ALPHAp+2*ipp),W0(i_ALPHAq+2*ipq), 325 & W0(i_R0),W0(i_IJK),W0(i_P),W0(i_Q),W0(i_PQ), 326 & W0(i_ff),W0(i_Rj), 327 & dbl_mb(i_ffscr),dbl_mb(i_rscr), 328 & MPP,MPQ,Lr,Lr3) 329 else 330 call hf2mkr(Cxyz,Dxyz,Axyz,Bxyz, 331 & W0(i_ALPHAq+2*ipq),W0(i_ALPHAp+2*ipp), 332 & W0(i_R0),W0(i_IJK),W0(i_Q),W0(i_P),W0(i_PQ), 333 & W0(i_ff),W0(i_Rj), 334 & dbl_mb(i_ffscr),dbl_mb(i_rscr), 335 & MPQ,MPP,Lr,Lr3) 336 end if 337 if(doscreen) then 338 if (.not.ma_chop_stack(l_ffscr)) call errquit( 339 E 'hf2: machopstack failed',0,MEM_ERR) 340 endif 341 end if 342 343c Compute the ERI. 344 345c Assign pointers to scratch space. 346 347 i_eri = i_P 348 i_t2 = i_eri + ninti 349 if (gencon) then 350 if (doit_PQ) then 351 i_E3 = i_t2 + MPQ*Lq3*max(MPP,NCP) 352 i_t1 = i_E3 353 i_t3 = i_E3 + MPQ*Lq3 354 i_t4 = i_t3 + MPQ 355 i_top = max(i_t1+MPQ*Lq3*NPB*NCA,i_E3+MPP,i_t4+NPD)-1 356 else 357 i_E3 = i_t2 + MPP*Lp3*max(MPQ,NCQ) 358 i_t1 = i_E3 359 i_t3 = i_E3 + MPP*Lp3 360 i_t4 = i_t3 + MPP 361 i_top = max(i_t1+MPP*Lp3*NPD*NCC,i_E3+MPQ,i_t4+NPB)-1 362 end if 363 else 364 i_E3 = i_t2 + MPR*Lpq3 365 i_top = i_E3 + max(MPP,MPQ)-1 366 end if 367 if( DryRun )then 368 MaxMem = max( MaxMem, i_top ) 369 else if (i_top.gt.maxW0) then 370 write(LuOut,*) 'HF2: Insufficient scratch space.' 371 write(LuOut,*) ' needed ',i_top 372 write(LuOut,*) ' allocated ',maxW0 373 write(LuOut,*) 'ALPHAp: ',i_ALPHAp 374 write(LuOut,*) 'IPAIRp: ',i_IPAIRp 375 write(LuOut,*) 'ALPHAq: ',i_ALPHAq 376 write(LuOut,*) 'IPAIRq: ',i_IPAIRq 377 write(LuOut,*) 'Ep : ',i_Ep 378 write(LuOut,*) 'Eq : ',i_Eq 379 write(LuOut,*) 'IJK : ',i_IJK 380 write(LuOut,*) 'R0 : ',i_R0 381 write(LuOut,*) 'eri : ',i_eri 382 write(LuOut,*) 't2 : ',i_t2 383 write(LuOut,*) 'E3 : ',i_E3 384 if (gencon) then 385 write(LuOut,*) 't1 : ',i_t1 386 write(LuOut,*) 't3 : ',i_t3 387 write(LuOut,*) 't4 : ',i_t4 388 end if 389 call errquit('hf2: Not enough scratch space for hferi',911, 390 & MEM_ERR) 391 else if (doit_PQ) then 392 if (GenCon) then 393 call hferi_gen(W0(i_Ep+3*ipp),W0(i_Eq+3*ipq), 394 & W0(i_IPAIRp+istep*ipp),W0(i_IPAIRq+istep*ipq), 395 & W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3), 396 & W0(i_t1),W0(i_t2),W0(i_t3),W0(i_t4), 397 & MPP,MPQ,NPP,NPQ,La,Lb,Lc,Ld,La2,Lb2,Lc2,Ld2, 398 & Lq,Lq3,Lr,Acoefs,Bcoefs,Ccoefs,Dcoefs, 399 & NPA,NPB,NPC,NPD,NCA,NCB,NCC,NCD, 400 & MXD,canAB,canCD,canPQ) 401 else 402 call hferi(W0(i_Ep+3*ipp),W0(i_Eq+3*ipq), 403 & W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3),W0(i_t2), 404 & MPP,MPQ,NPP,NPQ,Nints,La,Lb,Lc,Ld,Lr,MXD, 405 & canAB,canCD,canPQ) 406 end if 407 call daxpy(ninti,1.0D0,W0(i_eri),1,beri,1) 408 else 409 if (GenCon) then 410 call hferi_gen(W0(i_Eq+3*ipq),W0(i_Ep+3*ipp), 411 & W0(i_IPAIRq+istep*ipq),W0(i_IPAIRp+istep*ipp), 412 & W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3), 413 & W0(i_t1),W0(i_t2),W0(i_t3),W0(i_t4), 414 & MPQ,MPP,NPQ,NPP,Lc,Ld,La,Lb,Lc2,Ld2,La2,Lb2, 415 & Lp,Lp3,Lr,Ccoefs,Dcoefs,Acoefs,Bcoefs, 416 & NPC,NPD,NPA,NPB,NCC,NCD,NCA,NCB, 417 & MXD,canCD,canAB,canPQ) 418 else 419 call hferi(W0(i_Eq+3*ipq),W0(i_Ep+3*ipp), 420 & W0(i_R0),W0(i_IJK),W0(i_eri),W0(i_E3),W0(i_t2), 421 & MPQ,MPP,NPQ,NPP,Ninti,Lc,Ld,La,Lb,Lr,MXD, 422 & canCD,canAB,canPQ) 423 end if 424 j = 1 425 do i = 0,lszp-1 426 call daxpy(lszq,1.0d0,W0(i_eri+i),lszp,beri(j),1) 427 j = j+lszq 428 end do 429 end if 430 431 end do 432 end do 433 434c Return the maximum amount of scratch space required by a "dry run". 435 436 99 if( DryRun ) maxW0 = MaxMem 437c 438 end 439 440