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 Deck /* cc3_abmatt3zu */ 19 SUBROUTINE CC3_ABMATT3ZU(LISTX,IDLSTX,LISTZU,IDLSTZU, 20 & OMEGA2, 21 & OMEGA2EFF, 22 * ISYRES,WORK,LWORK) 23C 24 IMPLICIT NONE 25#include "priunit.h" 26C 27 CHARACTER*3 LISTX,LISTZU 28 INTEGER IDLSTX,IDLSTZU 29C 30 INTEGER ISYRES, LWORK 31C 32#if defined (SYS_CRAY) 33 REAL OMEGA2(*),OMEGA2EFF(*) 34 REAL WORK(LWORK) 35#else 36 DOUBLE PRECISION OMEGA2(*),OMEGA2EFF(*) 37 DOUBLE PRECISION WORK(LWORK) 38#endif 39C 40 CALL QENTER('ABT3ZU') 41C 42 CALL CC3_AMATT3ZU(LISTX,IDLSTX,LISTZU,IDLSTZU, 43 & OMEGA2,OMEGA2EFF, 44 * ISYRES,WORK,LWORK) 45C 46 CALL CC3_BMATT3ZU(LISTX,IDLSTX,LISTZU,IDLSTZU, 47 & OMEGA2,OMEGA2EFF, 48 * ISYRES,WORK,LWORK) 49C 50C---------- 51C End . 52C---------- 53C 54 CALL QEXIT('ABT3ZU') 55C 56 RETURN 57 END 58C /* Deck cc3_amatt3zu */ 59 SUBROUTINE CC3_AMATT3ZU(LISTX,IDLSTX,LISTZU,IDLSTZU, 60 & OMEGA2, 61 & OMEGA2EFF, 62 * ISYRES,WORK,LWORK) 63C 64 IMPLICIT NONE 65#include "priunit.h" 66#include "ccsdsym.h" 67#include "ccorb.h" 68#include "ccsdinp.h" 69C 70 INTEGER LUCKJD, LUDELD, LU3VI, LUTOC,LU3VI2, LUDKBC, LU3FOP 71C 72 CHARACTER*6 FNCKJD, FNDELD, FN3VI 73 CHARACTER*8 FNTOC,FN3VI2 74 CHARACTER*4 FNDKBC 75 CHARACTER*5 FN3FOP 76C 77 PARAMETER (FNCKJD='CKJDEL', FNTOC = 'CCSDT_OC' ) 78 PARAMETER (FNDELD = 'CKDELD' , FNDKBC = 'DKBC' ) 79 PARAMETER (FN3VI = 'CC3_VI' , FN3VI2 = 'CC3_VI12') 80 PARAMETER (FN3FOP='PTFOP' ) 81C 82 CHARACTER*3 LISTX,LISTZU, 83 INTEGER IDLSTX,IDLSTZU 84C 85 INTEGER ISYRES, LWORK 86 integer isym 87C 88#if defined (SYS_CRAY) 89 REAL OMEGA2(*),OMEGA2EFF(*) 90 REAL WORK(LWORK) 91#else 92 DOUBLE PRECISION OMEGA2(*),OMEGA2EFF(*) 93 DOUBLE PRECISION WORK(LWORK) 94#endif 95C 96 CALL QENTER('AT3ZU') 97C 98C-------------------------------- 99C Open files 100C-------------------------------- 101C 102 LUCKJD = -1 103 LUTOC = -1 104 LUDELD = -1 105 LUDKBC = -1 106 LU3VI2 = -1 107 LU3VI = -1 108 LU3FOP = -1 109C 110 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 111 CALL WOPEN2(LUTOC,FNTOC,64,0) 112 CALL WOPEN2(LUDELD,FNDELD,64,0) 113 CALL WOPEN2(LUDKBC,FNDKBC,64,0) 114 CALL WOPEN2(LU3VI2,FN3VI2,64,0) 115 CALL WOPEN2(LU3VI,FN3VI,64,0) 116 CALL WOPEN2(LU3FOP,FN3FOP,64,0) 117C 118 CALL CC3_AMATT3ZUVIR(LISTX,IDLSTX,LISTZU,IDLSTZU, 119 & OMEGA2, 120 & OMEGA2EFF, 121 * ISYRES,WORK,LWORK, 122 * LUCKJD,FNCKJD,LUTOC,FNTOC, 123 * LUDELD,FNDELD,LUDKBC,FNDKBC, 124 * LU3VI2,FN3VI2,LU3VI,FN3VI) 125C 126c 127 CALL CC3_AMATT3ZUOCC(LISTX,IDLSTX,LISTZU,IDLSTZU, 128 & OMEGA2, 129 & OMEGA2EFF, 130 * ISYRES,WORK,LWORK, 131 * LUCKJD,FNCKJD,LUTOC,FNTOC, 132 * LUDELD,FNDELD,LU3FOP,FN3FOP, 133 * LU3VI,FN3VI) 134c 135 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 136 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 137 CALL WCLOSE2(LUDELD,FNDELD,'KEEP') 138 CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP') 139 CALL WCLOSE2(LU3VI2,FN3VI2,'KEEP') 140 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 141 CALL WCLOSE2(LU3FOP,FN3FOP,'KEEP') 142C 143C---------- 144C End . 145C---------- 146C 147 CALL QEXIT('AT3ZU') 148C 149 RETURN 150 END 151C /* Deck cc3_amatt3zuvir */ 152 SUBROUTINE CC3_AMATT3ZUVIR(LISTX,IDLSTX,LISTZU,IDLSTZU, 153 & OMEGA2, 154 & OMEGA2EFF, 155 * ISYRES,WORK,LWORK, 156 * LUCKJD,FNCKJD,LUTOC,FNTOC, 157 * LUDELD,FNDELD,LUDKBC,FNDKBC, 158 * LU3VI2,FN3VI2,LU3VI,FN3VI) 159C 160* 161******************************************************************* 162* 163* This routine calculates the contribution to cubic response 164* (originating from both F matrix and B matrix): 165* 166* omega2eff = <mu2|[H^X,T3^ZU]|HF> 167* 168* which is contracted outside with the multipliers (zeroth-order for 169* F matrix or first-order for B matrix). 170* 171* 172* 173* T3^ZU is calculated using W^BD intermmediate. 174* 175* !!! ACTUALLY HERE WE CALCULATE ONLY THE PART OF THE TERM, WHICH ORIGINATES 176* FROM THE AY-MATRIX CONTRIBUTIONS TO T3^ZU 177* (to get the rest of contribution call cc3_bmatt3zu routine) !!! 178* 179* Thus here: 180* 181* W^BD = <mu3|[[Z,T1^U],T3^0]|HF> + <mu3|[[Z,T2^U],T2^0]|HF> 182* + <mu3|[Z,T3^U]|HF> 183* 184* (!!! TO GET THE TOTAL CONTRIBUTION FROM THE LAST TERM CALL CC3_AMATT3ZUOCC !!!) 185* 186* F. Pawlowski, P. Jorgensen, Aarhus Spring 2003. 187* 188************************************************************** 189 190C 191 IMPLICIT NONE 192#include "priunit.h" 193#include "ccsdsym.h" 194#include "dummy.h" 195#include "iratdef.h" 196#include "inftap.h" 197#include "ccinftap.h" 198#include "ccorb.h" 199#include "ccsdinp.h" 200#include "ccr1rsp.h" 201#include "ccexci.h" 202C 203 CHARACTER*(*) FNCKJD,FNTOC,FNDELD,FNDKBC,FN3VI2,FN3VI 204 INTEGER LUCKJD,LUTOC,LUDELD,LUDKBC,LU3VI2,LU3VI 205C 206 CHARACTER*12 FN3SRTR, FNCKJDRZ, FNDELDRZ, FNDKBCRZ 207 CHARACTER*13 FNCKJDRU,FNDELDRU,FNDKBCRU 208 INTEGER LU3SRTR, LUCKJDRZ, LUDELDRZ, LUDKBCRZ 209 INTEGER LUCKJDRU,LUDELDRU,LUDKBCRU 210C 211 CHARACTER CDUMMY*1 212C 213 PARAMETER(FN3SRTR = 'CCSDT_FBMAT1', FNCKJDRZ = 'CCSDT_FBMAT2', 214 * FNDELDRZ = 'CCSDT_FBMAT3', FNDKBCRZ = 'CCSDT_FBMAT4') 215C 216 PARAMETER(FNCKJDRU = 'CCSDT_FBMAT22',FNDELDRU = 'CCSDT_FBMAT33', 217 * FNDKBCRU = 'CCSDT_FBMAT44') 218C 219 CHARACTER*3 LISTRZ, LISTRU, LISTX, LISTZU 220 CHARACTER*8 LABELRZ,LABELRU 221 INTEGER IDLSTRZ,IDLSTRU,IDLSTX,IDLSTZU 222 INTEGER ISYMRZ, ISYMRU 223C 224 LOGICAL T2XNET2Z 225 LOGICAL LORXRZ,LORXRU 226C 227 INTEGER ISYRES, LWORK 228 INTEGER ISYM0,ISINT1,ISINT2,ISYMT2,ISYMT3 229 INTEGER KT2TP,KFOCKD,KLAMDP,KLAMDH,KEND00,LWRK00,KXIAJB 230 INTEGER IOPTTCME,IOPT,IOFF 231 INTEGER ISYMWZ,ISYMWZU,ISYFCKZ,ISYFCKU,ISYMZU 232 INTEGER KT2Z,KT2U,KFCKZ,KFCKU,KEND0,LWRK0 233 INTEGER KFCKZUV,KFCKZUO,KFCKUZV,KFCKUZO,KT1Z,KT1U,KEND1,LWRK1 234 INTEGER KLAMDPU,KLAMDHU,KLAMDPZ,KLAMDHZ 235 INTEGER ISINTRZ1,ISINTRZ2,ISINTRU1,ISINTRU2 236 INTEGER KRBJIAZU,KT3OG1,KT3OG2,KW3ZOGZ1,KW3UOGU1,KTROC,KTROC1 237 INTEGER KINTOC 238 INTEGER ISYMD,ISCKB1,ISCKB2,ISCKB2Z,ISCKB2U 239 INTEGER KT3VDG3,KT3VDG2,KEND2,LWRK2 240 INTEGER KT3VDG1,KEND3,LWRK3,KW3ZVDGZ1,KW3UVDGU1,KT3VDGZ3,KT3VDGU3 241 INTEGER KTRVI,KTRVI1,MAXZU 242 INTEGER ISYMB,ISYALJ,ISYALJ2,ISYMBD,ISCKIJ,ISCKD2,ISCKD2Z,ISCKD2U 243 INTEGER ISWMATZ,ISWMATU,ISWMATZU 244 INTEGER KSMAT,KSMAT3,KUMAT,KUMAT3,KDIAG,KT3MAT 245 INTEGER KDIAGWZ,KWMATZ,KINDSQWZ,KDIAGWU,KWMATU,KINDSQWU 246 INTEGER KINDSQ,KINDEX,KINDEX2 247 INTEGER MAXZZUU,KTMATW,KT3VBG1,KT3VBG2,KT3VBG3,KEND4,LWRK4 248 INTEGER KW3ZVDGZ2,KW3UVDGU2,KT3VBGZ3,KT3VBGU3,KINTVI 249 INTEGER KINDSQWZU,KDIAGWZU,KWMATZU 250 INTEGER LENSQ,LENSQWZ,LENSQWU,LENSQWZU 251 INTEGER AIBJCK_PERM 252c 253 INTEGER ISYMRX,ISYFCKX 254 INTEGER KLAMDPX,KLAMDHX,KFOCKX,KT1X,ISYMC,ISYMK,KOFF1,KOFF2 255 INTEGER ISINT1X,ISCKB1X 256c 257 integer kx3am 258c 259#if defined (SYS_CRAY) 260 REAL OMEGA2(*),OMEGA2EFF(*) 261 REAL WORK(LWORK) 262 REAL FREQRZ,FREQRU,FREQZU 263 REAL XNORMVAL,DDOT 264#else 265 DOUBLE PRECISION OMEGA2(*),OMEGA2EFF(*) 266 DOUBLE PRECISION WORK(LWORK) 267 DOUBLE PRECISION FREQRZ,FREQRU,FREQZU 268 DOUBLE PRECISION XNORMVAL,DDOT 269#endif 270C 271 CALL QENTER('AT3ZUV') 272c 273C 274C-------------------------- 275C Save and set flags. 276C-------------------------- 277C 278C Set symmetry flags. 279C 280C 281C ISYMT2 is symmetry of T2TP 282C ISINT2 is symmetry of integrals in triples equation (ISINT2=1) 283C ISINT1 is symmetry of integrals in contraction (ISINT1=1) 284C ISYMT3 is symmetry of S and U intermediate 285C ISYMW is symmetry of W intermediate 286C ISYRES is symmetry of result vector 287C 288C------------------------------------------------------------- 289C 290 IPRCC = IPRINT 291 ISYM0 = 1 292 ISINT1 = 1 293 ISINT2 = 1 294 ISYMT2 = 1 295 ISYMT3 = MULD2H(ISINT2,ISYMT2) 296C 297C-------------------------------- 298C Open temporary files 299C-------------------------------- 300C 301 LU3SRTR = -1 302 LUCKJDRZ = -1 303 LUCKJDRU = -1 304 LUDELDRZ = -1 305 LUDELDRU = -1 306 LUDKBCRZ = -1 307 LUDKBCRU = -1 308C 309 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 310 CALL WOPEN2(LUCKJDRZ,FNCKJDRZ,64,0) 311 CALL WOPEN2(LUCKJDRU,FNCKJDRU,64,0) 312 CALL WOPEN2(LUDELDRZ,FNDELDRZ,64,0) 313 CALL WOPEN2(LUDELDRU,FNDELDRU,64,0) 314 CALL WOPEN2(LUDKBCRZ,FNDKBCRZ,64,0) 315 CALL WOPEN2(LUDKBCRU,FNDKBCRU,64,0) 316C 317C---------------------------------------------- 318C Calculate the zeroth order stuff once 319C---------------------------------------------- 320C 321 KT2TP = 1 322 KFOCKD = KT2TP + NT2SQ(ISYMT2) 323 KLAMDP = KFOCKD + NORBTS 324 KLAMDH = KLAMDP + NLAMDT 325 KEND00 = KLAMDH + NLAMDT 326 LWRK00 = LWORK - KEND00 327C 328 KXIAJB = KEND00 329 KEND00 = KXIAJB + NT2AM(ISINT1) 330 LWRK00 = LWORK - KEND00 331C 332 IF (LWRK00 .LT. 0) THEN 333 CALL QUIT('Out of memory in CC3_AMATT3ZUVIR (zeroth allo.') 334 ENDIF 335C 336C------------------------ 337C Construct L(ia,jb). 338C------------------------ 339C 340 REWIND(LUIAJB) 341 CALL READI(LUIAJB,IRAT*NT2AM(ISINT1),WORK(KXIAJB)) 342 IOPTTCME = 1 343 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISINT1,IOPTTCME) 344C 345 IF ( IPRINT .GT. 55) THEN 346 XNORMVAL = DDOT(NT2AM(ISINT1),WORK(KXIAJB),1, 347 * WORK(KXIAJB),1) 348 WRITE(LUPRI,*) 'Norm of IAJB ',XNORMVAL 349 ENDIF 350C 351C---------------------------------------------------------------- 352C Read t1 and calculate the zero'th order Lambda matrices 353C---------------------------------------------------------------- 354C 355 CALL GET_LAMBDA0(WORK(KLAMDP),WORK(KLAMDH),WORK(KEND00),LWRK00) 356C 357C------------------------------------------- 358C Read in t2 , square it and reorder 359C------------------------------------------- 360C 361 IOPT = 2 362 CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KT2TP),'R0',0,ISYMT2, 363 * WORK(KEND00),LWRK00) 364 365 IF (IPRINT .GT. 55) THEN 366 XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 367 WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL 368 ENDIF 369C 370C-------------------------------------- 371C Read in orbital energies 372C-------------------------------------- 373C 374C 375 CALL GET_ORBEN(WORK(KFOCKD),WORK(KEND00),LWRK00) 376COMMENT COMMENT COMMENT 377COMMENT COMMENT COMMENT If we want to sum the T3 amplitudes 378COMMENT COMMENT COMMENT 379 if (.false.) then 380 kx3am = kend00 381 kend00 = kx3am + nt1amx*nt1amx*nt1amx 382 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 383 lwrk00 = lwork - kend00 384 if (lwrk00 .lt. 0) then 385 write(lupri,*) 'Memory available : ',lwork 386 write(lupri,*) 'Memory needed : ',kend00 387 call quit('Insufficient space (T3) in CC3_AMATT3ZUVIR') 388 END IF 389 endif 390COMMENT COMMENT COMMENT 391COMMENT COMMENT COMMENT 392COMMENT COMMENT COMMENT 393C 394C determining R1 labels 395C 396 ISYMRX = ISYLRT(IDLSTX) 397 !FREQuency not needed for LISTX 398 !LABELX = LRTLBL(IDLSTX) not needed for LISTX 399C 400 IF (LISTZU(1:3).EQ.'R2 ') THEN 401 LISTRZ = 'R1 ' 402 LABELRZ = LBLR2T(IDLSTZU,1) 403 ISYMRZ = ISYR2T(IDLSTZU,1) 404 FREQRZ = FRQR2T(IDLSTZU,1) 405 LORXRZ = LORXR2T(IDLSTZU,1) 406 IDLSTRZ = IR1TAMP(LABELRZ,LORXRZ,FREQRZ,ISYMRZ) 407C 408 IF (LORXRZ) CALL QUIT('NO ORBITAL RELAX. IN CC3_AMATT3ZUVIR') 409 410 LISTRU = 'R1 ' 411 LABELRU = LBLR2T(IDLSTZU,2) 412 ISYMRU = ISYR2T(IDLSTZU,2) 413 FREQRU = FRQR2T(IDLSTZU,2) 414 LORXRU = LORXR2T(IDLSTZU,2) 415 IDLSTRU = IR1TAMP(LABELRU,LORXRU,FREQRU,ISYMRU) 416C 417 IF (LORXRU) CALL QUIT('NO ORBITAL RELAX. IN CC3_AMATT3ZUVIR') 418C 419 ELSE 420 CALL QUIT('Unknown list for LISTZU in CC3_AMATT3ZUVIR.') 421 END IF 422C 423 FREQZU = FREQRZ + FREQRU 424C 425 ISYMWZ = MULD2H(ISYMT3,ISYMRZ) 426C 427 ISYMWZU = MULD2H(ISYMWZ,ISYMRU) 428C 429 IF (ISYRES.NE.ISYMWZU) THEN 430 WRITE(LUPRI,*) 'INCONSISTENT SYMMETRY SPECIFICATION' 431 WRITE(LUPRI,*) 432 * 'ISYRES =',ISYRES,'ISYMWZU =',ISYMWZU 433 CALL QUIT('CC3_AMATT3ZUVIR inconsistent symmetry specification') 434 END IF 435C 436 ISYFCKZ = ISYMRZ 437 ISYFCKU = ISYMRU 438 ISYMZU = MULD2H(ISYMRZ,ISYMRU) 439C 440C------------------------------------------------- 441C Read T1^Z and T2^Z 442C Calculate (ck|de)-tilde(Z) and (ck|lm)-tilde(Z) 443C Used to construct T3^Z 444C 445C Read T1^U and T2^U 446C Calculate (ck|de)-tilde(U) and (ck|lm)-tilde(U) 447C Used to construct T3^U 448C------------------------------------------------- 449C 450 ISYFCKX = MULD2H(ISYMOP,ISYMRX) 451C 452 KT2Z = KEND00 453 KT2U = KT2Z + NT2SQ(ISYMRZ) 454 KFCKZ = KT2U + NT2SQ(ISYMRU) 455 KFCKU = KFCKZ + N2BST(ISYFCKZ) 456 KEND0 = KFCKU + N2BST(ISYFCKU) 457 LWRK0 = LWORK - KEND0 458C 459 KFCKZUV = KEND0 460 KFCKZUO = KFCKZUV + N2BST(ISYMZU) 461 KFCKUZV = KFCKZUO + N2BST(ISYMZU) 462 KFCKUZO = KFCKUZV + N2BST(ISYMZU) 463 KEND0 = KFCKUZO + N2BST(ISYMZU) 464C 465 KLAMDPX = KEND0 466 KLAMDHX = KLAMDPX + NLAMDT 467 KEND0 = KLAMDHX + NLAMDT 468 LWRK0 = LWORK - KEND0 469C 470 KFOCKX = KEND0 471 KEND0 = KFOCKX + N2BST(ISYFCKX) 472 LWRK0 = LWORK - KEND0 473C 474 KT1Z = KEND0 475 KT1U = KT1Z + NT1AM(ISYMRZ) 476 KEND1 = KT1U + NT1AM(ISYMRU) 477C 478 KT1X = KEND1 479 KEND1 = KT1X + NT1AM(ISYMRX) 480 LWRK1 = LWORK - KEND1 481C 482 KLAMDPU = KEND1 483 KLAMDHU = KLAMDPU + NLAMDT 484 KLAMDPZ = KLAMDHU + NLAMDT 485 KLAMDHZ = KLAMDPZ + NLAMDT 486 KEND1 = KLAMDHZ + NLAMDT 487 LWRK1 = LWORK - KEND1 488C 489 IF (LWRK1 .LT. NT2SQ(ISYMRZ)) THEN 490 CALL QUIT('Out of memory in CC3_AMATT3ZUVIR (TOT_T3Y) ') 491 ENDIF 492C 493C--------------------------- 494C Read Lambda_X matrices 495C--------------------------- 496C 497 CALL GET_LAMBDAX(WORK(KLAMDPX),WORK(KLAMDHX),LISTX,IDLSTX, 498 * ISYMRX,WORK(KLAMDP),WORK(KLAMDH),WORK(KEND1), 499 * LWRK1) 500C 501C T1^Z and T2^Z amplitudes 502C 503 IOPT = 3 504 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1Z),WORK(KT2Z),LISTRZ,IDLSTRZ, 505 * ISYMRZ,WORK(KEND1),LWRK1) 506C 507 IF (IPRINT .GT. 55) THEN 508 XNORMVAL = DDOT(NT2SQ(ISYMRZ),WORK(KT2Z),1,WORK(KT2Z),1) 509 WRITE(LUPRI,*) 'Norm of T2X ',XNORMVAL 510 XNORMVAL = DDOT(NT1AM(ISYMRZ),WORK(KT1Z),1,WORK(KT1Z),1) 511 WRITE(LUPRI,*) 'Norm of T1X ',XNORMVAL 512 ENDIF 513C 514C T1^U and T2^U amplitudes 515C 516 IOPT = 3 517 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1U),WORK(KT2U),LISTRU,IDLSTRU, 518 * ISYMRU,WORK(KEND1),LWRK1) 519C 520 IF (IPRINT .GT. 55) THEN 521 XNORMVAL = DDOT(NT2SQ(ISYMRU),WORK(KT2U),1,WORK(KT2U),1) 522 WRITE(LUPRI,*) 'Norm of T2Y ',XNORMVAL 523 XNORMVAL = DDOT(NT1AM(ISYMRU),WORK(KT1U),1,WORK(KT1U),1) 524 WRITE(LUPRI,*) 'Norm of T1Y ',XNORMVAL 525 ENDIF 526C 527C------------------ 528C Read in T1^X 529C------------------ 530C 531 IOPT = 1 532 CALL GET_T1_T2(IOPT,.FALSE.,WORK(KT1X),DUMMY,LISTX,IDLSTX, 533 * ISYMRX,WORK(KEND1),LWRK1) 534C 535 IF (IPRINT .GT. 55) THEN 536 XNORMVAL = DDOT(NT1AM(ISYMRX),WORK(KT1X),1,WORK(KT1X),1) 537 WRITE(LUPRI,*) 'Norm of T1A ',XNORMVAL 538 ENDIF 539C 540C------------------------------------------ 541C Calculate the F^X matrix 542C------------------------------------------ 543C 544 CALL CC3LR_MFOCK(WORK(KEND1),WORK(KT1X),WORK(KXIAJB),ISYFCKX) 545C 546C Put the F_{kc} part into correct F_{pq} 547C 548 CALL DZERO(WORK(KFOCKX),N2BST(ISYFCKX)) 549C 550 DO ISYMC = 1, NSYM 551 ISYMK = MULD2H(ISYFCKX,ISYMC) 552 DO K = 1, NRHF(ISYMK) 553 DO C = 1, NVIR(ISYMC) 554 KOFF1 = KFOCKX - 1 555 * + IFCVIR(ISYMK,ISYMC) 556 * + NORB(ISYMK)*(C - 1) 557 * + K 558 KOFF2 = KEND1 - 1 559 * + IT1AM(ISYMC,ISYMK) 560 * + NVIR(ISYMC)*(K - 1) 561 * + C 562C 563 WORK(KOFF1) = WORK(KOFF2) 564C 565 ENDDO 566 ENDDO 567 ENDDO 568C 569 IF (IPRINT .GT. 55) THEN 570 CALL AROUND( 'In CC3_BMATT3ZU: Fock^A MO matrix' ) 571 CALL CC_PRFCKMO(WORK(KFOCKX),ISYFCKX) 572 ENDIF 573C 574 575C 576C Integrals H^Z in T3^Z 577C 578 579 ISINTRZ1 = MULD2H(ISINT1,ISYMRZ) 580 ISINTRZ2 = MULD2H(ISINT2,ISYMRZ) 581C 582 CALL CC3_BARINT(WORK(KT1Z),ISYMRZ,WORK(KLAMDP), 583 * WORK(KLAMDH),WORK(KEND1),LWRK1, 584 * LU3SRTR,FN3SRTR,LUCKJDRZ,FNCKJDRZ) 585C 586 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTRZ1,LU3SRTR,FN3SRTR, 587 * LUDELDRZ,FNDELDRZ,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 588 * IDUMMY,CDUMMY) 589C 590 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTRZ1, 591 * LUDELDRZ,FNDELDRZ,LUDKBCRZ,FNDKBCRZ) 592C 593C Integrals H^U in T3^U 594C 595 596 ISINTRU1 = MULD2H(ISINT1,ISYMRU) 597 ISINTRU2 = MULD2H(ISINT2,ISYMRU) 598C 599 CALL CC3_BARINT(WORK(KT1U),ISYMRU,WORK(KLAMDP), 600 * WORK(KLAMDH),WORK(KEND1),LWRK1, 601 * LU3SRTR,FN3SRTR,LUCKJDRU,FNCKJDRU) 602C 603 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTRU1,LU3SRTR,FN3SRTR, 604 * LUDELDRU,FNDELDRU,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 605 * IDUMMY,CDUMMY) 606C 607 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTRU1, 608 * LUDELDRU,FNDELDRU,LUDKBCRU,FNDKBCRU) 609C 610C------------------------------------------ 611C Calculate the F^Z matrix 612C------------------------------------------ 613C 614 CALL GET_FOCKX(WORK(KFCKZ),LABELRZ,IDLSTRZ,ISYMRZ,WORK(KLAMDP), 615 * ISYM0,WORK(KLAMDH),ISYM0,WORK(KEND1),LWRK1) 616C 617C------------------------------------------ 618C Calculate the F^U matrix 619C------------------------------------------ 620C 621 CALL GET_FOCKX(WORK(KFCKU),LABELRU,IDLSTRU,ISYMRU,WORK(KLAMDP), 622 * ISYM0,WORK(KLAMDH),ISYM0,WORK(KEND1),LWRK1) 623C 624C------------------------------------------ 625C Calculate the [Z,T1^U] matrix 626C Recall that we only need the occ-occ and vir-vir block. 627C------------------------------------------ 628C 629 CALL GET_LAMBDAX(WORK(KLAMDPU),WORK(KLAMDHU),LISTRU,IDLSTRU, 630 * ISYMRU,WORK(KLAMDP),WORK(KLAMDH),WORK(KEND1), 631 * LWRK1) 632 ! get vir-vir block Z_(c-,d) 633 CALL GET_FOCKX(WORK(KFCKZUV),LABELRZ,IDLSTRZ,ISYMRZ,WORK(KLAMDPU), 634 * ISYMRU,WORK(KLAMDH),ISYM0,WORK(KEND1),LWRK1) 635 ! get occ-occ block Z_(l,k-) 636 CALL GET_FOCKX(WORK(KFCKZUO),LABELRZ,IDLSTRZ,ISYMRZ,WORK(KLAMDP), 637 * ISYM0,WORK(KLAMDHU),ISYMRU,WORK(KEND1),LWRK1) 638C 639C------------------------------------------ 640C Calculate the [U,T1^Z] matrix 641C Recall that we only need the occ-occ and vir-vir block. 642C------------------------------------------ 643C 644 CALL GET_LAMBDAX(WORK(KLAMDPZ),WORK(KLAMDHZ),LISTRZ,IDLSTRZ, 645 * ISYMRZ,WORK(KLAMDP),WORK(KLAMDH),WORK(KEND1), 646 * LWRK1) 647 ! get vir-vir block U_(c-,d) 648 CALL GET_FOCKX(WORK(KFCKUZV),LABELRU,IDLSTRU,ISYMRU,WORK(KLAMDPZ), 649 * ISYMRZ,WORK(KLAMDH),ISYM0,WORK(KEND1),LWRK1) 650 ! get occ-occ block U_(l,k-) 651 CALL GET_FOCKX(WORK(KFCKUZO),LABELRU,IDLSTRU,ISYMRU,WORK(KLAMDP), 652 * ISYM0,WORK(KLAMDHZ),ISYMRZ,WORK(KEND1),LWRK1) 653C 654C----------------------------- 655C Read occupied integrals. 656C----------------------------- 657C 658C Memory allocation. 659C 660C--------------------------------------------------------- 661C Work allocation 662C--------------------------------------------------------- 663C 664 ISINT1X = MULD2H(ISINT1,ISYMRX) 665C 666 KRBJIAZU = KEND0 667 KT3OG1 = KRBJIAZU + NT2SQ(ISYRES) 668 KT3OG2 = KT3OG1 + NTRAOC(ISINT2) 669 KEND1 = KT3OG2+ NTRAOC(ISINT2) 670 LWRK1 = LWORK - KEND1 671C 672 KW3ZOGZ1 = KEND1 673 KEND1 = KW3ZOGZ1 + NTRAOC(ISINTRZ2) 674C 675 KW3UOGU1 = KEND1 676 KEND1 = KW3UOGU1 + NTRAOC(ISINTRU2) 677 LWRK1 = LWORK - KEND1 678C 679 KTROC = KEND1 680 KTROC1 = KTROC + NTRAOC(ISINT1X) !KTROC - int. in contraction 681 KEND1 = KTROC1 + NTRAOC(ISINT1X) !KTROC1 - int. in contraction 682 LWRK1 = LWORK - KEND1 683C 684 KINTOC = KEND1 685 KEND1 = KINTOC + NTOTOC(ISINT1) !KINTOC - temporary storage 686 LWRK1 = LWORK - KEND1 687C 688 IF (LWRK1 .LT. 0) THEN 689 WRITE(LUPRI,*) 'Memory available : ',LWORK 690 WRITE(LUPRI,*) 'Memory needed : ',KEND1 691 CALL QUIT('Insufficient space in CC3_AMATT3ZUVIR') 692 END IF 693C 694 CALL DZERO(WORK(KRBJIAZU),NT2SQ(ISYRES)) 695C 696C Occupied integrals. 697C 698C----------------------- 699C Read in integrals for SMAT etc. 700C----------------------- 701C 702 CALL INTOCC_T30(LUCKJD,FNCKJD,WORK(KLAMDP),ISINT2,WORK(KT3OG1), 703 * WORK(KT3OG2),WORK(KEND1),LWRK1) 704C 705C------------------------------------------ 706C Z transformed Occupied integrals. 707C----------------------------------------- 708C 709 CALL INTOCC_T3X(LUCKJDRZ,FNCKJDRZ,WORK(KLAMDP),ISINTRZ2, 710 * WORK(KW3ZOGZ1),WORK(KEND1),LWRK1) 711C 712C------------------------------------------ 713C U transformed Occupied integrals. 714C----------------------------------------- 715C 716 CALL INTOCC_T3X(LUCKJDRU,FNCKJDRU,WORK(KLAMDP),ISINTRU2, 717 * WORK(KW3UOGU1),WORK(KEND1),LWRK1) 718C 719C----------------------------------- 720C Occupied integrals in contraction 721C----------------------------------- 722C 723 IOFF = 1 724 IF (NTOTOC(ISINT1) .GT. 0) THEN 725 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT1)) 726 ENDIF 727C 728 !Write out norms of integrals. 729 IF (IPRINT .GT. 55) THEN 730 XNORMVAL = DDOT(NTOTOC(ISINT1),WORK(KINTOC),1,WORK(KINTOC),1) 731 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XNORMVAL 732 ENDIF 733C 734 !Transform (ia|j delta) integrals to (ia|j k) and sort as (i,j,k,a) 735 CALL CCLR_TROCC(WORK(KINTOC),WORK(KTROC),WORK(KLAMDHX),ISYMRX, 736 * WORK(KEND1),LWRK1) 737C 738 !sort (i,j,k,a) as (a,i,j,k) 739 CALL CCSDT_SRTOC2(WORK(KTROC),WORK(KTROC1),ISINT1X, 740 * WORK(KEND1),LWRK1) 741c 742C 743C---------------------------- 744C General loop structure. 745C---------------------------- 746C 747 DO ISYMD = 1,NSYM 748 749 ISCKB1 = MULD2H(ISINT1,ISYMD) 750 ISCKB1X = MULD2H(ISINT1X,ISYMD) 751 ISCKB2 = MULD2H(ISINT2,ISYMD) 752 ISCKB2Z = MULD2H(ISINTRZ2,ISYMD) 753 ISCKB2U = MULD2H(ISINTRU2,ISYMD) 754C 755 IF (IPRINT .GT. 55) THEN 756C 757 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISCKB2 :',ISCKB2 758 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISCKB2Z:',ISCKB2Z 759 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISCKB2U:',ISCKB2U 760 761 ENDIF 762C 763C-------------------------- 764C Memory allocation. 765C-------------------------- 766C 767 KT3VDG3 = KEND1 768 KT3VDG2 = KT3VDG3 + NCKATR(ISCKB2) 769 KEND2 = KT3VDG2 + NCKATR(ISCKB2) 770 LWRK2 = LWORK - KEND2 771 772 KT3VDG1 = KEND2 773 KEND3 = KT3VDG1 + NCKATR(ISCKB2) 774 LWRK3 = LWORK - KEND3 775C 776 KW3ZVDGZ1 = KEND3 777 KEND3 = KW3ZVDGZ1 + NCKATR(ISCKB2Z) 778 LWRK3 = LWORK - KEND3 779C 780 KW3UVDGU1 = KEND3 781 KEND3 = KW3UVDGU1 + NCKATR(ISCKB2U) 782 LWRK3 = LWORK - KEND3 783C 784 KT3VDGZ3 = KEND3 785 KEND3 = KT3VDGZ3 + NCKATR(ISCKB2Z) 786C 787 KT3VDGU3 = KEND3 788 KEND3 = KT3VDGU3 + NCKATR(ISCKB2U) 789C 790 KTRVI = KEND3 791 KTRVI1 = KTRVI + NCKATR(ISCKB1X) !KTRVI - int. in contraction 792 KEND3 = KTRVI1 + NCKATR(ISCKB1X) !KTRVI1 - int. in contraction 793 LWRK3 = LWORK - KEND3 794C 795 KINTVI = KEND3 796 MAXZU = MAX(NCKA(ISCKB2Z),NCKA(ISCKB2U)) 797 KEND3 = KINTVI + MAX(MAXZU,NCKA(ISCKB1)) 798 LWRK3 = LWORK - KEND3 799C 800 IF (LWRK3 .LT. 0) THEN 801 WRITE(LUPRI,*) 'Memory available : ',LWORK 802 WRITE(LUPRI,*) 'Memory needed : ',KEND3 803 CALL QUIT('Insufficient space in CC3_AMATT3ZUVIR') 804 END IF 805C 806C--------------------- 807C Sum over D 808C--------------------- 809C 810 DO D = 1,NVIR(ISYMD) 811C 812C----------------------------------------------- 813C Read virtual integrals used in first s3am. 814C----------------------------------------------- 815C 816 CALL INTVIR_T30_D(LUDKBC,FNDKBC,LUDELD,FNDELD,ISINT2, 817 * WORK(KT3VDG1),WORK(KT3VDG2), 818 * WORK(KT3VDG3),WORK(KLAMDH),ISYMD,D, 819 * WORK(KEND3),LWRK3) 820C 821C----------------------------------------------------------------------- 822C Read Z transformed virtual integrals (H^Z) used for W^Z 823C----------------------------------------------------------------------- 824C 825 IOFF = ICKBD(ISCKB2Z,ISYMD) + NCKATR(ISCKB2Z)*(D - 1) + 1 826 IF (NCKATR(ISCKB2Z) .GT. 0) THEN 827 CALL GETWA2(LUDKBCRZ,FNDKBCRZ,WORK(KW3ZVDGZ1),IOFF, 828 & NCKATR(ISCKB2Z)) 829 ENDIF 830C 831C----------------------------------------------------------------------- 832C Read U transformed virtual integrals (H^U) used for W^U 833C----------------------------------------------------------------------- 834C 835 IOFF = ICKBD(ISCKB2U,ISYMD) + NCKATR(ISCKB2U)*(D - 1) + 1 836 IF (NCKATR(ISCKB2U) .GT. 0) THEN 837 CALL GETWA2(LUDKBCRU,FNDKBCRU,WORK(KW3UVDGU1),IOFF, 838 & NCKATR(ISCKB2U)) 839 ENDIF 840C 841C-------------------------------------------------------------------- 842C Read virtual integrals [H,T1Z] where Z is LISTRZ (used in W^Z) 843C-------------------------------------------------------------------- 844C 845 IF (NCKA(ISCKB2Z) .GT. 0) THEN 846 IOFF = ICKAD(ISCKB2Z,ISYMD) + 847 & NCKA(ISCKB2Z)*(D - 1) + 1 848 CALL GETWA2(LUDELDRZ,FNDELDRZ,WORK(KINTVI),IOFF, 849 * NCKA(ISCKB2Z)) 850 ENDIF 851C 852 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KT3VDGZ3), 853 * WORK(KLAMDH),ISYMD,D,ISINTRZ2, 854 * WORK(KEND3),LWRK3) 855C 856C-------------------------------------------------------------------- 857C Read virtual integrals [H,T1U] where U is LISTRU (used in W^U) 858C-------------------------------------------------------------------- 859C 860 IF (NCKA(ISCKB2U) .GT. 0) THEN 861 IOFF = ICKAD(ISCKB2U,ISYMD) + 862 & NCKA(ISCKB2U)*(D - 1) + 1 863 CALL GETWA2(LUDELDRU,FNDELDRU,WORK(KINTVI),IOFF, 864 * NCKA(ISCKB2U)) 865 ENDIF 866C 867 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KT3VDGU3), 868 * WORK(KLAMDH),ISYMD,D,ISINTRU2, 869 * WORK(KEND3),LWRK3) 870C 871C----------------------------------------------------- 872C Virtual integrals in pojection (fixed D) 873C----------------------------------------------------- 874C 875 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 876 IF (NCKA(ISCKB1) .GT. 0) THEN 877 CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF, 878 * NCKA(ISCKB1)) 879 ENDIF 880 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVI),WORK(KLAMDPX), 881 * ISYMRX, 882 * ISYMD,D,ISYMOP,WORK(KEND3),LWRK3) 883C 884 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 885 IF (NCKA(ISCKB1) .GT. 0) THEN 886 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 887 * NCKA(ISCKB1)) 888 ENDIF 889 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI1),WORK(KLAMDPX), 890 * ISYMRX, 891 * ISYMD,D,ISYMOP,WORK(KEND3),LWRK3) 892C 893C-------------------------------------------------------- 894C Sum over ISYMB 895C-------------------------------------------------------- 896C 897 DO ISYMB = 1,NSYM 898 899 ISYALJ = MULD2H(ISYMB,ISYMT2) 900 ISYALJ2 = MULD2H(ISYMD,ISYMT2) 901 ISYMBD = MULD2H(ISYMB,ISYMD) 902 ISCKIJ = MULD2H(ISYMBD,ISYMT3) 903 ISCKD2 = MULD2H(ISINT2,ISYMB) 904 ISCKD2Z = MULD2H(ISINTRZ2,ISYMB) 905 ISCKD2U = MULD2H(ISINTRU2,ISYMB) 906 ISWMATZ = MULD2H(ISYMRZ,ISCKIJ) 907 ISWMATU = MULD2H(ISYMRU,ISCKIJ) 908 ISWMATZU = MULD2H(ISWMATZ,ISYMRU) 909 910 IF (IPRINT .GT. 55) THEN 911 912 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISYMD :',ISYMD 913 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISYMB :',ISYMB 914 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISYALJ :',ISYALJ 915 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISYALJ2:',ISYALJ2 916 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISYMBD :',ISYMBD 917 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISCKIJ :',ISCKIJ 918 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISWMATZ :',ISWMATZ 919 WRITE(LUPRI,*) 'In CC3_AMATT3ZUVIR: ISWMATU :',ISWMATU 920 921 ENDIF 922C 923 KSMAT = KEND3 924 KSMAT3 = KSMAT + NCKIJ(ISCKIJ) 925 KUMAT = KSMAT3 + NCKIJ(ISCKIJ) 926 KUMAT3 = KUMAT + NCKIJ(ISCKIJ) 927 KDIAG = KUMAT3 + NCKIJ(ISCKIJ) 928 KT3MAT = KDIAG + NCKIJ(ISCKIJ) 929 KEND3 = KT3MAT + NCKIJ(ISCKIJ) 930C 931 KDIAGWZ = KEND3 932 KWMATZ = KDIAGWZ + NCKIJ(ISWMATZ) 933 KINDSQWZ = KWMATZ + NCKIJ(ISWMATZ) 934 KEND3 = KINDSQWZ + (6*NCKIJ(ISWMATZ) - 1)/IRAT + 1 935C 936 KDIAGWU = KEND3 937 KWMATU = KDIAGWU + NCKIJ(ISWMATU) 938 KINDSQWU = KWMATU + NCKIJ(ISWMATU) 939 KEND3 = KINDSQWU + (6*NCKIJ(ISWMATU) - 1)/IRAT + 1 940C 941 KINDSQ = KEND3 942 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 943 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 944 KEND3 = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 945C 946 MAXZU = MAX(NCKIJ(ISWMATZ),NCKIJ(ISWMATU)) 947 MAXZZUU = MAX(MAXZU,NCKIJ(ISWMATZU)) 948C 949 KTMATW = KEND3 950 KT3VBG1 = KTMATW + MAX(MAXZZUU,NCKIJ(ISCKIJ)) 951 KT3VBG2 = KT3VBG1 + NCKATR(ISCKD2) 952 KT3VBG3 = KT3VBG2 + NCKATR(ISCKD2) 953 KEND4 = KT3VBG3 + NCKATR(ISCKD2) 954 LWRK4 = LWORK - KEND4 955C 956 KW3ZVDGZ2 = KEND4 957 KEND4 = KW3ZVDGZ2 + NCKATR(ISCKD2Z) 958C 959 KW3UVDGU2 = KEND4 960 KEND4 = KW3UVDGU2 + NCKATR(ISCKD2U) 961C 962 KT3VBGZ3 = KEND4 963 KEND4 = KT3VBGZ3 + NCKATR(ISCKD2Z) 964C 965 KT3VBGU3 = KEND4 966 KEND4 = KT3VBGU3 + NCKATR(ISCKD2U) 967C 968 KINTVI = KEND4 969 KEND4 = KINTVI + MAX(NCKA(ISCKD2Z),NCKA(ISCKD2U)) 970 LWRK4 = LWORK - KEND4 971C 972 KINDSQWZU = KEND4 973 KDIAGWZU = KINDSQWZU + (6*NCKIJ(ISWMATZU) - 1)/IRAT + 1 974 KWMATZU = KDIAGWZU + NCKIJ(ISWMATZU) 975 KEND4 = KWMATZU + NCKIJ(ISWMATZU) 976 LWRK4 = LWORK - KEND4 977C 978 IF (LWRK4 .LT. 0) THEN 979 WRITE(LUPRI,*) 'Memory available : ',LWORK 980 WRITE(LUPRI,*) 'Memory needed : ',KEND4 981 CALL QUIT('Insufficient space in CC3_AMATT3ZUVIR') 982 END IF 983C 984C--------------------------------------------- 985C Construct part of the diagonal. 986C--------------------------------------------- 987C 988 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 989 CALL CC3_DIAG(WORK(KDIAGWZ),WORK(KFOCKD),ISWMATZ) 990 CALL CC3_DIAG(WORK(KDIAGWU),WORK(KFOCKD),ISWMATU) 991 CALL CC3_DIAG(WORK(KDIAGWZU),WORK(KFOCKD),ISWMATZU) 992C 993 IF (IPRINT .GT. 55) THEN 994 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1, 995 * WORK(KDIAG),1) 996 WRITE(LUPRI,*) 'Norm of DIA ',XNORMVAL 997 ENDIF 998C 999C------------------------------------- 1000C Construct index arrays. 1001C---------------------------------- 1002C 1003 LENSQ = NCKIJ(ISCKIJ) 1004 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 1005 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 1006 CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2) 1007 LENSQWZ = NCKIJ(ISWMATZ) 1008 CALL CC3_INDSQ(WORK(KINDSQWZ),LENSQWZ,ISWMATZ) 1009 LENSQWU = NCKIJ(ISWMATU) 1010 CALL CC3_INDSQ(WORK(KINDSQWU),LENSQWU,ISWMATU) 1011 LENSQWZU = NCKIJ(ISWMATZU) 1012 CALL CC3_INDSQ(WORK(KINDSQWZU),LENSQWZU,ISWMATZU) 1013 1014 1015 DO B = 1,NVIR(ISYMB) 1016C 1017C------------------------------------ 1018C Initialise 1019C--------------------------------------- 1020C 1021 CALL DZERO(WORK(KWMATZ),NCKIJ(ISWMATZ)) 1022 CALL DZERO(WORK(KWMATU),NCKIJ(ISWMATU)) 1023 CALL DZERO(WORK(KWMATZU),NCKIJ(ISWMATZU)) 1024C 1025C----------------------------------------------- 1026C Read virtual integrals used in second s3am. 1027C----------------------------------------------- 1028C 1029 CALL INTVIR_T30_D(LUDKBC,FNDKBC,LUDELD,FNDELD,ISINT2, 1030 * WORK(KT3VBG1),WORK(KT3VBG2), 1031 * WORK(KT3VBG3),WORK(KLAMDH),ISYMB,B, 1032 * WORK(KEND4),LWRK4) 1033C 1034C-------------------------------------------------------------------- 1035C Read Z transformed virtual integrals (H^Z) used in W^Z 1036C-------------------------------------------------------------------- 1037C 1038 IOFF = ICKBD(ISCKD2Z,ISYMB) + 1039 & NCKATR(ISCKD2Z)*(B - 1) + 1 1040 IF (NCKATR(ISCKD2Z) .GT. 0) THEN 1041 CALL GETWA2(LUDKBCRZ,FNDKBCRZ,WORK(KW3ZVDGZ2),IOFF, 1042 & NCKATR(ISCKD2Z)) 1043 ENDIF 1044C 1045C-------------------------------------------------------------------- 1046C Read U transformed virtual integrals (H^U) used in W^U 1047C-------------------------------------------------------------------- 1048C 1049 IOFF = ICKBD(ISCKD2U,ISYMB) + 1050 & NCKATR(ISCKD2U)*(B - 1) + 1 1051 IF (NCKATR(ISCKD2U) .GT. 0) THEN 1052 CALL GETWA2(LUDKBCRU,FNDKBCRU,WORK(KW3UVDGU2),IOFF, 1053 & NCKATR(ISCKD2U)) 1054 ENDIF 1055C 1056C-------------------------------------------------------------------- 1057C Read virtual integrals [H,T1Z] where Z is LISTRZ (used in W^Z) 1058C-------------------------------------------------------------------- 1059C 1060 IF (NCKA(ISCKD2Z) .GT. 0) THEN 1061 IOFF = ICKAD(ISCKD2Z,ISYMB) + 1062 & NCKA(ISCKD2Z)*(B - 1) + 1 1063 CALL GETWA2(LUDELDRZ,FNDELDRZ,WORK(KINTVI),IOFF, 1064 * NCKA(ISCKD2Z)) 1065 ENDIF 1066C 1067 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KT3VBGZ3), 1068 * WORK(KLAMDH),ISYMB,B,ISINTRZ2, 1069 * WORK(KEND4),LWRK4) 1070C 1071C-------------------------------------------------------------------- 1072C Read virtual integrals [H,T1U] where U is LISTRU (used in W^U) 1073C-------------------------------------------------------------------- 1074C 1075 IF (NCKA(ISCKD2U) .GT. 0) THEN 1076 IOFF = ICKAD(ISCKD2U,ISYMB) + 1077 & NCKA(ISCKD2U)*(B - 1) + 1 1078 CALL GETWA2(LUDELDRU,FNDELDRU,WORK(KINTVI),IOFF, 1079 * NCKA(ISCKD2U)) 1080 ENDIF 1081C 1082 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KT3VBGU3), 1083 * WORK(KLAMDH),ISYMB,B,ISINTRU2, 1084 * WORK(KEND4),LWRK4) 1085C 1086C------------------------------------------------------------------------ 1087C Calculate the S(ci,bk,dj) matrix for T3 for B,D. 1088C------------------------------------------------------------------- 1089C 1090 CALL GET_T30_BD(ISYMT3,ISINT2,WORK(KT2TP),ISYMT2, 1091 * WORK(KT3MAT),WORK(KFOCKD),WORK(KDIAG), 1092 * WORK(KINDSQ),LENSQ,WORK(KSMAT), 1093 * WORK(KT3VDG1),WORK(KT3VDG2), 1094 * WORK(KT3OG1),WORK(KINDEX), 1095 * WORK(KSMAT3),WORK(KT3VBG1), 1096 * WORK(KT3VBG2),WORK(KINDEX2), 1097 * WORK(KUMAT),WORK(KT3VDG3), 1098 * WORK(KT3OG2),WORK(KUMAT3), 1099 * WORK(KT3VBG3),ISYMB,B,ISYMD,D, 1100 * ISCKIJ,WORK(KEND4),LWRK4) 1101 1102c 1103c call sum_pt3(work(KT3MAT),isymb,b,isymd,d, 1104c * 1,work(kx3am),3) 1105 1106C 1107C------------------------------------------------------ 1108C Calculate WBD^Z 1109C------------------------------------------------------ 1110C 1111 1112 ! Copy T3^0 from KT3MAT to KTMATW. KTMATW used later 1113 ! as help array and KTMATW overwritten. 1114c CALL DCOPY(NCKIJ(ISCKIJ),WORK(KT3MAT),1, 1115c * WORK(KTMATW),1) 1116C 1117 1118C 1119C------------------------------------------------------ 1120C Calculate the term <mu3|[Z,T30]|HF> occupied contribution 1121C added in W^BD(a,i-,k-,j-), where i- means one-index transformation of i 1122C (the three contributions are added: fixed BD, fixed aB, fixed Da 1123C added in W^BD). 1124C------------------------------------------------------ 1125C 1126 AIBJCK_PERM = 4 1127C write(lupri,*)'entering WXBD_O... ' 1128 CALL WXBD_O(AIBJCK_PERM,WORK(KT3MAT),ISCKIJ, 1129 * WORK(KFCKZ),ISYMRZ, 1130 * WORK(KWMATZ),ISWMATZ,WORK(KEND4),LWRK4) 1131C 1132C------------------------------------------------------ 1133C Calculate the term <mu3|[[Z,T2],T2]|HF> 1134C added in W^BD(a,i-,k-,j-). 1135C------------------------------------------------------ 1136C 1137C write(lupri,*)'entering WXBD_T2... ' 1138 CALL WXBD_T2(AIBJCK_PERM,B,ISYMB,D,ISYMD,WORK(KT2TP), 1139 * ISYMT2,WORK(KFCKZ), 1140 * ISYMRZ,WORK(KINDSQWZ),LENSQWZ, 1141 * WORK(KWMATZ),ISWMATZ,WORK(KEND4),LWRK4) 1142 1143C 1144C---------------------------------------------------------------- 1145C Calculate the terms <mu3|[H^Z,T2]|HF> + <mu3|[H,T2^Z]|HF> 1146C added in W^BD(a,i-,k-,j-). 1147C---------------------------------------------------------------- 1148C 1149C write(lupri,*)'entering WXBD_GROUND(1)... ' 1150 CALL WXBD_GROUND(AIBJCK_PERM, 1151 * WORK(KT2Z),ISYMRZ,WORK(KTMATW), 1152 * WORK(KT3VDG1),WORK(KT3VBG1),WORK(KT3VBG3), 1153 * WORK(KT3VDG3), 1154 * WORK(KT3OG1),ISINT2, 1155 * WORK(KWMATZ),WORK(KEND4),LWRK4, 1156 * WORK(KINDSQWZ),LENSQWZ, 1157 * ISYMB,B,ISYMD,D) 1158C 1159C write(lupri,*)'entering WXBD_GROUND(2)... ' 1160 CALL WXBD_GROUND(AIBJCK_PERM, 1161 * WORK(KT2TP),ISYMT2,WORK(KTMATW), 1162 * WORK(KW3ZVDGZ1),WORK(KW3ZVDGZ2), 1163 * WORK(KT3VBGZ3),WORK(KT3VDGZ3), 1164 * WORK(KW3ZOGZ1),ISINTRZ2, 1165 * WORK(KWMATZ),WORK(KEND4),LWRK4, 1166 * WORK(KINDSQWZ),LENSQWZ, 1167 * ISYMB,B,ISYMD,D) 1168C 1169C------------------------------------------------ 1170C Divide by the energy difference and 1171C remove the forbidden elements 1172C------------------------------------------------ 1173C 1174 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQRZ,ISWMATZ, 1175 * WORK(KWMATZ),WORK(KDIAGWZ),WORK(KFOCKD)) 1176 CALL T3_FORBIDDEN(WORK(KWMATZ),ISYMRZ,ISYMB,B, 1177 * ISYMD,D) 1178 1179c 1180c call sum_pt3(work(KWMATZ),isymb,b,isymd,d, 1181c * ISWMATZ,work(kx3am),5) 1182 1183 1184C 1185c CALL GET_T3X_BD(ISYMRZ,WORK(KWMATZ),ISWMATZ, 1186c * WORK(KTMATW),ISCKIJ,WORK(KFCKZ), 1187c * ISYMRZ,WORK(KT2TP),ISYMT2, 1188c * WORK(KT2Z),ISYMRZ,WORK(KT3VDG1), 1189c * WORK(KT3VBG1),WORK(KT3OG1), 1190c * ISINT2,WORK(KW3ZVDGZ1), 1191c * WORK(KW3ZVDGZ2),WORK(KW3ZOGZ1), 1192c * ISINTRZ2,FREQRZ,WORK(KDIAGWZ), 1193c * WORK(KFOCKD),WORK(KINDSQWZ),LENSQWZ, 1194c * B,ISYMB,D,ISYMD,WORK(KEND4),LWRK4) 1195 1196C Calculate WUZ intermmediate for <mu3|[U,wZ^{aBC}_{ijk}]|HF> 1197C where wZ^{aBC}_{ijk} is part of T3^Z with the virtual 1198C contribution (WXBD_V routine) taken out. 1199 1200c 1201 CALL WBD_O(WORK(KWMATZ),ISWMATZ,WORK(KFCKU),ISYMRU, 1202 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1203C 1204 CALL WBD_V(WORK(KWMATZ),ISWMATZ,WORK(KFCKU),ISYMRU, 1205 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1206C 1207C Calculate <mu3|[[U,T1Z],T30]|HF> 1208C 1209 CALL WBD_O(WORK(KT3MAT),ISCKIJ,WORK(KFCKUZO),ISYMZU, 1210 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1211C 1212 CALL WBD_V(WORK(KT3MAT),ISCKIJ,WORK(KFCKUZV),ISYMZU, 1213 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1214C 1215C Calculate <mu3|[[U,T2Z],T20]|HF> 1216C 1217 T2XNET2Z = .TRUE. 1218 CALL WBD_T2(T2XNET2Z,B,ISYMB,D,ISYMD, 1219 * WORK(KT2Z),ISYMRZ,WORK(KT2TP),ISYMT2, 1220 * WORK(KFCKU),ISYMRU, 1221 * WORK(KINDSQWZU),LENSQWZU,WORK(KWMATZU), 1222 * ISWMATZU,WORK(KEND4),LWRK4) 1223C 1224C------------------------------------------------ 1225C Divide by the energy difference and 1226C remove the forbidden elements (here only for debugging) 1227C------------------------------------------------ 1228C 1229c CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQZU,ISWMATZU, 1230c * WORK(KWMATZU),WORK(KDIAGWZU),WORK(KFOCKD)) 1231c CALL T3_FORBIDDEN(WORK(KWMATZU),ISYMZU,ISYMB,B, 1232c * ISYMD,D) 1233 1234 1235c call sum_pt3(work(KWMATZU),isymb,b,isymd,d, 1236c * ISWMATZU,work(kx3am),5) 1237 1238 1239C 1240C------------------------------------------------------ 1241C Calculate WBD^U 1242C------------------------------------------------------ 1243C 1244 1245 ! Copy T3^0 from KT3MAT to KTMATW. KTMATW used later 1246 ! as help array and KTMATW overwritten. 1247c CALL DCOPY(NCKIJ(ISCKIJ),WORK(KT3MAT),1, 1248c * WORK(KTMATW),1) 1249C 1250C 1251C------------------------------------------------------ 1252C Calculate the term <mu3|[U,T30]|HF> occupied contribution 1253C added in W^BD(a,i-,k-,j-), where i- means one-index transformation of i 1254C (the three contributions are added: fixed BD, fixed aB, fixed Da 1255C added in W^BD). 1256C------------------------------------------------------ 1257C 1258 AIBJCK_PERM = 4 ! transform all occupied indecies 1259c write(lupri,*)'entering WXBD_O... U' 1260 CALL WXBD_O(AIBJCK_PERM,WORK(KT3MAT),ISCKIJ, 1261 * WORK(KFCKU),ISYMRU, 1262 * WORK(KWMATU),ISWMATU,WORK(KEND4),LWRK4) 1263C 1264C------------------------------------------------------ 1265C Calculate the term <mu3|[[U,T2],T2]|HF> 1266C added in W^BD(a,i-,k-,j-). 1267C------------------------------------------------------ 1268C 1269c write(lupri,*)'entering WXBD_T2... U' 1270 CALL WXBD_T2(AIBJCK_PERM,B,ISYMB,D,ISYMD,WORK(KT2TP), 1271 * ISYMT2,WORK(KFCKU), 1272 * ISYMRU,WORK(KINDSQWU),LENSQWU, 1273 * WORK(KWMATU),ISWMATU,WORK(KEND4),LWRK4) 1274 1275C 1276C---------------------------------------------------------------- 1277C Calculate the terms <mu3|[H^U,T2]|HF> + <mu3|[H,T2^U]|HF> 1278C added in W^BD(a,i-,k-,j-). 1279C---------------------------------------------------------------- 1280C 1281c write(lupri,*)'entering WXBD_GROUND(1)... U' 1282 CALL WXBD_GROUND(AIBJCK_PERM, 1283 * WORK(KT2U),ISYMRU,WORK(KTMATW), 1284 * WORK(KT3VDG1),WORK(KT3VBG1),WORK(KT3VBG3), 1285 * WORK(KT3VDG3), 1286 * WORK(KT3OG1),ISINT2, 1287 * WORK(KWMATU),WORK(KEND4),LWRK4, 1288 * WORK(KINDSQWU),LENSQWU, 1289 * ISYMB,B,ISYMD,D) 1290C 1291c write(lupri,*)'entering WXBD_GROUND(2)... U' 1292 CALL WXBD_GROUND(AIBJCK_PERM, 1293 * WORK(KT2TP),ISYMT2,WORK(KTMATW), 1294 * WORK(KW3UVDGU1),WORK(KW3UVDGU2), 1295 * WORK(KT3VBGU3),WORK(KT3VDGU3), 1296 * WORK(KW3UOGU1),ISINTRU2, 1297 * WORK(KWMATU),WORK(KEND4),LWRK4, 1298 * WORK(KINDSQWU),LENSQWU, 1299 * ISYMB,B,ISYMD,D) 1300C 1301C------------------------------------------------ 1302C Divide by the energy difference and 1303C remove the forbidden elements 1304C------------------------------------------------ 1305C 1306 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQRU,ISWMATU, 1307 * WORK(KWMATU),WORK(KDIAGWU),WORK(KFOCKD)) 1308 CALL T3_FORBIDDEN(WORK(KWMATU),ISYMRU,ISYMB,B, 1309 * ISYMD,D) 1310 1311c 1312c call sum_pt3(work(KWMATU),isymb,b,isymd,d, 1313c * ISWMATU,work(kx3am),5) 1314 1315 1316c CALL GET_T3X_BD(ISYMRU,WORK(KWMATU),ISWMATU, 1317c * WORK(KTMATW),ISCKIJ,WORK(KFCKU), 1318c * ISYMRU,WORK(KT2TP),ISYMT2, 1319c * WORK(KT2U),ISYMRU,WORK(KT3VDG1), 1320c * WORK(KT3VBG1),WORK(KT3OG1), 1321c * ISINT2,WORK(KW3UVDGU1), 1322c * WORK(KW3UVDGU2),WORK(KW3UOGU1), 1323c * ISINTRU2,FREQRU,WORK(KDIAGWU), 1324c * WORK(KFOCKD),WORK(KINDSQWU),LENSQWU, 1325c * B,ISYMB,D,ISYMD,WORK(KEND4),LWRK4) 1326 1327C Calculate WZU intermmediate for <mu3|[Z,wY^{aBC}_{ijk}]|HF> 1328C where wY^{aBC}_{ijk} is part of T3^U with the virtual 1329C contribution (WXBD_V routine) taken out. 1330 1331 CALL WBD_O(WORK(KWMATU),ISWMATU,WORK(KFCKZ),ISYMRZ, 1332 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1333C 1334 CALL WBD_V(WORK(KWMATU),ISWMATU,WORK(KFCKZ),ISYMRZ, 1335 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1336C 1337C Calculate <mu3|[[Z,T1U],T30]|HF> 1338C 1339 CALL WBD_O(WORK(KT3MAT),ISCKIJ,WORK(KFCKZUO),ISYMZU, 1340 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1341C 1342 CALL WBD_V(WORK(KT3MAT),ISCKIJ,WORK(KFCKZUV),ISYMZU, 1343 * WORK(KWMATZU),ISWMATZU,WORK(KEND4),LWRK4) 1344C 1345C Calculate <mu3|[[Z,T2U],T20]|HF> 1346C 1347 T2XNET2Z = .TRUE. 1348 CALL WBD_T2(T2XNET2Z,B,ISYMB,D,ISYMD, 1349 * WORK(KT2U),ISYMRU,WORK(KT2TP),ISYMT2, 1350 * WORK(KFCKZ),ISYMRZ, 1351 * WORK(KINDSQWZU),LENSQWZU,WORK(KWMATZU), 1352 * ISWMATZU,WORK(KEND4),LWRK4) 1353C 1354C------------------------------------------------ 1355C Divide by the energy difference and 1356C remove the forbidden elements 1357C------------------------------------------------ 1358C 1359 CALL T3_FORBIDDEN(WORK(KWMATZU),ISYMZU,ISYMB,B, 1360 * ISYMD,D) 1361 1362c call sum_pt3(work(KWMATZU),isymb,b,isymd,d, 1363c * ISWMATZU,work(kx3am),4) 1364 1365 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQZU,ISWMATZU, 1366 * WORK(KWMATZU),WORK(KDIAGWZU),WORK(KFOCKD)) 1367 1368 1369C 1370C----------------------------------------------- 1371C Carry out the contractions... 1372C----------------------------------------------- 1373C 1374C 1375C------------------------------------------------------ 1376C Calculate the term <mu2|[H,WBD^ZU(3)]|HF> 1377C ( Fock matrix cont ) 1378C added in OMEGA2EFF 1379C------------------------------------------------------ 1380C 1381 CALL CC3_CY2F(OMEGA2EFF,ISYRES,WORK(KWMATZU),ISWMATZU, 1382 * WORK(KTMATW),WORK(KFOCKX),ISYFCKX, 1383 * WORK(KINDSQWZU), 1384 * LENSQWZU,WORK(KEND4),LWRK4, 1385 * ISYMB,B,ISYMD,D) 1386c 1387C------------------------------------------------------ 1388C Calculate the term <mu2|[H,WBD^ZU(3)]|HF> 1389C ( Occupied cont ) 1390C added in OMEGA2EFF 1391C------------------------------------------------------ 1392C 1393 CALL CC3_CY2O(OMEGA2EFF,ISYRES,WORK(KWMATZU),ISWMATZU, 1394 * WORK(KTMATW),WORK(KTROC), 1395 * WORK(KTROC1),ISINT1,WORK(KEND4),LWRK4, 1396 * WORK(KINDSQWZU),LENSQWZU, 1397 * ISYMB,B,ISYMD,D) 1398C 1399C------------------------------------------------------ 1400C Calculate the term <mu2|[H,WBD^ZU(3)]|HF> 1401C ( Virtual cont ) 1402C added in OMEGA2EFF 1403C------------------------------------------------------ 1404C 1405 CALL CC3_CY2V(OMEGA2EFF,ISYRES,WORK(KRBJIAZU), 1406 * WORK(KWMATZU), 1407 * ISWMATZU,WORK(KTMATW),WORK(KTRVI), 1408 * WORK(KTRVI1),ISINT1,WORK(KEND4),LWRK4, 1409 * WORK(KINDSQWZU),LENSQWZU, 1410 * ISYMB,B,ISYMD,D) 1411C 1412 END DO !B 1413 END DO !ISYMB 1414C 1415 END DO !D 1416 END DO !ISYMD 1417c 1418c 1419c write(lupri,*) 'KWMATZ in CC3_AAMAT : ' 1420c call print_pt3(work(kx3am),1,4) 1421c 1422 1423C 1424C 1425C------------------------------------------------------ 1426C Accumulate RBJIA from <mu2|[H,WBD^ZU(3)]|HF> ( Virtual cont ) 1427C in OMEGA2EFF 1428C------------------------------------------------------ 1429C 1430 CALL CC3_RBJIA(OMEGA2EFF,ISYRES,WORK(KRBJIAZU)) 1431C 1432 IF (IPRINT .GT. 55) THEN 1433 XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 1434 WRITE(LUPRI,*)'Norm of OMEGA2EFF final ',XNORMVAL 1435 ENDIF 1436C 1437C 1438C 1439c write(lupri,*) 'cont to t3x in CC3_AMATT3ZUVIR' 1440c call print_pt3(work(kx3am),1,4) 1441C write(lupri,*) 'The summed S terms : ' 1442C call print_pt3(work(kx3am),1,1) 1443C 1444C-------------------------------- 1445C Close files for "response" 1446C-------------------------------- 1447C 1448 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 1449 CALL WCLOSE2(LUCKJDRZ,FNCKJDRZ,'DELETE') 1450 CALL WCLOSE2(LUCKJDRU,FNCKJDRU,'DELETE') 1451 CALL WCLOSE2(LUDELDRZ,FNDELDRZ,'DELETE') 1452 CALL WCLOSE2(LUDELDRU,FNDELDRU,'DELETE') 1453 CALL WCLOSE2(LUDKBCRZ,FNDKBCRZ,'DELETE') 1454 CALL WCLOSE2(LUDKBCRU,FNDKBCRU,'DELETE') 1455C 1456C------------- 1457C End 1458C------------- 1459C 1460 CALL QEXIT('AT3ZUV') 1461C 1462 RETURN 1463 END 1464C /* Deck cc3_amatt3zuocc */ 1465 SUBROUTINE CC3_AMATT3ZUOCC(LISTX,IDLSTX,LISTZU,IDLSTZU, 1466 & OMEGA2, 1467 & OMEGA2EFF, 1468 * ISYRES,WORK,LWORK, 1469 * LUCKJD,FNCKJD,LUTOC,FNTOC, 1470 * LUDELD,FNDELD,LU3FOP,FN3FOP, 1471 * LU3VI,FN3VI) 1472* 1473******************************************************************* 1474* 1475* This routine calculates the contribution to cubic response 1476* (originating form both F matrix and B matrix): 1477* 1478* omega2eff = <mu2|[H^X,T3^ZU]|HF> 1479* 1480* which is contracted outside with the multipliers (zeroth-order for 1481* F matrix or first-order for B matrix). 1482* 1483* 1484* 1485* T3^ZU is calculated using W^JK intermmediate. 1486* 1487* !!! ACTUALLY HERE WE CALCULATE ONLY THE PART OF THE TERM, WHICH ORIGINATES 1488* FROM THE A^U-MATRIX CONTRIBUTION TO T3^ZU 1489* (to get the rest of contribution call cc3_bmatt3zu routine) !!! 1490* 1491* Thus here: 1492* 1493* W^BD = <mu3|[Z,T3^U]|HF> 1494* 1495* (!!! TO GET THE TOTAL CONTRIBUTION FROM THIS TERM AND THE OTHER 1496* A^U-MATRIX TERMS CALL CC3_AMATT3ZUVIR !!!) 1497* 1498* 1499* F. Pawlowski, P. Jorgensen, Aarhus Spring 2003. 1500* 1501************************************************************** 1502C 1503 IMPLICIT NONE 1504#include "priunit.h" 1505#include "ccsdsym.h" 1506#include "ccl1rsp.h" 1507#include "ccr1rsp.h" 1508#include "ccorb.h" 1509#include "dummy.h" 1510#include "iratdef.h" 1511#include "ccinftap.h" 1512#include "ccsdinp.h" 1513C 1514 CHARACTER*(*) FNCKJD,FNTOC,FNDELD,FN3FOP,FN3VI 1515 INTEGER LUCKJD,LUTOC,LUDELD,LU3FOP,LU3VI 1516C 1517 CHARACTER CDUMMY*1 1518C 1519 CHARACTER*3 LISTX, LISTRZU, LISTRZ, LISTRU 1520 CHARACTER*8 LABELRZ,LABELRU 1521 INTEGER IDLSTX,IDLSTRZU,IDLSTRZ,IDLSTRU 1522 INTEGER ISYMRZ, ISYMRU 1523C 1524 INTEGER ISYRES, LWORK 1525 INTEGER ISYM0,ISYMT2,ISYMT3 1526 INTEGER KT2TP,KFOCKD,KLAMDP,KLAMDH,KEND00,LWRK00,KXIAJB,KFOCK0 1527 INTEGER IOPTTCME,IOPT 1528 INTEGER ISYMWZ,ISYMWU,ISYMWZU 1529 INTEGER ISYFCKZ,ISYFCKU,ISYMZU,KFCKZ,KFCKU,KEND1,LWRK1 1530 INTEGER ISINT1,ISINT2,KT3BOG2,KT3OG1 1531 INTEGER KT3OG2,KT3VIJG1,KXGADCK,KXLADCK 1532 INTEGER ISYMTETAZ,ISYMTETAU,ISYMTETAZU 1533 INTEGER ISYMK,ISYML,ISYMKL,ISYT30KL,ISTETAZKL,ISTETAUKL,ISTETAZUKL 1534 INTEGER KT30KL,KTETAZKL,KTETAUKL,KTETAZUKL,KTMAT 1535c 1536 INTEGER ISYMRX,ISYFCKX,KLAMDPX,KLAMDHX,KFOCKX,KT1X,ISYMC 1537 INTEGER KOFF1,KOFF2,ISINT2X,KT3BOG1,KT3BOL2,KINTOC,IOFF 1538c 1539 integer kx3am 1540c 1541#if defined (SYS_CRAY) 1542 REAL OMEGA2(*),OMEGA2EFF(*) 1543 REAL WORK(LWORK) 1544 REAL FREQRZ,FREQRU,FREQZU 1545 REAL XNORMVAL,DDOT 1546#else 1547 DOUBLE PRECISION OMEGA2(*),OMEGA2EFF(*) 1548 DOUBLE PRECISION WORK(LWORK) 1549 DOUBLE PRECISION FREQRZ,FREQRU,FREQZU 1550 DOUBLE PRECISION XNORMVAL,DDOT 1551#endif 1552C 1553 CALL QENTER('AT3ZUOCC') 1554C 1555C-------------------------- 1556C Save and set flags. 1557C-------------------------- 1558C 1559 IPRCC = IPRINT 1560 ISYM0 = 1 1561 ISINT1 = 1 1562 ISINT2 = 1 1563 ISYMT2 = 1 1564 ISYMT3 = MULD2H(ISINT2,ISYMT2) 1565C 1566C---------------------------------------------- 1567C Calculate the zeroth order stuff once 1568C---------------------------------------------- 1569C 1570 KT2TP = 1 1571 KFOCKD = KT2TP + NT2SQ(ISYMT2) 1572 KLAMDP = KFOCKD + NORBTS 1573 KLAMDH = KLAMDP + NLAMDT 1574 KEND00 = KLAMDH + NLAMDT 1575 LWRK00 = LWORK - KEND00 1576C 1577 KXIAJB = KEND00 1578 KEND00 = KXIAJB + NT2AM(ISINT1) 1579 LWRK00 = LWORK - KEND00 1580C 1581 IF (LWRK00 .LT. 0) THEN 1582 CALL QUIT('Out of memory in CC3_AMATT3ZUOCC (zeroth allo.') 1583 ENDIF 1584C 1585C------------------------ 1586C Construct L(ia,jb). 1587C------------------------ 1588C 1589 REWIND(LUIAJB) 1590 CALL READI(LUIAJB,IRAT*NT2AM(ISINT1),WORK(KXIAJB)) 1591 IOPTTCME = 1 1592 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISINT1,IOPTTCME) 1593C 1594 IF ( IPRINT .GT. 55) THEN 1595 XNORMVAL = DDOT(NT2AM(ISINT1),WORK(KXIAJB),1, 1596 * WORK(KXIAJB),1) 1597 WRITE(LUPRI,*) 'Norm of IAJB ',XNORMVAL 1598 ENDIF 1599C 1600C---------------------------------------------------------------- 1601C Read t1 and calculate the zero'th order Lambda matrices 1602C---------------------------------------------------------------- 1603C 1604 CALL GET_LAMBDA0(WORK(KLAMDP),WORK(KLAMDH),WORK(KEND00),LWRK00) 1605C 1606C------------------------------------------- 1607C Read in t2 , square it and reorder 1608C------------------------------------------- 1609C 1610 IOPT = 2 1611 CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KT2TP),'R0',0,ISYMT2, 1612 * WORK(KEND00),LWRK00) 1613 1614 IF (IPRINT .GT. 55) THEN 1615 XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 1616 WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL 1617 ENDIF 1618C 1619C-------------------------------------- 1620C Read in orbital energies 1621C-------------------------------------- 1622C 1623C 1624 CALL GET_ORBEN(WORK(KFOCKD),WORK(KEND00),LWRK00) 1625COMMENT COMMENT COMMENT 1626COMMENT COMMENT COMMENT If we want to sum the T3 amplitudes 1627COMMENT COMMENT COMMENT 1628 if (.false.) then 1629 kx3am = kend00 1630 kend00 = kx3am + nt1amx*nt1amx*nt1amx 1631 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 1632 lwrk00 = lwork - kend00 1633 if (lwrk00 .lt. 0) then 1634 write(lupri,*) 'Memory available : ',lwork 1635 write(lupri,*) 'Memory needed : ',kend00 1636 call quit('Insufficient space (T3) in CC3_AMATT3ZUOCC') 1637 END IF 1638 endif 1639COMMENT COMMENT COMMENT 1640COMMENT COMMENT COMMENT 1641COMMENT COMMENT COMMENT 1642C 1643C determining R1 labels 1644C 1645 ISYMRX = ISYLRT(IDLSTX) 1646 !FREQuency not needed for LISTX 1647 !LABELX = LRTLBL(IDLSTX) not needed for LISTX 1648C 1649 IF (LISTZU(1:3).EQ.'R2 ') THEN 1650 LISTRZ = 'R1 ' 1651 LABELRZ = LBLR2T(IDLSTZU,1) 1652 ISYMRZ = ISYR2T(IDLSTZU,1) 1653 FREQRZ = FRQR2T(IDLSTZU,1) 1654 LORXRZ = LORXR2T(IDLSTZU,1) 1655 IDLSTRZ = IR1TAMP(LABELRZ,LORXRZ,FREQRZ,ISYMRZ) 1656C 1657 IF (LORXRZ) CALL QUIT('NO ORBITAL RELAX. IN CC3_AMATT3ZUOCC') 1658 1659 LISTRU = 'R1 ' 1660 LABELRU = LBLR2T(IDLSTZU,2) 1661 ISYMRU = ISYR2T(IDLSTZU,2) 1662 FREQRU = FRQR2T(IDLSTZU,2) 1663 LORXRU = LORXR2T(IDLSTZU,2) 1664 IDLSTRU = IR1TAMP(LABELRU,LORXRU,FREQRU,ISYMRU) 1665C 1666 IF (LORXRU) CALL QUIT('NO ORBITAL RELAX. IN CC3_AMATT3ZUOCC') 1667C 1668 ELSE 1669 CALL QUIT('Unknown list for LISTZU in CC3_AMATT3ZUOCC.') 1670 END IF 1671C 1672 FREQZU = FREQRZ + FREQRU 1673C 1674 ISYMWZ = MULD2H(ISYMT3,ISYMRZ) 1675 ISYMWU = MULD2H(ISYMT3,ISYMRU) 1676C 1677 ISYMWZU = MULD2H(ISYMWZ,ISYMRU) 1678C 1679 IF (ISYRES.NE.ISYMWZU) THEN 1680 WRITE(LUPRI,*) 'INCONSISTENT SYMMETRY SPECIFICATION' 1681 WRITE(LUPRI,*) 1682 * 'ISYRES =',ISYRES,'ISYMWZU =',ISYMWZU 1683 CALL QUIT('CC3_AMATT3ZUOCC inconsistent symmetry specification') 1684 END IF 1685C 1686 ISYFCKZ = ISYMRZ 1687 ISYFCKU = ISYMRU 1688 ISYMZU = MULD2H(ISYMRZ,ISYMRU) 1689 ISYFCKX = MULD2H(ISYMOP,ISYMRX) 1690C 1691C------------------------------------------------- 1692C Read property matrices and Lambda matrices 1693C------------------------------------------------- 1694C 1695 KFCKZ = KEND00 1696 KFCKU = KFCKZ + N2BST(ISYFCKZ) 1697 KEND1 = KFCKU + N2BST(ISYFCKU) 1698 LWRK1 = LWORK - KEND1 1699C 1700 KLAMDPX = KEND1 1701 KLAMDHX = KLAMDPX + NLAMDT 1702 KEND1 = KLAMDHX + NLAMDT 1703 LWRK1 = LWORK - KEND1 1704C 1705 KFOCKX = KEND1 1706 KEND1 = KFOCKX + N2BST(ISYFCKX) 1707 LWRK1 = LWORK - KEND1 1708C 1709 KT1X = KEND1 1710 KEND1 = KT1X + NT1AM(ISYMRX) 1711 LWRK1 = LWORK - KEND1 1712C 1713 IF (LWRK1 .LT. NT2SQ(ISYMRZ)) THEN 1714 CALL QUIT('Out of memory in CC3_AMATT3ZUOCC (TOT_T3Y) ') 1715 ENDIF 1716C 1717C--------------------------- 1718C Read Lambda_X matrices 1719C--------------------------- 1720C 1721 CALL GET_LAMBDAX(WORK(KLAMDPX),WORK(KLAMDHX),LISTX,IDLSTX, 1722 * ISYMRX,WORK(KLAMDP),WORK(KLAMDH),WORK(KEND1), 1723 * LWRK1) 1724C 1725C------------------ 1726C Read in T1^X 1727C------------------ 1728C 1729 IOPT = 1 1730 CALL GET_T1_T2(IOPT,.FALSE.,WORK(KT1X),DUMMY,LISTX,IDLSTX, 1731 * ISYMRX,WORK(KEND1),LWRK1) 1732C 1733C------------------------------------------ 1734C Calculate the F^X matrix 1735C------------------------------------------ 1736C 1737 CALL CC3LR_MFOCK(WORK(KEND1),WORK(KT1X),WORK(KXIAJB),ISYFCKX) 1738C 1739C Put the F_{kc} part into correct F_{pq} 1740C 1741 CALL DZERO(WORK(KFOCKX),N2BST(ISYFCKX)) 1742C 1743 DO ISYMC = 1, NSYM 1744 ISYMK = MULD2H(ISYFCKX,ISYMC) 1745 DO K = 1, NRHF(ISYMK) 1746 DO C = 1, NVIR(ISYMC) 1747 KOFF1 = KFOCKX - 1 1748 * + IFCVIR(ISYMK,ISYMC) 1749 * + NORB(ISYMK)*(C - 1) 1750 * + K 1751 KOFF2 = KEND1 - 1 1752 * + IT1AM(ISYMC,ISYMK) 1753 * + NVIR(ISYMC)*(K - 1) 1754 * + C 1755C 1756 WORK(KOFF1) = WORK(KOFF2) 1757C 1758 ENDDO 1759 ENDDO 1760 ENDDO 1761 1762C 1763C------------------------------------------ 1764C Calculate the F^Z matrix 1765C------------------------------------------ 1766C 1767 CALL GET_FOCKX(WORK(KFCKZ),LABELRZ,IDLSTRZ,ISYMRZ,WORK(KLAMDP), 1768 * ISYM0,WORK(KLAMDH),ISYM0,WORK(KEND1),LWRK1) 1769C 1770C------------------------------------------ 1771C Calculate the F^U matrix 1772C------------------------------------------ 1773C 1774 CALL GET_FOCKX(WORK(KFCKU),LABELRU,IDLSTRU,ISYMRU,WORK(KLAMDP), 1775 * ISYM0,WORK(KLAMDH),ISYM0,WORK(KEND1),LWRK1) 1776C 1777C----------------------------- 1778C Read integrals. 1779C----------------------------- 1780C 1781C 1782C----------------------------- 1783C Memory allocation. 1784C----------------------------- 1785C 1786C isint1, isint2 - symmetry of integrals in standard H, transformed 1787C with LambdaH_0 1788C isintr1 - symmetry of integrals in standard H, transformed 1789C with LambdaH_R1 1790 1791 ISINT1 = 1 1792 ISINT2 = 1 1793 ISINT2X = MULD2H(ISINT2,ISYMRX) 1794C 1795 KT3BOG1 = KEND1 1796 KT3BOG2 = KT3BOG1 + NTRAOC(ISINT2X) 1797 KT3BOL2 = KT3BOG2 + NTRAOC(ISINT2X) 1798 KT3OG1 = KT3BOL2 + NTRAOC(ISINT2X) 1799 KT3OG2 = KT3OG1 + NTRAOC(ISINT2) 1800 KEND1 = KT3OG2 + NTRAOC(ISINT2) 1801C 1802 KT3VIJG1 = KEND1 1803 KEND1 = KT3VIJG1 + NMAABCI(ISYM0) 1804 LWRK1 = LWORK - KEND1 1805C 1806 KXGADCK = KEND1 1807 KXLADCK = KXGADCK + NMAABCI(ISYMRX) 1808 KEND1 = KXLADCK + NMAABCI(ISYMRX) 1809 LWRK1 = LWORK - KEND1 1810C 1811 KINTOC = KEND1 1812 KEND1 = KINTOC + NTOTOC(ISINT2) 1813 LWRK1 = LWORK - KEND1 1814C 1815 IF (LWRK1 .LT. 0) THEN 1816 WRITE(LUPRI,*) 'Memory available : ',LWORK 1817 WRITE(LUPRI,*) 'Memory needed : ',KEND1 1818 CALL QUIT('Insufficient space in CC3_AMATT3ZUOCC (3)') 1819 END IF 1820C 1821C----------------------------------------------------------------- 1822C Construct occupied integrals which are required for occupied 1823C part of contraction (KT3BOG2) 1824C----------------------------------------------------------------- 1825C 1826c CALL INTOCC_T3BAR0(LUTOC,FNTOC,WORK(KLAMDH),ISYM0,WORK(KT3BOG1), 1827c * WORK(KT3BOL1),WORK(KT3BOG2),WORK(KT3BOL2), 1828c * WORK(KEND1),LWRK1) 1829 IOFF = 1 1830 IF (NTOTOC(ISINT2) .GT. 0) THEN 1831 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 1832 ENDIF 1833C 1834 !Transform (ia|j delta) integrals to (ia|j k-) and sort G(j,i,k-,a) 1835 CALL CCLR_TROCC(WORK(KINTOC),WORK(KT3BOG1),WORK(KLAMDHX),ISYMRX, 1836 * WORK(KEND1),LWRK1) 1837C 1838 !integrals g(ia|j k-) sorted as G(j,i,k-,a) resorted as G(k-,i,j,a) 1839 !KT3BOL2 are not used 1840 CALL CC3_SRTOCC(WORK(KT3BOG1),WORK(KT3BOL2),ISINT2X) 1841C 1842 ! (ia|jk) sorted as T3BOG1(kija) and now sorted as T3BOG2(ajik) 1843 CALL CCFOP_SORT(WORK(KT3BOG1),WORK(KT3BOG2),ISINT2X,1) 1844C 1845C----------------------------------------------------------------- 1846C Construct occupied integrals which are required to calculate 1847C t3_0 amplitudes 1848C----------------------------------------------------------------- 1849C 1850 CALL INTOCC_T30(LUCKJD,FNCKJD,WORK(KLAMDP),ISINT2,WORK(KT3OG1), 1851 * WORK(KT3OG2),WORK(KEND1),LWRK1) 1852C 1853C---------------------------------------------- 1854C Get virtual integrals for t30 amplitudes 1855C KT3VIJG1 : (ck|da) sorted as I(ad|ck) 1856C---------------------------------------------- 1857C 1858 CALL INTVIR_T30_IJ(WORK(KT3VIJG1),ISYM0,WORK(KLAMDH),LUDELD, 1859 * FNDELD,WORK(KEND1),LWRK1) 1860C 1861C------------------------------------------------------------------- 1862C Get virtual integrals for virtual part of contraction (KXGADCK) 1863C KXGADCK g(kcad) = (kc ! ad) sorted as I(adck) 1864C KXLADCK L(kcad) sorted as I(adck) 1865C------------------------------------------------------------------- 1866C 1867 CALL INTVIR_T3B0_JK(WORK(KXGADCK),WORK(KXLADCK),ISYM0, 1868 * WORK(KLAMDPX),ISYMRX,LU3VI,FN3VI, 1869 * LU3FOP,FN3FOP, 1870 * WORK(KEND1),LWRK1) 1871C 1872C---------------------------- 1873C Loop over K 1874C---------------------------- 1875C 1876 ISYMTETAZ = MULD2H(ISYM0,ISYMRZ) 1877 ISYMTETAU = MULD2H(ISYM0,ISYMRU) 1878 ISYMTETAZU = MULD2H(ISYM0,ISYMZU) 1879 DO ISYMK = 1,NSYM 1880 1881 DO K = 1,NRHF(ISYMK) 1882C 1883C---------------------------- 1884C Loop over L 1885C---------------------------- 1886C 1887 DO ISYML = 1,NSYM 1888C 1889 ISYMKL = MULD2H(ISYMK,ISYML) 1890 ISYT30KL = MULD2H(ISYMKL,ISYMT3) 1891 ISTETAZKL = MULD2H(ISYMKL,ISYMTETAZ) 1892 ISTETAUKL = MULD2H(ISYMKL,ISYMTETAU) 1893 ISTETAZUKL = MULD2H(ISYMKL,ISYMTETAZU) 1894C 1895 KT30KL = KEND1 1896 KTETAZKL = KT30KL + NMAABCI(ISYT30KL) 1897 KTETAUKL = KTETAZKL + NMAABCI(ISTETAZKL) 1898 KTETAZUKL = KTETAUKL + NMAABCI(ISTETAUKL) 1899 KEND1 = KTETAZUKL + NMAABCI(ISTETAZUKL) 1900 LWRK1 = LWORK - KEND1 1901C 1902 KTMAT = KEND1 1903 KEND1 = KTMAT + NMAABCI(ISTETAZUKL) 1904 LWRK1 = LWORK - KEND1 1905C 1906 IF (LWRK1 .LT. 0) THEN 1907 WRITE(LUPRI,*) 'Memory available : ',LWORK 1908 WRITE(LUPRI,*) 'Memory needed : ',KEND1 1909 CALL QUIT('Insufficient space in CC3_AMATT3ZUOCC (4)') 1910 END IF 1911C 1912 DO L = 1,NRHF(ISYML) 1913C 1914C 1915C------------------------------------------- 1916C Get T30^KL amplitudes 1917C------------------------------------------- 1918C 1919 CALL DZERO(WORK(KT30KL),NMAABCI(ISYT30KL)) 1920C 1921 CALL GET_T30_IJ_O(WORK(KT30KL),ISYT30KL,WORK(KT2TP), 1922 * ISYM0, 1923 * WORK(KT3OG2),ISYM0,ISYML,L,ISYMK,K, 1924 * WORK(KEND1),LWRK1) 1925C 1926 CALL GET_T30_IJ_V(WORK(KT30KL),ISYT30KL,WORK(KT2TP), 1927 * ISYM0,WORK(KT3VIJG1), 1928 * ISYM0,ISYML,L,ISYMK,K, 1929 * WORK(KEND1),LWRK1) 1930C 1931C-------------------------------------------------------------- 1932C Divide by orbital energy difference and remove 1933C forbidden elements 1934C-------------------------------------------------------------- 1935C 1936 CALL T3JK_DIA(WORK(KT30KL),ISYT30KL,ISYML,L,ISYMK,K, 1937 * WORK(KFOCKD)) 1938 CALL T3_FORBIDDEN_JK(WORK(KT30KL),ISYMT3,ISYML,L, 1939 * ISYMK,K) 1940C 1941c call sum_pt3_jk(work(kt30kl),isyml,l,isymk,k,isyt30kl, 1942c * work(kx3am),3) 1943C 1944 IF (IPRINT .GT. 55) THEN 1945 WRITE(LUPRI,*)'ISYML,L,ISYMK,K ', ISYML,L,ISYMK,K 1946 XNORMVAL = DDOT(NMAABCI(ISYT30KL),WORK(KT30KL),1, 1947 * WORK(KT30KL),1) 1948 WRITE(LUPRI,*)'NORM OF KT30KL IN CC3_AMATT3ZUOCC ', 1949 * XNORMVAL 1950 END IF 1951C 1952C--------------------------------- 1953C Calculate KTETAZKL 1954C--------------------------------- 1955C 1956 CALL DZERO(WORK(KTETAZKL),NMAABCI(ISTETAZKL)) 1957C 1958 CALL TETAX_JK_BC(WORK(KT30KL),ISYT30KL,WORK(KFCKZ), 1959 * ISYFCKZ,WORK(KTETAZKL),ISTETAZKL, 1960 * WORK(KEND1),LWRK1) 1961C 1962c call w3jk_dia(work(KTETAZKL),ISTETAZKL,isyml,l,isymk, 1963c * k,work(kfockd),freqr1) 1964c call t3_forbidden_jk(work(KTETAZKL),ISYMTETAZ,isyml,l, 1965c * isymk,k) 1966c 1967c call sum_pt3_jk(work(KTETAZKL),isyml,l,isymk,k, 1968c * ISTETAZKL, 1969c * work(kx3am),1) 1970c 1971c xnormval = ddot(nmaabci(ISTETAZKL),work(KTETAZKL),1, 1972c * work(KTETAZKL),1) 1973c write(lupri,*)'norm of KTETAZKL in CC3_AMATT3ZUOCC ', 1974c * xnormval 1975C 1976 CALL TETAX_JK_A(WORK(KT30KL),ISYT30KL,WORK(KFCKZ), 1977 * ISYFCKZ,WORK(KTETAZKL),ISTETAZKL, 1978 * WORK(KEND1),LWRK1) 1979C 1980 CALL W3JK_DIA(WORK(KTETAZKL),ISTETAZKL,ISYML,L,ISYMK, 1981 * K,WORK(KFOCKD),FREQRZ) 1982 CALL T3_FORBIDDEN_JK(WORK(KTETAZKL),ISYMTETAZ,ISYML,L, 1983 * ISYMK,K) 1984 1985c call sum_pt3_jk(work(KTETAZKL),isyml,l,isymk,k, 1986c * ISTETAZKL, 1987c * work(kx3am),4) 1988 1989C 1990 IF (IPRINT .GT. 55) THEN 1991 XNORMVAL = DDOT(NMAABCI(ISTETAZKL),WORK(KTETAZKL),1, 1992 * WORK(KTETAZKL),1) 1993 WRITE(LUPRI,*)'(FINAL) NORM OF KTETAZKL ', XNORMVAL 1994 END IF 1995C 1996C--------------------------------- 1997C Calculate KTETAUKL 1998C--------------------------------- 1999C 2000 CALL DZERO(WORK(KTETAUKL),NMAABCI(ISTETAUKL)) 2001C 2002 CALL TETAX_JK_BC(WORK(KT30KL),ISYT30KL,WORK(KFCKU), 2003 * ISYFCKU,WORK(KTETAUKL),ISTETAUKL, 2004 * WORK(KEND1),LWRK1) 2005C 2006C 2007 CALL TETAX_JK_A(WORK(KT30KL),ISYT30KL,WORK(KFCKU), 2008 * ISYFCKU,WORK(KTETAUKL),ISTETAUKL, 2009 * WORK(KEND1),LWRK1) 2010C 2011 CALL W3JK_DIA(WORK(KTETAUKL),ISTETAUKL,ISYML,L,ISYMK, 2012 * K,WORK(KFOCKD),FREQRU) 2013 CALL T3_FORBIDDEN_JK(WORK(KTETAUKL),ISYMTETAU,ISYML,L, 2014 * ISYMK,K) 2015 2016c call sum_pt3_jk(work(KTETAUKL),isyml,l,isymk,k, 2017c * ISTETAUKL, 2018c * work(kx3am),1) 2019 2020C 2021 IF (IPRINT .GT. 55) THEN 2022 XNORMVAL = DDOT(NMAABCI(ISTETAUKL),WORK(KTETAUKL),1, 2023 * WORK(KTETAUKL),1) 2024 WRITE(LUPRI,*)'(FINAL) NORM OF KTETAUKL ', XNORMVAL 2025C 2026 END IF 2027C 2028C------------------------------------------------------------- 2029C Calculate KTETAZUKL = <mu3|[U,KTETAZKL]|HF> 2030C------------------------------------------------------------- 2031C 2032 CALL DZERO(WORK(KTETAZUKL),NMAABCI(ISTETAZUKL)) 2033C 2034 ! virtual contribution: 2035 ! KTETAZUKL = sum_d U_ad KTETAZKL^{d-,b-,c-}_{iKL} 2036 CALL TETAX_JK_A(WORK(KTETAZKL),ISTETAZKL,WORK(KFCKU), 2037 * ISYFCKU,WORK(KTETAZUKL),ISTETAZUKL, 2038 * WORK(KEND1),LWRK1) 2039C 2040 ! occupied contribution: 2041 ! KTETAZUKL = sum_m U_mi KTETAZKL^{a-,b-,c-}_{mKL} 2042 CALL TETAX_JK_I(WORK(KTETAZKL),ISTETAZKL,WORK(KFCKU), 2043 * ISYFCKU,WORK(KTETAZUKL),ISTETAZUKL, 2044 * WORK(KEND1),LWRK1) 2045C 2046C------------------------------------------------------------- 2047C Calculate KTETAZUKL = <mu3|[Z,KTETAUKL]|HF> 2048C------------------------------------------------------------- 2049C 2050 ! virtual contribution: 2051 ! KTETAZUKL = sum_d Z_ad KTETAUKL^{d-,b-,c-}_{iKL} 2052 CALL TETAX_JK_A(WORK(KTETAUKL),ISTETAUKL,WORK(KFCKZ), 2053 * ISYFCKZ,WORK(KTETAZUKL),ISTETAZUKL, 2054 * WORK(KEND1),LWRK1) 2055C 2056 ! occupied contribution: 2057 ! KTETAZUKL = sum_m Z_mi KTETAUKL^{a-,b-,c-}_{mKL} 2058 CALL TETAX_JK_I(WORK(KTETAUKL),ISTETAUKL,WORK(KFCKZ), 2059 * ISYFCKZ,WORK(KTETAZUKL),ISTETAZUKL, 2060 * WORK(KEND1),LWRK1) 2061C 2062C--------------------------------------------------------------- 2063C Divide by orbital energy difference and remove 2064C forbidden elements 2065C--------------------------------------------------------------- 2066C 2067 CALL T3_FORBIDDEN_JK(WORK(KTETAZUKL),ISYMTETAZU,ISYML, 2068 * L,ISYMK,K) 2069 2070c call sum_pt3_jk(work(KTETAZUKL),isyml,l,isymk,k, 2071c * ISTETAZUKL, 2072c * work(kx3am),4) 2073 2074 CALL W3JK_DIA(WORK(KTETAZUKL),ISTETAZUKL,ISYML,L, 2075 * ISYMK,K,WORK(KFOCKD),FREQZU) 2076 2077 2078C 2079 IF (IPRINT .GT. 55) THEN 2080 XNORMVAL = DDOT(NMAABCI(ISTETAZUKL),WORK(KTETAZUKL), 2081 * 1,WORK(KTETAZUKL),1) 2082 WRITE(LUPRI,*)'(FINAL) NORM OF KTETAZUKL ', XNORMVAL 2083C 2084 END IF 2085C 2086C------------------------------------------------------- 2087C Project up to doubles space: 2088C omega2eff = - <mu2|[H,THETAZU^LK]|HF> 2089C------------------------------------------------------- 2090C 2091 ! Fock matrix contribution 2092 CALL CC3_CY2F_JK(OMEGA2EFF,ISYRES,WORK(KTETAZUKL), 2093 * ISTETAZUKL,WORK(KTMAT), 2094 * WORK(KFOCKX),ISYFCKX, 2095 * WORK(KEND1),LWRK1, 2096 * ISYML,L,ISYMK,K) 2097C 2098 ! occupied integrals contribution 2099 CALL CC3_CY2O_JK(OMEGA2EFF,ISYRES,WORK(KTETAZUKL), 2100 * ISTETAZUKL,WORK(KTMAT),WORK(KT3BOG2), 2101 * ISYM0,WORK(KEND1),LWRK1, 2102 * ISYML,L,ISYMK,K) 2103 ! virtual integrals contribution 2104 CALL CC3_CY2V_JK(OMEGA2EFF,ISYRES,WORK(KTETAZUKL), 2105 * ISTETAZUKL,WORK(KTMAT),WORK(KXGADCK), 2106 * ISYM0,WORK(KEND1),LWRK1, 2107 * ISYML,L,ISYMK,K) 2108C 2109 ENDDO ! L 2110 ENDDO ! ISYML 2111 ENDDO ! K 2112 ENDDO ! ISYMK 2113C 2114c write(lupri,*) 'T30KL in CC3_AMATT3ZUOCC' 2115c call print_pt3(work(kx3am),isym0,4) 2116C 2117C 2118C------------------------------------------------------ 2119C add: omega contributions to omegaeff 2120C------------------------------------------------------ 2121C 2122c write(lupri,*)'OMEGA2EFF in CC3_AASDOCC' 2123c call outpak(OMEGA2EFF,NT1AM(isyres),1,lupri) 2124 2125 IF (IPRINT .GT. 55) THEN 2126 XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 2127 WRITE(LUPRI,*)'Norm of OMEGA2EFF final ',XNORMVAL 2128 ENDIF 2129C 2130C------------- 2131C End 2132C------------- 2133C 2134 CALL QEXIT('AT3ZUOCC') 2135C 2136 RETURN 2137 END 2138C /* Deck cc3_bmatt3zu */ 2139 SUBROUTINE CC3_BMATT3ZU(LISTX,IDLSTX,LISTZU,IDLSTZU, 2140 * OMEGA2,OMEGA2EFF, 2141 * ISYRES,WORK,LWORK) 2142* 2143******************************************************************* 2144* 2145* This routine calculates the contribution to cubic response 2146* (originating form both F matrix and B matrix): 2147* 2148* omega2eff = <mu2|[H^X,T3^ZU]|HF> 2149* 2150* which is contracted outside with the multipliers (zeroth-order for 2151* F matrix or first-order for B matrix). 2152* 2153* 2154* 2155* T3^ZU is calculated using W^BD intermmediate. 2156* 2157* !!! ACTUALLY HERE WE CALCULATE ONLY THE PART OF THE TERM, WHICH ORIGINATES 2158* FROM THE B-MATRIX CONTRIBUTIONS TO T3^ZU + 2 EXTRA CONTRIBUTIONS FROM 2159* JACOBIAN (to get the rest of contribution call cc3_amatt3zu routine) !!! 2160* 2161* Thus here: 2162* 2163* W^BD = 1/2<mu3|[[[H,T1^Z],T1^U],T2^0]|HF> + <mu3|[H^Z,T2^U]|HF> !b-mat 2164* + <mu3|[[H,T1^ZU],T2^0]|HF> + <mu3|[H,T2^ZU]|HF> !jacobian 2165* 2166* Note (from the expression for T3^ZU) that the permutation operator P(Z,U) 2167* does NOT act on the last two terms. 2168* 2169* 2170* F. Pawlowski, P. Jorgensen, Aarhus Spring 2003. 2171* 2172******************************************************************* 2173C 2174 IMPLICIT NONE 2175C 2176#include "priunit.h" 2177#include "ccr1rsp.h" 2178#include "ccinftap.h" 2179#include "iratdef.h" 2180#include "ccsdsym.h" 2181#include "ccsdinp.h" 2182#include "dummy.h" 2183#include "ccorb.h" 2184C 2185 CHARACTER*6 FN3VI 2186 CHARACTER*8 FNTOC,FN3VI2 2187 CHARACTER*6 FNCKJD 2188 CHARACTER*4 FNDKBC 2189C 2190 PARAMETER (FNTOC = 'CCSDT_OC' ) 2191 PARAMETER (FN3VI = 'CC3_VI' , FN3VI2 = 'CC3_VI12') 2192 PARAMETER (FNCKJD = 'CKJDEL' ) 2193 PARAMETER (FNDKBC = 'DKBC' ) 2194 2195 CHARACTER*12 FN3SRTR, FNCKJDR, FNDELDR, FNDKBCR 2196 CHARACTER*13 FNCKJDR1, FNDKBCR1 2197 CHARACTER*13 FNCKJDR2, FNDKBCR2 2198 CHARACTER*13 FNCKJDR3, FNDKBCR3 2199C 2200 PARAMETER(FN3SRTR = 'CCSDT_FBMAT1', FNCKJDR = 'CCSDT_FBMAT2', 2201 * FNDELDR = 'CCSDT_FBMAT3', FNDKBCR = 'CCSDT_FBMAT4') 2202C 2203 PARAMETER(FNCKJDR1 = 'CCSDT_FBMAT21',FNDKBCR1 = 'CCSDT_FBMAT41') 2204 PARAMETER(FNCKJDR2 = 'CCSDT_FBMAT22',FNDKBCR2 = 'CCSDT_FBMAT42') 2205 PARAMETER(FNCKJDR3 = 'CCSDT_FBMAT23',FNDKBCR3 = 'CCSDT_FBMAT43') 2206C 2207 INTEGER LUTOC,LU3VI,LU3VI2,LUCKJD,LUDKBC 2208 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR 2209 INTEGER LUCKJDR1, LUDKBCR1 2210 INTEGER LUCKJDR2, LUDKBCR2 2211 INTEGER LUCKJDR3, LUDKBCR3 2212C 2213 CHARACTER*3 LISTZ, LISTU, LISTZU, LISTX 2214 CHARACTER*8 LABELZ, LABELU, LABELX, LABELZU 2215 INTEGER IDLSTZ,IDLSTU,IDLSTZU,IDLSTX 2216 INTEGER ISYMRZ,ISYMRU 2217C 2218 CHARACTER CDUMMY*1 2219C 2220 INTEGER ISYRES 2221 INTEGER LWORK 2222C 2223 INTEGER ISYM0,ISINT1,ISINT2,ISYMT2 2224 INTEGER KT2TP,KFOCKD,KLAMDP,KLAMDH,KXIAJB,KEND00,LWRK00 2225 INTEGER IOPTTCME,IOPT 2226 INTEGER KT2U,KT1Z,KEND1,LWRK1 2227 INTEGER ISINTR1Z,ISINTR2Z 2228 INTEGER ISINTZU 2229 INTEGER KRBJIA,KTROC,KTROC1,KTROC0Z,KTROC3,KINTOC 2230 INTEGER MAXOCC,KEND2,LWRK2 2231 INTEGER IOFF 2232 INTEGER ISYMD,ISCKB1,ISCKB2Z,ISCKBZU 2233 INTEGER KTRVI,KTRVI1,KTRVI0Z,KEND3,LWRK3,KTRVI7,KINTVI 2234 INTEGER ISYMB,ISYMZU,ISYMBD,ISWMAT,ISCKD2Z,ISCKDZU 2235 INTEGER KDIAGW,KWMAT,KTMAT,KEND4,KTRVI8Z,LWRK4,KTRVI9 2236 INTEGER LENSQW,KINDSQW 2237 INTEGER KEND0,LWRK0 2238 INTEGER ISINTR1U,ISINTR2U 2239 INTEGER KT2Z,KT1U,KTROC0U 2240 INTEGER ISCKB2U,KTRVI0U 2241 INTEGER ISCKD2U,KTRVI8U 2242 INTEGER MMAXOCC 2243c 2244 INTEGER ISYMRX,ISYMRZU,ISYFCKX,KT2ZU,KLAMDPX,KLAMDHX 2245 INTEGER KFOCKX,KT1ZU,KT1X,ISYMC,ISYMK,KOFF1,KOFF2,ISINTR1ZU 2246 INTEGER ISINTR2ZU,ISINT1X,KTROC0ZU,KT3OG1,MMMAXOCC,ISCKB1X 2247 INTEGER ISCKB2ZU,KTRVI0ZU,KT3VDG1,ISCKD2ZU,ISCKD1,KTRVI8ZU 2248 INTEGER KT3VBG1 2249c 2250 integer kx3am 2251c 2252C 2253#if defined (SYS_CRAY) 2254 REAL OMEGA2(*),OMEGA2EFF(*) 2255 REAL WORK(LWORK) 2256 REAL FREQZ,FREQU,FREQZU 2257 REAL XNORMVAL,DDOT 2258#else 2259 DOUBLE PRECISION OMEGA2(*),OMEGA2EFF(*) 2260 DOUBLE PRECISION WORK(LWORK) 2261 DOUBLE PRECISION FREQZ,FREQU,FREQZU 2262 DOUBLE PRECISION XNORMVAL,DDOT 2263#endif 2264C 2265 2266 CALL QENTER('CC3_BMATT3ZU') 2267C 2268C----------------------------- 2269C Lists handling and check 2270C----------------------------- 2271C 2272 ISYMRX = ISYLRT(IDLSTX) 2273 !FREQuency not needed for LISTX 2274 !LABELX = LRTLBL(IDLSTX) not needed for LISTX 2275C 2276C 2277 IF (LISTZU(1:3).EQ.'R2 ') THEN 2278 LISTZ = 'R1 ' 2279 LABELZ = LBLR2T(IDLSTZU,1) 2280 ISYMRZ = ISYR2T(IDLSTZU,1) 2281 FREQZ = FRQR2T(IDLSTZU,1) 2282 LORXRZ = LORXR2T(IDLSTZU,1) 2283 IDLSTZ = IR1TAMP(LABELRZ,LORXRZ,FREQRZ,ISYMRZ) 2284C 2285 IF (LORXRZ) CALL QUIT('NO ORBITAL RELAX. IN CC3_BMATT3ZU') 2286 2287 LISTU = 'R1 ' 2288 LABELU = LBLR2T(IDLSTZU,2) 2289 ISYMRU = ISYR2T(IDLSTZU,2) 2290 FREQU = FRQR2T(IDLSTZU,2) 2291 LORXRU = LORXR2T(IDLSTZU,2) 2292 IDLSTU = IR1TAMP(LABELRU,LORXRU,FREQRU,ISYMRU) 2293C 2294 IF (LORXRU) CALL QUIT('NO ORBITAL RELAX. IN CC3_BMATT3ZU') 2295C 2296 ELSE 2297 CALL QUIT('Unknown list for LISTZU in CC3_BMATT3ZU.') 2298 END IF 2299C 2300 ISYMZU = MULD2H(ISYMRZ,ISYMRU) 2301C 2302 IF (LISTZU(1:3).EQ.'R2 ') THEN 2303 ISYMRZU = ISYLRT(IDLSTZU) 2304 !FREQuency for second-order amplitudes not needed - we get it 2305 ! below as the sum of "first-order" frequencies 2306 !LABELZU = LRTLBL(IDLSTZU) not needed for LISTZU 2307 ELSE 2308 CALL QUIT('Unknown list for LISTZU in CC3_BMATT3ZU') 2309 END IF 2310C 2311 !Symmetry check 2312 IF (ISYMZU .NE. ISYMRZU) THEN 2313 WRITE(LUPRI,*)'ISYMZU = ', ISYMZU 2314 WRITE(LUPRI,*)'ISYMRZU = ', ISYMRZU 2315 WRITE(LUPRI,*)'ISYMZU should be the same as ISYMRZU !!! ' 2316 CALL QUIT('Symmetry mismatch in CC3_BMATT3ZU') 2317 END IF 2318C 2319C------------------------------------------- 2320C Construct FREQZU = (omega_Z + omega_U) 2321C------------------------------------------- 2322C 2323 FREQZU = FREQZ + FREQU 2324C 2325C-------------------------- 2326C Save and set flags. 2327C-------------------------- 2328C 2329 ISYM0 = 1 2330 ISINT1 = 1 2331 ISINT2 = 1 2332 ISYMT2 = 1 2333C 2334C------------------------------------------ 2335C Open files (integrals in contraction) 2336C------------------------------------------ 2337C 2338 LUTOC = -1 2339 LU3VI = -1 2340 LU3VI2 = -1 2341C 2342 CALL WOPEN2(LUTOC,FNTOC,64,0) 2343 CALL WOPEN2(LU3VI,FN3VI,64,0) 2344 CALL WOPEN2(LU3VI2,FN3VI2,64,0) 2345C 2346C------------------------------------------ 2347C Open files (H0 integrals for WBD) 2348C------------------------------------------ 2349C 2350 LUCKJD = -1 2351 LUDKBC = -1 2352C 2353 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 2354 CALL WOPEN2(LUDKBC,FNDKBC,64,0) 2355C 2356C----------------------------------------- 2357C Open temporary files (H^Z integrals) 2358C----------------------------------------- 2359C 2360 LU3SRTR = -1 2361 LUCKJDR = -1 2362 LUDELDR = -1 2363 LUDKBCR = -1 2364C 2365 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 2366 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 2367 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 2368 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 2369C 2370C------------------------------------------ 2371C Open temporary files (H^{Z,U} integrals) !double one-index transformed 2372C !with first-order amplitudes 2373C------------------------------------------ 2374C 2375 LUCKJDR1 = -1 2376 LUDKBCR1 = -1 2377C 2378 CALL WOPEN2(LUCKJDR1,FNCKJDR1,64,0) 2379 CALL WOPEN2(LUDKBCR1,FNDKBCR1,64,0) 2380C 2381C------------------------------------------ 2382C Open temporary files (H^U integrals) 2383C------------------------------------------ 2384C 2385 LUCKJDR2 = -1 2386 LUDKBCR2 = -1 2387C 2388 CALL WOPEN2(LUCKJDR2,FNCKJDR2,64,0) 2389 CALL WOPEN2(LUDKBCR2,FNDKBCR2,64,0) 2390C 2391C------------------------------------------ 2392C Open temporary files (H^ZU integrals) !one-index transformed 2393C !with second-order amplitudes 2394C------------------------------------------ 2395C 2396 LUCKJDR3 = -1 2397 LUDKBCR3 = -1 2398C 2399 CALL WOPEN2(LUCKJDR3,FNCKJDR3,64,0) 2400 CALL WOPEN2(LUDKBCR3,FNDKBCR3,64,0) 2401C 2402C---------------------------------------------- 2403C Calculate the zeroth order stuff once 2404C---------------------------------------------- 2405C 2406 KT2TP = 1 2407 KFOCKD = KT2TP + NT2SQ(ISYMT2) 2408 KLAMDP = KFOCKD + NORBTS 2409 KLAMDH = KLAMDP + NLAMDT 2410 KXIAJB = KLAMDH + NLAMDT 2411 KEND00 = KXIAJB + NT2AM(ISINT1) 2412 LWRK00 = LWORK - KEND00 2413C 2414 IF (LWRK00 .LT. 0) THEN 2415 WRITE(LUPRI,*)'Memory available: ',LWORK 2416 WRITE(LUPRI,*)'Memory needed: ',KEND00 2417 CALL QUIT('Out of memory in CC3_BMATT3ZU (zeroth allo.)') 2418 ENDIF 2419C 2420C------------------------ 2421C Construct L(ia,jb). 2422C------------------------ 2423C 2424 REWIND(LUIAJB) 2425 CALL READI(LUIAJB,IRAT*NT2AM(ISINT1),WORK(KXIAJB)) 2426 IOPTTCME = 1 2427 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISINT1,IOPTTCME) 2428C 2429 IF ( IPRINT .GT. 55) THEN 2430 XNORMVAL = DDOT(NT2AM(ISINT1),WORK(KXIAJB),1, 2431 * WORK(KXIAJB),1) 2432 WRITE(LUPRI,*) 'Norm of IAJB ',XNORMVAL 2433 ENDIF 2434C 2435C---------------------------------------------------------------- 2436C Read t1 and calculate the zero'th order Lambda matrices 2437C---------------------------------------------------------------- 2438C 2439 CALL GET_LAMBDA0(WORK(KLAMDP),WORK(KLAMDH),WORK(KEND00),LWRK00) 2440C 2441C------------------------------------------- 2442C Read in t2 , square it and reorder 2443C------------------------------------------- 2444C 2445 IOPT = 2 2446 CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KT2TP),'R0',0,ISYMT2, 2447 * WORK(KEND00),LWRK00) 2448 2449 IF (IPRINT .GT. 55) THEN 2450 XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 2451 WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL 2452 ENDIF 2453C 2454C-------------------------------------- 2455C Read in orbital energies 2456C-------------------------------------- 2457C 2458 CALL GET_ORBEN(WORK(KFOCKD),WORK(KEND00),LWRK00) 2459c 2460c If we want to sum the T3 amplitudes 2461c 2462 if (.false.) then 2463 kx3am = kend00 2464 kend00 = kx3am + nt1amx*nt1amx*nt1amx 2465 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 2466 lwrk00 = lwork - kend00 2467 if (lwrk00 .lt. 0) then 2468 write(lupri,*) 'Memory available : ',lwork 2469 write(lupri,*) 'Memory needed : ',kend00 2470 call quit('Insufficient space (T3) in CC3_BMATT3ZU (1)') 2471 END IF 2472 endif 2473C 2474C------------------------------------------------- 2475C Read T1^Z and T2^Z 2476C Read T1^U and T2^U 2477C Read T1^ZU and T2^ZU !second-order amplitudes 2478C 2479C Calculate (ck|de)-Utilde and (ck|lm)-Utilde 2480C Calculate (ck|de)-{Z,U}tilde and (ck|lm)-{Z,U}tilde 2481C Calculate (ck|de)-Ztilde and (ck|lm)-Ztilde 2482C 2483C Used to construct WBD intermmediate 2484C------------------------------------------------- 2485C 2486 ISYFCKX = MULD2H(ISYMOP,ISYMRX) 2487C 2488 KT2U = KEND00 2489 KEND0 = KT2U + NT2SQ(ISYMRU) 2490 LWRK0 = LWORK - KEND0 2491C 2492 KT2Z = KEND0 2493 KEND0 = KT2Z + NT2SQ(ISYMRZ) 2494 LWRK0 = LWORK - KEND0 2495C 2496 KT2ZU = KEND0 2497 KEND0 = KT2ZU + NT2SQ(ISYMRZU) 2498 LWRK0 = LWORK - KEND0 2499C 2500 KLAMDPX = KEND0 2501 KLAMDHX = KLAMDPX + NLAMDT 2502 KEND0 = KLAMDHX + NLAMDT 2503 LWRK0 = LWORK - KEND0 2504C 2505 KFOCKX = KEND0 2506 KEND0 = KFOCKX + N2BST(ISYFCKX) 2507 LWRK0 = LWORK - KEND0 2508C 2509 KT1U = KEND0 2510 KEND1 = KT1U + NT1AM(ISYMRU) 2511 LWRK1 = LWORK - KEND1 2512C 2513 KT1Z = KEND1 2514 KEND1 = KT1Z + NT1AM(ISYMRZ) 2515 LWRK1 = LWORK - KEND1 2516C 2517 KT1ZU = KEND1 2518 KEND1 = KT1ZU + NT1AM(ISYMRZU) 2519 LWRK1 = LWORK - KEND1 2520C 2521 KT1X = KEND1 2522 KEND1 = KT1X + NT1AM(ISYMRX) 2523 LWRK1 = LWORK - KEND1 2524C 2525 IF (LWRK1 .LT. NT2SQ(ISYMRZ)) THEN 2526 CALL QUIT('Out of memory in CC3_BMATT3ZU (TOT_T3Y) ') 2527 ENDIF 2528C 2529C--------------------------- 2530C Read Lambda_X matrices 2531C--------------------------- 2532C 2533 CALL GET_LAMBDAX(WORK(KLAMDPX),WORK(KLAMDHX),LISTX,IDLSTX, 2534 * ISYMRX,WORK(KLAMDP),WORK(KLAMDH),WORK(KEND1), 2535 * LWRK1) 2536C 2537C-------------------------- 2538C Read in T1^Z and T2^Z 2539C-------------------------- 2540C 2541 IOPT = 3 2542 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1Z),WORK(KT2Z),LISTZ,IDLSTZ, 2543 * ISYMRZ,WORK(KEND1),LWRK1) 2544C 2545 IF (IPRINT .GT. 55) THEN 2546 XNORMVAL = DDOT(NT1AM(ISYMRZ),WORK(KT1Z),1,WORK(KT1Z),1) 2547 WRITE(LUPRI,*) 'Norm of T1B ',XNORMVAL 2548 XNORMVAL = DDOT(NT2SQ(ISYMRZ),WORK(KT2Z),1,WORK(KT2Z),1) 2549 WRITE(LUPRI,*) 'Norm of T2B ',XNORMVAL 2550 ENDIF 2551C 2552C-------------------------- 2553C Read in T1^U and T2^U 2554C-------------------------- 2555C 2556 IOPT = 3 2557 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1U),WORK(KT2U),LISTU,IDLSTU, 2558 * ISYMRU,WORK(KEND1),LWRK1) 2559C 2560 IF (IPRINT .GT. 55) THEN 2561 XNORMVAL = DDOT(NT1AM(ISYMRU),WORK(KT1U),1,WORK(KT1U),1) 2562 WRITE(LUPRI,*) 'Norm of T1C ',XNORMVAL 2563 XNORMVAL = DDOT(NT2SQ(ISYMRU),WORK(KT2U),1,WORK(KT2U),1) 2564 WRITE(LUPRI,*) 'Norm of T2C ',XNORMVAL 2565 ENDIF 2566C 2567C------------------------------------------------------- 2568C Read in T1^ZU and T2^ZU !second-order amplitudes 2569C------------------------------------------------------- 2570C 2571 IOPT = 3 2572 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1ZU),WORK(KT2ZU),LISTZU,IDLSTZU, 2573 * ISYMRZU,WORK(KEND1),LWRK1) 2574C 2575 IF (IPRINT .GT. 55) THEN 2576 XNORMVAL = DDOT(NT1AM(ISYMRZU),WORK(KT1ZU),1,WORK(KT1ZU),1) 2577 WRITE(LUPRI,*) 'Norm of T1BC ',XNORMVAL 2578 XNORMVAL = DDOT(NT2SQ(ISYMRZU),WORK(KT2ZU),1,WORK(KT2ZU),1) 2579 WRITE(LUPRI,*) 'Norm of T2BC ',XNORMVAL 2580 ENDIF 2581C 2582C------------------ 2583C Read in T1^X 2584C------------------ 2585C 2586 IOPT = 1 2587 CALL GET_T1_T2(IOPT,.FALSE.,WORK(KT1X),DUMMY,LISTX,IDLSTX, 2588 * ISYMRX,WORK(KEND1),LWRK1) 2589C 2590 IF (IPRINT .GT. 55) THEN 2591 XNORMVAL = DDOT(NT1AM(ISYMRX),WORK(KT1X),1,WORK(KT1X),1) 2592 WRITE(LUPRI,*) 'Norm of T1A ',XNORMVAL 2593 ENDIF 2594C 2595C------------------------------------------ 2596C Calculate the F^X matrix 2597C------------------------------------------ 2598C 2599 CALL CC3LR_MFOCK(WORK(KEND1),WORK(KT1X),WORK(KXIAJB),ISYFCKX) 2600C 2601C Put the F_{kc} part into correct F_{pq} 2602C 2603 CALL DZERO(WORK(KFOCKX),N2BST(ISYFCKX)) 2604C 2605 DO ISYMC = 1, NSYM 2606 ISYMK = MULD2H(ISYFCKX,ISYMC) 2607 DO K = 1, NRHF(ISYMK) 2608 DO C = 1, NVIR(ISYMC) 2609 KOFF1 = KFOCKX - 1 2610 * + IFCVIR(ISYMK,ISYMC) 2611 * + NORB(ISYMK)*(C - 1) 2612 * + K 2613 KOFF2 = KEND1 - 1 2614 * + IT1AM(ISYMC,ISYMK) 2615 * + NVIR(ISYMC)*(K - 1) 2616 * + C 2617C 2618 WORK(KOFF1) = WORK(KOFF2) 2619C 2620 ENDDO 2621 ENDDO 2622 ENDDO 2623C 2624 IF (IPRINT .GT. 55) THEN 2625 CALL AROUND( 'In CC3_BMATT3ZU: Fock^A MO matrix' ) 2626 CALL CC_PRFCKMO(WORK(KFOCKX),ISYFCKX) 2627 ENDIF 2628C 2629C---------------------------------------------------- 2630C Integrals (ck|de)-tilde(Z) and (ck|lm)-tilde(Z) 2631C---------------------------------------------------- 2632C 2633 ISINTR1Z = MULD2H(ISINT1,ISYMRZ) 2634 ISINTR2Z = MULD2H(ISINT2,ISYMRZ) 2635C 2636 CALL CC3_BARINT(WORK(KT1Z),ISYMRZ,WORK(KLAMDP), 2637 * WORK(KLAMDH),WORK(KEND1),LWRK1, 2638 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 2639C 2640 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTR1Z,LU3SRTR,FN3SRTR, 2641 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 2642 * IDUMMY,CDUMMY) 2643C 2644 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTR1Z, 2645 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 2646C 2647C-------------------------------- 2648C Re-use some temporary files 2649C-------------------------------- 2650C 2651 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 2652 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 2653C 2654 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 2655 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 2656C 2657C------------------------------------------------------ 2658C Calculate the (ck|de)-{Z,U}tilde and (ck|lm)-{Z,U}tilde 2659C (double one-index transformed with first-order amplitudes) 2660C------------------------------------------------------ 2661C 2662 CALL CC3_3BARINT(ISYMRZ,LISTZ,IDLSTZ,ISYMRU,LISTU,IDLSTU, 2663 * IDUMMY,CDUMMY,IDUMMY,.FALSE., 2664 * WORK(KLAMDP),WORK(KLAMDH),WORK(KEND1),LWRK1, 2665 * LU3SRTR,FN3SRTR,LUCKJDR1,FNCKJDR1) 2666C 2667 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISYMZU,LU3SRTR,FN3SRTR, 2668 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 2669 * IDUMMY,CDUMMY) 2670C 2671 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISYMZU, 2672 * LUDELDR,FNDELDR,LUDKBCR1,FNDKBCR1) 2673C 2674C-------------------------------- 2675C Re-use some temporary files 2676C-------------------------------- 2677C 2678 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 2679 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 2680C 2681 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 2682 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 2683C 2684C---------------------------------------------------- 2685C Integrals (ck|de)-tilde(U) and (ck|lm)-tilde(U) 2686C---------------------------------------------------- 2687C 2688 ISINTR1U = MULD2H(ISINT1,ISYMRU) 2689 ISINTR2U = MULD2H(ISINT2,ISYMRU) 2690C 2691 CALL CC3_BARINT(WORK(KT1U),ISYMRU,WORK(KLAMDP), 2692 * WORK(KLAMDH),WORK(KEND1),LWRK1, 2693 * LU3SRTR,FN3SRTR,LUCKJDR2,FNCKJDR2) 2694C 2695 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTR1U,LU3SRTR,FN3SRTR, 2696 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 2697 * IDUMMY,CDUMMY) 2698C 2699 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTR1U, 2700 * LUDELDR,FNDELDR,LUDKBCR2,FNDKBCR2) 2701C 2702C-------------------------------- 2703C Re-use some temporary files 2704C-------------------------------- 2705C 2706 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 2707 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 2708C 2709 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 2710 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 2711C 2712C---------------------------------------------------------- 2713C Integrals (ck|de)-tilde(ZU) and (ck|lm)-tilde(ZU) 2714C (one-index transformed with second-order amplitudes) 2715C---------------------------------------------------------- 2716C 2717 ISINTR1ZU = MULD2H(ISINT1,ISYMRZU) 2718 ISINTR2ZU = MULD2H(ISINT2,ISYMRZU) 2719C 2720 CALL CC3_BARINT(WORK(KT1ZU),ISYMRZU,WORK(KLAMDP), 2721 * WORK(KLAMDH),WORK(KEND1),LWRK1, 2722 * LU3SRTR,FN3SRTR,LUCKJDR3,FNCKJDR3) 2723C 2724 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTR1ZU,LU3SRTR,FN3SRTR, 2725 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 2726 * IDUMMY,CDUMMY) 2727C 2728 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTR1ZU, 2729 * LUDELDR,FNDELDR,LUDKBCR3,FNDKBCR3) 2730C 2731C-------------------------------- 2732C Read occupied integrals 2733C------------------------------- 2734C 2735 ISINTZU = MULD2H(ISINT2,ISYMZU) 2736 ISINT1X = MULD2H(ISINT1,ISYMRX) 2737C 2738 !Use KEND0, because KT1Z, KT1U, KT1X, and KT1ZU are not needed any more 2739 KRBJIA = KEND0 2740 KTROC = KRBJIA + NT2SQ(ISYRES) 2741 KTROC1 = KTROC + NTRAOC(ISINT1X) !KTROC - int. in contraction 2742 KEND1 = KTROC1 + NTRAOC(ISINT1X) !KTROC1 - int. in contraction 2743 LWRK1 = LWORK - KEND1 2744C 2745 KTROC0Z = KEND1 2746 KEND1 = KTROC0Z + NTRAOC(ISINTR2Z)!KTROC0Z - int. in <mu3|[H^Z,T2^U]|HF> 2747C 2748 KTROC0U = KEND1 2749 KEND1 = KTROC0U + NTRAOC(ISINTR2U)!KTROC0U - int. in <mu3|[H^U,T2^Z]|HF>C 2750 KTROC0ZU = KEND1 2751 KEND1 = KTROC0ZU + NTRAOC(ISINTR2ZU)!KTROC0ZU <mu3|[H^ZU,T2^0]|HF> 2752C ==== 2753 KTROC3 = KEND1 2754 KEND1 = KTROC3 + NTRAOC(ISINTZU)!KTROC3- int. in <mu3|[H^{Z,U},T2]|HF> 2755C ======= 2756 KT3OG1 = KEND1 2757 KEND1 = KT3OG1 + NTRAOC(ISINT1) ! <mu3|[H^0,T2^ZU]|HF> 2758C === 2759 KINTOC = KEND1 2760 MAXOCC = MAX(NTOTOC(ISINTR2Z),NTOTOC(ISINTZU)) 2761 MMAXOCC = MAX(NTOTOC(ISINTR2U),MAXOCC) 2762 MMMAXOCC = MAX(NTOTOC(ISINTR2ZU),MMAXOCC) 2763 KEND2 = KINTOC + MAX(NTOTOC(ISINT1),MMMAXOCC)!KINTOC temporary storage 2764 LWRK2 = LWORK - KEND2 2765C 2766 IF (LWRK2 .LT. 0) THEN 2767 WRITE(LUPRI,*) 'Memory available : ',LWORK 2768 WRITE(LUPRI,*) 'Memory needed : ',KEND2 2769 CALL QUIT('Insufficient space in CC3_BMATT3ZU (2)') 2770 END IF 2771C 2772 CALL DZERO(WORK(KRBJIA),NT2SQ(ISYRES)) 2773C 2774C------------------------------------------------------------- 2775C occupied integrals for <mu3|[H^0,T2^ZU]|HF> 2776C------------------------------------------------------------- 2777C 2778 IOFF = 1 2779 IF (NTOTOC(ISINT1) .GT. 0) THEN 2780 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISINT1)) 2781 ENDIF 2782 2783 !Transform (ai| delta j) integrals to (ai|kj) and sort as I(k,i,j,a) 2784 CALL CC3_TROCC(WORK(KINTOC),WORK(KT3OG1),WORK(KLAMDP),WORK(KEND2), 2785 * LWRK2,ISINT1) 2786C 2787C------------------------------------------------------------- 2788C Z-transformed occupied integrals for <mu3|[H^Z,T2^U]|HF> 2789C------------------------------------------------------------- 2790C 2791 IOFF = 1 2792 IF (NTOTOC(ISINTR2Z) .GT. 0) THEN 2793 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF, 2794 * NTOTOC(ISINTR2Z)) 2795 ENDIF 2796C 2797 IF (IPRINT .GT. 55) THEN 2798 XNORMVAL = DDOT(NTOTOC(ISINTR2Z),WORK(KINTOC),1, 2799 * WORK(KINTOC),1) 2800 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (B transformed) ', 2801 * XNORMVAL 2802 ENDIF 2803C 2804 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0Z),WORK(KLAMDP), 2805 * WORK(KEND2),LWRK2,ISINTR2Z) 2806C 2807C------------------------------------------------------------- 2808C U-transformed occupied integrals for <mu3|[H^U,T2^Z]|HF> 2809C------------------------------------------------------------- 2810C 2811 IOFF = 1 2812 IF (NTOTOC(ISINTR2U) .GT. 0) THEN 2813 CALL GETWA2(LUCKJDR2,FNCKJDR2,WORK(KINTOC),IOFF, 2814 * NTOTOC(ISINTR2U)) 2815 ENDIF 2816C 2817 IF (IPRINT .GT. 55) THEN 2818 XNORMVAL = DDOT(NTOTOC(ISINTR2U),WORK(KINTOC),1, 2819 * WORK(KINTOC),1) 2820 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (C transformed) ', 2821 * XNORMVAL 2822 ENDIF 2823C 2824 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0U),WORK(KLAMDP), 2825 * WORK(KEND2),LWRK2,ISINTR2U) 2826C 2827C-------------------------------------------------------------------------- 2828C Z,U-transformed occupied integrals for <mu3|[H^{Z,U},T2]|HF> 2829C (double one-index transformed with first-order amplitudes) 2830C-------------------------------------------------------------------------- 2831C 2832 IOFF = 1 2833 IF (NTOTOC(ISINTZU) .GT. 0) THEN 2834 CALL GETWA2(LUCKJDR1,FNCKJDR1,WORK(KINTOC),IOFF, 2835 * NTOTOC(ISINTZU)) 2836 ENDIF 2837C 2838 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC3),WORK(KLAMDP), 2839 * WORK(KEND2),LWRK2,ISINTZU) 2840C 2841C--------------------------------------------------------------- 2842C ZU-transformed occupied integrals for <mu3|[H^ZU,T2^0]|HF> 2843C (one-index transformed with second-order amplitudes) 2844C--------------------------------------------------------------- 2845C 2846 IOFF = 1 2847 IF (NTOTOC(ISINTR2ZU) .GT. 0) THEN 2848 CALL GETWA2(LUCKJDR3,FNCKJDR3,WORK(KINTOC),IOFF, 2849 * NTOTOC(ISINTR2ZU)) 2850 ENDIF 2851C 2852 IF (IPRINT .GT. 55) THEN 2853 XNORMVAL = DDOT(NTOTOC(ISINTR2ZU),WORK(KINTOC),1, 2854 * WORK(KINTOC),1) 2855 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (BC transformed) ', 2856 * XNORMVAL 2857 ENDIF 2858C 2859 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0ZU),WORK(KLAMDP), 2860 * WORK(KEND2),LWRK2,ISINTR2ZU) 2861C 2862C----------------------------------- 2863C Occupied integrals in contraction 2864C----------------------------------- 2865C 2866 IOFF = 1 2867 IF (NTOTOC(ISINT1) .GT. 0) THEN 2868 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT1)) 2869 ENDIF 2870C 2871 !Write out norms of integrals. 2872 IF (IPRINT .GT. 55) THEN 2873 XNORMVAL = DDOT(NTOTOC(ISINT1),WORK(KINTOC),1,WORK(KINTOC),1) 2874 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XNORMVAL 2875 ENDIF 2876C 2877 !Transform (ia|j delta) integrals to (ia|j k)-tildeX and sort as (i,j,k,a) 2878 CALL CCLR_TROCC(WORK(KINTOC),WORK(KTROC),WORK(KLAMDHX),ISYMRX, 2879 * WORK(KEND2),LWRK2) 2880C 2881 !sort (i,j,k,a) as (a,i,j,k) 2882 CALL CCSDT_SRTOC2(WORK(KTROC),WORK(KTROC1),ISINT1X, 2883 * WORK(KEND2),LWRK2) 2884C 2885C--------------------- 2886C Start ISYMD-loop 2887C--------------------- 2888C 2889 DO ISYMD = 1,NSYM 2890C 2891 ISCKB1 = MULD2H(ISINT1,ISYMD) 2892 ISCKB1X = MULD2H(ISINT1X,ISYMD) 2893 ISCKB2Z = MULD2H(ISINTR2Z,ISYMD) 2894 ISCKB2U = MULD2H(ISINTR2U,ISYMD) 2895 ISCKB2ZU = MULD2H(ISINTR2ZU,ISYMD) 2896 ISCKBZU = MULD2H(ISINTZU,ISYMD) 2897C 2898C---------------------------------------- 2899C Read virtual integrals (fixed D) 2900C---------------------------------------- 2901C 2902 !Use KEND1, because KINTOC is not needed any more 2903 KTRVI = KEND1 2904 KTRVI1 = KTRVI + NCKATR(ISCKB1X) !KTRVI - int. in contraction 2905 KEND2 = KTRVI1 + NCKATR(ISCKB1X) !KTRVI1 - int. in contraction 2906 LWRK2 = LWORK - KEND2 2907C 2908 KTRVI0Z = KEND2 2909 KEND3 = KTRVI0Z + NCKATR(ISCKB2Z)!KTRVI0Z - int. in <mu3|[H^Z,T2^U]|HF> 2910 LWRK3 = LWORK - KEND3 ! === 2911C 2912 KTRVI0U = KEND3 2913 KEND3 = KTRVI0U + NCKATR(ISCKB2U)!KTRVI0U - int. in <mu3|[H^U,T2^Z]|HF> 2914 LWRK3 = LWORK - KEND3 ! === 2915C 2916 KTRVI0ZU = KEND3 2917 KEND3 = KTRVI0ZU + NCKATR(ISCKB2ZU)!KTRVI0ZU <mu3|[H^ZU,T2^0]|HF> 2918 LWRK3 = LWORK - KEND3 ! ==== 2919C 2920 KTRVI7 = KEND3 2921 KEND3 = KTRVI7 + NCKATR(ISCKBZU)!KTRVI7 int. in <mu3|[H^{Z,U},T2^0]|HF> 2922 LWRK3 = LWORK - KEND3 ! ====== 2923C 2924 KT3VDG1 = KEND3 2925 KEND3 = KT3VDG1 + NCKATR(ISCKB1)! <mu3|[H^0,T2^ZU]|HF> 2926C === 2927 KINTVI = KEND3 2928 KEND4 = KINTVI + NCKA(ISCKB1)!KINTVI - temporary storage 2929 LWRK4 = LWORK - KEND4 2930C 2931 IF (LWRK4 .LT. 0) THEN 2932 WRITE(LUPRI,*) 'Memory available : ',LWORK 2933 WRITE(LUPRI,*) 'Memory needed : ',KEND4 2934 CALL QUIT('Insufficient space in CC3_BMATT3ZU (3)') 2935 END IF 2936C 2937C-------------------- 2938C Start D-loop 2939C-------------------- 2940C 2941 DO D = 1,NVIR(ISYMD) 2942C 2943C---------------------------------------------------------------------------- 2944C virtual integrals for <mu3|[H^0,T2^ZU]|HF> (fixed D) 2945C---------------------------------------------------------------------------- 2946C 2947 IOFF = ICKBD(ISCKB1,ISYMD) + NCKATR(ISCKB1)*(D - 1) + 1 2948 IF (NCKATR(ISCKB1) .GT. 0) THEN 2949 CALL GETWA2(LUDKBC,FNDKBC,WORK(KT3VDG1),IOFF, 2950 * NCKATR(ISCKB1)) 2951 ENDIF 2952C 2953C---------------------------------------------------------------------------- 2954C Z-transformed virtual integrals for <mu3|[H^Z,T2^U]|HF> (fixed D) 2955C---------------------------------------------------------------------------- 2956C 2957 IOFF = ICKBD(ISCKB2Z,ISYMD) + NCKATR(ISCKB2Z)*(D - 1) + 1 2958 IF (NCKATR(ISCKB2Z) .GT. 0) THEN 2959 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI0Z),IOFF, 2960 & NCKATR(ISCKB2Z)) 2961 ENDIF 2962C 2963C---------------------------------------------------------------------------- 2964C U-transformed virtual integrals for <mu3|[H^U,T2^Z]|HF> (fixed D) 2965C---------------------------------------------------------------------------- 2966C 2967 IOFF = ICKBD(ISCKB2U,ISYMD) + NCKATR(ISCKB2U)*(D - 1) + 1 2968 IF (NCKATR(ISCKB2U) .GT. 0) THEN 2969 CALL GETWA2(LUDKBCR2,FNDKBCR2,WORK(KTRVI0U),IOFF, 2970 & NCKATR(ISCKB2U)) 2971 ENDIF 2972C 2973C----------------------------------------------------------------------------- 2974C Z,U-transformed virtual integrals for <mu3|[H^{Z,U},T2^0]|HF> 2975C (fixed D) 2976C (doubly one-index transformed with first-order amplitudes) 2977C----------------------------------------------------------------------------- 2978C 2979 IOFF = ICKBD(ISCKBZU,ISYMD) + NCKATR(ISCKBZU)*(D - 1) + 1 2980 IF (NCKATR(ISCKBZU) .GT. 0) THEN 2981 CALL GETWA2(LUDKBCR1,FNDKBCR1,WORK(KTRVI7),IOFF, 2982 & NCKATR(ISCKBZU)) 2983 ENDIF 2984C 2985C---------------------------------------------------------------------------- 2986C ZU-transformed virtual integrals for <mu3|[H^ZU,T2^0]|HF> (fixed D) 2987C (one-index transformed with second-order amplitudes) 2988C---------------------------------------------------------------------------- 2989C 2990 IOFF = ICKBD(ISCKB2ZU,ISYMD) + NCKATR(ISCKB2ZU)*(D - 1) + 1 2991 IF (NCKATR(ISCKB2ZU) .GT. 0) THEN 2992 CALL GETWA2(LUDKBCR3,FNDKBCR3,WORK(KTRVI0ZU),IOFF, 2993 & NCKATR(ISCKB2ZU)) 2994 ENDIF 2995C 2996C----------------------------------------------------- 2997C Virtual integrals in contraction (fixed D) 2998C----------------------------------------------------- 2999C 3000 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 3001 IF (NCKA(ISCKB1) .GT. 0) THEN 3002 CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF, 3003 * NCKA(ISCKB1)) 3004 ENDIF 3005 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVI),WORK(KLAMDPX), 3006 * ISYMRX, 3007 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3008C 3009 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 3010 IF (NCKA(ISCKB1) .GT. 0) THEN 3011 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 3012 * NCKA(ISCKB1)) 3013 ENDIF 3014 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVI1),WORK(KLAMDPX), 3015 * ISYMRX, 3016 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 3017C 3018C--------------------------- 3019C Start ISYMB-loop 3020C--------------------------- 3021C 3022 DO ISYMB = 1,NSYM 3023C 3024 ISYMZU = MULD2H(ISYMRZ,ISYMRU) 3025 ISYMBD = MULD2H(ISYMB,ISYMD) 3026 ISWMAT = MULD2H(ISYMZU,ISYMBD) 3027 ISCKD2Z = MULD2H(ISINTR2Z,ISYMB) 3028 ISCKD2U = MULD2H(ISINTR2U,ISYMB) 3029 ISCKD2ZU = MULD2H(ISINTR2ZU,ISYMB) 3030 ISCKDZU = MULD2H(ISINTZU,ISYMB) 3031 ISCKD1 = MULD2H(ISINT1,ISYMB) 3032C 3033 !Use KEND3, because KINTVI is not needed any more 3034 KDIAGW = KEND3 3035 KWMAT = KDIAGW + NCKIJ(ISWMAT) 3036 KINDSQW = KWMAT + NCKIJ(ISWMAT) 3037 KTMAT = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 3038 KEND4 = KTMAT + NCKIJ(ISWMAT) 3039C 3040 KTRVI8Z = KEND4 3041 KEND4 = KTRVI8Z + NCKATR(ISCKD2Z)!KTRVI8Z: <mu3|[H^Z,T2^U]|HF> 3042 LWRK4 = LWORK - KEND4 ! === 3043C 3044 KTRVI8U = KEND4 3045 KEND4 = KTRVI8U + NCKATR(ISCKD2U)!KTRVI8U: <mu3|[H^U,T2^Z]|HF> 3046 LWRK4 = LWORK - KEND4 ! === 3047C 3048 KTRVI8ZU = KEND4 3049 KEND4 = KTRVI8ZU + NCKATR(ISCKD2ZU)! <mu3|[H^ZU,T2^0]|HF> 3050 LWRK4 = LWORK - KEND4 ! ==== 3051C 3052 KT3VBG1 = KEND4 3053 KEND4 = KT3VBG1 + NCKATR(ISCKD1)! <mu3|[H^0,T2^ZU]|HF> 3054C 3055 KTRVI9 = KEND4 3056 KEND4 = KTRVI9 + NCKATR(ISCKDZU)!KTRVI9: <mu3|[H^{Z,U},T2^0]|HF> 3057 LWRK4 = LWORK - KEND4 ! ======= 3058C 3059 IF (LWRK4 .LT. 0) THEN 3060 WRITE(LUPRI,*) 'Memory available : ',LWORK 3061 WRITE(LUPRI,*) 'Memory needed : ',KEND4 3062 CALL QUIT('Insufficient space in CC3_BMATT3ZU (4)') 3063 END IF 3064C 3065C--------------------------------------------- 3066C Construct part of the diagonal. 3067C--------------------------------------------- 3068C 3069 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 3070C 3071 IF (IPRINT .GT. 55) THEN 3072 XNORMVAL = DDOT(NCKIJ(ISWMAT),WORK(KDIAGW),1, 3073 * WORK(KDIAGW),1) 3074 WRITE(LUPRI,*) 'Norm of DIA ',XNORMVAL 3075 ENDIF 3076C 3077C------------------------------------- 3078C Construct index arrays. 3079C------------------------------------- 3080C 3081 LENSQW = NCKIJ(ISWMAT) 3082 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 3083C 3084C-------------------------- 3085C Start B-loop 3086C-------------------------- 3087C 3088 DO B = 1,NVIR(ISYMB) 3089C 3090C--------------------------------------- 3091C Initialise 3092C--------------------------------------- 3093C 3094 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 3095C 3096C---------------------------------------------------------------------------- 3097C virtual integrals for <mu3|[H^0,T2^ZU]|HF> (fixed B) 3098C---------------------------------------------------------------------------- 3099C 3100 IOFF = ICKBD(ISCKD1,ISYMB) + NCKATR(ISCKD1)*(B - 1)+ 1 3101 IF (NCKATR(ISCKD1) .GT. 0) THEN 3102 CALL GETWA2(LUDKBC,FNDKBC,WORK(KT3VBG1),IOFF, 3103 * NCKATR(ISCKD1)) 3104 ENDIF 3105C 3106C---------------------------------------------------------------------------- 3107C Z-transformed virtual integrals for <mu3|[H^Z,T2^U]|HF> 3108C (fixed B) 3109C---------------------------------------------------------------------------- 3110C 3111 IOFF = ICKBD(ISCKD2Z,ISYMB) + NCKATR(ISCKD2Z)*(B-1) +1 3112 IF (NCKATR(ISCKD2Z) .GT. 0) THEN 3113 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI8Z),IOFF, 3114 * NCKATR(ISCKD2Z)) 3115 ENDIF 3116C 3117C---------------------------------------------------------------------------- 3118C U-transformed virtual integrals for <mu3|[H^U,T2^Z]|HF> 3119C (fixed B) 3120C---------------------------------------------------------------------------- 3121C 3122 IOFF = ICKBD(ISCKD2U,ISYMB) + NCKATR(ISCKD2U)*(B-1) +1 3123 IF (NCKATR(ISCKD2U) .GT. 0) THEN 3124 CALL GETWA2(LUDKBCR2,FNDKBCR2,WORK(KTRVI8U),IOFF, 3125 * NCKATR(ISCKD2U)) 3126 ENDIF 3127C 3128C----------------------------------------------------------------------------- 3129C Z,U-transformed virtual integrals for <mu3|[H^{Z,U},T2^0]|HF> 3130C (fixed B) 3131C (doubly one-index transformed with first-order amplitudes) 3132C----------------------------------------------------------------------------- 3133C 3134 IOFF = ICKBD(ISCKDZU,ISYMB) + NCKATR(ISCKDZU)*(B-1) +1 3135 IF (NCKATR(ISCKDZU) .GT. 0) THEN 3136 CALL GETWA2(LUDKBCR1,FNDKBCR1,WORK(KTRVI9),IOFF, 3137 * NCKATR(ISCKDZU)) 3138 ENDIF 3139C 3140C---------------------------------------------------------------------------- 3141C ZU-transformed virtual integrals for <mu3|[H^ZU,T2^0]|HF> 3142C (fixed B) 3143C (one-index transformed with second-order amplitudes) 3144C---------------------------------------------------------------------------- 3145C 3146 IOFF = ICKBD(ISCKD2ZU,ISYMB)+ NCKATR(ISCKD2ZU)*(B-1)+1 3147 IF (NCKATR(ISCKD2ZU) .GT. 0) THEN 3148 CALL GETWA2(LUDKBCR3,FNDKBCR3,WORK(KTRVI8ZU),IOFF, 3149 * NCKATR(ISCKD2ZU)) 3150 ENDIF 3151C 3152C---------------------------------------------- 3153C Calculate <mu3|[H^U,T2^Z]|HF> !b-matrix contribution 3154C---------------------------------------------- 3155C 3156 CALL WBD_GROUND(WORK(KT2Z),ISYMRZ,WORK(KTMAT), 3157 * WORK(KTRVI0U),WORK(KTRVI8U), 3158 * WORK(KTROC0U),ISINTR2U,WORK(KWMAT), 3159 * WORK(KEND4),LWRK4, 3160 * WORK(KINDSQW),LENSQW, 3161 * ISYMB,B,ISYMD,D) 3162C 3163C---------------------------------------------- 3164C Calculate <mu3|[H^Z,T2^U]|HF> !b-matrix contribution 3165C---------------------------------------------- 3166C 3167 CALL WBD_GROUND(WORK(KT2U),ISYMRU,WORK(KTMAT), 3168 * WORK(KTRVI0Z),WORK(KTRVI8Z), 3169 * WORK(KTROC0Z),ISINTR2Z,WORK(KWMAT), 3170 * WORK(KEND4),LWRK4, 3171 * WORK(KINDSQW),LENSQW, 3172 * ISYMB,B,ISYMD,D) 3173C 3174C---------------------------------------------------------- 3175C Calculate <mu3|[[[H,T1^Z],T1^U],T2^0]|HF> !b-matrix contrib. 3176C---------------------------------------------------------- 3177C 3178 CALL WBD_GROUND(WORK(KT2TP),ISYMT2,WORK(KTMAT), 3179 * WORK(KTRVI7),WORK(KTRVI9), 3180 * WORK(KTROC3),ISINTZU,WORK(KWMAT), 3181 * WORK(KEND4),LWRK4, 3182 * WORK(KINDSQW),LENSQW, 3183 * ISYMB,B,ISYMD,D) 3184C 3185C---------------------------------------------------- 3186C Calculate <mu3|[[H,T1^ZU],T2^0]|HF> !jacobian contribution 3187C---------------------------------------------------- 3188C 3189 CALL WBD_GROUND(WORK(KT2TP),ISYMT2,WORK(KTMAT), 3190 * WORK(KTRVI0ZU),WORK(KTRVI8ZU), 3191 * WORK(KTROC0ZU),ISINTR2ZU,WORK(KWMAT), 3192 * WORK(KEND4),LWRK4, 3193 * WORK(KINDSQW),LENSQW, 3194 * ISYMB,B,ISYMD,D) 3195C 3196C----------------------------------------------- 3197C Calculate <mu3|[H^0,T2^ZU]|HF> !jacobian contribution 3198C----------------------------------------------- 3199C 3200 CALL WBD_GROUND(WORK(KT2ZU),ISYMRZU,WORK(KTMAT), 3201 * WORK(KT3VDG1),WORK(KT3VBG1), 3202 * WORK(KT3OG1),ISINT1,WORK(KWMAT), 3203 * WORK(KEND4),LWRK4, 3204 * WORK(KINDSQW),LENSQW, 3205 * ISYMB,B,ISYMD,D) 3206C 3207C---------------------------------------------------- 3208C Divide by orbital energy difference 3209C and remove the forbidden elements 3210C---------------------------------------------------- 3211C 3212C 3213 CALL T3_FORBIDDEN(WORK(KWMAT),ISYMZU, 3214 * ISYMB,B,ISYMD,D) 3215 3216c call sum_pt3(work(kwmat),isymb,b,isymd,d, 3217c * iswmat,work(kx3am),4) 3218c write(lupri,*) 'Total norm of WBD : ', 3219c * ddot(nckij(iswmat),work(kwmat),1,work(kwmat),1) 3220 3221 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQZU, 3222 * ISWMAT,WORK(KWMAT), 3223 * WORK(KDIAGW),WORK(KFOCKD)) 3224 3225C 3226C----------------------------------------------- 3227C Carry out the contractions... 3228C----------------------------------------------- 3229C 3230C 3231C 3232C------------------------------------------------------ 3233C Calculate the term <mu2|[H,W^BD(3)]|HF> 3234C ( Fock matrix cont ) 3235C added in OMEGA2EFF 3236C------------------------------------------------------ 3237C 3238 CALL CC3_CY2F(OMEGA2EFF,ISYRES,WORK(KWMAT),ISWMAT, 3239 * WORK(KTMAT),WORK(KFOCKX),ISYFCKX, 3240 * WORK(KINDSQW), 3241 * LENSQW,WORK(KEND4),LWRK4, 3242 * ISYMB,B,ISYMD,D) 3243C 3244C------------------------------------------------------ 3245C Calculate the term <mu2|[H,W^BD(3)]|HF> 3246C ( Occupied cont ) 3247C added in OMEGA2EFF 3248C------------------------------------------------------ 3249C 3250 CALL CC3_CY2O(OMEGA2EFF,ISYRES,WORK(KWMAT),ISWMAT, 3251 * WORK(KTMAT),WORK(KTROC), 3252 * WORK(KTROC1),ISINT1,WORK(KEND4),LWRK4, 3253 * WORK(KINDSQW),LENSQW, 3254 * ISYMB,B,ISYMD,D) 3255C 3256C------------------------------------------------------ 3257C Calculate the term <mu2|[H,W^BD(3)]|HF> 3258C ( Virtual cont ) 3259C added in OMEGA2EFF 3260C------------------------------------------------------ 3261C 3262 CALL CC3_CY2V(OMEGA2EFF,ISYRES,WORK(KRBJIA), 3263 * WORK(KWMAT), 3264 * ISWMAT,WORK(KTMAT),WORK(KTRVI), 3265 * WORK(KTRVI1),ISINT1,WORK(KEND4),LWRK4, 3266 * WORK(KINDSQW),LENSQW, 3267 * ISYMB,B,ISYMD,D) 3268C 3269 END DO !B 3270 END DO !ISYMB 3271C 3272 END DO !D 3273 END DO !ISYMD 3274c 3275c 3276c write(lupri,*) 'T3XY in CC3_BMATT3ZU : ' 3277c call print_pt3(work(kx3am),1,4) 3278c 3279 3280C 3281C 3282C------------------------------------------------------ 3283C Accumulate RBJIA from <mu2|[H,W^BD(3)]|HF> ( Virtual cont ) 3284C in OMEGA2EFF 3285C------------------------------------------------------ 3286C 3287 CALL CC3_RBJIA(OMEGA2EFF,ISYRES,WORK(KRBJIA)) 3288c 3289c write(lupri,*)'OMEGA2EFF after CC3_BMATT3ZU' 3290c call outpak(OMEGA2EFF,NT1AM(1),1,lupri) 3291 3292C 3293 IF (IPRINT .GT. 55) THEN 3294 XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 3295 WRITE(LUPRI,*)'Norm of OMEGA2EFF final before added ',XNORMVAL 3296 ENDIF 3297C 3298 DO I = 1,NT2AM(ISYRES) 3299 OMEGA2EFF(I) = OMEGA2EFF(I) + OMEGA2(I) 3300 END DO 3301C 3302 IF (IPRINT .GT. 55) THEN 3303 XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 3304 WRITE(LUPRI,*)'Norm of OMEGA2EFF final, OMEGA2EFF + OMEGA2F ', 3305 * XNORMVAL 3306 ENDIF 3307C 3308C------------------------------------------- 3309C Close files (integrals in contraction) 3310C------------------------------------------- 3311C 3312 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 3313 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 3314 CALL WCLOSE2(LU3VI2,FN3VI2,'KEEP') 3315C 3316C------------------------------------------ 3317C Close files (H0 integrals for WBD) 3318C------------------------------------------ 3319C 3320 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 3321 CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP') 3322C 3323C------------------------------------------ 3324C Close temporary files (H^Z integrals) 3325C------------------------------------------ 3326C 3327 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 3328 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 3329 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 3330 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 3331C 3332C------------------------------------------ 3333C Close temporary files (H^{Z,U} integrals) 3334C------------------------------------------ 3335C 3336 CALL WCLOSE2(LUCKJDR1,FNCKJDR1,'DELETE') 3337 CALL WCLOSE2(LUDKBCR1,FNDKBCR1,'DELETE') 3338C 3339C------------------------------------------ 3340C Close temporary files (H^U integrals) 3341C------------------------------------------ 3342C 3343 CALL WCLOSE2(LUCKJDR2,FNCKJDR2,'DELETE') 3344 CALL WCLOSE2(LUDKBCR2,FNDKBCR2,'DELETE') 3345C 3346C------------------------------------------ 3347C Close temporary files (H^ZU integrals) 3348C------------------------------------------ 3349C 3350 CALL WCLOSE2(LUCKJDR3,FNCKJDR3,'DELETE') 3351 CALL WCLOSE2(LUDKBCR3,FNDKBCR3,'DELETE') 3352C 3353C------------- 3354C End 3355C------------- 3356C 3357 CALL QEXIT('CC3_BMATT3ZU ') 3358C 3359 RETURN 3360 END 3361C 3362