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_bmatsd */ 20 SUBROUTINE CC3_BMATSD(OMEGA1,OMEGA2,OMEGA1EFF,OMEGA2EFF,ISYRES, 21 * LISTB,IDLSTB,LISTC,IDLSTC, 22 * WORK,LWORK) 23* 24******************************************************************* 25* 26* This routine calculates the part of the right hand side vector 27* for the paritioned second-order amplitude equations. 28* 29* Calculate the B matrix contributions (using W intermmediate): 30* 31* W^BD = 1/2<mu3|[H^BC,T2^0]|HF> + <mu3|[H^B,T2^C]|HF> 32* 33* projected up into the Singles-and-Doubles space: 34* 35* omega1eff = omega1eff + <mu1|[H,W^BD]|HF> 36* omega2eff = omega2eff + <mu2|[H,W^BD]|HF> 37* 38* 39* At the end put omega to omegaeff: 40* 41* omega1eff = omega1eff + omega1 42* omega2eff = omega2eff + omega2 43* 44* 45* F. Pawlowski, P. Jorgensen, Aarhus Spring 2003. 46* 47******************************************************************* 48C 49 IMPLICIT NONE 50C 51#include "priunit.h" 52#include "ccr1rsp.h" 53#include "ccinftap.h" 54#include "iratdef.h" 55#include "ccsdsym.h" 56#include "ccsdinp.h" 57#include "dummy.h" 58#include "ccorb.h" 59#include "ccexci.h" 60C 61 CHARACTER*6 FN3VI 62 CHARACTER*8 FNTOC,FN3VI2 63C 64 PARAMETER (FNTOC = 'CCSDT_OC' ) 65 PARAMETER (FN3VI = 'CC3_VI' , FN3VI2 = 'CC3_VI12') 66 67 CHARACTER*12 FN3SRTR, FNCKJDR, FNDELDR, FNDKBCR 68 CHARACTER*13 FNCKJDR1, FNDKBCR1 69 CHARACTER*13 FNCKJDR2, FNDKBCR2 70C 71 PARAMETER(FN3SRTR = 'CCSDT_FBMAT1', FNCKJDR = 'CCSDT_FBMAT2', 72 * FNDELDR = 'CCSDT_FBMAT3', FNDKBCR = 'CCSDT_FBMAT4') 73C 74 PARAMETER(FNCKJDR1 = 'CCSDT_FBMAT21',FNDKBCR1 = 'CCSDT_FBMAT41') 75 PARAMETER(FNCKJDR2 = 'CCSDT_FBMAT22',FNDKBCR2 = 'CCSDT_FBMAT42') 76C 77 INTEGER LUTOC,LU3VI,LU3VI2 78 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR 79 INTEGER LUCKJDR1, LUDKBCR1 80 INTEGER LUCKJDR2, LUDKBCR2 81C 82 CHARACTER*3 LISTB, LISTC 83 CHARACTER*8 LABELB, LABELC 84 INTEGER IDLSTB,IDLSTC 85 INTEGER ISYMRB,ISYMRC 86C 87 CHARACTER CDUMMY*1 88C 89 INTEGER ISYRES 90 INTEGER LWORK 91C 92 INTEGER ISYM0,ISINT1,ISINT2,ISYMT2 93 INTEGER KT2TP,KFOCKD,KLAMDP,KLAMDH,KXIAJB,KFOCK0,KEND00,LWRK00 94 INTEGER KT1AM,KEND01,LWRK01 95 INTEGER IOPTTCME,IOPT 96 INTEGER KT2C,KT1B,KEND1,LWRK1 97 INTEGER ISINTR1B,ISINTR2B 98 INTEGER ISINTBC 99 INTEGER KRBJIA,KTROC,KTROC1,KTROC0Y,KTROC3,KINTOC 100 INTEGER MAXOCC,KEND2,LWRK2 101 INTEGER IOFF 102 INTEGER ISYMD,ISCKB1,ISCKB2Y,ISCKBBC 103 INTEGER KTRVI,KTRVI1,KTRVI0Y,KEND3,LWRK3,KTRVI7,KINTVI 104 INTEGER ISYMB,ISYMBC,ISYMBD,ISWMAT,ISCKD2Y,ISCKDBC 105 INTEGER KDIAGW,KWMAT,KTMAT,KEND4,KTRVI8Y,LWRK4,KTRVI9 106 INTEGER LENSQW,KINDSQW 107 INTEGER KEND0,LWRK0 108 INTEGER ISINTR1C,ISINTR2C 109 INTEGER KT2B,KT1C,KTROC0X 110 INTEGER ISCKB2X,KTRVI0X 111 INTEGER ISCKD2X,KTRVI8X 112 INTEGER MMAXOCC 113 INTEGER ILSTSYM 114c 115 integer kx3am 116c 117C 118#if defined (SYS_CRAY) 119 REAL OMEGA1(*),OMEGA1EFF(*),OMEGA2(*),OMEGA2EFF(*) 120 REAL WORK(LWORK) 121 REAL FREQB,FREQC,FREQBC 122 REAL XNORMVAL,DDOT 123#else 124 DOUBLE PRECISION OMEGA1(*),OMEGA1EFF(*),OMEGA2(*),OMEGA2EFF(*) 125 DOUBLE PRECISION WORK(LWORK) 126 DOUBLE PRECISION FREQB,FREQC,FREQBC 127 DOUBLE PRECISION XNORMVAL,DDOT 128#endif 129C 130 131 CALL QENTER('CC3_BMATSD') 132C 133C---------------------- 134C Lists check 135C---------------------- 136C 137 IF ( (LISTB(1:3).EQ.'R1 ') .AND. (LISTC(1:3).EQ.'R1 ') ) THEN 138 CONTINUE 139 ELSE IF ( (LISTB(1:3).EQ.'R1 ') .AND. (LISTC(1:3).EQ.'RE ') ) THEN 140 CONTINUE 141 ELSE 142 WRITE(LUPRI,*)'LISTB = ', LISTB 143 WRITE(LUPRI,*)'LISTC = ', LISTC 144 WRITE(LUPRI,*)'Only implemented when ' 145 * //'LISTB = R1 and LISTC = R1 or RE ' 146 CALL QUIT('Case not implemented in CC3_BMATSD') 147 END IF 148C 149C-------------------- 150C Lists handling 151C-------------------- 152C 153 IF (LISTB(1:3).EQ.'R1 ') THEN 154 ISYMRB = ISYLRT(IDLSTB) 155 FREQB = FRQLRT(IDLSTB) 156 LABELB = LRTLBL(IDLSTB) 157 ELSE 158 CALL QUIT('Unknown list for LISTB in CC3_BMATSD') 159 END IF 160C 161 IF (LISTC(1:3).EQ.'R1 ') THEN 162 ISYMRC = ISYLRT(IDLSTC) 163 FREQC = FRQLRT(IDLSTC) 164 LABELC = LRTLBL(IDLSTC) 165 ELSE IF (LISTC(1:3).EQ.'RE ') THEN 166 ISYMRC = ILSTSYM(LISTC,IDLSTC) 167 FREQC = EIGVAL(IDLSTC) 168 LABELC = '- none -' 169 ELSE 170 CALL QUIT('Unknown list for LISTC in CC3_BMATSD') 171 END IF 172C 173C------------------------------------------- 174C Construct FREQBC = (omega_B + omega_C) 175C------------------------------------------- 176C 177 FREQBC = FREQB + FREQC 178C 179C-------------------------- 180C Save and set flags. 181C-------------------------- 182C 183 ISYM0 = 1 184 ISINT1 = 1 185 ISINT2 = 1 186 ISYMT2 = 1 187C 188C------------------------------------------ 189C Open files (integrals in contraction) 190C------------------------------------------ 191C 192 LUTOC = -1 193 LU3VI = -1 194 LU3VI2 = -1 195C 196 CALL WOPEN2(LUTOC,FNTOC,64,0) 197 CALL WOPEN2(LU3VI,FN3VI,64,0) 198 CALL WOPEN2(LU3VI2,FN3VI2,64,0) 199C 200C----------------------------------------- 201C Open temporary files (H^B integrals) 202C----------------------------------------- 203C 204 LU3SRTR = -1 205 LUCKJDR = -1 206 LUDELDR = -1 207 LUDKBCR = -1 208C 209 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 210 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 211 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 212 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 213C 214C------------------------------------------ 215C Open temporary files (H^BC integrals) 216C------------------------------------------ 217C 218 LUCKJDR1 = -1 219 LUDKBCR1 = -1 220C 221 CALL WOPEN2(LUCKJDR1,FNCKJDR1,64,0) 222 CALL WOPEN2(LUDKBCR1,FNDKBCR1,64,0) 223C 224C------------------------------------------ 225C Open temporary files (H^C integrals) 226C------------------------------------------ 227C 228 LUCKJDR2 = -1 229 LUDKBCR2 = -1 230C 231 CALL WOPEN2(LUCKJDR2,FNCKJDR2,64,0) 232 CALL WOPEN2(LUDKBCR2,FNDKBCR2,64,0) 233C 234C---------------------------------------------- 235C Calculate the zeroth order stuff once 236C---------------------------------------------- 237C 238 KT2TP = 1 239 KFOCKD = KT2TP + NT2SQ(ISYMT2) 240 KLAMDP = KFOCKD + NORBTS 241 KLAMDH = KLAMDP + NLAMDT 242 KXIAJB = KLAMDH + NLAMDT 243 KFOCK0 = KXIAJB + NT2AM(ISINT1) 244 KEND00 = KFOCK0 + N2BST(ISYM0) 245 LWRK00 = LWORK - KEND00 246C 247 KT1AM = KEND00 248 KEND01 = KT1AM + NT1AM(ISYMT2) 249 LWRK01 = LWORK - KEND01 250C 251 IF (LWRK01 .LT. 0) THEN 252 WRITE(LUPRI,*)'Memory available: ',LWORK 253 WRITE(LUPRI,*)'Memory needed: ',KEND01 254 CALL QUIT('Out of memory in CC3_BMATSD (zeroth allo.)') 255 ENDIF 256C 257C------------------------ 258C Construct L(ia,jb). 259C------------------------ 260C 261 REWIND(LUIAJB) 262 CALL READI(LUIAJB,IRAT*NT2AM(ISINT1),WORK(KXIAJB)) 263 IOPTTCME = 1 264 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISINT1,IOPTTCME) 265C 266 IF ( IPRINT .GT. 55) THEN 267 XNORMVAL = DDOT(NT2AM(ISINT1),WORK(KXIAJB),1, 268 * WORK(KXIAJB),1) 269 WRITE(LUPRI,*) 'Norm of IAJB ',XNORMVAL 270 ENDIF 271C 272C---------------------------------------------------------------- 273C Read t1 and calculate the zero'th order Lambda matrices 274C---------------------------------------------------------------- 275C 276 CALL GET_LAMBDA0(WORK(KLAMDP),WORK(KLAMDH),WORK(KEND01),LWRK01) 277C 278C------------------------------------------- 279C Read in t2 , square it and reorder 280C------------------------------------------- 281C 282 IOPT = 2 283 CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KT2TP),'R0',0,ISYMT2, 284 * WORK(KEND01),LWRK01) 285 286 IF (IPRINT .GT. 55) THEN 287 XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 288 WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL 289 ENDIF 290C 291C------------------------------------------------------------------ 292C Read in Fock matrix in AO basis (from file) and transform to 293C Lambda_0 basis. 294C------------------------------------------------------------------ 295C 296 CALL GET_FOCK0(WORK(KFOCK0),WORK(KLAMDP),WORK(KLAMDH), 297 * WORK(KEND01),LWRK01) 298C 299C-------------------------------------- 300C Read in orbital energies 301C-------------------------------------- 302C 303 CALL GET_ORBEN(WORK(KFOCKD),WORK(KEND01),LWRK01) 304c 305c If we want to sum the T3 amplitudes 306c 307 if (.false.) then 308 kx3am = kend00 309 kend00 = kx3am + nt1amx*nt1amx*nt1amx 310 call dzero(work(kx3am),nt1amx*nt1amx*nt1amx) 311 lwrk00 = lwork - kend00 312 if (lwrk00 .lt. 0) then 313 write(lupri,*) 'Memory available : ',lwork 314 write(lupri,*) 'Memory needed : ',kend00 315 call quit('Insufficient space (T3) in CC3_BMATSD (1)') 316 END IF 317 endif 318C 319C------------------------------------------------- 320C Read T1^B and T2^C 321C Calculate (ck|de)-Btilde and (ck|lm)-Btilde 322C Calculate (ck|de)-BCtilde and (ck|lm)-BCtilde 323C Used to construct WBD intermmediate 324C------------------------------------------------- 325C 326 KT2C = KEND00 327 KEND0 = KT2C + NT2SQ(ISYMRC) 328 LWRK0 = LWORK - KEND0 329C 330 KT2B = KEND0 331 KEND0 = KT2B + NT2SQ(ISYMRB) 332 LWRK0 = LWORK - KEND0 333C 334 KT1C = KEND0 335 KEND0 = KT1C + NT1AM(ISYMRC) 336 LWRK0 = LWORK - KEND0 337C 338 KT1B = KEND0 339 KEND1 = KT1B + NT1AM(ISYMRB) 340 LWRK1 = LWORK - KEND1 341C 342 IF (LWRK1 .LT. NT2SQ(ISYMRB)) THEN 343 CALL QUIT('Out of memory in CC3_BMATSD (TOT_T3Y) ') 344 ENDIF 345C 346C-------------------------- 347C Read in T1^B and T2^B 348C-------------------------- 349C 350 IOPT = 3 351 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1B),WORK(KT2B),LISTB,IDLSTB, 352 * ISYMRB,WORK(KEND1),LWRK1) 353C 354 IF (IPRINT .GT. 55) THEN 355 XNORMVAL = DDOT(NT1AM(ISYMRB),WORK(KT1B),1,WORK(KT1B),1) 356 WRITE(LUPRI,*) 'Norm of T1B ',XNORMVAL 357 XNORMVAL = DDOT(NT2SQ(ISYMRB),WORK(KT2B),1,WORK(KT2B),1) 358 WRITE(LUPRI,*) 'Norm of T2B ',XNORMVAL 359 ENDIF 360C 361C-------------------------- 362C Read in T1^C and T2^C 363C-------------------------- 364C 365 IOPT = 3 366 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1C),WORK(KT2C),LISTC,IDLSTC, 367 * ISYMRC,WORK(KEND1),LWRK1) 368C 369 IF (IPRINT .GT. 55) THEN 370 XNORMVAL = DDOT(NT1AM(ISYMRC),WORK(KT1C),1,WORK(KT1C),1) 371 WRITE(LUPRI,*) 'Norm of T1C ',XNORMVAL 372 XNORMVAL = DDOT(NT2SQ(ISYMRC),WORK(KT2C),1,WORK(KT2C),1) 373 WRITE(LUPRI,*) 'Norm of T2C ',XNORMVAL 374 ENDIF 375C 376C---------------------------------------------------- 377C Integrals (ck|de)-tilde(B) and (ck|lm)-tilde(B) 378C---------------------------------------------------- 379C 380 ISINTR1B = MULD2H(ISINT1,ISYMRB) 381 ISINTR2B = MULD2H(ISINT2,ISYMRB) 382C 383 CALL CC3_BARINT(WORK(KT1B),ISYMRB,WORK(KLAMDP), 384 * WORK(KLAMDH),WORK(KEND1),LWRK1, 385 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 386C 387 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTR1B,LU3SRTR,FN3SRTR, 388 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 389 * IDUMMY,CDUMMY) 390C 391 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTR1B, 392 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 393C 394C-------------------------------- 395C Re-use some temporary files 396C-------------------------------- 397C 398 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 399 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 400C 401 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 402 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 403C 404C------------------------------------------------------ 405C Calculate the (ck|de)-BCtilde and (ck|lm)-BCtilde 406C------------------------------------------------------ 407C 408 ISYMBC = MULD2H(ISYMRB,ISYMRC) 409C 410 CALL CC3_3BARINT(ISYMRB,LISTB,IDLSTB,ISYMRC,LISTC,IDLSTC, 411 * IDUMMY,CDUMMY,IDUMMY,.FALSE., 412 * WORK(KLAMDP),WORK(KLAMDH),WORK(KEND1),LWRK1, 413 * LU3SRTR,FN3SRTR,LUCKJDR1,FNCKJDR1) 414C 415 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISYMBC,LU3SRTR,FN3SRTR, 416 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 417 * IDUMMY,CDUMMY) 418C 419 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISYMBC, 420 * LUDELDR,FNDELDR,LUDKBCR1,FNDKBCR1) 421C 422C-------------------------------- 423C Re-use some temporary files 424C-------------------------------- 425C 426 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 427 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 428C 429 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 430 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 431C 432C---------------------------------------------------- 433C Integrals (ck|de)-tilde(C) and (ck|lm)-tilde(C) 434C---------------------------------------------------- 435C 436 ISINTR1C = MULD2H(ISINT1,ISYMRC) 437 ISINTR2C = MULD2H(ISINT2,ISYMRC) 438C 439 CALL CC3_BARINT(WORK(KT1C),ISYMRC,WORK(KLAMDP), 440 * WORK(KLAMDH),WORK(KEND1),LWRK1, 441 * LU3SRTR,FN3SRTR,LUCKJDR2,FNCKJDR2) 442C 443 CALL CC3_SORT1(WORK(KEND1),LWRK1,2,ISINTR1C,LU3SRTR,FN3SRTR, 444 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 445 * IDUMMY,CDUMMY) 446C 447 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND1),LWRK1,ISINTR1C, 448 * LUDELDR,FNDELDR,LUDKBCR2,FNDKBCR2) 449C 450C-------------------------------- 451C Read occupied integrals 452C------------------------------- 453C 454 ISINTBC = MULD2H(ISINT2,ISYMBC) 455C 456 !Use KEND0, because KT1B is not needed any more 457 KRBJIA = KEND0 458 KTROC = KRBJIA + NT2SQ(ISYRES) 459 KTROC1 = KTROC + NTRAOC(ISINT1) !KTROC - int. in contraction 460 KEND1 = KTROC1 + NTRAOC(ISINT1) !KTROC1 - int. in contraction 461 LWRK1 = LWORK - KEND1 462C 463 KTROC0Y = KEND1 464 KEND1 = KTROC0Y + NTRAOC(ISINTR2B)!KTROC0Y - int. in <mu3|[H^B,T2^C]|HF> 465C 466 KTROC0X = KEND1 467 KEND1 = KTROC0X + NTRAOC(ISINTR2C)!KTROC0X - int. in <mu3|[H^C,T2^B]|HF> 468C === 469 KTROC3 = KEND1 470 KEND1 = KTROC3 + NTRAOC(ISINTBC)!KTROC3 - int. in <mu3|[H^BC,T2]|HF> 471C === 472 KINTOC = KEND1 473 MAXOCC = MAX(NTOTOC(ISINTR2B),NTOTOC(ISINTBC)) 474 MMAXOCC = MAX(NTOTOC(ISINTR2C),MAXOCC) 475 KEND2 = KINTOC + MAX(NTOTOC(ISINT1),MMAXOCC)!KINTOC - temporary storage 476 LWRK2 = LWORK - KEND2 477C 478 IF (LWRK2 .LT. 0) THEN 479 WRITE(LUPRI,*) 'Memory available : ',LWORK 480 WRITE(LUPRI,*) 'Memory needed : ',KEND2 481 CALL QUIT('Insufficient space in CC3_BMATSD (2)') 482 END IF 483C 484 CALL DZERO(WORK(KRBJIA),NT2SQ(ISYRES)) 485C 486C 487C------------------------------------------------------------- 488C B-transformed occupied integrals for <mu3|[H^B,T2^C]|HF> 489C------------------------------------------------------------- 490C 491 IOFF = 1 492 IF (NTOTOC(ISINTR2B) .GT. 0) THEN 493 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF, 494 * NTOTOC(ISINTR2B)) 495 ENDIF 496C 497 IF (IPRINT .GT. 55) THEN 498 XNORMVAL = DDOT(NTOTOC(ISINTR2B),WORK(KINTOC),1, 499 * WORK(KINTOC),1) 500 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (B transformed) ', 501 * XNORMVAL 502 ENDIF 503C 504 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0Y),WORK(KLAMDP), 505 * WORK(KEND2),LWRK2,ISINTR2B) 506C 507C------------------------------------------------------------- 508C C-transformed occupied integrals for <mu3|[H^C,T2^B]|HF> 509C------------------------------------------------------------- 510C 511 IOFF = 1 512 IF (NTOTOC(ISINTR2C) .GT. 0) THEN 513 CALL GETWA2(LUCKJDR2,FNCKJDR2,WORK(KINTOC),IOFF, 514 * NTOTOC(ISINTR2C)) 515 ENDIF 516C 517 IF (IPRINT .GT. 55) THEN 518 XNORMVAL = DDOT(NTOTOC(ISINTR2C),WORK(KINTOC),1, 519 * WORK(KINTOC),1) 520 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT (C transformed) ', 521 * XNORMVAL 522 ENDIF 523C 524 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0X),WORK(KLAMDP), 525 * WORK(KEND2),LWRK2,ISINTR2C) 526C 527C-------------------------------------------------------------------------- 528C BC-transformed occupied integrals for <mu3|[H^BC,T2]|HF> 529C-------------------------------------------------------------------------- 530C 531 IOFF = 1 532 IF (NTOTOC(ISINTBC) .GT. 0) THEN 533 CALL GETWA2(LUCKJDR1,FNCKJDR1,WORK(KINTOC),IOFF, 534 * NTOTOC(ISINTBC)) 535 ENDIF 536C 537 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC3),WORK(KLAMDP), 538 * WORK(KEND2),LWRK2,ISINTBC) 539C 540C----------------------------------- 541C Occupied integrals in contraction 542C----------------------------------- 543C 544 IOFF = 1 545 IF (NTOTOC(ISINT1) .GT. 0) THEN 546 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT1)) 547 ENDIF 548C 549 !Write out norms of integrals. 550 IF (IPRINT .GT. 55) THEN 551 XNORMVAL = DDOT(NTOTOC(ISINT1),WORK(KINTOC),1, 552 * WORK(KINTOC),1) 553 WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XNORMVAL 554 ENDIF 555C 556 !Transform (ia|j delta) integrals to (ia|j k) and sort as (i,j,k,a) 557 CALL CCSDT_TROCC(WORK(KINTOC),WORK(KTROC),WORK(KLAMDH), 558 * WORK(KEND2),LWRK2) 559C 560 !sort (i,j,k,a) as (a,i,j,k) 561 CALL CCSDT_SRTOC2(WORK(KTROC),WORK(KTROC1),ISINT1, 562 * WORK(KEND2),LWRK2) 563C 564C--------------------- 565C Start ISYMD-loop 566C--------------------- 567C 568 DO ISYMD = 1,NSYM 569C 570 ISCKB1 = MULD2H(ISINT1,ISYMD) 571 ISCKB2Y = MULD2H(ISINTR2B,ISYMD) 572 ISCKB2X = MULD2H(ISINTR2C,ISYMD) 573 ISCKBBC = MULD2H(ISINTBC,ISYMD) 574C 575C---------------------------------------- 576C Read virtual integrals (fixed D) 577C---------------------------------------- 578C 579 !Use KEND1, because KINTOC is not needed any more 580 KTRVI = KEND1 581 KTRVI1 = KTRVI + NCKATR(ISCKB1) !KTRVI - int. in contraction 582 KEND2 = KTRVI1 + NCKATR(ISCKB1) !KTRVI1 - int. in contraction 583 LWRK2 = LWORK - KEND2 584C 585 KTRVI0Y = KEND2 586 KEND3 = KTRVI0Y + NCKATR(ISCKB2Y)!KTRVI0Y - int. in <mu3|[H^B,T2^C]|HF> 587 LWRK3 = LWORK - KEND3 ! === 588C 589 KTRVI0X = KEND3 590 KEND3 = KTRVI0X + NCKATR(ISCKB2X)!KTRVI0X - int. in <mu3|[H^C,T2^B]|HF> 591 LWRK3 = LWORK - KEND3 ! === 592C 593 KTRVI7 = KEND3 594 KEND3 = KTRVI7 + NCKATR(ISCKBBC)!KTRVI7 - int. in <mu3|[H^BC,T2^0]|HF> 595 LWRK3 = LWORK - KEND3 ! ==== 596C 597 KINTVI = KEND3 598 KEND4 = KINTVI + NCKA(ISCKB1)!KINTVI - temporary storage 599 LWRK4 = LWORK - KEND4 600C 601 IF (LWRK4 .LT. 0) THEN 602 WRITE(LUPRI,*) 'Memory available : ',LWORK 603 WRITE(LUPRI,*) 'Memory needed : ',KEND4 604 CALL QUIT('Insufficient space in CC3_BMATSD (3)') 605 END IF 606C 607C-------------------- 608C Start D-loop 609C-------------------- 610C 611 DO D = 1,NVIR(ISYMD) 612C 613C---------------------------------------------------------------------------- 614C B-transformed virtual integrals for <mu3|[H^B,T2^C]|HF> (fixed D) 615C---------------------------------------------------------------------------- 616C 617 IOFF = ICKBD(ISCKB2Y,ISYMD) + NCKATR(ISCKB2Y)*(D - 1) + 1 618 IF (NCKATR(ISCKB2Y) .GT. 0) THEN 619 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI0Y),IOFF, 620 & NCKATR(ISCKB2Y)) 621 ENDIF 622C 623C---------------------------------------------------------------------------- 624C C-transformed virtual integrals for <mu3|[H^C,T2^B]|HF> (fixed D) 625C---------------------------------------------------------------------------- 626C 627 IOFF = ICKBD(ISCKB2X,ISYMD) + NCKATR(ISCKB2X)*(D - 1) + 1 628 IF (NCKATR(ISCKB2X) .GT. 0) THEN 629 CALL GETWA2(LUDKBCR2,FNDKBCR2,WORK(KTRVI0X),IOFF, 630 & NCKATR(ISCKB2X)) 631 ENDIF 632C 633C----------------------------------------------------------------------------- 634C B-transformed virtual integrals for <mu3|[H^BC,T2^0]|HF> (fixed D) 635C----------------------------------------------------------------------------- 636C 637 IOFF = ICKBD(ISCKBBC,ISYMD) + NCKATR(ISCKBBC)*(D - 1) + 1 638 IF (NCKATR(ISCKBBC) .GT. 0) THEN 639 CALL GETWA2(LUDKBCR1,FNDKBCR1,WORK(KTRVI7),IOFF, 640 & NCKATR(ISCKBBC)) 641 ENDIF 642C 643C----------------------------------------------------- 644C Virtual integrals in contraction (fixed D) 645C----------------------------------------------------- 646C 647 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 648 IF (NCKA(ISCKB1) .GT. 0) THEN 649 CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF, 650 * NCKA(ISCKB1)) 651 ENDIF 652 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI),WORK(KLAMDP), 653 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 654C 655 IOFF = ICKAD(ISCKB1,ISYMD) + NCKA(ISCKB1)*(D - 1) + 1 656 IF (NCKA(ISCKB1) .GT. 0) THEN 657 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 658 * NCKA(ISCKB1)) 659 ENDIF 660 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI1),WORK(KLAMDP), 661 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 662C 663C--------------------------- 664C Start ISYMB-loop 665C--------------------------- 666C 667 DO ISYMB = 1,NSYM 668C 669 ISYMBC = MULD2H(ISYMRB,ISYMRC) 670 ISYMBD = MULD2H(ISYMB,ISYMD) 671 ISWMAT = MULD2H(ISYMBC,ISYMBD) 672 ISCKD2Y = MULD2H(ISINTR2B,ISYMB) 673 ISCKD2X = MULD2H(ISINTR2C,ISYMB) 674 ISCKDBC = MULD2H(ISINTBC,ISYMB) 675C 676 !Use KEND3, because KINTVI is not needed any more 677 KDIAGW = KEND3 678 KWMAT = KDIAGW + NCKIJ(ISWMAT) 679 KINDSQW = KWMAT + NCKIJ(ISWMAT) 680 KTMAT = KINDSQW + (6*NCKIJ(ISWMAT) - 1)/IRAT + 1 681 KEND4 = KTMAT + NCKIJ(ISWMAT) 682C 683 KTRVI8Y = KEND4 684 KEND4 = KTRVI8Y + NCKATR(ISCKD2Y)!KTRVI8Y: <mu3|[H^B,T2^C]|HF> 685 LWRK4 = LWORK - KEND4 ! === 686C 687 KTRVI8X = KEND4 688 KEND4 = KTRVI8X + NCKATR(ISCKD2X)!KTRVI8X: <mu3|[H^C,T2^B]|HF> 689 LWRK4 = LWORK - KEND4 ! === 690C 691 KTRVI9 = KEND4 692 KEND4 = KTRVI9 + NCKATR(ISCKDBC)!KTRVI9: <mu3|[H^BC,T2^0]|HF> 693 LWRK4 = LWORK - KEND4 ! ==== 694C 695 IF (LWRK4 .LT. 0) THEN 696 WRITE(LUPRI,*) 'Memory available : ',LWORK 697 WRITE(LUPRI,*) 'Memory needed : ',KEND4 698 CALL QUIT('Insufficient space in CC3_BMATSD (4)') 699 END IF 700C 701C--------------------------------------------- 702C Construct part of the diagonal. 703C--------------------------------------------- 704C 705 CALL CC3_DIAG(WORK(KDIAGW),WORK(KFOCKD),ISWMAT) 706C 707 IF (IPRINT .GT. 55) THEN 708 XNORMVAL = DDOT(NCKIJ(ISWMAT),WORK(KDIAGW),1, 709 * WORK(KDIAGW),1) 710 WRITE(LUPRI,*) 'Norm of DIA ',XNORMVAL 711 ENDIF 712C 713C------------------------------------- 714C Construct index arrays. 715C------------------------------------- 716C 717 LENSQW = NCKIJ(ISWMAT) 718 CALL CC3_INDSQ(WORK(KINDSQW),LENSQW,ISWMAT) 719C 720C-------------------------- 721C Start B-loop 722C-------------------------- 723C 724 DO B = 1,NVIR(ISYMB) 725C 726C--------------------------------------- 727C Initialise 728C--------------------------------------- 729C 730 CALL DZERO(WORK(KWMAT),NCKIJ(ISWMAT)) 731C 732C---------------------------------------------------------------------------- 733C B-transformed virtual integrals for <mu3|[H^B,T2^C]|HF> 734C (fixed B) 735C---------------------------------------------------------------------------- 736C 737 IOFF = ICKBD(ISCKD2Y,ISYMB) + NCKATR(ISCKD2Y)*(B-1) +1 738 IF (NCKATR(ISCKD2Y) .GT. 0) THEN 739 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI8Y),IOFF, 740 * NCKATR(ISCKD2Y)) 741 ENDIF 742C 743C---------------------------------------------------------------------------- 744C C-transformed virtual integrals for <mu3|[H^C,T2^B]|HF> 745C (fixed B) 746C---------------------------------------------------------------------------- 747C 748 IOFF = ICKBD(ISCKD2X,ISYMB) + NCKATR(ISCKD2X)*(B-1) +1 749 IF (NCKATR(ISCKD2X) .GT. 0) THEN 750 CALL GETWA2(LUDKBCR2,FNDKBCR2,WORK(KTRVI8X),IOFF, 751 * NCKATR(ISCKD2X)) 752 ENDIF 753C 754C----------------------------------------------------------------------------- 755C B-transformed virtual integrals for <mu3|[H^BC,T2^0]|HF> 756C (fixed B) 757C----------------------------------------------------------------------------- 758C 759 IOFF = ICKBD(ISCKDBC,ISYMB) + NCKATR(ISCKDBC)*(B-1) +1 760 IF (NCKATR(ISCKDBC) .GT. 0) THEN 761 CALL GETWA2(LUDKBCR1,FNDKBCR1,WORK(KTRVI9),IOFF, 762 * NCKATR(ISCKDBC)) 763 ENDIF 764 765C 766C---------------------------------------------- 767C Calculate <mu3|[H^C,T2^B]|HF> 768C---------------------------------------------- 769C 770 CALL WBD_GROUND(WORK(KT2B),ISYMRB,WORK(KTMAT), 771 * WORK(KTRVI0X),WORK(KTRVI8X), 772 * WORK(KTROC0X),ISINTR2C,WORK(KWMAT), 773 * WORK(KEND4),LWRK4, 774 * WORK(KINDSQW),LENSQW, 775 * ISYMB,B,ISYMD,D) 776C 777C---------------------------------------------- 778C Calculate <mu3|[H^B,T2^C]|HF> 779C---------------------------------------------- 780C 781 CALL WBD_GROUND(WORK(KT2C),ISYMRC,WORK(KTMAT), 782 * WORK(KTRVI0Y),WORK(KTRVI8Y), 783 * WORK(KTROC0Y),ISINTR2B,WORK(KWMAT), 784 * WORK(KEND4),LWRK4, 785 * WORK(KINDSQW),LENSQW, 786 * ISYMB,B,ISYMD,D) 787C 788C---------------------------------------------- 789C Calculate <mu3|[H^BC,T2^0]|HF> 790C---------------------------------------------- 791C 792 CALL WBD_GROUND(WORK(KT2TP),ISYMT2,WORK(KTMAT), 793 * WORK(KTRVI7),WORK(KTRVI9), 794 * WORK(KTROC3),ISINTBC,WORK(KWMAT), 795 * WORK(KEND4),LWRK4, 796 * WORK(KINDSQW),LENSQW, 797 * ISYMB,B,ISYMD,D) 798C 799C---------------------------------------------------- 800C Divide by orbital energy difference 801C and remove the forbidden elements 802C---------------------------------------------------- 803C 804C 805 806 CALL T3_FORBIDDEN(WORK(KWMAT),ISYMBC, 807 * ISYMB,B,ISYMD,D) 808 809c call sum_pt3(work(kwmat),isymb,b,isymd,d, 810c * iswmat,work(kx3am),4) 811c write(lupri,*) 'Total norm of WBD : ', 812c * ddot(nckij(iswmat),work(kwmat),1,work(kwmat),1) 813 814 CALL WBD_DIA(B,ISYMB,D,ISYMD,FREQBC, 815 * ISWMAT,WORK(KWMAT), 816 * WORK(KDIAGW),WORK(KFOCKD)) 817 818C 819C----------------------------------------------- 820C Carry out the contractions... 821C----------------------------------------------- 822C 823C 824C-------------------------------------------------------- 825C Calculate the term <mu1|[H,W^BD(3)]|HF> 826C added in OMEGA1EFF 827C-------------------------------------------------------- 828C 829 CALL CC3_CY1(OMEGA1EFF,ISYRES,WORK(KWMAT),ISWMAT, 830 * WORK(KTMAT),WORK(KXIAJB), 831 * ISINT1,WORK(KINDSQW),LENSQW, 832 * WORK(KEND4),LWRK4, 833 * ISYMB,B,ISYMD,D) 834C 835C------------------------------------------------------ 836C Calculate the term <mu2|[H,W^BD(3)]|HF> 837C ( Fock matrix cont ) 838C added in OMEGA2EFF 839C------------------------------------------------------ 840C 841 CALL CC3_CY2F(OMEGA2EFF,ISYRES,WORK(KWMAT),ISWMAT, 842 * WORK(KTMAT),WORK(KFOCK0),ISYM0, 843 * WORK(KINDSQW), 844 * LENSQW,WORK(KEND4),LWRK4, 845 * ISYMB,B,ISYMD,D) 846c 847C------------------------------------------------------ 848C Calculate the term <mu2|[H,W^BD(3)]|HF> 849C ( Occupied cont ) 850C added in OMEGA2EFF 851C------------------------------------------------------ 852C 853 CALL CC3_CY2O(OMEGA2EFF,ISYRES,WORK(KWMAT),ISWMAT, 854 * WORK(KTMAT),WORK(KTROC), 855 * WORK(KTROC1),ISINT1,WORK(KEND4),LWRK4, 856 * WORK(KINDSQW),LENSQW, 857 * ISYMB,B,ISYMD,D) 858C 859C 860C------------------------------------------------------ 861C Calculate the term <mu2|[H,W^BD(3)]|HF> 862C ( Virtual cont ) 863C added in OMEGA2EFF 864C------------------------------------------------------ 865C 866 CALL CC3_CY2V(OMEGA2EFF,ISYRES,WORK(KRBJIA), 867 * WORK(KWMAT), 868 * ISWMAT,WORK(KTMAT),WORK(KTRVI), 869 * WORK(KTRVI1),ISINT1,WORK(KEND4),LWRK4, 870 * WORK(KINDSQW),LENSQW, 871 * ISYMB,B,ISYMD,D) 872C 873 END DO !B 874 END DO !ISYMB 875C 876 END DO !D 877 END DO !ISYMD 878c 879c 880c write(lupri,*) 'T3XY in CC3_BMATSD : ' 881c call print_pt3(work(kx3am),1,4) 882c 883 884C 885C 886C------------------------------------------------------ 887C Accumulate RBJIA from <mu2|[H,W^BD(3)]|HF> ( Virtual cont ) 888C in OMEGA2EFF 889C------------------------------------------------------ 890C 891 CALL CC3_RBJIA(OMEGA2EFF,ISYRES,WORK(KRBJIA)) 892c 893c write(lupri,*)'OMEGA1EFF after CC3_BMATSD' 894c call output(OMEGA1EFF,1,nvir(1),1,nrhf(1),nvir(1),nrhf(1),1,lupri) 895c 896c write(lupri,*)'OMEGA2EFF after CC3_BMATSD' 897c call outpak(OMEGA2EFF,NT1AM(1),1,lupri) 898 899C 900 IF (IPRINT .GT. 55) THEN 901 XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 902 WRITE(LUPRI,*)'Norm of OMEGA2EFF final before added ',XNORMVAL 903 ENDIF 904C 905 DO I = 1,NT2AM(ISYRES) 906 OMEGA2EFF(I) = OMEGA2EFF(I) + OMEGA2(I) 907 END DO 908C 909 IF (IPRINT .GT. 55) THEN 910 XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 911 WRITE(LUPRI,*)'Norm of OMEGA2EFF final, OMEGA2EFF + OMEGA2F ', 912 * XNORMVAL 913 ENDIF 914C 915 DO I = 1,NT1AM(ISYRES) 916 OMEGA1EFF(I) = OMEGA1EFF(I) + OMEGA1(I) 917 END DO 918C 919 IF (IPRINT .GT. 55) THEN 920 XNORMVAL = DDOT(NT1AM(ISYRES),OMEGA1EFF,1,OMEGA1EFF,1) 921 WRITE(LUPRI,*) 'Norm of OMEGA1EFF final, OMEGA1EFF + OMEGA1 ', 922 * XNORMVAL 923 ENDIF 924C 925C------------------------------------------- 926C Close files (integrals in contraction) 927C------------------------------------------- 928C 929 CALL WCLOSE2(LUTOC,FNTOC,'KEEP') 930 CALL WCLOSE2(LU3VI,FN3VI,'KEEP') 931 CALL WCLOSE2(LU3VI2,FN3VI2,'KEEP') 932C 933C------------------------------------------ 934C Close temporary files (H^B integrals) 935C------------------------------------------ 936C 937 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 938 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 939 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 940 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 941C 942C------------------------------------------ 943C Close temporary files (H^BC integrals) 944C------------------------------------------ 945C 946 CALL WCLOSE2(LUCKJDR1,FNCKJDR1,'DELETE') 947 CALL WCLOSE2(LUDKBCR1,FNDKBCR1,'DELETE') 948C 949C------------------------------------------ 950C Close temporary files (H^C integrals) 951C------------------------------------------ 952C 953 CALL WCLOSE2(LUCKJDR2,FNCKJDR2,'DELETE') 954 CALL WCLOSE2(LUDKBCR2,FNDKBCR2,'DELETE') 955C 956C------------- 957C End 958C------------- 959C 960 CALL QEXIT('CC3_BMATSD ') 961C 962 RETURN 963 END 964C 965