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