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 19*---------------------------------------------------------------------* 20c/* Deck CC_FTST */ 21*=====================================================================* 22 SUBROUTINE CC_FTST(WORK,LWORK) 23#if defined (IMPLICIT_NONE) 24 IMPLICIT NONE 25#else 26# include "implicit.h" 27#endif 28#include "priunit.h" 29#include "ccsdinp.h" 30#include "ccsdsym.h" 31#include "ccorb.h" 32 33* local parameters: 34 CHARACTER MSGDBG*(18) 35 PARAMETER (MSGDBG='[debug] CC_FTST> ') 36 37 LOGICAL LOCDBG 38 PARAMETER (LOCDBG = .FALSE.) 39 INTEGER MXBTRAN 40 PARAMETER (MXBTRAN = 2) 41 42 INTEGER LWORK 43#if defined (SYS_CRAY) 44 REAL WORK(LWORK) 45 REAL DDOT, RDUM 46#else 47 DOUBLE PRECISION WORK(LWORK) 48 DOUBLE PRECISION DDOT, RDUM 49#endif 50 51 CHARACTER*(3) LISTR, LISTL 52 CHARACTER*(8) FILFMA 53 CHARACTER*(10) MODEL 54 INTEGER IOPTRES 55 INTEGER IFTRAN(3,MXBTRAN), NFTRAN 56 INTEGER IDLSTL, IDLSTR, ISYML, ISYMR, ISYRES, IOPT 57 INTEGER KRESLT1, KRESLT2, KT1AMP, KT2AMP, KZETA1, KZETA2 58 INTEGER KTHETA1, KTHETA2, KEND1, LWRK1, IDUM 59 60* external function: 61 INTEGER IR1TAMP 62 INTEGER IL1ZETA 63 INTEGER ILSTSYM 64 65 66 CALL QENTER('CC_FTST') 67 68 69*---------------------------------------------------------------------* 70* call F matrix transformation: 71*---------------------------------------------------------------------* 72 LISTL = 'L0 ' 73 LISTR = 'R1 ' 74 IDLSTL = 0 75 ISYML = 1 76C IDLSTL = IL1ZETA('ZDIPLEN ',.FALSE.,0.0D0,ISYML) 77 IDLSTR = IR1TAMP('ZDIPLEN ',.FALSE.,0.0D0,ISYMR) 78 IFTRAN(1,1) = IDLSTL 79 IFTRAN(2,1) = IDLSTR 80 NFTRAN = 1 81 82 ISYRES = MULD2H(ISYML,ISYMR) 83 84 IOPTRES = 1 85 FILFMA = 'CC__FMAT' 86 87 CALL CC_FMATRIX(IFTRAN, NFTRAN, LISTL, LISTR, 88 & IOPTRES, FILFMA, IDUM, RDUM, 0, WORK, LWORK ) 89 90 IDUM = 0 91 WRITE (LUPRI,*) 'WORK(0):',WORK(IDUM) 92 DEBUG = .TRUE. 93 IPRINT = 51 94 CALL CC_FMATOLD(LISTL,IDLSTL,LISTR,IDLSTR,WORK,LWORK) 95 96 97c KTHETA1 = IBTRAN(3,1) 98c KTHETA2 = KTHETA1 + NT1AM(ISYMAB) 99 100c IF (NSYM.EQ.1 .AND. LOCDBG) THEN 101c KT1AMPB = KTHETA2 + NT2AM(ISYMAB) 102c KT2AMPB = KT1AMPB + NT1AM(ISYMB) 103c KT1AMPA = KT2AMPB + NT2AM(ISYMB) 104c KT2AMPA = KT1AMPA + NT1AM(ISYMA) 105c KRESLT1 = KT2AMPA + NT2AM(ISYMA) 106c KRESLT2 = KRESLT1 + NT1AM(ISYMAB) 107c KEND1 = KRESLT2 + NT2AM(ISYMAB) 108c LEND1 = LWORK - KEND1 109 110c IF (LEND1 .LT. 0) THEN 111c CALL QUIT('Insufficient work space in CC_FTST.') 112c END IF 113 114c IOPT = 3 115c Call CC_RDRSP(LISTA,IDLSTA,ISYMA,IOPT,MODEL, 116c & WORK(KT1AMPA),WORK(KT2AMPA)) 117c Call CC_RDRSP(LISTB,IDLSTB,ISYMB,IOPT,MODEL, 118c & WORK(KT1AMPB),WORK(KT2AMPB)) 119 120 ! zero doubles of B and/or C vector: 121C CALL DZERO(WORK(KT1AMPA),NT1AM(ISYMA)) 122C CALL DZERO(WORK(KT1AMPB),NT1AM(ISYMB)) 123C CALL DZERO(WORK(KT2AMPA),NT2AM(ISYMA)) 124C CALL DZERO(WORK(KT2AMPB),NT2AM(ISYMB)) 125c CALL DZERO(WORK(KRESLT1),NT1AM(ISYMAB)+NT2AM(ISYMAB)) 126c IPRINT = 5 127 128c CALL CC_FDB(NT1AM(ISYMAB),NT2AM(ISYMAB), 129c > WORK(KT1AMPA), WORK(KT1AMPB), WORK(KRESLT1), 130c > WORK(KEND1), LEND1) 131 132c IPRINT = 0 133 134c IF (.TRUE.) THEN 135c WRITE (LUPRI,*) 'LISTA, IDLSTA, ISYMA:',LISTA,IDLSTA,ISYMA 136c WRITE (LUPRI,*) 'LISTB, IDLSTB, ISYMB:',LISTB,IDLSTB,ISYMB 137c WRITE (LUPRI,*) 'ISYMAB:',ISYMAB 138c WRITE (LUPRI,*) 139c WRITE (LUPRI,*) 'finite difference Theta vector:' 140c Call CC_PRP(WORK(KRESLT1),WORK(KRESLT2),ISYMAB,1,1) 141c WRITE (LUPRI,*) 'analytical Theta vector:' 142c Call CC_PRP(WORK(KTHETA1),WORK(KTHETA2),ISYMAB,1,1) 143c END IF 144 145c Call DAXPY(NT1AM(ISYMAB),-1.0d0,WORK(KTHETA1),1, 146c & WORK(KRESLT1),1) 147c IF (.NOT.(CCS.OR.CCSTST)) THEN 148c Call DAXPY(NT2AM(ISYMAB),-1.0d0,WORK(KTHETA2),1, 149c & WORK(KRESLT2),1) 150c ELSE 151c Call DZERO(WORK(KRESLT2),NT2AM(ISYMAB)) 152c END IF 153 154c WRITE (LUPRI,*) 'Norm of difference between analytical THETA ' 155c > // 'vector and the numerical result:' 156c WRITE (LUPRI,*) 'singles excitation part:', 157c > DSQRT(DDOT(NT1AM(ISYMA),WORK(KRESLT1),1,WORK(KRESLT1),1)) 158c WRITE (LUPRI,*) 'double excitation part: ', 159c > DSQRT(DDOT(NT2AM(ISYMA),WORK(KRESLT2),1,WORK(KRESLT2),1)) 160 161c WRITE (LUPRI,*) 'difference vector:' 162c Call CC_PRP(WORK(KRESLT1),WORK(KRESLT2),ISYMAB,1,1) 163 164c CALL FLSHFO(LUPRI) 165 166c ELSE IF (NSYM.NE.1 .AND. LOCDBG) THEN 167c WRITE (LUPRI,*) 'CC_FTST> can not calculate finite difference B matrix' 168c WRITE (LUPRI,*) 'CC_FTST> with symmetry.' 169c END IF 170 171 CALL QEXIT('CC_FTST') 172 173 RETURN 174 END 175*=====================================================================* 176*=====================================================================* 177 SUBROUTINE CC_FMATOLD(LLIST,ILLNR,RLIST,IRLNR,WORK,LWORK) 178*--------------------------------------------------------------------* 179* 180* Purpose: transformation of trial vectors with the CC F-matrix 181* 182* Gamma(mu) = <L|[[H,R],Tau,mu]|CC> 183* 184* Left hand vector is read in according to LLIST,ILLNR. 185* Right vector according to RLIST,IRLNR. 186* 187* For LLIST .NE. 'L0' skip the HF term. 188* 189* Result vector is returned in the beginning of the work space 190* 191* Written by Ove Christiansen April-october 1996 192* Heavily revised to reduce operation count by C. Haettig July 1998 193* `Left' B intermediate revised in November 1998, C. Haettig 194*---------------------------------------------------------------------* 195#include "implicit.h" 196#include "priunit.h" 197#include "dummy.h" 198#include "maxash.h" 199#include "maxorb.h" 200#include "mxcent.h" 201#include "aovec.h" 202#include "iratdef.h" 203#include "ccorb.h" 204#include "ccisao.h" 205#include "ccsdsym.h" 206#include "ccsdinp.h" 207#include "ccfield.h" 208#include "cclr.h" 209#include "blocks.h" 210#include "ccsdio.h" 211#include "distcl.h" 212#include "cbieri.h" 213#include "eritap.h" 214#include "ccroper.h" 215#include "ccinftap.h" 216#include "r12int.h" 217#include "second.h" 218C 219 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, XMTWO= -2.0D0) 220 PARAMETER (XMHALF=-0.5D0,XHALF=0.5D0,THREE=3.0D0,XMONE=-1.0D0) 221C 222 LOGICAL NEWZWV 223 PARAMETER (NEWZWV = .TRUE.) 224C 225 INTEGER LUTR, LUPQR0, LUPQRR, LUBFI, LUBFD 226C 227 CHARACTER*(6) TRFIL, FILPQIM, FILXYM, FNBFI 228 CHARACTER*(8) FNBFD 229 PARAMETER (FILPQIM = 'CCPQ00', TRFIL = 'CC_TRA', FILXYM='CC_XYM') 230 PARAMETER (FNBFI = 'CCFBFI', FNBFD = 'CCBFDENS') 231C 232 INTEGER ISYM0 233 PARAMETER ( ISYM0 = 1 ) 234C 235 INTEGER LWORK 236 DIMENSION INDEXA(MXCORB_CC) 237 DIMENSION IADRPQ0(MXCORB_CC+IRAT), IADRPQR(MXCORB_CC+IRAT) 238 DIMENSION IADRBFI(MXCORB_CC), IADRBFD(MXCORB_CC) 239 DIMENSION WORK(*) 240 CHARACTER*6 CFIL,DFIL,CTFIL,DTFIL 241 CHARACTER*7 O3FIL,O3TFIL 242 CHARACTER LLIST*(*),RLIST*(*),MODEL*10,MODELX*10 243 LOGICAL FCKCON,ETRAN,L1TST,L2TST 244 LOGICAL T2TSAV,ORSAVE 245 LOGICAL LGAMMA, LGIM, LO3BF, LBFZETA 246 INTEGER ISYMR, ISYML, ISYRES, ISYMH, ISTARTBFI, ISTARTBFD 247 CHARACTER*6 FILPQR0, FILPQRR 248 INTEGER NBFRHF(8), IBFRHF(8,8) 249C 250 CALL QENTER('CC_FMATOLD') 251 252C----------------------------------- 253C initialize overall timing 254C----------------------------------- 255C 256 TIMALL = SECOND() 257 TIMIO = ZERO 258 TIMT2SQ = ZERO 259C 260C----------------------------------- 261C echo input 262C----------------------------------- 263C 264 IF (DEBUG) THEN 265 WRITE (LUPRI,*) 'entered CC_FMATOLD:' 266 WRITE (LUPRI,*) 'RLIST, IRLNR :',RLIST,IRLNR 267 WRITE (LUPRI,*) 'LLIST, ILLNR :',LLIST,ILLNR 268 WRITE (LUPRI,*) 'WORK(1):',WORK(1) 269 WRITE (LUPRI,*) 'LWORK:',LWORK 270 END IF 271C 272C--------------------------------------------- 273C Set symmetries, test flags and CC model: 274C--------------------------------------------- 275C 276 ISYML = ILSTSYM(LLIST,ILLNR) 277 ISYMR = ILSTSYM(RLIST,IRLNR) 278 ISYRES = MULD2H(ISYML,ISYMR) 279C 280C ---------------------------------------------------------- 281C initialize file addresses for BF density and intermediate: 282C ---------------------------------------------------------- 283C 284 ISTARTBFI = 1 285 ISTARTBFD = 1 286C 287C -------------------------------------- 288C precalculate symmetry array for BFRHF: 289C -------------------------------------- 290C 291 DO ISYM = 1, NSYM 292 ICOUNT = 0 293 DO ISYMAK = 1, NSYM 294 ISYBET = MULD2H(ISYMAK,ISYM) 295 IBFRHF(ISYMAK,ISYBET) = ICOUNT 296 ICOUNT = ICOUNT + NT1AO(ISYMAK)*NBAS(ISYBET) 297 END DO 298 NBFRHF(ISYM) = ICOUNT 299 END DO 300C 301 L1TST = .FALSE. 302 L2TST = .FALSE. 303C 304 MODEL = 'CCSD ' 305 IF (CCS) MODEL = 'CCS ' 306 IF (CC2) MODEL = 'CC2 ' 307C 308 IF (ISYRES .NE. MULD2H(ISYML,ISYMR)) THEN 309 WRITE(LUPRI,*) 'Symmetry mismatch in CC_FMATOLD' 310 CALL QUIT('Symmetry mismatch in CC_FMATOLD') 311 ENDIF 312C 313 IF (IPRINT .GT. 15) THEN 314 CALL AROUND(' START OF CC_FMATOLD ') 315 IF (DIRECT) WRITE(LUPRI,*) ' Atomic direct transformation' 316 ENDIF 317C 318 IF ( CCSDT ) THEN 319 WRITE(LUPRI,*) ' F-Matrix transformations not implemented ' 320 * //'for triples yet: Go and do it.' 321 CALL QUIT(' no triples F-matrix') 322 ENDIF 323C 324C--------------------------------------- 325C Set common block control variables 326C--------------------------------------- 327C 328 T2TSAV = T2TCOR 329 IF (CCS .OR. CC2) T2TCOR = .FALSE. 330C COMMENTED OUT! OBSOLETE SUBROUTINE 331C IF (CC2 .AND.(NSIDE .EQ. 2)) T2TCOR = T2TSAV 332C 333 ORSAVE = OMEGOR 334 IF (CCS .OR. CC2) THEN 335 OMEGOR = .FALSE. 336 ELSE 337 OMEGOR = .TRUE. 338 ENDIF 339C 340 IF (.NOT.DUMPCD) THEN 341 CALL QUIT('CC_FMATOLD requires DUMPCD=.TRUE.') 342 END IF 343C 344C-------------------------------------------------------------- 345C open files for result vector and for BFZeta intermediate: 346C-------------------------------------------------------------- 347C 348 IERRTR = 0 349 IERRBF = 0 350C 351 LUTR = -1 352 IF ( .NOT. CCS ) THEN 353 CALL WOPEN2(LUTR, TRFIL,64,0) 354 ENDIF 355C 356 LUBFI = -1 357 LUBFD = -1 358 IF (CCSD) THEN 359 CALL WOPEN2(LUBFI,FNBFI,64,0) 360 CALL WOPEN2(LUBFD,FNBFD,64,0) 361 END IF 362C 363*--------------------------------------------------------------------* 364* calculate P and Q intermediates, write them to file and 365* open the file for read within the loop over AO integrals 366*--------------------------------------------------------------------* 367 TIMEZ = SECOND() 368 369 IF (.NOT. (CCS.OR.CC2)) THEN 370 371 IF (LLIST(1:2).EQ.'L0') THEN 372 FILPQR0 = FILPQIM 373 ELSE 374 FILPQR0 = 'CCPQR0' 375 CALL CC_PQIM(LLIST,ILLNR,'R0',0,FILPQR0,WORK,LWORK) 376 END IF 377 378 FILPQRR = 'CCPQRR' 379 CALL CC_PQIM(LLIST,ILLNR,RLIST,IRLNR,FILPQRR,WORK,LWORK) 380 381C ----------------------------------------- 382C open file for P and Q intermediates: 383C ----------------------------------------- 384 LUPQR0 = -1 385 LUPQRR = -1 386 CALL WOPEN2(LUPQR0,FILPQR0,64,0) 387 CALL WOPEN2(LUPQRR,FILPQRR,64,0) 388 389 CALL GETWA2(LUPQR0,FILPQR0,IADRPQ0,1,(MXCORB_CC+IRAT)/IRAT) 390 CALL GETWA2(LUPQRR,FILPQRR,IADRPQR,1,(MXCORB_CC+IRAT)/IRAT) 391 392 END IF 393 394 TIMEZ = SECOND() - TIMEZ 395 396C 397C-------------------------------------------------------------------- 398C Allocate space for solution vector and right hand trial vector. 399C-------------------------------------------------------------------- 400C 401 IF (CCS) THEN 402 NRHO2 = 2 403 ELSE IF (CC2) THEN 404 NRHO2 = MAX(NT2AM(ISYRES),NT2AM(ISYMR),NT2AM(ISYML),NT2AM(1)) 405 ELSE 406 NRHO2 = MAX(NT2AM(ISYRES),NT2AM(ISYMR),NT2AM(ISYML),NT2AM(1), 407 * 2*NT2ORT(ISYRES),2*NT2ORT(ISYMR),NT2AO(ISYRES)) 408 ENDIF 409C 410 NC2AM = MAX(NT2SQ(ISYMR),NT2SQ(1),NT2AM(ISYMR)+2*NT2ORT(ISYMR), 411 * NT2AM(1)+2*NT2ORT(1),NT2AM(ISYRES)) 412 IF ( CC2 ) NC2AM = MAX(NT2SQ(ISYMR),NT2SQ(1)) 413 IF ( CCS ) NC2AM = 2 414C 415 KRHO1 = 1 416 KRHO2 = KRHO1 + NT1AM(ISYRES) 417 KC1AM = KRHO2 + NRHO2 418 KC2AM = KC1AM + NT1AM(ISYMR) 419 KL1AM = KC2AM + NC2AM 420 KL2AM = KL1AM + NT1AM(ISYML) 421 KL1R2 = KL2AM + NT2SQ(ISYML) 422 KL1R2PK = KL1R2 + N2BST(ISYRES) 423 KT1AM = KL1R2PK + NNBST(ISYRES) 424 KEND0 = KT1AM + NT1AM(ISYM0) 425 LWRK0 = LWORK - KEND0 426 IF (LWRK0 .LE. 0 ) THEN 427 CALL QUIT('Too little workspace in CC_FMAT.') 428 END IF 429C 430C------------------------------------------------------------------- 431C Read the lagrangian multiplier vector from disk and square up: 432C------------------------------------------------------------------- 433C 434 IF (.NOT. ( CCS .AND. LLIST(1:2).EQ.'L0') ) THEN 435 DTIME = SECOND() 436 IOPT = 3 437 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODEL, 438 * WORK(KL1AM),WORK(KRHO2)) 439 DTIME = SECOND() - DTIME 440 TIMIO = TIMIO + DTIME 441 END IF 442 443 IF (.NOT. CCS) THEN 444 DTIME = SECOND() 445 CALL CC_T2SQ(WORK(KRHO2),WORK(KL2AM),ISYML) 446 DTIME = SECOND() - DTIME 447 TIMT2SQ = TIMT2SQ + DTIME 448 END IF 449C 450C ------------------ 451C Test options. 452C ------------------ 453C 454 IF (L1TST .AND. (.NOT. CCS)) CALL DZERO(WORK(KL2AM),NT2SQ(ISYML)) 455 IF (L2TST) CALL DZERO(WORK(KL1AM),NT1AM(ISYML)) 456 457 IF ( DEBUG ) THEN 458 RHO1N = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM ),1) 459 WRITE(LUPRI,1) 'Norm of L1AM as read from file: ',RHO1N 460 ENDIF 461 IF ( DEBUG .AND. (.NOT. CCS)) THEN 462 RHO2N = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 463 WRITE(LUPRI,1) 'Norm of L2AM -after square: ',RHO2N 464 ENDIF 465C 466C----------------------------------------------------------------- 467C Read the CC reference amplitudes from disk, T2AMP0 in KRHO2 468C (double excitation part only needed if LLIST .NE. 'L0', 469C i.e. if X, Y and M intermediates have to be calculated) 470C----------------------------------------------------------------- 471C 472 IOPT = 1 473 IF (LLIST(1:2).NE.'L0') IOPT = 3 474 475 DTIME = SECOND() 476 CALL CC_RDRSP('R0',0,1,IOPT,MODELX,WORK(KT1AM),WORK(KRHO2)) 477 DTIME = SECOND() - DTIME 478 TIMIO = TIMIO + DTIME 479C 480C------------------------------------------------------------------- 481C initialize C1 and C2 vectors: 482C we start the calculations with C2 vector packed in WORK(KC2AM) 483C------------------------------------------------------------------- 484C 485 DTIME = SECOND() 486 IOPT = 3 487 CALL CC_RDRSP(RLIST,IRLNR,ISYMR,IOPT,MODELX, 488 * WORK(KC1AM),WORK(KC2AM)) 489 DTIME = SECOND() - DTIME 490 TIMIO = TIMIO + DTIME 491 492 IF (.NOT. CCS) CALL CCLR_DIASCL(WORK(KC2AM),TWO,ISYMR) 493 494 IF ( DEBUG ) THEN 495 RHO1N = DDOT(NT1AM(ISYMR),WORK(KC1AM),1,WORK(KC1AM),1) 496 WRITE(LUPRI,1) 'CC_FMAT: Norm of T1-response vector:',RHO1N 497 IF (.NOT. CCS) THEN 498 RHO2N = DDOT(NT2AM(ISYMR),WORK(KC2AM),1,WORK(KC2AM),1) 499 WRITE(LUPRI,1) 'CC_FMAT: Norm of T2-response vector:',RHO2N 500 ENDIF 501 ENDIF 502C 503C------------------------------------------------------------------- 504C initialize result vector: 505C (note that double-excitation result vector is kept on file...) 506C------------------------------------------------------------------- 507C 508 CALL DZERO(WORK(KRHO1),NT1AM(ISYRES)) 509C 510C------------------------------------------------------------ 511C Calculate <HF| [[H,tau,mu],R1]|HF> contribution. 512C Note inside FHF a Integral(ai,bj) is allocated. 513C Of course released afterwards. 514C------------------------------------------------------------ 515C 516 IF ((LLIST(1:2).EQ.'L0')) THEN 517 518 CALL CC_FHF(WORK(KC1AM),WORK(KRHO1),ISYMR,WORK(KEND0),LWRK0) 519 520 IF (DEBUG) THEN 521 WRITE(LUPRI,*) 522 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 523 WRITE(LUPRI,1) 'Norm of Rho1 -after FHF: ',RHO1N 524 ENDIF 525 526 ENDIF 527C 528C------------------------------------------ 529C If CCS all is done and we can return. 530C------------------------------------------ 531C 532 IF ( CCS .AND. (LLIST(1:2).EQ.'L0')) THEN 533 IF (IPRINT .GT. 50 ) THEN 534 CALL AROUND('END OF CC_FMAT :RHO ') 535 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,1,0) 536 ENDIF 537 RETURN 538 ENDIF 539C 540C---------------- 541C Open files. 542C---------------- 543C 544 IF (.NOT.(CCS .OR. CC2)) THEN 545 LUC = -1 546 LUD = -1 547C 548 CTFIL = 'CCLR_C' 549 DTFIL = 'CCLR_D' 550C 551 CALL WOPEN2(LUC,CTFIL,64,0) 552 CALL WOPEN2(LUD,DTFIL,64,0) 553C 554 LUCIM = -1 555 LUDIM = -1 556C 557 CFIL = 'PMAT_C' 558 DFIL = 'PMAT_D' 559C 560 CALL WOPEN2(LUCIM,CFIL,64,0) 561 CALL WOPEN2(LUDIM,DFIL,64,0) 562C 563 END IF 564C 565C----------------------------------------- 566C Open scratch file CCINT3O and O3TFIL 567C----------------------------------------- 568C 569 LUO3 = -1 570 O3FIL = 'CCINT3O' 571C 572 CALL WOPEN2(LUO3,O3FIL,64,0) 573C 574 LUO3T = -1 575 O3TFIL = 'CCO3INT' 576C 577 CALL WOPEN2(LUO3T,O3TFIL,64,0) 578C 579C------------------------ 580C Time Initialiation. 581C------------------------ 582C 583 TIMER1 = 0.0D00 584 TIMER2 = 0.0D00 585 TIMRDAO = 0.0D00 586C 587 TIMA = 0.0D00 588 TIMBF = 0.0D00 589 TIMC = 0.0D00 590 TIMD = 0.0D00 591 TIME = 0.0D00 592 TIMF = 0.0D00 593 TIMG = 0.0D00 594 TIMH = 0.0D00 595 TIMI = 0.0D00 596 TIMEI = 0.0D00 597 TIMEXI = 0.0D00 598 TIMEYI = 0.0D00 599 TIMEMI = 0.0D00 600 TIMEH = 0.0D00 601 TIMEFCK = 0.0D00 602 TIME1C1 = 0.0D00 603 TIM1C1F = 0.0D00 604 TIMCIO = 0.0D00 605 TIMDIO = 0.0D00 606 TIM22D = 0.0D00 607 TIMEC = 0.0D00 608 TIMEG = 0.0D00 609 TIM2EM = 0.0D00 610 TIME3O = 0.0D00 611C 612 TIMLAM = 0.0D00 613 TIMAOD = 0.0D00 614 TIMFCK = 0.0D00 615 TIMTRBT = 0.0D00 616 TIMT2AO = 0.0D00 617 TIMTCME = 0.0D00 618 TIMT2TR = 0.0D00 619 TIMT2SQ = 0.0D00 620 TIMT2BT = 0.0D00 621 TIMT2MO = 0.0D00 622 TIMFCKMO = 0.0D00 623 TIMC1T2 = 0.0D00 624 TIM12B = 0.0D00 625C 626 TIMT3I = 0.0D00 627 TIMO31 = 0.0D00 628 TIMO32 = 0.0D00 629 TIMO33 = 0.0D00 630C 631 TIMIO = 0.0D00 632C 633C----------------------------------------------------- 634C Work space allocation for general intermediates. 635C----------------------------------------------------- 636C 637 NE1TOT = MAX(NEMAT1(ISYMR),NEMAT1(1)) 638 NE2TOT = MAX(NMATIJ(ISYMR),NMATIJ(1)) 639 N2C2C2 = 0 640 IF ( T2TCOR ) N2C2C2 = NT2SQ(ISYMR) 641 IF ( CCS ) THEN 642 N2C2C2 = 2 643 ENDIF 644C 645 IF ( IPRINT .GT. 25) 646 * WRITE(LUPRI,*) ' In CC_FMAT work in start is:',LWORK 647C 648 K2C2C2 = KEND0 649 KLAMDP = K2C2C2 + N2C2C2 650 KLAMDH = KLAMDP + NLAMDT 651 KDENSI = KLAMDH + NLAMDT 652 KDNSPKI = KDENSI + N2BAST 653 KFOCK = KDNSPKI + NNBST(ISYMOP) 654 KEMAT1 = KFOCK + N2BST(ISYMOP) 655 KEMAT2 = KEMAT1 + NE1TOT 656 KEND1 = KEMAT2 + NE2TOT 657 LWRK1 = LWORK - KEND1 658C 659C---------------------------------------------------------------- 660C Extra Work space allocation for C1 dependent intermediates. 661C---------------------------------------------------------------- 662C 663 KLAMPC = KEND1 664 KLAMHC = KLAMPC + NGLMDT(ISYMR) 665 KDENSC = KLAMHC + NGLMDT(ISYMR) 666 KDNSPKC = KDENSC + N2BST(ISYMR) 667 KFOCKC = KDNSPKC + NNBST(ISYMR) 668 KFOCKT = KFOCKC + N2BST(ISYMR) 669 KFCKLR = KFOCKT + MAX(N2BST(ISYMR),N2BST(ISYMOP)) 670 KL1AOC = KFCKLR + N2BST(ISYRES) 671 KRHO1AO = KL1AOC + NGLMDT(ISYRES) 672 KEND1 = KRHO1AO + NT1AO(ISYRES) 673 LWRK1 = LWORK - KEND1 674C 675 NRHOR = NT2AOIJ(ISYMR) 676 NMGD = NT2AOIJ(ISYRES) 677C 678 IF (.NOT.((CC2.AND.(.NOT.NONHF)).OR.CCS)) THEN 679 KXMAT = KEND1 680 KXMATX = KXMAT + NMATIJ(ISYML) 681 KYMAT = KXMATX + NMATIJ(ISYRES) 682 KYMATX = KYMAT + NMATAB(ISYML) 683 KYDENX = KYMATX + NMATAB(ISYRES) 684 KYDENB = KYDENX + N2BST(ISYRES) 685 KMINT = KYDENB + N2BST(ISYRES) 686 KMINTX = KMINT + N3ORHF(ISYML) 687 KCHIM = KMINTX + N3ORHF(ISYRES) 688 KRHOR = KCHIM + NMATIJ(ISYRES) 689 KMGDL = KRHOR + NRHOR 690 KEND1 = KMGDL + NMGD 691 LWRK1 = LWORK - KEND1 692 ENDIF 693C 694 IF ( IPRINT .GT. 25) 695 * WRITE(LUPRI,*) ' In CC_FMAT 1. alloc. work. left:',LWRK1 696C 697 IF (LWRK1 .LE. 0) CALL QUIT('Too little workspace in CC_FMAT ') 698C 699C--------------------------------------- 700C Initialize t1am and intermediates. 701C--------------------------------------- 702C 703C CALL DZERO(WORK(KT1AM),NT1AM(1)) 704 CALL DZERO(WORK(KEMAT1),NE1TOT) 705 CALL DZERO(WORK(KEMAT2),NE2TOT) 706 CALL DZERO(WORK(KDENSI),N2BST(ISYMOP)) 707 CALL DZERO(WORK(KFOCK),N2BST(ISYMOP)) 708C 709C------------------------------------------- 710C Initialize C1 dependent intermediates. 711C------------------------------------------- 712C 713 CALL DZERO(WORK(KL1R2),N2BST(ISYRES)) 714 CALL DZERO(WORK(KFCKLR),N2BST(ISYRES)) 715 CALL DZERO(WORK(KFOCKC),N2BST(ISYMR)) 716 CALL DZERO(WORK(KDENSC),N2BST(ISYMR)) 717 CALL DZERO(WORK(KRHO1AO),NT1AO(ISYRES)) 718C 719C---------------------------------- 720C Calculate the lamda matrices: 721C---------------------------------- 722C 723 DTIME = SECOND() 724 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 725 * WORK(KEND1),LWRK1) 726 727 CALL CCLR_LAMTRA(WORK(KLAMDP),WORK(KLAMPC), 728 * WORK(KLAMDH),WORK(KLAMHC),WORK(KC1AM),ISYMR) 729 730 DTIME = SECOND() - DTIME 731 TIMLAM = DTIME + TIMLAM 732 733 IF (IPRINT .GT.45) THEN 734 CALL AROUND('Usual Lambda matrices ') 735 CALL CC_PRLAM(WORK(KLAMDP),WORK(KLAMDH),ISYM0) 736 CALL AROUND('C1 transformed Lambda matrices ') 737 CALL CC_PRLAM(WORK(KLAMPC),WORK(KLAMHC),ISYMR) 738 ENDIF 739C 740 IF (.NOT.(CCS.OR.(CC2.AND.(.NOT.NONHF)))) THEN 741C 742 IF ( LLIST(1:2) .EQ. 'L0' ) THEN 743C 744C --------------------------------------------------- 745C Read X, Y and M intermediates from 746C zeroth-order Zeta and T2AMP0 from file 747C --------------------------------------------------- 748C 749 DTIME = SECOND() 750 LUXYM = -1 751 CALL GPOPEN(LUXYM,FILXYM,'OLD',' ','UNFORMATTED',IDUMMY, 752 & .FALSE.) 753 754 REWIND(LUXYM,ERR=999) 755 756 IF (.NOT.(CCS.OR.(CC2.AND.(.NOT.NONHF)))) THEN 757 READ(LUXYM,ERR=999) (WORK(KYMAT-1+I),I=1,NMATAB(ISYML)) 758 READ(LUXYM,ERR=999) (WORK(KXMAT-1+I),I=1,NMATIJ(ISYML)) 759 END IF 760 761 IF (.NOT.(CCS.OR.(CC2.AND.(.NOT.NONHF)))) THEN 762 READ(LUXYM,ERR=999) (WORK(KMINT-1+I),I=1,N3ORHF(ISYML)) 763 END IF 764 765 CALL GPCLOSE(LUXYM,'KEEP') 766 DTIME = SECOND() - DTIME 767 TIMIO = TIMIO + DTIME 768 769 ELSE 770C 771C --------------------------------------------------- 772C Calculate X, Y and M intermediates, from 773C response Zeta and T2AMP0 774C --------------------------------------------------- 775C 776 TIMEYI = SECOND() 777 CALL CC_YI(WORK(KYMAT),WORK(KL2AM),ISYML,WORK(KRHO2),ISYM0, 778 & WORK(KEND1),LWRK1) 779 TIMEYI = SECOND() - TIMEYI 780 781 TIMEXI = SECOND() 782 CALL CC_XI(WORK(KXMAT),WORK(KL2AM),ISYML,WORK(KRHO2),ISYM0, 783 & WORK(KEND1),LWRK1) 784 TIMEXI = SECOND() - TIMEXI 785 786 IF (.NOT.(CCS.OR.(CC2.AND.(.NOT.NONHF)))) THEN 787 TIMEMI = SECOND() 788 CALL CC_MI(WORK(KMINT),WORK(KL2AM),ISYML,WORK(KRHO2), 789 & ISYM0,WORK(KEND1),LWRK1) 790 TIMEMI = SECOND() - TIMEMI 791 ENDIF 792 793 END IF 794C 795C ----------------------------------- 796C Calculate X, Y and M intermediates: 797C from Zeta and response T2AMP 798C ----------------------------------- 799C 800 DTIME = SECOND() 801 CALL CC_YI(WORK(KYMATX),WORK(KL2AM),ISYML, 802 * WORK(KC2AM),ISYMR,WORK(KEND1),LWRK1) 803 DTIME = SECOND() - DTIME 804 TIMEYI = TIMEYI + DTIME 805C 806 DTIME = SECOND() 807 CALL CC_XI(WORK(KXMATX),WORK(KL2AM),ISYML, 808 * WORK(KC2AM),ISYMR,WORK(KEND1),LWRK1) 809 DTIME = SECOND() - DTIME 810 TIMEXI = TIMEXI + DTIME 811C 812 IF (.NOT.(CCS.OR.(CC2.AND.(.NOT.NONHF)))) THEN 813 DTIME = SECOND() 814 CALL CC_MI(WORK(KMINTX),WORK(KL2AM),ISYML, 815 * WORK(KC2AM),ISYMR,WORK(KEND1),LWRK1) 816 DTIME = SECOND() - DTIME 817 TIMEMI = TIMEMI + DTIME 818 ENDIF 819C 820C ---------------------- 821C calculate Y densities: 822C ---------------------- 823C 824 DTIME = SECOND() 825 CALL CC_YD(WORK(KYDENB),WORK(KYMAT),ISYML,WORK(KLAMDH), 826 * WORK(KLAMPC),ISYMR,WORK(KEND1),LWRK1) 827C 828 CALL CC_YD(WORK(KYDENX),WORK(KYMATX),ISYRES,WORK(KLAMDH), 829 * WORK(KLAMDP),ISYMOP,WORK(KEND1),LWRK1) 830 DTIME = SECOND() - DTIME 831 TIMEYI = TIMEYI + DTIME 832C 833C ------------------------------ 834C calculate response Chi matrix: 835C (not backtransformed to AO) 836C ------------------------------ 837C 838 DO ISYMK = 1, NSYM 839 ISYMC = MULD2H(ISYMK,ISYMR) 840 ISYMI = MULD2H(ISYMC,ISYML) 841 842 KOFF1 = KC1AM + IT1AM(ISYMC,ISYMK) 843 KOFF2 = KL1AM + IT1AM(ISYMC,ISYMI) 844 KOFF3 = KCHIM + IMATIJ(ISYMK,ISYMI) 845 846 NRHFK = MAX(NRHF(ISYMK),1) 847 NVIRC = MAX(NVIR(ISYMC),1) 848 849 CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC), 850 * -ONE,WORK(KOFF1),NVIRC,WORK(KOFF2),NVIRC, 851 * ZERO,WORK(KOFF3),NRHFK) 852 853 END DO 854 CALL DAXPY(NMATIJ(ISYRES),-ONE,WORK(KXMATX),1,WORK(KCHIM),1) 855C 856C ------------------------------------------------------- 857C precalculate effective density for left B intermediate: 858C ------------------------------------------------------- 859C 860 IF (CCSD) THEN 861 862 DTIME = SECOND() 863 CALL CC_BFDENF(WORK(KL2AM),ISYML,WORK(KMINTX),ISYRES, 864 * WORK(KLAMDP),ISYMOP,WORK(KCHIM),ISYRES, 865 * WORK(KC1AM),ISYMR,WORK(KMGDL), 866 * WORK(KEND1),LWRK1) 867 TIMBF = TIMBF + SECOND() - DTIME 868 869 END IF 870 871 ENDIF 872C 873C-------------------------------- 874C Calculate L1R2 contraction: 875C-------------------------------- 876C 877 IOPT = 3 878 TIMC1T2 = SECOND() 879 IF (DEBUG) THEN 880 RHO1N = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM ),1) 881 WRITE(LUPRI,1) 'Norm of L1AM before CC_C1T2C: ',RHO1N 882 RHO1N = DDOT(NT2AM(ISYMR),RHO2,1,RHO2,1) 883 WRITE(LUPRI,1) 'Norm of c2am before C1T2C: ',RHO1N 884 RHO1N = DDOT(NGLMDT(ISYMOP),WORK(KLAMDH),1,WORK(KLAMDH),1) 885 WRITE(LUPRI,1) 'Norm of lamdaH before C1T2C: ',RHO1N 886 RHO1N = DDOT(NGLMDT(ISYMR ),WORK(KLAMPC),1,WORK(KLAMPC),1) 887 WRITE(LUPRI,1) 'Norm of lamdPC before C1T2C: ',RHO1N 888 RHO1N = DDOT(NGLMDT(ISYMR ),WORK(KLAMHC),1,WORK(KLAMHC),1) 889 WRITE(LUPRI,1) 'Norm of lamdHC before C1T2C: ',RHO1N 890 RHO1N = DDOT(NGLMDT(ISYMOP),WORK(KLAMDP),1,WORK(KLAMDP),1) 891 WRITE(LUPRI,1) 'Norm of lamdaP before C1T2C: ',RHO1N 892 ENDIF 893 CALL CC_C1T2C(WORK(KL1R2),WORK(KL1AM),WORK(KC2AM), 894 * WORK(KLAMDP),WORK(KLAMDH), 895 * WORK(KLAMPC),WORK(KLAMHC), 896 * WORK(KEND1),LWRK1,ISYML,ISYMR,IOPT) 897 TIMC1T2 = SECOND() - TIMC1T2 898 IF (DEBUG) THEN 899 RHO2N = DDOT(NT2AM(ISYMR),WORK(KC2AM),1,WORK(KC2AM),1) 900 RHO1N = DDOT(N2BST(ISYRES),WORK(KL1R2),1,WORK(KL1R2),1) 901 WRITE(LUPRI,1) 'Norm of RHO2 after CC_C1T2C: ',RHO2N 902 WRITE(LUPRI,1) 'Norm of FCKR2 after CC_C1T2C: ',RHO1N 903 ENDIF 904C 905C ------------------------------------------------------------ 906C Include the two LT21A contributions here by 2N^2 operations: 907C ------------------------------------------------------------ 908C 909 IF (.NOT.(CC2.OR.CCS)) THEN 910 911 CALL DAXPY(N2BST(ISYRES),ONE,WORK(KYDENX),1,WORK(KL1R2),1) 912 913 CALL DAXPY(N2BST(ISYRES),ONE,WORK(KYDENB),1,WORK(KL1R2),1) 914 915 END IF 916 917 CALL CC_DNSPK(WORK(KL1R2),WORK(KL1R2PK),ISYRES) 918C 919C------------------------------------------------------------------ 920C Calculate the density matrices: 921C IC=1 include core contribution for zeroth-order matrix DENSI 922C IC=0 no core contribution for C1 transformed matrix DENSC 923C------------------------------------------------------------------ 924C 925 DTIME = SECOND() 926C 927 IC = 1 928 CALL CC_AODENS(WORK(KLAMDP),WORK(KLAMDH),WORK(KDENSI),ISYM0, 929 * IC,WORK(KEND1),LWRK1) 930 CALL CC_DNSPK(WORK(KDENSI),WORK(KDNSPKI),ISYM0) 931C 932 IC = 0 933 CALL CC_AODENS(WORK(KLAMDP),WORK(KLAMHC),WORK(KDENSC),ISYMR, 934 * IC,WORK(KEND1),LWRK1) 935 CALL CC_DNSPK(WORK(KDENSC),WORK(KDNSPKC),ISYMR) 936C 937 DTIME = SECOND() - DTIME 938 TIMAOD = TIMAOD + DTIME 939C 940 IF (IPRINT .GT. 45) THEN 941 CALL AROUND('CC_FMAT : Usual Lamda density matrix') 942 CALL CC_PRAODEN(WORK(KDENSI),ISYM0) 943 CALL AROUND('CC_FMAT : C1 trans. Lamda density matrix') 944 CALL CC_PRAODEN(WORK(KDENSC),ISYMR) 945 ENDIF 946C 947C-------------------------------------------- 948C Calculate L1LamdahaC1 Lambda vector. 949C-------------------------------------------- 950C 951 CALL CCLL_LAMTRA(WORK(KL1AM),ISYML,WORK(KLAMPC), 952 * ISYMR,WORK(KL1AOC)) 953C 954 IF (IPRINT .GT. 45) THEN 955 CALL AROUND('L1 transformed Lambda C1 P matrix ') 956 CALL CC_PRLAM(WORK(KL1AOC),WORK(KLAMHC),ISYRES) 957 ENDIF 958C 959C-------------------------------------------------------------- 960C Prepare C2 vector and its transposed counterpart in core. 961C-------------------------------------------------------------- 962C 963 IF ( .NOT. CCS ) THEN 964 965 IOPT = 2 966 CALL CC_RDRSP(RLIST,IRLNR,ISYMR,IOPT,MODELX, 967 * DUMMY,WORK(KRHO2)) 968 CALL CCLR_DIASCL(WORK(KRHO2),TWO,ISYMR) 969 CALL CC_T2SQ(WORK(KRHO2),WORK(KC2AM),ISYMR) 970 971 IF (T2TCOR) THEN 972 DTIME = SECOND() 973 CALL DCOPY(NT2SQ(ISYMR),WORK(KC2AM),1,WORK(K2C2C2),1) 974 CALL CCSD_T2TP(WORK(K2C2C2),WORK(KEND1),LWRK1,ISYMR) 975 DTIME = SECOND() - DTIME 976 TIMT2TR = TIMT2TR + DTIME 977 ENDIF 978 979 ENDIF 980C 981C -------------------------------------------------------- 982C precalculate effective density for right B intermediate: 983C -------------------------------------------------------- 984C 985 IF (CCSD) THEN 986 987 IVEC = 1 988 IOPT = 3 989 DTIME = SECOND() 990 CALL CC_BFDEN(WORK(KC2AM),ISYMR,DUMMY,IDUMMY, 991 * WORK(KLAMDH),ISYMOP,WORK(KLAMDH),ISYMOP, 992 * WORK(KLAMHC),ISYMR,DUMMY,IDUMMY, 993 * FNBFD,LUBFD,IADRBFD,ISTARTBFD, 994 * IVEC, IOPT, .FALSE., WORK(KEND1), LWRK1) 995 TIMBF = TIMBF + SECOND() - DTIME 996 997 END IF 998C 999C------------------------------------------------ 1000C Read one-electron integrals in Fock-matrix. 1001C------------------------------------------------ 1002C 1003 TIMFCK = SECOND() 1004 CALL CCRHS_ONEAO(WORK(KFOCK),WORK(KEND1),LWRK1) 1005 DO IF = 1, NFIELD 1006 FF = EFIELD(IF) 1007 CALL CC_ONEP(WORK(KFOCK),WORK(KEND1),LWRK1,FF,1,LFIELD(IF)) 1008 END DO 1009 TIMFCK = SECOND() - TIMFCK 1010C 1011C------------------------------------------- 1012C initialize B intermediates with zeros: 1013C------------------------------------------- 1014C 1015 IF (.NOT. CCS) THEN 1016 CALL DZERO(WORK(KRHO2),NRHO2) 1017 CALL CC_WVEC(LUTR,TRFIL,NRHO2,NRHO2,1,WORK(KRHO2)) 1018 ENDIF 1019 IF (.NOT. (CC2.OR.CCS)) THEN 1020 CALL DZERO(WORK(KRHOR),NRHOR) 1021 ENDIF 1022C 1023C==================================================== 1024C Start the loop over distributions of integrals. 1025C==================================================== 1026C 1027 KENDS2 = KEND1 1028 LWRKS2 = LWRK1 1029C 1030 IF (DIRECT) THEN 1031 NTOSYM = 1 1032 DTIME = SECOND() 1033C 1034 IF (HERDIR) THEN 1035 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 1036 ELSE 1037 KCCFB1 = KEND1 1038 KINDXB = KCCFB1 + MXPRIM*MXCONT 1039 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 1040 LWRK1 = LWORK - KEND1 1041C 1042 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 1043 & KODPP1,KODPP2,KRDPP1,KRDPP2, 1044 & KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 1045 & WORK(KEND1),LWRK1,IPRERI) 1046 KEND1 = KFREE 1047 LWRK1 = LFREE 1048 END IF 1049 1050 DTIME = SECOND() - DTIME 1051 TIMER1 = TIMER1 + DTIME 1052 ELSE 1053 NTOSYM = NSYM 1054 ENDIF 1055C 1056 KENDSV = KEND1 1057 LWRKSV = LWRK1 1058C 1059 ICDEL1 = 0 1060 ICDEL2 = 0 1061C 1062 DO 100 ISYMD1 = 1,NTOSYM 1063C 1064 IF (DIRECT) THEN 1065 IF (HERDIR) THEN 1066 NTOT = MAXSHL 1067 ELSE 1068 NTOT = MXCALL 1069 ENDIF 1070 ELSE 1071 NTOT = NBAS(ISYMD1) 1072 ENDIF 1073C 1074 DO 110 ILLL = 1,NTOT 1075C 1076C------------------------------------------ 1077C If direct calculate the integrals. 1078C------------------------------------------ 1079C 1080 IF (DIRECT) THEN 1081C 1082 KEND1 = KENDSV 1083 LWRK1 = LWRKSV 1084C 1085 DTIME = SECOND() 1086 IF (HERDIR) THEN 1087 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 1088 & IPRINT) 1089 ELSE 1090 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 1091 & WORK(KODCL1),WORK(KODCL2), 1092 & WORK(KODBC1),WORK(KODBC2), 1093 & WORK(KRDBC1),WORK(KRDBC2), 1094 & WORK(KODPP1),WORK(KODPP2), 1095 & WORK(KRDPP1),WORK(KRDPP2), 1096 & WORK(KCCFB1),WORK(KINDXB), 1097 & WORK(KEND1), LWRK1,IPRERI) 1098 END IF 1099 DTIME = SECOND() - DTIME 1100 TIMER2 = TIMER2 + DTIME 1101C 1102 KRECNR = KEND1 1103 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 1104 LWRK1 = LWORK - KEND1 1105 IF (LWRK1 .LT. 0) THEN 1106 CALL QUIT('Insufficient core in CC_FMAT') 1107 END IF 1108C 1109 ELSE 1110 NUMDIS = 1 1111 ENDIF 1112C 1113C-------------------------------------------------- 1114C Loop over number of distributions in disk. 1115C-------------------------------------------------- 1116C 1117 DO 130 IDEL2 = 1,NUMDIS 1118C 1119 IF (DIRECT) THEN 1120 IDEL = INDEXA(IDEL2) 1121 IF (NOAUXB) THEN 1122 IDUM = 1 1123 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 1124 END IF 1125 ISYMD = ISAO(IDEL) 1126 ELSE 1127 IDEL = IBAS(ISYMD1) + ILLL 1128 ISYMD = ISYMD1 1129 ENDIF 1130C 1131 ISYDIS = MULD2H(ISYMD,ISYMOP) 1132 ISYAIK = MULD2H(ISYDIS,ISYMR) 1133C 1134C---------------------------------------------------------- 1135C Calculate adresses for c,cio,d,dio routines. 1136C---------------------------------------------------------- 1137C 1138 IT2DEL(IDEL) = ICDEL1 1139 ICDEL1 = ICDEL1 + NT2BCD(ISYDIS) 1140C 1141 IT2DLR(IDEL,1) = ICDEL2 1142 ICDEL2 = ICDEL2 + NT2BCD(ISYAIK) 1143C 1144C------------------------------------------ 1145C Work space allocation no. 2. 1146C------------------------------------------ 1147C 1148 KXINT = KEND1 1149 KEND2 = KXINT + NDISAO(ISYDIS) 1150 LWRK2 = LWORK - KEND2 1151C 1152 IF ( IPRINT .GT. 55) 1153 * WRITE(LUPRI,*) ' In CC_FMAT 2. alloc. work left:', 1154 * LWRK2 1155C 1156 IF (LWRK2 .LT. 0) THEN 1157 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 1158 CALL QUIT('Insufficient space in CC_FMAT-2') 1159 ENDIF 1160C 1161C----------------------------------------- 1162C Read in batch of integrals. 1163C----------------------------------------- 1164C 1165 DTIME = SECOND() 1166 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 1167 * WORK(KRECNR),DIRECT) 1168 DTIME = SECOND() - DTIME 1169 TIMRDAO = TIMRDAO + DTIME 1170C 1171C------------------------------------------------------------------ 1172C Calculate an AO-Fock matrix with C1. trans. density. 1173C------------------------------------------------------------------ 1174C 1175 ISYDEN = ISYMR 1176 DTIME = SECOND() 1177 IOPT = 1 1178 CALL CC_AOFOCK2(WORK(KXINT),WORK(KDENSC),WORK(KDNSPKC), 1179 * WORK(KFOCKC),WORK(KEND2),LWRK2, 1180 * IDEL,ISYDIS,ISYMD,ISYDEN,IOPT) 1181 DTIME = SECOND() - DTIME 1182 TIMFCK = TIMFCK + DTIME 1183C 1184C--------------------------------------------------------- 1185C Calculate AO-Fock matrix with L1R2 density. 1186C--------------------------------------------------------- 1187C 1188 ISYDEN = ISYRES 1189 DTIME = SECOND() 1190 IOPT = 1 1191 CALL CC_AOFOCK2(WORK(KXINT),WORK(KL1R2),WORK(KL1R2PK), 1192 * WORK(KFCKLR),WORK(KEND2),LWRK2, 1193 * IDEL,ISYDIS,ISYMD,ISYDEN,IOPT) 1194 DTIME = SECOND() - DTIME 1195 TIMFCK = TIMFCK + DTIME 1196C 1197C------------------------------------------------------------------- 1198C IF CCS calculate fock matrix since it is not on file. 1199C IF CCS jump to end of loop. 1200C------------------------------------------------------------------- 1201C 1202 IF ( CCS ) THEN 1203 ISYDEN = 1 1204 DTIME = SECOND() 1205 IOPT = 1 1206 CALL CC_AOFOCK2(WORK(KXINT),WORK(KDENSI), 1207 * WORK(KDNSPKI),WORK(KFOCK), 1208 * WORK(KEND2),LWRK2, 1209 * IDEL,ISYDIS,ISYMD,ISYDEN,IOPT) 1210 DTIME = SECOND() - DTIME 1211 TIMFCK = TIMFCK + DTIME 1212 ENDIF 1213C 1214C------------------------------- 1215C IF CCS then JUMP. 1216C------------------------------- 1217C 1218 IF (CCS) GOTO 130 1219C 1220C------------------------------------------------------- 1221C Work space allocation no. 3. 1222C------------------------------------------------------- 1223C 1224 ISAIJK = MULD2H(ISYMD,ISYMR) 1225 IAIJK2 = MULD2H(ISYMD,ISYMOP) 1226C 1227 KDSRHF = KEND2 1228 K3OINT = KDSRHF + NDSRHF(ISYMD) 1229 KO3T = K3OINT + NMAIJK(ISAIJK) 1230 KEND3 = KO3T + NMAIJK(IAIJK2) 1231C 1232 LWRK3 = LWORK - KEND3 1233 IF ( IPRINT .GT. 55) 1234 * WRITE(LUPRI,*) ' In CC_FMAT 3. alloc. work left:', 1235 * LWRK3 1236C 1237 IF (LWRK3 .LT. 0) THEN 1238 WRITE(LUPRI,*) 'Need : ',KEND3,'Available : ',LWORK 1239 CALL QUIT('Insufficient space in CC_FMAT-3 ') 1240 ENDIF 1241C 1242C-------------------------------------------------------- 1243C Transform one index in the integral batch. 1244C-------------------------------------------------------- 1245C 1246 DTIME = SECOND() 1247 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMDP), 1248 * ISYM0,WORK(KEND3),LWRK3,ISYDIS) 1249 TIMTRBT = TIMTRBT + SECOND() - DTIME 1250C 1251C------------------------------------------------------------------- 1252C calculate the BF intermediates: 1253C------------------------------------------------------------------- 1254C 1255 IF ( CCSD ) THEN 1256 ISYMGDR = MULD2H(ISYMD,ISYMR) 1257 NMGDR = NT2BGD(ISYMGDR) 1258C 1259 KMGDR = KEND3 1260 KBFRHF = KMGDR + NMGDR 1261 KEND4 = KBFRHF + NBFRHF(ISYMD) 1262 LWRK4 = LWORK - KEND4 1263C 1264 IF (LWRK4 .LT. 0) THEN 1265 CALL QUIT('Insufficient space in CC_FMAT-4 ') 1266 ENDIF 1267C 1268 IADR = IADRBFD(IDEL) 1269 CALL GETWA2(LUBFD,FNBFD,WORK(KMGDR),IADR,NMGDR) 1270C 1271 CALL CC_BFBSORT(WORK(KDSRHF),WORK(KBFRHF),ISYDIS) 1272C 1273 CALL CC_BFIB(WORK(KRHOR),WORK(KBFRHF),ISYDIS, 1274 * WORK(KMGDR),ISYMGDR,WORK(KEND4),LWRK4) 1275C 1276 CALL CC_BFIF(WORK(KBFRHF),ISYDIS,WORK(KMGDL),ISYRES, 1277 * LUBFI,FNBFI,IADRBFI,ISTARTBFI,IDEL, 1278 * WORK(KEND3),LWRK3) 1279 END IF 1280C 1281C------------------------------------------------------------------- 1282C Calculate integral batch with three occupied indices. 1283C------------------------------------------------------------------- 1284C 1285 DTIME = SECOND() 1286 CALL CC_INT3O(WORK(K3OINT),WORK(KDSRHF),WORK(KLAMHC), 1287 * ISYMR,WORK(KLAMDP),WORK(KEND3), 1288 * LWRK3,IDEL,ISYMD,LUO3,O3FIL) 1289 TIME3O = TIME3O + SECOND() - DTIME 1290C 1291C------------------------------------------------------------------- 1292C Calculate integral batch with three occupied indices. 1293C but C1 independent. 1294C------------------------------------------------------------------- 1295C 1296 IF (.NOT.(CC2.OR.CCS)) THEN 1297 DTIME = SECOND() 1298 CALL CC_INT3O(WORK(KO3T),WORK(KDSRHF),WORK(KLAMDH), 1299 * ISYMOP,WORK(KLAMDP),WORK(KEND3), 1300 * LWRK3,IDEL,ISYMD,LUO3T,O3TFIL) 1301 TIME3O = TIME3O + SECOND() - DTIME 1302 ENDIF 1303C 1304C-------------------------------------------------------------- 1305C Calculate intermediates needed for the 21-block. 1306C This section is for both perturbed an unperturbed. 1307C-------------------------------------------------------------- 1308C 1309 IF (.NOT. (CC2.OR.CCS)) THEN 1310C 1311C --------------------------------------------------- 1312C calculate first 21I and 21H terms with zeroth-order 1313C T2 amplitudes and first-order integrals (or LamdaP) 1314C --------------------------------------------------- 1315C 1316 ISYMTD = MULD2H(ISYMOP,ISYMD) 1317 ISYLTD = MULD2H(ISYMTD,ISYML) 1318C 1319 ISYMRD = MULD2H(ISYMR,ISYMD) 1320 ISYLRD = MULD2H(ISYMRD,ISYML) 1321C 1322 KSCTZI = KEND3 1323 KSCTVI = KSCTZI + NT2BCD(ISYLTD) 1324 KSCRZI = KSCTVI + NT2BCD(ISYLTD) 1325 KSCRVI = KSCRZI + NT2BCD(ISYLRD) 1326 KSCRWI = KSCRVI + NT2BCD(ISYLRD) 1327 KSCTWI = KSCRWI + NT2BCD(ISYLRD) 1328 KEND4 = KSCTWI + NT2BCD(ISYLTD) 1329 LWRK4 = LWORK - KEND4 1330 IF (LWRK4 .LT. 0) THEN 1331 CALL QUIT('Insufficient memory in CC_FMAT (4).') 1332 ENDIF 1333C 1334 CALL DZERO(WORK(KSCRWI),NT2BCD(ISYLRD)) 1335 CALL DZERO(WORK(KSCTWI),NT2BCD(ISYLTD)) 1336C 1337C --------------------------------------- 1338C read P0 and Q0 intermediates from file: 1339C --------------------------------------- 1340C 1341 DTIME = SECOND() 1342 CALL GETWA2(LUPQR0,FILPQR0,WORK(KSCTZI), 1343 & IADRPQ0(IDEL),NT2BCD(ISYLTD)) 1344 1345 CALL GETWA2(LUPQR0,FILPQR0,WORK(KSCTVI), 1346 & IADRPQ0(IDEL)+NT2BCD(ISYLTD),NT2BCD(ISYLTD)) 1347 TIMIO = TIMIO + SECOND() - DTIME 1348C 1349C --------------------------------------- 1350C read PR and QR intermediates from file: 1351C --------------------------------------- 1352C 1353 DTIME = SECOND() 1354 CALL GETWA2(LUPQRR,FILPQRR,WORK(KSCRZI), 1355 & IADRPQR(IDEL),NT2BCD(ISYLRD)) 1356 1357 CALL GETWA2(LUPQRR,FILPQRR,WORK(KSCRVI), 1358 & IADRPQR(IDEL)+NT2BCD(ISYLRD),NT2BCD(ISYLRD)) 1359 TIMIO = TIMIO + SECOND() - DTIME 1360C 1361C --------------------------------- 1362C Calculate the LT21I contribution: 1363C --------------------------------- 1364C 1365 IOPT = 2 1366 DTIME = SECOND() 1367 CALL CC_21I2(WORK(KRHO1AO), 1368 * WORK(KXINT),ISYDIS,DUMMY,0, 1369 * WORK(KSCTZI),WORK(KSCTVI),ISYLTD, 1370 * WORK(KSCRZI),WORK(KSCRVI),ISYLRD, 1371 * WORK(KLAMDP),WORK(KLAMDH),ISYMOP, 1372 * WORK(KLAMPC),ISYMR, 1373 * WORK(KEND4),LWRK4,IOPT, 1374 * .TRUE.,.FALSE.,.FALSE.) 1375 TIMEI = TIMEI + SECOND() - DTIME 1376C 1377C --------------------------------- 1378C Calculate the LT21H contribution: 1379C --------------------------------- 1380C 1381 ISYVWZ = ISYLTD 1382 ISYINT = ISYMR 1383 DTIME = SECOND() 1384 CALL CC_21H(WORK(KRHO1),ISYRES,WORK(KSCTVI), 1385 * WORK(KSCTWI),WORK(KSCTZI),ISYVWZ, 1386 * WORK(K3OINT),ISYINT, 1387 * WORK(KEND4),LWRK4,ISYMD,NEWZWV) 1388 TIMEH = TIMEH + SECOND() - DTIME 1389 WRITE (LUPRI,*) 'Norm(RHO1) = ', 1390 * DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1391C 1392C ---------------------------------------------------- 1393C Calculate the LT21H contribution with first-order 1394C T2 amplitudes and zeroth-order integrals 1395C ---------------------------------------------------- 1396C 1397 ISYVWZ = ISYLRD 1398 ISYINT = ISYM0 1399 DTIME = SECOND() 1400 CALL CC_21H(WORK(KRHO1),ISYRES,WORK(KSCRVI), 1401 * WORK(KSCRWI),WORK(KSCRZI),ISYVWZ, 1402 * WORK(KO3T),ISYINT, 1403 * WORK(KEND4),LWRK4,ISYMD,NEWZWV) 1404 TIMEH = TIMEH + SECOND() - DTIME 1405 WRITE (LUPRI,*) 'Norm(RHO1) = ', 1406 * DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1407 1408 ENDIF 1409C 1410C---------------------------------------------------------- 1411C Calculate the LT12B term in CC2-calculation. 1412C rho(ai,bj)=P(ai,bj)(sum(c)(l(c,j)*~cbia)) 1413C---------------------------------------------------------- 1414C 1415 IF (CC2) THEN 1416 1417 DTIME = SECOND() 1418 CALL CC_12B(WORK(KRHO2),WORK(KDSRHF),WORK(KL1AOC), 1419 * ISYRES,WORK(KLAMDH),WORK(KEND3),LWRK3, 1420 * IDEL,ISYMD) 1421 TIM12B = TIM12B + SECOND() - DTIME 1422 1423 ENDIF 1424C 1425C-------------------------------------------------------------------- 1426C Construct the partially transformed L2-amplitudes 1427C-------------------------------------------------------------------- 1428C 1429 ISYDVI = MULD2H(ISYMD,ISYMR) 1430 KSCRLM = KEND3 1431 KSCRLC = KSCRLM + NT2MMO(ISYMD,ISYML) 1432 KEND41 = KSCRLC + NT2MMO(ISYDVI,ISYML) 1433 LWRK41 = LWORK - KEND41 1434C 1435 DTIME = SECOND() 1436 ICON = 1 1437 ISYMLP = 1 1438 CALL CC_T2AO(WORK(KL2AM),WORK(KLAMDP),ISYMLP, 1439 * WORK(KSCRLM),WORK(KEND41),LWRK41, 1440 * IDEL,ISYMD,ISYML,ICON) 1441 TIMT2AO = TIMT2AO + SECOND() - DTIME 1442C 1443 DTIME = SECOND() 1444 ICON = 1 1445 ISYMLP = ISYMR 1446 CALL CC_T2AO(WORK(KL2AM),WORK(KLAMPC),ISYMLP, 1447 * WORK(KSCRLC),WORK(KEND41),LWRK41, 1448 * IDEL,ISYMD,ISYML,ICON) 1449 TIMT2AO = TIMT2AO + SECOND() - DTIME 1450C 1451C---------------------------------------------------------------- 1452C Calculate the 21C- and D contributions. 1453C If not CC2 skip 21C calculation and calculate this 1454C as part of the 21B terms. 1455C---------------------------------------------------------------- 1456C 1457 DTIME = SECOND() 1458 IOPT = 2 1459 IF ( CC2 ) THEN 1460 ICON = 2 1461 ELSE 1462 ICON = 1 1463 ENDIF 1464 1465 ISYMM = MULD2H(ISYMD,ISYML) 1466 ISYMMC = MULD2H(ISYDVI,ISYML) 1467 CALL CC_21DC(WORK(KRHO1),WORK(KL2AM),ISYML, 1468 * WORK(KSCRLM),ISYMM, 1469 * WORK(KSCRLC),ISYMMC,WORK(KXINT), 1470 * WORK(KLAMDH),1,WORK(KLAMPC),ISYMR, 1471 * WORK(KLAMHC),ISYMR,WORK(KLAMDP),1, 1472 * WORK(KEND41),LWRK41,IDEL,ISYMD,IOPT,ICON) 1473 TIMEC = TIMEC + SECOND() - DTIME 1474C 1475C-------------------------------------------------------------------- 1476C Construct the partially transformed C2-amplitudes, CM. 1477C-------------------------------------------------------------------- 1478C 1479 KSCRCM = KEND41 1480 KEND4 = KSCRCM + NT2MMO(ISYMD,ISYMR) 1481 LWRK4 = LWORK - KEND4 1482C 1483 DTIME = SECOND() 1484 ICON = 1 1485 ISYMLH = 1 1486 CALL CC_T2AO(WORK(KC2AM),WORK(KLAMDH),ISYMLH, 1487 * WORK(KSCRCM),WORK(KEND4),LWRK4, 1488 * IDEL,ISYMD,ISYMR,ICON) 1489 TIMT2AO = TIMT2AO + SECOND() - DTIME 1490C 1491C--------------------------------------- 1492C Transform CM to 2CM - CM. 1493C--------------------------------------- 1494C 1495 DTIME = SECOND() 1496 CALL CC_MTCME(WORK(KSCRCM),WORK(KEND4),LWRK4, 1497 * ISYMD,ISYMR) 1498 TIMTCME = TIMTCME + SECOND() - DTIME 1499C 1500C------------------------------------------------------- 1501C Calculate the C-tilde local intermediate. 1502C------------------------------------------------------- 1503C 1504 IF ( .NOT. (CC2.OR.CCS)) THEN 1505C 1506 FACTC = XMONE 1507 ICON = 3 1508 IOPTR12 = 0 1509 IOPTE = 0 1510C 1511 DTIME = SECOND() 1512 IF (.NOT. T2TCOR) THEN 1513 CALL CCRHS_C(WORK(KXINT),WORK(KDSRHF), 1514 * DUMMY, WORK(KC2AM),ISYMR, 1515 * WORK(KLAMDP),WORK(KLAMDP), 1516 * WORK(KLAMDH),WORK(KLAMPC),ISYMR, 1517 * WORK(KLAMHC),ISYMR, 1518 * DUMMY,DUMMY,WORK(KEND4), 1519 * LWRK4,IDEL,ISYMD,FACTC,ICON,IOPTR12, 1520 * IOPTE,LUC,CTFIL,IDUMMY,DUMMY,1) 1521 ELSE 1522 CALL CCRHS_C(WORK(KXINT),WORK(KDSRHF), 1523 * DUMMY, WORK(K2C2C2),ISYMR, 1524 * WORK(KLAMDP),WORK(KLAMDP), 1525 * WORK(KLAMDH),WORK(KLAMPC),ISYMR, 1526 * WORK(KLAMHC),ISYMR, 1527 * DUMMY,DUMMY,WORK(KEND4), 1528 * LWRK4,IDEL,ISYMD,FACTC,ICON,IOPTR12, 1529 * IOPTE,LUC,CTFIL,IDUMMY,DUMMY,1) 1530 ENDIF 1531 TIMC = TIMC + SECOND() - DTIME 1532 1533 ENDIF 1534C 1535C--------------------------------------- 1536C Transform C2 to 2C2 - C2. 1537C--------------------------------------- 1538C 1539 DTIME = SECOND() 1540 IF (T2TCOR) THEN 1541 CALL DSCAL(NT2SQ(ISYMR),TWO,WORK(KC2AM),1) 1542 CALL DAXPY(NT2SQ(ISYMR),-ONE,WORK(K2C2C2),1, 1543 * WORK(KC2AM),1) 1544 ELSE 1545 CALL CCRHS_T2TR(WORK(KC2AM),WORK(KEND4),LWRK4,ISYMR) 1546 ENDIF 1547 TIMT2TR = TIMT2TR + SECOND() - DTIME 1548C 1549C------------------------------------------------------- 1550C Calculate the D-tilde local intermediate. 1551C------------------------------------------------------- 1552C 1553 IF ( .NOT. (CC2.OR.CCS)) THEN 1554C 1555 ICON = 3 1556 FACTD = 1.0D00 1557 IOPTR12 = 0 1558 IOPTE = 0 1559C 1560 DTIME = SECOND() 1561 CALL CCRHS_D(WORK(KXINT),WORK(KDSRHF),DUMMY, 1562 * WORK(KC2AM),ISYMR, 1563 * WORK(KLAMDP),DUMMY,WORK(KLAMDH), 1564 * WORK(KLAMPC),ISYMR,WORK(KLAMHC),ISYMR, 1565 * DUMMY,DUMMY,WORK(KEND4),LWRK4,IDEL, 1566 * ISYMD,FACTD,ICON,IOPTR12,IOPTE, 1567 * LUD,DTFIL,IDUMMY,DUMMY,1) 1568 TIMD = TIMD + SECOND() - DTIME 1569C 1570 ENDIF 1571C 1572C---------------------------------------- 1573C Calculate E-intermediates. 1574C---------------------------------------- 1575C 1576 IF (.NOT.CCS) THEN 1577 DTIME = SECOND() 1578 CALL CCRHS_EI(WORK(KDSRHF),WORK(KEMAT1), 1579 * WORK(KEMAT2),WORK(KC2AM), 1580 * WORK(KSCRCM),WORK(KLAMDP), 1581 * WORK(KLAMDH),WORK(KEND4),LWRK4, 1582 * IDEL,ISYMD,ISYDIS,ISYMR) 1583 ENDIF 1584 TIMEI = TIMEI + SECOND() - DTIME 1585C 1586C--------------------------------------------- 1587C BackTransform C2 from 2C2 - C2. 1588C--------------------------------------------- 1589C 1590 DTIME = SECOND() 1591 IF (T2TCOR) THEN 1592 CALL DAXPY(NT2SQ(ISYMR),ONE,WORK(K2C2C2),1, 1593 * WORK(KC2AM),1) 1594 CALL DSCAL(NT2SQ(ISYMR),XHALF,WORK(KC2AM),1) 1595 ELSE 1596 CALL CCRHS_T2BT(WORK(KC2AM),WORK(KEND4),LWRK4,ISYMR) 1597 END IF 1598 TIMT2BT = TIMT2BT + SECOND() - DTIME 1599C 1600 IF (IPRINT .GT. 120) THEN 1601 CALL AROUND('END of LOOP CCLR:RHO ') 1602 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,1,1) 1603 ENDIF 1604C 1605 IF (IPRINT .GT. 20 .OR. DEBUG) THEN 1606 RHO1N=DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1607 IF (OMEGOR) THEN 1608 RHO2N = DDOT(2*NT2ORT(ISYRES),WORK(KRHO2),1, 1609 * WORK(KRHO2),1) 1610 ELSE 1611 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1, 1612 * WORK(KRHO2),1) 1613 ENDIF 1614 WRITE(LUPRI,1)'Norm of Rho1 - end of loop: ', 1615 & RHO1N 1616 WRITE(LUPRI,1)'Norm of Rho2 - end of loop: ', 1617 & RHO2N 1618 ENDIF 1619C 1620 130 CONTINUE 1621C 1622C-------------------------------------- 1623C write out result vector. 1624C-------------------------------------- 1625C 1626 IF (.NOT. CCS) THEN 1627 DTIME = SECOND() 1628 CALL CC_WVEC(LUTR,TRFIL,NRHO2,NRHO2,1,WORK(KRHO2)) 1629 DTIME = SECOND() - DTIME 1630 TIMIO = TIMIO + DTIME 1631 ENDIF 1632C 1633 110 CONTINUE 1634 100 CONTINUE 1635C 1636C------------------------ 1637C Recover work space. 1638C------------------------ 1639C 1640 KT2AMP = KENDS2 1641 KENDS2 = KT2AMP + MAX(NT2AM(ISYMOP),2*NT2ORT(ISYMOP)) 1642 LWRKS2 = LWORK - KENDS2 1643 1644 KEND1 = KENDS2 1645 LWRK1 = LWRKS2 1646C 1647C-------------------------------------------- 1648C Allocate space for the gamma matrix. 1649C-------------------------------------------- 1650C 1651 IF (NEWGAM) THEN 1652C 1653 KGAMMC = KEND1 1654 KEND1 = KGAMMC + NGAMMA(ISYMR) 1655 LWRK1 = LWORK - KEND1 1656C 1657 IF (LWRK1 .LT. 0) CALL QUIT('Insufficient memory in CC_FMAT') 1658C 1659 END IF 1660C 1661C---------------------------------- 1662C Usual Fock Matrix. 1663C Save AO fock matric in fockt. 1664C---------------------------------- 1665C 1666 ISYFAO = 1 1667 ISYMPA = 1 1668 ISYMHO = 1 1669C 1670 IF ( .NOT. CCS) THEN 1671 LFOCK = -1 1672 CALL GPOPEN(LFOCK,'CC_FCKH','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1673 * .FALSE.) 1674 REWIND(LFOCK) 1675 READ(LFOCK) (WORK(KFOCK + I - 1) , I = 1,N2BST(ISYMOP)) 1676 CALL GPCLOSE(LFOCK,'KEEP') 1677 ENDIF 1678C 1679 IF (IPRINT .GT.140) THEN 1680 CALL AROUND( 'Usual Fock AO matrix' ) 1681 CALL CC_PRFCKAO(WORK(KFOCK),ISYFAO) 1682 ENDIF 1683C 1684 CALL DCOPY(N2BST(ISYMOP),WORK(KFOCK),1,WORK(KFOCKT),1) 1685C 1686 DTIME = SECOND() 1687 CALL CC_FCKMO(WORK(KFOCK),WORK(KLAMDP),WORK(KLAMDH), 1688 * WORK(KEND1),LWRK1,ISYFAO,ISYMPA,ISYMHO) 1689 DTIME = SECOND() - DTIME 1690 TIMFCKMO = DTIME + TIMFCKMO 1691C 1692 IF (IPRINT .GT. 50) THEN 1693 CALL AROUND( 'Usual Fock MO matrix' ) 1694 CALL CC_PRFCKMO(WORK(KFOCK),ISYFAO) 1695 ENDIF 1696C 1697C----------------------------------------- 1698C Calculate the LT11A contribution. 1699C----------------------------------------- 1700C 1701 DTIME = SECOND() 1702 IF (DEBUG) THEN 1703 RHO1N = DDOT(N2BST(ISYRES),WORK(KFCKLR),1,WORK(KFCKLR),1) 1704 WRITE(LUPRI,1) 'Norm of FLR1 before CC_11A: ',RHO1N 1705 ENDIF 1706 1707 CALL CC_11A(WORK(KRHO1),WORK(KFCKLR),ISYRES,WORK(KLAMDH), 1708 * WORK(KLAMDP),WORK(KEND1),LWRK1) 1709 1710 IF (DEBUG) THEN 1711 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1712 WRITE(LUPRI,1) 'Norm of Rho1 -after CC_11A: ',RHO1N 1713 ENDIF 1714C 1715C-------------------------------------------------------------- 1716C transform CC_21I contribution to MO and add to RHO1: 1717C-------------------------------------------------------------- 1718C 1719 CALL CC_T1AM(WORK(KRHO1),ISYRES,WORK(KRHO1AO),ISYRES, 1720 * WORK(KLAMDH),ISYMOP,ONE) 1721C 1722C=============================================================== 1723C After Loop CCSD and CC2 Contributions with C2SQ in memory. 1724C=============================================================== 1725C 1726 IF ( .NOT. (CC2.OR.CCS)) THEN 1727C 1728C-------------------------------------------------------------- 1729C Transform the Omega2 vector to the MO basis. 1730C We thus have the B(L2) term. 1731C-------------------------------------------------------------- 1732C 1733 ICON = 1 1734 IOPTG = 0 1735 ISYMHC = 1 1736 ISYMBF = ISYRES 1737C 1738 LGAMMA = .FALSE. 1739 LGIM = .FALSE. 1740 LO3BF = .FALSE. 1741 LBFZETA = .TRUE. 1742 DTIME = SECOND() 1743 1744 CALL DZERO(WORK(KC2AM),NT2AM(ISYRES)) 1745 1746 CALL CC_BFIFSORT(WORK(KRHO2),ISYRES,LUBFI,FNBFI,IADRBFI, 1747 * WORK(KEND1),LWRK1) 1748 1749 CALL CC_T2MO3(DUMMY,DUMMY,ISYMOP,WORK(KRHO2), 1750 * WORK(KC2AM),DUMMY,DUMMY,DUMMY, 1751 * WORK(KLAMDH),ISYM0,WORK(KLAMDH),ISYM0, 1752 * WORK(KEND1),LWRK1,ISYMBF, 1753 * ICON,LGAMMA,IOPTG,LO3BF,LBFZETA) 1754 DTIME = SECOND() -DTIME 1755 TIMT2MO = TIMT2MO + DTIME 1756 NEWGAM = .TRUE. 1757C 1758 CALL DCOPY(NT2AM(ISYRES),WORK(KC2AM),1,WORK(KRHO2),1) 1759C 1760 IF (DEBUG) THEN 1761 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1762 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 1763 WRITE(LUPRI,1) 'Norm of Rho1 -after loop in MO: ',RHO1N 1764 WRITE(LUPRI,1) 'Norm of Rho2 -after loop in MO: ',RHO2N 1765 ENDIF 1766C 1767 IF (IPRINT .GT. 120) THEN 1768 CALL AROUND('After T2MO-1:BF(C1,C2) ') 1769 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,1,1) 1770 ENDIF 1771C 1772C-------------------------------------- 1773C Calculate the LT21G contribution. 1774C-------------------------------------- 1775C 1776 IF (DEBUG) THEN 1777 xn = DDOT(N3ORHF(ISYRES),WORK(KMINTX),1,WORK(KMINTX),1) 1778 WRITE(LUPRI,1) 'Norm of MX -before CC_21G: ',xn 1779 xn = DDOT(N3ORHF(ISYML),WORK(KMINT),1,WORK(KMINT),1) 1780 WRITE(LUPRI,1) 'Norm of M -before CC_21G: ',xn 1781 ENDIF 1782C 1783 TIMEG = SECOND() 1784 CALL CC_21G(WORK(KRHO1),WORK(KMINTX),ISYRES,WORK(KLAMDH), 1785 * WORK(KEND1),LWRK1,ISYMOP,LUO3T,O3TFIL) 1786 CALL CC_21G(WORK(KRHO1),WORK(KMINT),ISYML,WORK(KLAMDH), 1787 * WORK(KEND1),LWRK1,ISYMR,LUO3,O3FIL) 1788 TIMEG = SECOND() - TIMEG 1789C 1790 IF (DEBUG) THEN 1791 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1792 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 1793 WRITE(LUPRI,1) 'Norm of Rho1 -after CC_21G: ',RHO1N 1794 WRITE(LUPRI,1) 'Norm of Rho2 -after CC_21G: ',RHO2N 1795 ENDIF 1796C 1797C-------------------------------------------------------------- 1798C Transform the Omega2 vector to the MO basis. 1799C We thus have the B(C2) and (a,i-bar|bj) contributions. 1800C-------------------------------------------------------------- 1801C 1802 IF (DEBUG) THEN 1803 RHO2N = DDOT(NRHOR,WORK(KRHOR),1,WORK(KRHOR),1) 1804 WRITE(LUPRI,1) 'Norm of BF(C1,C2) intermediate: ',RHO2N 1805 RHO2N = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 1806 WRITE(LUPRI,1) 'Norm of L2: ',RHO2N 1807 ENDIF 1808C 1809 CALL DZERO(WORK(KGAMMC),NGAMMA(ISYMR)) 1810C 1811 ICON = 3 1812 IOPTG = 0 1813 LGAMMA = .TRUE. 1814 LGIM = .FALSE. 1815 LO3BF = .TRUE. 1816 LBFZETA= .FALSE. 1817 ISYMPC = 1 1818 ISYMBF = ISYMR 1819C 1820 DTIME = SECOND() 1821 CALL CC_T2MO3(WORK(KRHO1),WORK(KL2AM),ISYML, 1822 * WORK(KRHOR),DUMMY,WORK(KGAMMC),DUMMY,DUMMY, 1823 * WORK(KLAMDP),ISYM0,WORK(KLAMDP),ISYMPC, 1824 * WORK(KEND1),LWRK1,ISYMBF, 1825 * ICON,LGAMMA,IOPTG,LO3BF,LBFZETA) 1826 DTIME = SECOND() -DTIME 1827 TIMT2MO = TIMT2MO + DTIME 1828C 1829 IF ( DEBUG ) THEN 1830 XNGAM = DDOT(NGAMMA(ISYMR),WORK(KGAMMC),1,WORK(KGAMMC),1) 1831 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1832 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 1833 WRITE(LUPRI,1) 'Norm of GamC -after T2MO: ',XNGAM 1834 WRITE(LUPRI,1) 'Norm of Rho1 -after L2*BF(C1,C1): ',RHO1N 1835 WRITE(LUPRI,1) 'Norm of Rho2 -after L2*BF(C1,C1): ',RHO2N 1836 ENDIF 1837C 1838C---------------------------------- 1839C Read in integrals (ia|jb). 1840C---------------------------------- 1841C 1842 REWIND(LUIAJB) 1843 READ(LUIAJB) (WORK(KC2AM-1+J), J = 1,NT2AM(ISYMOP)) 1844C 1845C----------------------------------------------------------------- 1846C Calculate 22EM contributions: 1847C calculate 2(ia|jb) - (ja|ib) combination of the integrals 1848C and calculate coulomb and exchange contributions together 1849C----------------------------------------------------------------- 1850C 1851 DTIME = SECOND() 1852 1853 KXIM = KEND1 1854 KYIM = KXIM + NMATIJ(ISYRES) 1855 KEND1 = KYIM + NMATAB(ISYRES) 1856 LWRK1 = LWORK - KEND1 1857 1858 CALL DZERO(WORK(KXIM),NMATIJ(ISYRES)) 1859 CALL DCOPY(NMATAB(ISYRES),WORK(KYMATX),1,WORK(KYIM),1) 1860 CALL DSCAL(NMATAB(ISYRES),XHALF,WORK(KYIM),1) 1861 IOPTTCME = 1 1862 CALL CCSD_TCMEPK(WORK(KC2AM),ONE,ISYMOP,IOPTTCME) 1863 1864 CALL CC_22EC(WORK(KRHO2),WORK(KC2AM),WORK(KXIM), 1865 * WORK(KYIM),ISYRES,WORK(KEND1),LWRK1) 1866 1867 KEND1 = KXIM 1868 LWRK1 = LWORK - KEND1 1869C 1870 TIM2EM = SECOND() - DTIME 1871 1872 IF (DEBUG) THEN 1873 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1874 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 1875 WRITE(LUPRI,1) 'Norm of Rho1 -after Int*X + Int*Y: ',RHO1N 1876 WRITE(LUPRI,1) 'Norm of Rho2 -after Int*X + Int*Y: ',RHO2N 1877 ENDIF 1878C 1879 IF (IPRINT .GT. 120) THEN 1880 CALL AROUND('After Int*X and Int*Y E-terms) ') 1881 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,1,1) 1882 ENDIF 1883C 1884C-------------------------------- 1885C write out result vector. 1886C-------------------------------- 1887C 1888 DTIME = SECOND() 1889 CALL CC_WVEC(LUTR,TRFIL,NRHO2,NT2AM(ISYRES),1,WORK(KRHO2)) 1890 DTIME = SECOND() - DTIME 1891 TIMIO = TIMIO + DTIME 1892C 1893 ENDIF 1894C 1895 IF (.NOT.(CCS.OR.CC2)) THEN 1896C 1897C-------------------------------- 1898C Read Omega intermediate. 1899C-------------------------------- 1900C 1901 TIMEBF = SECOND() 1902 LUOM = -1 1903 CALL GPOPEN(LUOM,'CC_BFIM','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1904 & .FALSE.) 1905 REWIND(LUOM) 1906 READ(LUOM) (WORK(KT2AMP+I-1),I = 1,2*NT2ORT(1)) 1907 CALL GPCLOSE(LUOM,'KEEP') 1908C 1909C------------------------------------------ 1910C Calculate the LT21BF contribution. 1911C------------------------------------------ 1912C 1913 ICON = 3 1914C 1915 IF (DEBUG) THEN 1916 RHO2N=DDOT(2*NT2ORT(ISYMOP),WORK(KT2AMP),1,WORK(KT2AMP),1) 1917 WRITE(LUPRI,1) 'Norm of BF intermediate from disk: ',RHO2N 1918 RHO2N = DDOT(NGLMDT(ISYMR),WORK(KLAMPC),1,WORK(KLAMPC),1) 1919 WRITE(LUPRI,1) 'Norm of LMPC1: ',RHO2N 1920 ENDIF 1921C 1922 NEWGAM = .FALSE. 1923 CALL CC_T2MO(WORK(KRHO1),WORK(KL2AM),ISYML, 1924 * WORK(KT2AMP),PHONEY,DUMMY, 1925 * WORK(KLAMDP),WORK(KLAMPC),ISYMR,WORK(KEND1), 1926 * LWRK1,ISYMOP,ICON) 1927 NEWGAM = .TRUE. 1928 TIMEBF = SECOND() - TIMEBF 1929C 1930 IF (DEBUG) THEN 1931 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1932 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 1933 WRITE(LUPRI,1) 'Norm of Rho1 -after CC_T2MO 21BF: ',RHO1N 1934 WRITE(LUPRI,1) 'Norm of Rho2 -after CC_T2MO 21BF: ',RHO2N 1935 ENDIF 1936C 1937C----------------------------------- 1938C Readin C2 amplitude in RHO. 1939C----------------------------------- 1940C 1941 DTIME = SECOND() 1942 IOPT = 2 1943 CALL CC_RDRSP(RLIST,IRLNR,ISYMR,IOPT,MODELX, 1944 * DUMMY,WORK(KRHO2)) 1945 DTIME = SECOND() - DTIME 1946 TIMIO = TIMIO + DTIME 1947C 1948C-------------------------------- 1949C Square up C2 amplitudes. 1950C-------------------------------- 1951C 1952 DTIME = SECOND() 1953 CALL CCLR_DIASCL(WORK(KRHO2),TWO,ISYMR) 1954 CALL CC_T2SQ(WORK(KRHO2),WORK(KC2AM),ISYMR) 1955 DTIME = SECOND() - DTIME 1956 TIMT2SQ = TIMT2SQ + DTIME 1957C 1958 IF (IPRINT.GT.50) THEN 1959 CALL AROUND('CC_FMAT: (C1,C2) vector readin') 1960 CALL CC_PRSQ(WORK(KC1AM),WORK(KC2AM),ISYMR,1,1) 1961 ENDIF 1962C 1963C------------------------------ 1964C Read in result vector. 1965C------------------------------ 1966C 1967 DTIME = SECOND() 1968 CALL CC_RVEC(LUTR,TRFIL,NRHO2,NT2AM(ISYRES),1,WORK(KRHO2)) 1969 DTIME = SECOND() - DTIME 1970 TIMIO = TIMIO + DTIME 1971C 1972 ENDIF 1973C 1974 IF (IPRINT .GT. 15) THEN 1975 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 1976 WRITE(LUPRI,1) 'Norm of Rho1 -after loop in mo: ',RHO1N 1977 IF (.NOT. CCS) THEN 1978 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 1979 WRITE(LUPRI,1) 'Norm of Rho2 -after loop in mo: ',RHO2N 1980 ENDIF 1981 ENDIF 1982C 1983C------------------------------------------ 1984C Transform AO Fock matrix to MO basis. 1985C------------------------------------------ 1986C 1987C 1988C------------------------------------- 1989C C1 transformed Fock Matrix. 1990C first make 1991C F-dot: F`ai = SUM-k-(Lai,kk-bar). 1992C------------------------------------- 1993C 1994 ISYFCK = ISYMR 1995 ISYMPA = 1 1996 ISYMHO = 1 1997C 1998 IF (IPRINT .GT.140) THEN 1999 CALL AROUND( 'Fock AO matrix calc. from C1 transf. dens.' ) 2000 CALL CC_PRFCKAO(WORK(KFOCKC),ISYFCK) 2001 ENDIF 2002C 2003 DTIME = SECOND() 2004 CALL CC_FCKMO(WORK(KFOCKC),WORK(KLAMDP),WORK(KLAMDH), 2005 * WORK(KEND1),LWRK1,ISYFCK,ISYMPA,ISYMHO) 2006 DTIME = SECOND() - DTIME 2007 TIMFCKMO = DTIME + TIMFCKMO 2008C 2009 IF (IPRINT .GT. 50) THEN 2010 CALL AROUND( 'Fock MO matrix calc. from C1 transf. dens.' ) 2011 CALL CC_PRFCKMO(WORK(KFOCKC),ISYFCK) 2012 ENDIF 2013C 2014C----------------------------------------------------- 2015C Make Fock tilde: 2016C F~ai = SUM-k-(Lai,kk-bar) + Fa-bar,i + Fa,i-bar. 2017C----------------------------------------------------- 2018C 2019 KFCKAO = KEND1 2020 KEND2 = KFCKAO + MAX(N2BST(ISYMOP),N2BST(ISYMR)) 2021 LEND2 = LWORK - KEND2 2022C 2023 IF (IPRINT .GT.140) THEN 2024 CALL AROUND( 'Fock tilde -1.' ) 2025 CALL CC_PRFCKMO(WORK(KFOCKC),ISYMR) 2026 ENDIF 2027C 2028 ISYFAO = 1 2029 ISYMPA = ISYMR 2030 ISYMHO = 1 2031C 2032 CALL DCOPY(N2BST(ISYMOP),WORK(KFOCKT),1,WORK(KFCKAO),1) 2033 DTIME = SECOND() 2034 CALL CC_FCKMO(WORK(KFCKAO),WORK(KLAMPC),WORK(KLAMDH), 2035 * WORK(KEND2),LWRK2,ISYFAO,ISYMPA,ISYMHO) 2036 DTIME = SECOND() - DTIME 2037 TIMFCKMO = DTIME + TIMFCKMO 2038C 2039 CALL DAXPY(N2BST(ISYMR),ONE,WORK(KFCKAO),1,WORK(KFOCKC),1) 2040C 2041 IF (IPRINT .GT.140) THEN 2042 CALL AROUND( 'Fock tilde - 2.' ) 2043 CALL CC_PRFCKMO(WORK(KFOCKC),ISYMR) 2044 ENDIF 2045C 2046 ISYFAO = 1 2047 ISYMPA = 1 2048 ISYMHO = ISYMR 2049C 2050 CALL DCOPY(N2BST(ISYMOP),WORK(KFOCKT),1,WORK(KFCKAO),1) 2051 DTIME = SECOND() 2052 CALL CC_FCKMO(WORK(KFCKAO),WORK(KLAMDP),WORK(KLAMHC), 2053 * WORK(KEND2),LWRK2,ISYFAO,ISYMPA,ISYMHO) 2054 DTIME = SECOND() - DTIME 2055 TIMFCKMO = DTIME + TIMFCKMO 2056C 2057 CALL DAXPY(N2BST(ISYMR),ONE,WORK(KFCKAO),1,WORK(KFOCKC),1) 2058C 2059 IF (DEBUG) THEN 2060 XE1 = DDOT(N2BST(ISYMR),WORK(KFOCKC),1,WORK(KFOCKC),1) 2061 WRITE(LUPRI,1) 'Norm of FCKC1 after construction: ',XE1 2062 ENDIF 2063 IF (IPRINT .GT. 50) THEN 2064 CALL AROUND( 'Fock-Tilde MO matrix ' ) 2065 CALL CC_PRFCKMO(WORK(KFOCKC),ISYMR) 2066 ENDIF 2067C 2068C------------------------------------------------------------- 2069C Calculate simple fock contributions: 2070C rho2(ai,bj) = Paibj(2L1am(ai)FCKMO(jb)-L1am(aj)*FCKMO(ib) 2071C------------------------------------------------------------- 2072C 2073 IF (.NOT. CCS) THEN 2074 IF (IPRINT .GT. 15) THEN 2075 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2076 WRITE(LUPRI,1) 'Norm of Rho1 before L1FCK: ',RHO1N 2077 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2078 WRITE(LUPRI,1) 'Norm of Rho2 before L1FCK: ',RHO2N 2079 ENDIF 2080C 2081 CALL CC_L1FCK(WORK(KRHO2),WORK(KL1AM),WORK(KFOCKC),ISYML,ISYMR, 2082 * WORK(KEND2),LWRK2) 2083C 2084 IF (DEBUG) THEN 2085 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2086 WRITE(LUPRI,1) 'Norm of Rho1 -after L1FCK: ',RHO1N 2087 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2088 WRITE(LUPRI,1) 'Norm of Rho2 -after L1FCK: ',RHO2N 2089 WRITE(LUPRI,*) ' ' 2090 ENDIF 2091 ENDIF 2092C 2093C---------------------------------------------------- 2094C Calculate the 12C contribution: 2095C rho2(ai,bj) = P(ai.bj)(-sum(k)CTR1(a,k)*L(jbik) 2096C If .not.cc2 LT21B contributions to rho1. 2097C---------------------------------------------------- 2098C 2099 TIME12C = SECOND() 2100 IF (.NOT. CCS) THEN 2101 IOPT = 2 2102 CALL CC_21B12C(WORK(KRHO1),WORK(KRHO2),WORK(KL1AM),ISYML, 2103 * WORK(KLAMDH),ISYMOP,WORK(KXMAT),ISYML,ISYMR, 2104 * WORK(KEND2),LWRK2,LUO3,O3FIL,IOPT) 2105 ENDIF 2106 IF (.NOT.(CC2.OR.CCS)) THEN 2107 IOPT = 1 2108 CALL CC_21B12C(WORK(KRHO1),WORK(KRHO2),WORK(KL1AM),ISYML, 2109 * WORK(KLAMDH),ISYMOP,WORK(KXMATX),ISYRES,ISYMOP, 2110 * WORK(KEND2),LWRK2,LUO3T,O3TFIL,IOPT) 2111 ENDIF 2112 TIME12C = SECOND() - TIME12C 2113C 2114 IF (DEBUG .AND. (.NOT. CCS)) THEN 2115 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2116 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2117 WRITE(LUPRI,1) 'Norm of Rho1 -after CC_12C: ',RHO1N 2118 WRITE(LUPRI,1) 'Norm of Rho2 -after CC_12C: ',RHO2N 2119 ENDIF 2120C 2121C---------------------------------------- 2122C Calculate the LT21EFM contribution. 2123C---------------------------------------- 2124C 2125 IF ((NFIELD.GT.0).AND.(CC2.AND.NONHF)) THEN 2126 ISYMXY = MULD2H(ISYML,ISYMR) 2127 IOPT = 1 2128 CALL CC_LCC2FF(WORK(KRHO1),WORK(KLAMDP),WORK(KLAMDH), 2129 * WORK(KLAMPC),WORK(KLAMHC),ISYMOP, 2130 * WORK(KXMATX),WORK(KYMATX),ISYMXY, 2131 * WORK(KEND1),LWRK1,IOPT) 2132 ENDIF 2133 IF (.NOT.(CC2.OR.CCS)) THEN 2134C 2135 DTIME = SECOND() 2136 ISYFCK = ISYMR 2137 ISYMXY = ISYML 2138 CALL CC_21EFM(WORK(KRHO1),WORK(KFOCKC),ISYFCK,WORK(KXMAT), 2139 * WORK(KYMAT),ISYMXY, 2140 * WORK(KEND1),LWRK1) 2141 DTIME = SECOND() - DTIME 2142 TIMEEM = DTIME 2143 DTIME = SECOND() 2144C 2145 ISYFCK = ISYMOP 2146 ISYMXY = MULD2H(ISYML,ISYMR) 2147 CALL CC_21EFM(WORK(KRHO1),WORK(KFOCK),ISYFCK,WORK(KXMATX), 2148 * WORK(KYMATX),ISYMXY, 2149 * WORK(KEND1),LWRK1) 2150 DTIME = SECOND() - DTIME 2151 TIMEEM = TIMEEM + DTIME 2152C 2153 ENDIF 2154C 2155C----------------------------- 2156C write out result vector. 2157C----------------------------- 2158C 2159 IF (.NOT. CCS) THEN 2160 DTIME = SECOND() 2161 CALL CC_WVEC(LUTR,TRFIL,NRHO2,NT2AM(ISYRES),1,WORK(KRHO2)) 2162 DTIME = SECOND() - DTIME 2163 TIMIO = TIMIO + DTIME 2164 ENDIF 2165C 2166C 2167C========================================================== 2168C MO basis T2SQ section. 2169C Contract intermediates constructed in loop with T2SQ. 2170C========================================================== 2171C 2172C------------------------------------- 2173C Read in T2 amplitude in RHO2. 2174C------------------------------------- 2175C 2176 DTIME = SECOND() 2177C 2178 IF (.NOT. CCS) THEN 2179 IOPT = 2 2180 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KRHO2)) 2181C 2182 DTIME = SECOND() - DTIME 2183 TIMIO = TIMIO + DTIME 2184C 2185C----------------------------------- 2186C Square up T2 amplitudes. 2187C----------------------------------- 2188C 2189 DTIME = SECOND() 2190 CALL CC_T2SQ(WORK(KRHO2),WORK(KC2AM),1) 2191 DTIME = SECOND() - DTIME 2192 TIMT2SQ = TIMT2SQ + DTIME 2193C 2194 IF (IPRINT.GT.50) THEN 2195 CALL AROUND('CC_FMAT: (T1,T2) vector readin') 2196 CALL CC_PRSQ(WORK(KL1AM),WORK(KC2AM),1,0,1) 2197 ENDIF 2198C 2199C--------------------------------- 2200C Read in result vector. 2201C--------------------------------- 2202C 2203 DTIME = SECOND() 2204 CALL CC_RVEC(LUTR,TRFIL,NRHO2,NT2AM(ISYRES),1,WORK(KRHO2)) 2205 DTIME = SECOND() - DTIME 2206 TIMIO = TIMIO + DTIME 2207 ENDIF 2208C 2209C---------------------------------------------------- 2210C Calculate E-term: T2*E(C1,C2) 2211C---------------------------------------------------- 2212C 2213 IF ( IPRINT .GT. 50) THEN 2214 CALL AROUND( 'E-intermdiates calc. EI(C1,C2) - ao.') 2215 CALL CC_PREI(WORK(KEMAT1),WORK(KEMAT2),ISYMR,0) 2216 ENDIF 2217 IF ( DEBUG ) THEN 2218 XE1 = DDOT(NEMAT1(ISYMR),WORK(KEMAT1),1,WORK(KEMAT1),1) 2219 XE2 = DDOT(NMATIJ(ISYMR),WORK(KEMAT2),1,WORK(KEMAT2),1) 2220 WRITE(LUPRI,1) 'Norm of EI1 -ao: ',XE1 2221 WRITE(LUPRI,1) 'Norm of EI2 -ao: ',XE2 2222 XE1 = DDOT(N2BST(ISYMR),WORK(KFOCKC),1,WORK(KFOCKC),1) 2223 WRITE(LUPRI,1) 'Norm of FCKCC1: ',XE1 2224 ENDIF 2225C 2226 FCKCON = .TRUE. 2227 ETRAN = .TRUE. 2228 ISYMEI = ISYMR 2229 DTIME = SECOND() 2230 CALL CCRHS_EFCK(WORK(KEMAT1),WORK(KEMAT2),WORK(KLAMDH), 2231 * WORK(KFOCKC),WORK(KEND1),LWRK1,FCKCON, 2232 * ETRAN,ISYMR) 2233 DTIME = SECOND() - DTIME 2234 TIMEFCK = TIMEFCK + DTIME 2235C 2236 IF ( DEBUG ) THEN 2237 XE1 = DDOT(NMATAB(ISYMR),WORK(KEMAT1),1,WORK(KEMAT1),1) 2238 XE2 = DDOT(NMATIJ(ISYMR),WORK(KEMAT2),1,WORK(KEMAT2),1) 2239 XL1 = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 2240 WRITE(LUPRI,1) 'Norm of EI1 -mo: ',XE1 2241 WRITE(LUPRI,1) 'Norm of EI2 -mo: ',XE2 2242 ENDIF 2243C 2244 IF (IPRINT.GT.54) THEN 2245 CALL AROUND( 'E-intermdiates calc EI(C1,C2) -mo' ) 2246 CALL CC_PREI(WORK(KEMAT1),WORK(KEMAT2),ISYMR,1) 2247 ENDIF 2248C 2249 CALL CCLR_E1C1(WORK(KRHO1),WORK(KL1AM), 2250 * WORK(KEMAT1),WORK(KEMAT2), 2251 * WORK(KEND1),LWRK1,ISYML,ISYMR,'T') 2252C 2253 IF (.NOT. (CCS.OR.CC2)) THEN 2254C 2255 ISYVEC = ISYML 2256 ISYMIM = ISYMR 2257 DTIME = SECOND() 2258C 2259C------------------------------------------------------------ 2260C Prepare the E-intermediates for contraction with L2. 2261C------------------------------------------------------------ 2262C 2263 CALL CC_EITR(WORK(KEMAT1),WORK(KEMAT2),WORK(KEND1),LWRK1, 2264 * ISYMIM) 2265 IF ( DEBUG ) THEN 2266 XE1 = DDOT(NMATAB(ISYMR),WORK(KEMAT1),1,WORK(KEMAT1),1) 2267 XE2 = DDOT(NMATIJ(ISYMR),WORK(KEMAT2),1,WORK(KEMAT2),1) 2268 WRITE(LUPRI,1) 'Norm of EI1 -after EITR: ',XE1 2269 WRITE(LUPRI,1) 'Norm of EI2 -after EITR: ',XE2 2270 ENDIF 2271C 2272C 2273C--------------------------------------- 2274C L2*EI(T1,T2,R1,R2) contraction. 2275C--------------------------------------- 2276C 2277 CALL CCRHS_E(WORK(KRHO2),WORK(KL2AM), 2278 * WORK(KEMAT1),WORK(KEMAT2), 2279 * WORK(KEND1),LWRK1,ISYVEC,ISYMIM) 2280 DTIME = SECOND() - DTIME 2281 TIME = TIME + DTIME 2282C 2283 IF (IPRINT .GT. 120) THEN 2284 CALL AROUND(' AFTER CCRHS_E 1.cont. CCLR:RHO ') 2285 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,0,1) 2286 ENDIF 2287C 2288 ENDIF 2289C 2290 IF ((NFIELD .GT.0).AND.(CC2.AND.NONHF)) THEN 2291 IOPT = -1 2292 CALL CC_FCC2FF(WORK(KRHO2),WORK(KL2AM),ISYML, 2293 * WORK(KLAMDP),WORK(KLAMDH), 2294 * WORK(KLAMPC),WORK(KLAMHC),ISYMR, 2295 * WORK(KEND1),LWRK1,IOPT) 2296 ENDIF 2297C 2298 IF (DEBUG) THEN 2299 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2300 WRITE(LUPRI,1) 'Norm of Rho1 -after EI(C) terms: ',RHO1N 2301 IF (.NOT. CCS) THEN 2302 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2303 WRITE(LUPRI,1) 'Norm of Rho2 -after EI(C) terms: ',RHO2N 2304 ENDIF 2305 WRITE(LUPRI,*) ' ' 2306 ENDIF 2307C 2308C---------------------------------------- 2309C Calculate A-term from gamma matrix. 2310C---------------------------------------- 2311C 2312 IF (.NOT. (CCS .OR. CC2)) THEN 2313C 2314 IF ( DEBUG ) THEN 2315 XGAM = DDOT(NGAMMA(ISYMR),WORK(KGAMMC),1,WORK(KGAMMC),1) 2316 WRITE(LUPRI,1) 'Norm of Gamma matrix from T2MO: ',XGAM 2317 ENDIF 2318C 2319 DTIME = SECOND() 2320 IOPT = 2 2321 CALL CCRHS_A(WORK(KRHO2),WORK(KL2AM),WORK(KGAMMC), 2322 * WORK(KEND1),LWRK1, 2323 * ISYMR,ISYML,IOPT) 2324 DTIME = SECOND() - DTIME 2325 TIMA = TIMA + DTIME 2326C 2327 IF ( DEBUG ) THEN 2328 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2329 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2330 WRITE(LUPRI,1) 'Norm of Rho1 -after L*Gam(R,T): ',RHO1N 2331 WRITE(LUPRI,1) 'Norm of Rho2 -after L*Gam(R,T): ',RHO2N 2332 ENDIF 2333C 2334 ENDIF 2335C 2336C------------------------------------ 2337C Recover work from gamma matrix. 2338C------------------------------------ 2339C 2340 KEND1 = KENDS2 2341 LWRK1 = LWRKS2 2342C 2343C------------------------------------------- 2344C Calculate the D-term: L2*DI(C1,C2,T1). 2345C------------------------------------------- 2346C 2347 IF (.NOT. (CCS .OR. CC2)) THEN 2348C 2349 ISYDIM = ISYMR 2350 ISYVEC = ISYML 2351 IOPT = 2 2352C 2353 IF ( DEBUG ) THEN 2354 RHO2N = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 2355 WRITE(LUPRI,1) 'Norm of L2AM -before D-term(R,T): ',RHO2N 2356 ENDIF 2357C 2358 DTIME = SECOND() 2359 CALL CC_22D(WORK(KRHO2),WORK(KL2AM),ISYVEC,WORK(KLAMDH), 2360 * WORK(KEND1),LWRK1, 2361 * ISYDIM,LUD,DTFIL,1,IOPT) 2362 DTIME = SECOND() - DTIME 2363 TIMDIO = TIM22D + DTIME 2364C 2365 IF (IPRINT .GT. 120) THEN 2366 CALL AROUND(' AFTER (2T2-T2)*DI(C1) D-term RHO is ') 2367 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,0,1) 2368 ENDIF 2369C 2370 IF ( DEBUG ) THEN 2371 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2372 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2373 WRITE(LUPRI,1) 'Norm of Rho1 -after T2*DI(R,T): ',RHO1N 2374 WRITE(LUPRI,1) 'Norm of Rho2 -after T2*DI(R,T): ',RHO2N 2375 ENDIF 2376C 2377 ENDIF 2378C 2379C------------------------------------------- 2380C Calculate the C-term: L2*CI(C1,C2,T1). 2381C------------------------------------------- 2382C 2383 IF (.NOT. (CCS .OR. CC2)) THEN 2384C 2385 ISYCIM = ISYMR 2386 ISYVEC = ISYML 2387C 2388 IOPT = 2 2389C 2390C------------------- 2391C Prepare L2. 2392C------------------- 2393C 2394 CALL CCRHS_T2BT(WORK(KL2AM),WORK(KEND1),LWRK1,ISYML) 2395 CALL DSCAL(NT2SQ(ISYML),THREE,WORK(KL2AM),1) 2396 CALL CCSD_T2TP(WORK(KL2AM),WORK(KEND1),LWRK1,ISYML) 2397C 2398C------------------ 2399C Calculate. 2400C------------------ 2401C 2402 DTIME = SECOND() 2403 CALL CC_22C(WORK(KRHO2),WORK(KL2AM),ISYVEC,WORK(KLAMDH), 2404 * WORK(KEND1),LWRK1,ISYCIM, 2405 * LUC,CTFIL,1,IOPT) 2406 DTIME = SECOND() - DTIME 2407 TIMCIO = TIMCIO + DTIME 2408C 2409 IF (IPRINT .GT. 120) THEN 2410 CALL AROUND(' AFTER T2*CIM(C1,C2) C-term RHO is ') 2411 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,0,1) 2412 ENDIF 2413C 2414 IF ( DEBUG ) THEN 2415 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2416 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2417 WRITE(LUPRI,1) 'Norm of Rho1 -after L2*CI(R,T): ',RHO1N 2418 WRITE(LUPRI,1) 'Norm of Rho2 -after L2*CI(R,T): ',RHO2N 2419 WRITE(LUPRI,*) ' ' 2420 ENDIF 2421C 2422 ENDIF 2423C 2424C-------------------------------------- 2425C End loop over vectors in scratch. 2426C-------------------------------------- 2427C 2428 KEND1 = KENDS2 2429 LWRK1 = LWRKS2 2430C 2431C============================================================= 2432C 2433C End section. 2434C 2435C============================================================= 2436C 2437 IF (IPRINT .GT. 50 ) THEN 2438 CALL AROUND('END OF CC_FMAT :RHO ') 2439 NC2 = 1 2440 IF ( CCS ) NC2 = 0 2441 CALL CC_PRP(WORK(KRHO1),WORK(KRHO2),ISYRES,1,NC2) 2442 ENDIF 2443C 2444 IF ((IPRINT .GT.15).OR.DEBUG) THEN 2445 RHO1N = DDOT(NT1AM(ISYRES),WORK(KRHO1),1,WORK(KRHO1),1) 2446 RHO2N = DDOT(NT2AM(ISYRES),WORK(KRHO2),1,WORK(KRHO2),1) 2447 WRITE(LUPRI,1) 'Norm of Rho1 -end of CC_FMAT: ',RHO1N 2448 WRITE(LUPRI,1) 'Norm of Rho2 -end of CC_FMAT: ',RHO2N 2449 ENDIF 2450C 2451C------------- 2452C Timings. 2453C------------- 2454C 2455 TIMALL = SECOND() - TIMALL 2456C 2457 IF (IPRINT .GT. 4) THEN 2458 TIMCCSD = TIMA + TIMBF + TIMC + TIMD + TIME + 2459 * TIMF + TIMG + TIMH + TIMI + TIMEI + 2460 * TIMEFCK + TIME1C1 + TIM1C1F + TIMCIO + TIMEZ + 2461 * TIMDIO + TIMLAM + TIMAOD + TIMTRBT + TIMEXI + 2462 * TIMT2AO + TIMTCME + TIMT2TR + TIMT2BT + TIMEYI + 2463 * TIMT2MO + TIMFCKMO + TIMFCK + TIMT2SQ + TIMEMI + 2464 * TIMEH + TIMC1T2 + TIM12B + TIMEG + TIM2EM + 2465 * TIMEEM 2466 WRITE(LUPRI,'(/)') 2467 WRITE(LUPRI,9998) 'CC_FMAT in total ', TIMALL 2468 WRITE(LUPRI,9998) 'CC part ', TIMCCSD 2469 WRITE(LUPRI,9998) 'Int. calc. & read ', TIMER1 + TIMER2 2470 * + TIMRDAO 2471 WRITE(LUPRI,9998) 'IO of vectors/I.M. ', TIMIO 2472 WRITE(LUPRI,*) ' ' 2473 ENDIF 2474C 2475 IF (IPRINT .GT. 9) THEN 2476 WRITE(LUPRI,'(A)') 2477 IF ( .NOT. (CC2.OR.CCS)) THEN 2478 WRITE(LUPRI,9999) 'A ', TIMA 2479 WRITE(LUPRI,9999) 'BF ', TIMBF 2480 WRITE(LUPRI,9999) 'C ', TIMC 2481 WRITE(LUPRI,9999) 'CIO ', TIMCIO 2482 WRITE(LUPRI,9999) 'D ', TIMD 2483 WRITE(LUPRI,9999) 'DIO ', TIMDIO 2484 WRITE(LUPRI,9999) 'ZWV ', TIMEZ 2485 ENDIF 2486 IF ( .NOT. CCS ) THEN 2487 WRITE(LUPRI,9999) 'E ', TIME 2488 IF (CC2) WRITE(LUPRI,9999) 'F ', TIMF 2489 WRITE(LUPRI,9999) 'G ', TIMG 2490 WRITE(LUPRI,9999) '21G ', TIMEG 2491 WRITE(LUPRI,9999) '22EM ', TIM2EM 2492 WRITE(LUPRI,9999) 'LT21EFM ', TIMEEM 2493 WRITE(LUPRI,9999) '12C ', TIME12C 2494 WRITE(LUPRI,9999) 'H+21H ', TIMH+TIMEH 2495 WRITE(LUPRI,9999) 'I ', TIMI 2496 WRITE(LUPRI,9999) 'EI+21I ', TIMEI 2497 WRITE(LUPRI,9999) 'X,Y,M ', TIMEXI+TIMEYI+TIMEMI 2498 WRITE(LUPRI,9999) 'INT3O ', TIME3O 2499 WRITE(LUPRI,9999) 'CC_C1T2C ', TIMC1T2 2500 ENDIF 2501 WRITE(LUPRI,9999) 'EFCK ', TIMEFCK 2502 WRITE(LUPRI,9999) 'E1C1 ', TIME1C1 2503 WRITE(LUPRI,9999) '1C1F ', TIM1C1F 2504 WRITE(LUPRI,9999) 'LAMTRA ', TIMLAM 2505 WRITE(LUPRI,9999) 'AODENS ', TIMAOD 2506 WRITE(LUPRI,9999) 'FCK ', TIMFCK 2507 WRITE(LUPRI,9999) 'TRBT ', TIMTRBT 2508 WRITE(LUPRI,9999) '21DC ', TIMEC 2509 WRITE(LUPRI,9999) '12B ', TIM12B 2510 IF ( .NOT. CCS) THEN 2511 WRITE(LUPRI,9999) 'T2AO ', TIMT2AO 2512 WRITE(LUPRI,9999) 'TCME ', TIMTCME 2513 WRITE(LUPRI,9999) 'T2TR ', TIMT2TR 2514 WRITE(LUPRI,9999) 'T2BT ', TIMT2BT 2515 WRITE(LUPRI,9999) 'T2SQ ', TIMT2SQ 2516 WRITE(LUPRI,9999) 'T2MO ', TIMT2MO 2517 ENDIF 2518 WRITE(LUPRI,9999) 'FCKMO ', TIMFCKMO 2519C 2520 WRITE(LUPRI,'(A)') 2521 WRITE(LUPRI,9999) 'RDAO ', TIMRDAO 2522 WRITE(LUPRI,9999) 'ERIDI1 ', TIMER1 2523 WRITE(LUPRI,9999) 'ERIDI2 ', TIMER2 2524C 2525 IF (CCSDT) THEN 2526 WRITE(LUPRI,'(A)') 2527 WRITE(LUPRI,9999) 'T3INT ', TIMT3I 2528 WRITE(LUPRI,9999) 'OMEG-1 ', TIMO31 2529 WRITE(LUPRI,9999) 'OMEG-2 ', TIMO32 2530 WRITE(LUPRI,9999) 'OMEG-3 ', TIMO33 2531 ENDIF 2532 WRITE(LUPRI,*) ' ' 2533 ENDIF 2534 IF (IPRINT .GT. 15) THEN 2535 IF (DIRECT) WRITE(LUPRI,*) ' Atomic direct calculation' 2536 CALL AROUND(' END OF CC_FMAT ') 2537 ENDIF 2538C 2539 1 FORMAT(1x,A35,1X,E20.10) 25409998 FORMAT(1x,'Time used in',2x,A20,2x,': ',f10.2,' seconds') 25419999 FORMAT(1x,'Time used in',2x,A18,2x,': ',f10.2,' seconds') 2542C 2543C 2544C----------------- 2545C Close files. 2546C----------------- 2547C 2548 IF (.NOT.(CCS.OR.CC2)) THEN 2549 CALL WCLOSE2(LUC,CTFIL,'DELETE') 2550 CALL WCLOSE2(LUD,DTFIL,'DELETE') 2551C 2552 CALL WCLOSE2(LUCIM,CFIL,'KEEP') 2553 CALL WCLOSE2(LUDIM,DFIL,'KEEP') 2554C 2555 ENDIF 2556C 2557C------------------------------------------------ 2558C Close (and delete) files for intermediates: 2559C------------------------------------------------ 2560C 2561 IF (.NOT.(CCS.OR.CC2)) THEN 2562 IF (LLIST(1:2).EQ.'L0') THEN 2563 CALL WCLOSE2(LUPQR0,FILPQIM,'KEEP') 2564 ELSE 2565 CALL WCLOSE2(LUPQR0,FILPQR0,'DELETE') 2566 END IF 2567 CALL WCLOSE2(LUPQRR,FILPQRR,'DELETE') 2568 END IF 2569C 2570 CALL WCLOSE2(LUO3, O3FIL, 'DELETE') 2571 CALL WCLOSE2(LUO3T,O3TFIL,'DELETE') 2572C 2573 IF ( .NOT. CCS) THEN 2574 CALL WCLOSE2(LUTR, TRFIL,'DELETE') 2575 ENDIF 2576C 2577 IF (CCSD) THEN 2578 CALL WCLOSE2(LUBFI,FNBFI,'DELETE') 2579 CALL WCLOSE2(LUBFD,FNBFD,'DELETE') 2580 END IF 2581C 2582C------------------------------- 2583C Restore T2TCOR and OMEGOR. 2584C------------------------------- 2585C 2586 T2TCOR = T2TSAV 2587 OMEGOR = ORSAVE 2588C 2589 RETURN 2590 2591*---------------------------------------------------------------------* 2592* handle I/O error: 2593*---------------------------------------------------------------------* 2594999 CONTINUE 2595 WRITE (LUPRI,*) 'I/O error in CC_FMAT.' 2596 CALL QUIT('I/O error in CC_FMAT.') 2597 2598 CALL QEXIT('CC_FMATOLD') 2599 2600 END 2601*=====================================================================* 2602* END OF SUBROUTINE CC_FMAT 2603*=====================================================================* 2604 SUBROUTINE CCLR_FTST(WORK,LWORK) 2605C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2606C 2607C Written by Ove Christiansen April 1996, 2608C to calculate the F-matrix by finite difference and by 2609C construction with analytical transformations. 2610C This does not include the HF contribution to the F-matrix. 2611C (<HF|[[H,taumu],tauny]|HF>) 2612C 2613C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2614C 2615#include "implicit.h" 2616#include "priunit.h" 2617#include "maxorb.h" 2618#include "iratdef.h" 2619#include "ccorb.h" 2620#include "aovec.h" 2621C 2622 DIMENSION WORK(LWORK) 2623 CHARACTER*(2) LLIST 2624C 2625#include "ccsdinp.h" 2626#include "cclr.h" 2627#include "ccsdsym.h" 2628#include "ccsdio.h" 2629#include "leinf.h" 2630C 2631 CALL QENTER('CCLR_FTST') 2632C 2633 IF (IPRINT .GT. 10) THEN 2634 CALL AROUND(' START OF CCLR_FTST ') 2635 ENDIF 2636C 2637 IF (NSYM.GT.1) THEN 2638 WRITE(LUPRI,*) 'FTST does only work without symmetry' 2639 CALL QUIT('FTST does only work without symmetry') 2640 ENDIF 2641C 2642 LLIST = 'L1' 2643 ILLNR = 3 2644C 2645C-------------------- 2646C Initialization. 2647C-------------------- 2648C 2649 ISYMTR = 1 2650 IPRLE = IPRINT 2651 NC1VEC = NT1AM(ISYMTR) 2652 NC2VEC = NT2AM(ISYMTR) 2653 NTEMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 2654 NTEMP2 = NTEMP*(NC1VEC + NC2VEC) 2655C 2656C------------------------------------------------------------ 2657C Calculate F-matrix by Transformation with unit vectors. 2658C First elements in workspace contains the jacobian. 2659C------------------------------------------------------------ 2660C 2661 IF (IPRINT.GT.10) WRITE(LUPRI,*) 'CCLR_FTST Workspace:',LWORK 2662C 2663 KFANA = 1 2664 KEND1 = 1 + NTEMP2 2665 LWRK1 = LWORK - KEND1 2666 IF ( LWRK1 .LE. 0 ) CALL QUIT('Too little workspace in CCLR_FTST') 2667 CALL AROUND( 'Calculation of analytical F-matrix. ') 2668 CALL CCLR_FANA(LLIST,ILLNR,NC1VEC,NC2VEC,WORK(KFANA),LWORK) 2669C 2670C------------------------------------------------------- 2671C Calculate F-matrix by finite difference. 2672C First elements in workspace contains the F-matrix. 2673C------------------------------------------------------- 2674C 2675 KFDF = KEND1 2676 KEND2 = KFDF + NTEMP2 2677 LWRK2 = LWORK - KEND2 2678 IF ( LWRK2 .LE. 0 ) 2679 & CALL QUIT('Too little workspace in CCLR_FTST-2') 2680 CALL AROUND( 'Calculation of finite difference CC F-matrix ') 2681 CALL CCLR_FDF(LLIST,ILLNR,NC1VEC,NC2VEC,WORK(KEND1),LWRK1) 2682C 2683C----------------------------------------------------- 2684C calculate difference between the two F-matrices. 2685C----------------------------------------------------- 2686C 2687 IF (.TRUE.) THEN 2688 CALL DAXPY(NTEMP2,-1.0D00,WORK(KFDF),1,WORK(KFANA),1) 2689 IF ( IPRINT .GT. 40) THEN 2690 CALL AROUND( 'DIFFERENCE. CC F-matrix - 11 & 21 PART ' ) 2691 CALL OUTPUT(WORK(KFANA),1,NTEMP,1,NC1VEC,NTEMP,NC1VEC,1, 2692 & LUPRI) 2693 CALL AROUND( 'DIFFERENCE. CC F-matrix - 12 & 22 PART ' ) 2694 CALL OUTPUT(WORK(KFANA+NTEMP*NC1VEC),1,NTEMP,1,NC2VEC, 2695 * NTEMP,NC2VEC,1,LUPRI) 2696 ENDIF 2697 DIFNOR = DDOT(NTEMP2,WORK(KFANA),1,WORK(KFANA),1) 2698 WRITE(LUPRI,*) ' ' 2699 WRITE(LUPRI,*) ' The norm of the difference between fd and '// 2700 * 'ana F-matrix is ',SQRT(DIFNOR) 2701 WRITE(LUPRI,*) ' ' 2702 WRITE(LUPRI,*) ' END OF F-TEST' 2703 WRITE(LUPRI,*) ' ' 2704 ENDIF 2705C 2706 IF (IPRINT .GT. 10) THEN 2707 CALL AROUND(' END OF CCLR_FTST') 2708 ENDIF 2709C 2710 CALL QEXIT('CCLR_FTST') 2711C 2712 RETURN 2713 END 2714*=====================================================================* 2715 SUBROUTINE CCLR_FANA(LLIST,ILLNR,NC1VEC,NC2VEC,WORK,LWORK) 2716C---------------------------------------------------------------------- 2717C Test routine for calculating the CC F-matrix by Transformation of 2718C unit vectors. 2719C Ove Christiansen April 1996 2720C---------------------------------------------------------------------- 2721#include "implicit.h" 2722#include "priunit.h" 2723#include "maxorb.h" 2724#include "iratdef.h" 2725#include "ccorb.h" 2726#include "aovec.h" 2727#include "ccsdinp.h" 2728#include "cclr.h" 2729#include "ccsdsym.h" 2730#include "ccsdio.h" 2731#include "leinf.h" 2732C 2733 DIMENSION WORK(LWORK),ITADR(2) 2734 PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00 ) 2735 CHARACTER*10 MODEL 2736 CHARACTER*(*) LLIST 2737C 2738 CALL QENTER('CCLR_FANA') 2739C 2740 IF (IPRINT .GT. 10) THEN 2741 CALL AROUND(' START OF CCLR_FANA') 2742 ENDIF 2743C 2744 MODEL = 'CCSD ' 2745 IF (CCS) MODEL = 'CCS ' 2746 IF (CC2) MODEL = 'CC2 ' 2747C 2748C---------------------------- 2749C Work space allocations. 2750C---------------------------- 2751C 2752 NTAMP = NT1AMX + NT2AMX 2753 NTAMP2 = NTAMP*(NC1VEC + NC2VEC ) 2754 KF = 1 2755 KEND1 = 1 + NTAMP2 2756 LWRK1 = LWORK - KEND1 2757 IF (LWRK1 .LT. 0 ) 2758 & CALL QUIT('INSUFFICIENT WORK SPACE IN CCLR_FANA ') 2759C 2760C---------------------------------------------------------------- 2761C Make the transformation with unit vectors. 2762C Transformed vectors are returned as first elements in work. 2763C---------------------------------------------------------------- 2764C 2765 KC1 = KEND1 2766 KC2 = KC1 + NT1AMX 2767 KEND1 = KC2 + NT2AMX 2768 LWRK1 = LWORK - KEND1 2769 IF (LWRK1 .LT. 0 ) 2770 & CALL QUIT('INSUFFICIENT WORK SPACE IN CCLR_FANA ') 2771 2772 CALL DZERO(WORK(KC1),NTAMP) 2773 2774 DO IVEC = 1, NC1VEC+NC2VEC 2775 2776 WORK(KC1-1+IVEC) = 1.0d0 2777 2778 IOPT=3 2779 CALL CC_WRRSP('D0',0,1,IOPT,MODEL,DUMMY,WORK(KC1),WORK(KC2), 2780 & WORK(KEND1),LWRK1) 2781 2782 WORK(KC1-1+IVEC) = 0.0d0 2783 2784 CALL CC_FMATOLD(LLIST,ILLNR,'D0',1,WORK(KEND1),LWRK1) 2785 2786 KOFF3 = KF + NTAMP*(IVEC-1) 2787 CALL DCOPY(NTAMP,WORK(KEND1),1,WORK(KOFF3),1) 2788 2789 END DO 2790C 2791C--------------------- 2792C Calculate norms. 2793C--------------------- 2794C 2795 XNJ = DDOT(NTAMP2,WORK(KF),1,WORK(KF),1) 2796C 2797 WRITE(LUPRI,*) ' ' 2798 WRITE(LUPRI,*) ' NORM OF UNIT VEC. TRA. F.', SQRT(XNJ) 2799 WRITE(LUPRI,*) ' ' 2800C 2801 CALL CCLR_NORMS(X11,X21,XR1,X12,X22,XR2,X1R,X2R,XRR, 2802 & WORK(KF),NTAMP,NC1VEC,NC2VEC) 2803C 2804 WRITE(LUPRI,*) 'NORM OF 11 PART OF UNIT VECT. F. ', SQRT(X11) 2805 WRITE(LUPRI,*) 'NORM OF 21 PART OF UNIT VECT. F. ', SQRT(X21) 2806 WRITE(LUPRI,*) 'NORM OF 12 PART OF UNIT VECT. F. ', SQRT(X12) 2807 WRITE(LUPRI,*) 'NORM OF 22 PART OF UNIT VECT. F. ', SQRT(X22) 2808C 2809C------------------------------------------------ 2810C Print the columns of the F-Matrix obtained. 2811C------------------------------------------------ 2812C 2813 IF (IPRINT .GT. 4) THEN 2814C 2815 write(LUPRI,'(/,1X,A,/)') 2816 * 'F matrix without <HF| [[H,tau,mu],tau,ny]|HF> cont. ' 2817 CALL AROUND( 'CC analytical F-Matrix - F*C1 PART ' ) 2818 CALL OUTPUT(WORK(KF),1,NTAMP,1,NC1VEC,NTAMP,NC1VEC,1,LUPRI) 2819 CALL AROUND( 'CC analytical F-Matrix - F*C2 PART ' ) 2820 CALL OUTPUT(WORK(KF+NTAMP*NC1VEC),1,NTAMP,1,NC2VEC, 2821 * NTAMP,NC2VEC,1,LUPRI) 2822C 2823 WRITE(LUPRI,*) ' ' 2824 WRITE(LUPRI,*) ' NORM OF FIN. DIFF. F-Matrix.', SQRT(XNJ) 2825 WRITE(LUPRI,*) ' ' 2826 WRITE(LUPRI,*) ' NORM OF 11 PART OF FD. F-mat.: ', SQRT(X11) 2827 WRITE(LUPRI,*) ' NORM OF 21 PART OF FD. F-mat.: ', SQRT(X21) 2828 WRITE(LUPRI,*) ' NORM OF 12 PART OF FD. F-mat.: ', SQRT(X12) 2829 WRITE(LUPRI,*) ' NORM OF 22 PART OF FD. F-mat.: ', SQRT(X22) 2830 ENDIF 2831C 2832 IF (IPRINT .GT. 10) THEN 2833 CALL AROUND(' END OF CCLR_FANA ') 2834 ENDIF 2835C 2836 CALL QEXIT('CCLR_FANA') 2837C 2838 RETURN 2839 END 2840*=====================================================================* 2841 SUBROUTINE CCLR_FDF(LLIST,ILLNR,NC1VEC,NC2VEC,WORK,LWORK) 2842C 2843C---------------------------------------------------------------------- 2844C Test routine for calculating the CC F_matrix by 2845C finite difference on the cc-left-hand transformation. 2846C Ove Christiansen April 1996 2847C--------------------------------------------------------------------- 2848C 2849#include "implicit.h" 2850#include "priunit.h" 2851#include "dummy.h" 2852#include "maxorb.h" 2853#include "iratdef.h" 2854#include "ccorb.h" 2855#include "aovec.h" 2856#include "ccsdinp.h" 2857#include "cclr.h" 2858#include "ccsdsym.h" 2859#include "ccsdio.h" 2860#include "leinf.h" 2861C 2862 DIMENSION WORK(LWORK),ITADR(2) 2863 PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00, DELTA = 1.0D-07) 2864 CHARACTER MODEL*10, LLIST(*), APROXR12*3 2865 LOGICAL L1TST,L2TST 2866C 2867 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 2868C 2869 CALL QENTER('CCLR_FDF') 2870C 2871 MODEL = 'CCSD ' 2872 IF (CCS) MODEL = 'CCS ' 2873 IF (CC2) MODEL = 'CC2 ' 2874C 2875 IF (CCR12) CALL QUIT('CCLR_FDF for CCR12 not adapted; use '// 2876 & 'CC_FTSTNEW / CC_FDF') 2877C 2878 IF (IPRINT.GT.5) THEN 2879 CALL AROUND( 'IN CCLR_FDF : MAKING FINITE DIFF. CC F-Matrix') 2880 ENDIF 2881C 2882C---------------------- 2883C Set Test options. 2884C---------------------- 2885C 2886 L1TST = .FALSE. 2887 L2TST = .FALSE. 2888C 2889C---------------------------- 2890C Work space allocations. 2891C---------------------------- 2892C 2893 ISYMTR = 1 2894 ISYMOP = 1 2895C 2896 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 2897 NTAMP2 = NTAMP*(NC1VEC + NC2VEC ) 2898 KF = 1 2899 KRHO1 = KF + NTAMP2 2900 KRHO2 = KRHO1 + NT1AMX 2901 KC1AM = KRHO2 + MAX(NT2AMX,NT2AM(ISYMTR)) 2902 KC2AM = KC1AM + NT1AM(ISYMTR) 2903 KEND1 = KC2AM 2904 * + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR), 2905 * 2*NT2ORT(ISYMTR)) 2906 LWRK1 = LWORK - KEND1 2907C 2908 KRHO1D = KEND1 2909 KRHO2D = KRHO1D + NT1AMX 2910 KEND2 = KRHO2D 2911 * + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR), 2912 * 2*NT2ORT(ISYMTR)) 2913 LWRK2 = LWORK - KEND1 2914C 2915 IF (IPRINT .GT. 100 ) THEN 2916 WRITE(LUPRI,*) ' IN CCLR_FDF: KF = ',KF 2917 WRITE(LUPRI,*) ' IN CCLR_FDF: KRHO1 = ',KRHO1 2918 WRITE(LUPRI,*) ' IN CCLR_FDF: KRHO2 = ',KRHO2 2919 WRITE(LUPRI,*) ' IN CCLR_FDF: KC1AM = ',KC1AM 2920 WRITE(LUPRI,*) ' IN CCLR_FDF: KC2AM = ',KC2AM 2921 WRITE(LUPRI,*) ' IN CCLR_FDF: KRHO1D = ',KRHO1D 2922 WRITE(LUPRI,*) ' IN CCLR_FDF: KRHO2D = ',KRHO2D 2923 WRITE(LUPRI,*) ' IN CCLR_FDF: KEND2 = ',KEND2 2924 WRITE(LUPRI,*) ' IN CCLR_FDF: LWRK2 = ',LWRK2 2925 ENDIF 2926 IF (LWRK2.LT.0 ) THEN 2927 WRITE(LUPRI,*) 'Too little work space in cclr_fdf ' 2928 WRITE(LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 2929 WRITE(LUPRI,*) 'NEEDED (AT LEAST) = ',KEND2 2930 CALL QUIT('TOO LITTLE WORKSPACE IN CCLR_FDF ') 2931 ENDIF 2932 KF2 = KF + NC1VEC*NTAMP 2933C 2934C--------------------- 2935C Initializations. 2936C--------------------- 2937C 2938 CALL DZERO(WORK(KC1AM),NT1AMX) 2939 CALL DZERO(WORK(KC2AM),NT2AMX) 2940 CALL DZERO(WORK(KF),NTAMP2) 2941 IF (ABS(DELTA) .GT. 1.0D-15 ) THEN 2942 DELTAI = 1.0D00/DELTA 2943 ELSE 2944 DELTAI = 1 2945 ENDIF 2946 X11 = 0.0D00 2947 X12 = 0.0D00 2948 X21 = 0.0D00 2949 X22 = 0.0D00 2950 XNJ = 0.0D00 2951C 2952C------------------------------------------------ 2953C Read the CC reference amplitudes From disk. 2954C------------------------------------------------ 2955C 2956 IOPT = 3 2957 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KC1AM),WORK(KC2AM)) 2958C 2959C---------------------------------------------- 2960C Save the CC reference amplitudes on disk. 2961C---------------------------------------------- 2962C 2963 LUTAM = -1 2964 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 2965 * .FALSE.) 2966 REWIND(LUTAM) 2967 WRITE(LUTAM) (WORK(KC1AM + I -1 ), I = 1, NT1AMX) 2968 WRITE(LUTAM) (WORK(KC2AM + I -1 ), I = 1, NT2AMX) 2969 CALL GPCLOSE(LUTAM,'KEEP') 2970C 2971 IF (IPRINT.GT.125) THEN 2972 RHO1N = DDOT(NT1AMX,WORK(KC1AM),1,WORK(KC1AM),1) 2973 RHO2N = DDOT(NT2AMX,WORK(KC2AM),1,WORK(KC2AM),1) 2974 WRITE(LUPRI,*) 'Norm of T1AM: ',RHO1N 2975 WRITE(LUPRI,*) 'Norm of T2AM: ',RHO2N 2976 CALL CC_PRP(WORK(KC1AM),WORK(KC2AM),1,1,1) 2977 ENDIF 2978 RSPIM = .TRUE. 2979C 2980C------------------------------------------ 2981C Get the CC left amplitudes from disk. 2982C------------------------------------------ 2983C 2984 IOPT = 3 2985 ISYML = ILSTSYM(LLIST,ILLNR) 2986 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODEL, 2987 & WORK(KC1AM),WORK(KC2AM)) 2988C 2989C------------------------------------ 2990C For Test zero part of L vector. 2991C------------------------------------ 2992C 2993 IF ( L1TST ) THEN 2994 CALL DZERO(WORK(KC2AM),NT2AMX) 2995 ENDIF 2996 IF ( L2TST ) THEN 2997 CALL DZERO(WORK(KC1AM),NT1AMX) 2998 ENDIF 2999C 3000 IF (IPRINT.GT.125) THEN 3001 RHO1N = DDOT(NT1AMX,WORK(KC1AM),1,WORK(KC1AM),1) 3002 RHO2N = DDOT(NT2AMX,WORK(KC2AM),1,WORK(KC2AM),1) 3003 WRITE(LUPRI,*) 'Norm of T1AM: ',RHO1N 3004 WRITE(LUPRI,*) 'Norm of T2AM: ',RHO2N 3005 CALL CC_PRP(WORK(KC1AM),WORK(KC2AM),1,1,1) 3006 ENDIF 3007C 3008C------------------------------------------ 3009C Calculate reference omega=LHS vector. 3010C------------------------------------------ 3011C 3012 CALL DCOPY(NT1AMX,WORK(KC1AM),1,WORK(KRHO1D),1) 3013 CALL DCOPY(NT2AMX,WORK(KC2AM),1,WORK(KRHO2D),1) 3014 CALL CC_ATRR(0.0D0,1,-1,WORK(KRHO1D),LWRK1,.FALSE.,DUMMY, 3015 & APROXR12,.FALSE.) 3016C 3017C------------------------- 3018C Zero out components. 3019C------------------------- 3020C 3021 IF (LCOR .OR. LSEC) THEN 3022C 3023 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 3024C 3025 ENDIF 3026C 3027 IF (IPRINT.GT.2) THEN 3028 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 3029 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 3030 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ref' 3031 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ref' 3032 ENDIF 3033 IF (IPRINT.GT.125) THEN 3034 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 3035 ENDIF 3036C 3037 CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1),1) 3038 CALL DCOPY(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2),1) 3039C 3040C============================================= 3041C Calculate F-matrix by finite difference. 3042C============================================= 3043C 3044 DO 100 I = 1, NC1VEC 3045C 3046C---------------------------------------- 3047C Add finite displadement to t and 3048C calculate new intermediates. 3049C---------------------------------------- 3050C 3051 LUTAM = -1 3052 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 3053 * .FALSE.) 3054 READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX) 3055 READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX) 3056 CALL GPCLOSE(LUTAM,'KEEP') 3057C 3058 TI = SECOND() 3059 WORK(KC1AM +I -1) = WORK(KC1AM +I -1 ) + DELTA 3060 IF (LCOR .OR. LSEC) THEN 3061 CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR) 3062 ENDIF 3063C 3064 IOPT = 3 3065 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 3066 * WORK(KC2AM),WORK(KEND2),LWRK2) 3067C 3068 RSPIM = .TRUE. 3069 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 3070 * WORK(KC2AM),WORK(KEND2),LWRK2,'XXX') 3071C 3072C--------------------------------------------- 3073C Get the CC left amplitudes from disk. 3074C--------------------------------------------- 3075C 3076 IOPT = 3 3077 ISYML = ILSTSYM(LLIST,ILLNR) 3078 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODEL, 3079 & WORK(KC1AM),WORK(KC2AM)) 3080C 3081C--------------------------------------- 3082C For Test zero part of L vector. 3083C--------------------------------------- 3084C 3085 IF ( L1TST ) THEN 3086 CALL DZERO(WORK(KC2AM),NT2AMX) 3087 ENDIF 3088 IF ( L2TST ) THEN 3089 CALL DZERO(WORK(KC1AM),NT1AMX) 3090 ENDIF 3091C 3092C------------------ 3093C Transform. 3094C------------------ 3095C 3096 CALL DCOPY(NT1AMX,WORK(KC1AM),1,WORK(KRHO1D),1) 3097 CALL DCOPY(NT2AMX,WORK(KC2AM),1,WORK(KRHO2D),1) 3098 CALL CC_ATRR(0.0D0,1,-1,WORK(KRHO1D),LWRK1,.FALSE.,DUMMY, 3099 & APROXR12,.FALSE.) 3100C 3101 IF (LCOR .OR. LSEC) THEN 3102 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 3103 ENDIF 3104C 3105 IF (IPRINT.GT.2) THEN 3106 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 3107 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 3108 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ai=',I 3109 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ai=',I 3110 ENDIF 3111 IF (IPRINT.GT.125) THEN 3112 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 3113 ENDIF 3114 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 3115 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 3116 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 3117 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 3118 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 3119 * WORK(KF+NTAMP*(I-1)),1) 3120 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 3121 * WORK(KF+NTAMP*(I-1)+NT1AMX),1) 3122 X11 = X11 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 3123 X21 = X21 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 3124C 3125 TI = SECOND() - TI 3126 IF (IPRINT.GT.5 ) THEN 3127 WRITE(LUPRI,*) ' ' 3128 WRITE(LUPRI,*) 'FDF ROW NR. ',I,' DONE IN ',TI,' SEC.' 3129 ENDIF 3130C 3131 100 CONTINUE 3132C 3133C---------------------------------------------------------------- 3134C Loop over T2 amplitudes. Take care of diagonal t2 elements 3135C is in a different convention in the energy code. 3136C Factor 1/2 from right , and factor 2 from left. 3137C---------------------------------------------------------------- 3138C 3139 DO 200 NAI = 1, NT1AMX 3140 DO 300 NBJ = 1, NAI 3141 I = INDEX(NAI,NBJ) 3142C 3143 IF (I.LE.NC2VEC) THEN 3144C 3145C-------------------------------------------- 3146C Add finite displacement to t and 3147C calculate new intermediates. 3148C------------------------------------------- 3149C 3150 LUTAM = -1 3151 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED', 3152 * IDUMMY,.FALSE.) 3153 READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX) 3154 READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX) 3155 CALL GPCLOSE(LUTAM,'KEEP') 3156C 3157 TI = SECOND() 3158 DELT = DELTA 3159 IF (NAI.EQ.NBJ) DELT = 2*DELTA 3160 WORK(KC2AM + I -1) = WORK(KC2AM+I -1) + DELT 3161 IF (LCOR .OR. LSEC) THEN 3162 CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR) 3163 ENDIF 3164C 3165 IOPT = 3 3166 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 3167 * WORK(KC2AM),WORK(KEND2),LWRK2) 3168C 3169 RSPIM = .TRUE. 3170 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 3171 * WORK(KC2AM),WORK(KEND2),LWRK2,'XXX') 3172C 3173C----------------------------------------------- 3174C Get the CC left amplitudes from disk. 3175C----------------------------------------------- 3176C 3177 IOPT = 3 3178 ISYML = ILSTSYM(LLIST,ILLNR) 3179 CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODEL, 3180 & WORK(KC1AM),WORK(KC2AM)) 3181C 3182C----------------------------------------- 3183C For Test zero part of L vector. 3184C----------------------------------------- 3185C 3186 IF ( L1TST ) THEN 3187 CALL DZERO(WORK(KC2AM),NT2AMX) 3188 ENDIF 3189 IF ( L2TST ) THEN 3190 CALL DZERO(WORK(KC1AM),NT1AMX) 3191 ENDIF 3192C 3193 RHO1N = DDOT(NT1AMX,WORK(KC1AM),1,WORK(KC1AM),1) 3194 RHO2N = DDOT(NT2AMX,WORK(KC2AM),1,WORK(KC2AM),1) 3195 IF ( DEBUG ) THEN 3196 WRITE(LUPRI,*) 'Norm of L1AM-inp: ',RHO1N 3197 WRITE(LUPRI,*) 'Norm of L2AM-inp: ',RHO2N 3198 ENDIF 3199C 3200C-------------------- 3201C Transform. 3202C-------------------- 3203C 3204 CALL DCOPY(NT1AMX,WORK(KC1AM),1,WORK(KRHO1D),1) 3205 CALL DCOPY(NT2AMX,WORK(KC2AM),1,WORK(KRHO2D),1) 3206 CALL CC_ATRR(0.0D0,1,-1,WORK(KRHO1D),LWRK1,.FALSE.,DUMMY, 3207 & APROXR12,.FALSE.) 3208C 3209 IF (LCOR .OR. LSEC) THEN 3210 CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR) 3211 ENDIF 3212C 3213 IF (IPRINT.GT.2) THEN 3214 RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 3215 RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 3216 WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'aibj=',I 3217 WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'aibj=',I 3218 ENDIF 3219 IF (IPRINT.GT.125) THEN 3220 CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1) 3221 ENDIF 3222C 3223 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 3224 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 3225 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 3226 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 3227C 3228 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 3229 * WORK(KF2+NTAMP*(I-1)),1) 3230 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 3231 * WORK(KF2+NTAMP*(I-1)+NT1AMX),1) 3232C 3233 X12 = X12 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 3234 X22 = X22 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 3235 TI = SECOND() - TI 3236 IF (IPRINT.GT.5 ) THEN 3237 WRITE(LUPRI,*) ' ' 3238 WRITE(LUPRI,*) 'FDF ROW NR. ',I+NT1AMX, 3239 * ' DONE IN ',TI,' SEC.' 3240 ENDIF 3241C 3242 ENDIF 3243C 3244 300 CONTINUE 3245 200 CONTINUE 3246C 3247 WRITE(LUPRI,*) ' ' 3248 WRITE(LUPRI,*) '** FINITE DIFF WITH DELTA ',DELTA, '**' 3249 WRITE(LUPRI,*) ' ' 3250 IF ((IPRINT .GT. 4)) THEN 3251 write(LUPRI,'(/,1X,A,/)') 3252 * 'F matrix without <HF| [[H,tau,mu],tau,ny]|HF> cont. ' 3253 CALL AROUND( 'FINITE DIFF. CC F-Matrix - 11 & 21 PART ' ) 3254 CALL OUTPUT(WORK(KF),1,NTAMP,1,NC1VEC,NTAMP,NC1VEC,1,LUPRI) 3255 CALL AROUND( 'FINITE DIFF. CC F-Matrix - 12 & 22 PART ' ) 3256 CALL OUTPUT(WORK(KF+NTAMP*NC1VEC),1,NTAMP,1,NC2VEC, 3257 * NTAMP,NC2VEC,1,LUPRI) 3258 ENDIF 3259 IF (.TRUE.) THEN 3260 XNJ = X11 + X12 + X21 + X22 3261 WRITE(LUPRI,*) ' ' 3262 WRITE(LUPRI,*) ' NORM OF FIN. DIFF. F-Matrix.', SQRT(XNJ) 3263 WRITE(LUPRI,*) ' ' 3264 WRITE(LUPRI,*) ' NORM OF 11 PART OF FD. F-mat.: ', SQRT(X11) 3265 WRITE(LUPRI,*) ' NORM OF 21 PART OF FD. F-mat.: ', SQRT(X21) 3266 WRITE(LUPRI,*) ' NORM OF 12 PART OF FD. F-mat.: ', SQRT(X12) 3267 WRITE(LUPRI,*) ' NORM OF 22 PART OF FD. F-mat.: ', SQRT(X22) 3268 ENDIF 3269C 3270C------------------------------------------------- 3271C Restore the CC reference amplitudes on disk. 3272C------------------------------------------------- 3273C 3274 LUTAM = -1 3275 CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY, 3276 * .FALSE.) 3277 REWIND(LUTAM) 3278 READ(LUTAM) (WORK(KC1AM + I -1 ) , I = 1, NT1AMX) 3279 READ(LUTAM) (WORK(KC2AM + I -1 ) , I = 1, NT2AMX) 3280 CALL GPCLOSE(LUTAM,'DELETE') 3281C 3282 IOPT = 3 3283 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM), 3284 * WORK(KC2AM),WORK(KEND2),LWRK2) 3285C 3286 RSPIM = .TRUE. 3287 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 3288 & WORK(KC2AM), 3289 * WORK(KEND2),LWRK2,'XXX') 3290C 3291 IF (IPRINT .GT. 10) THEN 3292 CALL AROUND(' END OF CCLR_FDF ') 3293 ENDIF 3294C 3295 CALL QEXIT('CCLR_FDF') 3296C 3297 RETURN 3298 END 3299C---------------------------------------------------------------------- 3300*====================================================================== 3301 SUBROUTINE CC_PQIM(LLIST,ILLNR,RLIST,IRLNR,FILNAM,WORK,LWORK) 3302*---------------------------------------------------------------------- 3303* 3304* Purpose: Calculate the P and Q intermediates from the 3305* Lagrangian multipiers specified by LLIST,ILLNR 3306* and the amplitude vector specified by RLIST,IRLNR 3307* and write them to a file direct access file FILNAM 3308* 3309* P{aik,del} = sum_{dl} (2 t(dl,k;del) - t(dk,l;del)) Zeta(dl;ai) 3310* P{aik,del} = sum_{dl} t(dl,k;del) (2Zeta(di;al) + Zeta(dl;ai)) 3311* 3312* Christof Haettig, June 1998 3313* 3314*====================================================================== 3315#if defined (IMPLICIT_NONE) 3316 IMPLICIT NONE 3317#else 3318# include "implicit.h" 3319#endif 3320#include "priunit.h" 3321#include "ccorb.h" 3322#include "maxorb.h" 3323#include "ccsdsym.h" 3324#include "ccsdinp.h" 3325#include "iratdef.h" 3326#include "dummy.h" 3327#include "second.h" 3328 3329 LOGICAL LOCDBG 3330 PARAMETER (LOCDBG = .FALSE.) 3331 3332 INTEGER LUPQIM 3333 3334 INTEGER ILLNR, IRLNR, LWORK 3335 CHARACTER*(*) LLIST, RLIST, FILNAM 3336 3337#if defined (SYS_CRAY) 3338 REAL ONE, TWO, THREE, HALF, XNORM, DDOT, ZERO 3339 REAL WORK(*) 3340 REAL TIMET, TIMIO, TIMZWV, TIMTAM, DTIME, TIMTI 3341 REAL TIMZET,TIMSCL 3342#else 3343 DOUBLE PRECISION ONE, TWO, THREE, HALF, XNORM, DDOT, ZERO 3344 DOUBLE PRECISION WORK(*) 3345 DOUBLE PRECISION TIMET,TIMIO,TIMTI,TIMZWV,TIMTAM,DTIME 3346 DOUBLE PRECISION TIMZET,TIMSCL 3347#endif 3348 PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,HALF=0.5D0,ZERO=0.0D0) 3349 3350 CHARACTER*(10) MODEL 3351 INTEGER IADRPQ(MXCORB_CC+IRAT) 3352 INTEGER ISYTAMP, ISYZETA, ISYMD, ISYTIN, ISYAIK 3353 INTEGER KZETA, KTAMP, KEND1, LWRK1, KEND2, LWRK2 3354 INTEGER IERR, IADR, IOPT, IDEL, ILLL 3355 INTEGER KTINT, KTJNT, KPINT, KQINT, KLAMDH, KLAMDP, KT1AM 3356 3357 3358* external function: 3359 INTEGER ILSTSYM 3360 3361 CALL QENTER('CC_PQIM') 3362 3363*---------------------------------------------------------------------- 3364* set symmetries and allocate work space: 3365*---------------------------------------------------------------------- 3366 ISYTAMP = ILSTSYM(RLIST,IRLNR) 3367 ISYZETA = ILSTSYM(LLIST,ILLNR) 3368 3369 KZETA = 1 3370 KTAMP = KZETA + NT2SQ(ISYZETA) 3371 KLAMDH = KTAMP + MAX(NT2AM(ISYTAMP),NT2AM(ISYZETA)) 3372 KLAMDP = KLAMDH + NLAMDT 3373 KT1AM = KLAMDP + NLAMDT 3374 KEND1 = KT1AM + NT1AMX 3375 LWRK1 = LWORK - KEND1 3376 3377 IF ( (LWRK1 .LT. 0) .OR. ((LWORK-KTAMP).LT.NT2AM(ISYZETA)) ) THEN 3378 WRITE (LUPRI,*) 'Insufficient memory in CC_PQIM.' 3379 CALL QUIT('Insufficient memory in CC_PQIM.') 3380 END IF 3381 3382 TIMET = SECOND() 3383 TIMIO = ZERO 3384 TIMTI = ZERO 3385 TIMZWV = ZERO 3386 TIMTAM = ZERO 3387 TIMZET = ZERO 3388 TIMSCL = ZERO 3389 3390*---------------------------------------------------------------------- 3391* open file and initialize start address: 3392*---------------------------------------------------------------------- 3393 LUPQIM = -1 3394 CALL WOPEN2(LUPQIM,FILNAM,64,0) 3395* reserve at the beginning a record for the start addresses: 3396 IADR = (MXCORB_CC+IRAT)/IRAT + 1 3397 3398*---------------------------------------------------------------------- 3399* get XLAMD matrices from zero T1 amplitudes: 3400*---------------------------------------------------------------------- 3401 CALL DZERO(WORK(KT1AM),NT1AMX) 3402 3403 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 3404 & WORK(KEND1),LWRK1) 3405 3406*---------------------------------------------------------------------- 3407* read the zeta vector and amplitude vector into core: 3408*---------------------------------------------------------------------- 3409* first the packed zeta vector in TAMP: 3410 IOPT = 2 3411 DTIME = SECOND() 3412 CALL CC_RDRSP(LLIST,ILLNR,ISYZETA,IOPT,MODEL, 3413 & DUMMY,WORK(KTAMP)) 3414 TIMIO = TIMIO + SECOND() - DTIME 3415 3416* square up zeta vector and store now in ZETAA: 3417 DTIME = SECOND() 3418 CALL CC_T2SQ(WORK(KTAMP),WORK(KZETA),ISYZETA) 3419 TIMZET = TIMZET + SECOND() - DTIME 3420 3421* read now the amplitude vector into TAMP: 3422 IOPT = 2 3423 DTIME = SECOND() 3424 CALL CC_RDRSP(RLIST,IRLNR,ISYTAMP,IOPT,MODEL, 3425 & DUMMY,WORK(KTAMP)) 3426 TIMIO = TIMIO + SECOND() - DTIME 3427 3428* fix normalization for response vectors: 3429 IF (RLIST(1:2).NE.'R0') THEN 3430 DTIME = SECOND() 3431 CALL CCLR_DIASCL(WORK(KTAMP),TWO,ISYTAMP) 3432 TIMSCL = TIMSCL + SECOND() - DTIME 3433 END IF 3434 3435*---------------------------------------------------------------------- 3436* calculate the P intermediate in a loop over AO index delta: 3437*---------------------------------------------------------------------- 3438 DO ISYMD = 1, NSYM 3439 DO ILLL = 1, NBAS(ISYMD) 3440 IDEL = IBAS(ISYMD) + ILLL 3441 3442 ISYTIN = MULD2H(ISYMD,ISYTAMP) 3443 ISYAIK = MULD2H(ISYTIN,ISYZETA) 3444 3445 KTINT = KEND1 3446 KTJNT = KTINT + NT2BCD(ISYTIN) 3447 KPINT = KTJNT + NT2BCD(ISYTIN) 3448 KQINT = KPINT + NT2BCD(ISYAIK) 3449 KEND2 = KQINT + NT2BCD(ISYAIK) 3450 LWRK2 = LWORK - KEND2 3451 3452 IF (LWRK2 .LT. 0) THEN 3453 WRITE (LUPRI,*) 'Insufficient memory in CC_PQIM.' 3454 CALL QUIT('Insufficient memory in CC_PQIM.') 3455 END IF 3456C 3457C calculate delta batch of backtransformed amplitudes: 3458C 3459 DTIME = SECOND() 3460 CALL CC_TI(WORK(KTINT),ISYTIN,WORK(KTAMP),ISYTAMP, 3461 & WORK(KLAMDH),1,WORK(KEND2),LWRK2,IDEL,ISYMD) 3462 TIMTI = TIMTI + SECOND() - DTIME 3463C 3464C calculate 2 x Coul - Exc combination of backtransformed T: 3465C 3466 DTIME = SECOND() 3467 CALL DCOPY(NT2BCD(ISYTIN),WORK(KTINT),1,WORK(KTJNT),1) 3468 3469 CALL DSCAL(NT2BCD(ISYTIN),-ONE,WORK(KTJNT),1) 3470 3471 CALL CCLT_P21I(WORK(KTJNT),ISYTIN,WORK(KEND2),LWRK2, 3472 & IT2BCD, NT2BCD, IT1AM, NT1AM, NVIR) 3473 3474 CALL DAXPY(NT2BCD(ISYTIN),TWO,WORK(KTINT),1,WORK(KTJNT),1) 3475 TIMTAM = TIMTAM + SECOND() - DTIME 3476C 3477C calculate the P intermediate: 3478C 3479 IOPT = 1 3480 DTIME = SECOND() 3481 CALL CC_ZWVI(WORK(KPINT),WORK(KZETA),ISYZETA,WORK(KTJNT), 3482 & ISYTIN,WORK(KEND2),LWRK2,IOPT) 3483 TIMZWV = TIMZWV + SECOND() - DTIME 3484 3485 DTIME = SECOND() 3486 CALL DSCAL(NT2BCD(ISYAIK),HALF,WORK(KPINT),1) 3487 TIMSCL = TIMSCL + SECOND() - DTIME 3488C 3489C write the intermediate to file: 3490C 3491 IADRPQ(IDEL) = IADR 3492 3493 DTIME = SECOND() 3494 CALL PUTWA2(LUPQIM,FILNAM,WORK(KPINT),IADR,NT2BCD(ISYAIK)) 3495 TIMIO = TIMIO + SECOND() - DTIME 3496 3497 IADR = IADR + 2*NT2BCD(ISYAIK) 3498 3499 IF (LOCDBG) THEN 3500 WRITE (LUPRI,*) 'CC_PQIM> P interm. in AO for IDEL = ',IDEL 3501 WRITE (LUPRI,'(5F12.8)') (WORK(KPINT+I),I=0,NT2BCD(ISYAIK)) 3502 END IF 3503 3504 END DO 3505 END DO 3506 3507*---------------------------------------------------------------------- 3508* calculate (2 x Exc + Coul)/3 combination of ZETA: 3509* (we interrupt here the loop over AO to calculate the modified 3510* ZETA and accept that we recalculate the backtransformed 3511* amplitude since it the transformation of the amplitudes is 3512* less expansive than the transformation and restruction of ZETA 3513* for each delta. both the transformation of TAMP and the 3514* transformation/restruction of ZETA scale with N V^2 O^2. ) 3515*---------------------------------------------------------------------- 3516 DTIME = SECOND() 3517 CALL CCRHS_T2BT(WORK(KZETA),WORK(KEND2),LWRK2,ISYZETA) 3518 CALL CCSD_T2TP(WORK(KZETA),WORK(KEND2),LWRK2,ISYZETA) 3519 TIMZET = TIMZET + SECOND() - DTIME 3520 3521*---------------------------------------------------------------------- 3522* calculate the Q intermediate in a loop over AO index delta: 3523*---------------------------------------------------------------------- 3524 DO ISYMD = 1, NSYM 3525 DO ILLL = 1, NBAS(ISYMD) 3526 IDEL = IBAS(ISYMD) + ILLL 3527 3528 ISYTIN = MULD2H(ISYMD,ISYTAMP) 3529 ISYAIK = MULD2H(ISYTIN,ISYZETA) 3530 3531 KTINT = KEND1 3532 KTJNT = KTINT + NT2BCD(ISYTIN) 3533 KPINT = KTJNT + NT2BCD(ISYTIN) 3534 KQINT = KPINT + NT2BCD(ISYAIK) 3535 KEND2 = KQINT + NT2BCD(ISYAIK) 3536 LWRK2 = LWORK - KEND2 3537 3538 IF (LWRK2 .LT. 0) THEN 3539 WRITE (LUPRI,*) 'Insufficient memory in CC_PQIM.' 3540 CALL QUIT('Insufficient memory in CC_PQIM.') 3541 END IF 3542C 3543C calculate delta batch of backtransformed amplitudes: 3544C 3545 DTIME = SECOND() 3546 CALL CC_TI(WORK(KTINT),ISYTIN,WORK(KTAMP),ISYTAMP, 3547 & WORK(KLAMDH),1,WORK(KEND2),LWRK2,IDEL,ISYMD) 3548 TIMTI = TIMTI + SECOND() - DTIME 3549C 3550C calculate Q intermediate: 3551C 3552 IOPT = 2 3553 DTIME = SECOND() 3554 CALL CC_ZWVI(WORK(KQINT),WORK(KZETA),ISYZETA,WORK(KTINT), 3555 & ISYTIN,WORK(KEND2),LWRK2,IOPT) 3556 TIMZWV = TIMZWV + SECOND() - DTIME 3557 3558 DTIME = SECOND() 3559 CALL DSCAL(NT2BCD(ISYAIK),THREE*HALF,WORK(KQINT),1) 3560 TIMSCL = TIMSCL + SECOND() - DTIME 3561C 3562C write the intermediate to file: 3563C 3564 IADR = IADRPQ(IDEL) + NT2BCD(ISYAIK) 3565 3566 DTIME = SECOND() 3567 CALL PUTWA2(LUPQIM,FILNAM,WORK(KQINT),IADR,NT2BCD(ISYAIK)) 3568 TIMIO = TIMIO + SECOND() - DTIME 3569 3570 IF (LOCDBG) THEN 3571 WRITE (LUPRI,*) 'CC_PQIM> Q interm. in AO for IDEL = ',IDEL 3572 WRITE (LUPRI,'(5F12.8)') (WORK(KQINT+I),I=0,NT2BCD(ISYAIK)) 3573 END IF 3574 3575 END DO 3576 END DO 3577 3578*---------------------------------------------------------------------- 3579* save start addresses and close file and return: 3580*---------------------------------------------------------------------- 3581 CALL PUTWA2(LUPQIM,FILNAM,IADRPQ,1,(MXCORB_CC+IRAT)/IRAT) 3582 3583 CALL WCLOSE2(LUPQIM,FILNAM,'KEEP') 3584 3585 IF (IPRINT.GE.10) THEN 3586 TIMET = SECOND() - TIMET 3587 WRITE(LUPRI,*) ' Timings of CC_PQIM:' 3588 WRITE(LUPRI,1) 'I/O ', TIMIO 3589 WRITE(LUPRI,1) 'CC_TI ', TIMTI 3590 WRITE(LUPRI,1) '2C-E of T2 ', TIMTAM 3591 WRITE(LUPRI,1) '2E+C of L2 ', TIMZET 3592 WRITE(LUPRI,1) 'CC_ZWV ', TIMZWV 3593 WRITE(LUPRI,1) 'scaling etc. ', TIMSCL 3594 WRITE(LUPRI,1) 'total CC_PQIM: ', TIMET 3595 END IF 3596 3597 1 FORMAT(1x,'Time used for',2x,A18,2x,': ',f10.2,' seconds') 3598 3599 CALL QEXIT('CC_PQIM') 3600 3601 RETURN 3602 END 3603*====================================================================== 3604*====================================================================== 3605 SUBROUTINE CC_XYMIM(LLIST,ILLNR,RLIST,IRLNR,FILNAM,WORK,LWORK) 3606*---------------------------------------------------------------------- 3607* 3608* Purpose: Calculate the Y, X and M intermediates used in 3609* the left transformation and the F matrix 3610* and write them to file FILNAM 3611* 3612* Christof Haettig, June 1998 3613* 3614*====================================================================== 3615#if defined (IMPLICIT_NONE) 3616 IMPLICIT NONE 3617#else 3618# include "implicit.h" 3619#endif 3620#include "priunit.h" 3621#include "ccorb.h" 3622#include "ccsdsym.h" 3623#include "iratdef.h" 3624#include "dummy.h" 3625 3626 INTEGER LUXYM 3627 3628 INTEGER ILLNR, IRLNR, LWORK 3629 CHARACTER*(*) LLIST, RLIST, FILNAM 3630 3631#if defined (SYS_CRAY) 3632 REAL ONE, TWO, THREE, HALF 3633 REAL WORK(*) 3634#else 3635 DOUBLE PRECISION ONE, TWO, THREE, HALF 3636 DOUBLE PRECISION WORK(*) 3637#endif 3638 PARAMETER( ONE=1.0D0, TWO=2.0D0, THREE=3.0D0, HALF=0.5D0 ) 3639 3640 CHARACTER*(10) MODEL 3641 INTEGER ISYTAMP, ISYZETA, ISYXYM 3642 INTEGER KZETA, KTAMP, KYIM, KXIM, KMIM, KEND1, LWRK1, IOPT 3643 3644 3645* external function: 3646 INTEGER ILSTSYM 3647 3648 CALL QENTER('CC_XYMIM') 3649 3650*---------------------------------------------------------------------- 3651* set symmetries and allocate work space: 3652*---------------------------------------------------------------------- 3653 ISYTAMP = ILSTSYM(RLIST,IRLNR) 3654 ISYZETA = ILSTSYM(LLIST,ILLNR) 3655 ISYXYM = MULD2H(ISYTAMP,ISYZETA) 3656 3657 KZETA = 1 3658 KTAMP = KZETA + NT2SQ(ISYZETA) 3659 KYIM = KTAMP + NT2AM(ISYTAMP) 3660 KXIM = KYIM + NMATAB(ISYXYM) 3661 KMIM = KXIM + NMATIJ(ISYXYM) 3662 KEND1 = KMIM + N3ORHF(ISYXYM) 3663 LWRK1 = LWORK - KEND1 3664 3665 IF ( (LWRK1 .LT. 0) .OR. ((LWORK-KTAMP).LT.NT2AM(ISYZETA)) ) THEN 3666 WRITE (LUPRI,*) 'Insufficient memory in CC_XYMIM.' 3667 CALL QUIT('Insufficient memory in CC_XYMIM.') 3668 END IF 3669 3670*---------------------------------------------------------------------- 3671* read the zeta vector and amplitude vector into core: 3672*---------------------------------------------------------------------- 3673* first the packed zeta vector in TAMP: 3674 IOPT = 2 3675 CALL CC_RDRSP(LLIST,ILLNR,ISYZETA,IOPT,MODEL, 3676 & DUMMY,WORK(KTAMP)) 3677 3678* square up zeta vector and store now in ZETAA: 3679 CALL CC_T2SQ(WORK(KTAMP),WORK(KZETA),ISYZETA) 3680 3681* read now the amplitude vector into TAMP: 3682 IOPT = 2 3683 CALL CC_RDRSP(RLIST,IRLNR,ISYTAMP,IOPT,MODEL, 3684 & DUMMY,WORK(KTAMP)) 3685 3686* fix normalization for response vectors: 3687 IF (RLIST(1:2).NE.'R0') THEN 3688 CALL CCLR_DIASCL(WORK(KTAMP),TWO,ISYTAMP) 3689 END IF 3690 3691*---------------------------------------------------------------------- 3692* calculate the intermediates and save them on file: 3693*---------------------------------------------------------------------- 3694 3695 CALL CC_YI(WORK(KYIM),WORK(KZETA),ISYZETA,WORK(KTAMP),ISYTAMP, 3696 & WORK(KEND1),LWRK1) 3697 3698 CALL CC_XI(WORK(KXIM),WORK(KZETA),ISYZETA,WORK(KTAMP),ISYTAMP, 3699 & WORK(KEND1),LWRK1) 3700 3701 CALL CC_MI(WORK(KMIM),WORK(KZETA),ISYZETA,WORK(KTAMP),ISYTAMP, 3702 & WORK(KEND1),LWRK1) 3703 3704 LUXYM = -1 3705 CALL GPOPEN(LUXYM,FILNAM,'UNKNOWN',' ','UNFORMATTED',IDUMMY, 3706 & .FALSE.) 3707 3708 REWIND(LUXYM,ERR=999) 3709 3710 WRITE(LUXYM,ERR=999) (WORK(KYIM-1+I),I=1,NMATAB(ISYXYM)) 3711 WRITE(LUXYM,ERR=999) (WORK(KXIM-1+I),I=1,NMATIJ(ISYXYM)) 3712 WRITE(LUXYM,ERR=999) (WORK(KMIM-1+I),I=1,N3ORHF(ISYXYM)) 3713 3714 CALL GPCLOSE(LUXYM,'KEEP') 3715 3716 RETURN 3717 3718*---------------------------------------------------------------------- 3719* handle I/O errors: 3720*---------------------------------------------------------------------- 3721999 CONTINUE 3722 WRITE (LUPRI,*) 'I/O error for in CC_XYMIM:' 3723 WRITE (LUPRI,*) 'FILE:',FILNAM 3724 WRITE (LUPRI,*) 'UNIT:',LUXYM 3725 CALL QUIT('I/O error for in CC_XYMIM.') 3726 3727 CALL QEXIT('CC_XYMIM') 3728 3729 END 3730*====================================================================== 3731