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 rxr */ 20 SUBROUTINE RXR(R2AM,V2AM,SF,TF,Q11,NOCDIM,NSPAIR) 21C 22C This subroutine evaluates the matrix elements 23C 24C Q(ijkl) = Sum_ab (ia|r12|jb) * (ka|(K1 + K2)r12|lb) 25C 26C On input, the array R2AM contains the integrals (ia|r12|jb). 27C NOCDIM is the number of (nonfrozen) occupied orbitals and NSPAIR 28C is the number of pairs of (nonfrozen) occupied orbitals. 29C 30C On output, the arrays SF and TF contain the matrices Q(ijkl) for 31C singlet and triplet coupled pairs, respectively. 32C 33C V2AM (of length 2*NT2AMX) and Q11 (of length 2*NOCDIM**4) are used 34C for scratch space. 35C 36C The one-particle matrix X is read from disk (file name AUXFCK). 37C This is the primary+secondary exchange matrix in the orthonormal basis. 38C It is assumed that this matrix is diagonal in the secondary space. 39C 40C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 41C 42#include "implicit.h" 43#include "priunit.h" 44 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 45 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0) 46 PARAMETER (THRDIA = 1D-9) 47 DIMENSION R2AM(*), SF(*), TF(*), ISB(8) 48 REAL*8 Q11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), V2AM(*), 49 * RR, VV, XBD, XAC, SFAC 50 LOGICAL AVIRT, BVIRT, LDUM 51 INTEGER IDUM 52#include "ccorb.h" 53#include "ccsdsym.h" 54#include "ccsdinp.h" 55#include "r12int.h" 56 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 57 ISB(1) = 0 58 DO 100 ISYM = 2, NSYM 59 NBASF = NORB1(ISYM-1) + NORB2(ISYM-1) 60 NNBASF = NBASF * (NBASF + 1) / 2 61 ISB(ISYM) = ISB(ISYM-1) + NNBASF 62 100 CONTINUE 63 NBASF = NORB1(NSYM) + NORB2(NSYM) 64 NNBASF = ISB(NSYM) + NBASF * (NBASF + 1) / 2 65 LUMULB = -34 66 CALL GPOPEN(LUMULB,'AUXFCK','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 67 READ (LUMULB,'(4E30.20)') (SF(I), I = 1, NNBASF) 68 CALL GPCLOSE (LUMULB,'KEEP') 69 IF (.NOT. (R12NOB .OR. NORXR)) THEN 70 DO 200 ISYM = 1, NSYM 71 DO 210 I = NORB1(ISYM) + 2, NVIR(ISYM) 72 DO 220 J = NORB1(ISYM) + 1 , I - 1 73 IJ = ISB(ISYM) + INDEX(I,J) 74 IF (ABS(SF(IJ)) .GT. THRDIA) THEN 75 WRITE(LUPRI,'(/A/A,3I5,E20.10/)') 76 * '@ WARNING : Exchange matrix not diagonal', 77 * ' Nondiagonal element is :', 78 * ISYM,I,J,SF(IJ) 79 IF (ABS(SF(IJ)) .GT. SQRT(THRDIA)) 80 * CALL QUIT('Exchange matrix not diagonal') 81 END IF 82 220 CONTINUE 83 210 CONTINUE 84 200 CONTINUE 85 END IF 86C 87C Compute (ak|bl) = Sum_c X(ac)*(ck|bl) + Sum_d X(bd)*(ak|dl) 88C The sum over c runs over the auxiliary basis only 89C 90 DO 301 ISYMKA = 1, NSYM 91 ISYMLB = ISYMKA 92 DO 302 ISYMK = 1, NSYM 93 ISYMA = MULD2H(ISYMK,ISYMKA) 94 ISYMC = ISYMA 95 ISYMKC = ISYMKA 96 DO 303 ISYML = 1, NSYM 97 ISYMB = MULD2H(ISYML,ISYMLB) 98 ISYMD = ISYMB 99 ISYMLD = ISYMLB 100 DO 304 K = 1, NRHF(ISYMK) 101 DO 305 L = 1, NRHF(ISYML) 102 DO 306 A = 1, NVIR(ISYMA) 103 IF (R12CBS) THEN 104 NSTRC = 1 105 NENDC = NVIR(ISYMC) 106 ELSE 107 IF (A. LE. NORB1(ISYMA)) THEN 108 NSTRC = NORB1(ISYMC) + 1 109 NENDC = NVIR(ISYMC) 110 ELSE 111 NSTRC = A 112 NENDC = A 113 END IF 114 END IF 115 NAK = IT1AM(ISYMA,ISYMK) 116 * + NVIR(ISYMA)*(K-1) + A 117 DO 307 B = 1, NVIR(ISYMB) 118 IF (R12CBS) THEN 119 NSTRD = 1 120 NENDD = NVIR(ISYMD) 121 ELSE 122 IF (B. LE. NORB1(ISYMB)) THEN 123 NSTRD = NORB1(ISYMD) + 1 124 NENDD = NVIR(ISYMD) 125 ELSE 126 NSTRD = B 127 NENDD = B 128 END IF 129 END IF 130 NBL = IT1AM(ISYMB,ISYML) 131 * + NVIR(ISYMB)*(L-1) + B 132 NAKBL = IT2AM(ISYMKA,ISYMLB) 133 * + INDEX(NAK,NBL) 134C 135 V2AM(NAKBL) = D0 136C 137 DO 308 C = NSTRC, NENDC 138 NCK = IT1AM(ISYMC,ISYMK) 139 * + NVIR(ISYMC)*(K-1) + C 140 NCKBL = IT2AM(ISYMKC,ISYMLB) 141 * + INDEX(NCK,NBL) 142 RR = R2AM(NCKBL) 143 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 144 V2AM(NAKBL) = V2AM(NAKBL) + XAC * RR 145 308 CONTINUE 146 DO 309 D = NSTRD, NENDD 147 NDL = IT1AM(ISYMD,ISYML) 148 * + NVIR(ISYMD)*(L-1) + D 149 NAKDL = IT2AM(ISYMKA,ISYMLD) 150 * + INDEX(NAK,NDL) 151 RR = R2AM(NAKDL) 152 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 153 V2AM(NAKBL) = V2AM(NAKBL) + XBD * RR 154 309 CONTINUE 155C 156 307 CONTINUE 157 306 CONTINUE 158 305 CONTINUE 159 304 CONTINUE 160 303 CONTINUE 161 302 CONTINUE 162 301 CONTINUE 163C 164 DO 999 IFLAG = 1, 2 165C 166 DO 110 I = 1, NOCDIM**4 167 Q11(I,1,1,1) = D0 168 110 CONTINUE 169 DO 401 ISYMIA = 1, NSYM 170 ISYMJB = ISYMIA 171 DO 402 ISYMI = 1, NSYM 172 ISYMA = MULD2H(ISYMI,ISYMIA) 173 DO 403 ISYMJ = 1, NSYM 174 ISYMB = MULD2H(ISYMJ,ISYMJB) 175 DO 404 I = 1, NRHF(ISYMI) 176 KOFFI = IRHF(ISYMI) + I 177 DO 405 J = 1, NRHF(ISYMJ) 178 KOFFJ = IRHF(ISYMJ) + J 179 DO 406 A = 1, NVIR(ISYMA) 180 ISYMJA = MULD2H(ISYMJ,ISYMA) 181 DO 407 B = 1, NVIR(ISYMB) 182C 183 AVIRT = A .GT. NRHFSA(ISYMA) .AND. 184 * A .LE. NORB1(ISYMA) 185 BVIRT = B .GT. NRHFSA(ISYMB) .AND. 186 * B .LE. NORB1(ISYMB) 187 188 IF (R12CBS) THEN 189 SFAC = D1 190 IF (IFLAG .EQ. 1) THEN 191 IF (A .LE. NORB1(ISYMA) .OR. 192 * B .LE. NORB1(ISYMB)) GOTO 407 193 ELSE IF (IFLAG .EQ. 2) THEN 194 IF (A .LE. NRHFSA(ISYMA) .OR. 195 * B .LE. NRHFSA(ISYMB)) GOTO 407 196 ELSE IF (IFLAG .EQ. 3) THEN 197 IF (A .LE. NRHFSA(ISYMA) .OR. 198 * B .LE. NRHFSA(ISYMB) .OR. 199 * (AVIRT .AND. BVIRT)) GOTO 407 200 END IF 201 ELSE 202 IF (IFLAG .EQ. 2) THEN 203 IF (AVIRT .OR. BVIRT) GOTO 407 204 END IF 205 IF ((A .GT. NORB1(ISYMA) .AND. 206 * B .GT. NORB1(ISYMB)) .OR. 207 * (A .LE. NORB1(ISYMA) .AND. 208 * B .LE. NORB1(ISYMB))) THEN 209 SFAC = D1 210 ELSE 211 SFAC = - D1 212 END IF 213 IF (IFLAG .EQ. 3) THEN 214 IF ((A .GT. NORB1(ISYMA) .AND. 215 * B .GT. NORB1(ISYMB)) .OR. 216 * (A .LE. NRHFSA(ISYMA) .AND. 217 * B .LE. NRHFSA(ISYMB))) THEN 218 SFAC = D1 219 ELSE 220 SFAC = - D1 221 END IF 222 IF (AVIRT) THEN 223 IF (.NOT. BVIRT) GOTO 407 224 ELSE IF (BVIRT) THEN 225 IF (.NOT. AVIRT) GOTO 407 226 END IF 227 END IF 228 END IF 229C 230 NAI = IT1AM(ISYMA,ISYMI) 231 * + NVIR(ISYMA)*(I-1) + A 232 NBJ = IT1AM(ISYMB,ISYMJ) 233 * + NVIR(ISYMB)*(J-1) + B 234 NAIBJ = IT2AM(ISYMIA,ISYMJB) 235 * + INDEX(NAI,NBJ) 236 DO 408 ISYMKA = 1, NSYM 237 ISYMLB = ISYMKA 238 ISYMK = MULD2H(ISYMA,ISYMKA) 239 ISYML = MULD2H(ISYMB,ISYMLB) 240 DO 409 K = 1, NRHF(ISYMK) 241 KOFFK = IRHF(ISYMK) + K 242 DO 410 L = 1, NRHF(ISYML) 243 KOFFL = IRHF(ISYML) + L 244 NAK = IT1AM(ISYMA,ISYMK) 245 * + NVIR(ISYMA)*(K-1) + A 246 NBL = IT1AM(ISYMB,ISYML) 247 * + NVIR(ISYMB)*(L-1) + B 248 NAKBL = IT2AM(ISYMKA,ISYMLB) 249 * + INDEX(NAK,NBL) 250 RR = R2AM(NAIBJ) 251 VV = V2AM(NAKBL) 252 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 253 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 254 * SFAC * RR * VV 255 410 CONTINUE 256 409 CONTINUE 257 408 CONTINUE 258 407 CONTINUE 259 406 CONTINUE 260 405 CONTINUE 261 404 CONTINUE 262 403 CONTINUE 263 402 CONTINUE 264 401 CONTINUE 265C 266 IJKL = 0 267 FF = D1 / SQRT(D2) 268 DO 600 K = 1, NOCDIM 269 DO 601 L = 1, K 270 DO 602 I = 1, NOCDIM 271 DO 603 J = 1, I 272 IJKL = IJKL + 1 273 RR = Q11(I,J,K,L) + Q11(I,J,L,K) 274 VV = Q11(I,J,K,L) - Q11(I,J,L,K) 275 SF(IJKL) = RR 276 TF(IJKL) = VV 277 IF (I .EQ. J) SF(IJKL) = FF * SF(IJKL) 278 IF (K .EQ. L) SF(IJKL) = FF * SF(IJKL) 279 SF(IJKL) = SF(IJKL) * DP25 * DP5 280 TF(IJKL) = TF(IJKL) * DP75 * DP5 281 603 CONTINUE 282 602 CONTINUE 283 601 CONTINUE 284 600 CONTINUE 285C 286 CALL ERISFK(SF,NSPAIR,1) 287 CALL ERISFK(TF,NSPAIR,1) 288C 289 IF (IPRINT .GT. 3) THEN 290 GOTO (703,704,705), IFLAG 291 CALL QUIT('INVALID IFLAG IN RXR') 292 703 CALL AROUND('Singlet <RXR> matrix') 293 CALL OUTPUT(SF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 294 CALL AROUND('Triplet <RXR> matrix') 295 CALL OUTPUT(TF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 296 GOTO 706 297 704 CALL AROUND('Singlet <RXR@> matrix') 298 CALL OUTPUT(SF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 299 CALL AROUND('Triplet <RXR@> matrix') 300 CALL OUTPUT(TF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 301 GOTO 706 302 705 CALL AROUND('Singlet <RXR3> matrix') 303 CALL OUTPUT(SF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 304 CALL AROUND('Triplet <RXR3> matrix') 305 CALL OUTPUT(TF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 306 706 WRITE(LUPRI,*) 307 END IF 308C 309 LUMUL = -34 310 GOTO (701,702,707), IFLAG 311 CALL QUIT('INVALID IFLAG IN RXR') 312 701 CALL GPOPEN(LUMUL,'XMATSN','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 313 WRITE(LUMUL,'(4E30.20)') (SF(KL), KL = 1, NSPAIR*NSPAIR) 314 CALL GPCLOSE(LUMUL,'KEEP') 315 CALL GPOPEN(LUMUL,'XMATTN','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 316 WRITE(LUMUL,'(4E30.20)') (TF(KL), KL = 1, NSPAIR*NSPAIR) 317 CALL GPCLOSE(LUMUL,'KEEP') 318 GOTO 999 319 702 CALL GPOPEN(LUMUL,'XMATSP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 320 WRITE(LUMUL,'(4E30.20)') (SF(KL), KL = 1, NSPAIR*NSPAIR) 321 CALL GPCLOSE(LUMUL,'KEEP') 322 CALL GPOPEN (LUMUL,'XMATTP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 323 WRITE(LUMUL,'(4E30.20)') (TF(KL), KL = 1, NSPAIR*NSPAIR) 324 CALL GPCLOSE(LUMUL,'KEEP') 325cWK GOTO 999 326 707 CALL GPOPEN(LUMUL,'XMATSQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 327 WRITE(LUMUL,'(4E30.20)') (SF(KL), KL = 1, NSPAIR*NSPAIR) 328 CALL GPCLOSE(LUMUL,'KEEP') 329 CALL GPOPEN(LUMUL,'XMATTQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 330 WRITE(LUMUL,'(4E30.20)') (TF(KL), KL = 1, NSPAIR*NSPAIR) 331 CALL GPCLOSE(LUMUL,'KEEP') 332 999 CONTINUE 333 RETURN 334 END 335C /* Deck makegx */ 336 SUBROUTINE MAKEGX(V2AM,GX,F2AM,WORK,LWORK,IPBAS) 337C 338C This subroutine computes the X(p'k) matrix (exchange operator) 339C by summing two-electron integrals (ip'|ik) in the orthonormal 340C basis. The orbitals i and k belong to the primary basis, 341C whereas p' belongs to the secondary basis. 342C 343C On input, V2AM contains the (ip'|ik) integrals and F2AM 344C contains the (Ip'|Ik) integrals where the orbitals I 345C belong to the frozen core. The X(p'k) exchange matrix 346C is returned as GX. 347C 348C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 349C 350C Modified for natural virtual orbitals NATVIR (WK/UniKA/22-11-2005). 351C 352#include "implicit.h" 353#include "priunit.h" 354 DIMENSION V2AM(*), GX(*), F2AM(*), WORK(*) 355 integer iv1fro(8,8),it1vm(8,8),it2vm(8,8),iv2fro(8,8) 356#include "ccorb.h" 357#include "ccsdsym.h" 358#include "ccsdinp.h" 359#include "ccfro.h" 360#include "r12int.h" 361 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 362C 363 IF (NATVIR) THEN 364 CALL DZERO(GX,NT1AMX) 365 NNBASF = 0 366 DO ISYM = 1, NSYM 367 NONV = NVIR(ISYM) + NRHFS(ISYM) 368 NNBASF = NNBASF + NONV * (NONV + 1) / 2 369 END DO 370 LUNT = -1 371 IF (NNBASF .GT. LWORK) CALL QUIT('Not enough memory in MAKEGX') 372 CALL GPOPEN(LUNT,'AUXFCKN','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 373 READ(LUNT,'(4E30.20)') (WORK(I), I = 1, NNBASF) 374 CALL GPCLOSE(LUNT,'KEEP') 375 IOFF = 0 376 DO ISYMP = 1, NSYM 377 ISYMK = ISYMP 378 NPK = IT1AM(ISYMP,ISYMK) 379 DO K =1 + NRHFFR(ISYMK), NRHFS(ISYMK) 380 DO P = 1 + NRHFS(ISYMP), NVIR(ISYMP) + NRHFS(ISYMP) 381 NPK = NPK + 1 382 GX(NPK) = WORK(IOFF + P*(P-1)/2 + K) 383 END DO 384 END DO 385 NONV = NVIR(ISYMP) + NRHFS(ISYMP) 386 IOFF = IOFF + NONV * (NONV + 1) / 2 387 END DO 388c LUNT = -1 389c CALL GPOPEN(LUNT,'AUXFCK','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 390c READ(LUNT,'(4E30.20)') (WORK(I), I = 1, NNBASF) 391c CALL GPCLOSE(LUNT,'KEEP') 392c IOFF = 0 393c DO ISYMP = 1, NSYM 394c ISYMK = ISYMP 395c NPK = IT1AM(ISYMP,ISYMK) 396c DO K = NRHFFR(ISYMK)+1, NRHFSA(ISYMK) 397c DO P = 1, NVIR(ISYMP) 398c NPK = NPK + 1 399c GX(NPK) = WORK(IOFF + P*(P-1)/2 + K) 400c END DO 401c END DO 402c IOFF = IOFF + NVIR(ISYMP) * (NVIR(ISYMP) + 1) / 2 403c END DO 404 GOTO 100 405 END IF 406 407 IF (R12PRP) THEN 408 CALL DZERO(GX,NT1VMX) 409 DO ISYMK = 1,NSYM 410 NV1FRO(ISYMK) = 0 411 DO ISYMJ = 1,NSYM 412 ISYMI = MULD2H(ISYMK,ISYMJ) 413 NV1FRO(ISYMK) = NV1FRO(ISYMK) + NRHFFR(ISYMJ)* 414 & (NORB1(ISYMI)-NRHFFR(ISYMI)) 415 ENDDO 416 ENDDO 417 418 DO ISYMK = 1, NSYM 419 ICOUN1 = 0 420 ICOUN3 = 0 421 ICOUN4 = 0 422 DO ISYMJ = 1,NSYM 423 ISYMI = MULD2h(ISYMK,ISYMJ) 424 IT1VM(ISYMI,ISYMJ) = ICOUN1 425 ICOUN1 = ICOUN1 + (NORB1(ISYMJ)-NRHFFR(ISYMJ)) 426 & *NVIR(ISYMI) 427 IV2FRO(ISYMI,ISYMJ) = ICOUN3 428 IV1FRO(ISYMI,ISYMJ) = ICOUN4 429 ICOUN3 = ICOUN3 + NT1FRO(ISYMI)*NV1FRO(ISYMJ) 430 ICOUN4 = ICOUN4 + (NRHFFR(ISYMI)) 431 & *(NORB1(ISYMJ)-NRHFFR(ISYMJ)) 432 ENDDO 433 ENDDO 434 ELSE 435 CALL DZERO(GX,NT1AMX) 436 ENDIF 437C 438C Compute SUM_i ( i p' | i k ) 439C 440celena 441 IF (R12PRP) THEN 442 DO ISYMIP = 1, NSYM 443 ISYMIK = ISYMIP 444 DO ISYMI = 1, NSYM 445 ISYMP = MULD2H(ISYMI,ISYMIP) 446 ISYMK = ISYMP 447 DO I = 1, NRHF(ISYMI) 448 NPK = IT1vM(ISYMP,ISYMK) 449 DO K = NRHFFR(ISYMK)+1, NORB1(ISYMK) 450 DO P = 1, NVIR(ISYMP) 451 IF (.NOT. ONEAUX) THEN 452 NPK = NPK + 1 453 NPI = IT1AM(ISYMP,ISYMI) 454 * + NVIR(ISYMP)*(I-1) + P 455 NKI = IT1AM(ISYMK,ISYMI) 456 * + NVIR(ISYMK)*(I-1) + K 457 NPIKI = IT2AM(ISYMIP,ISYMIK) 458 * + INDEX(NPI,NKI) 459 GX(NPK) = GX(NPK) + V2AM(NPIKI) 460 ENDIF 461 ENDDO 462 ENDDO 463 ENDDO 464 ENDDO 465 ENDDO 466ccelena 467 ELSE 468 DO 101 ISYMIP = 1, NSYM 469 ISYMIK = ISYMIP 470 IKOFF = NH1AM(ISYMIK) * (NH1AM(ISYMIK) + 1) / 2 471 DO 102 ISYMI = 1, NSYM 472 ISYMP = MULD2H(ISYMI,ISYMIP) 473 ISYMK = ISYMP 474 DO 201 I = 1, NRHFA(ISYMI) 475 NPK = IT1AM(ISYMP,ISYMK) 476 DO 202 K = NRHFFR(ISYMK)+1, NRHFS(ISYMK) 477 DO 203 P = 1, NVIR(ISYMP) 478 NPK = NPK + 1 479 IF (ONEAUX) THEN 480 IF (P .LE. NORB1(ISYMP)) THEN 481 NPI = IH1AM(ISYMP,ISYMI) 482 * + NORB1(ISYMP)*(I-1) + P 483 NKI = IH1AM(ISYMK,ISYMI) 484 * + NORB1(ISYMK)*(I-1) + K 485 NPIKI = IH2AM(ISYMIP,ISYMIK) 486 * + INDEX(NPI,NKI) 487 ELSE 488 NPI = IG1AM(ISYMP,ISYMI) 489 * + NORB2(ISYMK)*(I-1) + P - NORB1(ISYMP) 490 NKI = IH1AM(ISYMK,ISYMI) 491 * + NORB1(ISYMK)*(I-1) + K 492 NPIKI = IKOFF + IH2AM(ISYMIP,ISYMIK) 493 * + NH1AM(ISYMIK) * (NPI - 1) + NKI 494 END IF 495 ELSE 496 NPI = IT1AM(ISYMP,ISYMI) 497 * + NVIR(ISYMP)*(I-1) + P 498 NKI = IH1AM(ISYMK,ISYMI) 499 * + NVIR(ISYMK)*(I-1) + K 500 NPIKI = IT2AM(ISYMIP,ISYMIK) 501 * + INDEX(NPI,NKI) 502 END IF 503 GX(NPK) = GX(NPK) + V2AM(NPIKI) 504 203 CONTINUE 505 202 CONTINUE 506 201 CONTINUE 507 102 CONTINUE 508 101 CONTINUE 509 ENDIF 510C 511 IF (FROIMP) THEN 512C 513C Add contribution from frozen orbitals 514C 515 IF (R12PRP) THEN 516 DO ISYMIP = 1, NSYM 517 ISYMIK = ISYMIP 518 DO ISYMI = 1, NSYM 519 ISYMP = MULD2H(ISYMI,ISYMIP) 520 ISYMK = ISYMP 521 DO I = 1, NRHFFR(ISYMI) 522 NPK = IT1VM(ISYMP,ISYMK) 523 DO K = 1, NORB1(ISYMK)-NRHFFR(ISYMK) 524 DO P = 1, NVIR(ISYMP) 525 NPK = NPK + 1 526 NKI = IV1FRO(ISYMI,ISYMK) 527 * + NRHFFR(ISYMI)*(K-1) + I 528 NPI = IT1FRO(ISYMP,ISYMI) 529 * + NVIR(ISYMP)*(I-1) + P 530 NPIKI = IV2FRO(ISYMIP,ISYMIK) 531 * + NT1FRO(ISYMIP)*(NKI-1) 532 * + NPI 533 GX(NPK) = GX(NPK) + F2AM(NPIKI) 534 ENDDO 535 ENDDO 536 ENDDO 537 ENDDO 538 ENDDO 539 ELSE 540 DO 301 ISYMIP = 1, NSYM 541 ISYMIK = ISYMIP 542 DO 302 ISYMI = 1, NSYM 543 ISYMP = MULD2H(ISYMI,ISYMIP) 544 ISYMK = ISYMP 545 DO 401 I = 1, NRHFFR(ISYMI) 546 NPK = IT1AM(ISYMP,ISYMK) 547 DO 402 K = 1, NRHF(ISYMK) 548 DO 403 P = 1, NVIR(ISYMP) 549 NPK = NPK + 1 550 NKI = IF1FRO(ISYMI,ISYMK) 551 * + NRHFFR(ISYMI)*(K-1) + I 552 NPI = IT1FRO(ISYMP,ISYMI) 553 * + NVIR(ISYMP)*(I-1) + P 554 NPIKI = IF2FRO(ISYMIP,ISYMIK) 555 * + NT1FRO(ISYMIP)*(NKI-1) 556 * + NPI 557 GX(NPK) = GX(NPK) + F2AM(NPIKI) 558 403 CONTINUE 559 402 CONTINUE 560 401 CONTINUE 561 302 CONTINUE 562 301 CONTINUE 563 ENDIF 564 END IF 565C 566 100 CONTINUE 567C 568 IF (R12OLD .OR. R12CBS) GOTO 700 569C 570C Zero primary block 571C 572 IF (.NOT. R12HYB .OR. (R12HYB .AND. IPBAS .EQ. 2)) THEN 573 DO 500 ISYM = 1, NSYM 574 IF (R12PRP) THEN 575 NPK = IT1vM(ISYM,ISYM) + 1 576 DO K = 1, NORB1(ISYM)-NRHFFR(ISYM) 577 CALL DZERO(GX(NPK),NORB1(ISYM)) 578 NPK = NPK + NVIR(ISYM) 579 ENDDO 580 ELSE 581 NPK = IT1AM(ISYM,ISYM) + 1 582 DO 501 K = 1, NRHF(ISYM) 583 CALL DZERO(GX(NPK),NORB1(ISYM)) 584 NPK = NPK + NVIR(ISYM) 585 501 CONTINUE 586 ENDIF 587 500 CONTINUE 588 END IF 589C 590 700 CONTINUE 591C 592C Zero secondary block 593C 594 IF (R12HYB .AND. IPBAS .EQ. 1) THEN 595 DO 600 ISYM = 1, NSYM 596 IF (R12PRP) THEN 597 NPK = IT1vM(ISYM,ISYM) + NORB1(ISYM)+ 1 598 DO K = 1, NORB1(ISYM)-NRHFFR(ISYM) 599 CALL DZERO(GX(NPK),NORB2(ISYM)) 600 NPK = NPK + NVIR(ISYM) 601 ENDDO 602 ELSE 603 NPK = IT1AM(ISYM,ISYM) + NORB1(ISYM) + 1 604 DO 601 K = 1, NRHF(ISYM) 605 CALL DZERO(GX(NPK),NORB2(ISYM)) 606 NPK = NPK + NVIR(ISYM) 607 601 CONTINUE 608 ENDIF 609 600 CONTINUE 610 END IF 611C 612C Print section 613C 614 IF (IPRINT .GT. 3) THEN 615 DO 300 ISYM = 1, NSYM 616 IF (NRHF(ISYM)*NVIR(ISYM) .NE. 0) 617 * WRITE(LUPRI,'(/A,I2)') ' X-MATRIX/ISYM =',ISYM 618 IF (R12PRP) THEN 619 NPK = IT1vM(ISYM,ISYM) + 1 620 CALL OUTPUT(GX(NPK),1,NVIR(ISYM),1,NORB1(ISYM), 621 * NVIR(ISYM),NORB1(ISYM),1,LUPRI) 622 ELSE 623 NPK = IT1AM(ISYM,ISYM) + 1 624 CALL OUTPUT(GX(NPK),1,NVIR(ISYM),1,NRHF(ISYM), 625 * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 626 ENDIF 627 300 CONTINUE 628 END IF 629 RETURN 630 END 631C /* Deck r12drv */ 632 SUBROUTINE R12DRV(V2AM,R2AM,U2AM,S2AM,T2AM, 633 * EV,RXRS,RXRT,WRK,LWRK, 634 * QQ2,QQ4,QQ6) 635C 636C Driver routine for the MP2-R12 calculation. 637C 638C V2AM = ( ip | 1/r12 | jq) 639C 640C R2AM = ( ip | r12 | jq) 641C 642C U2AM = ( ip | [T1+T2,r12] | jq) 643C 644C S2AM = ( ip | [r12,K1+K2] | jq) 645C 646C T2AM = ( ip | (K1+K2) r12 | jq) 647C 648C QQ2 = ( ik | (1/r12)*f(12) | jl ) 649C 650C QQ4 = ( ik | f(12)**2 | jl ) 651C 652C QQ6 = ( ik | [[f12,T1],f12] | jl ) 653C 654C p and q denote orbitals of both the primary and 655C secondary basis. Integrals where both p and q belong 656C to the secondary basis are not needed. The numerical 657C values at the corresponding places in the arrays are, 658C meaningless! 659C 660C EV is an array with orbital energies and RXRS and RXRT 661C contain the matrices computed in the subroutine RXR. 662C 663C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 664C 665#include "implicit.h" 666#include "priunit.h" 667#include "ccorb.h" 668#include "ccsdsym.h" 669#include "ccsdinp.h" 670#include "r12int.h" 671#include "maxorb.h" 672#include "infinp.h" 673 DIMENSION R2AM(*), V2AM(*), U2AM(*), EV(*), WRK(*), 674 * RXRS(*), RXRT(*), S2AM(*), T2AM(*), 675 * QQ2(*), QQ4(*), QQ6(*),ISB(8) 676 LOGICAL LGLOCAL,FRSWRT 677C 678 LGLOCAL = BOYORB.OR.PIPORB 679C 680 NOCDIM = NRHFT 681 NICDIM = NRHFTS 682 NACDIM = 0 683 DO ISYM = 1, NSYM 684 NACDIM = NACDIM + NRHFA(ISYM) 685 END DO 686 NAPAIR = NACDIM * (NACDIM + 1) / 2 687 NSPAIR = NOCDIM * (NOCDIM + 1) / 2 688 NTPAIR = NOCDIM * (NOCDIM - 1) / 2 689 NPAIR2 = NSPAIR ** 2 690 NODIM4 = 2 * NOCDIM ** 4 691 NIDIM4 = NICDIM ** 4 692 IF (LGLOCAL) THEN 693 NCT = 0 694 NRT = 0 695 DO ISYM = 1, NSYM 696 NCT = NCT + NRHFFR(ISYM) 697 NRT = NRT + NORB1(ISYM) 698 END DO 699 NBASM = (NRT-NCT)*(NRT-NCT+1)/2 700 ELSE 701 NBASM = NOCDIM*(NOCDIM+1)/2 702 END IF 703C 704 IES = 1 705 IET = IES + NSPAIR 706 IFS = IET + NSPAIR 707 IFT = IFS + NSPAIR 708 IEVS = IFT + NSPAIR 709 ININV = IEVS + NSPAIR 710 IQQ11 = ININV + NSPAIR*10 711C 712 IV11 = IQQ11 + NIDIM4 713 IU11 = IV11 + NODIM4 714 IB11 = IU11 + NODIM4 715 IW11 = IB11 + NODIM4 * NSPAIR 716 IQ11 = IW11 + NODIM4 717 IR11 = IQ11 + NODIM4 718C 719 IVS11 = IR11 + NODIM4 720 IUS11 = IVS11 + NPAIR2 721 IBS11 = IUS11 + NPAIR2 722 IWS11 = IBS11 + NPAIR2 * NSPAIR 723 IQS11 = IWS11 + NPAIR2 724 IRS11 = IQS11 + NPAIR2 725C 726 IVT11 = IRS11 + NPAIR2 727 IUT11 = IVT11 + NPAIR2 728 IBT11 = IUT11 + NPAIR2 729 IWT11 = IBT11 + NPAIR2 * NSPAIR 730 IQT11 = IWT11 + NPAIR2 731 IRT11 = IQT11 + NPAIR2 732C 733 IFLMAT= IRT11 + NPAIR2 734 ININ = IFLMAT+ NBASM 735 III = ININ + NOCDIM ** 2 736 IJJ = III + NSPAIR 737 ICS11 = IJJ + NSPAIR 738 ICT11 = ICS11 + NPAIR2 739 IOLWS = ICT11 + NPAIR2 740 IOLWT = IOLWS + NPAIR2 741 IOLYS = IOLWT + NPAIR2 742 IOLYT = IOLYS + NPAIR2 743 IOLZS = IOLYT + NPAIR2 744 IOLZT = IOLZS + NPAIR2 745 IEVL = IOLZT + NPAIR2 746C 747 IENDW = IEVL + NSPAIR 748C 749c Elena work space allocation for r12grad 750 IF (R12PRP) THEN 751 ISB(1) = 0 752 DO ISYM = 2, NSYM 753 NBASF = NORB1(ISYM-1) + NORB2(ISYM-1) 754 NNBASF = NBASF * (NBASF + 1) / 2 755 ISB(ISYM) = ISB(ISYM-1) + NNBASF 756 ENDDO 757 NBASF = NORB1(NSYM) + NORB2(NSYM) 758 NNBASF = ISB(NSYM) + NBASF * (NBASF + 1) / 2 759 760 IGRAD = IENDW + NNBASF 761 LWRK1 = LWRK - IGRAD 762 IF (IGRAD .GT. LWRK) THEN 763 WRITE(LUPRI,'(A,I20/A,I20)') 764 * ' WORK SPACE REQUIRED = ',IGRAD, 765 * ' WORK SPACE AVAILABLE = ',LWRK 766 CALL QUIT('INSUFFICIENT WORK SPACE IN R12DRV') 767 ENDIF 768c Elena 769C 770 ELSE 771 IF (IENDW .GT. LWRK) THEN 772 WRITE(LUPRI,'(A,I20/A,I20)') 773 * ' WORK SPACE REQUIRED = ',IENDW, 774 * ' WORK SPACE AVAILABLE = ',LWRK 775 CALL QUIT('INSUFFICIENT WORK SPACE IN R12DRV') 776 ENDIF 777 ENDIF 778 779 IF (R12ECO) THEN 780 CALL ECODRV(V2AM,R2AM,EV,WRK(IV11),WRK(IB11),WRK(IQ11), 781 * WRK(IVS11),WRK(IVT11),WRK(IWS11),WRK(IWT11), 782 * WRK(IQS11),WRK(IQT11), 783 * WRK(IBS11),WRK(IBT11),WRK(IRS11), 784 * WRK(IRT11),NICDIM,NOCDIM,NSPAIR,NTPAIR, 785 * WRK(IES),WRK(IET),WRK(IFS),WRK(IFT), 786 * WRK(IEVS),WRK(ININV),WRK(IENDW),WRK(ICS11), 787 * WRK(ICT11)) 788 RETURN 789 ENDIF 790 IF (.NOT. NOTONE .OR. R12OLD .OR. LGLOCAL) THEN 791 CALL AROUND('MP2-R12 ansatz := (1 - P1 - P2 + P1*P2) * R12') 792 CALL MAKEVR(V2AM,R2AM,U2AM,EV,WRK(IQQ11), 793 * WRK(IV11),WRK(IU11),WRK(IB11), 794 * WRK(IW11),WRK(IQ11),WRK(IR11), 795 * WRK(IVS11),WRK(IUS11),WRK(IBS11), 796 * WRK(IWS11),WRK(IQS11),WRK(IRS11), 797 * WRK(IVT11),WRK(IUT11),WRK(IBT11), 798 * WRK(IWT11),WRK(IQT11),WRK(IRT11), 799 * NICDIM,NOCDIM,NSPAIR,NTPAIR, 800 * WRK(IES),WRK(IET),WRK(IFS),WRK(IFT), 801 * WRK(IEVS),WRK(ININV),RXRS,RXRT,QQ2,QQ4,QQ6, 802 * R2AM,WRK(IFLMAT),WRK(ININ),WRK(III),WRK(IJJ), 803 * WRK(ICS11),WRK(ICT11),WRK(IOLWS),WRK(IOLWT), 804 * WRK(IOLYS),WRK(IOLYT),WRK(IOLZS),WRK(IOLZT), 805 * WRK(IEVL),LGLOCAL,WRK(IGRAD),WRK(IENDW),LWRK-IENDW, 806 * NACDIM,NAPAIR) 807 END IF 808 IF (R12OLD .OR. LGLOCAL .OR. NOTTWO) GOTO 999 809 IANSAZ = 2 810 FRSWRT = .TRUE. 811 CALL AROUND('MP2-R12 ansatz := (1 - O1 - O2 + O1*O2) * R12') 812 CALL MAKERV(V2AM,R2AM,U2AM,S2AM,T2AM,EV,WRK(IQQ11), 813 * WRK(IV11),WRK(IU11),WRK(IB11), 814 * WRK(IW11),WRK(IQ11),WRK(IR11), 815 * WRK(IVS11),WRK(IUS11),WRK(IBS11), 816 * WRK(IWS11),WRK(IQS11),WRK(IRS11), 817 * WRK(IVT11),WRK(IUT11),WRK(IBT11), 818 * WRK(IWT11),WRK(IQT11),WRK(IRT11), 819 * NICDIM,NOCDIM,NSPAIR,NTPAIR, 820 * WRK(IES),WRK(IET),WRK(IFS),WRK(IFT), 821 * WRK(IEVS),WRK(ININV),QQ2,QQ4,QQ6,WRK(ICS11),WRK(ICT11), 822 * IANSAZ,FRSWRT,WRK(IENDW),LWRK-IENDW, 823 * NACDIM,NAPAIR) 824 999 CONTINUE 825 IF (NOTTRE .OR. R12OLD .OR. LGLOCAL) RETURN 826 IANSAZ = 3 827 FRSWRT = .TRUE. 828 CALL AROUND('MP2-R12 ansatz := '// 829 * '(1 - V1*V2) * (1 - O1 - O2 + O1*O2) * R12') 830 CALL MAKERV(V2AM,R2AM,U2AM,S2AM,T2AM,EV,WRK(IQQ11), 831 * WRK(IV11),WRK(IU11),WRK(IB11), 832 * WRK(IW11),WRK(IQ11),WRK(IR11), 833 * WRK(IVS11),WRK(IUS11),WRK(IBS11), 834 * WRK(IWS11),WRK(IQS11),WRK(IRS11), 835 * WRK(IVT11),WRK(IUT11),WRK(IBT11), 836 * WRK(IWT11),WRK(IQT11),WRK(IRT11), 837 * NICDIM,NOCDIM,NSPAIR,NTPAIR, 838 * WRK(IES),WRK(IET),WRK(IFS),WRK(IFT), 839 * WRK(IEVS),WRK(ININV),QQ2,QQ4,QQ6,WRK(ICS11),WRK(ICT11), 840 * IANSAZ,FRSWRT,WRK(IENDW),LWRK-IENDW, 841 * NACDIM,NAPAIR) 842 RETURN 843 END 844C /* Deck makevr */ 845 SUBROUTINE MAKEVR(V2AM,R2AM,U2AM,EV,QQ11,V11,U11,B11,W11,Q11, 846 * R11, 847 * VS11,US11,BS11,WS11,QS11,RS11,VT11,UT11,BT11, 848 * WT11,QT11,RT11,NICDIM,NOCDIM,NSPAIR,NTPAIR, 849 * ES,ET,FS,FT,EVS,CNINV,RXRS,RXRT,QQ2,QQ4,QQ6, 850 * TIJAB,FLMAT,NIND,III,JJJ,CS11,CT11,OLWS,OLWT, 851 * OLYS,OLYT,OLZS,OLZT,EVL,LGLOCAL,GRAD, 852 * WORK,LWORK,NACDIM,NAPAIR) 853C 854C This subroutine computes the V(klmn), X(klmn), and B(klmn) 855C matrices and evaluates the MP2-R12 energies. 856C 857C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 858C Modified by Claire C.M. Samson (University of Karlsruhe, 28 April 2003) 859C for local orbitals. 860C 861#include "implicit.h" 862#include "priunit.h" 863 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 864 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0, 865 * THRES = 1.0D-10) 866 DIMENSION R2AM(*), V2AM(*), U2AM(*), EV(*) 867 DIMENSION TIJAB(*), FLMAT(*), WORK(*), GRAD(*) 868 DIMENSION NIND(NOCDIM,NOCDIM), III(NSPAIR), JJJ(NSPAIR) 869 REAL*8 R1, R2, V1, U1, U2, VR, UR, RR, 870 * V11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 871 * U11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 872 * B11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 873 * W11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 874 * Q11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 875 * R11(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 876 DIMENSION QQ2(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 877 * QQ4(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 878 * QQ6(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 879 DIMENSION VS11(NSPAIR,NSPAIR) 880 DIMENSION US11(NSPAIR,NSPAIR) 881 DIMENSION BS11(NSPAIR,NSPAIR) 882 DIMENSION WS11(NSPAIR,NSPAIR) 883 DIMENSION QS11(NSPAIR,NSPAIR) 884 DIMENSION RS11(NSPAIR,NSPAIR) 885 DIMENSION VT11(NSPAIR,NSPAIR),WT11(NSPAIR,NSPAIR) 886 DIMENSION UT11(NSPAIR,NSPAIR),BT11(NSPAIR,NSPAIR) 887 DIMENSION QT11(NSPAIR,NSPAIR),RT11(NSPAIR,NSPAIR) 888 DIMENSION RXRS(NSPAIR,NSPAIR),RXRT(NSPAIR,NSPAIR) 889 DIMENSION ES(NSPAIR),ET(NSPAIR),FS(NSPAIR),FT(NSPAIR) 890 DIMENSION EVS(NSPAIR),CNINV(NSPAIR,10) 891 DIMENSION QQ11(NICDIM,NICDIM,NICDIM,NICDIM) 892 DIMENSION CS11(NSPAIR,NSPAIR),CT11(NSPAIR,NSPAIR) 893 DIMENSION OLWS(NSPAIR,NSPAIR),OLWT(NSPAIR,NSPAIR) 894 DIMENSION OLYS(NSPAIR,NSPAIR),OLYT(NSPAIR,NSPAIR) 895 DIMENSION OLZS(NSPAIR,NSPAIR),OLZT(NSPAIR,NSPAIR) 896 DIMENSION EVL(NSPAIR) 897 INTEGER IDUM,LWORK 898 LOGICAL LDUM 899 LOGICAL LGLOCAL, LOCRST 900 CHARACTER*11 CFN 901 CHARACTER*3 APPROX(0:7,2), IPCC 902 INTEGER NKILJ(8) 903#include "r12int.h" 904#include "ccr12int.h" 905#include "ccorb.h" 906#include "ccsdsym.h" 907#include "ccsdinp.h" 908 DATA APPROX 909 &/"0 ","A ","A ","A' ","B ","B ","xxx","xxx", 910 & "0 ","A ","A ","A' ","B' ","B' ","B ","B "/ 911 912 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 913 1001 FORMAT('CSMAT.',I2.2,'.',I2.2) 914 1002 FORMAT('CTMAT.',I2.2,'.',I2.2) 915Ckr We need to have the next two unit numbers in consecutive order, and 916Ckr thus assign a very large unit number in order to avoid conflicts with 917Ckr already open records.\kr-231007\ 918Ckr 919 LUMULO = 90 920 LUMULN = 91 921 IF (R12CBS) THEN 922 FACX = - D1 923 ELSE 924 FACX = D1 925 END IF 926 DELTA = R12LEV 927 FF = D1 / SQRT(D2) 928 LU33 = -33 929 CALL GPOPEN(LU33,'AUXQ12','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 930 READ(LU33,'(4E30.20)') QQ11 931 CALL GPCLOSE(LU33,'KEEP') 932 DO NPIJ = 1, NSPAIR 933 IJ = 0 934 DO I = 1, NOCDIM 935 DO J = 1, I 936 IJ = IJ + 1 937 IF (IJ.EQ.NPIJ) THEN 938 III(NPIJ) = I 939 JJJ(NPIJ) = J 940 END IF 941 END DO 942 END DO 943 END DO 944 DO I = 1, NOCDIM 945 DO J = 1, NOCDIM 946 IF (I.GT.J) THEN 947 NIND(I,J) = I * (I-1) /2 + J 948 ELSE 949 NIND(I,J) = J * (J-1) /2 + I 950 END IF 951 END DO 952 END DO 953 IJ = 0 954 DO I = 1, NOCDIM 955 DO J = 1, I 956 IJ = IJ + 1 957 EVS(IJ) = EV(I) + EV(J) 958 END DO 959 END DO 960 IF (LGLOCAL) THEN 961 NCT = 0 962 NRT = 0 963 DO ISYM = 1, NSYM 964 NCT = NCT + NRHFFR(ISYM) 965 NRT = NRT + NORB1(ISYM) 966 END DO 967 LU99 = -99 968 CALL GPOPEN(LU99,'FLOCA','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 969 NBASM = (NRT-NCT)*(NRT-NCT+1)/2 970 DO I=1,NBASM 971 READ(LU99,'(D30.20)') FLMAT(I) 972 END DO 973 CALL GPCLOSE(LU99,'KEEP') 974 DO IJ = 1, NSPAIR 975 NI = III(IJ) 976 NJ = JJJ(IJ) 977 NFI = NIND(NI,NI) 978 NFJ = NIND(NJ,NJ) 979 EVL(IJ) = FLMAT(NFI) + FLMAT(NFJ) 980 END DO 981 ELSE 982 NFMAT = NOCDIM*(NOCDIM+1)/2 983 CALL DZERO(FLMAT,NFMAT) 984 DO IJ = 1, NSPAIR 985 NI = III(IJ) 986 NJ = JJJ(IJ) 987 NFI = NIND(NI,NI) 988 EVL(IJ) = EVS(IJ) 989 IF (NI.EQ.NJ) FLMAT(NFI) = EVL(IJ) * DP5 990 END DO 991 END IF 992 IF (IPRINT .GE. 4 .OR. LGLOCAL) THEN 993 CALL AROUND('Fock matrix in occupied space') 994 WRITE(LUPRI,*) 995 CALL OUTPAK(FLMAT,NOCDIM,1,LUPRI) 996 CALL AROUND('Pairs of orbital energies') 997 CALL OUTPUT(EVL,1,NSPAIR,1,1,NSPAIR,1,1,LUPRI) 998 END IF 999C 1000 DO 205 I = 1, NOCDIM**4 1001 R11(I,1,1,1) = D0 1002 Q11(I,1,1,1) = D0 1003 V11(I,1,1,1) = D0 1004 U11(I,1,1,1) = D0 1005 B11(I,1,1,1) = D0 1006 W11(I,1,1,1) = D0 1007 205 CONTINUE 1008 CALL DZERO(ES,NSPAIR) 1009 CALL DZERO(ET,NSPAIR) 1010C 1011C PRINT ORBITAL ENERGIES 1012C 1013 IF (IPRINT .LE. 5) GOTO 993 1014 DO 51 ISYM = 1, NSYM 1015 WRITE(LUPRI,'(/A,I2)') ' OCCUPIED ORBITALS OF SYMMETRY',ISYM 1016 DO 52 I = 1, NRHF(ISYM) 1017 KOFFI = IRHF(ISYM) + I 1018 WRITE(LUPRI,'(2I5,G20.10)') I,KOFFI,EV(KOFFI) 1019 52 CONTINUE 1020 WRITE(LUPRI,'(/A,I2)') ' PRIMARY ORBITALS OF SYMMETRY',ISYM 1021 DO 53 I = 1, NORB1(ISYM) 1022 KOFFI = IVIR(ISYM) + I 1023 WRITE(LUPRI,'(2I5,G20.10)') I,KOFFI,EV(KOFFI) 1024 53 CONTINUE 1025 WRITE(LUPRI,'(/A,I2)') ' SECONDARY ORBITALS OF SYMMETRY',ISYM 1026 DO 54 I = NORB1(ISYM) + 1, NVIR(ISYM) 1027 KOFFI = IVIR(ISYM) + I 1028 WRITE(LUPRI,'(2I5,G20.10)') I,KOFFI,EV(KOFFI) 1029 54 CONTINUE 1030 51 CONTINUE 1031 993 CONTINUE 1032C 1033C COMPUTE ITERATIVELY THE AMPLITUDES 1034C 1035 IF (LGLOCAL) THEN 1036 ICOUNT=1 1037 CALL DZERO(TIJAB,NH2AMX) 1038 900 ICONV=0 1039 E2 = D0 1040 DO 701 ISYMIA = 1, NSYM 1041 ISYMJB = ISYMIA 1042 DO 702 ISYMI = 1, NSYM 1043 ISYMA = MULD2H(ISYMI,ISYMIA) 1044 DO 703 ISYMJ = 1, NSYM 1045 ISYMB = MULD2H(ISYMJ,ISYMJB) 1046 DO 801 I = 1, NRHF(ISYMI) 1047 KOFFI = IRHF(ISYMI) + I 1048 DO 802 J = 1, NRHF(ISYMJ) 1049 KOFFJ = IRHF(ISYMJ) + J 1050 DO 803 A = NRHFS(ISYMA) + 1, NORB1(ISYMA) 1051 ISYMJA = MULD2H(ISYMJ,ISYMA) 1052 KOFFA = IVIR(ISYMA) + A 1053 DO 804 B = NRHFS(ISYMB) + 1, NORB1(ISYMB) 1054 ISYMIB = MULD2H(ISYMI,ISYMB) 1055 KOFFB = IVIR(ISYMB) + B 1056 IF (ONEAUX) THEN 1057 NAI = IH1AM(ISYMA,ISYMI) + 1058 * NORB1(ISYMA)*(I-1) + A 1059 NBJ = IH1AM(ISYMB,ISYMJ) + 1060 * NORB1(ISYMB)*(J-1) + B 1061 NAIBJ = IH2AM(ISYMIA,ISYMJB) + 1062 * INDEX(NAI,NBJ) 1063 NAJ = IH1AM(ISYMA,ISYMJ) + 1064 * NORB1(ISYMA)*(J-1) + A 1065 NBI = IH1AM(ISYMB,ISYMI) + 1066 * NORB1(ISYMB)*(I-1) + B 1067 NAJBI = IH2AM(ISYMJA,ISYMIB) + 1068 * INDEX(NAJ,NBI) 1069 ELSE 1070 NAI = IT1AM(ISYMA,ISYMI) + 1071 * NVIR(ISYMA)*(I-1) + A 1072 NBJ = IT1AM(ISYMB,ISYMJ) + 1073 * NVIR(ISYMB)*(J-1) + B 1074 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 1075 * INDEX(NAI,NBJ) 1076 NAJ = IT1AM(ISYMA,ISYMJ) + 1077 * NVIR(ISYMA)*(J-1) + A 1078 NBI = IT1AM(ISYMB,ISYMI) + 1079 * NVIR(ISYMB)*(I-1) + B 1080 NAJBI = IT2AM(ISYMJA,ISYMIB) + 1081 * INDEX(NAJ,NBI) 1082 END IF 1083 VAIBJ = V2AM(NAIBJ) 1084 VAJBI = V2AM(NAJBI) 1085 VV = VAIBJ 1086C 1087 TKI = D0 1088 TKJ = D0 1089 TKA = D0 1090 TKB = D0 1091 TINI = TIJAB(NAIBJ) 1092C 1093C LOOP FOR COEFFICIENTS (KJAB) 1094C 1095 ISYMKA = ISYMJB 1096 ISYMK = MULD2H(ISYMA,ISYMKA) 1097 DO 902 K = 1, NRHF(ISYMK) 1098 IF ((K.EQ.I).AND.(ISYMI.EQ.ISYMK)) 1099 * GOTO 902 1100 IF (ONEAUX) THEN 1101 NAK = IH1AM(ISYMA,ISYMK) + 1102 * NORB1(ISYMA)*(K-1) + A 1103 NAKBJ = IH2AM(ISYMKA,ISYMJB) + 1104 * INDEX(NAK,NBJ) 1105 ELSE 1106 NAK = IT1AM(ISYMA,ISYMK) + 1107 * NVIR(ISYMA)*(K-1) + A 1108 NAKBJ = IT2AM(ISYMKA,ISYMJB) + 1109 * INDEX(NAK,NBJ) 1110 END IF 1111 NIK = INDEX(I,K) 1112 XXX = FLMAT(NIK) * TIJAB(NAKBJ) 1113 TKI = TKI + XXX 1114 902 CONTINUE 1115C 1116C LOOP FOR COEFFICIENTS (IKAB) 1117C 1118 ISYMKB = ISYMIA 1119 ISYMK = MULD2H(ISYMB,ISYMKB) 1120 DO 904 K = 1, NRHF(ISYMK) 1121 IF ((K.EQ.J).AND.(ISYMK.EQ.ISYMJ)) 1122 * GOTO 904 1123 IF (ONEAUX) THEN 1124 NBK = IH1AM(ISYMB,ISYMK) + 1125 * NORB1(ISYMB)*(K-1) + B 1126 NAIBK = IH2AM(ISYMKB,ISYMIA) + 1127 * INDEX(NBK,NAI) 1128 ELSE 1129 NBK = IT1AM(ISYMB,ISYMK) + 1130 * NVIR(ISYMB)*(K-1) + B 1131 NAIBK = IT2AM(ISYMKB,ISYMIA) + 1132 * INDEX(NBK,NAI) 1133 END IF 1134 NJK = INDEX(J,K) 1135 XXX = FLMAT(NJK)*TIJAB(NAIBK) 1136 TKJ = TKJ + XXX 1137 904 CONTINUE 1138C 1139C LOOP FOR COEFFICIENTS (IJCB) 1140C 1141 ISYMIC = ISYMJB 1142 ISYMC = MULD2H(ISYMI,ISYMIC) 1143 DO 906 NC = NRHFS(ISYMC) + 1, NORB1(ISYMC) 1144 IF ((NC.EQ.A).AND.(ISYMA.EQ.ISYMC)) 1145 * GOTO 906 1146 IF (ONEAUX) THEN 1147 NCI = IH1AM(ISYMC,ISYMI) + 1148 * NORB1(ISYMC)*(I-1) + NC 1149 NCIBJ = IH2AM(ISYMIC,ISYMJB) + 1150 * INDEX(NCI,NBJ) 1151 ELSE 1152 NCI = IT1AM(ISYMC,ISYMI) + 1153 * NVIR(ISYMC)*(I-1) + NC 1154 NCIBJ = IT2AM(ISYMIC,ISYMJB) + 1155 * INDEX(NCI,NBJ) 1156 END IF 1157 NCA = INDEX(NC-NCT,A-NCT) 1158 XXX = FLMAT(NCA)*TIJAB(NCIBJ) 1159 TKA = TKA + XXX 1160 906 CONTINUE 1161C 1162C LOOP FOR COEFFICIENTS (IJAC) 1163C 1164 ISYMJC = ISYMIA 1165 ISYMC = MULD2H(ISYMJ,ISYMJC) 1166 DO 908 NC = NRHFS(ISYMC) + 1, NORB1(ISYMC) 1167 IF ((NC.EQ.B).AND.(ISYMB.EQ.ISYMC)) 1168 * GOTO 908 1169 IF (ONEAUX) THEN 1170 NCJ = IH1AM(ISYMC,ISYMJ) + 1171 * NORB1(ISYMC)*(J-1) + NC 1172 NAICJ = IH2AM(ISYMJC,ISYMIA) + 1173 * INDEX(NCJ,NAI) 1174 ELSE 1175 NCJ = IT1AM(ISYMC,ISYMJ) + 1176 * NVIR(ISYMC)*(J-1) + NC 1177 NAICJ = IT2AM(ISYMJC,ISYMIA) + 1178 * INDEX(NCJ,NAI) 1179 END IF 1180 NCB = INDEX(NC-NCT,B-NCT) 1181 XXX = FLMAT(NCB)*TIJAB(NAICJ) 1182 TKB = TKB + XXX 1183 908 CONTINUE 1184C 1185 NAA = INDEX(A-NCT,A-NCT) 1186 NBB = INDEX(B-NCT,B-NCT) 1187 NII = INDEX(I,I) 1188 NJJ = INDEX(J,J) 1189 DENOM = D1 / (FLMAT(NAA) + FLMAT(NBB) - 1190 * FLMAT(NII) - FLMAT(NJJ)) 1191 TIJAB(NAIBJ) = (TKI+TKJ-TKA-TKB-VV) 1192 * * DENOM 1193 VV = (D2 * VAIBJ - VAJBI) 1194 E2 = E2 + VV * TIJAB(NAIBJ) 1195 IF (ABS(TIJAB(NAIBJ)-TINI).GT.THRES) 1196 * ICONV = 1 1197 804 CONTINUE 1198 803 CONTINUE 1199 802 CONTINUE 1200 801 CONTINUE 1201 703 CONTINUE 1202 702 CONTINUE 1203 701 CONTINUE 1204 IF (ICOUNT .EQ. 1) WRITE(LUPRI,*) 1205 WRITE(LUPRI,'(A,I3,3X,F15.9)') ' ENERGY IN ITERATION', 1206 & ICOUNT,E2 1207 IF (ICOUNT.LE.50) THEN 1208 IF (ICONV.EQ.1) THEN 1209 ICOUNT= ICOUNT + 1 1210 GOTO 900 1211 ELSE 1212 WRITE(LUPRI,*) 1213 WRITE(LUPRI,'(A,I3)') ' FINAL ITERATION OF LOCAL MP2', 1214 & ICOUNT 1215 END IF 1216 ELSE 1217 WRITE(LUPRI,'(A)') 'MORE THAN 50 ITERATIONS' 1218 END IF 1219 END IF 1220C 1221C CONSTRUCT PRODUCTS (IA|JB)*(KA|LB) 1222C 1223 E2 = D0 1224 E2S = D0 1225 E2T = D0 1226 DO 101 ISYMIA = 1, NSYM 1227 ISYMJB = ISYMIA 1228 JBOFF = NH1AM(ISYMJB) * (NH1AM(ISYMJB) + 1) / 2 1229 DO 102 ISYMI = 1, NSYM 1230 ISYMA = MULD2H(ISYMI,ISYMIA) 1231 DO 103 ISYMJ = 1, NSYM 1232 ISYMB = MULD2H(ISYMJ,ISYMJB) 1233C 1234C COMPUTE MP2 ENERGY 1235C 1236 DO 201 I = 1, NRHFA(ISYMI) 1237 KOFFI = IRHF(ISYMI) + I 1238 KOFFIA = IRHFA(ISYMI) + I 1239 DO 202 J = 1, NRHFA(ISYMJ) 1240 KOFFJ = IRHF(ISYMJ) + J 1241 KOFFJA = IRHFA(ISYMJ) + J 1242 DO 203 A = NRHFSA(ISYMA) + 1, NORB1(ISYMA) 1243 ISYMJA = MULD2H(ISYMJ,ISYMA) 1244 KOFFA = IVIR(ISYMA) + A 1245 DO 204 B = NRHFSA(ISYMB) + 1, NORB1(ISYMB) 1246 ISYMIB = MULD2H(ISYMI,ISYMB) 1247 KOFFB = IVIR(ISYMB) + B 1248 IF (ONEAUX) THEN 1249 NAI = IH1AM(ISYMA,ISYMI) + 1250 * NORB1(ISYMA)*(I-1) + A 1251 NBJ = IH1AM(ISYMB,ISYMJ) + 1252 * NORB1(ISYMB)*(J-1) + B 1253 NAIBJ = IH2AM(ISYMIA,ISYMJB) + 1254 * INDEX(NAI,NBJ) 1255 NAJ = IH1AM(ISYMA,ISYMJ) + 1256 * NORB1(ISYMA)*(J-1) + A 1257 NBI = IH1AM(ISYMB,ISYMI) + 1258 * NORB1(ISYMB)*(I-1) + B 1259 NAJBI = IH2AM(ISYMJA,ISYMIB) + 1260 * INDEX(NAJ,NBI) 1261 ELSE 1262 NAI = IT1AM(ISYMA,ISYMI) + 1263 * NVIR(ISYMA)*(I-1) + A 1264 NBJ = IT1AM(ISYMB,ISYMJ) + 1265 * NVIR(ISYMB)*(J-1) + B 1266 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 1267 * INDEX(NAI,NBJ) 1268 NAJ = IT1AM(ISYMA,ISYMJ) + 1269 * NVIR(ISYMA)*(J-1) + A 1270 NBI = IT1AM(ISYMB,ISYMI) + 1271 * NVIR(ISYMB)*(I-1) + B 1272 NAJBI = IT2AM(ISYMJA,ISYMIB) + 1273 * INDEX(NAJ,NBI) 1274 END IF 1275 VAIBJ = V2AM(NAIBJ) 1276 VAJBI = V2AM(NAJBI) 1277C 1278 IF (LGLOCAL) THEN 1279C LOCALIZED MP2 1280 VV = (D2 * VAIBJ - VAJBI) 1281 DENOM = TIJAB(NAIBJ) 1282 VS = (VAIBJ + VAJBI) 1283 VT = (VAIBJ - VAJBI) 1284 ELSE 1285C CANONICAL MP2 1286 VV = VAIBJ * (D2 * VAIBJ - VAJBI) 1287 DENOM = D1 / (EV(KOFFI) + EV(KOFFJ) - 1288 * EV(KOFFA) - EV(KOFFB)) 1289 VS = (VAIBJ + VAJBI)**2 1290 VT = (VAIBJ - VAJBI)**2 1291 END IF 1292 E2 = E2 + VV * DENOM 1293 IJ = INDEX(KOFFIA,KOFFJA) 1294 ES(IJ) = ES(IJ) + VS * DENOM 1295 ET(IJ) = ET(IJ) + VT * DENOM 1296 204 CONTINUE 1297 203 CONTINUE 1298 202 CONTINUE 1299 201 CONTINUE 1300 103 CONTINUE 1301 102 CONTINUE 1302 101 CONTINUE 1303 1304C 1305 IF (LGLOCAL) THEN 1306C 1307C WRITE AMPLITUDES AND RESTORE R12 INTEGRALS 1308C 1309 LU43 = -43 1310 CALL GPOPEN(LU43,'TIJAB', 1311 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 1312 WRITE(LU43) (R2AM(I), I = 1,NH2AMX) 1313 CALL GPCLOSE(LU43,'KEEP') 1314 CALL GPOPEN(LU43,FR12R12, 1315 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 1316 READ(LU43) (R2AM(I), I = 1,NH2AMX) 1317 CALL GPCLOSE(LU43,'KEEP') 1318C 1319 END IF 1320C 1321 DO 1101 ISYMIA = 1, NSYM 1322 ISYMJB = ISYMIA 1323 JBOFF = NH1AM(ISYMJB) * (NH1AM(ISYMJB) + 1) / 2 1324 DO 1102 ISYMI = 1, NSYM 1325 ISYMA = MULD2H(ISYMI,ISYMIA) 1326 DO 1103 ISYMJ = 1, NSYM 1327 ISYMB = MULD2H(ISYMJ,ISYMJB) 1328C 1329 IF (ONEAUX) THEN 1330 DO 1104 ISYMKA = 1, NSYM 1331 ISYMLB = ISYMKA 1332 ISYMK = MULD2H(ISYMA,ISYMKA) 1333 ISYML = MULD2H(ISYMB,ISYMLB) 1334 DO 1105 I = 1, NRHF(ISYMI) 1335 KOFFI = IRHF(ISYMI) + I 1336 DO 1106 K = 1, NRHF(ISYMK) 1337 KOFFK = IRHF(ISYMK) + K 1338 DO 1107 A = 1, NORB1(ISYMA) 1339 NAI = IH1AM(ISYMA,ISYMI) + 1340 * NORB1(ISYMA)*(I-1) + A 1341 NAK = IH1AM(ISYMA,ISYMK) + 1342 * NORB1(ISYMA)*(K-1) + A 1343 DO 1108 J = 1, NRHF(ISYMJ) 1344 KOFFJ = IRHF(ISYMJ) + J 1345 DO 1109 L = 1, NRHF(ISYML) 1346 KOFFL = IRHF(ISYML) + L 1347 DO 1110 B = 1, NORB1(ISYMB) 1348 NBJ = IH1AM(ISYMB,ISYMJ) + 1349 * NORB1(ISYMB)*(J-1) + B 1350 NBL = IH1AM(ISYMB,ISYML) + 1351 * NORB1(ISYMB)*(L-1) + B 1352 NAIBJ = IH2AM(ISYMIA,ISYMJB) + 1353 * INDEX(NAI,NBJ) 1354 NAKBL = IH2AM(ISYMKA,ISYMLB) + 1355 * INDEX(NAK,NBL) 1356 R1 = R2AM(NAIBJ) 1357 V1 = V2AM(NAIBJ) 1358 U1 = U2AM(NAIBJ) 1359 R2 = R2AM(NAKBL) 1360 U2 = U2AM(NAKBL) 1361 VR = V1 * R2 1362 UR = U1 * R2 + U2 * R1 1363 RR = R1 * R2 1364 V11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1365 * V11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 1366 U11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1367 * U11(KOFFI,KOFFJ,KOFFK,KOFFL) + UR 1368 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1369 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 1370 1110 CONTINUE 1371 1109 CONTINUE 1372 1108 CONTINUE 1373 1107 CONTINUE 1374 1106 CONTINUE 1375 1105 CONTINUE 1376 1104 CONTINUE 1377 DO 2104 ISYMKA = 1, NSYM 1378 ISYMLB = ISYMKA 1379 LBOFF = NH1AM(ISYMLB) * (NH1AM(ISYMLB) + 1) / 2 1380 ISYMK = MULD2H(ISYMA,ISYMKA) 1381 ISYML = MULD2H(ISYMB,ISYMLB) 1382 DO 2105 I = 1, NRHF(ISYMI) 1383 KOFFI = IRHF(ISYMI) + I 1384 DO 2106 K = 1, NRHF(ISYMK) 1385 KOFFK = IRHF(ISYMK) + K 1386 DO 2107 A = 1, NORB2(ISYMA) 1387 NAI = IG1AM(ISYMA,ISYMI) + 1388 * NORB2(ISYMA)*(I-1) + A 1389 NAK = IG1AM(ISYMA,ISYMK) + 1390 * NORB2(ISYMA)*(K-1) + A 1391 DO 2108 J = 1, NRHF(ISYMJ) 1392 KOFFJ = IRHF(ISYMJ) + J 1393 DO 2109 L = 1, NRHF(ISYML) 1394 KOFFL = IRHF(ISYML) + L 1395 DO 2110 B = 1, NORB1(ISYMB) 1396 NBJ = IH1AM(ISYMB,ISYMJ) + 1397 * NORB1(ISYMB)*(J-1) + B 1398 NBL = IH1AM(ISYMB,ISYML) + 1399 * NORB1(ISYMB)*(L-1) + B 1400 NAIBJ = IH2AM(ISYMIA,ISYMJB) + 1401 * NH1AM(ISYMJB)*(NAI - 1) + NBJ + 1402 * JBOFF 1403 NAKBL = IH2AM(ISYMKA,ISYMLB) + 1404 * NH1AM(ISYMLB)*(NAK - 1) + NBL + 1405 * LBOFF 1406 R1 = R2AM(NAIBJ) 1407 V1 = V2AM(NAIBJ) 1408 U1 = U2AM(NAIBJ) 1409 R2 = R2AM(NAKBL) 1410 U2 = U2AM(NAKBL) 1411 VR = V1 * R2 1412 UR = U1 * R2 + U2 * R1 1413 RR = R1 * R2 1414 W11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1415 * W11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 1416 B11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1417 * B11(KOFFI,KOFFJ,KOFFK,KOFFL) + UR 1418 R11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1419 * R11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 1420 W11(KOFFJ,KOFFI,KOFFL,KOFFK) = 1421 * W11(KOFFJ,KOFFI,KOFFL,KOFFK) + VR 1422 B11(KOFFJ,KOFFI,KOFFL,KOFFK) = 1423 * B11(KOFFJ,KOFFI,KOFFL,KOFFK) + UR 1424 R11(KOFFJ,KOFFI,KOFFL,KOFFK) = 1425 * R11(KOFFJ,KOFFI,KOFFL,KOFFK) + RR 1426 2110 CONTINUE 1427 2109 CONTINUE 1428 2108 CONTINUE 1429 2107 CONTINUE 1430 2106 CONTINUE 1431 2105 CONTINUE 1432 2104 CONTINUE 1433 ELSE 1434 DO 104 ISYMKA = 1, NSYM 1435 ISYMLB = ISYMKA 1436 ISYMK = MULD2H(ISYMA,ISYMKA) 1437 ISYML = MULD2H(ISYMB,ISYMLB) 1438 DO 105 I = 1, NRHF(ISYMI) 1439 KOFFI = IRHF(ISYMI) + I 1440 DO 106 K = 1, NRHF(ISYMK) 1441 KOFFK = IRHF(ISYMK) + K 1442 DO 107 A = 1, NVIR(ISYMA) 1443 NAI = IT1AM(ISYMA,ISYMI) + 1444 * NVIR(ISYMA)*(I-1) + A 1445 NAK = IT1AM(ISYMA,ISYMK) + 1446 * NVIR(ISYMA)*(K-1) + A 1447 DO 108 J = 1, NRHF(ISYMJ) 1448 KOFFJ = IRHF(ISYMJ) + J 1449 DO 109 L = 1, NRHF(ISYML) 1450 KOFFL = IRHF(ISYML) + L 1451 DO 110 B = 1, NVIR(ISYMB) 1452 NBJ = IT1AM(ISYMB,ISYMJ) + 1453 * NVIR(ISYMB)*(J-1) + B 1454 NBL = IT1AM(ISYMB,ISYML) + 1455 * NVIR(ISYMB)*(L-1) + B 1456 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 1457 * INDEX(NAI,NBJ) 1458 NAKBL = IT2AM(ISYMKA,ISYMLB) + 1459 * INDEX(NAK,NBL) 1460 R1 = R2AM(NAIBJ) 1461 V1 = V2AM(NAIBJ) 1462 U1 = U2AM(NAIBJ) 1463 R2 = R2AM(NAKBL) 1464 U2 = U2AM(NAKBL) 1465 VR = V1 * R2 1466 UR = U1 * R2 + U2 * R1 1467 RR = R1 * R2 1468C 1469 IF (A .GT. NORB1(ISYMA)) GOTO 112 1470 IF (B .GT. NORB1(ISYMB)) GOTO 112 1471C 1472 V11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1473 * V11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 1474 U11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1475 * U11(KOFFI,KOFFJ,KOFFK,KOFFL) + UR 1476 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1477 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 1478 GOTO 111 1479 112 CONTINUE 1480C 1481 IF (A .GT. NORB1(ISYMA) .AND. 1482 * B .GT. NORB1(ISYMB)) GOTO 111 1483C 1484 W11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1485 * W11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 1486 B11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1487 * B11(KOFFI,KOFFJ,KOFFK,KOFFL) + UR 1488 R11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1489 * R11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 1490C 1491 111 CONTINUE 1492 110 CONTINUE 1493 109 CONTINUE 1494 108 CONTINUE 1495 107 CONTINUE 1496 106 CONTINUE 1497 105 CONTINUE 1498 104 CONTINUE 1499 ENDIF 1500 1103 CONTINUE 1501 1102 CONTINUE 1502 1101 CONTINUE 1503C 1504 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER PAIR ENERGIES:' 1505 DO 230 I=1,NAPAIR 1506 IF (LGLOCAL) THEN 1507 ES(I) = ES(I) * 0.5D0 1508 ET(I) = ET(I) * 1.5D0 1509 ELSE 1510 ES(I) = ES(I) * DP25 1511 ET(I) = ET(I) * DP75 1512 END IF 1513 E2S = E2S + ES(I) 1514 E2T = E2T + ET(I) 1515 WRITE(LUPRI,'(17X,I4,2F15.9)') I,ES(I),ET(I) 1516 230 CONTINUE 1517 WRITE(LUPRI,'(/A6,3F15.9)') ' MP2 =',E2,E2S,E2T 1518C 1519 IJ = 0 1520 DO 300 I = 1, NOCDIM 1521 DO 301 J = 1, I 1522 IJ = IJ + 1 1523 KL = 0 1524 DO 302 K = 1, NOCDIM 1525 DO 303 L = 1, K 1526 KL = KL + 1 1527C 1528 IF (R12EOR) THEN 1529 RR = - D2 *(QQ6(I,J,K,L) + QQ6(I,J,L,K)) 1530 RR = RR - (U11(I,J,K,L) + U11(I,J,L,K)) 1531 DX = D0 1532 ELSE 1533 RR = - (U11(I,J,K,L) + U11(I,J,L,K)) 1534 DX = D2 1535 END IF 1536 US11(KL,IJ) = RR 1537 IF (I .EQ. J) US11(KL,IJ) = FF * US11(KL,IJ) 1538 IF (K .EQ. L) US11(KL,IJ) = FF * US11(KL,IJ) 1539 IF (IJ .EQ. KL) US11(KL,IJ) = US11(KL,IJ) - DX 1540 IF (R12EOR) THEN 1541 RR = - D2 *(QQ6(I,J,K,L) - QQ6(I,J,L,K)) 1542 RR = RR - (U11(I,J,K,L) - U11(I,J,L,K)) 1543 DX = D0 1544 ELSE 1545 RR = - (U11(I,J,K,L) - U11(I,J,L,K)) 1546 DX = D2 1547 END IF 1548 UT11(KL,IJ) = RR 1549 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 1550 * UT11(KL,IJ) = UT11(KL,IJ) - DX 1551C 1552 IF (R12EOR) THEN 1553 RR = (QQ2(I,J,K,L) + QQ2(I,J,L,K)) 1554 RR = RR - (V11(I,J,K,L) + V11(I,J,L,K)) 1555 DX = D0 1556 ELSE 1557 RR = - (V11(I,J,K,L) + V11(I,J,L,K)) 1558 DX = D1 1559 END IF 1560 VS11(KL,IJ) = RR 1561 IF (I .EQ. J) VS11(KL,IJ) = FF * VS11(KL,IJ) 1562 IF (K .EQ. L) VS11(KL,IJ) = FF * VS11(KL,IJ) 1563 IF (IJ .EQ. KL) VS11(KL,IJ) = VS11(KL,IJ) + DX 1564 IF (R12EOR) THEN 1565 RR = (QQ2(I,J,K,L) - QQ2(I,J,L,K)) 1566 RR = RR - (V11(I,J,K,L) - V11(I,J,L,K)) 1567 DX = D0 1568 ELSE 1569 RR = - (V11(I,J,K,L) - V11(I,J,L,K)) 1570 DX = D1 1571 END IF 1572 VT11(KL,IJ) = RR 1573 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 1574 * VT11(KL,IJ) = VT11(KL,IJ) + DX 1575C 1576 IF (R12EOR) THEN 1577 RR = - D2 *(QQ6(I,J,K,L) + QQ6(I,J,L,K)) 1578 RR = RR + (U11(I,J,K,L) + U11(I,J,L,K)) * FACX 1579 RR = RR - (B11(I,J,K,L) + B11(I,J,L,K)) 1580 DX = D0 1581 ELSE 1582 RR = (U11(I,J,K,L) + U11(I,J,L,K)) * FACX 1583 * - (B11(I,J,K,L) + B11(I,J,L,K)) 1584 DX = D2 1585 END IF 1586 BS11(KL,IJ) = RR 1587 IF (I .EQ. J) BS11(KL,IJ) = FF * BS11(KL,IJ) 1588 IF (K .EQ. L) BS11(KL,IJ) = FF * BS11(KL,IJ) 1589 IF (IJ .EQ. KL) BS11(KL,IJ) = BS11(KL,IJ) - DX 1590 IF (R12EOR) THEN 1591 RR = - D2 *(QQ6(I,J,K,L) - QQ6(I,J,L,K)) 1592 RR = RR + (U11(I,J,K,L) - U11(I,J,L,K)) * FACX 1593 RR = RR - (B11(I,J,K,L) - B11(I,J,L,K)) 1594 DX = D0 1595 ELSE 1596 RR = (U11(I,J,K,L) - U11(I,J,L,K)) * FACX 1597 * - (B11(I,J,K,L) - B11(I,J,L,K)) 1598 DX = D2 1599 END IF 1600 BT11(KL,IJ) = RR 1601 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 1602 * BT11(KL,IJ) = BT11(KL,IJ) - DX 1603C 1604 IF (R12EOR) THEN 1605 RR = (QQ2(I,J,K,L) + QQ2(I,J,L,K)) 1606 RR = RR + (V11(I,J,K,L) + V11(I,J,L,K)) * FACX 1607 RR = RR - (W11(I,J,K,L) + W11(I,J,L,K)) 1608 DX = D0 1609 ELSE 1610 RR = (V11(I,J,K,L) + V11(I,J,L,K)) * FACX 1611 * - (W11(I,J,K,L) + W11(I,J,L,K)) 1612 DX = D1 1613 END IF 1614 WS11(KL,IJ) = RR 1615 IF (I .EQ. J) WS11(KL,IJ) = FF * WS11(KL,IJ) 1616 IF (K .EQ. L) WS11(KL,IJ) = FF * WS11(KL,IJ) 1617 IF (IJ .EQ. KL) WS11(KL,IJ) = WS11(KL,IJ) + DX 1618 IF (R12EOR) THEN 1619 RR = (QQ2(I,J,K,L) - QQ2(I,J,L,K)) 1620 RR = RR + (V11(I,J,K,L) - V11(I,J,L,K)) * FACX 1621 RR = RR - (W11(I,J,K,L) - W11(I,J,L,K)) 1622 DX = D0 1623 ELSE 1624 RR = (V11(I,J,K,L) - V11(I,J,L,K)) * FACX 1625 * - (W11(I,J,K,L) - W11(I,J,L,K)) 1626 DX = D1 1627 END IF 1628 WT11(KL,IJ) = RR 1629 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 1630 * WT11(KL,IJ) = WT11(KL,IJ) + DX 1631C 1632 IF (R12EOR) THEN 1633 QQIJKL = QQ4(I,J,K,L) 1634 QQIJLK = QQ4(I,J,L,K) 1635 ELSE 1636 QQIJKL = QQ11(IQ(I),IQ(J),IQ(K),IQ(L)) 1637 QQIJLK = QQ11(IQ(I),IQ(J),IQ(L),IQ(K)) 1638 END IF 1639 RR = QQIJKL + QQIJLK 1640 RR = RR - (Q11(I,J,K,L) + Q11(I,J,L,K)) 1641 QS11(KL,IJ) = RR 1642 IF (I .EQ. J) QS11(KL,IJ) = FF * QS11(KL,IJ) 1643 IF (K .EQ. L) QS11(KL,IJ) = FF * QS11(KL,IJ) 1644 RR = QQIJKL - QQIJLK 1645 RR = RR - (Q11(I,J,K,L) - Q11(I,J,L,K)) 1646 QT11(KL,IJ) = RR 1647C 1648 RR = QQIJKL + QQIJLK 1649 RR = RR + (Q11(I,J,K,L) + Q11(I,J,L,K)) * FACX 1650 * - (R11(I,J,K,L) + R11(I,J,L,K)) 1651ccn rr = + (Q11(I,J,K,L) + Q11(I,J,L,K)) * FACX 1652 RS11(KL,IJ) = RR 1653 IF (I .EQ. J) RS11(KL,IJ) = FF * RS11(KL,IJ) 1654 IF (K .EQ. L) RS11(KL,IJ) = FF * RS11(KL,IJ) 1655 RR = QQIJKL - QQIJLK 1656 RR = RR + (Q11(I,J,K,L) - Q11(I,J,L,K)) * FACX 1657 * - (R11(I,J,K,L) - R11(I,J,L,K)) 1658ccn rr = + (Q11(I,J,K,L) - Q11(I,J,L,K)) * FACX 1659 RT11(KL,IJ) = RR 1660C 1661 US11(KL,IJ) = US11(KL,IJ) * DP5 * DP25 1662 BS11(KL,IJ) = BS11(KL,IJ) * DP5 * DP25 1663 VS11(KL,IJ) = VS11(KL,IJ) * DP5 1664 WS11(KL,IJ) = WS11(KL,IJ) * DP5 1665 QS11(KL,IJ) = QS11(KL,IJ) * DP25 1666 RS11(KL,IJ) = RS11(KL,IJ) * DP25 1667C 1668 UT11(KL,IJ) = UT11(KL,IJ) * D1P5 * DP25 1669 BT11(KL,IJ) = BT11(KL,IJ) * D1P5 * DP25 1670 VT11(KL,IJ) = VT11(KL,IJ) * D1P5 1671 WT11(KL,IJ) = WT11(KL,IJ) * D1P5 1672 QT11(KL,IJ) = QT11(KL,IJ) * DP75 1673 RT11(KL,IJ) = RT11(KL,IJ) * DP75 1674c 1675 303 CONTINUE 1676 302 CONTINUE 1677 301 CONTINUE 1678 300 CONTINUE 1679C 1680 IF (R12OLD) THEN 1681 CALL DCOPY(NSPAIR*NSPAIR,VS11,1,WS11,1) 1682 CALL DCOPY(NSPAIR*NSPAIR,US11,1,BS11,1) 1683 CALL DCOPY(NSPAIR*NSPAIR,QS11,1,RS11,1) 1684 CALL DCOPY(NSPAIR*NSPAIR,VT11,1,WT11,1) 1685 CALL DCOPY(NSPAIR*NSPAIR,UT11,1,BT11,1) 1686 CALL DCOPY(NSPAIR*NSPAIR,QT11,1,RT11,1) 1687 ENDIF 1688C 1689 CALL VSHRNK(VS11,NSPAIR,NRHF,NRHFA,NSYM) 1690 CALL VSHRNK(VT11,NSPAIR,NRHF,NRHFA,NSYM) 1691 CALL VSHRNK(WS11,NSPAIR,NRHF,NRHFA,NSYM) 1692 CALL VSHRNK(WT11,NSPAIR,NRHF,NRHFA,NSYM) 1693C 1694 IF (IPRINT .LE. 3) GOTO 998 1695 CALL AROUND('OLD singlet <V> matrix') 1696 CALL OUTPUT(VS11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 1697 CALL AROUND('OLD singlet <B> matrix') 1698 CALL OUTPUT(US11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1699 CALL AROUND('OLD singlet <S> matrix') 1700 CALL OUTPUT(QS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1701 CALL AROUND('OLD triplet <V> matrix') 1702 CALL OUTPUT(VT11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 1703 CALL AROUND('OLD triplet <B> matrix') 1704 CALL OUTPUT(UT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1705 CALL AROUND('OLD triplet <S> matrix') 1706 CALL OUTPUT(QT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1707 IF (R12OLD) GOTO 998 1708 CALL AROUND('NEW singlet <V> matrix') 1709 CALL OUTPUT(WS11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 1710 CALL AROUND('NEW singlet <B> matrix') 1711 CALL OUTPUT(BS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1712 CALL AROUND('NEW singlet <S> matrix') 1713 CALL OUTPUT(RS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1714 CALL AROUND('NEW triplet <V> matrix') 1715 CALL OUTPUT(WT11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 1716 CALL AROUND('NEW triplet <B> matrix') 1717 CALL OUTPUT(BT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1718 CALL AROUND('NEW triplet <S> matrix') 1719 CALL OUTPUT(RT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1720 998 CONTINUE 1721 CALL GPOPEN (LUMULO,'AOR12O','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1722 CALL GPOPEN (LUMULN,'AOR12N','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1723 WRITE(LUMULO,'(4E30.20)') VS11 1724 WRITE(LUMULO,'(4E30.20)') US11 1725 WRITE(LUMULO,'(4E30.20)') QS11 1726 WRITE(LUMULO,'(4E30.20)') VT11 1727 WRITE(LUMULO,'(4E30.20)') UT11 1728 WRITE(LUMULO,'(4E30.20)') QT11 1729 WRITE(LUMULN,'(4E30.20)') WS11 1730 WRITE(LUMULN,'(4E30.20)') BS11 1731 WRITE(LUMULN,'(4E30.20)') RS11 1732 WRITE(LUMULN,'(4E30.20)') WT11 1733 WRITE(LUMULN,'(4E30.20)') BT11 1734 WRITE(LUMULN,'(4E30.20)') RT11 1735C 1736 LU43 = -43 1737 LU44 = -44 1738 LU45 = -45 1739 IF (R12OLD) THEN 1740 LUMULA = LUMULO 1741 LUMULE = LUMULO 1742 ELSE 1743 IF (R12XXL) THEN 1744 LUMULA = LUMULO 1745 ELSE 1746 LUMULA = LUMULN 1747 END IF 1748 LUMULE = LUMULN 1749 END IF 1750 DO 999 LUMULX = LUMULA, LUMULE 1751 IF (LUMULX .EQ. LUMULO) THEN 1752 CALL AROUND('Original MP2-R12 method') 1753 IF (.NOT. R12NOB .OR. NATVIR) THEN 1754 CALL GPOPEN(LU43,'QMATSO','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1755 READ(LU43,'(4E30.20)') US11 1756 CALL GPCLOSE(LU43,'KEEP') 1757 CALL GPOPEN(LU43,'QMATTO','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1758 READ(LU43,'(4E30.20)') UT11 1759 CALL GPCLOSE(LU43,'KEEP') 1760 END IF 1761 ELSE 1762 CALL AROUND('Auxiliary-basis MP2-R12 method') 1763 IF (.NOT. R12NOB .OR. NATVIR) THEN 1764 CALL GPOPEN(LU43,'QMATSN','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1765 READ(LU43,'(4E30.20)') US11 1766 CALL GPCLOSE(LU43,'KEEP') 1767 CALL GPOPEN(LU43,'QMATTN','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1768 READ(LU43,'(4E30.20)') UT11 1769 CALL GPCLOSE(LU43,'KEEP') 1770 IF (.NOT. NORXR) THEN 1771 CALL GPOPEN(LU43,'XMATSN','UNKNOWN',' ','FORMATTED', 1772 & IDUM,LDUM) 1773 READ(LU43,'(4E30.20)') RXRS 1774 CALL GPCLOSE(LU43,'KEEP') 1775 CALL GPOPEN(LU43,'XMATTN','UNKNOWN',' ','FORMATTED', 1776 & IDUM,LDUM) 1777 READ(LU43,'(4E30.20)') RXRT 1778 CALL GPCLOSE(LU43,'KEEP') 1779 END IF 1780 END IF 1781 END IF 1782 REWIND LUMULX 1783 READ(LUMULX,'(4E30.20)') WS11 1784 READ(LUMULX,'(4E30.20)') BS11 1785 READ(LUMULX,'(4E30.20)') RS11 1786 READ(LUMULX,'(4E30.20)') WT11 1787 READ(LUMULX,'(4E30.20)') BT11 1788 READ(LUMULX,'(4E30.20)') RT11 1789 LUM = LUMULX 1790 CALL GPCLOSE(LUM,'KEEP') 1791celena 1792 IF (R12PRP) THEN 1793 CALL GPOPEN(LU43,'CCR12_XP','UNKNOWN',' ','FORMATTED', 1794 & IDUM,LDUM) 1795 WRITE(LU43,'(I3)') 1 1796 WRITE(LU43,'(4E30.20)') RS11 1797 WRITE(LU43,'(4E30.20)') RT11 1798 CALL GPCLOSE(LU43,'KEEP') 1799 ENDIF 1800celena 1801 1802c Write out V-matrix, etc. for CC2-R12 model (HF/UniKA/02-05-2003). 1803c modified by C. Neiss: no not use singlet/triplet format any more, 1804c make intermediates compatible with correlation factor r_12 (and 1805c not 0.5*r_12 as in the MP2-R12 code)! 1806c (April 2005) 1807c 1808 IF (CCR12) THEN 1809 KSCR = 1 1810 KEND1 = KSCR + NGAMMA(1) !NGAMMA(1) in MP2-R12 >= NKILJ 1811c WRITE(LUPRI,*) 'NGAMMA(1) in CC_R12: ', NGAMMA(1) 1812 IF (NGAMMA(1) .gt. LWORK) THEN 1813 CALL QUIT('Insufficient WORK space in MAKEVR') 1814 END IF 1815 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,WS11,1) 1816 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,WT11,1) 1817 CALL DSCAL(NSPAIR*NSPAIR,1.0D0/3.0D0,WT11,1) 1818 CALL CCR12PCK(WORK(KSCR),1,WS11,WT11,NRHF,NRHFA,NKILJ) 1819c WRITE(LUPRI,*) 'VIJKL written in MP2R12:' 1820c CALL CC_PRPR12(WORK(KSCR),1,1,.FALSE.) 1821 CALL GPOPEN(LU43,FCCR12V,'UNKNOWN',' ','UNFORMATTED', 1822 & IDUM,LDUM) 1823 WRITE(LU43) 1 1824 WRITE(LU43) (WORK(KSCR-1+I), I=1, NKILJ(1)) 1825 CALL GPCLOSE(LU43,'KEEP') 1826 ! undo scaling 1827 CALL DSCAL(NSPAIR*NSPAIR,3.0D0,WT11,1) 1828 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,WS11,1) 1829 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,WT11,1) 1830c 1831c Write out X-matrix for CC2-R12 model (HF/UniKA/02-05-2003). 1832c 1833 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,RS11,1) 1834 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,RT11,1) 1835 CALL DSCAL(NSPAIR*NSPAIR,1.0D0/3.0D0,RT11,1) 1836 CALL CCR12PCK(WORK(KSCR),1,RS11,RT11,NRHF,NRHF,NKILJ) 1837 CALL GPOPEN(LU43,FCCR12X,'UNKNOWN',' ','UNFORMATTED', 1838 & IDUM,LDUM) 1839 WRITE(LU43) 1 1840 WRITE(LU43) (WORK(KSCR-1+I), I=1, NKILJ(1)) 1841 CALL GPCLOSE(LU43,'KEEP') 1842 ! undo scaling 1843 CALL DSCAL(NSPAIR*NSPAIR,3.0D0,RT11,1) 1844 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,RS11,1) 1845 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,RT11,1) 1846c 1847c Write out orbital energies for CC2-R12 model (HF/UniKA/02-05-2003). 1848c 1849 CALL GPOPEN(LU43,FCCR12E,'UNKNOWN',' ','FORMATTED', 1850 & IDUM,LDUM) 1851 WRITE(LU43,'(4E30.20)') EVL 1852 CALL GPCLOSE(LU43,'KEEP') 1853 END IF 1854c 1855c Open files for B and C (HF/UniKA/25-06-2003). 1856c 1857 IF (CCR12 .AND. (LUMULA .EQ. LUMULX)) THEN 1858 CALL GPOPEN(LU44,FCCR12B,'UNKNOWN',' ','UNFORMATTED', 1859 & IDUM,LDUM) 1860 CALL GPOPEN(LU45,FCCR12D,'UNKNOWN',' ','UNFORMATTED', 1861 & IDUM,LDUM) 1862 END IF 1863C 1864 DO 999 IPDD = 0, 7 1865chf 1866 IF (R12HYB .OR. NORXR) THEN 1867 IPCC = APPROX(IPDD,1) 1868 ELSE 1869 IPCC = APPROX(IPDD,2) 1870 END IF 1871chf 1872 IF (NATVIR .AND. IPDD .EQ. 3) THEN 1873 DO I = 1, NSPAIR 1874 DO J = 1, NSPAIR 1875 BS11(I,J) = BS11(I,J) - US11(I,J) 1876 BT11(I,J) = BT11(I,J) - UT11(I,J) 1877 END DO 1878 END DO 1879 END IF 1880 IF (IPDD .EQ. 4 .AND. .NOT. NATVIR .AND. .NOT. R12NOB) THEN 1881 DO I = 1, NSPAIR 1882 DO J = 1, NSPAIR 1883 BS11(I,J) = BS11(I,J) - US11(I,J) 1884 BT11(I,J) = BT11(I,J) - UT11(I,J) 1885 END DO 1886 END DO 1887 END IF 1888 IF (IPDD .EQ. 6) THEN 1889 IF (.NOT. R12NOB .AND..NOT. NORXR .AND. LUMULX .EQ. LUMULN) THEN 1890 DO I = 1, NSPAIR 1891 DO J = 1, NSPAIR 1892 BS11(I,J) = BS11(I,J) + RXRS(I,J) 1893 BT11(I,J) = BT11(I,J) + RXRT(I,J) 1894 END DO 1895 END DO 1896 END IF 1897 END IF 1898 IF (LUMULX .EQ. LUMULO .OR. NORXR) THEN 1899 IF (.NOT. R12XXL .AND. 1900 * (IPDD .NE. 2 .AND. IPDD .NE. 5)) GOTO 999 1901 ELSE 1902 IF (.NOT. R12XXL .AND. 1903 * (IPDD .NE. 2 .AND. IPDD .NE. 7)) GOTO 999 1904 END IF 1905 IF (IPDD .LE. 3 .AND. R12NOA) GOTO 999 1906 IF (IPDD .EQ. 3 .AND. R12NOP) GOTO 999 1907 IF (IPDD .GE. 4 .AND. R12NOB) GOTO 999 1908 IF (IPDD .GE. 6 .AND. (LUMULX .EQ. LUMULO .OR. NORXR)) GOTO 999 1909 IF (IPDD .EQ. 0) THEN 1910 XPDD = D0 1911 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER NONINVARIANT R12/0'// 1912 * ' PAIR ENERGIES WITH FIXED R12-COEFFICIENTS:' 1913 ELSE IF (IPDD .EQ. 1) THEN 1914 XPDD = D0 1915 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER NONINVARIANT R12/A'// 1916 * ' PAIR ENERGIES:' 1917 ELSE IF (IPDD .EQ. 2) THEN 1918 XPDD = D0 1919 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER R12/A PAIR ENERGIES:' 1920 ELSE IF (IPDD .EQ. 3) THEN 1921 XPDD = DP5 1922 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER R12/A'' PAIR ENERGIES:' 1923 ELSE IF (IPDD .EQ. 4) THEN 1924 XPDD = D0 1925 IF (R12HYB .OR. (LUMULX .EQ. LUMULO)) THEN 1926 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER NONINVARIANT R12/B'// 1927 * ' PAIR ENERGIES:' 1928 ELSE 1929 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER NONINVARIANT R12/B'''// 1930 * ' PAIR ENERGIES:' 1931 END IF 1932 ELSE IF (IPDD .EQ. 5) THEN 1933 XPDD = DP5 1934 IF (R12HYB .OR. (LUMULX .EQ. LUMULO)) THEN 1935 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER R12/B PAIR ENERGIES:' 1936 ELSE 1937 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER R12/B'' PAIR ENERGIES:' 1938 END IF 1939 ELSE IF (IPDD .EQ. 6) THEN 1940 XPDD = D0 1941 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER NONINVARIANT R12/B'// 1942 * ' PAIR ENERGIES:' 1943 ELSE IF (IPDD .EQ. 7) THEN 1944 XPDD = DP5 1945 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER R12/B PAIR ENERGIES:' 1946 END IF 1947 F2S = D0 1948 F2T = D0 1949 CSALL = -D1 1950 CTALL = -D1 1951c 1952c Write out B-matrix for CC2-R12 model (HF/UniKA/02-05-2003). 1953c 1954 IF (CCR12) THEN 1955 IF (LUMULX .EQ. LUMULO) THEN 1956 WRITE(LU44) 0,IPDD,IPCC 1957 ELSE 1958 WRITE(LU44) 1,IPDD,IPCC 1959 END IF 1960 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,BS11,1) 1961 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,BT11,1) 1962 CALL DSCAL(NSPAIR*NSPAIR,1.0D0/3.0D0,BT11,1) 1963 CALL CCR12PCK(WORK(KSCR),1,BS11,BT11,NRHF,NRHF,NKILJ) 1964 WRITE(LU44) (WORK(KSCR-1+I), I=1, NKILJ(1)) 1965 ! undo scaling 1966 CALL DSCAL(NSPAIR*NSPAIR,3.0D0,BT11,1) 1967 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,BS11,1) 1968 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,BT11,1) 1969 END IF 1970C 1971C CCMS 1972C 1973 IF (LGLOCAL .AND. (IPDD.EQ.3 .OR. IPDD.EQ.5 .OR. IPDD.EQ.7) ) THEN 1974 CALL R12FXL(OLWS,OLWT,FLMAT,RS11,RT11,NIND,III,JJJ,NOCDIM,NSPAIR) 1975 IF (IPRINT .GE. 5) THEN 1976 CALL AROUND('Singlet <W> matrix') 1977 CALL OUTPUT(OLWS,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1978 CALL AROUND('Triplet <W> matrix') 1979 CALL OUTPUT(OLWT,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1980 END IF 1981 IF (R12RST) THEN 1982 WRITE(CFN,1001) IPDD,LUMULX 1983 CALL GPOPEN(LU43,CFN,'UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1984 READ(LU43,'(4E30.20)') CS11 1985 CALL GPCLOSE(LU43,'KEEP') 1986 WRITE(CFN,1002) IPDD,LUMULX 1987 CALL GPOPEN(LU43,CFN,'UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1988 READ(LU43,'(4E30.20)') CT11 1989 CALL GPCLOSE(LU43,'KEEP') 1990 ITER = 1 1991 ELSE 1992 ITER = 0 1993 END IF 1994 800 CONTINUE 1995 IF (ITER .GT. 0) THEN 1996 CALL R12FCL(OLZS,OLZT,OLWS,OLWT,WS11,WT11,OLYS,OLYT,RS11,RT11, 1997 * CS11,CT11,FLMAT,NIND,III,JJJ,NOCDIM,NSPAIR,DELTA) 1998 IF (IPRINT .GE. 5) THEN 1999 CALL AROUND('Singlet <Z> matrix') 2000 CALL OUTPUT(OLZS,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2001 CALL AROUND('Triplet <Z> matrix') 2002 CALL OUTPUT(OLZT,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2003 END IF 2004 END IF 2005 F2SI = F2S 2006 F2TI = F2T 2007 F2S = D0 2008 F2T = D0 2009 IJ = 0 2010 DO 500 I = 1, NACDIM 2011 DO 501 J = 1, I 2012 IJ = IJ + 1 2013 IF (ITER .EQ. 0) THEN 2014 IF (R12SVD) THEN 2015 CALL R12MP2(WS11(1,IJ),WT11(1,IJ),WS11(1,IJ),WT11(1,IJ), 2016 * RS11,RT11,EVL(IJ),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 2017 * QS11,QT11,BS11,BT11,EVL,XPDD,QQ11,IJ,WSMIN,WTMIN,D0, 2018 * CS11(1,IJ),CT11(1,IJ)) 2019 ELSE IF (R12DIA) THEN 2020 CALL MP2R12(WS11(1,IJ),WT11(1,IJ),WS11(1,IJ),WT11(1,IJ), 2021 * RS11,RT11,EVL(IJ),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 2022 * QS11,QT11,BS11,BT11,EVL,XPDD,QQ11,IJ,WSMIN,WTMIN,D0, 2023 * CS11(1,IJ),CT11(1,IJ)) 2024 ELSE 2025 CALL QUIT('No solver selected (R12SVD or R12DIA)') 2026 END IF 2027 ELSE 2028 IF (R12SVD) THEN 2029 CALL R12MP2(OLZS(1,IJ),OLZT(1,IJ),WS11(1,IJ),WT11(1,IJ), 2030 * RS11,RT11,EVL(IJ),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 2031 * QS11,QT11,BS11,BT11,EVL,XPDD,QQ11,IJ,WSMIN,WTMIN,DELTA, 2032 * CS11(1,IJ),CT11(1,IJ)) 2033 ELSE IF (R12DIA) THEN 2034 CALL MP2R12(OLZS(1,IJ),OLZT(1,IJ),WS11(1,IJ),WT11(1,IJ), 2035 * RS11,RT11,EVL(IJ),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 2036 * QS11,QT11,BS11,BT11,EVL,XPDD,QQ11,IJ,WSMIN,WTMIN,DELTA, 2037 * CS11(1,IJ),CT11(1,IJ)) 2038 ELSE 2039 CALL QUIT('No solver selected (R12SVD or R12DIA)') 2040 END IF 2041 END IF 2042 F2S = F2S + FS(IJ) 2043 F2T = F2T + FT(IJ) 2044 501 CONTINUE 2045 500 CONTINUE 2046 IF (IPRINT .GE. 4) THEN 2047 CALL AROUND('Singlet <C> matrix') 2048 CALL OUTPUT(CS11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 2049 CALL AROUND('Triplet <C> matrix') 2050 CALL OUTPUT(CT11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 2051 END IF 2052 CSALLI = CSALL 2053 CTALLI = CTALL 2054 CSALL = DDOT(NSPAIR*NAPAIR,CS11,1,CS11,1) 2055 CTALL = DDOT(NSPAIR*NAPAIR,CT11,1,CT11,1) 2056 CSALL = SQRT(CSALL) 2057 CTALL = SQRT(CTALL) 2058C 2059 IF (ITER.LT.100) THEN 2060C IF (ABS(F2S-F2SI).GT.1.0d-10.OR.ABS(F2T-F2TI).GT.1.0d-10) THEN 2061 IF (ABS(CSALLI-CSALL).GT.1.0d-8. 2062 * OR.ABS(CTALLI-CTALL).GT.1.0d-8) THEN 2063 WRITE(LUPRI,'(I3,A,6F15.9)') ITER,'@@', 2064 * E2S,E2S+F2S,E2T,E2T+F2T, 2065 * CSALL, CTALL 2066 ITER = ITER + 1 2067 GO TO 800 2068 ELSE 2069 WRITE(LUPRI,'(I3,A,6F15.9/)') ITER,'@@', 2070 * E2S,E2S+F2S,E2T,E2T+F2T, 2071 * CSALL, CTALL 2072 DO IJ = 1, NAPAIR 2073 WRITE(LUPRI,'(I4,4F15.9,2G16.6)') IJ,ES(IJ),ES(IJ)+FS(IJ), 2074 * ET(IJ),ET(IJ)+FT(IJ), 2075 * WSMIN,WTMIN 2076 END DO 2077 WRITE(LUPRI,'(/A,I4,A)')' CONVERGED IN:',ITER,' ITERATIONS' 2078 END IF 2079 ELSE 2080 WRITE(LUPRI,'(A)') ' ITERATIONS FOR THE ENERGY EXCEEDED 100' 2081 END IF 2082 IF (IPRINT .GE. 4) THEN 2083 CALL AROUND('Singlet <C> matrix') 2084 CALL OUTPUT(CS11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 2085 CALL AROUND('Triplet <C> matrix') 2086 CALL OUTPUT(CT11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 2087 END IF 2088 WRITE(CFN,1001) IPDD,LUMULX 2089 CALL GPOPEN(LU43,CFN,'UNKNOWN',' ','FORMATTED',IDUM,LDUM) 2090 WRITE(LU43,'(4E30.20)') CS11 2091 CALL GPCLOSE(LU43,'KEEP') 2092 WRITE(CFN,1002) IPDD,LUMULX 2093 CALL GPOPEN(LU43,CFN,'UNKNOWN',' ','FORMATTED',IDUM,LDUM) 2094 WRITE(LU43,'(4E30.20)') CT11 2095 CALL GPCLOSE(LU43,'KEEP') 2096C 2097 ELSE 2098C 2099 CALL DZERO(FS,NSPAIR) 2100 CALL DZERO(FT,NSPAIR) 2101 CALL DZERO(CS11,NSPAIR*NSPAIR) 2102 CALL DZERO(CT11,NSPAIR*NSPAIR) 2103 F2S = D0 2104 F2T = D0 2105 IJ = 0 2106 JI = 0 2107 II = 0 2108 DO 511 ISYM = 1, NSYM 2109 DO 511 I = 1, NRHF(ISYM) 2110 II = II + 1 2111 JJ = 0 2112 DO 512 JSYM = 1, NSYM 2113 DO 512 J = 1, NRHF(JSYM) 2114 JJ = JJ + 1 2115 IF (JJ .GT. II) GOTO 512 2116 JI = JI + 1 2117 IF (I .GT. NRHFA(ISYM) .OR. 2118 * J .GT. NRHFA(JSYM)) GOTO 512 2119 IJ = IJ + 1 2120 CNINV(IJ,1) = WS11(JI,IJ) 2121 CNINV(IJ,2) = BS11(JI,JI) 2122 CNINV(IJ,9) = -BS11(JI,JI) / RS11(JI,JI) 2123 IF (IPDD .EQ. 0) THEN 2124 CNINV(IJ,3) = 1D0 2125 CNINV(IJ,4) = 2D0*WS11(JI,IJ) - BS11(JI,JI) 2126 ELSE 2127 IF (CNINV(IJ,2) .LT. -VCLTHR) THEN 2128 CNINV(IJ,3) = CNINV(IJ,1) / CNINV(IJ,2) 2129 ELSE 2130 WRITE(LUPRI,'(/A,E15.8,I4/)') 2131 & ' @ WARNING: NEGATIVE DIAGONAL ELEMENT OF SINGLET'// 2132 & ' B MATRIX',CNINV(IJ,2),IJ 2133 CNINV(IJ,3) = D0 2134 END IF 2135 CNINV(IJ,4) = CNINV(IJ,1) * CNINV(IJ,3) 2136 END IF 2137 IF (II .NE. JJ) THEN 2138 CNINV(IJ,5) = WT11(JI,IJ) 2139 CNINV(IJ,6) = BT11(JI,JI) 2140 CNINV(IJ,10) = -BT11(JI,JI) / RT11(JI,JI) 2141 IF (IPDD .EQ. 0) THEN 2142 CNINV(IJ,7) = 0.5D0 2143 CNINV(IJ,8) = WT11(JI,IJ) - 0.25D0*BT11(JI,JI) 2144 ELSE 2145 IF (CNINV(IJ,6) .LT. -VCLTHR) THEN 2146 CNINV(IJ,7) = CNINV(IJ,5) / CNINV(IJ,6) 2147 ELSE 2148 WRITE(LUPRI,'(/A,E15.8,I4/)') 2149 & ' @ WARNING: NEGATIVE DIAGONAL ELEMENT OF TRIPLET'// 2150 & ' B MATRIX',CNINV(IJ,6),IJ 2151 CNINV(IJ,7) = D0 2152 END IF 2153 CNINV(IJ,8) = CNINV(IJ,5) * CNINV(IJ,7) 2154 END IF 2155 ELSE 2156 CNINV(IJ,5) = D0 2157 CNINV(IJ,6) = D0 2158 CNINV(IJ,7) = D0 2159 CNINV(IJ,8) = D0 2160 CNINV(IJ,10) = D0 2161 END IF 2162 FS(IJ) = CNINV(IJ,4) 2163 FT(IJ) = CNINV(IJ,8) 2164 IF (IPDD .GT. 1 .AND. IPDD .NE. 4 .AND. IPDD .NE. 6) THEN 2165 IF (R12SVD) THEN 2166 CALL R12MP2(WS11(1,IJ),WT11(1,IJ),WS11(1,IJ),WT11(1,IJ), 2167 * RS11,RT11,EVS(JI),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 2168 * QS11,QT11,BS11,BT11,EVS,XPDD,QQ11,IJ,WSMIN,WTMIN,D0, 2169 * CS11(1,IJ),CT11(1,IJ)) 2170 ELSE IF (R12DIA) THEN 2171 CALL MP2R12(WS11(1,IJ),WT11(1,IJ),WS11(1,IJ),WT11(1,IJ), 2172 * RS11,RT11,EVS(JI),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 2173 * QS11,QT11,BS11,BT11,EVS,XPDD,QQ11,IJ,WSMIN,WTMIN,D0, 2174 * CS11(1,IJ),CT11(1,IJ)) 2175 ELSE 2176 CALL QUIT('No solver selected (R12SVD or R12DIA)') 2177 END IF 2178 WRITE(LUPRI,'(I4,4F15.9,2G16.6)') IJ,ES(IJ),ES(IJ)+FS(IJ), 2179 * ET(IJ),ET(IJ)+FT(IJ), 2180 * WSMIN,WTMIN 2181 ELSE 2182 WRITE(LUPRI,'(I4,4F15.9,2G16.6)') IJ,ES(IJ),ES(IJ)+FS(IJ), 2183 * ET(IJ),ET(IJ)+FT(IJ) 2184 CS11(IJ,IJ) = CNINV(IJ,3) 2185 CT11(IJ,IJ) = CNINV(IJ,7) 2186 END IF 2187 F2S = F2S + FS(IJ) 2188 F2T = F2T + FT(IJ) 2189 512 CONTINUE 2190 511 CONTINUE 2191C 2192 END IF 2193C 2194 WRITE(LUPRI,'(/4X,6F15.9)') E2S,E2S+F2S,E2T,E2T+F2T 2195c 2196c Write out amplitudes for CC2-R12 model (HF/UniKA/25-06-2003). 2197c 2198 IF (CCR12) THEN 2199 IF (LUMULX .EQ. LUMULO) THEN 2200 WRITE(LU45) 0,IPDD,IPCC 2201 ELSE 2202 WRITE(LU45) 1,IPDD,IPCC 2203 END IF 2204 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,CS11,1) 2205 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,CT11,1) 2206 CALL CCR12PCK(WORK(KSCR),1,CS11,CT11,NRHF,NRHFA,NKILJ) 2207 WRITE(LU45) (WORK(KSCR-1+I), I=1, NKILJ(1)) 2208 ! undo scaling: 2209 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,CS11,1) 2210 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,CT11,1) 2211 END IF 2212c 2213 WRITE(CFN,1001) IPDD,LUMULX 2214 CALL GPOPEN(LU43,CFN,'UNKNOWN',' ','FORMATTED',IDUM,LDUM) 2215 WRITE(LU43,'(4E30.20)') CS11 2216 CALL GPCLOSE(LU43,'KEEP') 2217 WRITE(CFN,1002) IPDD,LUMULX 2218 CALL GPOPEN(LU43,CFN,'UNKNOWN',' ','FORMATTED',IDUM,LDUM) 2219 WRITE(LU43,'(4E30.20)') CT11 2220 CALL GPCLOSE(LU43,'KEEP') 2221c 2222 IF (IPDD .EQ. 0) THEN 2223 WRITE(LUPRI,'(/A,F15.9//A/)') 2224 * ' Noninvariant MP2-R12/0 correlation energy =', 2225 * E2S+E2T+F2S+F2T, 2226 * ' IJ V(IJ) U(IJ) C(IJ) '// 2227 * ' V(IJ) U(IJ) C(IJ) ' 2228 IJ = 0 2229 DO I = 1, NACDIM 2230 DO J = 1, I 2231 IJ = IJ + 1 2232 WRITE(LUPRI,'(I4,8E12.4)') IJ,(CNINV(IJ,K),K=1,3), 2233 * (CNINV(IJ,K),K=5,7),CNINV(IJ,9), 2234 * CNINV(IJ,10) 2235 END DO 2236 END DO 2237 ELSE IF (IPDD .EQ. 1) THEN 2238 WRITE(LUPRI,'(/A,F15.9//A/)') 2239 * ' Noninvariant MP2-R12/A correlation energy =', 2240 * E2S+E2T+F2S+F2T, 2241 * ' IJ V(IJ) U(IJ) C(IJ) '// 2242 * ' V(IJ) U(IJ) C(IJ) ' 2243 IJ = 0 2244 DO 600 I = 1, NACDIM 2245 DO 601 J = 1, I 2246 IJ = IJ + 1 2247 WRITE(LUPRI,'(I4,8E12.4)') IJ,(CNINV(IJ,K),K=1,3), 2248 * (CNINV(IJ,K),K=5,7),CNINV(IJ,9), 2249 * CNINV(IJ,10) 2250 601 CONTINUE 2251 600 CONTINUE 2252 ELSE IF (IPDD .EQ. 2) THEN 2253 IF (IPRINT .GE. 4) THEN 2254c 2255c Print MP2-R12/A -Ansatz 1- amplitudes for 2256c 2257 CALL AROUND('Singlet <C> matrix (Ansatz 1, MP2-R12/A)') 2258 CALL OUTPUT(CS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2259 CALL AROUND('Triplet <C> matrix (Ansatz 1, MP2-R12/A)') 2260 CALL OUTPUT(CT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2261 ENDIF 2262c 2263c Write out amplitudes for MP2-R12-Gradient(Elena) 2264c 2265 LUNIT = -1 2266 IF (R12PRP .AND. LUMULX .NE. LUMULO) THEN 2267 LUNIT = -1 2268 CALL GPOPEN(LUNIT,'CCR12_C','UNKNOWN',' ','FORMATTED', 2269 & IDUM,LDUM) 2270 WRITE(LUNIT,'(4E30.20)') CS11 2271 WRITE(LUNIT,'(4E30.20)') CT11 2272 CALL GPCLOSE(LUNIT,'KEEP') 2273C 2274 CALL CC_R12GRAD(IPDD,EV,V2AM,GRAD,LWORK) 2275 IF (FROIMP) THEN 2276 CALL CC_R12GRADFR(IPDD,EV,V2AM,GRAD,LWORK) 2277 END IF 2278 END IF 2279c 2280cElena 2281 2282 WRITE(LUPRI,'(/A,F15.9 )') 2283 * ' MP2-R12/A correlation energy =', 2284 * E2S+E2T+F2S+F2T 2285 ELSE IF (IPDD .EQ. 3) THEN 2286 2287 LUNIT = -1 2288 IF (R12PRP .AND. LUMULX .NE. LUMULO) THEN 2289 IF (IPRINT .GE. 4) THEN 2290 CALL AROUND('Singlet <C> matrix (Ansatz 1, MP2-R12/A'')') 2291 CALL OUTPUT(CS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2292 CALL AROUND('Triplet <C> matrix (Ansatz 1, MP2-R12/A'')') 2293 CALL OUTPUT(CT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2294 END IF 2295 LUNIT = -1 2296 CALL GPOPEN(LUNIT,'CCR12_C','UNKNOWN',' ','FORMATTED', 2297 & IDUM,LDUM) 2298 WRITE(LUNIT,'(4E30.20)') CS11 2299 WRITE(LUNIT,'(4E30.20)') CT11 2300 CALL GPCLOSE(LUNIT,'KEEP') 2301c KEND1 = 1 2302c LWRK1 = LWORK - KEND1 2303 2304 CALL CC_R12GRAD(IPDD,EV,V2AM,GRAD,LWORK) 2305 IF (FROIMP) THEN 2306 CALL CC_R12GRADFR(IPDD,EV,V2AM,GRAD,LWORK) 2307 END IF 2308 END IF 2309c 2310 WRITE(LUPRI,'(/A,F15.9 )') 2311 * ' MP2-R12/A'' correlation energy =', 2312 * E2S+E2T+F2S+F2T 2313 ELSE IF (IPDD .EQ. 4) THEN 2314 IF (R12HYB .OR. (LUMULX .EQ. LUMULO)) THEN 2315 WRITE(LUPRI,'(/A,F15.9//A/)') 2316 * ' Noninvariant MP2-R12/B correlation energy =', 2317 * E2S+E2T+F2S+F2T, 2318 * ' IJ V(IJ) U(IJ) C(IJ) '// 2319 * ' V(IJ) U(IJ) C(IJ) ' 2320 ELSE 2321 WRITE(LUPRI,'(/A,F15.9//A/)') 2322 * ' Noninvariant MP2-R12/B'' correlation energy =', 2323 * E2S+E2T+F2S+F2T, 2324 * ' IJ V(IJ) U(IJ) C(IJ) '// 2325 * ' V(IJ) U(IJ) C(IJ) ' 2326 END IF 2327 IJ = 0 2328 DO 602 I = 1, NACDIM 2329 DO 603 J = 1, I 2330 IJ = IJ + 1 2331 WRITE(LUPRI,'(I4,6E12.4)') IJ,(CNINV(IJ,K),K=1,3), 2332 * (CNINV(IJ,K),K=5,7) 2333 603 CONTINUE 2334 602 CONTINUE 2335 ELSE IF (IPDD .EQ. 5) THEN 2336 IF (R12HYB .OR. (LUMULX .EQ. LUMULO)) THEN 2337 IF (IPRINT. GE. 4) THEN 2338 CALL AROUND('Singlet <C> matrix (Ansatz 1, MP2-R12/B)') 2339 CALL OUTPUT(CS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2340 CALL AROUND('Triplet <C> matrix (Ansatz 1, MP2-R12/B)') 2341 CALL OUTPUT(CT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 2342 END IF 2343 LUNIT = -1 2344 IF (R12PRP .AND. LUMULX .NE. LUMULO) THEN 2345 LUNIT = -1 2346 CALL GPOPEN(LUNIT,'CCR12_C','UNKNOWN',' ','FORMATTED', 2347 & IDUM,LDUM) 2348 WRITE(LUNIT,'(4E30.20)') CS11 2349 WRITE(LUNIT,'(4E30.20)') CT11 2350 CALL GPCLOSE(LUNIT,'KEEP') 2351 CALL CC_R12GRAD(IPDD,EV,V2AM,GRAD,LWORK) 2352 IF (FROIMP) THEN 2353 CALL CC_R12GRADFR(IPDD,EV,V2AM,GRAD,LWORK) 2354 END IF 2355 ENDIF 2356 WRITE(LUPRI,'(/A,F15.9 )') 2357 * ' MP2-R12/B correlation energy =', 2358 * E2S+E2T+F2S+F2T 2359 ELSE 2360 WRITE(LUPRI,'(/A,F15.9 )') 2361 * ' MP2-R12/B'' correlation energy =', 2362 * E2S+E2T+F2S+F2T 2363 END IF 2364 ELSE IF (IPDD .EQ. 6) THEN 2365 WRITE(LUPRI,'(/A,F15.9//A/)') 2366 * ' Noninvariant MP2-R12/B correlation energy =', 2367 * E2S+E2T+F2S+F2T, 2368 * ' IJ V(IJ) U(IJ) C(IJ) '// 2369 * ' V(IJ) U(IJ) C(IJ) ' 2370 IJ = 0 2371 DO 604 I = 1, NACDIM 2372 DO 605 J = 1, I 2373 IJ = IJ + 1 2374 WRITE(LUPRI,'(I4,6E12.4)') IJ,(CNINV(IJ,K),K=1,3), 2375 * (CNINV(IJ,K),K=5,7) 2376 605 CONTINUE 2377 604 CONTINUE 2378 ELSE IF (IPDD .EQ. 7) THEN 2379 WRITE(LUPRI,'(/A,F15.9 )') 2380 * ' MP2-R12/B correlation energy =', 2381 * E2S+E2T+F2S+F2T 2382 END IF 2383 999 CONTINUE 2384 IF (CCR12) THEN 2385 CALL GPCLOSE(LU44,'KEEP') 2386 CALL GPCLOSE(LU45,'KEEP') 2387 ENDIF 2388 RETURN 2389 END 2390C /* Deck mp2r12 */ 2391 SUBROUTINE MP2R12(VS,VT,VOS,VOT,SS,ST,EKL,NOCDIM,NSPAIR,FS,FT, 2392 * BB,UU,PP,QQ,SF,TF,EV,XF,RR,IJPAIR, 2393 * WS,WT,DELTA,CCS,CCT) 2394C 2395C This subroutine evaluates the MP2-R12 energy by 2396C solving sets of linear equations. 2397C 2398C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 2399C 2400#include "implicit.h" 2401#include "priunit.h" 2402 PARAMETER (QNIL = 1D-12, QTHR = 1D-10, D0 = 0D0, 2403 * D1 = 1D0, D2 = 2D0, DP5 = 0.5D0, D99 = 1D99) 2404 DIMENSION VS(NSPAIR), VT(NSPAIR), EV(NSPAIR) 2405 DIMENSION SS(NSPAIR,NSPAIR), ST(NSPAIR,NSPAIR) 2406 DIMENSION SF(NSPAIR,NSPAIR), TF(NSPAIR,NSPAIR) 2407 DIMENSION BB(*), PP(*), QQ(NSPAIR,NSPAIR), UU(NSPAIR,NSPAIR) 2408 DIMENSION RR(NSPAIR,NSPAIR) 2409 DIMENSION CCS(NSPAIR),CCT(NSPAIR) 2410 DIMENSION VOS(NSPAIR,NSPAIR),VOT(NSPAIR,NSPAIR) 2411C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2412C 2413 2414 N12 = NSPAIR * (NSPAIR + 1) / 2 2415 N2 = NSPAIR * NSPAIR 2416 WS = D99 2417 WT = D99 2418C 2419C SINGLET PAIRS 2420C 2421 CALL DZERO(UU,N2) 2422 IJ = 0 2423 DO 100 I = 1, NSPAIR 2424 UU(I,I) = D1 2425 DO 101 J = 1, I 2426 IJ = IJ + 1 2427 BB(IJ) = DP5 * (SS(I,J) + SS(J,I)) 2428 101 CONTINUE 2429 100 CONTINUE 2430 CALL JACO(BB,UU,NSPAIR,NSPAIR,NSPAIR,PP,QQ) 2431 CALL DZERO(RR,N2) 2432 II = 0 2433 DO 102 I = 1, NSPAIR 2434 II = II + I 2435 BBI = BB(II) 2436 IF (BBI .GT. QTHR) THEN 2437 BBI = D1 / SQRT(BBI) 2438 ELSE IF (ABS(BBI) .LT. QNIL) THEN 2439 BBI = D0 2440 CALL DZERO(UU(1,I),NSPAIR) 2441 ELSE 2442 IF (IJPAIR .EQ. 1) 2443 & WRITE(LUPRI,'(/A,E15.8/)') 2444 & ' @ WARNING: SMALL EIGENVALUE OF SINGLET'// 2445 & ' R12 OVERLAP',BBI 2446 BBI = D0 2447 CALL DZERO(UU(1,I),NSPAIR) 2448 END IF 2449 DO 103 L = 1, NSPAIR 2450 BU = BBI * UU(L,I) 2451 DO 104 K = 1, NSPAIR 2452 RR(K,L) = RR(K,L) + UU(K,I) * BU 2453 104 CONTINUE 2454 103 CONTINUE 2455 102 CONTINUE 2456 IJ = 0 2457 DO 105 I = 1, NSPAIR 2458 DO 106 J = 1, I 2459 IJ = IJ + 1 2460 PP(IJ) = - SF(I,J) - SF(J,I) 2461 * - XF * (D2 * EKL - EV(I) - EV(J)) 2462 * * (SS(I,J) + SS(J,I)) 2463 PP(IJ) = PP(IJ) * DP5 2464 106 CONTINUE 2465 PP(IJ) = PP(IJ) + DELTA 2466 105 CONTINUE 2467 CALL UTHU(PP,BB,RR,UU,NSPAIR,NSPAIR) 2468 CALL DZERO(UU,N2) 2469 DO 107 I = 1, NSPAIR 2470 UU(I,I) = D1 2471 107 CONTINUE 2472 CALL JACO(BB,UU,NSPAIR,NSPAIR,NSPAIR,PP,QQ) 2473 CALL DZERO(QQ,N2) 2474 II = 0 2475 DO 108 I = 1, NSPAIR 2476 II = II + I 2477 BBI = BB(II) 2478 IF (BBI .GT. QTHR) THEN 2479 WS = MIN(WS,BBI) 2480 BBI = D1 / BBI 2481 ELSE IF (ABS(BBI) .LT. QNIL) THEN 2482 BBI = D0 2483 CALL DZERO(UU(1,I),NSPAIR) 2484 ELSE 2485 WRITE(LUPRI,'(/A,E15.8,I4/)') 2486 & ' @ WARNING: NEGATIVE EIGENVALUE OF SINGLET'// 2487 & ' B MATRIX',BBI,IJPAIR 2488 BBI = D0 2489 CALL DZERO(UU(1,I),NSPAIR) 2490 END IF 2491 DO 109 L = 1, NSPAIR 2492 BU = BBI * UU(L,I) 2493 DO 110 K = 1, NSPAIR 2494 QQ(K,L) = QQ(K,L) + UU(K,I) * BU 2495 110 CONTINUE 2496 109 CONTINUE 2497 108 CONTINUE 2498 CALL MPAB(RR,NSPAIR,NSPAIR,NSPAIR,NSPAIR, 2499 * VS,NSPAIR,1,NSPAIR,1, 2500 * PP,NSPAIR,1) 2501 CALL MPAB(RR,NSPAIR,NSPAIR,NSPAIR,NSPAIR, 2502 * VOS,NSPAIR,1,NSPAIR,1, 2503 * BB,NSPAIR,1) 2504 FS = D0 2505 DO 111 I = 1, NSPAIR 2506 UU(I,1) = D0 2507 DO 112 J = 1, NSPAIR 2508 UU(I,1) = UU(I,1) - QQ(I,J) * PP(J) 2509 112 CONTINUE 2510 FS = FS + BB(I) * UU(I,1) 2511 111 CONTINUE 2512 CALL MPAB(RR,NSPAIR,NSPAIR,NSPAIR,NSPAIR, 2513 * UU,NSPAIR,1,NSPAIR,1, 2514 * CCS,NSPAIR,1) 2515C 2516C TRIPLET PAIRS 2517C 2518 CALL DZERO(UU,N2) 2519 IJ = 0 2520 DO 300 I = 1, NSPAIR 2521 UU(I,I) = D1 2522 DO 301 J = 1, I 2523 IJ = IJ + 1 2524 BB(IJ) = DP5 * (ST(I,J) + ST(J,I)) 2525 301 CONTINUE 2526 300 CONTINUE 2527 CALL JACO(BB,UU,NSPAIR,NSPAIR,NSPAIR,PP,QQ) 2528 CALL DZERO(RR,N2) 2529 II = 0 2530 DO 302 I = 1, NSPAIR 2531 II = II + I 2532 BBI = BB(II) 2533 IF (BBI .GT. QTHR) THEN 2534 BBI = D1 / SQRT(BBI) 2535 ELSE IF (ABS(BBI) .LT. QNIL) THEN 2536 BBI = D0 2537 CALL DZERO(UU(1,I),NSPAIR) 2538 ELSE 2539 IF (IJPAIR .EQ. 1) 2540 & WRITE(LUPRI,'(/A,E15.8/)') 2541 & ' @ WARNING: SMALL EIGENVALUE OF TRIPLET'// 2542 & ' R12 OVERLAP',BBI 2543 BBI = D0 2544 CALL DZERO(UU(1,I),NSPAIR) 2545 END IF 2546 DO 303 L = 1, NSPAIR 2547 BU = BBI * UU(L,I) 2548 DO 304 K = 1, NSPAIR 2549 RR(K,L) = RR(K,L) + UU(K,I) * BU 2550 304 CONTINUE 2551 303 CONTINUE 2552 302 CONTINUE 2553 IJ = 0 2554 DO 305 I = 1, NSPAIR 2555 DO 306 J = 1, I 2556 IJ = IJ + 1 2557 PP(IJ) = - TF(I,J) - TF(J,I) 2558 * - XF * (D2 * EKL - EV(I) - EV(J)) 2559 * * (ST(I,J) + ST(J,I)) 2560 PP(IJ) = PP(IJ) * DP5 2561 306 CONTINUE 2562 IF (ABS(PP(IJ)) .GT. 1D-12) PP(IJ) = PP(IJ) + DELTA 2563 305 CONTINUE 2564 CALL UTHU(PP,BB,RR,UU,NSPAIR,NSPAIR) 2565 CALL DZERO(UU,N2) 2566 DO 307 I = 1, NSPAIR 2567 UU(I,I) = D1 2568 307 CONTINUE 2569 CALL JACO(BB,UU,NSPAIR,NSPAIR,NSPAIR,PP,QQ) 2570 CALL DZERO(QQ,N2) 2571 II = 0 2572 DO 308 I = 1, NSPAIR 2573 II = II + I 2574 BBI = BB(II) 2575 IF (BBI .GT. QTHR) THEN 2576 WT = MIN(WT,BBI) 2577 BBI = D1 / BBI 2578 ELSE IF (ABS(BBI) .LT. QNIL) THEN 2579 BBI = D0 2580 CALL DZERO(UU(1,I),NSPAIR) 2581 ELSE 2582 WRITE(LUPRI,'(/A,E15.8,I4/)') 2583 & ' @ WARNING: NEGATIVE EIGENVALUE OF TRIPLET'// 2584 & ' B MATRIX',BBI,IJPAIR 2585 BBI = D0 2586 CALL DZERO(UU(1,I),NSPAIR) 2587 END IF 2588 DO 309 L = 1, NSPAIR 2589 BU = BBI * UU(L,I) 2590 DO 310 K = 1, NSPAIR 2591 QQ(K,L) = QQ(K,L) + UU(K,I) * BU 2592 310 CONTINUE 2593 309 CONTINUE 2594 308 CONTINUE 2595 CALL MPAB(RR,NSPAIR,NSPAIR,NSPAIR,NSPAIR, 2596 * VT,NSPAIR,1,NSPAIR,1, 2597 * PP,NSPAIR,1) 2598 CALL MPAB(RR,NSPAIR,NSPAIR,NSPAIR,NSPAIR, 2599 * VOT,NSPAIR,1,NSPAIR,1, 2600 * BB,NSPAIR,1) 2601 FT = D0 2602 DO 311 I = 1, NSPAIR 2603 UU(I,1) = D0 2604 DO 312 J = 1, NSPAIR 2605 UU(I,1) = UU(I,1) - QQ(I,J) * PP(J) 2606 312 CONTINUE 2607 FT = FT + BB(I) * UU(I,1) 2608 311 CONTINUE 2609 CALL MPAB(RR,NSPAIR,NSPAIR,NSPAIR,NSPAIR, 2610 * UU,NSPAIR,1,NSPAIR,1, 2611 * CCT,NSPAIR,1) 2612 RETURN 2613 END 2614C /* Deck iq */ 2615 INTEGER FUNCTION IQ(I) 2616#include "implicit.h" 2617#include "priunit.h" 2618#include "ccorb.h" 2619#include "ccsdsym.h" 2620#include "r12int.h" 2621 J = 0 2622 K = 0 2623 DO 100 ISYM = 1, NSYM 2624 K = K + NRHFFR(ISYM) 2625 DO 101 L = 1, NRHF(ISYM) 2626 J = J + 1 2627 K = K + 1 2628 IF (J .EQ. I) THEN 2629 IQ = K 2630 GOTO 99 2631 END IF 2632 101 CONTINUE 2633 100 CONTINUE 2634 99 RETURN 2635 END 2636C /* Deck testuu */ 2637 SUBROUTINE TESTUU(U2AM,R2AM) 2638C 2639C This subroutine generates on R2AM an array (ip|[T1+T2,r12)|jq) 2640C with (ia).le.(jb) from the asymmetric integrals (ip|[T1,r12)|jq) 2641C on U2AM without this triangular constraint. 2642C 2643C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 2644C 2645#include "implicit.h" 2646#include "priunit.h" 2647 DIMENSION U2AM(*), R2AM(*) 2648#include "ccorb.h" 2649#include "ccsdinp.h" 2650#include "ccsdsym.h" 2651#include "r12int.h" 2652 DO 101 ISYMIA = 1, NSYM 2653 ISYMJB = ISYMIA 2654 NAIBJ = IT2AM(ISYMIA,ISYMJB) 2655 KOFFU = IU2AM(ISYMIA,ISYMJB) 2656 NTOTAI = NT1AM(ISYMIA) 2657 DO 301 IA = 1, NTOTAI 2658 DO 302 JB = 1, IA 2659 NAIBJ = NAIBJ + 1 2660 R2AM(NAIBJ) = U2AM(KOFFU + (IA-1)*NTOTAI + JB) + 2661 & U2AM(KOFFU + (JB-1)*NTOTAI + IA) 2662 302 CONTINUE 2663 301 CONTINUE 2664 101 CONTINUE 2665 RETURN 2666 END 2667C /* Deck makeug */ 2668 SUBROUTINE MAKEUG(GX,UG,QQ,QV,WORK,LWORK,IPBAS) 2669C 2670C This subroutine computes the orbitals i' by means 2671C of the matrix product |i'> = Sum_p' X(p'i) |p'>. 2672C 2673C It then calls R12QQ to compute the integrals 2674C (i'k|r12**2|jl), which factorize into products 2675C of one-particle matrices. 2676C 2677C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 2678C 2679#include "implicit.h" 2680#include "priunit.h" 2681 DIMENSION GX(*), UG(*), WORK(*), QQ(*),QV(*) 2682 CHARACTER*7 GUNAM 2683 INTEGER IDUM 2684 LOGICAL LDUM 2685#include "ccorb.h" 2686#include "ccsdsym.h" 2687#include "ccsdinp.h" 2688#include "r12int.h" 2689C INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2690C 2691celena 2692 NORB1T = 0 2693 NORB2T = 0 2694 NRHFST = 0 2695 DO ISYM = 1, NSYM 2696 NORB1T = NORB1T + NORB1(ISYM) 2697 NORB2T = NORB2T + NORB2(ISYM) 2698 NRHFST = NRHFST + NRHFS(ISYM) 2699 ENDDO 2700 LUSIRG = -1 2701 CALL GPOPEN(LUSIRG,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,LDUM) 2702 REWIND LUSIRG 2703 2704 CALL MOLLAB(LABEL,LUSIRG,LUPRI) 2705 READ (LUSIRG) 2706C 2707 READ (LUSIRG) 2708 READ (LUSIRG) (WORK(i), I=1,NLAMDS) 2709C 2710 CALL GPCLOSE(LUSIRG,'KEEP') 2711 2712 2713 2714celena 2715 2716 NOCCT = NRHFT 2717 LUSIRG = -30 2718 KEND = NLAMDS + 2*NOCCT*NBAST + 14*NOCCT*NOCCT + 1 2719 IF (KEND .GT. LWORK) CALL STOPIT('MAKEUG',' ',LWORK,KEND) 2720 CALL GPOPEN(LUSIRG,'SIRIFC','OLD',' ','UNFORMATTED',IDUM,LDUM) 2721 REWIND LUSIRG 2722 2723 CALL MOLLAB(LABEL,LUSIRG,LUPRI) 2724 READ (LUSIRG) 2725C 2726 READ (LUSIRG) 2727 READ (LUSIRG) (WORK(I), I=1,NLAMDS) 2728C 2729 CALL GPCLOSE(LUSIRG,'KEEP') 2730C 2731 IWORK = 1 2732 IGX = 1 2733 IUG = 1 2734 IF (R12PRP .AND. R12HYB) THEN 2735 DO ISYM = 1, NSYM 2736 IWORK = IWORK + NBAS(ISYM)*NRHFS(ISYM) 2737 CALL MPAB(WORK(IWORK),NBAS(ISYM),NVIR(ISYM), 2738 * NBAS(ISYM),NVIR(ISYM), 2739 * GX(IGX),NVIR(ISYM),NORB1(ISYM), 2740 * NVIR(ISYM),NORB1(ISYM), 2741 * UG(IUG),NBAS(ISYM),NORB1(ISYM)) 2742 IF (IPRINT .GT. 3) THEN 2743 IF (NORB1(ISYM)*NBAS(ISYM) .NE. 0) 2744 * WRITE(LUPRI,'(/A,I2)') ' U(NORB1) x X/ISYM =',ISYM 2745 CALL OUTPUT(UG(IUG),1,NBAS(ISYM),1,NORB1(ISYM), 2746 * NBAS(ISYM),NORB1(ISYM),1,LUPRI) 2747 WRITE(LUPRI,'(/A,I2)') ' G(NORB1) x X/ISYM =',ISYM 2748 CALL OUTPUT(GX(IGX),1,NBAS(ISYM),1,NORB1(ISYM), 2749 * NBAS(ISYM),NORB1(ISYM),1,LUPRI) 2750 ENDIF 2751 IWORK = IWORK + NBAS(ISYM)*NVIRS(ISYM) 2752 IGX = IGX + NVIR(ISYM)*(NORB1(ISYM)-NRHFFR(ISYM)) 2753 IUG = IUG + NBAS(ISYM)*(NORB1(ISYM)-NRHFFR(ISYM)) 2754 ENDDO 2755 ELSE 2756 DO 100 ISYM = 1, NSYM 2757 IWORK = IWORK + NBAS(ISYM)*NRHFS(ISYM) 2758 CALL MPAB(WORK(IWORK),NBAS(ISYM),NVIR(ISYM), 2759 * NBAS(ISYM),NVIR(ISYM), 2760 * GX(IGX),NVIR(ISYM),NRHF(ISYM), 2761 * NVIR(ISYM),NRHF(ISYM), 2762 * UG(IUG),NBAS(ISYM),NRHF(ISYM)) 2763 IF (IPRINT .GT. 3) THEN 2764 IF (NRHF(ISYM)*NBAS(ISYM) .NE. 0) 2765 * WRITE(LUPRI,'(/A,I2)') ' U x X/ISYM =',ISYM 2766 CALL OUTPUT(UG(IUG),1,NBAS(ISYM),1,NRHF(ISYM), 2767 * NBAS(ISYM),NRHF(ISYM),1,LUPRI) 2768 ENDIF 2769 IWORK = IWORK + NBAS(ISYM)*NVIRS(ISYM) 2770 IGX = IGX + NVIR(ISYM)*NRHF(ISYM) 2771 IUG = IUG + NBAS(ISYM)*NRHF(ISYM) 2772 100 CONTINUE 2773 ENDIF 2774 NTOTGU = IUG - 1 2775C 2776C Add contribution from natural virtual orbitals (WK/UniKA/22-11-2005). 2777C 2778 IF (NATVIR) THEN 2779 NR12xVIR = 0 2780 DO ISYM = 1, NSYM 2781 NR12xVIR = NR12xVIR + (NORB1(ISYM)-NRHFSA(ISYM))*NRXR12(ISYM) 2782 END DO 2783 KFCKMIX = KEND 2784 KEND = KFCKMIX + NR12xVIR 2785 KFCKHF = KEND + NR12xVIR * NBAST 2786 KEND = KFCKHF 2787 IF (KEND .GT. LWORK) CALL STOPIT('MAKEUG',' ',LWORK,KEND) 2788 LUNIT = -1 2789 CALL GPOPEN(LUNIT,'R12FVIR','UNKNOWN',' ','UNFORMATTED',IDUMMY, 2790 & .FALSE.) 2791 READ (LUNIT) (WORK(KFCKMIX+I-1), I=1, NR12xVIR) 2792 CALL GPCLOSE(LUNIT,'KEEP') 2793 XALPH = 1.0D0 2794 IF (R12NOB) THEN 2795 XBETA = 0.0D0 2796 ELSE 2797 XBETA = 1.0D0 2798 END IF 2799 IUG = 1 2800 JUG = 1 2801 IWORK = 1 2802 DO ISYM = 1, NSYM 2803 IUG = IUG + NBAS(ISYM)*NRHFA(ISYM) 2804 IWORK = IWORK + NBAS(ISYM)*(NRHFS(ISYM)+NRHFSA(ISYM)) 2805 NREAVR = NORB1(ISYM)-NRHFSA(ISYM) 2806 IF (NRXR12(ISYM).GT.0) 2807 * CALL DGEMM('N','N',NBAS(ISYM),NRXR12(ISYM),NREAVR,XALPH, 2808 * WORK(IWORK),NBAS(ISYM),WORK(KFCKMIX),NREAVR, 2809 * XBETA,UG(IUG),NBAS(ISYM)) 2810 IF (IPRINT .GT. 4) THEN 2811 IF (NRXR12(ISYM)*NBAS(ISYM) .NE. 0) THEN 2812 WRITE(LUPRI,'(/A,I2)') ' MO MATRIX/ISYM =',ISYM 2813 CALL OUTPUT(WORK(IWORK),1,NBAS(ISYM),1,NREAVR, 2814 * NBAS(ISYM),NREAVR,1,LUPRI) 2815 WRITE(LUPRI,'(/A,I2)') ' FOCK MATRIX/ISYM =',ISYM 2816 CALL OUTPUT(WORK(KFCKMIX),1,NREAVR,1,NRXR12(ISYM), 2817 * NREAVR,NRXR12(ISYM),1,LUPRI) 2818 ENDIF 2819 WRITE(LUPRI,'(/A,I2)') ' U x (X + F)/ISYM =',ISYM 2820 CALL OUTPUT(UG(JUG),1,NBAS(ISYM),1,NRHF(ISYM), 2821 * NBAS(ISYM),NRHF(ISYM),1,LUPRI) 2822 JUG = JUG + NBAS(ISYM)*NRHF(ISYM) 2823 ENDIF 2824 IUG = IUG + NBAS(ISYM)*NRXR12(ISYM) 2825 IWORK = IWORK + NBAS(ISYM)*(NREAVR+NORB2(ISYM)) 2826 KFCKMIX = KFCKMIX + NREAVR*NRXR12(ISYM) 2827 END DO 2828 IF (IWORK-1 .NE. NLAMDS) 2829 * CALL STOPIT('MAKEUGa',' ',IWORK-1,NLAMDS) 2830 IF (NTOTGU .NE. IUG-1) 2831 * CALL STOPIT('MAKEUGb',' ',NTOTGU,IUG-1) 2832 END IF 2833C 2834 NTOTGU = IUG - 1 2835 WRITE(GUNAM,'(A6,I1)') 'GUMAT.',IPBAS 2836 LU43 = -43 2837 CALL GPOPEN(LU43,GUNAM,'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 2838 REWIND(LU43) 2839 WRITE(LU43) NTOTGU 2840 WRITE(LU43) (UG(I), I = 1, NTOTGU) 2841 CALL GPCLOSE(LU43,'KEEP') 2842C 2843 ISMOU = 1 2844 ISMO0 = 1 2845 ISMO1 = ISMO0 + NLAMDS 2846celena 2847 ISMO2 = ISMO1 + NORB1T*NBAST 2848celena 2849 CALL DZERO(WORK(ISMO1),NORB1T*NBAST) 2850 CALL DZERO(WORK(ISMO2),NORB1T*NBAST) 2851 DO 505 ISYM=1,NSYM 2852 ISMO0 = ISMO0 + NBAS(ISYM)*NRHFFR(ISYM) 2853 DO 510 I = 1, NRHF(ISYM) 2854 CALL DCOPY(NBAS(ISYM),WORK(ISMO0),1,WORK(ISMO1),1) 2855 CALL DCOPY(NBAS(ISYM),UG(ISMOU),1,WORK(ISMO2),1) 2856 ISMO1 = ISMO1 + NBAST 2857 ISMO2 = ISMO2 + NBAST 2858 ISMO0 = ISMO0 + NBAS(ISYM) 2859 ISMOU = ISMOU + NBAS(ISYM) 2860 510 CONTINUE 2861 ISMO0 = ISMO0 + NVIR(ISYM)*NBAS(ISYM) 2862 IF (R12PRP) ISMOU=ISMOU + (NORB1(ISYM)-NRHFS(ISYM))*NBAS(ISYM) 2863 ISMO1 = ISMO1 + NBAS(ISYM) 2864 ISMO2 = ISMO2 + NBAS(ISYM) 2865 505 CONTINUE 2866celena 2867 IF (R12PRP) THEN 2868 ISMO5 = 1 2869 ISMO3 = 1 + NLAMDS + NOCCT*NBAST 2870 ISMO4 = 1 2871 ISMO6 = ISMO3 + NORB1T*NBAST 2872 CALL DZERO(WORK(ISMO3),(NORB1T-NRHFST)*NBAST) 2873 DO ISYM=1,NSYM 2874 ISMO4 = ISMO4 + 2*NBAS(ISYM)*NRHFS(ISYM) 2875 ISMO5 = ISMO5 + NBAS(ISYM)*NRHF(ISYM) 2876 DO I = 1, NORB1(ISYM)-NRHFS(ISYM) 2877 CALL DCOPY(NBAS(ISYM),WORK(ISMO4),1,WORK(ISMO3),1) 2878 CALL DCOPY(NBAS(ISYM),UG(ISMO5),1,WORK(ISMO6),1) 2879 ISMO3 = ISMO3 + NBAST 2880 ISMO4 = ISMO4 + NBAS(ISYM) 2881 ISMO6 = ISMO6 + NBAST 2882 ISMO5 = ISMO5 + NBAS(ISYM) 2883 ENDDO 2884 ISMO3 = ISMO3 + NBAS(ISYM) 2885 ISMO6 = ISMO6 + NBAS(ISYM) 2886 ISMO4 = ISMO4 + NORB2(ISYM)*NBAS(ISYM) 2887 ENDDO 2888 ENDIF 2889celena 2890 2891 ISMO0 = 1 2892 ISMA1 = ISMO0 + 2*NOCCT*NBAST 2893 ISMO1 = ISMO0 + NLAMDS 2894 ISMO2 = ISMO1 + NORB1T*NBAST 2895 ISMS0 = ISMO2 + NORB1T*NBAST 2896 ISMX1 = ISMS0 + NOCCT*NOCCT*2 2897 ISMY1 = ISMX1 + NOCCT*NOCCT*2 2898 ISMZ1 = ISMY1 + NOCCT*NOCCT*2 2899 ISMX2 = ISMZ1 + NOCCT*NOCCT*2 2900 ISMY2 = ISMX2 + NOCCT*NOCCT*2 2901 ISMZ2 = ISMY2 + NOCCT*NOCCT*2 2902 IEND = ISMZ2 + NOCCT*NOCCT*2 2903celena 2904 IF (R12PRP .AND. R12HYB) THEN 2905 ISMRO0 = 1 2906 ISMRA1 = ISMRO0 + 2*NOCCT*NBAST 2907 ISMRO1 = ISMRO0 + NLAMDS 2908 ISMRO2 = ISMRO1 + NORB1T*NBAST 2909 ISMRS0 = ISMRO2 + NORB1T*NBAST 2910 ISMRX1 = ISMRS0 + 2*NORB1T*NORB1T 2911 ISMRY1 = ISMRX1 + 2*NORB1T*NORB1T 2912 ISMRZ1 = ISMRY1 + 2*NORB1T*NORB1T 2913 ISMRX2 = ISMRZ1 + 2*NORB1T*NORB1T 2914 ISMRY2 = ISMRX2 + 2*NORB1T*NORB1T 2915 ISMRZ2 = ISMRY2 + 2*NORB1T*NORB1T 2916 IEND = ISMRZ2 + 2*NORB1T*NORB1T 2917 ENDIF 2918celena 2919 LQQ = LWORK - IEND 2920 IF (IPRINT .GT. 3) THEN 2921 CALL AROUND('Active orbitals') 2922 CALL OUTPUT(WORK(ISMO1),1,NBAST,1,NOCCT,NBAST,NOCCT,1,LUPRI) 2923 CALL AROUND('Primed active orbitals') 2924 CALL OUTPUT(WORK(ISMO2),1,NBAST,1,NOCCT,NBAST,NOCCT,1,LUPRI) 2925 END IF 2926 IF (LQQ .LT. 0) CALL STOPIT('MAKEUG','R12QQ',LWORK,IEND) 2927 CALL R12QQ(QQ,WORK(ISMS0),WORK(ISMX1),WORK(ISMY1),WORK(ISMZ1), 2928 * WORK(ISMX2),WORK(ISMY2),WORK(ISMZ2), 2929 * WORK(ISMO2),WORK(ISMO1), 2930 * WORK(IEND),LQQ,NBAST,NOCCT) 2931 LUNIT = -1 2932celena 2933 IF (R12PRP .AND. R12HYB) THEN 2934 NVIR0T = NORB1T-NRHFST 2935 2936 CALL DZERO(WORK(ISMRS0),14*NORB1T*NORB1T) 2937 CALL DZERO(QV,NVIR0T*NRHFT**3) 2938 CALL R12QV(.FALSE.,.FALSE.,QV,WORK(ISMRS0),WORK(ISMRX1), 2939 * WORK(ISMRY1), 2940 * WORK(ISMRZ1), 2941 * WORK(ISMRX2),WORK(ISMRY2), 2942 * WORK(ISMRZ2), 2943 * WORK(ISMRO2),WORK(ISMRO1), 2944 * WORK(IEND),LQQ,NBAST,NORB1T,NOCCT,NVIR0T,.FALSE.) 2945 LUNIT = -1 2946 CALL GPOPEN(LUNIT,'AUXQSA12','UNKNOWN','SEQUENTIAL', 2947 & 'FORMATTED',IDUMMY,LDUMMY) 2948 WRITE(LUNIT,'(4E30.20)') (QV(J+1), J = 0, (NORB1T-NRHFST) 2949 & *NOCCT**3 - 1) 2950 CALL GPCLOSE(LUNIT,'KEEP') 2951 2952 CALL DZERO(WORK(ISMRS0),14*NORB1T*NORB1T) 2953 CALL DZERO(QV,NVIR0T*NRHFT**3) 2954 CALL R12QV(.FALSE.,.FALSE.,QV,WORK(ISMRS0),WORK(ISMRX1), 2955 * WORK(ISMRY1), 2956 * WORK(ISMRZ1), 2957 * WORK(ISMRX2),WORK(ISMRY2), 2958 * WORK(ISMRZ2), 2959 * WORK(ISMRO2),WORK(ISMRO1), 2960 * WORK(IEND),LQQ,NBAST,NORB1T,NOCCT,NVIR0T,.TRUE.) 2961 LUNIT = -1 2962 CALL GPOPEN(LUNIT,'AUXQSAJ12','UNKNOWN','SEQUENTIAL', 2963 & 'FORMATTED',IDUMMY,LDUMMY) 2964 WRITE(LUNIT,'(4E30.20)') (QV(J+1), J = 0, (NVIR0T) 2965 & *NRHFT**3 - 1) 2966 CALL GPCLOSE(LUNIT,'KEEP') 2967 2968 CALL DZERO(WORK(ISMRS0),14*NORB1T*NORB1T) 2969 CALL DZERO(QV,NVIR0T*NRHFT**3) 2970 CALL R12QV(.FALSE.,.TRUE.,QV,WORK(ISMRS0),WORK(ISMRX1), 2971 * WORK(ISMRY1), 2972 * WORK(ISMRZ1), 2973 * WORK(ISMRX2),WORK(ISMRY2), 2974 * WORK(ISMRZ2), 2975 * WORK(ISMRO2),WORK(ISMRO1), 2976 * WORK(IEND),LQQ,NBAST,NORB1T,NOCCT,NVIR0T,.FALSE.) 2977 LUNIT = -1 2978 CALL GPOPEN(LUNIT,'AUXQSAI12','UNKNOWN','SEQUENTIAL', 2979 & 'FORMATTED',IDUMMY,LDUMMY) 2980 WRITE(LUNIT,'(4E30.20)') (QV(J+1), J = 0, (NVIR0T) 2981 & *NRHFT**3 - 1) 2982 CALL GPCLOSE(LUNIT,'KEEP') 2983 2984 ENDIF 2985celena 2986 RETURN 2987 END 2988C /* Deck r12qq */ 2989 SUBROUTINE R12QQ(QQ,S0,X1,Y1,Z1,X2,Y2,Z2,UU,ZZ, 2990 * WORK,LWORK,NBAST,NOCCT) 2991C 2992C This subroutine computes the integrals (i'k|r12**2|jl). 2993C 2994C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 2995C 2996#include "implicit.h" 2997#include "priunit.h" 2998 DIMENSION QQ(*),X1(*),Y1(*),Z1(*),X2(*),Y2(*),Z2(*), 2999 * UU(*),ZZ(*),WORK(LWORK),S0(*) 3000 LOGICAL FOUND 3001 KEND = 1 + NBAST*NBAST 3002 LWRK = LWORK - KEND 3003 CALL RDPROP('CM000000',WORK,FOUND) 3004 CALL UTHZ(WORK,UU,ZZ,S0,WORK(KEND),LWRK,NBAST,NOCCT) 3005 CALL RDPROP('CM010000',WORK,FOUND) 3006 CALL UTHZ(WORK,UU,ZZ,X1,WORK(KEND),LWRK,NBAST,NOCCT) 3007 CALL RDPROP('CM000100',WORK,FOUND) 3008 CALL UTHZ(WORK,UU,ZZ,Y1,WORK(KEND),LWRK,NBAST,NOCCT) 3009 CALL RDPROP('CM000001',WORK,FOUND) 3010 CALL UTHZ(WORK,UU,ZZ,Z1,WORK(KEND),LWRK,NBAST,NOCCT) 3011 CALL RDPROP('CM020000',WORK,FOUND) 3012 CALL UTHZ(WORK,UU,ZZ,X2,WORK(KEND),LWRK,NBAST,NOCCT) 3013 CALL RDPROP('CM000200',WORK,FOUND) 3014 CALL UTHZ(WORK,UU,ZZ,Y2,WORK(KEND),LWRK,NBAST,NOCCT) 3015 CALL RDPROP('CM000002',WORK,FOUND) 3016 CALL UTHZ(WORK,UU,ZZ,Z2,WORK(KEND),LWRK,NBAST,NOCCT) 3017 CALL MAKEQQ(S0,X1,Y1,Z1,X2,Y2,Z2,QQ,NOCCT) 3018 RETURN 3019 END 3020C /* Deck uthz */ 3021 SUBROUTINE UTHZ(AA,UU,ZZ,BB,WORK,LWORK,NBAST,NOCCT) 3022C 3023C B(1) = Z(Transpose) * A * Z 3024C B(2) = U(Transpose) * A * Z 3025C 3026C On input, AA is the upper triangle of the symmetric 3027C matrix A. B is returned through the array BB as 3028C square matrix. 3029C 3030C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 3031C 3032#include "implicit.h" 3033#include "priunit.h" 3034 DIMENSION AA(*),UU(NBAST,NOCCT),ZZ(NBAST,NOCCT), 3035 * WORK(NBAST,NOCCT),BB(NOCCT,NOCCT,2) 3036 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3037 IF (NBAST*NOCCT .GT. LWORK) 3038 * CALL STOPIT('UTHZ',' ',LWORK,NBAST*NOCCT) 3039 CALL DZERO(WORK,NBAST*NOCCT) 3040 DO 103 K = 1, NBAST 3041 DO 102 J = 1, NOCCT 3042 ZZKJ = ZZ(K,J) 3043 DO 101 I = 1, NBAST 3044 IK = INDEX(I,K) 3045 WORK(I,J) = WORK(I,J) + AA(IK) * ZZKJ 3046 101 CONTINUE 3047 102 CONTINUE 3048 103 CONTINUE 3049 DO 202 I = 1, NOCCT 3050 DO 201 J = 1, NOCCT 3051 BB(I,J,1) = DDOT(NBAST,ZZ(1,I),1,WORK(1,J),1) 3052 BB(I,J,2) = DDOT(NBAST,UU(1,I),1,WORK(1,J),1) 3053 201 CONTINUE 3054 202 CONTINUE 3055 RETURN 3056 END 3057C /* Deck makeqq */ 3058 SUBROUTINE MAKEQQ(S0,X1,Y1,Z1,X2,Y2,Z2,QQ,NOCCT) 3059C 3060C This subroutine computes the integrals (i'k|r12**2|jl). 3061C 3062C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 3063C 3064#include "implicit.h" 3065 PARAMETER (D2 = 2D0) 3066#include "priunit.h" 3067 DIMENSION X1(NOCCT,NOCCT,2),Y1(NOCCT,NOCCT,2), 3068 * Z1(NOCCT,NOCCT,2),X2(NOCCT,NOCCT,2), 3069 * Y2(NOCCT,NOCCT,2),Z2(NOCCT,NOCCT,2), 3070 * S0(NOCCT,NOCCT,2), 3071 * QQ(NOCCT,NOCCT,NOCCT,NOCCT) 3072 DO 100 I = 1, NOCCT 3073 DO 200 J = 1, NOCCT 3074 DO 300 K = 1, NOCCT 3075 DO 400 L = 1, NOCCT 3076 QQ(I,J,K,L) = 3077 * (X2(I,K,2) + Y2(I,K,2) + Z2(I,K,2)) * S0(J,L,1) + 3078 * (X2(J,L,1) + Y2(J,L,1) + Z2(J,L,1)) * S0(I,K,2) - 3079 * D2 * ( X1(I,K,2) * X1(J,L,1) + 3080 * Y1(I,K,2) * Y1(J,L,1) + 3081 * Z1(I,K,2) * Z1(J,L,1) ) 3082 400 CONTINUE 3083 300 CONTINUE 3084 200 CONTINUE 3085 100 CONTINUE 3086 RETURN 3087 END 3088C /* Deck qcal */ 3089 SUBROUTINE QCAL(QQ,R2AM,U2AM,WORK,LWORK,IFLAG,IANSAZ) 3090C 3091C Driver routine for the computation of the Q term 3092C of R12/B theory. 3093C 3094C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 3095C 3096#include "implicit.h" 3097#include "priunit.h" 3098#include "ccorb.h" 3099#include "ccsdsym.h" 3100#include "r12int.h" 3101 DIMENSION R2AM(*), U2AM(*), WORK(*), QQ(*) 3102C 3103 NOCDIM = NRHFT 3104 NSPAIR = NOCDIM * (NOCDIM + 1) / 2 3105C 3106 IR11 = 1 3107 IQ11 = IR11 + 2 * NOCDIM**4 3108 IRS11 = IQ11 + 2 * NOCDIM**4 3109 IQS11 = IRS11 + NSPAIR**2 3110 IRT11 = IQS11 + NSPAIR**2 3111 IQT11 = IRT11 + NSPAIR**2 3112 IEND = IQT11 + NSPAIR**2 3113 IF (IEND .GT. LWORK) 3114 *CALL STOPIT('QCAL',' ',LWORK,IEND) 3115 CALL QCAL1(R2AM,U2AM,QQ,WORK(IR11),WORK(IQ11), 3116 * WORK(IRS11),WORK(IQS11),WORK(IRT11), 3117 * WORK(IQT11),NOCDIM,NSPAIR,IFLAG,IANSAZ) 3118 RETURN 3119 END 3120C /* Deck qcal1 */ 3121 SUBROUTINE QCAL1(R2AM,U2AM,QQ11,R11,Q11,RS11,RT11,QS11,QT11, 3122 * NOCDIM,NSPAIR,IFLAG,IANSAZ) 3123C 3124C This subroutine computes the Q term of R12/B: 3125C 3126C (i'k|r12**2|jl) = Sum_pq (i'p'|r12|jq')(kp'|r12|lq') 3127C 3128C Note that i' refers to the orbital |i'> = Sum_p' X(p'i) |p'> 3129C 3130C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 3131C 3132#include "implicit.h" 3133#include "priunit.h" 3134 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 3135 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0) 3136 DIMENSION R2AM(*), U2AM(*) 3137 REAL*8 Q11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 3138 * R11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 3139 * U1, R2, RR, QQIJKL, QQIJLK, UIJKL, UJIKL, UIJLK, UJILK 3140 DIMENSION QS11(NSPAIR,NSPAIR), RS11(NSPAIR,NSPAIR) 3141 DIMENSION QT11(NSPAIR,NSPAIR),RT11(NSPAIR,NSPAIR) 3142 DIMENSION QQ11(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 3143 INTEGER IDUM 3144 LOGICAL LDUM 3145#include "r12int.h" 3146#include "ccorb.h" 3147#include "ccsdsym.h" 3148#include "ccsdinp.h" 3149 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 3150C 3151 IF (R12NOB) THEN 3152 CALL DZERO(RS11,NSPAIR**2) 3153 CALL DZERO(RT11,NSPAIR**2) 3154 CALL DZERO(QS11,NSPAIR**2) 3155 CALL DZERO(QT11,NSPAIR**2) 3156 GOTO 998 3157 ENDIF 3158C 3159 FF = D1 / SQRT(D2) 3160 DO 90 I = 1, NOCDIM**4 3161 R11(I,1,1,1) = D0 3162 Q11(I,1,1,1) = D0 3163 90 CONTINUE 3164C 3165 IF (R12CBS) THEN 3166 FACX = - D1 3167 ELSE 3168 FACX = D1 3169 END IF 3170C 3171C CONSTRUCT PRODUCTS (IA|JB)*(KA|LB) 3172C 3173 DO 101 ISYMIA = 1, NSYM 3174 ISYMJB = ISYMIA 3175 JBOFF = NH1AM(ISYMJB) * (NH1AM(ISYMJB) + 1) / 2 3176 DO 102 ISYMI = 1, NSYM 3177 ISYMA = MULD2H(ISYMI,ISYMIA) 3178 DO 103 ISYMJ = 1, NSYM 3179 ISYMB = MULD2H(ISYMJ,ISYMJB) 3180 DO 104 ISYMKA = 1, NSYM 3181 ISYMLB = ISYMKA 3182 LBOFF = NH1AM(ISYMLB) * (NH1AM(ISYMLB) + 1) / 2 3183 ISYMK = MULD2H(ISYMA,ISYMKA) 3184 ISYML = MULD2H(ISYMB,ISYMLB) 3185 DO 105 I = 1, NRHF(ISYMI) 3186 KOFFI = IRHF(ISYMI) + I 3187 DO 106 K = 1, NRHF(ISYMK) 3188 KOFFK = IRHF(ISYMK) + K 3189 DO 107 A = 1, NVIR(ISYMA) 3190 IF (ONEAUX) THEN 3191 IF (A .LE. NORB1(ISYMA)) THEN 3192 NAI = IH1AM(ISYMA,ISYMI) + 3193 * NORB1(ISYMA)*(I-1) + A 3194 NAK = IH1AM(ISYMA,ISYMK) + 3195 * NORB1(ISYMA)*(K-1) + A 3196 ELSE 3197 NAI = IG1AM(ISYMA,ISYMI) + 3198 * NORB2(ISYMA)*(I-1) + A - 3199 * NORB1(ISYMA) 3200 NAK = IG1AM(ISYMA,ISYMK) + 3201 * NORB2(ISYMA)*(K-1) + A - 3202 * NORB1(ISYMA) 3203 END IF 3204 ELSE 3205 NAI = IT1AM(ISYMA,ISYMI) + 3206 * NVIR(ISYMA)*(I-1) + A 3207 NAK = IT1AM(ISYMA,ISYMK) + 3208 * NVIR(ISYMA)*(K-1) + A 3209 END IF 3210 DO 108 J = 1, NRHF(ISYMJ) 3211 KOFFJ = IRHF(ISYMJ) + J 3212 DO 109 L = 1, NRHF(ISYML) 3213 KOFFL = IRHF(ISYML) + L 3214 IF (ONEAUX) THEN 3215 NENDB = NORB1(ISYMB) 3216 ELSE 3217 NENDB = NVIR(ISYMB) 3218 END IF 3219 DO 110 B = 1, NENDB 3220 IF (ONEAUX) THEN 3221 NBJ = IH1AM(ISYMB,ISYMJ) + 3222 * NORB1(ISYMB)*(J-1) + B 3223 NBL = IH1AM(ISYMB,ISYML) + 3224 * NORB1(ISYMB)*(L-1) + B 3225 IF (A .LE. NORB1(ISYMA)) THEN 3226 NAIBJ = IH2AM(ISYMIA,ISYMJB) 3227 * + INDEX(NAI,NBJ) 3228 NAKBL = IH2AM(ISYMKA,ISYMLB) 3229 * + INDEX(NAK,NBL) 3230 ELSE 3231 NAIBJ = IH2AM(ISYMIA,ISYMJB) 3232 * + NH1AM(ISYMJB)*(NAI-1) 3233 * + NBJ + JBOFF 3234 NAKBL = IH2AM(ISYMKA,ISYMLB) 3235 * + NH1AM(ISYMLB)*(NAK-1) 3236 * + NBL + LBOFF 3237 END IF 3238 ELSE 3239 NBJ = IT1AM(ISYMB,ISYMJ) + 3240 * NVIR(ISYMB)*(J-1) + B 3241 NBL = IT1AM(ISYMB,ISYML) + 3242 * NVIR(ISYMB)*(L-1) + B 3243 NAIBJ = IU2AM(ISYMIA,ISYMJB) 3244 * + NT1AM(ISYMJB)*(NAI-1) 3245 * + NBJ 3246 NAKBL = IT2AM(ISYMKA,ISYMLB) + 3247 * INDEX(NAK,NBL) 3248 END IF 3249 U1 = U2AM(NAIBJ) 3250 R2 = R2AM(NAKBL) 3251 RR = U1 * R2 3252C 3253 GOTO (113,114), IFLAG 3254 CALL QUIT('INVALID IFLAG IN QCAL1') 3255C 3256 113 CONTINUE 3257 IF (A .GT. NORB1(ISYMA)) GOTO 112 3258 IF (B .GT. NORB1(ISYMB)) GOTO 112 3259C 3260 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 3261 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 3262C 3263 GOTO 111 3264 112 CONTINUE 3265C 3266 IF (A .GT. NORB1(ISYMA) .AND. 3267 * B .GT. NORB1(ISYMB)) GOTO 111 3268C 3269 R11(KOFFI,KOFFJ,KOFFK,KOFFL) = 3270 * R11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 3271 IF (ONEAUX) THEN 3272 R11(KOFFJ,KOFFI,KOFFL,KOFFK) = 3273 * R11(KOFFJ,KOFFI,KOFFL,KOFFK) + RR 3274 END IF 3275C 3276 GOTO 111 3277C 3278 114 CONTINUE 3279 IF (R12CBS) THEN 3280 NQUEA = 0 3281 NQUEB = 0 3282 ELSE 3283 NQUEA = NORB1(ISYMA) 3284 NQUEB = NORB1(ISYMB) 3285 END IF 3286 IF (A .LE. NRHFSA(ISYMA) .AND. 3287 * B .LE. NRHFSA(ISYMB)) THEN 3288C 3289C Sum |ij><ij| 3290C 3291 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 3292 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) - RR 3293C 3294 END IF 3295 IF (A .LE. NRHFSA(ISYMA) .AND. 3296 * B .GT. NQUEB) THEN 3297C 3298C Sum |iq'><iq'| 3299C 3300 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 3301 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 3302C 3303 END IF 3304 IF (A .GT. NQUEA .AND. 3305 * B .LE. NRHFSA(ISYMB)) THEN 3306C 3307C Sum |p'j><p'j| 3308C 3309 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 3310 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 3311 IF (ONEAUX.AND.A.GT.NORB1(ISYMA)) 3312 * Q11(KOFFJ,KOFFI,KOFFL,KOFFK) = 3313 * Q11(KOFFJ,KOFFI,KOFFL,KOFFK) + RR 3314C 3315 END IF 3316 IF (A .GT. NRHFSA(ISYMA) .AND. 3317 * A .LE. NORB1(ISYMA) .AND. 3318 * B .GT. NRHFSA(ISYMB) .AND. 3319 * B .LE. NORB1(ISYMB)) THEN 3320C 3321C Sum |ab><ab| 3322C 3323cWK IF (IANSAZ .EQ. 3) 3324cWK * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 3325cWK * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 3326C 3327 END IF 3328 111 CONTINUE 3329 110 CONTINUE 3330 109 CONTINUE 3331 108 CONTINUE 3332 107 CONTINUE 3333 106 CONTINUE 3334 105 CONTINUE 3335 104 CONTINUE 3336 103 CONTINUE 3337 102 CONTINUE 3338 101 CONTINUE 3339C 3340 IJ = 0 3341 DO 300 I = 1, NOCDIM 3342 DO 301 J = 1, I 3343 IJ = IJ + 1 3344 KL = 0 3345 DO 302 K = 1, NOCDIM 3346 DO 303 L = 1, K 3347 KL = KL + 1 3348C 3349 UIJKL = QQ11(I,J,K,L) 3350 UJILK = QQ11(J,I,L,K) 3351 UIJLK = QQ11(I,J,L,K) 3352 UJIKL = QQ11(J,I,K,L) 3353C 3354 QQIJKL = UIJKL + UJILK 3355 QQIJLK = UIJLK + UJIKL 3356 3357 RR = (QQIJKL + QQIJLK) 3358 * - (Q11(I,J,K,L) + Q11(I,J,L,K)) 3359 IF (.NOT. ONEAUX) 3360 * RR = RR - (Q11(J,I,L,K) + Q11(J,I,K,L)) 3361 QS11(KL,IJ) = RR 3362 IF (I .EQ. J) QS11(KL,IJ) = FF * QS11(KL,IJ) 3363 IF (K .EQ. L) QS11(KL,IJ) = FF * QS11(KL,IJ) 3364 RR = (QQIJKL - QQIJLK) 3365 * - (Q11(I,J,K,L) - Q11(I,J,L,K)) 3366 IF (.NOT. ONEAUX) 3367 * RR = RR - (Q11(J,I,L,K) - Q11(J,I,K,L)) 3368 QT11(KL,IJ) = RR 3369 QS11(KL,IJ) = QS11(KL,IJ) * DP25 * DP5 3370 QT11(KL,IJ) = QT11(KL,IJ) * DP75 * DP5 3371C 3372 IF (IFLAG .EQ. 2) GOTO 303 3373 RR = (QQIJKL + QQIJLK) 3374 * + (Q11(I,J,K,L) + Q11(I,J,L,K)) * FACX 3375 * - (R11(I,J,K,L) + R11(I,J,L,K)) 3376 IF (.NOT. ONEAUX) RR = RR 3377 * + (Q11(J,I,L,K) + Q11(J,I,K,L)) * FACX 3378 * - (R11(J,I,L,K) + R11(J,I,K,L)) 3379 RS11(KL,IJ) = RR 3380 IF (I .EQ. J) RS11(KL,IJ) = FF * RS11(KL,IJ) 3381 IF (K .EQ. L) RS11(KL,IJ) = FF * RS11(KL,IJ) 3382 RR = (QQIJKL - QQIJLK) 3383 * + (Q11(I,J,K,L) - Q11(I,J,L,K)) * FACX 3384 * - (R11(I,J,K,L) - R11(I,J,L,K)) 3385 IF (.NOT. ONEAUX) RR = RR 3386 * + (Q11(J,I,L,K) - Q11(J,I,K,L)) * FACX 3387 * - (R11(J,I,L,K) - R11(J,I,K,L)) 3388 RT11(KL,IJ) = RR 3389 RS11(KL,IJ) = RS11(KL,IJ) * DP25 * DP5 3390 RT11(KL,IJ) = RT11(KL,IJ) * DP75 * DP5 3391 303 CONTINUE 3392 302 CONTINUE 3393 301 CONTINUE 3394 300 CONTINUE 3395C 3396 CALL ERISFK(QS11,NSPAIR,1) 3397 CALL ERISFK(QT11,NSPAIR,1) 3398 IF (IFLAG .EQ. 1) THEN 3399 CALL ERISFK(RS11,NSPAIR,1) 3400 CALL ERISFK(RT11,NSPAIR,1) 3401 END IF 3402C 3403 IF (R12OLD) THEN 3404 CALL DCOPY(NSPAIR*NSPAIR,QS11,1,RS11,1) 3405 CALL DCOPY(NSPAIR*NSPAIR,QT11,1,RT11,1) 3406 ENDIF 3407C 3408 IF (IPRINT .LE. 3) GOTO 998 3409 GOTO (991,992), IFLAG 3410 CALL QUIT('INVALID IFLAG IN QCAL1') 3411 991 CALL AROUND('OLD singlet <Q> matrix') 3412 CALL OUTPUT(QS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 3413 CALL AROUND('OLD triplet <Q> matrix') 3414 CALL OUTPUT(QT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 3415 IF (R12OLD) GOTO 998 3416 CALL AROUND('NEW singlet <Q> matrix') 3417 CALL OUTPUT(RS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 3418 CALL AROUND('NEW triplet <Q> matrix') 3419 CALL OUTPUT(RT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 3420 GOTO 998 3421 992 CALL AROUND('Singlet <Q@> matrix') 3422 CALL OUTPUT(QS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 3423 CALL AROUND('Triplet <Q@> matrix') 3424 CALL OUTPUT(QT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 3425C 3426 998 LUMUL = -43 3427 GOTO (996,997), IFLAG 3428 CALL QUIT('INVALID IFLAG IN QCAL1') 3429 996 CALL GPOPEN(LUMUL,'QMATSO','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3430 WRITE(LUMUL,'(4E30.20)') QS11 3431 CALL GPCLOSE(LUMUL,'KEEP') 3432 CALL GPOPEN(LUMUL,'QMATTO','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3433 WRITE(LUMUL,'(4E30.20)') QT11 3434 CALL GPCLOSE(LUMUL,'KEEP') 3435 CALL GPOPEN(LUMUL,'QMATSN','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3436 WRITE(LUMUL,'(4E30.20)') RS11 3437 CALL GPCLOSE(LUMUL,'KEEP') 3438 CALL GPOPEN(LUMUL,'QMATTN','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3439 WRITE(LUMUL,'(4E30.20)') RT11 3440 CALL GPCLOSE(LUMUL,'KEEP') 3441 RETURN 3442 997 IF (IANSAZ .EQ. 2) THEN 3443 CALL GPOPEN(LUMUL,'QMATSP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3444 WRITE(LUMUL,'(4E30.20)') QS11 3445 CALL GPCLOSE(LUMUL,'KEEP') 3446 CALL GPOPEN(LUMUL,'QMATTP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3447 WRITE(LUMUL,'(4E30.20)') QT11 3448 CALL GPCLOSE(LUMUL,'KEEP') 3449 ELSE 3450 CALL GPOPEN(LUMUL,'QMATSQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3451 WRITE(LUMUL,'(4E30.20)') QS11 3452 CALL GPCLOSE(LUMUL,'KEEP') 3453 CALL GPOPEN(LUMUL,'QMATTQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 3454 WRITE(LUMUL,'(4E30.20)') QT11 3455 CALL GPCLOSE(LUMUL,'KEEP') 3456 END IF 3457 RETURN 3458 END 3459C /* Deck ccsd_r12 */ 3460 SUBROUTINE CCSD_R12(WORK,LWORK,WXRK,LWXRK,CCR12RSP) 3461C 3462C R12 driver. 3463C 3464C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 3465C 3466#include "implicit.h" 3467#include "maxorb.h" 3468#include "priunit.h" 3469#include "ccorb.h" 3470#include "iratdef.h" 3471#include "ccsdinp.h" 3472#include "ccsdsym.h" 3473#include "ccfro.h" 3474#include "ccsdio.h" 3475#include "cclr.h" 3476#include "ccfdgeo.h" 3477#include "cbirea.h" 3478#include "r12int.h" 3479#include "ccr12int.h" 3480 PARAMETER (D0 = 0D0) 3481 CHARACTER*6 TIMLAB, LABINT 3482 CHARACTER*10 APPROX 3483 DIMENSION WORK(LWORK), WXRK(LWXRK) 3484 INTEGER IDUM, LU43 3485 LOGICAL PROFIL, LDUM, TEMP_HERDIR, TEMP_DIRECT, TEMP_ONEAUX, LHTF, 3486 & MKVAJKL, CCR12RSP,PROJVIR 3487 CALL QENTER('CCSD_R12') 3488 TEMP_DIRECT = DIRECT 3489 DIRECT = .TRUE. 3490 MKVAJKL=.FALSE. 3491 LUSIRG = -30 3492 LU43 = -43 3493 PROFIL = .FALSE. 3494 R12CAL = R12CAL .OR. LMULBS 3495 IF (R12CAL) THEN 3496 R12OLD = .TRUE. 3497 DO ISYM = 1, 8 3498 R12OLD = R12OLD .AND. NORB2(ISYM) .EQ. 0 3499 IF (NVIRFR(ISYM) .NE. 0) 3500 & CALL QUIT('DO NOT FREEZE VIRTUALS IN R12 CALCULATIONS.') 3501 END DO 3502 NORXR = NORXR .OR. R12OLD .OR. R12HYB 3503 NORXR = NORXR .OR. R12OLD 3504 ONEAUX = (R12HYB .OR. R12NOB) .AND. .NOT. R12OLD 3505 & .AND. .NOT. R12EOR 3506 & .AND. .NOT. CCR12RSP 3507celena 3508 IF (R12PRP) ONEAUX = .FALSE. 3509celena 3510 R12HYB = R12HYB .OR. ONEAUX 3511 LABEL = 'FULLBAS ' 3512 write(lupri,*) 'Before INIT1 in ccsd_r12' 3513 call flshfo(lupri) ! asm 3514 CALL CCSD_INIT1(WORK,LWORK) 3515 write(lupri,*) 'After INIT1 in ccsd_r12' 3516 call flshfo(lupri) ! asm 3517 NORB1T = 0 3518 NORB2T = 0 3519 NNBASF = 0 3520 NSPAIR = NRHFT * (NRHFT + 1) / 2 3521 NOCXBS = 0 3522 NVCXBS = 0 3523 NBASTX = 0 3524 DO 230 ISYM = 1, NSYM 3525 NORB1T = NORB1T + NORB1(ISYM) 3526 NORB2T = NORB2T + NORB2(ISYM) 3527 NBASF = NORB1(ISYM) + NORB2(ISYM) 3528 NNBASF = NNBASF + NBASF * (NBASF + 1) / 2 3529 NOCXBS = NOCXBS + NRHF(ISYM) * NBAS(ISYM) 3530 NVCXBS = NVCXBS + NORB1(ISYM) * NBAS(ISYM) 3531 NBASTX = NBASTX + NBAS(ISYM) 3532 230 CONTINUE 3533C 3534 WRITE(LUPRI,'(1X,A,8I8)') 3535 * ' NRHF :',(NRHF(ISYM),ISYM=1,NSYM) 3536 WRITE(LUPRI,'(1X,A,8I8)') 3537 * ' NRHFS :',(NRHFS(ISYM),ISYM=1,NSYM) 3538 WRITE(LUPRI,'(1X,A,8I8)') 3539 * ' NORB1 :',(NORB1(ISYM),ISYM=1,NSYM) 3540 WRITE(LUPRI,'(1X,A,8I8)') 3541 * ' NORB2 :',(NORB2(ISYM),ISYM=1,NSYM) 3542 WRITE(LUPRI,'(/1X,A,I12)') 3543 * 'Number of primary orbitals : ',NORB1T 3544 WRITE(LUPRI,'(1X,A,I12)') 3545 * 'Number of secondary orbitals : ',NORB2T 3546 WRITE(LUPRI,'(1X,A,I12)') 3547 * 'Total number of orbitals : ',NORB2T + NORB1T 3548 WRITE(LUPRI,'(1X,A,I12)') 3549 * 'Total number of basis functions : ',NBASTX 3550 IF (ONEAUX) THEN 3551 WRITE(LUPRI,'(1X,A,I12)') 3552 * 'Number of (i|p) integrals : ',NH1AMX 3553 WRITE(LUPRI,'(1X,A,I12)') 3554 * 'Number of (i|p'') integrals : ',NT1AMX 3555 WRITE(LUPRI,'(1X,A,I12)') 3556 * 'Number of (i|u'') integrals : ',NOCXBS 3557 WRITE(LUPRI,'(1X,A,I12)') 3558 * 'Number of (Ip''|jL) integrals : ',NF2FRO(1) 3559 WRITE(LUPRI,'(1X,A,I12)') 3560 * 'Number of (ip''|jq) integrals : ',NH2AMX 3561 ELSE 3562 WRITE(LUPRI,'(1X,A,I12)') 3563 * 'Number of (i|p'') integrals : ',NT1AMX 3564 WRITE(LUPRI,'(1X,A,I12)') 3565 * 'Number of (i|u'') integrals : ',NOCXBS 3566 WRITE(LUPRI,'(1X,A,I12)') 3567 * 'Number of (Ip''|1/r12|jL) integrals : ',NF2FRO(1) 3568 WRITE(LUPRI,'(1X,A,I12)') 3569 * 'Number of (ip''|1/r12|jq'') integrals : ',NT2AMX 3570 WRITE(LUPRI,'(1X,A,I12)') 3571 * 'Number of (ip''|r12|jq'') integrals : ',NT2AMX 3572 WRITE(LUPRI,'(1X,A,I12)') 3573 * 'Number of (ip''|[T1,r12]|jq'') integrals : ',NU2AMX 3574 END IF 3575 WRITE(LUPRI,'(1X,A,G13.6)') 3576 * 'Threshold for R12 contributions =',VCLTHR 3577 WRITE(LUPRI,'(1X,A,G13.6/)') 3578 * 'Threshold for singular value decomposition =',SVDTHR 3579 IF (R12EOR) THEN 3580 IF (SLATER) THEN 3581 WRITE(LUPRI,'(1X,A/1X,A)') 3582 * 'Calculation with Gaussian-type geminals (GTGs)', 3583 * 'Expansion coefficients and exponents:' 3584 ELSE 3585 WRITE(LUPRI,'(1X,A/1X,A)') 3586 * 'Calculation with Gaussian-damped linear-r12 terms', 3587 * 'Expansion coefficients and exponents:' 3588 END IF 3589 DO I = 1, NTOGAM 3590 WRITE(LUPRI,'(I5,2F25.15)') I, GAMAB(I), GAMAX(I) 3591 END DO 3592 ENDIF 3593 IF (R12DIA) THEN 3594 WRITE(LUPRI,'(1X,A/A/)') 3595 * 'The B matrix is transformed into an orthogonal basis', 3596 * ' and inverted by diagonalization' 3597 ELSE IF (R12SVD) THEN 3598 WRITE(LUPRI,'(1X,A/)') 3599 * 'The B matrix is inverted by singular value decomposition' 3600 ELSE 3601 CALL QUIT('No solver selected') 3602 END IF 3603 IF (R12RST) WRITE(LUPRI,'(1X,A)') 3604 * 'The iterative local-MP2-R12 solver is restarted' 3605 IF (R12LEV .NE. D0) WRITE(LUPRI,'(1X,A,G13.6)') 3606 * 'Level-shift in local-MP2-R12 solver =',R12LEV 3607 WRITE(LUPRI,'(1X,78A1)') ('*',I=1,78) 3608 IF (R12HYB) THEN 3609 WRITE(LUPRI,'(1X,A,4(/1X,A))') 3610 * 'The MP2-R12/A and MP2-R12/B energies'// 3611 * ' are computed as described in:', 3612 * 'W. Klopper,'// 3613 * ' J. Chem. Phys. 120, 10890-10895 (2004).', 3614 * 'To avoid the hybride scheme, the keyword ''.NO HYB''', 3615 * 'must be specified in the Dalton input file.', 3616 * '----------------------------------------'// 3617 * '--------------------------------------' 3618 ELSE 3619 WRITE(LUPRI,'(1X,A,2(/1X,A))') 3620 * 'The MP2-R12/A and MP2-R12/B energies'// 3621 * ' are computed as described in:', 3622 * 'W. Klopper and C. C. M. Samson,'// 3623 * ' J. Chem. Phys. 116, 6397-6410 (2002).', 3624 * '----------------------------------------'// 3625 * '--------------------------------------' 3626 ENDIF 3627 IF (R12NOA) THEN 3628 WRITE(APPROX,'(A10)') ' B. ' 3629 ELSE IF (R12NOB) THEN 3630 WRITE(APPROX,'(A10)') ' A. ' 3631 ELSE 3632 WRITE(APPROX,'(A10)') 's A and B.' 3633 END IF 3634 IF (NOTTWO .AND. NOTTRE) THEN 3635 WRITE(LUPRI,'(A)') 3636 * ' Ansatz 1 is used and results are'// 3637 * ' printed for approximation'//APPROX 3638 ELSE IF (NOTONE .AND. NOTTRE) THEN 3639 WRITE(LUPRI,'(A)') 3640 * ' Ansatz 2 is used and results are'// 3641 * ' printed for approximation'//APPROX 3642 ELSE IF (NOTONE .AND. NOTTWO) THEN 3643 WRITE(LUPRI,'(A)') 3644 * ' Ansatz 3 is used and results are'// 3645 * ' printed for approximation'//APPROX 3646 ELSE IF (NOTTRE) THEN 3647 WRITE(LUPRI,'(A)') 3648 * ' Ansaetze 1 and 2 are used and results are'// 3649 * ' printed for approximation'//APPROX 3650 ELSE IF (NOTTWO) THEN 3651 WRITE(LUPRI,'(A)') 3652 * ' Ansaetze 1 and 3 are used and results are'// 3653 * ' printed for approximation'//APPROX 3654 ELSE IF (NOTONE) THEN 3655 WRITE(LUPRI,'(A)') 3656 * ' Ansaetze 2 and 3 are used and results are'// 3657 * ' printed for approximation'//APPROX 3658 ELSE 3659 WRITE(LUPRI,'(A)') 3660 * ' Ansaetze 1, 2, 3 are used and results are'// 3661 * ' printed for approximation'//APPROX 3662 END IF 3663 IF (R12EOR) WRITE(LUPRI,'(A,3(/A))') 3664 * '----------------------------------------'// 3665 * '--------------------------------------', 3666 * ' The Gaussian-damped R12 integrals'// 3667 * ' are computed as described in:', 3668 * ' C. C. M. Samson, W. Klopper and T. Helgaker,', 3669 * ' Comp. Phys. Commun. 149, 1-10 (2002).' 3670 WRITE(LUPRI,'(1X,78A1/)') ('*',I=1,78) 3671 CALL FLSHFO(LUPRI) 3672C 3673 KFOCKD = 1 3674 KFOCKQ = KFOCKD + NORBTS 3675 KSF = KFOCKQ + 2 * NRHFT ** 4 3676 KTF = KSF + MAX(NSPAIR ** 2, NNBASF) 3677 KT1AM = KTF + MAX(NSPAIR ** 2, NNBASF) 3678 KGXAM = KT1AM + NT1AMX * 10 3679 IF (R12PRP) THEN 3680 KQQAM = KGXAM + NT1VMX 3681 ELSE 3682 KQQAM = KGXAM + NT1AMX 3683 ENDIF 3684 IF (R12EOR) THEN 3685 KQQ2M = KQQAM + NRHFT ** 4 3686 KQQ4M = KQQ2M + NRHFT ** 4 3687 KQQ6M = KQQ4M + NRHFT ** 4 3688 KCMOM = KQQ6M + NRHFT ** 4 3689 KQVAM = KCMOM + 2*NLAMDS 3690 ELSE 3691 KQQ2M = KQQAM 3692 KQQ4M = KQQ2M 3693 KQQ6M = KQQ4M 3694 KCMOM = KQQ6M 3695 KQVAM = KCMOM + NRHFT**4 3696 ENDIF 3697celena 3698 KUGAM = KQVAM + NBAST*NRHFT**3 3699 3700celena 3701 IF (R12PRP) THEN 3702 KFRIN = KUGAM + NVCXBS 3703 ELSE 3704 KFRIN = KUGAM + NOCXBS 3705 ENDIF 3706C 3707celena 3708 KFRIN1= KFRIN + NFROVR(1) 3709 KS2AM = KFRIN1 + NFROVF(1) 3710celena 3711 IF (ONEAUX) THEN 3712 KT2AM = KS2AM + NH2AMX 3713 KU2AM = KT2AM 3714 KR2AM = KT2AM + NH2AMX 3715cWK IF (R12XXL) THEN 3716 KV2AM = KR2AM + NH2AMX 3717 KW2AM = KV2AM + NH2AMX 3718cWK ELSE 3719c KV2AM = KR2AM 3720c KW2AM = KV2AM 3721cWK END IF 3722 KEND1 = KW2AM + NH2AMX 3723 LENT2AM = NH2AMX 3724 ELSE 3725 KT2AM = KS2AM + NT2AMX 3726 KU2AM = KT2AM 3727 KR2AM = KT2AM + NT2AMX 3728 KV2AM = MAX(KR2AM + NT2AMX, KU2AM + NU2AMX) 3729 KW2AM = KV2AM + NT2AMX 3730 KEND1 = KW2AM + NT2AMX 3731 LENT2AM = NT2AMX 3732 END IF 3733 WRITE(LUPRI,'(A,I15,A/)')' Core memory required by CCSD_R12 :', 3734 & KEND1,' (double words)' 3735 CALL FLSHFO(LUPRI) 3736 LWRK1 = LWORK - KEND1 3737 IF (KEND1 .GT. LWORK) CALL STOPIT('CCSD_R12','R12 WORK SPACE', 3738 & LWORK,KEND1) 3739 IF (R12EOR) THEN 3740C 3741 LUSIFC = -1 3742 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 3743 & .FALSE.) 3744 REWIND LUSIFC 3745 CALL MOLLAB(LABEL,LUSIFC,LUPRI) 3746 READ (LUSIFC) 3747 READ (LUSIFC) 3748 READ (LUSIFC) (WORK(KCMOM+I-1),I=1,NLAMDS) 3749 CALL DCOPY(NLAMDS,WORK(KCMOM),1,WORK(KCMOM+NLAMDS),1) 3750 KCMO1 = KCMOM 3751 DO ISYM = 1, NSYM 3752 KCMO2 = KCMO1 + NRHFS(ISYM)*NBAS(ISYM) 3753 LENMO = NBAS(ISYM)*NRHFS(ISYM) 3754 CALL DCOPY(LENMO,WORK(KCMO1),1,WORK(KCMO2),1) 3755 KCMO1 = KCMO1 + 3756 & (NRHFS(ISYM)+NORB1(ISYM)+NORB2(ISYM))*NBAS(ISYM) 3757 ENDDO 3758 REWIND LUSIFC 3759 CALL MOLLAB(LABEL,LUSIFC,LUPRI) 3760 READ (LUSIFC) 3761 READ (LUSIFC) 3762 WRITE (LUSIFC) (WORK(KCMOM+I-1),I=1,NLAMDS) 3763 CALL GPCLOSE(LUSIFC,'KEEP') 3764C 3765 R12EIN = .TRUE. 3766C 3767C Correlation factor: f12 = r12 * exp(-gamma*r12**2) 3768C 3769C COMPUTE (IA|[[F12,T1],F12]|JB) INTEGRALS 3770C 3771 CALL TIMER('START ',TIMSTR,TIMEND) 3772 R12TRA = .TRUE. 3773 R12INT = .FALSE. 3774 R12SQR = .FALSE. 3775 U12INT = .FALSE. 3776 U21INT = .FALSE. 3777 V12INT = .FALSE. 3778 MBSMAX = 4 3779 INTGAC = 6 3780 CALL DZERO(WORK(KT1AM),NT1AMX) 3781 LHTF = .FALSE. 3782 write(Lupri,*) 'ccsd_r12.1' 3783 call flshfo(lupri) !asm 3784 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 3785 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 3786 write(Lupri,*) 'ccsd_r12.1' 3787 call flshfo(lupri) !asm 3788 CALL CCSD_IKJL(WORK(KQQ6M),NRHFT,WXRK(KS2AM)) 3789 WRITE(LUPRI,'(1X,A)') 3790 * 'Computation of (ik|[[f12,T1],f12]|jl) integrals done' 3791 WRITE(TIMLAB,1013) '#6 ',MBSMAX 3792 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 3793 CALL FLSHFO(LUPRI) 3794C 3795C COMPUTE (IA|F12**2|JB) INTEGRALS 3796C 3797 CALL TIMER('START ',TIMSTR,TIMEND) 3798 R12TRA = .TRUE. 3799 R12INT = .FALSE. 3800 R12SQR = .FALSE. 3801 U12INT = .FALSE. 3802 U21INT = .FALSE. 3803 V12INT = .FALSE. 3804 MBSMAX = 4 3805 INTGAC = 4 3806 CALL DZERO(WORK(KT1AM),NT1AMX) 3807 LHTF = .FALSE. 3808 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 3809 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 3810 CALL CCSD_IKJL(WORK(KQQ4M),NRHFT,WXRK(KS2AM)) 3811 WRITE(LUPRI,'(1X,A)') 3812 * 'Computation of (ik|f12**2|jl) integrals done' 3813 WRITE(TIMLAB,1013) '#4 ',MBSMAX 3814 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 3815 CALL FLSHFO(LUPRI) 3816C 3817C COMPUTE (IA|(1/R12)*F12|JB) INTEGRALS 3818C 3819 CALL TIMER('START ',TIMSTR,TIMEND) 3820 R12TRA = .TRUE. 3821 R12INT = .FALSE. 3822 R12SQR = .FALSE. 3823 U12INT = .FALSE. 3824 U21INT = .FALSE. 3825 V12INT = .FALSE. 3826 CCR12SM = .TRUE. 3827 MBSMAX = 4 3828 INTGAC = 2 3829 CALL DZERO(WORK(KT1AM),NT1AMX) 3830 LHTF = CCR12 3831 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 3832 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 3833 CALL CCSD_IKJL(WORK(KQQ2M),NRHFT,WXRK(KS2AM)) 3834 WRITE(LUPRI,'(1X,A)') 3835 * 'Computation of (ik|f12/r12|jl) integrals done' 3836 WRITE(TIMLAB,1013) '#2 ',MBSMAX 3837 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 3838 CALL FLSHFO(LUPRI) 3839 R12EIN = .FALSE. 3840 CCR12SM = .FALSE. 3841 LUSIFC = -1 3842 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 3843 & .FALSE.) 3844 REWIND LUSIFC 3845 CALL MOLLAB(LABEL,LUSIFC,LUPRI) 3846 READ (LUSIFC) 3847 READ (LUSIFC) 3848 WRITE (LUSIFC) (WORK(NLAMDS+KCMOM+I-1),I=1,NLAMDS) 3849 CALL GPCLOSE(LUSIFC,'KEEP') 3850 ELSE 3851 CALL TIMER('START ',TIMSTR,TIMEND) 3852 END IF 3853 CALL FLSHFO(LUPRI) 3854C 3855C COMPUTE (IA|R12|JB) INTEGRALS 3856C 3857 R12TRA = .TRUE. 3858 IF (R12EOR) THEN 3859 R12EIN = .TRUE. 3860 R12INT = .FALSE. 3861 INTGAC = 3 3862 ELSE 3863 R12INT = .TRUE. 3864 END IF 3865 R12SQR = .FALSE. 3866 U12INT = .FALSE. 3867 U21INT = .FALSE. 3868 V12INT = .FALSE. 3869 IF (R12OLD) THEN 3870 MBSMAX = 4 3871 ELSE IF ((R12NOB.AND.(.NOT.CCR12RSP)) .OR. 3872 & (NORXR.AND.(.NOT.CCR12RSP))) THEN 3873 MBSMAX = 5 3874 ELSE 3875 MBSMAX = 6 3876 END IF 3877 CALL DZERO(WORK(KT1AM),NT1AMX) 3878C........LHTF must be variable (WK/14-06-2004). 3879 LHTF = CCR12.AND.(CCR12RSP.OR.IANR12.EQ.2.OR.IANR12.EQ.3) 3880 write(lupri,*) 'ccsd_r12.11' 3881 call flshfo(lupri) ! asm 3882 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM),LHTF, 3883 & CCR12RSP,MKVAJKL,WORK(KEND1),LWRK1) 3884 write(lupri,*) 'ccsd_r12.12' 3885 call flshfo(lupri) ! asm 3886c------------------------------------------------------ 3887chf get r12 integrals on file 3888c------------------------------------------------------ 3889 CALL GPOPEN(LU43,FR12R12, 3890 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 3891 WRITE(LU43) (WXRK(KS2AM+I-1), I = 1,NH2AMX) 3892 CALL GPCLOSE(LU43,'KEEP') 3893 IF (R12EOR) THEN 3894 IF (ONEAUX) THEN 3895 WRITE(LUPRI,'(1X,A)') 3896 * 'Computation of (ip''|f12|jq) integrals done' 3897 ELSE 3898 WRITE(LUPRI,'(1X,A)') 3899 * 'Computation of (ip''|f12|jq'') integrals done' 3900 END IF 3901 WRITE(TIMLAB,1013) '#3 ',MBSMAX 3902 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 3903 ELSE 3904 IF (ONEAUX) THEN 3905 WRITE(LUPRI,'(1X,A)') 3906 * 'Computation of (ip''|r12|jq) integrals done' 3907 ELSE 3908 WRITE(LUPRI,'(1X,A)') 3909 * 'Computation of (ip''|r12|jq'') integrals done' 3910 END IF 3911 WRITE(TIMLAB,1013) 'R ',MBSMAX 3912 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 3913 END IF 3914 CALL FLSHFO(LUPRI) 3915C 3916 IF (CC2R12INT .OR. CCSDR12INT .OR. R12PRP) THEN 3917 IF (IANR12.EQ.1.OR.IANR12.EQ.3) THEN 3918 PROJVIR = (IANR12.EQ.3) 3919 CALL R12BATF(WXRK(KS2AM),FNBACK,DUMMY,.FALSE., 3920 & PROJVIR,WORK(KEND1),LWRK1) 3921 END IF 3922 IF (IANR12.EQ.2.OR.IANR12.EQ.3) THEN 3923 CALL CC_R12BATF2(WORK(KEND1),LWRK1) 3924 END IF 3925 WRITE(LUPRI,'(1X,A)') 3926 * 'Backtransformation of (ip''|r12|jq'') integrals done' 3927 CALL TIMER('R12BATF',TIMSTR,TIMEND) 3928CCN 3929 CALL CC_R12MKSAK(WORK(KEND1),LWRK1) 3930 3931 IF (CCR12RSP) THEN 3932 CALL CC_R12VXINT2(WXRK(KS2AM),WORK(KEND1),LWRK1) 3933 END IF 3934CCN 3935 END IF 3936celena 3937 IF (IANR12.EQ.1.AND.R12PRP) THEN 3938 MKVAJKL = .TRUE. 3939 FNVAJKL = 'CCR12XAJKL' 3940 CALL DZERO(WORK(KT1AM),NT1AMX) 3941 LHTF = .FALSE. 3942 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 3943 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 3944 MKVAJKL = .FALSE. 3945 ENDIF 3946celena 3947C 3948C COMPUTE (IA|1/R12|JB) INTEGRALS 3949C 3950 R12TRA = .TRUE. 3951 R12EIN = .FALSE. 3952 R12INT = .FALSE. 3953 R12SQR = .FALSE. 3954 U12INT = .FALSE. 3955 U21INT = .FALSE. 3956 V12INT = .TRUE. 3957 IF (PROFIL .OR. R12ECO) THEN 3958 MBSMAX = 6 3959 ELSE IF (R12OLD) THEN 3960 MBSMAX = 4 3961 ELSE 3962 MBSMAX = 5 3963 END IF 3964 CALL DZERO(WORK(KT1AM),NT1AMX) 3965 IF (CC2R12INT .OR. (R12PRP.AND.IANR12.EQ.1)) THEN 3966 MKVAJKL = .TRUE. 3967 IF (R12PRP.AND. IANR12.EQ.1) FNVAJKL = 'CCR12VAJKL' 3968 END IF 3969C........LHTF must be variable (WK/14-06-2004). 3970 LHTF = CCR12.AND.(IANR12.EQ.2.OR.IANR12.EQ.3) 3971 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM),LHTF, 3972 & CCR12RSP,MKVAJKL,WORK(KEND1),LWRK1) 3973 MKVAJKL = .FALSE. 3974 CALL GPOPEN(LU43,FR12V12, 3975 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 3976 WRITE(LU43) (WXRK(KS2AM+I-1), I = 1,NH2AMX) 3977 CALL GPCLOSE(LU43,'KEEP') 3978C 3979 IF (R12PRP.AND. IANR12.EQ.1) THEN 3980 CALL R12BATF(WXRK(KS2AM),FV12BACK,DUMMY,.FALSE.,.FALSE., 3981 & WORK(KEND1),LWRK1) 3982 END IF 3983C 3984 IF (R12ECO) GOTO 1094 3985 IF (FROIMP) THEN 3986 call dzero(WORK(KFRIN),NFROVR(1)) 3987 call dzero(WORK(KFRIN1),NFROVF(1)) 3988 IF (IANR12.EQ.1.AND.R12PRP) THEN 3989 LUNIT = -1 3990 CALL GPOPEN(LUNIT,'INCJDA', 3991 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 3992 REWIND(LUNIT) 3993 READ(LUNIT) (WORK(KFRIN+I-1), I = 1,NFROVR(1)) 3994 CALL GPCLOSE(LUNIT,'DELETE') 3995 3996 LUNIT = -1 3997 CALL GPOPEN(LUNIT,'INCJDI', 3998 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 3999 REWIND(LUNIT) 4000 READ(LUNIT) (WORK(KFRIN1+I-1), I = 1,NFROVF(1)) 4001 CALL GPCLOSE(LUNIT,'DELETE') 4002 ELSE 4003 CALL GPOPEN(LU43,'INCJDK', 4004 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 4005 REWIND(LU43) 4006 READ(LU43) (WORK(KFRIN+I-1), I = 1,NF2FRO(1)) 4007 CALL GPCLOSE(LU43,'KEEP') 4008 4009 ENDIF 4010 4011 END IF 4012 CALL MAKEGX(WXRK(KS2AM),WORK(KGXAM),WORK(KFRIN),WORK(KEND1), 4013 & LWRK1,2) 4014celena 4015 CALL MAKEUG(WORK(KGXAM),WORK(KUGAM), 4016 & WORK(KQQAM),WORK(KQVAM),WORK(KEND1),LWRK1,2) 4017 IF (R12HYB) THEN 4018 CALL MAKEGX(WXRK(KS2AM),WORK(KGXAM),WORK(KFRIN),WORK(KEND1), 4019 & LWRK1,1) 4020 CALL MAKEUG(WORK(KGXAM),WORK(KUGAM), 4021 & WORK(KQQAM),WORK(KQVAM),WORK(KEND1),LWRK1,1) 4022celena 4023 END IF 4024 4025 IF (IANR12.EQ.1.AND.R12PRP) CALL CC_R12MKSFR(WXRK(KS2AM), 4026 & WORK(KFRIN1),WORK(KEND1),LWRK1) 4027 1094 CONTINUE 4028 CALL GPOPEN(LU43,FR12R12, 4029 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 4030 READ(LU43) (WXRK(KS2AM+I-1), I = 1,NH2AMX) 4031 CALL GPCLOSE(LU43,'KEEP') 4032 IF (ONEAUX) THEN 4033 WRITE(LUPRI,'(1X,A)') 4034 * 'Computation of (ip''|1/r12|jq) integrals done' 4035 ELSE 4036 WRITE(LUPRI,'(1X,A)') 4037 * 'Computation of (ip''|1/r12|jq'') integrals done' 4038 END IF 4039 IF (R12EOR) THEN 4040 WRITE(TIMLAB,1013) '#1 ',MBSMAX 4041 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4042 ELSE 4043 WRITE(TIMLAB,1013) '1/R',MBSMAX 4044 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4045 END IF 4046 IF (PROFIL .OR. R12ECO) THEN 4047 GOTO 1096 4048 ELSE IF (R12OLD .OR. R12NOB .OR. NORXR) THEN 4049 GOTO 1095 4050 ELSE 4051 CALL RXR(WXRK(KS2AM),WXRK(KT2AM), 4052 * WORK(KSF),WORK(KTF), 4053 * WORK(KFOCKQ),NRHFT,NSPAIR) 4054 WRITE(LUPRI,'(1X,A)') 4055 * 'Computation of (ip''|r12|jq'')X(p''r'')'// 4056 * '(kr''|r12|lq'') integrals done' 4057 CALL TIMER('RXR ',TIMSTR,TIMEND) 4058 END IF 4059C 4060 1095 CONTINUE 4061C 4062 IF (IANR12.EQ.1.AND.R12PRP) THEN 4063 MKVAJKL = .TRUE. 4064chf FNBACK = 'CCV12_BACK' 4065 FNVAJKL = 'CCR12VIJAL' 4066 R12TRA = .TRUE. 4067 V12INT = .FALSE. 4068 IF (R12EOR) THEN 4069 R12EIN = .TRUE. 4070 R12INT = .FALSE. 4071 INTGAC = 3 4072 ELSE 4073 R12INT = .TRUE. 4074 END IF 4075 CALL DZERO(WORK(KT1AM),NT1AMX) 4076 LHTF = .FALSE. 4077 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 4078 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4079 4080 MKVAJKL = .FALSE. 4081 END IF 4082C 4083C COMPUTE (I'A|R12|JB) INTEGRALS 4084C 4085 COMBSS = .NOT. R12HYB 4086 TEMP_ONEAUX = ONEAUX 4087 R12TRA = .TRUE. 4088 IF (R12EOR) THEN 4089 ONEAUX = .FALSE. 4090 R12EIN = .TRUE. 4091 R12INT = .FALSE. 4092 ELSE 4093 R12INT = .TRUE. 4094 END IF 4095 R12SQR = .TRUE. 4096 U12INT = .FALSE. 4097 U21INT = .FALSE. 4098 V12INT = .FALSE. 4099 IF (R12EOR) THEN 4100 R12EIN = .TRUE. 4101 IF (.NOT. R12NOB) THEN 4102 LUSIFC = -1 4103 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 4104 & IDUMMY,.FALSE.) 4105 REWIND LUSIFC 4106 CALL MOLLAB(LABEL,LUSIFC,LUPRI) 4107 READ (LUSIFC) 4108 READ (LUSIFC) 4109 READ (LUSIFC) (WORK(KCMOM+I-1),I=1,NLAMDS) 4110 CALL DCOPY(NLAMDS,WORK(KCMOM),1,WORK(KCMOM+NLAMDS),1) 4111 KCMO1 = KCMOM 4112 DO ISYM = 1, NSYM 4113 KCMO2 = KCMO1 + NRHFS(ISYM)*NBAS(ISYM) 4114 LENMO = NBAS(ISYM)*NRHFS(ISYM) 4115 CALL DCOPY(LENMO,WORK(KCMO1),1,WORK(KCMO2),1) 4116 KCMO1 = KCMO1 + 4117 & (NRHFS(ISYM)+NORB1(ISYM)+NORB2(ISYM))*NBAS(ISYM) 4118 ENDDO 4119 REWIND LUSIFC 4120 CALL MOLLAB(LABEL,LUSIFC,LUPRI) 4121 READ (LUSIFC) 4122 READ (LUSIFC) 4123 WRITE (LUSIFC) (WORK(KCMOM+I-1),I=1,NLAMDS) 4124 CALL GPCLOSE(LUSIFC,'KEEP') 4125 INTGAC = 4 4126 IF (R12OLD .OR. R12HYB) THEN 4127 MBSMAX = 4 4128 ELSE 4129 MBSMAX = 5 4130 END IF 4131 CALL DZERO(WORK(KT1AM),NT1AMX) 4132 LHTF = .FALSE. 4133 CALL CCSD_IAJB(WXRK(KU2AM),WORK(KT1AM), 4134 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4135 CALL CCSD_IKJL(WORK(KQQAM),NRHFT,WXRK(KU2AM)) 4136 LUSIFC = -1 4137 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 4138 & IDUMMY,.FALSE.) 4139 REWIND LUSIFC 4140 CALL MOLLAB(LABEL,LUSIFC,LUPRI) 4141 READ (LUSIFC) 4142 READ (LUSIFC) 4143 WRITE (LUSIFC) (WORK(NLAMDS+KCMOM+I-1),I=1,NLAMDS) 4144 CALL GPCLOSE(LUSIFC,'KEEP') 4145 WRITE(TIMLAB,1013) '#4 ',MBSMAX 4146 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4147 END IF 4148 INTGAC = 3 4149 END IF 4150 IF (R12OLD) THEN 4151 IF (R12NOB .AND. .NOT. NATVIR) GOTO 1097 4152 MBSMAX = 4 4153 ELSE IF (R12NOB .OR. R12HYB) THEN 4154 MBSMAX = 5 4155 ELSE 4156 MBSMAX = 6 4157 END IF 4158 ONEAUX = TEMP_ONEAUX 4159celena 4160 IF (R12PRP) ONEAUX = .FALSE. 4161celena 4162 IF (IANR12.EQ.1.AND.R12PRP) THEN 4163 MBSMAX = 5 4164 MKVAJKL = .TRUE. 4165 FNVAJKL= 'CCR12QAJKL' 4166 R12INT = .TRUE. 4167 R12SQR = .FALSE. 4168 U12INT = .FALSE. 4169 U21INT = .FALSE. 4170 V12INT = .FALSE. 4171 CALL DZERO(WORK(KT1AM),NT1AMX) 4172 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 4173 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4174 R12SQR = .TRUE. 4175 MKVAJKL = .TRUE. 4176 FNVAJKL = 'CCR12QIJAL' 4177 LHTF = .FALSE. 4178 CALL DZERO(WORK(KT1AM),NT1AMX) 4179 CALL CCSD_IAJB(WXRK(KU2AM),WORK(KT1AM), 4180 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4181 MBSMAX = 5 4182 CALL R12BATF(WXRK(KU2AM),FQ12BACK,DUMMY,.FALSE.,.FALSE., 4183 & WORK(KEND1),LWRK1) 4184 FNVAJKL = 'CCR12UIJAL' 4185 R12SQR = .FALSE. 4186 CALL DZERO(WORK(KT1AM),NT1AMX) 4187 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 4188 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4189 CALL R12BATF(WXRK(KU2AM),FU12BACK,DUMMY,.FALSE.,.FALSE., 4190 & WORK(KEND1),LWRK1) 4191 FNVAJKL = 'CCR12UAJKL' 4192 CALL DZERO(WORK(KT1AM),NT1AMX) 4193 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 4194 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4195 R12SQR = .TRUE. 4196 MKVAJKL = .FALSE. 4197 ELSE 4198 CALL DZERO(WORK(KT1AM),NT1AMX) 4199 LHTF = .FALSE. 4200 CALL CCSD_IAJB(WXRK(KU2AM),WORK(KT1AM), 4201 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4202 ENDIF 4203 CALL QCAL(WORK(KQQAM),WXRK(KS2AM),WXRK(KU2AM), 4204 * WORK(KEND1),LWRK1,1,1) 4205 IF (.NOT. R12OLD) THEN 4206 CALL QCAL(WORK(KQQAM),WXRK(KS2AM),WXRK(KU2AM), 4207 * WORK(KEND1),LWRK1,2,2) 4208 CALL QCAL(WORK(KQQAM),WXRK(KS2AM),WXRK(KU2AM), 4209 * WORK(KEND1),LWRK1,2,3) 4210 IF (R12HYB) THEN 4211 IF (R12EOR) THEN 4212 IF (ONEAUX) THEN 4213 WRITE(LUPRI,'(1X,A)') 4214 * 'Computation of (i* p''|f12|jq) integrals done' 4215 ELSE 4216 WRITE(LUPRI,'(1X,A)') 4217 * 'Computation of (i* p''|f12|jq'') integrals done' 4218 END IF 4219 WRITE(TIMLAB,1013) '#3 ',MBSMAX 4220 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4221 ELSE 4222 IF (ONEAUX) THEN 4223 WRITE(LUPRI,'(1X,A)') 4224 * 'Computation of (i* p''|r12|jq) integrals done' 4225 ELSE 4226 WRITE(LUPRI,'(1X,A)') 4227 * 'Computation of (i* p''|r12|jq'') integrals done' 4228 END IF 4229 WRITE(TIMLAB,1013) 'R ',MBSMAX 4230 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4231 END IF 4232 COMBSS = .TRUE. 4233 CALL DZERO(WORK(KT1AM),NT1AMX) 4234 LHTF = .FALSE. 4235 CALL CCSD_IAJB(WXRK(KU2AM),WORK(KT1AM), 4236 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4237 END IF 4238 CALL COMKR(WXRK(KS2AM),WXRK(KU2AM),WORK(KV2AM), 4239 * WORK(KW2AM),WORK(KEND1),LWRK1) 4240 END IF 4241 IF (R12HYB) THEN 4242 NLBINT = 3 4243 LABINT(1:NLBINT) = 'i''p' 4244 ELSE 4245 NLBINT = 4 4246 LABINT(1:NLBINT) = 'i''p''' 4247 END IF 4248 IF (R12EOR) THEN 4249 IF (ONEAUX) THEN 4250 WRITE(LUPRI,'(1X,A)') 4251 * 'Computation of ('//LABINT(1:NLBINT)// 4252 * '|f12|jq) integrals done' 4253 ELSE 4254 WRITE(LUPRI,'(1X,A)') 4255 * 'Computation of ('//LABINT(1:NLBINT)// 4256 * '|f12|jq'') integrals done' 4257 END IF 4258 WRITE(TIMLAB,1013) '#3 ',MBSMAX 4259 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4260 ELSE 4261 IF (ONEAUX) THEN 4262 WRITE(LUPRI,'(1X,A)') 4263 * 'Computation of ('//LABINT(1:NLBINT)// 4264 * '|r12|jq) integrals done' 4265 ELSE 4266 WRITE(LUPRI,'(1X,A)') 4267 * 'Computation of ('//LABINT(1:NLBINT)// 4268 * '|r12|jq'') integrals done' 4269 END IF 4270 WRITE(TIMLAB,1013) 'R ',MBSMAX 4271 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4272 END IF 4273 CALL FLSHFO(LUPRI) 4274C 4275C COMPUTE (IA|[T1,R12]|JB) INTEGRALS ASM 4276C 4277 1097 CONTINUE 4278 R12TRA = .TRUE. 4279 IF (R12EOR) THEN 4280 R12EIN = .TRUE. 4281 INTGAC = 5 4282 END IF 4283 R12SQR = .FALSE. 4284 R12INT = .FALSE. 4285 U12INT = .TRUE. 4286 U21INT = .FALSE. 4287 V12INT = .FALSE. 4288 IF (R12OLD) THEN 4289 MBSMAX = 4 4290 ELSE 4291 MBSMAX = 5 4292 END IF 4293C 4294 TEMP_HERDIR = HERDIR 4295 TEMP_ONEAUX = ONEAUX 4296 IF (R12EOR) ONEAUX = .FALSE. 4297 IF (R12PRP .OR. ONEAUX) THEN 4298 HERDIR = .TRUE. 4299 U21INT = .TRUE. 4300 END IF 4301 IF (IANR12.EQ.1.AND.R12PRP) THEN 4302 MKVAJKL = .TRUE. 4303chf FNBACK = 'CCR12_BACK' 4304 FNVAJKL = 'CCR12BAJKL' 4305 END IF 4306 CALL DZERO(WORK(KT1AM),NT1AMX) 4307 IF (ONEAUX) THEN 4308 LHTF = .FALSE. 4309 CALL CCSD_IAJB(WXRK(KS2AM),WORK(KT1AM), 4310 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4311 ELSE 4312 LHTF = .FALSE. 4313 CALL CCSD_IAJB(WXRK(KU2AM),WORK(KT1AM), 4314 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4315 CALL TESTUU(WXRK(KU2AM),WXRK(KS2AM)) 4316 END IF 4317 HERDIR = TEMP_HERDIR 4318 ONEAUX = TEMP_ONEAUX 4319 MKVAJKL = .FALSE. 4320 U21INT = .FALSE. 4321C 4322 IF (IANR12.EQ.1.AND.R12PRP) THEN 4323chf FNBACK = 'CCT12_BACK' 4324 CALL R12BATF(WXRK(KS2AM),FT12BACK,DUMMY,.FALSE.,.FALSE., 4325 & WORK(KEND1),LWRK1) 4326 R12TRA = .TRUE. 4327 IF (R12EOR) THEN 4328 R12EIN = .TRUE. 4329 R12INT = .FALSE. 4330 INTGAC = 3 4331 ELSE 4332 R12INT = .TRUE. 4333 END IF 4334 R12SQR = .FALSE. 4335 U12INT = .FALSE. 4336 U21INT = .FALSE. 4337 V12INT = .FALSE. 4338 IF (R12OLD) THEN 4339 MBSMAX = 4 4340 ELSE 4341 MBSMAX = 5 4342 END IF 4343 MKVAJKL = .TRUE. 4344chf FNBACK = 'CCT12_BACK' 4345 FNVAJKL = 'CCR12BIJAL' 4346 CALL DZERO(WORK(KT1AM),NT1AMX) 4347 LHTF = .FALSE. 4348 CALL CCSD_IAJB(WXRK(KT2AM),WORK(KT1AM), 4349 & LHTF,.FALSE.,MKVAJKL,WORK(KEND1),LWRK1) 4350 MKVAJKL = .FALSE. 4351 END IF 4352C 4353 IF (R12EOR) THEN 4354 IF (ONEAUX) THEN 4355 WRITE(LUPRI,'(1X,A)') 4356 * 'Computation of (ip''|[T,f12]|jq) integrals done' 4357 ELSE 4358 WRITE(LUPRI,'(1X,A)') 4359 * 'Computation of (ip''|[T,f12]|jq'') integrals done' 4360 END IF 4361 WRITE(TIMLAB,1013) '#5 ',MBSMAX 4362 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4363 ELSE 4364 IF (ONEAUX) THEN 4365 WRITE(LUPRI,'(1X,A)') 4366 * 'Computation of (ip''|[T,r12]|jq) integrals done' 4367 ELSE 4368 WRITE(LUPRI,'(1X,A)') 4369 * 'Computation of (ip''|[T,r12]|jq'') integrals done' 4370 END IF 4371 WRITE(TIMLAB,1013) 'T,R',MBSMAX 4372 CALL TIMER(TIMLAB,TIMSTR,TIMEND) 4373 END IF 4374C 4375 1096 CONTINUE 4376 4377 IF (CC2R12INT.AND.(IANR12.EQ.2.OR.IANR12.EQ.3)) THEN 4378 ! Resort integrals on WORK(KS2AM) to the order needed in 4379 ! CC-R12, result returned on WXRK(KT2AM) 4380 CALL CC_R12GETKIAJB(WORK(KS2AM),WXRK(KT2AM),ONEAUX, 4381 & NELEM,NH2AMX,IH2AM,NH1AM,IH1AM,NORB1,NRHF,NRHFSA) 4382 4383 CALL GPOPEN(LU43,FTR12, 4384 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 4385 WRITE(LU43) (WXRK(KT2AM+I-1), I = 1,NELEM) 4386 CALL GPCLOSE(LU43,'KEEP') 4387 4388 ! Resort integrals on WORK(KV2AM) to the order needed in 4389 ! CC-R12, result returned on WXRK(KT2AM) 4390 CALL CC_R12GETKIAJB(WORK(KV2AM),WXRK(KT2AM),ONEAUX, 4391 & NELEM,NH2AMX,IH2AM,NH1AM,IH1AM,NORB1,NRHF,NRHFSA) 4392 CALL GPOPEN(LU43,FKR12,'UNKNOWN','SEQUENTIAL', 4393 & 'UNFORMATTED',IDUM,LDUM) 4394 WRITE(LU43) (WXRK(KT2AM+I-1), I = 1,NELEM) 4395 CALL GPCLOSE(LU43,'KEEP') 4396 END IF 4397C 4398 IF (CCSDR12INT.AND.(IANR12.EQ.2.OR.IANR12.EQ.3)) THEN 4399 !write integrals r_(bp')^(mn) to file 4400 CALL CC_R12GETRAMNP(WORK(KS2AM),'R12RMNAP',ONEAUX,NORB1, 4401 & NORB2,NRHF,NRHFSA,IH1AM,IH2AM, 4402 & WXRK(KT2AM),LENT2AM) 4403 END IF 4404C 4405 CALL GPOPEN(LU43,FR12V12, 4406 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 4407 READ(LU43) (WXRK(KT2AM+I-1), I = 1,NH2AMX) 4408 CALL GPCLOSE(LU43,'KEEP') 4409 CALL GPOPEN(LU43,FR12R12, 4410 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 4411 READ(LU43) (WXRK(KR2AM+I-1), I = 1,NH2AMX) 4412 CALL GPCLOSE(LU43,'KEEP') 4413 CALL GPOPEN(LUSIRG,'SIRIFC', 4414 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 4415 REWIND LUSIRG 4416 CALL MOLLAB('SIR IPH ',LUSIRG,LUPRI) 4417 READ (LUSIRG) POTNUC,EMY,EACTIV,EMCSCF,ISTATE,ISPIN, 4418 * NACTEL,LSYM,MS2 4419 ESCF = EMCSCF 4420 CALL MOLLAB(LABEL,LUSIRG,LUPRI) 4421 READ (LUSIRG) 4422 READ (LUSIRG) (WORK(KFOCKD+I-1), I=1,NORBTS) 4423 CALL GPCLOSE(LUSIRG,'KEEP') 4424 IF (FROIMP .OR. FROEXP) 4425 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK1) 4426 CALL FOCK_REORDER(WORK(KFOCKD),WORK(KEND1),LWRK1) 4427C 4428C COMPUTE MP2-R12 ENERGY 4429C 4430 IF (PROFIL) THEN 4431 CALL PRFDRV(WXRK(KT2AM),WXRK(KR2AM),WORK(KFOCKD), 4432 * WORK(KEND1),LWRK1) 4433 CALL QUIT('PROFILES ALL DONE') 4434 ELSE 4435 CALL R12DRV(WXRK(KT2AM),WXRK(KR2AM),WXRK(KS2AM), 4436 * WORK(KV2AM),WORK(KW2AM), 4437 * WORK(KFOCKD),WORK(KSF),WORK(KTF), 4438 * WORK(KEND1),LWRK1, 4439 * WORK(KQQ2M),WORK(KQQ4M),WORK(KQQ6M)) 4440 END IF 4441 CALL TIMER('E(R12)',TIMSTR,TIMEND) 4442 4443 IF (CC2R12INT.AND.(IANR12.EQ.2.OR.IANR12.EQ.3)) THEN 4444 ! Resort R12 integrals into WORK, sort them to the order 4445 ! needed in CC-R12, result returned on WXRK, and write them 4446 ! back to file 4447 CALL GPOPEN(LU43,FR12R12, 4448 & 'UNKNOWN',' ','UNFORMATTED',IDUM,LDUM) 4449 READ(LU43) (WORK(I), I = 1,NH2AMX) 4450 CALL CC_R12GETKIAJB(WORK,WXRK(KR2AM),ONEAUX,NELEM, 4451 & NH2AMX,IH2AM,NH1AM,IH1AM,NORB1,NRHF,NRHFSA) 4452 REWIND(LU43) 4453 WRITE(LU43) (WXRK(KR2AM + I - 1), I = 1,NELEM) 4454 CALL GPCLOSE(LU43,'KEEP') 4455 END IF 4456 ENDIF 4457 R12TRA = .FALSE. 4458 R12INT = .FALSE. 4459 R12SQR = .FALSE. 4460 U12INT = .FALSE. 4461 U21INT = .FALSE. 4462 R12EIN = .FALSE. 4463 V12INT = .TRUE. 4464 INTGAC = 0 4465 MBSMAX = 4 4466 ONEAUX = .FALSE. 4467 DIRECT = TEMP_DIRECT 4468 CALL QEXIT('CCSD_R12') 4469 RETURN 4470 1013 FORMAT(A3,'[',I1,']') 4471 END 4472C /* Deck ccsd_ikjl */ 4473 SUBROUTINE CCSD_IKJL(QQ,NOCDIM,XX) 4474C 4475C This subroutine collects the (ik|jl) MO-integrals. 4476C 4477C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 4478C 4479#include "implicit.h" 4480#include "priunit.h" 4481 DIMENSION XX(*),QQ(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 4482#include "r12int.h" 4483#include "ccorb.h" 4484#include "ccsdsym.h" 4485#include "ccsdinp.h" 4486 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4487 CALL DZERO(QQ,NOCDIM**4) 4488 DO ISYMIK = 1, NSYM 4489 ISYMJL = ISYMIK 4490 DO ISYMI = 1, NSYM 4491 ISYMK = MULD2H(ISYMI,ISYMIK) 4492 DO ISYMJ = 1, NSYM 4493 ISYML = MULD2H(ISYMJ,ISYMJL) 4494 DO I = 1, NRHF(ISYMI) 4495 KOFFI = IRHF(ISYMI) + I 4496 DO K = NRHFFR(ISYMK) + 1, NRHFS(ISYMK) 4497 KOFFK = IRHF(ISYMK) + K - NRHFFR(ISYMK) 4498 DO J = 1, NRHF(ISYMJ) 4499 KOFFJ = IRHF(ISYMJ) + J 4500 DO L = NRHFFR(ISYML) + 1, NRHFS(ISYML) 4501 KOFFL = IRHF(ISYML) + L - NRHFFR(ISYML) 4502 NKI = IT1AM(ISYMK,ISYMI) + 4503 * NVIR(ISYMK)*(I-1) + K 4504 NLJ = IT1AM(ISYML,ISYMJ) + 4505 * NVIR(ISYML)*(J-1) + L 4506 IF (ONEAUX) THEN 4507 IF (R12SQR) THEN 4508 NKILJ = IU2AM(ISYMIK,ISYMJL) 4509 * + NT1AM(ISYMJL)*(NKI - 1) 4510 * + NLJ 4511 ELSE 4512 NKILJ = IH2AM(ISYMIK,ISYMJL) 4513 * + INDEX(NKI,NLJ) 4514 END IF 4515 ELSE 4516 IF (R12SQR) THEN 4517 NKILJ = IU2AM(ISYMIK,ISYMJL) 4518 * + NT1AM(ISYMJL)*(NKI - 1) 4519 * + NLJ 4520 ELSE 4521 NKILJ = IT2AM(ISYMIK,ISYMJL) 4522 * + INDEX(NKI,NLJ) 4523 END IF 4524 END IF 4525 QQ(KOFFI,KOFFJ,KOFFK,KOFFL)=XX(NKILJ) 4526 ENDDO 4527 ENDDO 4528 ENDDO 4529 ENDDO 4530 ENDDO 4531 ENDDO 4532 ENDDO 4533 IF (IPRINT .GE. 5) THEN 4534 NP = NOCDIM**2 4535 CALL AROUND('Output from CCSD_IKJL') 4536 CALL OUTPUT(QQ,1,NP,1,NP,NP,NP,1,LUPRI) 4537 END IF 4538 RETURN 4539 END 4540C /* Deck vclean */ 4541 SUBROUTINE VCLEAN(A,N,THR) 4542C 4543C This subroutine removes small matrix elements. 4544C 4545C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 4546C 4547#include "implicit.h" 4548#include "priunit.h" 4549 PARAMETER (D0 = 0D0) 4550 DIMENSION A(N,N) 4551 DO 100 I = 1, N 4552 DO 100 J = 1, N 4553 IF (ABS(A(I,J)) .LT. THR) A(I,J) = 0D0 4554 100 CONTINUE 4555 RETURN 4556 END 4557C /* Deck r12mp2 */ 4558 SUBROUTINE R12MP2(VS,VT,VOS,VOT,SS,ST,EKL,NOCDIM,NSPAIR,FS,FT, 4559 * BB,UU,WW,QQ,SF,TF,EV,XF,VV,IJPAIR,WS,WT,DELTA,CCS,CCT) 4560C 4561C This subroutine evaluates the MP2-R12 energy by 4562C solving sets of linear equations using 4563C singular value decomposition. 4564C 4565C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 4566C 4567#include "implicit.h" 4568#include "priunit.h" 4569#include "r12int.h" 4570 PARAMETER (D0 = 0D0, DP5 = 5D-1, D2 = 2D0, D99 = 1D99) 4571 DIMENSION VS(NSPAIR), VT(NSPAIR), EV(NSPAIR) 4572 DIMENSION SS(NSPAIR,NSPAIR), ST(NSPAIR,NSPAIR) 4573 DIMENSION SF(NSPAIR,NSPAIR), TF(NSPAIR,NSPAIR) 4574 DIMENSION BB(NSPAIR,NSPAIR), QQ(NSPAIR,NSPAIR) 4575 DIMENSION VV(NSPAIR,NSPAIR), WW(NSPAIR) 4576 DIMENSION UU(NSPAIR,NSPAIR) 4577 DIMENSION CCS(NSPAIR),CCT(NSPAIR) 4578 DIMENSION VOS(NSPAIR,NSPAIR),VOT(NSPAIR,NSPAIR) 4579 LOGICAL MATV, MATU 4580 MATV = .TRUE. 4581 MATU = .TRUE. 4582 IERR = 0 4583 WS = D99 4584 WT = D99 4585C 4586C SINGLET PAIRS 4587C 4588 IF (R12ECO) THEN 4589 DO J = 1, NSPAIR 4590 DO I = 1, NSPAIR 4591 BB(I,J) = SF(I,J) - EKL * SS(I,J) 4592 END DO 4593 END DO 4594 ELSE 4595 DO J = 1, NSPAIR 4596 DO I = 1, NSPAIR 4597 BB(I,J) = - SF(I,J) - SF(J,I) 4598 * - XF * (D2 * EKL - EV(I) - EV(J)) 4599 * * (SS(I,J) + SS(J,I)) 4600 BB(I,J) = BB(I,J) * DP5 4601 END DO 4602 BB(J,J) = BB(J,J) + DELTA 4603 END DO 4604 END IF 4605 CALL DZERO(WW,NSPAIR) 4606 CALL DZERO(UU,NSPAIR*NSPAIR) 4607 CALL DZERO(VV,NSPAIR*NSPAIR) 4608 CALL DZERO(QQ,NSPAIR*NSPAIR) 4609 CALL SVD(NSPAIR,NSPAIR,NSPAIR,BB,WW,MATU,UU,MATV,VV,IERR,QQ) 4610 DO I = 1, NSPAIR 4611 IF (ABS(WW(I)) .LT. SVDTHR) WW(I) = D0 4612 IF (WW(I) .NE. D0) THEN 4613 WS = MIN(WS, WW(I)) 4614 IF (WW(I) .LT. VCLTHR) WW(I) = D0 4615 END IF 4616 END DO 4617 IF (IERR .NE. 0) THEN 4618 WRITE(LUPRI,'(//A,2I10/)') ' SVD ERROR / IJPAIR, IERR =', 4619 * IJPAIR, IERR 4620 CALL QUIT('SVD ERROR -- SINGLET') 4621 END IF 4622 CALL SVBKSB(UU,WW,VV,NSPAIR,NSPAIR,NSPAIR,NSPAIR,VS,BB,QQ) 4623 DO I = 1, NSPAIR 4624 CCS(I)= - BB(I,1) 4625 END DO 4626 FS = -DDOT(NSPAIR,VOS,1,BB,1) 4627C 4628C TRIPLET PAIRS 4629C 4630 IF (R12ECO) THEN 4631 DO J = 1, NSPAIR 4632 DO I = 1, NSPAIR 4633 BB(I,J) = TF(I,J) - EKL * ST(I,J) 4634 END DO 4635 END DO 4636 ELSE 4637 DO J = 1, NSPAIR 4638 DO I = 1, NSPAIR 4639 BB(I,J) = - TF(I,J) - TF(J,I) 4640 * - XF * (D2 * EKL - EV(I) - EV(J)) 4641 * * (ST(I,J) + ST(J,I)) 4642 BB(I,J) = BB(I,J) * DP5 4643 END DO 4644 IF (ABS(BB(J,J)) .GT. 1D-12) BB(J,J) = BB(J,J) + DELTA 4645 END DO 4646 END IF 4647 CALL DZERO(WW,NSPAIR) 4648 CALL DZERO(UU,NSPAIR*NSPAIR) 4649 CALL DZERO(VV,NSPAIR*NSPAIR) 4650 CALL DZERO(QQ,NSPAIR*NSPAIR) 4651 CALL SVD(NSPAIR,NSPAIR,NSPAIR,BB,WW,MATU,UU,MATV,VV,IERR,QQ) 4652 DO I = 1, NSPAIR 4653 IF (ABS(WW(I)) .LT. SVDTHR) WW(I) = D0 4654 IF (WW(I) .NE. D0) THEN 4655 WT = MIN(WT, WW(I)) 4656 IF (WW(I) .LT. VCLTHR) WW(I) = D0 4657 END IF 4658 END DO 4659 IF (IERR .NE. 0) THEN 4660 WRITE(LUPRI,'(//A,2I10/)') ' SVD ERROR / IJPAIR, IERR =', 4661 * IJPAIR, IERR 4662 CALL QUIT('SVD ERROR -- TRIPLET') 4663 END IF 4664 CALL SVBKSB(UU,WW,VV,NSPAIR,NSPAIR,NSPAIR,NSPAIR,VT,BB,QQ) 4665 DO I =1,NSPAIR 4666 CCT(I) = - BB(I,1) 4667 END DO 4668 FT = -DDOT(NSPAIR,VOT,1,BB,1) 4669 RETURN 4670 END 4671C /* Deck r12fxl */ 4672 SUBROUTINE R12FXL(XS,XT,FLM,SS,ST,NIND,II,JJ,NOCDIM,NSPAIR) 4673C Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003). 4674#include "implicit.h" 4675#include "priunit.h" 4676 PARAMETER ( DP5 = 0.5D0, D1 = 1.0D0, D1M = -1.0D0, D2 = 2.0D0 ) 4677 PARAMETER ( SYMTHR = 1.0D-10 ) 4678 DIMENSION FLM(*) 4679 DIMENSION II(NSPAIR), JJ(NSPAIR), NIND(NOCDIM,NOCDIM) 4680 DIMENSION WS(NSPAIR,NSPAIR), WT(NSPAIR,NSPAIR) 4681 DIMENSION XS(NSPAIR,NSPAIR), XT(NSPAIR,NSPAIR) 4682 DIMENSION SS(NSPAIR,NSPAIR), ST(NSPAIR,NSPAIR) 4683C 4684 DSQ2 = SQRT(D2) 4685 DSQ2I = D1 / DSQ2 4686 NSP2 = NSPAIR * NSPAIR 4687 CALL DZERO(XS,NSP2) 4688 CALL DZERO(XT,NSP2) 4689 DO NPIJ = 1, NSPAIR 4690 NI = II(NPIJ) 4691 NJ = JJ(NPIJ) 4692 DO NPKL = 1, NSPAIR 4693 NK = II(NPKL) 4694 NL = JJ(NPKL) 4695 DO NO = 1, NOCDIM 4696 NIO = NIND(NI,NO) 4697 NJO = NIND(NJ,NO) 4698 NKO = NIND(NK,NO) 4699 NLO = NIND(NL,NO) 4700 IF (NO.NE.NI) THEN 4701 IF (NO.LT.NJ) THEN 4702 FAC = D1M 4703 ELSE 4704 FAC = D1 4705 END IF 4706 IF (NO.EQ.NJ) THEN 4707 FA = DSQ2 4708 ELSE 4709 FA = D1 4710 END IF 4711 IF (NI.EQ.NJ) FA = FA * DSQ2I 4712 XS(NPIJ,NPKL) = XS(NPIJ,NPKL) + FA * 4713 & SS(NJO,NPKL)*FLM(NIO) 4714 XT(NPIJ,NPKL) = XT(NPIJ,NPKL) + 4715 & FAC * ST(NJO,NPKL)*FLM(NIO) 4716 END IF 4717 IF (NO.NE.NJ) THEN 4718 IF (NO.GT.NI) THEN 4719 FAC = D1M 4720 ELSE 4721 FAC = D1 4722 END IF 4723 IF (NO.EQ.NI) THEN 4724 FA = DSQ2 4725 ELSE 4726 FA = D1 4727 END IF 4728 IF (NI.EQ.NJ) FA = FA * DSQ2I 4729 XS(NPIJ,NPKL) = XS(NPIJ,NPKL) + FA * 4730 & SS(NIO,NPKL)*FLM(NJO) 4731 XT(NPIJ,NPKL) = XT(NPIJ,NPKL) + 4732 & FAC * ST(NIO,NPKL)*FLM(NJO) 4733 END IF 4734 IF (NO.NE.NK) THEN 4735 IF (NO.LT.NL) THEN 4736 FAC = D1M 4737 ELSE 4738 FAC = D1 4739 END IF 4740 IF (NO.EQ.NL) THEN 4741 FA = DSQ2 4742 ELSE 4743 FA = D1 4744 END IF 4745 IF (NK.EQ.NL) FA = FA * DSQ2I 4746 XS(NPIJ,NPKL) = XS(NPIJ,NPKL) + FA * 4747 & SS(NPIJ,NLO)*FLM(NKO) 4748 XT(NPIJ,NPKL) = XT(NPIJ,NPKL) + 4749 & FAC * ST(NPIJ,NLO)*FLM(NKO) 4750 END IF 4751 IF (NO.NE.NL) THEN 4752 IF (NO.GT.NK) THEN 4753 FAC = D1M 4754 ELSE 4755 FAC = D1 4756 END IF 4757 IF (NO.EQ.NK) THEN 4758 FA = DSQ2 4759 ELSE 4760 FA = D1 4761 END IF 4762 IF (NK.EQ.NL) FA = FA * DSQ2I 4763 XS(NPIJ,NPKL) = XS(NPIJ,NPKL) + FA * 4764 & SS(NPIJ,NKO)*FLM(NLO) 4765 XT(NPIJ,NPKL) = XT(NPIJ,NPKL) + 4766 & FAC * ST(NPIJ,NKO)*FLM(NLO) 4767 END IF 4768 END DO 4769 XS(NPIJ,NPKL) = DP5 * XS(NPIJ,NPKL) 4770 XT(NPIJ,NPKL) = DP5 * XT(NPIJ,NPKL) 4771 END DO 4772 END DO 4773 DO NPIJ = 1, NSPAIR 4774 DO NPKL = 1, NSPAIR 4775 IF (ABS(XS(NPIJ,NPKL)-XS(NPKL,NPIJ)) .GT. SYMTHR) THEN 4776 WRITE (LUPRI,'(A,2I5,2E20.10)') 4777 & 'Singlet matrix <XF+FX> not symmetric',NPIJ,NPKL, 4778 & XS(NPIJ,NPKL),XS(NPKL,NPIJ) 4779 CALL QUIT('Singlet matrix <XF+FX> not symmetric') 4780 END IF 4781 IF (ABS(XT(NPIJ,NPKL)-XT(NPKL,NPIJ)) .GT. SYMTHR) THEN 4782 WRITE (LUPRI,'(A,2I5,2E20.10)') 4783 & 'Triplet matrix <XF+FX> not symmetric',NPIJ,NPKL, 4784 & XT(NPIJ,NPKL),XT(NPKL,NPIJ) 4785 CALL QUIT('Triplet matrix <XF+FX> not symmetric') 4786 END IF 4787 END DO 4788 END DO 4789 RETURN 4790 END 4791C /* Deck r12fcl */ 4792 SUBROUTINE R12FCL(ZS,ZT,WS,WT,VS,VT,YS,YT,SS,ST,CCS,CCT,FLM,NIND, 4793 & II,JJ,NOCDIM,NSPAIR,DELTA) 4794C Written by Claire C.M. Samson (University of Karlsruhe, 28 April 2003). 4795#include "implicit.h" 4796#include "priunit.h" 4797 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, D1M = -1.0D0 , D2 = 2D0 ) 4798 DIMENSION FLM(*) 4799 DIMENSION II(NSPAIR), JJ(NSPAIR), NIND(NOCDIM,NOCDIM) 4800 DIMENSION ZS(NSPAIR,NSPAIR), ZT(NSPAIR,NSPAIR) 4801 DIMENSION WS(NSPAIR,NSPAIR), WT(NSPAIR,NSPAIR) 4802 DIMENSION VS(NSPAIR,NSPAIR), VT(NSPAIR,NSPAIR) 4803 DIMENSION YS(NSPAIR,NSPAIR), YT(NSPAIR,NSPAIR) 4804 DIMENSION SS(NSPAIR,NSPAIR), ST(NSPAIR,NSPAIR) 4805 DIMENSION CCS(NSPAIR,NSPAIR), CCT(NSPAIR,NSPAIR) 4806C 4807 DSQ2 = SQRT(D2) 4808 DSQ2I = D1 / DSQ2 4809 NSP2 = NSPAIR * NSPAIR 4810 CALL DZERO(YS,NSP2) 4811 CALL DZERO(YT,NSP2) 4812 DO NPIJ = 1, NSPAIR 4813 NI = II(NPIJ) 4814 NJ = JJ(NPIJ) 4815 DO NPKL = 1, NSPAIR 4816 DO NO = 1, NOCDIM 4817 NIO = NIND(NI,NO) 4818 NJO = NIND(NJ,NO) 4819 IF (NO.NE.NI) THEN 4820 IF (NO.LT.NJ) THEN 4821 FAC = D1M 4822 ELSE 4823 FAC = D1 4824 END IF 4825 IF (NO.EQ.NJ) THEN 4826 FA = DSQ2 4827 ELSE 4828 FA = D1 4829 END IF 4830 YS(NPKL,NPIJ) = YS(NPKL,NPIJ) + FA * 4831 * CCS(NPKL,NJO)*FLM(NIO) 4832 YT(NPKL,NPIJ) = YT(NPKL,NPIJ) + 4833 * FAC * CCT(NPKL,NJO)*FLM(NIO) 4834 END IF 4835 IF (NO.NE.NJ) THEN 4836 IF (NO.GT.NI) THEN 4837 FAC = D1M 4838 ELSE 4839 FAC = D1 4840 END IF 4841 IF (NO.EQ.NI) THEN 4842 FA = DSQ2 4843 ELSE 4844 FA = D1 4845 END IF 4846 YS(NPKL,NPIJ) = YS(NPKL,NPIJ) + FA * 4847 * CCS(NPKL,NIO)*FLM(NJO) 4848 YT(NPKL,NPIJ) = YT(NPKL,NPIJ) + 4849 * FAC * CCT(NPKL,NIO)*FLM(NJO) 4850 END IF 4851 END DO 4852 IF (NI.EQ.NJ) YS(NPKL,NPIJ) = YS(NPKL,NPIJ) * DSQ2I 4853 END DO 4854 END DO 4855 DO NPIJ = 1, NSPAIR 4856 DO NPKL = 1, NSPAIR 4857 AAA = D0 4858 BBB = D0 4859 CCC = D0 4860 DDD = D0 4861 DO NPMN = 1,NSPAIR 4862 AAA = AAA + SS(NPIJ,NPMN) * YS(NPMN,NPKL) 4863 BBB = BBB + WS(NPIJ,NPMN) * CCS(NPMN,NPKL) 4864 CCC = CCC + ST(NPIJ,NPMN) * YT(NPMN,NPKL) 4865 DDD = DDD + WT(NPIJ,NPMN) * CCT(NPMN,NPKL) 4866 END DO 4867 ZS(NPIJ,NPKL)=VS(NPIJ,NPKL)-DELTA*CCS(NPIJ,NPKL)-AAA+BBB 4868 ZT(NPIJ,NPKL)=VT(NPIJ,NPKL)-DELTA*CCT(NPIJ,NPKL)-CCC+DDD 4869 END DO 4870 END DO 4871 RETURN 4872 END 4873C /* Deck svbksb */ 4874 SUBROUTINE SVBKSB(U,W,V,M,N,MP,NP,B,X,TMP) 4875C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 4876#include "implicit.h" 4877 PARAMETER (D0 = 0D0) 4878 DIMENSION U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP),TMP(N) 4879 DO 12 J=1,N 4880 S=D0 4881 IF (W(J).NE.D0) THEN 4882 DO 11 I=1,M 4883 S=S+U(I,J)*B(I) 4884 11 CONTINUE 4885 S=S/W(J) 4886 ENDIF 4887 TMP(J)=S 4888 12 CONTINUE 4889 DO 14 J=1,N 4890 S=0.0D0 4891 DO 13 JJ=1,N 4892 S=S+V(J,JJ)*TMP(JJ) 4893 13 CONTINUE 4894 X(J)=S 4895 14 CONTINUE 4896 RETURN 4897 END 4898C /* Deck ecodrv */ 4899 SUBROUTINE ECODRV(V2AM,R2AM,EV,V11,B11,Q11,VS11,VT11,WS11,WT11, 4900 * QS11,QT11,BS11,BT11,RS11,RT11, 4901 * NICDIM,NOCDIM,NSPAIR,NTPAIR,ES,ET,FS,FT,EVS, 4902 * CNINV,SF,CS11,CT11) 4903C 4904C This subroutine computes the V(klmn), X(klmn), and B(klmn) 4905C matrices and evaluates the MP2-R12 energies. 4906C 4907C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 4908C 4909#include "implicit.h" 4910#include "priunit.h" 4911 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 4912 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0, DELTA=0.0d0) 4913 PARAMETER (THRDIA = 1D-9) 4914 DIMENSION V2AM(*), R2AM(*), EV(*), SF(*) 4915 REAL*8 R1, R2, V1, VR, RR, BB, 4916 * V11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 4917 * B11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 4918 * Q11(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 4919 DIMENSION VS11(NSPAIR,NSPAIR), VT11(NSPAIR,NSPAIR) 4920 DIMENSION WS11(NSPAIR,NSPAIR), WT11(NSPAIR,NSPAIR) 4921 DIMENSION QS11(NSPAIR,NSPAIR), QT11(NSPAIR,NSPAIR) 4922 DIMENSION BS11(NSPAIR,NSPAIR), BT11(NSPAIR,NSPAIR) 4923 DIMENSION RS11(NSPAIR,NSPAIR), RT11(NSPAIR,NSPAIR) 4924 DIMENSION ES(NSPAIR),ET(NSPAIR),FS(NSPAIR),FT(NSPAIR) 4925 DIMENSION EVS(NSPAIR),CNINV(NSPAIR,8) 4926 DIMENSION CS11(NSPAIR,NSPAIR),CT11(NSPAIR,NSPAIR) 4927 DIMENSION ISB(8) 4928 INTEGER IDUM 4929 LOGICAL LDUM 4930#include "r12int.h" 4931#include "ccorb.h" 4932#include "ccsdsym.h" 4933#include "ccsdinp.h" 4934 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 4935 FF = D1 / SQRT(D2) 4936 ISB(1) = 0 4937 DO 100 ISYM = 2, NSYM 4938 NBASF = NORB1(ISYM-1) + NORB2(ISYM-1) 4939 NNBASF = NBASF * (NBASF + 1) / 2 4940 ISB(ISYM) = ISB(ISYM-1) + NNBASF 4941 100 CONTINUE 4942 NBASF = NORB1(NSYM) + NORB2(NSYM) 4943 NNBASF = ISB(NSYM) + NBASF * (NBASF + 1) / 2 4944 LUMULB = -34 4945 CALL GPOPEN(LUMULB,'AUXFCK','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 4946 READ (LUMULB,'(4E30.20)') (SF(I), I = 1, NNBASF) 4947 CALL GPCLOSE(LUMULB,'KEEP') 4948 DO 200 ISYM = 1, NSYM 4949 DO 210 I = NORB1(ISYM) + 2, NVIR(ISYM) 4950 DO 220 J = NORB1(ISYM) + 1 , I - 1 4951 IJ = ISB(ISYM) + INDEX(I,J) 4952 IF (ABS(SF(IJ)) .GT. THRDIA) THEN 4953 WRITE(LUPRI,'(/A/A,3I5,E20.10/)') 4954 * '@ WARNING : Fock matrix not diagonal', 4955 * ' Nondiagonal element is :', 4956 * ISYM,I,J,SF(IJ) 4957 IF (ABS(SF(IJ)) .GT. SQRT(THRDIA)) 4958 * CALL QUIT('Fock matrix not diagonal') 4959 END IF 4960 220 CONTINUE 4961 210 CONTINUE 4962 200 CONTINUE 4963C 4964 DO 205 I = 1, NOCDIM**4 4965 Q11(I,1,1,1) = D0 4966 V11(I,1,1,1) = D0 4967 B11(I,1,1,1) = D0 4968 205 CONTINUE 4969 CALL DZERO(ES,NSPAIR) 4970 CALL DZERO(ET,NSPAIR) 4971C 4972C CONSTRUCT PRODUCTS (IA|JB)*(KA|LB) 4973C 4974 E2 = D0 4975 E2S = D0 4976 E2T = D0 4977 DO 101 ISYMIA = 1, NSYM 4978 ISYMJB = ISYMIA 4979 DO 102 ISYMI = 1, NSYM 4980 ISYMA = MULD2H(ISYMI,ISYMIA) 4981 DO 103 ISYMJ = 1, NSYM 4982 ISYMB = MULD2H(ISYMJ,ISYMJB) 4983C 4984C COMPUTE MP2 ENERGY 4985C 4986 DO 201 I = 1, NRHF(ISYMI) 4987 KOFFI = IRHF(ISYMI) + I 4988 DO 202 J = 1, NRHF(ISYMJ) 4989 KOFFJ = IRHF(ISYMJ) + J 4990 DO 203 A = NRHFS(ISYMA) + 1, NORB1(ISYMA) 4991 ISYMJA = MULD2H(ISYMJ,ISYMA) 4992 KOFFA = IVIR(ISYMA) + A 4993 DO 204 B = NRHFS(ISYMB) + 1, NORB1(ISYMB) 4994 ISYMIB = MULD2H(ISYMI,ISYMB) 4995 KOFFB = IVIR(ISYMB) + B 4996 NAI = IT1AM(ISYMA,ISYMI) + 4997 * NVIR(ISYMA)*(I-1) + A 4998 NBJ = IT1AM(ISYMB,ISYMJ) + 4999 * NVIR(ISYMB)*(J-1) + B 5000 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 5001 * INDEX(NAI,NBJ) 5002 NAJ = IT1AM(ISYMA,ISYMJ) + 5003 * NVIR(ISYMA)*(J-1) + A 5004 NBI = IT1AM(ISYMB,ISYMI) + 5005 * NVIR(ISYMB)*(I-1) + B 5006 NAJBI = IT2AM(ISYMJA,ISYMIB) + 5007 * INDEX(NAJ,NBI) 5008 VAIBJ = V2AM(NAIBJ) 5009 VAJBI = V2AM(NAJBI) 5010 VV = VAIBJ * (D2 * VAIBJ - VAJBI) 5011 DENOM = D1 / (EV(KOFFI) + EV(KOFFJ) - 5012 * EV(KOFFA) - EV(KOFFB)) 5013 E2 = E2 + VV * DENOM 5014 IJ = INDEX(KOFFI,KOFFJ) 5015 VS = (VAIBJ + VAJBI)**2 5016 VT = (VAIBJ - VAJBI)**2 5017 ES(IJ) = ES(IJ) + VS * DENOM 5018 ET(IJ) = ET(IJ) + VT * DENOM 5019 204 CONTINUE 5020 203 CONTINUE 5021 202 CONTINUE 5022 201 CONTINUE 5023C 5024 DO 104 ISYMKA = 1, NSYM 5025 ISYMLB = ISYMKA 5026 ISYMK = MULD2H(ISYMA,ISYMKA) 5027 ISYML = MULD2H(ISYMB,ISYMLB) 5028 DO 105 I = 1, NRHF(ISYMI) 5029 KOFFI = IRHF(ISYMI) + I 5030 DO 106 K = 1, NRHF(ISYMK) 5031 KOFFK = IRHF(ISYMK) + K 5032 DO 107 A = NORB1(ISYMA)+1, NVIR(ISYMA) 5033 NAI = IT1AM(ISYMA,ISYMI) + 5034 * NVIR(ISYMA)*(I-1) + A 5035 NAK = IT1AM(ISYMA,ISYMK) + 5036 * NVIR(ISYMA)*(K-1) + A 5037 DO 108 J = 1, NRHF(ISYMJ) 5038 KOFFJ = IRHF(ISYMJ) + J 5039 DO 109 L = 1, NRHF(ISYML) 5040 KOFFL = IRHF(ISYML) + L 5041 DO 110 B = NORB1(ISYMB)+1, NVIR(ISYMB) 5042 NBJ = IT1AM(ISYMB,ISYMJ) + 5043 * NVIR(ISYMB)*(J-1) + B 5044 NBL = IT1AM(ISYMB,ISYML) + 5045 * NVIR(ISYMB)*(L-1) + B 5046 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 5047 * INDEX(NAI,NBJ) 5048 NAKBL = IT2AM(ISYMKA,ISYMLB) + 5049 * INDEX(NAK,NBL) 5050 V1 = V2AM(NAIBJ) 5051 R1 = R2AM(NAIBJ) 5052 R2 = R2AM(NAKBL) 5053 VR = V1 * R2 5054 RR = R1 * R2 5055 BB = SF(ISB(ISYMA) + INDEX(A,A)) 5056 * + SF(ISB(ISYMB) + INDEX(B,B)) 5057 BB = R1 * R2 * BB 5058C 5059 V11(KOFFI,KOFFJ,KOFFK,KOFFL) = 5060 * V11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 5061 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 5062 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 5063 B11(KOFFI,KOFFJ,KOFFK,KOFFL) = 5064 * B11(KOFFI,KOFFJ,KOFFK,KOFFL) + BB 5065C 5066 111 CONTINUE 5067 110 CONTINUE 5068 109 CONTINUE 5069 108 CONTINUE 5070 107 CONTINUE 5071 106 CONTINUE 5072 105 CONTINUE 5073 104 CONTINUE 5074 103 CONTINUE 5075 102 CONTINUE 5076 101 CONTINUE 5077C 5078 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER PAIR ENERGIES:' 5079 DO 230 I=1,NSPAIR 5080 ES(I) = ES(I) * DP25 5081 ET(I) = ET(I) * DP75 5082 E2S = E2S + ES(I) 5083 E2T = E2T + ET(I) 5084 WRITE(LUPRI,'(17X,I4,2F15.9)') I,ES(I),ET(I) 5085 230 CONTINUE 5086 WRITE(LUPRI,'(/A6,3F15.9)') ' MP2 =',E2,E2S,E2T 5087C 5088 IJ = 0 5089 DO 300 I = 1, NOCDIM 5090 DO 301 J = 1, I 5091 IJ = IJ + 1 5092 KL = 0 5093 DO 302 K = 1, NOCDIM 5094 DO 303 L = 1, K 5095 KL = KL + 1 5096C 5097 RR = V11(I,J,K,L) + V11(I,J,L,K) 5098 VS11(KL,IJ) = RR 5099 IF (I .EQ. J) VS11(KL,IJ) = FF * VS11(KL,IJ) 5100 IF (K .EQ. L) VS11(KL,IJ) = FF * VS11(KL,IJ) 5101 RR = V11(I,J,K,L) - V11(I,J,L,K) 5102 VT11(KL,IJ) = RR 5103C 5104 BB = B11(I,J,K,L) + B11(I,J,L,K) 5105 BS11(KL,IJ) = BB 5106 IF (I .EQ. J) BS11(KL,IJ) = FF * BS11(KL,IJ) 5107 IF (K .EQ. L) BS11(KL,IJ) = FF * BS11(KL,IJ) 5108 BB = B11(I,J,K,L) - B11(I,J,L,K) 5109 BT11(KL,IJ) = BB 5110C 5111 RR = Q11(I,J,K,L) + Q11(I,J,L,K) 5112 RS11(KL,IJ) = RR 5113 IF (I .EQ. J) RS11(KL,IJ) = FF * RS11(KL,IJ) 5114 IF (K .EQ. L) RS11(KL,IJ) = FF * RS11(KL,IJ) 5115 RR = Q11(I,J,K,L) - Q11(I,J,L,K) 5116 RT11(KL,IJ) = RR 5117C 5118 VT11(KL,IJ) = VT11(KL,IJ) * D3 5119 BT11(KL,IJ) = BT11(KL,IJ) * D3 5120 RT11(KL,IJ) = RT11(KL,IJ) * D3 5121C 5122 303 CONTINUE 5123 302 CONTINUE 5124 301 CONTINUE 5125 300 CONTINUE 5126C 5127 IF (IPRINT .LE. 3) GOTO 998 5128 CALL AROUND('Singlet <Vab> matrix') 5129 CALL OUTPUT(VS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 5130 CALL AROUND('Singlet <inglet <V> matrixab> matrix') 5131 CALL OUTPUT(BS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 5132 CALL AROUND('Singlet <Sab> matrix') 5133 CALL OUTPUT(RS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 5134 CALL AROUND('Triplet <Vab> matrix') 5135 CALL OUTPUT(VT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 5136 CALL AROUND('Triplet <Bab> matrix') 5137 CALL OUTPUT(BT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 5138 CALL AROUND('Triplet <Sab> matrix') 5139 CALL OUTPUT(RT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 5140 998 CONTINUE 5141C 5142 CALL DZERO(FS,NSPAIR) 5143 CALL DZERO(FT,NSPAIR) 5144 F2S = D0 5145 F2T = D0 5146 IJ = 0 5147 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER MP2-EC PAIR ENERGIES:' 5148 DO 500 I = 1, NOCDIM 5149 DO 501 J = 1, I 5150 IJ = IJ + 1 5151 IF (R12SVD) THEN 5152 CALL R12MP2(VS11(1,IJ),VT11(1,IJ),VS11(1,IJ),VT11(1,IJ), 5153 * RS11,RT11,EV(I)+EV(J),NOCDIM,NSPAIR,FS(IJ),FT(IJ),WS11,WT11, 5154 * QS11,QT11,BS11,BT11,EVS,D1,V11,IJ,WSMIN,WTMIN,DELTA, 5155 * CS11(1,IJ),CT11(1,IJ)) 5156 ELSE IF (R12DIA) THEN 5157 CALL MP2R12(VS11(1,IJ),VT11(1,IJ),VS11(1,IJ),VT11(1,IJ), 5158 * RS11,RT11,EV(I)+EV(J),NOCDIM,NSPAIR,FS(IJ),FT(IJ),WS11,WT11, 5159 * QS11,QT11,BS11,BT11,EVS,D1,V11,IJ,WSMIN,WTMIN,DELTA, 5160 * CS11(1,IJ),CT11(1,IJ)) 5161 ENDIF 5162 WRITE(LUPRI,'(I4,4F15.9,2G16.6)') IJ,ES(IJ),ES(IJ)+FS(IJ), 5163 * ET(IJ),ET(IJ)+FT(IJ), 5164 * WSMIN,WTMIN 5165 F2S = F2S + FS(IJ) 5166 F2T = F2T + FT(IJ) 5167 501 CONTINUE 5168 500 CONTINUE 5169 WRITE(LUPRI,'(/A,F15.9 )') 5170 * ' MP2-EC correlation energy =', 5171 * E2S+E2T+F2S+F2T 5172 RETURN 5173 END 5174C /* Deck prfdrv */ 5175 SUBROUTINE PRFDRV(V2AM,R2AM,EV,WRK,LWRK) 5176C 5177C Driver for computation of basis-set profiles. 5178C 5179C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 5180C 5181#include "implicit.h" 5182#include "priunit.h" 5183#include "ccorb.h" 5184#include "ccsdsym.h" 5185#include "r12int.h" 5186 DIMENSION R2AM(*), V2AM(*), EV(*), WRK(*) 5187 NOCDIM = NRHFT 5188 NICDIM = NRHFTS 5189 NSPAIR = NOCDIM * (NOCDIM + 1) / 2 5190 NTPAIR = NOCDIM * (NOCDIM - 1) / 2 5191 NPAIR2 = NSPAIR ** 2 5192 NODIM4 = NOCDIM ** 4 5193 NIDIM4 = NICDIM ** 4 5194 IVS = 1 5195 IVT = IVS + NOCDIM * NOCDIM 5196 IENDW = IVT + NOCDIM * NOCDIM 5197 IF (IENDW .GT. LWRK) THEN 5198 WRITE(LUPRI,'(A,I20/A,I20)') 5199 * ' WORK SPACE REQUIRED = ',IENDW, 5200 * ' WORK SPACE AVAILABLE = ',LWRK 5201 CALL QUIT('INSUFFICIENT WORK SPACE IN R12DRV') 5202 END IF 5203 CALL PRFEVR(V2AM,R2AM,WRK(IVS),WRK(IVT),NICDIM,NOCDIM,NSPAIR, 5204 * NTPAIR) 5205 RETURN 5206 END 5207C /* Deck prfevr */ 5208 SUBROUTINE PRFEVR(V2AM,R2AM,VS,VT,NICDIM,NOCDIM,NSPAIR,NTPAIR) 5209C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 5210#include "implicit.h" 5211#include "priunit.h" 5212 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 5213 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0) 5214 DIMENSION R2AM(*), V2AM(*) 5215 DIMENSION VS(NOCDIM,NOCDIM), VT(NOCDIM,NOCDIM) 5216 INTEGER IDUM 5217 LOGICAL LDUM 5218#include "ccorb.h" 5219#include "ccsdsym.h" 5220#include "r12int.h" 5221C 5222 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 5223C 5224 CALL DZERO(VS,NOCDIM**2) 5225 CALL DZERO(VT,NOCDIM**2) 5226C 5227C CONSTRUCT PRODUCTS (IA|JB)*(IA||JB) 5228C 5229 DO 101 ISYMIA = 1, NSYM 5230 ISYMJB = ISYMIA 5231 DO 102 ISYMI = 1, NSYM 5232 ISYMA = MULD2H(ISYMI,ISYMIA) 5233 DO 103 ISYMJ = 1, NSYM 5234 ISYMB = MULD2H(ISYMJ,ISYMJB) 5235 DO 105 I = 1, NRHF(ISYMI) 5236 KOFFI = IRHF(ISYMI) + I 5237 DO 107 A = 1 + NORB1(ISYMA), NVIR(ISYMA) 5238 KOFFA = IVIR(ISYMA) + A 5239 NAI = IT1AM(ISYMA,ISYMI) + 5240 * NVIR(ISYMA)*(I-1) + A 5241 DO 108 J = 1, NRHF(ISYMJ) 5242 KOFFJ = IRHF(ISYMJ) + J 5243 NAJ = IT1AM(ISYMA,ISYMJ) + 5244 * NVIR(ISYMA)*(J-1) + A 5245 ISYMJA = MULD2H(ISYMJ,ISYMA) 5246 DO 110 B = 1 + NORB1(ISYMB), NVIR(ISYMB) 5247 ISYMIB = MULD2H(ISYMI,ISYMB) 5248 KOFFB = IVIR(ISYMB) + B 5249 NBJ = IT1AM(ISYMB,ISYMJ) + 5250 * NVIR(ISYMB)*(J-1) + B 5251 NBI = IT1AM(ISYMB,ISYMI) + 5252 * NVIR(ISYMB)*(I-1) + B 5253 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 5254 * INDEX(NAI,NBJ) 5255 NAJBI = IT2AM(ISYMJA,ISYMIB) + 5256 * INDEX(NAJ,NBI) 5257 VR = V2AM(NAIBJ) * R2AM(NAIBJ) 5258 VP = V2AM(NAIBJ) * R2AM(NAJBI) 5259 VS(KOFFI,KOFFJ) = VS(KOFFI,KOFFJ) + VR + VP 5260 VT(KOFFI,KOFFJ) = VT(KOFFI,KOFFJ) + VR - VP 5261 110 CONTINUE 5262 108 CONTINUE 5263 107 CONTINUE 5264 105 CONTINUE 5265 103 CONTINUE 5266 102 CONTINUE 5267 101 CONTINUE 5268C 5269 LU33 = -33 5270 LU34 = -34 5271 LU35 = -35 5272 CALL GPOPEN(LU33,'PROFILE.1','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 5273 CALL GPOPEN(LU34,'PROFILE.2','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 5274 CALL GPOPEN(LU35,'PROFILE.3','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 5275 DO 201 I = 1, NOCDIM 5276 DO 202 J = 1, NOCDIM 5277 WRITE(LU33,'(2I4,2F20.14)') I,J,VS(I,J) 5278 WRITE(LU34,'(2I4,2F20.14)') I,J,VT(I,J) 5279 WRITE(LU35,'(2I4,2F20.14)') I,J,(VS(I,J)+VT(I,J))*DP5 5280 202 CONTINUE 5281 WRITE(LU33,*) 5282 WRITE(LU34,*) 5283 WRITE(LU35,*) 5284 201 CONTINUE 5285 CALL GPCLOSE(33,'KEEP') 5286 CALL GPCLOSE(34,'KEEP') 5287 CALL GPCLOSE(35,'KEEP') 5288 RETURN 5289 END 5290C /* Deck mpab */ 5291 SUBROUTINE MPAB(A,NROWA,NCOLA,NRDIMA,NCDIMA, 5292 1 B,NROWB,NCOLB,NRDIMB,NCDIMB, 5293 2 C,NRDIMC,NCDIMC) 5294C----------------------------------------------------------- 5295C MATRIX PRODUCT A TIMES B = C. 5296C Written by George D. Purvis 1983 5297C Revised 6-Nov-1984 Hans Joergen Aa. Jensen 5298C------------------------------------------------------------- 5299#include "implicit.h" 5300 DIMENSION A(NRDIMA,NCDIMA),B(NRDIMB,NCDIMB),C(NRDIMC,NCDIMC) 5301 PARAMETER (ZERO=0.D00) 5302C 5303 IF (NCOLA.NE.NROWB) THEN 5304 WRITE (*,9000) NCOLA,NROWB 5305 CALL QUIT('ERROR, inconsistent matrix dimensions in MPAB') 5306 ENDIF 5307 9000 FORMAT(/' MPAB error: NCOLA .ne. NROWB, values =',2I10) 5308C 5309 IF (NROWB .EQ. 0) RETURN 5310 DO 40 J = 1,NCOLB 5311 DO 10 I = 1,NROWA 5312 10 C(I,J) = A(I,1)*B(1,J) 5313 DO 30 K = 2,NROWB 5314 IF (B(K,J).EQ.ZERO) GO TO 30 5315 BKJ = B(K,J) 5316 DO 20 I = 1,NROWA 5317 20 C(I,J) = BKJ*A(I,K) + C(I,J) 5318 30 CONTINUE 5319 40 CONTINUE 5320C 5321 RETURN 5322 END 5323C=====================================================================* 5324C /* Deck vshrnk */ 5325 SUBROUTINE VSHRNK(VS11,NSPAIR,NRHF,NRHFA,NSYM) 5326#include "implicit.h" 5327 DIMENSION VS11(NSPAIR,NSPAIR),NRHF(NSYM),NRHFA(NSYM) 5328 IJ = 0 5329 JI = 0 5330 II = 0 5331 DO ISYM = 1, NSYM 5332 DO I = 1, NRHF(ISYM) 5333 II = II + 1 5334 JJ = 0 5335 DO JSYM = 1, NSYM 5336 DO J = 1, NRHF(JSYM) 5337 JJ = JJ + 1 5338 IF (JJ .LE. II) THEN 5339 IJ = IJ + 1 5340 IF (I .LE. NRHFA(ISYM) .AND. 5341 * J .LE. NRHFA(JSYM)) THEN 5342 JI = JI + 1 5343 IF (IJ .NE. JI) 5344 * CALL DCOPY(NSPAIR,VS11(1,IJ),1,VS11(1,JI),1) 5345 END IF 5346 END IF 5347 END DO 5348 END DO 5349 END DO 5350 END DO 5351 IF (JI .LT. NSPAIR) 5352 *CALL DZERO(VS11(1,JI+1),NSPAIR*(NSPAIR-JI)) 5353 RETURN 5354 END 5355