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_lhtr_l3 */ 20 SUBROUTINE CC3_FMAT(LISTL,IDLSTL,LISTB,IDLSTB,IOPTRES, 21 * OMEGA1,OMEGA2,OMEGA1EFF,OMEGA2EFF, 22 * IDOTS,DOTPROD,LISTDP,ITRAN, 23 * NFTRAN,MXVEC,WORK,LWORK, 24 * LUCKJD,FNCKJD,LUDKBC,FNDKBC,LUTOC,FNTOC, 25 * LU3VI,FN3VI,LUDKBC3,FNDKBC3, 26 * LU3FOPX,FN3FOPX,LU3FOP2X,FN3FOP2X) 27C 28C Written by K. Hald, Spring 2002. 29C 30C 31 IMPLICIT NONE 32C 33#include "priunit.h" 34#include "dummy.h" 35#include "ccsdsym.h" 36#include "inftap.h" 37#include "ccsdinp.h" 38#include "ccorb.h" 39#include "iratdef.h" 40#include "ccinftap.h" 41#include "second.h" 42C 43 INTEGER IDLSTL, IDLSTB, IOPTRES, ITRAN, NFTRAN, MXVEC, LWORK 44 INTEGER LUCKJD, LUDKBC, LUTOC, LU3VI, LUDKBC3, LU3FOPX, LU3FOP2X 45 INTEGER IDOTS(MXVEC,NFTRAN) 46 INTEGER LUDKBC4, LU4V, ISYMTR, ISINT1, ISINT2, ISYMIM 47 INTEGER ISYMT1, ISYMT2, ISYML, ISYML1, ISYML2, KT1AM 48 INTEGER KL1AM, KL2TP, KLAMDP, KLAMDH, IOPT 49 INTEGER KFOCKD, KFCKBA, KEND0, LWRK0, KEND1, LWRK1 50 INTEGER KEND2, LWRK2, KEND3, LWRK3, KEND4, LWRK4 51 INTEGER LUFCK, ISYMC, ISYMK, KOFF1, KOFF2, KTROC, KTROC0, KTROC1 52 INTEGER KTROC2, KXIAJB, KINTOC, LENGTH, ISYOPE, IOPTTCME 53 INTEGER IOFF, ISYMD, ISAIJ1, ISYCKB, ISCKB1, ISCKB2, KTRVI 54 INTEGER KTRVI1, KTRVI2, KTRVI0, KTRVI3, KTRVI4 55 INTEGER KTRVI5, KTRVI6, KINTVI, ISYMB, ISYALJ, ISAIJ2, ISYMBD 56 INTEGER ISCKIJ, KSMAT, KQMAT, KDIAG, KINDSQ, KINDEX, KTMAT 57 INTEGER KRMAT1, KRMAT2, LENSQ, ISYRES 58 INTEGER ISYMB1, ISYMB2, KB1AM, KB2TP, KT2TP 59 INTEGER KVVVVO, KOIOOOO, KOIOVVO, KOIOOVV 60 INTEGER KVVVVB, KBIOOOO, KBIOVVO, KBIOOVV, LU4VB 61 INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR, LUDKBCR4 62 INTEGER ISINTR1 63 INTEGER KTROC3, KTROC4, ISCKBR, KTRVI7, KTRVI8 64c 65 integer kt3am 66c 67C Functions : 68 INTEGER ILSTSYM 69C 70#if defined (SYS_CRAY) 71 REAL OMEGA1(*), OMEGA2(*), OMEGA1EFF(*), OMEGA2EFF(*) 72 REAL DOTPROD(MXVEC,NFTRAN), WORK(LWORK) 73 REAL TITRAN, TISORT, TISMAT, TIQMAT, TICONT, TIOME1 74 REAL HALF, ONE, TWO 75 REAL DTIME, DDOT, XNORM 76#else 77 DOUBLE PRECISION OMEGA1(*), OMEGA2(*), OMEGA1EFF(*), OMEGA2EFF(*) 78 DOUBLE PRECISION DOTPROD(MXVEC,NFTRAN), WORK(LWORK) 79 DOUBLE PRECISION TITRAN, TISORT, TISMAT, TIQMAT, TICONT, TIOME1 80 DOUBLE PRECISION HALF, ONE, TWO 81 DOUBLE PRECISION DTIME, DDOT, XNORM 82#endif 83C 84 CHARACTER*(*) FNCKJD, FNDKBC, FNTOC, FN3VI, FNDKBC3, FN3FOPX 85 CHARACTER*(*) FN3FOP2X 86 CHARACTER*1 CDUMMY 87 CHARACTER*3 LISTL, LISTB, LISTDP 88 CHARACTER*10 MODEL 89 CHARACTER*13 FNDKBC4, FN4V, FN4VB, FN3SRTR, FNCKJDR 90 CHARACTER*13 FNDELDR, FNDKBCR, FNDKBCR4 91C 92 PARAMETER(HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 93C 94 CALL QENTER('CC3_FMAT') 95C 96C---------------------------------------------------- 97C Initialise character strings and open files 98C---------------------------------------------------- 99C 100 CDUMMY = ' ' 101C 102 LUDKBC4 = -1 103 LU4V = -1 104 LU4VB = -1 105 LU3SRTR = -1 106 LUCKJDR = -1 107 LUDELDR = -1 108 LUDKBCR = -1 109 LUDKBCR4 = -1 110C 111 FNDKBC4 = 'CC3_FMAT_TMP1' 112 FN4V = 'CC3_FMAT_TMP2' 113 FN4VB = 'CC3_FMAT_TMP3' 114 FN3SRTR = 'CC3_FMAT_TMP4' 115 FNCKJDR = 'CC3_FMAT_TMP5' 116 FNDELDR = 'CC3_FMAT_TMP6' 117 FNDKBCR = 'CC3_FMAT_TMP7' 118 FNDKBCR4 = 'CC3_FMAT_TMP8' 119C 120 CALL WOPEN2(LUDKBC4,FNDKBC4,64,0) 121 CALL WOPEN2(LU4V,FN4V,64,0) 122 CALL WOPEN2(LU4VB,FN4VB,64,0) 123 CALL WOPEN2(LU3SRTR,FN3SRTR,64,0) 124 CALL WOPEN2(LUCKJDR,FNCKJDR,64,0) 125 CALL WOPEN2(LUDELDR,FNDELDR,64,0) 126 CALL WOPEN2(LUDKBCR,FNDKBCR,64,0) 127 CALL WOPEN2(LUDKBCR4,FNDKBCR4,64,0) 128C 129C------------------------------------------------------------- 130C Set symmetry flags. 131C 132C omega = int1*T2*int2 133C isymres is symmetry of result(omega) 134C isint1 is symmetry of integrals in contraction.(int1) 135C isint2 is symmetry of integrals in the triples equation.(int2) 136C isymim is symmetry of S and Q intermediates.(t2*int2) 137C (sym is for all index of S and Q (cbd,klj) 138C thus cklj=b*d*isymim) 139C------------------------------------------------------------- 140C 141 IPRCC = IPRINT 142 ISYMTR = ISYMOP 143 ISINT1 = ISYMOP 144 ISINT2 = ISYMOP 145 ISYMT1 = 1 146 ISYMT2 = 1 147 ISYML = ILSTSYM(LISTL,IDLSTL) 148 ISYML1 = ISYML 149 ISYML2 = ISYML 150 ISYMB = ILSTSYM(LISTB,IDLSTB) 151 ISYMB1 = ISYMB 152 ISYMB2 = ISYMB 153 ISINTR1 = MULD2H(ISINT1,ISYMB1) 154 ISYMIM = MULD2H(ISYML,ISYMOP) 155 ISYRES = MULD2H(ISYMIM,ISYMB) 156C 157C-------------------- 158C Time variables. 159C-------------------- 160C 161 TITRAN = 0.0D0 162 TISORT = 0.0D0 163 TISMAT = 0.0D0 164 TIQMAT = 0.0D0 165 TICONT = 0.0D0 166 TIOME1 = 0.0D0 167C 168C----------------------------------- 169C Work space allocation 170C----------------------------------- 171C 172 KFOCKD = 1 173 KFCKBA = KFOCKD + NORBTS 174 KT1AM = KFCKBA + N2BST(ISYMOP) 175 KL1AM = KT1AM + NT1AM(ISYMT1) 176 KL2TP = KL1AM + NT1AM(ISYML1) 177 KB1AM = KL2TP + NT2SQ(ISYML2) 178 KB2TP = KB1AM + NT1AM(ISYMB1) 179 KT2TP = KB2TP + NT2SQ(ISYMB2) 180 KLAMDP = KT2TP + NT2SQ(ISYMT2) 181 KLAMDH = KLAMDP + NLAMDT 182 KEND0 = KLAMDH + NLAMDT 183 LWRK0 = LWORK - KEND0 184C 185 if (.false.) then 186 kt3am = kend0 187 kend0 = kt3am + nt1amx*nt1amx*nt1amx 188 lwrk0 = lwork - kend0 189 call dzero(work(kt3am),nt1amx*nt1amx*nt1amx) 190 endif 191C 192 IF (LWRK0 .LT. 0) THEN 193 WRITE(LUPRI,*) 'Memory available : ',LWORK 194 WRITE(LUPRI,*) 'Memory needed : ',KEND0 195 CALL QUIT('Insufficient space in CC3_FMAT') 196 END IF 197C 198C--------------------------------------------------------- 199C Read the zero'th order amplitudes and calculate 200C the zero'th order Lambda matrices 201C--------------------------------------------------------- 202C 203 IOPT = 1 204 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 205C 206 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 207 & WORK(KEND0),LWRK0) 208C 209C----------------------------------------------------------- 210C Resort DKBC -> DKBC4 211C----------------------------------------------------------- 212C 213 CALL CC3_TCME(WORK(KLAMDP),ISINT1,WORK(KEND0),LWRK0,IDUMMY,CDUMMY, 214 * LUDKBC,FNDKBC,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 215 * IDUMMY,CDUMMY,LUDKBC4,FNDKBC4,2) 216C 217C-------------------------------------- 218C Reorder the t2-amplitudes i T2TP. 219C-------------------------------------- 220C 221 IF (LWORK .LT. NT2SQ(1)) THEN 222 CALL QUIT('Not enough memory to construct T2TP in CC3_FMAT') 223 ENDIF 224C 225 IOPT = 2 226 CALL CC_RDRSP('R0',0,1,IOPT,MODEL, 227 * DUMMY,WORK(KT2TP)) 228 CALL CC_T2SQ(WORK(KT2TP),WORK(KEND0),ISYMT2) 229 CALL CC3_T2TP(WORK(KT2TP),WORK(KEND0),ISYMT2) 230C 231 IF (IPRINT .GT. 55) THEN 232 XNORM = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1) 233 WRITE(LUPRI,*) 'Norm of T2TP ',XNORM 234 ENDIF 235C 236C-------------------------------------- 237C Reorder the l2-amplitudes i L2TP. 238C-------------------------------------- 239C 240 IF (LWORK .LT. NT2SQ(ISYML2)) THEN 241 CALL QUIT('Not enough memory to construct L2TP in CC3_FMAT') 242 ENDIF 243C 244 IOPT = 3 245 CALL CC_RDRSP(LISTL,IDLSTL,ISYML2,IOPT,MODEL, 246 * WORK(KL1AM),WORK(KL2TP)) 247 CALL CC_T2SQ(WORK(KL2TP),WORK(KEND0),ISYML2) 248 CALL CC3_T2TP(WORK(KL2TP),WORK(KEND0),ISYML2) 249C 250 IF (IPRINT .GT. 55) THEN 251 XNORM = DDOT(NT2SQ(ISYML2),WORK(KL2TP),1,WORK(KL2TP),1) 252 WRITE(LUPRI,*) 'Norm of L2TP ',XNORM 253 ENDIF 254C 255C-------------------------------------- 256C Reorder the B2-amplitudes i B2TP. 257C-------------------------------------- 258C 259 IF (LWORK .LT. NT2SQ(ISYMB2)) THEN 260 CALL QUIT('Not enough memory to construct B2TP in CC3_FMAT') 261 ENDIF 262C 263 IOPT = 3 264 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB2,IOPT,MODEL, 265 & WORK(KB1AM),WORK(KB2TP)) 266 CALL CCLR_DIASCL(WORK(KB2TP),TWO,ISYMB2) 267 CALL CC_T2SQ(WORK(KB2TP),WORK(KEND0),ISYMB2) 268 CALL CC3_T2TP(WORK(KB2TP),WORK(KEND0),ISYMB2) 269C 270 IF (IPRINT .GT. 55) THEN 271 XNORM = DDOT(NT2SQ(ISYMB2),WORK(KB2TP),1,WORK(KB2TP),1) 272 WRITE(LUPRI,*) 'Norm of B2TP ',XNORM 273 ENDIF 274C 275C------------------------------------------------------ 276C Calculate the (ck|de)-tilde and (ck|lm)-tilde 277C------------------------------------------------------ 278C 279 CALL CC3_BARINT(WORK(KB1AM),ISYMB1,WORK(KLAMDP), 280 * WORK(KLAMDH),WORK(KEND0),LWRK0, 281 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 282C 283 CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISINTR1,LU3SRTR,FN3SRTR, 284 * LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY, 285 * IDUMMY,CDUMMY) 286C 287 CALL CC3_SINT(WORK(KLAMDH),WORK(KEND0),LWRK0,ISINTR1, 288 * LUDELDR,FNDELDR,LUDKBCR,FNDKBCR) 289C 290 CALL CC3_TCME(WORK(KLAMDP),ISINTR1,WORK(KEND0),LWRK0, 291 * IDUMMY,CDUMMY,LUDKBCR,FNDKBCR, 292 * IDUMMY,CDUMMY,IDUMMY,CDUMMY, 293 * IDUMMY,CDUMMY,LUDKBCR4,FNDKBCR4,2) 294C 295C-------------------------------------------------------------- 296C Calculate the normal g^0 integrals for 297C OOOO, OVVO and OOVV integrals. VVVV is stored on file 298C-------------------------------------------------------------- 299C 300 KOIOOOO = KEND0 301 KOIOVVO = KOIOOOO + N3ORHF(ISYMOP) 302 KOIOOVV = KOIOVVO + NT2SQ(ISYMOP) 303 KEND0 = KOIOOVV + NT2SQ(ISYMOP) 304 LWRK0 = LWORK - KEND0 305C 306 IF (LWRK0 .LT. 0) THEN 307 CALL QUIT('Out of memory in CC3_FMAT (g^0[2o2v] kind)') 308 ENDIF 309C 310 CALL DZERO(WORK(KOIOVVO),NT2SQ(ISYMOP)) 311 CALL DZERO(WORK(KOIOOVV),NT2SQ(ISYMOP)) 312C 313 CALL CC3_INT(WORK(KOIOOOO),WORK(KOIOOVV),WORK(KOIOVVO), 314 * ISYMOP,LU4V,FN4V, 315 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 316 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 317 * WORK(KEND0),LWRK0) 318C 319C--------------------------------------------------------------- 320C Calculate the g^B integrals for 321C OOOO, OVVO and OOVV integrals. VVVV is stored on file 322C--------------------------------------------------------------- 323C 324 KBIOOOO = KEND0 325 KBIOVVO = KBIOOOO + N3ORHF(ISINTR1) 326 KBIOOVV = KBIOVVO + NT2SQ(ISINTR1) 327 KEND0 = KBIOOVV + NT2SQ(ISINTR1) 328 LWRK0 = LWORK - KEND0 329C 330 IF (LWRK0 .LT. 0) THEN 331 CALL QUIT('Out of memory in CC3_FMAT (g^B[2o2v] kind)') 332 ENDIF 333C 334 CALL DZERO(WORK(KBIOOOO),N3ORHF(ISINTR1)) 335 CALL DZERO(WORK(KBIOVVO),NT2SQ(ISINTR1)) 336 CALL DZERO(WORK(KBIOOVV),NT2SQ(ISINTR1)) 337C 338 CALL CC3_INT(WORK(KBIOOOO),WORK(KBIOOVV),WORK(KBIOVVO), 339 * ISINTR1,LU4VB,FN4VB, 340 * .TRUE.,LISTB,IDLSTB,ISYMB1, 341 * .FALSE.,'DUMMY',IDUMMY,IDUMMY, 342 * WORK(KEND0),LWRK0) 343C 344C--------------------------------------------------------- 345C Read canonical orbital energies and MO coefficients. 346C--------------------------------------------------------- 347C 348 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 349 & .FALSE.) 350 REWIND LUSIFC 351C 352 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 353 READ (LUSIFC) 354 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 355C 356 CALL GPCLOSE(LUSIFC,'KEEP') 357C 358C--------------------------------------------- 359C Delete frozen orbitals in Fock diagonal. 360C--------------------------------------------- 361C 362 IF (FROIMP .OR. FROEXP) 363 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND0),LWRK0) 364C 365C----------------------------------------------------- 366C Construct the transformed Fock matrix 367C----------------------------------------------------- 368C 369 LUFCK = -1 370C This AO Fock matrix is constructed from the T1 transformed density 371 CALL GPOPEN(LUFCK,'CC_FCKH','UNKNOWN',' ','UNFORMATTED', 372 * IDUMMY,.FALSE.) 373C This AO Fock matrix is constructed from the CMO transformed density 374C CALL GPOPEN(LUFCK,'CC_FCKREF','UNKNOWN',' ','UNFORMATTED', 375C * IDUMMY,.FALSE.) 376 REWIND(LUFCK) 377 READ(LUFCK)(WORK(KFCKBA + I-1),I = 1,N2BST(ISYMOP)) 378 CALL GPCLOSE(LUFCK,'KEEP' ) 379C 380 IF (IPRINT .GT. 140) THEN 381 CALL AROUND( 'Usual Fock AO matrix' ) 382 CALL CC_PRFCKAO(WORK(KFCKBA),ISYMOP) 383 ENDIF 384C 385 CALL CC_FCKMO(WORK(KFCKBA),WORK(KLAMDP),WORK(KLAMDH), 386 * WORK(KEND0),LWRK0,1,1,1) 387C 388 IF (IPRINT .GT. 50) THEN 389 CALL AROUND( 'In CC3_FMAT: Triples Fock MO matrix' ) 390 CALL CC_PRFCKMO(WORK(KFCKBA),ISYMOP) 391 ENDIF 392C 393C Sort the fock matrix 394C 395C 396 CALL DCOPY(N2BST(ISINT1),WORK(KFCKBA),1,WORK(KEND0),1) 397C 398 DO ISYMC = 1,NSYM 399C 400 ISYMK = MULD2H(ISYMC,ISINT1) 401C 402 DO K = 1,NRHF(ISYMK) 403C 404 DO C = 1,NVIR(ISYMC) 405C 406 KOFF1 = KEND0 + IFCVIR(ISYMK,ISYMC) + 407 * NORB(ISYMK)*(C - 1) + K - 1 408 KOFF2 = KFCKBA + IT1AM(ISYMC,ISYMK) 409 * + NVIR(ISYMC)*(K - 1) + C - 1 410C 411 WORK(KOFF2) = WORK(KOFF1) 412C 413 ENDDO 414 ENDDO 415 ENDDO 416C 417 IF (IPRINT .GT. 50) THEN 418 CALL AROUND('In CC3_FMAT: Triples Fock MO matrix (sort)') 419 CALL CC_PRFCKMO(WORK(KFCKBA),ISYMOP) 420 ENDIF 421C 422C----------------------------- 423C Read occupied integrals. 424C----------------------------- 425C 426C Memory allocation. 427C 428 KTROC = KEND0 429 KTROC0 = KTROC + NTRAOC(ISINT2) 430 KTROC1 = KTROC0 + NTRAOC(ISINT2) 431 KTROC2 = KTROC1 + NTRAOC(ISINT2) 432 KTROC3 = KTROC2 + NTRAOC(ISINT2) 433 KTROC4 = KTROC3 + NTRAOC(ISINTR1) 434 KXIAJB = KTROC4 + NTRAOC(ISINTR1) 435 KEND1 = KXIAJB + NT2AM(ISYMOP) 436 LWRK1 = LWORK - KEND1 437C 438 KINTOC = KEND1 439 KEND2 = KINTOC + MAX(NTOTOC(ISYMOP),NTOTOC(ISINT2)) 440 LWRK2 = LWORK - KEND2 441C 442 IF (LWRK2 .LT. 0) THEN 443 WRITE(LUPRI,*) 'Memory available : ',LWORK 444 WRITE(LUPRI,*) 'Memory needed : ',KEND2 445 CALL QUIT('Insufficient space in CC3_LHTR_L3') 446 END IF 447C 448C------------------------ 449C Construct L(ia,jb). 450C------------------------ 451C 452 LENGTH = IRAT*NT2AM(ISYMOP) 453C 454 REWIND(LUIAJB) 455 CALL READI(LUIAJB,LENGTH,WORK(KXIAJB)) 456C 457 ISYOPE = ISYMOP 458 IOPTTCME = 1 459 CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISYOPE,IOPTTCME) 460C 461 IF ( IPRINT .GT. 55) THEN 462 XNORM = DDOT(NT2AM(ISYMOP),WORK(KXIAJB),1, 463 * WORK(KXIAJB),1) 464 WRITE(LUPRI,*) 'Norm of IAJB ',XNORM 465 ENDIF 466C 467C-------------------------------------------------------------------------- 468C Read in B-transformed integrals used in contractions and transform. 469C-------------------------------------------------------------------------- 470C 471 IOFF = 1 472 IF (NTOTOC(ISINTR1) .GT. 0) THEN 473 CALL GETWA2(LUCKJDR,FNCKJDR,WORK(KINTOC),IOFF,NTOTOC(ISINTR1)) 474 ENDIF 475C 476 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC3),WORK(KLAMDP), 477 * WORK(KEND2),LWRK2,ISINTR1) 478C 479 CALL CCFOP_SORT(WORK(KTROC3),WORK(KTROC4),ISINTR1,1) 480C 481 CALL CC3_LSORT1(WORK(KTROC3),ISINTR1,WORK(KEND2),LWRK2,5) 482C 483C------------------------------------------------------------ 484C Read in integrals used in contractions and transform. 485C------------------------------------------------------------ 486C 487 IOFF = 1 488 IF (NTOTOC(ISINT2) .GT. 0) THEN 489 CALL GETWA2(LUCKJD,FNCKJD,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 490 ENDIF 491C 492 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC),WORK(KLAMDP), 493 * WORK(KEND2),LWRK2,ISINT2) 494C 495 CALL CCFOP_SORT(WORK(KTROC),WORK(KTROC1),ISINT2,1) 496C 497 CALL CC3_LSORT1(WORK(KTROC),ISINT2,WORK(KEND2),LWRK2,5) 498C 499C------------------------------------------------------------ 500C Read in integrals used in L3 and transform. 501C------------------------------------------------------------ 502C 503 IOFF = 1 504 IF (NTOTOC(ISINT2) .GT. 0) THEN 505 CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISINT2)) 506 ENDIF 507C 508 CALL CC3_TROCC(WORK(KINTOC),WORK(KTROC0),WORK(KLAMDH), 509 * WORK(KEND2),LWRK2,ISINT2) 510C 511C----------------------------------------------------------- 512C Construct 2*C-E of the integrals. 513C----------------------------------------------------------- 514C 515 CALL CCSDT_TCMEOCC(WORK(KTROC0),WORK(KTROC2),ISINT2) 516C 517C------------------------------- 518C Write out norms of arrays. 519C------------------------------- 520C 521 IF (IPRINT .GT. 55) THEN 522 XNORM = DDOT(NTOTOC(ISINT2),WORK(KINTOC),1, 523 * WORK(KINTOC),1) 524 WRITE(LUPRI,*) 'Norm of CKJDEL-INT ',XNORM 525 ENDIF 526C 527 IF (IPRINT .GT. 55) THEN 528 XNORM = DDOT(NTRAOC(ISINT2),WORK(KTROC0),1, 529 * WORK(KTROC0),1) 530 WRITE(LUPRI,*) 'Norm of TROC0 CC3_FMAT : ',XNORM 531 ENDIF 532C 533 IF (IPRINT .GT. 55) THEN 534 XNORM = DDOT(NTRAOC(ISINT2),WORK(KTROC2),1, 535 * WORK(KTROC2),1) 536 WRITE(LUPRI,*) 'Norm of TROC2 CC3_FMAT : ',XNORM 537 ENDIF 538C 539C---------------------------- 540C General loop structure. 541C---------------------------- 542C 543 DO ISYMD = 1,NSYM 544C 545 ISAIJ1 = MULD2H(ISYMD,ISYRES) 546 ISYCKB = MULD2H(ISYMD,ISYMOP) 547 ISCKB1 = MULD2H(ISINT1,ISYMD) 548 ISCKB2 = MULD2H(ISINT2,ISYMD) 549 ISCKBR = MULD2H(ISINTR1,ISYMD) 550C 551 IF (IPRINT .GT. 55) THEN 552C 553 WRITE(LUPRI,*) 'In CC3_FMAT: ISAIJ1 :',ISAIJ1 554 WRITE(LUPRI,*) 'In CC3_FMAT: ISYCKB :',ISYCKB 555 WRITE(LUPRI,*) 'In CC3_FMAT: ISCKB2 :',ISCKB2 556 WRITE(LUPRI,*) 'In CC3_FMAT: ISCKBR :',ISCKBR 557C 558 ENDIF 559C 560C-------------------------- 561C Memory allocation. 562C-------------------------- 563C 564 KTRVI = KEND1 565 KTRVI1 = KTRVI + NCKATR(ISCKB1) 566 KTRVI7 = KTRVI1 + NCKATR(ISCKB1) 567 KTRVI8 = KTRVI7 + NCKATR(ISCKBR) 568 KTRVI2 = KTRVI8 + NCKATR(ISCKBR) 569 KRMAT1 = KTRVI2 + NCKATR(ISCKB2) 570 KEND2 = KRMAT1 + NCKI(ISAIJ1) 571 LWRK2 = LWORK - KEND2 572C 573 KTRVI0 = KEND2 574 KTRVI3 = KTRVI0 + NCKATR(ISCKB2) 575 KTRVI4 = KTRVI3 + NCKATR(ISCKB2) 576 KTRVI5 = KTRVI4 + NCKATR(ISCKB2) 577 KTRVI6 = KTRVI5 + NCKATR(ISCKB2) 578 KVVVVO = KTRVI6 + NCKATR(ISCKB2) 579 KVVVVB = KVVVVO + NMAABC(ISCKB2) 580 KEND3 = KVVVVB + NMAABC(ISCKBR) 581 LWRK3 = LWORK - KEND3 582C 583 KINTVI = KEND3 584 KEND4 = KINTVI + MAX(NCKA(ISYMD),NCKA(ISCKB2)) 585 LWRK4 = LWORK - KEND4 586C 587 IF (LWRK4 .LT. 0) THEN 588 WRITE(LUPRI,*) 'Memory available : ',LWORK 589 WRITE(LUPRI,*) 'Memory needed : ',KEND4 590 CALL QUIT('Insufficient space in CC3_LHTR_L3') 591 END IF 592C 593C--------------------- 594C Sum over D 595C--------------------- 596C 597 DO D = 1,NVIR(ISYMD) 598C 599C------------------------------------ 600C Initialize the R1 matrix. 601C------------------------------------ 602C 603 CALL DZERO(WORK(KRMAT1),NCKI(ISAIJ1)) 604C 605C-------------------------------------------- 606C Read in g^0_{vvvv} for a given D 607C-------------------------------------------- 608C 609 IF (NMAABC(ISCKB2) .GT. 0) THEN 610 IOFF = I3VVIR(ISCKB2,ISYMD) 611 * + NMAABC(ISCKB2)*(D-1) 612 * + 1 613 CALL GETWA2(LU4V,FN4V,WORK(KVVVVO),IOFF,NMAABC(ISCKB2)) 614 ENDIF 615C 616C-------------------------------------------- 617C Read in g^B_{vvvv} for a given D 618C-------------------------------------------- 619C 620 IF (NMAABC(ISCKBR) .GT. 0) THEN 621 IOFF = I3VVIR(ISCKBR,ISYMD) 622 * + NMAABC(ISCKBR)*(D-1) 623 * + 1 624 CALL GETWA2(LU4VB,FN4VB,WORK(KVVVVB),IOFF,NMAABC(ISCKBR)) 625 ENDIF 626C 627C-------------------------------------------------------------- 628C Read B-transformed integrals used in contraction 629C-------------------------------------------------------------- 630C 631 IF (NCKATR(ISCKBR) .GT. 0) THEN 632 IOFF = ICKBD(ISCKBR,ISYMD) + NCKATR(ISCKBR)*(D - 1) + 1 633 CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI7),IOFF, 634 & NCKATR(ISCKBR)) 635 ENDIF 636C 637 IF (LWRK4 .LT. NCKATR(ISCKBR)) THEN 638 CALL QUIT('Insufficient space for allocation in '// 639 & 'CC3_FMAT (TRVI7)') 640 END IF 641C 642 DTIME = SECOND() 643 CALL CCSDT_SRVIR3(WORK(KTRVI7),WORK(KEND4),ISYMD,D,ISINTR1) 644C 645 IF (NCKATR(ISCKBR) .GT. 0) THEN 646 IOFF = ICKBD(ISCKBR,ISYMD) + NCKATR(ISCKBR)*(D - 1) + 1 647 CALL GETWA2(LUDKBCR4,FNDKBCR4,WORK(KTRVI8),IOFF, 648 & NCKATR(ISCKBR)) 649 ENDIF 650C 651 IF (LWRK4 .LT. NCKATR(ISCKBR)) THEN 652 CALL QUIT('Insufficient space for allocation in '// 653 & 'CC3_FMAT (TRVI8)') 654 END IF 655C 656 DTIME = SECOND() 657 CALL CCSDT_SRVIR3(WORK(KTRVI8),WORK(KEND4),ISYMD,D,ISINTR1) 658C 659C------------------------------------------------------------ 660C Read and transform integrals used in contraction. 661C------------------------------------------------------------ 662C 663 IF (NCKATR(ISCKB1) .GT. 0) THEN 664 IOFF = ICKBD(ISCKB1,ISYMD) + NCKATR(ISCKB1)*(D - 1) + 1 665 CALL GETWA2(LUDKBC,FNDKBC,WORK(KTRVI),IOFF, 666 & NCKATR(ISCKB1)) 667 ENDIF 668C 669 IF (LWRK4 .LT. NCKATR(ISCKB1)) THEN 670 CALL QUIT('Insufficient space for allocation in '// 671 & 'CC3_L3 (TRVI)') 672 END IF 673C 674 DTIME = SECOND() 675 CALL CCSDT_SRVIR3(WORK(KTRVI),WORK(KEND4),ISYMD,D,ISINT1) 676C 677 DTIME = SECOND() - DTIME 678 TISORT = TISORT + DTIME 679C 680 IF (NCKATR(ISCKB1) .GT. 0) THEN 681 IOFF = ICKBD(ISCKB1,ISYMD) + NCKATR(ISCKB1)*(D - 1) + 1 682 CALL GETWA2(LUDKBC4,FNDKBC4,WORK(KTRVI1),IOFF, 683 & NCKATR(ISCKB1)) 684 ENDIF 685C 686 IF (LWRK4 .LT. NCKATR(ISCKB1)) THEN 687 CALL QUIT('Insufficient space for allocation in '// 688 & 'CC3_L3 (TRVI1)') 689 END IF 690C 691 DTIME = SECOND() 692 CALL CCSDT_SRVIR3(WORK(KTRVI1),WORK(KEND4),ISYMD,D,ISINT1) 693C 694 DTIME = SECOND() - DTIME 695 TISORT = TISORT + DTIME 696C 697C----------------------------------------------- 698C Integrals used in L3am. 699C----------------------------------------------- 700C 701 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 702 IF (NCKATR(ISCKB2) .GT. 0) THEN 703 CALL GETWA2(LUDKBC3,FNDKBC3,WORK(KTRVI0),IOFF, 704 & NCKATR(ISCKB2)) 705 ENDIF 706C 707 CALL CCSDT_SRVIR3(WORK(KTRVI0),WORK(KEND4),ISYMD,D,ISINT2) 708C 709C------------------------------------------------------ 710C Read 2*C-E of integral used for L3 711C------------------------------------------------------ 712C 713 IOFF = ICKBD(ISCKB2,ISYMD) + NCKATR(ISCKB2)*(D - 1) + 1 714 IF (NCKATR(ISCKB2) .GT. 0) THEN 715 CALL GETWA2(LU3FOP2X,FN3FOP2X,WORK(KTRVI4),IOFF, 716 & NCKATR(ISCKB2)) 717 ENDIF 718C 719C----------------------------------------------------------- 720C Sort the integrals for s3am and for t3-bar 721C----------------------------------------------------------- 722C 723 DTIME = SECOND() 724 CALL CCSDT_SRTVIR(WORK(KTRVI0),WORK(KTRVI2),WORK(KEND4), 725 * LWRK4,ISYMD,ISINT2) 726C 727 CALL CCSDT_SRTVIR(WORK(KTRVI4),WORK(KTRVI5),WORK(KEND4), 728 * LWRK4,ISYMD,ISINT2) 729C 730 DTIME = SECOND() - DTIME 731 TISORT = TISORT + DTIME 732C 733 IF (IPRINT .GT. 55) THEN 734 XNORM = DDOT(NCKATR(ISCKB2),WORK(KTRVI0),1, 735 * WORK(KTRVI0),1) 736 WRITE(LUPRI,*) 'Norm of TRVI0 ',XNORM 737 ENDIF 738C 739 IF (IPRINT .GT. 55) THEN 740 XNORM = DDOT(NCKATR(ISCKB2),WORK(KTRVI2),1, 741 * WORK(KTRVI2),1) 742 WRITE(LUPRI,*) 'Norm of TRVI2 ',XNORM 743 ENDIF 744C 745 IF (IPRINT .GT. 55) THEN 746 XNORM = DDOT(NCKATR(ISCKB2),WORK(KTRVI4),1, 747 * WORK(KTRVI4),1) 748 WRITE(LUPRI,*) 'Norm of TRVI4 ',XNORM 749 ENDIF 750C 751 IF (IPRINT .GT. 55) THEN 752 XNORM = DDOT(NCKATR(ISCKB2),WORK(KTRVI5),1, 753 * WORK(KTRVI5),1) 754 WRITE(LUPRI,*) 'Norm of TRVI5 ',XNORM 755 ENDIF 756C 757C------------------------------------------------------ 758C Calculate virtual integrals used in q3am. 759C------------------------------------------------------ 760C 761 IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1 762 IF (NCKA(ISYCKB) .GT. 0) THEN 763 CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF, 764 * NCKA(ISYCKB)) 765 ENDIF 766C 767 CALL CCSDT_TRVIR(WORK(KINTVI),WORK(KTRVI3),WORK(KLAMDP), 768 * ISYMD,D,ISYMOP,WORK(KEND4),LWRK4) 769C 770 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 771 CALL QUIT('Insufficient space for allocation in '// 772 * 'CC3_FMAT (TRVI3)') 773 END IF 774C 775C Can use kend3 since dont need the integrals anymore 776 DTIME = SECOND() 777 CALL CCSDT_SRVIR3(WORK(KTRVI3),WORK(KEND3),ISYMD,D,ISINT2) 778C 779C--------------------------------------------------------------- 780C Read virtual integrals used in q3am/u3am for t3-bar. 781C--------------------------------------------------------------- 782C 783 IOFF = ICKBD(ISYCKB,ISYMD) + NCKATR(ISYCKB)*(D - 1) + 1 784 IF (NCKATR(ISYCKB) .GT. 0) THEN 785 CALL GETWA2(LU3FOPX,FN3FOPX,WORK(KTRVI6),IOFF, 786 * NCKATR(ISYCKB)) 787 ENDIF 788C 789 IF (LWRK3 .LT. NCKATR(ISYCKB)) THEN 790 CALL QUIT('Insufficient space for allocation in '// 791 & 'CC3_FMAT (TRVI6)') 792 END IF 793C 794C Can use kend3 since dont need the integrals anymore 795 DTIME = SECOND() 796 CALL CCSDT_SRVIR3(WORK(KTRVI6),WORK(KEND4),ISYMD,D,ISINT2) 797C 798 IF (IPRINT .GT. 55) THEN 799 XNORM = DDOT(NCKATR(ISCKB2),WORK(KTRVI6),1, 800 * WORK(KTRVI6),1) 801 WRITE(LUPRI,*) 'Norm of TRVI6 ',XNORM 802 ENDIF 803C 804C--------------------- 805C Calculate. 806C--------------------- 807C 808 DO ISYMB = 1,NSYM 809C 810 ISYALJ = MULD2H(ISYMB,ISYML2) 811 ISAIJ2 = MULD2H(ISYMB,ISYRES) 812 ISYMBD = MULD2H(ISYMB,ISYMD) 813 ISCKIJ = MULD2H(ISYMBD,ISYMIM) 814C 815 IF (IPRINT .GT. 55) THEN 816C 817 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYMD :',ISYMD 818 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYMB :',ISYMB 819 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYALJ:',ISYALJ 820 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISAIJ2:',ISAIJ2 821 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISYMBD:',ISYMBD 822 WRITE(LUPRI,*) 'In CC3_LHTR_L3: ISCKIJ:',ISCKIJ 823C 824 ENDIF 825C 826C Can use kend3 since we do not need the integrals anymore. 827 KSMAT = KEND3 828 KQMAT = KSMAT + NCKIJ(ISCKIJ) 829 KDIAG = KQMAT + NCKIJ(ISCKIJ) 830 KINDSQ = KDIAG + NCKIJ(ISCKIJ) 831 KINDEX = KINDSQ + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1 832 KTMAT = KINDEX + (NCKI(ISYALJ) - 1)/IRAT + 1 833 KRMAT2 = KTMAT + NCKIJ(ISCKIJ) 834 KEND4 = KRMAT2 + NCKI(ISAIJ2) 835 LWRK4 = LWORK - KEND4 836C 837 IF (LWRK4 .LT. 0) THEN 838 WRITE(LUPRI,*) 'Memory available : ',LWORK 839 WRITE(LUPRI,*) 'Memory needed : ',KEND4 840 CALL QUIT('Insufficient space in CC3_LHTR_L3 (inner)') 841 END IF 842C 843C--------------------------------------------- 844C Construct part of the diagonal. 845C--------------------------------------------- 846C 847 CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ) 848C 849 IF (IPRINT .GT. 55) THEN 850 XNORM = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1, 851 * WORK(KDIAG),1) 852 WRITE(LUPRI,*) 'Norm of DIA ',XNORM 853 ENDIF 854 855C 856C------------------------------------- 857C Construct index arrays. 858C------------------------------------- 859C 860 LENSQ = NCKIJ(ISCKIJ) 861 CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ) 862 CALL CC3_INDEX(WORK(KINDEX),ISYALJ) 863C 864 DO B = 1,NVIR(ISYMB) 865C 866C----------------------------------------- 867C Initialize the R2 matrix. 868C----------------------------------------- 869C 870 CALL DZERO(WORK(KRMAT2),NCKI(ISAIJ2)) 871C 872C--------------------------------------------------------------------------- 873C Calculate the S(ci,bk,dj) matrix for for B,D for L3 874C Scale the smat with a half since have a factor 875C of two in the routine. 876C--------------------------------------------------------------------------- 877C 878 DTIME = SECOND() 879C 880 CALL DZERO(WORK(KSMAT),NCKIJ(ISCKIJ)) 881C 882 CALL CCFOP_SMAT(0.0D0,WORK(KL1AM),ISYML1,WORK(KL2TP), 883 * ISYML2,WORK(KTMAT), 884 * WORK(KFCKBA),WORK(KXIAJB),ISINT1, 885 * WORK(KTRVI0),WORK(KTRVI2), 886 * WORK(KTRVI4),WORK(KTRVI5), 887 * WORK(KTROC0),WORK(KTROC2), 888 * ISINT2,WORK(KFOCKD),WORK(KDIAG), 889 * WORK(KSMAT),WORK(KEND4),LWRK4, 890 * WORK(KINDEX),WORK(KINDSQ),LENSQ, 891 * ISYMB,B,ISYMD,D) 892C 893 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KSMAT),1) 894C 895 CALL T3_FORBIDDEN(WORK(KSMAT),ISYMIM,ISYMB,B,ISYMD,D) 896COMMENT COMMENT 897COMMENT COMMENT 898C call sum_pt3(work(ksmat),isymb,b,isymd,d, 899C * isckij,work(kt3am),1) 900COMMENT COMMENT 901COMMENT COMMENT 902C 903 DTIME = SECOND() - DTIME 904 TISMAT = TISMAT + DTIME 905C 906 IF (IPRINT .GT. 55) THEN 907 XNORM = DDOT(NCKIJ(ISCKIJ),WORK(KSMAT),1, 908 * WORK(KSMAT),1) 909 WRITE(LUPRI,*) 'Norm of SMAT-L3 ',XNORM 910 ENDIF 911C 912C------------------------------------------------------------------- 913C Calculate Q(ci,jk) for fixed b,d for t3-bar. 914C------------------------------------------------------------------- 915C 916 DTIME = SECOND() 917C 918 CALL DZERO(WORK(KQMAT),NCKIJ(ISCKIJ)) 919C 920 ! g integrals 921c call dzero(WORK(KTRVI3),NCKATR(ISCKB2)) 922 ! L integrals 923c call dzero(WORK(KTRVI6),NCKATR(ISCKB2)) 924C 925 ! g integrals 926c call dzero(WORK(KTROC0),NTRAOC(ISINT2)) 927c ! L integrals 928c call dzero(WORK(KTROC2),NTRAOC(ISINT2)) 929c 930c call dzero(WORK(KXIAJB),NT2AM(ISYMOP)) 931c call dzero(WORK(KFCKBA),NT1AM(1)) 932 CALL CCFOP_QMAT(0.0D0,WORK(KL1AM),ISYML1, 933 * WORK(KL2TP),ISYML2, 934 * WORK(KTMAT),WORK(KFCKBA), 935 * WORK(KXIAJB),ISINT1,WORK(KTRVI3), 936 * WORK(KTRVI6),WORK(KTROC0), 937 * WORK(KTROC2),ISINT2,WORK(KFOCKD), 938 * WORK(KDIAG),WORK(KQMAT), 939 * WORK(KEND4),LWRK4,WORK(KINDSQ), 940 * LENSQ,ISYMB,B,ISYMD,D) 941C 942 CALL DSCAL(NCKIJ(ISCKIJ),HALF,WORK(KQMAT),1) 943C 944 CALL T3_FORBIDDEN(WORK(KQMAT),ISYMIM,ISYMB,B,ISYMD,D) 945COMMENT COMMENT COMMENT 946COMMENT COMMENT COMMENT 947c call sum_pt3(work(kqmat),isymb,b,isymd,d, 948c * isckij,work(kt3am),2) 949COMMENT COMMENT COMMENT 950COMMENT COMMENT COMMENT 951C 952 DTIME = SECOND() - DTIME 953 TIQMAT = TIQMAT + DTIME 954C 955 IF (IPRINT .GT. 55) THEN 956 XNORM = DDOT(NCKIJ(ISCKIJ),WORK(KQMAT),1, 957 * WORK(KQMAT),1) 958 WRITE(LUPRI,*) 'Norm of QMAT-L3 ',XNORM 959 ENDIF 960C 961C----------------------------------------------------------------------- 962C Calculate the contributions to omega2 963C----------------------------------------------------------------------- 964C 965 CALL CCFOP_CONVIR(WORK(KRMAT2),WORK(KSMAT), 966 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 967 * WORK(KTRVI7),WORK(KTRVI8),ISINTR1, 968 * WORK(KEND4),LWRK4,WORK(KINDSQ), 969 * LENSQ,ISYMB,B,ISYMD,D) 970C 971 CALL CCFOP_CONOCC(OMEGA2,WORK(KRMAT1), 972 * WORK(KRMAT2),WORK(KSMAT), 973 * WORK(KTMAT),ISYMIM, 974 * WORK(KTROC3),WORK(KTROC4),ISINTR1, 975 * WORK(KEND4),LWRK4,WORK(KINDSQ), 976 * LENSQ,ISYMB,B,ISYMD,D) 977C 978C------------------------------------------------------------------------ 979C Calculate the L3 contribution to omega1 980C------------------------------------------------------------------------ 981C 982 CALL DSCAL(NCKIJ(ISCKIJ),-ONE,WORK(KSMAT),1) 983 CALL DSCAL(NCKIJ(ISCKIJ),-ONE,WORK(KQMAT),1) 984C 985C t^{B}_{1} <L3|[H^0,T^{C}_{2}],tau_{1}]|HF> 986 CALL CC3_L3_OMEGA1(OMEGA1,ISYRES,WORK(KSMAT), 987 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 988 * WORK(KOIOOOO),WORK(KOIOVVO), 989 * WORK(KOIOOVV),WORK(KVVVVO),ISINT1, 990 * WORK(KB2TP),ISYMB2,WORK(KEND4), 991 * LWRK4,LENSQ,WORK(KINDSQ), 992 * ISYMB,B,ISYMD,D) 993C 994 IF (IPRINT .GT. 55) THEN 995 XNORM = DDOT(NT1AM(ISYRES),OMEGA1,1,OMEGA1,1) 996 WRITE(LUPRI,*) 'Norm CC3_FMAT after CC3_L3_OMEGA1', 997 * XNORM 998 ENDIF 999C 1000C t^{B}_{1} <L3|[H^B,T^{0}_{2}],tau_{1}]|HF> 1001 CALL CC3_L3_OMEGA1(OMEGA1,ISYRES,WORK(KSMAT), 1002 * WORK(KQMAT),WORK(KTMAT),ISYMIM, 1003 * WORK(KBIOOOO),WORK(KBIOVVO), 1004 * WORK(KBIOOVV),WORK(KVVVVB),ISINTR1, 1005 * WORK(KT2TP),ISYMT2,WORK(KEND4), 1006 * LWRK4,LENSQ,WORK(KINDSQ), 1007 * ISYMB,B,ISYMD,D) 1008C 1009 IF (IPRINT .GT. 55) THEN 1010 XNORM = DDOT(NT1AM(ISYRES),OMEGA1,1,OMEGA1,1) 1011 WRITE(LUPRI,*) 'Norm CC3_FMAT after CC3_L3_OMEGA1', 1012 * XNORM 1013 ENDIF 1014C 1015C------------------------------------------- 1016C Accumulate R2 in Omega2 1017C------------------------------------------- 1018C 1019 CALL CC3_RACC(OMEGA2,WORK(KRMAT2),ISYMB,B, 1020 * ISYRES) 1021C 1022 IF (IPRINT .GT. 55) THEN 1023 XNORM = DDOT(NT2AM(ISYRES),OMEGA2,1,OMEGA2,1) 1024 WRITE(LUPRI,*) 'Norm of Rho22-after CC3_RACC',XNORM 1025 ENDIF 1026C 1027 IF (IPRINT .GT. 220) THEN 1028 CALL AROUND('After CC3_RACC: ') 1029 CALL CC_PRP(DUMMY,OMEGA2,ISYRES,0,1) 1030 ENDIF 1031C 1032 ENDDO ! B 1033 ENDDO ! ISYMB 1034C 1035C--------------------------------------- 1036C Accumulate R1 in Omega2. 1037C--------------------------------------- 1038C 1039 CALL CC3_RACC(OMEGA2,WORK(KRMAT1),ISYMD,D,ISYRES) 1040C 1041 IF (IPRINT .GT. 55) THEN 1042 XNORM = DDOT(NT2AM(ISYRES),OMEGA2,1,OMEGA2,1) 1043 WRITE(LUPRI,*) 'Norm of Rho22-after CC3_RACC-2',XNORM 1044 ENDIF 1045C 1046 IF (IPRINT .GT. 220) THEN 1047 CALL AROUND('After CC3_RACC-2: ') 1048 CALL CC_PRP(DUMMY,OMEGA2,ISYRES,0,1) 1049 ENDIF 1050C 1051 ENDDO ! D 1052 ENDDO ! ISYMD 1053C 1054COMMENT COMMENT COMMENT COMMENT COMMENT 1055COMMENT COMMENT COMMENT COMMENT COMMENT 1056COMMENT COMMENT COMMENT COMMENT COMMENT 1057c write(lupri,*) 'L3 Amplitudes in cc3_fmat : ' 1058C call dscal(nt1amx*nt1amx*nt1amx,-1.0D0,work(kt3am),1) 1059c call print_pt3(work(kt3am),isckij,3) 1060COMMENT COMMENT COMMENT COMMENT COMMENT 1061COMMENT COMMENT COMMENT COMMENT COMMENT 1062COMMENT COMMENT COMMENT COMMENT COMMENT 1063C------------------------------- 1064C Close and delete files 1065C------------------------------- 1066C 1067 CALL WCLOSE2(LUDKBC4,FNDKBC4,'DELETE') 1068 CALL WCLOSE2(LU4V,FN4V,'DELETE') 1069 CALL WCLOSE2(LU4VB,FN4VB,'DELETE') 1070 CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE') 1071 CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE') 1072 CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE') 1073 CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE') 1074 CALL WCLOSE2(LUDKBCR4,FNDKBCR4,'DELETE') 1075C 1076COMMENT COMMENT COMMENT COMMENT 1077C write(lupri,*) 'OMEGA1EFF in cc3_fmat before added' 1078C do i = 1,nt1am(isyres) 1079C write(lupri,*) i,OMEGA1EFF(i) 1080C end do 1081COMMENT COMMENT COMMENT COMMENT 1082C write(lupri,*) 'OMEGA2EFF in cc3_fmat before added' 1083C call outpak(OMEGA2EFF,NT2AM(isyres),1,lupri) 1084COMMENT COMMENT COMMENT COMMENT 1085 1086C DO I = 1,NT1AM(ISYRES) 1087C OMEGA1EFF(I) = OMEGA1EFF(I) + OMEGA1(I) 1088C END DO 1089C 1090C DO I = 1,NT2AM(ISYRES) 1091C OMEGA2EFF(I) = OMEGA2EFF(I) + OMEGA2(I) 1092C END DO 1093COMMENT COMMENT COMMENT COMMENT 1094C write(lupri,*) 'OMEGA1EFF in cc3_fmat after added' 1095C do i = 1,nt1am(isyres) 1096C write(lupri,*) i,OMEGA1EFF(i) 1097C end do 1098COMMENT COMMENT COMMENT COMMENT 1099C write(lupri,*) 'OMEGA2EFF in cc3_fmat after added' 1100C call outpak(OMEGA2EFF,NT2AM(isyres),1,lupri) 1101COMMENT COMMENT COMMENT COMMENT 1102 1103C 1104C 1105C------------------- 1106C Print timings. 1107C------------------- 1108C 1109 IF (IPRINT .GT. 9) THEN 1110 WRITE(LUPRI,*) 1111 WRITE(LUPRI,*) 1112 WRITE(LUPRI,1) 'CC3_TRAN : ',TITRAN 1113 WRITE(LUPRI,1) 'CC3_SORT : ',TISORT 1114 WRITE(LUPRI,1) 'CC3_SMAT : ',TISMAT 1115 WRITE(LUPRI,1) 'CC3_QMAT : ',TIQMAT 1116 WRITE(LUPRI,1) 'CC3_OME1 : ',TIOME1 1117 WRITE(LUPRI,*) 1118 END IF 1119C 1120C------------- 1121C End 1122C------------- 1123C 1124 CALL QEXIT('CC3_FMAT') 1125C 1126 RETURN 1127C 1128 1 FORMAT(7X,'Time used in',2X,A12,F12.2,' seconds') 1129C 1130 END 1131C /* Deck cc3_int */ 1132 SUBROUTINE CC3_INT(XOOOO,XOOVV,XOVVO,ISYINT,LU4V,FN4V, 1133 * DO_LAMB,LISTB,IDLSTB,ISYMB, 1134 * DO_LAMC,LISTC,IDLSTC,ISYMC, 1135 * WORK,LWORK) 1136C 1137C Written by Kasper Hald, Summer 2002 1138C 1139C Purpose: 1140C 1141C Calculate g_{oooo}, g_{ovvo} and g_{oovv} 1142C 1143C IF DO_LAMB = TRUE we calculate the T^{B}_{1} transformed integrals. 1144C IF DO_LAMB = TRUE and DO_LAMC = TRUE calculated T^{B}_{1} T^{C}_{1} 1145C (doubly) transformed integrals. 1146C Otherwise we calculate the ordinary (lambda^{(0)}) integrals. 1147C 1148C DO_LAMC can only be TRUE if DO_LAMB is TRUE !!! 1149C 1150C IF LVVVV = .TRUE. calculate VVVV integrals and store them on LU4V file. 1151C (LVVVV is sitting in a common block in the ccsdinp.h file and defined 1152C through the input).) 1153C 1154 IMPLICIT NONE 1155C 1156#include "priunit.h" 1157#include "dummy.h" 1158#include "maxorb.h" 1159#include "maxash.h" 1160#include "mxcent.h" 1161#include "aovec.h" 1162#include "iratdef.h" 1163#include "ccorb.h" 1164#include "ccisao.h" 1165#include "r12int.h" 1166#include "blocks.h" 1167#include "ccsdinp.h" 1168#include "ccsdsym.h" 1169#include "distcl.h" 1170#include "cbieri.h" 1171#include "eritap.h" 1172#include "second.h" 1173C 1174 INTEGER ISYINT, LU4V, IDLSTB, ISYMB, LWORK 1175 INTEGER KLAMDP, KLAMDH, KEND1, LWRK1, KEND2, LWRK2, KENDS2, LWRKS2 1176 INTEGER KENDSV, LWRKSV, KEND3, LWRK3 1177 INTEGER LUIN3O, LUIN3O2, LU3V, LU3VB 1178 INTEGER KCCFB1, KINDXB, KODCL1, KODCL2, KODBC1, KODBC2 1179 INTEGER KRDBC1, KRDBC2, KODPP1, KODPP2, KRDPP1, KRDPP2 1180 INTEGER KFREE, LFREE, NTOSYM, ISYMD1, ILLL, KRECNR 1181 INTEGER IDEL2, IDEL, ISYMD, ISYDIS, KXINT, ISYMTI, ISYMTI2 1182 INTEGER KDSRHF, K3OINT, NTOT, NUMDIS, IOPT 1183 INTEGER KLAMPB, KLAMHB, KT1AM, ISYMLP, K3OINT2 1184 INTEGER KOOOO, ISYMB2 1185 INTEGER INDEXA(MXCORB_CC) 1186 INTEGER ISYMC2,ISYMC,IDLSTC,ISYMBC,MAXBC,KLAMPC,KLAMHC,IDUM 1187C 1188#if defined (SYS_CRAY) 1189 REAL XOOOO(*), XOOVV(*), XOVVO(*), WORK(LWORK) 1190 REAL TIMTOT, TIRDAO, TIME1O, TIME3O 1191 REAL TIMEINTDEL, TIME2O2V, TIMHE1, TIMHE2 1192 REAL DTIME, ZERO, ONE 1193 real ddot ,xnormval 1194#else 1195 DOUBLE PRECISION XOOOO(*), XOOVV(*), XOVVO(*),WORK(LWORK) 1196 DOUBLE PRECISION TIMTOT, TIRDAO, TIME1O, TIME3O 1197 DOUBLE PRECISION TIMEINTDEL, TIME2O2V, TIMHE1, TIMHE2 1198 DOUBLE PRECISION DTIME, ZERO, ONE 1199 DOUBLE PRECISION ddot ,xnormval 1200#endif 1201C 1202 LOGICAL DO_LAMB,DO_LAMC 1203C 1204 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1205C 1206 CHARACTER*(*) FN4V, LISTB, LISTC 1207 CHARACTER*12 FNIN3O, FNIN3O2, FN3V, FN3VB 1208 CHARACTER*10 MODEL 1209C 1210 CALL QENTER('CC3_INT') 1211C 1212C DO_LAMC can only be TRUE if DO_LAMB is TRUE !!! 1213C 1214 IF ( (DO_LAMC) .AND. (.NOT.DO_LAMB) ) THEN 1215 WRITE(LUPRI,*)'DO_LAMC can only be TRUE if DO_LAMB is TRUE !!!' 1216 WRITE(LUPRI,*)'DO_LAMC = ', DO_LAMC 1217 WRITE(LUPRI,*)'DO_LAMB = ', DO_LAMB 1218 CALL QUIT('Incorrect logic in CC3_INT') 1219 END IF 1220 1221 TIMTOT = SECOND() 1222 TIRDAO = ZERO 1223 TIME1O = ZERO 1224 TIME3O = ZERO 1225 TIMEINTDEL = ZERO 1226 TIME2O2V = ZERO 1227 TIME2O2V = ZERO 1228 TIMHE1 = ZERO 1229 TIMHE2 = ZERO 1230C 1231 IF (.NOT. DO_LAMB) THEN 1232 ISYMB2 = ISYMOP 1233 ELSE 1234 ISYMB2 = ISYMB 1235 ENDIF 1236C 1237 IF (.NOT. DO_LAMC) THEN 1238 ISYMC2 = ISYMOP 1239 ELSE 1240 ISYMC2 = ISYMC 1241 ENDIF 1242C 1243 ISYMBC = MULD2H(ISYMB2,ISYMC2) 1244 IF (ISYMBC .NE. ISYINT) THEN 1245 WRITE(LUPRI,*)'ISYMBC and ISYINT should be the same' 1246 WRITE(LUPRI,*)'ISYMBC = ',ISYMBC 1247 WRITE(LUPRI,*)'ISYINT = ',ISYINT 1248 CALL QUIT('Symmetry mismatch in CC3_INT') 1249 END IF 1250C 1251C----------------------------- 1252C Work-space allocation 1. 1253C----------------------------- 1254C 1255 MAXBC = MAX(NT1AM(ISYMB2),NT1AM(ISYMC2)) 1256 KLAMDP = 1 1257 KLAMDH = KLAMDP + NLAMDT 1258 KLAMPB = KLAMDH + NLAMDT 1259 KLAMHB = KLAMPB + NGLMDT(ISYMB2) 1260 KLAMPC = KLAMHB + NGLMDT(ISYMB2) 1261 KLAMHC = KLAMPC + NGLMDT(ISYMC2) 1262 KT1AM = KLAMHC + NGLMDT(ISYMC2) 1263 KEND1 = KT1AM + MAX(NT1AM(ISYMOP),MAXBC) 1264 LWRK1 = LWORK - KEND1 1265C 1266 IF (LWRK1 .LT. 0) THEN 1267 CALL QUIT('Insufficient space in CC3_OVVO_INT') 1268 ENDIF 1269C 1270C------------------------------- 1271C Open scratch files 1272C------------------------------- 1273C 1274 LUIN3O = -1 1275 LU3V = -1 1276C 1277 FNIN3O = 'CC3_INT_TMP1' 1278 FN3V = 'CC3_INT_TMP2' 1279C 1280 CALL WOPEN2(LUIN3O,FNIN3O,64,0) 1281 IF (LVVVV) THEN 1282 CALL WOPEN2(LU3V,FN3V,64,0) 1283 END IF 1284C 1285 IF (DO_LAMB) THEN 1286 LUIN3O2 = -1 1287 FNIN3O2 = 'CC3_INT_TMP3' 1288 CALL WOPEN2(LUIN3O2,FNIN3O2,64,0) 1289 ENDIF 1290C 1291C-------------------------------------------- 1292C Calculate the zero'th lamda matrices. 1293C-------------------------------------------- 1294C 1295 IOPT = 1 1296 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 1297C 1298 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 1299 & LWRK1) 1300C 1301C-------------------------------------------- 1302C Calculate the "B" lamda matrices. 1303C-------------------------------------------- 1304C 1305 IF (DO_LAMB) THEN 1306 IOPT = 1 1307 CALL CC_RDRSP(LISTB,IDLSTB,ISYMB2,IOPT,MODEL, 1308 * WORK(KT1AM),DUMMY) 1309C 1310 1311 CALL CCLR_LAMTRA(WORK(KLAMDP),WORK(KLAMPB),WORK(KLAMDH), 1312 * WORK(KLAMHB),WORK(KT1AM),ISYMB2) 1313 ELSE 1314 CALL DCOPY(NLAMDT,WORK(KLAMDP),1,WORK(KLAMPB),1) 1315 CALL DCOPY(NLAMDT,WORK(KLAMDH),1,WORK(KLAMHB),1) 1316 ENDIF 1317C 1318C-------------------------------------------- 1319C Calculate the "C" lamda matrices. 1320C-------------------------------------------- 1321C 1322 IF (DO_LAMC) THEN 1323 IOPT = 1 1324 CALL CC_RDRSP(LISTC,IDLSTC,ISYMC2,IOPT,MODEL, 1325 * WORK(KT1AM),DUMMY) 1326C 1327 1328 CALL CCLR_LAMTRA(WORK(KLAMDP),WORK(KLAMPC),WORK(KLAMDH), 1329 * WORK(KLAMHC),WORK(KT1AM),ISYMC2) 1330 ELSE 1331 CALL DCOPY(NLAMDT,WORK(KLAMDP),1,WORK(KLAMPC),1) 1332 CALL DCOPY(NLAMDT,WORK(KLAMDH),1,WORK(KLAMHC),1) 1333 ENDIF 1334C 1335C------------------------------------------ 1336C Regain work space from T1-amplitudes. 1337C------------------------------------------ 1338C 1339 KEND1 = KT1AM 1340 LWRK1 = LWORK - KEND1 1341C 1342C----------------------------------- 1343C Start the loop over integrals. 1344C----------------------------------- 1345C 1346 KENDS2 = KEND1 1347 LWRKS2 = LWRK1 1348C 1349 IF (DIRECT) THEN 1350 DTIME = SECOND() 1351 IF (HERDIR) THEN 1352 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 1353 ELSE 1354 KCCFB1 = KEND1 1355 KINDXB = KCCFB1 + MXPRIM*MXCONT 1356 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 1357 LWRK1 = LWORK - KEND1 1358 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 1359 * KODPP1,KODPP2,KRDPP1,KRDPP2, 1360 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 1361 * WORK(KEND1),LWRK1,IPRERI) 1362 KEND1 = KFREE 1363 LWRK1 = LFREE 1364 ENDIF 1365 DTIME = SECOND() - DTIME 1366 TIMHE1 = TIMHE1 + DTIME 1367 NTOSYM = 1 1368 ELSE 1369 NTOSYM = NSYM 1370 ENDIF 1371C 1372 KENDSV = KEND1 1373 LWRKSV = LWRK1 1374C 1375 DO ISYMD1 = 1,NTOSYM 1376C 1377 IF (DIRECT) THEN 1378 IF (HERDIR) THEN 1379 NTOT = MAXSHL 1380 ELSE 1381 NTOT = MXCALL 1382 ENDIF 1383 ELSE 1384 NTOT = NBAS(ISYMD1) 1385 ENDIF 1386C 1387 DO ILLL = 1,NTOT 1388C 1389C----------------------------------------------------------------- 1390C If direct calculate the integrals and transposed CTR2. 1391C----------------------------------------------------------------- 1392C 1393 IF (DIRECT) THEN 1394C 1395 KEND1 = KENDSV 1396 LWRK1 = LWRKSV 1397C 1398 DTIME = SECOND() 1399 IF (HERDIR) THEN 1400 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 1401 & IPRINT) 1402 ELSE 1403 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 1404 * WORK(KODCL1),WORK(KODCL2), 1405 * WORK(KODBC1),WORK(KODBC2), 1406 * WORK(KRDBC1),WORK(KRDBC2), 1407 * WORK(KODPP1),WORK(KODPP2), 1408 * WORK(KRDPP1),WORK(KRDPP2), 1409 * WORK(KCCFB1),WORK(KINDXB), 1410 * WORK(KEND1),LWRK1,IPRERI) 1411 ENDIF 1412 DTIME = SECOND() - DTIME 1413 TIMHE2 = TIMHE2 + DTIME 1414C 1415 KRECNR = KEND1 1416 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 1417 LWRK1 = LWORK - KEND1 1418 IF (LWRK1 .LT. 0) THEN 1419 CALL QUIT('Insufficient core in CC3_INT') 1420 END IF 1421C 1422 ELSE 1423 NUMDIS = 1 1424 ENDIF 1425C 1426C----------------------------------------------------- 1427C Loop over number of distributions in disk. 1428C----------------------------------------------------- 1429C 1430 DO IDEL2 = 1,NUMDIS 1431C 1432 IF (DIRECT) THEN 1433 IDEL = INDEXA(IDEL2) 1434 IF (NOAUXB) THEN 1435 IDUM = 1 1436 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 1437 END IF 1438 ISYMD = ISAO(IDEL) 1439 ELSE 1440 IDEL = IBAS(ISYMD1) + ILLL 1441 ISYMD = ISYMD1 1442 ENDIF 1443C 1444 ISYDIS = MULD2H(ISYMD,ISYMOP) 1445C 1446C------------------------------------------ 1447C Work space allocation no. 2. 1448C------------------------------------------ 1449C 1450 KXINT = KEND1 1451 KEND2 = KXINT + NDISAO(ISYDIS) 1452 LWRK2 = LWORK - KEND2 1453C 1454 IF (LWRK2 .LT. 0) THEN 1455 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 1456 CALL QUIT('Insufficient space in CC3_INT') 1457 ENDIF 1458C 1459C-------------------------------------------- 1460C Read AO integral distribution. 1461C-------------------------------------------- 1462C 1463 DTIME = SECOND() 1464 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 1465 * WORK(KRECNR),DIRECT) 1466 DTIME = SECOND() - DTIME 1467 TIRDAO = TIRDAO + DTIME 1468C 1469C------------------------------------------ 1470C Work space allocation no. 3. 1471C------------------------------------------ 1472C 1473 ISYMTI = MULD2H(ISYMD,ISYMC2) 1474 ISYMTI2 = MULD2H(ISYMD,ISYMB2) 1475C 1476 KDSRHF = KEND2 1477 K3OINT = KDSRHF + NDSRHF(ISYMD) 1478 KEND3 = K3OINT + NMAIJK(ISYMTI) 1479 IF (DO_LAMB) THEN 1480 K3OINT2 = KEND3 1481 KEND3 = K3OINT2 + NMAIJK(ISYMTI2) 1482 ENDIF 1483 LWRK3 = LWORK - KEND3 1484C 1485 IF (LWRK3 .LT. 0) THEN 1486 WRITE(LUPRI,*) 'Need : ',KEND3,'Available : ',LWORK 1487 CALL QUIT('Insufficient space in CC3_INT') 1488 ENDIF 1489C 1490C----------------------------------------------------------------------- 1491C Transform one index in the integral batch. 1492C (alpha beta | j delta) meaning that we only need lam^0 1493C----------------------------------------------------------------------- 1494C 1495 DTIME = SECOND() 1496 ISYMLP = 1 1497 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMDP), 1498 * ISYMLP,WORK(KEND3),LWRK3,ISYDIS) 1499 DTIME = SECOND() - DTIME 1500 TIME1O = TIME1O + DTIME 1501C 1502C------------------------------------------------------------------- 1503C Calculate integral batch with three occupied indices. 1504C------------------------------------------------------------------- 1505C 1506 DTIME = SECOND() 1507 CALL CC_INT3O(WORK(K3OINT),WORK(KDSRHF),WORK(KLAMHC), 1508 * ISYMC2,WORK(KLAMDP),WORK(KEND3),LWRK3, 1509 * IDEL,ISYMD,LUIN3O,FNIN3O) 1510 DTIME = SECOND() - DTIME 1511 TIME3O = TIME3O + DTIME 1512C 1513 IF (DO_LAMB) THEN 1514 DTIME = SECOND() 1515 CALL CC_INT3O(WORK(K3OINT2),WORK(KDSRHF),WORK(KLAMHB), 1516 * ISYMB2,WORK(KLAMDP),WORK(KEND3),LWRK3, 1517 * IDEL,ISYMD,LUIN3O2,FNIN3O2) 1518c xnormval = ddot(NMAIJK(ISYMTI2),WORK(K3OINT2),1, 1519c * WORK(K3OINT2),1) 1520c write(lupri,*)'IDEL,ISYMD,luin3o2',IDEL,ISYMD 1521c write(lupri,*)'WORK(K3OINT2)',xnormval 1522 1523 DTIME = SECOND() - DTIME 1524 TIME3O = TIME3O + DTIME 1525 ENDIF 1526C 1527C-------------------------------------------------------------------- 1528C Calculate integrals for CC3. 1529C Integrals with 3 virtual indices are stored on disc 1530C-------------------------------------------------------------------- 1531C 1532 IF (LVVVV) THEN 1533 DTIME = SECOND() 1534 IF (.NOT. DO_LAMB) THEN 1535 CALL CC3_INTDEL(WORK(KXINT),1,LU3V,FN3V,WORK(KLAMDP), 1536 * 1,WORK(KLAMDH),1,1, 1537 * WORK(KEND3),LWRK3,IDEL,ISYMD) 1538 ELSE 1539 CALL CC3_INTDEL2(WORK(KXINT),ISYMOP,LU3V,FN3V, 1540 * WORK(KLAMPC),ISYMC2, 1541 * WORK(KLAMDH),ISYMOP, 1542 * WORK(KLAMPB),ISYMB2, 1543 * WORK(KLAMHB),ISYMB2, 1544 * ISYINT,WORK(KEND3),LWRK3,IDEL,ISYMD) 1545 ENDIF 1546 DTIME = SECOND() - DTIME 1547 TIMEINTDEL = TIMEINTDEL + DTIME 1548 END IF 1549C 1550 DTIME = SECOND() 1551 CALL CC3_2O2V(WORK(KXINT),ISYMOP,WORK(KDSRHF),ISYMOP, 1552 * XOVVO,XOOVV,WORK(KLAMDP),ISYMOP, 1553 * WORK(KLAMDH),ISYMOP, 1554 * WORK(KLAMPB),ISYMB2, 1555 * WORK(KLAMHC),ISYMC2,ISYINT, 1556 * WORK(KEND3),LWRK3,IDEL,ISYMD) 1557C 1558 IF (DO_LAMB) THEN 1559 CALL CC3_2O2V(WORK(KXINT),ISYMOP,WORK(KDSRHF),ISYMOP, 1560 * XOVVO,XOOVV,WORK(KLAMDP),ISYMOP, 1561 * WORK(KLAMDH),ISYMOP, 1562 * WORK(KLAMPC),ISYMC2, 1563 * WORK(KLAMHB),ISYMB2,ISYINT, 1564 * WORK(KEND3),LWRK3,IDEL,ISYMD) 1565 ENDIF 1566C 1567 DTIME = SECOND() - DTIME 1568 TIME2O2V = TIME2O2V + DTIME 1569C 1570C--------------------------------------- 1571C END ALL LOOPS 1572C--------------------------------------- 1573C 1574 ENDDO ! IDEL2 1575 ENDDO ! ILLL 1576 ENDDO ! ISYMD1 1577C 1578C------------------------ 1579C Recover work space. 1580C------------------------ 1581C 1582 KEND1 = KENDS2 1583 LWRK1 = LWRKS2 1584C 1585 CALL DZERO(XOOOO,N3ORHF(ISYINT)) 1586C 1587C for DO_LAMB lamH^B is in KLAMHB otherwise lamH^0 1588C 1589 IF (LVVVV) THEN 1590 IOPT = 3 1591 ELSE 1592 IOPT = 1 1593 END IF 1594 CALL CC3_INTSTORE(LUIN3O,FNIN3O,XOOOO,ISYINT, 1595 * WORK(KLAMHB),ISYMB2,WORK(KLAMDH),ISYMOP, 1596 * LU3V,FN3V,LU4V,FN4V,ISYINT, 1597 * WORK(KEND1),LWRK1,IOPT) 1598C 1599 IF (DO_LAMB) THEN 1600 KOOOO = KEND1 1601 KEND1 = KOOOO + N3ORHF(ISYINT) 1602 LWRK1 = LWORK - KEND1 1603C 1604 IF (LWRK1 .LT. 0) THEN 1605 CALL QUIT('Out of memory (koooo) in CC3_INT') 1606 ENDIF 1607C 1608 CALL DZERO(WORK(KOOOO),N3ORHF(ISYINT)) 1609C 1610 IOPT = 1 1611 CALL CC3_INTSTORE(LUIN3O2,FNIN3O2,WORK(KOOOO),ISYINT, 1612 * WORK(KLAMHC),ISYMC2,WORK(KLAMDH),ISYMOP, 1613 * LU3V,FN3V,LU4V,FN4V,ISYINT, 1614 * WORK(KEND1),LWRK1,IOPT) 1615C 1616 CALL DAXPY(N3ORHF(ISYINT),ONE,WORK(KOOOO),1,XOOOO,1) 1617 ENDIF 1618C 1619 CALL CC3_SORT4O(XOOOO,ISYINT,WORK(KEND1),LWRK1) 1620C 1621C------------------------------- 1622C Delete intermediate files. 1623C------------------------------- 1624C 1625 CALL WCLOSE2(LUIN3O,FNIN3O,'DELETE') 1626 IF (LVVVV) THEN 1627 CALL WCLOSE2(LU3V,FN3V,'DELETE') 1628 END IF 1629 IF (DO_LAMB) THEN 1630 CALL WCLOSE2(LUIN3O2,FNIN3O2,'DELETE') 1631 ENDIF 1632C 1633 TIMTOT = SECOND() - TIMTOT 1634C 1635C------------------------------- 1636C Write out program timings. 1637C------------------------------- 1638C 1639 IF (IPRINT .GT. 3) THEN 1640 WRITE(LUPRI,*) ' ' 1641 WRITE(LUPRI,*) 'Time used in CC3_INT :',TIMTOT 1642 WRITE(LUPRI,*) ' ' 1643 ENDIF 1644 IF (IPRINT .GT. 9) THEN 1645 WRITE(LUPRI,*) 'Time used as follows:' 1646 WRITE(LUPRI,*) ' ' 1647 WRITE(LUPRI,*) 'Time used in CCINT3O:',TIME3O 1648 ENDIF 1649C 1650 1 FORMAT(1x,A35,1X,E20.10) 1651C 1652C 1653 CALL QEXIT('CC3_INT') 1654 RETURN 1655 END 1656C /* Deck cc3_intdel2 */ 1657 SUBROUTINE CC3_INTDEL2(AOINT,ISYMAO,LUINT,FNINT, 1658 * XLAMP0,ISYMLP0,XLAMH0,ISYMLH0, 1659 * XLAMPB,ISYMLPB,XLAMHB,ISYMLHB, 1660 * ISYINT,WORK,LWORK,IDEL,ISYMD) 1661C 1662C Written by K. Hald, Summer 2002. 1663C 1664C Calculate integrals that are needed for the CC3 left hand side, 1665C and store on file. 1666C 1667C VVV,delta (V=vir.) are needed ... 1668C equals (v-bar v|v delta) + (vv|v-bar delta) 1669C 1670C 1671 IMPLICIT NONE 1672C 1673 INTEGER ISYMAO, LUINT, ISYMLP0, ISYMLH0, ISYMLPB, ISYMLHB 1674 INTEGER ISYINT, LWORK, IDEL, ISYMD 1675 INTEGER ISYABG, ISYTMP, ISYABC, KVVVV, KEND1, KEND2, LWRK1, LWRK2 1676 INTEGER ISYMG, ISYMC, ISALBE, ISYMAB, KINT, KSCR1, KSCR2 1677 INTEGER KOFF1, KOFF2, KOFF3, ISYMB, ISYMBE, ISYMAL, ISYMA 1678 INTEGER NBASAL, NBASBE, NVIRA, NAB, NBASG, IOFF 1679C 1680#if defined (SYS_CRAY) 1681 REAL AOINT(*), XLAMP0(*), XLAMH0(*) 1682 REAL XLAMPB(*), XLAMHB(*) 1683 REAL WORK(LWORK), ZERO, ONE 1684#else 1685 DOUBLE PRECISION AOINT(*), XLAMP0(*), XLAMH0(*) 1686 DOUBLE PRECISION XLAMPB(*), XLAMHB(*) 1687 DOUBLE PRECISION WORK(LWORK), ZERO, ONE 1688#endif 1689C 1690#include "priunit.h" 1691#include "ccsdsym.h" 1692#include "ccorb.h" 1693C 1694 CHARACTER FNINT*(*) 1695C 1696 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1697C 1698 CALL QENTER('CC3_INTDEL') 1699C 1700C------------------------------------------- 1701C Work space allocation. 1702C------------------------------------------- 1703C 1704 ISYABG = MULD2H(ISYMAO,ISYMD) 1705C 1706 ISYTMP = MULD2H(ISYINT,ISYMD) 1707 ISYABC = MULD2H(ISYTMP,ISYMLH0) 1708C 1709 KVVVV = 1 1710 KEND1 = KVVVV + NMAABC(ISYABC) 1711 LWRK1 = LWORK - KEND1 1712C 1713 IF (LWRK1 .LT. 0) THEN 1714 CALL QUIT('Out of memory in CC3_INTDEL') 1715 ENDIF 1716C 1717 CALL DZERO(WORK(KVVVV),NMAABC(ISYABC)) 1718C 1719C----------------------------------------------------- 1720C Transform AO-integrals to g_{v-barvv,delta} 1721C----------------------------------------------------- 1722C 1723 DO ISYMG = 1, NSYM 1724 ISYMC = MULD2H(ISYMG,ISYMLP0) 1725 ISALBE = MULD2H(ISYABG,ISYMG) 1726 ISYMAB = MULD2H(ISYABC,ISYMC) 1727 ISYTMP = MULD2H(ISYMAB,ISYMLH0) 1728C 1729 KINT = KEND1 1730 KSCR1 = KINT + NMATAB(ISYMAB)*NBAS(ISYMG) 1731 KSCR2 = KSCR1 + N2BST(ISALBE) 1732 KEND2 = KSCR2 + NEMAT1(ISYTMP) 1733 LWRK2 = LWORK - KEND2 1734COMMENT 1735COMMENT allocate to much space for kscr2 at the moment 1736COMMENT 1737C 1738 IF (LWRK2 .LT. 0) THEN 1739 CALL QUIT('Out of memory in CC3_INTDEL (2)') 1740 ENDIF 1741C 1742 DO G = 1, NBAS(ISYMG) 1743C 1744 KOFF1 = IDSAOG(ISYMG,ISYMD) + NNBST(ISALBE)*(G-1) + 1 1745 CALL CCSD_SYMSQ(AOINT(KOFF1),ISALBE,WORK(KSCR1)) 1746C 1747 DO ISYMB = 1,NSYM 1748C 1749 ISYMBE = MULD2H(ISYMB,ISYMLH0) 1750 ISYMAL = MULD2H(ISYMBE,ISALBE) 1751 ISYMA = MULD2H(ISYMAL,ISYMLPB) 1752C 1753 KOFF1 = KSCR1 1754 * + IAODIS(ISYMAL,ISYMBE) 1755 KOFF2 = IGLMVI(ISYMBE,ISYMB) + 1 1756 KOFF3 = KSCR2 1757C 1758 NBASAL = MAX(NBAS(ISYMAL),1) 1759 NBASBE = MAX(NBAS(ISYMBE),1) 1760C 1761 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMB),NBAS(ISYMBE), 1762 * ONE,WORK(KOFF1),NBASAL,XLAMH0(KOFF2),NBASBE, 1763 * ZERO,WORK(KOFF3),NBASAL) 1764C 1765 KOFF1 = IGLMVI(ISYMAL,ISYMA) + 1 1766 KOFF2 = KSCR2 1767 KOFF3 = KINT 1768 * + NMATAB(ISYMAB)*(G - 1) 1769 * + IMATAB(ISYMA,ISYMB) 1770C 1771 NBASAL = MAX(NBAS(ISYMAL),1) 1772 NVIRA = MAX(NVIR(ISYMA),1) 1773C 1774 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NBAS(ISYMAL), 1775 * ONE,XLAMPB(KOFF1),NBASAL,WORK(KOFF2),NBASAL, 1776 * ZERO,WORK(KOFF3),NVIRA) 1777C 1778 ENDDO 1779C 1780 ENDDO 1781C 1782 KOFF2 = IGLMVI(ISYMG,ISYMC) + 1 1783 KOFF3 = KVVVV 1784 * + IMAABC(ISYMAB,ISYMC) 1785C 1786 NAB = MAX(NMATAB(ISYMAB),1) 1787 NBASG = MAX(NBAS(ISYMG),1) 1788C 1789 CALL DGEMM('N','N',NMATAB(ISYMAB),NVIR(ISYMC),NBAS(ISYMG), 1790 * ONE,WORK(KINT),NAB,XLAMP0(KOFF2),NBASG, 1791 * ONE,WORK(KOFF3),NAB) 1792C 1793 ENDDO 1794C 1795C----------------------------------------------------- 1796C Transform AO-integrals to g_{vvv-bar,delta} 1797C----------------------------------------------------- 1798C 1799 DO ISYMG = 1, NSYM 1800 ISYMC = MULD2H(ISYMG,ISYMLPB) 1801 ISALBE = MULD2H(ISYABG,ISYMG) 1802 ISYMAB = MULD2H(ISYABC,ISYMC) 1803 ISYTMP = MULD2H(ISYMAB,ISYMLH0) 1804C 1805 KINT = KEND1 1806 KSCR1 = KINT + NMATAB(ISYMAB)*NBAS(ISYMG) 1807 KSCR2 = KSCR1 + N2BST(ISALBE) 1808 KEND2 = KSCR2 + NEMAT1(ISYTMP) 1809 LWRK2 = LWORK - KEND2 1810COMMENT 1811COMMENT allocate to much space for kscr2 at the moment 1812COMMENT 1813C 1814 IF (LWRK2 .LT. 0) THEN 1815 CALL QUIT('Out of memory in CC3_INTDEL (2)') 1816 ENDIF 1817C 1818 DO G = 1, NBAS(ISYMG) 1819C 1820 KOFF1 = IDSAOG(ISYMG,ISYMD) + NNBST(ISALBE)*(G-1) + 1 1821 CALL CCSD_SYMSQ(AOINT(KOFF1),ISALBE,WORK(KSCR1)) 1822C 1823 DO ISYMB = 1,NSYM 1824C 1825 ISYMBE = MULD2H(ISYMB,ISYMLH0) 1826 ISYMAL = MULD2H(ISYMBE,ISALBE) 1827 ISYMA = MULD2H(ISYMAL,ISYMLP0) 1828C 1829 KOFF1 = KSCR1 1830 * + IAODIS(ISYMAL,ISYMBE) 1831 KOFF2 = IGLMVI(ISYMBE,ISYMB) + 1 1832 KOFF3 = KSCR2 1833C 1834 NBASAL = MAX(NBAS(ISYMAL),1) 1835 NBASBE = MAX(NBAS(ISYMBE),1) 1836C 1837 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMB),NBAS(ISYMBE), 1838 * ONE,WORK(KOFF1),NBASAL,XLAMH0(KOFF2),NBASBE, 1839 * ZERO,WORK(KOFF3),NBASAL) 1840C 1841 KOFF1 = IGLMVI(ISYMAL,ISYMA) + 1 1842 KOFF2 = KSCR2 1843 KOFF3 = KINT 1844 * + NMATAB(ISYMAB)*(G - 1) 1845 * + IMATAB(ISYMA,ISYMB) 1846C 1847 NBASAL = MAX(NBAS(ISYMAL),1) 1848 NVIRA = MAX(NVIR(ISYMA),1) 1849C 1850 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NBAS(ISYMAL), 1851 * ONE,XLAMP0(KOFF1),NBASAL,WORK(KOFF2),NBASAL, 1852 * ZERO,WORK(KOFF3),NVIRA) 1853C 1854 ENDDO 1855C 1856 ENDDO 1857C 1858 KOFF2 = IGLMVI(ISYMG,ISYMC) + 1 1859 KOFF3 = KVVVV 1860 * + IMAABC(ISYMAB,ISYMC) 1861C 1862 NAB = MAX(NMATAB(ISYMAB),1) 1863 NBASG = MAX(NBAS(ISYMG),1) 1864C 1865 CALL DGEMM('N','N',NMATAB(ISYMAB),NVIR(ISYMC),NBAS(ISYMG), 1866 * ONE,WORK(KINT),NAB,XLAMPB(KOFF2),NBASG, 1867 * ONE,WORK(KOFF3),NAB) 1868C 1869 ENDDO 1870C 1871C---------------------------------------- 1872C Save the g_{vvv,delta} to disc. 1873C---------------------------------------- 1874C 1875 IF (NMAABC(ISYABC) .GT. 0) THEN 1876 KOFF1 = IDEL - IBAS(ISYMD) 1877 IOFF = I3VDEL(ISYABC,ISYMD) + NMAABC(ISYABC)*(KOFF1-1) + 1 1878 CALL PUTWA2(LUINT,FNINT,WORK(KVVVV),IOFF,NMAABC(ISYABC)) 1879 ENDIF 1880C 1881C-------------------------- 1882C End. 1883C-------------------------- 1884C 1885 CALL QEXIT('CC3_INTDEL') 1886C 1887 RETURN 1888 END 1889C /* Deck cc3_intdel2 */ 1890 SUBROUTINE CC3_BARINT(B1AM,ISYMB,XLAMDP,XLAMDH,WORK,LWORK, 1891 * LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR) 1892C 1893C Written by Kasper Hald, Summer 2002. 1894C 1895C Interface to calculate the (ck|de)-bar and (ck|lm)-bar integrals 1896C 1897C NSIMTR = 1 !!!!!!! 1898C 1899 IMPLICIT NONE 1900C 1901#include "priunit.h" 1902#include "dummy.h" 1903#include "maxorb.h" 1904#include "maxash.h" 1905#include "mxcent.h" 1906#include "aovec.h" 1907#include "iratdef.h" 1908#include "ccorb.h" 1909#include "ccisao.h" 1910#include "r12int.h" 1911#include "blocks.h" 1912#include "ccsdinp.h" 1913#include "ccsdsym.h" 1914#include "distcl.h" 1915#include "cbieri.h" 1916#include "eritap.h" 1917C#include "cclr.h" 1918C 1919 INTEGER ISYMB, LWORK, LU3SRTR, LUCKJDR 1920 INTEGER KEND1, LWRK1, KEND2, LWRK2, KENDSV, LWRKSV 1921 INTEGER KCCFB1, KINDXB, KODCL1, KODCL2, KODBC1, KODBC2 1922 INTEGER KRDBC1, KRDBC2, KODPP1, KODPP2, KRDPP1, KRDPP2 1923 INTEGER KFREE, LFREE, NTOSYM, ISYMD1, NTOT, ILLL 1924 INTEGER KRECNR, NUMDIS, IDEL2, IDEL, ISYMD, ISYDIS 1925 INTEGER ISYAIK, KXINT, KLAMPB, KLAMHB, ISYCKJ, KCKJ 1926 INTEGER INDEXA(MXCORB_CC),IDUM 1927C 1928#if defined (SYS_CRAY) 1929#else 1930 DOUBLE PRECISION B1AM(*), XLAMDP(*), XLAMDH(*), WORK(LWORK) 1931#endif 1932C 1933 CHARACTER*(*) FN3SRTR, FNCKJDR 1934 CHARACTER*1 CDUMMY 1935C 1936 CALL QENTER('CC3_BARINT') 1937C 1938 CDUMMY = ' ' 1939C 1940 KLAMPB = 1 1941 KLAMHB = KLAMPB + NLAMDT 1942 KEND1 = KLAMHB + NLAMDT 1943 LWRK1 = LWORK - KEND1 1944C 1945C-------------------------------------- 1946C Calculate lambda-B 1947C-------------------------------------- 1948C 1949 CALL CCLR_LAMTRA(XLAMDP,WORK(KLAMPB),XLAMDH, 1950 * WORK(KLAMHB),B1AM,ISYMB) 1951C 1952C-------------------------------- 1953C Calculate integrals 1954C-------------------------------- 1955C 1956 IF (DIRECT) THEN 1957 IF (HERDIR) THEN 1958 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 1959 ELSE 1960 KCCFB1 = KEND1 1961 KINDXB = KCCFB1 + MXPRIM*MXCONT 1962 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 1963 LWRK1 = LWORK - KEND1 1964 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 1965 & KODPP1,KODPP2,KRDPP1,KRDPP2, 1966 & KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 1967 & WORK(KEND1),LWRK1,IPRERI) 1968 KEND1 = KFREE 1969 LWRK1 = LFREE 1970 ENDIF 1971 NTOSYM = 1 1972 ELSE 1973 NTOSYM = NSYM 1974 ENDIF 1975C 1976 KENDSV = KEND1 1977 LWRKSV = LWRK1 1978C 1979 DO ISYMD1 = 1,NTOSYM 1980C 1981 IF (DIRECT) THEN 1982 IF (HERDIR) THEN 1983 NTOT = MAXSHL 1984 ELSE 1985 NTOT = MXCALL 1986 ENDIF 1987 ELSE 1988 NTOT = NBAS(ISYMD1) 1989 ENDIF 1990C 1991 DO ILLL = 1,NTOT 1992C 1993C------------------------------------------ 1994C If direct calculate the integrals. 1995C------------------------------------------ 1996C 1997 IF (DIRECT) THEN 1998C 1999 KEND1 = KENDSV 2000 LWRK1 = LWRKSV 2001C 2002 IF (HERDIR) THEN 2003 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 2004 & IPRERI) 2005 ELSE 2006 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 2007 & WORK(KODCL1),WORK(KODCL2), 2008 & WORK(KODBC1),WORK(KODBC2), 2009 & WORK(KRDBC1),WORK(KRDBC2), 2010 & WORK(KODPP1),WORK(KODPP2), 2011 & WORK(KRDPP1),WORK(KRDPP2), 2012 & WORK(KCCFB1),WORK(KINDXB), 2013 & WORK(KEND1), LWRK1,IPRERI) 2014 ENDIF 2015C 2016 KRECNR = KEND1 2017 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 2018 LWRK1 = LWORK - KEND1 2019 IF (LWRK1 .LT. 0) THEN 2020 CALL QUIT('Insufficient core in CC3_BARINT') 2021 END IF 2022C 2023 ELSE 2024 NUMDIS = 1 2025 ENDIF 2026C 2027C-------------------------------------------------- 2028C Loop over number of distributions in disk. 2029C-------------------------------------------------- 2030C 2031 DO IDEL2 = 1,NUMDIS 2032C 2033 IF (DIRECT) THEN 2034 IDEL = INDEXA(IDEL2) 2035 IF (NOAUXB) THEN 2036 IDUM = 1 2037 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 2038 END IF 2039 ISYMD = ISAO(IDEL) 2040 ELSE 2041 IDEL = IBAS(ISYMD1) + ILLL 2042 ISYMD = ISYMD1 2043 ENDIF 2044C 2045 ISYDIS = MULD2H(ISYMD,ISYMOP) 2046 ISYAIK = MULD2H(ISYDIS,ISYMOP) 2047 ISYCKJ = ISYDIS 2048C 2049C------------------------------------------ 2050C Work space allocation no. 2. 2051C------------------------------------------ 2052C 2053 KXINT = KEND1 2054 KCKJ = KXINT + NDISAO(ISYDIS) 2055 KEND2 = KCKJ + NCKI(ISYCKJ) 2056 LWRK2 = LWORK - KEND2 2057C 2058 IF (LWRK2 .LT. 0) THEN 2059 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 2060 CALL QUIT('Insufficient space in CC3_FMAT ') 2061 ENDIF 2062C 2063 CALL DZERO(WORK(KCKJ),NCKI(ISYCKJ)) 2064C 2065C----------------------------------------- 2066C Read in batch of integrals. 2067C----------------------------------------- 2068C 2069 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 2070 * WORK(KRECNR),DIRECT) 2071C 2072C----------------------------------------------------------- 2073C Calculate transformed integrals used in t3am. 2074C----------------------------------------------------------- 2075C 2076 CALL CC3_T3INT(WORK(KXINT),XLAMDP,XLAMDH, 2077 * B1AM,ISYMB,WORK(KEND2),LWRK2, 2078 * IDEL,ISYMD,2,LU3SRTR,FN3SRTR, 2079 * LUCKJDR,FNCKJDR) 2080C 2081 ENDDO ! IDEL2 2082 ENDDO ! ILLL 2083 ENDDO ! ISYMD1 2084C 2085 CALL QEXIT('CC3_BARINT') 2086 RETURN 2087 END 2088C /* Deck wbd_ground */ 2089 SUBROUTINE WBD_GROUND(T2TP,ISYMT2,TMAT,TRVIR,TRVIR2,TROCC,ISYINT, 2090 * WMAT,WORK,LWORK,INDSQ,LENSQ, 2091 * ISYMB,B,ISYMC,C) 2092C 2093C Kasper Hald, Summer 2002 2094C 2095C WBD(ai,k,j) = WBD(ai,k,j) 2096C + t(ai,ej) (Dk|Be) + t(ai,ej) (Bk|De) 2097C - t(ai,Bl) (Dk|lj) - t(ai,Dl) (Bk|lj) 2098C 2099C ISYINT is the symmetry of the integrals. 2100C ISYMT2 is the symmetry of T2TP. 2101C 2102 IMPLICIT NONE 2103C 2104#include "priunit.h" 2105#include "ccorb.h" 2106#include "ccsdinp.h" 2107#include "ccsdsym.h" 2108C#include "cclr.h" 2109C 2110 INTEGER ISYMT2, ISYINT, LWORK, LENSQ, ISYMB, ISYMC 2111 INTEGER INDSQ(LENSQ,6) 2112 INTEGER ISYMBC, ISYRES, JSAIKJ, ISYMDK, LENGTH, ISYMK 2113 INTEGER ISYMD, ISYAIJ, KOFF1, KOFF2, KOFF3, NTOAIJ, NVIRD 2114 INTEGER ISYAIL, ISYLKJ, ISYMJ, ISYMLK, ISYML, ISYMAI 2115 INTEGER ISYAIK, NTOTAI, NRHFL 2116C 2117 integer isymdkb 2118C 2119#if defined (SYS_CRAY) 2120 REAL T2TP(*), TRVIR(*), TRVIR2(*), TROCC(*) 2121 REAL TMAT(*), WMAT(*), WORK(LWORK) 2122 REAL ZERO, ONE, DDOT, XNORM 2123#else 2124 DOUBLE PRECISION T2TP(*), TRVIR(*), TRVIR2(*), TROCC(*) 2125 DOUBLE PRECISION TMAT(*), WMAT(*), WORK(LWORK) 2126 DOUBLE PRECISION ZERO, ONE, DDOT, XNORM 2127#endif 2128C 2129 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2130C 2131 CALL QENTER('WBD_GROUND') 2132 2133C 2134C-------------------------------- 2135C First virtual contribution. 2136C-------------------------------- 2137C 2138 ISYMBC = MULD2H(ISYMB,ISYMC) 2139 ISYRES = MULD2H(ISYINT,ISYMT2) 2140 JSAIKJ = MULD2H(ISYMBC,ISYRES) 2141 ISYMDK = MULD2H(ISYMBC,ISYINT) 2142C 2143 LENGTH = NCKIJ(JSAIKJ) 2144C 2145 IF (LWORK .LT. LENGTH) THEN 2146 CALL QUIT('Insufficient core in WBD_GROUND (1)') 2147 ENDIF 2148C 2149 DO ISYMK = 1,NSYM 2150C 2151 ISYMD = MULD2H(ISYMK,ISYMDK) 2152 ISYAIJ = MULD2H(ISYMK,JSAIKJ) 2153C 2154 KOFF1 = IT2SP(ISYAIJ,ISYMD) + 1 2155 KOFF2 = ICKATR(ISYMDK,ISYMB) + NT1AM(ISYMDK)*(B - 1) 2156 * + IT1AM(ISYMD,ISYMK) + 1 2157 KOFF3 = ISAIKJ(ISYAIJ,ISYMK) + 1 2158C 2159 NTOAIJ = MAX(1,NCKI(ISYAIJ)) 2160 NVIRD = MAX(NVIR(ISYMD),1) 2161C 2162 CALL DGEMM('N','N',NCKI(ISYAIJ),NRHF(ISYMK), 2163 * NVIR(ISYMD),ONE,T2TP(KOFF1),NTOAIJ, 2164 * TRVIR(KOFF2),NVIRD,ZERO, 2165 * WORK(KOFF3),NTOAIJ) 2166C 2167 ENDDO 2168C 2169 2170 DO I = 1,LENGTH 2171 WMAT(I) = WMAT(I) + WORK(INDSQ(I,3)) 2172 ENDDO 2173 2174C 2175 IF (IPRINT .GT. 55) THEN 2176 XNORM = DDOT(NCKIJ(JSAIKJ),WMAT,1,WMAT,1) 2177 WRITE(LUPRI,*) 'In WBD_GROUND: 1. Norm of WMAT ',XNORM 2178 ENDIF 2179C 2180C--------------------------------- 2181C Second virtual contribution. 2182C--------------------------------- 2183C 2184 ISYMBC = MULD2H(ISYMB,ISYMC) 2185 ISYRES = MULD2H(ISYINT,ISYMT2) 2186 JSAIKJ = MULD2H(ISYMBC,ISYRES) 2187 ISYMDK = MULD2H(ISYMBC,ISYINT) 2188C 2189 LENGTH = NCKIJ(JSAIKJ) 2190C 2191 IF (LWORK .LT. LENGTH) THEN 2192 CALL QUIT('Insufficient core in WBD_GROUND (2)') 2193 ENDIF 2194C 2195 DO ISYMK = 1,NSYM 2196C 2197 ISYMD = MULD2H(ISYMK,ISYMDK) 2198 ISYAIJ = MULD2H(ISYMK,JSAIKJ) 2199C 2200 KOFF1 = IT2SP(ISYAIJ,ISYMD) + 1 2201 KOFF2 = ICKATR(ISYMDK,ISYMC) + NT1AM(ISYMDK)*(C - 1) 2202 * + IT1AM(ISYMD,ISYMK) + 1 2203 KOFF3 = ISAIKJ(ISYAIJ,ISYMK) + 1 2204C 2205 NTOAIJ = MAX(1,NCKI(ISYAIJ)) 2206 NVIRD = MAX(NVIR(ISYMD),1) 2207C 2208 CALL DGEMM('N','N',NCKI(ISYAIJ),NRHF(ISYMK), 2209 * NVIR(ISYMD),ONE,T2TP(KOFF1),NTOAIJ, 2210 * TRVIR2(KOFF2),NVIRD,ZERO, 2211 * WORK(KOFF3),NTOAIJ) 2212C 2213 ENDDO 2214C 2215 DO I = 1,LENGTH 2216 WMAT(I) = WMAT(I) + WORK(I) 2217 ENDDO 2218C 2219 IF (IPRINT .GT. 55) THEN 2220 XNORM = DDOT(NCKIJ(JSAIKJ),WMAT,1,WMAT,1) 2221 WRITE(LUPRI,*) 'In WBD_GROUND: 2. Norm of WMAT ',XNORM 2222 ENDIF 2223C 2224C--------------------------------- 2225C First occupied contribution. 2226C--------------------------------- 2227C 2228 ISYAIL = MULD2H(ISYMB,ISYMT2) 2229 ISYLKJ = MULD2H(ISYMC,ISYINT) 2230C 2231 DO ISYMJ = 1,NSYM 2232C 2233 ISYMLK = MULD2H(ISYMJ,ISYLKJ) 2234C 2235 DO J = 1,NRHF(ISYMJ) 2236C 2237 DO ISYMK = 1,NSYM 2238C 2239 ISYML = MULD2H(ISYMK,ISYMLK) 2240 ISYMAI = MULD2H(ISYAIL,ISYML) 2241 ISYAIK = MULD2H(ISYMAI,ISYMK) 2242C 2243 KOFF1 = IT2SP(ISYAIL,ISYMB) 2244 * + NCKI(ISYAIL)*(B - 1) 2245 * + ICKI(ISYMAI,ISYML) + 1 2246 KOFF2 = ISJIKA(ISYLKJ,ISYMC) 2247 * + NMAJIK(ISYLKJ)*(C - 1) 2248 * + ISJIK(ISYMLK,ISYMJ) 2249 * + NMATIJ(ISYMLK)*(J - 1) 2250 * + IMATIJ(ISYML,ISYMK) + 1 2251 KOFF3 = ISAIKJ(ISYAIK,ISYMJ) 2252 * + NCKI(ISYAIK)*(J - 1) 2253 * + ICKI(ISYMAI,ISYMK) + 1 2254C 2255 NTOTAI = MAX(1,NT1AM(ISYMAI)) 2256 NRHFL = MAX(1,NRHF(ISYML)) 2257C 2258 CALL DGEMM('N','N',NT1AM(ISYMAI),NRHF(ISYMK), 2259 * NRHF(ISYML),-ONE,T2TP(KOFF1),NTOTAI, 2260 * TROCC(KOFF2),NRHFL,ONE,WMAT(KOFF3), 2261 * NTOTAI) 2262C 2263 ENDDO 2264 ENDDO 2265 ENDDO 2266C 2267 IF (IPRINT .GT. 55) THEN 2268 XNORM = DDOT(NCKIJ(JSAIKJ),WMAT,1,WMAT,1) 2269 WRITE(LUPRI,*) 'In WBD_GROUND: 3. Norm of WMAT ',XNORM 2270 ENDIF 2271C 2272C---------------------------------- 2273C Second occupied contribution. 2274C---------------------------------- 2275C 2276 ISYAIL = MULD2H(ISYMC,ISYMT2) 2277 ISYLKJ = MULD2H(ISYMB,ISYINT) 2278C 2279 DO ISYMJ = 1,NSYM 2280C 2281 ISYMLK = MULD2H(ISYMJ,ISYLKJ) 2282C 2283 DO J = 1,NRHF(ISYMJ) 2284C 2285 DO ISYMK = 1,NSYM 2286C 2287 ISYML = MULD2H(ISYMK,ISYMLK) 2288 ISYMAI = MULD2H(ISYAIL,ISYML) 2289 ISYAIK = MULD2H(ISYMAI,ISYMK) 2290C 2291 KOFF1 = IT2SP(ISYAIL,ISYMC) 2292 * + NCKI(ISYAIL)*(C - 1) 2293 * + ICKI(ISYMAI,ISYML) + 1 2294 KOFF2 = ISJIKA(ISYLKJ,ISYMB) 2295 * + NMAJIK(ISYLKJ)*(B - 1) 2296 * + ISJIK(ISYMLK,ISYMJ) 2297 * + NMATIJ(ISYMLK)*(J - 1) 2298 * + IMATIJ(ISYML,ISYMK) + 1 2299 KOFF3 = ISAIKJ(ISYAIK,ISYMJ) 2300 * + NCKI(ISYAIK)*(J - 1) 2301 * + ICKI(ISYMAI,ISYMK) + 1 2302C 2303 NTOTAI = MAX(1,NT1AM(ISYMAI)) 2304 NRHFL = MAX(1,NRHF(ISYML)) 2305C 2306 CALL DGEMM('N','N',NT1AM(ISYMAI),NRHF(ISYMK), 2307 * NRHF(ISYML),-ONE,T2TP(KOFF1),NTOTAI, 2308 * TROCC(KOFF2),NRHFL,ZERO,TMAT(KOFF3), 2309 * NTOTAI) 2310C 2311 ENDDO 2312 ENDDO 2313 ENDDO 2314C 2315 DO I = 1,LENGTH 2316 WMAT(I) = WMAT(I) + TMAT(INDSQ(I,3)) 2317 ENDDO 2318C 2319 IF (IPRINT .GT. 55) THEN 2320 XNORM = DDOT(NCKIJ(JSAIKJ),WMAT,1,WMAT,1) 2321 WRITE(LUPRI,*) 'In WBD_GROUND: 4. Norm of WMAT ',XNORM 2322 ENDIF 2323C 2324C--------------------------------------- 2325C End. 2326C--------------------------------------- 2327C 2328 CALL QEXIT('WBD_GROUND') 2329C 2330 RETURN 2331 END 2332C /* Deck cclr_trvir */ 2333 SUBROUTINE CCLR_TRVIR(XINT,TRINT,XLAMDP,ISYMLP,ISYMD,D,ISYINT, 2334 * WORK,LWORK) 2335C 2336C Henrik Koch and Alfredo Sanchez. Dec 1994 2337C 2338C Transform (ck|alpha d) integrals to (ck|a d) 2339C 2340C Ove Christiansen 11-1-1996: ISYINT is symmetry of (ck|alpha d) 2341C 2342C Kasper Hald, 01062002 : General symmetry of Lambda 2343C 2344#include "implicit.h" 2345C 2346 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 2347 DIMENSION XINT(*),TRINT(*),XLAMDP(*),WORK(LWORK) 2348C 2349#include "priunit.h" 2350#include "ccorb.h" 2351#include "ccsdsym.h" 2352C 2353 CALL QENTER('CCLR_TRVIR') 2354C 2355 ISYDIS = MULD2H(ISYMD,ISYINT) 2356C 2357 DO 100 ISYMA = 1,NSYM 2358C 2359COMMENT COMMENT 2360COMMENT COMMENT 2361 ISYMAL = MULD2H(ISYMA,ISYMLP) 2362C ISYMCK = MULD2H(ISYMA,ISYDIS) 2363 ISYMCK = MULD2H(ISYMAL,ISYDIS) 2364COMMENT COMMENT 2365COMMENT COMMENT 2366C 2367 NTOTCK = MAX(NT1AM(ISYMCK),1) 2368COMMENT COMMENT 2369C NBASA = MAX(NBAS(ISYMA),1) 2370 NBASA = MAX(NBAS(ISYMAL),1) 2371COMMENT COMMENT 2372C 2373COMMENT COMMENT COMMENT 2374C KOFF1 = ICKALP(ISYMCK,ISYMA) + 1 2375 KOFF1 = ICKALP(ISYMCK,ISYMAL) + 1 2376C KOFF2 = ILMVIR(ISYMA) + 1 2377 KOFF2 = IGLMVI(ISYMAL,ISYMA) 2378 * + 1 2379COMMENT COMMENT COMMENT 2380 KOFF3 = ICKATR(ISYMCK,ISYMA) + 1 2381C 2382COMMENT COMMENT 2383C CALL DGEMM('N','N',NT1AM(ISYMCK),NVIR(ISYMA),NBAS(ISYMA), 2384 CALL DGEMM('N','N',NT1AM(ISYMCK),NVIR(ISYMA),NBAS(ISYMAL), 2385COMMENT COMMENT 2386 * ONE,XINT(KOFF1),NTOTCK,XLAMDP(KOFF2),NBASA, 2387 * ZERO,TRINT(KOFF3),NTOTCK) 2388C 2389 100 CONTINUE 2390C 2391 CALL QEXIT('CCLR_TRVIR') 2392C 2393 RETURN 2394 END 2395C /* Deck cclr_trocc */ 2396 SUBROUTINE CCLR_TROCC(XINT,TRINT,XLAMDH,ISYMLH,WORK,LWORK) 2397C 2398C Henrik Koch and Alfredo Sanchez. Dec 1994 2399C 2400C Transform (ia|j delta) integrals to (ia|j k) and sort as (ij,k,a) 2401C 2402C Kasper Hald, Summer 2002 -> generalisation to nonsymmetric lambda 2403C 2404 IMPLICIT NONE 2405C 2406#include "priunit.h" 2407#include "ccorb.h" 2408#include "ccsdsym.h" 2409C 2410 INTEGER ISYMLH, LWORK 2411 INTEGER ISYMK, ISYINT, ISYAIJ, LENGTH, ISYMD 2412 INTEGER NTOIAJ, NBASD, KOFF1, KOFF2, ISYMJ, ISYMAI 2413 INTEGER ISYMI, ISYMA, ISYMJI, ISYJIK, NTOJIK 2414C 2415#if defined (SYS_CRAY) 2416#else 2417 DOUBLE PRECISION XINT(*), TRINT(*), XLAMDH(*), WORK(LWORK) 2418 DOUBLE PRECISION ZERO, ONE, DDOT 2419#endif 2420C 2421 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 2422C 2423C 2424 CALL QENTER('CCLR_TROCC') 2425C 2426 DO ISYMD = 1,NSYM 2427C 2428COMMENT COMMENT 2429C ISYMK = ISYMD 2430C ISYAIJ = MULD2H(ISYMD,ISYMOP) 2431 ISYMK = MULD2H(ISYMD,ISYMLH) 2432C ISYINT = MULD2H(ISYMLH,ISYMOP) 2433C ISYAIJ = MULD2H(ISYMK,ISYINT) 2434 ISYAIJ = MULD2H(ISYMD,ISYMOP) 2435COMMENT COMMENT 2436C 2437 LENGTH = NCKI(ISYAIJ) * NRHF(ISYMK) 2438 IF (LENGTH .GT. LWORK) THEN 2439 CALL QUIT('Not enough space in CCLR_TROCC') 2440 END IF 2441C 2442 NTOIAJ = MAX(NCKI(ISYAIJ),1) 2443 NBASD = MAX(NBAS(ISYMD),1) 2444C 2445 KOFF1 = ICKID(ISYAIJ,ISYMD) + 1 2446COMMENT COMMENT 2447COMMENT COMMENT 2448C KOFF2 = ILMRHF(ISYMD) + 1 2449 KOFF2 = IGLMRH(ISYMD,ISYMK) + 1 2450COMMENT COMMENT 2451COMMENT COMMENT 2452C 2453 CALL DGEMM('N','N',NCKI(ISYAIJ),NRHF(ISYMK),NBAS(ISYMD), 2454 * ONE,XINT(KOFF1),NTOIAJ,XLAMDH(KOFF2),NBASD, 2455 * ZERO,WORK,NTOIAJ) 2456C 2457C Sort 2458C 2459 DO ISYMJ = 1,NSYM 2460C 2461 ISYMAI = MULD2H(ISYAIJ,ISYMJ) 2462C 2463 DO ISYMI = 1,NSYM 2464C 2465 ISYMA = MULD2H(ISYMAI,ISYMI) 2466 ISYMJI = MULD2H(ISYMJ,ISYMI) 2467 ISYJIK = MULD2H(ISYMJI,ISYMK) 2468C 2469 DO K = 1,NRHF(ISYMK) 2470C 2471 DO J = 1,NRHF(ISYMJ) 2472C 2473 DO I = 1,NRHF(ISYMI) 2474C 2475 NTOJIK = NMAJIK(ISYJIK) 2476C 2477 KOFF1 = NCKI(ISYAIJ)*(K - 1) 2478 * + ISAIK(ISYMAI,ISYMJ) 2479 * + NT1AM(ISYMAI)*(J - 1) 2480 * + IT1AM(ISYMA,ISYMI) 2481 * + NVIR(ISYMA)*(I - 1) + 1 2482 KOFF2 = ISJIKA(ISYJIK,ISYMA) 2483 * + ISJIK(ISYMJI,ISYMK) 2484 * + NMATIJ(ISYMJI)*(K - 1) 2485 * + IMATIJ(ISYMJ,ISYMI) 2486 * + NRHF(ISYMJ)*(I - 1) + J 2487C 2488 CALL DCOPY(NVIR(ISYMA),WORK(KOFF1),1, 2489 * TRINT(KOFF2),NTOJIK) 2490C 2491 ENDDO ! I 2492 ENDDO ! J 2493 ENDDO ! K 2494 ENDDO ! ISYMI 2495 ENDDO ! ISYMJ 2496C 2497 ENDDO ! ISYMD 2498C 2499 CALL QEXIT('CCLR_TROCC') 2500C 2501 RETURN 2502 END 2503