1!... Copyright (c) 2015 by the authors of Dalton (see below). 2!... All Rights Reserved. 3!... 4!... The source code in this file is part of 5!... "Dalton, a molecular electronic structure program, 6!... Release DALTON2015 (2015), see http://daltonprogram.org" 7!... 8!... This source code is provided under a written licence and may be 9!... used, copied, transmitted, or stored only in accord with that 10!... written licence. 11!... 12!... In particular, no part of the source code or compiled modules may 13!... be distributed outside the research group of the licence holder. 14!... This means also that persons (e.g. post-docs) leaving the research 15!... group of the licence holder may not take any part of Dalton, 16!... including modified files, with him/her, unless that person has 17!... obtained his/her own licence. 18!... 19!... For further information, including how to get a licence, see: 20!... http://daltonprogram.org 21! 22! 23C 24*=====================================================================* 25 SUBROUTINE CCEOM_XIETA( IXETRAN, NXETRAN, IOPTRES, IORDER, LISTL, 26 & FILXI, IXDOTS, XCONS, 27 & FILETA, IEDOTS, ECONS, 28 & MXVEC, WORK, LWORK ) 29*---------------------------------------------------------------------* 30* 31* Purpose: batched loop over Xi and Eta vector calculations 32* for EOM/CCCI 33* 34* for more detailed documentation see: CC_XIETA1 35* 36* Written by Sonia and Filip. 2015, 37* Based on CC_XIETA by Christof Haettig 38* 39*=====================================================================* 40#if defined (IMPLICIT_NONE) 41 IMPLICIT NONE 42#else 43# include "implicit.h" 44#endif 45#include "priunit.h" 46#include "maxorb.h" 47#include "cclists.h" 48#include "ccsdinp.h" 49#include "ccnoddy.h" 50Cholesky 51#include "ccdeco.h" 52Cholesky 53!sonia 54#include "second.h" 55 56 LOGICAL LOCDBG 57 PARAMETER (LOCDBG=.false.) 58 59 INTEGER MAXXETRAN, MXVEC 60 PARAMETER (MAXXETRAN = 100) 61 62 CHARACTER*(*) LISTL, FILXI, FILETA 63 INTEGER IOPTRES, IORDER 64 INTEGER NXETRAN, LWORK 65 INTEGER IXETRAN(MXDIM_XEVEC,NXETRAN) 66 INTEGER IXDOTS(MXVEC,NXETRAN), IEDOTS(MXVEC,NXETRAN) 67 68#if defined (SYS_CRAY) 69 REAL WORK(LWORK) 70 REAL XCONS(MXVEC,NXETRAN), ECONS(MXVEC,NXETRAN) 71#else 72 DOUBLE PRECISION WORK(LWORK) 73 DOUBLE PRECISION XCONS(MXVEC,NXETRAN), ECONS(MXVEC,NXETRAN) 74#endif 75 76 INTEGER NTRAN, IFIRST, IBATCH, NBATCH, IDUM 77 INTEGER ISTART, IEND, KEND, LEND 78 79 CALL QENTER('CCEOM_XIETA') 80 81*---------------------------------------------------------------------* 82* Main section: CCEOM_XIETA1, driven by a loop over Xi/Eta vectors 83* ------------- 84* singles and doubles models: 85* calculation of the single and double excitation 86* parts of the Xi and Eta vectors and respective 87* dot products (IOPTRES=5). 88*---------------------------------------------------------------------* 89 90 NBATCH = (NXETRAN+MAXXETRAN-1)/MAXXETRAN 91 92 IF (LOCDBG) THEN 93 WRITE (LUPRI,*) 'Batching over Xi/Eta vector calculations:' 94 WRITE (LUPRI,*) 'nb. of batches needed:', NBATCH 95 CALL FLSHFO(LUPRI) 96 END IF 97 98 DO IBATCH = 1, NBATCH 99 IFIRST = (IBATCH-1) * MAXXETRAN + 1 100 NTRAN = MIN(NXETRAN-(IFIRST-1),MAXXETRAN) 101 102 ISTART = 1 103 IEND = ISTART + NTRAN 104 KEND = IEND + NTRAN 105 LEND = LWORK - KEND 106 if (locdbg) then 107 WRITE (LUPRI,*) 'istart,iend,kend,lend,lwork,ntran', 108 & istart,iend,kend,lend,lwork,ntran 109 end if 110 111 IF (LEND .LT. 0) THEN 112 WRITE (LUPRI,*) 'Insufficient work space in CCEOM_XIETA.' 113 WRITE (LUPRI,*) 'Available :',LWORK,' words,' 114 WRITE (LUPRI,*) 'Need at least:',KEND, ' words.' 115 CALL QUIT('Insufficient work space in CCEOM_XIETA.') 116 END IF 117 118 IF (LOCDBG) THEN 119 WRITE (LUPRI,*) 'Batch No.:',IBATCH 120 WRITE (LUPRI,*) 'start at :',IFIRST 121 WRITE (LUPRI,*) '# transf.:',NTRAN 122 WRITE (LUPRI,*) 'kend :',KEND 123 IDUM = 0 124 WRITE (LUPRI,*) 'work(0) :',WORK(IDUM) 125 END IF 126 127 CALL CCEOM_XIETA1(IXETRAN(1,IFIRST),NTRAN,IOPTRES, 128 & IORDER,LISTL,MXVEC, 129 & FILXI, IXDOTS(1,IFIRST),XCONS(1,IFIRST), 130 & FILETA,IEDOTS(1,IFIRST),ECONS(1,IFIRST), 131 & WORK(ISTART), WORK(IEND), 132 & WORK(KEND), LEND ) 133 134 IF (LOCDBG) THEN 135 WRITE (LUPRI,*) 'returned from CCEOM_XIETA1:' 136 IDUM = 0 137 WRITE (LUPRI,*) 'work(0) :',WORK(IDUM) 138 END IF 139 140 END DO 141 142 CALL QEXIT('CCEOM_XIETA') 143 144 RETURN 145 END 146 147*---------------------------------------------------------------------* 148* END OF SUBROUTINE CCEOM_XIETA * 149*---------------------------------------------------------------------* 150c /* deck ccEOM_xieta1 */ 151*=====================================================================* 152 SUBROUTINE CCEOM_XIETA1( IXETRAN, NXETRAN, IOPTRES, 153 & IORDER, LISTL, MXVEC, 154 & FILXI, IXDOTS, XCONS, 155 & FILETA, IEDOTS, ECONS, 156 & ISTART, IEND, 157 & WORK, LWORK ) 158*---------------------------------------------------------------------* 159* NODDY VERSION OF EOM XIETA 160* Purpose: Calculation of the Jacobian 161* left transformation with a generalized Hamiltonian 162* describing a (first-order) perturbation including the 163* contributions of derivative integrals, orbital relaxation 164* and reorthonormalization (orbital connection). 165* 166* Used to calculate right-hand-side vectors for: 167* -- general one-electron perturbations T response "O1" 168* -- geometrical first-derivatives T response "O1" 169* -- magnetic first-derivatives T response "O1" 170* -- general one-electron perturbations T-bar response "X1" 171* -- geometrical first-derivatives T-bar response "X1" 172* -- magnetic first-derivatives T-bar response "X1" 173* -- Cauchy expansion of one-electron T response "CO1" 174* 175* <mu|exp(-T^0) Hbar^(1) |CC> and/or <L|exp(-T^0)[Hbar^(1),mu]|CC> 176* 177* Hbar^(1) = hbar^(1)_pq E_pq + 1/2 gbar^(1)_pqrs e_pqrs 178* 179* hbar^(1)_pq = h^[1]_pq + 180* h^[0]_alp,q LQ^p_alp,p + h^[0]_p,bet LQ^h_bet,q 181* gbar_pqrs defined analogously 182* 183* IXETRAN(1,*) -- operator indices in LBLOPR array 184* IXETRAN(2,*) -- indices for left vectors for ETA result vector 185* IXETRAN(3,*) -- indices for output Xi vector on FILXI list 186* IXETRAN(4,*) -- indices for output Eta vector on FILETA list 187* IXETRAN(5,*) -- indices for the 1. orbital relaxation vectors 188* (for unrelaxed use 0) 189* IXETRAN(6,*) -- indices for the 2. orbital relaxation vectors 190* (for unrelaxed use 0, only partially implemented) 191* IXETRAN(7,*) -- indices for the 3. orbital relaxation vectors 192* (not yet used....) 193* IXETRAN(8,*) -- indices for the 4. orbital relaxation vectors 194* (not yet used....) 195* 196* NXETRAN -- number of requested transformations/vectors 197* 198* FILXI -- list type of the Xi result vectors (IOPTRES=3,4) or 199* of the vectors which are dotted on the Xi vectors 200* (IOPTRES=5) 201* FILETA -- list type of the Eta result vectors (IOPTRES=3,4) or 202* of the vectors which are dotted on the Eta vectors 203* (IOPTRES=5) 204* 205* IXCONS -- indices of vectors to be dotted on the Xi vectors 206* IECONS -- indices of vectors to be dotted on the Eta vectors 207* 208* XCONS -- contains for IOPTRES=5 the Xi dot products on return 209* ECONS -- contains for IOPTRES=5 the Eta dot products on return 210* 211* IOPTRES = 0 : all result vectors are written to direct access 212* files. FILXI and FILETA are used as file names, 213* the start addresses of the vectors are returned 214* in IXETRAN(3,*) and IXETRAN(4,*) 215* 216* IOPTRES = 3 : each result vector is written to its own file by 217* a CALL to CC_WRRSP using FILXI/FILETA as list type 218* and IXETRAN(3,*)/IXETRAN(4,*) as list index 219* 220* IOPTRES = 4 : each result vector is added to a vector on file by 221* a CALL to CC_WARSP using FILXI/FILETA as list type 222* and IXETRAN(3,*)/IXETRAN(4,*) as list index 223* 224* IOPTRES = 5 : the result vectors are dotted on an array of vectors, 225* the type of the arrays is given by FILXI for the 226* Xi result vectors and by FILETA for the Eta result 227* vectors and the indices are taken from the matrices 228* IXDOTS and IEDOTS, respectively. the resulting 229* scalar results are returned in the matrices XCONS 230* and ECONS. 231* 232* IT2DEL0,IT2DELB -- integer scratch arrays of dimensions 233* MXCORB_CC and MXCORB_CC x NXETRAN 234* 235* operator labels: 236* HAM0 : unperturbed Hamiltonian (1-el & 2-el part) 237* (for test purposes) 238* 1DHAMxxx : geometrical first derivatives (1-el & 2-el part) 239* 240* ELSE the label is intepreted as a one-electron operator 241* and the reorthonormalization terms are skipped 242* (see subroutine CC_GET_RMAT for details about the 243* treatment of the connection matrix) 244* 245* Written by Christof Haettig, May 1998 -- Jan 1999. 246* 247*=====================================================================* 248#if defined (IMPLICIT_NONE) 249 IMPLICIT NONE 250#else 251# include "implicit.h" 252#endif 253#include "priunit.h" 254#include "aovec.h" 255#include "maxorb.h" 256#include "blocks.h" 257#include "mxcent.h" 258#include "maxaqn.h" 259#include "nuclei.h" 260#include "symmet.h" 261#include "inftap.h" 262#include "ccorb.h" 263#include "ccfield.h" 264#include "ccsdsym.h" 265#include "ccsdinp.h" 266#include "iratdef.h" 267#include "distcl.h" 268#include "eribuf.h" 269#include "ccisao.h" 270#include "ccroper.h" 271#include "ccfro.h" 272#include "cco1rsp.h" 273#include "ccx1rsp.h" 274#include "ccr1rsp.h" 275#include "cclists.h" 276#include "chrxyz.h" 277#include "dummy.h" 278#include "ccsections.h" 279#include "qm3.h" 280!#include "qmmm.h" 281#include "ccexci.h" 282#include "ccnoddy.h" 283#include "cch2d.h" 284#include "ccrc1rsp.h" 285#include "r12int.h" 286#include "ccr12int.h" 287 288* local parameters: 289 LOGICAL LOCDBG 290 PARAMETER (LOCDBG = .false.) 291 INTEGER ISYM0 292 PARAMETER (ISYM0 = 1) ! symmetry of the reference state 293#if defined (SYS_CRAY) 294 REAL DDOT, ONE, TWO, THREE, HALF, ZERO 295#else 296 DOUBLE PRECISION DDOT, ONE, TWO, THREE, HALF, ZERO 297#endif 298 PARAMETER (ONE = 1.0d0, TWO = 2.0d0, THREE = 3.0d0, HALF = 0.5d0) 299 PARAMETER (ZERO = 0.0d0) 300 301 CHARACTER*10 MODEL 302 CHARACTER*(*) LISTL, FILXI, FILETA 303 CHARACTER*8 LABELH, LAB1, LABEL1, LABEL2 304 INTEGER LWORK, NXETRAN, IOPTRES, MAXSIM, MXVEC, IORDER 305 INTEGER IXETRAN(MXDIM_XEVEC,NXETRAN) 306 INTEGER IXDOTS(MXVEC,NXETRAN), IEDOTS(MXVEC,NXETRAN) 307 INTEGER ISTART(NXETRAN), IEND(NXETRAN) 308!Sonia 309 INTEGER IADRF_ETA,IADRF_XI,IATOM, IBATCH, IDLSTL,IDLSTL_OLD 310 INTEGER IFILE,IOPER,IOPER_OLD,IOPT,IOPTW,IOPTWE,IOPTWR12 311 integer IOPTYP,IRELAX,IRELAX2_OLD,IRELAX2,IRELAX_OLD,ISYCTR 312 integer ISYETA,isyhop,isyres, itran,itst 313 integer IVEC,KEND0,KEND1,KEND1SV,KEND2,KETA1,KETA2 314 logical LLrelax,lltwoel,lrelax,ltwoel,LZEROIMDONE,newtype 315 logical skipeta,skipxi 316 integer lueta,luxi,LWRK1,LWRK1SV ,LWRK2, kzeta1,kzeta2 317 integer mscratch0,mscratch1,mwork, nbatch 318 CHARACTER*10 MODELW 319 INTEGER :: ISPNETA, ISPNOP, ISPNTR 320 INTEGER :: KEND3, LEND3 321 322#if defined (SYS_CRAY) 323 REAL WORK(LWORK) 324 REAL AVERAGE, FREQ, FREQLST, KTEST 325 REAL XCONS(MXVEC,NXETRAN), ECONS(MXVEC,NXETRAN) 326#else 327 DOUBLE PRECISION WORK(LWORK) 328 DOUBLE PRECISION AVERAGE, FREQ, FREQLST, KTEST 329 DOUBLE PRECISION XCONS(MXVEC,NXETRAN), ECONS(MXVEC,NXETRAN) 330#endif 331 332* external functions: 333 INTEGER ILSTSYM 334 INTEGER IROPER 335 336 CALL QENTER('CCEOM_XIETA1') 337 338*---------------------------------------------------------------------* 339* begin: check wavefunction model and flush print unit 340*---------------------------------------------------------------------* 341 342 IF ( .NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN 343 WRITE (LUPRI,'(/1x,a)') 'CCEOM_XIETA Called for a CC ' 344 & //'method not implemented in CCEOM_XIETA...' 345 CALL QUIT('Unknown CC method in CCEOM_XIETA.') 346 END IF 347 348 CALL FLSHFO(LUPRI) 349 350 ! save the present value of the DIRECT flag 351 !DIRSAV = DIRECT 352 353 IF (LOCDBG) THEN 354 ITST = 0 355 WRITE (LUPRI,'(/1x,a,i15)') 'Work space in CCEOM_XIETA:',LWORK 356 WRITE (LUPRI,*) 'FILXI = ',FILXI(1:3) 357 WRITE (LUPRI,*) 'FILETA = ',FILETA(1:3) 358 WRITE (LUPRI,*) 'NXETRAN, MXVEC:',NXETRAN,MXVEC 359 WRITE (LUPRI,*) 'WORK(0) = ',WORK(ITST) 360 CALL FLSHFO(LUPRI) 361 END IF 362 363*---------------------------------------------------------------------* 364* check return option for the result vectors and initialize output: 365*---------------------------------------------------------------------* 366 LUXI = -1 367 LUETA = -1 368 IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN 369 CALL WOPEN2(LUXI, FILXI, 64,0) 370 CALL WOPEN2(LUETA,FILETA,64,0) 371 IADRF_XI = 1 372 IADRF_ETA = 1 373 IF (CCSDT) CALL QUIT('Problem in CCEOM_XIETA.') 374 ELSE IF (IOPTRES.EQ.2) THEN 375 CALL QUIT('IOPTRES=2 option not implemented in CC_XIETA.') 376 ELSE IF (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4) THEN 377 CONTINUE 378 ELSE IF (IOPTRES.EQ.5) THEN 379 IF (MXVEC*NXETRAN .LE. 0) THEN 380 WRITE (LUPRI,*) 381 & 'WARNING: CCEOM_XIETA called, but nothing to do!' 382 RETURN 383 END IF 384 ELSE 385 CALL QUIT('Illegal value IOPTRES in CCEOM_XIETA.') 386 END IF 387 388 IOPTWR12 = 0 389 IF (CCS) THEN 390 MODELW = 'CCS ' 391 IOPTW = 1 392 ELSE IF (CC2) THEN 393 MODELW = 'CC2 ' 394 IF (CCR12) THEN 395 MODELW = 'CC2-R12 ' 396 IOPTWR12 = 32 397 END IF 398 IOPTW = 3 399 ELSE IF (CCSD) THEN 400 MODELW = 'CCSD ' 401 IF (CCR12) THEN 402 MODELW = 'CCSD-R12 ' 403 IOPTWR12 = 32 404 END IF 405 IOPTW = 3 406 ELSE IF (CC3) THEN 407 MODELW = 'CC3 ' 408 IF (CCR12) THEN 409 MODELW = 'CC3-R12 ' 410 IOPTWR12 = 32 411CCN CALL QUIT('No CC3-R12 yet in CC_XIETA!') 412 END IF 413 IOPTW = 3 414 IOPTWE = 24 415 ELSE 416 CALL QUIT('Unknown CC model in CCEOM_XIETA.') 417 END IF 418*---------------------------------------------------------------------* 419* prepare batching, make sure that for derivatives all transf. 420* in a batch belong to coordinates of the same atom, 421* do not mix derivatives with other perturbations 422* -> IOPTYP = 0 : 1el perturbations 423 424 NBATCH = 1 425 MWORK = 0 426 IATOM = 0 427 IOPTYP = 0 428 ISTART(NBATCH) = 1 429 DO ITRAN = 1, NXETRAN 430 431 IOPER = IXETRAN(1,ITRAN) 432 LABELH = LBLOPR(IOPER) 433 434 NEWTYPE = (IOPTYP.NE.0) 435 IATOM = 0 436 IOPTYP = 0 437! IF ( (NEWTYPE .AND. ITRAN.GT.1) .OR. 438! & ((MWORK+MSCRATCH1+MSCRATCH0).GT.LWORK) ) THEN 439! IEND(NBATCH) = ITRAN-1 440! NBATCH = NBATCH + 1 441! ISTART(NBATCH) = ITRAN 442! MWORK = 0 443! END IF 444 445! MWORK = MWORK + MSCRATCH1 446! IF ((MWORK+MSCRATCH0) .GT. LWORK) THEN 447! WRITE (LUPRI,*) 'Insufficient work space in CC_XIETA.' 448! WRITE (LUPRI,*) 'Available:',LWORK,' words.' 449! WRITE (LUPRI,*) 'Need at least:',MWORK+MSCRATCH0,' words.' 450! CALL FLSHFO(LUPRI) 451! CALL QUIT('Insufficient work space in CC_XIETA.') 452! END IF 453 454 END DO 455 IEND(NBATCH) = NXETRAN 456 457 IF (LOCDBG) THEN 458 WRITE (LUPRI,*) 'CCEOM_XIETA> batching statistics:' 459 WRITE (LUPRI,*) 'CCEOM_XIETA> NBATCH : ',NBATCH 460 WRITE (LUPRI,*) 'CCEOM_XIETA> ISTART : ', 461 & (ISTART(IBATCH),IBATCH=1,NBATCH) 462 WRITE (LUPRI,*) 'CCEOM_XIETA> IEND : ', 463 & (IEND(IBATCH), IBATCH=1,NBATCH) 464 WRITE (LUPRI,*) 'WORK(0) = ',WORK(ITST) 465 END IF 466 467*=====================================================================* 468* loop over the batches of transformations: 469*=====================================================================* 470 471 472 DO IBATCH = 1, NBATCH 473 474*---------------------------------------------------------------------* 475* loop over transformations, allocate work space and 476* initialze R, G, F, BF & BFZ intermediates: 477*---------------------------------------------------------------------* 478 LLRELAX = .FALSE. 479 LLTWOEL = .FALSE. 480 LZEROIMDONE = .FALSE. 481 482 KEND0 = 1 !Sonia ??? here??? 483 KEND1 = KEND0 484 LWRK1 = LWORK - KEND1 485 486 IOPER_OLD = -1 487 IRELAX_OLD = -1 488 IDLSTL_OLD = -1 489 490 IRELAX2 = -1 491 IRELAX2_OLD = -1 492 493*---------------------------------------------------------------------* 494* start a new loop over the transformations and calculate now the 495* complete vectors from the intermediates: 496*---------------------------------------------------------------------* 497 KEND1SV = KEND1 498 LWRK1SV = LWRK1 499 500 DO ITRAN = ISTART(IBATCH), IEND(IBATCH) 501 502 IOPER = IXETRAN(1,ITRAN) ! operator index 503 IDLSTL = IXETRAN(2,ITRAN) ! ZETA vector index 504 IRELAX = IXETRAN(5,ITRAN) ! flag for relax. contrib. 505 506 ISYHOP = ISYOPR(IOPER) ! symmetry of hamiltonian 507 LABELH = LBLOPR(IOPER) ! operator label 508 LTWOEL = LPDBSOP(IOPER) ! two-electron contribution 509 ISPNOP = 1 ! Only singlet operators for now 510 511 SKIPETA = ( IXETRAN(4,ITRAN) .EQ. -1 ) 512 LRELAX = LTWOEL .OR. (IRELAX.GE.1) 513 514 SKIPXI = .true. 515C 516 IF (.NOT. SKIPETA) THEN 517 ISYCTR = ILSTSYM(LISTL,IDLSTL) ! sym. of ZETA vector 518 ISYETA = MULD2H(ISYHOP,ISYCTR) ! sym. of ETA result vector 519 IF (LISTL(1:2).EQ.'LE') THEN 520 ISPNTR = IMULTE(IDLSTL) 521 ELSE 522 ISPNTR = 1 523 ENDIF 524 ISPNETA = ISPNTR ! Only singlet operators work for now 525 END IF 526 527 528 KZETA1 = KEND1 529 KETA1 = KZETA1 + NT1AM(ISYCTR) 530 KEND2 = KETA1 + NT1AM(ISYETA) 531 IF (CC2 .OR. CCSD .OR. CCSDT) THEN 532 !write(lupri,*)' EOM_XIETA' 533 !call flshfo(lupri) 534 KETA2 = KEND2 535 KEND2 = KETA2 + NT2AM(ISYETA) 536 IF (ISPNETA.EQ.3) THEN ! Eta is triplet: twice as long! 537 KEND2 = KEND2 + NT2AM(ISYETA) 538 END IF 539 END IF 540 LWRK2 = LWORK - KEND2 541 542 ! Perform self-test: Calculate Xi and compare!! 543C IF (ISPNTR.EQ.3) THEN 544 IF (.FALSE.) THEN 545 CALL CC_XIC13(ISYMOP,LABELH,WORK(KETA1),FILETA,IDLSTL, 546 & .TRUE.,DUMMY,WORK(KEND2),LWRK2) 547C 548 KEND3 = KEND2 + MXVEC 549 LEND3 = LWRK2 - MXVEC 550 IOPT = 5 551C 552 CALL CCDOTRSP(IEDOTS,WORK(KEND2),IOPT,LISTL,ITRAN,MXVEC, 553 & NXETRAN, 554 & WORK(KETA1),WORK(KETA2),ISYCTR, 555 & WORK(KEND3),LEND3) 556 557 IVEC = 1 558 DO WHILE (IEDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 559 WRITE (LUPRI,*) 560 & ' 2 ECONS:',ITRAN,IVEC,WORK(KEND2-1+IVEC),IOPT 561 IVEC = IVEC + 1 562 END DO 563 564 END IF 565 !write(lupri,*)'SONIA: KETA1:', KETA1, ' KETA2:', KETA2 566 !write(lupri,*)'SONIA: KEND2:', KEND2 567 !call dzero(WORK(KETA1),NT1AM(ISYETA)+NT2AM(ISYETA)) 568 569 IF (LWRK2 .LT. 0) THEN 570 CALL QUIT('Insufficient work space in CC_XIETA. (CCETA2)') 571 END IF 572 573 IF (ISPNTR.EQ.3.AND.ISPNOP.EQ.1) THEN 574 575 CALL CC_ETAC13(ISYHOP,LABELH,WORK(KETA1),LISTL,IDLSTL, 576 & .TRUE., DUMMY,WORK(KEND2),LWRK2) 577 ELSE 578 call CCCI_ETAC(ISYHOP,LABELH,WORK(KETA1),LISTL,IDLSTL, 579 & 0, DUMMY,WORK(KEND2),LWRK2) 580 END IF 581 582 IF (LOCDBG) THEN 583 WRITE (LUPRI,*) 'CCEOM_XIETA> IOPER, IRELAX = ',IOPER,IRELAX 584 WRITE (LUPRI,*) 'CCEOM_XIETA> LTWOEL,LRELAX = ',LTWOEL,LRELAX 585 WRITE (LUPRI,*) 'CCEOM_XIETA> LABELH,ISYHOP = ',LABELH,ISYHOP 586 WRITE (LUPRI,*) 'CCEOM_XIETA> SKIPXI = ',SKIPXI 587 WRITE (LUPRI,*) 'WORK(0) = ',WORK(ITST) 588 CALL FLSHFO(LUPRI) 589 END IF 590 591 KEND1 = KEND1SV 592 593*=====================================================================* 594* calculate now the ETA vector: 595*=====================================================================* 596 IF (LOCDBG) THEN 597 WRITE (LUPRI,*) 'SKIPETA:',SKIPETA 598 WRITE (LUPRI,*) 'IDLSTL :',IDLSTL 599 WRITE (LUPRI,*) 'IXETRAN(4,ITRAN):',IXETRAN(2,ITRAN) 600 END IF 601 602 IF (.NOT. SKIPETA) THEN 603 604 ISYRES = ISYETA ! symmetry of result vector 605 606 IF (LWRK2 .LT. 0) THEN 607 CALL QUIT('Insufficient work space in CC_XIETA. (CCETA2)') 608 END IF 609 610 IOPT = 1 611 !useless 612 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL, 613 & WORK(KZETA1),DUMMY) 614 615 IF (IOPTRES.EQ.0 .OR. IOPTRES.EQ.1) THEN 616 IXETRAN(4,ITRAN) = IADRF_ETA 617 CALL PUTWA2(LUETA,FILETA,WORK(KETA1),IADRF_ETA, 618 & NT1AM(ISYRES)) 619 IADRF_ETA = IADRF_ETA + NT1AM(ISYRES) 620 IF (.NOT.CCS) THEN 621 CALL PUTWA2(LUETA,FILETA,WORK(KETA2),IADRF_ETA, 622 & NT2AM(ISYRES)) 623 IADRF_ETA = IADRF_ETA + NT2AM(ISYRES) 624 END IF 625 IF (CCSDT) CALL QUIT('Problem in CC_XIETA') 626 ELSE IF (IOPTRES.EQ.3) THEN 627 IFILE = IXETRAN(4,ITRAN) 628 IF (ILSTSYM(FILETA,IFILE).NE.ISYRES) THEN 629 CALL QUIT('Symmetry mismatch for Eta vector in CC_XIETA.') 630 END IF 631 CALL CC_WRRSP(FILETA,IFILE,ISYRES,IOPTW,MODELW,DUMMY, 632 & WORK(KETA1),WORK(KETA2),WORK(KEND2),LWRK2) 633 ELSE IF (IOPTRES.EQ.4) THEN 634 IFILE = IXETRAN(4,ITRAN) 635 IF (ILSTSYM(FILETA,IFILE).NE.ISYRES) THEN 636 CALL QUIT('Symmetry mismatch for Eta vector in CC_XIETA.') 637 END IF 638 CALL CC_WARSP(FILETA,IFILE,ISYRES,IOPTW,MODELW,DUMMY, 639 & WORK(KETA1),WORK(KETA2),WORK(KEND2),LWRK2) 640 ELSE IF (IOPTRES.EQ.5) THEN 641 IF (LOCDBG) THEN 642 IVEC = 1 643 WRITE(LUPRI,*) 'ECONS TRIPLES CONTRIBUTION:' 644 DO WHILE (IEDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 645 WRITE (LUPRI,*) 646 & '1 ECONS:',IVEC,ITRAN,ECONS(IVEC,ITRAN),IOPTW 647 IVEC = IVEC + 1 648 END DO 649 END IF 650C 651 !shall we do it? EOM, SONIA? 652 IF ((.NOT.CCS).AND.(ISPNETA.EQ.1)) 653 & CALL CCLR_DIASCL(WORK(KETA2),TWO,ISYRES) 654 iopt = ioptw 655 if (ispneta.eq.3) iopt = 5 656 ! 657 CALL CCDOTRSP(IEDOTS,ECONS,IOPT,FILETA,ITRAN,NXETRAN,MXVEC, 658 & WORK(KETA1),WORK(KETA2),ISYRES, 659 & WORK(KEND2),LWRK2) 660 IF (LOCDBG) THEN 661 IVEC = 1 662 DO WHILE (IEDOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC) 663 WRITE (LUPRI,*) 664 & ' 2 ECONS:',IVEC,ITRAN,ECONS(IVEC,ITRAN),IOPT 665 IVEC = IVEC + 1 666 END DO 667 END IF 668C 669 ELSE 670 CALL QUIT('Illegal value for IOPTRES in CCEOM_XIETA.') 671 END IF 672C 673 IF (LOCDBG) THEN 674 WRITE (LUPRI,*) 'Final result of CCEOM_XIETA:' 675 WRITE (LUPRI,*) 'operator label: ',LBLOPR(IOPER) 676 WRITE (LUPRI,*) 'output file type: ',FILETA 677 WRITE (LUPRI,*) 'index of output file:',IFILE 678 WRITE (LUPRI,*) 'two-electron contr.: ',LTWOEL 679 WRITE (LUPRI,*) 'relax/reorth. contr.:',LRELAX 680 I = 1 681 IF (CCS) I = 0 682 CALL CC_PRP(WORK(KETA1),WORK(KETA2),ISYRES,1,I) 683 CALL FLSHFO(LUPRI) 684 END IF 685 686 END IF 687*=====================================================================* 688* end of loop over transformations within this batch: 689*=====================================================================* 690 691 END DO ! ITRAN 692 693*---------------------------------------------------------------------* 694* close the loop over batches 695*---------------------------------------------------------------------* 696 END DO ! IBATCH 697 698 IF (IOPTRES.EQ.0) THEN 699 CALL WCLOSE2(LUXI, FILXI, 'KEEP') 700 CALL WCLOSE2(LUETA,FILETA,'KEEP') 701 ELSE IF (IOPTRES.EQ.1) THEN 702 CALL WCLOSE2(LUXI, FILXI, 'DELETE') 703 CALL WCLOSE2(LUETA,FILETA,'DELETE') 704 CALL QUIT('IOPTRES=1 not yet implemented in CCEOM_XIETA.') 705 END IF 706 707 IF (LOCDBG) THEN 708 WRITE (LUPRI,*) 'CCEOM_XIETA ended successfully (?) ... '// 709 & 'return now.' 710 CALL FLSHFO(LUPRI) 711 END IF 712 713 ! restore the DIRECT flag 714 !DIRECT = DIRSAV 715 716 CALL FLSHFO(LUPRI) 717 CALL QEXIT('CCEOM_XIETA1') 718 RETURN 719 END 720*=====================================================================* 721* END OF SUBROUTINE CCEOM_XIETA1 * 722*=====================================================================* 723 724c*DECK CCCI_ETAC 725 SUBROUTINE CCCI_ETAC(ISYMC,LBLC, !Operator stuff 726 & ETAC, !Result vector 727 & LIST,ILSTNR, !Left vector 728 & IOPTCC2, 729 & XINT,WORK,LWORK) 730C 731C----------------------------------------------------------------------- 732C Purpose: Calculate the CCCI/EOMCC etaC(L0) or 733C etaC(L) vector (Modified version of the CCSD etaC(l0/l1) code 734C 735C Important note: Requires work space of dimension of 736C NT2AM + NT2SQ in addition to ETAC, so please take care. 737C 738C eta(tau_nu)= (<HF| + Sum(mu)L(0 or 1)<mu|) 739C exp(-t)[C,tau_nu]exp(T)|HF> 740C 741C LIST= 'L0' for zeroth order left amplitudes. 742C ISYML should be ISYMOP in this case. 743C 744C 'L1' for first order left amplitudes, read in from file 745C In this case the vector is found according to its list 746C number ILSTNR. 747C 748C For L1 HF contribution is skipped. 749C 750C IOPTCC2 = 1 -- transform for CC2 the Fock matrix entering the 751C E term contribution with CMO vector instead with 752C Lambda matrices 753C 754C IOPTL0 = 1 -- The L0 vectors are passed in memory 755C 756C C property integrals read according to LBLC 757C 758C SLV98,OC: Allow for input of integrals if 759C LBLC.eq.'GIVE INT' 760C 761C Sonia & Filip, Maj 2015 762C 763C----------------------------------------------------------------------- 764C 765#include "implicit.h" 766#include "priunit.h" 767#include "dummy.h" 768#include "maxorb.h" 769#include "ccorb.h" 770#include "iratdef.h" 771#include "cclr.h" 772#include "ccexci.h" 773#include "ccsdsym.h" 774#include "ccsdio.h" 775#include "ccsdinp.h" 776C 777!Sonia & Filip: find a better place 778#include "ccsections.h" 779! 780 PARAMETER( TWO = 2.0D00, XHALF = 0.5D00 ) 781 DIMENSION ETAC(*),WORK(LWORK),XINT(*) 782 CHARACTER LBLC*(*),LIST*(*),MODEL*10 783 INTEGER IOPTCC2 784 LOGICAL FCKCON,ETRAN, LOCDBG, LEOM 785 PARAMETER( LOCDBG = .false., LEOM = .TRUE. ) 786C 787! CALL AROUND( 'IN CCCI_ETAC: Constructing ETAC(L1) vector ') 788 IF ( (IPRINT .GT. 10).or.locdbg ) THEN 789 CALL AROUND( 'IN CCCI_ETAC: Constructing ETAC(L1) vector ') 790 ENDIF 791C 792C-------------------------------- 793C find symmetry of D operator. 794C-------------------------------- 795C 796 ISYML = ILSTSYM(LIST,ILSTNR) 797C 798 ISYRES = MULD2H(ISYML,ISYMC) 799 IF (( LIST .EQ. 'L0').AND.(ISYML.NE.1)) THEN 800 CALL QUIT('Misuse of CCCI_ETAC') 801 ENDIF 802C 803 TIMEC = SECOND() 804C 805 MODEL = 'CCSD ' 806 IF (CCS) MODEL = 'CCS ' 807 IF (CC2) MODEL = 'CC2 ' 808C 809C-------------------- 810C Allocate space. 811C-------------------- 812C 813 KCTMO = 1 814 KT1AM = KCTMO + N2BST(ISYMC) 815 KLAMDP = KT1AM + NT1AM(ISYMOP) 816 KLAMDH = KLAMDP + NLAMDT 817 KEND1 = KLAMDH + NLAMDT 818C 819 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 820 KCMO = KEND1 821 KFCKHF = KCMO + NLAMDT 822 KEND1 = KFCKHF + N2BST(ISYMC) 823 END IF 824C 825 LEND1 = LWORK - KEND1 826C 827 IF ( .NOT. CCS) THEN 828 829 !if (EOMCCSD.or.CCCISD) then 830 !Integrals in MO basis 831 !highly redundant, but I need a second check 832 !write(lupri,*)'CCCI_ETAC: EOM alloc for X_mo integrals' 833 KINTIA = KEND1 834 KEND1 = KINTIA + NT1AM(ISYMC) 835 LEND1 = LWORK - KEND1 836 CALL DZERO(WORK(KintIA),NT1AM(isymc)) 837 !end if 838C 839 KL1AM = KEND1 840 KL2AM = KL1AM + NT1AM(ISYML) 841 KEND2 = KL2AM + NT2SQ(ISYML) 842 LEND2 = LWORK - KEND2 843 KT2AM = KEND2 844 KEND21= KT2AM + MAX(NT2AM(ISYML),NT2AM(1)) 845 LEND21= LWORK - KEND21 846C 847 ELSE 848C 849 KL1AM = KEND1 850 KEND2 = KL1AM + NT1AM(ISYML) 851 LEND2 = LEND1 852 KEND21= KEND1 853 LEND21= LEND1 854C 855 ENDIF 856C 857 IF (LEND21.LT. 0 ) CALL QUIT('TOO LITTLE WORKSPACE IN CCCI_ETAC') 858C 859C----------------------- 860C get T1 amplitudes. 861C----------------------- 862C 863 CALL DZERO(WORK(KT1AM),NT1AM(1)) 864 IF ( .NOT. CCS) THEN 865 IOPT = 1 866 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 867 !skod = DDOT(NT1AM(1),WORK(KT1AM),1,WORK(KT1AM),1) 868 !WRITE(LUPRI,1) 'Norm of T0_1: ',skod 869 ENDIF 870C 871 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 872 * WORK(KEND21),LEND21) 873C 874 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 875 LUSIFC = -1 876 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 877 * IDUMMY,.FALSE.) 878 REWIND(LUSIFC) 879 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 880 READ(LUSIFC) 881 READ(LUSIFC) 882 READ(LUSIFC) (WORK(KCMO+I-1),I=1,NLAMDS) 883 CALL GPCLOSE(LUSIFC,'KEEP') 884 CALL CMO_REORDER(WORK(KCMO),WORK(KEND21),LEND21) 885 END IF 886C 887C------------------------------- 888C get AO property integrals. 889C------------------------------- 890C 891 CALL DZERO(WORK(KCTMO),N2BST(ISYMC)) 892 FF = 1.0D0 893C SLV98,OC give integrals option 894 IF (LBLC.EQ.'GIVE INT') THEN 895 CALL DCOPY(N2BST(ISYMC),XINT(1),1,WORK(KCTMO),1) 896 ELSE 897 FF = 1.0D0 898 CALL CC_ONEP(WORK(KCTMO),WORK(KEND21),LEND21,FF,ISYMC,LBLC) 899 ENDIF 900C 901 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 902 CALL DCOPY(N2BST(ISYMC),WORK(KCTMO),1,WORK(KFCKHF),1) 903 CALL CC_FCKMO(WORK(KFCKHF),WORK(KCMO),WORK(KCMO), 904 * WORK(KEND21),LEND21,ISYMC,1,1) 905 END IF 906C 907C----------------------------------------------- 908C Make MO T1-transformed property integrals. 909C----------------------------------------------- 910C 911 !X_pq 912 CALL CC_FCKMO(WORK(KCTMO),WORK(KLAMDP),WORK(KLAMDH), 913 * WORK(KEND21),LEND21,ISYMC,1,1) 914C 915C---------------------------------------------------------- 916C Calculate 2Cia (stored ia) Hartree-Fock contribution. 917C---------------------------------------------------------- 918C 919 CALL DZERO(ETAC,NT1AM(ISYRES)) 920 921 DO 100 ISYMI = 1,NSYM 922C 923 ISYMA = MULD2H(ISYMI,ISYMC) 924C 925 DO 110 A = 1,NVIR(ISYMA) 926C 927 DO 120 I = 1,NRHF(ISYMI) 928C 929 KOFF1 = KINTIA + IT1AM(ISYMA,ISYMI) 930 & + NVIR(ISYMA)*(I - 1) + A-1 931 KOFF2 = KCTMO + IFCVIR(ISYMI,ISYMA) 932 * + NORB(ISYMI)*(A - 1) + I - 1 933C 934 WORK(KOFF1) = WORK(KOFF2) 935C 936 120 CONTINUE 937 110 CONTINUE 938C 939 100 CONTINUE 940C 941 IF (LIST .EQ. 'L0') THEN 942 CALL DAXPY(NT1AM(ISYRES),TWO,WORK(KINTIA),1,ETAC(1),1) 943 ENDIF 944C---------------------------------------- 945 IF ( DEBUG.or.locdbg ) THEN 946 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 947 WRITE(LUPRI,*) ' ' 948 WRITE(LUPRI,1) 'Norm of ETAC - 2Cia contribution:',ETA1 949 ENDIF 950 951 !return 952C 953C------------------------ 954C IF CCS then return. 955C------------------------ 956C 957 IF ( CCS .AND. (LIST .EQ. 'L0')) RETURN 958C 959C---------------------------------------------- 960C Read zero'th order multipliers. 961C Tbar2 are SQUARED 962C---------------------------------------------- 963C 964 IOPT = 3 965 CALL CC_RDRSP(LIST,ILSTNR,ISYML,IOPT,MODEL, 966 * WORK(KL1AM),WORK(KT2AM)) 967C skod = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 968C skod = skod+DDOT(NT2AM(ISYML),WORK(KT2AM),1,WORK(KT2AM),1) 969 !WRITE(LUPRI,1) 'Norm of Tbar (full): ',skod 970 IF (.NOT. CCS) CALL CC_T2SQ(WORK(KT2AM),WORK(KL2AM),ISYML) 971C 972C-------------------------------- 973C Put T2 (packed) amplitudes in etac2. 974C-------------------------------- 975C 976 IF (.NOT. CCS) THEN 977 !read in T2AM (packed) 978 IOPT = 2 979 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KT2AM)) 980C skod = DDOT(NT2AM(1),WORK(KT2AM),1,WORK(KT2AM),1) 981 ! WRITE(LUPRI,1) 'Norm of T2(0) (full): ',skod 982 ENDIF 983C 984C-------------------------------- 985C Make X and Y intermediates. 986C-------------------------------- 987C 988 IF (.NOT. CCS) THEN 989 KXMAT = KEND21 990 KYMAT = KXMAT + NMATIJ(ISYML) 991 KEND3 = KYMAT + NMATAB(ISYML) 992 LEND3 = LWORK - KEND3 993 IF (LEND3.LT. 0 ) 994 & CALL QUIT(' TOO LITTLE WORKSPACE IN CC_ETAC-2') 995C 996 IF ( DEBUG.or.LOCDBG ) THEN 997 XYI = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 998 WRITE(LUPRI,1) 'CC_ETAC: L1AM vector: ',XYI 999 XYI = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 1000 WRITE(LUPRI,1) 'CC_ETAC: L2AM vector: ',XYI 1001 XXI = DDOT(NT2AM(ISYMOP),WORK(KT2AM),1,WORK(KT2AM),1) 1002 WRITE(LUPRI,1) 'T2AM vector : ',XXI 1003 ENDIF 1004 CALL CC_XI(WORK(KXMAT),WORK(KL2AM),ISYML,WORK(KT2AM),1, 1005 * WORK(KEND3),LEND3) 1006 CALL CC_YI(WORK(KYMAT),WORK(KL2AM),ISYML,WORK(KT2AM),1, 1007 * WORK(KEND3),LEND3) 1008 IF ( DEBUG.or.LOCDBG ) THEN 1009 XYI = DDOT(NMATAB(ISYML),WORK(KYMAT),1,WORK(KYMAT),1) 1010 WRITE(LUPRI,1) 'CC_ETAC: YI intermediate is: ',XYI 1011 XXI = DDOT(NMATIJ(ISYML),WORK(KXMAT),1,WORK(KXMAT),1) 1012 WRITE(LUPRI,1) 'CC_ETAC: XI intermediate is: ',XXI 1013 ENDIF 1014 ELSE 1015 KEND3 = KEND2 1016 LEND3 = LEND2 1017 ENDIF 1018C 1019C---------------------------------------------- 1020C Calculate X and Y contributions to etac1. 1021C etac1 = -sum(e)Cie*Yae - sum(l)Cla*Xli 1022C---------------------------------------------- 1023C 1024 IF ( (.NOT.CCS) .AND. (.NOT.(CC2.AND.IOPTCC2.EQ.1)) ) THEN 1025C DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1026 !WRITE(LUPRI,*) 'Integral norm before CC_21EFM out norm:',DN 1027 CALL CC_21EFM(ETAC,WORK(KCTMO),ISYMC,WORK(KXMAT), 1028 * WORK(KYMAT),ISYML,WORK(KEND3),LEND3) 1029C DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1030 !WRITE(LUPRI,*) 'Integral norm after CC_21EFM out norm:',DN 1031C 1032 IF ( DEBUG.or.locdbg ) THEN 1033 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1034 WRITE(LUPRI,1) 'Norm of eta1-after X&Y cont: ',ETA1 1035 ENDIF 1036 ENDIF 1037 1038 IF (LEOM) THEN 1039C 1040C--------------------------------------------------------------- 1041C EOM contribution to ETAC: 1042C ~ ~ 1043C etac = 2sum(ei) L_{ckei}*(C_{ei} + t_{ei,fn}*C_{nf}) 1044C--------------------------------------------------------------- 1045 KCINT = KEND21 1046 KETAC = KCINT 1047 KINTAI = KCINT + MAX(NT1AM(ISYMC),NT1AM(ISYRES)) 1048C Extract C_{ai} 1049 DO ISYMI = 1, NSYM 1050 ISYMA = MULD2H(ISYMI,ISYMC) 1051 DO I = 1, NRHF(ISYMI) 1052 KOFF1 = KINTAI + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) 1053 KOFF2 = KCTMO + IFCRHF(ISYMA,ISYMI) + 1054 & NORB(ISYMA)*(I-1) + NRHF(ISYMA) 1055 CALL DCOPY(NVIR(ISYMA),WORK(KOFF2),1,WORK(KOFF1),1) 1056 END DO 1057 END DO 1058 !get ttilde (overwrites T2am) 1059 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,1,1) 1060C ~ 1061C Add t_{em,fn} *C_{nf} to C_{em} 1062 CALL CCG_LXD(WORK(KCINT),ISYMC,WORK(KINTIA),ISYMC, 1063 & WORK(KT2AM),1,0) 1064 CALL DAXPY(NT1AM(ISYMC),1.D0,WORK(KCINT),1,WORK(KINTAI),1) 1065C Calculate the term as L_{ai,em} * C'_{em} 1066 CALL CCG_LXD(WORK(KETAC),ISYRES,WORK(KINTAI),ISYMC, 1067 & WORK(KL2AM),ISYML,1) 1068 IF ( DEBUG.or.locdbg ) THEN 1069 ETA1 = DDOT(NT1AM(ISYRES),WORK(KETAC),1,WORK(KETAC),1) 1070 WRITE(LUPRI,1) 'Norm of alone Tbar_ck,ei*Cei: ',ETA1 1071 ENDIF 1072 !removed factor 2 to get FCI limit! 1073 CALL DAXPY(NT1AM(ISYRES),1.0D0,WORK(KETAC),1,ETAC(1),1) 1074 END IF 1075C 1076C------------------------------------------------ 1077C Workspace for T2AM and X and Y is now free. 1078C etac2 = P(ab,ij)(2l(ai)*Cjb - l(aj)*c(ib)) 1079C------------------------------------------------ 1080C 1081 IF (.NOT. CCS) THEN 1082 CALL DZERO(ETAC(1+NT1AM(ISYRES)),NT2AM(ISYRES)) 1083 1084C DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1085 !WRITE(LUPRI,*) 'Integral norm before CC_L1FCK out norm:',DN 1086 1087 CALL CC_L1FCK(ETAC(1+NT1AM(ISYRES)),WORK(KL1AM),WORK(KCTMO), 1088 * ISYML,ISYMC,WORK(KEND2),LEND2) 1089 1090C DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1091! WRITE(LUPRI,*) 'Integral norm after CC_L1FCK out norm:',DN 1092! * ISYML,ISYMC,WORK(KEND2),LEND2) 1093C 1094 IF ( DEBUG.or.locdbg ) THEN 1095 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1096 ETA2 = DDOT(NT2AM(ISYRES),ETAC(1+NT1AM(ISYRES)),1, 1097 * ETAC(1+NT1AM(ISYRES)),1) 1098 WRITE(LUPRI,1) 'Norm of eta1-after L1c cont: ',ETA1 1099 WRITE(LUPRI,1) 'Norm of eta2-after L1c cont: ',ETA2 1100 ENDIF 1101 ENDIF 1102C 1103 KEI1 = KEND2 1104 KEI2 = KEI1 + NEMAT1(ISYMC) 1105 KEND3 = KEI2 + NMATIJ(ISYMC) 1106 LEND3 = LWORK - KEND3 1107 IF (LEND3.LT. 0 ) CALL QUIT(' TOO LITTLE WORKSPACE IN CC_ETAC-3') 1108C 1109C-------------------------------- 1110C Put A into E matrix format. 1111C-------------------------------- 1112C 1113 FCKCON = .TRUE. 1114 ETRAN = .FALSE. 1115 CALL CCRHS_EFCK(WORK(KEI1),WORK(KEI2),WORK(KLAMDH), 1116 * WORK(KCTMO),WORK(KEND3),LEND3,FCKCON, 1117 * ETRAN,ISYMC) 1118C 1119C-------------------------------------------- 1120C etac1 = sum(b)Lbi*Cba - sum(j)Laj*Cij. 1121C-------------------------------------------- 1122C 1123 IF ( DEBUG.or.locdbg ) THEN 1124 XE1 = DDOT(NMATAB(ISYMC),WORK(KEI1),1,WORK(KEI1),1) 1125 XE2 = DDOT(NMATIJ(ISYMC),WORK(KEI2),1,WORK(KEI2),1) 1126 WRITE(LUPRI,1) 'Norm of EI1 -after EFCK: ',XE1 1127 WRITE(LUPRI,1) 'Norm of EI2 -after EFCK: ',XE2 1128 XE2 = DDOT(NMATIJ(ISYMC),WORK(KINTIJ),1,WORK(KINTIJ),1) 1129 WRITE(LUPRI,1) 'Norm of Cij integrals -after EFCK: ',XE2 1130 ETA1 = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 1131 WRITE(LUPRI,1) 'Norm of L1AM before CCLR_E1C1: ',ETA1 1132 ENDIF 1133C 1134 CALL CCLR_E1C1(ETAC,WORK(KL1AM),WORK(KEI1),WORK(KEI2), 1135 * WORK(KEND3),LEND3,ISYML,ISYMC,'T') 1136C 1137 IF ( DEBUG.or.locdbg ) THEN 1138 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1139 WRITE(LUPRI,1) 'Norm of eta1 - after CCLR_E1C1: ',ETA1 1140 ENDIF 1141C 1142C--------------------------------------------------------------- 1143C etac2 = P(ab,ij)(sum(e)2L(aiej)*Ceb - sym(k)L(aibk)*C(jk)) 1144C--------------------------------------------------------------- 1145C 1146 IF (.NOT. CCS) THEN 1147C 1148 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 1149 FCKCON = .TRUE. 1150 ETRAN = .FALSE. 1151 CALL CCRHS_EFCK(WORK(KEI1),WORK(KEI2),WORK(KCMO), 1152 * WORK(KFCKHF),WORK(KEND3),LEND3,FCKCON,ETRAN,ISYMC) 1153 END IF 1154 1155 CALL CC_EITR(WORK(KEI1),WORK(KEI2),WORK(KEND3),LEND3, 1156 * ISYMC) 1157C 1158 CALL CCRHS_E(ETAC(1+NT1AM(ISYRES)),WORK(KL2AM), 1159 * WORK(KEI1),WORK(KEI2),WORK(KEND3), 1160 * LEND3,ISYML,ISYMC) 1161C 1162 IF (IPRINT .GT. 40 ) THEN 1163 CALL AROUND( 'In CCCI_ETAC: EtaC vector ' ) 1164 CALL CC_PRP(ETAC(1),ETAC(1+NT1AM(ISYRES)),ISYMC,1,1) 1165 ENDIF 1166C 1167 IF (DEBUG .OR. ( IPRINT .GT. 20 )) THEN 1168 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1169 ETA2 = DDOT(NT2AM(ISYRES),ETAC(1+NT1AM(ISYRES)),1, 1170 * ETAC(1+NT1AM(ISYRES)),1) 1171 WRITE(LUPRI,1) 'Norm of eta1 - end of CCCI_ETAC: ',ETA1 1172 WRITE(LUPRI,1) 'Norm of eta2 - end of CCCI_ETAC: ',ETA2 1173 CALL AROUND( 'END OF CC_ETAC ') 1174 ENDIF 1175 ENDIF 1176 1177 !end if 1178 IF (IPRINT .GT. 5 ) THEN 1179 TIMEC = SECOND() - TIMEC 1180 WRITE(LUPRI,9999) 'CCCI_ETA^C ', TIMEC 1181 ENDIF 1182C 1183 1 FORMAT(1x,A35,1X,E20.10) 11849999 FORMAT(1x,'Time used in',2x,A18,2x,': ',f10.2,' seconds') 1185C 1186 END 1187c*DECK CCCI_ETAC2 1188 SUBROUTINE CCCI_ETAC2(ISYMC,LBLC, !Operator stuff 1189 & ETAC, !Result vector 1190 & LIST,ILSTNR, !Left vector 1191 & IOPTCC2, 1192 & ioptl0,xl0, 1193 & XINT,WORK,LWORK) 1194C 1195C----------------------------------------------------------------------- 1196C Purpose: Calculate the CCCI/EOMCC etaC(L0) or 1197C etaC(L) vector (Modified version of the CCSD etaC(l0/l1) code 1198C 1199C Important note: Requires work space of dimension of 1200C NT2AM + NT2SQ in addition to ETAC, so please take care. 1201C 1202C eta(tau_nu)= (<HF| + Sum(mu)L(0 or 1)<mu|) 1203C exp(-t)[C,tau_nu]exp(T)|HF> 1204C 1205C LIST= 'L0' for zeroth order left amplitudes. 1206C ISYML should be ISYMOP in this case. 1207C 1208C 'L1' for first order left amplitudes, read in from file 1209C In this case the vector is found according to its list 1210C number ILSTNR. 1211C 1212C For L1 HF contribution is skipped. 1213C 1214C IOPTCC2 = 1 -- transform for CC2 the Fock matrix entering the 1215C E term contribution with CMO vector instead with 1216C Lambda matrices 1217C 1218C IOPTL0 = 1 -- The L0 vectors are passed in memory and Hartree-Fock part skipped. 1219C IOPTL0 = 2 -- only Hartree-Fock part calculated. 1220C 1221C C property integrals read according to LBLC 1222C 1223C SLV98,OC: Allow for input of integrals if 1224C LBLC.eq.'GIVE INT' 1225C 1226C Sonia & Filip, Maj 2015 1227C 1228C----------------------------------------------------------------------- 1229C 1230#include "implicit.h" 1231#include "priunit.h" 1232#include "dummy.h" 1233#include "maxorb.h" 1234#include "ccorb.h" 1235#include "iratdef.h" 1236#include "cclr.h" 1237#include "ccexci.h" 1238#include "ccsdsym.h" 1239#include "ccsdio.h" 1240#include "ccsdinp.h" 1241C 1242!Sonia & Filip: find a better place 1243#include "ccsections.h" 1244#include "fdeom.h" 1245! 1246 PARAMETER( TWO = 2.0D00, XHALF = 0.5D00 ) 1247 DIMENSION ETAC(*),WORK(LWORK),XINT(*), Xl0(*) 1248 CHARACTER LBLC*(*),LIST*(*),MODEL*10 1249 INTEGER IOPTCC2, ioptl0 1250 LOGICAL FCKCON,ETRAN, LOCDBG 1251 PARAMETER( LOCDBG = .true. ) 1252C 1253 IF ( (IPRINT .GT. 10).or.locdbg ) THEN 1254 CALL AROUND( 'IN CCCI_ETAC2: Constructing ETAC(L1) vector ') 1255 ENDIF 1256C 1257 if ( (ioptl0.ne.1) .and. (ioptl0.ne.2) ) then 1258 call quit('Unknown value of ioptl0 in CCCI_ETAC2') 1259 end if 1260C 1261 if ( (ioptl0.eq.2) .and. (LIST(1:2).ne.'L0') ) then 1262 call quit('Wrong setup in CCCI_ETAC2') 1263 end if 1264C 1265C-------------------------------- 1266C find symmetry of D operator. 1267C-------------------------------- 1268C 1269 ISYML = ILSTSYM(LIST,ILSTNR) 1270C 1271 ISYRES = MULD2H(ISYML,ISYMC) 1272 IF (( LIST .EQ. 'L0').AND.(ISYML.NE.1)) THEN 1273 CALL QUIT('Misuse of CCCI_ETAC') 1274 ENDIF 1275C 1276 TIMEC = SECOND() 1277C 1278 MODEL = 'CCSD ' 1279 IF (CCS) MODEL = 'CCS ' 1280 IF (CC2) MODEL = 'CC2 ' 1281C 1282C-------------------- 1283C Allocate space. 1284C-------------------- 1285C 1286 KCTMO = 1 1287 KT1AM = KCTMO + N2BST(ISYMC) 1288 KLAMDP = KT1AM + NT1AM(ISYMOP) 1289 KLAMDH = KLAMDP + NLAMDT 1290 KEND1 = KLAMDH + NLAMDT 1291C 1292 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 1293 KCMO = KEND1 1294 KFCKHF = KCMO + NLAMDT 1295 KEND1 = KFCKHF + N2BST(ISYMC) 1296 END IF 1297C 1298 LEND1 = LWORK - KEND1 1299C 1300 IF ( .NOT. CCS) THEN 1301 1302 !if (EOMCCSD.or.CCCISD.or.FDEOM) then 1303 !Integrals in MO basis 1304 !highly redundant, but I need a second check 1305 write(lupri,*)'CCCI_ETAC2: EOM alloc for X_mo integrals' 1306 KINTAI = KEND1 1307 KINTAB = KINTAI + NT1AM(ISYMC) 1308 KINTIJ = KINTAB + NMATAB(ISYMC) 1309 KINTIA = KINTIJ + NMATIJ(ISYMC) 1310 KEND1 = KINTIA + NT1AM(ISYMC) 1311 LEND1 = LWORK - KEND1 1312 CALL DZERO(WORK(KINTAB),NMATAB(isymc)) 1313 CALL DZERO(WORK(KintIJ),NMATIJ(isymc)) 1314 CALL DZERO(WORK(KintIA),NT1AM(isymc)) 1315 CALL DZERO(WORK(KintAI),NT1AM(isymc)) 1316 !end if 1317C 1318 KL1AM = KEND1 1319 KL2AM = KL1AM + NT1AM(ISYML) 1320 KEND2 = KL2AM + NT2SQ(ISYML) 1321 LEND2 = LWORK - KEND2 1322 KT2AM = KEND2 1323 KEND21= KT2AM + MAX(NT2AM(ISYML),NT2AM(1)) 1324 LEND21= LWORK - KEND21 1325C 1326 ELSE 1327C 1328 KL1AM = KEND1 1329 KEND2 = KL1AM + NT1AM(ISYML) 1330 LEND2 = LEND1 1331 KEND21= KEND1 1332 LEND21= LEND1 1333C 1334 ENDIF 1335C 1336 IF (LEND21.LT. 0 ) CALL QUIT('TOO LITTLE WORKSPACE IN CCCI_ETAC2') 1337C 1338C----------------------- 1339C get T1 amplitudes. 1340C----------------------- 1341C 1342 CALL DZERO(WORK(KT1AM),NT1AM(1)) 1343 IF ( .NOT. CCS) THEN 1344 IOPT = 1 1345 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 1346 skod = DDOT(NT1AM(1),WORK(KT1AM),1,WORK(KT1AM),1) 1347 WRITE(LUPRI,1) 'Norm of T0_1: ',skod 1348 ENDIF 1349C 1350 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 1351 * WORK(KEND21),LEND21) 1352C 1353 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 1354 LUSIFC = -1 1355 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 1356 * IDUMMY,.FALSE.) 1357 REWIND(LUSIFC) 1358 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1359 READ(LUSIFC) 1360 READ(LUSIFC) 1361 READ(LUSIFC) (WORK(KCMO+I-1),I=1,NLAMDS) 1362 CALL GPCLOSE(LUSIFC,'KEEP') 1363 CALL CMO_REORDER(WORK(KCMO),WORK(KEND21),LEND21) 1364 END IF 1365C 1366C------------------------------- 1367C get AO property integrals. 1368C------------------------------- 1369C 1370 CALL DZERO(WORK(KCTMO),N2BST(ISYMC)) 1371 FF = 1.0D0 1372C SLV98,OC give integrals option 1373 IF (LBLC.EQ.'GIVE INT') THEN 1374 CALL DCOPY(N2BST(ISYMC),XINT(1),1,WORK(KCTMO),1) 1375 ELSE 1376 FF = 1.0D0 1377 CALL CC_ONEP(WORK(KCTMO),WORK(KEND21),LEND21,FF,ISYMC,LBLC) 1378 ENDIF 1379C 1380 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 1381 CALL DCOPY(N2BST(ISYMC),WORK(KCTMO),1,WORK(KFCKHF),1) 1382 CALL CC_FCKMO(WORK(KFCKHF),WORK(KCMO),WORK(KCMO), 1383 * WORK(KEND21),LEND21,ISYMC,1,1) 1384 END IF 1385C 1386C----------------------------------------------- 1387C Make MO T1-transformed property integrals. 1388C----------------------------------------------- 1389C 1390 !if (eomccsd.or.cccisd.or.fdeom) then 1391 1392 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1393 WRITE(LUPRI,*) 'IN ONEP: FOCK out norm:',DN 1394 1395 write(lupri,*) 'X_mo integrals of EOM from CCDINTMO' 1396 ISYM = ISYMC 1397 CALL CCDINTMO(WORK(KINTIJ),WORK(KINTIA),WORK(KINTAB), 1398 * WORK(KINTAI),WORK(KCTMO),WORK(KLAMDP), 1399 * WORK(KLAMDH),WORK(KEND21),LEND21,ISYM) 1400 IF ( DEBUG.or.locdbg ) THEN 1401 ETA1 = DDOT(NT1AM(ISYMC),WORK(KINTIA),1,WORK(KINTIA),1) 1402 WRITE(LUPRI,*) ' ' 1403 WRITE(LUPRI,1) '4*Norm of C_ie for EOM:',ETA1*4.0d0 1404 ENDIF 1405 IF ((IPRINT .GT. 50).or.locdbg) THEN 1406C 1407 CALL AROUND('One electron integrals in MO-basis') 1408 DO 10 ISYM = 1,NSYM 1409 WRITE(LUPRI,333) 'Symmetry block number:', ISYM 1410 333 FORMAT(/3X,A22,2X,I1) 1411 KOFFIJ = KINTIJ + IMATIJ(ISYM,ISYM) 1412 KOFFAI = KINTAI + IT1AM(ISYM,ISYM) 1413 KOFFIA = KINTIA + IT1AM(ISYM,ISYM) 1414 KOFFAB = KINTAB + IMATAB(ISYM,ISYM) 1415! CALL AROUND('occ.occ part') 1416! CALL OUTPUT(WORK(KOFFIJ),1,NRHF(ISYM),1,NRHF(ISYM), 1417! * NRHF(ISYM),NRHF(ISYM),1,LUPRI) 1418! CALL AROUND('vir.occ part') 1419! CALL OUTPUT(WORK(KOFFAI),1,NVIR(ISYM),1,NRHF(ISYM), 1420! * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 1421! CALL AROUND('occ.vir part') 1422! CALL OUTPUT(WORK(KOFFIA),1,NVIR(ISYM),1,NRHF(ISYM), 1423! * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 1424! CALL AROUND('vir.vir part') 1425! CALL OUTPUT(WORK(KOFFAB),1,NVIR(ISYM),1,NVIR(ISYM), 1426! * NVIR(ISYM),NVIR(ISYM),1,LUPRI) 1427C 1428 10 CONTINUE 1429 HIJNO = DDOT(NMATIJ(isymc),WORK(KINTIJ),1,WORK(KINTIJ),1) 1430 HAINO = DDOT(NT1AM(isymc),WORK(KINTAI),1,WORK(KINTAI),1) 1431 HIANO = DDOT(NT1AM(isymc),WORK(KINTIA),1,WORK(KINTIA),1) 1432 HABNO = DDOT(NMATAB(isymc),WORK(KINTAB),1,WORK(KINTAB),1) 1433 WRITE(LUPRI,*) ' ' 1434 WRITE(LUPRI,*) 'Norm of one electron integrals in MO-basis:' 1435 WRITE(LUPRI,*) HIJNO, HAINO, HIANO, HABNO 1436 WRITE(LUPRI,*) 'Integrals sum: ', HIJNO+HAINO+HIANO+HABNO 1437 ENDIF 1438 !end if 1439 1440 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1441 WRITE(LUPRI,*) 'IN ONEP: FOCK out norm:',DN 1442 1443 !X_pq 1444 CALL CC_FCKMO(WORK(KCTMO),WORK(KLAMDP),WORK(KLAMDH), 1445 * WORK(KEND21),LEND21,ISYMC,1,1) 1446 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1447 WRITE(LUPRI,*) 'Integral norm after FCKMO out norm:',DN 1448C 1449C---------------------------------------------------------- 1450C Calculate 2Cia (stored ia) Hartree-Fock contribution. 1451C---------------------------------------------------------- 1452C 1453 CALL DZERO(ETAC,NT1AM(ISYRES)) 1454C 1455 IF ((LIST .EQ. 'L0').and.(ioptl0.ne.1)) THEN 1456 DO 100 ISYMI = 1,NSYM 1457C 1458 ISYMA = MULD2H(ISYMI,ISYMC) 1459C 1460 DO 110 A = 1,NVIR(ISYMA) 1461C 1462 DO 120 I = 1,NRHF(ISYMI) 1463C 1464 KOFF1 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 1465 KOFF2 = KCTMO + IFCVIR(ISYMI,ISYMA) 1466 * + NORB(ISYMI)*(A - 1) + I - 1 1467C 1468 ETAC(KOFF1) = TWO*WORK(KOFF2) 1469C 1470 120 CONTINUE 1471 110 CONTINUE 1472 100 CONTINUE 1473C 1474 ENDIF 1475! 1476! For ioptl0 = 2 we want only the Hartree-Fock contribution (which has just been 1477! (calculated above), so skip the rest: 1478 if (ioptl0.eq.2) then 1479 go to 8888 1480 end if 1481C---------------------------------------- 1482 IF ( DEBUG.or.locdbg ) THEN 1483 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1484 WRITE(LUPRI,*) ' ' 1485 WRITE(LUPRI,1) 'Norm of ETAC - 2Cia contribution:',ETA1 1486 ENDIF 1487 1488 !return 1489C 1490C------------------------ 1491C IF CCS then return. 1492C------------------------ 1493C 1494 IF ( CCS .AND. (LIST .EQ. 'L0')) RETURN 1495C 1496C---------------------------------------------- 1497C Read zero'th order multipliers. 1498C Tbar2 are SQUARED 1499C---------------------------------------------- 1500C 1501 if (ioptl0.eq.1) then 1502 call dcopy(NT1AM(ISYML),XL0(1),1,WORK(KL1AM),1) 1503 call dcopy(NT2AM(ISYML),XL0(1+NT1AM(ISYML)),1,WORK(KT2AM),1) 1504 else 1505 IOPT = 3 1506 CALL CC_RDRSP(LIST,ILSTNR,ISYML,IOPT,MODEL, 1507 * WORK(KL1AM),WORK(KT2AM)) 1508 end if 1509 skod = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 1510 skod = skod+DDOT(NT2AM(ISYML),WORK(KT2AM),1,WORK(KT2AM),1) 1511 WRITE(LUPRI,1) 'Norm of Tbar (full): ',skod 1512 IF (.NOT. CCS) CALL CC_T2SQ(WORK(KT2AM),WORK(KL2AM),ISYML) 1513C 1514C-------------------------------- 1515C Put T2 (packed) amplitudes in etac2. 1516C-------------------------------- 1517C 1518 IF (.NOT. CCS) THEN 1519 !read in T2AM (packed) 1520 IOPT = 2 1521 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KT2AM)) 1522 skod = DDOT(NT2AM(1),WORK(KT2AM),1,WORK(KT2AM),1) 1523 WRITE(LUPRI,1) 'Norm of T2(0) (full): ',skod 1524 ENDIF 1525C 1526C-------------------------------- 1527C Make X and Y intermediates. 1528C-------------------------------- 1529C 1530 IF (.NOT. CCS) THEN 1531 KXMAT = KEND21 1532 KYMAT = KXMAT + NMATIJ(ISYML) 1533 KEND3 = KYMAT + NMATAB(ISYML) 1534 LEND3 = LWORK - KEND3 1535 IF (LEND3.LT. 0 ) 1536 & CALL QUIT(' TOO LITTLE WORKSPACE IN CC_ETAC-2') 1537C 1538 !IF ( DEBUG ) THEN 1539 XYI = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 1540 WRITE(LUPRI,1) 'CC_ETAC: L1AM vector: ',XYI 1541 XYI = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 1542 WRITE(LUPRI,1) 'CC_ETAC: L2AM vector: ',XYI 1543 XXI = DDOT(NT2AM(ISYMOP),WORK(KT2AM),1,WORK(KT2AM),1) 1544 WRITE(LUPRI,1) 'T2AM vector : ',XXI 1545 !ENDIF 1546 CALL CC_XI(WORK(KXMAT),WORK(KL2AM),ISYML,WORK(KT2AM),1, 1547 * WORK(KEND3),LEND3) 1548 CALL CC_YI(WORK(KYMAT),WORK(KL2AM),ISYML,WORK(KT2AM),1, 1549 * WORK(KEND3),LEND3) 1550 !IF ( DEBUG ) THEN 1551 XYI = DDOT(NMATAB(ISYML),WORK(KYMAT),1,WORK(KYMAT),1) 1552 WRITE(LUPRI,1) 'CC_ETAC: YI intermediate is: ',XYI 1553 XXI = DDOT(NMATIJ(ISYML),WORK(KXMAT),1,WORK(KXMAT),1) 1554 WRITE(LUPRI,1) 'CC_ETAC: XI intermediate is: ',XXI 1555 !ENDIF 1556 ELSE 1557 KEND3 = KEND2 1558 LEND3 = LEND2 1559 ENDIF 1560C 1561C---------------------------------------------- 1562C Calculate X and Y contributions to etac1. 1563C etac1 = -sum(e)Cie*Yae - sum(l)Cla*Xli 1564C---------------------------------------------- 1565C 1566 IF ( (.NOT.CCS) .AND. (.NOT.(CC2.AND.IOPTCC2.EQ.1)) ) THEN 1567 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1568 WRITE(LUPRI,*) 'Integral norm before CC_21EFM out norm:',DN 1569 CALL CC_21EFM(ETAC,WORK(KCTMO),ISYMC,WORK(KXMAT), 1570 * WORK(KYMAT),ISYML,WORK(KEND3),LEND3) 1571 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1572 WRITE(LUPRI,*) 'Integral norm after CC_21EFM out norm:',DN 1573C 1574 IF ( DEBUG.or.locdbg ) THEN 1575 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1576 WRITE(LUPRI,1) 'Norm of eta1-after X&Y cont: ',ETA1 1577 ENDIF 1578 ENDIF 1579 1580C 1581C------------------------------------------------ 1582C Workspace for T2AM and X and Y is now free. 1583C etac2 = P(ab,ij)(2l(ai)*Cjb - l(aj)*c(ib)) 1584C------------------------------------------------ 1585C 1586 IF (.NOT. CCS) THEN 1587 CALL DZERO(ETAC(1+NT1AM(ISYRES)),NT2AM(ISYRES)) 1588 1589 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1590 WRITE(LUPRI,*) 'Integral norm before CC_L1FCK out norm:',DN 1591 1592 CALL CC_L1FCK(ETAC(1+NT1AM(ISYRES)),WORK(KL1AM),WORK(KCTMO), 1593 * ISYML,ISYMC,WORK(KEND2),LEND2) 1594 1595 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1596 WRITE(LUPRI,*) 'Integral norm after CC_L1FCK out norm:',DN 1597! * ISYML,ISYMC,WORK(KEND2),LEND2) 1598C 1599 IF ( DEBUG.or.locdbg ) THEN 1600 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1601 ETA2 = DDOT(NT2AM(ISYRES),ETAC(1+NT1AM(ISYRES)),1, 1602 * ETAC(1+NT1AM(ISYRES)),1) 1603 WRITE(LUPRI,1) 'Norm of eta1-after L1c cont: ',ETA1 1604 WRITE(LUPRI,1) 'Norm of eta2-after L1c cont: ',ETA2 1605 ENDIF 1606 ENDIF 1607C 1608 KEI1 = KEND2 1609 KEI2 = KEI1 + NEMAT1(ISYMC) 1610 KEND3 = KEI2 + NMATIJ(ISYMC) 1611 LEND3 = LWORK - KEND3 1612 IF (LEND3.LT. 0 ) CALL QUIT(' TOO LITTLE WORKSPACE IN CC_ETAC-3') 1613C 1614C-------------------------------- 1615C Put A into E matrix format. 1616C-------------------------------- 1617C 1618 FCKCON = .TRUE. 1619 ETRAN = .FALSE. 1620 DN = DDOT(N2BST(ISYMC),WORK(KCTMO),1,WORK(KCTMO),1) 1621 WRITE(LUPRI,*) 'Integral norm before CCRHS_EFCK out norm:',DN 1622 CALL CCRHS_EFCK(WORK(KEI1),WORK(KEI2),WORK(KLAMDH), 1623 * WORK(KCTMO),WORK(KEND3),LEND3,FCKCON, 1624 * ETRAN,ISYMC) 1625C 1626C-------------------------------------------- 1627C etac1 = sum(b)Lbi*Cba - sum(j)Laj*Cij. 1628C-------------------------------------------- 1629C 1630 IF ( DEBUG.or.locdbg ) THEN 1631 XE1 = DDOT(NMATAB(ISYMC),WORK(KEI1),1,WORK(KEI1),1) 1632 XE2 = DDOT(NMATIJ(ISYMC),WORK(KEI2),1,WORK(KEI2),1) 1633 WRITE(LUPRI,1) 'Norm of EI1 -after EFCK: ',XE1 1634 WRITE(LUPRI,1) 'Norm of EI2 -after EFCK: ',XE2 1635 XE2 = DDOT(NMATIJ(ISYMC),WORK(KINTIJ),1,WORK(KINTIJ),1) 1636 WRITE(LUPRI,1) 'Norm of Cij integrals -after EFCK: ',XE2 1637 ETA1 = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 1638 WRITE(LUPRI,1) 'Norm of L1AM before CCLR_E1C1: ',ETA1 1639 ENDIF 1640C 1641c test 1642c kei11= kend3 1643c kei21= kei11+ NMATAB(ISYMC) 1644c kend3 = kei21+ NMATIJ(ISYMC) 1645c lend3 = lwork -kend3 1646c call dzero(work(kei11),NMATAB(ISYMC)) 1647c call dzero(work(kei21),NMATIJ(ISYMC)) 1648c call dcopy(NMATAB(ISYMC),work(kei1),1,work(kei11),1) 1649c call dcopy(NMATIJ(ISYMC),work(kei2),1,work(kei21),1) 1650c CALL CCLR_E1C1(ETAC,WORK(KL1AM),WORK(KEI11),WORK(KEI21), 1651c * WORK(KEND3),LEND3,ISYML,ISYMC,'T') 1652c test 1653C 1654 CALL CCLR_E1C1(ETAC,WORK(KL1AM),WORK(KEI1),WORK(KEI2), 1655 * WORK(KEND3),LEND3,ISYML,ISYMC,'T') 1656C 1657 IF ( DEBUG.or.locdbg ) THEN 1658 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1659 WRITE(LUPRI,1) 'Norm of eta1 - after CCLR_E1C1: ',ETA1 1660 ENDIF 1661C 1662C--------------------------------------------------------------- 1663C etac2 = P(ab,ij)(sum(e)2L(aiej)*Ceb - sym(k)L(aibk)*C(jk)) 1664C--------------------------------------------------------------- 1665C 1666 IF (.NOT. CCS) THEN 1667C 1668 IF (CC2 .AND. IOPTCC2.EQ.1) THEN 1669 FCKCON = .TRUE. 1670 ETRAN = .FALSE. 1671 CALL CCRHS_EFCK(WORK(KEI1),WORK(KEI2),WORK(KCMO), 1672 * WORK(KFCKHF),WORK(KEND3),LEND3,FCKCON,ETRAN,ISYMC) 1673 END IF 1674 1675 CALL CC_EITR(WORK(KEI1),WORK(KEI2),WORK(KEND3),LEND3, 1676 * ISYMC) 1677C 1678 CALL CCRHS_E(ETAC(1+NT1AM(ISYRES)),WORK(KL2AM), 1679 * WORK(KEI1),WORK(KEI2),WORK(KEND3), 1680 * LEND3,ISYML,ISYMC) 1681C 1682 IF (IPRINT .GT. 40 ) THEN 1683 CALL AROUND( 'In CC_ETAC: EtaC vector ' ) 1684 CALL CC_PRP(ETAC(1),ETAC(1+NT1AM(ISYRES)),ISYMC,1,1) 1685 ENDIF 1686C 1687 IF (DEBUG .OR. ( IPRINT .GT. 20 )) THEN 1688 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1689 ETA2 = DDOT(NT2AM(ISYRES),ETAC(1+NT1AM(ISYRES)),1, 1690 * ETAC(1+NT1AM(ISYRES)),1) 1691 WRITE(LUPRI,1) 'Norm of eta1 - end of CC_ETAC: ',ETA1 1692 WRITE(LUPRI,1) 'Norm of eta2 - end of CC_ETAC: ',ETA2 1693 CALL AROUND( 'END OF CC_ETAC ') 1694 ENDIF 1695 ENDIF 1696C 1697 !if (cccisd.or.eomccsd.or.fdeom) then 1698 1699 !restart from zero to make sure not to mess up 1700 KETAC = KEND21 1701 KT2AM = KETAC + NT1AM(ISYRES) 1702 KT2SQ = KT2AM + NT2AM(1) 1703 KL2AM = KT2SQ + NT2SQ(1) 1704 KL2SQ = KL2AM + NT2AM(ISYML) 1705 KZint = KL2SQ + NT2SQ(ISYML) 1706 KEND3 = KZint + NT2SQ(isyml) 1707 LEND3 = LWORK - KEND3 1708 IF (LEND3 .LT. 0 ) 1709 & CALL QUIT('INSUFFICIENT WORK SPACE IN CC_EOMETA') 1710 1711C 1712C--------------------------------------------------------------- 1713C etac1 = 2sum(ei)L(ckei)*Cei 1714C requires packed L2, so I reread it into T2AM 1715C can also used it square, if last argument in called routine 1716C is set = 1 1717C--------------------------------------------------------------- 1718Css 1719 if (ioptl0.eq.1) then 1720 call dcopy(NT2AM(ISYML),XL0(1+NT1AM(ISYML)),1,WORK(KL2AM),1) 1721 else 1722 IOPT = 2 1723 CALL CC_RDRSP(LIST,ILSTNR,ISYML,IOPT,MODEL, 1724 * WORK(KEND3),WORK(KL2AM)) 1725 end if 1726 !CALL DSCAL(NT1AM(isymc),TWO,WORK(KINTAI),1) 1727 ! 1728! CALL CCG_LXD(ETAC(1),ISYRES,WORK(KINTAI),ISYMC, 1729! & WORK(KL2AM),ISYML,0) 1730 1731 !call dzero(WORK(KETAC),NT1AM(ISYRES)) 1732 !warning: result vector is zeroed inside! 1733 CALL CCG_LXD(WORK(KETAC),ISYRES,WORK(KINTAI),ISYMC, 1734 & WORK(KL2AM),ISYML,0) 1735 IF ( DEBUG.or.locdbg ) THEN 1736 ETA1 = DDOT(NT1AM(ISYRES),WORK(KETAC),1,WORK(KETAC),1) 1737 WRITE(LUPRI,1) 'Norm of alone Tbar_ck,ei*Cei: ',ETA1 1738 ENDIF 1739 !removed factor 2 to get FCI limit! 1740 CALL DAXPY(NT1AM(ISYRES),1.0D0,WORK(KETAC),1,ETAC(1),1) 1741 ! add to total tensor 1742 IF ( DEBUG.or.locdbg ) THEN 1743 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1744 WRITE(LUPRI,1)'Norm Eta^C after +Tbar_ck,ei*Cei:',ETA1 1745 ENDIF 1746C 1747C--------------------------------------------------------------- 1748C etac2 = 2(sum(nf)*Cnf(sum(bj)tbar(ck,bj)*ttilde(fn,bj)) 1749C the sum(bj)tbar(ck,bj)*ttilde(fn,bj) term is part of the 1750C d_aijb density!!! 1751C--------------------------------------------------------------- 1752C 1753 !I reread T2AM 1754 IOPT = 2 1755 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KT2AM)) 1756 ISYOPE = 1 1757 IOPTTCME = 1 1758 !get ttilde (overwrites T2am) 1759 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME) 1760 write(lupri,*) ' norm of ttilde ', 1761 & ddot(NT2AM(1),WORK(KT2AM),1,WORK(KT2AM),1) 1762 !square it up 1763 CALL CC_T2SQ(WORK(KT2AM),WORK(KT2SQ),1) 1764 !Square up L2 as well 1765 IF (.NOT. CCS) CALL CC_T2SQ(WORK(KL2AM),WORK(KL2SQ),ISYML) 1766 1767 !zero-ed inside the subroutine called 1768 call CC_ZIsq(work(KZint),work(kl2sq),ISYML, 1769 & work(kT2sq),1,WORK(kend3),LEND3) 1770 write(lupri,*) ' Norm of Z intermediate ', 1771 & ddot(NT2SQ(isyml),WORK(KZINT),1,WORK(KZINT),1) 1772 1773 !if (iprint.ge.45) then 1774 ! write(lupri,*) 'CC_EOMETA: Zint intermediate' 1775 ! call cc_prsq(work(kend3),WORK(KZint),1,0,1) 1776 !end if 1777 1778 call dzero(WORK(KETAC),NT1AM(ISYRES)) 1779 CALL CCG_LXD(WORK(KETAC),ISYRES,WORK(KINTIA),ISYMC, 1780 & WORK(KZINT),ISYML,1) 1781 1782 IF ( DEBUG.or.locdbg ) THEN 1783 ETA1 = DDOT(NT1AM(ISYRES),WORK(KETAC),1,WORK(KETAC),1) 1784 !ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1785 WRITE(LUPRI,1) 'Norm of eta1 (alone) Zint*C: ',ETA1 1786 ENDIF 1787 !removed factor 2 to get FCI limit 1788 CALL DAXPY(NT1AM(ISYRES),1.0D0,WORK(KETAC),1,ETAC(1),1) 1789 ! 1790 IF ( DEBUG.or.locdbg ) THEN 1791 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 1792 WRITE(LUPRI,1) 'Norm of eta after 2*Zint*Cei: ',ETA1 1793 ENDIF 1794 1795 !end if 1796 IF (IPRINT .GT. 5 ) THEN 1797 TIMEC = SECOND() - TIMEC 1798 WRITE(LUPRI,9999) 'CCCI_ETA^C ', TIMEC 1799 ENDIF 1800C 1801 1 FORMAT(1x,A35,1X,E20.10) 18029999 FORMAT(1x,'Time used in',2x,A18,2x,': ',f10.2,' seconds') 1803C 18048888 CONTINUE 1805 END 1806 1807C /* Deck cc_zisq */ 1808 SUBROUTINE CC_ZIsq(ZMAT,CLTR2,ISYMLV,T2AM,ISYMRV, 1809 & WORK,LWORK) 1810C 1811C Written by Sonia Coriani Halkier 27/04 - 2015 1812C 1813C Version: 1.0 1814C 1815C Purpose: To calculate the Z-intermediate entering the EOM 1816C etaC. 1817C 1818C It is assumed that both CLTR2 and T2am are squared up 1819C 1820C General symmetry of both CLTR2 and T2AM for F-mat and etaC. 1821C 1822#include "implicit.h" 1823 LOGICAL LOCDBG 1824 PARAMETER (LOCDBG = .false.) 1825 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1826 DIMENSION ZMAT(*),CLTR2(*),T2AM(*),WORK(LWORK) 1827#include "priunit.h" 1828#include "ccorb.h" 1829#include "ccsdsym.h" 1830#include "cclr.h" 1831C 1832 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1833C 1834 ISYRES = MULD2H(ISYMLV,ISYMRV) 1835C 1836C----------------------------- 1837C Initialize result array. 1838C----------------------------- 1839C 1840 CALL DZERO(ZMAT,NT2SQ(ISYRES)) 1841C 1842 DO 100 ISYMCK = 1,NSYM 1843C 1844 ISYMBJ = MULD2H(ISYMCK,ISYMLV) !left vector (tbar) 1845 ISYMFN = MULD2H(ISYMBJ,ISYMRV) !right vector (ttilde) 1846C 1847C------------------------------------- 1848C Contract the two matrices. 1849C------------------------------------------- 1850C 1851 KOFF1 = IT2SQ(ISYMCK,ISYMBJ) + 1 1852 KOFF2 = IT2SQ(ISYMBJ,ISYMFN) + 1 1853 KOFF3 = IT2SQ(ISYMCK,ISYMFN) + 1 1854C 1855 NTOTCK = MAX(NT1AM(ISYMCK),1) 1856 NTOTBJ = MAX(NT1AM(ISYMBJ),1) 1857 NTOTFN = MAX(NT1AM(ISYMFN),1) 1858C 1859 CALL DGEMM('N','N',NTOTCK,NTOTFN,NTOTBJ, 1860 & ONE,CLTR2(KOFF1),NTOTCK,T2AM(KOFF2),NTOTFN, 1861 & ONE,ZMAT(KOFF3),NTOTCK) 1862C 1863 120 CONTINUE 1864 110 CONTINUE 1865 100 CONTINUE 1866C 1867 IF (LOCDBG) THEN 1868 WRITE(LUPRI,*) '>>>>>CC_ZI: Norm of Z-intermediate:', 1869 & DDOT(NT2SQ(ISYRES),ZMAT,1,ZMAT,1) 1870 END IF 1871 1872! DO 100 ISYMCK = 1,NSYM 1873C 1874! ISYMBJ = MULD2H(ISYMCK,ISYMLV) !left vector (tbar) 1875! ISYMFN = MULD2H(ISYMBJ,ISYMRV) !right vector (ttilde) 1876C 1877C------------------------------------- 1878C Contract the two matrices. 1879C------------------------------------------- 1880C 1881! KOFF1 = IMATAB(ISYME,ISYMF) + 1 1882C 1883! NTOTCK = MAX(NT1AM(ISYMCK),1) 1884! NTOTBJ = MAX(NT1AM(ISYMBJ),1) 1885! NTOTFN = MAX(NT1AM(ISYMFN),1) 1886C 1887! CALL DGEMM('N','N',NTOTCK,NTOTFN,NT1AM(ISYMBJ), 1888! & ONE,CLTR2,NTOTCK,KT2AM,NTOTFN, 1889! & ONE,ZMAT,NTOTBJ) 1890C 1891! 120 CONTINUE 1892! 110 CONTINUE 1893! 100 CONTINUE 1894C 1895 RETURN 1896 END 1897