1C /* Deck cc_t2sq3a */ 2 SUBROUTINE CC_T2SQ3A(T2AM,T2SQ,ISYM,FACT) 3! 4!-------------------------------------------------------- 5! Rasmus Faber 2017: 6! 7! Adds "FACT" times the content of the antisymmetric, 8! packed array T2AM to the square array T2SQ 9! 10! Based on CC_T2SQ3 by Kasper Hald 11!-------------------------------------------------------- 12! 13 IMPLICIT NONE 14!if defined (SYS_CRAY) 15! REAL T2AM(*), T2SQ(*) 16!else 17 DOUBLE PRECISION, PARAMETER :: ONEM = -1.0D0 18 DOUBLE PRECISION, INTENT(IN) :: T2AM(*), FACT 19 DOUBLE PRECISION, INTENT(INOUT) :: T2SQ(*) 20!endif 21 INTEGER, INTENT(IN) :: ISYM 22 INTEGER KOFF1, KOFF2, ISYMBJ, KOFF, ISYMAI, NAMP, NAI 23 INTEGER NBJ 24#include "priunit.h" 25#include "ccorb.h" 26#include "ccsdsym.h" 27! 28 CALL QENTER('CC_T2SQ3') 29! 30 IF (ISYM.EQ.1) THEN 31 KOFF1 = 1 32 KOFF2 = 1 33 DO 100 ISYMBJ = 1,NSYM 34 IF (NT1AM(ISYMBJ) .GT. 0) THEN 35 CALL SQMATR3A(NT1AM(ISYMBJ),T2AM(KOFF1),T2SQ(KOFF2),FACT) 36 KOFF1 = KOFF1 + NT1AM(ISYMBJ)*(NT1AM(ISYMBJ)+1)/2 37 KOFF2 = KOFF2 + NT1AM(ISYMBJ)*NT1AM(ISYMBJ) 38 ENDIF 39 100 CONTINUE 40! 41 ELSE 42! 43 KOFF = 1 44 DO 200 ISYMBJ = 1,NSYM 45 ISYMAI = MULD2H(ISYM,ISYMBJ) 46! 47 IF (ISYMBJ.LT.ISYMAI) CYCLE 48! 49 NAMP = NT1AM(ISYMAI)*NT1AM(ISYMBJ) 50 KOFF = IT2AM(ISYMAI,ISYMBJ) + 1 51! 52 IF (NAMP .GT. 0) THEN 53 KOFF1 = IT2SQ(ISYMAI,ISYMBJ) + 1 54 CALL DAXPY(NAMP,FACT,T2AM(KOFF),1,T2SQ(KOFF1),1) 55 NAI = MAX(NT1AM(ISYMAI),1) 56 NBJ = MAX(NT1AM(ISYMBJ),1) 57 KOFF2 = IT2SQ(ISYMBJ,ISYMAI) + 1 58 CALL TRM3A(T2AM(KOFF),NT1AM(ISYMAI),NT1AM(ISYMBJ), 59 * T2SQ(KOFF2),ONEM*FACT) 60! 61 ENDIF 62! 63 200 CONTINUE 64! 65 ENDIF 66! 67 CALL QEXIT('CC_T2SQ3') 68! 69 RETURN 70 CONTAINS 71! 72 SUBROUTINE TRM3A(A,M,N,B,FACTM) 73! 74!--------------------------------------------------------------- 75! 76! Adds FACT times the transpose of A to the array B 77! 78! Based on trm3 by Kasper Hald 79! Based on TRM by Ove Christiansen. 80!--------------------------------------------------------------- 81! 82 INTEGER, INTENT(IN) :: M, N 83 DOUBLE PRECISION, INTENT(IN) :: A(M,N), FACTM 84 DOUBLE PRECISION, INTENT(INOUT) :: B(N,M) 85! 86 DOUBLE PRECISION :: XMONE 87 INTEGER :: I 88! 89 DO 100 I = 1, M 90! 91C CALL DAXPY(N,FACT,A(I,1),M,B(1,I),1) 92 DO J = 1, N 93 B(J,I) = B(J,I) + FACTM*A(I,J) 94 END DO 95! 96 100 CONTINUE 97! 98 RETURN 99 END SUBROUTINE 100C 101 PURE SUBROUTINE SQMATR3A(NDIM,PKMAT,SQMAT,FACT) 102! 103!----------------------------------------------------- 104! 105! This subroutine adds the packed antisymmetric matrix 106! PKMAT*FACT to the square matrix SQMAT 107! 108! Based on SQMATR3 by Kasper Hald. 109! Based on SQMATR by Henrik Koch. 110!----------------------------------------------------- 111! 112 INTEGER, INTENT(IN) :: NDIM 113 DOUBLE PRECISION, INTENT(IN) :: FACT, PKMAT(*) 114 DOUBLE PRECISION, INTENT(INOUT) :: SQMAT(NDIM,NDIM) 115! 116 INTEGER I, J, IJ 117 DOUBLE PRECISION ZERO,XMONE 118! 119 PARAMETER(XMONE = -1.0D00) 120! 121 IJ = 0 122 DO 100 I = 1,NDIM 123 DO 110 J = 1,I-1 124 IJ = IJ + 1 125 SQMAT(I,J) = -FACT*PKMAT(IJ) + SQMAT(I,J) 126 SQMAT(J,I) = FACT*PKMAT(IJ) + SQMAT(J,I) 127 110 CONTINUE 128 IJ = IJ + 1 129 SQMAT(I,I) = 0.0D0 130 100 CONTINUE 131! 132 RETURN 133 END SUBROUTINE 134 135 END 136C /* Deck cc_12c3 */ 137 SUBROUTINE CC_12C3(RHO2,CTR1,ISYMV,XLAMDH,ISYMH, 138 * ISINT,WORK,LWORK, 139 * LUO3,O3FIL) 140C 141C Rasmus Faber - 2017 142C 143C Based on CC_21B12C by Asger Halkier & Henrik Koch 13/9 - 1995. 144C 145C Note: the file-read performed here is also done in 146C CC_21G: we might want to merge those routines. 147C 148 IMPLICIT NONE 149 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, ONE = 1.0D0, 150 & ONEM = -1.0D0 151 CHARACTER(len=*), PARAMETER :: myname = 'CC_12C3' 152 DOUBLE PRECISION, INTENT(INOUT) :: RHO2(*) 153 DOUBLE PRECISION, INTENT(IN) :: CTR1(*),XLAMDH(*) 154 DOUBLE PRECISION, INTENT(OUT) :: WORK(LWORK) 155 CHARACTER, INTENT(IN) :: O3FIL*(*) 156 INTEGER, INTENT(IN) :: ISYMV, ISYMH, ISINT, LWORK, LUO3 157#include "priunit.h" 158#include "ccorb.h" 159#include "ccsdsym.h" 160#include "ccsdinp.h" 161#include "cclr.h" 162C 163 INTEGER :: ISYMA, ISYMB, ISYMD, ISYMI, ISYMJ, ISYMK, 164 & ISYMAI, ISYMBJ, ISYMIK, ISYIKJ, 165 & ISYMLI, ISYRES 166 INTEGER :: NTOT, NTOTA, NTOTD, NTOTI, NTOIKJ 167 INTEGER :: KSCR1, KSCR3, KOFF1, KOFF2, KOFF3, KEND1 168 INTEGER :: LWRK1, IOFF, NBJ 169C 170 CALL QENTER(myname) 171C 172 ISYMLI = MULD2H(ISINT,ISYMH) 173 ISYRES = MULD2H(ISYMV,ISYMLI) 174C 175 DO 100 ISYMD = 1,NSYM 176C 177 IF (NBAS(ISYMD) .EQ. 0) GOTO 100 178C 179 ISYIKJ = MULD2H(ISINT,ISYMD) 180 ISYMB = MULD2H(ISYMH,ISYMD) 181C 182C--------------------------------- 183C Allocation of work space. 184C--------------------------------- 185C 186 KSCR1 = 1 187 KSCR3 = KSCR1 + NMAIJK(ISYIKJ)*NVIR(ISYMB) 188 KEND1 = KSCR3 + NMAIJK(ISYIKJ)*NBAS(ISYMD) 189 LWRK1 = LWORK - KEND1 190C 191 IF (LWRK1 .LT. 0) THEN 192 CALL QUIT('Insufficient work space area in '//myname) 193 ENDIF 194C 195C------------------------------------ 196C Read integrals from disc. 197C------------------------------------ 198C 199 NTOT = NMAIJK(ISYIKJ)*NBAS(ISYMD) 200 IOFF = I3ODEL(ISYIKJ,ISYMD) + 1 201C 202 CALL GETWA2(LUO3,O3FIL,WORK(KSCR3),IOFF,NTOT) 203C 204C--------------------------- 205C Transform AO index. 206C--------------------------- 207C 208 KOFF1 = ILMVIR(ISYMB) + 1 209C 210 NTOIKJ = MAX(NMAIJK(ISYIKJ),1) 211 NTOTD = MAX(NBAS(ISYMD),1) 212C 213 CALL DGEMM('N','N',NMAIJK(ISYIKJ),NVIR(ISYMB),NBAS(ISYMD), 214 * ONE,WORK(KSCR3),NTOIKJ,XLAMDH(KOFF1),NTOTD, 215 * ZERO,WORK(KSCR1),NTOIKJ) 216C 217 DO 110 B = 1,NVIR(ISYMB) 218C 219 DO 170 ISYMJ = 1,NSYM 220C 221 ISYMBJ = MULD2H(ISYMJ,ISYMB) 222 ISYMAI = MULD2H(ISYMBJ,ISYRES) 223 ISYMIK = MULD2H(ISYMJ,ISYIKJ) 224C 225 DO J = 1, NRHF(ISYMJ) 226C 227C----------------------------------------------------------- 228C Contract integrals with trial vector CTR1. 229C----------------------------------------------------------- 230C 231 DO 190 ISYMK = 1,NSYM 232C 233 ISYMI = MULD2H(ISYMK,ISYMIK) 234 ISYMA = MULD2H(ISYMK,ISYMV) 235 NBJ = IT1AM(ISYMB,ISYMJ) 236 * + NVIR(ISYMB)*(J-1) + B 237C 238 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 239 KOFF2 = KSCR1 + NMAIJK(ISYIKJ)*(B-1) 240 * + IMAIJK(ISYMIK,ISYMJ) 241 * + NMATIJ(ISYMIK)*(J-1) 242 * + IMATIJ(ISYMI,ISYMK) 243 KOFF3 = IT2SQ(ISYMAI,ISYMBJ) 244 * + NT1AM(ISYMAI)*(NBJ-1) 245 * + IT1AM(ISYMA,ISYMI) + 1 246C 247 NTOTA = MAX(NVIR(ISYMA),1) 248 NTOTI = MAX(NRHF(ISYMI),1) 249C 250 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI), 251 * NRHF(ISYMK),ONEM,CTR1(KOFF1),NTOTA, 252 * WORK(KOFF2),NTOTI,ONE,RHO2(KOFF3), 253 * NTOTA) 254C 255 190 CONTINUE 256 END DO 257 170 CONTINUE 258 110 CONTINUE 259 100 CONTINUE 260C 261 CALL QEXIT(myname) 262C 263 RETURN 264 END 265C /* Deck cc_xd3 */ 266 SUBROUTINE CC_XD3(XOUT,XMAT,ISYMX,XLAMDH,XLAMDP,ISYMLH, 267 * WORK,LWORK) 268 269C 270C Rasmus Faber 2017 271C 272C Purpose: Transform X-matrix to AO-basis! 273C 274C Based on CC_YD by Asger Halkier 8/12 - 1995. 275C 276 IMPLICIT NONE 277 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, ONE = 1.0D0, 278 & ONEM = -1.0D0 279 CHARACTER(LEN=*), PARAMETER :: myname = 'CC_XD3' 280 DOUBLE PRECISION, INTENT(IN) :: XMAT(*), XLAMDH(*), XLAMDP(*) 281 DOUBLE PRECISION, INTENT(INOUT) :: XOUT(*), WORK(LWORK) 282 INTEGER, INTENT(IN) :: ISYMX, ISYMLH, LWORK 283 284 INTEGER :: ISYMM, ISYMN, ISYMAL, ISYMBE, ISYMXD 285 INTEGER :: NTOTAL, NTOTBE, NTOTM 286 INTEGER :: KOFF1, KOFF2, LWRKCH 287 288C#include "priunit.h" 289#include "ccorb.h" 290#include "ccsdsym.h" 291#include "cclr.h" 292C 293 ISYMXD = MULD2H(ISYMLH,ISYMX) 294C 295C------------------------------------ 296C Transform X-matrix to AO-basis. 297C------------------------------------ 298C 299 DO 100 ISYMAL = 1,NSYM 300C 301 ISYMM = ISYMAL 302 ISYMBE = MULD2H(ISYMAL,ISYMXD) 303 ISYMN = MULD2H(ISYMBE,ISYMLH) 304C 305 LWRKCH = LWORK - NBAS(ISYMAL)*NVIR(ISYMN) 306C 307 IF (LWRKCH .LT. 0) THEN 308 CALL QUIT('Insufficient work space in '//myname) 309 ENDIF 310C 311 KOFF1 = ILMRHF(ISYMM) + 1 312 KOFF2 = IMATIJ(ISYMM,ISYMN) + 1 313C 314 NTOTAL = MAX(NBAS(ISYMAL),1) 315 NTOTM = MAX(NRHF(ISYMM),1) 316C 317C lambdp (alpha,m) * x(m,n) => work(alpha,n) 318 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMN),NRHF(ISYMM), 319 * ONE,XLAMDP(KOFF1),NTOTAL,XMAT(KOFF2),NTOTM, 320 * ZERO,WORK,NTOTAL) 321C 322 KOFF1 = IGLMRH(ISYMBE,ISYMN) + 1 323 KOFF2 = IAODIS(ISYMAL,ISYMBE) + 1 324C 325 NTOTAL = MAX(NBAS(ISYMAL),1) 326 NTOTBE = MAX(NBAS(ISYMBE),1) 327C work(alpha,n) * lambdh(beta,n) => xout(alpha,beta) 328 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NRHF(ISYMN), 329 * ONEM,WORK,NTOTAL,XLAMDH(KOFF1),NTOTBE, 330 * ONE,XOUT(KOFF2),NTOTAL) 331C 332 100 CONTINUE 333C 334 RETURN 335 END 336C 337C 338 SUBROUTINE CC_T2SQTRANSP(X2SQ,ISYMTR) 339C 340C Transpose the squared doubles amplitude-like array, 341C x2 (ai,bj) -> (bj,ai) 342C Needed in the triplet case, when the squared array 343C typically contain a linear combination of (+) and (-) 344C components 345 IMPLICIT NONE 346 347#include "ccorb.h" 348#include "ccsdsym.h" 349 350 INTEGER, INTENT(IN) :: ISYMTR 351 DOUBLE PRECISION, INTENT(INOUT) :: X2SQ(*) 352 INTEGER :: ISYMAI, ISYMBJ, IXPOS, IXPOS2 353C 354 IF (ISYMTR.EQ.1) THEN 355 356 DO ISYMBJ = 1, NSYM 357 ISYMAI = ISYMBJ 358 IXPOS = IT2SQ(ISYMAI,ISYMBJ) + 1 359 360 CALL TRANSP(X2SQ(IXPOS),NT1AM(ISYMAI)) 361 362 END DO 363 364 ELSE 365 366 DO ISYMBJ = 1, NSYM 367 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 368 IF (ISYMAI .GT. ISYMBJ) CYCLE 369 IXPOS = IT2SQ(ISYMAI,ISYMBJ) + 1 370 IXPOS2 = IT2SQ(ISYMBJ,ISYMAI) + 1 371 CALL TRANSPA(X2SQ(IXPOS),X2SQ(IXPOS2), 372 & NT1AM(ISYMAI),NT1AM(ISYMBJ)) 373 END DO 374 375 END IF 376C 377 CONTAINS 378C 379 PURE SUBROUTINE TRANSP(A,N) 380C 381 DOUBLE PRECISION, INTENT(INOUT) :: A(N,N) 382 INTEGER, INTENT(IN) :: N 383 DOUBLE PRECISION :: TMP 384 INTEGER I, J 385C 386 DO I = 1, N 387 DO J = 1, I-1 388 TMP = A(I,J) 389 A(I,J) = A(J,I) 390 A(J,I) = TMP 391 END DO 392 END DO 393C 394 END SUBROUTINE 395C 396 PURE SUBROUTINE TRANSPA(A,B,M,N) 397C 398 DOUBLE PRECISION, INTENT(INOUT) :: A(M,N), B(N,M) 399 INTEGER, INTENT(IN) :: N, M 400 INTEGER :: I, J 401 DOUBLE PRECISION TMP 402C 403 DO I = 1, M 404 DO J = 1, N 405 TMP = A(I,J) 406 A(I,J) = B(J,I) 407 B(J,I) = TMP 408 END DO 409 END DO 410C 411 END SUBROUTINE 412 413 END SUBROUTINE 414C 415 SUBROUTINE CC_T2SQSYMSCAL(X2SQ,ISYMTR,FACT) 416C 417C Scale the symmetric part of the squared doubles array by FACT: 418C x2 (ai,bj) <- (1+FACT)/2 (ai,bj) - (1-FACT)/2 (bj,ai) 419C Needed in the triplet case, when the squared array 420C typically contain a linear combination of (+) and (-) 421C components 422 IMPLICIT NONE 423 424#include "ccorb.h" 425#include "ccsdsym.h" 426 427 INTEGER, INTENT(IN) :: ISYMTR 428 DOUBLE PRECISION, INTENT(INOUT) :: X2SQ(*) 429 DOUBLE PRECISION, INTENT(IN) :: FACT 430 DOUBLE PRECISION, PARAMETER :: HALF = 0.5D0, ONE = 1.0D0 431 INTEGER :: ISYMAI, ISYMBJ, IXPOS, IXPOS2 432 DOUBLE PRECISION :: FN, FT 433C 434 FN = HALF*(ONE+FACT) 435 FT = HALF*(ONE-FACT) 436C 437 IF (ISYMTR.EQ.1) THEN 438 439 DO ISYMBJ = 1, NSYM 440 ISYMAI = ISYMBJ 441 IXPOS = IT2SQ(ISYMAI,ISYMBJ) + 1 442 443 CALL SCAL1(X2SQ(IXPOS),NT1AM(ISYMAI)) 444 445 END DO 446 447 ELSE 448 449 DO ISYMBJ = 1, NSYM 450 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 451 IF (ISYMAI .GT. ISYMBJ) CYCLE 452 IXPOS = IT2SQ(ISYMAI,ISYMBJ) + 1 453 IXPOS2 = IT2SQ(ISYMBJ,ISYMAI) + 1 454 CALL SCAL2(X2SQ(IXPOS),X2SQ(IXPOS2), 455 & NT1AM(ISYMAI),NT1AM(ISYMBJ)) 456 END DO 457 458 END IF 459C 460 CONTAINS 461C 462 PURE SUBROUTINE SCAL1(A,N) 463C 464 DOUBLE PRECISION, INTENT(INOUT) :: A(N,N) 465 INTEGER, INTENT(IN) :: N 466 DOUBLE PRECISION :: TMP1, TMP2 467 INTEGER I, J 468C 469 DO I = 1, N 470 DO J = 1, I-1 471 TMP1 = A(I,J) 472 TMP2 = A(J,I) 473 A(I,J) = FN*TMP1 - FT*TMP2 474 A(J,I) = FN*TMP2 - FT*TMP1 475 END DO 476 END DO 477C 478 END SUBROUTINE 479C 480 PURE SUBROUTINE SCAL2(A,B,M,N) 481C 482 DOUBLE PRECISION, INTENT(INOUT) :: A(M,N), B(N,M) 483 INTEGER, INTENT(IN) :: N, M 484 INTEGER :: I, J 485 DOUBLE PRECISION TMP1, TMP2 486C 487 DO I = 1, M 488 DO J = 1, N 489 TMP1 = A(I,J) 490 TMP2 = B(J,I) 491 A(I,J) = FN*TMP1 - FT*TMP2 492 B(J,I) = FN*TMP2 - FT*TMP1 493 END DO 494 END DO 495C 496 END SUBROUTINE 497 498 END SUBROUTINE 499C 500 SUBROUTINE CCSD_TCMEPK3(T2AM,ISYOPE) 501C 502C Swap I and J in the antisymmetric doubles array 503C T2AM(ai,bj) -> T2AM(aj,bi) 504C Since only the lower half is stored, 505C in practice A and B is swapped instead: 506C T2AM(ai,bj) -> - T2AM(bi,aj) 507C only for i != j 508C 509C Based on CCSD_TCMEPK by 510C Henrik Koch and Alfredo Sanchez. Dec 1994 511C 512#include "implicit.h" 513C 514 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, ONEM = -1.0D0) 515C 516 DIMENSION T2AM(*) 517C 518#include "priunit.h" 519#include "ccorb.h" 520#include "ccsdinp.h" 521#include "ccsdsym.h" 522C 523 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 524C 525 CALL QENTER('CCSD_TCMEPK3') 526C 527 DO 100 ISYMIJ = 1,NSYM 528C 529 ISYMAB = MULD2H(ISYMIJ,ISYOPE) 530C 531 DO 110 ISYMJ = 1,NSYM 532C 533 ISYMI = MULD2H(ISYMJ,ISYMIJ) 534C 535 IF (ISYMI .GT. ISYMJ) GOTO 110 536C 537 DO 120 ISYMB = 1,NSYM 538C 539 ISYMA = MULD2H(ISYMB,ISYMAB) 540C 541 IF (ISYMA .GT. ISYMB) GOTO 120 542C 543 ISYMAI = MULD2H(ISYMA,ISYMI) 544 ISYMBJ = MULD2H(ISYMB,ISYMJ) 545 ISYMBI = MULD2H(ISYMB,ISYMI) 546 ISYMAJ = MULD2H(ISYMA,ISYMJ) 547C 548C 549 DO 130 J = 1,NRHF(ISYMJ) 550C 551 IF (ISYMI .EQ. ISYMJ) THEN 552 NRHFI = J - 1 553 ELSE 554 NRHFI = NRHF(ISYMI) 555 ENDIF 556C 557 DO 140 I = 1,NRHFI 558C 559 DO 150 B = 1,NVIR(ISYMB) 560C 561 IF (ISYMB .EQ. ISYMA) THEN 562 NVIRA = B 563 ELSE 564 NVIRA = NVIR(ISYMA) 565 ENDIF 566C 567 NBI = IT1AM(ISYMB,ISYMI) 568 * + NVIR(ISYMB)*(I - 1) + B 569 NBJ = IT1AM(ISYMB,ISYMJ) 570 * + NVIR(ISYMB)*(J - 1) + B 571C 572 DO 160 A = 1,NVIRA 573C 574 NAI = IT1AM(ISYMA,ISYMI) 575 * + NVIR(ISYMA)*(I - 1) + A 576 NAJ = IT1AM(ISYMA,ISYMJ) 577 * + NVIR(ISYMA)*(J - 1) + A 578C 579 FACT = ONE 580 IF (ISYMAI.EQ.ISYMBJ) THEN 581 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 582 * + INDEX(NAI,NBJ) 583 ELSE IF (ISYMAI.LT.ISYMBJ) THEN 584 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 585 * + NT1AM(ISYMAI)*(NBJ-1)+NAI 586 ELSE IF (ISYMBJ.LT.ISYMAI) THEN 587 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 588 * + NT1AM(ISYMBJ)*(NAI-1)+NBJ 589 FACT = FACT*ONEM 590 END IF 591C 592 IF (ISYMAJ.EQ.ISYMBI) THEN 593 NAJBI = IT2AM(ISYMAJ,ISYMBI) 594 * + INDEX(NAJ,NBI) 595 FACT = FACT*ONEM 596 ELSE IF (ISYMAJ.LT.ISYMBI) THEN 597 NAJBI = IT2AM(ISYMAJ,ISYMBI) 598 * + NT1AM(ISYMAJ)*(NBI-1)+NAJ 599 ELSE IF (ISYMBI.LT.ISYMAJ) THEN 600 NAJBI = IT2AM(ISYMAJ,ISYMBI) 601 * + NT1AM(ISYMBI)*(NAJ-1)+NBI 602 FACT = FACT*ONEM 603 END IF 604C 605 XAIBJ = T2AM(NAJBI) 606 XAJBI = T2AM(NAIBJ) 607C 608 T2AM(NAJBI) = FACT*XAJBI 609 T2AM(NAIBJ) = FACT*XAIBJ 610C 611 160 CONTINUE 612 150 CONTINUE 613 140 CONTINUE 614 130 CONTINUE 615 120 CONTINUE 616 110 CONTINUE 617 100 CONTINUE 618C 619C--------------------------------------- 620C Scale diagonal elements of result. 621C--------------------------------------- 622C 623 IF (ISYOPE .NE. 1) THEN 624 CALL QEXIT('CCSD_TCMEPK3') 625 RETURN 626 ENDIF 627C 628 DO 200 ISYMAI = 1,NSYM 629 DO 210 NAI = 1,NT1AM(ISYMAI) 630 NAIAI = IT2AM(ISYMAI,ISYMAI) + INDEX(NAI,NAI) 631 T2AM(NAIAI) = 0.0D0 632 210 CONTINUE 633 200 CONTINUE 634C 635 CALL QEXIT('CCSD_TCMEPK3') 636C 637 RETURN 638 END 639C /* Deck cc_zwvi3 */ 640 SUBROUTINE CC_ZWVI3(ZINT,CTR2,ISYMC2 ,TINT,ISYTIN, 641 * WORK,LWORK) 642C 643C Written by Asger Halkier 26/10 - 1995 644C 645C Version: 1.0 646C 647C Purpose: To calculate the intermediates entering some of the 648C terms in the 2.1-block. 649C Ove Christiansen 1-10-1996: 650C general symmetry of ctr2 (isymc2) 651C 652#include "implicit.h" 653 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 654 DIMENSION ZINT(*), CTR2(*), TINT(*), WORK(LWORK) 655#include "priunit.h" 656#include "ccorb.h" 657#include "ccsdsym.h" 658#include "cclr.h" 659C 660C----------------------------- 661C Initialize result array. 662C----------------------------- 663C 664 ISYZIN = MULD2H(ISYMC2,ISYTIN) 665C 666 CALL DZERO(ZINT,NT2BCD(ISYZIN)) 667C 668C------------------------ 669C Do the calculation. 670C------------------------ 671C 672 DO 100 ISYMDK = 1,NSYM 673C 674 ISYMEI = MULD2H(ISYMDK,ISYMC2) 675 ISYMJ = MULD2H(ISYMDK,ISYTIN) 676C 677 KOFF1 = IT2SQ(ISYMDK,ISYMEI) + 1 678 KOFF2 = IT2BCD(ISYMDK,ISYMJ) + 1 679 KOFF3 = IT2BCD(ISYMEI,ISYMJ) + 1 680C 681 NTOTEI = MAX(NT1AM(ISYMEI),1) 682 NTOTDK = MAX(NT1AM(ISYMDK),1) 683C 684 CALL DGEMM('T','N',NT1AM(ISYMEI),NRHF(ISYMJ),NT1AM(ISYMDK), 685 * ONE,CTR2(KOFF1),NTOTDK,TINT(KOFF2),NTOTDK,ZERO, 686 * ZINT(KOFF3),NTOTEI) 687C 688 100 CONTINUE 689 690 RETURN 691 END 692 693C /* Deck cc_t2ao3 */ 694 SUBROUTINE CC_T2AO3(T2AM,XLAMDH,ISYMLH,SCRM,WORK,LWORK, 695 * IDEL,ISYMD,ISYMTR) 696C 697C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 698C Written by Rasmus Faber, 2017 699C Based on CC_T2AO by H. Kock, A. Sanchez, O. Christiansen and 700C A. Halkier. 701C PURPOSE: 702C Tdjci -> Tci,j (delta) 703C i.e. transform the other index than CC_T2AO, needed in case 704C of non-symmetric T2AM 705C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 706C 707 implicit none 708#include "priunit.h" 709#include "iratdef.h" 710 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, ONE = 1.0D0 711 CHARACTER(LEN=*), PARAMETER :: myname = 'CC_T2AO3' 712C 713 INTEGER, INTENT(IN) :: ISYMLH, LWORK, IDEL, ISYMD, ISYMTR 714 DOUBLE PRECISION, INTENT(IN) :: T2AM(*),XLAMDH(*) 715 DOUBLE PRECISION, INTENT(OUT) :: SCRM(*) 716 DOUBLE PRECISION, INTENT(INOUT) :: WORK(LWORK) 717C 718#include "ccorb.h" 719#include "ccsdsym.h" 720C 721 INTEGER :: ISYDVI, ISYMM, ISYMCI, ISYMJ, ISYMDJ 722 INTEGER :: KOFF1, KOFF2, KOFF3 723 INTEGER :: NTOTCI, NCI, NTOTDVI 724 725C 726C----------------------------------------------------- 727C Calculate the transformed t2-amplitude and save. 728C----------------------------------------------------- 729C 730 ISYDVI = MULD2H(ISYMLH,ISYMD) 731 ISYMM = MULD2H(ISYMTR,ISYDVI) 732 CALL DZERO(SCRM,NT2BCD(ISYMM)) 733C 734 IF ( LWORK .LT. NVIR(ISYDVI)) THEN 735 CALL QUIT('Insufficient core in '//myname) 736 ENDIF 737C 738 CALL DZERO(WORK,NVIR(ISYDVI)) 739C 740 KOFF1 = IGLMVI(ISYMD,ISYDVI) + IDEL - IBAS(ISYMD) 741C 742 CALL DCOPY(NVIR(ISYDVI),XLAMDH(KOFF1),NBAS(ISYMD),WORK,1) 743C 744 DO 100 ISYMCI = 1,NSYM 745C 746 ISYMDJ = MULD2H(ISYMTR,ISYMCI) 747 ISYMJ = MULD2H(ISYMDJ,ISYDVI) 748C 749 NTOTCI = MAX(1,NT1AM(ISYMCI)) 750 NTOTDVI = MAX(1,NVIR(ISYDVI)) 751C 752 DO 110 NCI = 1, NT1AM(ISYMCI) 753C 754 KOFF2 = IT2SQ(ISYMDJ,ISYMCI) 755 * + NT1AM(ISYMDJ)*(NCI - 1) 756 & + IT1AM(ISYDVI,ISYMJ) + 1 757 KOFF3 = IT2BCD(ISYMCI,ISYMJ) 758 * + NCI 759C 760 CALL DGEMV('T',NVIR(ISYDVI),NRHF(ISYMJ),ONE, 761 * T2AM(KOFF2),NTOTDVI,WORK,1,ZERO, 762 * SCRM(KOFF3),NTOTCI) 763C 764 110 CONTINUE 765 100 CONTINUE 766C 767 RETURN 768 END 769C /* Deck cc_t2mot */ 770 SUBROUTINE CC_T2MOT(RHO1,CTR2,ISYMC2,OMEGA2,RHO2,GAMMA,XLAMDP, 771 * XLAMPC,ISYMPC,WORK,LWORK,ISYMBF,ICON) 772C 773C Rasmus Faber 774C 775C Based on cc_t2mo by Henrik Koch and Alfredo Sanchez. 776C 777C Transform the Omega2 vector from the AO basis to the MO 778C basis. Triplet version. 779C 780C 781 implicit none 782C 783#include "priunit.h" 784#include "maxorb.h" 785 CHARACTER(len=*), PARAMETER :: myname = "CC_T2MOT" 786 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0 787 DOUBLE PRECISION, PARAMETER :: FOURTH = 0.25D0 788 DOUBLE PRECISION, INTENT(IN):: CTR2(*), OMEGA2(*), 789 & GAMMA(*), XLAMDP(*), XLAMPC(*) 790 DOUBLE PRECISION, INTENT(OUT):: RHO1(*), RHO2(*), WORK(*) 791#include "ccorb.h" 792#include "ccsdsym.h" 793#include "symsq.h" 794C 795 INTEGER, INTENT(IN) :: ISYMC2, ISYMPC, ISYMBF, ICON, LWORK 796 797 integer :: index 798 INTEGER :: ISYMI, ISYMJ, ISYAL, ISYBE, ISALBE, ISYALI, ISYBEJ, 799 & ISYMA, ISYMB, ISYMAI, ISYMBJ, ISYMO1, ISYMO2, 800 & ISYMIJ, ISYMAB 801 INTEGER :: NVA, NRA, NVB, NRB, NIJP, NIJM, NABP, NABIJP, NABIJM, 802 & NAB, NAI, NBJ, NBASA, NBASB, NVIRA, NAIBJ 803 INTEGER :: KSCR1, KSCR2, KSCR3, KSCR4, KSCR5 804 INTEGER :: KEND1, LWRK1 805 INTEGER :: KOFF1, KOFF2 806 DOUBLE PRECISION :: FACP, FACM 807C 808 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 809C 810 ISYMO2 = MULD2H(ISYMBF,ISYMPC) 811 ISYMO1 = MULD2H(ISYMO2,ISYMC2) 812C 813 IF ((ICON.EQ.1).OR.(ICON.EQ.2)) THEN 814 CALL DZERO(RHO2,NT2SQ(ISYMO2)) 815 ENDIF 816C 817 DO 100 ISYMJ = 1,NSYM 818 DO 110 ISYMI = 1,NSYM 819C 820 ISYMIJ = MULD2H(ISYMI,ISYMJ) 821 ISALBE = MULD2H(ISYMIJ,ISYMBF) 822 ISYMAB = MULD2H(ISYMIJ,ISYMO2) 823C 824 DO 120 ISYBE = 1,NSYM 825C 826 ISYAL = MULD2H(ISYBE,ISALBE) 827 ISYALI = MULD2H(ISYAL,ISYMI) 828 ISYBEJ = MULD2H(ISYBE,ISYMJ) 829C 830C----------------------------------------------- 831C Dynamic allocation of work space. 832C----------------------------------------------- 833C 834 ISYMA = MULD2H(ISYAL,ISYMPC) 835 NVA = MAX(NVIR(ISYMA),NVIR(ISYAL)) 836 NRA = MAX(NRHF(ISYMA),NRHF(ISYAL)) 837 ISYMB = MULD2H(ISYBE,ISYMPC) 838 NVB = MAX(NVIR(ISYMB),NVIR(ISYBE),NRHF(ISYBE)) 839 NRB = MAX(NRHF(ISYMB),NRHF(ISYBE)) 840C 841 KSCR1 = 1 842 IF (ICON.NE.4) THEN 843 KSCR2 = KSCR1 + NBAS(ISYAL)*NBAS(ISYBE) 844 KSCR3 = KSCR2 + NBAS(ISYAL)*NVB 845 KEND1 = KSCR3 + NVA*NVB 846 ELSE 847 KSCR4 = KSCR1 + NBAS(ISYAL)*NBAS(ISYBE) 848 KEND1 = KSCR4 + NBAS(ISYAL)*NRB 849 END IF 850 LWRK1 = LWORK - KEND1 851C 852 IF (LWRK1 .LT. 0) THEN 853 CALL QUIT('Not enough space in '//myname) 854 END IF 855C 856 DO 130 J = 1,NRHF(ISYMJ) 857 DO 140 I = 1,NRHF(ISYMI) 858C 859C------------------------------------------ 860C Squareup the AB block. 861C------------------------------------------ 862C Plus version is stored for i<j with an implicit 863C symmetry i<j = - i>j 864C 865 IF (ISYMI .EQ. ISYMJ) THEN 866 NIJP = IMIJP(ISYMI,ISYMJ) + INDEX(I,J) 867 FACP = ONE 868 IF (I .EQ. J) FACP = ZERO 869 IF (I .GT. J) FACP = -ONE 870 ELSE IF (ISYMI .LT. ISYMJ) THEN 871 NIJP = IMIJP(ISYMI,ISYMJ) 872 * + NRHF(ISYMI)*(J - 1) + I 873 FACP = ONE 874 ELSE 875 NIJP = IMIJP(ISYMJ,ISYMI) 876 * + NRHF(ISYMJ)*(I - 1) + J 877 FACP = -ONE 878 ENDIF 879C Minus version: i,j stored as square array 880 NIJM = IMATIJ(ISYMI,ISYMJ) 881 & + NRHF(ISYMI)*(J-1) + I 882C 883 DO 166 B = 1, NBAS(ISYBE) 884 DO 167 A = 1, NBAS(ISYAL) 885C 886 IF (ISYAL .EQ. ISYBE) THEN 887 NABP = IAODPK(ISYAL,ISYBE) 888 * + INDEX(A,B) 889 FACM = ONE 890 IF (A .EQ. B) FACM = ZERO 891 IF (A .GT. B) FACM = -ONE 892 ELSE IF (ISYAL .LT. ISYBE) THEN 893 NABP = IAODPK(ISYAL,ISYBE) 894 * + NBAS(ISYAL)*(B - 1) + A 895 FACM = ONE 896 ELSE 897 NABP = IAODPK(ISYBE,ISYAL) 898 * + NBAS(ISYBE)*(A - 1) + B 899 FACM = -ONE 900 ENDIF 901C 902 NABIJP = IT2ORT(ISALBE,ISYMIJ) 903 * + NNBST(ISALBE)*(NIJP - 1) + NABP 904C 905 NABIJM = NT2ORT(ISYMBF) 906 * + IT2ORT3(ISALBE,ISYMIJ) 907 * + NNBST(ISALBE)*(NIJM - 1) + NABP 908C 909 NAB = KSCR1 + NBAS(ISYAL)*(B - 1) + A - 1 910C 911 WORK(NAB) = FOURTH* 912 & (FACP*OMEGA2(NABIJP) + 913 & FACM*OMEGA2(NABIJM)) 914C 915 167 CONTINUE 916 166 CONTINUE 917C 918C------------------------------------------------------------ 919C Transform the AB block to virtual space. 920C------------------------------------------------------------ 921C 922 IF ((ICON.EQ.1).OR.(ICON.EQ.2)) THEN 923C 924 ISYMA = MULD2H(ISYAL,ISYMPC) 925 ISYMB = ISYBE 926 ISYMAI = MULD2H(ISYMA,ISYMI) 927 ISYMBJ = MULD2H(ISYMB,ISYMJ) 928C 929 NBASA = MAX(NBAS(ISYAL),1) 930 NBASB = MAX(NBAS(ISYBE),1) 931 NVIRA = MAX(NVIR(ISYMA),1) 932C 933 KOFF1 = ILMVIR(ISYBE) + 1 934C 935 CALL DGEMM('N','N',NBAS(ISYAL),NVIR(ISYMB), 936 * NBAS(ISYBE),ONE,WORK(KSCR1),NBASA, 937 * XLAMDP(KOFF1),NBASB,ZERO,WORK(KSCR2), 938 * NBASA) 939C 940 KOFF2 = IGLMVI(ISYAL,ISYMA) + 1 941C 942 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB), 943 * NBAS(ISYAL),ONE,XLAMPC(KOFF2),NBASA, 944 * WORK(KSCR2),NBASA,ZERO,WORK(KSCR3), 945 * NVIRA) 946C 947C-------------------------------------------- 948C Store the omega2 vector. 949C-------------------------------------------- 950C 951 DO 170 B = 1,NVIR(ISYMB) 952 NBJ = IT1AM(ISYMB,ISYMJ) 953 * + NVIR(ISYMB)*(J-1) + B 954 DO 180 A = 1,NVIR(ISYMA) 955C 956 NAI = IT1AM(ISYMA,ISYMI) 957 * + NVIR(ISYMA)*(I-1) + A 958 NAB = KSCR3 + NVIR(ISYMA)*(B - 1) + A - 1 959C 960 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) + 961 & NT1AM(ISYMAI)*(NBJ-1) + NAI 962C 963 RHO2(NAIBJ) = RHO2(NAIBJ) + WORK(NAB) 964C 965 180 CONTINUE 966 170 CONTINUE 967C 968 ENDIF 969C 970C-------------------------------------- 971C CCLR contribution. 972C-------------------------------------- 973C 974 IF (ICON .EQ. 2 ) THEN 975C 976 CALL DZERO(WORK(KSCR2),NVA*NVB) 977 ISYMA = ISYAL 978 ISYMB = MULD2H(ISYBE,ISYMPC) 979 ISYMAI = MULD2H(ISYMA,ISYMI) 980 ISYMBJ = MULD2H(ISYMB,ISYMJ) 981C 982 NBASA = MAX(NBAS(ISYAL),1) 983 NBASB = MAX(NBAS(ISYBE),1) 984 NVIRA = MAX(NVIR(ISYMA),1) 985C 986 KOFF1 = IGLMVI(ISYBE,ISYMB) + 1 987C 988 CALL DGEMM('N','N',NBAS(ISYAL),NVIR(ISYMB), 989 * NBAS(ISYBE),ONE,WORK(KSCR1),NBASA, 990 * XLAMPC(KOFF1),NBASB,ZERO,WORK(KSCR2), 991 * NBASA) 992C 993 KOFF2 = ILMVIR(ISYAL) + 1 994C 995 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB), 996 * NBAS(ISYAL),ONE,XLAMDP(KOFF2),NBASA, 997 * WORK(KSCR2),NBASA,ZERO,WORK(KSCR3), 998 * NVIRA) 999C 1000C-------------------------------------------- 1001C Store the omega2 vector. 1002C-------------------------------------------- 1003 DO 171 B = 1,NVIR(ISYMB) 1004 NBJ = IT1AM(ISYMB,ISYMJ) 1005 * + NVIR(ISYMB)*(J-1) + B 1006 DO 181 A = 1,NVIR(ISYMA) 1007C 1008 NAI = IT1AM(ISYMA,ISYMI) 1009 * + NVIR(ISYMA)*(I-1) + A 1010 NAB = KSCR3 + NVIR(ISYMA)*(B - 1) + A - 1 1011C 1012 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) + 1013 & NT1AM(ISYMAI)*(NBJ-1) + NAI 1014C 1015 RHO2(NAIBJ) = RHO2(NAIBJ) + WORK(NAB) 1016C 1017 181 CONTINUE 1018 171 CONTINUE 1019C 1020C 1021 ENDIF 1022C 1023C 1024 999 CONTINUE 1025 140 CONTINUE 1026 130 CONTINUE 1027 120 CONTINUE 1028 110 CONTINUE 1029 100 CONTINUE 1030C 1031 RETURN 1032 END 1033 1034 SUBROUTINE CC_TRIP_EXTRACT(X2SQ,X2AMP,ISYMTR,ANTI) 1035C 1036C Rasmus Faber, Aug 2017 1037C 1038C Extracts the components of the "+" triplet doubles vector from the 1039C square array "X2SQ" and puts it in X2AMP ordered AI<BJ 1040C If anti is true, extract (-) component instead. 1041C 1042 IMPLICIT NONE 1043#include "ccorb.h" 1044#include "ccsdsym.h" 1045C 1046 DOUBLE PRECISION, INTENT(IN) :: X2SQ(*) 1047 DOUBLE PRECISION, INTENT(OUT) :: X2AMP(*) 1048 INTEGER, INTENT(IN) :: ISYMTR 1049 LOGICAL, INTENT(IN) :: ANTI 1050C 1051 INTEGER :: ISYMAI, ISYMBJ 1052 INTEGER :: NAI, NBJ, NAIBJS, NAIBJT, NBJAIS 1053 DOUBLE PRECISION :: FACTOR, XAIBJ, XBJAI 1054C 1055 FACTOR = 1.0D0 1056 IF (ANTI) FACTOR = -1.0D0 1057C 1058 DO ISYMBJ = 1, NSYM 1059 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 1060 1061 IF (ISYMTR.EQ.1) THEN 1062 DO NBJ = 1, NT1AM(ISYMBJ) 1063 DO NAI = 1, NBJ 1064C 1065 NAIBJS = IT2SQ(ISYMAI,ISYMBJ) 1066 & + NT1AM(ISYMAI)*(NBJ-1) + NAI 1067 NBJAIS = IT2SQ(ISYMBJ,ISYMAI) 1068 % + NT1AM(ISYMBJ)*(NAI-1) + NBJ 1069 NAIBJT = IT2AM(ISYMAI,ISYMBJ) 1070 & + NBJ*(NBJ-1)/2 + NAI 1071C 1072 XAIBJ = X2SQ(NAIBJS) 1073 XBJAI = X2SQ(NBJAIS) 1074 X2AMP(NAIBJT) = XAIBJ + FACTOR*XBJAI 1075 END DO 1076 END DO 1077 ELSE IF ( ISYMAI .LT. ISYMBJ) THEN 1078 DO NBJ = 1, NT1AM(ISYMBJ) 1079 DO NAI = 1, NT1AM(ISYMAI) 1080 NAIBJS = IT2SQ(ISYMAI,ISYMBJ) 1081 & + NT1AM(ISYMAI)*(NBJ-1) + NAI 1082 NBJAIS = IT2SQ(ISYMBJ,ISYMAI) 1083 % + NT1AM(ISYMBJ)*(NAI-1) + NBJ 1084 NAIBJT = IT2AM(ISYMAI,ISYMBJ) 1085 & + NT1AM(ISYMAI)*(NBJ-1) + NAI 1086C 1087 XAIBJ = X2SQ(NAIBJS) 1088 XBJAI = X2SQ(NBJAIS) 1089 X2AMP(NAIBJT) = XAIBJ + FACTOR*XBJAI 1090 END DO 1091 END DO 1092 END IF 1093 1094 END DO 1095C 1096 END SUBROUTINE 1097C 1098C /* Deck cc_22am3 */ 1099 SUBROUTINE CC_22AM3(RHO2,XIAJB,XMINT,ISYMIM,WORK,LWORK) 1100C 1101C Written by Rasmus Faber - 2017 1102C Based on cc_22am by Asger Halkier 1103C 1104C Purpose: Contract the M intermediate with the (ma|nb) integrals 1105C to obtain the E contribution to the left transformed doubles 1106C vector. 1107C 1108C 1109 IMPLICIT NONE 1110C 1111 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, ONE = 1.0D0, 1112 & HALF = 0.5D0 1113 DOUBLE PRECISION :: RHO2(*), XIAJB(*), XMINT(*), WORK(LWORK) 1114 INTENT(IN) :: XIAJB 1115 INTENT(INOUT) :: XMINT 1116 INTENT(OUT) :: RHO2, WORK 1117 INTEGER, INTENT(IN) :: ISYMIM, LWORK 1118#include "priunit.h" 1119#include "ccorb.h" 1120#include "ccsdsym.h" 1121#include "cclr.h" 1122C 1123 INTEGER :: INDEX 1124 INTEGER :: ISYMJ, ISYMB, ISYMA, ISYMI, ISYMM, ISYMN, 1125 & ISYRES, ISYMIN, ISYMBJ, ISYMAI, 1126 & ISYMAM, ISYMAN, ISYMMI, ISYMBN 1127 INTEGER :: KSCR1, KSCR2 1128 INTEGER :: LWRK1, LWRK2 1129 INTEGER :: KEND1, KEND2 1130 INTEGER :: KOFF1, KOFF2, KOFF3 1131 INTEGER :: NBJ, NBN, NAM, NAI, NAMBN, NAIBJ 1132 INTEGER :: NTOTA, NTOTM 1133C 1134 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1135C 1136 ISYRES = ISYMIM 1137C 1138 CALL TRANSP_MI(XMINT,ISYRES,WORK(1)) 1139C 1140 DO 100 ISYMJ = 1,NSYM 1141C 1142 ISYMIN = MULD2H(ISYMJ,ISYRES) 1143C 1144 IF (NRHF(ISYMJ) .EQ. 0) GOTO 100 1145C 1146 DO 110 J = 1,NRHF(ISYMJ) 1147C 1148 DO 120 ISYMB = 1,NSYM 1149C 1150 ISYMBJ = MULD2H(ISYMB,ISYMJ) 1151 ISYMAI = MULD2H(ISYMBJ,ISYRES) 1152C 1153 IF (NVIR(ISYMB) .EQ. 0) GOTO 120 1154C 1155C---------------------------------------- 1156C Work space allocation one. 1157C---------------------------------------- 1158C 1159 KSCR1 = 1 1160 KEND1 = KSCR1 + NT1AM(ISYMAI) 1161 LWRK1 = LWORK - KEND1 1162C 1163 IF (LWRK1 .LT. 0) THEN 1164 CALL QUIT('Insufficient work space in CC_22AM3') 1165 ENDIF 1166C 1167 DO 130 B = 1,NVIR(ISYMB) 1168C 1169 CALL DZERO(WORK(KSCR1),NT1AM(ISYMAI)) 1170C 1171 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 1172C 1173 DO 140 ISYMN = 1,NSYM 1174C 1175 ISYMBN = MULD2H(ISYMN,ISYMB) 1176 ISYMAM = MULD2H(ISYMBN,ISYMOP) 1177 ISYMMI = MULD2H(ISYMN,ISYMIN) 1178C 1179C---------------------------------------------- 1180C Work space allocation two. 1181C---------------------------------------------- 1182C 1183 KSCR2 = KEND1 1184 KEND2 = KSCR2 + NT1AM(ISYMAM) 1185 LWRK2 = LWORK - KEND2 1186C 1187 IF (LWRK2 .LT. 0) THEN 1188 CALL QUIT('Insufficient work space in CC_22AM3') 1189 ENDIF 1190C 1191 DO 150 N = 1,NRHF(ISYMN) 1192C 1193 NBN = IT1AM(ISYMB,ISYMN) 1194 & + NVIR(ISYMB)*(N - 1) + B 1195C 1196C------------------------------------------------------- 1197C Copy submatrix out of integrals. 1198C------------------------------------------------------- 1199C 1200 DO 160 NAM = 1,NT1AM(ISYMAM) 1201C 1202 NAMBN = IT2AM(ISYMAM,ISYMBN) 1203 & + INDEX(NAM,NBN) 1204C 1205 WORK(KSCR2 + NAM - 1) = XIAJB(NAMBN) 1206C 1207 160 CONTINUE 1208C 1209C-------------------------------------------------------------------- 1210C Contraction of integrals with M-intermediate. 1211C-------------------------------------------------------------------- 1212C 1213 DO 170 ISYMA = 1,NSYM 1214C 1215 ISYMM = MULD2H(ISYMA,ISYMAM) 1216 ISYMI = MULD2H(ISYMA,ISYMAI) 1217C 1218 KOFF1 = KSCR2 + IT1AM(ISYMA,ISYMM) 1219 KOFF2 = I3ORHF(ISYMIN,ISYMJ) 1220 & + NMAIJK(ISYMIN)*(J - 1) 1221 & + IMAIJK(ISYMMI,ISYMN) 1222 & + NMATIJ(ISYMMI)*(N - 1) 1223 & + IMATIJ(ISYMM,ISYMI) + 1 1224 KOFF3 = KSCR1 + IT1AM(ISYMA,ISYMI) 1225C 1226 NTOTA = MAX(NVIR(ISYMA),1) 1227 NTOTM = MAX(NRHF(ISYMM),1) 1228C 1229 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 1230 & NRHF(ISYMM),ONE,WORK(KOFF1), 1231 & NTOTA,XMINT(KOFF2),NTOTM,ONE, 1232 & WORK(KOFF3),NTOTA) 1233C 1234 170 CONTINUE 1235C 1236 150 CONTINUE 1237 140 CONTINUE 1238C 1239C----------------------------------------------- 1240C Storage in RHO2 result vector. 1241C----------------------------------------------- 1242C 1243 DO NAI = 1, NT1AM(ISYMAI) 1244 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) 1245 & + NT1AM(ISYMAI)*(NBJ-1) + NAI 1246 RHO2(NAIBJ) = RHO2(NAIBJ) 1247 & + HALF*WORK(KSCR1+NAI-1) 1248 END DO 1249C 1250 130 CONTINUE 1251 120 CONTINUE 1252 110 CONTINUE 1253 100 CONTINUE 1254C 1255 RETURN 1256C 1257 CONTAINS 1258 SUBROUTINE TRANSP_MI(MINT,ISYMINT,WORK) 1259C 1260C Transposes the M intermediate: 1261C M(M,I,N,J) -> M(N,J,M,I) 1262C 1263 DOUBLE PRECISION, INTENT(INOUT) :: MINT(*), WORK(*) 1264 INTEGER, INTENT(IN) :: ISYMINT 1265 1266 INTEGER :: I, J, N 1267 INTEGER :: ISYMI, ISYMJ, ISYMM, ISYMN, ISYMIN, ISYNJM, 1268 & ISYMMI, ISYMNJ, ISYMMN 1269 INTEGER :: KOFF1, KOFF2 1270 INTEGER :: NTOTI, NTOTM 1271 1272 DO ISYMJ = 1, NSYM 1273 ISYMIN = MULD2H(ISYMJ,ISYMINT) 1274 DO ISYMI = 1, ISYMJ 1275 ISYMMN = MULD2H(ISYMI,ISYMIN) 1276 ISYNJM = MULD2H(ISYMI,ISYMINT) 1277 NTOTI = NRHF(ISYMI) 1278 DO J = 1, NRHF(ISYMJ) 1279 IF (ISYMI.EQ.ISYMJ) NTOTI = J - 1 1280 DO I = 1, NTOTI 1281 DO ISYMN = 1, NSYM 1282 ISYMM = MULD2H(ISYMN,ISYMMN) 1283 ISYMMI = MULD2H(ISYMIN,ISYMN) 1284 ISYMNJ = MULD2H(ISYNJM,ISYMM) 1285 DO N = 1, NRHF(ISYMN) 1286 KOFF1 = I3ORHF(ISYMIN,ISYMJ) 1287 & + NMAIJK(ISYMIN)*(J-1) 1288 & + IMAIJK(ISYMMI,ISYMN) 1289 & + NMATIJ(ISYMMI)*(N-1) 1290 & + IMATIJ(ISYMM,ISYMI) 1291 & + NRHF(ISYMM)*(I-1) + 1 1292 1293 KOFF2 = I3ORHF(ISYNJM,ISYMI) 1294 & + NMAIJK(ISYNJM)*(I-1) 1295 & + IMAIJK(ISYMNJ,ISYMM) 1296 & + IMATIJ(ISYMN,ISYMJ) 1297 & + NRHF(ISYMN)*(J-1) + N 1298 1299 CALL DCOPY(NRHF(ISYMM), 1300 & MINT(KOFF2),NMATIJ(ISYMNJ), 1301 & WORK(1),1) 1302 1303 CALL DCOPY(NRHF(ISYMM), 1304 & MINT(KOFF1),1, 1305 & MINT(KOFF2),NMATIJ(ISYMNJ)) 1306 1307 CALL DCOPY(NRHF(ISYMM), 1308 & WORK(1),1, 1309 & MINT(KOFF1),1) 1310 END DO ! N 1311 END DO ! ISYMN 1312 END DO ! I 1313! 1314! Handle I = J case 1315 IF (ISYMI.EQ.ISYMJ) THEN 1316 I = J 1317 DO ISYMN = 1, NSYM 1318 ISYMM = MULD2H(ISYMMN,ISYMN) 1319 IF (ISYMM.GT.ISYMN) CYCLE 1320 NTOTM = NRHF(ISYMM) 1321 ISYMMI = MULD2H(ISYMIN,ISYMN) 1322 ISYMNJ = MULD2H(ISYNJM,ISYMM) 1323 DO N = 1, NRHF(ISYMN) 1324 IF (ISYMM.EQ.ISYMN) NTOTM = N - 1 1325 KOFF1 = I3ORHF(ISYMIN,ISYMJ) 1326 & + NMAIJK(ISYMIN)*(J-1) 1327 & + IMAIJK(ISYMMI,ISYMN) 1328 & + NMATIJ(ISYMMI)*(N-1) 1329 & + IMATIJ(ISYMM,ISYMI) 1330 & + NRHF(ISYMM)*(I-1) + 1 1331 1332 KOFF2 = I3ORHF(ISYNJM,ISYMI) 1333 & + NMAIJK(ISYNJM)*(I-1) 1334 & + IMAIJK(ISYMNJ,ISYMM) 1335 & + IMATIJ(ISYMN,ISYMJ) 1336 & + NRHF(ISYMN)*(J-1) + N 1337 1338 CALL DCOPY(NTOTM, 1339 & MINT(KOFF2),NMATIJ(ISYMNJ), 1340 & WORK(1),1) 1341 1342 CALL DCOPY(NTOTM, 1343 & MINT(KOFF1),1, 1344 & MINT(KOFF2),NMATIJ(ISYMNJ)) 1345 1346 CALL DCOPY(NTOTM, 1347 & WORK(1),1, 1348 & MINT(KOFF1),1) 1349 END DO 1350 END DO 1351 END IF 1352 END DO ! J 1353 END DO ! ISYMI 1354 END DO ! ISYMJ 1355 1356 END SUBROUTINE 1357 END 1358C 1359 SUBROUTINE CCRHS_ASQ(OMEGA2,T2AM,GAMMA,WORK,LWORK,ISYGAM,ISYVEC, 1360 * IOPT) 1361C 1362C Rasmus Faber 2017 1363C Based on CCRHS_A by Henrik Koch & Ove Christiansen 1364C 1365C Purpose: Calculates the doubles A term 1366C Contract the "Gamma" intermediate with the quadratic 1367C trial-vector "T2AM". 1368C Use 1369C IOPT = 1 for right transformation 1370C IOPT = 2 for left transformation 1371C 1372C This version uses a quadratic storage on the output vector, OMEGA2 1373C 1374 IMPLICIT NONE 1375C 1376 DOUBLE PRECISION, PARAMETER :: ZERO=0.0D0, ONE=1.0D0, HALF=0.5D0 1377 CHARACTER(LEN=*), PARAMETER :: myname = "CCRHS_ASQ" 1378C 1379 DOUBLE PRECISION, INTENT(INOUT) :: OMEGA2(*) 1380 DOUBLE PRECISION, INTENT(IN) :: T2AM(*), GAMMA(*) 1381 DOUBLE PRECISION, INTENT(OUT) :: WORK(LWORK) 1382 INTEGER, INTENT(IN) :: LWORK, ISYGAM, ISYVEC, IOPT 1383#include "priunit.h" 1384#include "ccorb.h" 1385#include "ccsdsym.h" 1386C 1387 INTEGER :: INDEX 1388 INTEGER :: ISYMA, ISYMB, ISYMI, ISYMJ, ISYMK, ISYML, 1389 & ISYMAI, ISYMAK, ISYMBJ, ISYMBL, ISYMKI, ISYMLJ, ISAIBJ 1390 INTEGER :: NAI, NBJ, NBL, NLJ, NKI, NAIBJ, NKILJ, NSTO 1391 INTEGER :: NRHFK, NVIRA 1392 INTEGER :: KSCR1, KSCR2 1393 INTEGER :: KEND1, KEND2 1394 INTEGER :: LWRK1, LWRK2 1395 INTEGER :: KOFF1, KOFF2, KOFF3 1396C 1397 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1398C 1399C---------------------------- 1400C Calculate contribution. 1401C---------------------------- 1402C 1403 ISAIBJ = MULD2H(ISYGAM,ISYVEC) 1404C 1405 DO 100 ISYMLJ = 1,NSYM 1406C 1407 ISYMKI = MULD2H(ISYMLJ,ISYGAM) 1408C 1409 KSCR1 = 1 1410 KEND1 = KSCR1 + NMATIJ(ISYMKI) 1411 LWRK1 = LWORK - KEND1 1412C 1413 IF (LWRK1 .LT. 0) THEN 1414 CALL QUIT('Insufficient space for allocation in '//myname) 1415 END IF 1416C 1417 DO 110 ISYMJ = 1,NSYM 1418C 1419 ISYML = MULD2H(ISYMJ,ISYMLJ) 1420C 1421 DO 120 J = 1,NRHF(ISYMJ) 1422C 1423 DO 130 L = 1,NRHF(ISYML) 1424C 1425 IF (IOPT .EQ. 1) THEN 1426C 1427 NLJ = IMATIJ(ISYML,ISYMJ) 1428 * + NRHF(ISYML)*(J - 1) + L 1429C 1430 ELSE IF (IOPT .EQ. 2) THEN 1431C 1432 NLJ = IMATIJ(ISYMJ,ISYML) 1433 * + NRHF(ISYMJ)*(L - 1) + J 1434C 1435 ENDIF 1436C 1437 DO 140 ISYMK = 1,NSYM 1438C 1439 ISYMI = MULD2H(ISYMK,ISYMKI) 1440C 1441 DO 150 I = 1,NRHF(ISYMI) 1442C 1443 DO 160 K = 1,NRHF(ISYMK) 1444C 1445 IF (IOPT .EQ. 1) THEN 1446C 1447 NKI = IMATIJ(ISYMK,ISYMI) 1448 * + NRHF(ISYMK)*(I - 1) + K 1449C 1450 ELSE IF (IOPT .EQ. 2) THEN 1451C 1452 NKI = IMATIJ(ISYMI,ISYMK) 1453 * + NRHF(ISYMI)*(K - 1) + I 1454C 1455 ENDIF 1456C 1457 IF (ISYMKI .EQ. ISYMLJ) THEN 1458 NKILJ = IGAMMA(ISYMKI,ISYMLJ) 1459 * + INDEX(NKI,NLJ) 1460 ELSE 1461 IF (ISYMKI .LT. ISYMLJ) THEN 1462 NKILJ = IGAMMA(ISYMKI,ISYMLJ) 1463 * + NMATIJ(ISYMKI)*(NLJ - 1) + NKI 1464 ELSE 1465 NKILJ = IGAMMA(ISYMLJ,ISYMKI) 1466 * + NMATIJ(ISYMLJ)*(NKI - 1) + NLJ 1467 ENDIF 1468 ENDIF 1469C 1470 NSTO = IMATIJ(ISYMK,ISYMI) 1471 * + NRHF(ISYMK)*(I - 1) + K 1472C 1473 WORK(KSCR1 + NSTO - 1) = GAMMA(NKILJ) 1474C 1475 160 CONTINUE 1476 150 CONTINUE 1477 140 CONTINUE 1478C 1479 DO 170 ISYMB = 1,NSYM 1480C 1481 ISYMBJ = MULD2H(ISYMB,ISYMJ) 1482 ISYMAI = MULD2H(ISYMBJ,ISAIBJ) 1483 ISYMBL = MULD2H(ISYMB,ISYML) 1484 ISYMAK = MULD2H(ISYVEC,ISYMBL) 1485C 1486 KSCR2 = KEND1 1487 KEND2 = KSCR2 + NT1AM(ISYMAI) 1488 LWRK2 = LWORK - KEND2 1489C 1490 IF (LWRK2 .LT. 0) THEN 1491 CALL QUIT('Insufficient space in '//myname) 1492 END IF 1493C 1494C IF (ISYMAI .GT. ISYMBJ) GOTO 170 1495C 1496 DO 180 B = 1,NVIR(ISYMB) 1497C 1498 NBJ = IT1AM(ISYMB,ISYMJ) 1499 * + NVIR(ISYMB)*(J - 1) + B 1500 NBL = IT1AM(ISYMB,ISYML) 1501 * + NVIR(ISYMB)*(L - 1) + B 1502C 1503 CALL DZERO(WORK(KSCR2),NT1AM(ISYMAI)) 1504C 1505 DO 190 ISYMI = 1,NSYM 1506C 1507 ISYMK = MULD2H(ISYMI,ISYMKI) 1508 ISYMA = MULD2H(ISYMK,ISYMAK) 1509C 1510 NVIRA = MAX(NVIR(ISYMA),1) 1511 NRHFK = MAX(NRHF(ISYMK),1) 1512C 1513 KOFF1 = IT2SQ(ISYMAK,ISYMBL) 1514 * + NT1AM(ISYMAK)*(NBL - 1) 1515 * + IT1AM(ISYMA,ISYMK) + 1 1516 KOFF2 = KSCR1 + IMATIJ(ISYMK,ISYMI) 1517 KOFF3 = KSCR2 + IT1AM(ISYMA,ISYMI) 1518C 1519 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 1520 * NRHF(ISYMK),ONE,T2AM(KOFF1), 1521 * NVIRA,WORK(KOFF2),NRHFK,ZERO, 1522 * WORK(KOFF3),NVIRA) 1523C 1524 190 CONTINUE 1525C 1526 DO 200 NAI = 1, NT1AM(ISYMAI) 1527C 1528 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) 1529 & + NT1AM(ISYMAI)*(NBJ-1) + NAI 1530C 1531 OMEGA2(NAIBJ) = OMEGA2(NAIBJ) 1532 * + HALF*WORK(KSCR2 + NAI - 1) 1533C 1534 200 CONTINUE 1535C 1536 180 CONTINUE 1537 170 CONTINUE 1538C 1539 130 CONTINUE 1540 120 CONTINUE 1541 110 CONTINUE 1542 100 CONTINUE 1543C 1544 RETURN 1545 END 1546C /* Deck cc_22e3 */ 1547 SUBROUTINE CC_22E3(RHO2,T2AM,XMAT,YMAT,ISYMXY,WORK,LWORK) 1548C 1549C Rasmus Faber - 2017 1550C Based on CC_22EC by Asger Halkier 21/11 - 1995 1551C 1552C Purpose: To calculate the triplet doubles F-contribution 1553C to RHO2. 1554C 1555C 1556C 1557 IMPLICIT NONE 1558C 1559 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, ONE = 1.0D0, 1560 & ONEM = -ONE 1561 DOUBLE PRECISION, INTENT(INOUT) :: RHO2(*) 1562 DOUBLE PRECISION, INTENT(IN) :: T2AM(*), XMAT(*), YMAT(*) 1563 DOUBLE PRECISION, INTENT(OUT) :: WORK(LWORK) 1564 INTEGER, INTENT(IN) :: LWORK, ISYMXY 1565#include "priunit.h" 1566#include "ccorb.h" 1567#include "ccsdsym.h" 1568#include "cclr.h" 1569C 1570 INTEGER :: INDEX 1571 INTEGER :: ISYMB, ISYMJ, ISYMAI, ISYMAN, ISYMAT, ISYMAX, ISYMAY, 1572 & ISYMBJ, ISYMFI, ISYMFY, ISYMIX, ISYMIY, ISYMNX, ISYRES 1573 INTEGER :: KOFF1, KOFF2, KOFF3, KOFF4, KOFF5, KOFF6 1574 INTEGER :: NBJ, NAN, NANBJ 1575 INTEGER :: NTOTA, NTOTF, NTOTN 1576 INTEGER :: KSCR1, KSCR2, KEND1, LWRK1 1577C 1578 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1579C 1580 ISYRES = MULD2H(ISYMOP,ISYMXY) 1581 ISYMAT = ISYRES 1582C 1583 DO 100 ISYMJ = 1,NSYM 1584C 1585 IF (NRHF(ISYMJ) .EQ. 0) GOTO 100 1586C 1587 DO 110 J = 1,NRHF(ISYMJ) 1588C 1589 DO 120 ISYMB = 1,NSYM 1590C 1591 ISYMBJ = MULD2H(ISYMB,ISYMJ) 1592 ISYMAI = MULD2H(ISYMBJ,ISYRES) 1593 ISYMAN = MULD2H(ISYMBJ,ISYMOP) 1594 ISYMFI = ISYMAN 1595C 1596 IF (NVIR(ISYMB) .EQ. 0) GOTO 120 1597C 1598C------------------------------------ 1599C Work space allocation. 1600C------------------------------------ 1601C 1602 KSCR1 = 1 1603 KSCR2 = KSCR1 + NT1AM(ISYMAI) 1604 KEND1 = KSCR2 + NT1AM(ISYMAN) 1605 LWRK1 = LWORK - KEND1 1606C 1607 IF (LWRK1 .LT. 0) THEN 1608 CALL QUIT('Insufficient work space in CC_22E3') 1609 ENDIF 1610C 1611 DO 130 B = 1,NVIR(ISYMB) 1612C 1613 CALL DZERO(WORK(KSCR1),NT1AM(ISYMAI)) 1614C 1615 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 1616C 1617C------------------------------------------------- 1618C Copy submatrix out of integrals. 1619C------------------------------------------------- 1620C 1621 DO 140 NAN = 1,NT1AM(ISYMAN) 1622C 1623 NANBJ = IT2AM(ISYMAN,ISYMBJ) + INDEX(NAN,NBJ) 1624C 1625 WORK(KSCR2 + NAN - 1) = T2AM(NANBJ) 1626C 1627 140 CONTINUE 1628C 1629C----------------------------------------------------------------- 1630C Contraction of integrals with X- and Y-matrices. 1631C----------------------------------------------------------------- 1632C 1633 DO 150 ISYMAX = 1,NSYM 1634C 1635 ISYMNX = MULD2H(ISYMAX,ISYMAN) 1636 ISYMIX = MULD2H(ISYMNX,ISYMAT) 1637 ISYMFY = ISYMAX 1638 ISYMIY = MULD2H(ISYMFY,ISYMFI) 1639 ISYMAY = MULD2H(ISYMFY,ISYMAT) 1640C 1641 KOFF1 = KSCR2 + IT1AM(ISYMAX,ISYMNX) 1642 KOFF2 = IMATIJ(ISYMNX,ISYMIX) + 1 1643 KOFF3 = KSCR1 + IT1AM(ISYMAX,ISYMIX) 1644C 1645 NTOTA = MAX(NVIR(ISYMAX),1) 1646 NTOTN = MAX(NRHF(ISYMNX),1) 1647C 1648 CALL DGEMM('N','N',NVIR(ISYMAX),NRHF(ISYMIX), 1649 & NRHF(ISYMNX),ONE,WORK(KOFF1),NTOTA, 1650 & XMAT(KOFF2),NTOTN,ONE,WORK(KOFF3), 1651 & NTOTA) 1652C 1653 KOFF4 = IMATAB(ISYMFY,ISYMAY) + 1 1654 KOFF5 = KSCR2 + IT1AM(ISYMFY,ISYMIY) 1655 KOFF6 = KSCR1 + IT1AM(ISYMAY,ISYMIY) 1656C 1657 NTOTF = MAX(NVIR(ISYMFY),1) 1658 NTOTA = MAX(NVIR(ISYMAY),1) 1659C 1660 CALL DGEMM('T','N',NVIR(ISYMAY),NRHF(ISYMIY), 1661 & NVIR(ISYMFY),ONE,YMAT(KOFF4),NTOTF, 1662 & WORK(KOFF5),NTOTF,ONE,WORK(KOFF6), 1663 & NTOTA) 1664C 1665 150 CONTINUE 1666C 1667C----------------------------------------------- 1668C Storage in RHO2 result vector. 1669C----------------------------------------------- 1670C Try to rearrange the above DGEMMs to 1671C do it there directly 1672 KOFF1 = IT2SQ(ISYMAI,ISYMBJ) 1673 & + NT1AM(ISYMAI)*(NBJ-1) + 1 1674 CALL DAXPY(NT1AM(ISYMAI),ONEM,WORK(KSCR1),1, 1675 & RHO2(KOFF1),1) 1676 1677 130 CONTINUE 1678 120 CONTINUE 1679 110 CONTINUE 1680 100 CONTINUE 1681C 1682 RETURN 1683 END 1684C /* Deck cc_22cd3 */ 1685 SUBROUTINE CC_22CD3(RHO2,CTR2,ISYVEC,XLAMDH,WORK,LWORK, 1686 * ISYDIM,LUF,FIL,IV,IOPT,FACT) 1687C 1688C Written by Rasmus Faber, 2017 1689C Based on CC_22D by Asger Halkier 1690C 1691C Purpose: To calculate the triplet 22C or 22D contribution to RHO2. 1692C If 22D: CTR2 needs to be transposed before this routine. 1693C If 22C: RHO2 must be transposed. 1694C 1695C if iopt = 1 the D intermediate is assumed 1696C to be as in energy calc. 1697C 1698C if iopt ne. 1 we use the intermediate 1699C on lud with address given according to 1700C transformed vector nr iv (iv is not 1 1701C if several vectors are transformed 1702C simultaneously.) 1703C 1704C in Lambda vector calc: iv=1,iopt=1 1705C 1706 IMPLICIT NONE 1707 CHARACTER(LEN=*), PARAMETER :: myname = 'CC_22CD3' 1708 DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0, ONE = 1.0D0 1709 DOUBLE PRECISION, INTENT(INOUT) :: RHO2(*), WORK(LWORK) 1710 DOUBLE PRECISION, INTENT(IN) :: CTR2(*), XLAMDH(*), FACT 1711 CHARACTER, INTENT(IN) :: FIL*(*) 1712 INTEGER, INTENT(IN) :: ISYVEC, LWORK, ISYDIM, LUF, IV, IOPT 1713C 1714#include "priunit.h" 1715#include "ccorb.h" 1716#include "ccsdsym.h" 1717#include "cclr.h" 1718#include "maxorb.h" 1719#include "ccsdio.h" 1720C 1721 INTEGER :: ISYMB, ISYMD, ISYMJ, ISYMAI, ISYMBJ, ISYMEM, 1722 & ISYAIJ, ISYJEM, ISYRES 1723 INTEGER :: IDEL, ID, IOFF, IOFFD, LEN 1724 INTEGER :: NTOTD, NTOTJ, NTOTAI, NTOTEM 1725 INTEGER :: KOFF1, KOFF2, KOFF3, KOFF4, KOFF5 1726 INTEGER :: KPINT, KSCRN, KLAMD, KEND1 1727 INTEGER :: LWRK1 1728 1729 1730C 1731 ISYRES = MULD2H(ISYVEC,ISYDIM) 1732C 1733 DO 100 ISYMD = 1,NSYM 1734C 1735 ISYMB = ISYMD 1736 ISYJEM = MULD2H(ISYMB,ISYDIM) 1737 ISYAIJ = MULD2H(ISYMB,ISYRES) 1738C 1739C------------------------------ 1740C Work space allocation. 1741C------------------------------ 1742C 1743 KPINT = 1 1744 KSCRN = KPINT + NT2BCD(ISYJEM) 1745 KLAMD = KSCRN + NT2BCD(ISYAIJ) 1746 KEND1 = KLAMD + NVIR(ISYMB) 1747 LWRK1 = LWORK - KEND1 1748C 1749 IF (LWRK1 .LT. 0) THEN 1750 CALL QUIT('Insufficient work space in '//myname) 1751 ENDIF 1752C 1753 LEN = NT2BCD(ISYJEM) 1754 IF (LEN .LE. 0) CYCLE ! Nothing to do for this sym 1755 1756 DO 110 IDEL = 1,NBAS(ISYMD) 1757C 1758C----------------------------------------- 1759C Copy row "IDEL" out of XLAMDH. 1760C----------------------------------------- 1761C 1762 IOFFD = ILMVIR(ISYMB) + IDEL 1763 NTOTD = MAX(NBAS(ISYMD),1) 1764 CALL DCOPY(NVIR(ISYMB),XLAMDH(IOFFD),NTOTD, 1765 & WORK(KLAMD),1) 1766 1767C 1768C--------------------------------------------------- 1769C Read P-intermediate submatrix from disc. 1770C--------------------------------------------------- 1771C 1772 ID = IDEL + IBAS(ISYMD) 1773C 1774 IF (IOPT .EQ. 1) THEN 1775 IOFF = IT2DEL(ID) + 1 1776 ELSE 1777 IOFF = IT2DLR(ID,IV) + 1 1778 ENDIF 1779C 1780 CALL GETWA2(LUF,FIL,WORK(KPINT),IOFF,LEN) 1781C 1782 DO 120 ISYMJ = 1,NSYM 1783C 1784C-------------------------------------------------------- 1785C Contraction of intermediate and trial vector. 1786C-------------------------------------------------------- 1787C 1788 ISYMAI = MULD2H(ISYMJ,ISYAIJ) 1789 ISYMEM = MULD2H(ISYMAI,ISYVEC) 1790 ISYMBJ = MULD2H(ISYMJ,ISYMB) 1791C 1792 KOFF1 = IT2SQ(ISYMEM,ISYMAI) + 1 1793 KOFF2 = KPINT + IT2BCT(ISYMJ,ISYMEM) 1794 KOFF3 = KSCRN + IT2BCD(ISYMAI,ISYMJ) 1795C 1796 NTOTAI = MAX(NT1AM(ISYMAI),1) 1797 NTOTEM = MAX(NT1AM(ISYMEM),1) 1798 NTOTJ = MAX(NRHF(ISYMJ),1) 1799C 1800 CALL DGEMM('T','T',NT1AM(ISYMAI),NRHF(ISYMJ), 1801 & NT1AM(ISYMEM),ONE,CTR2(KOFF1),NTOTEM, 1802 & WORK(KOFF2),NTOTJ,ZERO,WORK(KOFF3),NTOTAI) 1803C 1804C----------------------------------------------------------------- 1805C Final scaling with XLAMDH-element and storage in RHO2. 1806C----------------------------------------------------------------- 1807C 1808 DO J = 1, NRHF(ISYMJ) 1809 1810 KOFF4 = KOFF3 + NT1AM(ISYMAI)*(J-1) 1811 1812 KOFF5 = IT2SQ(ISYMAI,ISYMBJ) + 1813 & NT1AM(ISYMAI)* (IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1)) 1814 & + 1 1815 1816 CALL DGER(NT1AM(ISYMAI),NVIR(ISYMB),FACT, 1817 & WORK(KOFF4),1,WORK(KLAMD),1, 1818 & RHO2(KOFF5),NTOTAI) 1819 1820 END DO 1821C 1822 120 CONTINUE 1823 1824 110 CONTINUE 1825 100 CONTINUE 1826C 1827 RETURN 1828 END 1829C /* Deck cc_L1FCK3 */ 1830 SUBROUTINE CC_L1FCK3(RHO2,C1AM,FCKMO,ISYMV,ISYMF,WORK,LWORK) 1831C 1832C Rasmus Faber 2017, 1833C based on CC_L1FCK by Asger Halkier 1834C 1835C Purpose: To calculate the disconnected (first) term 1836C of the H contribution to RHO2. Triplet version 1837C 1838C rho2(ai,bj) = L1am(ai)*FCKMO(jb) 1839C 1840 IMPLICIT NONE 1841 DOUBLE PRECISION, PARAMETER ::ONE = 1.0D0, TWO = 2.0D0 1842 CHARACTER(LEN=*), PARAMETER :: myname = 'CC_L1FCK3' 1843C 1844 DOUBLE PRECISION, INTENT(INOUT) :: RHO2(*), WORK(LWORK) 1845 DOUBLE PRECISION, INTENT(IN) :: C1AM(*), FCKMO(*) 1846 INTEGER, INTENT(IN) :: ISYMV, ISYMF, LWORK 1847#include "priunit.h" 1848#include "ccorb.h" 1849#include "ccsdsym.h" 1850#include "cclr.h" 1851 INTEGER :: ISYMC, ISYMK, ISYMAI, ISYMBJ, ISYRES 1852 INTEGER :: KOFF1, KOFF2, KOFF3, NTOTAI 1853C 1854 CALL QENTER(myname) 1855C 1856 ISYRES = MULD2H(ISYMV,ISYMF) 1857C 1858 IF (LWORK .LT. NT1AM(ISYMF)) THEN 1859 CALL QUIT('Insufficient work-space-area in '//myname) 1860 ENDIF 1861C 1862C----------------------------------- 1863C Extract the fock matrix F(jb). 1864C----------------------------------- 1865C 1866 DO 50 ISYMC = 1,NSYM 1867C 1868 ISYMK = MULD2H(ISYMC,ISYMF) 1869C 1870 DO 60 K = 1,NRHF(ISYMK) 1871C 1872 DO 70 C = 1,NVIR(ISYMC) 1873C 1874 KOFF1 = IFCVIR(ISYMK,ISYMC) + NORB(ISYMK)*(C - 1) + K 1875 KOFF2 = IT1AM(ISYMC,ISYMK) + NVIR(ISYMC)*(K - 1) + C 1876C 1877 WORK(KOFF2) = FCKMO(KOFF1) 1878C 1879 70 CONTINUE 1880 60 CONTINUE 1881 50 CONTINUE 1882C 1883C-------------------------------------------------------------- 1884C Calculate the contraction of C1AM with FCKMO. 1885C-------------------------------------------------------------- 1886C 1887 ISYMAI = ISYMV 1888 ISYMBJ = ISYMF 1889 NTOTAI = MAX(1,NT1AM(ISYMAI)) 1890 1891 KOFF3 = IT2SQ(ISYMAI,ISYMBJ) + 1 1892 1893 CALL DGER(NT1AM(ISYMAI),NT1AM(ISYMBJ),ONE,C1AM,1,WORK,1, 1894 & RHO2(KOFF3),NTOTAI) 1895C 1896C All done 1897 CALL QEXIT(myname) 1898C 1899 RETURN 1900 END 1901 1902C /* Deck cc_L1FCK3P */ 1903 SUBROUTINE CC_L1FCK3P(RHO2,C1AM,FCKMO,ISYMV,ISYMF) 1904C 1905C Rasmus Faber 2017, 1906C based on CC_L1FCK by Asger Halkier 1907C 1908C Purpose: To calculate the disconnected (first) term 1909C of the H contribution to RHO2. 1910C Triplet version, with packed RHO2 1911C 1912C rho2+(ai,bj) = L1am(ai)*FCKMO(jb) + FCKMO(ai)*L1am(bj) 1913C rho2-(ai,bj) = L1am(ai)*FCKMO(jb) - FCKMO(ai)*L1am(bj) 1914C 1915 IMPLICIT NONE 1916 CHARACTER(LEN=*), PARAMETER :: myname = 'CC_L1FCK3P' 1917C 1918 DOUBLE PRECISION, INTENT(INOUT) :: RHO2(*) 1919 DOUBLE PRECISION, INTENT(IN) :: C1AM(*), FCKMO(*) 1920 INTEGER, INTENT(IN) :: ISYMV, ISYMF 1921#include "priunit.h" 1922#include "ccorb.h" 1923#include "ccsdsym.h" 1924#include "cclr.h" 1925 INTEGER :: ISYMAI, ISYMBJ, ISYRES 1926 INTEGER :: NAI, NBJ 1927 INTEGER :: KP, KM, KRHOP, KRHOM 1928 INTEGER :: KOFFP, KOFFM, KOFF3, NTOTAI 1929 DOUBLE PRECISION :: TEMP, TEMP1, TEMP2 1930C 1931 CALL QENTER(myname) 1932C 1933 ISYRES = MULD2H(ISYMV,ISYMF) 1934 KRHOP = 1 1935 KRHOM = NT2AM(ISYRES) + 1 1936C 1937C-------------------------------------------------------------- 1938C Calculate the contraction of C1AM with FCKMO. 1939C-------------------------------------------------------------- 1940C 1941 ISYMAI = ISYMV 1942 ISYMBJ = ISYMF 1943 1944 IF (ISYMV.EQ.ISYMF) THEN 1945 ISYMAI = ISYMV 1946 ISYMBJ = ISYMAI 1947 DO NBJ = 1, NT1AM(ISYMBJ) 1948 KOFFP = KRHOP + IT2AM(ISYMAI,ISYMBJ) + NBJ*(NBJ-1)/2 1949 KOFFM = KRHOM + IT2AM(ISYMAI,ISYMBJ) + NBJ*(NBJ-1)/2 1950 DO NAI = 1, NBJ-1 1951 KP = KOFFP + (NAI-1) 1952 KM = KOFFM + (NAI-1) 1953 TEMP1 = C1AM(NAI)*FCKMO(NBJ) 1954 TEMP2 = FCKMO(NAI)*C1AM(NBJ) 1955 RHO2(KP) = RHO2(KP) + TEMP1 + TEMP2 1956 RHO2(KM) = RHO2(KM) + TEMP1 - TEMP2 1957 END DO 1958 RHO2(KOFFP+NBJ-1) = 0.0D0 1959 RHO2(KOFFM+NBJ-1) = 0.0D0 1960 END DO 1961 ELSE IF (ISYMV.LT.ISYMF) THEN ! only first term is non-zero 1962 ISYMAI = ISYMV 1963 ISYMBJ = ISYMF 1964 DO NBJ = 1, NT1AM(ISYMBJ) 1965 KOFFP = KRHOP + IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1) 1966 KOFFM = KRHOM + IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1) 1967 DO NAI = 1, NT1AM(ISYMAI) 1968 KP = KOFFP + (NAI-1) 1969 KM = KOFFM + (NAI-1) 1970 TEMP = C1AM(NAI)*FCKMO(NBJ) 1971 RHO2(KP) = RHO2(KP) + TEMP 1972 RHO2(KM) = RHO2(KM) + TEMP 1973 END DO 1974 END DO 1975 ELSE ! Only the second term is stored 1976 ISYMAI = ISYMF 1977 ISYMBJ = ISYMV 1978 DO NBJ = 1, NT1AM(ISYMBJ) 1979 KOFFP = KRHOP + IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1) 1980 KOFFM = KRHOM + IT2AM(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ-1) 1981 DO NAI = 1, NT1AM(ISYMAI) 1982 KP = KOFFP + (NAI-1) 1983 KM = KOFFM + (NAI-1) 1984 TEMP = C1AM(NBJ)*FCKMO(NAI) 1985 RHO2(KP) = RHO2(KP) + TEMP 1986 RHO2(KM) = RHO2(KM) - TEMP 1987 END DO 1988 END DO 1989 END IF 1990C 1991C All done 1992 CALL QEXIT(myname) 1993C 1994 RETURN 1995 END 1996