1C 2C Triplet versions of Xi and Eta routines 3C 4c*DECK CC_ETAC13 5 SUBROUTINE CC_ETAC13(ISYMC,LBLC, !Operator stuff 6 & ETAC, !Result vector 7 & LIST,ILSTNR, !Left vector 8 & LEOM, 9 & XINT,WORK,LWORK) 10C 11C----------------------------------------------------------------------- 12C Purpose: Calculate the etaC(L) vector, when L is a triplet vector 13C and C is a singlet operator. 14C This means that the result vector will be a triplet vector. 15C (Modified version of the CCSD etaC(l0/l1) code 16C 17C eta(tau_nu)= (<HF| + Sum(mu)L(1)<mu|) 18C exp(-t)[C,tau_nu]exp(T)|HF> 19C 20C Input: 21C 22C 23C LIST= 'L0' for zeroth order left amplitudes. 24C ISYML should be ISYMOP in this case. 25C 26C 'L1' for first order left amplitudes, read in from file 27C In this case the vector is found according to its list 28C number ILSTNR. 29C 30C For L1 HF contribution is skipped. 31C 32C C property integrals read according to LBLC 33C 34C SLV98,OC: Allow for input of integrals if 35C LBLC.eq.'GIVE INT' 36C 37C Sonia & Filip, Maj 2015 38C 39C----------------------------------------------------------------------- 40C 41 implicit none 42 43 character(len=*), parameter :: myname = 'CC_ETAC13' 44#include "priunit.h" 45#include "dummy.h" 46#include "maxorb.h" 47#include "ccorb.h" 48#include "iratdef.h" 49#include "cclr.h" 50#include "ccexci.h" 51#include "ccsdsym.h" 52#include "ccsdio.h" 53#include "ccsdinp.h" 54C 55!Sonia & Filip: find a better place 56#include "ccsections.h" 57!sonia 58#include "second.h" 59 60 DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0, HALF = 0.5D0, 61 & ZERO = 0.0D0, ONEM = -1.0D0 62 LOGICAL, PARAMETER :: LOCDBG = .false. 63C 64C Input: 65C 66 INTEGER, INTENT(IN) :: ISYMC, ! Symmetry of C 67 & ILSTNR, ! List number for input 68 & LWORK ! 69 CHARACTER(LEN=*), INTENT(IN) :: LBLC, ! 70 & LIST ! See above 71 DOUBLE PRECISION, INTENT(OUT):: ETAC(*), ! Output array 72 & WORK(LWORK) ! Work array 73 DOUBLE PRECISION, INTENT(INOUT) :: XINT(*) 74 LOGICAL, INTENT(IN) :: LEOM ! Whether to include extra EOM-CC 75 ! Contributions 76 77 CHARACTER MODEL*10 78 LOGICAL :: FCKCON, ETRAN 79 INTEGER :: KT1AM, KT2AM, KL1AM, KL2AM, KCTMO, KLAMDP, KLAMDH, 80 & KEI1, KEI2, KEND1, KEND2, KEND21, KEND3, KEND4, 81 & KINTAI, KINTIA, KXMAT, KYMAT, KETAC, KCINT 82 INTEGER :: KOFF1, KOFF2, KOFFP, KOFFM 83 INTEGER :: ISYML, ISYRES, ISYMA, ISYMI 84 INTEGER :: LEND1, LEND2, LEND21, LEND3, LEND4 85 INTEGER :: IOPT 86 DOUBLE PRECISION :: FF, XXI, XYI, ETA1, ETA2 87 DOUBLE PRECISION :: TIMEC 88C EXTERNAL FUNCTIONS 89 DOUBLE PRECISION :: DDOT 90 INTEGER :: ILSTSYM 91C 92 CALL AROUND( 'Constructing ETA^{'// LBLC 93 & //'}('// LIST //') vector ') 94C IF ( (IPRINT .GT. 10).or.locdbg ) THEN 95C CALL AROUND( 'IN CCCI_ETAC: Constructing ETAC(LE) vector ') 96C ENDIF 97C 98C-------------------------------- 99C find symmetry of D operator. 100C-------------------------------- 101C 102 ISYML = ILSTSYM(LIST,ILSTNR) 103C 104 ISYRES = MULD2H(ISYML,ISYMC) 105 IF (( LIST .EQ. 'L0')) THEN 106 CALL QUIT('Misuse of '//myname) 107 ENDIF 108C 109 TIMEC = SECOND() 110C 111 MODEL = 'CCSD ' 112 IF (CCS) MODEL = 'CCS ' 113 IF (CC2) MODEL = 'CC2 ' 114C 115C-------------------- 116C Allocate space. 117C-------------------- 118C 119 KCTMO = 1 120 KT1AM = KCTMO + N2BST(ISYMC) 121 KLAMDP = KT1AM + NT1AM(1) 122 KLAMDH = KLAMDP + NLAMDT 123 KEND1 = KLAMDH + NLAMDT 124C 125 LEND1 = LWORK - KEND1 126C 127 IF ( .NOT. CCS) THEN 128 129 KINTIA = KEND1 130 KEND1 = KINTIA + NT1AM(ISYMC) 131 LEND1 = LWORK - KEND1 132 CALL DZERO(WORK(KintIA),NT1AM(isymc)) 133C 134 KL1AM = KEND1 135 KL2AM = KL1AM + NT1AM(ISYML) 136 KEND2 = KL2AM + NT2SQ(ISYML) 137 LEND2 = LWORK - KEND2 138 KT2AM = KEND2 139 KEND21= KT2AM + MAX(NT2AM(1),NT2AM(ISYML)) 140 LEND21= LWORK - KEND21 141C 142 ELSE 143C 144 KL1AM = KEND1 145 KEND2 = KL1AM + NT1AM(ISYML) 146 LEND2 = LEND1 147 KEND21= KEND1 148 LEND21= LEND1 149C 150 ENDIF 151 KEI1 = KEND21 152 KEI2 = KEI1 + NEMAT1(ISYMC) 153 KEND3 = KEI2 + NMATIJ(ISYMC) 154 LEND3 = LWORK - KEND3 155 IF (LEND3.LT. 0 ) CALL QUIT(' TOO LITTLE WORKSPACE IN '//myname) 156C 157C----------------------- 158C get T1 amplitudes. 159C----------------------- 160C 161 CALL DZERO(WORK(KT1AM),NT1AM(1)) 162 IF ( .NOT. CCS) THEN 163 IOPT = 1 164 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 165 ENDIF 166C 167 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 168 * WORK(KEND21),LEND21) 169C 170C------------------------------- 171C get AO property integrals. 172C------------------------------- 173C 174 CALL DZERO(WORK(KCTMO),N2BST(ISYMC)) 175 FF = 1.0D0 176C SLV98,OC give integrals option 177 IF (LBLC.EQ.'GIVE INT') THEN 178 CALL DCOPY(N2BST(ISYMC),XINT(1),1,WORK(KCTMO),1) 179 ELSE 180 FF = 1.0D0 181 CALL CC_ONEP(WORK(KCTMO),WORK(KEND21),LEND21,FF,ISYMC,LBLC) 182 ENDIF 183C 184C----------------------------------------------- 185C Make MO T1-transformed property integrals. 186C----------------------------------------------- 187C 188 CALL CC_FCKMO(WORK(KCTMO),WORK(KLAMDP),WORK(KLAMDH), 189 * WORK(KEND21),LEND21,ISYMC,1,1) 190C 191C---------------------------------------------------------- 192C Extract Cia (stored ia) and reorder ai 193C---------------------------------------------------------- 194C 195 DO 100 ISYMI = 1,NSYM 196C 197 ISYMA = MULD2H(ISYMI,ISYMC) 198C 199 DO 110 A = 1,NVIR(ISYMA) 200C 201 DO 120 I = 1,NRHF(ISYMI) 202C 203 KOFF1 = KINTIA + IT1AM(ISYMA,ISYMI) 204 & + NVIR(ISYMA)*(I - 1) + A-1 205 KOFF2 = KCTMO + IFCVIR(ISYMI,ISYMA) 206 * + NORB(ISYMI)*(A - 1) + I - 1 207C 208 WORK(KOFF1) = WORK(KOFF2) 209C 210 120 CONTINUE 211 110 CONTINUE 212C 213 100 CONTINUE 214C 215C Initialize ETAC 216 CALL DZERO(ETAC,NT1AM(ISYRES)+2*NT2AM(ISYRES)) 217 218C---------------------------------------------- 219C Read L1, and L2(+) multipliers/left vectors. 220C---------------------------------------------- 221C 222 IOPT = 65 223 CALL CC_RDRSP(LIST,ILSTNR,ISYML,IOPT,MODEL, 224 * WORK(KL1AM),WORK(KT2AM)) 225C 226C-------------------------------- 227C Put C into E matrix format. 228C-------------------------------- 229C 230 FCKCON = .TRUE. 231 ETRAN = .FALSE. 232 CALL CCRHS_EFCK(WORK(KEI1),WORK(KEI2),WORK(KLAMDH), 233 * WORK(KCTMO),WORK(KEND3),LEND3,FCKCON, 234 * ETRAN,ISYMC) 235C 236C-------------------------------------------- 237C etac1 = sum(b)Lbi*Cba - sum(j)Laj*Cij. 238C-------------------------------------------- 239C 240 CALL CCLR_E1C1(ETAC,WORK(KL1AM),WORK(KEI1),WORK(KEI2), 241 * WORK(KEND3),LEND3,ISYML,ISYMC,'T') 242C ~ 243C Square L2(+) onto L2 244 IF (.NOT. CCS) THEN 245 CALL CCRHS3_R2IJ(WORK(KT2AM),WORK(KEND3),LEND3,ISYML) 246 CALL CC_T2SQ(WORK(KT2AM),WORK(KL2AM),ISYML) 247C--------------------------------- 248C Calculate the doubles (+)G term 249C--------------------------------- 250C Transpose i,j and a,b blocks of C 251 CALL CC_EITR(WORK(KEI1),WORK(KEI2),WORK(KEND3),LEND3,ISYMC) 252C Scale C by 1/2, as this term has a factor of 1/2 253 CALL DSCAL(NMATAB(ISYMC),HALF,WORK(KEI1),1) 254 CALL DSCAL(NMATIJ(ISYMC),HALF,WORK(KEI2),1) 255C Calculate actual term 256 KOFFP = NT1AM(ISYRES) + 1 257 CALL CCRHS_E(ETAC(KOFFP),WORK(KL2AM),WORK(KEI1),WORK(KEI2), 258 & WORK(KEND3),LEND3,ISYML,ISYMC) 259C Restore factor of C 260 CALL DSCAL(NMATAB(ISYMC),TWO,WORK(KEI1),1) 261 CALL DSCAL(NMATIJ(ISYMC),TWO,WORK(KEI2),1) 262 END IF 263C 264C---------------------------------------------- 265C Read L2(-) multipliers/left vector. 266C---------------------------------------------- 267C The routine needs to be modified to allow this 268 IOPT = 128 269 CALL CC_RDRSP(LIST,ILSTNR,ISYML,IOPT,MODEL, 270 * WORK(KL1AM),WORK(KT2AM)) 271C 272C We need L(+) - L(-) for X and Y in the following. 273 CALL CC_T2SQ3A(WORK(KT2AM),WORK(KL2AM),ISYML,ONEM) 274C 275C------------------------------------------ 276C Put T2 (packed) amplitudes in etac2. 277C------------------------------------------ 278C 279 IF (.NOT. CCS) THEN 280 !read in T2AM (packed) 281 IOPT = 2 282 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KT2AM)) 283 ENDIF 284C 285C-------------------------------- 286C Make X and Y intermediates. 287C-------------------------------- 288CRF Triplet requires (L(+) - L(-)) 289 IF (.NOT. CCS) THEN 290 KXMAT = KEND3 291 KYMAT = KXMAT + NMATIJ(ISYML) 292 KEND4 = KYMAT + NMATAB(ISYML) 293 LEND4 = LWORK - KEND4 294 IF (LEND4.LT. 0 ) 295 & CALL QUIT(' TOO LITTLE WORKSPACE IN CC_ETAC-2') 296C 297 IF ( DEBUG.or.LOCDBG ) THEN 298 XYI = DDOT(NT1AM(ISYML),WORK(KL1AM),1,WORK(KL1AM),1) 299 WRITE(LUPRI,1) 'CC_ETAC: L1AM vector: ',XYI 300 XYI = DDOT(NT2SQ(ISYML),WORK(KL2AM),1,WORK(KL2AM),1) 301 WRITE(LUPRI,1) 'CC_ETAC: L2AM vector: ',XYI 302 XXI = DDOT(NT2AM(ISYMOP),WORK(KT2AM),1,WORK(KT2AM),1) 303 WRITE(LUPRI,1) 'T2AM vector : ',XXI 304 ENDIF 305 CALL CC_XI(WORK(KXMAT),WORK(KL2AM),ISYML,WORK(KT2AM),1, 306 * WORK(KEND4),LEND4) 307 CALL CC_YI(WORK(KYMAT),WORK(KL2AM),ISYML,WORK(KT2AM),1, 308 * WORK(KEND4),LEND4) 309 IF ( DEBUG.or.LOCDBG ) THEN 310 XYI = DDOT(NMATAB(ISYML),WORK(KYMAT),1,WORK(KYMAT),1) 311 WRITE(LUPRI,1) 'CC_ETAC: YI intermediate is: ',XYI 312 XXI = DDOT(NMATIJ(ISYML),WORK(KXMAT),1,WORK(KXMAT),1) 313 WRITE(LUPRI,1) 'CC_ETAC: XI intermediate is: ',XXI 314 ENDIF 315 ELSE 316 KEND4 = KEND3 317 LEND4 = LEND3 318 ENDIF 319C 320C---------------------------------------------- 321C Calculate X and Y contributions to etac1. 322C etac1 = -sum(e)Cie*Yae - sum(l)Cla*Xli 323C---------------------------------------------- 324C 325 IF ( .NOT.CCS ) THEN 326C 327 CALL CC_21EFM(ETAC,WORK(KCTMO),ISYMC,WORK(KXMAT), 328 * WORK(KYMAT),ISYML,WORK(KEND4),LEND4) 329 IF ( DEBUG.or.locdbg ) THEN 330 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 331 WRITE(LUPRI,1) 'Norm of eta1-after X&Y cont: ',ETA1 332 ENDIF 333 ENDIF 334C We have L(+)-L(-), we need L(+)+L(-) 335 CALL CC_T2SQTRANSP(WORK(KL2AM),ISYML) 336C 337 IF (LEOM) THEN 338C 339C--------------------------------------------------------------- 340C EOM contribution to ETAC: ~ 341C etac = 2sum(ei)(L(+) + L(-))_{ckei}*(C_{ei} + t_{ei,fn}*C_{nf}) 342C--------------------------------------------------------------- 343 KCINT = KEND3 344 KETAC = KCINT 345 KINTAI = KCINT + MAX(NT1AM(ISYMC),NT1AM(ISYRES)) 346C Extract C_{ai} 347 DO ISYMI = 1, NSYM 348 ISYMA = MULD2H(ISYMI,ISYMC) 349 DO I = 1, NRHF(ISYMI) 350 KOFF1 = KINTAI + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) 351 KOFF2 = KCTMO + IFCRHF(ISYMA,ISYMI) + 352 & NORB(ISYMA)*(I-1) + NRHF(ISYMA) 353 CALL DCOPY(NVIR(ISYMA),WORK(KOFF2),1,WORK(KOFF1),1) 354 END DO 355 END DO 356 !get ttilde (overwrites T2am) 357 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,1,1) 358C ~ 359C Add t_{em,fn} *C_{nf} to C_{em} 360 CALL CCG_LXD(WORK(KCINT),ISYMC,WORK(KINTIA),ISYMC, 361 & WORK(KT2AM),1,0) 362 CALL DAXPY(NT1AM(ISYMC),1.D0,WORK(KCINT),1,WORK(KINTAI),1) 363C Calculate the term as L_{ai,em} * C'_{em} 364 CALL CCG_LXD(WORK(KETAC),ISYRES,WORK(KINTAI),ISYMC, 365 & WORK(KL2AM),ISYML,1) 366 IF ( DEBUG.or.locdbg ) THEN 367 ETA1 = DDOT(NT1AM(ISYRES),WORK(KETAC),1,WORK(KETAC),1) 368 WRITE(LUPRI,1) 'Norm of alone Tbar_ck,ei*Cei: ',ETA1 369 ENDIF 370 !removed factor 2 to get FCI limit! 371 CALL DAXPY(NT1AM(ISYRES),1.0D0,WORK(KETAC),1,ETAC(1),1) 372 END IF 373C 374C--------------------------------- 375C Calculate the doubles (-)G term 376C--------------------------------- 377C 378 KOFFM = KOFFP + NT2AM(ISYRES) 379 CALL CC_T2SQSYMSCAL(WORK(KL2AM),ISYML,ZERO) 380C 381 CALL CCRHS_E3(DUMMY,.FALSE.,WORK(KL2AM),WORK(KEI1),WORK(KEI2), 382 & WORK(KEND3),LEND3,ISYML,ISYMC, 383 & ETAC(KOFFM),.TRUE.) 384C 385C------------------------------------------------ 386C Workspace for T2AM and X and Y is now free. 387C etac2 = P(ab,ij)(2l(ai)*Cjb - l(aj)*c(ib)) 388C------------------------------------------------ 389C 390 IF (.NOT. CCS) THEN 391C 392 CALL CC_L1FCK3P(ETAC(1+NT1AM(ISYRES)),WORK(KL1AM),WORK(KINTIA), 393 * ISYML,ISYMC) 394C 395 IF ( DEBUG.or.locdbg ) THEN 396 ETA1 = DDOT(NT1AM(ISYRES),ETAC(1),1,ETAC(1),1) 397 ETA2 = DDOT(2*NT2AM(ISYRES),ETAC(1+NT1AM(ISYRES)),1, 398 * ETAC(1+NT1AM(ISYRES)),1) 399 WRITE(LUPRI,1) 'Norm of eta1-after L1c cont: ',ETA1 400 WRITE(LUPRI,1) 'Norm of eta2-after L1c cont: ',ETA2 401 ENDIF 402 ENDIF 403C 404C Permute indices of plus vector 405 CALL CCRHS3_IJ(ETAC(KOFFP),WORK(KEND2),LEND2,ISYRES) 406C For now, explicitly zero (-) diagonal! 407 IF (ISYRES.EQ.1) CALL CCLR_DIASCL(ETAC(KOFFM),ZERO,ISYRES) 408C 409 IF (IPRINT .GT. 5 ) THEN 410 TIMEC = SECOND() - TIMEC 411 WRITE(LUPRI,9999) 'CCCI_ETA^C ', TIMEC 412 ENDIF 413C 414 1 FORMAT(1x,A35,1X,E20.10) 4159999 FORMAT(1x,'Time used in',2x,A18,2x,': ',f10.2,' seconds') 416C 417 END 418C 419 SUBROUTINE CC_XIC13(ISYMC,LBLC, !Operator stuff 420 & XIC, !Result vector 421 & LIST,ILSTNR, !Left vector 422 & LEOM, 423 & XINT,WORK,LWORK) 424C 425C----------------------------------------------------------------------- 426C Purpose: Calculate the xiC(R) vector, when R is a triplet vector 427C and C is a singlet operator. 428C This means that the result vector will be a triplet vector. 429C (Modified version of the CCSD etaC(l0/l1) code 430C 431C Input: 432C 433C ISYMC Symmetry of C 434C LBLC Label of C integrals 435C XIC Result vector (output) 436C LIST List type of R vector 437C ILSTNR List number of R vector 438C LEOM Whether to include EOM disconnected terms (logical) 439C XINT Can be used to pass C integrals in memory 440C WORK 441C LWORK 442C 443C R. Faber 2017 444C 445C----------------------------------------------------------------------- 446C 447 implicit none 448 449 character(len=*), parameter :: myname = 'CC_XI13' 450#include "priunit.h" 451#include "dummy.h" 452#include "maxorb.h" 453#include "ccorb.h" 454#include "iratdef.h" 455#include "cclr.h" 456#include "ccexci.h" 457#include "ccsdsym.h" 458#include "ccsdio.h" 459#include "ccsdinp.h" 460C 461!Sonia & Filip: find a better place 462#include "ccsections.h" 463#include "second.h" 464 465 DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0, HALF = 0.5D0, 466 & ZERO = 0.0D0, ONEM = -1.0D0 467 LOGICAL, PARAMETER :: LOCDBG = .false. 468C 469C Input: 470C 471 INTEGER, INTENT(IN) :: ISYMC, ! Symmetry of C 472 & ILSTNR, ! List number for input 473 & LWORK ! 474 CHARACTER(LEN=*), INTENT(IN) :: LBLC, ! 475 & LIST ! See above 476 DOUBLE PRECISION, INTENT(OUT):: XIC(*), ! Output array 477 & WORK(LWORK) ! Work array 478 DOUBLE PRECISION, INTENT(INOUT) :: XINT(*) 479 LOGICAL, INTENT(IN) :: LEOM ! Whether to include extra EOM-CC 480 ! Contributions 481 482 CHARACTER MODEL*10 483 LOGICAL :: FCKCON, ETRAN 484 INTEGER :: KT1AM, KT2AM, KL1AM, KL2AM, KCTMO, KLAMDP, KLAMDH, 485 & KEI1, KEI2, KEND1, KEND2, KEND21, KEND3, 486 & KINTAI, KINTIA, KXMAT, KYMAT, KETAC, KCINT 487 INTEGER :: KOFF1, KOFF2, KOFFP, KOFFM 488 INTEGER :: ISYMR, ISYRES, ISYMA, ISYMI 489 INTEGER :: LEND1, LEND2, LEND21, LEND3 490 INTEGER :: IOPT 491 DOUBLE PRECISION :: FF, XXI, XYI, ETA1, ETA2 492 DOUBLE PRECISION :: TIMEC 493C EXTERNAL FUNCTIONS 494 DOUBLE PRECISION :: DDOT 495 INTEGER :: ILSTSYM 496C 497 CALL AROUND( 'Constructing XI^{'// LBLC 498 & //'}('// LIST //') vector ') 499C IF ( (IPRINT .GT. 10).or.locdbg ) THEN 500C CALL AROUND( 'IN CCCI_ETAC: Constructing ETAC(LE) vector ') 501C ENDIF 502C 503C-------------------------------- 504C find symmetry of D operator. 505C-------------------------------- 506C 507 ISYMR = ILSTSYM(LIST,ILSTNR) 508C 509 ISYRES = MULD2H(ISYMR,ISYMC) 510 IF (( LIST .EQ. 'R0')) THEN 511 CALL QUIT('Misuse of '//myname) 512 ENDIF 513C 514 TIMEC = SECOND() 515C 516 MODEL = 'CCSD ' 517 IF (CCS) MODEL = 'CCS ' 518 IF (CC2) MODEL = 'CC2 ' 519C 520C-------------------- 521C Allocate space. 522C-------------------- 523C 524 KCTMO = 1 525 KT1AM = KCTMO + N2BST(ISYMC) 526 KLAMDP = KT1AM + NT1AM(ISYMOP) 527 KLAMDH = KLAMDP + NLAMDT 528 KEND1 = KLAMDH + NLAMDT 529C 530 LEND1 = LWORK - KEND1 531 IF ( .NOT. CCS) THEN 532 533 KINTIA = KEND1 534 KEND1 = KINTIA + NT1AM(ISYMC) 535 LEND1 = LWORK - KEND1 536 CALL DZERO(WORK(KintIA),NT1AM(isymc)) 537C 538 KL1AM = KEND1 539 KL2AM = KL1AM + NT1AM(ISYMR) 540 KEND2 = KL2AM + NT2SQ(ISYMR) 541 LEND2 = LWORK - KEND2 542 KT2AM = KEND2 543 KEND21= KT2AM + MAX(NT2AM(1),NT2AM(ISYMR)) 544 LEND21= LWORK - KEND21 545C 546 ELSE 547C 548 KL1AM = KEND1 549 KEND2 = KL1AM + NT1AM(ISYMR) 550 LEND2 = LEND1 551 KEND21= KEND1 552 LEND21= LEND1 553C 554 ENDIF 555 KEI1 = KEND21 556 KEI2 = KEI1 + NEMAT1(ISYMC) 557 KEND3 = KEI2 + NMATIJ(ISYMC) 558 LEND3 = LWORK - KEND3 559 IF (LEND3.LT. 0 ) CALL QUIT(' TOO LITTLE WORKSPACE IN '//myname) 560 IF (LEND21.LT. 0 ) CALL QUIT(' TOO LITTLE WORKSPACE IN '//myname) 561C 562 IOPT = 3 563C CALL CC_RDRSP(LIST,ILSTNR,ISYMR,IOPT,MODEL, 564C * ETAC,ETAC(1+NT1AM(ISYMR))) 565C ff = ddot(NT1AM(isyml)+2*nt2am(isyml),etac,1,etac,1) 566C write(lupri,*) 'Norm of total', ff 567 568C 569C----------------------- 570C get T1 amplitudes. 571C----------------------- 572C 573 CALL DZERO(WORK(KT1AM),NT1AM(1)) 574 IF ( .NOT. CCS) THEN 575 IOPT = 1 576 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 577 ENDIF 578C 579 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 580 * WORK(KEND21),LEND21) 581C 582C------------------------------- 583C get AO property integrals. 584C------------------------------- 585C 586 CALL DZERO(WORK(KCTMO),N2BST(ISYMC)) 587 FF = 1.0D0 588C SLV98,OC give integrals option 589 IF (LBLC.EQ.'GIVE INT') THEN 590 CALL DCOPY(N2BST(ISYMC),XINT(1),1,WORK(KCTMO),1) 591 ELSE 592 FF = 1.0D0 593 CALL CC_ONEP(WORK(KCTMO),WORK(KEND21),LEND21,FF,ISYMC,LBLC) 594 ENDIF 595C 596C----------------------------------------------- 597C Make MO T1-transformed property integrals. 598C----------------------------------------------- 599C 600 CALL CC_FCKMO(WORK(KCTMO),WORK(KLAMDP),WORK(KLAMDH), 601 * WORK(KEND21),LEND21,ISYMC,1,1) 602C 603C---------------------------------------------------------- 604C Extract Cia (stored ia) and reorder ai 605C---------------------------------------------------------- 606C 607 DO 100 ISYMI = 1,NSYM 608C 609 ISYMA = MULD2H(ISYMI,ISYMC) 610C 611 DO 110 A = 1,NVIR(ISYMA) 612C 613 DO 120 I = 1,NRHF(ISYMI) 614C 615 KOFF1 = KINTIA + IT1AM(ISYMA,ISYMI) 616 & + NVIR(ISYMA)*(I - 1) + A-1 617 KOFF2 = KCTMO + IFCVIR(ISYMI,ISYMA) 618 * + NORB(ISYMI)*(A - 1) + I - 1 619C 620 WORK(KOFF1) = WORK(KOFF2) 621C 622 120 CONTINUE 623 110 CONTINUE 624C 625 100 CONTINUE 626C 627C Initialize XIC 628 CALL DZERO(XIC,NT1AM(ISYRES)+2*NT2AM(ISYRES)) 629 630C---------------------------------------------- 631C Read R1, and R2(+) multipliers/left vectors. 632C---------------------------------------------- 633C 634 IOPT = 65 635 CALL CC_RDRSP(LIST,ILSTNR,ISYMR,IOPT,MODEL, 636 * WORK(KL1AM),WORK(KT2AM)) 637C 638C-------------------------------- 639C Put C into E matrix format. 640C-------------------------------- 641C 642 FCKCON = .TRUE. 643 ETRAN = .FALSE. 644 CALL CCRHS_EFCK(WORK(KEI1),WORK(KEI2),WORK(KLAMDH), 645 * WORK(KCTMO),WORK(KEND3),LEND3,FCKCON, 646 * ETRAN,ISYMC) 647C ~ 648C Square L2(+) onto L2 649 IF (.NOT. CCS) THEN 650 CALL CCRHS3_R2IJ(WORK(KT2AM),WORK(KEND3),LEND3,ISYMR) 651 CALL CC_T2SQ(WORK(KT2AM),WORK(KL2AM),ISYMR) 652C--------------------------------- 653C Calculate the doubles (+)G term 654C--------------------------------- 655C Scale C by 1/2, as this term has a factor of 1/2 656 CALL DSCAL(NMATAB(ISYMC),HALF,WORK(KEI1),1) 657 CALL DSCAL(NMATIJ(ISYMC),HALF,WORK(KEI2),1) 658C Calculate actual term 659 KOFFP = NT1AM(ISYRES) + 1 660 CALL CCRHS_E(XIC(KOFFP),WORK(KL2AM),WORK(KEI1),WORK(KEI2), 661 & WORK(KEND3),LEND3,ISYMR,ISYMC) 662C Restore factor of C 663 CALL DSCAL(NMATAB(ISYMC),TWO,WORK(KEI1),1) 664 CALL DSCAL(NMATIJ(ISYMC),TWO,WORK(KEI2),1) 665C 666C---------------------------------------------- 667C Read R2(-) multipliers/left vector. 668C---------------------------------------------- 669C The routine needs to be modified to allow this 670 IOPT = 128 671 CALL CC_RDRSP(LIST,ILSTNR,ISYMR,IOPT,MODEL, 672 * WORK(KL1AM),WORK(KT2AM)) 673C 674C We need R(+) + R(-) in the following. 675 CALL CC_T2SQ3A(WORK(KT2AM),WORK(KL2AM),ISYMR,1.0D0) 676C 677C Doubles contribution to Xi1: 678C Xi1 = sum_{em} C_{me} ((+)R+(-)R)_{ai,em} 679 CALL CCG_LXD(XIC,ISYRES,WORK(KINTIA),ISYMC, 680 & WORK(KL2AM),ISYMR,1) 681C 682C--------------------------------- 683C Calculate the doubles (-)G term 684C--------------------------------- 685C 686 KOFFM = KOFFP + NT2AM(ISYRES) 687 CALL CC_T2SQSYMSCAL(WORK(KL2AM),ISYMR,ZERO) 688 CALL CCRHS_E3(DUMMY,.FALSE.,WORK(KL2AM),WORK(KEI1),WORK(KEI2), 689 & WORK(KEND3),LEND3,ISYMR,ISYMC, 690 & XIC(KOFFM),.TRUE.) 691C 692 END IF 693C 694C-------------------------------------------- 695C xi1 = sum(b)Lbi*Cba - sum(j)Laj*Cij. 696C-------------------------------------------- 697C 698 CALL CCLR_E1C1(XIC,WORK(KL1AM),WORK(KEI1),WORK(KEI2), 699 * WORK(KEND3),LEND3,ISYMR,ISYMC,'N') 700C 701C----------------------------------------- 702C Put T2 (packed) amplitudes in KT2AM. 703C----------------------------------------- 704C 705 IF (.NOT. CCS) THEN 706 !read in T2AM (packed) 707 IOPT = 2 708 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KT2AM)) 709 ENDIF 710C 711 IF (LEOM) THEN 712C 713C--------------------------------------------------------------- 714C EOM contribution to XIC: ~ 715C (+/-) XiC2 = R_{ai}*(C_{bj} + t_{bj,fn}*C_{nf}) 716C--------------------------------------------------------------- 717 KCINT = KEND21 718 KINTAI = KCINT + NT1AM(ISYMC) 719C Extract C_{ai} 720 DO ISYMI = 1, NSYM 721 ISYMA = MULD2H(ISYMI,ISYMC) 722 DO I = 1, NRHF(ISYMI) 723 KOFF1 = KINTAI + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) 724 KOFF2 = KCTMO + IFCRHF(ISYMA,ISYMI) + 725 & NORB(ISYMA)*(I-1) + NRHF(ISYMA) 726 CALL DCOPY(NVIR(ISYMA),WORK(KOFF2),1,WORK(KOFF1),1) 727 END DO 728 END DO 729 !get ttilde (overwrites T2am) 730 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,1,1) 731C ~ 732C Add t_{em,fn} *C_{nf} to C_{em} 733 CALL CCG_LXD(WORK(KCINT),ISYMC,WORK(KINTIA),ISYMC, 734 & WORK(KT2AM),1,0) 735 CALL DAXPY(NT1AM(ISYMC),1.D0,WORK(KCINT),1,WORK(KINTAI),1) 736C Calculate the term as R_{ai} * C'_{bj} 737 CALL CC_L1FCK3P(XIC(1+NT1AM(ISYRES)),WORK(KL1AM),WORK(KINTAI), 738 * ISYMR,ISYMC) 739 IF ( DEBUG.or.locdbg ) THEN 740 ETA1 = DDOT(NT1AM(ISYRES),WORK(KETAC),1,WORK(KETAC),1) 741 WRITE(LUPRI,1) 'Norm of alone Tbar_ck,ei*Cei: ',ETA1 742 ENDIF 743 END IF 744C 745C------------------------------------------------ 746C Workspace for T2AM and X and Y is now free. 747C etac2 = P(ab,ij)(2l(ai)*Cjb - l(aj)*c(ib)) 748C------------------------------------------------ 749C 750 IF (.NOT. CCS) THEN 751C 752C 753 IF ( DEBUG.or.locdbg ) THEN 754 ETA1 = DDOT(NT1AM(ISYRES),XIC(1),1,XIC(1),1) 755 ETA2 = DDOT(2*NT2AM(ISYRES),XIC(1+NT1AM(ISYRES)),1, 756 * XIC(1+NT1AM(ISYRES)),1) 757 WRITE(LUPRI,1) 'Norm of eta1-after L1c cont: ',ETA1 758 WRITE(LUPRI,1) 'Norm of eta2-after L1c cont: ',ETA2 759 ENDIF 760 ENDIF 761C 762C Permute indices of plus vector 763 CALL CCRHS3_IJ(XIC(KOFFP),WORK(KEND2),LEND2,ISYRES) 764C For now, explicitly zero (-) diagonal! 765 CALL CCLR_DIASCL(XIC(KOFFM),ZERO,ISYRES) 766 767 768 IF (IPRINT .GT. 5 ) THEN 769 TIMEC = SECOND() - TIMEC 770 WRITE(LUPRI,9999) 'CCCI_ETA^C ', TIMEC 771 ENDIF 772C 773 1 FORMAT(1x,A35,1X,E20.10) 7749999 FORMAT(1x,'Time used in',2x,A18,2x,': ',f10.2,' seconds') 775C 776 END 777 778