1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck cc3_xi */ 20 SUBROUTINE CCSDT_FBMAT(IBTRAN,IDOTS,DOTPROD,NBTRAN,MXVEC, 21 * LISTL,LISTB,LISTC,FNCKJD,FNDKBC,FNDELD, 22 * FNTOC,FN3VI,FN3VI2,WORK,LWORK) 23C 24C Written by K. Hald, Summer 2002. 25C 26C Calculate <L2|[[H,R_1],R_3]|HF> 27C 28C Initially construct 29C 30C t3(ai,Bj,Ck) = S^BC(aikj) + U^BC(aikj) + S^CB(aijk) + U^CB(aijk) 31C 32C W^BC(aikj) = sum_d X(ad)*t3(di,Bj,Ck) - sum_l X(li)*t3(al,Bj,Ck) 33C - sum_ld X(ld)*t2(ai,Cl)*t2(Bj,dk) 34C - sum_ld X(ld)*t2(ai,Bl)*t2(Ck,dj) 35C + ( -<mu3 | [[H,T1^Y],T2] | HF> -<mu3| [H,T2^Y] |HF> ) 36C 37************************************************************************** 38* 39* Spring 2004, Filip Pawlowski, Aarhus: Modified to treat also the 40* Cauchy moments. 41************************************************************************** 42 IMPLICIT NONE 43C 44#include "priunit.h" 45#include "dummy.h" 46#include "iratdef.h" 47#include "ccsdsym.h" 48#include "inftap.h" 49#include "ccinftap.h" 50#include "ccorb.h" 51#include "ccsdinp.h" 52#include "ccr1rsp.h" 53#include "ccexci.h" 54#include "ccrc1rsp.h" 55C 56 LOGICAL LCAUCHY 57C 58 CHARACTER*(*) FNCKJD, FNDKBC, FNDELD, FNTOC, FN3VI, FN3VI2 59 CHARACTER*12 FN3SRTR, FNDELDR 60 CHARACTER*16 FNCKJDR, FNDKBCR 61 CHARACTER*10 MODEL 62 CHARACTER*8 LABELB, LABELC 63 CHARACTER*3 LISTL, LISTB, LISTC 64 CHARACTER CDUMMY*1 65C 66 PARAMETER(FN3SRTR = 'CCSDT_FBMAT1', FNDELDR = 'CCSDT_FBMAT3') 67C 68 INTEGER NBTRAN, MXVEC, LWORK 69 INTEGER IBTRAN(3,NBTRAN), IDOTS(MXVEC,NBTRAN) 70 INTEGER ISINT1, ISINT2, ISYMT2, ISYMT3, ISYMW 71 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR, LUDKBCR4 72 INTEGER KT1AM, KT2TP, KRBJIA, KFOCKD, KEND0, LWRK0, IOPT 73 INTEGER IDLSTB, ISYMB, KT1B, KT2B 74 INTEGER ISINTR1, ISINTR2, KLAMDP, KLAMDH, KEND1, LWRK1 75 INTEGER KTROC, KTROC1, KTROC0, KTROC02, KTROC0Y 76 INTEGER KINTOC, KEND2, LWRK2 77 INTEGER IOFF, ISYMD, ISCKB1, ISCKB2, ISCKB2Y 78 INTEGER KTRVI, KTRVI1, KTRVI2, KTRVI3, KTRVI0, KTRVI0Y 79 INTEGER KEND3, LWRK3, KEND4, LWRK4, KINTVI, LENSQ, LENSQW 80 INTEGER ISYALJ, ISYALJ2, ISYMBD, ISCKIJ, ISCKD2, ISCKD2Y 81 INTEGER ISWMAT, KSMAT, KSMAT3, KUMAT, KUMAT3, KDIAG, KDIAGW 82 INTEGER KWMAT, KINDSQ, KINDSQW, KINDEX, KINDEX2, KTMAT 83 INTEGER KTRVI8, KTRVI8Y, KTRVI9, KTRVI10, KEND5, LWRK5 84 INTEGER LUCKJD, LUTOC, LU3VI, LU3VI2, LUDKBC, LUDELD 85 INTEGER KOMEG2, KSAVE, KEND00, LWRK00, KEND01, LWRK01 86 INTEGER ITRAN, IDLSTC, ISYMRB, ISYMRC, ISINT1C, ISYRES 87 INTEGER KFCKB, LUFCK, IVEC, IDLSTL, ISYML 88 INTEGER KLAMPC, KLAMHC, KT1C, ISYFCKB, ISYFCKC, KFCKC 89 INTEGER IRREP, IERR, KXIAJB, IOPTTCME, ISYMK, KOFF1, KOFF2 90 INTEGER ISYMC, IRREP2, ISCKB1C, ISYCKB 91C 92 integer kx3am 93C 94 INTEGER IDLSTBM 95 INTEGER NCAU,MCAU,KOFFOCC,KOFFINTD,KOFFINTB 96C 97C Functions : 98C 99 INTEGER ILSTSYM,ILRCAMP 100C 101#if defined (SYS_CRAY) 102 REAL DOTPROD(MXVEC,NBTRAN), WORK(LWORK) 103 REAL FREQB, FREQC, TCON 104 REAL XNORMVAL, DDOT, HALF, TWO, ZERO 105#else 106 DOUBLE PRECISION DOTPROD(MXVEC,NBTRAN), WORK(LWORK) 107 DOUBLE PRECISION FREQB, FREQC, TCON 108 DOUBLE PRECISION XNORMVAL, DDOT, HALF, TWO, ZERO 109#endif 110C 111 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, TWO = 2.0D0) 112C 113 CALL QENTER('FBMAT') 114C 115C-------------------------- 116C Save and set flags. 117C-------------------------- 118C 119C Set symmetry flags. 120C 121C 122C ISYMT2 is symmetry of T2TP 123C ISINT2 is symmetry of integrals in triples equation (ISINT2=1) 124C ISINT1 is symmetry of integrals in contraction (ISINT1=1) 125C ISYMT3 is symmetry of S and U intermediate 126C ISYMW is symmetry of W intermediate 127C ISYRES is symmetry of result vector (here the same as ISYFKY) 128C ISYMW = ISYFKY*ISYMT3 129C ISYRES = ISYMT3*ISYFKY*ISINT1 130C 131C------------------------------------------------------------- 132C 133 IPRCC = IPRINT 134 ISINT1 = 1 135 ISINT2 = 1 136 ISYMT2 = 1 137 ISYMT3 = MULD2H(ISINT2,ISYMT2) 138C 139C-------------------------------- 140C Open files 141C-------------------------------- 142C 143 LUCKJD = -1 144 LUTOC = -1 145 LU3VI = -1 146 LU3VI2 = -1 147 LUDKBC = -1 148 LUDELD = -1 149C 150 CALL WOPEN2(LUCKJD,FNCKJD,64,0) 151 CALL WOPEN2(LUTOC,FNTOC,64,0) 152 CALL WOPEN2(LU3VI,FN3VI,64,0) 153 CALL WOPEN2(LU3VI2,FN3VI2,64,0) 154 CALL WOPEN2(LUDKBC,FNDKBC,64,0) 155 CALL WOPEN2(LUDELD,FNDELD,64,0) 156C 157C---------------------------------------------- 158C Calculate the zeroth order stuff once 159C---------------------------------------------- 160C 161 KT2TP = 1 162 KFOCKD = KT2TP + NT2SQ(ISYMT2) 163 KLAMDP = KFOCKD + NORBTS 164 KLAMDH = KLAMDP + NLAMDT 165 KXIAJB = KLAMDH + NLAMDT 166 KEND00 = KXIAJB + NT2AM(ISINT1) 167 LWRK00 = LWORK - KEND00 168C 169 KT1AM = KEND00 170 KEND01 = KT1AM + NT1AM(ISYMT2) 171 LWRK01 = LWORK - KEND01 172C 173 IF (LWRK01 .LT. 0) THEN 174 CALL QUIT('Out of memory in CCSDT_FBMAT (zeroth allo.') 175 ENDIF 176C 177C------------------------ 178C Construct L(ia,jb). 179C------------------------ 180C 181 REWIND(LUIAJB) 182 CALL READI(LUIAJB,IRAT*NT2AM(ISINT1),WORK(KXIAJB)) 183 IOPTTCME = 1 184 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISINT1,IOPTTCME) 185C 186 IF ( IPRINT .GT. 55) THEN 187 XNORMVAL = DDOT(NT2AM(ISINT1),WORK(KXIAJB),1, 188 * WORK(KXIAJB),1) 189 WRITE(LUPRI,*) 'Norm of IAJB ',XNORMVAL 190 ENDIF 191C 192C---------------------------------------------------------------- 193C Read t1 and calculate the zero'th order Lambda matrices 194C---------------------------------------------------------------- 195C 196 IOPT = 1 197 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 198C 199 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 200 & WORK(KEND01),LWRK01) 201C 202C------------------------------------------- 203C Read in t2 , square it and reorder 204C------------------------------------------- 205C 206 IOPT = 2 207 CALL CC_RDRSP('R0',0,ISYMT2,IOPT,MODEL,DUMMY,WORK(KEND00)) 208 CALL CC_T2SQ(WORK(KEND00),WORK(KT2TP),ISYMT2) 209 IF (LWRK00 .LT. NT2SQ(ISYMT2)) THEN 210 CALL QUIT('Not enough memory to construct T2TP in CCSDT_FBMAT') 211 ENDIF 212C 213 CALL DCOPY(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KEND00),1) 214 CALL CC3_T2TP(WORK(KT2TP),WORK(KEND00),ISYMT2) 215C 216 IF (IPRINT .GT. 55) THEN 217 XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 218 WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL 219 ENDIF 220C 221C-------------------------------------- 222C Read in orbital energies 223C-------------------------------------- 224C 225 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 226 & .FALSE.) 227 REWIND LUSIFC 228C 229 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 230 READ (LUSIFC) 231 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 232C 233 CALL GPCLOSE(LUSIFC,'KEEP') 234C 235C--------------------------------------------- 236C Delete frozen orbitals in Fock diagonal. 237C--------------------------------------------- 238C 239 IF (FROIMP .OR. FROEXP) 240 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND00),LWRK00) 241C 242C If we want to sum the T3 amplitudes 243C 244 if (.false.) then 245 kx3am = kend00 246 kend00 = kx3am + nt1amx*nt1amx*nt1amx 247 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 248 lwrk00 = lwork - kend00 249 if (lwrk00 .lt. 0) then 250 write(lupri,*) 'Memory available : ',lwork 251 write(lupri,*) 'Memory needed : ',kend00 252 call quit('Insufficient space (T3) in CCSDT_FBMAT') 253 END IF 254 endif 255C 256C---------------------------------------------------- 257C Loop over left and right amplitude vectors: 258C---------------------------------------------------- 259C 260 DO ITRAN = 1, NBTRAN 261C 262 IDLSTB = IBTRAN(1,ITRAN) 263 IDLSTC = IBTRAN(2,ITRAN) 264C 265 IF ( (LISTB(1:3).EQ.'RE ') .AND. (LISTC(1:3).EQ.'RE ') ) THEN 266 WRITE(LUPRI,*)'LISTB = ',LISTB 267 WRITE(LUPRI,*)'LISTC = ',LISTC 268 WRITE(LUPRI,*)'LISTB and LISTC equal RE not implemented !!!' 269 CALL QUIT('Not implemented in CCSDT_FBMAT.') 270 END IF 271C 272 IF ( (LISTB(1:3).EQ.'RC ') .AND. (LISTC(1:3).NE.'RC ') ) THEN 273 WRITE(LUPRI,*)'LISTB = ',LISTB 274 WRITE(LUPRI,*)'LISTC = ',LISTC 275 WRITE(LUPRI,*)'LISTB and LISTC have to be RC simultanously!' 276 CALL QUIT('Not implemented in CCSDT_FBMAT.') 277 END IF 278C 279 IF ( (LISTC(1:3).EQ.'RC ') .AND. (LISTB(1:3).NE.'RC ') ) THEN 280 WRITE(LUPRI,*)'LISTB = ',LISTB 281 WRITE(LUPRI,*)'LISTC = ',LISTC 282 WRITE(LUPRI,*)'LISTB and LISTC have to be RC simultanously!' 283 CALL QUIT('Not implemented in CCSDT_FBMAT.') 284 END IF 285C 286 !initialize LCAUCHY and NCAU 287 LCAUCHY = .FALSE. 288 NCAU = 0 289C 290 IF (LISTB(1:3).EQ.'R1 ') THEN 291cfp here was the problem 292c originally I had: 293c ISYMRB = ISYLRT(IDLSTB) 294c but now I put: 295 ISYMRB = ILSTSYM(LISTB,IDLSTB) 296cfp 297 FREQB = FRQLRT(IDLSTB) 298 LABELB = LRTLBL(IDLSTB) 299 ELSE IF (LISTB(1:3).EQ.'RE ') THEN 300 ISYMRB = ILSTSYM(LISTB,IDLSTB) 301 FREQB = EIGVAL(IDLSTB) 302 LABELB = '- none -' 303 ELSE IF (LISTB(1:3).EQ.'RC ') THEN 304* ISYMRB = ILSTSYM(LISTB,IDLSTB) 305 ISYMRB = ISYLRC(IDLSTB) 306 FREQB = ZERO 307 LABELB = LRCLBL(IDLSTB) 308 NCAU = ILRCAU(IDLSTB) 309 LCAUCHY = .TRUE. 310 ELSE 311 WRITE(LUPRI,*)'LISTB = ',LISTB 312 CALL QUIT('Unknown list for LISTB in CCSDT_FBMAT.') 313 END IF 314C 315 IF ( (LISTC(1:3).EQ.'R1 ') .OR. (LISTC(1:3).EQ.'R2 ') ) THEN 316cfp here was the problem 317c originally I had: 318c ISYMRC = ISYLRT(IDLSTC) 319c but now I put: 320 ISYMRC = ILSTSYM(LISTC,IDLSTC) 321cfp 322 FREQC = FRQLRT(IDLSTC) 323 LABELC = LRTLBL(IDLSTC) 324 ELSE IF (LISTC(1:3).EQ.'RE ') THEN 325 ISYMRC = ILSTSYM(LISTC,IDLSTC) 326 FREQC = EIGVAL(IDLSTC) 327 LABELC = '- none -' 328 ELSE IF (LISTC(1:3).EQ.'ER1') THEN 329 ISYMRC = ILSTSYM(LISTC,IDLSTC) 330 FREQC = DUMMY 331 LABELC = '- none -' 332 ELSE IF (LISTC(1:3).EQ.'RC ') THEN 333* ISYMRC = ILSTSYM(LISTC,IDLSTC) 334 ISYMRC = ISYLRC(IDLSTC) 335 FREQC = ZERO 336 LABELC = LRCLBL(IDLSTC) 337* NCAU = ILRCAU(IDLSTC) 338 ELSE 339 WRITE(LUPRI,*)'LISTC = ',LISTC 340 CALL QUIT('Unknown list for LISTC in CCSDT_FBMAT.') 341 END IF 342C 343 ISYMW = MULD2H(ISYMT3,ISYMRB) 344 ISINT1C = MULD2H(ISINT1,ISYMRC) 345 ISYRES = MULD2H(ISYMW,ISINT1C) 346C 347 ISYFCKB = MULD2H(ISYMOP,ISYMRB) 348 ISYFCKC = MULD2H(ISYMOP,ISYMRC) 349C 350C------------------------------------------------- 351C Read T1^B and T2^B 352C Calculate (ck|de)-tilde(B) and (ck|lm)-tilde(B) 353C Used to construct T3^B 354C------------------------------------------------- 355C 356 KT2B = KEND00 357 KFCKB = KT2B + (NCAU+1)*NT2SQ(ISYMRB) 358 KFCKC = KFCKB + N2BST(ISYFCKB) 359 KEND0 = KFCKC + N2BST(ISYFCKC) 360 LWRK0 = LWORK - KEND0 361C 362 KT1B = KEND0 363 KEND1 = KT1B + (NCAU+1)*NT1AM(ISYMRB) 364 LWRK1 = LWORK - KEND1 365C 366 IF (LWRK1 .LT. NT2SQ(ISYMRB)) THEN 367 WRITE(LUPRI,*)'Memory available :',LWRK1 368 WRITE(LUPRI,*)'Memory needed :',NT2SQ(ISYMRB) 369 CALL QUIT('Out of memory in CCSDT_FBMAT (T2) ') 370 ENDIF 371 372C 373C-------------------------------------------------------------------------- 374C 1) get KT1B and KT2B vectors 375C 2) get T1B-transformed integrals needed in construction of W3 376C ( <mu3|[[H^0,T1B],T2^0]|HF> term). 377C 378C For non-Cacuhy calculations (LCAUCHY = .FALSE., NCAU = 0) the 379C loop below is dummy. 380C 381C For Cacuhy calculations (LCAUCHY = .TRUE., NCAU >= 0) 382C KT1B and KT2B correspond to singles and doubles part of Cauchy 383C vvectors. 384C-------------------------------------------------------------------------- 385C 386 DO MCAU = 0, NCAU 387C 388 !Open temporary files 389 LU3SRTR = -1 390 LUDELDR = -1 391 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 392 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 393C 394 IF (LCAUCHY) THEN 395 !Get the list for current MCAU 396 IDLSTBM = ILRCAMP(LABELB,MCAU,ISYMRB) 397 !Consistency check 398 IF (MCAU.EQ.NCAU .AND. IDLSTBM.NE.IDLSTB) THEN 399 CALL QUIT('Inconsistency in Cauchy loop in CCSDT_FBMAT') 400 END IF 401 ELSE 402 IDLSTBM = IDLSTB 403 END IF 404C 405 !Read in C1 and C2 Cauchy vectors for each MCAU 406 IOPT = 3 407 KOFF1 = KT1B + MCAU*NT1AM(ISYMRB) 408 KOFF2 = KT2B + MCAU*NT2SQ(ISYMRB) 409 CALL CC_RDRSP(LISTB,IDLSTBM,ISYMRB,IOPT,MODEL, 410 * WORK(KOFF1),WORK(KOFF2)) 411 CALL CCLR_DIASCL(WORK(KOFF2),TWO,ISYMRB) 412 CALL CC_T2SQ(WORK(KOFF2),WORK(KEND1),ISYMRB) 413 CALL CC3_T2TP(WORK(KOFF2),WORK(KEND1),ISYMRB) 414C 415 IF (IPRINT .GT. 55) THEN 416 XNORMVAL = DDOT(NT1AM(ISYMRB),WORK(KOFF1),1,WORK(KOFF1),1) 417 WRITE(LUPRI,*) 'Norm of T1B ',XNORMVAL 418 XNORMVAL = DDOT(NT2SQ(ISYMRB),WORK(KOFF2),1,WORK(KOFF2),1) 419 WRITE(LUPRI,*) 'Norm of T2B ',XNORMVAL 420 ENDIF 421C 422 ISINTR1 = MULD2H(ISINT1,ISYMRB) 423 ISINTR2 = MULD2H(ISINT2,ISYMRB) 424C 425 !Generate file names for T1-transformed integrals 426 IF (LCAUCHY) THEN 427 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 428 ELSE 429 FNCKJDR = 'CCSDT_____FBMAT2' 430 FNDKBCR = 'CCSDT_____FBMAT4' 431 END IF 432C 433 !Open the files for T1-transformed integrals 434 LUCKJDR = -1 435 LUDKBCR = -1 436 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 437 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 438C 439 !Construct T1-transformed integrals and store on file 440 CALL CC3_BARINT(WORK(KOFF1),ISYMRB,WORK(KLAMDP), 441 * WORK(KLAMDH),WORK(KEND1),LWRK1, 442 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 443C 444 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTR1,LU3SRTR,FN3SRTR, 445 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 446 * IDUMMY,CDUMMY) 447C 448 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTR1, 449 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 450C 451 !Close the files for T1-transformed integrals 452 CALL WCLOSE2(LUCKJDR,FNCKJDR,'KEEP') 453 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 454C 455 !Close and delete temporary files 456 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 457 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 458C 459 END DO !MCAU 460C 461C------------------------------------------ 462C Calculate the F^B matrix 463C------------------------------------------ 464C 465 !Use KEND0 since KT1B is not needed any longer 466 467 IF ((LISTB(1:3).EQ.'R1 ').OR.(LISTB(1:3).EQ.'RC ')) THEN 468C 469 CALL CCPRPAO(LABELB,.TRUE.,WORK(KFCKB),ISYFCKB,IRREP,IERR, 470 * WORK(KEND0),LWRK0) 471C 472 IF (IERR.GT.0) THEN 473 CALL QUIT('CCSDT_FBMAT : error reading operator '//LABELB) 474 ELSE IF (IERR.LT.0) THEN 475 CALL DZERO(WORK(KFCKB),N2BST(ISYMRB)) 476 END IF 477C 478 CALL CC_FCKMO(WORK(KFCKB),WORK(KLAMDP),WORK(KLAMDH), 479 * WORK(KEND0),LWRK0,ISYMRB,ISYMOP,ISYMOP) 480C 481 IF (IPRINT .GT. 50) THEN 482 CALL AROUND( 'In CCSDT_FBMAT: Fock^B MO matrix' ) 483 CALL CC_PRFCKMO(WORK(KFCKB),ISYFCKB) 484 ENDIF 485C 486 END IF 487C 488C--------------------------------------------- 489C Create lam^C 490C--------------------------------------------- 491C 492 KLAMPC = KEND0 493 KLAMHC = KLAMPC + NLAMDT 494 KEND0 = KLAMHC + NLAMDT 495 LWRK0 = LWORK - KEND0 496C 497 KT1C = KEND0 498 KEND1 = KT1C + NT1AM(ISYMRC) 499 LWRK1 = LWORK - KEND1 500C 501 IF (LWRK1.LT.NT1AM(ISYFCKC)) THEN 502 WRITE(LUPRI,*)'Memory available: ', LWRK1 503 WRITE(LUPRI,*)'Memory needed : ', NT1AM(ISYFCKC) 504 CALL QUIT('Insufficient memory in CCSDT_FBMAT (FCKC)') 505 END IF 506C 507 IOPT = 1 508 CALL CC_RDRSP(LISTC,IDLSTC,ISYMRC,IOPT,MODEL, 509 * WORK(KT1C),DUMMY) 510C 511 CALL CCLR_LAMTRA(WORK(KLAMDP),WORK(KLAMPC),WORK(KLAMDH), 512 * WORK(KLAMHC),WORK(KT1C),ISYMRC) 513C 514C------------------------------------------ 515C Calculate the F^C matrix 516C------------------------------------------ 517C 518 CALL CC3LR_MFOCK(WORK(KEND1),WORK(KT1C),WORK(KXIAJB),ISYFCKC) 519C 520C Put the F_{kc} part into correct F_{pq} 521C 522 CALL DZERO(WORK(KFCKC),N2BST(ISYFCKC)) 523C 524 DO ISYMC = 1, NSYM 525 ISYMK = MULD2H(ISYFCKC,ISYMC) 526 DO K = 1, NRHF(ISYMK) 527 DO C = 1, NVIR(ISYMC) 528 KOFF1 = KFCKC - 1 529 * + IFCVIR(ISYMK,ISYMC) 530 * + NORB(ISYMK)*(C - 1) 531 * + K 532 KOFF2 = KEND1 - 1 533 * + IT1AM(ISYMC,ISYMK) 534 * + NVIR(ISYMC)*(K - 1) 535 * + C 536C 537 WORK(KOFF1) = WORK(KOFF2) 538C 539 ENDDO 540 ENDDO 541 ENDDO 542C 543 IF (IPRINT .GT. 50) THEN 544 CALL AROUND( 'In CCSDT_FBMAT: Fock^C MO matrix' ) 545 CALL CC_PRFCKMO(WORK(KFCKC),ISYFCKC) 546 ENDIF 547C 548C 549C------------------------------------- 550C Summation over vectors! 551C------------------------------------- 552C 553 IVEC = 1 554 DO WHILE (IDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 555C 556 IDLSTL = IDOTS(IVEC,ITRAN) 557 ISYML = ILSTSYM(LISTL,IDLSTL) 558C 559C------------------------ 560C Symmetry check 561C------------------------ 562C 563 IF (ISYML .NE. ISYRES) THEN 564 WRITE(LUPRI,*)'Symmetry mismatch in CCSDT_FBMAT ' 565 WRITE(LUPRI,*)'ISYML = ', ISYML, 'not equal ISYRES =',ISYRES 566 CALL QUIT('Symmetry mismatch in CCSDT_FBMAT ') 567 END IF 568C 569C----------------------------- 570C Read occupied integrals. 571C----------------------------- 572C 573 KOMEG2 = KEND0 574 KRBJIA = KOMEG2 + NT2AM(ISYRES) 575 KTROC = KRBJIA + NT2SQ(ISYRES) 576 KTROC1 = KTROC + NTRAOC(ISINT1C) 577 KTROC0 = KTROC1 + NTRAOC(ISINT1C) 578 KTROC02= KTROC0 + NTRAOC(ISINT2) 579 KEND1 = KTROC02+ NTRAOC(ISINT2) 580 LWRK1 = LWORK - KEND1 581C 582 KTROC0Y = KEND1 583 KEND1 = KTROC0Y + (NCAU+1)*NTRAOC(ISINTR2) 584C 585 KINTOC = KEND1 586 KEND2 = KINTOC 587 * + MAX(NTOTOC(ISINTR2),NTOTOC(ISINT2),NTOTOC(ISINT1), 588 * NTOTOC(ISYMOP)) 589 LWRK2 = LWORK - KEND2 590C 591 IF (LWRK2 .LT. 0) THEN 592 WRITE(LUPRI,*) 'Memory available : ',LWORK 593 WRITE(LUPRI,*) 'Memory needed : ',KEND2 594 CALL QUIT('Insufficient space in CCSDT_FBMAT (intoc)') 595 END IF 596C 597C---------------------------------------- 598C Initialise the result vectors 599C---------------------------------------- 600C 601 CALL DZERO(WORK(KRBJIA),NT2SQ(ISYRES)) 602 CALL DZERO(WORK(KOMEG2),NT2AM(ISYRES)) 603C 604C------------------------ 605C Occupied integrals. 606C 607C Read in integrals for SMAT etc. 608C----------------------- 609C 610 IOFF = 1 611 IF (NTOTOC(ISINT2) .GT. 0) THEN 612 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 613 ENDIF 614C 615C---------------------------------- 616C Write out norms of Integrals. 617C---------------------------------- 618C 619 IF (IPRINT .GT. 55) THEN 620 XNORMVAL = DDOT(NTOTOC(ISINT2),WORK(KINTOC),1, 621 * WORK(KINTOC),1) 622 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XNORMVAL 623 ENDIF 624C 625C---------------------------------------------------------------------- 626C Transform (ai|j delta) integrals to (ai|j k) and sort as (ij,k,a) 627C---------------------------------------------------------------------- 628C 629 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0),WORK(KLAMDP), 630 * WORK(KEND2),LWRK2,ISINT2) 631C 632C---------------------------------------------------------------------- 633C (ai|j k) sorted as (ij,k,a) 634C---------------------------------------------------------------------- 635C 636 CALL CCFOP_SORT(WORK(KTROC0),WORK(KTROC02),ISINT2,1) 637C 638C------------------------------- 639C Write out norms of arrays. 640C------------------------------- 641C 642 IF (IPRINT .GT. 55) THEN 643 XNORMVAL = DDOT(NTOTOC(ISINT2),WORK(KINTOC),1, 644 * WORK(KINTOC),1) 645 WRITE(LUPRI,*) 'Norm of CKJDEL-INT ',XNORMVAL 646 ENDIF 647C 648 IF (IPRINT .GT. 55) THEN 649 XNORMVAL = DDOT(NTRAOC(ISINT2),WORK(KTROC0),1, 650 * WORK(KTROC0),1) 651 WRITE(LUPRI,*) 'Norm of TROC0 ',XNORMVAL 652 ENDIF 653C 654 IF (IPRINT .GT. 55) THEN 655 XNORMVAL = DDOT(NTRAOC(ISINT2),WORK(KTROC02),1, 656 * WORK(KTROC02),1) 657 WRITE(LUPRI,*) 'Norm of TROC02 ',XNORMVAL 658 ENDIF 659C 660C--------------------------------------------------------------------- 661C Read in T1B-transformed occupied integrals used to construct W3 662C--------------------------------------------------------------------- 663C 664 DO MCAU = 0,NCAU 665C 666 !Generate file name for T1B-transformed integrals 667 IF (LCAUCHY) THEN 668 !(FNDKBCR is not needed here) 669 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 670 ELSE 671 FNCKJDR = 'CCSDT_____FBMAT2' 672 END IF 673C 674 !Open the file for T1B-transformed integrals 675 LUCKJDR = -1 676 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 677C 678 !Get the offset for the integrals depending on MCAU 679 KOFF1 = KTROC0Y + MCAU*NTRAOC(ISINTR2) 680C 681 IOFF = 1 682 IF (NTOTOC(ISINTR2) .GT. 0) THEN 683 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF, 684 * NTOTOC(ISINTR2)) 685 ENDIF 686C 687 IF (IPRINT .GT. 55) THEN 688 XNORMVAL = DDOT(NTOTOC(ISINTR2),WORK(KINTOC),1, 689 * WORK(KINTOC),1) 690 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (B transformed) ', 691 * XNORMVAL 692 ENDIF 693C 694 CALL CC3_TROCC(WORK(KINTOC),WORK(KOFF1),WORK(KLAMDP), 695 * WORK(KEND2),LWRK2,ISINTR2) 696 !Close and keep the file for T1B-transformed occupied integrals 697 CALL WCLOSE2(LUCKJDR,FNCKJDR,'KEEP') !it was deleted before, but apparently it is needed (filip) 698C 699 END DO !MCAU 700C 701C------------------------------------------- 702C Occupied integrals used in contraction 703C------------------------------------------- 704C 705 IOFF = 1 706 IF (NTOTOC(ISYMOP) .GT. 0) THEN 707 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYMOP)) 708 ENDIF 709C 710C---------------------------------- 711C Write out norms of Integrals. 712C---------------------------------- 713C 714 IF (IPRINT .GT. 55) THEN 715 XNORMVAL = DDOT(NTOTOC(ISYMOP),WORK(KINTOC),1, 716 * WORK(KINTOC),1) 717 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XNORMVAL 718 ENDIF 719C 720C---------------------------------------------------------------------- 721C Transform (ia|j delta) integrals to (ia|j k) and sort as (i,j,k,a) 722C---------------------------------------------------------------------- 723C 724 CALL CCLR_TROCC(WORK(KINTOC),WORK(KTROC), 725 * WORK(KLAMHC),ISYMRC, 726 * WORK(KEND2),LWRK2) 727C 728C---------------------------------------------------------------------- 729C sort (i,j,k,a) as (a,i,j,k) 730C---------------------------------------------------------------------- 731C 732 CALL CCSDT_SRTOC2(WORK(KTROC),WORK(KTROC1),ISINT1C, 733 * WORK(KEND2),LWRK2) 734C 735C---------------------------- 736C General loop structure. 737C---------------------------- 738C 739 DO ISYMD = 1,NSYM 740C 741 ISYCKB = MULD2H(ISYMOP,ISYMD) 742 ISCKB2 = MULD2H(ISINT2,ISYMD) 743 ISCKB2Y = MULD2H(ISINTR2,ISYMD) 744 ISCKB1 = MULD2H(ISINT1,ISYMD) 745 ISCKB1C = MULD2H(ISINT1C,ISYMD) 746C 747 IF (IPRINT .GT. 55) THEN 748C 749 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISYCKB :',ISYCKB 750 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISCKB2 :',ISCKB2 751 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISCKB2Y:',ISCKB2Y 752 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISCKB1 :',ISCKB1 753 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISCKB1C:',ISCKB1C 754C 755 ENDIF 756 757C 758C-------------------------- 759C Memory allocation. 760C-------------------------- 761C 762 KTRVI = KEND1 763 KTRVI1 = KTRVI + NCKATR(ISCKB1C) 764 KTRVI3 = KTRVI1 + NCKATR(ISCKB1C) 765 KTRVI2 = KTRVI3 + NCKATR(ISCKB2) 766 KEND2 = KTRVI2 + NCKATR(ISCKB2) 767 LWRK2 = LWORK - KEND2 768C 769 KTRVI0 = KEND2 770 KEND3 = KTRVI0 + NCKATR(ISCKB2) 771 LWRK3 = LWORK - KEND3 772C 773 KTRVI0Y = KEND3 774 KEND3 = KTRVI0Y + (NCAU+1)*NCKATR(ISCKB2Y) 775 LWRK3 = LWORK - KEND3 776C 777 KINTVI = KEND3 778 KEND4 = KINTVI 779 * + MAX(NCKA(ISCKB2),NCKA(ISCKB1),NCKA(ISYCKB)) 780 LWRK4 = LWORK - KEND4 781C 782 IF (LWRK4 .LT. 0) THEN 783 WRITE(LUPRI,*) 'Memory available : ',LWORK 784 WRITE(LUPRI,*) 'Memory needed : ',KEND4 785 CALL QUIT('Insufficient space in CC3_XI') 786 END IF 787C 788C--------------------- 789C Sum over D 790C--------------------- 791C 792 DO D = 1,NVIR(ISYMD) 793C 794C------------------------------------------------------------ 795C Read and transform integrals used in contraction 796C with W intermediate. 797C------------------------------------------------------------ 798C 799 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 800 IF (NCKA(ISYCKB) .GT. 0) THEN 801 CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF, 802 & NCKA(ISYCKB)) 803 ENDIF 804 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVI), 805 * WORK(KLAMPC),ISYMRC, 806 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 807C 808 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 809 IF (NCKA(ISYCKB) .GT. 0) THEN 810 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 811 & NCKA(ISYCKB)) 812 ENDIF 813 CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVI1), 814 * WORK(KLAMPC),ISYMRC, 815 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 816C 817C----------------------------------------------- 818C Read virtual integrals used in first s3am. 819C----------------------------------------------- 820C 821 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 822 IF (NCKATR(ISCKB2) .GT. 0) THEN 823 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI0),IOFF, 824 & NCKATR(ISCKB2)) 825 ENDIF 826C 827C----------------------------------------------------------- 828C Sort the integrals for s3am 829C----------------------------------------------------------- 830C 831 CALL CCSDT_SRTVIR(WORK(KTRVI0),WORK(KTRVI2),WORK(KEND4), 832 * LWRK4,ISYMD,ISINT2) 833C 834 IF (IPRINT .GT. 55) THEN 835 XNORMVAL = DDOT(NCKATR(ISCKB2),WORK(KTRVI0),1, 836 * WORK(KTRVI0),1) 837 WRITE(LUPRI,*) 'Norm of TRVI0 ',XNORMVAL 838 ENDIF 839C 840 IF (IPRINT .GT. 55) THEN 841 XNORMVAL = DDOT(NCKATR(ISCKB2),WORK(KTRVI2),1, 842 * WORK(KTRVI2),1) 843 WRITE(LUPRI,*) 'Norm of TRVI2 ',XNORMVAL 844 ENDIF 845C 846C------------------------------------------------------------------------ 847C Read T1B-transformed virt ints (fixed D) used to construct W3 848C------------------------------------------------------------------------ 849C 850 DO MCAU = 0, NCAU 851C 852 !Generate file names for T1B-transformed integrals 853 IF (LCAUCHY) THEN 854 !(FNCKJDR is not needed here) 855 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 856 ELSE 857 FNDKBCR = 'CCSDT_____FBMAT4' 858 END IF 859C 860 !Open the files for T1B-transformed integrals 861 LUDKBCR = -1 862 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 863C 864 !Get the offset for a given MCAU 865 KOFF1 = KTRVI0Y + MCAU*NCKATR(ISCKB2Y) 866C 867 !Read in the integrals 868 IOFF = ICKBD(ISCKB2Y,ISYMD)+NCKATR(ISCKB2Y)*(D-1)+1 869 IF (NCKATR(ISCKB2) .GT. 0) THEN 870 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KOFF1),IOFF, 871 & NCKATR(ISCKB2Y)) 872 ENDIF 873C 874 !Close the file for T1B-transformed virtual integrals 875 !(and keep it: it will be needed in B-loop) 876 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 877C 878 END DO !MCAU 879C 880C------------------------------------------------------ 881C Read virtual integrals used in first U. 882C------------------------------------------------------ 883C 884 IOFF = ICKAD(ISCKB2,ISYMD) + NCKA(ISCKB2)*(D - 1) + 1 885 IF (NCKA(ISCKB2) .GT. 0) THEN 886 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 887 & NCKA(ISCKB2)) 888 ENDIF 889C 890 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI3),WORK(KLAMDH), 891 * ISYMD,D,ISINT2,WORK(KEND4),LWRK4) 892C 893C-------------------------------------------------------- 894C Sum over ISYMB 895C-------------------------------------------------------- 896C 897 DO ISYMB = 1,NSYM 898C 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 ISCKD2Y = MULD2H(ISINTR2,ISYMB) 905 ISWMAT = MULD2H(ISYMRB,ISCKIJ) 906C 907 IF (IPRINT .GT. 55) THEN 908C 909 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISYMD :',ISYMD 910 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISYMB :',ISYMB 911 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISYALJ :',ISYALJ 912 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISYALJ2:',ISYALJ2 913 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISYMBD :',ISYMBD 914 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISCKIJ :',ISCKIJ 915 WRITE(LUPRI,*) 'In CCSDT_FBMAT: ISWMAT :',ISWMAT 916C 917 ENDIF 918C 919C Can use kend3 since we do not need the integrals anymore. 920 KSMAT = KEND3 921 KSMAT3 = KSMAT + NCKIJ(ISCKIJ) 922 KUMAT = KSMAT3 + NCKIJ(ISCKIJ) 923 KUMAT3 = KUMAT + NCKIJ(ISCKIJ) 924 KDIAG = KUMAT3 + NCKIJ(ISCKIJ) 925 KDIAGW = KDIAG + NCKIJ(ISCKIJ) 926 KWMAT = KDIAGW + NCKIJ(ISWMAT) 927 KINDSQW = KWMAT + NCKIJ(ISWMAT) 928 KINDSQ = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 929 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 930 KINDEX2 = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 931 KTMAT = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1 932 KTRVI8 = KTMAT + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWMAT)) 933 KTRVI9 = KTRVI8 + NCKATR(ISCKD2) 934 KTRVI10 = KTRVI9 + NCKATR(ISCKD2) 935 KEND4 = KTRVI10 + NCKATR(ISCKD2) 936 LWRK4 = LWORK - KEND4 937C 938 KTRVI8Y = KEND4 939 KEND4 = KTRVI8Y + (NCAU+1)*NCKATR(ISCKD2Y) 940C 941C 942 KINTVI = KEND4 943 KEND5 = KINTVI + NCKA(ISCKD2) 944 LWRK5 = LWORK - KEND5 945C 946 IF (LWRK5 .LT. 0) THEN 947 WRITE(LUPRI,*) 'Memory available : ',LWORK 948 WRITE(LUPRI,*) 'Memory needed : ',KEND5 949 CALL QUIT('Insufficient space in CC3_XI') 950 END IF 951C 952C--------------------------------------------- 953C Construct part of the diagonal. 954C--------------------------------------------- 955C 956 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 957 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 958C 959 IF (IPRINT .GT. 55) THEN 960 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1, 961 * WORK(KDIAG),1) 962 WRITE(LUPRI,*) 'Norm of DIA ',XNORMVAL 963 ENDIF 964C 965C------------------------------------- 966C Construct index arrays. 967C------------------------------------- 968C 969 LENSQ = NCKIJ(ISCKIJ) 970 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 971 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 972 CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2) 973 LENSQW = NCKIJ(ISWMAT) 974 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 975C 976 DO B = 1,NVIR(ISYMB) 977C 978C--------------------------------------- 979C Initialise 980C--------------------------------------- 981C 982 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 983C 984C----------------------------------------------- 985C Read virtual integrals used in second s3am. 986C----------------------------------------------- 987C 988 IOFF = ICKBD(ISCKD2,ISYMB) + 989 & NCKATR(ISCKD2)*(B - 1) + 1 990 IF (NCKATR(ISCKD2) .GT. 0) THEN 991 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI8),IOFF, 992 & NCKATR(ISCKD2)) 993 ENDIF 994C 995C-------------------------------------------------------------------- 996C Read T1B-transformed vir int (fixed B) used for W3 997C-------------------------------------------------------------------- 998C 999 DO MCAU = 0, NCAU 1000C 1001 !Generate file names for T1B-transformed integrals 1002 IF (LCAUCHY) THEN 1003 !(FNCKJDR is not needed here) 1004 CALL GET_FILENAME(FNCKJDR,FNDKBCR,MCAU) 1005 ELSE 1006 FNDKBCR = 'CCSDT_____FBMAT4' 1007 END IF 1008C 1009 !Open the files for T1B-transformed integrals 1010 LUDKBCR = -1 1011 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 1012C 1013 !Get the offset for a given MCAU 1014 KOFF1 = KTRVI8Y + MCAU*NCKATR(ISCKD2Y) 1015C 1016 !Read in the integrals 1017 IOFF = ICKBD(ISCKD2Y,ISYMB) + 1018 & NCKATR(ISCKD2Y)*(B - 1) + 1 1019 IF (NCKATR(ISCKD2) .GT. 0) THEN 1020 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KOFF1),IOFF, 1021 & NCKATR(ISCKD2Y)) 1022 ENDIF 1023C 1024 !Close and keep the file for C1-transformed virt int 1025 !...or delte it if you don't need it anymore 1026 1027* IF ( (ISYMD.EQ.NSYM) .AND. (ISYMB.EQ.NSYM) 1028* * .AND. (D.EQ.NVIR(ISYMD)) 1029* * .AND. (B.EQ.NVIR(ISYMB))) THEN 1030* 1031* CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 1032* ELSE 1033* CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 1034* END IF 1035 1036 !It actually seems thath the file for C1-transformed virt int 1037 !is needed all the way through (filip) 1038 CALL WCLOSE2(LUDKBCR,FNDKBCR,'KEEP') 1039C 1040 END DO !MCAU 1041C 1042C----------------------------------------------------------- 1043C Sort the integrals for s3am 1044C----------------------------------------------------------- 1045C 1046 CALL CCSDT_SRTVIR(WORK(KTRVI8), 1047 * WORK(KTRVI9),WORK(KEND4), 1048 * LWRK4,ISYMB,ISINT2) 1049C 1050C 1051C---------------------------------------------------------- 1052C Read virtual integrals used in second U 1053C---------------------------------------------------------- 1054C 1055 IOFF = ICKAD(ISCKD2,ISYMB) 1056 * + NCKA(ISCKD2)*(B - 1) + 1 1057 IF (NCKA(ISCKD2) .GT. 0) THEN 1058 CALL GETWA2(LUDELD,FNDELD,WORK(KINTVI),IOFF, 1059 * NCKA(ISCKD2)) 1060 ENDIF 1061C 1062 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI10), 1063 * WORK(KLAMDH),ISYMB,B,ISINT2, 1064 & WORK(KEND5),LWRK5) 1065C 1066C------------------------------------------------------------------------ 1067C Calculate the S(ci,bk,dj) matrix for T3 for B,D. 1068C------------------------------------------------------------------- 1069C 1070 CALL CC3_SMAT(0.0D0,WORK(KT2TP),ISYMT2, 1071 * WORK(KTMAT),WORK(KTRVI0), 1072 * WORK(KTRVI2),WORK(KTROC0),ISINT2, 1073 * WORK(KFOCKD),WORK(KDIAG), 1074 * WORK(KSMAT),WORK(KEND4),LWRK4, 1075 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 1076 * ISYMB,B,ISYMD,D) 1077C 1078 CALL T3_FORBIDDEN(WORK(KSMAT),ISYMT3,ISYMB,B, 1079 * ISYMD,D) 1080C 1081* call sum_pt3(work(ksmat),isymb,b,isymd,d, 1082* * isckij,work(kx3am),1) 1083C 1084 IF (IPRINT .GT. 55) THEN 1085 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT),1, 1086 * WORK(KSMAT),1) 1087 WRITE(LUPRI,*) 'Norm of SMAT ',XNORMVAL 1088 ENDIF 1089C 1090C------------------------------------------------------------------- 1091C Calculate the S(ci,bk,dj) matrix for T3 for D,B. 1092C------------------------------------------------------------------- 1093C 1094 CALL CC3_SMAT(0.0D0,WORK(KT2TP),ISYMT2, 1095 * WORK(KTMAT),WORK(KTRVI8), 1096 * WORK(KTRVI9),WORK(KTROC0),ISINT2, 1097 * WORK(KFOCKD),WORK(KDIAG), 1098 * WORK(KSMAT3),WORK(KEND4),LWRK4, 1099 * WORK(KINDEX2),WORK(KINDSQ),LENSQ, 1100 * ISYMD,D,ISYMB,B) 1101C 1102 CALL T3_FORBIDDEN(WORK(KSMAT3),ISYMT3, 1103 * ISYMD,D,ISYMB,B) 1104C 1105* call sum_pt3(work(ksmat3),isymd,d,isymb,b, 1106* * isckij,work(kx3am),1) 1107C 1108 IF (IPRINT .GT. 55) THEN 1109 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT3),1, 1110 * WORK(KSMAT3),1) 1111 WRITE(LUPRI,*) 'Norm of SMAT3 ',XNORMVAL 1112 ENDIF 1113C 1114C--------------------------------------------------------------------------- 1115C Calculate U(ci,jk) for fixed b,d. 1116C-------------------------------------------------- 1117C 1118 CALL CC3_UMAT(0.0D0,WORK(KT2TP),ISYMT2, 1119 * WORK(KTRVI3),WORK(KTROC02), 1120 * ISINT2,WORK(KFOCKD), 1121 * WORK(KDIAG),WORK(KUMAT),WORK(KTMAT), 1122 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1123 * ISYMB,B,ISYMD,D) 1124C 1125 CALL T3_FORBIDDEN(WORK(KUMAT),ISYMT3, 1126 * ISYMB,B,ISYMD,D) 1127C 1128* call sum_pt3(work(kumat),isymb,b,isymd,d, 1129* * isckij,work(kx3am),3) 1130C 1131 IF (IPRINT .GT. 55) THEN 1132 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT),1, 1133 * WORK(KUMAT),1) 1134 WRITE(LUPRI,*) 'Norm of UMAT ',XNORMVAL 1135 ENDIF 1136C 1137C-------------------------------------------------- 1138C Calculate U(ci,jk) for fixed d,b. 1139C-------------------------------------------------- 1140C 1141 CALL CC3_UMAT(0.0D0,WORK(KT2TP),ISYMT2, 1142 * WORK(KTRVI10),WORK(KTROC02), 1143 * ISINT2,WORK(KFOCKD), 1144 * WORK(KDIAG),WORK(KUMAT3),WORK(KTMAT), 1145 * WORK(KEND4),LWRK4,WORK(KINDSQ),LENSQ, 1146 * ISYMD,D,ISYMB,B) 1147C 1148 CALL T3_FORBIDDEN(WORK(KUMAT3),ISYMT3, 1149 * ISYMD,D,ISYMB,B) 1150C 1151* call sum_pt3(work(kumat3),isymd,d,isymb,b, 1152* * isckij,work(kx3am),3) 1153C 1154 IF (IPRINT .GT. 55) THEN 1155 XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KUMAT),1, 1156 * WORK(KUMAT),1) 1157 WRITE(LUPRI,*) 'Norm of UMAT3 ',XNORMVAL 1158 ENDIF 1159C 1160C-------------------------------------------------- 1161C Sum up S and U intermediates to get T3 BD amplitudes 1162C-------------------------------------------------- 1163C 1164 CALL CC3_T3BD(ISCKIJ,WORK(KSMAT),WORK(KSMAT3), 1165 * WORK(KUMAT),WORK(KUMAT3),WORK(KTMAT), 1166 * WORK(KINDSQ),LENSQ) 1167C 1168C------------------------------------------------------ 1169C Based on T3 BD amplitudes calculate W BD intermediates 1170C------------------------------------------------------ 1171C 1172 IF ((LISTB(1:3).EQ.'R1 ').OR. 1173 * (LISTB(1:3).EQ.'RC ')) THEN 1174C 1175C------------------------------------------------------ 1176C Calculate the term <mu3|[Y,T3]|HF> virtual contribution 1177C added in W^BD(KWMAT) 1178C------------------------------------------------------ 1179C 1180 CALL WBD_V(WORK(KTMAT),ISCKIJ,WORK(KFCKB),ISYMRB, 1181 * WORK(KWMAT), 1182 * ISWMAT,WORK(KEND4),LWRK4) 1183C 1184C------------------------------------------------------ 1185C Calculate the term <mu3|[Y,T3]|HF> occupied contribution 1186C added in W^BD(KWMAT) 1187C------------------------------------------------------ 1188C 1189 CALL WBD_O(WORK(KTMAT),ISCKIJ,WORK(KFCKB),ISYMRB, 1190 * WORK(KWMAT), 1191 * ISWMAT,WORK(KEND4),LWRK4) 1192C 1193C------------------------------------------------------ 1194C Calculate the term <mu3|[[Y,T2],T2]|HF> 1195C added in W^BD(KWMAT) 1196C------------------------------------------------------ 1197C 1198 CALL WBD_T2(.FALSE.,B,ISYMB,D,ISYMD, 1199 * WORK(KT2TP),ISYMT2, 1200 * WORK(KT2TP),ISYMT2, 1201 * WORK(KFCKB),ISYMRB, 1202 * WORK(KINDSQW),LENSQW,WORK(KWMAT), 1203 * ISWMAT,WORK(KEND4),LWRK4) 1204C 1205 END IF 1206C 1207C------------------------------------------------------------------- 1208C Calculate the terms <mu3|[H^0,T2B]|HF> and <mu3|[H^B,T2^0]|HF> 1209C added in W^BD(KWMAT) 1210C------------------------------------------------------------------- 1211C 1212C 1213 DO MCAU = 0, NCAU !for non-Cauchy calc this is dummy loop 1214C 1215 !Get offset for KT2B vec for a given MCAU 1216 KOFF2 = KT2B + MCAU*NT2SQ(ISYMRB) 1217 1218 !Calculate the term <mu3|[H^0,T2B]|HF> 1219 CALL WBD_GROUND(WORK(KOFF2),ISYMRB,WORK(KTMAT), 1220 * WORK(KTRVI0),WORK(KTRVI8), 1221 * WORK(KTROC0),ISINT2, 1222 * WORK(KWMAT), 1223 * WORK(KEND4),LWRK4, 1224 * WORK(KINDSQW),LENSQW, 1225 * ISYMB,B,ISYMD,D) 1226C 1227 !Get offset for T1B-trans occ int for a given MCAU 1228 KOFFOCC = KTROC0Y + MCAU*NTRAOC(ISINTR2) 1229C 1230 !Get offset for T1B-trans virt int with fixed D for MCAU 1231 KOFFINTD = KTRVI0Y + MCAU*NCKATR(ISCKB2Y) 1232 !Get offset for T1B-trans virt int with fixed B for MCAU 1233 KOFFINTB = KTRVI8Y + MCAU*NCKATR(ISCKD2Y) 1234C 1235 !Calculate the term <mu3|[H^B,T2^0]|HF> 1236 CALL WBD_GROUND(WORK(KT2TP),ISYMT2,WORK(KTMAT), 1237 * WORK(KOFFINTD),WORK(KOFFINTB), 1238 * WORK(KOFFOCC),ISINTR2, 1239 * WORK(KWMAT), 1240 * WORK(KEND4),LWRK4, 1241 * WORK(KINDSQW),LENSQW, 1242 * ISYMB,B,ISYMD,D) 1243 1244c call sum_pt3(work(kwmat),isymb,b,isymd,d, 1245c * isckij,work(kx3am),4) 1246 1247 !Divide by the energy difference and 1248 !remove the forbidden elements 1249 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQB, 1250 * ISWMAT,WORK(KWMAT), 1251 * WORK(KDIAGW),WORK(KFOCKD)) 1252C 1253 CALL T3_FORBIDDEN(WORK(KWMAT),ISYMRB, 1254 * ISYMB,B,ISYMD,D) 1255C 1256 END DO !MCAU 1257C 1258 !AFTER ALL THE SIGN SHOULD NOT BE TURNED !!!!! 1259 1260* IF (LCAUCHY) THEN 1261* !trun the sign C_T <- -C_T 1262* CALL DSCAL(NCKIJ(ISWMAT),-1.0D0,WORK(KWMAT),1) 1263* END IF !LCAUCHY 1264 1265* call sum_pt3(work(kwmat),isymb,b,isymd,d, 1266* * isckij,work(kx3am),4) 1267 1268C 1269C---------------------------------------------------------------------- 1270C Calculate the term <mu2|[H,W^BD(3)]|HF> ( Fock matrix cont ) 1271C added in WORK(KOMEG2) 1272C---------------------------------------------------------------------- 1273C 1274 CALL CC3_CY2F(WORK(KOMEG2),ISYRES,WORK(KWMAT), 1275 * ISWMAT,WORK(KTMAT),WORK(KFCKC), 1276 * ISYFCKC,WORK(KINDSQW),LENSQW, 1277 * WORK(KEND4),LWRK4, 1278 * ISYMB,B,ISYMD,D) 1279C 1280C------------------------------------------------------ 1281C Calculate the term <mu2|[H,W^BD(3)]|HF> ( Occupied cont ) 1282C added in WORK(KOMEG2) 1283C------------------------------------------------------ 1284C 1285 CALL CC3_CY2O(WORK(KOMEG2),ISYRES,WORK(KWMAT), 1286 * ISWMAT,WORK(KTMAT),WORK(KTROC), 1287 * WORK(KTROC1),ISINT1C, 1288 * WORK(KEND4),LWRK4, 1289 * WORK(KINDSQW),LENSQW, 1290 * ISYMB,B,ISYMD,D) 1291C 1292C------------------------------------------------------ 1293C Calculate the term <mu2|[H,W^BD(3)]|HF> ( Virtual cont ) 1294C added in WORK(KOMEG2) 1295C------------------------------------------------------ 1296C 1297 CALL CC3_CY2V(WORK(KOMEG2),ISYRES,WORK(KRBJIA), 1298 * WORK(KWMAT),ISWMAT,WORK(KTMAT), 1299 * WORK(KTRVI),WORK(KTRVI1),ISINT1C, 1300 * WORK(KEND4),LWRK4,WORK(KINDSQW), 1301 * LENSQW,ISYMB,B,ISYMD,D) 1302C 1303 END DO ! B 1304 END DO ! ISYMB 1305C 1306 END DO ! D 1307 END DO ! ISYMD 1308C 1309* write(lupri,*) 'The summed W terms in ccsdt_fbmat, ncau = ',ncau 1310* call print_pt3(work(kx3am),1,4) 1311* write(lupri,*) 'The summed S terms : ' 1312* call print_pt3(work(kx3am),1,1) 1313C 1314C 1315C------------------------------------------------------ 1316C Accumulate RBJIA from <mu2|[H,W^BD(3)]|HF> ( Vccupied cont ) 1317C in WORK(KOMEG2) 1318C------------------------------------------------------ 1319C 1320 CALL CC3_RBJIA(WORK(KOMEG2),ISYRES,WORK(KRBJIA)) 1321C 1322 IF (IPRINT .GT. 55) THEN 1323 XNORMVAL = DDOT(NT2AM(ISYRES),WORK(KOMEG2),1,WORK(KOMEG2),1) 1324 WRITE(LUPRI,*) 'Norm of final WORK(KOMEG2) ',XNORMVAL 1325 ENDIF 1326C 1327C-------------------------------------------- 1328C Multiply the result vector with L2 1329C and put the result into DOTPROD 1330C-------------------------------------------- 1331C 1332 IOPT = 2 1333 CALL CC_RDRSP(LISTL,IDLSTL,ISYML,IOPT,MODEL,DUMMY,WORK(KEND1)) 1334 CALL CCLR_DIASCL(WORK(KEND1),HALF,ISYML) 1335C 1336 TCON = DDOT(NT2AM(ISYRES),WORK(KOMEG2),1,WORK(KEND1),1) 1337C 1338 DOTPROD(IVEC,ITRAN) = DOTPROD(IVEC,ITRAN) 1339 * + TCON 1340C 1341 IF (IPRINT .GT. 11) THEN 1342 write(lupri,*) 'LISTB, LISTC ', LISTB, LISTC 1343 WRITE(LUPRI,*) 'IVEC, ITRAN, Contribution: ',IVEC,ITRAN,TCON 1344 ENDIF 1345C 1346C------------------------------------- 1347C End loop over vectors. 1348C End loop over left and right 1349C------------------------------------- 1350C 1351 IVEC = IVEC + 1 1352C 1353 ENDDO ! DO WHILE 1354C 1355 ENDDO ! ITRAN 1356C 1357C-------------------------------- 1358C Close files 1359C-------------------------------- 1360C 1361 CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP') 1362 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 1363 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 1364 CALL WCLOSE2(LU3VI2,FN3VI2,'KEEP') 1365 CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP') 1366 CALL WCLOSE2(LUDELD,FNDELD,'KEEP') 1367C 1368C--------------------------------------------------------------------- 1369C It might have happened that 'CC3_CAUINT_V*' or 'CCSDT_____FBMAT4' 1370C files have not been deleted up there. Do it now! 1371C--------------------------------------------------------------------- 1372C 1373 IF (LCAUCHY) THEN 1374 CALL DELETE_FILES('CC3_CAUINT_V*') 1375 ELSE 1376 CALL DELETE_FILES('CCSDT_____FBMAT4') 1377 END IF 1378C 1379C------------- 1380C End 1381C------------- 1382C 1383 CALL QEXIT('FBMAT') 1384C 1385 RETURN 1386 END 1387