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 19*=====================================================================* 20 SUBROUTINE CCSDT_ADEN_NODDY(LISTL,IDLSTL,LISTR,IDLSTR, 21 * XLAMDP0,XLAMDH0,FOCK0, 22 * DIJ,DAB,DO_DIA,DIA,DO_XMMAT,XMMAT, 23 * WORK,LWORK, 24 * LUDELD,FNDELD,LUCKJD,FNCKJD,LUDKBC, 25 * FNDKBC,LUTOC,FNTOC,LU3VI,FN3VI, 26 * LUDKBC3,FNDKBC3,LU3FOPX,FN3FOPX, 27 * LU3FOP2X,FN3FOP2X) 28*---------------------------------------------------------------------* 29* Purpose: compute density corresponding to the triples * 30* contributions to a L A{O} R contraction * 31* * 32* Written by Christof Haettig, Jan 2003, Friedrichstal * 33*---------------------------------------------------------------------* 34 IMPLICIT NONE 35#include "priunit.h" 36#include "ccorb.h" 37#include "ccsdsym.h" 38#include "dummy.h" 39#include "ccr2rsp.h" 40#include "ccr1rsp.h" 41#include "ccer1rsp.h" 42#include "ccnoddy.h" 43 44 INTEGER ISYM0 45 PARAMETER (ISYM0 = 1) 46 47 LOGICAL LOCDBG, HAVETB3T3, PRINTDENS 48 PARAMETER (LOCDBG=.FALSE., HAVETB3T3=.TRUE., PRINTDENS=.FALSE.) 49 50 CHARACTER*3 LISTL, LISTR 51 CHARACTER*(*) FNDELD, FNCKJD, FNDKBC, FNTOC, FN3VI, FNDKBC3 52 CHARACTER*(*) FN3FOPX, FN3FOP2X 53 LOGICAL DO_DIA, DO_XMMAT 54 INTEGER LUDELD,LUCKJD,LUDKBC,LUTOC,LU3VI,LUDKBC3,LU3FOPX,LU3FOP2X 55 INTEGER LWORK, IDLSTL, IDLSTR 56 57#if defined (SYS_CRAY) 58 REAL XLAMDP0(*), XLAMDH0(*), FOCK0(*) 59 REAL DIJ(*), DAB(*), DIA(*), XMMAT(*) 60 REAL WORK(LWORK) 61 REAL HALF, ONE, ZERO, TWO, THREE, DDOT, XNORM 62 REAL FREQLST, FREQY, FREQX, FREQR 63#else 64 DOUBLE PRECISION XLAMDP0(*), XLAMDH0(*), FOCK0(*) 65 DOUBLE PRECISION DIJ(*), DAB(*), DIA(*), XMMAT(*) 66 DOUBLE PRECISION WORK(LWORK) 67 DOUBLE PRECISION HALF, ONE, TWO, ZERO, THREE, DDOT, XNORM 68 DOUBLE PRECISION FREQLST, FREQY, FREQX, FREQR, FREQL 69#endif 70 PARAMETER( ZERO=0.0D0, HALF=0.5D0, ONE=1.0D0, TWO = 2.0D0, 71 & THREE=3.0D0 ) 72 73 CHARACTER MODEL*10, LISTX*3, LISTY*3, LABELX*8, LABELY*8 74 LOGICAL CUBIC, LORXX, LORXY, HAVEX, HAVEY, HAVEXY, 75 & SPLIT_WBTB, SPLIT_WT 76 INTEGER KEND1, LWRK1, KWBAR, KWYINT, KTHETAY, KTBAR3, KTAMP3, 77 & NCK, NBJ, NAI, KAIBJCK, KCKAIBJ, KBJCKAI, IOFF, ISYMR, 78 & KWTILDE, KWBREVE, KTHEOCC, KTVIRT1, KTVIRT2, KEND2, 79 & LWRK2, KTHVIRT, KSWAP, KDABR, KDIJR, KDIAR, IOPT, 80 & KT0AMP, KTAMP2, KT0FMI, KTYFMI, KT0MI, KTYMI, 81 & KWDE, KWED, KSCR, KFNEMDL, KEMDLFN, KT3AMP0, KDLEI, 82 & NAN, NFN, NFI, NEM, NDL, KDLEMFN, KDLAN, KEMFI, NEI, 83 & KL2AMP, KDLFIAN, KDLFNAI, KDLFN, ISYML, 84 & KDAB0, KDAB0R, KDIJ0, KDIJ0R, KTAMP1, KXMMAT, KIMFN, 85 & KTHBAR, KWXINT, KTHETAX, KWXYVINT, KWXYOINT, KXPROP, 86 & KYPROP, KXYMAT, KFOCKD, KTHEOV, KTVIRT3, KSCRVVVO, 87 & ISYMX, ISYMY, IDLSTX, IDLSTY, LUSIFC, LUTEMP, 88 & KSCR1, KINT1T0, KINT2T0, KXIAJB, KYIAJB, 89 & KLAMPB, KLAMHB, KFOCKB, 90 & KINT1S0, KINT2S0, KINT1SB, KINT2SB, KDUM 91 92 INTEGER ILSTSYM, IR1TAMP 93 94*---------------------------------------------------------------------* 95* Begin: 96*---------------------------------------------------------------------* 97 CALL QENTER('CCTADNOD') 98 99 IF (NSYM.GT.1) CALL QUIT('No symmetry in CCSDT_ADEN_NODDY!') 100 101 IF (LISTR(1:3).EQ.'R1 ') THEN 102 ! we compute the first-order amplitude response vectors 103 CUBIC = .FALSE. 104 LORXX = LORXLRT(IDLSTR) 105 LORXY = .FALSE. 106 107 ELSE IF (LISTR(1:3).EQ.'RE ') THEN 108 ! we compute eigenvectors 109 CUBIC = .FALSE. 110 LORXX = .FALSE. 111 LORXY = .FALSE. 112 ELSE IF (LISTR(1:3).EQ.'R2 ') THEN 113 ! we compute the second-order amplitude response vectors 114 CUBIC = .TRUE. 115 116 LISTX = 'R1 ' 117 LABELX = LBLR2T(IDLSTR,1) 118 ISYMX = ISYR2T(IDLSTR,1) 119 FREQX = FRQR2T(IDLSTR,1) 120 LORXX = LORXR2T(IDLSTR,1) 121 IDLSTX = IR1TAMP(LABELX,LORXX,FREQX,ISYMX) 122 123 LISTY = 'R1 ' 124 LABELY = LBLR2T(IDLSTR,2) 125 ISYMY = ISYR2T(IDLSTR,2) 126 FREQY = FRQR2T(IDLSTR,2) 127 LORXY = LORXR2T(IDLSTR,2) 128 IDLSTY = IR1TAMP(LABELY,LORXY,FREQY,ISYMY) 129 130 ELSE IF (LISTR(1:3).EQ.'ER1') THEN 131 ! we compute first-order response of excited states 132 IF (.FALSE.) THEN 133 ! use quadratic response code 134 ! (for internal check of noddy code) 135 CUBIC = .FALSE. 136 LORXX = LORXER1(IDLSTR) 137 LORXY = .FALSE. 138 ELSE 139 ! use cubic response code 140 ! (this will be most similar to the real code) 141 CUBIC = .TRUE. 142 143 LISTX = 'R1 ' 144 LABELX = LBLER1(IDLSTR) 145 ISYMX = ISYOER1(IDLSTR) 146 FREQX = FRQER1(IDLSTR) 147 LORXX = LORXER1(IDLSTR) 148 IDLSTX = IR1TAMP(LABELX,LORXX,FREQX,ISYMX) 149 150 LISTY = 'RE ' 151 ISYMY = ISYSER1(IDLSTR) 152 FREQY = EIGER1(IDLSTR) 153 LORXY = .FALSE. 154 IDLSTY = ISTER1(IDLSTR) 155 END IF 156 ELSE 157 PRINT *,'LISTR:',LISTR 158 CALL QUIT('Unknown LISTR in CCSDT_ADEN_NODDY.') 159 END IF 160 161 IF (LORXX.OR.LORXY) THEN 162 CALL QUIT('Orbital relaxation not allowed in CCSDT_ADEN_NODDY.') 163 END IF 164 165 HAVEX = CUBIC .AND. (LISTX(1:3).EQ.'R1 ') 166 HAVEY = CUBIC .AND. (LISTY(1:3).EQ.'R1 ') 167 HAVEXY = HAVEX .OR. HAVEY 168 169 SPLIT_WBTB = CUBIC 170 SPLIT_WT = (LISTR(1:3).NE.'RE ') 171 172 ISYML = ILSTSYM(LISTL,IDLSTL) 173 ISYMR = ILSTSYM(LISTR,IDLSTR) 174 FREQR = FREQLST(LISTR,IDLSTR) 175 176*---------------------------------------------------------------------* 177* Allocate arrays kept until the end of the routine: 178*---------------------------------------------------------------------* 179 KEND1 = 1 180 181 KDAB0 = KEND1 182 KDIJ0 = KDAB0 + NVIRT*NVIRT 183 KEND1 = KDIJ0 + NRHFT*NRHFT 184 185 KFOCKD = KEND1 186 KEND1 = KFOCKD + NORBT 187 188 KWBAR = KEND1 189 KEND1 = KWBAR + NT1AMX*NT1AMX*NT1AMX 190 191 IF (SPLIT_WBTB) THEN 192 KTHBAR = KEND1 193 KEND1 = KTHBAR + NT1AMX*NT1AMX*NT1AMX 194 END IF 195 196 KWYINT = KEND1 197 KEND1 = KWYINT + NT1AMX*NT1AMX*NT1AMX 198 199 IF (SPLIT_WT) THEN 200 KTHETAY = KWYINT + NT1AMX*NT1AMX*NT1AMX 201 KEND1 = KTHETAY + NT1AMX*NT1AMX*NT1AMX 202 END IF 203 204 IF (CUBIC) THEN 205 KWXINT = KEND1 206 KTHETAX = KWXINT + NT1AMX*NT1AMX*NT1AMX 207 KEND1 = KTHETAX + NT1AMX*NT1AMX*NT1AMX 208 209 KWXYVINT = KEND1 210 KWXYOINT = KWXYVINT + NT1AMX*NT1AMX*NT1AMX 211 KEND1 = KWXYOINT + NT1AMX*NT1AMX*NT1AMX 212 213 KXPROP = KEND1 214 KYPROP = KXPROP + NBAST*NBAST 215 KXYMAT = KYPROP + NBAST*NBAST 216 KEND1 = KXYMAT + NBAST*NBAST 217 ELSE 218 KWXYVINT = KWYINT 219 KWXYOINT = KWYINT 220 END IF 221 222 KT3AMP0 = KEND1 223 KEND1 = KT3AMP0 + NT1AMX*NT1AMX*NT1AMX 224 225 KL2AMP = KEND1 226 KEND1 = KL2AMP + NT1AMX*NT1AMX 227 228 KTAMP1 = KEND1 229 KEND1 = KTAMP1 + NT1AMX 230 231 KT0AMP = KEND1 232 KTAMP2 = KT0AMP + NT1AMX*NT1AMX 233 KEND1 = KTAMP2 + NT1AMX*NT1AMX 234 235 IF (HAVETB3T3) THEN 236 KTBAR3 = KEND1 237 KTAMP3 = KTBAR3 + NT1AMX*NT1AMX*NT1AMX 238 KDABR = KTAMP3 + NT1AMX*NT1AMX*NT1AMX 239 KDIJR = KDABR + NVIRT*NVIRT 240 KDIAR = KDIJR + NRHFT*NRHFT 241 KDAB0R = KDIAR + NVIRT*NRHFT 242 KDIJ0R = KDAB0R + NVIRT*NVIRT 243 KEND1 = KDIJ0R + NRHFT*NRHFT 244 245 IF (DO_XMMAT) THEN 246 KXMMAT = KEND1 247 KEND1 = KXMMAT + NRHFT*NRHFT*NT1AMX 248 END IF 249 END IF 250 251 LWRK1 = LWORK - KEND1 252 IF (LWRK1.LT.0)CALL QUIT('Out of memory in CCSDT_ADEN_NODDY (1)') 253 254 IF (NSYM.GT.1) CALL QUIT('No symmetry in CCSDT_ADEN_NODDY!') 255 256*---------------------------------------------------------------------* 257* Read SCF orbital energies from file: 258*---------------------------------------------------------------------* 259 LUSIFC = -1 260 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 261 & .FALSE.) 262 REWIND LUSIFC 263 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 264 READ (LUSIFC) 265 READ (LUSIFC) (WORK(KFOCKD-1+I), I=1,NORBT) 266 CALL GPCLOSE(LUSIFC,'KEEP') 267 268*---------------------------------------------------------------------* 269* read doubles amplitudes (zeroth-order and response) from file: 270*---------------------------------------------------------------------* 271 IF (LWRK1.LT.NT2AMX) 272 & CALL QUIT('Out of memory in CCSDT_ADEN_NODDY (3a)') 273 274 ! read T^0 doubles amplitudes from file and square up 275 IOPT = 2 276 Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KEND1)) 277 CALL CC_T2SQ(WORK(KEND1),WORK(KT0AMP),ISYM0) 278 279 ISYMR = ILSTSYM(LISTR,IDLSTR) 280 281 ! read T^X doubles amplitudes from file and square up 282 IOPT = 3 283 Call CC_RDRSP(LISTR,IDLSTR,ISYMR,IOPT,MODEL, 284 & WORK(KTAMP1),WORK(KEND1)) 285 Call CCLR_DIASCL(WORK(KEND1),TWO,ISYMR) 286 CALL CC_T2SQ(WORK(KEND1),WORK(KTAMP2),ISYMR) 287 288 ISYML = ILSTSYM(LISTL,IDLSTL) 289 290 ! read left doubles amplitudes from file and square up 291 IOPT = 2 292 Call CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL,DUMMY,WORK(KEND1)) 293 CALL CC_T2SQ(WORK(KEND1),WORK(KL2AMP),ISYML) 294 295*---------------------------------------------------------------------* 296* Get W-bar (and Theta-bar) intermediate into memory: 297*---------------------------------------------------------------------* 298 IF (LISTL.EQ.'L0 ') THEN 299 300 LUTEMP = -1 301 CALL GPOPEN(LUTEMP,FILNODL30,'UNKNOWN',' ','UNFORMATTED', 302 & IDUMMY,.FALSE.) 303 READ(LUTEMP) (WORK(KWBAR+I-1), I=1,NT1AMX*NT1AMX*NT1AMX) 304 CALL GPCLOSE(LUTEMP,'KEEP') 305 306 IF (HAVETB3T3) THEN 307 CALL DCOPY(NT1AMX*NT1AMX*NT1AMX,WORK(KWBAR),1,WORK(KTBAR3),1) 308 END IF 309 310 ! cheat the code by putting 1/3 Tbar^(0) into W-bar 311 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,ONE/THREE,WORK(KWBAR),1) 312 313 IF (SPLIT_WBTB) THEN 314 ! Theta-bar is set to zero for this case 315 CALL DZERO(WORK(KTHBAR),NT1AMX*NT1AMX*NT1AMX) 316 END IF 317 318 ELSE IF (LISTL.EQ.'L1 ') THEN 319 320 CALL CCSDT_GWBIC(LISTL,IDLSTL,WORK(KWBAR),WORK(KTHBAR), 321 & SPLIT_WBTB,XLAMDP0,XLAMDH0,FOCK0, 322 & LUTOC,FNTOC,LU3VI,FN3VI,LUDKBC3,FNDKBC3, 323 & LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X, 324 & WORK(KEND1),LWRK1) 325 326 ! turn sign of Wbar intermediate 327 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,WORK(KWBAR),1) 328 329 IF (SPLIT_WBTB) THEN 330 ! turn also sign of Theta-bar intermediate 331 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,WORK(KTHBAR),1) 332 333 ! add THBAR intermediate to WBAR to obtain the full WBAR 334 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,ONE,WORK(KTHBAR),1, 335 & WORK(KWBAR),1) 336 END IF 337 338 IF (LOCDBG) THEN 339 XNORM = DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KWBAR),1,WORK(KWBAR),1) 340 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> Wbar^Y:',XNORM 341 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> LISTL,IDLSTL:',LISTL,IDLSTL 342 CALL PRINT_PT3_NODDY(WORK(KWBAR)) 343 END IF 344 345 IF (HAVETB3T3) THEN 346 ! accumulate Tbar^Y amplitudes in WORK and print 347 DO NCK = 1, NT1AMX 348 DO NBJ = 1, NT1AMX 349 DO NAI = 1, NT1AMX 350 KAIBJCK = ((NCK-1)*NT1AMX+NBJ-1)*NT1AMX+NAI 351 KCKAIBJ = ((NBJ-1)*NT1AMX+NAI-1)*NT1AMX+NCK 352 KBJCKAI = ((NAI-1)*NT1AMX+NCK-1)*NT1AMX+NBJ 353 WORK(KTBAR3-1+KAIBJCK) = WORK(KWBAR-1+KAIBJCK) + 354 & WORK(KWBAR-1+KCKAIBJ) + WORK(KWBAR-1+KBJCKAI) 355 END DO 356 END DO 357 END DO 358 END IF 359 360 ELSEIF (LISTL(1:3).EQ.'LE '.OR.LISTL(1:3).EQ.'M1 '.OR. 361 & LISTL(1:3).EQ.'E0 '.OR.LISTL(1:3).EQ.'N2 ' ) THEN 362 363 KSCR1 = KEND1 364 KEND2 = KSCR1 + NT1AMX 365 366 KINT1T0 = KEND2 367 KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT 368 KEND2 = KINT2T0 + NRHFT*NRHFT*NT1AMX 369 370 KXIAJB = KEND2 371 KYIAJB = KXIAJB + NT1AMX*NT1AMX 372 KEND2 = KYIAJB + NT1AMX*NT1AMX 373 374 LWRK2 = LWORK - KEND2 375 IF (LWRK2 .LT. 0) 376 & CALL QUIT('Out of memory in CCSDT_ADEN_NODDY (2)') 377 378C --------------------------------------- 379C Read precalculated integrals from file: 380C XINT1T0 = (KC|BD) 381C XINT2T0 = (KC|LJ) 382C XIAJB = 2(IA|JB) - (IB|JA) 383C YIAJB = (IA|JB) 384C --------------------------------------- 385 CALL CCSDT_READ_NODDY(.FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 386 & .TRUE.,WORK(KXIAJB),WORK(KYIAJB), 387 & .FALSE.,DUMMY,DUMMY, 388 & .TRUE.,WORK(KINT1T0),WORK(KINT2T0), 389 & .FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 390 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 391 392 CALL DZERO(WORK(KWBAR),NT1AMX*NT1AMX*NT1AMX) 393 CALL CCSDT_TBAR31_NODDY(WORK(KWBAR),FREQL,LISTL,IDLSTL, 394 & XLAMDP0,XLAMDH0,FOCK0, 395 & WORK(KFOCKD),WORK(KSCR1), 396 & WORK(KXIAJB),WORK(KINT1T0), 397 & WORK(KINT2T0),WORK(KEND2),LWRK2) 398 399 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KWBAR),1) 400 401 IF (HAVETB3T3) THEN 402 CALL DCOPY(NT1AMX*NT1AMX*NT1AMX,WORK(KWBAR),1,WORK(KTBAR3),1) 403 END IF 404 405 IF (SPLIT_WBTB) THEN 406 ! Theta-bar is set to zero for this case 407 CALL DZERO(WORK(KTHBAR),NT1AMX*NT1AMX*NT1AMX) 408 END IF 409 410 ! account for the fact that in the real code a combination 411 ! of three W intermediates give the complete Tbar3 412 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,ONE/THREE,WORK(KWBAR),1) 413 414 ELSE 415 CALL QUIT('Unknown or illegal list in CCDOTRSP_NODDY.') 416 END IF 417 418 IF (LOCDBG .AND. HAVETB3T3) THEN 419 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> Tbar^Y:' 420 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> LISTL,IDLSTL:',LISTL,IDLSTL 421 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> norm^2:', 422 & DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KTBAR3),1,WORK(KTBAR3),1) 423 CALL PRINT_PT3_NODDY(WORK(KTBAR3)) 424 END IF 425 426*---------------------------------------------------------------------* 427* Compute intermediates needed to represent the first- 428* or second-order response amplitudes: 429*---------------------------------------------------------------------* 430 IF (.NOT. CUBIC) THEN 431 432 IF (LISTR.EQ.'R1 ') THEN 433 ! Get w^y and theta^y intermediate into memory: 434 ! (returns in addition the zeroth-order amplitudes in T3AMP0, 435 ! and for double checking T3^Y on TAMP3) 436 CALL CCSDT_GWTIC(LISTR,IDLSTR,WORK(KWYINT),WORK(KTHETAY), 437 & WORK(KT3AMP0),HAVETB3T3,WORK(KTAMP3),LOCDBG, 438 & XLAMDP0,XLAMDH0,FOCK0, 439 & LUDELD,FNDELD,LUDKBC,FNDKBC,LUCKJD,FNCKJD, 440 & WORK(KEND1),LWRK1) 441 442 ! turn sign of t3, w and theta intermediates 443 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,WORK(KT3AMP0),1) 444 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,WORK(KWYINT),1) 445 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,WORK(KTHETAY),1) 446 IF (HAVETB3T3) 447 & CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,WORK(KTAMP3),1) 448 449 ELSE 450 451 KEND2 = KEND1 452 453 KINT1S0 = KEND2 454 KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT 455 KEND2 = KINT2S0 + NRHFT*NRHFT*NT1AMX 456 457 IF (LISTR.EQ.'RE ') THEN 458 KLAMPB = KEND2 459 KLAMHB = KLAMPB + NLAMDT 460 KFOCKB = KLAMHB + NLAMDT 461 KEND2 = KFOCKB + N2BST(1) 462 463 KINT1SB = KEND2 464 KINT2SB = KINT1SB + NT1AMX*NVIRT*NVIRT 465 KEND2 = KINT2SB + NRHFT*NRHFT*NT1AMX 466 ELSE IF (LISTR.EQ.'R2 '.OR. LISTR.EQ.'ER1') THEN 467 KSCR1 = KEND2 468 KEND2 = KSCR1 + NT1AMX 469 END IF 470 471 472 LWRK2 = LWORK - KEND2 473 IF (LWRK2 .LT. 0) 474 & CALL QUIT('Out of memory in CCSDT_ADEN_NODDY (3)') 475 476C --------------------------------------- 477C Read precalculated integrals from file: 478C XINT1S0 = (CK|BD) 479C XINT2S0 = (CK|LJ) 480C --------------------------------------- 481 CALL CCSDT_READ_NODDY(.FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 482 & .FALSE.,DUMMY,DUMMY, 483 & .TRUE.,WORK(KINT1S0),WORK(KINT2S0), 484 & .FALSE.,DUMMY,DUMMY, 485 & .FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 486 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 487 488 489 CALL DZERO(WORK(KWYINT),NT1AMX*NT1AMX*NT1AMX) 490 491 IF (LISTR.EQ.'RE '.OR.LISTR.EQ.'RC ') THEN 492 KDUM = KEND2 493 CALL CCSDT_T31_NODDY(WORK(KWYINT),LISTR,IDLSTR,FREQR,.FALSE., 494 & .FALSE.,WORK(KINT1S0),WORK(KINT2S0), 495 & .FALSE.,DUMMY,DUMMY, 496 & .FALSE.,DUMMY,DUMMY, 497 & WORK(KINT1SB),WORK(KINT2SB), 498 & WORK(KLAMPB),WORK(KLAMHB),WORK(KFOCKB), 499 & XLAMDP0,XLAMDH0,FOCK0, 500 & WORK(KDUM),WORK(KFOCKD), 501 & WORK(KEND2),LWRK2) 502 ELSE IF (LISTR.EQ.'R2 '.OR. LISTR.EQ.'ER1') THEN 503 CALL CCSDT_T32_NODDY(WORK(KWYINT),LISTR,IDLSTR,FREQR, 504 & WORK(KINT1S0),WORK(KINT2S0), 505 & XLAMDP0,XLAMDH0,FOCK0, 506 & WORK(KFOCKD),DUMMY,DUMMY, 507 & WORK(KSCR1),WORK(KEND2),LWRK2) 508 ELSE 509 CALL QUIT('Unknown or illegal LISTR in CCSDT_ADEN_NODDY.') 510 END IF 511 512 IF (HAVETB3T3) THEN 513 CALL DCOPY(NT1AMX*NT1AMX*NT1AMX,WORK(KWYINT),1,WORK(KTAMP3),1) 514 END IF 515 516 IF (SPLIT_WT) THEN 517 ! Theta-Y is set to zero for this case 518 CALL DZERO(WORK(KTHETAY),NT1AMX*NT1AMX*NT1AMX) 519 END IF 520 521 ! account for the fact that in the real code a combination 522 ! of three W intermediates give the complete Tbar3 523 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,ONE/THREE,WORK(KWYINT),1) 524 525 LUTEMP = -1 526 CALL GPOPEN(LUTEMP,FILNODT30,'UNKNOWN',' ','UNFORMATTED', 527 & IDUMMY,.FALSE.) 528 READ(LUTEMP) (WORK(KT3AMP0+I-1), I=1,NT1AMX*NT1AMX*NT1AMX) 529 CALL GPCLOSE(LUTEMP,'KEEP') 530 531 IF (LOCDBG .AND. HAVETB3T3) THEN 532 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> T^Y:' 533 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> LISTR,IDLSTR:',LISTR,IDLSTR 534 WRITE(LUPRI,*) 'CCSDT_ADEN_NODDY> norm^2:', 535 & DDOT(NT1AMX*NT1AMX*NT1AMX,WORK(KTAMP3),1,WORK(KTAMP3),1) 536 CALL PRINT_PT3_NODDY(WORK(KTAMP3)) 537 END IF 538 539 END IF 540 541 ELSE ! IF (CUBIC) THEN 542 543 ! Get for both perturbations the intermediates: 544 ! w^y and theta^y with same meaning as above 545 ! YPROP contains matrix elements of y 546 ! In addition some second-order intermediates are set up: 547 ! XYMAT contains matrix elements of [X,T1^Y] + [Y,T1^X] 548 ! WXYVINT contains theta(e-bar m-bar, f n-bar, d l-bar) 549 ! WXYOINT contains theta(e m-bar-bar, f n-bar-bar, d l-bar-bar) 550 ! (returns also the zeroth-order amplitudes in T3AMP0, 551 ! and for double checking the complete T3^XY on TAMP3) 552 CALL CCSDT_T32INT_NODDY(LISTR,IDLSTR,ISYMR,FREQR, 553 * LISTX,IDLSTX,ISYMX, 554 * WORK(KWXINT),WORK(KTHETAX),HAVEX,WORK(KXPROP), 555 * LISTY,IDLSTY,ISYMY, 556 * WORK(KWYINT),WORK(KTHETAY),HAVEY,WORK(KYPROP), 557 * WORK(KT3AMP0),HAVEXY,WORK(KXYMAT), 558 * WORK(KTAMP1),WORK(KTAMP2), 559 * WORK(KTAMP3),WORK(KWXYOINT),WORK(KWXYVINT), 560 * XLAMDP0,XLAMDH0,FOCK0,WORK(KFOCKD), 561 * LUDELD,FNDELD,LUDKBC,FNDKBC,LUCKJD,FNCKJD, 562 * LOCDBG,HAVETB3T3, 563 * WORK(KEND1),LWRK1) 564 565 ! turn sign of t3 for compatibility with routines below 566 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-ONE,WORK(KT3AMP0),1) 567 END IF 568 569*---------------------------------------------------------------------* 570* if we have Tbar3 and T3 in core evaluate reference density 571* using the simple formulas: 572* D_ab = +1/2 sum_dlemn tbar(dl,em,an) ty(dl,em,bn) 573* D_ij = -1/2 sum_dlemf ty(dl,em,fi) tbar(dl,em,fj) 574* D_ia = -sum_dln sum_efm tbar(dl,em,fn) ty(em,fi) t0(dl,an) 575* -sum_dln sum_efm tbar(dl,em,fn) t0(em,fi) ty(dl,an) 576* +sum_dlem tbar(dl,em) [ty(dl,em,ai)-ty(dl,ei,am)] 577*---------------------------------------------------------------------* 578 IF (HAVETB3T3) THEN 579 CALL DZERO(WORK(KDABR),NVIRT*NVIRT) 580 CALL DZERO(WORK(KDIJR),NRHFT*NRHFT) 581 CALL DZERO(WORK(KDIAR),NVIRT*NRHFT) 582 CALL DZERO(WORK(KDAB0R),NVIRT*NVIRT) 583 CALL DZERO(WORK(KDIJ0R),NRHFT*NRHFT) 584 585 DO N = 1, NRHFT 586 IOFF = NT1AMX*NT1AMX*NVIRT*(N-1) 587 CALL DGEMM('T','N',NVIRT,NVIRT, NT1AMX*NT1AMX, 588 & +HALF,WORK(KTBAR3+IOFF),NT1AMX*NT1AMX, 589 & WORK(KTAMP3+IOFF),NT1AMX*NT1AMX, 590 & ONE,WORK(KDABR),NVIRT) 591 592 CALL DGEMM('T','N',NVIRT,NVIRT, NT1AMX*NT1AMX, 593 & +HALF,WORK(KTBAR3+IOFF),NT1AMX*NT1AMX, 594 & WORK(KT3AMP0+IOFF),NT1AMX*NT1AMX, 595 & ONE,WORK(KDAB0R),NVIRT) 596 END DO 597 598 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NT1AMX*NVIRT, 599 & -HALF,WORK(KTAMP3),NT1AMX*NT1AMX*NVIRT, 600 & WORK(KTBAR3),NT1AMX*NT1AMX*NVIRT, 601 & ZERO,WORK(KDIJR),NRHFT) 602 603 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NT1AMX*NVIRT, 604 & -HALF,WORK(KT3AMP0),NT1AMX*NT1AMX*NVIRT, 605 & WORK(KTBAR3),NT1AMX*NT1AMX*NVIRT, 606 & ZERO,WORK(KDIJ0R),NRHFT) 607 608 DO I = 1, NRHFT 609 DO A = 1, NVIRT 610 DO N = 1, NRHFT 611 DO F = 1, NVIRT 612 613 NAN = NVIRT*(N-1) + A 614 NAI = NVIRT*(I-1) + A 615 NFN = NVIRT*(N-1) + F 616 NFI = NVIRT*(I-1) + F 617 618 DO NDL = 1, NT1AMX 619 KDLFIAN = ((NAN-1)*NT1AMX + NFI-1) * NT1AMX + NDL 620 KDLFNAI = ((NAI-1)*NT1AMX + NFN-1) * NT1AMX + NDL 621 KDLFN = (NFN-1)*NT1AMX + NDL 622 623 WORK(KDIAR-1+NAI) = WORK(KDIAR-1+NAI) + 624 & ( WORK(KTAMP3-1+KDLFNAI) - WORK(KTAMP3-1+KDLFIAN) ) * 625 & WORK(KL2AMP-1+KDLFN) 626 END DO 627 628 DO NEM = 1, NT1AMX 629 DO NDL = 1, NT1AMX 630 KDLEMFN = ((NFN-1)*NT1AMX + NEM-1) * NT1AMX + NDL 631 KFNEMDL = ((NDL-1)*NT1AMX + NEM-1) * NT1AMX + NFN 632 KEMDLFN = ((NFN-1)*NT1AMX + NDL-1) * NT1AMX + NEM 633 KDLAN = (NAN-1)*NT1AMX + NDL 634 KEMFI = (NFI-1)*NT1AMX + NEM 635 636 WORK(KDIAR-1+NAI) = WORK(KDIAR-1+NAI) - 637 & WORK(KTBAR3-1+KDLEMFN) * 638 & ( WORK(KT0AMP-1+KEMFI) * WORK(KTAMP2-1+KDLAN) + 639 & WORK(KTAMP2-1+KEMFI) * WORK(KT0AMP-1+KDLAN) ) 640 641 END DO 642 END DO 643 644 END DO 645 END DO 646 END DO 647 END DO 648 649 IF (DO_XMMAT) THEN 650 651 DO NFN = 1, NT1AMX 652 DO M = 1, NRHFT 653 DO I = 1, NRHFT 654 655 KIMFN = ((NFN-1)*NRHFT + M-1) * NRHFT + I 656 WORK(KXMMAT-1+KIMFN) = ZERO 657 658 DO E = 1, NVIRT 659 NEI = NVIRT*(I-1) + E 660 NEM = NVIRT*(M-1) + E 661 662 DO NDL = 1, NT1AMX 663 KDLEMFN = ((NFN-1)*NT1AMX + NEM-1) * NT1AMX + NDL 664 KDLEI = (NEI-1) * NT1AMX + NDL 665 666 WORK(KXMMAT-1+KIMFN) = WORK(KXMMAT-1+KIMFN) - 667 & WORK(KTBAR3-1+KDLEMFN) * WORK(KTAMP2-1+KDLEI) 668 END DO 669 END DO 670 671 END DO 672 END DO 673 END DO 674 675 END IF 676 677 IF (LOCDBG) THEN 678 WRITE(LUPRI,*) 'The virtual-virtual block of the reference:' 679 CALL OUTPUT(WORK(KDABR),1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,LUPRI) 680 WRITE(LUPRI,*)'The occupied-occupied block of the reference:' 681 CALL OUTPUT(WORK(KDIJR),1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,LUPRI) 682 WRITE(LUPRI,*)'The occupied-virtual block of the reference:' 683 CALL OUTPUT(WORK(KDIAR),1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,LUPRI) 684 END IF 685 END IF 686 687*---------------------------------------------------------------------* 688* compute the densities: 689* first do some contributions in a loop over one virtual, then 690* do the remaining contributions in a loop over two occupieds 691*---------------------------------------------------------------------* 692 CALL DZERO(WORK(KDAB0),NVIRT*NVIRT) 693 CALL DZERO(WORK(KDIJ0),NRHFT*NRHFT) 694 CALL DZERO(DAB,NVIRT*NVIRT) 695 CALL DZERO(DIJ,NRHFT*NRHFT) 696 IF (DO_DIA) CALL DZERO(DIA,NVIRT*NRHFT) 697 698 KWTILDE = KEND1 699 KTHEOCC = KWTILDE + NT1AMX*NT1AMX 700 KEND2 = KTHEOCC + NT1AMX*NT1AMX 701 702 IF (SPLIT_WT) THEN 703 KWBREVE = KEND2 704 KTVIRT1 = KWBREVE + NT1AMX*NT1AMX 705 KTVIRT2 = KTVIRT1 + NT1AMX*NT1AMX 706 KEND2 = KTVIRT2 + NT1AMX*NT1AMX 707 END IF 708 709 IF (CUBIC) THEN 710 KTHEOV = KEND2 711 KTVIRT3 = KTHEOV + NT1AMX*NT1AMX 712 KEND2 = KTVIRT3 + NT1AMX*NT1AMX 713 END IF 714 715 KT0FMI = KEND2 716 KTYFMI = KT0FMI + NT1AMX*NRHFT 717 KT0MI = KTYFMI + NT1AMX*NRHFT 718 KTYMI = KT0MI + NRHFT*NRHFT 719 KWDE = KTYMI + NRHFT*NRHFT 720 KWED = KWDE + NT1AMX*NRHFT 721 KSCR = KWED + NT1AMX*NRHFT 722 KEND2 = KSCR + NT1AMX*NRHFT 723 724 LWRK2 = LWORK - KEND2 725 IF (LWRK2.LT.0) 726 & CALL QUIT('Out of memory in CCSDT_ADEN_NODDY (3b)') 727 728 IF (.NOT.CUBIC) KWXYOINT = KWYINT 729 730 CALL CCSDT_ADENVIR_NODDY(DAB,DIJ,DO_DIA,DIA, 731 & WORK(KDAB0),WORK(KDIJ0), 732 & SPLIT_WT,WORK(KWBAR),WORK(KWXYOINT), 733 & WORK(KTHETAY),WORK(KT3AMP0), 734 & CUBIC,WORK(KTHBAR),WORK(KWXYVINT), 735 & WORK(KTHETAX),FREQR, 736 & WORK(KXPROP),WORK(KYPROP), 737 & WORK(KTHEOV),WORK(KFOCKD), 738 & WORK(KWTILDE),WORK(KWBREVE), 739 & WORK(KTHEOCC),WORK(KTVIRT1), 740 & WORK(KTVIRT2),WORK(KTVIRT3), 741 & WORK(KT0AMP),WORK(KTAMP2),WORK(KL2AMP), 742 & WORK(KT0FMI),WORK(KTYFMI), 743 & WORK(KT0MI), WORK(KTYMI), 744 & WORK(KWDE), WORK(KWED), WORK(KSCR)) 745 746 IF (LOCDBG) THEN 747 write(lupri,*)'DAB in noddy after VIR' 748 call output(DAB,1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,lupri) 749 write(lupri,*)'DIJ in noddy after VIR' 750 call output(DIJ,1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,lupri) 751 if (do_dia) then 752 write(lupri,*)'DIA in noddy after VIR' 753 call output(DIA,1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,lupri) 754 end if 755 END IF 756 757 IF (SPLIT_WT .OR. SPLIT_WBTB) THEN 758 759 KWTILDE = KEND1 760 KTHVIRT = KWTILDE + NVIRT*NVIRT*NT1AMX 761 KSWAP = KTHVIRT + NVIRT*NVIRT*NT1AMX 762 KEND2 = KSWAP + NVIRT 763 764 IF (CUBIC) THEN 765 KSCRVVVO = KEND2 766 KEND2 = KSCRVVVO + NVIRT*NVIRT*NT1AMX 767 END IF 768 769 LWRK2 = LWORK - KEND2 770 IF (LWRK2.LT.0) 771 & CALL QUIT('Out of memory in CCSDT_ADEN_NODDY (3c)') 772 773 CALL CCSDT_ADENOCC_NODDY(DAB,DIJ,DO_DIA,DIA,WORK(KWBAR), 774 & WORK(KTHETAY),WORK(KL2AMP), 775 & CUBIC,WORK(KTHBAR),FREQR,WORK(KTHETAX), 776 & WORK(KWYINT),WORK(KWXINT), 777 & WORK(KFOCKD),WORK(KXPROP),WORK(KYPROP), 778 & WORK(KT3AMP0),WORK(KXYMAT), 779 & WORK(KWTILDE),WORK(KTHVIRT), 780 & WORK(KSWAP),WORK(KSCRVVVO)) 781 782 IF (LOCDBG) THEN 783 write(lupri,*)'DAB in noddy after OCC' 784 call output(DAB,1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,lupri) 785 write(lupri,*)'DIJ in noddy after OCC' 786 call output(DIJ,1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,lupri) 787 if (do_dia) then 788 write(lupri,*)'DIA in noddy after OCC' 789 call output(DIA,1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,lupri) 790 end if 791 END IF 792 793 END IF 794 795*---------------------------------------------------------------------* 796* add the one-index transformed zeroth-order density to DIA: 797*---------------------------------------------------------------------* 798 IF (DO_DIA) THEN 799 800 ! read T^Y singles amplitudes from file 801 IOPT = 1 802 Call CC_RDRSP(LISTR,IDLSTR,ISYMR,IOPT,MODEL,WORK(KTAMP1),DUMMY) 803 804 CALL DGEMM('T','N',NVIRT,NRHFT,NVIRT, 805 & -ONE,WORK(KDAB0),NVIRT,WORK(KTAMP1),NVIRT, 806 & ONE,DIA,NVIRT) 807 808 CALL DGEMM('N','T',NVIRT,NRHFT,NRHFT, 809 & ONE,WORK(KTAMP1),NVIRT,WORK(KDIJ0),NRHFT, 810 & ONE,DIA,NVIRT) 811 812 813 IF (LOCDBG) THEN 814 write(lupri,*)'DAB(total) in noddy after OCC' 815 call output(DAB,1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,lupri) 816 write(lupri,*)'DIJ(total) in noddy after OCC' 817 call output(DIJ,1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,lupri) 818 write(lupri,*)'DIA(total) in noddy after OCC' 819 call output(DIA,1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,lupri) 820 END IF 821 822 IF (HAVETB3T3) THEN 823 CALL DGEMM('T','N',NVIRT,NRHFT,NVIRT, 824 & -ONE,WORK(KDAB0R),NVIRT,WORK(KTAMP1),NVIRT, 825 & ONE,WORK(KDIAR),NVIRT) 826 827 CALL DGEMM('N','T',NVIRT,NRHFT,NRHFT, 828 & ONE,WORK(KTAMP1),NVIRT,WORK(KDIJ0R),NRHFT, 829 & ONE,WORK(KDIAR),NVIRT) 830 END IF 831 832 END IF 833*---------------------------------------------------------------------* 834* quick hack solution for the M intermediates (use noddy-noddy): 835*---------------------------------------------------------------------* 836 IF (DO_XMMAT) THEN 837 IF (.NOT. HAVETB3T3) CALL QUIT('DO_XMMAT needs HAVETB3T3!') 838 839 CALL DCOPY(NRHFT*NRHFT*NT1AMX,WORK(KXMMAT),1,XMMAT,1) 840 END IF 841*---------------------------------------------------------------------* 842* print debug output and return: 843*---------------------------------------------------------------------* 844 845 IF (PRINTDENS) THEN 846 WRITE(LUPRI,*) 'A densities for:' 847 WRITE(LUPRI,*) 'LISTL,IDLSTL:',LISTL,IDLSTL 848 WRITE(LUPRI,*) 'LISTR,IDLSTR:',LISTR,IDLSTR 849 850 WRITE(LUPRI,*) 'The virtual-virtual block DAB:' 851 CALL OUTPUT(DAB,1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,LUPRI) 852 IF (HAVETB3T3) THEN 853 WRITE(LUPRI,*) 'The virtual-virtual block of the reference:' 854 CALL OUTPUT(WORK(KDABR),1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,LUPRI) 855 CALL DAXPY(NVIRT*NVIRT,-ONE,DAB,1,WORK(KDABR),1) 856 WRITE(LUPRI,*)'DAB norm^2(difference):', 857 & DDOT(NVIRT*NVIRT,WORK(KDABR),1,WORK(KDABR),1) 858 END IF 859 860 WRITE(LUPRI,*) 'The occupied-occupied block DIJ:' 861 CALL OUTPUT(DIJ,1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,LUPRI) 862 IF (HAVETB3T3) THEN 863 WRITE(LUPRI,*)'The occupied-occupied block of the reference:' 864 CALL OUTPUT(WORK(KDIJR),1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,LUPRI) 865 CALL DAXPY(NRHFT*NRHFT,-ONE,DIJ,1,WORK(KDIJR),1) 866 WRITE(LUPRI,*)'DIJ norm^2(difference):', 867 & DDOT(NRHFT*NRHFT,WORK(KDIJR),1,WORK(KDIJR),1) 868 END IF 869 870 IF (DO_DIA) THEN 871 WRITE(LUPRI,*) 'The occupied-virtual block DIA:' 872 CALL OUTPUT(DIA,1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,LUPRI) 873 IF (HAVETB3T3) THEN 874 WRITE(LUPRI,*)'The occupied-virtual block of the reference:' 875 CALL OUTPUT(WORK(KDIAR),1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,LUPRI) 876 CALL DAXPY(NVIRT*NRHFT,-ONE,DIA,1,WORK(KDIAR),1) 877 WRITE(LUPRI,*)'DIA norm^2(difference):', 878 & DDOT(NVIRT*NRHFT,WORK(KDIAR),1,WORK(KDIAR),1) 879 END IF 880 881 WRITE(LUPRI,*) 'The virtual-virtual block DAB0:' 882 CALL OUTPUT(WORK(KDAB0),1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,LUPRI) 883 IF (HAVETB3T3) THEN 884 WRITE(LUPRI,*) 'The virtual-virtual block of the reference:' 885 CALL OUTPUT(WORK(KDAB0R),1,NVIRT,1,NVIRT,NVIRT,NVIRT,1,LUPRI) 886 CALL DAXPY(NVIRT*NVIRT,-ONE,WORK(KDAB0),1,WORK(KDAB0R),1) 887 WRITE(LUPRI,*)'DAB0 norm^2(difference):', 888 & DDOT(NVIRT*NVIRT,WORK(KDAB0R),1,WORK(KDAB0R),1) 889 END IF 890 891 WRITE(LUPRI,*) 'The occupied-occupied block DIJ0:' 892 CALL OUTPUT(WORK(KDIJ0),1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,LUPRI) 893 IF (HAVETB3T3) THEN 894 WRITE(LUPRI,*)'The occupied-occupied block of the reference:' 895 CALL OUTPUT(WORK(KDIJ0R),1,NRHFT,1,NRHFT,NRHFT,NRHFT,1,LUPRI) 896 CALL DAXPY(NRHFT*NRHFT,-ONE,WORK(KDIJ0),1,WORK(KDIJ0R),1) 897 WRITE(LUPRI,*)'DIJ0 norm^2(difference):', 898 & DDOT(NRHFT*NRHFT,WORK(KDIJ0R),1,WORK(KDIJ0R),1) 899 END IF 900 END IF 901 END IF 902 903 CALL QEXIT('CCTADNOD') 904 RETURN 905 END 906*---------------------------------------------------------------------* 907* END OF SUBROUTINE CCSDT_ADEN_NODDY * 908*---------------------------------------------------------------------* 909*=====================================================================* 910 SUBROUTINE CCSDT_ADENVIR_NODDY(DAB,DIJ,DO_DIA,DIA,DAB0,DIJ0, 911 & SPLIT_WT,WBAR,WYINT,THETAY,T0AMP, 912 & CUBIC,THBAR,WXYV,THETAX,FREQXY, 913 & XPROP,YPROP,THEOV,FOCKD, 914 & WTILDE,WBREVE,THEOCC, 915 & TVIRT1,TVIRT2,TVIRT3, 916 & TAMP02,TAMPR2,L2AMP,TXFMI,TYFMI, 917 & TXMI,TYMI,WDE,WED,SCR) 918*---------------------------------------------------------------------* 919* Purpose: set up loop over virtual, sort Wbar, W, theta as they * 920* would be sorted in the real code and compute the * 921* contributions to the densities * 922* * 923* Written by Christof Haettig, Jan 2003, Friedrichstal * 924*---------------------------------------------------------------------* 925 IMPLICIT NONE 926#include "priunit.h" 927#include "ccorb.h" 928#include "ccsdsym.h" 929 930 INTEGER ISYMWB, ISYMTH 931 PARAMETER (ISYMWB = 1, ISYMTH = 1) 932 933 LOGICAL CUBIC, DO_DIA, SPLIT_WT 934 935#if defined (SYS_CRAY) 936 REAL DIJ(NRHFT,NRHFT), DAB(NVIRT,NVIRT),DIA(NVIRT,NRHFT), 937 & DIJ0(*), DAB0(*) 938 REAL WBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 939 REAL THBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 940 REAL WYINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 941 REAL WXYV(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 942 REAL THETAY(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 943 REAL THETAX(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 944 REAL T0AMP(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 945 REAL WTILDE(NVIRT,NRHFT,NVIRT,NRHFT) 946 REAL WBREVE(NVIRT,NRHFT,NVIRT,NRHFT) 947 REAL THEOCC(NVIRT,NRHFT,NVIRT,NRHFT) 948 REAL THEOV(NVIRT,NRHFT,NVIRT,NRHFT) 949 REAL TVIRT1(NVIRT,NRHFT,NVIRT,NRHFT) 950 REAL TVIRT2(NVIRT,NRHFT,NVIRT,NRHFT) 951 REAL TVIRT3(NVIRT,NRHFT,NVIRT,NRHFT) 952 REAL TAMP02(NT1AMX,NVIRT,NRHFT) 953 REAL TAMPR2(NT1AMX,NVIRT,NRHFT) 954 REAL L2AMP(NVIRT,NRHFT,NVIRT,NRHFT) 955 REAL TXFMI(NVIRT,NRHFT,NRHFT) 956 REAL TYFMI(NVIRT,NRHFT,NRHFT) 957 REAL TXMI(NRHFT,NRHFT), TYMI(NRHFT,NRHFT) 958 REAL WDE(NVIRT,NRHFT,NRHFT) 959 REAL WED(NVIRT,NRHFT,NRHFT) 960 REAL SCR(NT1AMX*NRHFT), FOCKD(*) 961 REAL XPROP(NORBT,NORBT),YPROP(NORBT,NORBT) 962 REAL HALF, ONE, ZERO, DDOT, FREQXY 963#else 964 DOUBLE PRECISION DIJ(NRHFT,NRHFT),DAB(NVIRT,NVIRT), 965 & DIA(NVIRT,NRHFT),DIJ0(*),DAB0(*) 966 DOUBLE PRECISION WBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 967 DOUBLE PRECISION THBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 968 DOUBLE PRECISION WYINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 969 DOUBLE PRECISION WXYV(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 970 DOUBLE PRECISION THETAY(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 971 DOUBLE PRECISION THETAX(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 972 DOUBLE PRECISION T0AMP(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 973 DOUBLE PRECISION WTILDE(NVIRT,NRHFT,NVIRT,NRHFT) 974 DOUBLE PRECISION WBREVE(NVIRT,NRHFT,NVIRT,NRHFT) 975 DOUBLE PRECISION THEOCC(NVIRT,NRHFT,NVIRT,NRHFT) 976 DOUBLE PRECISION THEOV(NVIRT,NRHFT,NVIRT,NRHFT) 977 DOUBLE PRECISION TVIRT1(NVIRT,NRHFT,NVIRT,NRHFT) 978 DOUBLE PRECISION TVIRT2(NVIRT,NRHFT,NVIRT,NRHFT) 979 DOUBLE PRECISION TVIRT3(NVIRT,NRHFT,NVIRT,NRHFT) 980 DOUBLE PRECISION TAMP02(NT1AMX,NVIRT,NRHFT) 981 DOUBLE PRECISION TAMPR2(NT1AMX,NVIRT,NRHFT) 982 DOUBLE PRECISION L2AMP(NVIRT,NRHFT,NVIRT,NRHFT) 983 DOUBLE PRECISION TXFMI(NVIRT,NRHFT,NRHFT) 984 DOUBLE PRECISION TYFMI(NVIRT,NRHFT,NRHFT) 985 DOUBLE PRECISION TXMI(NRHFT,NRHFT), TYMI(NRHFT,NRHFT) 986 DOUBLE PRECISION WDE(NVIRT,NRHFT,NRHFT) 987 DOUBLE PRECISION WED(NVIRT,NRHFT,NRHFT) 988 DOUBLE PRECISION SCR(NT1AMX*NRHFT), FOCKD(*) 989 DOUBLE PRECISION XPROP(NORBT,NORBT),YPROP(NORBT,NORBT) 990 DOUBLE PRECISION HALF, ONE, ZERO, DDOT, FREQXY 991#endif 992 PARAMETER( HALF = 0.5D0, ONE = 1.0D0, ZERO = 0.0D0 ) 993 994 INTEGER NDM, NDL 995 996*---------------------------------------------------------------------* 997* start loop over D and E to compute contributions to D(IA) using 998* more or less the same algorithms as in CC3_XI_DEN_IA 999*---------------------------------------------------------------------* 1000 DO D = 1, NVIRT 1001 DO L = 1, NRHFT 1002 DO E = 1, NVIRT 1003 1004c ---------------------------------------------- 1005c get WDE(fnm) = W^DE(fnml) and WED = W^ED(fnml) 1006c ---------------------------------------------- 1007 DO M = 1, NRHFT 1008 CALL DCOPY(NT1AMX,WBAR(1,1,E,M,D,L),1,WDE(1,1,M),1) 1009 CALL DCOPY(NT1AMX,WBAR(1,1,D,M,E,L),1,WED(1,1,M),1) 1010 END DO 1011 1012c ---------------------------------------------- 1013c get TXFMI(fmi) = t^x(fm,Ei) and TYFMI 1014c and TXMI(mi) = t^x(Dm,Ei) and TYMI 1015c ---------------------------------------------- 1016 DO I = 1, NRHFT 1017 CALL DCOPY(NT1AMX,TAMP02(1,E,I),1,TXFMI(1,1,I),1) 1018 CALL DCOPY(NT1AMX,TAMPR2(1,E,I),1,TYFMI(1,1,I),1) 1019 DO M = 1, NRHFT 1020 NDM = NVIRT*(M-1) + D 1021 TXMI(M,I) = TAMP02(NDM,E,I) 1022 TYMI(M,I) = TAMPR2(NDM,E,I) 1023 END DO 1024 END DO 1025 1026c ------------------------------------------ 1027c D(IA) -= W^DE(fmnl) t^x(fm,Ei) t^y(an,DL) 1028c scr(n,i) 1029c ------------------------------------------ 1030 IF (DO_DIA) THEN 1031 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX, 1032 & ONE,WDE,NT1AMX,TXFMI,NT1AMX, 1033 & ZERO,SCR,NRHFT) 1034 CALL DGEMM('N','N',NVIRT,NRHFT,NRHFT, 1035 & -ONE,TAMPR2(1,D,L),NVIRT,SCR,NRHFT, 1036 & ONE,DIA,NVIRT) 1037 END IF 1038 1039c ------------------------------------------ 1040c D(IA) -= W^DE(fmnl) t^y(fm,Ei) t^x(an,DL) 1041c scr(n,i) 1042c ------------------------------------------ 1043 IF (DO_DIA) THEN 1044 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX, 1045 & ONE,WDE,NT1AMX,TYFMI,NT1AMX, 1046 & ZERO,SCR,NRHFT) 1047 CALL DGEMM('N','N',NVIRT,NRHFT,NRHFT, 1048 & -ONE,TAMP02(1,D,L),NVIRT,SCR,NRHFT, 1049 & ONE,DIA,NVIRT) 1050 END IF 1051 1052c -------------------------------------------- 1053c resort: WDE(fmi) --> WDE(fim) 1054c TXFMI(fmi) --> TXFMI(fim) 1055c TYFMI(fmi) --> TYFMI(fim) 1056c -------------------------------------------- 1057 DO M = 1, NRHFT 1058 DO I = M+1, NRHFT 1059 CALL DCOPY(NVIRT,WDE(1,M,I),1,SCR,1) 1060 CALL DCOPY(NVIRT,WDE(1,I,M),1,WDE(1,M,I),1) 1061 CALL DCOPY(NVIRT,SCR,1,WDE(1,I,M),1) 1062 1063 CALL DCOPY(NVIRT,TXFMI(1,M,I),1,SCR,1) 1064 CALL DCOPY(NVIRT,TXFMI(1,I,M),1,TXFMI(1,M,I),1) 1065 CALL DCOPY(NVIRT,SCR,1,TXFMI(1,I,M),1) 1066 1067 CALL DCOPY(NVIRT,TYFMI(1,M,I),1,SCR,1) 1068 CALL DCOPY(NVIRT,TYFMI(1,I,M),1,TYFMI(1,M,I),1) 1069 CALL DCOPY(NVIRT,SCR,1,TYFMI(1,I,M),1) 1070 END DO 1071 END DO 1072 1073c ------------------------------------------ 1074c D(IA) -= W^DE(fnml) t^x(fi,Em) t^y(an,Dl) 1075c scr(n,i) 1076c ------------------------------------------ 1077 IF (DO_DIA) THEN 1078 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX, 1079 & ONE,WDE,NT1AMX,TXFMI,NT1AMX, 1080 & ZERO,SCR,NRHFT) 1081 CALL DGEMM('N','N',NVIRT,NRHFT,NRHFT, 1082 & -ONE,TAMPR2(1,D,L),NVIRT,SCR,NRHFT, 1083 & ONE,DIA,NVIRT) 1084 END IF 1085 1086 1087c ------------------------------------------ 1088c D(IA) -= W^DE(fnml) t^y(fi,Em) t^x(an,Dl) 1089c scr(n,i) 1090c ------------------------------------------ 1091 IF (DO_DIA) THEN 1092 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX, 1093 & ONE,WDE,NT1AMX,TYFMI,NT1AMX, 1094 & ZERO,SCR,NRHFT) 1095 CALL DGEMM('N','N',NVIRT,NRHFT,NRHFT, 1096 & -ONE,TAMP02(1,D,L),NVIRT,SCR,NRHFT, 1097 & ONE,DIA,NVIRT) 1098 END IF 1099 1100 1101c ------------------------------------------ 1102c D(IA) -= W^ED(fnml) t^x(Dm,Ei) t^y(fn,al) 1103c scr(fni) 1104c ------------------------------------------ 1105 IF (DO_DIA) THEN 1106 CALL DGEMM('N','N',NT1AMX,NRHFT,NRHFT, 1107 & ONE,WED,NT1AMX,TXMI,NRHFT, 1108 & ZERO,SCR,NT1AMX) 1109 CALL DGEMM('T','N',NVIRT,NRHFT,NT1AMX, 1110 & -ONE,TAMPR2(1,1,L),NT1AMX,SCR,NT1AMX, 1111 & ONE,DIA,NVIRT) 1112 END IF 1113 1114c ------------------------------------------ 1115c D(IA) -= W^ED(fnml) t^y(Dm,Ei) t^x(fn,al) 1116c scr(fni) 1117c ------------------------------------------ 1118 IF (DO_DIA) THEN 1119 CALL DGEMM('N','N',NT1AMX,NRHFT,NRHFT, 1120 & ONE,WED,NT1AMX,TYMI,NRHFT, 1121 & ZERO,SCR,NT1AMX) 1122 CALL DGEMM('T','N',NVIRT,NRHFT,NT1AMX, 1123 & -ONE,TAMP02(1,1,L),NT1AMX,SCR,NT1AMX, 1124 & ONE,DIA,NVIRT) 1125 END IF 1126 1127 END DO 1128 END DO 1129 END DO 1130 1131*---------------------------------------------------------------------* 1132* start loop over D and L to compute contributions to D(AB), D(IJ) 1133*---------------------------------------------------------------------* 1134 DO D = 1, NVIRT 1135 DO L = 1, NRHFT 1136 1137c ---------------------------------------------------- 1138c build Wtilde^DL(em,an) = W^Da(emnl) + 1/2 W^De(anml) 1139c ---------------------------------------------------- 1140 CALL DCOPY(NT1AMX*NT1AMX,WBAR(1,1,1,1,D,L),1,WTILDE,1) 1141 CALL CC_T2MOD(WTILDE,ISYMWB,HALF) 1142 1143c ---------------------------------------------------- 1144c resort Wbreve^DL(em,fj) = W^Df(emlj) 1145c ---------------------------------------------------- 1146 IF (SPLIT_WT) THEN 1147 DO J = 1, NRHFT 1148 CALL DCOPY(NT1AMX*NVIRT,WBAR(1,1,1,L,D,J),1, 1149 & WBREVE(1,1,1,J),1) 1150 END DO 1151 END IF 1152 1153c ========================================================== 1154c build TVIRT2^DL intermediate: 1155c ========================================================== 1156 IF (SPLIT_WT) THEN 1157 IF (.NOT. CUBIC) THEN 1158 1159c ------------------------------------------------------------ 1160c TVIRT2^DL(em,fi) = theta(e-bar m,fl,di)+theta(f-bar l,em,di) 1161c ------------------------------------------------------------ 1162 DO I = 1, NRHFT 1163 DO F = 1, NVIRT 1164 DO M = 1, NRHFT 1165 DO E = 1, NVIRT 1166 TVIRT2(E,M,F,I)=THETAY(E,M,F,L,D,I)+THETAY(F,L,E,M,D,I) 1167 END DO 1168 END DO 1169 END DO 1170 END DO 1171 1172 ELSE 1173c -------------------------------------------------------- 1174c TVIRT2^DL(em,fi) = theta(f-bar-bar l, e m, d i) 1175c +theta(f-bar l, e-bar m, d i) 1176c +theta(f l, e-bar-bar m, d i) 1177c +theta(f-bar l-bar, e m-bar, d i-bar) 1178c (uses TVIRT1 array as scratch area!) 1179c -------------------------------------------------------- 1180 CALL DZERO(TVIRT2,NT1AMX*NT1AMX) 1181 1182 ! build TVIRT1^DL(em,fi) = thetay(em,fl,di)+thetay(fl,em,di) 1183 DO I = 1, NRHFT 1184 DO F = 1, NVIRT 1185 DO M = 1, NRHFT 1186 DO E = 1, NVIRT 1187 TVIRT1(E,M,F,I)=THETAY(E,M,F,L,D,I)+THETAY(F,L,E,M,D,I) 1188 END DO 1189 END DO 1190 END DO 1191 END DO 1192 1193 ! one-index transform. e and f with XPROP 1194 CALL CCSDT_TX1_NODDY(TVIRT2,TVIRT1,XPROP,.FALSE.,.TRUE., 1195 & NRHFT,NVIRT,NORBT) 1196 1197 ! build TVIRT1^DL(em,fi) = thetax(em,fl,di)+thetax(fl,em,di) 1198 DO I = 1, NRHFT 1199 DO F = 1, NVIRT 1200 DO M = 1, NRHFT 1201 DO E = 1, NVIRT 1202 TVIRT1(E,M,F,I)=THETAX(E,M,F,L,D,I)+THETAX(F,L,E,M,D,I) 1203 END DO 1204 END DO 1205 END DO 1206 END DO 1207 1208 ! one-index transform. e and f with XPROP 1209 CALL CCSDT_TX1_NODDY(TVIRT2,TVIRT1,YPROP,.FALSE.,.TRUE., 1210 & NRHFT,NVIRT,NORBT) 1211 1212 1213 ! divide by orbital energies and remove forbid. elements: 1214 CALL CCSDT_DIVFCL_NODDY(TVIRT2,D,L,FREQXY,FOCKD,SCR, 1215 & NT1AMX,NVIRT,NRHFT,NORBT) 1216 CALL DSCAL(NT1AMX*NT1AMX,-1.0D0,TVIRT2,1) 1217 1218 DO I = 1, NRHFT 1219 DO F = 1, NVIRT 1220 DO M = 1, NRHFT 1221 DO E = 1, NVIRT 1222 TVIRT2(E,M,F,I) = TVIRT2(E,M,F,I) + WXYV(F,L,E,M,D,I) 1223 END DO 1224 END DO 1225 END DO 1226 END DO 1227 1228 END IF 1229 END IF 1230 1231c ========================================================== 1232c build TVIRT1^DL / TVIRT3^DL intermediates: 1233c ========================================================== 1234 IF (SPLIT_WT) THEN 1235 IF (.NOT. CUBIC) THEN 1236c ---------------------------------------------------- 1237c TVIRT1^DL(fi,em) = theta(f-bar i,em,dl) 1238c ---------------------------------------------------- 1239 CALL DCOPY(NT1AMX*NT1AMX,THETAY(1,1,1,1,D,L),1,TVIRT1,1) 1240 1241 ELSE 1242c ----------------------------------------------------------- 1243c build TVIRT1^DL(em,fi) = theta(e m, f-bar-bar i, d l) 1244c +theta(e-bar m, f-bar i, d l) 1245c and TVIRT3^DL(em,fi) = theta(e-bar m, f-bar i, d l) 1246c ----------------------------------------------------------- 1247 CALL DZERO(TVIRT1,NT1AMX*NT1AMX) 1248 CALL DZERO(TVIRT3,NT1AMX*NT1AMX) 1249 1250 DO I = 1, NRHFT 1251 DO F = 1, NVIRT 1252 DO M = 1, NRHFT 1253 DO E = 1, NVIRT 1254 DO C = 1, NVIRT 1255 TVIRT3(E,M,F,I) = TVIRT3(E,M,F,I) 1256 & - THETAY(E,M,C,I,D,L) * XPROP(NRHFT+F,NRHFT+C) 1257 & - THETAY(F,I,C,M,D,L) * XPROP(NRHFT+E,NRHFT+C) 1258 & - THETAX(E,M,C,I,D,L) * YPROP(NRHFT+F,NRHFT+C) 1259 & - THETAX(F,I,C,M,D,L) * YPROP(NRHFT+E,NRHFT+C) 1260 1261 TVIRT1(E,M,F,I) = TVIRT1(E,M,F,I) 1262 & - THETAY(C,I,E,M,D,L) * XPROP(NRHFT+F,NRHFT+C) 1263 & - THETAX(C,I,E,M,D,L) * YPROP(NRHFT+F,NRHFT+C) 1264 END DO 1265 END DO 1266 END DO 1267 END DO 1268 END DO 1269 1270 CALL DAXPY(NT1AMX*NT1AMX,ONE,TVIRT3,1,TVIRT1,1) 1271 1272 ! divide by orbital energies and remove forbid. elements: 1273 CALL CCSDT_DIVFCL_NODDY(TVIRT1,D,L,FREQXY,FOCKD,SCR, 1274 & NT1AMX,NVIRT,NRHFT,NORBT) 1275 CALL DSCAL(NT1AMX*NT1AMX,-1.0D0,TVIRT1,1) 1276 1277 ! divide by orbital energies and remove forbid. elements: 1278 CALL CCSDT_DIVFCL_NODDY(TVIRT3,D,L,FREQXY,FOCKD,SCR, 1279 & NT1AMX,NVIRT,NRHFT,NORBT) 1280 CALL DSCAL(NT1AMX*NT1AMX,-1.0D0,TVIRT3,1) 1281 1282 END IF 1283 END IF 1284 1285c ============================================================ 1286c build THEOCC^DL(em,fn) = w(em,fn,dl)+w(dl,em,fn)+w(fn,dl,em) 1287c ============================================================ 1288 DO E = 1, NVIRT 1289 DO F = 1, NVIRT 1290 1291 DO M = 1, NRHFT 1292 DO N = 1, NRHFT 1293 1294 THEOCC(E,M,F,N) = WYINT(E,M,F,N,D,L)+ 1295 & WYINT(D,L,E,M,F,N)+WYINT(F,N,D,L,E,M) 1296 1297 END DO 1298 END DO 1299 1300 IF (CUBIC) THEN 1301 ! include contributions from 1302 ! WXYV(fn,em,dl) = theta(D l-bar, e m-bar, f-bar n-bar) 1303 ! and from 1304 ! WXYV(em,fn,dl) = theta(D l-bar, f n-bar, e-bar m-bar) 1305 DO M = 1, NRHFT 1306 DO N = 1, NRHFT 1307 THEOCC(E,M,F,N) = THEOCC(E,M,F,N)+ 1308 & WXYV(F,N,E,M,D,L) + WXYV(E,M,F,N,D,L) 1309 END DO 1310 END DO 1311 END IF 1312 1313 END DO 1314 END DO 1315 1316 IF (CUBIC) THEN 1317c ----------------------------------------------------------- 1318c fetch THEOV^DL(em,fn) = theta(D l-bar,f n-bar,e-bar m-bar) 1319c ----------------------------------------------------------- 1320 DO E = 1, NVIRT 1321 DO F = 1, NVIRT 1322 DO M = 1, NRHFT 1323 DO N = 1, NRHFT 1324 1325 THEOV(E,M,F,N) = WXYV(E,M,F,N,D,L) 1326 1327 END DO 1328 END DO 1329 END DO 1330 END DO 1331 END IF 1332 1333*---------------------------------------------------------------------* 1334* now we have the following intermediates in O^2 N^2 arrays: 1335* 1336* WTILDE^DL(em,an) = W^Da(emnl) + 1/2 W^De(anml) 1337* WBREVE^DL(em,fj) = W^Df(emlj) 1338* 1339* for quadratic response: 1340* TVIRT2^DL(em,fi) = theta(e-bar m,fl,di)+theta(f-bar l,em,di) 1341* TVIRT1^DL(fi,em) = theta(f-bar i,em,dl) 1342* THEOCC^DL(em,fn) = w(em,fn,dl)+w(dl,em,fn)+w(fn,dl,em) 1343* 1344* for cubic response: 1345* TVIRT2^DL(em,fi) = theta(e-bar-bar m, f l,di) 1346* + theta(e-bar m, f-bar l,di) 1347* + theta(e m, f-bar-bar l,di) 1348* + theta(e m-bar, f-bar l-bar, d i-bar) 1349* TVIRT1^DL(em,fi) = theta(e m, f-bar-bar i,dl) 1350* + theta(e-bar m, f-bar i,dl) 1351* TVIRT3^DL(em,fi) = theta(e-bar m, f-bar i,dl) 1352* THEOCC^DL(em,fn) = wxy(em,fn,dl)+wxy(dl,em,fn)+wxy(fn,dl,em) 1353* + theta(f-bar n-bar, e m-bar, D l-bar) 1354* + theta(e-bar m-bar, f n-bar, D l-bar) 1355* THEOV^DL(em,fn) = theta(e-bar m-bar, f n-bar, D l-bar) 1356* 1357* with these we can compute the contributions to D(ij) and D(ab) 1358* for quadratic response with 3 calls to DGEMM, and 1359* for cubic response with 7 calls to DGEMM, 1360*---------------------------------------------------------------------* 1361 1362c -------------------------------------------- 1363c D(ab) += WTILDE^DL(em,an) * THEOCC^DL(em,bn) 1364c -------------------------------------------- 1365 DO N = 1, NRHFT 1366 CALL DGEMM('T','N',NVIRT,NVIRT,NT1AMX, 1367 & ONE,WTILDE(1,1,1,N),NT1AMX, 1368 & THEOCC(1,1,1,N),NT1AMX, 1369 & ONE,DAB,NVIRT) 1370 END DO 1371 1372c -------------------------------------------- 1373c D0(ab) += WTILDE^DL(em,an) * T0^DL(em,bn) 1374c -------------------------------------------- 1375 DO N = 1, NRHFT 1376 CALL DGEMM('T','N',NVIRT,NVIRT,NT1AMX, 1377 & ONE,WTILDE(1,1,1,N),NT1AMX, 1378 & T0AMP(1,1,1,N,D,L),NT1AMX, 1379 & ONE,DAB0,NVIRT) 1380 END DO 1381 1382c --------------------------------------------------------- 1383c compute the contribution from <L_2|[E_ia,T_3^xy]|HF> 1384c to D(ia) which can be done in the loop over virtuals: 1385c 1386c D(ia) += [THEOCC^DL(em,ai)-THEOCC^DL(ei,am)] * L(em,DL) 1387c 1388c for cubic response include in addition: 1389c D(DL) += THEOV^DL(ai,em) * L(ai,em) 1390c D(Di) -= THEOV^DL(am,ei) * L(am,eL) 1391c ----------------------------------------------------------- 1392 IF (DO_DIA) THEN 1393 DO I = 1, NRHFT 1394 DO A = 1, NVIRT 1395 DO M = 1, NRHFT 1396 DIA(A,I) = DIA(A,I) + 1397 & DDOT(NVIRT,THEOCC(1,M,A,I),1,L2AMP(1,M,D,L),1) - 1398 & DDOT(NVIRT,THEOCC(1,I,A,M),1,L2AMP(1,M,D,L),1) 1399 END DO 1400 END DO 1401 END DO 1402 1403 IF (CUBIC) THEN 1404 DIA(D,L) = DIA(D,L) + DDOT(NT1AMX*NT1AMX,THEOV,1,L2AMP,1) 1405 1406 CALL DGEMV('T',NT1AMX*NVIRT,NRHFT, 1407 * -ONE,THEOV,NT1AMX*NVIRT,L2AMP(1,1,1,L),1, 1408 * ONE,DIA(D,1),NVIRT) 1409 END IF 1410 END IF 1411 1412c ---------------------------------------------------------- 1413c add TVIRT1 to THEOCC: THEOCC^DL(em,fi) += TVIRT1^DL(fi,em) 1414c (skipped for cubic response) 1415c ---------------------------------------------------------- 1416 IF (SPLIT_WT .AND. (.NOT. CUBIC)) THEN 1417 DO M = 1, NRHFT 1418 DO E = 1, NVIRT 1419 DO I = 1, NRHFT 1420 DO F = 1, NVIRT 1421 THEOCC(E,M,F,I) = THEOCC(E,M,F,I)+TVIRT1(F,I,E,M) 1422 END DO 1423 END DO 1424 END DO 1425 END DO 1426 END IF 1427 1428c ----------------------------------------------------------- 1429c D(ij)+=WTILDE^DL(em,fj)*[THEOCC^DL(em,fi)+TVIRT1^DL(fi,em)] 1430c ----------------------------------------------------------- 1431 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NVIRT, 1432 & -ONE,THEOCC,NT1AMX*NVIRT, 1433 & WTILDE,NT1AMX*NVIRT, 1434 & ONE,DIJ,NRHFT) 1435 1436c ----------------------------------------------------------- 1437c D(ij) -= WBREVE^DL(em,fj) * TVIRT2^DL(em,fi) 1438c ----------------------------------------------------------- 1439 IF (SPLIT_WT) THEN 1440 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NVIRT, 1441 & -ONE,TVIRT2,NT1AMX*NVIRT, 1442 & WBREVE,NT1AMX*NVIRT, 1443 & ONE,DIJ,NRHFT) 1444 END IF 1445 1446c ----------------------------------------------------------- 1447c D0(ij) -= WTILDE^DL(em,fj) * T0AMP^DL(em,fi) 1448c ----------------------------------------------------------- 1449 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NVIRT, 1450 & -ONE,T0AMP(1,1,1,1,D,L),NT1AMX*NVIRT, 1451 & WTILDE,NT1AMX*NVIRT, 1452 & ONE,DIJ0,NRHFT) 1453 1454*---------------------------------------------------------------------* 1455* add contributions without counterpart in quadratic response: 1456*---------------------------------------------------------------------* 1457 IF (CUBIC) THEN 1458 1459c ----------------------------------------------------------- 1460c D(ij) -= W^DL(em,fj) * TVIRT1^DL(em,fi) 1461c ----------------------------------------------------------- 1462 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NVIRT, 1463 & -ONE,TVIRT1,NT1AMX*NVIRT, 1464 & WBAR(1,1,1,1,D,L),NT1AMX*NVIRT, 1465 & ONE,DIJ,NRHFT) 1466 1467c --------------------------------------------------------- 1468c resort WBREVE^DL(em,fj) = Wbar^De(fj,ml) = WBAR(fj,em,DL) 1469c (i.e. transpose the pair indeces em, fj) 1470c --------------------------------------------------------- 1471 DO F = 1, NVIRT 1472 DO J = 1, NRHFT 1473 CALL DCOPY(NT1AMX,WBAR(F,J,1,1,D,L),NT1AMX, 1474 & WBREVE(1,1,F,J),1) 1475 END DO 1476 END DO 1477 1478c ----------------------------------------------------------- 1479c D(ij) -= 1/2 WBREVE^DL(em,fj) * [THEOV^DL(em,fi) + ...] 1480c ----------------------------------------------------------- 1481 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NVIRT, 1482 & -HALF,THEOV,NT1AMX*NVIRT, 1483 & WBREVE,NT1AMX*NVIRT, 1484 & ONE,DIJ,NRHFT) 1485 1486c ----------------------------------------------------------- 1487c get in WTILDE^DL(em,fj) = wbar^D(ef;l-bar m j) 1488c ----------------------------------------------------------- 1489 DO J = 1, NRHFT 1490 DO F = 1, NVIRT 1491 DO M = 1, NRHFT 1492 DO E = 1, NVIRT 1493 WTILDE(E,M,F,J) = WBAR(D,L,E,M,F,J) - THBAR(D,L,E,M,F,J) 1494 END DO 1495 END DO 1496 END DO 1497 END DO 1498 1499c ----------------------------------------------------------- 1500c D(ij) -= W^ef(DL,jm) * Theta^DL(e-bar m, f-bar i) 1501c ----------------------------------------------------------- 1502 CALL DGEMM('T','N',NRHFT,NRHFT,NT1AMX*NVIRT, 1503 & -ONE,TVIRT3,NT1AMX*NVIRT, 1504 & WTILDE,NT1AMX*NVIRT, 1505 & ONE,DIJ,NRHFT) 1506 1507c ----------------------------------------------------------- 1508c D(ab) += [ 1/2 W^De(an,ml) + wbar^D(ea; l-bar m n) ] 1509c * THEOV^DL(em,bn) 1510c ----------------------------------------------------------- 1511 CALL DAXPY(NT1AMX*NT1AMX,HALF,WBREVE,1,WTILDE,1) 1512 DO N = 1, NRHFT 1513 CALL DGEMM('T','N',NVIRT,NVIRT,NT1AMX, 1514 & ONE,WTILDE(1,1,1,N),NT1AMX, 1515 & THEOV(1,1,1,N),NT1AMX, 1516 & ONE,DAB,NVIRT) 1517 END DO 1518 1519 END IF 1520 END DO 1521 END DO 1522 1523 RETURN 1524 END 1525*---------------------------------------------------------------------* 1526* END OF SUBROUTINE CCSDT_ADENVIR_NODDY * 1527*---------------------------------------------------------------------* 1528*=====================================================================* 1529 SUBROUTINE CCSDT_TX1_NODDY(TX,TAMP,XMAT,LOCC,LVIRT, 1530 & NRHFT,NVIRT,NORBT) 1531*---------------------------------------------------------------------* 1532c TX(em,fn) += sum_k T(ek,fn) X_km - sum_c T(cm,fn) X_ec 1533c +sum_k T(em,fk) X_kn - sum_c T(em,cn) X_fc 1534c LOCC = .TRUE. : include transformation of occupied index 1535c LVIRT = .TRUE. : include transformation of virtual index 1536*---------------------------------------------------------------------* 1537 IMPLICIT NONE 1538 1539 LOGICAL LOCC, LVIRT 1540 INTEGER NVIRT,NRHFT,NORBT,E,F,M,N,C,K 1541 1542#if defined (SYS_CRAY) 1543 REAL TAMP(NVIRT,NRHFT,NVIRT,NRHFT) 1544 REAL TX(NVIRT,NRHFT,NVIRT,NRHFT) 1545 REAL XMAT(NORBT,NORBT) 1546#else 1547 DOUBLE PRECISION TAMP(NVIRT,NRHFT,NVIRT,NRHFT) 1548 DOUBLE PRECISION TX(NVIRT,NRHFT,NVIRT,NRHFT) 1549 DOUBLE PRECISION XMAT(NORBT,NORBT) 1550#endif 1551 1552 DO F = 1, NVIRT 1553 DO N = 1, NRHFT 1554 1555 IF (LVIRT) THEN 1556 DO E = 1, NVIRT 1557 DO M = 1, NRHFT 1558 DO C = 1, NVIRT 1559 TX(E,M,F,N) = TX(E,M,F,N) - 1560 & ( TAMP(C,M,F,N)*XMAT(NRHFT+E,NRHFT+C) + 1561 & TAMP(E,M,C,N)*XMAT(NRHFT+F,NRHFT+C) ) 1562 END DO 1563 END DO 1564 END DO 1565 END IF 1566 1567 IF (LOCC) THEN 1568 DO E = 1, NVIRT 1569 DO M = 1, NRHFT 1570 DO K = 1, NRHFT 1571 TX(E,M,F,N) = TX(E,M,F,N) + 1572 & ( TAMP(E,K,F,N)*XMAT(K,M) + TAMP(E,M,F,K)*XMAT(K,N) ) 1573 END DO 1574 END DO 1575 END DO 1576 END IF 1577 1578 END DO 1579 END DO 1580 1581 RETURN 1582 END 1583*---------------------------------------------------------------------* 1584*=====================================================================* 1585 SUBROUTINE CCSDT_DIVFCL_NODDY(TAMP,NC,NK,FREQ,FOCKD,SCR1, 1586 & NT1AMX,NVIRT,NRHFT,NORBT) 1587*---------------------------------------------------------------------* 1588c 1) T^CK(ai,bj) = T^CK(ai,bj) / (eps_aibjck - freq) 1589c 2) remove forbidden elements 1590*---------------------------------------------------------------------* 1591 IMPLICIT NONE 1592 1593 INTEGER NVIRT,NRHFT,NORBT,NT1AMX,NAI,NBJ,NCK,NK,NC 1594 INTEGER NA,NB,NAK,NBK,NI,NJ,NCI,NCJ 1595 1596#if defined (SYS_CRAY) 1597 REAL TAMP(NT1AMX,NT1AMX), FOCKD(NORBT), SCR1(NT1AMX) 1598 REAL FREQ 1599#else 1600 DOUBLE PRECISION TAMP(NT1AMX,NT1AMX), FOCKD(NORBT), SCR1(NT1AMX) 1601 DOUBLE PRECISION FREQ 1602#endif 1603 1604 DO NI = 1,NRHFT 1605 DO NA = 1,NVIRT 1606 NAI = NVIRT*(NI-1) + NA 1607 SCR1(NAI) = FOCKD(NRHFT+NA) - FOCKD(NI) 1608 END DO 1609 END DO 1610 1611 NCK = NVIRT*(NK-1) + NC 1612 1613 DO NBJ = 1, NT1AMX 1614 DO NAI = 1, NT1AMX 1615 TAMP(NAI,NBJ) = TAMP(NAI,NBJ) / 1616 & (SCR1(NAI) + SCR1(NBJ) + SCR1(NCK) - FREQ) 1617 END DO 1618 END DO 1619 1620 DO NB = 1, NVIRT 1621 NBK = NVIRT*(NK-1) + NB 1622 DO NA = 1, NVIRT 1623 NAK = NVIRT*(NK-1) + NA 1624 TAMP(NAK,NBK) = 0.0d0 1625 ENDDO 1626 ENDDO 1627 DO NJ = 1, NRHFT 1628 NCJ = NVIRT*(NJ-1) + NC 1629 DO NI = 1, NRHFT 1630 NCI = NVIRT*(NI-1) + NC 1631 TAMP(NCI,NCJ) = 0.0d0 1632 ENDDO 1633 ENDDO 1634 1635 RETURN 1636 END 1637*---------------------------------------------------------------------* 1638*=====================================================================* 1639 SUBROUTINE CCSDT_ADENOCC_NODDY(DAB,DIJ,DO_DIA,DIA, 1640 & WBAR,THETA,L2AMP, 1641 & CUBIC,THBAR,FREQR2,THETAX, 1642 & WYINT,WXINT,FOCKD,XPROP,YPROP, 1643 & T0AMP,XYMAT, 1644 & WTILDE,THVIRT,SWAP,SCR) 1645*---------------------------------------------------------------------* 1646* Purpose: set up loop over occupied, sort Wbar, W, theta as they * 1647* would be sorted in the real code and compute the * 1648* contributions to the densities * 1649* * 1650* Written by Christof Haettig, Jan 2003, Friedrichstal * 1651*---------------------------------------------------------------------* 1652 IMPLICIT NONE 1653#include "priunit.h" 1654#include "ccorb.h" 1655#include "ccsdsym.h" 1656 1657 LOGICAL CUBIC, DO_DIA 1658 1659#if defined (SYS_CRAY) 1660 REAL DIJ(*), DAB(NVIRT,NVIRT), DIA(NVIRT,NRHFT) 1661 REAL T0AMP(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1662 REAL WBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1663 REAL THBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1664 REAL THETA(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1665 REAL THETAX(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1666 REAL WYINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1667 REAL WXINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1668 REAL WTILDE(NVIRT,NVIRT,NVIRT,NRHFT) 1669 REAL THVIRT(NVIRT,NVIRT,NVIRT,NRHFT) 1670 REAL SCR(NVIRT,NVIRT,NVIRT,NRHFT) 1671 REAL L2AMP(NVIRT,NRHFT,NVIRT,NRHFT) 1672 REAL FREQR2, XPROP(NORBT,NORBT), YPROP(NORBT,NORBT), FOCKD(*) 1673 REAL SWAP(NVIRT), XYMAT(NORBT,NORBT) 1674 REAL HALF, ONE, DDOT 1675#else 1676 DOUBLE PRECISION DIJ(*), DAB(NVIRT,NVIRT), DIA(NVIRT,NRHFT) 1677 DOUBLE PRECISION T0AMP(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1678 DOUBLE PRECISION WBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1679 DOUBLE PRECISION THBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1680 DOUBLE PRECISION THETA(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1681 DOUBLE PRECISION THETAX(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1682 DOUBLE PRECISION WYINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1683 DOUBLE PRECISION WXINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 1684 DOUBLE PRECISION WTILDE(NVIRT,NVIRT,NVIRT,NRHFT) 1685 DOUBLE PRECISION THVIRT(NVIRT,NVIRT,NVIRT,NRHFT) 1686 DOUBLE PRECISION SCR(NVIRT,NVIRT,NVIRT,NRHFT) 1687 DOUBLE PRECISION L2AMP(NVIRT,NRHFT,NVIRT,NRHFT) 1688 DOUBLE PRECISION FREQR2, XPROP(NORBT,NORBT), YPROP(NORBT,NORBT) 1689 DOUBLE PRECISION SWAP(NVIRT), FOCKD(*), XYMAT(NORBT,NORBT) 1690 DOUBLE PRECISION HALF, ONE, DDOT 1691#endif 1692 PARAMETER( HALF = 0.5D0, ONE = 1.0D0 ) 1693 1694*---------------------------------------------------------------------* 1695* start loop over L and M 1696*---------------------------------------------------------------------* 1697 DO L = 1, NRHFT 1698 DO M = 1, NRHFT 1699 1700c ---------------------------------------------------- 1701c resort WTILDE^LM(de,an) = W^DE(anml) = W^LM(naed) 1702c ---------------------------------------------------- 1703 DO N = 1, NRHFT 1704 DO A = 1, NVIRT 1705 DO E = 1, NVIRT 1706 DO D = 1, NVIRT 1707 WTILDE(D,E,A,N) = WBAR(A,N,E,M,D,L) 1708 END DO 1709 END DO 1710 END DO 1711 END DO 1712 1713 IF (.NOT. CUBIC) THEN 1714c ------------------------------------------------ 1715c resort THVIRT^LM(de,an) = theta(d-bar l,em,an) + 1716c theta(dl,e-bar m,an) 1717c ------------------------------------------------ 1718 DO N = 1, NRHFT 1719 DO A = 1, NVIRT 1720 DO E = 1, NVIRT 1721 DO D = 1, NVIRT 1722 THVIRT(D,E,A,N)=THETA(D,L,E,M,A,N)+THETA(E,M,D,L,A,N) 1723 END DO 1724 END DO 1725 END DO 1726 END DO 1727 1728 ELSE 1729c ---------------------------------------------------- 1730c for cubic response transform THVIRT further and add 1731c add contributions to obtain in THVIRT^LM the 1732c theta intermediate with double bar on all virtuals 1733c ---------------------------------------------------- 1734 CALL DZERO(THVIRT,NVIRT*NVIRT*NT1AMX) 1735 1736 ! get SCR(de,an) = theta(d-bar^y l,em,an) + 1737 ! theta(dl,e-bar^y m,an) + theta(dl,em,a-bar^y n) 1738 DO N = 1, NRHFT 1739 DO A = 1, NVIRT 1740 DO E = 1, NVIRT 1741 DO D = 1, NVIRT 1742 SCR(D,E,A,N) = THETA(D,L,E,M,A,N) + 1743 & THETA(E,M,A,N,D,L) + THETA(A,N,D,L,E,M) 1744 END DO 1745 END DO 1746 END DO 1747 END DO 1748 1749 ! transform virtual indeces with X 1750 call CCSDT_TXVIR_NODDY(THVIRT,SCR,XPROP,0, 1751 & NVIRT,NRHFT,NORBT) 1752 1753 ! get SCR(de,an) = theta(d-bar^x l,em,an) + 1754 ! theta(dl,e-bar^x m,an) + theta(dl,em,a-bar^x n) 1755 DO N = 1, NRHFT 1756 DO A = 1, NVIRT 1757 DO E = 1, NVIRT 1758 DO D = 1, NVIRT 1759 SCR(D,E,A,N) = THETAX(D,L,E,M,A,N) + 1760 & THETAX(E,M,A,N,D,L) + THETAX(A,N,D,L,E,M) 1761 END DO 1762 END DO 1763 END DO 1764 END DO 1765 1766 ! transform virtual indeces with Y 1767 call CCSDT_TXVIR_NODDY(THVIRT,SCR,YPROP,0, 1768 & NVIRT,NRHFT,NORBT) 1769 1770 ! divide by orbital energies and remove forbidden entries 1771 CALL CCSDT_DIVVIR_NODDY(THVIRT,L,M,FREQR2, 1772 & FOCKD,FOCKD(NRHFT+1),NVIRT,NRHFT) 1773 1774 CALL DSCAL(NVIRT*NVIRT*NT1AMX,-1.0D0,THVIRT,1) 1775 END IF 1776 1777c -------------------------------------------- 1778c D(ij) += WTILDE^LM(de,fj) * THVIRT^LM(de,fi) 1779c -------------------------------------------- 1780 CALL DGEMM('T','N',NRHFT,NRHFT,NVIRT*NVIRT*NVIRT, 1781 & -HALF,THVIRT,NVIRT*NVIRT*NVIRT, 1782 & WTILDE,NVIRT*NVIRT*NVIRT, 1783 & ONE,DIJ,NRHFT) 1784 1785 IF (.NOT. CUBIC) THEN 1786c -------------------------------------------- 1787c add THVIRT^LM(de,an) += theta(dl,em,a-bar n) 1788c -------------------------------------------- 1789 DO N = 1, NRHFT 1790 DO A = 1, NVIRT 1791 DO E = 1, NVIRT 1792 DO D = 1, NVIRT 1793 THVIRT(D,E,A,N) = THVIRT(D,E,A,N) + THETA(A,N,E,M,D,L) 1794 END DO 1795 END DO 1796 END DO 1797 END DO 1798 END IF 1799 1800c -------------------------------------------- 1801c D(ab) += WTILDE^LM(de,an) * THVIRT^LM(de,bn) 1802c -------------------------------------------- 1803 DO N = 1, NRHFT 1804 CALL DGEMM('T','N',NVIRT,NVIRT,NVIRT*NVIRT, 1805 & HALF,WTILDE(1,1,1,N),NVIRT*NVIRT, 1806 & THVIRT(1,1,1,N),NVIRT*NVIRT, 1807 & ONE,DAB,NVIRT) 1808 END DO 1809 1810c ----------------------------------------------------- 1811c D(ia) += [THVIRT^LM(de,ai)-THVIRT^LM(da,ei)]*L(dL,eM) 1812c ----------------------------------------------------- 1813 IF (DO_DIA) THEN 1814 DO I = 1, NRHFT 1815 DO A = 1, NVIRT 1816 DO E = 1, NVIRT 1817 DIA(A,I) = DIA(A,I) + 1818 & DDOT(NVIRT,THVIRT(1,E,A,I),1,L2AMP(1,L,E,M),1) - 1819 & DDOT(NVIRT,THVIRT(1,A,E,I),1,L2AMP(1,L,E,M),1) 1820 END DO 1821 END DO 1822 END DO 1823 END IF 1824 1825c -------------------------------------------- 1826c resort: WTILDE^LM(de,an) --> WTILDE(da,en) 1827c THVIRT^LM(de,an) --> THVIRT(da,en) 1828c -------------------------------------------- 1829 DO N = 1, NRHFT 1830 DO A = 1, NVIRT 1831 DO E = A+1, NVIRT 1832 CALL DCOPY(NVIRT,WTILDE(1,E,A,N),1,SWAP,1) 1833 CALL DCOPY(NVIRT,WTILDE(1,A,E,N),1,WTILDE(1,E,A,N),1) 1834 CALL DCOPY(NVIRT,SWAP,1,WTILDE(1,A,E,N),1) 1835 1836 CALL DCOPY(NVIRT,THVIRT(1,E,A,N),1,SWAP,1) 1837 CALL DCOPY(NVIRT,THVIRT(1,A,E,N),1,THVIRT(1,E,A,N),1) 1838 CALL DCOPY(NVIRT,SWAP,1,THVIRT(1,A,E,N),1) 1839 END DO 1840 END DO 1841 END DO 1842 1843c -------------------------------------------- 1844c D(ab) += WTILDE^LM(da,en) * THVIRT^LM(db,en) 1845c -------------------------------------------- 1846 DO N = 1, NRHFT 1847 CALL DGEMM('T','N',NVIRT,NVIRT,NVIRT*NVIRT, 1848 & ONE,WTILDE(1,1,1,N),NVIRT*NVIRT, 1849 & THVIRT(1,1,1,N),NVIRT*NVIRT, 1850 & ONE,DAB,NVIRT) 1851 END DO 1852 1853 IF (CUBIC) THEN 1854c ------------------------------------------------------ 1855c resort WTILDE^LM(de,an) = theta-bar(d-bar L, e M, a n) 1856c ------------------------------------------------------ 1857 DO N = 1, NRHFT 1858 DO A = 1, NVIRT 1859 DO E = 1, NVIRT 1860 DO D = 1, NVIRT 1861 WTILDE(D,E,A,N) = THBAR(D,L,E,M,A,N) 1862 END DO 1863 END DO 1864 END DO 1865 END DO 1866 1867c ------------------------------------------------------ 1868c set up THVIRT^LM(de,an) = theta(e-bar M, a-bar n, d L) 1869c ------------------------------------------------------ 1870 CALL DZERO(THVIRT,NVIRT*NVIRT*NT1AMX) 1871 1872 DO N = 1, NRHFT 1873 DO A = 1, NVIRT 1874 DO E = 1, NVIRT 1875 DO D = 1, NVIRT 1876 DO C = 1, NVIRT 1877 THVIRT(D,E,A,N) = THVIRT(D,E,A,N) - 1878 & ( THETA(A,N,C,M,D,L) * XPROP(NRHFT+E,NRHFT+C) + 1879 & THETA(E,M,C,N,D,L) * XPROP(NRHFT+A,NRHFT+C) + 1880 & THETAX(A,N,C,M,D,L) * YPROP(NRHFT+E,NRHFT+C) + 1881 & THETAX(E,M,C,N,D,L) * YPROP(NRHFT+A,NRHFT+C) ) 1882 END DO 1883 END DO 1884 END DO 1885 END DO 1886 END DO 1887 1888 ! divide by orbital energies and remove forbidden entries 1889 CALL CCSDT_DIVVIR_NODDY(THVIRT,L,M,FREQR2, 1890 & FOCKD,FOCKD(NRHFT+1),NVIRT,NRHFT) 1891 1892 CALL DSCAL(NVIRT*NVIRT*NT1AMX,-1.0D0,THVIRT,1) 1893 1894 1895c ---------------------------------------------------- 1896c D(ij) += theta-bar(d-bar L,eM,aj) * THVIRT^LM(de,fi) 1897c ---------------------------------------------------- 1898 CALL DGEMM('T','N',NRHFT,NRHFT,NVIRT*NVIRT*NVIRT, 1899 & -ONE,THVIRT,NVIRT*NVIRT*NVIRT, 1900 & WTILDE,NVIRT*NVIRT*NVIRT, 1901 & ONE,DIJ,NRHFT) 1902 1903 1904c ------------------------------------------------------ 1905c set up THVIRT^LM(da,en) = theta(d L, a M, e-bar n-bar) 1906c ------------------------------------------------------ 1907 CALL DZERO(THVIRT,NVIRT*NVIRT*NT1AMX) 1908 1909 DO N = 1, NRHFT 1910 DO A = 1, NVIRT 1911 DO E = 1, NVIRT 1912 DO D = 1, NVIRT 1913 DO C = 1, NVIRT 1914 THVIRT(D,A,E,N) = THVIRT(D,A,E,N) 1915 & - WYINT(C,N,A,M,D,L) * XPROP(NRHFT+E,NRHFT+C) 1916 & - WXINT(C,N,A,M,D,L) * YPROP(NRHFT+E,NRHFT+C) 1917 & + T0AMP(C,N,A,M,D,L) * XYMAT(NRHFT+E,NRHFT+C) 1918 END DO 1919 DO K = 1, NRHFT 1920 THVIRT(D,A,E,N) = THVIRT(D,A,E,N) 1921 & + THETA(E,K,A,M,D,L) * XPROP(K,N) 1922 & + THETAX(E,K,A,M,D,L) * YPROP(K,N) 1923 END DO 1924 END DO 1925 END DO 1926 END DO 1927 END DO 1928 1929 ! divide by orbital energies and remove forbidden entries 1930 CALL CCSDT_DIVVIR_NODDY(THVIRT,L,M,FREQR2, 1931 & FOCKD,FOCKD(NRHFT+1),NVIRT,NRHFT) 1932 1933 CALL DSCAL(NVIRT*NVIRT*NT1AMX,-1.0D0,THVIRT,1) 1934 1935c -------------------------------------------------- 1936c D(ab) += theta-bar(d-bar L, a M, e n) * 1937c theta(d L, b M, e-bar n-bar) 1938c -------------------------------------------------- 1939 DO N = 1, NRHFT 1940 DO E = 1, NVIRT 1941 CALL DGEMM('T','N',NVIRT,NVIRT,NVIRT, 1942 & ONE,WTILDE(1,1,E,N),NVIRT, 1943 & THVIRT(1,1,E,N),NVIRT, 1944 & ONE,DAB,NVIRT) 1945 END DO 1946 END DO 1947 1948 1949c -------------------------------------------------------- 1950c set up THVIRT^LM(de,an) = theta(d L, e-bar M, a n-bar) 1951c -------------------------------------------------------- 1952 CALL DZERO(THVIRT,NVIRT*NVIRT*NT1AMX) 1953 1954 DO N = 1, NRHFT 1955 DO A = 1, NVIRT 1956 DO E = 1, NVIRT 1957 DO D = 1, NVIRT 1958 DO C = 1, NVIRT 1959 THVIRT(D,E,A,N) = THVIRT(D,E,A,N) - 1960 & WYINT(A,N,C,M,D,L) * XPROP(NRHFT+E,NRHFT+C) - 1961 & WXINT(A,N,C,M,D,L) * YPROP(NRHFT+E,NRHFT+C) 1962 END DO 1963 DO K = 1, NRHFT 1964 THVIRT(D,E,A,N) = THVIRT(D,E,A,N) + 1965 & THETA(E,M,A,K,D,L) * XPROP(K,N) + 1966 & THETAX(E,M,A,K,D,L) * YPROP(K,N) 1967 END DO 1968 END DO 1969 END DO 1970 END DO 1971 END DO 1972 1973 ! divide by orbital energies and remove forbidden entries 1974 CALL CCSDT_DIVVIR_NODDY(THVIRT,L,M,FREQR2, 1975 & FOCKD,FOCKD(NRHFT+1),NVIRT,NRHFT) 1976 1977 CALL DSCAL(NVIRT*NVIRT*NT1AMX,-1.0D0,THVIRT,1) 1978 1979c -------------------------------------------------- 1980c D(ab) += theta-bar(d-bar L, e M, a n) * 1981c theta(d L, e-bar M, b n-bar) 1982c -------------------------------------------------- 1983 DO N = 1, NRHFT 1984 CALL DGEMM('T','N',NVIRT,NVIRT,NVIRT*NVIRT, 1985 & ONE,WTILDE(1,1,1,N),NVIRT*NVIRT, 1986 & THVIRT(1,1,1,N),NVIRT*NVIRT, 1987 & ONE,DAB,NVIRT) 1988 END DO 1989 1990c ------------------------------------------------------ 1991c resort WTILDE^LM(ae,dn) = theta-bar(d-bar n, e M, a L) 1992c ------------------------------------------------------ 1993 DO N = 1, NRHFT 1994 DO A = 1, NVIRT 1995 DO E = 1, NVIRT 1996 DO D = 1, NVIRT 1997 WTILDE(A,E,D,N) = THBAR(D,N,E,M,A,L) 1998 END DO 1999 END DO 2000 END DO 2001 END DO 2002 2003c -------------------------------------------------- 2004c D(ab) += theta-bar(a L, e M, d-bar n) * 2005c theta(b L, e-bar M, d n-bar) 2006c -------------------------------------------------- 2007 CALL DGEMM('N','T',NVIRT,NVIRT,NVIRT*NVIRT*NRHFT, 2008 & ONE,WTILDE,NVIRT, 2009 & THVIRT,NVIRT, 2010 & ONE,DAB,NVIRT) 2011 2012 END IF 2013 2014 END DO 2015 END DO 2016 2017 RETURN 2018 END 2019*---------------------------------------------------------------------* 2020* END OF SUBROUTINE CCSDT_ADENOCC_NODDY * 2021*---------------------------------------------------------------------* 2022*=====================================================================* 2023 SUBROUTINE CCSDT_TXVIR_NODDY(TX,TAMP,XPROP,IOPT, 2024 & NVIRT,NRHFT,NORBT) 2025*---------------------------------------------------------------------* 2026c 2027c IOPT = 1 : TX(de,fn) -= sum_c X_dc*T(ce,fn) (transf. 1 index) 2028c 2029c IOPT = 2 : TX(de,fn) -= sum_c X_ec*T(dc,fn) (transf. 2 index) 2030c 2031c IOPT = 3 : TX(de,fn) -= sum_c X_fc*T(de,cn) (transf. 3 index) 2032c 2033c IOPT = 0 : transform all three indeces 2034c 2035*---------------------------------------------------------------------* 2036 IMPLICIT NONE 2037 2038 INTEGER NVIRT, NRHFT, NORBT, IOPT, D, E, F, N, C 2039 2040#if defined (SYS_CRAY) 2041 REAL TAMP(NVIRT,NVIRT,NVIRT,NRHFT) 2042 REAL TX(NVIRT,NVIRT,NVIRT,NRHFT) 2043 REAL XPROP(NORBT,NORBT) 2044#else 2045 DOUBLE PRECISION TAMP(NVIRT,NVIRT,NVIRT,NRHFT) 2046 DOUBLE PRECISION TX(NVIRT,NVIRT,NVIRT,NRHFT) 2047 DOUBLE PRECISION XPROP(NORBT,NORBT) 2048#endif 2049 2050 DO N = 1, NRHFT 2051 DO F = 1, NVIRT 2052 2053 IF (IOPT.EQ.0 .OR. IOPT.EQ.1) THEN 2054 DO D = 1, NVIRT 2055 DO E = 1, NVIRT 2056 DO C = 1, NVIRT 2057 TX(D,E,F,N) = TX(D,E,F,N) - 2058 & XPROP(NRHFT+D,NRHFT+C) * TAMP(C,E,F,N) 2059 END DO 2060 END DO 2061 END DO 2062 END IF 2063 2064 IF (IOPT.EQ.0 .OR. IOPT.EQ.2) THEN 2065 DO D = 1, NVIRT 2066 DO E = 1, NVIRT 2067 DO C = 1, NVIRT 2068 TX(D,E,F,N) = TX(D,E,F,N) - 2069 & XPROP(NRHFT+E,NRHFT+C) * TAMP(D,C,F,N) 2070 END DO 2071 END DO 2072 END DO 2073 END IF 2074 2075 IF (IOPT.EQ.0 .OR. IOPT.EQ.3) THEN 2076 DO D = 1, NVIRT 2077 DO E = 1, NVIRT 2078 DO C = 1, NVIRT 2079 TX(D,E,F,N) = TX(D,E,F,N) - 2080 & XPROP(NRHFT+F,NRHFT+C) * TAMP(D,E,C,N) 2081 END DO 2082 END DO 2083 END DO 2084 END IF 2085 2086 END DO 2087 END DO 2088 2089 RETURN 2090 END 2091*---------------------------------------------------------------------* 2092*=====================================================================* 2093 SUBROUTINE CCSDT_DIVVIR_NODDY(TAMP,L,M,FREQ,FCKDOCC,FCKDVIR, 2094 & NVIRT,NRHFT) 2095*---------------------------------------------------------------------* 2096c 1) T^LM(defn) = T^LM(defn) / (eps_dl + eps_em + eps_fn - freq) 2097c 2) remove forbidden elements 2098*---------------------------------------------------------------------* 2099 IMPLICIT NONE 2100 2101 INTEGER NVIRT,NRHFT,L,M,N,D,E,F 2102 2103#if defined (SYS_CRAY) 2104 REAL TAMP(NVIRT,NVIRT,NVIRT,NRHFT) 2105 REAL FCKDOCC(NRHFT), FCKDVIR(NVIRT), FREQ, FLMN 2106#else 2107 DOUBLE PRECISION TAMP(NVIRT,NVIRT,NVIRT,NRHFT) 2108 DOUBLE PRECISION FCKDOCC(NRHFT), FCKDVIR(NVIRT), FREQ, FLMN 2109#endif 2110 2111 DO N = 1, NRHFT 2112 FLMN = FCKDOCC(L) + FCKDOCC(M) + FCKDOCC(N) + FREQ 2113 DO F = 1, NVIRT 2114 DO E = 1, NVIRT 2115 DO D = 1, NVIRT 2116 TAMP(D,E,F,N) = TAMP(D,E,F,N) / 2117 & (FCKDVIR(D) + FCKDVIR(E) + FCKDVIR(F) - FLMN) 2118 END DO 2119 END DO 2120 END DO 2121 END DO 2122 2123 IF (L.EQ.M) THEN 2124 DO F = 1, NVIRT 2125 DO E = 1, NVIRT 2126 DO D = 1, NVIRT 2127 TAMP(D,E,F,L) = 0.0d0 2128 END DO 2129 END DO 2130 END DO 2131 END IF 2132 2133 DO N = 1, NRHFT 2134 DO D = 1, NVIRT 2135 TAMP(D,D,D,N) = 0.0d0 2136 END DO 2137 END DO 2138 2139 RETURN 2140 END 2141*---------------------------------------------------------------------* 2142*=====================================================================* 2143 SUBROUTINE CCSDT_WXYVINT_NODDY(WXYV,THETAY,THETAX,YPROP,XPROP, 2144 & NVIRT,NRHFT,NORBT) 2145*---------------------------------------------------------------------* 2146c build wxy(em,fn,dl) = wxy(em,fn,dl) + 2147c thetay(ek,fn,dl) * xprop(km) + thetax(ek,fn,dl) * yprop(km) + 2148c thetay(em,fk,dl) * xprop(kn) + thetax(em,fk,dl) * yprop(kn) + 2149c thetay(em,fn,dk) * xprop(kl) + thetax(em,fn,dk) * yprop(kl) 2150*---------------------------------------------------------------------* 2151 IMPLICIT NONE 2152 2153 INTEGER NRHFT,NVIRT,NORBT, E, M, F, N, D, L, K 2154 2155#if defined (SYS_CRAY) 2156 REAL WXYV(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2157 REAL THETAY(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2158 REAL THETAX(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2159 REAL XPROP(NORBT,NORBT), YPROP(NORBT,NORBT) 2160#else 2161 DOUBLE PRECISION WXYV(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2162 DOUBLE PRECISION THETAY(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2163 DOUBLE PRECISION THETAX(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2164 DOUBLE PRECISION XPROP(NORBT,NORBT), YPROP(NORBT,NORBT) 2165#endif 2166 2167 DO D = 1, NVIRT 2168 DO L = 1, NRHFT 2169 DO F = 1, NVIRT 2170 DO N = 1, NRHFT 2171 DO E = 1, NVIRT 2172 DO M = 1, NRHFT 2173 2174 DO K = 1, NRHFT 2175 WXYV(E,M,F,N,D,L) = WXYV(E,M,F,N,D,L) + 2176 & THETAY(E,K,F,N,D,L) * XPROP(K,M) + 2177 & THETAX(E,K,F,N,D,L) * YPROP(K,M) + 2178 & THETAY(E,M,F,K,D,L) * XPROP(K,N) + 2179 & THETAX(E,M,F,K,D,L) * YPROP(K,N) + 2180 & THETAY(E,M,F,N,D,K) * XPROP(K,L) + 2181 & THETAX(E,M,F,N,D,K) * YPROP(K,L) 2182 END DO 2183 2184 END DO 2185 END DO 2186 END DO 2187 END DO 2188 END DO 2189 END DO 2190 2191 RETURN 2192 END 2193*---------------------------------------------------------------------* 2194*=====================================================================* 2195 SUBROUTINE CCSDT_T32INT_NODDY(LISTR,IDLSTR,ISYMR,FREQR, 2196 * LISTX,IDLSTX,ISYMX, 2197 * WXINT,THETAX,HAVEX,XPROP, 2198 * LISTY,IDLSTY,ISYMY, 2199 * WYINT,THETAY,HAVEY,YPROP, 2200 * T3AMP0,HAVEXY,XYMAT, 2201 * T1AMXY,T2AMXY,TAMXY3,WXYOINT,WXYVINT, 2202 * XLAMDP0,XLAMDH0,FOCK0,FOCKD, 2203 * LUDELD,FNDELD,LUDKBC,FNDKBC,LUCKJD,FNCKJD, 2204 * LOCDBG,HAVETB3T3, 2205 * WORK,LWORK) 2206*---------------------------------------------------------------------* 2207* calculate a number of intermediates needed in CCSDT_ADEN_NODDY 2208* to represent the second-order amplitudes T3^XY 2209*---------------------------------------------------------------------* 2210 IMPLICIT NONE 2211#include "priunit.h" 2212#include "ccorb.h" 2213#include "ccsdsym.h" 2214#include "dummy.h" 2215 2216 INTEGER ISYM0 2217 PARAMETER (ISYM0 = 1) 2218 2219 LOGICAL LOCDBG, HAVETB3T3, HAVEX, HAVEY, HAVEXY 2220 CHARACTER LISTX*3, LISTY*3, LISTR*3, 2221 & FNDELD*(*), FNDKBC*(*), FNCKJD*(*) 2222 INTEGER IDLSTX, ISYMX, IDLSTY, ISYMY, IDLSTR, ISYMR, LWORK 2223 INTEGER LUDELD, LUDKBC, LUCKJD 2224 2225#if defined (SYS_CRAY) 2226 REAL WXINT(*), THETAX(*), XPROP(*) 2227 REAL WYINT(*), THETAY(*), YPROP(*) 2228 REAL T3AMP0(*), XYMAT(*), WORK(*) 2229 REAL TAMXY3(*), WXYOINT(*), WXYVINT(*) 2230 REAL T1AMXY(*), T2AMXY(*) 2231 REAL XLAMDP0(*), XLAMDH0(*), FOCK0(*), FOCKD(*) 2232 REAL FREQR, DDOT 2233#else 2234 DOUBLE PRECISION WXINT(*), THETAX(*), XPROP(*) 2235 DOUBLE PRECISION WYINT(*), THETAY(*), YPROP(*) 2236 DOUBLE PRECISION T3AMP0(*), XYMAT(*), WORK(*) 2237 DOUBLE PRECISION TAMXY3(*), WXYOINT(*), WXYVINT(*) 2238 DOUBLE PRECISION T1AMXY(*), T2AMXY(*) 2239 DOUBLE PRECISION XLAMDP0(*), XLAMDH0(*), FOCK0(*), FOCKD(*) 2240 DOUBLE PRECISION FREQR, DDOT 2241#endif 2242 2243 CHARACTER*10 MODEL 2244 INTEGER KEND1, KLAMPX, KLAMHX, KLAMPY, KLAMHY, LWRK1, KEND2, 2245 & KINT1SXY, KINT2SXY, KINT1SX, KINT2SX, KINT1SY, KINT2SY, 2246 & KT2AMP0, KINT1S0, KINT2S0, KSCR1, KTAMPX3, KTAMPY3, 2247 & LWRK2, IOPT, LUSIFC, KTHEOCC, KWXYO, KWXYV, 2248 & KT2AMPX, KT2AMPY, KTHVIRT, KTXY3A, KTXY3B 2249 2250 CALL QENTER('CCT32INT') 2251*---------------------------------------------------------------------* 2252* allocate memory for local arrays: 2253*---------------------------------------------------------------------* 2254 KSCR1 = 1 2255 KEND1 = KSCR1 + NT1AMX 2256 2257 KLAMPX = KEND1 2258 KLAMHX = KLAMPX + NLAMDT 2259 KEND1 = KLAMHX + NLAMDT 2260 2261 KLAMPY = KEND1 2262 KLAMHY = KLAMPY + NLAMDT 2263 KEND1 = KLAMHY + NLAMDT 2264 2265 KT2AMP0 = KEND1 2266 KEND1 = KT2AMP0 + NT1AMX*NT1AMX 2267 2268 KT2AMPY = KEND1 2269 KT2AMPX = KT2AMPY + NT1AMX*NT1AMX 2270 KEND1 = KT2AMPX + NT1AMX*NT1AMX 2271 2272 KTXY3A = KEND1 2273 KTXY3B = KTXY3A + NT1AMX*NT1AMX*NT1AMX 2274 KEND1 = KTXY3B + NT1AMX*NT1AMX*NT1AMX 2275 2276 IF (HAVETB3T3) THEN 2277 KTAMPX3 = KEND1 2278 KTAMPY3 = KTAMPX3 + NT1AMX*NT1AMX*NT1AMX 2279 KEND1 = KTAMPY3 + NT1AMX*NT1AMX*NT1AMX 2280 END IF 2281 2282 LWRK1 = LWORK - KEND1 2283 IF (LWRK1 .LT. 0) CALL QUIT('Out of memory in CCSDT_T32INT_NODDY') 2284 2285*---------------------------------------------------------------------* 2286* Read SCF orbital energies from file: 2287*---------------------------------------------------------------------* 2288 LUSIFC = -1 2289 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 2290 & .FALSE.) 2291 REWIND LUSIFC 2292 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 2293 READ (LUSIFC) 2294 READ (LUSIFC) (FOCKD(I), I=1,NORBT) 2295 CALL GPCLOSE(LUSIFC,'KEEP') 2296 2297*---------------------------------------------------------------------* 2298* Get all the intermediates needed to compute the contributions 2299* from the A{x} t^y and A{y} t^x to T^xy, i.e. 2300* the w^x, theta^x, w^y, theta^y intermediates and the property 2301* integrals and the one-index transformed property integrals 2302*---------------------------------------------------------------------* 2303 CALL CCSDT_T32_AAPREP_NODDY( 2304 * LISTX,IDLSTX,ISYMX,WXINT,THETAX,XPROP, 2305 * WORK(KLAMPX),WORK(KLAMHX),WORK(KTAMPX3), 2306 * LISTY,IDLSTY,ISYMY,WYINT,THETAY,YPROP, 2307 * WORK(KLAMPY),WORK(KLAMHY),WORK(KTAMPY3), 2308 * T3AMP0,XYMAT,XLAMDP0,XLAMDH0,FOCK0, 2309 * LUDELD,FNDELD,LUDKBC,FNDKBC,LUCKJD,FNCKJD, 2310 * LOCDBG,HAVETB3T3,WORK(KEND1),LWRK1) 2311 2312*---------------------------------------------------------------------* 2313* Compute the "easy" contributions to T^XY_3 from the B matrix, 2314* which also in the real code will be accesible explicitly 2315* (or as W-like intermediate) 2316* <mu_3| [H^AB,T^0_2] + [H^A,T^B_2] + [H^B,T^A_2] |HF> 2317*---------------------------------------------------------------------* 2318 KINT1SXY = KEND1 2319 KINT2SXY = KINT1SXY + NT1AMX*NVIRT*NVIRT 2320 KEND2 = KINT2SXY + NRHFT*NRHFT*NT1AMX 2321 2322 KINT1SX = KEND2 2323 KINT2SX = KINT1SX + NT1AMX*NVIRT*NVIRT 2324 KEND2 = KINT2SX + NRHFT*NRHFT*NT1AMX 2325 2326 KINT1SY = KEND2 2327 KINT2SY = KINT1SY + NT1AMX*NVIRT*NVIRT 2328 KEND2 = KINT2SY + NRHFT*NRHFT*NT1AMX 2329 2330 LWRK2 = LWORK - KEND2 2331 IF (LWRK2 .LT. NT2AMX) THEN 2332 CALL QUIT('Insufficient space in CCSDT_T32INT_NODDY') 2333 ENDIF 2334 2335 ! read T^0 doubles amplitudes from file and square up 2336 IOPT = 2 2337 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KEND2)) 2338 CALL CC_T2SQ(WORK(KEND2),WORK(KT2AMP0),ISYM0) 2339 2340 2341 CALL CCSDT_INTS1_NODDY(.TRUE.,WORK(KINT1SX),WORK(KINT2SX), 2342 & .FALSE.,DUMMY,DUMMY, 2343 & XLAMDP0,XLAMDH0, 2344 & WORK(KLAMPX),WORK(KLAMHX), 2345 & WORK(KEND2),LWRK2) 2346 2347 CALL CCSDT_INTS1_NODDY(.TRUE.,WORK(KINT1SY),WORK(KINT2SY), 2348 & .FALSE.,DUMMY,DUMMY, 2349 & XLAMDP0,XLAMDH0, 2350 & WORK(KLAMPY),WORK(KLAMHY), 2351 & WORK(KEND2),LWRK2) 2352 2353 CALL CCSDT_INTS2_NODDY(WORK(KINT1SXY),WORK(KINT2SXY), 2354 & XLAMDP0,XLAMDH0, 2355 & WORK(KLAMPX),WORK(KLAMHX), 2356 & WORK(KLAMPY),WORK(KLAMHY), 2357 & WORK(KEND2),LWRK2) 2358 2359 CALL DZERO(WORK(KTXY3B),NT1AMX*NT1AMX*NT1AMX) 2360 2361 ! note: this routine will overwrite WORK(KT2AMP0) 2362 CALL CCSDT_B3AM(WORK(KTXY3B), 2363 & WORK(KINT1SXY),WORK(KINT2SXY),FOCKD, 2364 & LISTX,IDLSTX,WORK(KINT1SX),WORK(KINT2SX), 2365 & LISTY,IDLSTY,WORK(KINT1SY),WORK(KINT2SY), 2366 & WORK(KSCR1),WORK(KT2AMP0),WORK(KEND2),LWRK2) 2367 2368 2369 ! devide by orbital energy denominators and fix the sign 2370 ! such that it corresponds to that of T^XY in CCSDT_T32_NODDY 2371 CALL CCSDT_3AM(WORK(KTXY3B),FREQR,WORK(KSCR1),FOCKD, 2372 & .FALSE.,DUMMY,.FALSE.,WORK(KEND1)) 2373 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KTXY3B),1) 2374 2375 CALL CCSDT_CLEAN_T3(WORK(KTXY3B),NT1AMX,NVIRT,NRHFT) 2376 2377 2378*---------------------------------------------------------------------* 2379* Precompute in a similar way the "easy" contributions to T^XY_3 2380* comming from the Jacobian matrix: 2381* <mu_3| [[H,T^AB_1],T^0_2] + [H,T^AB_2] |HF> 2382*---------------------------------------------------------------------* 2383 KINT1S0 = KEND1 2384 KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT 2385 KEND2 = KINT2S0 + NRHFT*NRHFT*NT1AMX 2386 2387 KINT1SXY = KEND2 2388 KINT2SXY = KINT1SXY + NT1AMX*NVIRT*NVIRT 2389 KEND2 = KINT2SXY + NRHFT*NRHFT*NT1AMX 2390 2391 LWRK2 = LWORK - KEND2 2392 IF (LWRK2 .LT. NT2AMX) THEN 2393 CALL QUIT('Insufficient space in CCSDT_T32INT_NODDY') 2394 ENDIF 2395 2396 ! read T^0 doubles amplitudes from file and square up 2397 IOPT = 2 2398 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KEND2)) 2399 CALL CC_T2SQ(WORK(KEND2),WORK(KT2AMP0),ISYM0) 2400 2401 CALL CCSDT_READ_NODDY(.FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 2402 & .FALSE.,DUMMY,DUMMY, 2403 & .TRUE.,WORK(KINT1S0),WORK(KINT2S0), 2404 & .FALSE.,DUMMY,DUMMY, 2405 & .FALSE.,DUMMY,DUMMY,DUMMY,DUMMY, 2406 & NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX) 2407 2408 CALL DZERO(WORK(KTXY3A),NT1AMX*NT1AMX*NT1AMX) 2409 2410 CALL CCSDT_A3AM(WORK(KTXY3A),T1AMXY,T2AMXY, 2411 & DUMMY,ISYMR,WORK(KT2AMP0),T3AMP0, 2412 & WORK(KINT1S0), WORK(KINT2S0), 2413 & WORK(KINT1SXY),WORK(KINT2SXY), 2414 & FOCKD,XLAMDP0,XLAMDH0, 2415 & DUMMY,DUMMY,WORK(KSCR1),WORK(KEND2),LWRK2) 2416 2417 2418 ! devide by orbital energy denominators and fix the sign 2419 ! such that it corresponds to that of T^XY in CCSDT_T32_NODDY 2420 CALL CCSDT_3AM(WORK(KTXY3A),FREQR,WORK(KSCR1),FOCKD, 2421 & .FALSE.,DUMMY,.FALSE.,WORK(KEND1)) 2422 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KTXY3A),1) 2423 2424 CALL CCSDT_CLEAN_T3(WORK(KTXY3A),NT1AMX,NVIRT,NRHFT) 2425 2426 2427*---------------------------------------------------------------------* 2428* precalculate modified ^{XY}W intermediates 2429*---------------------------------------------------------------------* 2430 KTHEOCC = KEND1 2431 KWXYO = KTHEOCC + NT1AMX*NRHFT*NRHFT 2432 KWXYV = KWXYO + NT1AMX*NRHFT*NRHFT 2433 KEND2 = KWXYV + NT1AMX*NRHFT*NRHFT 2434 2435 LWRK2 = LWORK - KEND2 2436 IF (LWRK2.LT.NT2AMX) 2437 & CALL QUIT('Out of memory in CCSDT_T32INT_NODDY') 2438 2439 ! read T^X doubles amplitudes from file and square up 2440 IOPT = 2 2441 Call CC_RDRSP(LISTX,IDLSTX,ISYMX,IOPT,MODEL,DUMMY,WORK(KEND2)) 2442 Call CCLR_DIASCL(WORK(KEND2),2.0D0,ISYMX) 2443 CALL CC_T2SQ(WORK(KEND2),WORK(KT2AMPX),ISYMX) 2444 2445 ! read T^X doubles amplitudes from file and square up 2446 IOPT = 2 2447 Call CC_RDRSP(LISTY,IDLSTY,ISYMY,IOPT,MODEL,DUMMY,WORK(KEND2)) 2448 Call CCLR_DIASCL(WORK(KEND2),2.0D0,ISYMY) 2449 CALL CC_T2SQ(WORK(KEND2),WORK(KT2AMPY),ISYMY) 2450 2451 IF (HAVETB3T3) THEN 2452 CALL DZERO(TAMXY3,NT1AMX*NT1AMX*NT1AMX) 2453 END IF 2454 2455 CALL CCSDT_O32VIR_NODDY(WXINT,XPROP,WORK(KT2AMPX),HAVEX, 2456 & WYINT,YPROP,WORK(KT2AMPY),HAVEY, 2457 & T3AMP0,WORK(KT2AMP0), 2458 & XYMAT,HAVEXY, 2459 & WORK(KTHEOCC),WORK(KWXYV),WORK(KWXYO), 2460 & TAMXY3,HAVETB3T3, 2461 & WXYOINT,.TRUE.,WXYVINT,.TRUE.) 2462 2463 2464 2465 CALL CCSDT_CLEAN_T3(WXYOINT,NT1AMX,NVIRT,NRHFT) 2466 CALL CCSDT_CLEAN_T3(WXYVINT,NT1AMX,NVIRT,NRHFT) 2467 2468c ------------------------------------------------------------- 2469c add some more contributions to WXYVINT to obtain full 2470c theta(e-bar m-bar, f n-bar, d l-bar) 2471c ------------------------------------------------------------- 2472 CALL CCSDT_WXYVINT_NODDY(WXYVINT,THETAY,THETAX,YPROP,XPROP, 2473 & NVIRT,NRHFT,NORBT) 2474 CALL CCSDT_CLEAN_T3(WXYVINT,NT1AMX,NVIRT,NRHFT) 2475 2476 2477c ------------------------------------------------------------- 2478c devide by orbital energy denominators and fix the sign 2479c such that it corresponds to that of T^XY in CCSDT_T32_NODDY 2480c ------------------------------------------------------------- 2481 CALL CCSDT_3AM(WXYOINT,FREQR,WORK(KSCR1),FOCKD, 2482 & .FALSE.,DUMMY,.FALSE.,WORK(KEND1)) 2483 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WXYOINT,1) 2484 2485 CALL CCSDT_3AM(WXYVINT,FREQR,WORK(KSCR1),FOCKD, 2486 & .FALSE.,DUMMY,.FALSE.,WORK(KEND1)) 2487 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WXYVINT,1) 2488 2489 2490c ---------------------------------------------------------------- 2491c add to WXYOINT 1/3 of the contribution from the B T^X T^Y and 2492c A T^XY transformation such, that we can later generate from 2493c this intermediate the terms with double-bars on occupied indeces 2494c ---------------------------------------------------------------- 2495 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,1.0D0/3.0D0, 2496 & WORK(KTXY3B),1,WXYOINT,1) 2497 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,1.0D0/3.0D0, 2498 & WORK(KTXY3A),1,WXYOINT,1) 2499 2500 2501*---------------------------------------------------------------------* 2502* if the HAVETB3T3 test option is set, complete the T^XY_3 vector: 2503*---------------------------------------------------------------------* 2504 IF (HAVETB3T3) THEN 2505 KTHVIRT = KEND1 2506 KWXYV = KTHVIRT + NT1AMX*NVIRT*NVIRT 2507 KWXYO = KWXYV + NT1AMX*NVIRT*NVIRT 2508 KEND2 = KWXYO + NT1AMX*NVIRT*NVIRT 2509 2510 LWRK2 = LWORK - KEND2 2511 IF (LWRK2.LT.0) 2512 & CALL QUIT('Out of memory in CCSDT_T32INT_NODDY') 2513 2514 CALL CCSDT_O32OCC_NODDY(THETAX,XPROP,HAVEX, 2515 & THETAY,YPROP,HAVEY, 2516 & TAMXY3,WORK(KTHVIRT), 2517 & WORK(KWXYV),WORK(KWXYO),.TRUE.,.TRUE., 2518 & DUMMY,.FALSE.,DUMMY,.FALSE.) 2519 CALL CCSDT_CLEAN_T3(TAMXY3,NT1AMX,NVIRT,NRHFT) 2520 2521 2522 ! devide by orbital energy denominators and fix the sign 2523 ! such that it corresponds to that of T^XY in CCSDT_T32_NODDY 2524 CALL CCSDT_3AM(TAMXY3,FREQR,WORK(KSCR1),FOCKD, 2525 & .FALSE.,DUMMY,.FALSE.,WORK(KEND1)) 2526 CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,TAMXY3,1) 2527 2528 2529 ! add contributions from B T^X T^Y transformation 2530 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,1.0D0,WORK(KTXY3B),1,TAMXY3,1) 2531 2532 2533 ! add contributions from A T^XY transformation 2534 CALL DAXPY(NT1AMX*NT1AMX*NT1AMX,1.0D0,WORK(KTXY3A),1,TAMXY3,1) 2535 2536 2537C WRITE(LUPRI,*) 'CCSDT_T32INT_NODDY> T^XY:',LISTR,IDLSTR,FREQR 2538C CALL CCSDT_CLEAN_T3(TAMXY3,NT1AMX,NVIRT,NRHFT) 2539C CALL PRINT_PT3_NODDY(TAMXY3) 2540 END IF 2541 2542 CALL QEXIT('CCT32INT') 2543 RETURN 2544 END 2545*=====================================================================* 2546*---------------------------------------------------------------------* 2547C /* Deck ccsdt_gwbic */ 2548*=====================================================================* 2549 SUBROUTINE CCSDT_GWBIC(LISTL,IDLSTL,WBAR,THETA,SPLIT_WT, 2550 & XLAMDP0,XLAMDH0,FOCK0, 2551 & LUTOC,FNTOC,LU3VI,FN3VI,LUDKBC3,FNDKBC3, 2552 & LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X, 2553 & WORK,LWORK) 2554*---------------------------------------------------------------------* 2555* 2556* Purpose: set up WBAR as N^6 array in core for use in noddy code 2557* 2558* W3BAR^Y = ( <L2|[Y,tau3]|HF> + <L3|[Y^,tau3]|HF> 2559* + <L2|[H^Y,tau3]|HF> 2560* + <L1^Y + L2^Y|[H^,tau3]|HF> )/ (w+epsiln(tau3)) 2561* 2562* 2563* Christof Haettig, spring 2003 based CC3_XI_DEN. 2564* 2565*=====================================================================* 2566C 2567 IMPLICIT NONE 2568C 2569#include "priunit.h" 2570#include "dummy.h" 2571#include "iratdef.h" 2572#include "ccsdsym.h" 2573#include "ccorb.h" 2574#include "ccsdinp.h" 2575#include "ccinftap.h" 2576#include "inftap.h" 2577#include "cc3t3d.h" 2578#include "ccl1rsp.h" 2579#include "ccr1rsp.h" 2580 2581C 2582 LOGICAL LOCDBG 2583 PARAMETER (LOCDBG = .FALSE.) 2584C 2585 INTEGER ISYM0 2586 PARAMETER(ISYM0 = 1) 2587 CHARACTER LISTL*3, LISTR*3, MODEL*10, LABELY*8 2588 2589 CHARACTER*5 FN3FOP 2590 CHARACTER*6 FN3FOP2 2591 CHARACTER*8 FN3VI2 2592 2593 PARAMETER (FN3VI2 = 'CC3_VI12') 2594 PARAMETER (FN3FOP = 'PTFOP' , FN3FOP2= 'PTFOP2') 2595 2596 CHARACTER*(*) FNTOC,FN3VI,FNDKBC3,FN3FOPX,FN3FOP2X 2597 2598 LOGICAL LORX, SPLIT_WT 2599 INTEGER LWORK, IDLSTL, ISYMY, IDLSTR, LU3VI2, LU3FOPX, LU3FOP2X, 2600 & LUTOC, LU3VI, LUDKBC3, LU3FOP, LU3FOP2, 2601 & KFOCKD,KFCKBA,KL1AM,KL2TP,KEND0,LWRK0, 2602 & KL1BY, KL2TPBY, KFOCKY, KEND1, LWRK1, 2603 & IOPT, LENGTH, IOFF, IRREP, ISYM, IERR, 2604 & ISYMD, ISYML, ISYMDL, ISYMC, ISYMK, KOFF1, KOFF2, KXIAJB, 2605 & KTROC01, KTROC21, KTROC03, KTROC23, KT3OC0, KT3OC02, 2606 & KLAMPY, KLAMHY, KFCKYCK, KINTOC, KT1Y, 2607 & KTROCCG0,KTROCCL0,KTROCCGY,KTROCCLY, KEND2,LWRK2, 2608 & ISYCKB, ISYCKBY, 2609 & KTRVI4, KTRVI5, KTRVI6, KTRVI7, KEND3, LWRK3, 2610 & KTRVI11,KTRVI12,KTRVI13, 2611 & KTRVI14,KTRVI15,KTRVI16,KTRVI17,KTRVI18,KTRVI19,KTRVI20, 2612 & KVGDKCB0,KVGDKBC0,KVLDKCB0,KVLDKBC0, KINTVI, 2613 & KVGDKCBY,KVGDKBCY,KVLDKCBY,KVLDKBCY, KEND4, LWRK4, 2614 & ISYALJ, ISYALJ2, ISYALJBY, ISYALJDY, ISYMBD, ISCKIJ, 2615 & ISWMAT, ISYCKD, ISAIBJ, ISYMJ, ISYMAI, ISYAIL, 2616 & KSMAT2, KUMAT2, KDIAG, KDIAGW, KINDSQ, KINDSQW, KINDEX, 2617 & KINDEX2, KINDEXBY, KINDEXDY, KTMAT,KWMAT, KSMAT4,KUMAT4, 2618 & KEND5, LWRK5, LENSQ, LENSQW, ISYMB, ISYMBJ 2619 2620#if defined (SYS_CRAY) 2621 REAL XLAMDP0(*), XLAMDH0(*), FOCK0(*) 2622 REAL FREQ, HALF 2623 REAL WORK(LWORK) 2624 REAL WBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2625 REAL THETA(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2626#else 2627 DOUBLE PRECISION XLAMDP0(*), XLAMDH0(*), FOCK0(*) 2628 DOUBLE PRECISION FREQ, HALF 2629 DOUBLE PRECISION WORK(LWORK) 2630 DOUBLE PRECISION WBAR(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2631 DOUBLE PRECISION THETA(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 2632#endif 2633 PARAMETER(HALF = 0.5D0) 2634 2635* external functions: 2636 INTEGER IR1TAMP 2637 2638*---------------------------------------------------------------------* 2639* some initializations: 2640*---------------------------------------------------------------------* 2641 CALL QENTER('CCSDT_GWBIC') 2642 2643 CALL DZERO(WBAR,NT1AMX*NT1AMX*NT1AMX) 2644 2645 IF (LISTL(1:3).EQ.'L1 ') THEN 2646 ! get symmetry, frequency and integral label from common blocks 2647 ! defined in ccl1rsp.h 2648 ISYMY = ISYLRZ(IDLSTL) 2649 FREQ = FRQLRZ(IDLSTL) 2650 LABELY = LRZLBL(IDLSTL) 2651c 2652 LORX = LORXLRZ(IDLSTL) 2653 IF (LORX) CALL QUIT('NO ORBITAL RELAX. IN CCSDT_GWBIC') 2654 ELSE 2655 CALL QUIT('Illegal left vector type in CCSDT_GWBIC:'//LISTL) 2656 END IF 2657 2658 ! find right response vector: 2659 IF (LISTL(1:3).EQ.'L1 ') THEN 2660 LISTR = 'R1 ' 2661 IDLSTR = IR1TAMP(LABELY,LORX,FREQ,ISYMY) 2662 ELSE 2663 CALL QUIT('Unknown list in CCSDT_GWBIC') 2664 END IF 2665 2666*---------------------------------------------------------------------* 2667* open files: 2668*---------------------------------------------------------------------* 2669 LU3VI2 = -1 2670 LU3FOP = -1 2671 LU3FOP2 = -1 2672 2673 CALL WOPEN2(LU3VI2,FN3VI2,64,0) 2674 CALL WOPEN2(LU3FOP,FN3FOP,64,0) 2675 CALL WOPEN2(LU3FOP2,FN3FOP2,64,0) 2676 2677*---------------------------------------------------------------------* 2678* initial allocations, orbital energy, fock matrix and L2 : 2679*---------------------------------------------------------------------* 2680 KFOCKD = 1 2681 KFCKBA = KFOCKD + NORBTS 2682 KL1AM = KFCKBA + NT1AMX 2683 KL2TP = KL1AM + NT1AM(ISYM0) 2684 KEND0 = KL2TP + NT2SQ(ISYM0) 2685 LWRK0 = LWORK - KEND0 2686 2687 KL1BY = KEND0 2688 KL2TPBY = KL1BY + NT1AM(ISYMY) 2689 KFOCKY = KL2TPBY + NT2SQ(ISYMY) 2690 KEND0 = KFOCKY + N2BST(ISYMY) 2691 2692 KEND1 = KEND0 2693 LWRK1 = LWORK - KEND1 2694 2695C---------------------------------------------------------------------* 2696C Matrix with property integrals in MO basis: 2697C---------------------------------------------------------------------* 2698 ! read property integrals from file: 2699 CALL CCPRPAO(LABELY,.TRUE.,WORK(KFOCKY),IRREP,ISYM,IERR, 2700 & WORK(KEND1),LWRK1) 2701 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISYMY)) THEN 2702 CALL QUIT('CCSDT_GWBIC: error reading oper. '//LABELY) 2703 ELSE IF (IERR.LT.0) THEN 2704 CALL DZERO(WORK(KFOCKY),N2BST(ISYMY)) 2705 END IF 2706 2707 ! transform property integrals to Lambda-MO basis 2708 CALL CC_FCKMO(WORK(KFOCKY),XLAMDP0,XLAMDH0, 2709 & WORK(KEND1),LWRK1,ISYMY,1,1) 2710 2711C------------------------------------- 2712C Read L2 amplitudes 2713C------------------------------------- 2714 IOPT = 3 2715 CALL CC_RDRSP('L0 ',0,ISYM0,IOPT,MODEL, 2716 * WORK(KL1AM),WORK(KL2TP)) 2717 CALL CC_T2SQ(WORK(KL2TP),WORK(KEND1),ISYM0) 2718 CALL CC3_T2TP(WORK(KL2TP),WORK(KEND1),ISYM0) 2719 2720C------------------------------------- 2721C Read L1BY and L2BY multipliers 2722C------------------------------------- 2723 IOPT = 3 2724 CALL CC_RDRSP(LISTL,IDLSTL,ISYMY,IOPT,MODEL, 2725 & WORK(KL1BY),WORK(KL2TPBY)) 2726 CALL CC_T2SQ(WORK(KL2TPBY),WORK(KEND1),ISYMY) 2727 CALL CC3_T2TP(WORK(KL2TPBY),WORK(KEND1),ISYMY) 2728 2729C------------------------------------- 2730C Read canonical orbital energies: 2731C------------------------------------- 2732 LUSIFC = -1 2733 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 2734 & .FALSE.) 2735 REWIND LUSIFC 2736 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 2737 READ (LUSIFC) 2738 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 2739 CALL GPCLOSE(LUSIFC,'KEEP') 2740 2741 ! Delete frozen orbitals in Fock diagonal. 2742 IF (FROIMP .OR. FROEXP) 2743 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK0) 2744 2745C----------------------------------------------------- 2746C Sort the fock matrix 2747C----------------------------------------------------- 2748 DO ISYMC = 1,NSYM 2749 ISYMK = MULD2H(ISYMC,ISYM0) 2750 DO K = 1,NRHF(ISYMK) 2751 DO C = 1,NVIR(ISYMC) 2752 KOFF1 = IFCVIR(ISYMK,ISYMC) + 2753 * NORB(ISYMK)*(C - 1) + K 2754 KOFF2 = IT1AM(ISYMC,ISYMK) 2755 * + NVIR(ISYMC)*(K - 1) + C 2756 WORK(KFCKBA-1+KOFF2) = FOCK0(KOFF1) 2757 ENDDO 2758 ENDDO 2759 ENDDO 2760C 2761C----------------------------- 2762C Memory allocation. 2763C----------------------------- 2764 KXIAJB = KEND1 2765 KEND1 = KXIAJB + NT2AM(ISYM0) 2766 2767 KTROC01 = KEND1 2768 KTROC21 = KTROC01 + NTRAOC(ISYM0) 2769 KTROC03 = KTROC21 + NTRAOC(ISYM0) 2770 KTROC23 = KTROC03 + NTRAOC(ISYM0) 2771 KT3OC0 = KTROC23 + NTRAOC(ISYM0) 2772 KT3OC02 = KT3OC0 + NTRAOC(ISYM0) 2773 KLAMPY = KT3OC02 + NTRAOC(ISYM0) 2774 KLAMHY = KLAMPY + NLAMDT 2775 KEND1 = KLAMHY + NLAMDT 2776 LWRK1 = LWORK - KEND1 2777 2778 KFCKYCK = KEND1 2779 KTROCCG0 = KFCKYCK + NT1AM(ISYMY) 2780 KTROCCL0 = KTROCCG0 + NTRAOC(ISYM0) 2781 KTROCCGY = KTROCCL0 + NTRAOC(ISYM0) 2782 KTROCCLY = KTROCCGY + NTRAOC(ISYMY) 2783 KEND1 = KTROCCLY + NTRAOC(ISYMY) 2784 LWRK1 = LWORK - KEND1 2785 2786 KINTOC = KEND1 2787 KT1Y = KINTOC + MAX(NTOTOC(ISYM0),NTOTOC(ISYM0)) 2788 KEND2 = KT1Y + NT1AM(ISYMY) 2789 LWRK2 = LWORK - KEND2 2790 2791 IF (LWRK2 .LT. 0) THEN 2792 WRITE(LUPRI,*) 'Memory available : ',LWORK 2793 WRITE(LUPRI,*) 'Memory needed : ',KEND2 2794 CALL QUIT('Insufficient space in CCSDT_GWBIC') 2795 END IF 2796 2797C------------------------ 2798C Construct L(ia,jb). 2799C------------------------ 2800 LENGTH = IRAT*NT2AM(ISYM0) 2801 REWIND(LUIAJB) 2802 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 2803 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYM0,1) 2804 2805C------------------------------------------------------------ 2806C Create Lambda-p^Y and Lambda-h^Y: 2807C------------------------------------------------------------ 2808 IOPT = 1 2809 CALL CC_RDRSP(LISTR,IDLSTR,ISYMY,IOPT,MODEL,WORK(KT1Y),DUMMY) 2810 2811 CALL CCLR_LAMTRA(XLAMDP0,WORK(KLAMPY),XLAMDH0, 2812 * WORK(KLAMHY),WORK(KT1Y),ISYMY) 2813 2814C---------------------------------------------------------------------- 2815C Calculate the F^Y matrix (kc elements evaluated and stored as ck) 2816C---------------------------------------------------------------------- 2817 CALL CC3LR_MFOCK(WORK(KFCKYCK),WORK(KT1Y),WORK(KXIAJB),ISYMY) 2818 2819C----------------------------------------------------------- 2820C Construct 2*C-E of the integrals. 2821C Have integral for both (ij,k,a) and (a,k,j,i) 2822C----------------------------------------------------------- 2823 IOFF = 1 2824 IF (NTOTOC(ISYM0) .GT. 0) THEN 2825cch 2826C write(lupri,*) 'read from fntoc (1)' 2827cch 2828 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYM0)) 2829 ENDIF 2830 2831 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC01),XLAMDH0, 2832 * WORK(KEND2),LWRK2,ISYM0) 2833 2834 CALL CCSDT_TCMEOCC(WORK(KTROC01),WORK(KTROC21),ISYM0) 2835 2836 CALL CCFOP_SORT(WORK(KTROC01),WORK(KTROC03),ISYM0,1) 2837 2838 CALL CCFOP_SORT(WORK(KTROC21),WORK(KTROC23),ISYM0,1) 2839 2840C ------------------------------- 2841C Read in integrals for WMAT etc. 2842C ------------------------------- 2843 IOFF = 1 2844 IF (NTOTOC(ISYMOP) .GT. 0) THEN 2845cch 2846C write(lupri,*) 'read from fntoc (2)' 2847cch 2848 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYMOP)) 2849 ENDIF 2850 2851C ----------------- 2852C LAMPY TRANSFORMED 2853C ----------------- 2854 CALL CCLR_TROCC(WORK(KINTOC),WORK(KTROCCGY), 2855 * WORK(KLAMHY),ISYMY,WORK(KEND2),LWRK2) 2856 CALL CC3_SRTOCC(WORK(KTROCCGY),WORK(KTROCCLY),ISYMY) 2857 2858C ----------------- 2859C LAMP0 TRANSFORMED 2860C ----------------- 2861 CALL CCLR_TROCC(WORK(KINTOC),WORK(KTROCCG0), 2862 * XLAMDH0,ISYM0,WORK(KEND2),LWRK2) 2863 CALL CC3_SRTOCC(WORK(KTROCCG0),WORK(KTROCCL0),ISYM0) 2864C 2865C---------------------------- 2866C Loop over D 2867C---------------------------- 2868C 2869 DO ISYMD = 1,NSYM 2870 2871 ISYCKB = MULD2H(ISYMD,ISYM0) 2872 ISYCKBY = MULD2H(ISYMD,ISYMY) 2873 2874 DO D = 1,NVIR(ISYMD) 2875 2876C ------------------ 2877C Memory allocation. 2878C ------------------ 2879 KTRVI4 = KEND1 2880 KTRVI5 = KTRVI4 + NCKATR(ISYCKB) 2881 KTRVI7 = KTRVI5 + NCKATR(ISYCKB) 2882 KEND3 = KTRVI7 + NCKATR(ISYCKB) 2883 LWRK3 = LWORK - KEND3 2884 2885 KTRVI14 = KEND3 2886 KTRVI15 = KTRVI14 + NCKATR(ISYCKB) 2887 KTRVI18 = KTRVI15 + NCKATR(ISYCKB) 2888 KTRVI19 = KTRVI18 + NCKATR(ISYCKB) 2889 KEND3 = KTRVI19 + NCKATR(ISYCKB) 2890 LWRK3 = LWORK - KEND3 2891 2892 KVGDKCB0 = KEND3 2893 KVGDKBC0 = KVGDKCB0 + NCKATR(ISYCKB) 2894 KVLDKCB0 = KVGDKBC0 + NCKATR(ISYCKB) 2895 KVLDKBC0 = KVLDKCB0 + NCKATR(ISYCKB) 2896 KEND3 = KVLDKBC0 + NCKATR(ISYCKB) 2897 LWRK3 = LWORK - KEND3 2898 2899 KVGDKCBY = KEND3 2900 KVGDKBCY = KVGDKCBY + NCKATR(ISYCKBY) 2901 KVLDKCBY = KVGDKBCY + NCKATR(ISYCKBY) 2902 KVLDKBCY = KVLDKCBY + NCKATR(ISYCKBY) 2903 KEND3 = KVLDKBCY + NCKATR(ISYCKBY) 2904 LWRK3 = LWORK - KEND3 2905 2906 KINTVI = KEND3 2907 KTRVI6 = KINTVI + MAX(NCKA(ISYCKB),NCKA(ISYCKBY)) 2908 KEND4 = KTRVI6 + NCKATR(ISYCKB) 2909 LWRK4 = LWORK - KEND4 2910 2911 IF (LWRK4 .LT. 0) THEN 2912 WRITE(LUPRI,*) 'Memory available : ',LWORK 2913 WRITE(LUPRI,*) 'Memory needed : ',KEND4 2914 CALL QUIT('Insufficient space in CCSDT_GWBIC') 2915 END IF 2916C 2917C ------------------------------------------- 2918C Read 2*C-E of integral used for t3-bar 2919C ------------------------------------------- 2920 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 2921 IF (NCKATR(ISYCKB) .GT. 0) THEN 2922cch 2923C write(lupri,*) 'read from fn3fop2x (3)' 2924cch 2925 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI4),IOFF, 2926 & NCKATR(ISYCKB)) 2927 ENDIF 2928C 2929C ------------------------------------------------- 2930C Integrals used for t3-bar for cc3 2931C ------------------------------------------------- 2932 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 2933 IF (NCKATR(ISYCKB) .GT. 0) THEN 2934cch 2935C write(lupri,*) 'read from fndkbc3 (4)' 2936cch 2937 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI14),IOFF, 2938 & NCKATR(ISYCKB)) 2939 ENDIF 2940 CALL CCSDT_SRVIR3(WORK(KTRVI14),WORK(KEND4), 2941 * ISYMD,D,ISYM0) 2942 CALL CCSDT_SRTVIR(WORK(KTRVI14),WORK(KTRVI15),WORK(KEND4) 2943 * ,LWRK4,ISYMD,ISYM0) 2944C 2945C ------------------------------------------------ 2946C Sort the integrals for t3-bar 2947C ------------------------------------------------ 2948 CALL CCSDT_SRTVIR(WORK(KTRVI4),WORK(KTRVI5),WORK(KEND4), 2949 * LWRK4,ISYMD,ISYM0) 2950C 2951C ---------------------------------------------------- 2952C Read virtual integrals used in q3am/u3am for t3-bar. 2953C ---------------------------------------------------- 2954 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 2955 IF (NCKA(ISYCKB) .GT. 0) THEN 2956cch 2957C write(lupri,*) 'read from fn3vi (4)' 2958cch 2959 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 2960 * NCKA(ISYCKB)) 2961 ENDIF 2962 2963 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI19),XLAMDP0, 2964 * ISYMD,D,ISYM0,WORK(KEND4),LWRK4) 2965 2966 IF (LWRK4 .LT. NCKATR(ISYCKB)) THEN 2967 CALL QUIT('Insufficient space for allocation in '// 2968 * 'CCSDT_GWBIC (CC3 TRVI)') 2969 END IF 2970 2971 CALL CCSDT_SRTVIR(WORK(KTRVI19),WORK(KTRVI18),WORK(KEND4) 2972 * ,LWRK4,ISYMD,ISYM0) 2973 2974 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 2975 IF (NCKATR(ISYCKB) .GT. 0) THEN 2976cch 2977C write(lupri,*) 'gwic> 1: ioff,len:',ioff,nckatr(isyckb) 2978cch 2979 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI6),IOFF, 2980 * NCKATR(ISYCKB)) 2981 ENDIF 2982 2983 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 2984 CALL QUIT('Insufficient space for allocation in '// 2985 & 'CCSDT_GWBIC (2)') 2986 END IF 2987 2988 CALL DCOPY(NCKATR(ISYCKB),WORK(KTRVI6),1,WORK(KTRVI7),1) 2989C 2990C------------------------------------------------------------ 2991C Read in, transform and sort integrals used in 2992C construction of W intermediate. 2993C------------------------------------------------------------ 2994 2995C ------------------------------------------- 2996C Read in g(c1^h k1^p | delta b2^h) integrals 2997C ------------------------------------------- 2998 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 2999 IF (NCKA(ISYCKB) .GT. 0) THEN 3000cch 3001C write(lupri,*) 'read from fn3vi (5)' 3002cch 3003 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 3004 & NCKA(ISYCKB)) 3005 ENDIF 3006 3007C ------------------------------------------------- 3008C LAMPY TRANSFORMED: 3009C transform g(c1^h k1^p | delta b2^h) integrals 3010C to g(c1^h k1^p | d2^pY b2^h) with lampY 3011C and sort as g(d2^pY k1^p | c1^h b2^h) 3012C ------------------------------------------------- 3013 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVGDKCBY), 3014 * WORK(KLAMPY),ISYMY, 3015 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3016 CALL CCSDT_SRVIR3(WORK(KVGDKCBY),WORK(KEND4),ISYMD,D,ISYMY) 3017 3018C ------------------------------------------------- 3019C LAMP0 TRANSFORMED: 3020C transform g(c1^h k1^p | delta b2^h) integrals 3021C to g(c1^h k1^p | d2^p b2^h) with lamp0 3022C and sort as g(d2^p k1^p | c1^h b2^h) 3023C ------------------------------------------------- 3024 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVGDKCB0), 3025 * XLAMDP0,ISYM0, 3026 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3027 CALL CCSDT_SRVIR3(WORK(KVGDKCB0),WORK(KEND4),ISYMD,D,ISYM0) 3028 3029C ------------------------------------------------- 3030C Read in g(b2^h k1^p | delta c1^h) integrals and 3031C transform and sort --- see above 3032C ------------------------------------------------- 3033 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 3034 IF (NCKA(ISYCKB) .GT. 0) THEN 3035cch 3036C write(lupri,*) 'read from fn3vi2 (6)' 3037cch 3038 CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF, 3039 & NCKA(ISYCKB)) 3040 ENDIF 3041 3042C 3043C ------------------ 3044C LAMPY TRANSFORMED: 3045C ------------------ 3046 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVGDKBCY), 3047 * WORK(KLAMPY),ISYMY, 3048 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3049 CALL CCSDT_SRVIR3(WORK(KVGDKBCY),WORK(KEND4),ISYMD,D,ISYMY) 3050 3051C ------------------ 3052C LAMP0 TRANSFORMED: 3053C ------------------ 3054 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVGDKBC0),XLAMDP0,ISYM0, 3055 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3056 CALL CCSDT_SRVIR3(WORK(KVGDKBC0),WORK(KEND4),ISYMD,D,ISYM0) 3057 3058C -------------------------------------------- 3059C Read in L(c1^h k1^p | delta b2^h) integrals: 3060C -------------------------------------------- 3061 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 3062 IF (NCKA(ISYCKB) .GT. 0) THEN 3063cch 3064C write(lupri,*) 'read from fn3fop (7):',lu3fop,fn3fop 3065cch 3066 CALL GETWA2(LU3FOP,FN3FOP,WORK(KINTVI),IOFF,NCKA(ISYCKB)) 3067 ENDIF 3068 3069C --------------------------------------------- 3070C LAMPY TRANSFORMED: 3071C transform L(c1^h k1^p | delta b2^h) integrals 3072C to L(c1^h k1^p | d2^pY b2^h) 3073C and sort as L(d2^pY k1^p | c1^h b2^h) 3074C --------------------------------------------- 3075 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVLDKCBY), 3076 * WORK(KLAMPY),ISYMY, 3077 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3078 CALL CCSDT_SRVIR3(WORK(KVLDKCBY),WORK(KEND4),ISYMD,D,ISYMY) 3079 3080C --------------------------------------------- 3081C LAMP0 TRANSFORMED 3082C transform L(c1^h k1^p | delta b2^h) integrals 3083C to L(c1^h k1^p | d2^p b2^h) 3084C and sort as L(d2^p k1^p | c1^h b2^h) 3085C --------------------------------------------- 3086 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVLDKCB0), 3087 * XLAMDP0,ISYM0, 3088 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3089 CALL CCSDT_SRVIR3(WORK(KVLDKCB0),WORK(KEND4),ISYMD,D,ISYM0) 3090 3091C ----------------------------------------------- 3092C Read in L(b2^h k1^p | delta c1^h) integrals and 3093C transform and sort --- see above 3094C ----------------------------------------------- 3095cch 3096C write(lupri,*) 'read from fn3fop2 (8)' 3097cch 3098 CALL GETWA2(LU3FOP2,FN3FOP2,WORK(KINTVI),IOFF, 3099 & NCKA(ISYCKB)) 3100 3101C ------------------ 3102C LAMPY TRANSFORMED: 3103C ------------------ 3104 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVLDKBCY), 3105 * WORK(KLAMPY),ISYMY, 3106 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3107 CALL CCSDT_SRVIR3(WORK(KVLDKBCY),WORK(KEND4),ISYMD,D,ISYMY) 3108 3109C ------------------ 3110C LAMP0 TRANSFORMED: 3111C ------------------ 3112 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KVLDKBC0), 3113 * XLAMDP0,ISYM0, 3114 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3115 CALL CCSDT_SRVIR3(WORK(KVLDKBC0),WORK(KEND4),ISYMD,D,ISYM0) 3116 3117 3118 DO ISYMB = 1,NSYM 3119 3120 ISYALJ = MULD2H(ISYMB,ISYM0) 3121 ISYALJ2 = MULD2H(ISYMD,ISYM0) 3122 ISYALJBY= MULD2H(ISYMB,ISYMY) 3123 ISYALJDY= MULD2H(ISYMD,ISYMY) 3124 ISYMBD = MULD2H(ISYMD,ISYMB) 3125 ISCKIJ = MULD2H(ISYMBD,ISYM0) 3126 ISWMAT = MULD2H(ISCKIJ,ISYMY) 3127 ISYCKD = MULD2H(ISYM0,ISYMB) 3128 3129C Can use kend3 since we do not need the integrals anymore. 3130 KSMAT2 = KEND3 3131 KUMAT2 = KSMAT2 + NCKIJ(ISCKIJ) 3132 KDIAG = KUMAT2 + NCKIJ(ISCKIJ) 3133 KDIAGW = KDIAG + NCKIJ(ISCKIJ) 3134 KINDSQ = KDIAGW + NCKIJ(ISWMAT) 3135 KINDSQW = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 3136 KINDEX = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 3137 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 3138 KINDEXBY = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 3139 KINDEXDY = KINDEXBY + (NCKI(ISYALJBY) - 1)/IRAT + 1 3140 KTMAT = KINDEXDY + (NCKI(ISYALJDY) - 1)/IRAT + 1 3141 KWMAT = KTMAT + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWMAT)) 3142 KEND4 = KWMAT + NCKIJ(ISWMAT) 3143 LWRK4 = LWORK - KEND4 3144 3145 KTRVI16 = KEND4 3146 KTRVI17 = KTRVI16 + NCKATR(ISYCKD) 3147 KTRVI20 = KTRVI17 + NCKATR(ISYCKD) 3148 KEND4 = KTRVI20 + NCKATR(ISYCKD) 3149 LWRK4 = LWORK - KEND4 3150 3151 KSMAT4 = KEND4 3152 KUMAT4 = KSMAT4 + NCKIJ(ISCKIJ) 3153 KTRVI11 = KUMAT4 + NCKIJ(ISCKIJ) 3154 KTRVI12 = KTRVI11 + NCKATR(ISYCKD) 3155 KTRVI13 = KTRVI12 + NCKATR(ISYCKD) 3156 KEND4 = KTRVI13 + NCKATR(ISYCKD) 3157 LWRK4 = LWORK - KEND4 3158 3159 KINTVI = KEND4 3160 KEND5 = KINTVI + MAX(NCKA(ISYMB),NCKA(ISYCKD)) 3161 LWRK5 = LWORK - KEND5 3162 3163 IF (LWRK5 .LT. 0) THEN 3164 WRITE(LUPRI,*) 'Memory available : ',LWORK 3165 WRITE(LUPRI,*) 'Memory needed : ',KEND5 3166 CALL QUIT('Insufficient space in CCSDT_GWBIC') 3167 END IF 3168 3169C ------------------------------- 3170C Construct part of the diagonal. 3171C ------------------------------- 3172 CALL CC3_DIAG(WORK(KDIAG), WORK(KFOCKD),ISCKIJ) 3173 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 3174 3175C ----------------------- 3176C Construct index arrays. 3177C ----------------------- 3178 LENSQ = NCKIJ(ISCKIJ) 3179 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 3180 LENSQW = NCKIJ(ISWMAT) 3181 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 3182 CALL CC3_INDEX(WORK(KINDEX), ISYALJ) 3183 CALL CC3_INDEX(WORK(KINDEX2), ISYALJ2) 3184 CALL CC3_INDEX(WORK(KINDEXBY),ISYALJBY) 3185 CALL CC3_INDEX(WORK(KINDEXDY),ISYALJDY) 3186 3187 DO B = 1,NVIR(ISYMB) 3188C 3189 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 3190 3191C ------------------------------------ 3192C Read and transform integrals used in 3193C first S-bar and U-bar 3194C ------------------------------------ 3195 IOFF = ICKBD(ISYCKD,ISYMB) 3196 * + NCKATR(ISYCKD)*(B - 1) + 1 3197 IF (NCKATR(ISYCKD) .GT. 0) THEN 3198cch 3199C write(lupri,*) 'read from fndkbc3 (9)' 3200cch 3201 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI16),IOFF, 3202 * NCKATR(ISYCKD)) 3203 ENDIF 3204 CALL CCSDT_SRVIR3(WORK(KTRVI16),WORK(KEND5), 3205 * ISYMB,B,ISYM0) 3206 CALL CCSDT_SRTVIR(WORK(KTRVI16),WORK(KTRVI17), 3207 * WORK(KEND5),LWRK5,ISYMB,ISYM0) 3208C 3209C ------------------------------------ 3210C Read and transform integrals used in 3211C second S-bar and U-bar 3212C ------------------------------------ 3213 IOFF = ICKBD(ISYCKD,ISYMB) 3214 * + NCKATR(ISYCKD)*(B-1) + 1 3215 IF (NCKATR(ISYCKD) .GT. 0) THEN 3216cch 3217C write(lupri,*) 'read from fn3fop2x (10)' 3218cch 3219 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI11), 3220 * IOFF,NCKATR(ISYCKD)) 3221 ENDIF 3222 3223 CALL CCSDT_SRTVIR(WORK(KTRVI11),WORK(KTRVI12), 3224 * WORK(KEND5),LWRK5,ISYMB, 3225 * ISYM0) 3226 3227 IOFF = ICKBD(ISYCKD,ISYMB) 3228 * + NCKATR(ISYCKD)*(B - 1) + 1 3229 IF (NCKATR(ISYCKD) .GT. 0) THEN 3230cch 3231C write(lupri,*) 'gwic> 2: ioff,len:',ioff,nckatr(isyckd) 3232cch 3233 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI13),IOFF, 3234 * NCKATR(ISYCKD)) 3235 ENDIF 3236 3237 IOFF = ICKAD(ISYCKD,ISYMB) 3238 * + NCKA(ISYCKD)*(B - 1) + 1 3239 IF (NCKA(ISYCKD) .GT. 0) THEN 3240cch 3241C write(lupri,*) 'read from fn3vi (10)' 3242cch 3243 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 3244 * NCKA(ISYCKD)) 3245 ENDIF 3246 3247 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI20), 3248 * XLAMDP0,ISYMB,B,ISYM0, 3249 * WORK(KEND5),LWRK5) 3250 3251C ---------------------------------------------------- 3252C Calculate the S(ci,bk,dj) matrix for B,D for T3-BAR: 3253C ---------------------------------------------------- 3254 CALL DZERO(WORK(KSMAT2),NCKIJ(ISCKIJ)) 3255 3256 CALL CCFOP_SMAT(0.0D0,WORK(KL1AM),ISYM0,WORK(KL2TP), 3257 * ISYM0,WORK(KTMAT), 3258 * WORK(KFCKBA),WORK(KXIAJB),ISYM0, 3259 * WORK(KTRVI14),WORK(KTRVI15), 3260 * WORK(KTRVI4),WORK(KTRVI5), 3261 * WORK(KTROC01),WORK(KTROC21), 3262 * ISYM0,WORK(KFOCKD),WORK(KDIAG), 3263 * WORK(KSMAT2),WORK(KEND5),LWRK5, 3264 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 3265 * ISYMB,B,ISYMD,D) 3266 3267 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT2),1) 3268 3269 CALL T3_FORBIDDEN(WORK(KSMAT2),ISYM0,ISYMB,B,ISYMD,D) 3270C 3271C ---------------------------------------------------- 3272C Calculate the S(ci,bk,dj) matrix for D,B for T3-BAR: 3273C ---------------------------------------------------- 3274 CALL DZERO(WORK(KSMAT4),NCKIJ(ISCKIJ)) 3275 3276 CALL CCFOP_SMAT(0.0D0,WORK(KL1AM),ISYM0,WORK(KL2TP), 3277 * ISYM0,WORK(KTMAT),WORK(KFCKBA), 3278 * WORK(KXIAJB),ISYM0, 3279 * WORK(KTRVI16),WORK(KTRVI17), 3280 * WORK(KTRVI11),WORK(KTRVI12), 3281 * WORK(KTROC01),WORK(KTROC21), 3282 * ISYM0,WORK(KFOCKD),WORK(KDIAG), 3283 * WORK(KSMAT4),WORK(KEND5),LWRK5, 3284 * WORK(KINDEX2),WORK(KINDSQ), 3285 * LENSQ,ISYMD,D,ISYMB,B) 3286 3287 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT4),1) 3288 CALL T3_FORBIDDEN(WORK(KSMAT4),ISYM0,ISYMD,D,ISYMB,B) 3289C 3290C ------------------------------------------------ 3291C Calculate U(ci,jk) for fixed b,d for t3-bar. 3292C ------------------------------------------------ 3293 CALL DZERO(WORK(KUMAT2),NCKIJ(ISCKIJ)) 3294 3295 CALL CCFOP_UMAT(0.0D0,WORK(KL1AM),ISYM0,WORK(KL2TP), 3296 * ISYM0, 3297 * WORK(KXIAJB),ISYM0,WORK(KFCKBA), 3298 * WORK(KTRVI19),WORK(KTRVI7), 3299 * WORK(KTROC03),WORK(KTROC23),ISYM0, 3300 * WORK(KFOCKD),WORK(KDIAG), 3301 * WORK(KUMAT2), 3302 * WORK(KTMAT),WORK(KEND5),LWRK5, 3303 * WORK(KINDSQ),LENSQ,ISYMB,B,ISYMD,D) 3304 3305 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KUMAT2),1) 3306 CALL T3_FORBIDDEN(WORK(KUMAT2),ISYM0,ISYMB,B,ISYMD,D) 3307C 3308C ------------------------------------------------ 3309C Calculate U(ci,jk) for fixed d,b for t3-bar. 3310C ------------------------------------------------ 3311 CALL DZERO(WORK(KUMAT4),NCKIJ(ISCKIJ)) 3312 3313 CALL CCFOP_UMAT(0.0D0,WORK(KL1AM),ISYM0,WORK(KL2TP), 3314 * ISYM0,WORK(KXIAJB),ISYM0, 3315 * WORK(KFCKBA),WORK(KTRVI20), 3316 * WORK(KTRVI13),WORK(KTROC03), 3317 * WORK(KTROC23),ISYM0, 3318 * WORK(KFOCKD),WORK(KDIAG), 3319 * WORK(KUMAT4),WORK(KTMAT), 3320 * WORK(KEND5),LWRK5,WORK(KINDSQ), 3321 * LENSQ,ISYMD,D,ISYMB,B) 3322 3323 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KUMAT4),1) 3324 CALL T3_FORBIDDEN(WORK(KUMAT4),ISYM0,ISYMD,D,ISYMB,B) 3325C 3326C ------------------------------------------- 3327C Sum up the S-bar and U-bar to get a real T3-bar 3328C ------------------------------------------- 3329 CALL CC3_T3BD(ISCKIJ,WORK(KSMAT2),WORK(KSMAT4), 3330 * WORK(KUMAT2),WORK(KUMAT4), 3331 * WORK(KTMAT),WORK(KINDSQ),LENSQ) 3332C 3333C ---------------------------- 3334C Calculate <L3|[Y^,tau3]|HF>: 3335C ---------------------------- 3336 CALL WBARBD_V(WORK(KTMAT),ISCKIJ,WORK(KFOCKY),ISYMY, 3337 * WORK(KWMAT),ISWMAT,WORK(KEND5),LWRK5) 3338 3339 IF (SPLIT_WT) THEN 3340C --------------------------------------------------- 3341C 1)Store this contribution of WMAT in THETA(ai,bj,dl): 3342C but first divide by orbital energies and clean up 3343C 2)Store zeroth-order T3 amplitudes in T0AMP: 3344C --------------------------------------------------- 3345 CALL WBD_DIA(B,ISYMB,D,ISYMD,-FREQ, 3346 * ISWMAT,WORK(KWMAT), 3347 * WORK(KDIAGW),WORK(KFOCKD)) 3348 CALL T3_FORBIDDEN(WORK(KWMAT),ISYMY,ISYMB,B,ISYMD,D) 3349 3350 DO ISYML = 1, NSYM 3351 ISYMDL = MULD2H(ISYMD,ISYML) 3352 ISAIBJ = MULD2H(ISYMY,ISYMDL) 3353 DO L = 1, NRHF(ISYML) 3354 DO ISYMJ = 1, NSYM 3355 ISYMBJ = MULD2H(ISYMJ,ISYMB) 3356 ISYMAI = MULD2H(ISAIBJ,ISYMBJ) 3357 ISYAIL = MULD2H(ISYMAI,ISYML) 3358 DO J = 1, NRHF(ISYMJ) 3359 3360 KOFF1 = ISAIKJ(ISYAIL,ISYMJ)+ NCKI(ISYAIL)*(J-1) 3361 * + ISAIK(ISYMAI, ISYML)+NT1AM(ISYMAI)*(L-1) 3362 3363 CALL DCOPY(NT1AMX,WORK(KWMAT+KOFF1),1, 3364 * THETA(1,1,B,J,D,L),1) 3365 END DO 3366 END DO 3367 END DO 3368 END DO 3369 ! Reinitialize WMAT to avoid a double count: 3370 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 3371 END IF 3372 3373 CALL WBARBD_O(WORK(KTMAT),ISCKIJ,WORK(KFOCKY),ISYMY, 3374 * WORK(KWMAT),ISWMAT,WORK(KEND5),LWRK5) 3375 3376C --------------------------- 3377C Calculate <L2|[Y,tau3]|HF>: 3378C --------------------------- 3379 CALL WBARBD_T2(B,ISYMB,D,ISYMD,WORK(KL2TP),ISYM0, 3380 * WORK(KFOCKY),ISYMY,WORK(KWMAT),ISWMAT) 3381 3382C ----------------------------- 3383C calculate <L2|[H^Y,tau3]|HF>: 3384C ----------------------------- 3385 CALL WBARBD_TMAT(WORK(KL2TP),ISYM0,WORK(KWMAT), 3386 * WORK(KTMAT),ISWMAT,WORK(KFCKYCK), 3387 * ISYMY,WORK(KVLDKBCY), 3388 * WORK(KVLDKCBY), 3389 * WORK(KVGDKBCY),WORK(KVGDKCBY), 3390 * WORK(KTROCCLY),WORK(KTROCCGY), 3391 * ISYMY,WORK(KEND5),LWRK5, 3392 * WORK(KINDEX),WORK(KINDEX2), 3393 * WORK(KINDSQW),LENSQW,ISYMB,B,ISYMD,D) 3394 3395C ----------------------------- 3396C Calculate <L2Y|[H^,tau3]|HF>: 3397C ----------------------------- 3398 CALL WBARBD_TMAT(WORK(KL2TPBY),ISYMY,WORK(KWMAT), 3399 * WORK(KTMAT),ISWMAT,WORK(KFCKBA), 3400 * ISYM0,WORK(KVLDKBC0), 3401 * WORK(KVLDKCB0), 3402 * WORK(KVGDKBC0),WORK(KVGDKCB0), 3403 * WORK(KTROCCL0),WORK(KTROCCG0), 3404 * ISYM0,WORK(KEND5),LWRK5, 3405 * WORK(KINDEXBY),WORK(KINDEXDY), 3406 * WORK(KINDSQW),LENSQW,ISYMB,B,ISYMD,D) 3407 3408C ----------------------------- 3409C Calculate <L1Y|[H^,tau3]|HF>: 3410C ----------------------------- 3411 CALL WBARBD_L1(WORK(KL1BY),ISYMY,WORK(KTMAT), 3412 * WORK(KXIAJB),ISYM0, 3413 * WORK(KWMAT),WORK(KEND5),LWRK5, 3414 * WORK(KINDSQW),LENSQW,ISYMB,B,ISYMD,D) 3415 3416C ----------------------------------- 3417C Divide by the energy difference and 3418C remove the forbidden elements 3419C ----------------------------------- 3420 CALL WBD_DIA(B,ISYMB,D,ISYMD,-FREQ, 3421 * ISWMAT,WORK(KWMAT), 3422 * WORK(KDIAGW),WORK(KFOCKD)) 3423C 3424 CALL T3_FORBIDDEN(WORK(KWMAT),ISYMY,ISYMB,B,ISYMD,D) 3425C 3426C ----------------------------- 3427C Store WMAT as WBAR(ai,bj,dl): 3428C ----------------------------- 3429 DO ISYML = 1, NSYM 3430 ISYMDL = MULD2H(ISYMD,ISYML) 3431 ISAIBJ = MULD2H(ISYMY,ISYMDL) 3432 DO L = 1, NRHF(ISYML) 3433 DO ISYMJ = 1, NSYM 3434 ISYMBJ = MULD2H(ISYMJ,ISYMB) 3435 ISYMAI = MULD2H(ISAIBJ,ISYMBJ) 3436 ISYAIL = MULD2H(ISYMAI,ISYML) 3437 DO J = 1, NRHF(ISYMJ) 3438 3439 KOFF1 = ISAIKJ(ISYAIL,ISYMJ)+ NCKI(ISYAIL)*(J-1) 3440 * + ISAIK(ISYMAI, ISYML)+NT1AM(ISYMAI)*(L-1) 3441 3442 CALL DCOPY(NT1AMX,WORK(KWMAT+KOFF1),1, 3443 * WBAR(1,1,B,J,D,L),1) 3444 END DO 3445 END DO 3446 END DO 3447 END DO 3448 3449 ENDDO ! B 3450 ENDDO ! ISYMB 3451 ENDDO ! D 3452 ENDDO ! ISYMD 3453 3454*---------------------------------------------------------------------* 3455* Close files and return: 3456*---------------------------------------------------------------------* 3457 CALL WCLOSE2(LU3VI2,FN3VI2,'KEEP') 3458 CALL WCLOSE2(LU3FOP,FN3FOP,'KEEP') 3459 CALL WCLOSE2(LU3FOP2,FN3FOP2,'KEEP') 3460 3461 CALL QEXIT('CCSDT_GWBIC') 3462 3463 RETURN 3464 END 3465*=====================================================================* 3466*=====================================================================* 3467 SUBROUTINE CCSDT_GWTIC(LISTR,IDLSTR,WINT,THETA,T0AMP, 3468 & HAVET3Y,TAMP3Y,PRINT_T3, 3469 & XLAMDP0,XLAMDH0,FOCK0, 3470 & LUDELD,FNDELD,LUDKBC,FNDKBC,LUCKJD,FNCKJD, 3471 & WORK,LWORK) 3472*---------------------------------------------------------------------* 3473* 3474* Purpose: collect w^Y and theta^y intermediates in memory 3475* for use in noddy code 3476* 3477* Written by Christof Haettig, Jan 2003, Friedrichstal 3478* 3479*=====================================================================* 3480 IMPLICIT NONE 3481C 3482#include "priunit.h" 3483#include "dummy.h" 3484#include "iratdef.h" 3485#include "ccsdsym.h" 3486#include "ccorb.h" 3487#include "ccsdinp.h" 3488#include "ccr1rsp.h" 3489#include "ccinftap.h" 3490#include "inftap.h" 3491C 3492 LOGICAL LOCDBG, TOT_T3Y, HAVET3Y, PRINT_T3 3493 PARAMETER(LOCDBG = .FALSE., TOT_T3Y = .TRUE.) 3494C 3495 INTEGER ISYM0 3496 PARAMETER(ISYM0 = 1) 3497C 3498 CHARACTER*11 FN3SRTR, FNCKJDR, FNDELDR, FNDKBCR, FNDKBCR4 3499 PARAMETER(FN3SRTR = 'CCSDT_ETA_1', FNCKJDR = 'CCSDT_ETA_2', 3500 * FNDELDR = 'CCSDT_ETA_3', FNDKBCR = 'CCSDT_ETA_4', 3501 * FNDKBCR4 = 'CCSDT_ETA_5') 3502 3503 CHARACTER*(*) FNDKBC, FNDELD, FNCKJD 3504 3505 CHARACTER LISTR*3 3506 3507 CHARACTER MODEL*10, LABELY*8, CDUMMY*1 3508 3509 INTEGER LWORK, IDLSTR, ISTAMY 3510 3511 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR, LUDKBCR4, 3512 & LUCKJD, LUDKBC, LUDELD, LUFOCK, 3513 & KEND1, KFOCKY, KFOCKD, KFCKBA, KEND0, KT2TP, KT1B, KT2B, 3514 & LWRK0, KTROC0, KTROC02,LWRK1,KINTOC, KEND2, LWRK2, 3515 & KTRVI0, KTRVI1, KTRVI2, KEND3, LWRK3, KTRVI0Y, KINTVI, 3516 & KSMAT, KSMAT3, KUMAT, KUMAT3, KDIAG, KDIAGW, KINDSQW, 3517 & KINDSQ, KINDEX, KINDEX2, KTMAT, KTRVI8, KTRVI9, KTRVI10, 3518 & KEND4, LWRK4, KTRVI8Y, KWMAT, KTROC0Y, 3519 & IOPT, IRREP, IMAT, IERR, LENSQ, LENSQW, IOFF, 3520 & ISYMC, ISYMK, KOFF1, KOFF2, ISYMD, ISYCKB, ISCKBY, 3521 & ISYMB, ISYALJ, ISYALJ2, ISYMBD, ISCKIJ, ISYCKD, ISCKDY, 3522 & ISWMAT, ISYML, ISYMDL, ISAIBJ, ISYMJ, ISYMBJ, ISYMAI, 3523 & ISYAIL, KOFFA, KEND5, LWRK5 3524 3525#if defined (SYS_CRAY) 3526 REAL XLAMDP0(*), XLAMDH0(*), FOCK0(*) 3527 REAL WORK(LWORK), FREQY, DDOT 3528 REAL ZERO, ONE, TWO, HALF 3529 REAL WINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3530 REAL THETA(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3531 REAL T0AMP(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3532 REAL TAMP3Y(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3533#else 3534 DOUBLE PRECISION XLAMDP0(*), XLAMDH0(*), FOCK0(*) 3535 DOUBLE PRECISION WORK(LWORK), FREQY, DDOT 3536 DOUBLE PRECISION ZERO, ONE, TWO, HALF 3537 DOUBLE PRECISION WINT(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3538 DOUBLE PRECISION THETA(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3539 DOUBLE PRECISION T0AMP(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3540 DOUBLE PRECISION TAMP3Y(NVIRT,NRHFT,NVIRT,NRHFT,NVIRT,NRHFT) 3541#endif 3542C 3543 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 3544C 3545C------------------------------------------------------------ 3546C some initializations: 3547C------------------------------------------------------------ 3548 CALL QENTER('CCSDT_GWTIC') 3549 3550 IF (LISTR(1:3).EQ.'R1 ') THEN 3551 ISTAMY = ISYLRT(IDLSTR) 3552 FREQY = FRQLRT(IDLSTR) 3553 LABELY = LRTLBL(IDLSTR) 3554c 3555 IF (LORXLRT(IDLSTR)) 3556 & CALL QUIT('NO ORBITAL RELAXED PERTURBATION IN CCSDT_GWTIC') 3557 ELSE 3558 ! ups, probably higher-order response, not yet implemented here 3559 CALL QUIT('Unknown type of right vector in CCSDT_GWTIC.') 3560 END IF 3561C 3562C--------------------- 3563C Open some files: 3564C--------------------- 3565 IF (TOT_T3Y) THEN 3566 LU3SRTR = -1 3567 LUCKJDR = -1 3568 LUDELDR = -1 3569 LUDKBCR = -1 3570 LUDKBCR4 = -1 3571 3572 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 3573 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 3574 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 3575 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 3576 CALL WOPEN2(LUDKBCR4,FNDKBCR4,64,0) 3577 ENDIF 3578C 3579C----------------------------------------------- 3580C initial allocations and preparation of T2: 3581C----------------------------------------------- 3582 KEND0 = 1 3583 3584 IF (LISTR(1:3).EQ.'R1 ') THEN 3585 KFOCKY = KEND0 3586 KEND0 = KFOCKY + N2BST(ISTAMY) 3587 END IF 3588 3589 KFOCKD = KEND0 3590 KFCKBA = KFOCKD + NORBTS 3591 KEND0 = KFCKBA + NT1AMX 3592 3593 KT2TP = KEND0 3594 KEND0 = KT2TP + NT2SQ(ISYM0) 3595 3596 IF (TOT_T3Y) THEN 3597 KT1B = KEND0 3598 KT2B = KT1B + NT1AM(ISTAMY) 3599 KEND0 = KT2B + NT2SQ(ISTAMY) 3600 END IF 3601 3602 LWRK0 = LWORK - KEND0 3603 3604 IF (LWRK0.LT.NT2SQ(ISYM0) .OR. LWRK0.LT.NT2SQ(ISTAMY)) THEN 3605 CALL QUIT('Not enough memory to construct T2TP (CCSDT_GWTIC)') 3606 ENDIF 3607 3608 3609 IOPT = 2 3610 CALL CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KT2TP)) 3611 3612 CALL CC_T2SQ(WORK(KT2TP),WORK(KEND0),ISYM0) 3613 3614 CALL CC3_T2TP(WORK(KT2TP),WORK(KEND0),ISYM0) 3615 3616 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of T2TP ', 3617 * DDOT(NT2SQ(ISYM0),WORK(KT2TP),1,WORK(KT2TP),1) 3618 3619 IF (TOT_T3Y) THEN 3620 IOPT = 3 3621 CALL CC_RDRSP(LISTR,IDLSTR,ISTAMY,IOPT,MODEL, 3622 * WORK(KT1B),WORK(KT2B)) 3623 CALL CCLR_DIASCL(WORK(KT2B),TWO,ISTAMY) 3624 CALL CC_T2SQ(WORK(KT2B),WORK(KEND0),ISTAMY) 3625 CALL CC3_T2TP(WORK(KT2B),WORK(KEND0),ISTAMY) 3626 END IF 3627 3628C---------------------------------------------------- 3629C Calculate (ck|de)-tilde(Y) and (ck|lm)-tilde(Y) 3630C---------------------------------------------------- 3631 IF (TOT_T3Y) THEN 3632 CALL CC3_BARINT(WORK(KT1B),ISTAMY,XLAMDP0, 3633 * XLAMDH0,WORK(KEND0),LWRK0, 3634 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 3635C 3636 CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISTAMY,LU3SRTR,FN3SRTR, 3637 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 3638 * IDUMMY,CDUMMY) 3639C 3640 CALL CC3_SINT(XLAMDH0,WORK(KEND0),LWRK0,ISTAMY, 3641 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 3642C 3643 CALL CC3_TCME(XLAMDP0,ISTAMY,WORK(KEND0),LWRK0, 3644 * IDUMMY,CDUMMY,LUDKBCR,FNDKBCR, 3645 * IDUMMY,CDUMMY,IDUMMY,CDUMMY, 3646 * IDUMMY,CDUMMY,LUDKBCR4,FNDKBCR4,2) 3647 ENDIF 3648C 3649C------------------------------------- 3650C Read canonical orbital energies: 3651C------------------------------------- 3652C 3653 LUSIFC = -1 3654 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 3655 & .FALSE.) 3656 REWIND LUSIFC 3657 3658 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 3659 READ (LUSIFC) 3660 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 3661 3662 CALL GPCLOSE(LUSIFC,'KEEP') 3663 3664 IF (FROIMP .OR. FROEXP) 3665 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND0),LWRK0) 3666C 3667C----------------------------------------------------- 3668C Sort the fock matrix 3669C----------------------------------------------------- 3670C 3671 DO ISYMC = 1,NSYM 3672 ISYMK = MULD2H(ISYMC,ISYM0) 3673 DO K = 1,NRHF(ISYMK) 3674 DO C = 1,NVIR(ISYMC) 3675 KOFF1 = IFCVIR(ISYMK,ISYMC) + 3676 * NORB(ISYMK)*(C - 1) + K 3677 KOFF2 = IT1AM(ISYMC,ISYMK) 3678 * + NVIR(ISYMC)*(K - 1) + C 3679 WORK(KFCKBA-1+KOFF2) = FOCK0(KOFF1) 3680 ENDDO 3681 ENDDO 3682 ENDDO 3683 3684 IF (LOCDBG) THEN 3685 CALL AROUND('In CCSDT_GWTIC: Triples Fock MO matrix (sort)') 3686 CALL CC_PRFCKMO(WORK(KFCKBA),ISYM0) 3687 ENDIF 3688C 3689C----------------------------------------------------------- 3690C Prepare one-electron operators needed to compute the 3691C amplitude response vectors: 3692C----------------------------------------------------------- 3693C 3694 IF (LISTR(1:3).EQ.'R1 ') THEN 3695 ! read property integrals from file: 3696 CALL CCPRPAO(LABELY,.TRUE.,WORK(KFOCKY),IRREP,IMAT,IERR, 3697 * WORK(KEND0),LWRK0) 3698 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISTAMY)) THEN 3699 WRITE(LUPRI,*) 'ISTAMY:',ISTAMY 3700 WRITE(LUPRI,*) 'IRREP :',IRREP 3701 WRITE(LUPRI,*) 'IERR :',IERR 3702 CALL QUIT('CCSDT_GWTIC: error reading operator '//LABELY) 3703 ELSE IF (IERR.LT.0) THEN 3704 CALL DZERO(WORK(KFOCKY),N2BST(ISTAMY)) 3705 END IF 3706 3707 ! transform property integrals to Lambda-MO basis 3708 CALL CC_FCKMO(WORK(KFOCKY),XLAMDP0,XLAMDH0, 3709 & WORK(KEND0),LWRK0,ISTAMY,1,1) 3710 END IF 3711 3712C 3713C----------------------------- 3714C Memory allocation. 3715C----------------------------- 3716C 3717 KTROC0 = KEND0 3718 KTROC02 = KTROC0 + NTRAOC(ISYM0) 3719 KEND1 = KTROC02 + NTRAOC(ISYM0) 3720 IF (TOT_T3Y) THEN 3721 KTROC0Y = KEND1 3722 KEND1 = KTROC0Y + NTRAOC(ISTAMY) 3723 ENDIF 3724 LWRK1 = LWORK - KEND1 3725 3726 KINTOC = KEND1 3727 KEND2 = KINTOC + MAX(NTOTOC(ISYM0),NTOTOC(ISYM0)) 3728 LWRK2 = LWORK - KEND2 3729 3730 IF (LWRK2 .LT. 0) THEN 3731 WRITE(LUPRI,*) 'Memory available : ',LWORK 3732 WRITE(LUPRI,*) 'Memory needed : ',KEND2 3733 CALL QUIT('Insufficient space in CCSDT_GWTIC') 3734 END IF 3735C 3736C------------------------ 3737C Occupied integrals. 3738C------------------------ 3739C 3740 IOFF = 1 3741 IF (NTOTOC(ISYM0) .GT. 0) THEN 3742 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISYM0)) 3743 ENDIF 3744 3745 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of OCC-INT ', 3746 * DDOT(NTOTOC(ISYM0),WORK(KINTOC),1,WORK(KINTOC),1) 3747C 3748C---------------------------------------------------------------------- 3749C Transform (ia|j delta) integrals to (ia|j k) and sort as (ij,k,a) 3750C---------------------------------------------------------------------- 3751C 3752 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0),XLAMDP0, 3753 * WORK(KEND2),LWRK2,ISYM0) 3754 3755 CALL CCFOP_SORT(WORK(KTROC0),WORK(KTROC02),ISYM0,1) 3756 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of TROC0 ', 3757 * DDOT(NTRAOC(ISYM0),WORK(KTROC0),1,WORK(KTROC0),1) 3758C 3759C------------------------------------------ 3760C Y transformed Occupied integrals. 3761C----------------------------------------- 3762C 3763 IF (TOT_T3Y) THEN 3764 IOFF = 1 3765 IF (NTOTOC(ISTAMY) .GT. 0) THEN 3766 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF, 3767 * NTOTOC(ISTAMY)) 3768 ENDIF 3769C 3770 IF (LOCDBG) 3771 * WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (Y transformed) ', 3772 * DDOT(NTOTOC(ISTAMY),WORK(KINTOC),1,WORK(KINTOC),1) 3773C 3774 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0Y),XLAMDP0, 3775 * WORK(KEND2),LWRK2,ISTAMY) 3776 ENDIF 3777 3778C 3779C---------------------------- 3780C Loop over D 3781C---------------------------- 3782C 3783 DO ISYMD = 1,NSYM 3784 3785 ISYCKB = MULD2H(ISYMD,ISYM0) 3786 ISCKBY = MULD2H(ISTAMY,ISYMD) 3787 IF (LOCDBG) WRITE(LUPRI,*) 'In CCSDT_ETA_CON: ISYCKB :',ISYCKB 3788 3789 DO D = 1,NVIR(ISYMD) 3790C 3791C ------------------ 3792C Memory allocation. 3793C ------------------ 3794 KTRVI0 = KEND1 3795 KTRVI1 = KTRVI0 + NCKATR(ISYCKB) 3796 KTRVI2 = KTRVI1 + NCKATR(ISYCKB) 3797 KEND3 = KTRVI2 + NCKATR(ISYCKB) 3798 LWRK3 = LWORK - KEND3 3799 3800 IF (TOT_T3Y) THEN 3801 KTRVI0Y = KEND3 3802 KEND3 = KTRVI0Y + NCKATR(ISCKBY) 3803 LWRK3 = LWORK - KEND3 3804 ENDIF 3805 3806 KINTVI = KEND3 3807 KEND4 = KINTVI + MAX(NCKA(ISYMD),NCKA(ISYCKB)) 3808 LWRK4 = LWORK - KEND4 3809 3810 IF (LWRK4 .LT. 0) THEN 3811 WRITE(LUPRI,*) 'Memory available : ',LWORK 3812 WRITE(LUPRI,*) 'Memory needed : ',KEND4 3813 CALL QUIT('Insufficient space in CCSDT_GWTIC') 3814 END IF 3815C 3816C ------------------------------------ 3817C Integrals used in s3am. 3818C ------------------------------------ 3819 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 3820 IF (NCKATR(ISYCKB) .GT. 0) THEN 3821 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI0),IOFF, 3822 & NCKATR(ISYCKB)) 3823 ENDIF 3824C 3825C ------------------------------------------------------------ 3826C Read B transformed virtual integrals used for W for TOT_T3Y 3827C ------------------------------------------------------------ 3828 IF (TOT_T3Y) THEN 3829 IOFF = ICKBD(ISCKBY,ISYMD) + NCKATR(ISCKBY)*(D - 1) + 1 3830 IF (NCKATR(ISCKBY) .GT. 0) THEN 3831 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI0Y),IOFF, 3832 & NCKATR(ISCKBY)) 3833 ENDIF 3834 ENDIF 3835 3836C ------------------------------------------------ 3837C Sort the integrals for s3am: 3838C ------------------------------------------------ 3839 CALL CCSDT_SRTVIR(WORK(KTRVI0),WORK(KTRVI2),WORK(KEND4), 3840 * LWRK4,ISYMD,ISYM0) 3841 3842C ------------------------------------------- 3843C Read virtual integrals used in contraction. 3844C ------------------------------------------- 3845 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 3846 IF (NCKA(ISYCKB) .GT. 0) THEN 3847 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 3848 * NCKA(ISYCKB)) 3849 ENDIF 3850 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI1),XLAMDH0, 3851 * ISYMD,D,ISYM0,WORK(KEND4),LWRK4) 3852 3853 3854 DO ISYMB = 1,NSYM 3855 3856 ISYALJ = MULD2H(ISYMB,ISYM0) 3857 ISYALJ2 = MULD2H(ISYMD,ISYM0) 3858 ISYMBD = MULD2H(ISYMB,ISYMD) 3859 ISCKIJ = MULD2H(ISYMBD,ISYM0) 3860 ISYCKD = MULD2H(ISYM0,ISYMB) 3861 ISCKDY = MULD2H(ISTAMY,ISYMB) 3862 ISWMAT = MULD2H(ISCKIJ,ISTAMY) 3863 3864 IF (LOCDBG) THEN 3865 WRITE(LUPRI,*) 'In CCSDT_GWTIC: ISYMD :',ISYMD 3866 WRITE(LUPRI,*) 'In CCSDT_GWTIC: ISYMB :',ISYMB 3867 WRITE(LUPRI,*) 'In CCSDT_GWTIC: ISYALJ:',ISYALJ 3868 WRITE(LUPRI,*) 'In CCSDT_GWTIC: ISYMBD:',ISYMBD 3869 WRITE(LUPRI,*) 'In CCSDT_GWTIC: ISCKIJ:',ISCKIJ 3870 ENDIF 3871C 3872C Can use kend3 since we do not need the integrals anymore. 3873 KSMAT = KEND3 3874 KSMAT3 = KSMAT + NCKIJ(ISCKIJ) 3875 KUMAT = KSMAT3 + NCKIJ(ISCKIJ) 3876 KUMAT3 = KUMAT + NCKIJ(ISCKIJ) 3877 KDIAG = KUMAT3 + NCKIJ(ISCKIJ) 3878 KDIAGW = KDIAG + NCKIJ(ISCKIJ) 3879 KINDSQW = KDIAGW + NCKIJ(ISWMAT) 3880 KINDSQ = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 3881 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 3882 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 3883 KTMAT = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 3884 KTRVI8 = KTMAT + NCKIJ(ISCKIJ) 3885 KTRVI9 = KTRVI8 + NCKATR(ISYCKD) 3886 KTRVI10 = KTRVI9 + NCKATR(ISYCKD) 3887 KEND4 = KTRVI10 + NCKATR(ISYCKD) 3888 LWRK4 = LWORK - KEND4 3889 3890 IF (TOT_T3Y) THEN 3891 KTRVI8Y = KEND4 3892 KEND4 = KTRVI8Y + NCKATR(ISCKDY) 3893 ENDIF 3894 3895 KWMAT = KEND4 3896 KEND4 = KWMAT + NCKIJ(ISWMAT) 3897 LWRK4 = LWORK - KEND4 3898 3899 KINTVI = KEND4 3900 KEND5 = KINTVI + MAX(NCKA(ISYMB),NCKA(ISYCKD)) 3901 LWRK5 = LWORK - KEND5 3902 3903 IF (LWRK5 .LT. 0) THEN 3904 WRITE(LUPRI,*) 'Memory available : ',LWORK 3905 WRITE(LUPRI,*) 'Memory needed : ',KEND5 3906 CALL QUIT('Insufficient space in CCSDT_GWTIC') 3907 END IF 3908C 3909C ------------------------------- 3910C Construct part of the diagonal. 3911C ------------------------------- 3912 CALL CC3_DIAG(WORK(KDIAG), WORK(KFOCKD),ISCKIJ) 3913 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 3914 3915 IF (LOCDBG) THEN 3916 WRITE(LUPRI,*) 'Norm of DIA ', 3917 * DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1,WORK(KDIAG),1) 3918 WRITE(LUPRI,*) 'Norm of DIA_W', 3919 * DDOT(NCKIJ(ISWMAT),WORK(KDIAGW),1,WORK(KDIAGW),1) 3920 END IF 3921C 3922C ----------------------- 3923C Construct index arrays. 3924C ----------------------- 3925 LENSQ = NCKIJ(ISCKIJ) 3926 LENSQW = NCKIJ(ISWMAT) 3927 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 3928 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 3929 CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2) 3930 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 3931 3932 DO B = 1,NVIR(ISYMB) 3933 IF (LOCDBG) write(lupri,*) 'For B, D : ',B,D 3934 3935C --------------------------------------------- 3936C Initialize WMAT: 3937C --------------------------------------------- 3938 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 3939 3940C --------------------------------------------- 3941C Read and transform integrals used in second S 3942C --------------------------------------------- 3943 IOFF = ICKBD(ISYCKD,ISYMB) 3944 * + NCKATR(ISYCKD)*(B - 1) + 1 3945 IF (NCKATR(ISYCKD) .GT. 0) THEN 3946 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI8),IOFF, 3947 * NCKATR(ISYCKD)) 3948 ENDIF 3949 3950 CALL CCSDT_SRTVIR(WORK(KTRVI8),WORK(KTRVI9), 3951 * WORK(KEND4),LWRK4,ISYMB,ISYM0) 3952C 3953C ---------------------------------------------- 3954C Read B transformed virtual integrals used in W 3955C ---------------------------------------------- 3956 IF (TOT_T3Y) THEN 3957 IOFF = ICKBD(ISCKDY,ISYMB) + 3958 & NCKATR(ISCKDY)*(B - 1) + 1 3959 IF (NCKATR(ISCKDY) .GT. 0) THEN 3960 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI8Y),IOFF, 3961 & NCKATR(ISCKDY)) 3962 ENDIF 3963 ENDIF 3964C 3965C ----------------------------------------- 3966C Read virtual integrals used in second U 3967C ----------------------------------------- 3968 IOFF = ICKAD(ISYCKD,ISYMB) 3969 * + NCKA(ISYCKD)*(B - 1) + 1 3970 IF (NCKA(ISYCKD) .GT. 0) THEN 3971 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 3972 * NCKA(ISYCKD)) 3973 ENDIF 3974 3975 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI10), 3976 * XLAMDH0,ISYMB,B,ISYM0, 3977 * WORK(KEND5),LWRK5) 3978C 3979C -------------------------------------------------- 3980C Calculate the S(ci,bk,dj) matrix for T3 for B,D. 3981C -------------------------------------------------- 3982 CALL CC3_SMAT(0.0D0,WORK(KT2TP),ISYM0,WORK(KTMAT), 3983 * WORK(KTRVI0),WORK(KTRVI2), 3984 * WORK(KTROC0),ISYM0, 3985 * WORK(KFOCKD),WORK(KDIAG), 3986 * WORK(KSMAT),WORK(KEND4),LWRK4, 3987 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 3988 * ISYMB,B,ISYMD,D) 3989 3990 CALL T3_FORBIDDEN(WORK(KSMAT),ISYM0,ISYMB,B,ISYMD,D) 3991C 3992 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of SMAT :', 3993 * DDOT(NCKIJ(ISCKIJ),WORK(KSMAT),1,WORK(KSMAT),1) 3994C 3995C -------------------------------------------------- 3996C Calculate the S(ci,bk,dj) matrix for T3 for D,B. 3997C -------------------------------------------------- 3998 CALL CC3_SMAT(0.0D0,WORK(KT2TP),ISYM0,WORK(KTMAT), 3999 * WORK(KTRVI8),WORK(KTRVI9), 4000 * WORK(KTROC0),ISYM0, 4001 * WORK(KFOCKD),WORK(KDIAG), 4002 * WORK(KSMAT3),WORK(KEND4),LWRK4, 4003 * WORK(KINDEX2),WORK(KINDSQ),LENSQ, 4004 * ISYMD,D,ISYMB,B) 4005 4006 CALL T3_FORBIDDEN(WORK(KSMAT3),ISYM0,ISYMD,D,ISYMB,B) 4007C 4008 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of SMAT3:', 4009 * DDOT(NCKIJ(ISCKIJ),WORK(KSMAT3),1,WORK(KSMAT3),1) 4010C 4011C --------------------------------- 4012C Calculate U(ci,jk) for fixed b,d. 4013C --------------------------------- 4014 CALL CC3_UMAT(0.0D0,WORK(KT2TP),ISYM0,WORK(KTRVI1), 4015 * WORK(KTROC02),ISYM0,WORK(KFOCKD), 4016 * WORK(KDIAG),WORK(KUMAT),WORK(KTMAT), 4017 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 4018 * ISYMB,B,ISYMD,D) 4019 4020 CALL T3_FORBIDDEN(WORK(KUMAT),ISYM0,ISYMB,B,ISYMD,D) 4021C 4022 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of UMAT :', 4023 * DDOT(NCKIJ(ISCKIJ),WORK(KUMAT),1,WORK(KUMAT),1) 4024C 4025C --------------------------------- 4026C Calculate U(ci,jk) for fixed d,b. 4027C --------------------------------- 4028 CALL CC3_UMAT(0.0D0,WORK(KT2TP),ISYM0,WORK(KTRVI10), 4029 * WORK(KTROC02),ISYM0,WORK(KFOCKD), 4030 * WORK(KDIAG),WORK(KUMAT3),WORK(KTMAT), 4031 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 4032 * ISYMD,D,ISYMB,B) 4033 4034 CALL T3_FORBIDDEN(WORK(KUMAT3),ISYM0,ISYMD,D,ISYMB,B) 4035 4036 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of UMAT3:', 4037 * DDOT(NCKIJ(ISCKIJ),WORK(KUMAT3),1,WORK(KUMAT3),1) 4038C 4039C ----------------------------------- 4040C Sum up the S and U to get a real T3 4041C ----------------------------------- 4042 CALL CC3_T3BD(ISCKIJ,WORK(KSMAT),WORK(KSMAT3), 4043 * WORK(KUMAT),WORK(KUMAT3), 4044 * WORK(KTMAT),WORK(KINDSQ),LENSQ) 4045C 4046C --------------------------------- 4047C Use the T3 in TMAT to calculate W 4048C --------------------------------- 4049 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of TMAT:', 4050 * DDOT(NCKIJ(ISCKIJ),WORK(KTMAT),1,WORK(KTMAT),1) 4051 4052 CALL WBD_V(WORK(KTMAT),ISCKIJ,WORK(KFOCKY),ISTAMY, 4053 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 4054 4055 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_V):', 4056 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 4057 4058 IF (.TRUE.) THEN 4059C --------------------------------------------------- 4060C 1)Store this contribution of WMAT in THETA(ai,bj,dl): 4061C but first divide by orbital energies and clean up 4062C 2)Store zeroth-order T3 amplitudes in T0AMP: 4063C --------------------------------------------------- 4064 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQY, 4065 * ISWMAT,WORK(KWMAT), 4066 * WORK(KDIAGW),WORK(KFOCKD)) 4067 CALL T3_FORBIDDEN(WORK(KWMAT),ISTAMY,ISYMB,B,ISYMD,D) 4068 4069 DO ISYML = 1, NSYM 4070 ISYMDL = MULD2H(ISYMD,ISYML) 4071 ISAIBJ = MULD2H(ISTAMY,ISYMDL) 4072 DO L = 1, NRHF(ISYML) 4073 DO ISYMJ = 1, NSYM 4074 ISYMBJ = MULD2H(ISYMJ,ISYMB) 4075 ISYMAI = MULD2H(ISAIBJ,ISYMBJ) 4076 ISYAIL = MULD2H(ISYMAI,ISYML) 4077 DO J = 1, NRHF(ISYMJ) 4078 4079 KOFFA = ISAIKJ(ISYAIL,ISYMJ)+ NCKI(ISYAIL)*(J-1) 4080 * + ISAIK(ISYMAI, ISYML)+NT1AM(ISYMAI)*(L-1) 4081 4082 CALL DCOPY(NT1AMX,WORK(KWMAT+KOFFA),1, 4083 & THETA(1,1,B,J,D,L),1) 4084 4085 CALL DCOPY(NT1AMX,WORK(KTMAT+KOFFA),1, 4086 & T0AMP(1,1,B,J,D,L),1) 4087 4088 END DO 4089 END DO 4090 END DO 4091 END DO 4092 ! Reinitialize WMAT to avoid a double count: 4093 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 4094 END IF 4095 4096 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_V)', 4097 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 4098 4099 CALL WBD_O(WORK(KTMAT),ISCKIJ,WORK(KFOCKY),ISTAMY, 4100 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 4101 4102 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_O)', 4103 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 4104C 4105 CALL WBD_T2(.FALSE.,B,ISYMB,D,ISYMD, 4106 * WORK(KT2TP),ISYM0,WORK(KT2TP),ISYM0, 4107 * WORK(KFOCKY),ISTAMY,WORK(KINDSQW),LENSQW, 4108 * WORK(KWMAT),ISWMAT,WORK(KEND4),LWRK4) 4109C 4110 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_T2):', 4111 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 4112 4113 IF (TOT_T3Y) THEN 4114 CALL WBD_GROUND(WORK(KT2B),ISTAMY,WORK(KTMAT), 4115 * WORK(KTRVI0),WORK(KTRVI8), 4116 * WORK(KTROC0),1,WORK(KWMAT), 4117 * WORK(KEND4),LWRK4, 4118 * WORK(KINDSQW),LENSQW, 4119 * ISYMB,B,ISYMD,D) 4120c 4121 CALL WBD_GROUND(WORK(KT2TP),ISYM0,WORK(KTMAT), 4122 * WORK(KTRVI0Y),WORK(KTRVI8Y), 4123 * WORK(KTROC0Y),ISTAMY,WORK(KWMAT), 4124 * WORK(KEND4),LWRK4, 4125 * WORK(KINDSQW),LENSQW, 4126 * ISYMB,B,ISYMD,D) 4127 ENDIF 4128 4129 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQY, 4130 * ISWMAT,WORK(KWMAT), 4131 * WORK(KDIAGW),WORK(KFOCKD)) 4132 4133 CALL T3_FORBIDDEN(WORK(KWMAT),ISTAMY,ISYMB,B,ISYMD,D) 4134 4135 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of WMAT (WBD_DIA)', 4136 * DDOT(NCKIJ(ISWMAT),WORK(KWMAT),1,WORK(KWMAT),1) 4137 4138C ----------------------------- 4139C Store WMAT in WINT(ai,bj,dl): 4140C ----------------------------- 4141 DO ISYML = 1, NSYM 4142 ISYMDL = MULD2H(ISYMD,ISYML) 4143 ISAIBJ = MULD2H(ISTAMY,ISYMDL) 4144 DO L = 1, NRHF(ISYML) 4145 DO ISYMJ = 1, NSYM 4146 ISYMBJ = MULD2H(ISYMJ,ISYMB) 4147 ISYMAI = MULD2H(ISAIBJ,ISYMBJ) 4148 ISYAIL = MULD2H(ISYMAI,ISYML) 4149 DO J = 1, NRHF(ISYMJ) 4150 4151 KOFFA = ISAIKJ(ISYAIL,ISYMJ)+ NCKI(ISYAIL)*(J-1) 4152 * + ISAIK(ISYMAI, ISYML)+NT1AM(ISYMAI)*(L-1) 4153 4154 CALL DCOPY(NT1AMX,WORK(KWMAT+KOFFA),1, 4155 & WINT(1,1,B,J,D,L),1) 4156 4157 END DO 4158 END DO 4159 END DO 4160 END DO 4161 4162 ENDDO ! B 4163 ENDDO ! ISYMB 4164 4165 ENDDO ! D 4166 ENDDO ! ISYMD 4167c 4168 4169C------------- 4170C End 4171C------------- 4172 IF (TOT_T3Y) THEN 4173 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 4174 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 4175 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 4176 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 4177 CALL WCLOSE2(LUDKBCR4,FNDKBCR4,'DELETE') 4178 ENDIF 4179 4180 IF (PRINT_T3) THEN 4181 WRITE(LUPRI,*) 'CCSDT_GWTIC> W^Y:' 4182 WRITE(LUPRI,*) 'CCSDT_GWTIC> LISTR,IDLSTR:',LISTR,IDLSTR 4183 CALL PRINT_PT3_NODDY(WINT) 4184 WRITE(LUPRI,*) 'CCSDT_GWTIC> Theta^Y:' 4185 WRITE(LUPRI,*) 'CCSDT_GWTIC> LISTR,IDLSTR:',LISTR,IDLSTR 4186 CALL PRINT_PT3_NODDY(THETA) 4187 END IF 4188 4189 IF (HAVET3Y) THEN 4190 ! accumulate T^Y amplitudes in TAMP3Y and print 4191 DO K = 1, NRHFT 4192 DO C = 1, NVIRT 4193 DO J = 1, NRHFT 4194 DO B = 1, NVIRT 4195 DO I = 1, NRHFT 4196 DO A = 1, NVIRT 4197 TAMP3Y(A,I,B,J,C,K) = 4198 & WINT(A,I,B,J,C,K) + THETA(A,I,B,J,C,K) + 4199 & WINT(C,K,A,I,B,J) + THETA(C,K,A,I,B,J) + 4200 & WINT(B,J,C,K,A,I) + THETA(B,J,C,K,A,I) 4201 END DO 4202 END DO 4203 END DO 4204 END DO 4205 END DO 4206 END DO 4207 IF (PRINT_T3) THEN 4208 WRITE(LUPRI,*) 'CCSDT_GWTIC> T3^Y:' 4209 WRITE(LUPRI,*) 'CCSDT_GWTIC> LISTR,IDLSTR:',LISTR,IDLSTR 4210 CALL PRINT_PT3_NODDY(TAMP3Y) 4211 END IF 4212 END IF 4213 4214 CALL QEXIT('CCSDT_GWTIC') 4215 4216 RETURN 4217 END 4218*=====================================================================* 4219