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_bfmat */ 20 SUBROUTINE CC3_BFMAT(LISTL,IDLSTL,LISTR,IDLSTR, 21 * OMEGA1,OMEGA2,OMEGA1EFF,OMEGA2EFF, 22 * ISYRES, 23 * WORK,LWORK, 24 * LUTOC,FNTOC, 25 * LU3VI,FN3VI,LUDKBC3,FNDKBC3, 26 * LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X, 27 * LU3VI2,FN3VI2,LU3FOP,FN3FOP,LU3FOP2,FN3FOP2) 28* 29********************************************************************** 30* 31* Calculate following contributions to fmat: 32* 33* <L3|[H^0,T^{B}_{2}],tau_{1}]|HF> 34* <L3|[H^B,T^{0}_{2}],tau_{1}]|HF> 35* <L3|[H^B,\tau_nu_2]|HF> 36* 37* L0 or L1 list may be used for the multipliers; thus fmat or bmat 38* contributions may be calculated, respectively. 39* 40* If we have L0 then WB3X = .false. and SKIPGEI = .true. 41* 42* If we have L1 (or other similar lists---like LE,M1,N2---that require 43* W intermediate) then WB3X = .true. and SKIPGEI = .false. 44* 45* 46* F. Pawlowski and P. Jorgensen, Spring 2003. 47* 48* (based on the old cc3_fmat routine written by K. Hald) 49* 50* April-2004, Aarhus, FP: VVVV integrals removed, flags LVVVV and SKIPGEI 51* introduced. 52* 53********************************************************************** 54* 55 IMPLICIT NONE 56C 57#include "priunit.h" 58#include "ccorb.h" 59#include "ccl1rsp.h" 60#include "ccsdsym.h" 61#include "ccr1rsp.h" 62#include "dummy.h" 63#include "iratdef.h" 64#include "ccinftap.h" 65#include "cclrmrsp.h" 66#include "ccexci.h" 67#include "ccn2rsp.h" 68#include "ccsdinp.h" 69C 70 LOGICAL LSKIPL1R,SKIPGEI 71C 72 INTEGER ISYM0 73 PARAMETER(ISYM0 = 1) 74C 75 CHARACTER LISTL0*3, LISTL*3,LISTR*3,LISTL1R*3,LABELL1*8 76 CHARACTER*(*) FNTOC, FN3VI, FNDKBC3, FN3FOPX, FN3FOP2X 77 CHARACTER*(*) FN3VI2,FN3FOP,FN3FOP2 78C 79 CHARACTER*13 FNDELDR,FNDKBCR,FNDKBCR4,FNCKJDR,FN4V,FN4VB,FN3SRTR 80C 81 INTEGER IDLSTL0,IDLSTL,IDLSTR,IDLSTL1R 82 INTEGER LUTOC, LU3VI, LUDKBC3, LU3FOPX, LU3FOP2X 83 INTEGER LU3VI2,LU3FOP,LU3FOP2 84 INTEGER LUDELDR,LUDKBCR,LUDKBCR4,LUCKJDR,LU4V,LU4VB,LU3SRTR 85C 86 CHARACTER*6 FNGEI,FNFEI 87 CHARACTER*5 FNN1 88 PARAMETER( FNGEI = 'N1_GEI' , FNFEI = 'N1_FEI' , FNN1 = 'N1MAT' ) 89 INTEGER LUGEI,LUFEI,LUN1 90 91 CHARACTER*7 FNGEIB,FNFEIB 92 CHARACTER*6 FNN1B 93 PARAMETER(FNGEIB='N1_GEIB' , FNFEIB='N1_FEIB' , FNN1B='N1MATB' ) 94 INTEGER LUGEIB,LUFEIB,LUN1B 95C 96 CHARACTER CDUMMY*1 97C 98 LOGICAL LOCDBG,LORXL1 99 PARAMETER (LOCDBG = .FALSE.) 100C 101 INTEGER ISYRES,LWORK 102 INTEGER ISYML1,ISYML1R,ISYMR1 103 INTEGER ISINT2,KRBJIA,KLAMP0,KLAMH0,KFOCKD,KFOCK0CK,KT2TP,KL1AM 104 INTEGER KL2TP,KEND1,LWRK1 105 INTEGER KL1L1,KL2L1,KFOCKL1 106 INTEGER KT1R1,KT2R1,KFOCK0 107 INTEGER IOPT 108 INTEGER ISINT1R1,ISINT2R1 109 INTEGER KOIOOOO,KOIOVVO,KOIOOVV,KBIOOOO,KBIOVVO,KBIOOVV 110 INTEGER ISINT2L1R,ISYFCKL1R,KXIAJB,KT3BOG1,KT3BOL1,KT3BOG2 111 INTEGER KT3BOL2,KLAMPL1R,KLAMHL1R 112 INTEGER KFOCKL1RCK,KW3BXOGX1,KW3BXOLX1 113 INTEGER KW3BXOG1,KW3BXOL1,KT1L1R 114 INTEGER LENGTH 115 INTEGER ISINT1,ISINT1L1R 116 INTEGER IOFF 117 INTEGER ISYMD,ISYCKBD0,ISYCKBDL1R,ISYCKBDR1 118 INTEGER KVVVVO,KVVVVB 119 INTEGER KT3BVDL1,KT3BVDL2,KT3BVDL3,KEND3,LWRK3 120 INTEGER KT3BVDG1,KT3BVDG2,KT3BVDG3 121 INTEGER KW3BXVDG1,KW3BXVDG2,KW3BXVDL1,KW3BXVDL2,KW3BXVDGX1 122 INTEGER KW3BXVDGX2,KW3BXVDLX1,KW3BXVDLX2 123 INTEGER KTRVIR,KTRVIR1,KEND4,LWRK4 124 INTEGER ISYMB,ISYALJB0,ISYALJD0,ISYALJBL1,ISYALJDL1,ISYMBD 125 INTEGER ISCKIJ,ISWBMAT,ISYCKD 126 INTEGER KSMAT2,KUMAT2,KDIAG,KINDSQ 127 INTEGER KDIAGWB,KINDSQWB 128 INTEGER KINDEX,KINDEX2 129 INTEGER KINDEXBL1,KINDEXDL1,KTMAT,KW3BMAT 130 INTEGER KT3BVBG1,KT3BVBG2,KT3BVBG3,KSMAT4,KUMAT4,KT3BVBL1 131 INTEGER KT3BVBL2,KT3BVBL3,KEND5,LWRK5 132 INTEGER KINTOC,KTROCR,KTROCR1 133 INTEGER ISYML0 134 135 INTEGER LENSQ,LENSQWB 136 INTEGER ISTBD,KTB,LENSQTB,KINDSQTB,ISTB 137C 138 INTEGER IR1TAMP 139 INTEGER ILSTSYM 140C 141 INTEGER KEND0,LWRK0 142C 143 INTEGER ISYMN1,ISYMN2,ISYMN1B,ISYMN2B,KN2MAT,KN2MATB,KINDSQN 144 INTEGER KINDSQNB,LENSQN,LENSQNB,ISGEI,ISFEI,ISGEIB,ISFEIB,KGEI 145 INTEGER KFEI,KGEIB,KFEIB,IADR,KLAMPB,KLAMHB 146C 147 INTEGER LENGTHGEI,LENGTHGEIB 148c 149 integer kx3am 150c 151 LOGICAL WB3X 152C 153#if defined (SYS_CRAY) 154 REAL OMEGA1(*),OMEGA2(*),OMEGA1EFF(*),OMEGA2EFF(*) 155 REAL WORK(LWORK) 156 REAL FREQL1,FREQL1R 157 REAL DDOT,XNORMVAL 158#else 159 DOUBLE PRECISION OMEGA1(*),OMEGA2(*),OMEGA1EFF(*),OMEGA2EFF(*) 160 DOUBLE PRECISION WORK(LWORK) 161 DOUBLE PRECISION FREQL1,FREQL1R 162 DOUBLE PRECISION DDOT,XNORMVAL 163#endif 164C 165 CALL QENTER('BFMAT') 166C 167c write(lupri,*)'BEFORE ' 168c write(lupri,*)'omega1eff before CC3_BFMAT, LISTL ', LISTL 169c call PRINT_MATAI(OMEGA1EFF,ISYRES) 170c xnormval = ddot(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 171c write(lupri,*)'norm omega2eff before CC3_BFMAT ', xnormval 172c CALL OUTPAK(OMEGA2EFF,NT1AMX,1,LUPRI) 173C 174C--------------------------------------------------------------- 175C Initialise character strings, logical flags and open files 176C--------------------------------------------------------------- 177C 178 CDUMMY = ' ' 179C 180 SKIPGEI = .FALSE. 181C 182 LU4V = -1 183 LU4VB = -1 184 LU3SRTR = -1 185 LUCKJDR = -1 186 LUDELDR = -1 187 LUDKBCR = -1 188 LUDKBCR4 = -1 189C 190 FN4V = 'CC3_FMAT_TMP2' 191 FN4VB = 'CC3_FMAT_TMP3' 192 FN3SRTR = 'CC3_FMAT_TMP4' 193 FNCKJDR = 'CC3_FMAT_TMP5' 194 FNDELDR = 'CC3_FMAT_TMP6' 195 FNDKBCR = 'CC3_FMAT_TMP7' 196 FNDKBCR4 = 'CC3_FMAT_TMP8' 197C 198 IF (LVVVV) THEN 199 CALL WOPEN2(LU4V,FN4V,64,0) 200 CALL WOPEN2(LU4VB,FN4VB,64,0) 201 END IF 202C 203 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 204 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 205 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 206 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 207 CALL WOPEN2(LUDKBCR4,FNDKBCR4,64,0) 208C 209C------------------------------------------------------------ 210C lists handling 211C------------------------------------------------------------ 212C 213c LISTL0 = 'L0 ' 214c IDLSTL0 = 0 215c ISYML0 = ISYM0 216C 217 IF (LISTL(1:3).EQ.'L1 ') THEN 218 219 ! get symmetry, frequency and integral label from common blocks 220 ! defined in ccl1rsp.h 221 ISYML1 = ISYLRZ(IDLSTL) 222 FREQL1 = FRQLRZ(IDLSTL) 223 LABELL1 = LRZLBL(IDLSTL) 224 LORXL1 = LORXLRZ(IDLSTL) 225 226 IF (LORXL1) CALL QUIT('NO ORBITAL RELAX. IN CC3_BFMAT') 227 228 LISTL1R = 'R1 ' 229 IDLSTL1R = IR1TAMP(LABELL1,LORXL1,FREQL1,ISYML1) 230 ! get symmetry and frequency from common blocks 231 ! defined in ccl1rsp.h 232 ISYML1R = ISYLRT(IDLSTL1R) 233 FREQL1R = FRQLRT(IDLSTL1R) 234C 235 LISTL0 = 'L0 ' 236 IDLSTL0 = 0 237 ISYML0 = ISYM0 238C 239 IF (ISYML1 .NE. ISYML1R) THEN 240 WRITE(LUPRI,*)'ISYML1: ', ISYML1 241 WRITE(LUPRI,*)'ISYML1R: ', ISYML1R 242 CALL QUIT('Symmetry mismatch in CC3_BFMAT') 243 END IF 244C 245 IF (FREQL1R .NE. FREQL1) THEN 246 WRITE(LUPRI,*)'FREQL1R: ', FREQL1R 247 WRITE(LUPRI,*)'FREQL1: ', FREQL1 248 CALL QUIT('Frequency mismatch in CC3_BFMAT(L1)') 249 END IF 250C 251 ELSE IF (LISTL(1:3).EQ.'M1 ') THEN 252 ISYML1 = ILSTSYM(LISTL,IDLSTL) 253 FREQL1 = FRQLRM(IDLSTL) 254 LABELL1 = '- none -' 255C 256 ! find corresponding right eigenvector 257 LISTL1R = 'RE ' 258 IDLSTL1R = ILRM(IDLSTL) 259 ISYML1R = ISYML1 260 FREQL1R = EIGVAL(IDLSTL1R) 261C 262 LISTL0 = 'L0 ' 263 IDLSTL0 = 0 264 ISYML0 = ISYM0 265C 266 IF (FREQL1R .NE. FREQL1) THEN 267 WRITE(LUPRI,*)'FREQL1R: ', FREQL1R 268 WRITE(LUPRI,*)'FREQL1: ', FREQL1 269 CALL QUIT('Frequency mismatch in CC3_BFMAT(M1)') 270 END IF 271C 272 ELSE IF (LISTL(1:3).EQ.'N2 ') THEN 273 ISYML1 = ILSTSYM(LISTL,IDLSTL) 274 FREQL1 = FRQIN2(IDLSTL) + FRQFN2(IDLSTL) 275 LABELL1 = '- none -' 276C 277 ! find corresponding right eigenvector 278 LISTL1R = 'RE ' 279 IDLSTL1R = IFN2(IDLSTL) 280 ISYML1R = ILSTSYM(LISTL1R,IDLSTL1R) 281 FREQL1R = FRQFN2(IDLSTL) 282C 283 !LITSL0 corresponding to LISTL 284 LISTL0 = 'LE ' 285 IDLSTL0 = IIN2(IDLSTL) 286 ISYML0 = ILSTSYM(LISTL0,IDLSTL0) 287C 288 ELSE IF (LISTL(1:3).EQ.'LE ') THEN 289 ISYML1 = ILSTSYM(LISTL,IDLSTL) 290 FREQL1 = -EIGVAL(IDLSTL) 291 LABELL1 = '- none -' 292C 293 !we don't have any "right" vector entering a right hand side 294 LISTL1R = '---' 295 IDLSTL1R = -99 296 297C !LISTL0 not used for LE (...but does not hurt ot have them defined) 298 LISTL0 = 'L0 ' 299 IDLSTL0 = 0 300 ISYML0 = ISYM0 301C 302 ELSE IF (LISTL(1:3).EQ.'L0 ') THEN 303 LISTL0 = 'L0 ' 304 IDLSTL0 = 0 305 ISYML0 = ISYM0 306C 307 SKIPGEI = .TRUE. 308 ELSE 309 CALL QUIT('Unknown left list in CC3_BFMAT') 310 END IF 311 312C 313 IF (.NOT.LVVVV) THEN 314 !Open files for N1MAT intermediates 315 LUGEI = -1 316 LUFEI = -1 317 LUN1 = -1 318 CALL WOPEN2(LUFEI,FNFEI,64,0) 319 CALL WOPEN2(LUN1,FNN1,64,0) 320 IF (.NOT.SKIPGEI) THEN 321 CALL WOPEN2(LUGEI,FNGEI,64,0) 322 END IF 323 !Open files for N1MAT^B intermediates 324 LUGEIB = -1 325 LUFEIB = -1 326 LUN1B = -1 327 CALL WOPEN2(LUFEIB,FNFEIB,64,0) 328 CALL WOPEN2(LUN1B,FNN1B,64,0) 329 IF (.NOT.SKIPGEI) THEN 330 CALL WOPEN2(LUGEIB,FNGEIB,64,0) 331 END IF 332 END IF 333C 334c IF (LISTR(1:3).EQ.'R1 ') THEN 335c ! get symmetry, frequency and integral label for right list 336c ! from common blocks defined in ccr1rsp.h 337c ISYMR1 = ISYLRT(IDLSTR) 338c ELSE 339c CALL QUIT('Unknown right list in CC3_BFMAT') 340c END IF 341 ISYMR1 = ILSTSYM(LISTR,IDLSTR) 342C 343 ISINT1 = ISYM0 344 ISINT2 = ISYM0 345C 346 IF (.NOT.LVVVV) THEN 347 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 348 * .OR. (LISTL(1:3).EQ.'M1 ') 349 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 350 !Symmetries for N1 and N2 intermediates 351 ISYMN1 = MULD2H(ISYML1,ISYM0) 352 ISYMN2 = MULD2H(ISYML1,ISYM0) 353 !Symmetries for N1^B and N2^B intermediates 354 ISYMN1B = MULD2H(ISYML1,ISYMR1) 355 ISYMN2B = MULD2H(ISYML1,ISYMR1) 356 ELSE IF (LISTL(1:3).EQ.'L0 ') THEN 357 !Symmetries for N1 and N2 intermediates 358 ISYMN1 = MULD2H(ISYML0,ISYM0) 359 ISYMN2 = MULD2H(ISYML0,ISYM0) 360 !Symmetries for N1^B and N2^B intermediates 361 ISYMN1B = MULD2H(ISYML0,ISYMR1) 362 ISYMN2B = MULD2H(ISYML0,ISYMR1) 363 ELSE 364 CALL QUIT('Unknown left list in CC3_BFMAT(x)') 365 END IF 366 END IF 367C 368 IF (LVVVV) THEN 369 KRBJIA = 1 370 ELSE 371 KN2MAT = 1 372 KN2MATB = KN2MAT + NCKIJ(ISYMN2) 373 KRBJIA = KN2MATB + NCKIJ(ISYMN2B) 374 END IF 375 KLAMP0 = KRBJIA + NT2SQ(ISYRES) 376 KLAMH0 = KLAMP0 + NLAMDT 377 KFOCKD = KLAMH0 + NLAMDT 378 KFOCK0CK = KFOCKD + NORBTS 379 KT2TP = KFOCK0CK + NT1AMX 380 KL1AM = KT2TP + NT2SQ(ISYM0) 381 KL2TP = KL1AM + NT1AM(ISYML0) 382 KEND0 = KL2TP + NT2SQ(ISYML0) 383 LWRK0 = LWORK - KEND0 384C 385 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 386 * .OR. (LISTL(1:3).EQ.'M1 ') 387 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 388 KL1L1 = KEND0 389 KL2L1 = KL1L1 + NT1AM(ISYML1) 390 KFOCKL1 = KL2L1 + NT2SQ(ISYML1) 391 KEND0 = KFOCKL1 + N2BST(ISYML1) 392 LWRK0 = LWORK - KEND0 393 END IF 394C 395 KT1R1 = KEND0 396 KT2R1 = KT1R1 + NT1AM(ISYMR1) 397 KFOCK0 = KT2R1 + NT2SQ(ISYMR1) 398 KEND0 = KFOCK0 + N2BST(ISYM0) 399 LWRK0 = LWORK - KEND0 400C 401 IF (.NOT.LVVVV) THEN 402 KINDSQN = KEND0 403 KINDSQNB = KINDSQN + (6*NCKIJ(ISYMN2) - 1)/IRAT + 1 404 KEND0 = KINDSQNB + (6*NCKIJ(ISYMN2B) - 1)/IRAT + 1 405 LWRK0 = LWORK - KEND0 406 END IF 407C 408 IF (LWRK0 .LT. 0) THEN 409 WRITE(LUPRI,*) 'Memory available : ',LWORK 410 WRITE(LUPRI,*) 'Memory needed : ',KEND0 411 CALL QUIT('Insufficient space in CC3_BFMAT') 412 END IF 413C 414 CALL DZERO(WORK(KRBJIA),NT2SQ(ISYRES)) 415C 416 IF (.NOT.LVVVV) THEN 417 CALL DZERO(WORK(KN2MAT),NCKIJ(ISYMN2)) 418 CALL DZERO(WORK(KN2MATB),NCKIJ(ISYMN2B)) 419C 420 !index array for N2 421 LENSQN = NCKIJ(ISYMN2) 422 CALL CC3_INDSQ(WORK(KINDSQN),LENSQN,ISYMN2) 423 !index array for N2^B 424 LENSQNB = NCKIJ(ISYMN2B) 425 CALL CC3_INDSQ(WORK(KINDSQNB),LENSQNB,ISYMN2B) 426 END IF 427C 428C------------------------------------- 429C Read in lamdap and lamdh 430C------------------------------------- 431C 432 CALL GET_LAMBDA0(WORK(KLAMP0),WORK(KLAMH0),WORK(KEND0),LWRK0) 433C 434C--------------------------------------------------------------------- 435C Read zeroth-order AO Fock matrix from file and trasform it to 436C lambda basis 437C--------------------------------------------------------------------- 438C 439 CALL GET_FOCK0(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0),WORK(KEND0), 440 * LWRK0) 441C 442C--------------------------------------------------------------------- 443C Read the matrix the property integrals and trasform it to lambda 444C basis for L1 list 445C--------------------------------------------------------------------- 446C 447 IF (LISTL(1:3).EQ.'L1 ') THEN 448 CALL GET_FOCKX(WORK(KFOCKL1),LABELL1,IDLSTL,ISYML1, 449 * WORK(KLAMP0),ISYM0, 450 * WORK(KLAMH0),ISYM0,WORK(KEND0),LWRK0) 451 END IF 452C 453C------------------------------------- 454C Read T2 zeroth-order amplitudes 455C------------------------------------- 456C 457 IOPT = 2 458 CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KT2TP),'R0',0,ISYM0, 459 * WORK(KEND0),LWRK0) 460C 461 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of T2TP ', 462 * DDOT(NT2SQ(ISYM0),WORK(KT2TP),1,WORK(KT2TP),1) 463C 464C-------------------------------------------- 465C Read L1 and L2 zeroth-order multipliers 466C-------------------------------------------- 467C 468 IOPT = 3 469 CALL GET_T1_T2(IOPT,.FALSE.,WORK(KL1AM),WORK(KL2TP),LISTL0, 470 * IDLSTL0,ISYML0,WORK(KEND0),LWRK0) 471C 472 IF (LOCDBG) WRITE(LUPRI,*) 'Norm of L2TP ', 473 * DDOT(NT2SQ(ISYML0),WORK(KL2TP),1,WORK(KL2TP),1) 474C 475C--------------------------------------------- 476C Read L1L1 and L2L1 multipliers (L1 list) 477C--------------------------------------------- 478C 479 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 480 * .OR. (LISTL(1:3).EQ.'M1 ') 481 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 482 IOPT = 3 483 CALL GET_T1_T2(IOPT,.FALSE.,WORK(KL1L1),WORK(KL2L1),LISTL, 484 * IDLSTL,ISYML1,WORK(KEND0),LWRK0) 485 END IF 486C 487C------------------------------------- 488C Read T1R1 and L2R1 amplitudes 489C------------------------------------- 490C 491 IOPT = 3 492 CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1R1),WORK(KT2R1),LISTR, 493 * IDLSTR,ISYMR1,WORK(KEND0),LWRK0) 494C 495C---------------------------------------- 496C Integrals [H,T1Y] where Y is LISTR 497C---------------------------------------- 498C 499 ISINT1R1 = MULD2H(ISINT1,ISYMR1) 500 ISINT2R1 = MULD2H(ISINT2,ISYMR1) 501C 502 CALL CC3_BARINT(WORK(KT1R1),ISYMR1,WORK(KLAMP0), 503 * WORK(KLAMH0),WORK(KEND0),LWRK0, 504 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 505C 506 CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISINT1R1,LU3SRTR,FN3SRTR, 507 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 508 * IDUMMY,CDUMMY) 509C 510 CALL CC3_SINT(WORK(KLAMH0),WORK(KEND0),LWRK0,ISINT1R1, 511 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 512C 513 !with option IOPT = 2 it just calculates LUDKBCR4 needed 514 !in contractions 515 IOPT = 2 516 CALL CC3_TCME(DUMMY,ISINT1R1,WORK(KEND0),LWRK0,IDUMMY,CDUMMY, 517 * LUDKBCR,FNDKBCR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 518 * IDUMMY,CDUMMY,LUDKBCR4,FNDKBCR4,IOPT) 519C 520C--------------------------------------------------------------- 521C Read canonical orbital energies and delete frozen orbitals 522C in Fock diagonal, if required 523C--------------------------------------------------------------- 524C 525 CALL GET_ORBEN(WORK(KFOCKD),WORK(KEND0),LWRK0) 526C 527C-------------------------------------------- 528C Sort the Fock matrix to get F(ck) block 529C-------------------------------------------- 530C 531 CALL SORT_FOCKCK(WORK(KFOCK0CK),WORK(KFOCK0),ISYM0) 532C 533C---------------------------------------- 534C If we want to sum the T3 amplitudes 535C---------------------------------------- 536C 537 if (.false.) then 538 kx3am = kend0 539 kend0 = kx3am + nrhft*nrhft*nrhft*nvirt*nvirt*nvirt 540 call dzero(work(kx3am),nrhft*nrhft*nrhft*nvirt*nvirt*nvirt) 541 lwrk0 = lwork - kend0 542 if (lwrk0 .lt. 0) then 543 write(lupri,*) 'Memory available : ',lwork 544 write(lupri,*) 'Memory needed : ',kend0 545 call quit('Insufficient space (T3) in CC3_BFMAT') 546 END IF 547 endif 548C 549C write(lupri,*) 'WBMAT after dzero' 550C call print_pt3(work(kx3am),ISYML1,4) 551C 552C-------------------------------------------------------------- 553C Calculate the normal g^0 integrals for 554C OOOO, OVVO and OOVV integrals used in contraction 555C (VVVV is stored on file LU4V and read in in D-loop) 556C-------------------------------------------------------------- 557C 558 KOIOOOO = KEND0 559 KOIOVVO = KOIOOOO + N3ORHF(ISYM0) 560 KOIOOVV = KOIOVVO + NT2SQ(ISYM0) 561 KEND0 = KOIOOVV + NT2SQ(ISYM0) 562 LWRK0 = LWORK - KEND0 563C 564 IF (LWRK0 .LT. 0) THEN 565 CALL QUIT('Out of memory in CC3_BFMAT (g^0[2o2v] kind)') 566 ENDIF 567C 568 CALL DZERO(WORK(KOIOVVO),NT2SQ(ISYMOP)) 569 CALL DZERO(WORK(KOIOOVV),NT2SQ(ISYMOP)) 570C 571 CALL CC3_INT(WORK(KOIOOOO),WORK(KOIOOVV),WORK(KOIOVVO), 572 * ISYMOP,LU4V,FN4V, 573 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 574 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 575 * WORK(KEND0),LWRK0) 576C 577C--------------------------------------------------------------- 578C Calculate the g^B integrals for 579C OOOO^B, OVVO^B and OOVV^B integrals used in contraction 580C (VVVV^B is stored on file LU4VB and read in in D-loop) 581C--------------------------------------------------------------- 582C 583 KBIOOOO = KEND0 584 KBIOVVO = KBIOOOO + N3ORHF(ISINT1R1) 585 KBIOOVV = KBIOVVO + NT2SQ(ISINT1R1) 586 KEND0 = KBIOOVV + NT2SQ(ISINT1R1) 587 LWRK0 = LWORK - KEND0 588C 589 IF (LWRK0 .LT. 0) THEN 590 CALL QUIT('Out of memory in CC3_BFMAT (g^B[2o2v] kind)') 591 ENDIF 592C 593 CALL DZERO(WORK(KBIOOOO),N3ORHF(ISINT1R1)) 594 CALL DZERO(WORK(KBIOVVO),NT2SQ(ISINT1R1)) 595 CALL DZERO(WORK(KBIOOVV),NT2SQ(ISINT1R1)) 596C 597 CALL CC3_INT(WORK(KBIOOOO),WORK(KBIOOVV),WORK(KBIOVVO), 598 * ISINT1R1,LU4VB,FN4VB, 599 * .TRUE.,LISTR,IDLSTR,ISYMR1, 600 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 601 * WORK(KEND0),LWRK0) 602C 603C----------------------------- 604C Memory allocation. 605C----------------------------- 606C 607 IF ( (LISTL(1:3).EQ.'L1 ') 608 * .OR. (LISTL(1:3).EQ.'M1 ') 609 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 610 ISINT2L1R = MULD2H(ISYML1R,ISINT2) 611 ISYFCKL1R = MULD2H(ISYMOP,ISYML1R) 612 END IF 613C 614 KXIAJB = KEND0 615 KEND0 = KXIAJB + NT2AM(ISYM0) 616 617 KT3BOG1 = KEND0 618 KT3BOL1 = KT3BOG1 + NTRAOC(ISYM0) 619 KT3BOG2 = KT3BOL1 + NTRAOC(ISYM0) 620 KT3BOL2 = KT3BOG2 + NTRAOC(ISYM0) 621 KEND0 = KT3BOL2 + NTRAOC(ISYM0) 622 LWRK0 = LWORK - KEND0 623C 624 IF ( (LISTL(1:3).EQ.'L1 ') 625 * .OR. (LISTL(1:3).EQ.'M1 ') 626 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 627 KFOCKL1RCK = KEND0 628 KW3BXOGX1 = KFOCKL1RCK + NT1AM(ISYFCKL1R) 629 KW3BXOLX1 = KW3BXOGX1 + NTRAOC(ISINT2L1R) 630 KEND0 = KW3BXOLX1 + NTRAOC(ISINT2L1R) 631 LWRK0 = LWORK - KEND0 632 END IF 633C 634 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 635 * .OR. (LISTL(1:3).EQ.'M1 ') 636 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 637 KW3BXOG1 = KEND0 638 KW3BXOL1 = KW3BXOG1 + NTRAOC(ISYM0) 639 KEND0 = KW3BXOL1 + NTRAOC(ISYM0) 640 LWRK0 = LWORK - KEND0 641 END IF 642C 643 IF ( (LISTL(1:3).EQ.'L1 ') 644 * .OR. (LISTL(1:3).EQ.'M1 ') 645 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 646 KT1L1R = KEND0 647 KEND0 = KT1L1R + NT1AM(ISYML1R) 648 LWRK0 = LWORK - KEND0 649C 650 KLAMPL1R = KEND0 651 KLAMHL1R = KLAMPL1R + NLAMDT 652 KEND0 = KLAMHL1R + NLAMDT 653 LWRK0 = LWORK - KEND0 654 END IF 655C 656 KINTOC = KEND0 657 KTROCR = KINTOC + NTOTOC(ISINT1R1) 658 KTROCR1 = KTROCR + NTRAOC(ISINT1R1) 659 KEND0 = KTROCR1 + NTRAOC(ISINT1R1) 660 LWRK0 = LWORK - KEND0 661C 662 IF (LWRK0 .LT. 0) THEN 663 WRITE(LUPRI,*) 'Memory available : ',LWORK 664 WRITE(LUPRI,*) 'Memory needed : ',KEND0 665 CALL QUIT('Insufficient space in CC3_BFMAT') 666 END IF 667 668C 669C------------------------ 670C Construct L(ia,jb). 671C------------------------ 672C 673 LENGTH = IRAT*NT2AM(ISYM0) 674 675 REWIND(LUIAJB) 676 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 677 678 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYM0,1) 679C 680C--------------------------------------------------- 681C Prepare to construct the occupied integrals... 682C--------------------------------------------------- 683C 684C isint1 - symmetry of integrals in standard H, transformed 685C with LambdaH_0 686C ISINT1L1R - symmetry of integrals in standard H, transformed 687C with LambdaH_L1R 688 689 ISINT1 = 1 690C 691C-------------------------- 692C Get Lambda for right list depended on left LISTL list 693C-------------------------- 694C 695 IF ( (LISTL(1:3).EQ.'L1 ') 696 * .OR. (LISTL(1:3).EQ.'M1 ') 697 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 698C 699 ISINT1L1R = MULD2H(ISINT1,ISYML1R) 700C 701 CALL GET_LAMBDAX(WORK(KLAMPL1R),WORK(KLAMHL1R),LISTL1R, 702 * IDLSTL1R, 703 * ISYML1R, 704 * WORK(KLAMP0),WORK(KLAMH0),WORK(KEND0),LWRK0) 705 706C 707C------------------------------------------------------------------ 708C Calculate the F^L1R matrix (kc elements evaluated and stored 709C as ck) 710C------------------------------------------------------------------ 711C 712 IOPT = 1 713 CALL GET_T1_T2(IOPT,.FALSE.,WORK(KT1L1R),DUMMY,LISTL1R, 714 * IDLSTL1R,ISYML1R,WORK(KEND0),LWRK0) 715 CALL CC3LR_MFOCK(WORK(KFOCKL1RCK),WORK(KT1L1R),WORK(KXIAJB), 716 * ISYFCKL1R) 717C 718 END IF 719C 720C----------------------------------------------------------------- 721C Construct occupied integrals which are required to calculate 722C t3bar_0 multipliers 723C----------------------------------------------------------------- 724C 725 CALL INTOCC_T3BAR0(LUTOC,FNTOC,WORK(KLAMH0),ISYM0,WORK(KT3BOG1), 726 * WORK(KT3BOL1),WORK(KT3BOG2),WORK(KT3BOL2), 727 * WORK(KEND0),LWRK0) 728C 729C----------------------------------------------------------------- 730C Construct occupied integrals which are required to calculate 731C t3bar_Y multipliers 732C----------------------------------------------------------------- 733C 734 IF ( (LISTL(1:3).EQ.'L1 ') 735 * .OR. (LISTL(1:3).EQ.'M1 ') 736 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 737 LSKIPL1R = .FALSE. 738 CALL INTOCC_T3BARX(LSKIPL1R, 739 * LUTOC,FNTOC,ISYMOP,WORK(KLAMH0),ISYM0, 740 * ISINT1, 741 * WORK(KLAMHL1R),ISYML1R,ISINT1L1R, 742 * WORK(KW3BXOG1), 743 * WORK(KW3BXOL1),WORK(KW3BXOGX1), 744 * WORK(KW3BXOLX1), 745 * WORK(KEND0),LWRK0) 746 ELSE IF (LISTL(1:3).EQ.'LE ') THEN 747 LSKIPL1R = .TRUE. 748 CALL INTOCC_T3BARX(LSKIPL1R, 749 * LUTOC,FNTOC,ISYMOP,WORK(KLAMH0),ISYM0, 750 * ISINT1, 751 * DUMMY,IDUMMY,IDUMMY, 752 * WORK(KW3BXOG1), 753 * WORK(KW3BXOL1),DUMMY, 754 * DUMMY, 755 * WORK(KEND0),LWRK0) 756 END IF 757C 758C-------------------------------------------------------------------- 759C Read in R1-transformed occupied integrals used in contractions, 760C transform with lambda_P^0 and sort 761C-------------------------------------------------------------------- 762C 763 IOFF = 1 764 IF (NTOTOC(ISINT1R1) .GT. 0) THEN 765 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF,NTOTOC(ISINT1R1)) 766 ENDIF 767C 768 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROCR),WORK(KLAMP0), 769 * WORK(KEND0),LWRK0,ISINT1R1) 770C 771 CALL CCFOP_SORT(WORK(KTROCR),WORK(KTROCR1),ISINT1R1,1) 772C 773 CALL CC3_LSORT1(WORK(KTROCR),ISINT1R1,WORK(KEND0),LWRK0,5) 774C 775 DO ISYMD = 1,NSYM 776C 777 ISYCKBD0 = MULD2H(ISYMD,ISYM0) 778 ISYCKBDR1 = MULD2H(ISYMD,ISINT2R1) 779C 780 IF (.NOT.LVVVV) THEN 781 !Symmetry of arrays needed to construct N1MAT 782 ISGEI = MULD2H(ISYMN1,ISYMD) 783 ISFEI = MULD2H(ISYMN1,ISYMD) 784 !Symmetry of arrays needed to construct N1MAT^B 785 ISGEIB = MULD2H(ISYMN1B,ISYMD) 786 ISFEIB = MULD2H(ISYMN1B,ISYMD) 787 END IF 788C 789 IF ( (LISTL(1:3).EQ.'L1 ') 790 * .OR. (LISTL(1:3).EQ.'M1 ') 791 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 792 ISYCKBDL1R = MULD2H(ISYMD,ISYML1R) 793 END IF 794C 795 IF (LVVVV) THEN 796 KVVVVO = KEND0 797 KVVVVB = KVVVVO + NMAABC(ISYCKBD0) 798 KEND1 = KVVVVB + NMAABC(ISYCKBDR1) 799 LWRK1 = LWORK - KEND1 800 ELSE 801 IF (SKIPGEI) THEN !we want to have the dummy allocations 802 !to skip GEI, but keep the code simple. 803 LENGTHGEI = 1 804 LENGTHGEIB = 1 805 ELSE 806 LENGTHGEI = NCKATR(ISGEI) 807 LENGTHGEIB = NCKATR(ISGEIB) 808 END IF 809 !Arrays needed to construct N1MAT 810 KGEI = KEND0 811 KFEI = KGEI + LENGTHGEI 812 KEND1 = KFEI + NCKATR(ISFEI) 813 LWRK1 = LWORK - KEND1 814C 815 !Arrays needed to construct N1MAT^B 816 KGEIB = KEND1 817 KFEIB = KGEIB + LENGTHGEIB 818 KEND1 = KFEIB + NCKATR(ISFEIB) 819 LWRK1 = LWORK - KEND1 820C 821 END IF 822C 823 IF (LWRK1 .LT. 0) THEN 824 WRITE(LUPRI,*) 'Memory available : ',LWORK 825 WRITE(LUPRI,*) 'Memory needed : ',KEND1 826 CALL QUIT('Insufficient space (isymd-loop) in CC3_BFMAT') 827 END IF 828C 829 DO D = 1,NVIR(ISYMD) 830C 831 KT3BVDL1 = KEND1 832 KT3BVDL2 = KT3BVDL1 + NCKATR(ISYCKBD0) 833 KT3BVDL3 = KT3BVDL2 + NCKATR(ISYCKBD0) 834 KEND3 = KT3BVDL3 + NCKATR(ISYCKBD0) 835 LWRK3 = LWORK - KEND3 836 837 KT3BVDG1 = KEND3 838 KT3BVDG2 = KT3BVDG1 + NCKATR(ISYCKBD0) 839 KT3BVDG3 = KT3BVDG2 + NCKATR(ISYCKBD0) 840 KEND3 = KT3BVDG3 + NCKATR(ISYCKBD0) 841 LWRK3 = LWORK - KEND3 842 843 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 844 * .OR. (LISTL(1:3).EQ.'M1 ') 845 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 846 KW3BXVDG1 = KEND3 847 KW3BXVDG2 = KW3BXVDG1 + NCKATR(ISYCKBD0) 848 KW3BXVDL1 = KW3BXVDG2 + NCKATR(ISYCKBD0) 849 KW3BXVDL2 = KW3BXVDL1 + NCKATR(ISYCKBD0) 850 KEND3 = KW3BXVDL2 + NCKATR(ISYCKBD0) 851 LWRK3 = LWORK - KEND3 852 END IF 853C 854 IF ( (LISTL(1:3).EQ.'L1 ') 855 * .OR. (LISTL(1:3).EQ.'M1 ') 856 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 857 KW3BXVDGX1 = KEND3 858 KW3BXVDGX2 = KW3BXVDGX1 + NCKATR(ISYCKBDL1R) 859 KW3BXVDLX1 = KW3BXVDGX2 + NCKATR(ISYCKBDL1R) 860 KW3BXVDLX2 = KW3BXVDLX1 + NCKATR(ISYCKBDL1R) 861 KEND3 = KW3BXVDLX2 + NCKATR(ISYCKBDL1R) 862 LWRK3 = LWORK - KEND3 863 END IF 864C 865 KTRVIR = KEND3 866 KTRVIR1 = KTRVIR + NCKATR(ISYCKBDR1) 867 KEND4 = KTRVIR1 + NCKATR(ISYCKBDR1) 868 LWRK4 = LWORK - KEND4 869 870 IF (LWRK4 .LT. 0) THEN 871 WRITE(LUPRI,*) 'Memory available : ',LWORK 872 WRITE(LUPRI,*) 'Memory needed : ',KEND4 873 CALL QUIT('Insufficient space (d-loop) in CC3_BFMAT') 874 END IF 875C 876 IF (.NOT.LVVVV) THEN 877 CALL DZERO(WORK(KFEI),NCKATR(ISFEI)) 878 CALL DZERO(WORK(KFEIB),NCKATR(ISFEIB)) 879 CALL DZERO(WORK(KGEI),LENGTHGEI) 880 CALL DZERO(WORK(KGEIB),LENGTHGEIB) 881 END IF 882 883 884 IF (LVVVV) THEN 885C 886C----------------------------------------------------------------- 887C Read in g^0_{vvvv} (used in contraction) for a given D 888C----------------------------------------------------------------- 889C 890 IF (NMAABC(ISYCKBD0) .GT. 0) THEN 891 IOFF = I3VVIR(ISYCKBD0,ISYMD) 892 * + NMAABC(ISYCKBD0)*(D-1) 893 * + 1 894 CALL GETWA2(LU4V,FN4V,WORK(KVVVVO),IOFF,NMAABC(ISYCKBD0)) 895 ENDIF 896C 897C----------------------------------------------------------------- 898C Read in g^B_{vvvv} (used in contraction) for a given D 899C----------------------------------------------------------------- 900C 901 IF (NMAABC(ISYCKBDR1) .GT. 0) THEN 902 IOFF = I3VVIR(ISYCKBDR1,ISYMD) 903 * + NMAABC(ISYCKBDR1)*(D-1) 904 * + 1 905 CALL GETWA2(LU4VB,FN4VB,WORK(KVVVVB),IOFF, 906 * NMAABC(ISYCKBDR1)) 907 ENDIF 908C 909 END IF !LVVVV 910C 911C----------------------------------------------------------------------- 912C Construct virtual integrals (for fixed D) which are required 913C to calculate t3bar_0 multipliers 914C----------------------------------------------------------------------- 915C 916 CALL INTVIR_T3BAR0_D(LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X, 917 * LUDKBC3,FNDKBC3,LU3VI,FN3VI,ISYM0, 918 * WORK(KT3BVDL1),WORK(KT3BVDG1), 919 * WORK(KT3BVDG2),WORK(KT3BVDL2), 920 * WORK(KT3BVDG3),WORK(KT3BVDL3), 921 * WORK(KLAMP0),ISYMD,D,WORK(KEND4),LWRK4) 922C 923C----------------------------------------------------------------------- 924C Construct virtual integrals (for fixed D) which are required 925C to calculate t3bar_X multipliers 926C----------------------------------------------------------------------- 927C 928 IF ( (LISTL(1:3).EQ.'L1 ') 929 * .OR. (LISTL(1:3).EQ.'M1 ') 930 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 931 LSKIPL1R = .FALSE. 932 CALL INTVIR_T3BARX_D(LSKIPL1R, 933 * ISYMOP,LU3VI,FN3VI,LU3VI2,FN3VI2, 934 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2, 935 * WORK(KW3BXVDGX1),WORK(KW3BXVDG1), 936 * WORK(KW3BXVDGX2),WORK(KW3BXVDG2), 937 * WORK(KW3BXVDLX1),WORK(KW3BXVDL1), 938 * WORK(KW3BXVDLX2),WORK(KW3BXVDL2), 939 * WORK(KLAMPL1R),ISYML1R,WORK(KLAMP0), 940 * ISYM0,ISYMD,D,WORK(KEND4),LWRK4) 941 ELSE IF (LISTL(1:3).EQ.'LE ') THEN 942 LSKIPL1R = .TRUE. 943 CALL INTVIR_T3BARX_D(LSKIPL1R, 944 * ISYMOP,LU3VI,FN3VI,LU3VI2,FN3VI2, 945 * LU3FOP,FN3FOP,LU3FOP2,FN3FOP2, 946 * DUMMY,WORK(KW3BXVDG1), 947 * DUMMY,WORK(KW3BXVDG2), 948 * DUMMY,WORK(KW3BXVDL1), 949 * DUMMY,WORK(KW3BXVDL2), 950 * DUMMY,IDUMMY,WORK(KLAMP0), 951 * ISYM0,ISYMD,D,WORK(KEND4),LWRK4) 952 END IF 953C 954C--------------------------------------------------------------------- 955C Read and sort R1-transformed integrals used in contraction. 956C--------------------------------------------------------------------- 957C 958 IF (NCKATR(ISYCKBDR1) .GT. 0) THEN 959 IOFF = ICKBD(ISYCKBDR1,ISYMD) + NCKATR(ISYCKBDR1)*(D - 1) 960 * + 1 961 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVIR),IOFF, 962 & NCKATR(ISYCKBDR1)) 963 ENDIF 964C 965 IF (LWRK4 .LT. NCKATR(ISYCKBDR1)) THEN 966 CALL QUIT('Insufficient space for allocation in '// 967 & 'CC3_BFMAT (TRVI)') 968 END IF 969C 970 CALL CCSDT_SRVIR3(WORK(KTRVIR),WORK(KEND4),ISYMD,D,ISINT1R1) 971C 972C 973 IF (NCKATR(ISYCKBDR1) .GT. 0) THEN 974 IOFF = ICKBD(ISYCKBDR1,ISYMD) + NCKATR(ISYCKBDR1)*(D - 1) 975 * + 1 976 CALL GETWA2(LUDKBCR4,FNDKBCR4,WORK(KTRVIR1),IOFF, 977 & NCKATR(ISYCKBDR1)) 978 ENDIF 979C 980 IF (LWRK4 .LT. NCKATR(ISYCKBDR1)) THEN 981 CALL QUIT('Insufficient space for allocation in '// 982 & 'CC3_BFMAT (TRVI1)') 983 END IF 984C 985 CALL CCSDT_SRVIR3(WORK(KTRVIR1),WORK(KEND4),ISYMD,D, 986 * ISINT1R1) 987C 988 989C 990 DO ISYMB = 1,NSYM 991C 992 ISYALJB0 = MULD2H(ISYMB,ISYML0) 993 ISYALJD0 = MULD2H(ISYMD,ISYML0) 994 ISYMBD = MULD2H(ISYMD,ISYMB) 995 ISCKIJ = MULD2H(ISYMBD,ISYM0) 996 ISYCKD = MULD2H(ISYM0,ISYMB) 997 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 998 * .OR. (LISTL(1:3).EQ.'M1 ') 999 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1000 ISYALJBL1 = MULD2H(ISYMB,ISYML1) 1001 ISYALJDL1 = MULD2H(ISYMD,ISYML1) 1002 ISWBMAT = MULD2H(ISCKIJ,ISYML1) 1003 END IF 1004C 1005 KSMAT2 = KEND4 1006 KUMAT2 = KSMAT2 + NCKIJ(ISCKIJ) 1007 KDIAG = KUMAT2 + NCKIJ(ISCKIJ) 1008 KINDSQ = KDIAG + NCKIJ(ISCKIJ) 1009 KEND5 = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 1010 1011 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 1012 * .OR. (LISTL(1:3).EQ.'M1 ') 1013 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1014 KDIAGWB = KEND5 1015 KINDSQWB = KDIAGWB + NCKIJ(ISWBMAT) 1016 KEND5 = KINDSQWB + (6*NCKIJ(ISWBMAT) - 1)/IRAT + 1 1017 END IF 1018 1019 KINDEX = KEND5 1020 KINDEX2 = KINDEX + (NCKI(ISYALJB0) - 1)/IRAT + 1 1021 KEND5 = KINDEX2 + (NCKI(ISYALJD0) - 1)/IRAT + 1 1022C 1023 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 1024 * .OR. (LISTL(1:3).EQ.'M1 ') 1025 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1026 KINDEXBL1 = KEND5 1027 KINDEXDL1 = KINDEXBL1 + (NCKI(ISYALJBL1)-1)/IRAT + 1 1028 KTMAT = KINDEXDL1 + (NCKI(ISYALJDL1)-1)/IRAT + 1 1029 KW3BMAT = KTMAT + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWBMAT)) 1030 KEND5 = KW3BMAT + NCKIJ(ISWBMAT) 1031 LWRK5 = LWORK - KEND5 1032 ELSE 1033 KTMAT = KEND5 1034 KEND5 = KTMAT + NCKIJ(ISCKIJ) 1035 END IF 1036 1037 KT3BVBG1 = KEND5 1038 KT3BVBG2 = KT3BVBG1 + NCKATR(ISYCKD) 1039 KT3BVBG3 = KT3BVBG2 + NCKATR(ISYCKD) 1040 KEND5 = KT3BVBG3 + NCKATR(ISYCKD) 1041 LWRK5 = LWORK - KEND5 1042 1043 KSMAT4 = KEND5 1044 KUMAT4 = KSMAT4 + NCKIJ(ISCKIJ) 1045 KT3BVBL1 = KUMAT4 + NCKIJ(ISCKIJ) 1046 KT3BVBL2 = KT3BVBL1 + NCKATR(ISYCKD) 1047 KT3BVBL3 = KT3BVBL2 + NCKATR(ISYCKD) 1048 KEND5 = KT3BVBL3 + NCKATR(ISYCKD) 1049 LWRK5 = LWORK - KEND5 1050C 1051 IF (LWRK5 .LT. 0) THEN 1052 WRITE(LUPRI,*) 'Memory available : ',LWORK 1053 WRITE(LUPRI,*) 'Memory needed : ',KEND5 1054 CALL QUIT('Insufficient space(isymb-loop) 1055 * in CC3_BFMAT') 1056 END IF 1057C 1058C-------------------------------------------- 1059C Construct part of the diagonal 1060C-------------------------------------------- 1061C 1062 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 1063 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 1064 * .OR. (LISTL(1:3).EQ.'M1 ') 1065 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1066 CALL CC3_DIAG(WORK(KDIAGWB),WORK(KFOCKD),ISWBMAT) 1067 END IF 1068C 1069C------------------------------------ 1070C Construct index arrays 1071C------------------------------------ 1072C 1073 LENSQ = NCKIJ(ISCKIJ) 1074 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 1075 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 1076 * .OR. (LISTL(1:3).EQ.'M1 ') 1077 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1078 LENSQWB = NCKIJ(ISWBMAT) 1079 CALL CC3_INDSQ(WORK(KINDSQWB),LENSQWB,ISWBMAT) 1080 END IF 1081C 1082 CALL CC3_INDEX(WORK(KINDEX),ISYALJB0) 1083 CALL CC3_INDEX(WORK(KINDEX2),ISYALJD0) 1084 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 1085 * .OR. (LISTL(1:3).EQ.'M1 ') 1086 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1087 CALL CC3_INDEX(WORK(KINDEXBL1),ISYALJBL1) 1088 CALL CC3_INDEX(WORK(KINDEXDL1),ISYALJDL1) 1089 END IF 1090C 1091 DO B = 1,NVIR(ISYMB) 1092C 1093 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'LE ') 1094 * .OR. (LISTL(1:3).EQ.'M1 ') 1095 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1096 CALL DZERO(WORK(KW3BMAT),NCKIJ(ISWBMAT)) 1097 END IF 1098C 1099C----------------------------------------------------------------------- 1100C Construct virtual integrals (for fixed B) which are 1101C required to calculate t3bar_0 multipliers 1102C (the same routine as in d-loop is used) 1103C----------------------------------------------------------------------- 1104C 1105 CALL INTVIR_T3BAR0_D(LU3FOPX,FN3FOPX,LU3FOP2X, 1106 * FN3FOP2X,LUDKBC3,FNDKBC3, 1107 * LU3VI,FN3VI,ISYM0,WORK(KT3BVBL1), 1108 * WORK(KT3BVBG1),WORK(KT3BVBG2), 1109 * WORK(KT3BVBL2),WORK(KT3BVBG3), 1110 * WORK(KT3BVBL3),WORK(KLAMP0), 1111 * ISYMB,B,WORK(KEND5),LWRK5) 1112C 1113C--------------------------------------------------------- 1114C Get T3bar_BD multipliers (using S and U) 1115C--------------------------------------------------------- 1116C 1117 IF ( (LISTL(1:3).EQ.'L1 ') 1118 * .OR. (LISTL(1:3).EQ.'L0 ') ) THEN 1119 CALL GET_T3BAR0_BD(ISYM0,WORK(KL1AM),ISYML0, 1120 * WORK(KL2TP),ISYML0,WORK(KTMAT), 1121 * WORK(KFOCK0CK),WORK(KFOCKD), 1122 * WORK(KDIAG),WORK(KXIAJB),ISYM0, 1123 * ISYM0,WORK(KINDSQ),LENSQ, 1124 * WORK(KSMAT2),WORK(KT3BVDG1), 1125 * WORK(KT3BVDG2),WORK(KT3BVDL1), 1126 * WORK(KT3BVDL2),WORK(KT3BOG1), 1127 * WORK(KT3BOL1),WORK(KINDEX), 1128 * WORK(KSMAT4),WORK(KT3BVBG1), 1129 * WORK(KT3BVBG2),WORK(KT3BVBL1), 1130 * WORK(KT3BVBL2),WORK(KINDEX2), 1131 * WORK(KUMAT2),WORK(KT3BVDG3), 1132 * WORK(KT3BVDL3),WORK(KT3BOG2), 1133 * WORK(KT3BOL2),WORK(KUMAT4), 1134 * WORK(KT3BVBG3),WORK(KT3BVBL3), 1135 * ISYMB,B,ISYMD,D,ISCKIJ, 1136 * WORK(KEND5),LWRK5) 1137 END IF 1138c 1139* call sum_pt3(work(KTMAT),isymb,b,isymd,d, 1140* * ISYM0,work(kx3am),4) 1141C 1142C---------------------------------------------------- 1143C Get T3barX_BD multipliers (using W) 1144C---------------------------------------------------- 1145C 1146 IF (LISTL(1:3).EQ.'L1 ') THEN 1147 CALL GET_T3BARX_BD(.FALSE., 1148 * WORK(KTMAT),ISCKIJ,WORK(KFOCKL1), 1149 * ISYML1,WORK(KW3BMAT),ISWBMAT, 1150 * WORK(KL2TP),ISYML0,WORK(KFOCKL1RCK), 1151 * ISYFCKL1R,WORK(KW3BXVDLX2), 1152 * WORK(KW3BXVDLX1),WORK(KW3BXVDGX2), 1153 * WORK(KW3BXVDGX1),WORK(KW3BXOLX1), 1154 * WORK(KW3BXOGX1),ISINT2L1R, 1155 * WORK(KINDEX),WORK(KINDEX2), 1156 * WORK(KINDSQWB),LENSQWB,WORK(KL2L1), 1157 * ISYML1,WORK(KFOCK0CK),ISYM0, 1158 * WORK(KW3BXVDL2),WORK(KW3BXVDL1), 1159 * WORK(KW3BXVDG2),WORK(KW3BXVDG1), 1160 * WORK(KW3BXOL1),WORK(KW3BXOG1), 1161 * ISINT2,WORK(KINDEXBL1), 1162 * WORK(KINDEXDL1),WORK(KL1L1),ISYML1, 1163 * WORK(KXIAJB),ISINT1,-FREQL1, 1164 * WORK(KDIAGWB),WORK(KFOCKD),B,ISYMB, 1165 * D,ISYMD,ISYML1,WORK(KEND5),LWRK5) 1166c 1167* call sum_pt3(work(KW3BMAT),isymb,b,isymd,d, 1168* * ISWBMAT,work(kx3am),4) 1169 ELSE IF ( (LISTL(1:3).EQ.'M1 ') 1170 * .OR. (LISTL(1:3).EQ.'N2 ') ) THEN 1171 CALL GET_M3BAR_BD(WORK(KTMAT),WORK(KW3BMAT), 1172 * ISWBMAT,WORK(KL2TP),ISYML0,WORK(KFOCKL1RCK), 1173 * ISYFCKL1R,WORK(KW3BXVDLX2), 1174 * WORK(KW3BXVDLX1),WORK(KW3BXVDGX2), 1175 * WORK(KW3BXVDGX1),WORK(KW3BXOLX1), 1176 * WORK(KW3BXOGX1),ISINT2L1R, 1177 * WORK(KINDEX),WORK(KINDEX2), 1178 * WORK(KINDSQWB),LENSQWB,WORK(KL2L1), 1179 * ISYML1,WORK(KFOCK0CK),ISYM0, 1180 * WORK(KW3BXVDL2),WORK(KW3BXVDL1), 1181 * WORK(KW3BXVDG2),WORK(KW3BXVDG1), 1182 * WORK(KW3BXOL1),WORK(KW3BXOG1), 1183 * ISINT2,WORK(KINDEXBL1), 1184 * WORK(KINDEXDL1),WORK(KL1L1),ISYML1, 1185 * WORK(KXIAJB),ISINT1,-FREQL1, 1186 * WORK(KDIAGWB),WORK(KFOCKD),B,ISYMB, 1187 * D,ISYMD,ISYML1,WORK(KEND5),LWRK5) 1188 1189c 1190c call sum_pt3(work(KW3BMAT),isymb,b,isymd,d, 1191c * ISWBMAT,work(kx3am),4) 1192 ELSE IF (LISTL(1:3).EQ.'LE ') THEN 1193 CALL GET_L3BAR_BD(WORK(KTMAT),WORK(KW3BMAT), 1194 * ISWBMAT, 1195 * WORK(KINDSQWB),LENSQWB,WORK(KL2L1), 1196 * ISYML1,WORK(KFOCK0CK),ISYM0, 1197 * WORK(KW3BXVDL2),WORK(KW3BXVDL1), 1198 * WORK(KW3BXVDG2),WORK(KW3BXVDG1), 1199 * WORK(KW3BXOL1),WORK(KW3BXOG1), 1200 * ISINT2,WORK(KINDEXBL1), 1201 * WORK(KINDEXDL1),WORK(KL1L1),ISYML1, 1202 * WORK(KXIAJB),ISINT1,-FREQL1, 1203 * WORK(KDIAGWB),WORK(KFOCKD),B,ISYMB, 1204 * D,ISYMD,ISYML1,WORK(KEND5),LWRK5) 1205c 1206c call sum_pt3(work(KW3BMAT),isymb,b,isymd,d, 1207c * ISWBMAT,work(kx3am),4) 1208 END IF 1209C 1210C 1211C------------------------------------------------------------------- 1212C Set up variables depending on the LISTL (L1 or L0) 1213C------------------------------------------------------------------- 1214C 1215C If WB3X = .TRUE. (i.e, LISTL = L1), then KTB contains wbar_X 1216C else 1217C we get tbar_0 in KTB 1218C 1219C Note Wbar_X and Tbar_0 is stored in the same way for fixed B and D 1220C 1221C Inside the routines CC3_W3_OMEGA1, CC3_W3_CY2V, CC3_W3_CY2O terms 1222C are selected depending on WB3X 1223 1224 IF ( (LISTL(1:3).EQ.'L1 ') .OR. (LISTL(1:3).EQ.'M1 ') 1225 * .OR. (LISTL(1:3).EQ.'N2 ') 1226 * .OR. (LISTL(1:3).EQ.'LE ') ) THEN 1227 WB3X = .TRUE. 1228 ISTBD = ISWBMAT 1229 KTB = KW3BMAT 1230 LENSQTB = LENSQWB 1231 KINDSQTB = KINDSQWB 1232 ELSE 1233 WB3X = .FALSE. 1234 ISTBD = ISCKIJ 1235 !KTMAT will be destroyed so an extra array for 1236 !zeroth-order multipliers is needed 1237 KTB = KEND5 1238 KEND5 = KTB + NCKIJ(ISCKIJ) 1239 LWRK5 = LWORK - KEND5 1240 IF (LWRK5 .LT. 0) THEN 1241 WRITE(LUPRI,*) 'Memory available : ',LWORK 1242 WRITE(LUPRI,*) 'Memory needed : ',KEND5 1243 CALL QUIT('Too little space(setup) 1244 * in CC3_BFMAT') 1245 END IF 1246 CALL DCOPY(NCKIJ(ISCKIJ),WORK(KTMAT),1,WORK(KTB),1) 1247 1248 LENSQTB = LENSQ 1249 KINDSQTB = KINDSQ 1250 END IF 1251 ISTB = MULD2H(ISTBD,ISYMBD) 1252C 1253 !Check the consistency of logical flags 1254 IF (WB3X.AND.(.NOT.SKIPGEI)) THEN 1255 CONTINUE 1256 ELSE IF ((.NOT.WB3X).AND.SKIPGEI) THEN 1257 CONTINUE 1258 ELSE 1259 WRITE(LUPRI,*)'WB3X = ',WB3X 1260 WRITE(LUPRI,*)'SKIPGEI = ',SKIPGEI 1261 WRITE(LUPRI,*)'WB3X and SKIPGEI should be opposite' 1262 CALL QUIT('Inconsistent flags in CC3_BFMAT') 1263 END IF 1264C 1265C------------------------------------------------- 1266C Calculate <L3|[H^B,\tau_nu_2]|HF> 1267C------------------------------------------------- 1268C 1269 CALL CC3_W3_CY2V(OMEGA2EFF,ISYRES,WORK(KRBJIA), 1270 * WORK(KTB),ISTBD,WORK(KTMAT), 1271 * WORK(KTRVIR),WORK(KTRVIR1),ISINT1R1, 1272 * WORK(KEND5),LWRK5,WORK(KINDSQTB), 1273 * LENSQTB,ISYMB,B,ISYMD,D,WB3X) 1274C 1275 CALL CC3_W3_CY2O(OMEGA2EFF,ISYRES,WORK(KTB), 1276 * ISTBD, 1277 * WORK(KTMAT),WORK(KTROCR), 1278 * WORK(KTROCR1),ISINT1R1, 1279 * WORK(KEND5),LWRK5,WORK(KINDSQTB), 1280 * LENSQTB,ISYMB,B,ISYMD,D,WB3X) 1281C 1282C 1283 IF (LVVVV) THEN 1284 !Calculate <L3|[H^0,T^{B}_{2}],tau_{1}]|HF> 1285 CALL CC3_W3_OMEGA1(OMEGA1EFF,ISYRES,WORK(KTB), 1286 * WORK(KTMAT),ISTB, 1287 * WORK(KOIOOOO),WORK(KOIOVVO), 1288 * WORK(KOIOOVV),WORK(KVVVVO), 1289 * ISYM0,WORK(KT2R1),ISYMR1, 1290 * WORK(KEND5),LWRK5, 1291 * LENSQTB,WORK(KINDSQTB), 1292 * ISYMB,B,ISYMD,D,WB3X) 1293C 1294 !Calculate <L3|[H^B,T^{0}_{2}],tau_{1}]|HF> 1295 CALL CC3_W3_OMEGA1(OMEGA1EFF,ISYRES,WORK(KTB), 1296 * WORK(KTMAT),ISTB, 1297 * WORK(KBIOOOO),WORK(KBIOVVO), 1298 * WORK(KBIOOVV),WORK(KVVVVB), 1299 * ISINT1R1,WORK(KT2TP),ISYM0, 1300 * WORK(KEND5),LWRK5, 1301 * LENSQTB,WORK(KINDSQTB), 1302 * ISYMB,B,ISYMD,D,WB3X) 1303 ELSE 1304C 1305 CALL DSCAL(NCKIJ(ISTBD),-1.0D0,WORK(KTB),1) 1306 1307 !<L3|[H^0,T^{B}_{2}],tau_{1}]|HF> in terms of N1 and N2 1308 CALL WT2_N1N2(WORK(KTB),ISTB, 1309 * WORK(KT2R1),ISYMR1, 1310 * WORK(KGEIB),WORK(KFEIB), 1311 * ISYMN1B, 1312 * WORK(KN2MATB),ISYMN2B, 1313 * B,ISYMB,D,ISYMD, 1314 * WORK(KINDSQTB),LENSQTB, 1315 * WORK(KINDSQNB),LENSQNB, 1316 * WORK(KEND5),LWRK5, 1317 * WB3X) 1318C 1319 !<L3|[H^B,T^{0}_{2}],tau_{1}]|HF>in terms of N1 and N2 1320 CALL WT2_N1N2(WORK(KTB),ISTB, 1321 * WORK(KT2TP),ISYM0, 1322 * WORK(KGEI),WORK(KFEI), 1323 * ISYMN1, 1324 * WORK(KN2MAT),ISYMN2, 1325 * B,ISYMB,D,ISYMD, 1326 * WORK(KINDSQTB),LENSQTB, 1327 * WORK(KINDSQN),LENSQN, 1328 * WORK(KEND5),LWRK5, 1329 * WB3X) 1330 1331 END IF 1332 1333C 1334* call sum_pt3(work(ktb),isymb,b,isymd,d, 1335* * ISYM0,work(kx3am),5) 1336 ENDDO ! B 1337 ENDDO ! ISYMB 1338C 1339 IF (.NOT.LVVVV) THEN 1340C 1341C ---------------------------------------------------------- 1342C Put KGEI(ge,i)^F and KFEI(fe,i)^G (which are intermediates 1343C for N1MAT(fge,i) ) to files (for fixed F=D and G=D). 1344C ---------------------------------------------------------- 1345 1346 IF (.NOT.SKIPGEI) THEN 1347 !Put KGEI to file as (gei,F) (fixed F corresponds to D) 1348 IADR = ICKBD(ISGEI,ISYMD) + NCKATR(ISGEI)*(D-1) + 1 1349 CALL PUTWA2(LUGEI,FNGEI,WORK(KGEI),IADR,NCKATR(ISGEI)) 1350 END IF 1351C 1352 !Put KFEI to file as (fei,G) (fixed G corresponds to D) 1353 IADR = ICKBD(ISFEI,ISYMD) + NCKATR(ISFEI)*(D-1) + 1 1354 CALL PUTWA2(LUFEI,FNFEI,WORK(KFEI),IADR,NCKATR(ISFEI)) 1355C 1356 IF (.NOT.SKIPGEI) THEN 1357 !Put KGEIB to file as (gei,F) (fixed F corresponds to D) 1358 IADR = ICKBD(ISGEIB,ISYMD) + NCKATR(ISGEIB)*(D-1) + 1 1359 CALL PUTWA2(LUGEIB,FNGEIB,WORK(KGEIB),IADR,NCKATR(ISGEIB)) 1360 END IF 1361C 1362 !Put KFEIB to file as (fei,G) (fixed G corresponds to D) 1363 IADR = ICKBD(ISFEIB,ISYMD) + NCKATR(ISFEIB)*(D-1) + 1 1364 CALL PUTWA2(LUFEIB,FNFEIB,WORK(KFEIB),IADR,NCKATR(ISFEIB)) 1365C 1366 END IF 1367C 1368 ENDDO ! D 1369 ENDDO ! ISYMD 1370C 1371C------------------------------------------------------ 1372C Accumulate RBJIA from the virtual contribution 1373C in OMEGA2 1374C------------------------------------------------------ 1375C 1376 CALL CC3_RBJIA(OMEGA2EFF,ISYRES,WORK(KRBJIA)) 1377C 1378C 1379 IF (.NOT.LVVVV) THEN 1380C 1381 KLAMPB = KEND0 1382 KLAMHB = KLAMPB + NGLMDT(ISYMR1) 1383 KEND0 = KLAMHB + NGLMDT(ISYMR1) 1384 LWRK0 = LWORK - KEND0 1385 IF (LWRK0 .LT. 0) THEN 1386 WRITE(LUPRI,*) 'Memory available : ',LWORK 1387 WRITE(LUPRI,*) 'Memory needed : ',KEND0 1388 CALL QUIT('Insufficient space in CC3_BFMAT(final)') 1389 END IF 1390C 1391 !Lambda^B (particle) needed in backtransformation of N1 1392 CALL GET_LAMBDAX(WORK(KLAMPB),WORK(KLAMHB),LISTR,IDLSTR,ISYMR1, 1393 * WORK(KLAMP0),WORK(KLAMH0), 1394 * WORK(KEND0),LWRK0) 1395C 1396 !Read (gei,F) and (fei,G) intermediates from files 1397 !add them and put the result to a file as (fge,I) 1398 CALL N1_RESORT(ISYMN1,LUN1,FNN1,LUGEI,FNGEI,LUFEI,FNFEI, 1399 * WORK(KEND0),LWRK0,SKIPGEI) 1400C 1401 !Calculate <T3|[[H^B,T2],tau_ai]|HF> except VVVV contribution 1402 CALL N1N2_G(LUN1,FNN1, 1403 * ISYMN1, 1404 * WORK(KN2MAT),ISYMN2, 1405 * WORK(KBIOVVO),WORK(KBIOOVV), 1406 * WORK(KBIOOOO),ISINT1R1, 1407 * OMEGA1EFF,ISYRES, 1408 * WORK(KINDSQN),LENSQN, 1409 * WORK(KEND0),LWRK0) 1410C 1411 !Calculate VVVV contribution to <T3|[[H^B,T2],tau_ai]|HF> 1412 IOPT = 1 !T1 one-index backtransformation 1413 CALL N1_GV4(IOPT, 1414 * LUN1,FNN1, 1415 * ISYMN1, 1416 * WORK(KLAMPB),ISYMR1, 1417 * WORK(KLAMP0),1, 1418 * WORK(KLAMH0),1, 1419 * WORK(KLAMH0),1, 1420 * OMEGA1EFF,ISYRES, 1421 * WORK(KEND0),LWRK0) 1422C 1423C REPEAT THINGS FOR N1MATB 1424C 1425 !Read (gei,F) and (fei,G) intermediates from files 1426 !add them and put the result to a file as (fge,I) 1427 CALL N1_RESORT(ISYMN1B,LUN1B,FNN1B,LUGEIB,FNGEIB,LUFEIB,FNFEIB, 1428 * WORK(KEND0),LWRK0,SKIPGEI) 1429C 1430 !Calculate <T3|[[H^0,T2^B],tau_ai]|HF> except VVVV contribution 1431 CALL N1N2_G(LUN1B,FNN1B, 1432 * ISYMN1B, 1433 * WORK(KN2MATB),ISYMN2B, 1434 * WORK(KOIOVVO),WORK(KOIOOVV), 1435 * WORK(KOIOOOO),ISYM0, 1436 * OMEGA1EFF,ISYRES, 1437 * WORK(KINDSQNB),LENSQNB, 1438 * WORK(KEND0),LWRK0) 1439C 1440 !Calculate VVVV contribution to <T3|[[H^0,T2^B],tau_ai]|HF> 1441 IOPT = 0 !normal Lambda matrices used in backtransformation 1442 CALL N1_GV4(IOPT, 1443 * LUN1B,FNN1B, 1444 * ISYMN1B, 1445 * WORK(KLAMP0),1, 1446 * WORK(KLAMP0),1, 1447 * WORK(KLAMH0),1, 1448 * WORK(KLAMH0),1, 1449 * OMEGA1EFF,ISYRES, 1450 * WORK(KEND0),LWRK0) 1451C 1452 END IF 1453C 1454 DO I = 1,NT1AM(ISYRES) 1455 OMEGA1(I) = OMEGA1EFF(I) + OMEGA1(I) 1456 END DO 1457 1458 DO I = 1,NT2AM(ISYRES) 1459 OMEGA2(I) = OMEGA2EFF(I) + OMEGA2(I) 1460 END DO 1461C 1462c write(lupri,*)'AFTER ' 1463c write(lupri,*)'omega1eff after CC3_BFMAT,LISTL ', LISTL 1464c call PRINT_MATAI(OMEGA1EFF,ISYRES) 1465c xnormval = ddot(NT1AM(ISYRES),OMEGA1EFF,1,OMEGA1EFF,1) 1466c write(lupri,*)'norm omega1eff after CC3_BFMAT ', xnormval 1467c xnormval = ddot(NT2AM(ISYRES),OMEGA2EFF,1,OMEGA2EFF,1) 1468c write(lupri,*)'norm omega2eff after CC3_BFMAT ' 1469c CALL OUTPAK(OMEGA2EFF,NT1AMX,1,LUPRI) 1470C 1471c write(lupri,*) 't3barx in CC3_BFMAT' 1472c call print_pt3(work(kx3am),ISYM0,4) 1473C 1474C 1475C 1476C------------------------------- 1477C Close and delete files 1478C------------------------------- 1479C 1480 IF (LVVVV) THEN 1481 CALL WCLOSE2(LU4V,FN4V,'DELETE') 1482 CALL WCLOSE2(LU4VB,FN4VB,'DELETE') 1483 END IF 1484C 1485 IF (.NOT.LVVVV) THEN 1486 !Close files for N1MATB intermediates 1487 CALL WCLOSE2(LUFEIB,FNFEIB,'DELETE') 1488 CALL WCLOSE2(LUN1B,FNN1B,'DELETE') 1489 IF (.NOT.SKIPGEI) THEN 1490 CALL WCLOSE2(LUGEIB,FNGEIB,'DELETE') 1491 END IF 1492 !Close files for N1MAT intermediates 1493 CALL WCLOSE2(LUFEI,FNFEI,'DELETE') 1494 CALL WCLOSE2(LUN1,FNN1,'DELETE') 1495 IF (.NOT.SKIPGEI) THEN 1496 CALL WCLOSE2(LUGEI,FNGEI,'DELETE') 1497 END IF 1498 END IF 1499C 1500 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 1501 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 1502 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 1503 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 1504 CALL WCLOSE2(LUDKBCR4,FNDKBCR4,'DELETE') 1505 1506C------------- 1507C End 1508C------------- 1509C 1510 CALL QEXIT('BFMAT') 1511C 1512 RETURN 1513 END 1514 1515 1516