1! 2!... Copyright (c) 2015 by the authors of Dalton (see below). 3!... All Rights Reserved. 4!... 5!... The source code in this file is part of 6!... "Dalton, a molecular electronic structure program, 7!... Release DALTON2016 (2015), see http://daltonprogram.org" 8!... 9!... This source code is provided under a written licence and may be 10!... used, copied, transmitted, or stored only in accord with that 11!... written licence. 12!... 13!... In particular, no part of the source code or compiled modules may 14!... be distributed outside the research group of the licence holder. 15!... This means also that persons (e.g. post-docs) leaving the research 16!... group of the licence holder may not take any part of Dalton, 17!... including modified files, with him/her, unless that person has 18!... obtained his/her own licence. 19!... 20!... For further information, including how to get a licence, see: 21!... http://daltonprogram.org 22! 23! 24C 25C /* Deck cc_lhtr */ 26 SUBROUTINE CC_LHTR3(ECURR, 27 * FRHO1,LUFR1,FRHO2,LUFR2, 28 * FC1AM,LUFC1,FC2AM,LUFC2, 29 * RHO1,RHO2,CTR1,CTR2,WORK,LWORK, 30 * NSIMTR,IVEC,ITR,LRHO1) 31C 32C Written by Rasmus Faber, summer 2017 33C 34C Based on CC_LHTR by Asger Halkier & Henrik Koch. 35C 36C Purpose: 37C 38C Calculate left hand side transformation of a triplet trial vector 39C in an AO-integral direct fashion. 40C 41C (IVEC is first number for first vector on files, FRHO1,FC1AM,FC2AM, 42C FRHO12,FC12AM; 43C ITR is first vector on FRHO2) 44C 45C Current models are: CCSD..? 46C 47C 48C 49#include "implicit.h" 50#include "priunit.h" 51#include "dummy.h" 52#include "maxorb.h" 53#include "maxash.h" 54#include "mxcent.h" 55#include "aovec.h" 56#include "iratdef.h" 57#include "ccorb.h" 58C 59C #include "infind.h" replaced by: #include <ccisao.h> 60C C. Haettig 13.08.2004 61C 62#include "ccisao.h" 63#include "blocks.h" 64#include "ccsdinp.h" 65#include "ccsections.h" 66#include "ccfield.h" 67#include "ccsdsym.h" 68#include "ccsdio.h" 69#include "ccinftap.h" 70#include "distcl.h" 71#include "cbieri.h" 72#include "cclr.h" 73#include "eritap.h" 74#include "ccslvinf.h" 75#include "qm3.h" 76!#include "qmmm.h" 77#include "ccnoddy.h" 78#include "r12int.h" 79#include "ccr12int.h" 80C 81 CHARACTER(LEN=*), PARAMETER :: myname = 'CC_LHTR3' 82C 83 CHARACTER*10 MODEL 84 CHARACTER*8 FRHO2, FRHO1, FC2AM, FC1AM, FN3SRT, FNTOC, FRHO12, 85 & FC12AM 86 CHARACTER*7 FN3FOP2X 87 CHARACTER*6 CFIL, DFIL, PQFIL, FN3V, FNCKJD, FN3VI, FN3FOPX 88 CHARACTER*5 FNDKBC3 89 CHARACTER*4 FN4V, FNDKBC 90 CHARACTER*3 APROXR12 91 CHARACTER*1 LR 92C 93 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0, 94 & FOUR = 4.0D0) 95 PARAMETER (FOURTH = 0.25D0, THREE = 3.0D0) 96 PARAMETER (ONEM = -1.0D0, TWOM = -TWO) 97 DIMENSION INDEXA(MXCORB_CC) 98 DIMENSION RHO1(LRHO1,NSIMTR), CTR1(LRHO1,NSIMTR) 99 DIMENSION RHO2(*), CTR2(*) 100 DIMENSION WORK(LWORK) 101 DIMENSION IADRPQ(MXCORB_CC,MAXSIM) 102C 103 LOGICAL ETRAN,FCKCON,CCSDR12,SKIPMI 104C 105 CALL QENTER(myname) 106C 107 THIRD = ONE/THREE 108C 109C 110 !set CCSDR12 111 CCSDR12 = .FALSE. 112C 113 SKIPMI = .FALSE. 114C 115C--------------------------------- 116C Initialze timing parameters. 117C--------------------------------- 118C 119 TIMTOT = ZERO 120 TIMTOT = SECOND() 121 TIMHE1 = ZERO 122 TIMHE2 = ZERO 123 TIRDAO = ZERO 124 TIMIO = ZERO 125 TIMEYI = ZERO 126 TIMEXI = ZERO 127 TIMEMI = ZERO 128 TIMETI = ZERO 129 TIMEZ1 = ZERO 130 TIMEZ2 = ZERO 131 TIMFCK = ZERO 132 TIMEC = ZERO 133 TIMEA = ZERO 134 TIMEH = ZERO 135 TIMEI = ZERO 136 TIME3O = ZERO 137 TIME1O = ZERO 138 TIMETP = ZERO 139 TIMEBF = ZERO 140 TIM2BF = ZERO 141 TIMEEM = ZERO 142 TIMEB = ZERO 143 TIMEG = ZERO 144 TIMEC2 = ZERO 145 TIM12B = ZERO 146 TIM12A = ZERO 147 TIM11B = ZERO 148 TIMEAM = ZERO 149 TIM2EM = ZERO 150 TIM22E = ZERO 151 TIM22A = ZERO 152 TIM22C = ZERO 153 TIM22D = ZERO 154 TIM2AO = ZERO 155 TIM2MO = ZERO 156C 157C----------------------------- 158C Work-space allocation 1. 159C----------------------------- 160C 161 ISYRES = MULD2H(ISYMTR,ISYMOP) 162C 163 KLAMDP = 1 164 KLAMDH = KLAMDP + NLAMDT 165 KC1T2 = KLAMDH + NLAMDT 166 KCT1AO = KC1T2 + N2BST(ISYRES)*NSIMTR 167 KDENSI = KCT1AO + NGLMRH(ISYMTR)*NSIMTR 168 KFOCK = KDENSI + N2BST(1) 169 KFOCKG = KFOCK + N2BST(ISYMOP) 170 IF (CC2.AND.(.NOT.NONHF)) THEN 171 KEND0 = KFOCKG + N2BST(ISYMTR)*NSIMTR 172 ELSE 173 KYMAT = KFOCKG + N2BST(ISYMTR)*NSIMTR 174 KYDEN = KYMAT + NMATAB(ISYRES)*NSIMTR 175 KXMAT = KYDEN + N2BST(ISYRES)*NSIMTR 176 KEND0 = KXMAT + NMATIJ(ISYRES)*NSIMTR 177 ENDIF 178C 179 IF (CCR12) THEN 180 KRHOR12 = KEND0 181 KEND0 = KRHOR12 + NTR12SQ(ISYMTR) 182 END IF 183C 184 IF (CC3) THEN 185 KOVVO = KEND0 186 KOOVV = KOVVO + NT2SQ(1) 187 KEND0 = KOOVV + NT2SQ(1) 188 ENDIF 189C 190 IF (CC2) THEN 191 IF (NONHF) THEN 192 KFCKHF = KEND0 193 KEND0 = KFCKHF + N2BST(1) 194 END IF 195 KT2AM = KEND0 196 KT1AM = KT2AM + MAX(NT2AMX,NT2AM(ISYMOP)) 197 KEND1 = KT1AM + MAX(NT1AMX,NT1AM(ISYMOP)) 198 ELSE 199 KT2AM = KEND0 200 KT1AM = KT2AM + MAX(NT2AMX,NT2AM(ISYMOP),NT2AM(ISYRES)) 201 KEND1 = KT1AM + MAX(NT1AMX,NT1AM(ISYMOP)) 202 ENDIF 203 LWRK1 = LWORK - KEND1 204C 205 IF (LWRK1 .LT. 0) THEN 206 CALL QUIT('Insufficient space in '//myname) 207 ENDIF 208C 209C------------------------------- 210C Open scratch file CCINT3O. 211C------------------------------- 212C 213 LUIN30 = -1 214 CALL WOPEN2(LUIN30,'CCINT3O',64,0) 215C 216C------------------------------------------------------------- 217C Open files with C- & D-intermediates needed in 22-block. 218C------------------------------------------------------------- 219C 220 LUCIM = -1 221 LUDIM = -1 222C 223 CFIL = 'PMAT_C' 224 DFIL = 'PMAT_D' 225C 226 CALL WOPEN2(LUCIM,CFIL,64,0) 227 CALL WOPEN2(LUDIM,DFIL,64,0) 228C 229C------------------------------------------------------------- 230C Open files with P- & Q-intermediates needed in 21-block. 231C------------------------------------------------------------- 232C 233 IF (.TRUE.) THEN 234 235 LUPQ = -1 236 PQFIL = 'CCPQIM' 237C 238 CALL WOPEN2(LUPQ,PQFIL,64,0) 239 240 END IF 241C 242C------------------------------------------- 243C Read zero'th order cluster amplitudes. 244C------------------------------------------- 245C 246 DTIME = SECOND() 247 248 IOPT = 3 249 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM)) 250 251 DTIME = SECOND() - DTIME 252 TIMIO = TIMIO + DTIME 253C 254C---------------------------------------------------------------------- 255C Initialize result-arrays and zero out single-vectors in CCD-calc. 256C---------------------------------------------------------------------- 257C 258 NRHO2 = MAX(NT2SQ(ISYRES),NT2ORT(ISYRES)+NT2ORT3(ISYRES)) 259 IF (CC2) NRHO2 = NT2SQ(ISYRES) 260C 261 CALL DZERO(RHO1(1,1),NT1AM(ISYRES)*NSIMTR) 262 CALL DZERO(RHO2,NRHO2) 263C 264 IF (CCD) THEN 265 CALL DZERO(WORK(KT1AM),NT1AMX) 266 CALL DZERO(CTR1(1,1),NT1AM(ISYMTR)*NSIMTR) 267 ENDIF 268C 269 CALL DZERO(WORK(KFOCKG),N2BST(ISYMTR)*NSIMTR) 270C 271C---------------------------------- 272C Calculate the lamda matrices. 273C---------------------------------- 274C 275 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 276 & LWRK1) 277C 278C------------------------------------------ 279C Regain work space from T1-amplitudes. 280C------------------------------------------ 281C 282 KEND1 = KT1AM 283 LWRK1 = LWORK - KEND1 284C 285C---------------------------------- 286C Calculate the density matrix. 287C---------------------------------- 288C 289 ISYMH = 1 290 IC = 1 291 CALL CC_AODENS(WORK(KLAMDP),WORK(KLAMDH),WORK(KDENSI),ISYMH, 292 & IC,WORK(KEND1),LWRK1) 293C 294C-------------------------------------------------------- 295C initialize start adress for P & Q intermediates: 296C-------------------------------------------------------- 297C 298 IADRSTRT = 1 299C 300C------------------------------------------------------------ 301C Loop over trial vectors for constructing intermediates. 302C------------------------------------------------------------ 303C 304 DO 95 IV = 1,NSIMTR 305C 306C-------------------------------------------------------- 307C Calculate contraction of CTR1 with lamda-matrix. 308C-------------------------------------------------------- 309C 310 KOFFAO = KCT1AO + NGLMRH(ISYMTR)*(IV - 1) 311C 312 CALL CCLT_CT1AO(CTR1(1,IV),WORK(KLAMDP),WORK(KOFFAO)) 313C 314C----------------------------------------------- 315C Calculate contraction of CTR1 and T2AM. 316C----------------------------------------------- 317C 318 KOFFCT = KC1T2 + N2BST(ISYRES)*(IV - 1) 319C 320 IOPT = -1 321 CALL CC_C1T2C(WORK(KOFFCT),CTR1(1,IV),WORK(KT2AM), 322 * WORK(KLAMDP),WORK(KLAMDH),WORK(KLAMDP), 323 * WORK(KLAMDH),WORK(KEND1),LWRK1,ISYMTR, 324 * ISYMOP,IOPT) 325C 326C-------------------------------------- 327C Read CTR2+ from disc into RHO2. 328C-------------------------------------- 329C 330 DTIME = SECOND() 331 CALL CC_RVEC3(LUFC2,FC2AM,2*NT2AM(ISYMTR),NT2AM(ISYMTR), 332 * IVEC+IV-1,0,RHO2) 333 DTIME = SECOND() - DTIME 334 TIMIO = TIMIO + DTIME 335C 336C----------------------- 337C Square up CTR2. 338C----------------------- 339C 340 CALL CCRHS3_R2IJ(RHO2,WORK(KEND1),LWRK1,ISYMTR) 341 CALL CC_T2SQ(RHO2,CTR2,ISYMTR) 342C 343C-------------------------------------- 344C Read CTR2- from disc into RHO2. 345C-------------------------------------- 346C 347 DTIME = SECOND() 348 CALL CC_RVEC3(LUFC2,FC2AM,2*NT2AM(ISYMTR),NT2AM(ISYMTR), 349 * IVEC+IV-1,NT2AM(ISYMTR),RHO2) 350 DTIME = SECOND() - DTIME 351 TIMIO = TIMIO + DTIME 352C 353 CALL CC_T2SQ3A(RHO2,CTR2,ISYMTR,ONEM) 354C 355C 356 IF ( DEBUG ) THEN 357 RHO1N = DDOT(NT1AM(ISYMTR),CTR1(1,IV),1,CTR1(1,IV),1) 358 RHO2N = DDOT(NT2SQ(ISYMTR),CTR2,1,CTR2,1) 359 WRITE(LUPRI,1) 'Norm of CTR1 -Read in before loop: ',RHO1N 360 WRITE(LUPRI,1) 'Norm of CTR2 -Read in before loop: ',RHO2N 361 ENDIF 362C 363 IF (.NOT. (CC2 .AND.(.NOT.NONHF))) THEN 364C 365C-------------------------------------------------- 366C Calculate Y-intermediate of CTR2 and T2AM. 367C-------------------------------------------------- 368C 369 KOFFYI = KYMAT + NMATAB(ISYRES)*(IV - 1) 370 KOFFYD = KYDEN + N2BST(ISYRES)*(IV - 1) ! Not needed 371C 372 TIMEYI = SECOND() 373 CALL CC_YI(WORK(KOFFYI),CTR2,ISYMTR,WORK(KT2AM),ISYMOP, 374 * WORK(KEND1),LWRK1) 375 CALL CC_YD(WORK(KOFFYD),WORK(KOFFYI),ISYRES,WORK(KLAMDH), 376 * WORK(KLAMDP),ISYMOP,WORK(KEND1),LWRK1) 377 TIMEYI = SECOND() - TIMEYI 378 379 CALL DAXPY(N2BST(ISYRES),ONE, 380 & WORK(KOFFYD),1,WORK(KOFFCT),1) 381C 382C----------------------------------------------------- 383C Calculate X-intermediate of CTR2 and T2AM. 384C----------------------------------------------------- 385C 386 KOFFXI = KXMAT + NMATIJ(ISYRES)*(IV - 1) 387C 388 TIMEXI = SECOND() 389 CALL CC_XI(WORK(KOFFXI),CTR2,ISYMTR,WORK(KT2AM),ISYMOP, 390 * WORK(KEND1),LWRK1) 391 CALL CC_XD3(WORK(KOFFCT),WORK(KOFFXI),ISYRES,WORK(KLAMDH), 392 & WORK(KLAMDP),1,WORK(KEND1),LWRK1) 393 TIMEXI = SECOND() - TIMEXI 394 END IF 395C 396C------------------------------------------------------------ 397C Precalculate P and Q intermediates 398C and restore CTR2 which is overwritten in CC_PQI: 399C------------------------------------------------------------ 400C 401 IF (.NOT. CC2) THEN 402C 403 IF ( DEBUG ) THEN 404 RHO1N = DDOT(NT1AM(ISYMTR),CTR1(1,IV),1,CTR1(1,IV),1) 405 RHO2N = DDOT(NT2SQ(ISYMTR),CTR2,1,CTR2,1) 406 WRITE(LUPRI,1) 'Norm of CTR1 - before CC_PQI: ',RHO1N 407 WRITE(LUPRI,1) 'Norm of CTR2 - before CC_PQI: ',RHO2N 408 RHO2N = DDOT(NT2AM(ISYMOP),WORK(KT2AM),1,WORK(KT2AM),1) 409 WRITE(LUPRI,1) 'Norm of T2AM - before CC_PQI: ',RHO2N 410 ENDIF 411C 412C Calculate the triplet P-matrix 413C 414 CALL CC_PQI3(CTR2,ISYMTR,WORK(KT2AM),ISYMOP,PQFIL,LUPQ, 415 * IADRPQ,IADRSTRT,IV,WORK(KLAMDH),0, 416 & WORK(KEND1),LWRK1) 417C 418C We need to rearrange CTR2 in order to calculate triplet Q matrix: 419C First get -(+)L - (-)L from current (+)L - (-)L by transposing, 420C then scaling 421 422 CALL CC_T2SQTRANSP(CTR2,ISYMTR) 423C 424 CALL DSCAL(NT2SQ(ISYMTR),ONEM,CTR2,1) 425C 426C RHO2 still contains (-)L, reorder it so (-)L(ai<bj) to -(-)L(aj<bi) 427C 428 CALL CCSD_TCMEPK3(RHO2,ISYMTR) 429C 430C Create L'(ai,bj) = -(+)L(ai,bj)-(-)L(ai,bj) + 2(-)L(aj,bi) 431 CALL CC_T2SQ3A(RHO2,CTR2,ISYMTR,TWO) 432C 433 CALL CC_PQI3(CTR2,ISYMTR,WORK(KT2AM),ISYMOP,PQFIL,LUPQ, 434 * IADRPQ,IADRSTRT,IV,WORK(KLAMDH),1, 435 & WORK(KEND1),LWRK1) 436C 437C Create L'(ai,bj) = -(+)L(ai,bj)-(-)L(ai,bj) - 2(-)L(aj,bi) 438C 439 CALL CC_T2SQ3A(RHO2,CTR2,ISYMTR,-FOUR) 440C 441 CALL CC_PQI3(CTR2,ISYMTR,WORK(KT2AM),ISYMOP,PQFIL,LUPQ, 442 * IADRPQ,IADRSTRT,IV,WORK(KLAMDH),2, 443 & WORK(KEND1),LWRK1) 444 445C 446 IF (MINSCR) THEN 447C 448C We better clean up the mess.... 449C Restore (+)L + (-)L on CTR2 450 CALL CC_T2SQ3A(RHO2,CTR2,ISYMTR,TWO) 451 CALL DSCAL(NT2SQ(ISYMTR),ONEM,CTR2,1) 452 453 END IF 454C DTIME = SECOND() 455C CALL CC_RVEC(LUFC2,FC2AM,NT2AM(ISYMTR),NT2AM(ISYMTR), 456C * IVEC+IV-1,RHO2) 457C DTIME = SECOND() - DTIME 458C TIMIO = TIMIO + DTIME 459C 460C CALL CC_T2SQ(RHO2,CTR2,ISYMTR) 461C 462 IF ( DEBUG ) THEN 463 RHO1N = DDOT(NT1AM(ISYMTR),CTR1(1,IV),1,CTR1(1,IV),1) 464 RHO2N = DDOT(NT2SQ(ISYMTR),CTR2,1,CTR2,1) 465 WRITE(LUPRI,1) 'Norm of CTR1 - before CC_PQI: ',RHO1N 466 WRITE(LUPRI,1) 'Norm of CTR2 - before CC_PQI: ',RHO2N 467 RHO2N = DDOT(NT2AM(ISYMOP),WORK(KT2AM),1,WORK(KT2AM),1) 468 WRITE(LUPRI,1) 'Norm of T2AM - before CC_PQI: ',RHO2N 469 ENDIF 470 471 ENDIF 472C 473 95 CONTINUE 474C 475C 476C--------------------------------------------------------------- 477C Calculate transposed CTR2 if t2tcor and reinitialize RHO2. 478C--------------------------------------------------------------- 479C 480 IF (MINSCR) THEN 481C 482C In the integral loop, we need the + combinations of (+) and (-) 483C We currently have the - combination, the (-) part is in RHO2 484C 485C CALL CC_T2SQ3A(RHO2,CTR2,ISYMTR,TWO) 486 CALL DZERO(RHO2,NRHO2) 487 ENDIF 488C 489C------------------------------------------------------------- 490C Regain work space from T2-amplitudes in CC2-calculation. 491C------------------------------------------------------------- 492C 493 IF (CC2) THEN 494C 495 KEND1 = KT2AM 496 LWRK1 = LWORK - KEND1 497C 498 ENDIF 499C 500C----------------------------------- 501C Start the loop over integrals. 502C----------------------------------- 503C 504 KENDS2 = KEND1 505 LWRKS2 = LWRK1 506C 507 IF (DIRECT) THEN 508 DTIME = SECOND() 509 IF (HERDIR) THEN 510 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 511 ELSE 512 KCCFB1 = KEND1 513 KINDXB = KCCFB1 + MXPRIM*MXCONT 514 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 515 LWRK1 = LWORK - KEND1 516 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 517 * KODPP1,KODPP2,KRDPP1,KRDPP2, 518 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 519 * WORK(KEND1),LWRK1,IPRERI) 520 KEND1 = KFREE 521 LWRK1 = LFREE 522 ENDIF 523 DTIME = SECOND() - DTIME 524 TIMHE1 = TIMHE1 + DTIME 525 NTOSYM = 1 526 ELSE 527 NTOSYM = NSYM 528 ENDIF 529C 530 KENDSV = KEND1 531 LWRKSV = LWRK1 532C 533 ICDEL1 = 0 534 DO 100 ISYMD1 = 1,NTOSYM 535C 536 IF (DIRECT) THEN 537 IF (HERDIR) THEN 538 NTOT = MAXSHL 539 ELSE 540 NTOT = MXCALL 541 ENDIF 542 ELSE 543 NTOT = NBAS(ISYMD1) 544 ENDIF 545C 546 DO 110 ILLL = 1,NTOT 547C 548C----------------------------------------------------------------- 549C If direct calculate the integrals and transposed CTR2. 550C----------------------------------------------------------------- 551C 552 IF (DIRECT) THEN 553C 554 KEND1 = KENDSV 555 LWRK1 = LWRKSV 556C 557 DTIME = SECOND() 558 IF (HERDIR) THEN 559 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 560 & IPRINT) 561 ELSE 562 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 563 * WORK(KODCL1),WORK(KODCL2), 564 * WORK(KODBC1),WORK(KODBC2), 565 * WORK(KRDBC1),WORK(KRDBC2), 566 * WORK(KODPP1),WORK(KODPP2), 567 * WORK(KRDPP1),WORK(KRDPP2), 568 * WORK(KCCFB1),WORK(KINDXB), 569 * WORK(KEND1),LWRK1,IPRERI) 570 ENDIF 571 DTIME = SECOND() - DTIME 572 TIMHE2 = TIMHE2 + DTIME 573C 574 KRECNR = KEND1 575 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 576 LWRK1 = LWORK - KEND1 577 IF (LWRK1 .LT. 0) THEN 578 CALL QUIT('Insufficient core in CC_LHTR') 579 END IF 580C 581 IF ((MINSCR .AND. T2TCOR) .AND. (.NOT. CC2)) THEN 582C 583 KCTR2T = KEND1 584 KEND1 = KCTR2T + NT2SQ(ISYMTR) 585 LWRK1 = LWORK - KEND1 586 IF (LWRK1 .LT. 0) THEN 587 CALL QUIT('Insufficient core in CC_LHTR') 588 END IF 589C 590 JSYM = ISYMTR 591 CALL DCOPY(NT2SQ(ISYMTR),CTR2,1,WORK(KCTR2T),1) 592 AUTIME = SECOND() 593 CALL CCSD_T2TP(WORK(KCTR2T),WORK(KEND1),LWRK1,JSYM) 594 AUTIME = SECOND() - AUTIME 595 TIMETP = TIMETP + AUTIME 596C 597 END IF 598C 599 ELSE 600 NUMDIS = 1 601 ENDIF 602C 603C---------------------------------------------------- 604C Loop over number of trial vectors nsimtr. 605C---------------------------------------------------- 606C 607 DO 115 IV = 1, NSIMTR 608C 609 IF (.NOT. MINSCR) THEN 610 611 CALL QUIT( 612 & 'Triplet LHTR without MINSCR not implemented' 613 & ) 614C 615C------------------------------------------------------------- 616C Read CTR2 from disc into RHO2 and square up. 617C------------------------------------------------------------- 618C 619 DTIME = SECOND() 620 CALL CC_RVEC(LUFC2,FC2AM,NT2AM(ISYMTR),NT2AM(ISYMTR), 621 * IVEC+IV-1,RHO2) 622 DTIME = SECOND() - DTIME 623 TIMIO = TIMIO + DTIME 624C 625 CALL CC_T2SQ(RHO2,CTR2,ISYMTR) 626C 627C---------------------------------------------- 628C Read result vector from disc. 629C---------------------------------------------- 630C 631 DTIME = SECOND() 632 CALL CC_RVEC(LUFR2,FRHO2,NRHO2,NRHO2, 633 * IV+ITR-1,RHO2) 634 DTIME = SECOND() - DTIME 635 TIMIO = TIMIO + DTIME 636C 637C----------------------------------------------------- 638C Calculate transposed CTR2 if t2tcor. 639C----------------------------------------------------- 640C 641 IF (T2TCOR) THEN 642C 643 KCTR2T = KEND1 644 KEND1 = KCTR2T + NT2SQ(ISYMTR) 645 LWRK1 = LWORK - KEND1 646 IF (LWRK1 .LT. 0) THEN 647 CALL QUIT('Insufficient core in CC_LHTR') 648 END IF 649C 650 JSYM = ISYMTR 651 CALL DCOPY(NT2SQ(ISYMTR),CTR2,1,WORK(KCTR2T),1) 652 AUTIME = SECOND() 653 CALL CCSD_T2TP(WORK(KCTR2T),WORK(KEND1), 654 * LWRK1,JSYM) 655 AUTIME = SECOND() - AUTIME 656 TIMETP = TIMETP + AUTIME 657C 658 ENDIF 659 ENDIF 660C 661C-------------------------------------------------------------- 662C Calculate adresses for CTR-dependent intermediates. 663C-------------------------------------------------------------- 664C 665 KOFFAO = KCT1AO + NGLMRH(ISYMTR)*(IV - 1) 666 KOFFCT = KC1T2 + N2BST(ISYRES)*(IV - 1) 667 KOFFYI = KYMAT + NMATAB(ISYRES)*(IV - 1) 668 KOFFYD = KYDEN + N2BST(ISYRES)*(IV - 1) 669 KOFFFG = KFOCKG + N2BST(ISYMTR)*(IV - 1) 670C 671C----------------------------------------------------- 672C Loop over number of distributions in disk. 673C----------------------------------------------------- 674C 675 DO 120 IDEL2 = 1,NUMDIS 676C 677 IF (DIRECT) THEN 678 IDEL = INDEXA(IDEL2) 679C IF (NOAUXB) THEN 680c IDUM = 1 681C CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 682C END IF 683 ISYMD = ISAO(IDEL) 684 ELSE 685 IDEL = IBAS(ISYMD1) + ILLL 686 ISYMD = ISYMD1 687 ENDIF 688C 689 ISYDIS = MULD2H(ISYMD,ISYMOP) 690C 691C------------------------------------------ 692C Work space allocation no. 2. 693C------------------------------------------ 694C 695 KXINT = KEND1 696 KEND2 = KXINT + NDISAO(ISYDIS) 697 LWRK2 = LWORK - KEND2 698C 699 IF (LWRK2 .LT. 0) THEN 700 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 701 CALL QUIT('Insufficient space in CC_LHTR') 702 ENDIF 703C 704C-------------------------------------------- 705C Read AO integral distribution. 706C-------------------------------------------- 707C 708 AUTIME = SECOND() 709 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 710 * WORK(KRECNR),DIRECT) 711 AUTIME = SECOND() - AUTIME 712 TIRDAO = TIRDAO + AUTIME 713C 714C------------------------------------------- 715C Calculate the AO-fock-matrix. 716C Keep Exchange part only 717C------------------------------------------- 718C 719 AUTIME = SECOND() 720 ISYDEN = ISYMTR 721 CALL CC_AOFOCK3(WORK(KXINT),WORK(KOFFCT),WORK(KOFFFG), 722 * WORK(KEND2),LWRK2,IDEL,ISYMD,ISYDEN) 723 AUTIME = SECOND() - AUTIME 724 TIMFCK = TIMFCK + AUTIME 725C 726C------------------------------------------ 727C Work space allocation no. 3. 728C------------------------------------------ 729C 730 ISYMTI = MULD2H(ISYMD,ISYMOP) 731 ISY21I = MULD2H(ISYMD,ISYRES) 732C 733 KDSRHF = KEND2 734 K3OINT = KDSRHF + NDSRHF(ISYMD) 735 IF (CC2) THEN 736 KEND3 = K3OINT + NMAIJK(ISYMTI) 737 ELSE 738 KSCRTI = K3OINT + NMAIJK(ISYMTI) ! is KSCRTI ever used? 739 KSCRZI = KSCRTI + NT2BCD(ISYMTI) 740 KSCRWI = KSCRZI + NT2BCD(ISY21I) 741 KSCRVI = KSCRWI + NT2BCD(ISY21I) 742 KEND3 = KSCRVI + NT2BCD(ISY21I) 743 ENDIF 744 LWRK3 = LWORK - KEND3 745C 746 IF (LWRK3 .LT. 0) THEN 747 WRITE(LUPRI,*) 'Need : ',KEND3,'Available : ',LWORK 748 CALL QUIT('Insufficient space in CC_LHTR') 749 ENDIF 750C 751C-------------------------------------------------------- 752C Transform one index in the integral batch. 753C-------------------------------------------------------- 754C 755 AUTIME = SECOND() 756 ISYMLP = 1 757 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMDP), 758 * ISYMLP,WORK(KEND3),LWRK3,ISYDIS) 759 AUTIME = SECOND() - AUTIME 760 TIME1O = TIME1O + AUTIME 761C 762C------------------------------------------------------------------- 763C Calculate integral batch with three occupied indices. 764C------------------------------------------------------------------- 765C 766 AUTIME = SECOND() 767 CALL CC_INT3O(WORK(K3OINT),WORK(KDSRHF),WORK(KLAMDH), 768 * ISYMOP,WORK(KLAMDP),WORK(KEND3),LWRK3, 769 * IDEL,ISYMD,LUIN30,'CCINT3O') 770 AUTIME = SECOND() - AUTIME 771 TIME3O = TIME3O + AUTIME 772C 773C-------------------------------------------------------------- 774C Calculate intermediates needed for the 21-block. 775C-------------------------------------------------------------- 776C 777 IF (.NOT. CC2) THEN 778C 779 IADR = IADRPQ(IDEL,IV) 780 CALL GETWA2(LUPQ,PQFIL,WORK(KSCRZI), 781 * IADR,NT2BCD(ISY21I)) 782C 783 IADR = IADRPQ(IDEL,IV) + NT2BCD(ISY21I) 784 CALL GETWA2(LUPQ,PQFIL,WORK(KSCRVI), 785 * IADR,NT2BCD(ISY21I)) 786C 787 IADR = IADRPQ(IDEL,IV) + 2*NT2BCD(ISY21I) 788 CALL GETWA2(LUPQ,PQFIL,WORK(KSCRWI), 789 * IADR,NT2BCD(ISY21I)) 790C 791 IF ( DEBUG ) THEN 792 PN=DDOT(NT2BCD(ISY21I),WORK(KSCRZI),1,WORK(KSCRZI),1) 793 QN=DDOT(NT2BCD(ISY21I),WORK(KSCRVI),1,WORK(KSCRVI),1) 794 WRITE(LUPRI,*) 'IV,IDEL,IADR,P,Q',IV,IDEL,IADR,PN,QN 795 ENDIF 796 797C 798C-------------------------------------------------- 799C Calculate the LT21I contribution. 800C-------------------------------------------------- 801C 802 803 IOPT = 1 804 ISYLRD = MULD2H(ISYMD,ISYRES) 805 AUTIME = SECOND() 806 CALL CC_21I2(RHO1(1,IV),WORK(KXINT),ISYDIS,DUMMY,0, 807 * WORK(KSCRZI),WORK(KSCRVI),ISYLRD, 808 * DUMMY, DUMMY, 0, 809 * WORK(KLAMDP),WORK(KLAMDH),ISYMOP, 810 * WORK(KLAMDP),ISYMOP, 811 * WORK(KEND3),LWRK3,IOPT,.FALSE.,.FALSE., 812 * .FALSE.) 813 AUTIME = SECOND() - AUTIME 814 TIMEI = TIMEI + AUTIME 815 816 IF ( DEBUG ) THEN 817 RHO1N=DDOT(NT1AM(ISYMTR),RHO1(1,IV),1,RHO1(1,IV),1) 818 WRITE(LUPRI,*) 'Norm of RHO1 after CC_21I2:',IV, 819 & RHO1N,IDEL 820 ENDIF 821 822C 823C-------------------------------------------------- 824C Calculate the LT21H contribution. 825C------------------------------------------------- 826C 827 CALL DZERO(WORK(KSCRVI),NT2BCD(ISY21I)) 828 AUTIME = SECOND() 829 CALL CC_21H(RHO1(1,IV),ISYRES,WORK(KSCRWI), 830 * WORK(KSCRVI),WORK(KSCRZI),ISYLRD, 831 * WORK(K3OINT),ISYMOP,WORK(KEND3), 832 * LWRK3,ISYMD) 833 AUTIME = SECOND() - AUTIME 834 TIMEH = TIMEH + AUTIME 835 836 IF ( DEBUG ) THEN 837 RHO1N=DDOT(NT1AM(ISYMTR),RHO1(1,IV),1,RHO1(1,IV),1) 838 WRITE(LUPRI,*) 'Norm of RHO1 after CC_21H:',IV,RHO1N 839 ENDIF 840 841 ENDIF! (.NOT. CC) 842C 843C------------------------------------------ 844C Work space allocation no. 4. 845C------------------------------------------ 846C 847 ISSCRM = MULD2H(ISYMD,ISYMTR) 848C 849 KSCRM = KEND3 ! KSCR{T,V,W,Z} is not uused beyond this! 850 ! Reset to KSCRTI? 851 KSCRM2 = KSCRM + NT2BCD(ISSCRM) 852 KEND4 = KSCRM2 + NT2BCD(ISSCRM) 853 LWRK4 = LWORK - KEND4 854C 855 IF (LWRK4 .LT. 0) THEN 856 WRITE(LUPRI,*) 'Need : ',KEND4,'Available : ',LWORK 857 CALL QUIT('Insufficient space in'//myname) 858 ENDIF 859C 860C-------------------------------------------------------------- 861C Construct the partially transformed CTR2-vector. 862C-------------------------------------------------------------- 863C 864 AUTIME = SECOND() 865 ICON = 1 866 ISYMLP = 1 867 CALL CC_T2AO(CTR2,WORK(KLAMDP),ISYMLP,WORK(KSCRM), 868 * WORK(KEND4),LWRK4,IDEL,ISYMD, 869 * ISYMTR,ICON) 870 871 CALL CC_T2AO3(CTR2,WORK(KLAMDP),ISYMLP,WORK(KSCRM2), 872 * WORK(KEND4),LWRK4,IDEL,ISYMD, 873 * ISYMTR) 874C 875 AUTIME = SECOND() - AUTIME 876 TIM2AO = TIM2AO + AUTIME 877C 878C------------------------------------------------------- 879C Calculate the LT21C- and D contributions. 880C------------------------------------------------------- 881C 882 AUTIME = SECOND() 883 IOPT = 1 884 IF ( CC2 ) THEN 885 ICON = 2 886 ELSE 887 ICON = 1 888 ENDIF 889C 890 ISYMM = MULD2H(ISYMD,ISYMTR) 891 CALL CC_21DC(RHO1(1,IV),CTR2,ISYMTR,WORK(KSCRM),ISYMM, 892 * WORK(KSCRM),ISYMM,WORK(KXINT), 893 * WORK(KLAMDH),ISYMOP,WORK(KLAMDP),ISYMOP, 894 * WORK(KLAMDH),ISYMOP,WORK(KLAMDP),ISYMOP, 895 * WORK(KEND4),LWRK4,IDEL,ISYMD,IOPT,ICON) 896 AUTIME = SECOND() - AUTIME 897 TIMEC = TIMEC + AUTIME 898 899 IF ( DEBUG ) THEN 900 RHO1N=DDOT(NT1AM(ISYMTR),RHO1(1,IV),1,RHO1(1,IV),1) 901 WRITE(LUPRI,*) 'Norm of RHO1 after CC_21DC:', 902 & IV,RHO1N 903 ENDIF 904 905C 906C-------------------------------------------------------- 907C Calculate the LT12B & LT22B contributions. 908C-------------------------------------------------------- 909C 910 IF (.NOT. CC2) THEN 911C 912 AUTIME = SECOND() 913 IOPT = 2 914 ISYMM = MULD2H(ISYMD,ISYMTR) 915 CALL CC_BF3(WORK(KXINT),RHO2,WORK(KLAMDP),1, 916 * WORK(KOFFAO),ISYMTR,WORK(KLAMDP),1, 917 * WORK(KSCRM),ISYMM,WORK(KSCRM2),ISYMM, 918 * WORK(KEND4),LWRK4,IDEL,ISYMD,IOPT) 919 AUTIME = SECOND() - AUTIME 920 TIM2BF = TIM2BF + AUTIME 921C 922 ENDIF 923C 924 120 CONTINUE 925C 926 IF (.NOT. MINSCR) THEN 927C 928C----------------------------------------- 929C Write out result vector. 930C----------------------------------------- 931C 932 DTIME = SECOND() 933 CALL CC_WVEC(LUFR2,FRHO2,NRHO2,NRHO2,IV + ITR -1, 934 * RHO2) 935 DTIME = SECOND() - DTIME 936 TIMIO = TIMIO + DTIME 937C 938 ENDIF 939C 940 115 CONTINUE 941 110 CONTINUE 942 100 CONTINUE 943C 944C------------------------ 945C Recover work space. 946C------------------------ 947C 948 KEND1 = KENDS2 949 LWRK1 = LWRKS2 950C 951C----------------------------- 952C Loop over trial vectors. 953C----------------------------- 954C 955 DO 125 IV = 1,NSIMTR 956C 957C 958 IF (.NOT. MINSCR) THEN 959C 960C Read in CTR2+ and CTR. and 961C construct CTR+ - CTR- intermediate 962C T2 vectors are overwritten further down 963C and needs to be read in once more if 964C not first iteration of the loop. 965C Alternatively we could calculate and 966C write out the M intermediate before the 967C loop over integrals. 968 ELSE 969C 970C--------------------------------------- 971C We have in memory CTR2+ + CTR2- 972C in the following we need 973C CTR2+ - CTR- 974C--------------------------------------- 975C 976 CALL CC_T2SQTRANSP(CTR2,ISYMTR) 977 978 ENDIF 979C 980 IF ( DEBUG ) THEN 981 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 982 RHO2N = DDOT(NRHO2,RHO2,1,RHO2,1) 983 WRITE(LUPRI,1) 'Norm of RHO1 loop over vect. 1. ', RHO1N 984 WRITE(LUPRI,1) 'Norm of RHO2 loop over vect. 1. ', RHO2N 985 ENDIF 986C 987C----------------------------------------------------------- 988C Calculate adresses for CTR-dependent intermediates. 989C Skip large section for CC2. 990C----------------------------------------------------------- 991C 992 KOFFFG = KFOCKG + N2BST(ISYMTR)*(IV - 1) 993 KOFFYI = KYMAT + NMATAB(ISYRES)*(IV - 1) 994 KOFFXI = KXMAT + NMATIJ(ISYRES)*(IV - 1) 995C 996 IF (.NOT. CC2) THEN 997C 998C------------------------------ 999C Work space allocation. 1000C------------------------------ 1001C 1002 IF (.NOT.SKIPMI) THEN 1003 KMINT = KT2AM + MAX(NT2SQ(ISYRES),NT2AM(1)) 1004 KEND1 = KMINT + N3ORHF(ISYRES) 1005 LWRK1 = LWORK - KEND1 1006C 1007 IF (LWRK1 .LT. 0) THEN 1008 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1009 CALL QUIT('Insufficient memory for M-intermediate '// 1010 & 'in CC_LHTR') 1011 ENDIF 1012C 1013C----------------------------------------------------------- 1014C Calculate transposed M-intermediate of CTR2 and T2AM. 1015C-------------------------------------------------- 1016C 1017 TIMEMI = SECOND() 1018 CALL CC_MI(WORK(KMINT),CTR2,ISYMTR,WORK(KT2AM),ISYMOP, 1019 * WORK(KEND1),LWRK1) 1020 TIMEMI = SECOND() - TIMEMI 1021C 1022 END IF 1023C 1024 IF ( DEBUG ) THEN 1025 WRITE(LUPRI,1) 'Norm of CTR2 : ', 1026 & DDOT(NT2SQ(ISYMTR),CTR2,1,CTR2,1) 1027 WRITE(LUPRI,1) 'Norm of T2AM : ', 1028 & DDOT(NT2AMX,WORK(KT2AM),1,WORK(KT2AM),1) 1029 WRITE(LUPRI,1) 'Norm of M-INT : ', 1030 & DDOT(N3ORHF(ISYRES),WORK(KMINT),1,WORK(KMINT),1) 1031 ENDIF 1032C 1033C----------------------------------------- 1034C Calculate the LT21G contribution. 1035C----------------------------------------- 1036C 1037 TIMEG = SECOND() 1038 CALL CC_21G(RHO1(1,IV),WORK(KMINT),ISYRES,WORK(KLAMDH), 1039 * WORK(KEND1),LWRK1,ISYMOP,LUIN30,'CCINT3O') 1040 TIMEG = SECOND() - TIMEG 1041C 1042 CALL CC_T2SQTRANSP(CTR2,ISYMTR) 1043C 1044C---------------------------------------------- 1045C Transform the RHO2 vector to MO basis. 1046C---------------------------------------------- 1047C 1048 AUTIME = SECOND() 1049 CALL CC_T2MOT(PHONEY,FAKE,ISYMOP, 1050 * RHO2,WORK(KT2AM),DUMMY,WORK(KLAMDH), 1051 * WORK(KLAMDH),1,WORK(KEND1),LWRK1,ISYRES,1) 1052 TIM2MO = SECOND() - AUTIME 1053C 1054C----------------------------------------------------- 1055C Copy the MO vector back to the result vector. 1056C----------------------------------------------------- 1057C 1058 CALL DCOPY(NT2SQ(ISYRES),WORK(KT2AM),1,RHO2,1) 1059C 1060 IF (IPRINT .GT. 120) THEN 1061 CALL AROUND('Transformed vector after B-TERM') 1062 CALL CC_PRP(RHO1(1,IV),RHO2,ISYRES,0,1) 1063 ENDIF 1064C 1065 IF ( DEBUG ) THEN 1066 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1067 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1068 WRITE(LUPRI,1) 'Norm of RHO1 after B-TERM: ', RHO1N 1069 WRITE(LUPRI,1) 'Norm of RHO2 after B-TERM: ', RHO2N 1070 ENDIF 1071C 1072C------------------------------------------------------------------- 1073C Read integrals ( ma | nb ) stored as ( am | bn ) from disc. 1074C Ove: CCSD_IAJB is assumed open through the complete coupled 1075C cluster calculation. 1076C------------------------------------------------------------------- 1077C 1078 REWIND(LUIAJB) 1079 READ(LUIAJB) (WORK(KT2AM + J - 1), J = 1,NT2AM(ISYMOP)) 1080C 1081C------------------------------------------ 1082C Calculate the LT22AM contribution. 1083C This seems to be called the E contribution 1084C in the article..? 1085C------------------------------------------ 1086C 1087 AUTIME = SECOND() 1088 CALL CC_22AM3(RHO2,WORK(KT2AM),WORK(KMINT),ISYRES, 1089 * WORK(KEND1),LWRK1) 1090 TIMEAM = SECOND() - AUTIME 1091C 1092C------------------------------------------ 1093C Read Gamma-intermediate from disc. 1094C------------------------------------------ 1095C 1096 KGAMMI = KENDS2 1097 KENDGI = KGAMMI + NGAMMA(ISYMOP) 1098 LWRKGI = LWORK - KENDGI 1099C 1100 IF (LWRKGI .LT. 0) THEN 1101 CALL QUIT('Insufficient work space in CC_LHTR') 1102 ENDIF 1103C 1104 AUTIME = SECOND() 1105 LUGI = -1 1106 CALL GPOPEN(LUGI,'CC_GAMIM','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1107 & .FALSE.) 1108 REWIND(LUGI) 1109 READ(LUGI) (WORK(KGAMMI + J - 1), J = 1,NGAMMA(ISYMOP)) 1110 CALL GPCLOSE(LUGI,'KEEP') 1111C 1112C----------------------------------------- 1113C Calculate the LT22A contribution. 1114C----------------------------------------- 1115C 1116 ISYG = ISYMOP 1117 ISYV = ISYMTR 1118 IOPT = 2 1119C 1120C 1121 CALL CCRHS_ASQ(RHO2,CTR2,WORK(KGAMMI),WORK(KENDGI),LWRKGI, 1122 & ISYG,ISYV,IOPT) 1123 TIM22A = SECOND() - AUTIME 1124C 1125C----------------------------------------------------------------------- 1126C The above terms carries a factor of 1/2 on the (+) contribution 1127C relative to the (-) contribution. 1128C Take care of that by scaling the symmetric part py 1/2 1129C----------------------------------------------------------------------- 1130C 1131 CALL CC_T2SQSYMSCAL(RHO2,ISYMTR,HALF) 1132C 1133C------------------------------------------ 1134C Calculate the LT22EM contribution. 1135C This seems to be called the F term 1136C the article? 1137C------------------------------------------ 1138C 1139 AUTIME = SECOND() 1140 CALL CC_22E3(RHO2,WORK(KT2AM),WORK(KOFFXI),WORK(KOFFYI), 1141 * ISYRES,WORK(KEND1),LWRK1) 1142 TIM2EM = SECOND() - AUTIME 1143C 1144 IF (IPRINT .GT. 50) THEN 1145 CALL AROUND('Transformed vector after EM-TERM') 1146 CALL CC_PRP(RHO1(1,IV),RHO2,ISYRES,0,1) 1147 ENDIF 1148C 1149 IF ( DEBUG ) THEN 1150 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1151 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1152 WRITE(LUPRI,1) 'Norm of RHO1 after EM-TERM: ', RHO1N 1153 WRITE(LUPRI,1) 'Norm of RHO2 after EM-TERM: ', RHO2N 1154 ENDIF 1155C 1156C--------------------------------------------- 1157C Regain work space from T2-amplitudes. 1158C--------------------------------------------- 1159C 1160 KEND1 = KT2AM 1161 LWRK1 = LWORK - KEND1 1162C 1163C------------------------------------------ 1164C Ove: Read in AO Fock. 1165C Transform AO Fock matrix to MO basis. 1166C------------------------------------------ 1167C 1168 LFOCK = -1 1169 CALL GPOPEN(LFOCK,'CC_FCKH','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1170 & .FALSE.) 1171 REWIND(LFOCK) 1172 READ(LFOCK) (WORK(KFOCK + I - 1) , I = 1,N2BST(ISYMOP)) 1173 CALL GPCLOSE(LFOCK,'KEEP') 1174C 1175 IHELP = 1 1176C 1177 CALL CC_FCKMO(WORK(KFOCK),WORK(KLAMDP),WORK(KLAMDH), 1178 * WORK(KEND1),LWRK1,IHELP,IHELP,IHELP) 1179C 1180C------------------------------------------------ 1181C Calculate the LT11B and LT11C contribution. 1182C------------------------------------------------ 1183C 1184 CALL CCLT_11BC(RHO1(1,IV),CTR1(1,IV),WORK(KFOCK),WORK(KEND1), 1185 * LWRK1) 1186 TIM11B = SECOND() - AUTIME 1187C 1188C-------------------------------------- 1189C Calculate the LT12A (H1) contribution. 1190C-------------------------------------- 1191C 1192 AUTIME = SECOND() 1193 CALL CC_L1FCK3(RHO2,CTR1(1,IV),WORK(KFOCK),ISYMTR,ISYMOP, 1194 * WORK(KEND1),LWRK1) 1195 TIM12A = SECOND() - AUTIME 1196C 1197C----------------------------------------- 1198C LT12C / H3 contribution 1199C----------------------------------------- 1200C 1201 CALL CC_12C3(RHO2,CTR1(1,IV),ISYMTR,WORK(KLAMDH),1, 1202 & 1,WORK(KEND1),LWRK1,LUIN30,'CCINT3O') 1203C 1204C--------------------------------------------- 1205C Save RHO2 on disc to gain work space. 1206C--------------------------------------------- 1207C 1208 LURHO2 = -1 1209 CALL GPOPEN(LURHO2,'LSRHO2','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1210 & .FALSE.) 1211 REWIND(LURHO2) 1212 WRITE(LURHO2) (RHO2(I), I = 1,NT2SQ(ISYRES)) 1213C 1214C----------------------------------------------------------- 1215C Read Omega(albe,ij) written to disc by energy code. 1216C----------------------------------------------------------- 1217C 1218 TIMEBF = SECOND() 1219 LUOM = -1 1220 CALL GPOPEN(LUOM,'CC_BFIM','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1221 & .FALSE.) 1222 REWIND(LUOM) 1223 READ(LUOM) (RHO2(I),I = 1,2*NT2ORT(1) ) 1224 CALL GPCLOSE(LUOM,'KEEP') 1225C 1226C Need again C2(+) - C2(-) 1227 CALL CC_T2SQTRANSP(CTR2,ISYMTR) 1228C------------------------------------------ 1229C Calculate the LT21BF contribution. 1230C------------------------------------------ 1231C 1232 ISYM = 1 1233 ICON = 3 1234C 1235 NEWGAM = .FALSE. 1236 CALL CC_T2MO(RHO1(1,IV),CTR2,ISYMTR,RHO2,PHONEY,DUMMY, 1237 * WORK(KLAMDP),WORK(KLAMDP),ISYM,WORK(KEND1), 1238 * LWRK1,ISYMOP,ICON) 1239 NEWGAM = .TRUE. 1240 TIMEBF = SECOND() - TIMEBF 1241C 1242C----------------------------------- 1243C Restore RHO2 result vector. 1244C----------------------------------- 1245C 1246 REWIND(LURHO2) 1247 READ(LURHO2) (RHO2(I), I = 1,NT2SQ(ISYRES)) 1248 CALL GPCLOSE(LURHO2,'DELETE') 1249C 1250 IF (IPRINT .GT. 50) THEN 1251 CALL AROUND('Transformed vectors after 21BF-TERM') 1252 CALL CC_PRP(RHO1(1,IV),RHO2,ISYRES,1,1) 1253 ENDIF 1254C 1255 IF ( DEBUG ) THEN 1256 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1257 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1258 WRITE(LUPRI,1) 'Norm of RHO1 after 21BF-TERM: ', RHO1N 1259 WRITE(LUPRI,1) 'Norm of RHO2 after 21BF-TERM: ', RHO2N 1260 ENDIF 1261C 1262C----------------------------------------- 1263C Calculate the LT22D contribution. 1264C----------------------------------------- 1265C 1266 AUTIME = SECOND() 1267 IOPT = 1 1268 ISYDIM = 1 1269 FACT = HALF 1270 CALL CC_22CD3(RHO2,CTR2,ISYMTR,WORK(KLAMDH),WORK(KEND1), 1271 * LWRK1,ISYDIM,LUDIM,DFIL,1,IOPT,FACT) 1272 TIM22D = SECOND() - AUTIME 1273 CALL CC_T2SQTRANSP(CTR2,ISYMTR) 1274C 1275C----------------------------------------- 1276C Calculate the LT22C contribution. 1277C----------------------------------------- 1278C 1279 AUTIME = SECOND() 1280 IOPT = 1 1281 ISYCIM = 1 1282 FACT = -HALF 1283 1284 CALL CC_T2SQTRANSP(RHO2,ISYMTR) 1285 CALL CC_22CD3(RHO2,CTR2,ISYMTR,WORK(KLAMDH),WORK(KEND1), 1286 * LWRK1,ISYCIM,LUCIM,CFIL,1,IOPT,FACT) 1287 CALL CC_T2SQTRANSP(RHO2,ISYMTR) 1288C 1289 TIM22C = SECOND() - AUTIME 1290C 1291 IF ( DEBUG ) THEN 1292 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1293 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1294 WRITE(LUPRI,1) 'Norm of RHO1 after CC_22C: ', RHO1N 1295 WRITE(LUPRI,1) 'Norm of RHO2 after CC_22C: ', RHO2N 1296 ENDIF 1297C---------------------------------------------------------------- 1298C Extract the plus combination from the mixed doubles vector, 1299C---------------------------------------------------------------- 1300C use WORK(KEND1) as scratch 1301 CALL CC_TRIP_EXTRACT(RHO2,WORK(KEND1),ISYRES,.FALSE.) 1302C Write out (+) vector 1303 CALL CC_WVEC3(LUFR2,FRHO2,NT2AM(ISYMTR)+NT2AMA(ISYMTR), 1304 & NT2AM(ISYMTR),IV+ITR-1,0,WORK(KEND1)) 1305C 1306C Remove (+) parts from RHO2 and CTR, now they contain (-) parts only 1307C 1308 CALL CC_T2SQSYMSCAL(RHO2,ISYMTR,ZERO) 1309 CALL CC_T2SQSYMSCAL(CTR2,ISYMTR,ZERO) 1310C 1311C------------------------------------------ 1312C Calculate the LT22C- contribution. 1313C------------------------------------------ 1314C 1315 CALL CCSD_T2TP(CTR2,WORK(KEND1),LWRK1,ISYMTR) 1316 CALL CCSD_T2TP(RHO2,WORK(KEND1),LWRK1,ISYMTR) 1317C 1318 IOPT = 1 1319 ISYCIM = 1 1320 FACT = ONE 1321C 1322 CALL CC_22CD3(RHO2,CTR2,ISYMTR,WORK(KLAMDH),WORK(KEND1), 1323 * LWRK1,ISYCIM,LUCIM,CFIL,1,IOPT,FACT) 1324C 1325 CALL CCSD_T2TP(CTR2,WORK(KEND1),LWRK1,ISYMTR) 1326 CALL CCSD_T2TP(RHO2,WORK(KEND1),LWRK1,ISYMTR) 1327C 1328C------------------------------------------ 1329C Pack the RHO2(-) vector as (ai<bj) 1330C------------------------------------------ 1331C 1332 CALL CC_TRIP_EXTRACT(RHO2,WORK(KEND1),ISYRES,.TRUE.) 1333 CALL DCOPY(NT2AM(ISYMTR),WORK(KEND1),1,RHO2,1) 1334C 1335C--------------------------------------- 1336C Read E-intermediates from disc. 1337C--------------------------------------- 1338C 1339 KE1INT = KEND1 1340 KE2INT = KE1INT + NMATAB(ISYMOP) 1341 KENDTE = KE2INT + NMATIJ(ISYMOP) 1342 LWRKTE = LWORK - KENDTE 1343C 1344 IF (LWRKTE .LT. 0) THEN 1345 CALL QUIT('Insufficient work space in CC_LHTR') 1346 ENDIF 1347C 1348 AUTIME = SECOND() 1349 LUE1 = -1 1350 CALL GPOPEN(LUE1,'CC_E1IM','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1351 & .FALSE.) 1352 REWIND(LUE1) 1353 READ(LUE1) (WORK(KE1INT + J - 1), J = 1,NMATAB(ISYMOP)) 1354 CALL GPCLOSE(LUE1,'KEEP') 1355C 1356 LUE2 = -1 1357 CALL GPOPEN(LUE2,'CC_E2IM','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1358 & .FALSE.) 1359 REWIND(LUE2) 1360 READ(LUE2) (WORK(KE2INT + J - 1), J = 1,NMATIJ(ISYMOP)) 1361 CALL GPCLOSE(LUE2,'KEEP') 1362 1363C 1364C-------------------------------------------------------------- 1365C Prepare the E-intermediates for contraction with CTR2. 1366C-------------------------------------------------------------- 1367C 1368 CALL CC_EITR(WORK(KE1INT),WORK(KE2INT),WORK(KENDTE),LWRKTE, 1369 & ISYMOP) 1370C 1371C----------------------------------------- 1372C Calculate the LT22E (-) contribution. 1373C----------------------------------------- 1374C 1375 CALL CCRHS_E3(DUMMY,.FALSE.,CTR2,WORK(KE1INT),WORK(KE2INT), 1376 & WORK(KENDTE),LWRKTE,ISYMTR,ISYMOP, 1377 & RHO2,.TRUE.) 1378 TIM22E = SECOND() - AUTIME 1379C 1380 IF (IPRINT .GT. 50) THEN 1381 CALL AROUND('Transformed vector after E-TERM') 1382 CALL CC_PRP(RHO1(1,IV),RHO2,ISYRES,0,1) 1383 ENDIF 1384C 1385 IF ( DEBUG ) THEN 1386 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1387 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1388 WRITE(LUPRI,1) 'Norm of RHO1 after E-TERM: ', RHO1N 1389 WRITE(LUPRI,1) 'Norm of RHO2 after E-TERM: ', RHO2N 1390 ENDIF 1391C 1392C For now, explicitly zero (-) diagonal! 1393 IF (ISYMTR.EQ.1) CALL CCLR_DIASCL(RHO2,ZERO,ISYMTR) 1394C------------------------------------------ 1395C RHO2(-) done: write out (-) vector 1396C------------------------------------------ 1397C 1398 CALL CC_WVEC3(LUFR2,FRHO2,NT2AM(ISYMTR)+NT2AMA(ISYMTR), 1399 & NT2AM(ISYMTR),IV+ITR-1,NT2AM(ISYMTR),RHO2) 1400C 1401C-------------------------------------- 1402C Read CTR2+ from disc into RHO2. 1403C-------------------------------------- 1404C 1405 DTIME = SECOND() 1406 CALL CC_RVEC3(LUFC2,FC2AM,2*NT2AM(ISYMTR),NT2AM(ISYMTR), 1407 * IVEC+IV-1,0,RHO2) 1408 DTIME = SECOND() - DTIME 1409 TIMIO = TIMIO + DTIME 1410C 1411C-------------------------- 1412C Square up CTR2(+). 1413C-------------------------- 1414C 1415 CALL CCRHS3_R2IJ(RHO2,WORK(KENDTE),LWRKTE,ISYMTR) 1416 CALL CC_T2SQ(RHO2,CTR2,ISYMTR) 1417C----------------------------- 1418C Read RHO2 (+) back in 1419C----------------------------- 1420 CALL CC_RVEC3(LUFR2,FRHO2,NT2AM(ISYMTR)+NT2AMA(ISYMTR), 1421 & NT2AM(ISYMTR),IV+ITR-1,0,RHO2) 1422C 1423C--------------------------------------------- 1424C Calculate the LT22E (+) contribution. 1425C--------------------------------------------- 1426C This term carries a factor of 1/2: Scale E1, E2 1427 CALL DSCAL(NMATAB(1)+NMATIJ(1),HALF,WORK(KE1INT),1) 1428C 1429 AUTIME = SECOND() 1430 CALL CCRHS_E(RHO2,CTR2,WORK(KE1INT),WORK(KE2INT), 1431 & WORK(KENDTE),LWRKTE,ISYMTR,ISYMOP) 1432 TIM22E = TIM22E + SECOND() - AUTIME 1433 1434C Permute indices 1435 CALL CCRHS3_IJ(RHO2,WORK(KEND1), 1436 & LWRK1,ISYRES) 1437C Write out (+) vector 1438 CALL CC_WVEC3(LUFR2,FRHO2,NT2AM(ISYMTR)+NT2AMA(ISYMTR), 1439 & NT2AM(ISYMTR),IV+ITR-1,0,RHO2) 1440 1441 ENDIF 1442C 1443C-------------------------------------- 1444C Calculate the LT11A contribution. 1445C-------------------------------------- 1446C 1447 AUTIME = SECOND() 1448 CALL CC_11A(RHO1(1,IV),WORK(KOFFFG),ISYRES,WORK(KLAMDH), 1449 * WORK(KLAMDP),WORK(KEND1),LWRK1) 1450C 1451 IF ( DEBUG ) THEN 1452 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1453 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1454 WRITE(LUPRI,1) 'Norm of RHO1 after CC_11A: ', RHO1N 1455C WRITE(LUPRI,1) 'Norm of RHO2 after CC_11A: ', RHO2N 1456 ENDIF 1457C 1458 IF ( DEBUG ) THEN 1459 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1460 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1461 WRITE(LUPRI,1) 'Norm of RHO1 after CC_12C ', RHO1N 1462 WRITE(LUPRI,1) 'Norm of RHO2 after CC_12C ', RHO2N 1463 ENDIF 1464C 1465C---------------------------------------- 1466C Calculate the LT21EFM contribution. 1467C---------------------------------------- 1468C 1469 IF (.NOT. CC2) THEN 1470C 1471 KOFFYI = KYMAT + NMATAB(ISYRES)*(IV - 1) 1472 KOFFXI = KXMAT + NMATIJ(ISYRES)*(IV - 1) 1473 TIMEEM = SECOND() 1474 CALL CC_21EFM(RHO1(1,IV),WORK(KFOCK),ISYMOP,WORK(KOFFXI), 1475 * WORK(KOFFYI),ISYRES,WORK(KEND1),LWRK1) 1476 TIMEEM = SECOND() - TIMEEM 1477C 1478 ENDIF 1479C 1480C-------------------------------------------------------------------- 1481C Print out result vectors - zero out single vectors in CCD-calc. 1482C-------------------------------------------------------------------- 1483C 1484 IF (CCD) CALL DZERO(RHO1(1,IV),NT1AM(ISYRES)) 1485 1486 DTIME = SECOND() 1487 CALL CC_WVEC(LUFR1,FRHO1,NT1AM(ISYMTR),NT1AM(ISYMTR), 1488 * IV + IVEC -1,RHO1(1,IV)) 1489C 1490 IF (IPRINT .GT. 50) THEN 1491 CALL AROUND('Transformed vectors coming out of CC_LHTR') 1492 CALL CC_PRP(RHO1(1,IV),RHO2,ISYRES,1,1) 1493 ENDIF 1494C 1495 IF ( DEBUG ) THEN 1496 RHO1N = DDOT(NT1AM(ISYRES),RHO1(1,IV),1,RHO1(1,IV),1) 1497 RHO2N = DDOT(NT2AM(ISYRES),RHO2,1,RHO2,1) 1498 WRITE(LUPRI,1) 'Norm of RHO1 coming out of CC_LHTR:', RHO1N 1499 WRITE(LUPRI,1) 'Norm of RHO2 coming out of CC_LHTR:', RHO2N 1500 ENDIF 1501C 1502 1503 1504C 1505 DTIME = SECOND() - DTIME 1506 TIMIO = TIMIO + DTIME 1507C 1508 125 CONTINUE 1509 1510C 1511C------------------------------ 1512C Close intermediate files. 1513C------------------------------ 1514C 1515 CALL WCLOSE2(LUCIM,CFIL,'KEEP') 1516 CALL WCLOSE2(LUDIM,DFIL,'KEEP') 1517C 1518C---------------------------------------- 1519C Close intermediate files for P & Q. 1520C---------------------------------------- 1521C 1522 IF (.TRUE.) THEN 1523 CALL WCLOSE2(LUPQ,PQFIL,'KEEP') 1524 END IF 1525C 1526C------------------------------- 1527C Delete intermediate files. 1528C------------------------------- 1529C 1530 CALL WCLOSE2(LUIN30,'CCINT3O','DELETE') 1531C 1532 TIMTOT = SECOND() - TIMTOT 1533C 1534C------------------------------- 1535C Write out program timings. 1536C------------------------------- 1537C 1538 IF (IPRINT .GT. 3) THEN 1539 WRITE(LUPRI,*) ' ' 1540 WRITE(LUPRI,*) 'Time used in CC_LHTR :',TIMTOT 1541 WRITE(LUPRI,*) ' ' 1542 ENDIF 1543 IF (IPRINT .GT. 9) THEN 1544 WRITE(LUPRI,*) 'Time used as follows:' 1545 WRITE(LUPRI,*) ' ' 1546 WRITE(LUPRI,*) 'Time used in HERMIT1:',TIMHE1 1547 WRITE(LUPRI,*) 'Time used in HERMIT2:',TIMHE2 1548 WRITE(LUPRI,*) 'Time used in CCRDAO :',TIRDAO 1549 WRITE(LUPRI,*) 'Time used in Vect.IO:',TIMIO 1550 WRITE(LUPRI,*) 'Time used in CCLT_YI:',TIMEYI 1551 WRITE(LUPRI,*) 'Time used in CCLT_XI:',TIMEXI 1552 WRITE(LUPRI,*) 'Time used in CCLT_MI:',TIMEMI 1553 WRITE(LUPRI,*) 'Time used in CCLT_TI:',TIMETI 1554 WRITE(LUPRI,*) 'Time used in CCLT_Z1:',TIMEZ1 1555 WRITE(LUPRI,*) 'Time used in CCLT_Z2:',TIMEZ2 1556 WRITE(LUPRI,*) 'Time used in CCLT_C :',TIMEC 1557 WRITE(LUPRI,*) 'Time used in CCLT_A :',TIMEA 1558 WRITE(LUPRI,*) 'Time used in CCLT_I :',TIMEI 1559 WRITE(LUPRI,*) 'Time used in CCINT3O:',TIME3O 1560 WRITE(LUPRI,*) 'Time used in CCTRBT :',TIME1O 1561 WRITE(LUPRI,*) 'Time used in AOFOCK :',TIMFCK 1562 WRITE(LUPRI,*) 'Time used in CCT2TP :',TIMETP 1563 WRITE(LUPRI,*) 'Time used in CCT2AO :',TIM2AO 1564 WRITE(LUPRI,*) 'Time used in CCT2MO :',TIM2MO 1565 WRITE(LUPRI,*) 'Time used in CCLT_H :',TIMEH 1566 WRITE(LUPRI,*) 'Time used in LT_21BF:',TIMEBF 1567 WRITE(LUPRI,*) 'Time used in LT_22BF:',TIM2BF 1568 WRITE(LUPRI,*) 'Time used in CCLT_EM:',TIMEEM 1569 WRITE(LUPRI,*) 'Time used in CCLT_G :',TIMEG 1570 WRITE(LUPRI,*) 'Time used in 12C:',TIMEB 1571 WRITE(LUPRI,*) 'Time used in CC2_FCK:',TIMEC2 1572 WRITE(LUPRI,*) 'Time used in LT_12B :',TIM12B 1573 WRITE(LUPRI,*) 'Time used in LT_12A :',TIM12A 1574 WRITE(LUPRI,*) 'Time used in 11BLOCK:',TIM11B 1575 WRITE(LUPRI,*) 'Time used in LT_22AM:',TIMEAM 1576 WRITE(LUPRI,*) 'Time used in LT_22EM:',TIM2EM 1577 WRITE(LUPRI,*) 'Time used in LT_22E :',TIM22E 1578 WRITE(LUPRI,*) 'Time used in LT_22A :',TIM22A 1579 WRITE(LUPRI,*) 'Time used in LT_22C :',TIM22C 1580 WRITE(LUPRI,*) 'Time used in LT_22D :',TIM22D 1581 ENDIF 1582C 1583 1 FORMAT(1x,A35,1X,E20.10) 1584C 1585C 1586 CALL QEXIT(myname) 1587 RETURN 1588 END 1589