1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck hr2drv */ 20 SUBROUTINE HR2DRV(HERINT,INDHER,COOR12,COOR34,EXP12,EXP34,FAC12, 21 & FAC34,NTUV,WORK,LWORK,JMAX,MAXDER,NUABCD,IPQ0X, 22 & IPQ0Y,IPQ0Z,NOINT,ONECEN,NUC1,NUC2,NUC12,NUC3, 23 & NUC4,NUC34,THRESH,IPRINT,SIGNT,IODDHR) 24C 25C T. Helgaker, Sep. 91 26C 27C gnrinf.h needed for PANAS correction factor, and CHIVAL,ERFEXP, 28C SRINTS for partition of r_12 (MCDFT) 29C 30#include "implicit.h" 31#include "priunit.h" 32#include "maxaqn.h" 33 LOGICAL NOINT, ONECEN 34 DIMENSION HERINT(NUABCD,NTUV), INDHER(*), IODDHR(*), 35 & COOR12(*), COOR34(*), EXP12(*), EXP34(*), FAC12(*), 36 & FAC34(*), SIGNT(*), WORK(LWORK) 37#include "gnrinf.h" 38#include "twosta.h" 39 LRJ000 = (JMAX + 1)*NUABCD 40 KRJ000 = 1 41 KPQX = KRJ000 + LRJ000 42 KPQY = KPQX + NUABCD 43 KPQZ = KPQY + NUABCD 44 KLAST = KPQZ + NUABCD 45 IF (KLAST .GT. LWORK) CALL STOPIT('HR2DRV',' ',KLAST,LWORK) 46 MWRJ00 = MAX(MWRJ00,LRJ000) 47 LWTOT = LWTOT + KLAST 48 MWTOT = MAX(MWTOT,LWTOT) 49 LWRK = LWORK - KLAST + 1 50 CALL HR2DR1(HERINT,INDHER,COOR12,COOR34,EXP12,EXP34,FAC12,FAC34, 51 & WORK(KRJ000),WORK(KPQX),WORK(KPQY),WORK(KPQZ), 52 & WORK(KLAST),LWRK,NTUV,JMAX,MAXDER,NUABCD,IPQ0X,IPQ0Y, 53 & IPQ0Z,NOINT,ONECEN,NUC1,NUC2,NUC12,NUC3,NUC4,NUC34, 54 & THRESH,IPRINT,SIGNT,IODDHR,PANAS,CHIVAL,ERFEXP, 55 & VLAMBDA,COMLAM,SRINTS) 56 LWTOT = LWTOT - KLAST 57 RETURN 58 END 59C /* Deck hr2dr1 */ 60 SUBROUTINE HR2DR1(HERINT,INDHER,COOR12,COOR34,EXP12,EXP34,FAC12, 61 & FAC34,RJ000,PQX,PQY,PQZ,WORK,LWORK,NTUV,JMAX, 62 & MAXDER,NUABCD,IPQ0X,IPQ0Y,IPQ0Z,NOINT,ONECEN, 63 & NUC1,NUC2,NUC12,NUC3,NUC4,NUC34,THRESH,IPRINT, 64 & SIGNT,IODDHR,PANAS,CHIVAL,ERFEXP, 65 & VLAMBDA,COMLAM,SRINTS) 66C Modified for r12 integrals (WK/UniKA/15-11-2002). 67#include "implicit.h" 68#include "priunit.h" 69#include "r12int.h" 70#include "drw2el.h" 71#include "maxaqn.h" 72 PARAMETER (D1 = 1.D0) 73 LOGICAL NOINT, ONECEN, SRINTS, COMLAM 74 LOGICAL ERFEXP(0:*) 75#include "twosta.h" 76C HERINT is declared 3-dimensional for r12 integrals (WK/UniKA/14-11-2002). 77 DIMENSION HERINT(NUABCD,NTUV,2), INDHER(*), IODDHR(*), 78 & RJ000(NUABCD,*), WORK(LWORK), 79 & PQX(NUABCD), PQY(NUABCD), PQZ(NUABCD), 80 & COOR12(*), COOR34(*), EXP12(*), EXP34(*), FAC12(*), 81 & FAC34(*), SIGNT(*) 82C 83 IF (R12EIN) THEN 84C 85C Gaussian-damped R12 integrals (WK/UniKA/15-11-2002). 86C ==================================================== 87C 88 LRJ000 = NUABCD * (JMAX + 1) 89 LHRINT = NUABCD * NTUV 90 KFACNT = 1 91 KALPHA = KFACNT 92 KHARGE = KALPHA + NUABCD 93 KHARGF = KHARGE + NUABCD 94 KHALPH = KHARGF + NUABCD 95 KHBETA = KHALPH + NUABCD 96 KHPRE1 = KHBETA + NUABCD 97 KHPRE2 = KHPRE1 + NUABCD 98 KHPRE3 = KHPRE2 + NUABCD 99 KHEXPP = KHPRE3 + NUABCD 100 KHEXPQ = KHEXPP + NUABCD 101 KRO000 = KHEXPQ + NUABCD 102 KHRIN3 = KRO000 + LRJ000 103 KHRIN4 = KHRIN3 + LHRINT 104 KHRIN5 = KHRIN4 + LHRINT 105 LSTEIN = KHRIN5 + LHRINT 106 LWKHER = LWORK - LSTEIN + 1 107 IF (LWKHER .LE. 0) CALL STOPIT('HR2DR1',' ',LSTEIN,LWORK) 108 IF (TKTIME) TIMSTR = SECOND() 109 CALL R00G(RJ000,COOR12,COOR34,EXP12,EXP34,FAC12,FAC34,PQX,PQY, 110 & PQZ,JMAX,NOINT,NUABCD,NUC1,NUC2,NUC12,NUC3,NUC4, 111 & NUC34,THRESH,ONECEN,IPRINT,IPQ0X,IPQ0Y,IPQ0Z,SIGNT, 112 & WORK(KFACNT),WORK(KHEXPP),WORK(KHEXPQ)) 113 CALL ERIIFC(NTUV,NUABCD,IPQ0X,IPQ0Y,IPQ0Z,JMAX) 114 CALL ERIHER(HERINT(1,1,1),HERINT(1,1,2),RJ000,PQX,INDHER, 115 & IODDHR,WORK(KFACNT),WORK(KALPHA),WORK(KHARGE), 116 & WORK(KHARGF),WORK(KHALPH),WORK(KHBETA), 117 & WORK(KHPRE1),WORK(KHPRE2),WORK(KHPRE3), 118 & WORK(KRO000),WORK(KHRIN3),WORK(KHRIN4), 119 & WORK(KHRIN5),WORK(KHEXPP),WORK(KHEXPQ), 120 & WORK(LSTEIN),LWKHER,IPRINT) 121 IF (TKTIME) THEN 122 TIMEND = SECOND() 123 TIME = TIMEND - TIMSTR 124 THERIX(JMAX) = THERIX(JMAX) + TIME 125 THERI = THERI + TIME 126 END IF 127C 128 ELSE 129C 130C Incomplete gamma function 131C ========================= 132C 133 IF (TKTIME) TIMSTR = SECOND() 134 CALL R000(RJ000,COOR12,COOR34,EXP12,EXP34,FAC12,FAC34,PQX,PQY, 135 & PQZ,JMAX,NOINT,NUABCD,NUC1,NUC2,NUC12,NUC3,NUC4, 136 & NUC34,THRESH,ONECEN,IPRINT,IPQ0X,IPQ0Y,IPQ0Z,SIGNT, 137 & PANAS,CHIVAL,ERFEXP,VLAMBDA,COMLAM,SRINTS,WORK,LWORK) 138 IF (TKTIME) THEN 139 TIMEND = SECOND() 140 TIME = TIMEND - TIMSTR 141 TR000X(JMAX) = TR000X(JMAX) + TIME 142 TR000 = TR000 + TIME 143 TIMSTR = TIMEND 144 END IF 145C 146C Hermite integrals 147C ================= 148C 149 IF (.NOT.NOINT) THEN 150C 151C Calculate integrals 152C 153 CALL HERI(HERINT,WORK,LWORK,RJ000,PQX,PQY,PQZ,INDHER,JMAX, 154 & MAXDER,NUABCD,NTUV,IPQ0X,IPQ0Y,IPQ0Z,IODDHR,IPRINT) 155C 156 IF (R12INT .OR. U12INT .OR. BPH2OO) THEN 157C 158C Calculate Hermite integrals Q(t,u,v) over r12 from 159C R(t,u,v) integrals over 1/r12 (WK/UniKA/14-11-2002). 160C 161 KWKAMA = 1 162 KWKAMB = KWKAMA + NUABCD 163 KWKLST = KWKAMB + NUABCD 164 IF (KWKLST .GT. LWORK) 165 & CALL STOPIT('HR2DR1','ABEQ52',KWKLST,LWORK) 166 CALL ABEQ52(HERINT(1,1,2),HERINT(1,1,1),WORK(KWKAMA), 167 & WORK(KWKAMB),PQX,PQY,PQZ,INDHER,JMAX, 168 & EXP12,EXP34,NUC12,NUC34,NUABCD, 169 & NTUV,IPQ0X,IPQ0Y,IPQ0Z,IODDHR,IPRINT) 170 END IF 171 IF (TKTIME) THEN 172 TIMEND = SECOND() 173 TIME = TIMEND - TIMSTR 174 THERIX(JMAX) = THERIX(JMAX) + TIME 175 THERI = THERI + TIME 176 TIMSTR = TIMEND 177 END IF 178 END IF 179 END IF 180C 181 RETURN 182 END 183C /* Deck r000 */ 184 SUBROUTINE R000(RJ000,COOR12,COOR34,EXP12,EXP34,FAC12,FAC34,PQX, 185 & PQY,PQZ,JMAX,NOINT,NUABCD,NUC1,NUC2,NUC12,NUC3, 186 & NUC4,NUC34,THRESH,ONECEN,IPRINT,IPQ0X,IPQ0Y,IPQ0Z, 187 & SIGNT,PANAS,CHIVAL,ERFEXP,VLAMBDA,COMLAM,SRINTS, 188 & WORK,LWORK) 189C 190C TUH 84 191C 192#include "implicit.h" 193#include "priunit.h" 194#include "iratdef.h" 195#include "maxaqn.h" 196#include "twosta.h" 197#include "r12int.h" 198 LOGICAL ONECEN, NOINT, SRINTS, COMLAM 199 LOGICAL ERFEXP(0:*) 200 DIMENSION RJ000(*), PQX(*), PQY(*), PQZ(*), COOR12(*), COOR34(*), 201 & EXP12 (*), EXP34(*), FAC12 (*), FAC34(*), SIGNT(3), 202 & WORK(LWORK) 203C----------------------------------------------------------------------- 204C ============================================= 205C Parameters used for memory-allocations 206C ============================================= 207 LXJ000 = (JMAX+1)*NUABCD 208 IF(ERFEXP(0)) THEN 209C Allocate extra memory for exponential factors 210C in modified two-electron operator. \jkp 211 NTALPHX = 2 212 ELSE 213 NTALPHX = 1 214 ENDIF 215 IF (SRINTS) THEN 216 NTALPHX = NTALPHX + 1 217 END IF 218C ============================================= 219C Allocations 220C ============================================= 221 KWVALU = 1 222 KALPHA = KWVALU + NTALPHX*NUABCD 223 KALPHJ = KALPHA + NTALPHX*NUABCD 224 KINDAD = KALPHJ + NTALPHX*NUABCD 225 KEJ000 = KINDAD + (NUABCD + 1)/IRAT 226 IF(ERFEXP(0)) THEN 227 KSJ000 = KEJ000 + LXJ000 228 ELSE 229 KSJ000 = KEJ000 230 ENDIF 231 IF(SRINTS) THEN 232 KLAST = KSJ000 + LXJ000 233 ELSE 234 KLAST = KSJ000 235 ENDIF 236 IF (KLAST .GT. LWORK) CALL STOPIT('R000',' ',KLAST,LWORK) 237 LWRK = LWORK - KLAST + 1 238C ============================================= 239 CALL R0001(RJ000,COOR12,COOR34,EXP12,EXP34,FAC12,FAC34,PQX,PQY, 240 & PQZ,JMAX,NOINT,NUABCD,NUC1,NUC2,NUC12,NUC3,NUC4,NUC34, 241 & THRESH,ONECEN,IPRINT,IPQ0X,IPQ0Y,IPQ0Z,SIGNT, 242 & WORK(KWVALU),WORK(KALPHJ),WORK(KALPHA),WORK(KINDAD), 243 & PANAS,CHIVAL,ERFEXP,VLAMBDA,COMLAM,NTALPHX, 244 & WORK(KEJ000),WORK(KSJ000),SRINTS,WORK(KLAST),LWRK) 245 RETURN 246 END 247C /* Deck r0001 */ 248 SUBROUTINE R0001(RJ000,COOR12,COOR34,EXP12,EXP34,FAC12,FAC34,PQX, 249 & PQY,PQZ,JMAX,NOINT,NUABCD,NUC1,NUC2,NUC12,NUC3, 250 & NUC4,NUC34,THRESH,ONECEN,IPRINT,IPQ0X,IPQ0Y,IPQ0Z, 251 & SIGNT,WVALU,TALPHA,SCLFAC,INDADR,PANAS,CHIVAL, 252 & ERFEXP,VLAMBDA,COMLAM,NTALPHX,EJ000,SJ000, 253 & SRINTS,WORK,LWORK) 254C 255C TUH 84 256C 257#include "implicit.h" 258#include "priunit.h" 259#include "iratdef.h" 260#include "maxaqn.h" 261#include "aovec.h" 262#include "codata.h" 263#include "drw2el.h" 264 PARAMETER (D0=0.0D0,D1=1.0D0,D2=2.0D0,D3=3.0D0,D5=5.0D0, 265 & DP25=0.25D0,D10=10.0D0,D18=18.0D0) 266 LOGICAL ONECEN, NOINT, SRINTS, COMLAM 267 LOGICAL ERFEXP(0:*) 268 DIMENSION RJ000(NUABCD,0:JMAX),EJ000(NUABCD,0:JMAX), 269 & SJ000(NUABCD,0:JMAX),PQX(NUABCD),PQY(NUABCD), 270 & PQZ(NUABCD),COOR12(NUC1*NUC2,3),COOR34(NUC3*NUC4,3), 271 & EXP12 (*), EXP34(*), FAC12 (*), FAC34(*), SIGNT(3), 272 & WVALU(NUABCD,NTALPHX), TALPHA(NUABCD,NTALPHX), 273 & INDADR(NUABCD),SCLFAC(NUABCD,NTALPHX),WORK(LWORK) 274#include "twosta.h" 275#include "dftcom.h" 276#include "subdir.h" 277C----------------------------------------------------------------------- 278C 279C ***************** 280C ***** RJ000 ***** 281C ***************** 282C 283C 284C Special case: One-center Integrals 285C ================================== 286C 287C Note: There should be no testing for small integrals since 288C this may in the case of one-center integrals introduce 289C numerical instabilities for large exponents. 290C 291 NOINT = .FALSE. 292 IF (ERFEXP(0)) THEN 293 IEJ = 2 294 IF (SRINTS) IEJ = 3 295 CALL DZERO(EJ000(1,0),(JMAX+1)*NUABCD) 296 END IF 297 IF (AD2DAR) DRWFAC=DP25*ALPHA2*DARFAC 298 IF (DO2DAR) DRWFAC=DP25*ALPHA2 299 IF (S4CENT) DRWFAC=-DP25/PI 300 IF (ONECEN) THEN 301 IODS = 0 302 CALL DZERO(PQX,NUABCD) 303 CALL DZERO(PQY,NUABCD) 304 CALL DZERO(PQZ,NUABCD) 305 IPQ0X = 1 306 IPQ0Y = 1 307 IPQ0Z = 1 308 DO 100 IOD12 = 1, NUC12 309 EXPP = EXP12(IOD12) 310 FAC12I = FAC12(IOD12) 311 IF (DO2DAR .OR. S4CENT) THEN 312C 313C Compute two-electron Dirac delta function 314C (Wim Klopper, University of Utrecht, March 4, 1999) 315C Merged into Dalton1.1 by A. Halkier 29/11/99. 316C 317 DO IOD34 = 1, NUC34 318 IODS = IODS + 1 319 EXPQ = EXP34(IOD34) 320 EXPPQ = EXPP + EXPQ 321 EXPPQI = D1/EXPPQ 322 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 323 TALPHA(IODS,1) = - D2*EXPP*EXPQ*EXPPQI 324 SCLFAC(IODS,1) = FACTOR*TALPHA(IODS,1) 325 RJ000(IODS,0) = DRWFAC*SCLFAC(IODS,1) 326 END DO 327 ELSE IF (AD2DAR) THEN 328C 329C Add two-electron Darwin integrals with perturbation 330C parameter DARFAC onto repulsion integrals. 331C (Wim Klopper, University of Utrecht, March 4, 1999) 332C Merged into Dalton1.1 by A. Halkier 29/11/99. 333C 334 DO IOD34 = 1, NUC34 335 IODS = IODS + 1 336 EXPQ = EXP34(IOD34) 337 EXPPQ = EXPP + EXPQ 338 EXPPQI = D1/EXPPQ 339 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 340 TALPHA(IODS,1) = - D2*EXPP*EXPQ*EXPPQI 341 SCLFAC(IODS,1) = FACTOR*TALPHA(IODS,1) 342 RJ000(IODS,0) = FACTOR+DRWFAC*SCLFAC(IODS,1) 343 END DO 344 ELSE IF (HFXMU .NE. D0) THEN 345 DO IOD34 = 1, NUC34 346 IODS = IODS + 1 347 EXPQ = EXP34(IOD34) 348 EXPPQ = EXPP + EXPQ 349 EXPPQI = D1/EXPPQ 350 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 351 ALPHA = EXPP*EXPQ*EXPPQI 352 BETA = HFXMU**2/(ALPHA + HFXMU**2) 353 TALPHA(IODS,1) = - D2*ALPHA*BETA 354 RJ000(IODS,0) = SQRT(BETA)*FACTOR 355 END DO 356 ELSE IF (PANAS .EQ. D0 .AND. CHIVAL .EQ. -1.D0) THEN 357C 358C Regular two-electron integrals 359C 360 DO IOD34 = 1, NUC34 361 IODS = IODS + 1 362 EXPQ = EXP34(IOD34) 363 EXPPQ = EXPP + EXPQ 364 EXPPQI = D1/EXPPQ 365 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 366 TALPHA(IODS,1) = - D2*EXPP*EXPQ*EXPPQI 367 RJ000(IODS,0) = FACTOR 368 END DO 369 ELSE IF (SRINTS) THEN 370C 371C Short-range two-electron integrals (MCDFT) 372C 373 CHISQ = CHIVAL**2 374 DO IOD34 = 1, NUC34 375 IODS = IODS + 1 376c --- Regular 377 EXPQ = EXP34(IOD34) 378 EXPPQ = EXPP + EXPQ 379 EXPPQI = D1/EXPPQ 380 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 381 TALPHA(IODS,1) = - D2*EXPP*EXPQ*EXPPQI 382 RJ000(IODS,0) = FACTOR 383c --- Modified 384 EXPPI = D1/EXPP 385 EXPQI = D1/EXPQ 386 EXPPQ = EXPPI*EXPQI 387 EXPPQI = EXPPI + EXPQI + D1/CHISQ 388 EXPPQII= D1/EXPPQI 389 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQ*EXPPQII) 390 TALPHA(IODS,2) = - D2*EXPPQII 391 SJ000(IODS,0) = FACTOR 392 IF(COMLAM) THEN 393 SJ000(IODS,0) = SJ000(IODS,0) * VLAMBDA 394 ENDIF 395c --- Old ERFEXP function 2*chival/sqrt(pi)*exp(-chisq/3*r**2) 396 IF(ERFEXP(1)) THEN 397 AI = D3/CHISQ 398 AI = D3/CHISQ 399 EXPFAC = EXPPI + EXPQI + AI 400 EXPFACI = D1/EXPFAC 401 TALPHA(IODS,IEJ) = - D2*EXPFACI 402 EJ000(IODS,0) = - FAC12I*FAC34(IOD34)*CHIVAL* 403 & SQRT((AI*EXPFACI)**3/(EXPP*EXPQ)) 404 ENDIF 405C --- New ERFEXP function (10*chival)/(9*sqrt(pi))*exp(-3*chisq/5 *r**2) 406 IF(ERFEXP(2)) THEN 407 AI = D5/(D3*CHISQ) 408 EXPFAC = EXPPI + EXPQI + AI 409 EXPFACI = D1/EXPFAC 410 TALPHA(IODS,IEJ) = - D2*EXPFACI 411 EJ000(IODS,0) = 412 & - FAC12I*FAC34(IOD34)*(D10/D18)*CHIVAL* 413 & SQRT((AI*EXPFACI)**3/(EXPP*EXPQ)) 414 ENDIF 415c --- 416 END DO 417 ELSE IF (PANAS .EQ. D0) THEN 418C 419C Long-range two-electron integrals (MCDFT) 420C 421 CHISQ = CHIVAL**2 422 DO IOD34 = 1, NUC34 423 IODS = IODS + 1 424 EXPQ = EXP34(IOD34) 425 EXPPI = D1/EXPP 426 EXPQI = D1/EXPQ 427 EXPPQ = EXPPI*EXPQI 428 EXPPQI = EXPPI + EXPQI + D1/CHISQ 429 EXPPQII= D1/EXPPQI 430 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQ*EXPPQII) 431 TALPHA(IODS,1) = - D2*EXPPQII 432 RJ000(IODS,0) = FACTOR 433 IF(COMLAM) THEN 434 RJ000(IODS,0) = RJ000(IODS,0) * VLAMBDA 435 ENDIF 436 IF(ERFEXP(1)) THEN 437 AI = D3/CHISQ 438 EXPFAC = EXPPI + EXPQI + AI 439 EXPFACI = D1/EXPFAC 440 TALPHA(IODS,IEJ) = - D2*EXPFACI 441 EJ000(IODS,0) = - FAC12I*FAC34(IOD34)*CHIVAL* 442 & SQRT((AI*EXPFACI)**3/(EXPP*EXPQ)) 443 ENDIF 444 IF(ERFEXP(2)) THEN 445 AI = D5/(D3*CHISQ) 446 EXPFAC = EXPPI + EXPQI + AI 447 EXPFACI = D1/EXPFAC 448 TALPHA(IODS,IEJ) = - D2*EXPFACI 449 EJ000(IODS,0) = 450 & - FAC12I*FAC34(IOD34)*(D10/D18)*CHIVAL* 451 & SQRT((AI*EXPFACI)**3/(EXPP*EXPQ)) 452 ENDIF 453 END DO 454 ELSE 455C 456C Two-electron integrals regularized according to I.Panas 457C 458 DO IOD34 = 1, NUC34 459 IODS = IODS + 1 460 EXPQ = EXP34(IOD34) 461 EXPPQ = EXPP + EXPQ 462 EXPMOD = EXPPQ/D2 463 EXPPQP = EXPMOD + SQRT(EXPMOD**2 + EXPP*EXPQ*PANAS) 464 FACTOR = FAC12I*FAC34(IOD34)*SQRT(D1/EXPPQP) 465 TALPHA(IODS,1) = - D2*EXPP*EXPQ/EXPPQP 466 RJ000(IODS,0) = FACTOR 467 END DO 468 END IF 469 100 CONTINUE 470 IF (JMAX .GT. 0) THEN 471 IF (DO2DAR .OR. S4CENT) THEN 472 DO J = 1, JMAX 473 DO I = 1, NUABCD 474 SCLFAC(I,1) = TALPHA(I,1)*SCLFAC(I,1) 475 RJ000(I,J) = DRWFAC*SCLFAC(I,1) 476 END DO 477 END DO 478 ELSE IF (AD2DAR) THEN 479 DO J = 1, JMAX 480 FAC = D1/(2.0D0*J + 1.0D0) 481 DO I = 1, NUABCD 482 RJ000(I,J) = (FAC+DRWFAC*TALPHA(I,1))*SCLFAC(I,1) 483 SCLFAC(I,1) = TALPHA(I,1)*SCLFAC(I,1) 484 END DO 485 END DO 486 ELSE 487 CALL DCOPY(NUABCD,TALPHA(1,1),1,SCLFAC(1,1),1) 488 IF(SRINTS) CALL DCOPY(NUABCD,TALPHA(1,2),1,SCLFAC(1,2),1) 489 IF(ERFEXP(0)) 490 & CALL DCOPY(NUABCD,TALPHA(1,IEJ),1,SCLFAC(1,IEJ),1) 491 IF (SRINTS) THEN 492C ... Calculate SJ000 for short-range integrals 493 DO J = 1, JMAX 494 FAC = D1/(2.0D0*J + 1.0D0) 495 DO I = 1, NUABCD 496 SJ000(I,J) = FAC*SCLFAC(I,2)*SJ000(I,0) 497 SCLFAC(I,2) = TALPHA(I,2)*SCLFAC(I,2) 498 END DO 499 END DO 500 END IF 501 IF (ERFEXP(0)) THEN 502C ... Calculate EJ000 for long-range integrals 503 DO J = 1, JMAX 504 DO I = 1, NUABCD 505 EJ000(I,J) = EJ000(I,0)*SCLFAC(I,IEJ) 506 SCLFAC(I,IEJ) = TALPHA(I,IEJ)*SCLFAC(I,IEJ) 507 END DO 508 END DO 509 END IF 510C ... Calculate regular two-electron integrals 511 DO J = 1, JMAX 512 FAC = D1/(2.0D0*J + 1.0D0) 513 DO I = 1, NUABCD 514 RJ000(I,J) = FAC*SCLFAC(I,1)*RJ000(I,0) 515 SCLFAC(I,1) = TALPHA(I,1)*SCLFAC(I,1) 516 END DO 517 END DO 518 END IF 519 END IF ! JMAX .GT. 0 520C 521C General case: Multicenter Integrals 522C =================================== 523C 524 ELSE 525 IF (.NOT.DPATH1) THEN 526 SGN12X = - SIGNT(1) 527 SGN12Y = - SIGNT(2) 528 SGN12Z = - SIGNT(3) 529 SGN34X = - D1 530 SGN34Y = - D1 531 SGN34Z = - D1 532 ELSE 533 SGN12X = D1 534 SGN12Y = D1 535 SGN12Z = D1 536 SGN34X = SIGNT(1) 537 SGN34Y = SIGNT(2) 538 SGN34Z = SIGNT(3) 539 END IF 540C 541 IODS = 0 542 NODS = 0 543 DO 300 IOD12 = 1, NUC12 544 EXPP = EXP12(IOD12) 545 FAC12I = FAC12(IOD12) 546 PX = SGN12X*COOR12(IOD12,1) 547 PY = SGN12Y*COOR12(IOD12,2) 548 PZ = SGN12Z*COOR12(IOD12,3) 549 IF (HFXMU .NE. D0) THEN 550 DO IOD34 = 1, NUC34 551 IODS = IODS + 1 552 EXPQ = EXP34(IOD34) 553 EXPPQI = D1/(EXPP + EXPQ) 554 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 555 ALPHA = EXPP*EXPQ*EXPPQI 556 BETA = HFXMU**2/(ALPHA + HFXMU**2) 557 SCLFAC(IODS,1) = SQRT(BETA)*FACTOR 558 NODS = NODS + 1 559 TALPHA(IODS,1) = - D2*ALPHA*BETA 560 PQXI = PX - SGN34X*COOR34(IOD34,1) 561 PQYI = PY - SGN34Y*COOR34(IOD34,2) 562 PQZI = PZ - SGN34Z*COOR34(IOD34,3) 563 PQX(IODS) = PQXI 564 PQY(IODS) = PQYI 565 PQZ(IODS) = PQZI 566 WVALU(NODS,1)=ALPHA*BETA* 567 & (PQXI*PQXI+PQYI*PQYI+PQZI*PQZI) 568 INDADR(NODS) = IODS 569 END DO 570 ELSE IF (PANAS .EQ. D0 .AND. CHIVAL .EQ. -1.D0) THEN 571C 572C Regular two-electron integrals 573C 574 DO IOD34 = 1, NUC34 575 IODS = IODS + 1 576 EXPQ = EXP34(IOD34) 577 EXPPQI = D1/(EXPP + EXPQ) 578 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 579 SCLFAC(IODS,1) = FACTOR 580 NODS = NODS + 1 581 ALPHA = EXPP*EXPQ*EXPPQI 582 TALPHA(IODS,1) = - D2*ALPHA 583 PQXI = PX - SGN34X*COOR34(IOD34,1) 584 PQYI = PY - SGN34Y*COOR34(IOD34,2) 585 PQZI = PZ - SGN34Z*COOR34(IOD34,3) 586 PQX(IODS) = PQXI 587 PQY(IODS) = PQYI 588 PQZ(IODS) = PQZI 589 WVALU(NODS,1) = ALPHA*(PQXI*PQXI+PQYI*PQYI+PQZI*PQZI) 590 INDADR(NODS) = IODS 591 END DO 592 ELSE IF (SRINTS) THEN 593C 594C Short-range two-electron integrals (MCDFT) 595C 596 CHISQ = CHIVAL**2 597 DO IOD34 = 1, NUC34 598 IODS = IODS + 1 599c --- Regular 600 EXPQ = EXP34(IOD34) 601 EXPPQI = D1/(EXPP + EXPQ) 602 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQI) 603 SCLFAC(IODS,1) = FACTOR 604c --- Modified 605 EXPPI = D1/EXPP 606 EXPQI = D1/EXPQ 607 EXPPQ = EXPPI*EXPQI 608 EXPPQI2= EXPPI + EXPQI + D1/CHISQ 609 EXPPQII= D1/EXPPQI2 610 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQ*EXPPQII) 611 SCLFAC(IODS,2) = FACTOR 612 IF(COMLAM) THEN 613 SCLFAC(IODS,2) = SCLFAC(IODS,2) * VLAMBDA 614 ENDIF 615 NODS = NODS + 1 616 ALPHA = EXPP*EXPQ*EXPPQI 617 TALPHA(IODS,1) = - D2*ALPHA 618 TALPHA(IODS,2) = - D2*EXPPQII 619 PQXI = PX - SGN34X*COOR34(IOD34,1) 620 PQYI = PY - SGN34Y*COOR34(IOD34,2) 621 PQZI = PZ - SGN34Z*COOR34(IOD34,3) 622 PQX(IODS) = PQXI 623 PQY(IODS) = PQYI 624 PQZ(IODS) = PQZI 625 PQXYZI = PQXI*PQXI+PQYI*PQYI+PQZI*PQZI 626 IF(ERFEXP(1)) THEN 627 AI = D3/CHISQ 628 EXPFAC = EXPPI + EXPQI + AI 629 EXPFACI = D1/EXPFAC 630 FACTOR = FAC12I*FAC34(IOD34)* 631 & CHIVAL*SQRT(D1/(EXPP*EXPQ))*(SQRT(AI*EXPFACI))**3 632 EJ000(IODS,0) = - FACTOR * DEXP(- EXPFACI*PQXYZI) 633 TALPHA(IODS,IEJ) = - D2*EXPFACI 634 ENDIF 635 IF(ERFEXP(2)) THEN 636 AI = D5/(CHISQ*D3) 637 EXPFAC = EXPPI + EXPQI + AI 638 EXPFACI = D1/EXPFAC 639 FACTOR = FAC12I*FAC34(IOD34)*(D10/D18)* 640 & CHIVAL*SQRT(D1/(EXPP*EXPQ))*(SQRT(AI*EXPFACI))**3 641 EJ000(IODS,0) = - FACTOR * DEXP(- EXPFACI*PQXYZI) 642 TALPHA(IODS,IEJ) = - D2*EXPFACI 643 ENDIF 644 WVALU(NODS,1) = ALPHA*PQXYZI 645 WVALU(NODS,2) = EXPPQII*PQXYZI 646 INDADR(NODS) = IODS 647 END DO 648 ELSE IF (PANAS .EQ. D0) THEN 649C 650C Long-range two-electron integrals (MCSCF-DFT) 651C 652 CHISQ = CHIVAL**2 653 DO IOD34 = 1, NUC34 654 IODS = IODS + 1 655 EXPQ = EXP34(IOD34) 656 EXPPI = D1/EXPP 657 EXPQI = D1/EXPQ 658 EXPPQ = EXPPI*EXPQI 659 EXPPQI = EXPPI + EXPQI + D1/CHISQ 660 EXPPQII= D1/EXPPQI 661 FACTOR = FAC12I*FAC34(IOD34)*SQRT(EXPPQ*EXPPQII) 662 SCLFAC(IODS,1) = FACTOR 663 IF(COMLAM) THEN 664 SCLFAC(IODS,1) = SCLFAC(IODS,1) * VLAMBDA 665 ENDIF 666 NODS = NODS + 1 667 TALPHA(IODS,1) = - D2*EXPPQII 668 PQXI = PX - SGN34X*COOR34(IOD34,1) 669 PQYI = PY - SGN34Y*COOR34(IOD34,2) 670 PQZI = PZ - SGN34Z*COOR34(IOD34,3) 671 PQX(IODS) = PQXI 672 PQY(IODS) = PQYI 673 PQZ(IODS) = PQZI 674 PQXYZI = PQXI*PQXI+PQYI*PQYI+PQZI*PQZI 675 WVALU(NODS,1) = EXPPQII*PQXYZI 676 IF (ERFEXP(1)) THEN 677 AI = D3/CHISQ 678 EXPFAC = EXPPI + EXPQI + AI 679 EXPFACI = D1/EXPFAC 680 TALPHA(IODS,IEJ) = - D2*EXPFACI 681 FACTOR = FAC12I*FAC34(IOD34)*CHIVAL* 682 & SQRT(D1/(EXPP*EXPQ))*(SQRT(AI*EXPFACI))**3 683 EJ000(IODS,0) = - FACTOR* DEXP(- EXPFACI*PQXYZI) 684 ENDIF 685 IF(ERFEXP(2)) THEN 686 AI = D5/(CHISQ*D3) 687 EXPFAC = EXPPI + EXPQI + AI 688 EXPFACI = D1/EXPFAC 689 FACTOR = FAC12I*FAC34(IOD34)*(D10/D18)* 690 & CHIVAL*SQRT(D1/(EXPP*EXPQ))*(SQRT(AI*EXPFACI))**3 691 EJ000(IODS,0) = - FACTOR * DEXP(- EXPFACI*PQXYZI) 692 TALPHA(IODS,IEJ) = - D2*EXPFACI 693 ENDIF 694 INDADR(NODS) = IODS 695 END DO 696 ELSE 697C 698C Two-electron integrals regularized according to I.Panas 699C 700 DO IOD34 = 1, NUC34 701 IODS = IODS + 1 702 EXPQ = EXP34(IOD34) 703 EXPPQ = EXPP + EXPQ 704 EXPMOD = EXPPQ/D2 705 EXPPQP = EXPMOD + SQRT(EXPMOD**2 + EXPP*EXPQ*PANAS) 706 FACTOR = FAC12I*FAC34(IOD34)*SQRT(D1/EXPPQP) 707 SCLFAC(IODS,1) = FACTOR 708 NODS = NODS + 1 709 ALPHA = EXPP*EXPQ/EXPPQP 710 TALPHA(IODS,1) = - D2*ALPHA 711 PQXI = PX - SGN34X*COOR34(IOD34,1) 712 PQYI = PY - SGN34Y*COOR34(IOD34,2) 713 PQZI = PZ - SGN34Z*COOR34(IOD34,3) 714 PQX(IODS) = PQXI 715 PQY(IODS) = PQYI 716 PQZ(IODS) = PQZI 717 WVALU(NODS,1) = ALPHA*(PQXI*PQXI+PQYI*PQYI+PQZI*PQZI) 718 INDADR(NODS) = IODS 719 END DO 720 END IF 721 300 CONTINUE 722 IF (.NOT.NOINT) THEN 723 IPQ0X = 1 724 IPQ0Y = 1 725 IPQ0Z = 1 726 IF (DASUM(NUABCD,PQX,1) .GT. THRESH) IPQ0X = 0 727 IF (DASUM(NUABCD,PQY,1) .GT. THRESH) IPQ0Y = 0 728 IF (DASUM(NUABCD,PQZ,1) .GT. THRESH) IPQ0Z = 0 729C 730 CALL DZERO(RJ000(1,0),(JMAX+1)*NUABCD) 731 IF (SRINTS) CALL DZERO(SJ000(1,0),(JMAX+1)*NUABCD) 732C 733C Calculate gamma function 734C ======================== 735C 736C Allocations 737C 738 KINDAD = 1 739 KWVALS = KINDAD + (3*NODS + 1)/IRAT 740 KFJW = KWVALS + 3*NODS 741 KREXPW = KFJW + NODS*(JMAX + 1) 742 KLAST = KREXPW + NODS 743 IF (KLAST .GT. LWORK) CALL STOPIT('R0001',' ',KLAST,LWORK) 744 IF (.NOT. (DO2DAR .OR. S4CENT)) THEN 745 CALL GETGAM(NODS,INDADR,WVALU(1,1),RJ000,JMAX,NUABCD, 746 & WORK(KFJW),WORK(KINDAD),WORK(KWVALS), 747 & WORK(KREXPW),IPRINT) 748 IF (SRINTS) 749 & CALL GETGAM(NODS,INDADR,WVALU(1,2),SJ000,JMAX,NUABCD, 750 & WORK(KFJW),WORK(KINDAD),WORK(KWVALS), 751 & WORK(KREXPW),IPRINT) 752 IF (ERFEXP(0)) 753 & CALL DCOPY(NUABCD,TALPHA(1,IEJ),1,SCLFAC(1,IEJ),1) 754 END IF 755C 756C Scale gamma function 757C 758 IF (DO2DAR .OR. AD2DAR .OR. S4CENT) THEN 759 DO I = 1, NUABCD 760 WVALU(I,1) = DRWFAC*TALPHA(I,1)*EXP(-WVALU(I,1)) 761 END DO 762 IF (DO2DAR .OR. S4CENT) THEN 763 DO J = 0, JMAX 764 DO I = 1, NUABCD 765 RJ000(I,J) = WVALU(I,1) 766 END DO 767 END DO 768 ELSE IF (AD2DAR) THEN 769 DO J = 0, JMAX 770 DO I = 1, NUABCD 771 RJ000(I,J) = RJ000(I,J) + WVALU(I,1) 772 END DO 773 END DO 774 END IF 775 ENDIF 776 IF (SRINTS) THEN 777C ... Calculate SJ000 for short-range integrals 778 DO J = 0, JMAX 779 DO I = 1, NUABCD 780 SJ000(I,J) = SCLFAC(I,2)*SJ000(I,J) 781 SCLFAC(I,2) = TALPHA(I,2)*SCLFAC(I,2) 782 END DO 783 END DO 784 END IF 785 IF (ERFEXP(0)) THEN 786C ... Calculate EJ000 for long-range integrals 787 DO J = 1, JMAX 788 DO I = 1, NUABCD 789 EJ000(I,J) = EJ000(I,0)*SCLFAC(I,IEJ) 790 SCLFAC(I,IEJ) = TALPHA(I,IEJ)*SCLFAC(I,IEJ) 791 END DO 792 END DO 793 END IF 794C ... Calculate regular two-electron integrals 795 DO J = 0, JMAX 796 DO I = 1, NUABCD 797 RJ000(I,J) = SCLFAC(I,1)*RJ000(I,J) 798 SCLFAC(I,1) = TALPHA(I,1)*SCLFAC(I,1) 799 END DO 800 END DO 801 END IF 802C 803 END IF 804C 805C For modified operators, R000(I,J) can now be assigned correctly 806C 807 IF(SRINTS) THEN 808 IF (ERFEXP(0)) 809 & CALL DAXPY(NUABCD*(JMAX+1), D1,EJ000(1,0),1,SJ000(1,0),1) 810 CALL DAXPY(NUABCD*(JMAX+1),-D1,SJ000(1,0),1,RJ000(1,0),1) 811 ELSE 812 IF (ERFEXP(0)) 813 & CALL DAXPY(NUABCD*(JMAX+1),D1,EJ000(1,0),1,RJ000(1,0),1) 814 ENDIF 815C 816C ************************* 817C ***** PRINT SECTION ***** 818C ************************* 819C 820 IF (IPRINT .GT. 10) THEN 821 WRITE (LUPRI, 1000) 822 WRITE (LUPRI, 1010) JMAX 823 WRITE (LUPRI, 1020) NUC12, NUC34 824 WRITE (LUPRI, 1030) THRESH 825 WRITE (LUPRI, 1040) ONECEN 826 WRITE (LUPRI, 1045) NOINT 827 WRITE (LUPRI, 1050) NUABCD 828 IF (IPRINT .GT. 20) THEN 829 WRITE (LUPRI, 1100) 830 ISTART = 0 831 DO J = 0, JMAX 832 WRITE (LUPRI, 1110) J 833 IADR = 0 834 DO I = 1, NUC12 835 WRITE (LUPRI, 1120) I 836 WRITE (LUPRI, 1130) (RJ000(IADR + K,J), K = 1, NUC34) 837 IADR = IADR + NUC34 838 END DO 839 END DO 840 END IF 841 END IF 842 1000 FORMAT (//,' ********** SUBROUTINE R000 **********') 843 1010 FORMAT (//,' JMAX ',I7) 844 1020 FORMAT ( ' NUC ',2I7) 845 1030 FORMAT ( ' THRESH ',1P,D12.4) 846 1040 FORMAT ( ' ONECEN ',L7) 847 1045 FORMAT ( ' NOINT ',L7) 848 1050 FORMAT ( ' NUABCD ',I7) 849 1100 FORMAT(//,' ***** RJ000 - INTEGRALS *****') 850 1110 FORMAT( /,' J ',I7) 851 1120 FORMAT( /,' NUC12 ',I7,/) 852 1130 FORMAT(1P,6D12.4) 853 RETURN 854 END 855C /* Deck getgam */ 856 SUBROUTINE GETGAM(NODS,INDADR,WVALU,RJ000,JMAX,NUABCD,FJWS,INDADS, 857 & WVALS,REXPW,IPRINT) 858#include "implicit.h" 859#include "priunit.h" 860#include "maxaqn.h" 861#include "pi.h" 862 PARAMETER (HALF = 0.5D0) 863 PARAMETER (D1 = 1.D0, D10 = 10.D0) 864 PARAMETER (D2 = 2.D0, D4 = 4.D0, D12 = 12.D0, TENTH = 0.1D0) 865 PARAMETER (COEF2 = HALF, COEF3 = - D1/6.D0, COEF4 = D1/24.D0, 866 & COEF5 = - D1/120.D0, COEF6 = D1/720.D0) 867 PARAMETER (GFAC0 = D2*0.4999489092 D0, 868 & GFAC1 = -D2*0.2473631686 D0, 869 & GFAC2 = D2*0.321180909 D0, 870 & GFAC3 = -D2*0.3811559346 D0) 871 PARAMETER (SQRPIH = SQRTPI/D2) 872 PARAMETER (PID4 = PI/D4, PID4I = D4/PI) 873#include "gamcom.h" 874 DIMENSION FJWS(NODS,0:JMAX), INDADR(*), 875 & WVALU(*), RJ000(NUABCD,0:JMAX), INDADS(NODS,3), 876 & WVALS(NODS,3), REXPW(NODS) 877C----------------------------------------------------------------------- 878C 879 NODS1 = 0 880 NODS2 = 0 881 NODS3 = 0 882 D2JP36 = 2.0D0*JMAX + 36.0D0 883 DO 100 I = 1, NODS 884 WVAL = WVALU(I) 885 IF (WVAL .LT. D12) THEN 886 NODS1 = NODS1 + 1 887 INDADS(NODS1,1) = INDADR(I) 888 WVALS(NODS1,1) = WVAL 889 ELSE IF (WVAL .LE. D2JP36) THEN 890 NODS2 = NODS2 + 1 891 INDADS(NODS2,2) = INDADR(I) 892 WVALS(NODS2,2) = WVAL 893 ELSE 894 NODS3 = NODS3 + 1 895 INDADS(NODS3,3) = INDADR(I) 896 WVALS(NODS3,3) = WVAL 897 END IF 898 100 CONTINUE 899C 900C WVAL < 12 901C 902 IF (NODS1 .GT. 0) THEN 903 ISTRT0 = 1 + 121*JMAX 904 DO 200 I = 1, NODS1 905 WVAL = WVALS(I,1) 906 IPNT = NINT(D10*WVAL) 907 WDIF = WVAL - TENTH*IPNT 908 ISTART = ISTRT0 + IPNT 909 FJWS(I,JMAX) = (((((COEF6*TABFJW(ISTART + 726)*WDIF 910 & + COEF5*TABFJW(ISTART + 605))*WDIF 911 & + COEF4*TABFJW(ISTART + 484))*WDIF 912 & + COEF3*TABFJW(ISTART + 363))*WDIF 913 & + COEF2*TABFJW(ISTART + 242))*WDIF 914 & - TABFJW(ISTART + 121))*WDIF 915 & + TABFJW(ISTART) 916 200 CONTINUE 917 IF (JMAX .GT. 0) THEN 918 DO 300 I = 1, NODS1 919 REXPW(I) = HALF*EXP(-WVALS(I,1)) 920 300 CONTINUE 921 DO 310 J = JMAX - 1, 0, -1 922 FCT = D2/(2.0D0*J + 1.0D0) 923 DO 320 I = 1, NODS1 924 FJWS(I,J) = FCT*(WVALS(I,1)*FJWS(I,J+1) + REXPW(I)) 925 320 CONTINUE 926 310 CONTINUE 927 DO 330 J = 1, JMAX 928 DO 340 I = 1, NODS1 929 RJ000(INDADS(I,1),J) = FJWS(I,J) 930 340 CONTINUE 931 330 CONTINUE 932 END IF 933 DO 350 I = 1, NODS1 934 RJ000(INDADS(I,1),0) = FJWS(I,0) 935 350 CONTINUE 936 END IF 937C 938C Near asymptotic region 939C 940 IF (NODS2 .GT. 0) THEN 941 DO 400 I = 1, NODS2 942 WVAL = WVALS(I,2) 943 REXPW(I) = HALF*EXP(-WVAL) 944 WVALS(I,2) = D1/WVAL 945 400 CONTINUE 946 DO 410 I = 1, NODS2 947 RWVAL = WVALS(I,2) 948 GVAL = GFAC0 + RWVAL*(GFAC1 + RWVAL*(GFAC2 + RWVAL*GFAC3)) 949 FJWS(I,0) = SQRPIH*SQRT(RWVAL) - REXPW(I)*GVAL*RWVAL 950 410 CONTINUE 951 DO 420 I = 1, NODS2 952 RJ000(INDADS(I,2),0) = FJWS(I,0) 953 420 CONTINUE 954 DO 430 J = 1, JMAX 955 FCT = J - HALF 956 DO 440 I = 1, NODS2 957 FJWS(I,J) = (FCT*FJWS(I,J-1) - REXPW(I))*WVALS(I,2) 958 440 CONTINUE 959 DO 450 I = 1, NODS2 960 RJ000(INDADS(I,2),J) = FJWS(I,J) 961 450 CONTINUE 962 430 CONTINUE 963 END IF 964C 965C Asymptotic region 966C 967 IF (NODS3 .GT. 0) THEN 968 DO 500 I = 1, NODS3 969 WVALS(I,3) = PID4/WVALS(I,3) 970 500 CONTINUE 971 DO 510 I = 1, NODS3 972 FJWS(I,0) = SQRT(WVALS(I,3)) 973 510 CONTINUE 974 DO 520 I = 1, NODS3 975 RJ000(INDADS(I,3),0) = FJWS(I,0) 976 520 CONTINUE 977 DO 530 J = 1, JMAX 978 FACTOR = PID4I*(J - HALF) 979 DO 540 I = 1, NODS3 980 FJWS(I,J) = FACTOR*FJWS(I,J-1)*WVALS(I,3) 981 540 CONTINUE 982 DO 550 I = 1, NODS3 983 RJ000(INDADS(I,3),J) = FJWS(I,J) 984 550 CONTINUE 985 530 CONTINUE 986 END IF 987C 988 RETURN 989 END 990C /* Deck heri */ 991 SUBROUTINE HERI(HERINT,WORK,LWORK,RJ000,PQX,PQY,PQZ,INDHER,JMAX, 992 & MAXDER,NUABCD,NTUV,IPQ0X,IPQ0Y,IPQ0Z,IODDHR, 993 & IPRINT) 994C 995C tuh fall 1984 996C 997C Modified Jul 28 88 to avoid multiplying zero with 998C undetermined numbers - tuh 999C 1000C Modified for triangular addressing March 92 - tuh 1001C 1002#include "implicit.h" 1003#include "priunit.h" 1004#include "maxaqn.h" 1005 INTEGER T, U, V, TUV 1006 DIMENSION WORK(LWORK), RJ000(NUABCD,0:JMAX), IODDHR(*), 1007 & PQX(NUABCD), PQY(NUABCD), PQZ(NUABCD), 1008 & INDHER(0:JTOP,0:JTOP,0:JTOP), HERINT(NUABCD,NTUV) 1009#include "twosta.h" 1010#include "hertop.h" 1011C 1012C The final R(T,U,V) integrals are arranged as follows: 1013C 1014C R(000) 1015C R(100) R(010) R(001) 1016C R(200) R(110) R(101) R(020) R(011) R(002) 1017C R(300) R(210) R(201) R(120) R(111) R(102) R(030) R(021) R(012) 1018C R(003) 1019C Special case JMAX = 0 1020C ===================== 1021C 1022 IF (JMAX .EQ. 0) THEN 1023 CALL DCOPY(NUABCD,RJ000,1,HERINT,1) 1024 ELSE 1025C 1026C Allocate work space 1027C =================== 1028C 1029 KHRWRK = 1 1030 KLAST = KHRWRK + NTUV*NUABCD 1031 IF (KLAST .GT. LWORK) CALL STOPIT('HERI',' ',KLAST,LWORK) 1032 MWHRSQ = MAX(MWHRSQ,NTUV*NUABCD) 1033 LWTOT = LWTOT + KLAST 1034 MWTOT = MAX(MWTOT,LWTOT) 1035C 1036C Recursion loop for Hermite integrals 1037C ==================================== 1038C 1039 IPQX = IPQ0X + 1 1040 IPQY = IPQ0Y + 1 1041 IPQZ = IPQ0Z + 1 1042 CALL DZERO(HERINT,NUABCD*NTUV) 1043 DO 200 JVAL = 1, JMAX 1044 IF (MOD(JMAX-JVAL,2).EQ.0) THEN 1045 CALL HRECUR(HERINT,WORK(KHRWRK),JVAL,RJ000,PQX,PQY,PQZ, 1046 & INDHER,JMAX,MAXDER,NUABCD,NTUV,IPQX,IPQY, 1047 & IPQZ) 1048 ELSE 1049 CALL HRECUR(WORK(KHRWRK),HERINT,JVAL,RJ000,PQX,PQY,PQZ, 1050 & INDHER,JMAX,MAXDER,NUABCD,NTUV,IPQX,IPQY, 1051 & IPQZ) 1052 END IF 1053 200 CONTINUE 1054 LWTOT = LWTOT - KLAST 1055 END IF 1056C 1057C Print section 1058C ============= 1059C 1060 IF (IPRINT .GE. 10) THEN 1061 CALL TITLER('Output from HERI','*',103) 1062 WRITE (LUPRI,'(2X,A,I5)') 'JMAX ', JMAX 1063 WRITE (LUPRI,'(2X,A,I5)') 'NUABCD', NUABCD 1064 WRITE (LUPRI,'(2X,A,I5)') 'NTUV ', NTUV 1065 IF (IPRINT .GE. 20) THEN 1066 CALL HEADER('Hermite integrals R(t,u,v)',1) 1067 DO 300 J = 0, JMAX 1068 DO 320 T = J, 0, -1 1069 DO 330 U = J - T, 0, -1 1070 V = J - T - U 1071 TUV = INDHER(T,U,V) 1072 IF (IODDHR(TUV) .EQ. 0) THEN 1073 WRITE (LUPRI,'(2X,3(A,I1),A,2X,5F12.8/, 1074 & (12X,5F12.8))') 1075 & 'R(',T,',',U,',',V,')', (HERINT(I,TUV),I=1,NUABCD) 1076 WRITE (LUPRI,'()') 1077 END IF 1078 330 CONTINUE 1079 320 CONTINUE 1080 300 CONTINUE 1081 END IF 1082 END IF 1083 RETURN 1084 END 1085C /* Deck hrecur */ 1086 SUBROUTINE HRECUR(CUR,OLD,JVAL,RJ000,PQX,PQY,PQZ,INDHER,JMAX, 1087 & MAXDER,NUABCD,NTUV,IPQX,IPQY,IPQZ) 1088C 1089#include "implicit.h" 1090 INTEGER T, U, V, TUV 1091 LOGICAL PQXGT0, PQYGT0, PQZGT0 1092 DIMENSION CUR(NUABCD,NTUV), OLD(NUABCD,NTUV), 1093 & INDHER(0:JTOP,0:JTOP,0:JTOP), 1094 & PQX(NUABCD), PQY(NUABCD), PQZ(NUABCD), 1095 & RJ000(NUABCD,0:JMAX) 1096#include "doxyz.h" 1097#include "hertop.h" 1098C 1099 1100C 1101 PQXGT0 = IPQX .EQ. 1 1102 PQYGT0 = IPQY .EQ. 1 1103 PQZGT0 = IPQZ .EQ. 1 1104C 1105C JVAL = 1 1106C ======== 1107C 1108 IF (JVAL .EQ. 1) THEN 1109 CALL DCOPY(NUABCD,RJ000(1,JMAX-1),1,CUR(1,1),1) 1110 IF (PQXGT0) THEN 1111 DO 110 I = 1, NUABCD 1112 CUR(I,2) = PQX(I)*RJ000(I,JMAX) 1113 110 CONTINUE 1114 END IF 1115 IF (PQYGT0) THEN 1116 DO 120 I = 1, NUABCD 1117 CUR(I,3) = PQY(I)*RJ000(I,JMAX) 1118 120 CONTINUE 1119 END IF 1120 IF (PQZGT0) THEN 1121 DO 130 I = 1, NUABCD 1122 CUR(I,4) = PQZ(I)*RJ000(I,JMAX) 1123 130 CONTINUE 1124 END IF 1125C 1126C JVAL > 1 1127C ======== 1128C 1129 ELSE 1130 MAXT = JMAX 1131 MAXU = JMAX 1132 MAXV = JMAX 1133 IF (.NOT.DOX) MAXT = JMAX - MAXDER 1134 IF (.NOT.DOY) MAXU = JMAX - MAXDER 1135 IF (.NOT.DOZ) MAXV = JMAX - MAXDER 1136C 1137C R(0,0,0) 1138C 1139 CALL DCOPY(NUABCD,RJ000(1,JMAX-JVAL),1,CUR,1) 1140C 1141C R(T,0,0) 1142C 1143 IF (PQXGT0) THEN 1144 DO 200 I = 1, NUABCD 1145 CUR(I,2) = PQX(I)*OLD(I,1) 1146 200 CONTINUE 1147 DO 300 T = 2, MIN(MAXT,JVAL) 1148 TMIN1 = T - 1.0D0 1149 TUV = INDHER(T ,0,0) 1150 M1T = INDHER(T-1,0,0) 1151 M2T = INDHER(T-2,0,0) 1152 DO 310 I = 1, NUABCD 1153 CUR(I,TUV) = PQX(I)*OLD(I,M1T) + TMIN1*OLD(I,M2T) 1154 310 CONTINUE 1155 300 CONTINUE 1156 ELSE 1157 DO 400 T = 2, MIN(MAXT,JVAL), 2 1158 TMIN1 = T - 1.0D0 1159 TUV = INDHER(T ,0,0) 1160 M2T = INDHER(T-2,0,0) 1161 DO 410 I = 1, NUABCD 1162 CUR(I,TUV) = TMIN1*OLD(I,M2T) 1163 410 CONTINUE 1164 400 CONTINUE 1165 END IF 1166C 1167C R(T,U,0) 1168C 1169 IF (PQYGT0) THEN 1170 DO 500 T = 0, MIN(MAXT,JVAL - 1), IPQX 1171 TUV = INDHER(T,1,0) 1172 M1U = INDHER(T,0,0) 1173 DO 510 I = 1, NUABCD 1174 CUR(I,TUV) = PQY(I)*OLD(I,M1U) 1175 510 CONTINUE 1176 500 CONTINUE 1177 DO 600 U = 2, MIN(MAXU,JVAL) 1178 UMIN1 = U - 1.0D0 1179 DO 610 T = 0, MIN(MAXT,JVAL - U), IPQX 1180 TUV = INDHER(T,U ,0) 1181 M1U = INDHER(T,U-1,0) 1182 M2U = INDHER(T,U-2,0) 1183 DO 620 I = 1, NUABCD 1184 CUR(I,TUV) = PQY(I)*OLD(I,M1U) + UMIN1*OLD(I,M2U) 1185 620 CONTINUE 1186 610 CONTINUE 1187 600 CONTINUE 1188 ELSE 1189 DO 700 U = 2, MIN(MAXU,JVAL), 2 1190 UMIN1 = U - 1.0D0 1191 DO 710 T = 0, MIN(MAXT,JVAL - U), IPQX 1192 TUV = INDHER(T,U ,0) 1193 M2U = INDHER(T,U-2,0) 1194 DO 720 I = 1, NUABCD 1195 CUR(I,TUV) = UMIN1*OLD(I,M2U) 1196 720 CONTINUE 1197 710 CONTINUE 1198 700 CONTINUE 1199 END IF 1200C 1201C R(T,U,V) 1202C 1203 IF (PQZGT0) THEN 1204 IUMAX = JVAL - 1 1205 DO 800 U = 0, MIN(MAXU,IUMAX), IPQY 1206 DO 810 T = 0, MIN(MAXT,IUMAX - U), IPQX 1207 TUV = INDHER(T,U,1) 1208 M1V = INDHER(T,U,0) 1209 DO 820 I = 1, NUABCD 1210 CUR(I,TUV) = PQZ(I)*OLD(I,M1V) 1211 820 CONTINUE 1212 810 CONTINUE 1213 800 CONTINUE 1214 DO 900 V = 2, MIN(MAXV,JVAL) 1215 VMIN1 = V - 1.0D0 1216 IUMAX = JVAL - V 1217 DO 910 U = 0, MIN(MAXU,IUMAX), IPQY 1218 DO 920 T = 0, MIN(MAXT,IUMAX - U), IPQX 1219 TUV = INDHER(T,U,V ) 1220 M1V = INDHER(T,U,V-1) 1221 M2V = INDHER(T,U,V-2) 1222 DO 930 I = 1, NUABCD 1223 CUR(I,TUV) = PQZ(I)*OLD(I,M1V)+VMIN1*OLD(I,M2V) 1224 930 CONTINUE 1225 920 CONTINUE 1226 910 CONTINUE 1227 900 CONTINUE 1228 ELSE 1229 DO 1000 V = 2, MIN(MAXV,JVAL), 2 1230 VMIN1 = V - 1.0D0 1231 IUMAX = JVAL - V 1232 DO 1010 U = 0, MIN(MAXU,IUMAX), IPQY 1233 DO 1020 T = 0, MIN(MAXT,IUMAX - U), IPQX 1234 TUV = INDHER(T,U,V ) 1235 M2V = INDHER(T,U,V-2) 1236 DO 1030 I = 1, NUABCD 1237 CUR(I,TUV) = VMIN1*OLD(I,M2V) 1238 1030 CONTINUE 1239 1020 CONTINUE 1240 1010 CONTINUE 1241 1000 CONTINUE 1242 END IF 1243 END IF 1244 RETURN 1245 END 1246C /* Deck herswp */ 1247 SUBROUTINE HERSWP(HERINT,NTUV,NUABCD,WORK,LWORK,JMAX,NUCAB, 1248 & NUCCD,IPRINT) 1249C 1250C TUH 87 1251C 1252#include "implicit.h" 1253#include "priunit.h" 1254 DIMENSION HERINT(NUABCD,NTUV), WORK(LWORK) 1255 NHRINT = NTUV*NUABCD 1256 IF (NHRINT .GT. LWORK) CALL STOPIT('HERSWP',' ',LWORK,NHRINT) 1257 ISTR10 = 1 1258 ISTR2 = 1 1259 DO 100 I = 1, NUCAB 1260 ISTR1 = ISTR10 1261 DO 200 J = 1, NUCCD 1262 CALL DCOPY(NTUV,HERINT(ISTR2,1),NUABCD,WORK(ISTR1),NUABCD) 1263 ISTR1 = ISTR1 + NUCAB 1264 ISTR2 = ISTR2 + 1 1265 200 CONTINUE 1266 ISTR10 = ISTR10 + 1 1267 100 CONTINUE 1268 CALL DCOPY(NHRINT,WORK,1,HERINT,1) 1269 IF (IPRINT .GE. 25) THEN 1270 CALL HEADER('OUTPUT FROM HERSWP',-1) 1271 DO 500 ITUV = 1, NTUV 1272 WRITE (LUPRI,'(/,A,I5/)') ' NR ',ITUV 1273 WRITE (LUPRI,'(6F12.8)') (HERINT(I,ITUV),I=1,NUABCD) 1274 500 CONTINUE 1275 END IF 1276 RETURN 1277 END 1278C /* Deck r000x */ 1279#if defined (VAR_OLDGAM) 1280 SUBROUTINE R000X(RJ000,COOR12,COOR34,EXP12,EXP34,FAC12,FAC34,PQX, 1281 & PQY,PQZ,JMAX,NOINT,NUABCD,NUC1,NUC2,NUC12,NUC3, 1282 & NUC4,NUC34,THRESH,ONECEN,IPRINT,IPQ0X,IPQ0Y, 1283 & IPQ0Z) 1284C 1285C TUH 84 1286C 1287#include "implicit.h" 1288#include "priunit.h" 1289#include "maxaqn.h" 1290#include "aovec.h" 1291 PARAMETER (D0 = 0.D0, D1 = 1.D0, D2 = 2.D0) 1292 LOGICAL ONECEN, NOINT 1293 DIMENSION RJ000(NUABCD,0:JMAX), 1294 & PQX(NUABCD), PQY(NUABCD), PQZ(NUABCD), 1295 & COOR12(NUC1*NUC2,3), COOR34(NUC3*NUC4,3), 1296 & EXP12 (NUC1*NUC2), EXP34 (NUC3*NUC4), 1297 & FAC12 (NUC1*NUC2), FAC34 (NUC3*NUC4) 1298#include "gamcom.h" 1299#include "twosta.h" 1300#include "subdir.h" 1301 IF (.NOT.DPATH1) THEN 1302 SIGN = - D1 1303 ELSE 1304 SIGN = D1 1305 END IF 1306 JMAX0 = JMAX 1307 NOINT = .TRUE. 1308 IPQ0X = 1 1309 IPQ0Y = 1 1310 IPQ0Z = 1 1311C 1312C ***************** 1313C ***** RJ000 ***** 1314C ***************** 1315C 1316C Special case: One-center Integrals 1317C ================================== 1318C 1319 IF (ONECEN) THEN 1320 IODS = 0 1321 DO 200 IOD12 = 1, NUC12 1322 EXPP = EXP12(IOD12) 1323 FAC12I = FAC12(IOD12) 1324 DO 210 IOD34 = 1, NUC34 1325 IODS = IODS + 1 1326 EXPQ = EXP34(IOD34) 1327 FAC34I = FAC34(IOD34) 1328 EXPPQI = EXPP + EXPQ 1329 FACTOR = FAC12I*FAC34I/SQRT(EXPPQI) 1330 IF (ABS(FACTOR) .GT. THRESH) THEN 1331 NOINT = .FALSE. 1332 ALPHA = EXPP*EXPQ/EXPPQI 1333 PQX(IODS) = D0 1334 PQY(IODS) = D0 1335 PQZ(IODS) = D0 1336 TALPHA = ALPHA + ALPHA 1337 FJ0INV = D1 1338 DO 220 I = 0, JMAX 1339 RJ000(IODS,I) = FACTOR/FJ0INV 1340 FACTOR = - TALPHA*FACTOR 1341 FJ0INV = FJ0INV + D2 1342 220 CONTINUE 1343 ELSE 1344 DO 225 I = 0 ,JMAX 1345 RJ000(IODS,I) = D0 1346 225 CONTINUE 1347 END IF 1348 210 CONTINUE 1349 200 CONTINUE 1350C 1351C General case: Multicenter Integrals 1352C =================================== 1353C 1354 ELSE 1355 IODS = 0 1356 DO 300 IOD12 = 1, NUC12 1357 EXPP = EXP12(IOD12) 1358 FAC12I = FAC12(IOD12) 1359 PX = COOR12(IOD12,1) 1360 PY = COOR12(IOD12,2) 1361 PZ = COOR12(IOD12,3) 1362 DO 310 IOD34 = 1, NUC34 1363 IODS = IODS + 1 1364 EXPQ = EXP34(IOD34) 1365 FAC34I = FAC34(IOD34) 1366 EXPPQI = EXPP + EXPQ 1367 FACTOR = FAC12I*FAC34I/SQRT(EXPPQI) 1368 IF (ABS(FACTOR) .GT. THRESH) THEN 1369 NOINT = .FALSE. 1370 ALPHA = EXPP*EXPQ/EXPPQI 1371 PQXI = PX - COOR34(IOD34,1) 1372 PQYI = PY - COOR34(IOD34,2) 1373 PQZI = PZ - COOR34(IOD34,3) 1374 IF (ABS(PQXI) .GT. THRESH) IPQ0X = 0 1375 IF (ABS(PQYI) .GT. THRESH) IPQ0Y = 0 1376 IF (ABS(PQZI) .GT. THRESH) IPQ0Z = 0 1377 PQX(IODS) = SIGN*PQXI 1378 PQY(IODS) = SIGN*PQYI 1379 PQZ(IODS) = SIGN*PQZI 1380 WVAL = ALPHA*(PQXI*PQXI + PQYI*PQYI + PQZI*PQZI) 1381 CALL GAMFUN 1382 TALPHA = ALPHA + ALPHA 1383 DO 320 I = 0, JMAX 1384 RJ000(IODS,I) = FACTOR*FJW(I) 1385 FACTOR = - TALPHA*FACTOR 1386 320 CONTINUE 1387 ELSE 1388 DO 325 I = 0 ,JMAX 1389 RJ000(IODS,I) = D0 1390 325 CONTINUE 1391 END IF 1392 310 CONTINUE 1393 300 CONTINUE 1394 END IF 1395C 1396C ************************* 1397C ***** PRINT SECTION ***** 1398C ************************* 1399C 1400 IF (IPRINT .GT. 10) THEN 1401 WRITE (LUPRI, 1000) 1402 WRITE (LUPRI, 1010) JMAX 1403 WRITE (LUPRI, 1020) NUC12, NUC34 1404 WRITE (LUPRI, 1030) THRESH 1405 WRITE (LUPRI, 1040) ONECEN 1406 WRITE (LUPRI, 1045) NOINT 1407 WRITE (LUPRI, 1050) NUABCD 1408 IF (IPRINT .GT. 20) THEN 1409 WRITE (LUPRI, 1100) 1410 ISTART = 0 1411 DO 2000 J = 0, JMAX 1412 WRITE (LUPRI, 1110) J 1413 IADR = 0 1414 DO 2100 I = 1, NUC12 1415 WRITE (LUPRI, 1120) I 1416 WRITE (LUPRI, 1130) (RJ000(IADR + K,J), K = 1, NUC34) 1417 IADR = IADR + NUC34 1418 2100 CONTINUE 1419 2000 CONTINUE 1420 END IF 1421 END IF 1422 1000 FORMAT (//,' ********** SUBROUTINE R000 **********') 1423 1010 FORMAT (//,' JMAX ',I7) 1424 1020 FORMAT ( ' NUC ',2I7) 1425 1030 FORMAT ( ' THRESH ',1P,D12.4) 1426 1040 FORMAT ( ' ONECEN ',L7) 1427 1045 FORMAT ( ' NOINT ',L7) 1428 1050 FORMAT ( ' NUABCD ',I7) 1429 1100 FORMAT(//,' ***** RJ000 - INTEGRALS *****') 1430 1110 FORMAT( /,' J ',I7) 1431 1120 FORMAT( /,' NUC12 ',I7,/) 1432 1130 FORMAT(1P,6D12.4) 1433 RETURN 1434 END 1435#endif 1436