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 makerv*/ 20 SUBROUTINE MAKERV(V2AM,R2AM,U2AM,S2AM,T2AM, 21 * EV,QQ11,V11,U11,B11,W11,Q11,R11, 22 * VS11,US11,BS11,WS11, 23 * QS11,RS11,VT11,UT11,BT11,WT11,QT11,RT11, 24 * NICDIM,NOCDIM,NSPAIR,NTPAIR,ES,ET,FS,FT,EVS, 25 * CNINV,QQ2,QQ4,QQ6,CS11,CT11,IANSAZ,FRSTWR, 26 * WORK,LWORK,NACDIM,NAPAIR) 27C 28C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 29C Ansatz 3 added on 15 May 2004. 30C 31#include "implicit.h" 32#include "priunit.h" 33 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 34 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0, DELTA=0.0D0) 35 REAL*8 SU, US, SV, VU, VR, RR, UR, RU, UT, TU, QQIJKL, QQIJLK, 36 * EIJAB, EKLAB, EIJ, EMN, EKL, EIJKL, EAB, EMNAB, UR3, 37 * SAIBJ, RAKBL, UAKBL, SAKBL, RAIBJ, VAIBJ, UAIBJ, 38 * TAIBJ, TAKBL, 39 * V11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 40 * W11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 41 * B11(NOCDIM,NOCDIM,NOCDIM,NOCDIM,NAPAIR), 42 * U11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 43 * Q11(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 44 * R11(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 45 DIMENSION QQ2(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 46 * QQ4(NOCDIM,NOCDIM,NOCDIM,NOCDIM), 47 * QQ6(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 48 DIMENSION R2AM(*), V2AM(*), U2AM(*), S2AM(*), T2AM(*), EV(*) 49 DIMENSION WORK(*) 50 DIMENSION VS11(NSPAIR,NSPAIR) 51 DIMENSION US11(NSPAIR,NSPAIR) 52 DIMENSION WS11(NSPAIR,NSPAIR) 53 DIMENSION BS11(NSPAIR,NSPAIR,NSPAIR) 54 DIMENSION WT11(NSPAIR,NSPAIR) 55 DIMENSION BT11(NSPAIR,NSPAIR,NSPAIR) 56 DIMENSION QS11(NSPAIR,NSPAIR) 57 DIMENSION RS11(NSPAIR,NSPAIR) 58 DIMENSION VT11(NSPAIR,NSPAIR),UT11(NSPAIR,NSPAIR) 59 DIMENSION QT11(NSPAIR,NSPAIR),RT11(NSPAIR,NSPAIR) 60 DIMENSION ES(NSPAIR),ET(NSPAIR),FS(NSPAIR),FT(NSPAIR) 61 DIMENSION EVS(NSPAIR),CNINV(NSPAIR,10) 62 DIMENSION QQ11(NICDIM,NICDIM,NICDIM,NICDIM) 63 DIMENSION CS11(NSPAIR,NSPAIR),CT11(NSPAIR,NSPAIR) 64 LOGICAL DOA, DOB, DOAP, DOBP, ASTRIX, DONINV, LDUM 65 INTEGER IDUM, JDUM 66 CHARACTER*3 APPROX(0:14,2), IPCC 67 CHARACTER*8 EOFLAB 68 LOGICAL FRSTWR 69 INTEGER NKILJ(8) 70#include "r12int.h" 71#include "ccorb.h" 72#include "ccsdsym.h" 73#include "ccsdinp.h" 74#include "ccr12int.h" 75 DATA APPROX 76 &/'0 ','A ','A ','A'' ','A* ','A* ','A*''','B ', 77 & 'B ','xxx','xxx','B* ','B* ','xxx','xxx', 78 & '0 ','A ','A ','A'' ','A* ','A* ','A*''','B'' ', 79 & 'B'' ','B ','B ','B*''','B*''','B* ','B* '/ 80 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 81C 82 LU43 = -43 83 LU44 = -44 84 LU45 = -45 85 LU46 = -46 86 CALL GPOPEN(LU46,FCCR12K,'UNKNOWN',' ','FORMATTED', 87 & IDUM,LDUM) 88 FF = D1 / SQRT(D2) 89C 90 CALL AROUND('Auxiliary-basis MP2-R12 method') 91 DO 999 IPDD = 0, 14 92 IF (R12HYB) THEN 93 IF (IANSAZ .EQ. 3 .AND. .NOT. (IPDD .EQ. 7 .OR. IPDD .EQ. 8)) 94 & GOTO 999 95 IPCC = APPROX(IPDD,1) 96 ELSE 97 IF (IANSAZ .EQ. 3 .AND. .NOT. (IPDD .GE. 7 .AND. IPDD .LE. 10)) 98 & GOTO 999 99 IPCC = APPROX(IPDD,2) 100 END IF 101 IF (.NOT. R12XXL) THEN 102 IF (NORXR) THEN 103 IF (IPDD .NE. 2 .AND. IPDD .NE. 8) GOTO 999 104 ELSE 105 IF (IPDD .NE. 2 .AND. IPDD .NE. 10) GOTO 999 106 END IF 107 END IF 108C 109 DOA = IPDD .LE. 6 110 DOB = .NOT. DOA 111 DONINV = IPDD .LE. 1 .OR. IPDD .EQ. 4 .OR. 112 * IPDD .EQ. 7 .OR. IPDD .EQ. 9 .OR. 113 * IPDD .EQ. 11 .OR. IPDD .EQ. 13 114 DOAP = MOD(IPDD, 3) .EQ. 0 .AND. DOA 115 DOBP = IPDD .EQ. 9 .OR. IPDD .EQ. 10 .OR. IPDD .EQ. 13 116 * .OR. IPDD .EQ. 14 117 ASTRIX = (IPDD .GE. 4 .AND. IPDD .LE. 6) .OR. 118 * (IPDD .GE. 11 .AND. IPDD .LE. 14) 119 IF (DOA .AND. R12NOA) GOTO 999 120 IF (DOB .AND. R12NOB) GOTO 999 121 IF (DOAP .AND. R12NOP) GOTO 999 122 IF (DOBP .AND. NORXR) GOTO 999 123 IF (IPDD .EQ. 2 .OR. IPDD .EQ. 5) THEN 124 XPDD = D0 125 ELSE 126 XPDD = DP5 127 END IF 128 129 LU33 = -33 130 CALL GPOPEN(LU33,'AUXQ12','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 131 READ(LU33,'(4E30.20)') QQ11 132 CALL GPCLOSE(LU33,'KEEP') 133C 134 DO 205 I = 1, NOCDIM**4 135 R11(I,1,1,1) = D0 136 Q11(I,1,1,1) = D0 137 V11(I,1,1,1) = D0 138 W11(I,1,1,1) = D0 139 U11(I,1,1,1) = D0 140 205 CONTINUE 141 DO 206 I = 1, NAPAIR*NOCDIM**4 142 B11(I,1,1,1,1) = D0 143 206 CONTINUE 144 CALL DZERO(ES,NSPAIR) 145 CALL DZERO(ET,NSPAIR) 146C 147C CONSTRUCT PRODUCTS (IA|JB)*(KA|LB) 148C 149 E2 = D0 150 E2S = D0 151 E2T = D0 152 DO 101 ISYMIA = 1, NSYM 153 ISYMJB = ISYMIA 154 JBOFF = NH1AM(ISYMJB) * (NH1AM(ISYMJB) + 1) / 2 155 DO 102 ISYMI = 1, NSYM 156 ISYMA = MULD2H(ISYMI,ISYMIA) 157 DO 103 ISYMJ = 1, NSYM 158 ISYMB = MULD2H(ISYMJ,ISYMJB) 159C 160C COMPUTE MP2 ENERGY 161C 162 DO 201 I = 1, NRHFA(ISYMI) 163 KOFFI = IRHF(ISYMI) + I 164 KOFFIA = IRHFA(ISYMI) + I 165 DO 202 J = 1, NRHFA(ISYMJ) 166 KOFFJ = IRHF(ISYMJ) + J 167 KOFFJA = IRHFA(ISYMJ) + J 168 DO 203 A = NRHFSA(ISYMA) + 1, NORB1(ISYMA) 169 ISYMJA = MULD2H(ISYMJ,ISYMA) 170 KOFFA = IVIR(ISYMA) + A 171 DO 204 B = NRHFSA(ISYMB) + 1, NORB1(ISYMB) 172 ISYMIB = MULD2H(ISYMI,ISYMB) 173 KOFFB = IVIR(ISYMB) + B 174 IF (ONEAUX) THEN 175 NAI = IH1AM(ISYMA,ISYMI) + 176 * NORB1(ISYMA)*(I-1) + A 177 NBJ = IH1AM(ISYMB,ISYMJ) + 178 * NORB1(ISYMB)*(J-1) + B 179 NAIBJ = IH2AM(ISYMIA,ISYMJB) + 180 * INDEX(NAI,NBJ) 181 NAJ = IH1AM(ISYMA,ISYMJ) + 182 * NORB1(ISYMA)*(J-1) + A 183 NBI = IH1AM(ISYMB,ISYMI) + 184 * NORB1(ISYMB)*(I-1) + B 185 NAJBI = IH2AM(ISYMJA,ISYMIB) + 186 * INDEX(NAJ,NBI) 187 ELSE 188 NAI = IT1AM(ISYMA,ISYMI) + 189 * NVIR(ISYMA)*(I-1) + A 190 NBJ = IT1AM(ISYMB,ISYMJ) + 191 * NVIR(ISYMB)*(J-1) + B 192 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 193 * INDEX(NAI,NBJ) 194 NAJ = IT1AM(ISYMA,ISYMJ) + 195 * NVIR(ISYMA)*(J-1) + A 196 NBI = IT1AM(ISYMB,ISYMI) + 197 * NVIR(ISYMB)*(I-1) + B 198 NAJBI = IT2AM(ISYMJA,ISYMIB) + 199 * INDEX(NAJ,NBI) 200 END IF 201 VAIBJ = V2AM(NAIBJ) 202 VAJBI = V2AM(NAJBI) 203 VV = VAIBJ * (D2 * VAIBJ - VAJBI) 204 DENOM = D1 / (EV(KOFFI) + EV(KOFFJ) - 205 * EV(KOFFA) - EV(KOFFB)) 206 E2 = E2 + VV * DENOM 207 IJ = INDEX(KOFFIA,KOFFJA) 208 VS = (VAIBJ + VAJBI)**2 209 VT = (VAIBJ - VAJBI)**2 210 ES(IJ) = ES(IJ) + VS * DENOM 211 ET(IJ) = ET(IJ) + VT * DENOM 212 204 CONTINUE 213 203 CONTINUE 214 202 CONTINUE 215 201 CONTINUE 216C 217 DO 104 ISYMKA = 1, NSYM 218 ISYMLB = ISYMKA 219 LBOFF = NH1AM(ISYMLB) * (NH1AM(ISYMLB) + 1) / 2 220 ISYMK = MULD2H(ISYMA,ISYMKA) 221 ISYML = MULD2H(ISYMB,ISYMLB) 222 DO 105 I = 1, NRHF(ISYMI) 223 KOFFI = IRHF(ISYMI) + I 224 DO 106 K = 1, NRHF(ISYMK) 225 KOFFK = IRHF(ISYMK) + K 226 DO 107 A = 1, NVIR(ISYMA) 227 KOFFA = IVIR(ISYMA) + A 228 IF (ONEAUX) THEN 229 IF (A. LE. NORB1(ISYMA)) THEN 230 NAI = IH1AM(ISYMA,ISYMI) + 231 * NORB1(ISYMA)*(I-1) + A 232 NAK = IH1AM(ISYMA,ISYMK) + 233 * NORB1(ISYMA)*(K-1) + A 234 ELSE 235 NA1 = A - NORB1(ISYMA) 236 NAI = IG1AM(ISYMA,ISYMI) + 237 * NORB2(ISYMA)*(I-1) + NA1 238 NAK = IG1AM(ISYMA,ISYMK) + 239 * NORB2(ISYMA)*(K-1) + NA1 240 END IF 241 ELSE 242 NAI = IT1AM(ISYMA,ISYMI) + 243 * NVIR(ISYMA)*(I-1) + A 244 NAK = IT1AM(ISYMA,ISYMK) + 245 * NVIR(ISYMA)*(K-1) + A 246 END IF 247 DO 108 J = 1, NRHF(ISYMJ) 248 KOFFJ = IRHF(ISYMJ) + J 249 DO 109 L = 1, NRHF(ISYML) 250 KOFFL = IRHF(ISYML) + L 251 EIJKL = EV(KOFFI) + EV(KOFFJ) 252 * - EV(KOFFK) - EV(KOFFL) 253 IF (ONEAUX) THEN 254 NBEND = NORB1(ISYMB) 255 ELSE 256 NBEND = NVIR(ISYMB) 257 END IF 258 DO 110 B = 1, NBEND 259 KOFFB = IVIR(ISYMB) + B 260 IF (ONEAUX) THEN 261 NBJ = IH1AM(ISYMB,ISYMJ) + 262 * NORB1(ISYMB)*(J-1) + B 263 NBL = IH1AM(ISYMB,ISYML) + 264 * NORB1(ISYMB)*(L-1) + B 265 IF (A .GT. NORB1(ISYMA)) THEN 266 NAIBJ = IH2AM(ISYMIA,ISYMJB) + 267 * JBOFF + NBJ + 268 * NH1AM(ISYMJB)*(NAI-1) 269 NAKBL = IH2AM(ISYMKA,ISYMLB) + 270 * LBOFF + NBL + 271 * NH1AM(ISYMLB)*(NAK-1) 272 ELSE 273 NAIBJ = IH2AM(ISYMIA,ISYMJB) + 274 * INDEX(NAI,NBJ) 275 NAKBL = IH2AM(ISYMKA,ISYMLB) + 276 * INDEX(NAK,NBL) 277 END IF 278 ELSE 279 NBJ = IT1AM(ISYMB,ISYMJ) + 280 * NVIR(ISYMB)*(J-1) + B 281 NBL = IT1AM(ISYMB,ISYML) + 282 * NVIR(ISYMB)*(L-1) + B 283 NAIBJ = IT2AM(ISYMIA,ISYMJB) + 284 * INDEX(NAI,NBJ) 285 NAKBL = IT2AM(ISYMKA,ISYMLB) + 286 * INDEX(NAK,NBL) 287 END IF 288C 289 EIJAB = EV(KOFFI) + EV(KOFFJ) 290 * - EV(KOFFA) - EV(KOFFB) 291 EKLAB = EV(KOFFK) + EV(KOFFL) 292 * - EV(KOFFA) - EV(KOFFB) 293C 294 RAIBJ = R2AM(NAIBJ) 295 VAIBJ = V2AM(NAIBJ) 296 UAIBJ = U2AM(NAIBJ) 297 SAIBJ = S2AM(NAIBJ) 298 TAIBJ = T2AM(NAIBJ) 299C 300 RAKBL = R2AM(NAKBL) 301 UAKBL = U2AM(NAKBL) 302 SAKBL = S2AM(NAKBL) 303 TAKBL = T2AM(NAKBL) 304C 305 IF (ASTRIX) THEN 306C 307C SU := (Ek+El-Ea-Eb) * <ab|r12|kl> 308C US := (Ei+Ej-Ea-Eb) * <ab|r12|ij> 309C 310 SU = EKLAB * RAKBL 311 US = EIJAB * RAIBJ 312 ELSE 313C 314C SU := - <ab|[T1+T2,r12]|kl> + <ab|[K1+K2,r12]|kl> 315C US := - <ab|[T1+T2,r12]|ij> + <ab|[K1+K2,r12]|kl> 316C 317 SU = UAKBL - SAKBL 318 US = UAIBJ - SAIBJ 319 ENDIF 320C 321C SV := - <ab|(F1+F2-Ei-Ej)r12|kl> = 322C - <ab|[F1+F2,r12]|kl> 323C + (Ei+Ej-Ek-El) * <ab|r12|kl> = 324C = C^(ij)_kl,ab 325C 326 SV = SU + EIJKL * RAKBL 327C 328 IF (IANSAZ .EQ. 3) THEN 329 SV = SV - EIJAB * RAKBL 330 END IF 331C 332 VU = VAIBJ * SV 333 VR = VAIBJ * RAKBL 334 RR = RAIBJ * RAKBL 335 UR = UAIBJ * RAKBL + RAIBJ * UAKBL 336 IF (IANSAZ .EQ. 3) THEN 337 UR3 = (UAIBJ - SAIBJ) * RAKBL 338 * + RAIBJ * (UAKBL - SAKBL) 339 UR3 = UR3 + UAIBJ * RAKBL 340 * + RAIBJ * UAKBL 341 IF (DOBP) THEN 342 UR3 = UR3 - SAIBJ * RAKBL 343 * - RAIBJ * SAKBL 344 ELSE 345 UR3 = UR3 - TAIBJ * RAKBL 346 * - RAIBJ * TAKBL 347 END IF 348 UR3 = UR3 - (EIJAB + EKLAB) * RR 349 END IF 350 IF (R12CBS) THEN 351 NQUEA = 0 352 IF (ONEAUX) THEN 353 NQUEB = NORB1(ISYMB) 354 ELSE 355 NQUEB = 0 356 END IF 357 ELSE 358 NQUEA = NORB1(ISYMA) 359 NQUEB = NORB1(ISYMB) 360 END IF 361 362 IF (A .LE. NRHFSA(ISYMA) .AND. 363 * B .LE. NRHFSA(ISYMB)) THEN 364C 365C Sum |ij><ij| 366 367 V11(KOFFI,KOFFJ,KOFFK,KOFFL) = 368 * V11(KOFFI,KOFFJ,KOFFK,KOFFL) - VR 369 U11(KOFFI,KOFFJ,KOFFK,KOFFL) = 370 * U11(KOFFI,KOFFJ,KOFFK,KOFFL) - UR 371 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 372 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) - RR 373C 374 END IF 375C 376 IF (A .LE. NRHFSA(ISYMA) .AND. 377 * B .GT. NQUEB) THEN 378C 379C Sum |iq'><iq'| 380C 381 V11(KOFFI,KOFFJ,KOFFK,KOFFL) = 382 * V11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 383 U11(KOFFI,KOFFJ,KOFFK,KOFFL) = 384 * U11(KOFFI,KOFFJ,KOFFK,KOFFL) + UR 385 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 386 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 387C 388 END IF 389C 390 IF (A .GT. NQUEA .AND. 391 * B .LE. NRHFSA(ISYMB)) THEN 392C 393C Sum |p'j><p'j| 394C 395 V11(KOFFI,KOFFJ,KOFFK,KOFFL) = 396 * V11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 397 U11(KOFFI,KOFFJ,KOFFK,KOFFL) = 398 * U11(KOFFI,KOFFJ,KOFFK,KOFFL) + UR 399 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 400 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 401 IF (ONEAUX) THEN 402 V11(KOFFJ,KOFFI,KOFFL,KOFFK) = 403 * V11(KOFFJ,KOFFI,KOFFL,KOFFK) + VR 404 U11(KOFFJ,KOFFI,KOFFL,KOFFK) = 405 * U11(KOFFJ,KOFFI,KOFFL,KOFFK) + UR 406 Q11(KOFFJ,KOFFI,KOFFL,KOFFK) = 407 * Q11(KOFFJ,KOFFI,KOFFL,KOFFK) + RR 408 END IF 409C 410 END IF 411C 412 IF (A. GT. NRHFSA(ISYMA) .AND. 413 * A .LE. NORB1(ISYMA) .AND. 414 * B .GT. NRHFSA(ISYMB) .AND. 415 * B .LE. NORB1(ISYMB)) THEN 416C 417C Sum |ab><ab| 418C 419 IF (IANSAZ .EQ. 3) THEN 420 V11(KOFFI,KOFFJ,KOFFK,KOFFL) = 421 * V11(KOFFI,KOFFJ,KOFFK,KOFFL) + VR 422 U11(KOFFI,KOFFJ,KOFFK,KOFFL) = 423 * U11(KOFFI,KOFFJ,KOFFK,KOFFL) + 424 * UR3 425 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 426 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + RR 427 END IF 428C 429C Write out C^(ij)_kl,ab for 430C approximation A, A', B, and *, resp., 431C for CC2-R12 model (WK/UniKA/30-12-2003). 432C 433C 434 IF (FRSTWR 435 & .AND. A. LE. NORB1(ISYMA) 436 & .AND. B. LE. NORB1(ISYMB)) 437 & WRITE (LU46,'(6I5,2E30.20/ 438 & 30X,2E30.20)') 439 & KOFFI,KOFFJ,KOFFK,KOFFL, 440 & KOFFA,KOFFB, 441 & - UAKBL, 442 & -(UAKBL + EIJKL * RAKBL), 443 & -(UAKBL - SAKBL + EIJKL * RAKBL), 444 & -(EKLAB * RAKBL + EIJKL * RAKBL) 445C 446 EIJ = D1 / 447 * (EV(KOFFI) + EV(KOFFJ) - 448 * EV(KOFFA) - EV(KOFFB)) 449! 450! hjaaj: sometimes EIJ = Infinity or -Infinity 451! apparently W11 with these values are not used 452! However, ifort and -i8 stops with "invalid float from "-Infinity + Infinity" 453! My solution: reset value to 1.D100, then it will be obvious if the value is used. 454 IF (abs(eij) .gt. 1.D100) 455 & eij = 1.D100 456! 457C 458 W11(KOFFI,KOFFJ,KOFFK,KOFFL) = 459 * W11(KOFFI,KOFFJ,KOFFK,KOFFL) - 460 * VU * EIJ 461C 462 MN = 0 463 NM = 0 464 M = 0 465 DO 111 MSYMM = 1, NSYM 466 DO 111 MM = 1, NRHF(MSYMM) 467 M = M + 1 468 N = 0 469 DO 112 NSYMM = 1, NSYM 470 DO 112 NN = 1, NRHF(NSYMM) 471 N = N + 1 472 IF (N .GT. M) GOTO 111 473 NM = NM + 1 474 IF (MM .GT. NRHFA(MSYMM) 475 * .OR. 476 * NN .GT. NRHFA(NSYMM)) 477 * GOTO 112 478 MN = MN + 1 479 EMNAB = EV(M) + EV(N) 480 * - EV(KOFFA) - EV(KOFFB) 481 EMN = D1 / EMNAB 482 EIJ = EV(M) + EV(N) 483 * - EV(KOFFI) - EV(KOFFJ) 484 EKL = EV(M) + EV(N) 485 * - EV(KOFFK) - EV(KOFFL) 486C 487C RU := (Em+En-Ea-Eb) * <ab|r12|kl> 488C UR := (Em+En-Ea-Eb) * <ab|r12|ij> 489 RU = SU + EKL * RAKBL 490 UR = US + EIJ * RAIBJ 491C 492 IF (IANSAZ .EQ. 3) THEN 493 RU = RU - EMNAB * RAKBL 494 UR = UR - EMNAB * RAIBJ 495 END IF 496C 497 IF (DOA) THEN 498C 499C Eq. (68) 500C 501 TU = UAKBL 502 UT = UAIBJ 503 ELSE IF (DOBP) THEN 504C 505C Eq. (70) 506C 507 TU = UAKBL - SAKBL 508 UT = UAIBJ - SAIBJ 509 ELSE 510C 511C Eq. (70) without last sum over r' 512C 513 TU = UAKBL - TAKBL 514 UT = UAIBJ - TAIBJ 515 END IF 516 IF (DOAP .OR. DOB) THEN 517C 518C Contribution from Eq. (69) 519C 520C TU = C^(mn)_kl,ab 521C 522 TU = TU + EKL * RAKBL 523C 524C UT = C^(mn)_ij,ab 525C 526 UT = UT + EIJ * RAIBJ 527 END IF 528C 529 IF (IANSAZ .EQ. 3) THEN 530 TU = TU - EMNAB * RAKBL 531 UT = UT - EMNAB * RAIBJ 532 END IF 533C 534 UU = RU * UT + TU * UR 535 B11(KOFFI,KOFFJ, 536 * KOFFK,KOFFL,MN) = 537 * B11(KOFFI,KOFFJ, 538 * KOFFK,KOFFL,MN) - 539 * UU * EMN 540C 541 112 CONTINUE 542 111 CONTINUE 543C 544 END IF 545C 546 110 CONTINUE 547 109 CONTINUE 548 108 CONTINUE 549 107 CONTINUE 550 106 CONTINUE 551 105 CONTINUE 552 104 CONTINUE 553 103 CONTINUE 554 102 CONTINUE 555 101 CONTINUE 556 557 DO 230 I=1,NAPAIR 558 ES(I) = ES(I) * DP25 559 ET(I) = ET(I) * DP75 560 E2S = E2S + ES(I) 561 E2T = E2T + ET(I) 562 230 CONTINUE 563C 564 IJ = 0 565 DO 300 I = 1, NOCDIM 566 DO 301 J = 1, I 567 IJ = IJ + 1 568 KL = 0 569 DO 302 K = 1, NOCDIM 570 DO 303 L = 1, K 571 KL = KL + 1 572C 573 IF (R12EOR) THEN 574 RR = -D2 *(QQ6(I,J,K,L) + QQ6(I,J,L,K)) 575 RR = RR - (U11(I,J,K,L) + U11(I,J,L,K)) 576 DX = D0 577 ELSE 578 RR = - (U11(I,J,K,L) + U11(I,J,L,K)) 579 DX = D2 580 END IF 581 US11(KL,IJ) = RR 582 IF (I .EQ. J) US11(KL,IJ) = FF * US11(KL,IJ) 583 IF (K .EQ. L) US11(KL,IJ) = FF * US11(KL,IJ) 584 IF (IJ .EQ. KL) US11(KL,IJ) = US11(KL,IJ) - DX 585 IF (R12EOR) THEN 586 RR = -D2 *(QQ6(I,J,K,L) - QQ6(I,J,L,K)) 587 RR = RR - (U11(I,J,K,L) - U11(I,J,L,K)) 588 DX = D0 589 ELSE 590 RR = - (U11(I,J,K,L) - U11(I,J,L,K)) 591 DX = D2 592 END IF 593 UT11(KL,IJ) = RR 594 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 595 * UT11(KL,IJ) = UT11(KL,IJ) - DX 596C 597 IF (R12EOR) THEN 598 RR = (QQ2(I,J,K,L) + QQ2(I,J,L,K)) 599 RR = RR - (V11(I,J,K,L) + V11(I,J,L,K)) 600 DX = D0 601 ELSE 602 RR = - (V11(I,J,K,L) + V11(I,J,L,K)) 603 DX = D1 604 END IF 605 VS11(KL,IJ) = RR 606 IF (I .EQ. J) VS11(KL,IJ) = FF * VS11(KL,IJ) 607 IF (K .EQ. L) VS11(KL,IJ) = FF * VS11(KL,IJ) 608 IF (IJ .EQ. KL) VS11(KL,IJ) = VS11(KL,IJ) + DX 609 IF (R12EOR) THEN 610 RR = (QQ2(I,J,K,L) - QQ2(I,J,L,K)) 611 RR = RR - (V11(I,J,K,L) - V11(I,J,L,K)) 612 DX = D0 613 ELSE 614 RR = - (V11(I,J,K,L) - V11(I,J,L,K)) 615 DX = D1 616 END IF 617 VT11(KL,IJ) = RR 618 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 619 * VT11(KL,IJ) = VT11(KL,IJ) + DX 620C 621 IF (R12EOR) THEN 622 RR = (QQ2(I,J,K,L) + QQ2(I,J,L,K)) 623 RR = RR - (V11(I,J,K,L) + V11(I,J,L,K)) 624 DX = D0 625 ELSE 626 RR = - (V11(I,J,K,L) + V11(I,J,L,K)) 627 DX = D1 628 END IF 629 RR = RR + (W11(I,J,K,L) + W11(I,J,L,K)) 630 WS11(KL,IJ) = RR 631 IF (I .EQ. J) WS11(KL,IJ) = FF * WS11(KL,IJ) 632 IF (K .EQ. L) WS11(KL,IJ) = FF * WS11(KL,IJ) 633 IF (IJ .EQ. KL) WS11(KL,IJ) = WS11(KL,IJ) + DX 634 IF (R12EOR) THEN 635 RR = (QQ2(I,J,K,L) - QQ2(I,J,L,K)) 636 RR = RR - (V11(I,J,K,L) - V11(I,J,L,K)) 637 DX = D0 638 ELSE 639 RR = - (V11(I,J,K,L) - V11(I,J,L,K)) 640 DX = D1 641 END IF 642 RR = RR + (W11(I,J,K,L) - W11(I,J,L,K)) 643 WT11(KL,IJ) = RR 644 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 645 * WT11(KL,IJ) = WT11(KL,IJ) + DX 646C 647 DO 210 MN = 1, NAPAIR 648 IF (R12EOR) THEN 649 RR = -D2 *(QQ6(I,J,K,L) + QQ6(I,J,L,K)) 650 RR = RR - (U11(I,J,K,L) + U11(I,J,L,K)) 651 DX = D0 652 ELSE 653 RR = - (U11(I,J,K,L) + U11(I,J,L,K)) 654 DX = D2 655 END IF 656 RR = RR + (B11(I,J,K,L,MN) + B11(I,J,L,K,MN)) 657 BS11(KL,IJ,MN) = RR 658 IF (I .EQ. J) BS11(KL,IJ,MN) = FF * BS11(KL,IJ,MN) 659 IF (K .EQ. L) BS11(KL,IJ,MN) = FF * BS11(KL,IJ,MN) 660 IF (IJ .EQ. KL) 661 * BS11(KL,IJ,MN) = BS11(KL,IJ,MN) - DX 662 IF (R12EOR) THEN 663 RR = -D2 *(QQ6(I,J,K,L) - QQ6(I,J,L,K)) 664 RR = RR - (U11(I,J,K,L) - U11(I,J,L,K)) 665 DX = D0 666 ELSE 667 RR = - (U11(I,J,K,L) - U11(I,J,L,K)) 668 DX = D2 669 END IF 670 RR = RR + (B11(I,J,K,L,MN) - B11(I,J,L,K,MN)) 671 BT11(KL,IJ,MN) = RR 672 IF (IJ .EQ. KL .AND. I .NE. J .AND. K. NE. L) 673 * BT11(KL,IJ,MN) = BT11(KL,IJ,MN) - DX 674 210 CONTINUE 675C 676 IF (R12EOR) THEN 677 QQIJKL = QQ4(I,J,K,L) 678 QQIJLK = QQ4(I,J,L,K) 679 ELSE 680 QQIJKL = QQ11(IQ(I),IQ(J),IQ(K),IQ(L)) 681 QQIJLK = QQ11(IQ(I),IQ(J),IQ(L),IQ(K)) 682 END IF 683 RR = (QQIJKL + QQIJLK) 684 * - (Q11(I,J,K,L) + Q11(I,J,L,K)) 685 QS11(KL,IJ) = RR 686 IF (I .EQ. J) QS11(KL,IJ) = FF * QS11(KL,IJ) 687 IF (K .EQ. L) QS11(KL,IJ) = FF * QS11(KL,IJ) 688 RR = (QQIJKL - QQIJLK) 689 * - (Q11(I,J,K,L) - Q11(I,J,L,K)) 690 QT11(KL,IJ) = RR 691C 692 RR = (QQIJKL + QQIJLK) 693 * - (Q11(I,J,K,L) + Q11(I,J,L,K)) 694 * + (R11(I,J,K,L) + R11(I,J,L,K)) 695 RS11(KL,IJ) = RR 696 IF (I .EQ. J) RS11(KL,IJ) = FF * RS11(KL,IJ) 697 IF (K .EQ. L) RS11(KL,IJ) = FF * RS11(KL,IJ) 698 RR = (QQIJKL - QQIJLK) 699 * - (Q11(I,J,K,L) - Q11(I,J,L,K)) 700 * + (R11(I,J,K,L) - R11(I,J,L,K)) 701 RT11(KL,IJ) = RR 702C 703 US11(KL,IJ) = US11(KL,IJ) * DP5 * DP25 704 VS11(KL,IJ) = VS11(KL,IJ) * DP5 705 WS11(KL,IJ) = WS11(KL,IJ) * DP5 706 QS11(KL,IJ) = QS11(KL,IJ) * DP25 707 RS11(KL,IJ) = RS11(KL,IJ) * DP25 708C 709 UT11(KL,IJ) = UT11(KL,IJ) * D1P5 * DP25 710 VT11(KL,IJ) = VT11(KL,IJ) * D1P5 711 WT11(KL,IJ) = WT11(KL,IJ) * D1P5 712 QT11(KL,IJ) = QT11(KL,IJ) * DP75 713 RT11(KL,IJ) = RT11(KL,IJ) * DP75 714C 715 DO 211 MN = 1, NAPAIR 716 BS11(KL,IJ,MN) = BS11(KL,IJ,MN) * DP5 * DP25 717 BT11(KL,IJ,MN) = BT11(KL,IJ,MN) * D1P5 * DP25 718 211 CONTINUE 719c 720 303 CONTINUE 721 302 CONTINUE 722 301 CONTINUE 723 300 CONTINUE 724C 725 CALL VSHRNK(VS11,NSPAIR,NRHF,NRHFA,NSYM) 726 CALL VSHRNK(VT11,NSPAIR,NRHF,NRHFA,NSYM) 727 CALL VSHRNK(WS11,NSPAIR,NRHF,NRHFA,NSYM) 728 CALL VSHRNK(WT11,NSPAIR,NRHF,NRHFA,NSYM) 729C 730 IF (IPRINT .LE. 3) GOTO 998 731 WRITE(LUPRI,*) 'IPDD, IANSAZ = ', IPDD, IANSAZ 732 CALL AROUND('Singlet <W@> matrix') 733 CALL OUTPUT(WS11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 734 WRITE(LUPRI,*) 'Norm of Singlet <W@> matrix =', 735 * DDOT(NSPAIR*NAPAIR,WS11,1,WS11,1) 736 CALL AROUND('Singlet <V@> matrix') 737 CALL OUTPUT(VS11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 738 WRITE(LUPRI,*) 'Norm of Singlet <V@> matrix =', 739 * DDOT(NSPAIR*NAPAIR,VS11,1,VS11,1) 740 CALL AROUND('Singlet <B@> matrix') 741 CALL OUTPUT(US11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 742 WRITE(LUPRI,*) 'Norm of Singlet <B@> matrix =', 743 * DDOT(NSPAIR*NSPAIR,US11,1,US11,1) 744 CALL AROUND('Singlet <S@> matrix') 745 CALL OUTPUT(RS11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 746 WRITE(LUPRI,*) 'Norm of Singlet <S@> matrix =', 747 * DDOT(NSPAIR*NSPAIR,RS11,1,RS11,1) 748 CALL AROUND('Triplet <W@> matrix') 749 CALL OUTPUT(WT11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 750 WRITE(LUPRI,*) 'Norm of triplet <W@> matrix =', 751 * DDOT(NSPAIR*NAPAIR,WT11,1,WT11,1) 752 CALL AROUND('Triplet <V@> matrix') 753 CALL OUTPUT(VT11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 754 WRITE(LUPRI,*) 'Norm of Triplet <V@> matrix =', 755 * DDOT(NSPAIR*NAPAIR,VT11,1,VT11,1) 756 CALL AROUND('Triplet <B@> matrix') 757 CALL OUTPUT(UT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 758 WRITE(LUPRI,*) 'Norm of Triplet <B@> matrix =', 759 * DDOT(NSPAIR*NSPAIR,UT11,1,UT11,1) 760 CALL AROUND('Triplet <S@> matrix') 761 CALL OUTPUT(RT11,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 762 WRITE(LUPRI,*) 'Norm of Triplet <S@> matrix =', 763 * DDOT(NSPAIR*NSPAIR,RT11,1,RT11,1) 764 CALL AROUND('Singlet <B#> matrices') 765 DO MN = 1, NAPAIR 766 WRITE(LUPRI,'(A,I4)') ' PAIR =', MN 767 CALL OUTPUT(BS11(1,1,MN),1, 768 & NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 769 WRITE(LUPRI,*) 'Norm of Singlet <B#> matrix =', 770 * DDOT(NSPAIR*NSPAIR,BS11(1,1,MN),1,BS11(1,1,MN),1) 771 END DO 772 CALL AROUND('Triplet <B#> matrices') 773 DO MN = 1, NAPAIR 774 WRITE(LUPRI,'(A,I4)') ' PAIR =', MN 775 CALL OUTPUT(BT11(1,1,MN),1, 776 & NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 777 WRITE(LUPRI,*) 'Norm of Triplet <B#> matrix =', 778 * DDOT(NSPAIR*NSPAIR,BT11(1,1,MN),1,BT11(1,1,MN),1) 779 END DO 780C 781 998 IJ = 0 782 DO 490 I = 1, NOCDIM 783 DO 491 J = 1, I 784 IJ = IJ + 1 785 EVS(IJ) = EV(I) + EV(J) 786 491 CONTINUE 787 490 CONTINUE 788c 789c Write out orbital energies for CC2-R12 model 790c 791 CALL GPOPEN(LU43,FCCR12E,'UNKNOWN',' ','FORMATTED', 792 & IDUM,LDUM) 793 WRITE(LU43,'(4E30.20)') (EVS(I), I=1,NSPAIR) 794 CALL GPCLOSE(LU43,'KEEP') 795C 796 IF (FRSTWR) THEN 797C 798C Write out V-matrix, etc. for CC2-R12 model (WK/UniKA/30-12-2003). 799C modified by C. Neiss: no not use singlet/triplet format any more, 800C make intermediates compatible with correlation factor r_12 (and 801C not 0.5*r_12 as in the MP2-R12 code)! 802C (April 2005) 803C 804 KSCR = 1 805 KEND1 = KSCR + NGAMMA(1) 806 IF (NGAMMA(1) .gt. LWORK) THEN 807 CALL QUIT('Insufficient WORK space in MAKEVR') 808 END IF 809 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,VS11,1) 810 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,VT11,1) 811 CALL DSCAL(NSPAIR*NSPAIR,1.0D0/3.0D0,VT11,1) 812 CALL CCR12PCK(WORK(KSCR),1,VS11,VT11,NRHF,NRHFA,NKILJ) 813 814 CALL GPOPEN(LU43,FCCR12V,'UNKNOWN',' ','UNFORMATTED', 815 & IDUM,LDUM) 816 1244 READ(LU43,END=1243) IDUM 817 READ(LU43) 818 GOTO 1244 819 1243 CONTINUE 820#ifdef VAR_MFDS 821! backspace if multifile dataset 822 BACKSPACE (LU43) 823#endif 824 WRITE(LU43) IANSAZ 825 WRITE(LU43) (WORK(KSCR-1+I), I=1, NKILJ(1)) 826 CALL GPCLOSE(LU43,'KEEP') 827chf 828c write(lupri,*)'Norm V out of MP2-R12:', 829c & ddot(nkilj(1),work(kscr),1,work(kscr),1) 830chf 831 ! undo scaling 832 CALL DSCAL(NSPAIR*NSPAIR,3.0D0,VT11,1) 833 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,VS11,1) 834 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,VT11,1) 835C 836C Write out X-matrix for CC2-R12 model (WK/UniKA/30-12-2003). 837C 838 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,RS11,1) 839 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,RT11,1) 840 CALL DSCAL(NSPAIR*NSPAIR,1.0D0/3.0D0,RT11,1) 841 CALL CCR12PCK(WORK(KSCR),1,RS11,RT11,NRHF,NRHF,NKILJ) 842 CALL GPOPEN(LU43,FCCR12X,'UNKNOWN',' ','UNFORMATTED', 843 & IDUM,LDUM) 844 IF (NOTONE) THEN 845 1246 READ(LU43,END=1245) IDUM 846 READ(LU43) 847 GOTO 1246 848 1245 CONTINUE 849#ifdef VAR_MFDS 850! backspace if multifile dataset 851 BACKSPACE (LU43) 852#endif 853 END IF 854 WRITE(LU43) IANSAZ 855 WRITE(LU43) (WORK(KSCR-1+I), I=1, NKILJ(1)) 856 CALL GPCLOSE(LU43,'KEEP') 857 ! undo scaling 858 CALL DSCAL(NSPAIR*NSPAIR,3.0D0,RT11,1) 859 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,RS11,1) 860 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,RT11,1) 861C 862C Open files for B and C (WK/UniKA/30-12-2003). 863C 864 CALL GPOPEN(LU44,FCCR12B,'UNKNOWN',' ','UNFORMATTED', 865 & IDUM,LDUM) 866 IF (NOTONE) THEN 867 1248 READ(LU44,END=1247) IDUM,JDUM 868 READ(LU44) 869 GOTO 1248 870 1247 CONTINUE 871#ifdef VAR_MFDS 872! backspace if multifile dataset 873 BACKSPACE (LU44) 874#endif 875 END IF 876 CALL GPOPEN(LU45,FCCR12D,'UNKNOWN',' ','UNFORMATTED', 877 & IDUM,LDUM) 878 IF (NOTONE) THEN 879 1250 READ(LU45,END=1249) IDUM,JDUM 880 READ(LU45) 881 GOTO 1250 882 1249 CONTINUE 883#ifdef VAR_MFDS 884! backspace if multifile dataset 885 BACKSPACE (LU45) 886#endif 887 END IF 888 FRSTWR = .FALSE. 889 END IF 890C 891 IF (DOB) THEN 892 LU43 = -43 893 IF (IANSAZ .EQ. 2) THEN 894 CALL GPOPEN(LU43,'QMATSP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 895 READ(LU43,'(4E30.20)') QS11 896 CALL GPCLOSE(LU43,'KEEP') 897 CALL GPOPEN(LU43,'QMATTP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 898 READ(LU43,'(4E30.20)') QT11 899 CALL GPCLOSE(LU43,'KEEP') 900 ELSE 901 CALL GPOPEN(LU43,'QMATSQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 902 READ(LU43,'(4E30.20)') QS11 903 CALL GPCLOSE(LU43,'KEEP') 904 CALL GPOPEN(LU43,'QMATTQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 905 READ(LU43,'(4E30.20)') QT11 906 CALL GPCLOSE(LU43,'KEEP') 907 END IF 908 DO 701 MN = 1, NAPAIR 909 DO 702 J = 1, NSPAIR 910 DO 703 I = 1, NSPAIR 911 BS11(I,J,MN) = BS11(I,J,MN) - QS11(I,J) 912 BT11(I,J,MN) = BT11(I,J,MN) - QT11(I,J) 913 703 CONTINUE 914 702 CONTINUE 915 701 CONTINUE 916 DO J = 1, NSPAIR 917 DO I = 1, NSPAIR 918 US11(I,J) = US11(I,J) - QS11(I,J) 919 UT11(I,J) = UT11(I,J) - QT11(I,J) 920 END DO 921 END DO 922 IF (DOBP) THEN 923 LU43 = -43 924 IF (IANSAZ .EQ. 2) THEN 925 CALL GPOPEN(LU43,'XMATSP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 926 READ(LU43,'(4E30.20)') QS11 927 CALL GPCLOSE(LU43,'KEEP') 928 CALL GPOPEN(LU43,'XMATTP','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 929 READ(LU43,'(4E30.20)') QT11 930 CALL GPCLOSE(LU43,'KEEP') 931 ELSE 932 CALL GPOPEN(LU43,'XMATSQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 933 READ(LU43,'(4E30.20)') QS11 934 CALL GPCLOSE(LU43,'KEEP') 935 CALL GPOPEN(LU43,'XMATTQ','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 936 READ(LU43,'(4E30.20)') QT11 937 CALL GPCLOSE(LU43,'KEEP') 938 END IF 939 DO 704 MN = 1, NAPAIR 940 DO 705 J = 1, NSPAIR 941 DO 706 I = 1, NSPAIR 942 BS11(I,J,MN) = BS11(I,J,MN) + QS11(I,J) 943 BT11(I,J,MN) = BT11(I,J,MN) + QT11(I,J) 944 706 CONTINUE 945 705 CONTINUE 946 704 CONTINUE 947 DO J = 1, NSPAIR 948 DO I = 1, NSPAIR 949 US11(I,J) = US11(I,J) + QS11(I,J) 950 UT11(I,J) = UT11(I,J) + QT11(I,J) 951 END DO 952 END DO 953 END IF 954 END IF 955 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,US11,1) 956 CALL DSCAL(NSPAIR*NSPAIR,4.0D0,UT11,1) 957 CALL DSCAL(NSPAIR*NSPAIR,1.0D0/3.0D0,UT11,1) 958 CALL CCR12PCK(WORK(KSCR),1,US11,UT11,NRHF,NRHF,NKILJ) 959 WRITE(LU44) IANSAZ, IPDD, IPCC 960 WRITE(LU44) (WORK(KSCR-1+I), I=1, NKILJ(1)) 961 ! undo scaling 962 CALL DSCAL(NSPAIR*NSPAIR,3.0D0,UT11,1) 963 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,US11,1) 964 CALL DSCAL(NSPAIR*NSPAIR,0.25D0,UT11,1) 965C 966 IF (IPDD .EQ. 0) THEN 967 WRITE(LUPRI,'(/A/)')' SECOND-ORDER R12/0 PAIR ENERGIES:' 968 INVAR = 0 969 ELSE IF (IPDD .EQ. 1 .OR. IPDD .EQ. 4 .OR. IPDD .EQ. 7 .OR. 970 & IPDD .EQ. 9 .OR. IPDD .EQ. 11 .OR. IPDD .EQ. 13) THEN 971 WRITE(LUPRI,'(/A/)') ' SECOND-ORDER NONINVARIANT R12/'//IPCC// 972 * ' PAIR ENERGIES:' 973 INVAR = 0 974 ELSE 975 WRITE(LUPRI,'(/A/)')' SECOND-ORDER R12/'//IPCC//' PAIR ENERGIES:' 976 INVAR = 1 977 END IF 978 CALL DZERO(FS,NSPAIR) 979 CALL DZERO(FT,NSPAIR) 980 CALL DZERO(CS11,NSPAIR*NSPAIR) 981 CALL DZERO(CT11,NSPAIR*NSPAIR) 982 F2S = D0 983 F2T = D0 984 IJ = 0 985 JI = 0 986 II = 0 987 DO 500 ISYM = 1, NSYM 988 DO 500 I = 1, NRHF(ISYM) 989 II = II + 1 990 JJ = 0 991 DO 501 JSYM = 1, NSYM 992 DO 501 J = 1, NRHF(JSYM) 993 JJ = JJ + 1 994 IF (JJ .GT. II) GOTO 501 995 JI = JI + 1 996 IF (I .GT. NRHFA(ISYM) .OR. 997 * J .GT. NRHFA(JSYM)) GOTO 501 998 IJ = IJ + 1 999 CNINV(IJ,1) = WS11(JI,IJ) 1000 CNINV(IJ,2) = BS11(JI,JI,IJ) 1001 CNINV(IJ,9) = -BS11(JI,JI,IJ) / RS11(JI,JI) 1002 IF (IPDD .EQ. 0) THEN 1003 CNINV(IJ,3) = 1D0 1004 CNINV(IJ,4) = 2D0*WS11(JI,IJ) - BS11(JI,JI,IJ) 1005 ELSE 1006 IF (CNINV(IJ,2) .LT. -VCLTHR) THEN 1007 CNINV(IJ,3) = CNINV(IJ,1) / CNINV(IJ,2) 1008 ELSE 1009 WRITE(LUPRI,'(/A,E15.8,I4/)') 1010 & ' @ WARNING: NEGATIVE DIAGONAL ELEMENT OF SINGLET'// 1011 & ' B MATRIX',CNINV(IJ,2),IJ 1012 CNINV(IJ,3) = D0 1013 END IF 1014 CNINV(IJ,4) = CNINV(IJ,1) * CNINV(IJ,3) 1015 END IF 1016 IF (II .NE. JJ) THEN 1017 CNINV(IJ,5) = WT11(JI,IJ) 1018 CNINV(IJ,6) = BT11(JI,JI,IJ) 1019 CNINV(IJ,10) = -BT11(JI,JI,IJ) / RT11(JI,JI) 1020 IF (IPDD .EQ. 0) THEN 1021 CNINV(IJ,7) = 0.5D0 1022 CNINV(IJ,8) = WT11(JI,IJ) - 0.25D0*BT11(JI,JI,IJ) 1023 ELSE 1024 IF (CNINV(IJ,6) .LT. -VCLTHR) THEN 1025 CNINV(IJ,7) = CNINV(IJ,5) / CNINV(IJ,6) 1026 ELSE 1027 WRITE(LUPRI,'(/A,E15.8,I4/)') 1028 & ' @ WARNING: NEGATIVE DIAGONAL ELEMENT OF TRIPLET'// 1029 & ' B MATRIX',CNINV(IJ,6),IJ 1030 CNINV(IJ,7) = D0 1031 END IF 1032 CNINV(IJ,8) = CNINV(IJ,5) * CNINV(IJ,7) 1033 END IF 1034 ELSE 1035 CNINV(IJ,5) = D0 1036 CNINV(IJ,6) = D0 1037 CNINV(IJ,7) = D0 1038 CNINV(IJ,8) = D0 1039 CNINV(IJ,10) = D0 1040 END IF 1041 FS(IJ) = CNINV(IJ,4) 1042 FT(IJ) = CNINV(IJ,8) 1043 IF (.NOT. DONINV) THEN 1044 IF (R12SVD) THEN 1045 CALL R12MP2(WS11(1,IJ),WT11(1,IJ),WS11(1,IJ),WT11(1,IJ), 1046 * RS11,RT11,EVS(JI),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 1047 * QS11,QT11,BS11(1,1,IJ),BT11(1,1,IJ),EVS,XPDD,QQ11,IJ, 1048 * WSMIN,WTMIN,DELTA,CS11(1,IJ),CT11(1,IJ)) 1049 ELSE IF (R12DIA) THEN 1050 CALL MP2R12(WS11(1,IJ),WT11(1,IJ),WS11(1,IJ),WT11(1,IJ), 1051 * RS11,RT11,EVS(JI),NOCDIM,NSPAIR,FS(IJ),FT(IJ),VS11,VT11, 1052 * QS11,QT11,BS11(1,1,IJ),BT11(1,1,IJ),EVS,XPDD,QQ11,IJ, 1053 * WSMIN,WTMIN,DELTA,CS11(1,IJ),CT11(1,IJ)) 1054 ELSE 1055 CALL QUIT('No solver selected (R12SVD or R12DIA)') 1056 END IF 1057 WRITE(LUPRI,'(I4,4F15.9,2G16.6)') IJ,ES(IJ),ES(IJ)+FS(IJ), 1058 * ET(IJ),ET(IJ)+FT(IJ), 1059 * WSMIN,WTMIN 1060 ELSE 1061 WRITE(LUPRI,'(I4,4F15.9,2G16.6)') IJ,ES(IJ),ES(IJ)+FS(IJ), 1062 * ET(IJ),ET(IJ)+FT(IJ) 1063 CS11(IJ,IJ) = CNINV(IJ,3) 1064 CT11(IJ,IJ) = CNINV(IJ,7) 1065 ENDIF 1066 F2S = F2S + FS(IJ) 1067 F2T = F2T + FT(IJ) 1068 501 CONTINUE 1069 500 CONTINUE 1070 WRITE(LUPRI,'(/4X,6F15.9)') E2S,E2S+F2S,E2T,E2T+F2T 1071 IF (INVAR .EQ. 0) THEN 1072 WRITE(LUPRI,'(/A,F15.9//A/)') 1073 * ' Noninvariant MP2-R12/'//IPCC//' correlation energy =', 1074 * E2S+E2T+F2S+F2T, 1075 * ' IJ V(IJ) U(IJ) C(IJ) '// 1076 * ' V(IJ) U(IJ) C(IJ) ' 1077 IJ = 0 1078 DO 600 I = 1, NACDIM 1079 DO 601 J = 1, I 1080 IJ = IJ + 1 1081 WRITE(LUPRI,'(I4,8E12.4)') IJ,(CNINV(IJ,K),K=1,3), 1082 * (CNINV(IJ,K),K=5,7), 1083 * CNINV(IJ,9),CNINV(IJ,10) 1084 601 CONTINUE 1085 600 CONTINUE 1086 ELSE 1087 WRITE(LUPRI,'(/A,F15.9 )') 1088 * ' MP2-R12/'//IPCC//' correlation energy =', 1089 * E2S+E2T+F2S+F2T 1090 END IF 1091c 1092c Write out amplitudes for CC2-R12 model (WK/UniKA/30-12-2003). 1093c 1094 IF (IPRINT .GT. 3) THEN 1095 CALL AROUND('Singlet <C> matrix') 1096 CALL OUTPUT(CS11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 1097 CALL AROUND('Triplet <C> matrix') 1098 CALL OUTPUT(CT11,1,NSPAIR,1,NAPAIR,NSPAIR,NSPAIR,1,LUPRI) 1099 END IF 1100c 1101 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,CS11,1) 1102 CALL DSCAL(NSPAIR*NSPAIR,0.5D0,CT11,1) 1103 CALL CCR12PCK(WORK(KSCR),1,CS11,CT11,NRHF,NRHFA,NKILJ) 1104 WRITE(LU45) IANSAZ,IPDD,IPCC 1105 WRITE(LU45) (WORK(KSCR-1+I), I=1, NKILJ(1)) 1106chf 1107c write(lupri,*)'Norm MP2-R12 amplitudes c: ', 1108c & ddot(nkilj(1),work(kscr),1,work(kscr),1) 1109chf 1110 ! undo scaling: 1111 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,CS11,1) 1112 CALL DSCAL(NSPAIR*NSPAIR,2.0D0,CT11,1) 1113 999 CONTINUE 1114 CALL GPCLOSE(LU44,'KEEP') 1115 CALL GPCLOSE(LU45,'KEEP') 1116 CALL GPCLOSE(LU46,'KEEP') 1117 RETURN 1118 END 1119C /* Deck ryr*/ 1120 SUBROUTINE RYR(R2AM,SF,TF,Q11,NOCDIM,NSPAIR) 1121C 1122C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 1123C 1124#include "implicit.h" 1125#include "priunit.h" 1126 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 1127 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0) 1128 DIMENSION R2AM(*), SF(*), TF(*), 1129 * ISB(8), Q11(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 1130 INTEGER IDUM 1131 LOGICAL LDUM 1132#include "ccorb.h" 1133#include "ccsdsym.h" 1134#include "ccsdinp.h" 1135#include "r12int.h" 1136 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1137 CALL DZERO(Q11,NOCDIM**4) 1138 ISB(1) = 0 1139 DO 100 ISYM = 2, NSYM 1140 NBASF = NORB1(ISYM-1) + NORB2(ISYM-1) 1141 NNBASF = NBASF * (NBASF + 1) / 2 1142 ISB(ISYM) = ISB(ISYM-1) + NNBASF 1143 100 CONTINUE 1144 NBASF = NORB1(NSYM) + NORB2(NSYM) 1145 NNBASF = ISB(NSYM) + NBASF * (NBASF + 1) / 2 1146 LUMULB = -34 1147 CALL GPOPEN(LUMULB,'AUXFCK','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1148 READ (LUMULB,'(4E30.20)') (SF(I), I = 1, NNBASF) 1149 CALL GPCLOSE(LUMULB,'KEEP') 1150 CALL GPOPEN(LUMULB,'AUXOVL','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1151 READ (LUMULB,'(4E30.20)') (TF(I), I = 1, NNBASF) 1152 CALL GPCLOSE(LUMULB,'KEEP') 1153 DO 200 ISYM = 1, NSYM 1154 IJ = ISB(ISYM) 1155 DO 210 I = 1, NORB1(ISYM) + NORB2(ISYM) 1156 DO 220 J = 1, I 1157 IJ = IJ + 1 1158 IF ((I .LE. NORB1(ISYM) .AND. J .GT. NORB1(ISYM)) .OR. 1159 * (J .LE. NORB1(ISYM) .AND. I .GT. NORB1(ISYM))) THEN 1160 SF(IJ) = - SF(IJ) 1161 TF(IJ) = - TF(IJ) 1162 END IF 1163 220 CONTINUE 1164 210 CONTINUE 1165 200 CONTINUE 1166 DO 301 ISYMIA = 1, NSYM 1167 ISYMJB = ISYMIA 1168 DO 302 ISYMI = 1, NSYM 1169 ISYMA = MULD2H(ISYMI,ISYMIA) 1170 ISYMC = ISYMA 1171 DO 303 ISYMJ = 1, NSYM 1172 ISYMB = MULD2H(ISYMJ,ISYMJB) 1173 ISYMD = ISYMB 1174 DO 304 ISYMKC = 1, NSYM 1175 ISYMLD = ISYMKC 1176 ISYMK = MULD2H(ISYMC,ISYMKC) 1177 ISYML = MULD2H(ISYMD,ISYMLD) 1178 DO 305 I = 1, NRHF(ISYMI) 1179 KOFFI = IRHF(ISYMI) + I 1180 DO 306 A = 1, NVIR(ISYMA) 1181 NAI = IT1AM(ISYMA,ISYMI) 1182 * + NVIR(ISYMA)*(I-1) + A 1183 DO 307 J = 1, NRHF(ISYMJ) 1184 KOFFJ = IRHF(ISYMJ) + J 1185 DO 308 B = 1, NVIR(ISYMB) 1186 NBJ = IT1AM(ISYMB,ISYMJ) 1187 * + NVIR(ISYMB)*(J-1) + B 1188 NAIBJ = IT2AM(ISYMIA,ISYMJB) 1189 * + INDEX(NAI,NBJ) 1190 DO 309 K = 1, NRHF(ISYMK) 1191 KOFFK = IRHF(ISYMK) + K 1192 DO 310 L = 1, NRHF(ISYML) 1193 KOFFL = IRHF(ISYML) + L 1194C 1195 IF (A .LE. NORB1(ISYMA) .AND. 1196 * B .LE. NORB1(ISYMB)) THEN 1197C 1198 DO 311 C = 1, NVIR(ISYMC) 1199 NCK = IT1AM(ISYMC,ISYMK) 1200 * + NVIR(ISYMC)*(K-1) + C 1201 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1202 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1203 DO 312 D = 1, NVIR(ISYMD) 1204 NDL = IT1AM(ISYMD,ISYML) 1205 * + NVIR(ISYMD)*(L-1) + D 1206 NCKDL = IT2AM(ISYMKC,ISYMLD) 1207 * + INDEX(NCK,NDL) 1208 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1209 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1210 XSX = XAC * SBD + SAC * XBD 1211 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1212 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1213 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1214 312 CONTINUE 1215 311 CONTINUE 1216C 1217 ELSE IF (A .LE. NORB1(ISYMA) .AND. 1218 * B .GT. NORB1(ISYMB)) THEN 1219C 1220 DO 313 C = 1, NVIR(ISYMC) 1221 NCK = IT1AM(ISYMC,ISYMK) 1222 * + NVIR(ISYMC)*(K-1) + C 1223 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1224 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1225 D = B 1226 NDL = IT1AM(ISYMD,ISYML) 1227 * + NVIR(ISYMD)*(L-1) + D 1228 NCKDL = IT2AM(ISYMKC,ISYMLD) 1229 * + INDEX(NCK,NDL) 1230 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1231 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1232 XSX = XAC * SBD + SAC * XBD 1233 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1234 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1235 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1236 DO 314 D = 1, NORB1(ISYMD) 1237 NDL = IT1AM(ISYMD,ISYML) 1238 * + NVIR(ISYMD)*(L-1) + D 1239 NCKDL = IT2AM(ISYMKC,ISYMLD) 1240 * + INDEX(NCK,NDL) 1241 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1242 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1243 XSX = XAC * SBD + SAC * XBD 1244 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1245 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1246 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1247 314 CONTINUE 1248 313 CONTINUE 1249C 1250 ELSE IF (A .GT. NORB1(ISYMA) .AND. 1251 * B .LE. NORB1(ISYMB)) THEN 1252C 1253 DO 315 D = 1, NVIR(ISYMD) 1254 NDL = IT1AM(ISYMD,ISYML) 1255 * + NVIR(ISYMD)*(L-1) + D 1256 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1257 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1258 C = A 1259 NCK = IT1AM(ISYMC,ISYMK) 1260 * + NVIR(ISYMC)*(K-1) + C 1261 NCKDL = IT2AM(ISYMKC,ISYMLD) 1262 * + INDEX(NCK,NDL) 1263 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1264 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1265 XSX = XAC * SBD + SAC * XBD 1266 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1267 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1268 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1269 DO 316 C = 1, NORB1(ISYMC) 1270 NCK = IT1AM(ISYMC,ISYMK) 1271 * + NVIR(ISYMC)*(K-1) + C 1272 NCKDL = IT2AM(ISYMKC,ISYMLD) 1273 * + INDEX(NCK,NDL) 1274 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1275 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1276 XSX = XAC * SBD + SAC * XBD 1277 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1278 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1279 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1280 316 CONTINUE 1281 315 CONTINUE 1282C 1283 ELSE IF (A .GT. NORB1(ISYMA) .AND. 1284 * B .GT. NORB1(ISYMB)) THEN 1285C 1286 DO 317 C = 1, NORB1(ISYMC) 1287 NCK = IT1AM(ISYMC,ISYMK) 1288 * + NVIR(ISYMC)*(K-1) + C 1289 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1290 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1291 DO 318 D = 1, NORB1(ISYMD) 1292 NDL = IT1AM(ISYMD,ISYML) 1293 * + NVIR(ISYMD)*(L-1) + D 1294 NCKDL = IT2AM(ISYMKC,ISYMLD) 1295 * + INDEX(NCK,NDL) 1296 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1297 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1298 XSX = XAC * SBD + SAC * XBD 1299 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1300 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1301 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1302 318 CONTINUE 1303 317 CONTINUE 1304 C = A 1305 NCK = IT1AM(ISYMC,ISYMK) 1306 * + NVIR(ISYMC)*(K-1) + C 1307 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1308 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1309 DO 319 D = 1, NORB1(ISYMD) 1310 NDL = IT1AM(ISYMD,ISYML) 1311 * + NVIR(ISYMD)*(L-1) + D 1312 NCKDL = IT2AM(ISYMKC,ISYMLD) 1313 * + INDEX(NCK,NDL) 1314 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1315 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1316 XSX = XAC * SBD + SAC * XBD 1317 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1318 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1319 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1320 319 CONTINUE 1321 D = B 1322 NDL = IT1AM(ISYMD,ISYML) 1323 * + NVIR(ISYMD)*(L-1) + D 1324 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1325 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1326 DO 320 C = 1, NORB1(ISYMC) 1327 NCK = IT1AM(ISYMC,ISYMK) 1328 * + NVIR(ISYMC)*(K-1) + C 1329 NCKDL = IT2AM(ISYMKC,ISYMLD) 1330 * + INDEX(NCK,NDL) 1331 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1332 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1333 XSX = XAC * SBD + SAC * XBD 1334 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1335 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1336 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1337 320 CONTINUE 1338 C = A 1339 D = B 1340 NCK = IT1AM(ISYMC,ISYMK) 1341 * + NVIR(ISYMC)*(K-1) + C 1342 NDL = IT1AM(ISYMD,ISYML) 1343 * + NVIR(ISYMD)*(L-1) + D 1344 NCKDL = IT2AM(ISYMKC,ISYMLD) 1345 * + INDEX(NCK,NDL) 1346 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1347 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1348 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1349 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1350 XSX = XAC * SBD + SAC * XBD 1351 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1352 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1353 * R2AM(NAIBJ) * XSX * R2AM(NCKDL) 1354C 1355 END IF 1356C 1357 310 CONTINUE 1358 309 CONTINUE 1359 308 CONTINUE 1360 307 CONTINUE 1361 306 CONTINUE 1362 305 CONTINUE 1363 304 CONTINUE 1364 303 CONTINUE 1365 302 CONTINUE 1366 301 CONTINUE 1367 IJKL = 0 1368 FF = D1 / SQRT(D2) 1369 DO 600 K = 1, NOCDIM 1370 DO 601 L = 1, K 1371 DO 602 I = 1, NOCDIM 1372 DO 603 J = 1, I 1373 IJKL = IJKL + 1 1374 SF(IJKL) = Q11(I,J,K,L) + Q11(I,J,L,K) 1375 TF(IJKL) = Q11(I,J,K,L) - Q11(I,J,L,K) 1376 IF (I .EQ. J) SF(IJKL) = FF * SF(IJKL) 1377 IF (K .EQ. L) SF(IJKL) = FF * SF(IJKL) 1378 SF(IJKL) = SF(IJKL) * DP25 1379 TF(IJKL) = TF(IJKL) * DP75 1380 603 CONTINUE 1381 602 CONTINUE 1382 601 CONTINUE 1383 600 CONTINUE 1384 IF (IPRINT .GT. 3) THEN 1385 CALL AROUND('Singlet <RXR> matrix') 1386 CALL OUTPUT(SF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1387 CALL AROUND('Triplet <RXR> matrix') 1388 CALL OUTPUT(TF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1389 END IF 1390 RETURN 1391 END 1392C /* Deck rzr*/ 1393 SUBROUTINE RZR(R2AM,V2AM,SF,TF,Q11,NOCDIM,NSPAIR) 1394C 1395C This subroutine is not used! 1396C 1397C It evaluates the sum 1398C 1399C Q(ijkl) = Sum_abcd (ia|r12|jb) * [X(ac)*S(bd) + S(ac)*X(db)] * (kc|r12|ld) 1400C 1401C The sum over a,b,c,d runs over both the primary and the secondary basis. 1402C 1403C On input, the array R2AM contains the integrals (ia|r12|jb). 1404C NOCDIM is the number of (nonfrozen) occupied orbitals and NSPAIR 1405C is the number of pairs of (nonfrozen) occupied orbitals. 1406C 1407C On output, the arrays SF and TF contain the matrices Q(ijkl) for 1408C singlet and triplet coupled pairs, respectively. 1409C 1410C V2AM (of length NT2AMX) and Q11 (of length NOCDIM**4) are used 1411C for scratch space. 1412C 1413C The one-particle matrices X (AUXFCK) and S (AUXOVL) are read from disk. 1414C These are the primary+secondary exchange and overlap matrices in the 1415C the orthonormal bases, respectively. It is assumed that these matrices 1416C are diagonal in the secondary-secondary block. This is tested for. 1417C 1418C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 1419C 1420#include "implicit.h" 1421#include "priunit.h" 1422 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 1423 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0) 1424 PARAMETER (THRDIA = 1D-9) 1425 DIMENSION R2AM(*), V2AM(*), SF(*), TF(*), 1426 * ISB(8), Q11(NOCDIM,NOCDIM,NOCDIM,NOCDIM) 1427 LOGICAL LDUM 1428 INTEGER IDUM 1429#include "ccorb.h" 1430#include "ccsdsym.h" 1431#include "ccsdinp.h" 1432#include "r12int.h" 1433 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1434 CALL DZERO(Q11,NOCDIM**4) 1435 ISB(1) = 0 1436 DO 100 ISYM = 2, NSYM 1437 NBASF = NORB1(ISYM-1) + NORB2(ISYM-1) 1438 NNBASF = NBASF * (NBASF + 1) / 2 1439 ISB(ISYM) = ISB(ISYM-1) + NNBASF 1440 100 CONTINUE 1441 NBASF = NORB1(NSYM) + NORB2(NSYM) 1442 NNBASF = ISB(NSYM) + NBASF * (NBASF + 1) / 2 1443 LUMULB = -34 1444 CALL GPOPEN(LUMULB,'AUXFCK','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1445 READ (LUMULB,'(4E30.20)') (SF(I), I = 1, NNBASF) 1446 CALL GPCLOSE(LUMULB,'KEEP') 1447 CALL GPOPEN(LUMULB,'AUXOVL','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1448 READ (LUMULB,'(4E30.20)') (TF(I), I = 1, NNBASF) 1449 CALL GPCLOSE(LUMULB,'KEEP') 1450 DO 200 ISYM = 1, NSYM 1451 IJ = ISB(ISYM) 1452 DO 210 I = 1, NORB1(ISYM) + NORB2(ISYM) 1453 DO 220 J = 1, I 1454 IJ = IJ + 1 1455 IF ((I .LE. NORB1(ISYM) .AND. J .GT. NORB1(ISYM)) .OR. 1456 * (J .LE. NORB1(ISYM) .AND. I .GT. NORB1(ISYM))) THEN 1457 SF(IJ) = - SF(IJ) 1458 TF(IJ) = - TF(IJ) 1459 ELSE IF ((I .GT. NORB1(ISYM) .AND. J .GT. NORB1(ISYM))) THEN 1460 IF (I .NE. J .AND. .NOT. (R12NOB .OR. NORXR)) THEN 1461 IF (ABS(SF(IJ)) .GT. THRDIA) THEN 1462 WRITE(LUPRI,'(/A/A,3I5,E20.10/)') 1463 * '@ WARNING : Exchange matrix not diagonal', 1464 * ' Nondiagonal element is :', 1465 * ISYM,I,J,SF(IJ) 1466 IF (ABS(SF(IJ)) .GT. SQRT(THRDIA)) 1467 * CALL QUIT('Exchange matrix not diagonal') 1468 END IF 1469 IF (ABS(TF(IJ)) .GT. THRDIA) THEN 1470 WRITE(LUPRI,'(/A/A,3I5,E20.10/)') 1471 * '@ WARNING : Overlap matrix not diagonal', 1472 * ' Nondiagonal element is :', 1473 * ISYM,I,J,TF(IJ) 1474 IF (ABS(TF(IJ)) .GT. SQRT(THRDIA)) 1475 * CALL QUIT('Overlap matrix not diagonal') 1476 END IF 1477 END IF 1478 END IF 1479 220 CONTINUE 1480 210 CONTINUE 1481 200 CONTINUE 1482C 1483 DO 301 ISYMKA = 1, NSYM 1484 ISYMLB = ISYMKA 1485 DO 302 ISYMK = 1, NSYM 1486 ISYMA = MULD2H(ISYMK,ISYMKA) 1487 ISYMC = ISYMA 1488 ISYMKC = ISYMKA 1489 DO 303 ISYML = 1, NSYM 1490 ISYMB = MULD2H(ISYML,ISYMLB) 1491 ISYMD = ISYMB 1492 ISYMLD = ISYMLB 1493 DO 304 K = 1, NRHF(ISYMK) 1494 DO 305 L = 1, NRHF(ISYML) 1495 DO 306 A = 1, NVIR(ISYMA) 1496 NAK = IT1AM(ISYMA,ISYMK) 1497 * + NVIR(ISYMA)*(K-1) + A 1498 DO 307 B = 1, NVIR(ISYMB) 1499 NBL = IT1AM(ISYMB,ISYML) 1500 * + NVIR(ISYMB)*(L-1) + B 1501 NAKBL = IT2AM(ISYMKA,ISYMLB) 1502 * + INDEX(NAK,NBL) 1503C 1504 V2AM(NAKBL) = D0 1505C 1506 IF (A .LE. NORB1(ISYMA) .AND. 1507 * B .LE. NORB1(ISYMB)) THEN 1508C 1509 DO 308 C = 1, NVIR(ISYMC) 1510 NCK = IT1AM(ISYMC,ISYMK) 1511 * + NVIR(ISYMC)*(K-1) + C 1512 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1513 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1514 DO 309 D = 1, NVIR(ISYMD) 1515 NDL = IT1AM(ISYMD,ISYML) 1516 * + NVIR(ISYMD)*(L-1) + D 1517 NCKDL = IT2AM(ISYMKC,ISYMLD) 1518 * + INDEX(NCK,NDL) 1519 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1520 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1521 XSX = XAC * SBD + SAC * XBD 1522 V2AM(NAKBL) = V2AM(NAKBL) 1523 * + XSX * R2AM(NCKDL) 1524 309 CONTINUE 1525 308 CONTINUE 1526C 1527 ELSE IF (A .LE. NORB1(ISYMA) .AND. 1528 * B .GT. NORB1(ISYMB)) THEN 1529C 1530 DO 310 C = 1, NVIR(ISYMC) 1531 NCK = IT1AM(ISYMC,ISYMK) 1532 * + NVIR(ISYMC)*(K-1) + C 1533 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1534 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1535 D = B 1536 NDL = NBL 1537 NCKDL = IT2AM(ISYMKC,ISYMLD) 1538 * + INDEX(NCK,NDL) 1539 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1540 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1541 XSX = XAC * SBD + SAC * XBD 1542 V2AM(NAKBL) = V2AM(NAKBL) 1543 * + XSX * R2AM(NCKDL) 1544 DO 311 D = 1, NORB1(ISYMD) 1545 NDL = IT1AM(ISYMD,ISYML) 1546 * + NVIR(ISYMD)*(L-1) + D 1547 NCKDL = IT2AM(ISYMKC,ISYMLD) 1548 * + INDEX(NCK,NDL) 1549 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1550 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1551 XSX = XAC * SBD + SAC * XBD 1552 V2AM(NAKBL) = V2AM(NAKBL) 1553 * + XSX * R2AM(NCKDL) 1554 311 CONTINUE 1555 310 CONTINUE 1556C 1557 ELSE IF (A .GT. NORB1(ISYMA) .AND. 1558 * B .LE. NORB1(ISYMB)) THEN 1559C 1560 DO 312 D = 1, NVIR(ISYMD) 1561 NDL = IT1AM(ISYMD,ISYML) 1562 * + NVIR(ISYMD)*(L-1) + D 1563 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1564 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1565 C = A 1566 NCK = NAK 1567 NCKDL = IT2AM(ISYMKC,ISYMLD) 1568 * + INDEX(NCK,NDL) 1569 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1570 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1571 XSX = XAC * SBD + SAC * XBD 1572 V2AM(NAKBL) = V2AM(NAKBL) 1573 * + XSX * R2AM(NCKDL) 1574 DO 313 C = 1, NORB1(ISYMC) 1575 NCK = IT1AM(ISYMC,ISYMK) 1576 * + NVIR(ISYMC)*(K-1) + C 1577 NCKDL = IT2AM(ISYMKC,ISYMLD) 1578 * + INDEX(NCK,NDL) 1579 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1580 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1581 XSX = XAC * SBD + SAC * XBD 1582 V2AM(NAKBL) = V2AM(NAKBL) 1583 * + XSX * R2AM(NCKDL) 1584 313 CONTINUE 1585 312 CONTINUE 1586C 1587 ELSE IF (A .GT. NORB1(ISYMA) .AND. 1588 * B .GT. NORB1(ISYMB)) THEN 1589C 1590 DO 314 C = 1, NORB1(ISYMC) 1591 NCK = IT1AM(ISYMC,ISYMK) 1592 * + NVIR(ISYMC)*(K-1) + C 1593 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1594 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1595 DO 315 D = 1, NORB1(ISYMD) 1596 NDL = IT1AM(ISYMD,ISYML) 1597 * + NVIR(ISYMD)*(L-1) + D 1598 NCKDL = IT2AM(ISYMKC,ISYMLD) 1599 * + INDEX(NCK,NDL) 1600 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1601 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1602 XSX = XAC * SBD + SAC * XBD 1603 V2AM(NAKBL) = V2AM(NAKBL) 1604 * + XSX * R2AM(NCKDL) 1605 315 CONTINUE 1606 314 CONTINUE 1607 C = A 1608 NCK = NAK 1609 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1610 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1611 DO 316 D = 1, NORB1(ISYMD) 1612 NDL = IT1AM(ISYMD,ISYML) 1613 * + NVIR(ISYMD)*(L-1) + D 1614 NCKDL = IT2AM(ISYMKC,ISYMLD) 1615 * + INDEX(NCK,NDL) 1616 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1617 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1618 XSX = XAC * SBD + SAC * XBD 1619 V2AM(NAKBL) = V2AM(NAKBL) 1620 * + XSX * R2AM(NCKDL) 1621 316 CONTINUE 1622 D = B 1623 NDL = NBL 1624 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1625 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1626 DO 317 C = 1, NORB1(ISYMC) 1627 NCK = IT1AM(ISYMC,ISYMK) 1628 * + NVIR(ISYMC)*(K-1) + C 1629 NCKDL = IT2AM(ISYMKC,ISYMLD) 1630 * + INDEX(NCK,NDL) 1631 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1632 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1633 XSX = XAC * SBD + SAC * XBD 1634 V2AM(NAKBL) = V2AM(NAKBL) 1635 * + XSX * R2AM(NCKDL) 1636 317 CONTINUE 1637 C = A 1638 D = B 1639 NCK = NAK 1640 NDL = NBL 1641 NCKDL = NAKBL 1642 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1643 SAC = TF(ISB(ISYMA)+INDEX(A,C)) 1644 SBD = TF(ISB(ISYMB)+INDEX(B,D)) 1645 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1646 XSX = XAC * SBD + SAC * XBD 1647 V2AM(NAKBL) = V2AM(NAKBL) 1648 * + XSX * R2AM(NCKDL) 1649C 1650 END IF 1651C 1652 307 CONTINUE 1653 306 CONTINUE 1654 305 CONTINUE 1655 304 CONTINUE 1656 303 CONTINUE 1657 302 CONTINUE 1658 301 CONTINUE 1659C 1660 DO 401 ISYMIA = 1, NSYM 1661 ISYMJB = ISYMIA 1662 DO 402 ISYMI = 1, NSYM 1663 ISYMA = MULD2H(ISYMI,ISYMIA) 1664 DO 403 ISYMJ = 1, NSYM 1665 ISYMB = MULD2H(ISYMJ,ISYMJB) 1666 DO 404 I = 1, NRHF(ISYMI) 1667 KOFFI = IRHF(ISYMI) + I 1668 DO 405 J = 1, NRHF(ISYMJ) 1669 KOFFJ = IRHF(ISYMJ) + J 1670 DO 406 A = 1, NVIR(ISYMA) 1671 ISYMJA = MULD2H(ISYMJ,ISYMA) 1672 DO 407 B = 1, NVIR(ISYMB) 1673 NAI = IT1AM(ISYMA,ISYMI) 1674 * + NVIR(ISYMA)*(I-1) + A 1675 NBJ = IT1AM(ISYMB,ISYMJ) 1676 * + NVIR(ISYMB)*(J-1) + B 1677 NAIBJ = IT2AM(ISYMIA,ISYMJB) 1678 * + INDEX(NAI,NBJ) 1679 DO 408 ISYMKA = 1, NSYM 1680 ISYMLB = ISYMKA 1681 ISYMK = MULD2H(ISYMA,ISYMKA) 1682 ISYML = MULD2H(ISYMB,ISYMLB) 1683 DO 409 K = 1, NRHF(ISYMK) 1684 KOFFK = IRHF(ISYMK) + K 1685 DO 410 L = 1, NRHF(ISYML) 1686 KOFFL = IRHF(ISYML) + L 1687 NAK = IT1AM(ISYMA,ISYMK) 1688 * + NVIR(ISYMA)*(K-1) + A 1689 NBL = IT1AM(ISYMB,ISYML) 1690 * + NVIR(ISYMB)*(L-1) + B 1691 NAKBL = IT2AM(ISYMKA,ISYMLB) 1692 * + INDEX(NAK,NBL) 1693 Q11(KOFFI,KOFFJ,KOFFK,KOFFL) = 1694 * Q11(KOFFI,KOFFJ,KOFFK,KOFFL) + 1695 * R2AM(NAIBJ) * V2AM(NAKBL) 1696 410 CONTINUE 1697 409 CONTINUE 1698 408 CONTINUE 1699 407 CONTINUE 1700 406 CONTINUE 1701 405 CONTINUE 1702 404 CONTINUE 1703 403 CONTINUE 1704 402 CONTINUE 1705 401 CONTINUE 1706 IJKL = 0 1707 FF = D1 / SQRT(D2) 1708 DO 600 K = 1, NOCDIM 1709 DO 601 L = 1, K 1710 DO 602 I = 1, NOCDIM 1711 DO 603 J = 1, I 1712 IJKL = IJKL + 1 1713 SF(IJKL) = Q11(I,J,K,L) + Q11(I,J,L,K) 1714 TF(IJKL) = Q11(I,J,K,L) - Q11(I,J,L,K) 1715 IF (I .EQ. J) SF(IJKL) = FF * SF(IJKL) 1716 IF (K .EQ. L) SF(IJKL) = FF * SF(IJKL) 1717 SF(IJKL) = SF(IJKL) * DP25 1718 TF(IJKL) = TF(IJKL) * DP75 1719 603 CONTINUE 1720 602 CONTINUE 1721 601 CONTINUE 1722 600 CONTINUE 1723C 1724 IF (IPRINT .GT. 3) THEN 1725 CALL AROUND('Singlet <RXR> matrix') 1726 CALL OUTPUT(SF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1727 CALL AROUND('Triplet <RXR> matrix') 1728 CALL OUTPUT(TF,1,NSPAIR,1,NSPAIR,NSPAIR,NSPAIR,1,LUPRI) 1729 WRITE(LUPRI,*) 1730 END IF 1731 RETURN 1732 END 1733C /* Deck comkr*/ 1734 SUBROUTINE COMKR(R2AM,V2AM,S2AM,T2AM,SF,LSF) 1735C 1736C This subroutine computes the integrals <pq|[r12,K1+K2]|ij>, 1737C utilizing the secondary basis RI approximation. 1738C 1739C Written by Wim Klopper (University of Karlsruhe, 31 October 2002). 1740C 1741#include "implicit.h" 1742#include "priunit.h" 1743 PARAMETER (D0 = 0D0, D1 = 1D0, D2 = 2D0, D3 = 3D0, D1P5 = 1.5D0, 1744 * DP5 = 0.5D0, DP25 = 0.25D0, DP75 = 0.75D0) 1745 PARAMETER (THRDIA = 1D-9) 1746 DIMENSION R2AM(*), T2AM(*), V2AM(*), S2AM(*), SF(*), ISB(8) 1747 REAL*8 RR, VV, XBD, XAC, DSM 1748 LOGICAL LDUM 1749 INTEGER IDUM 1750#include "ccorb.h" 1751#include "ccsdsym.h" 1752#include "ccsdinp.h" 1753#include "r12int.h" 1754 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1755 ISB(1) = 0 1756 DO 100 ISYM = 2, NSYM 1757 NBASF = NORB1(ISYM-1) + NORB2(ISYM-1) 1758 NNBASF = NBASF * (NBASF + 1) / 2 1759 ISB(ISYM) = ISB(ISYM-1) + NNBASF 1760 100 CONTINUE 1761 NBASF = NORB1(NSYM) + NORB2(NSYM) 1762 NNBASF = ISB(NSYM) + NBASF * (NBASF + 1) / 2 1763 IF (NNBASF .GT. LSF) 1764 * CALL QUIT('NOT ENOUGH SPACE IN COMKR') 1765 LUMULB = -34 1766 CALL GPOPEN(LUMULB,'AUXFCK','UNKNOWN',' ','FORMATTED',IDUM,LDUM) 1767 READ (LUMULB,'(4E30.20)') (SF(I), I = 1, NNBASF) 1768 CALL GPCLOSE(LUMULB,'KEEP') 1769 IF (IPRINT. GE. 10) CALL AROUND('Exchange matrices') 1770 DO 200 ISYM = 1, NSYM 1771 IF (IPRINT. GE. 10) THEN 1772 WRITE (LUPRI,'(A,I2)') ' Symmetry =', ISYM 1773 CALL OUTPAK(SF(ISB(ISYM)+1),NVIR(ISYM),1,LUPRI) 1774 END IF 1775 IF (.NOT. (R12NOB .OR. NORXR)) THEN 1776 DO 210 I = NORB1(ISYM) + 2, NVIR(ISYM) 1777 DO 220 J = NORB1(ISYM) + 1, I - 1 1778 IJ = ISB(ISYM) + INDEX(I,J) 1779 IF (ABS(SF(IJ)) .GT. THRDIA) THEN 1780 WRITE(LUPRI,'(/A/A,3I5,E20.10/)') 1781 * '@ WARNING : Exchange matrix not diagonal', 1782 * ' Nondiagonal element is :', 1783 * ISYM,I,J,SF(IJ) 1784 IF (ABS(SF(IJ)) .GT. SQRT(THRDIA)) 1785 * CALL QUIT('Exchange matrix not diagonal') 1786 END IF 1787 220 CONTINUE 1788 210 CONTINUE 1789 END IF 1790 200 CONTINUE 1791C 1792 IF (ONEAUX) THEN 1793 DO II = 1, NH2AMX 1794 T2AM(II) = V2AM(II) 1795 END DO 1796 DO 311 ISYMKA = 1, NSYM 1797 ISYMLB = ISYMKA 1798 LBOFF = NH1AM(ISYMLB) * (NH1AM(ISYMLB) + 1) / 2 1799 KAOFF = LBOFF 1800 DO 312 ISYMK = 1, NSYM 1801 ISYMA = MULD2H(ISYMK,ISYMKA) 1802 ISYMC = ISYMA 1803 ISYMKC = ISYMKA 1804 DO 313 ISYML = 1, NSYM 1805 ISYMB = MULD2H(ISYML,ISYMLB) 1806 ISYMD = ISYMB 1807 ISYMLD = ISYMLB 1808 DO 314 K = 1, NRHF(ISYMK) 1809 DO 315 L = 1, NRHF(ISYML) 1810 DO 316 A = 1, NORB1(ISYMA) 1811 IF (R12CBS) THEN 1812 NSTRC = 1 1813 ELSE 1814 NSTRC = NORB1(ISYMC) + 1 1815 END IF 1816 NENDC = NVIR(ISYMC) 1817 NAK = IH1AM(ISYMA,ISYMK) 1818 * + NORB1(ISYMA)*(K-1) + A 1819 DO 317 B = 1, NORB1(ISYMB) 1820 IF (R12CBS) THEN 1821 NSTRD = 1 1822 ELSE 1823 NSTRD = NORB1(ISYMD) + 1 1824 END IF 1825 NENDD = NVIR(ISYMD) 1826 NBL = IH1AM(ISYMB,ISYML) 1827 * + NORB1(ISYMB)*(L-1) + B 1828 IF (NBL .GT. NAK) GOTO 317 1829 NAKBL = IH2AM(ISYMKA,ISYMLB) 1830 * + INDEX(NAK,NBL) 1831 DSM = D0 1832 DO 318 C = NSTRC, NENDC 1833 IF (C .LE. NORB1(ISYMC)) THEN 1834 NCK = IH1AM(ISYMC,ISYMK) 1835 * + NORB1(ISYMC)*(K-1) + C 1836 NCKBL = IH2AM(ISYMKC,ISYMLB) 1837 * + INDEX(NCK,NBL) 1838 ELSE 1839 NCK = IG1AM(ISYMC,ISYMK) 1840 * + NORB2(ISYMC)*(K-1) 1841 * - NORB1(ISYMC) + C 1842 NCKBL = IH2AM(ISYMKC,ISYMLB) 1843 * + NH1AM(ISYMLB)*(NCK - 1) 1844 * + NBL + LBOFF 1845 END IF 1846 RR = R2AM(NCKBL) 1847 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1848 DSM = DSM + XAC * RR 1849 318 CONTINUE 1850 S2AM(NAKBL) = T2AM(NAKBL) - DSM 1851 DSM = D0 1852 DO 319 D = NSTRD, NENDD 1853 IF (D .LE. NORB1(ISYMD)) THEN 1854 NDL = IH1AM(ISYMD,ISYML) 1855 * + NORB1(ISYMD)*(L-1) + D 1856 NDLAK = IH2AM(ISYMLD,ISYMKA) 1857 * + INDEX(NDL,NAK) 1858 ELSE 1859 NDL = IG1AM(ISYMD,ISYML) 1860 * + NORB2(ISYMD)*(L-1) 1861 * - NORB1(ISYMD) + D 1862 NDLAK = IH2AM(ISYMLD,ISYMKA) 1863 * + NH1AM(ISYMKA)*(NDL - 1) 1864 * + NAK + KAOFF 1865 END IF 1866 RR = R2AM(NDLAK) 1867 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1868 DSM = DSM + XBD * RR 1869 319 CONTINUE 1870 S2AM(NAKBL) = S2AM(NAKBL) - DSM 1871 317 CONTINUE 1872 316 CONTINUE 1873 315 CONTINUE 1874 314 CONTINUE 1875 313 CONTINUE 1876 312 CONTINUE 1877 311 CONTINUE 1878 ELSE 1879 DO 401 ISYMIA = 1, NSYM 1880 ISYMJB = ISYMIA 1881 NAIBJ = IT2AM(ISYMIA,ISYMJB) 1882 KOFFU = IU2AM(ISYMIA,ISYMJB) 1883 NTOTAI = NT1AM(ISYMIA) 1884 DO 402 IA = 1, NTOTAI 1885 DO 403 JB = 1, IA 1886 NAIBJ = NAIBJ + 1 1887 T2AM(NAIBJ) = V2AM(KOFFU + (IA-1)*NTOTAI + JB) + 1888 & V2AM(KOFFU + (JB-1)*NTOTAI + IA) 1889 403 CONTINUE 1890 402 CONTINUE 1891 401 CONTINUE 1892 DO 301 ISYMKA = 1, NSYM 1893 ISYMLB = ISYMKA 1894 DO 302 ISYMK = 1, NSYM 1895 ISYMA = MULD2H(ISYMK,ISYMKA) 1896 ISYMC = ISYMA 1897 ISYMKC = ISYMKA 1898 DO 303 ISYML = 1, NSYM 1899 ISYMB = MULD2H(ISYML,ISYMLB) 1900 ISYMD = ISYMB 1901 ISYMLD = ISYMLB 1902 DO 304 K = 1, NRHF(ISYMK) 1903 DO 305 L = 1, NRHF(ISYML) 1904 DO 306 A = 1, NORB1(ISYMA) 1905 IF (R12CBS) THEN 1906 NSTRC = 1 1907 ELSE 1908 NSTRC = NORB1(ISYMC) + 1 1909 END IF 1910 NENDC = NVIR(ISYMC) 1911 NAK = IT1AM(ISYMA,ISYMK) 1912 * + NVIR(ISYMA)*(K-1) + A 1913 DO 307 B = 1, NORB1(ISYMB) 1914 IF (R12CBS) THEN 1915 NSTRD = 1 1916 ELSE 1917 NSTRD = NORB1(ISYMD) + 1 1918 END IF 1919 NENDD = NVIR(ISYMD) 1920 NBL = IT1AM(ISYMB,ISYML) 1921 * + NVIR(ISYMB)*(L-1) + B 1922 IF (NBL .GT. NAK) GOTO 307 1923 NAKBL = IT2AM(ISYMKA,ISYMLB) 1924 * + INDEX(NAK,NBL) 1925 DSM = D0 1926 DO 308 C = NSTRC, NENDC 1927 NCK = IT1AM(ISYMC,ISYMK) 1928 * + NVIR(ISYMC)*(K-1) + C 1929 NCKBL = IT2AM(ISYMKC,ISYMLB) 1930 * + INDEX(NCK,NBL) 1931 RR = R2AM(NCKBL) 1932 XAC = SF(ISB(ISYMA)+INDEX(A,C)) 1933 DSM = DSM + XAC * RR 1934 308 CONTINUE 1935 S2AM(NAKBL) = T2AM(NAKBL) - DSM 1936 DSM = D0 1937 DO 309 D = NSTRD, NENDD 1938 NDL = IT1AM(ISYMD,ISYML) 1939 * + NVIR(ISYMD)*(L-1) + D 1940 NAKDL = IT2AM(ISYMKA,ISYMLD) 1941 * + INDEX(NAK,NDL) 1942 RR = R2AM(NAKDL) 1943 XBD = SF(ISB(ISYMB)+INDEX(B,D)) 1944 DSM = DSM + XBD * RR 1945 309 CONTINUE 1946 S2AM(NAKBL) = S2AM(NAKBL) - DSM 1947 307 CONTINUE 1948 306 CONTINUE 1949 305 CONTINUE 1950 304 CONTINUE 1951 303 CONTINUE 1952 302 CONTINUE 1953 301 CONTINUE 1954 END IF 1955C 1956 IF (IPRINT .GE. 10) THEN 1957 CALL AROUND('<pq|r12*(K1+K2)|ij> integrals') 1958 IF (ONEAUX) THEN 1959 DO ISYM = 1, NSYM 1960 WRITE(LUPRI,'(A,I2)') ' Symmetry =',ISYM 1961 NTOT = NH1AM(ISYM) 1962 KOFF = IH2AM(ISYM,ISYM) + 1 1963 CALL OUTPAK(T2AM(KOFF),NTOT,1,LUPRI) 1964 END DO 1965 ELSE 1966 DO ISYM = 1, NSYM 1967 WRITE(LUPRI,'(A,I2)') ' Symmetry =',ISYM 1968 NTOT = NT1AM(ISYM) 1969 KOFF = IT2AM(ISYM,ISYM) + 1 1970 CALL OUTPAK(T2AM(KOFF),NTOT,1,LUPRI) 1971 END DO 1972 END IF 1973 CALL AROUND('-<pq|(K1+K2)*r12|ij> integrals') 1974 IF (ONEAUX) THEN 1975 DO ISYM = 1, NSYM 1976 WRITE(LUPRI,'(A,I2)') ' Symmetry =',ISYM 1977 NTOT = NH1AM(ISYM) 1978 KOFF = IH2AM(ISYM,ISYM) + 1 1979 CALL OUTPAK(S2AM(KOFF),NTOT,1,LUPRI) 1980 END DO 1981 ELSE 1982 DO ISYM = 1, NSYM 1983 WRITE(LUPRI,'(A,I2)') ' Symmetry =',ISYM 1984 NTOT = NT1AM(ISYM) 1985 KOFF = IT2AM(ISYM,ISYM) + 1 1986 CALL OUTPAK(S2AM(KOFF),NTOT,1,LUPRI) 1987 END DO 1988 END IF 1989 END IF 1990 RETURN 1991 END 1992