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