1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck cc_yi */ 20 SUBROUTINE CC_YI(YMAT,CLTR2,ISYMLV,T2AM,ISYMRV,WORK,LWORK) 21C 22C Written by Asger Halkier 16/10 - 1995 23C 24C Version: 1.0 25C 26C Purpose: To calculate the Y-intermediat entering the E'-term 27C of both the 2.1-block and 2.2-block of Jacobian, 28C and in the F-matrix and etaC. 29C 30C It is assumed that CLTR2 is squared up, and that T2AM is not. 31C 32C Ove Christiansen 20-6-1996: 33C 34C General symmetry of both CLTR2 and T2AM for F-mat. and etaC. 35C 36#include "implicit.h" 37 LOGICAL LOCDBG 38 PARAMETER (LOCDBG = .FALSE.) 39 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 40 DIMENSION YMAT(*),CLTR2(*),T2AM(*),WORK(LWORK) 41#include "priunit.h" 42#include "ccorb.h" 43#include "ccsdsym.h" 44#include "cclr.h" 45C 46 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 47C 48 ISYRES = MULD2H(ISYMLV,ISYMRV) 49C 50C----------------------------- 51C Initialize result array. 52C----------------------------- 53C 54 CALL DZERO(YMAT,NMATAB(ISYRES)) 55C 56 DO 100 ISYMCM = 1,NSYM 57C 58 ISYMFN = MULD2H(ISYMCM,ISYMLV) 59 ISYMEN = MULD2H(ISYMCM,ISYMRV) 60C 61 DO 110 ISYMN = 1,NSYM 62C 63 ISYME = MULD2H(ISYMN,ISYMEN) 64 ISYMF = MULD2H(ISYMN,ISYMFN) 65C 66C------------------------------------- 67C Work space allocation one. 68C------------------------------------- 69C 70 KT2SM = 1 71 KC2SM = KT2SM + NT1AM(ISYMCM)*NVIR(ISYME) 72 KEND1 = KC2SM + NT1AM(ISYMCM)*NVIR(ISYMF) 73 LWRK1 = LWORK - KEND1 74C 75 IF (LWRK1 .LT. 0) THEN 76 CALL QUIT('Insufficient work space in CC_YI') 77 ENDIF 78C 79 CALL DZERO(WORK,KEND1) 80C 81 DO 120 N = 1,NRHF(ISYMN) 82C 83C------------------------------------------------ 84C Copy submatrix out of packed T2AM. 85C------------------------------------------------ 86C 87 DO 130 E = 1,NVIR(ISYME) 88C 89 NEN = IT1AM(ISYME,ISYMN) + NVIR(ISYME)*(N - 1) + E 90C 91 IF (ISYMCM .EQ. ISYMEN) THEN 92C 93 DO 140 NCM = 1,NT1AM(ISYMCM) 94C 95 NCMEN = IT2AM(ISYMCM,ISYMEN) + INDEX(NCM,NEN) 96 NECM = KT2SM + NVIR(ISYME)*(NCM - 1) + E - 1 97C 98 WORK(NECM) = T2AM(NCMEN) 99C 100 140 CONTINUE 101C 102 ELSE IF (ISYMCM .LT. ISYMEN) THEN 103C 104 DO 150 NCM = 1,NT1AM(ISYMCM) 105C 106 NCMEN = IT2AM(ISYMCM,ISYMEN) 107 & + NT1AM(ISYMCM)*(NEN - 1) + NCM 108 NECM = KT2SM + NVIR(ISYME)*(NCM - 1) + E - 1 109C 110 WORK(NECM) = T2AM(NCMEN) 111C 112 150 CONTINUE 113C 114 ELSE IF (ISYMCM .GT. ISYMEN) THEN 115C 116 DO 160 NCM = 1,NT1AM(ISYMCM) 117C 118 NCMEN = IT2AM(ISYMEN,ISYMCM) 119 & + NT1AM(ISYMEN)*(NCM - 1) + NEN 120 NECM = KT2SM + NVIR(ISYME)*(NCM - 1) + E - 1 121C 122 WORK(NECM) = T2AM(NCMEN) 123C 124 160 CONTINUE 125C 126 ENDIF 127C 128 130 CONTINUE 129C 130C----------------------------------------- 131C Copy submatrix out of CLTR2. 132C----------------------------------------- 133C 134 DO 170 F = 1,NVIR(ISYMF) 135C 136 NFN = IT1AM(ISYMF,ISYMN) + NVIR(ISYMF)*(N - 1) + F 137C 138 DO 180 NCM = 1,NT1AM(ISYMCM) 139C 140 NCMFN = IT2SQ(ISYMCM,ISYMFN) 141 & + NT1AM(ISYMCM)*(NFN - 1) + NCM 142 NCMF = KC2SM + NT1AM(ISYMCM)*(F - 1) + NCM - 1 143C 144 WORK(NCMF) = CLTR2(NCMFN) 145C 146 180 CONTINUE 147 170 CONTINUE 148C 149C------------------------------------------- 150C Contract the two submatrices. 151C------------------------------------------- 152C 153 KOFF1 = IMATAB(ISYME,ISYMF) + 1 154C 155 NTOTE = MAX(NVIR(ISYME),1) 156 NTOTCM = MAX(NT1AM(ISYMCM),1) 157C 158 CALL DGEMM('N','N',NVIR(ISYME),NVIR(ISYMF),NT1AM(ISYMCM), 159 & ONE,WORK(KT2SM),NTOTE,WORK(KC2SM),NTOTCM, 160 & ONE,YMAT(KOFF1),NTOTE) 161C 162 120 CONTINUE 163 110 CONTINUE 164 100 CONTINUE 165C 166 IF (LOCDBG) THEN 167 WRITE(LUPRI,*) '-----CC_YI: Norm of Y-intermediate:', 168 & DDOT(NMATAB(ISYRES),YMAT,1,YMAT,1) 169 END IF 170C 171 RETURN 172 END 173C /* Deck cc_xi */ 174 SUBROUTINE CC_XI(XMAT,CLTR2,ISYMLV,T2AM,ISYMRV,WORK,LWORK) 175C 176C Written by Asger Halkier 16/10 - 1995 177C 178C Version: 1.0 179C 180C Purpose: To calculate the X-intermediat entering the E'-term 181C of both the 2.1-block and 2.2-block of Jacobian, 182C and in the F-matrix and etaC. 183C 184C Ove Christiansen 20-6-1996: 185C 186C General symmetry of both CLTR2 and T2AM for F-mat. and etaC. 187C 188C 189#include "implicit.h" 190 LOGICAL LOCDBG 191 PARAMETER (LOCDBG=.FALSE.) 192 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 193 DIMENSION XMAT(*),CLTR2(*),T2AM(*),WORK(LWORK) 194#include "priunit.h" 195#include "ccorb.h" 196#include "ccsdsym.h" 197#include "cclr.h" 198C 199 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 200C 201C----------------------------- 202C Initialize result array. 203C----------------------------- 204C 205 ISYRES = MULD2H(ISYMLV,ISYMRV) 206C 207 CALL DZERO(XMAT,NMATIJ(ISYRES)) 208C 209 DO 100 ISYMEM = 1,NSYM 210C 211 ISYMFL = MULD2H(ISYMEM,ISYMRV) 212 ISYMFJ = MULD2H(ISYMEM,ISYMLV) 213C 214 DO 110 ISYMF = 1,NSYM 215C 216 ISYML = MULD2H(ISYMF,ISYMFL) 217 ISYMJ = MULD2H(ISYMF,ISYMFJ) 218C 219C------------------------------------- 220C Work space allocation one. 221C------------------------------------- 222C 223 KT2SM = 1 224 KC2SM = KT2SM + NT1AM(ISYMEM)*NRHF(ISYML) 225 KEND1 = KC2SM + NT1AM(ISYMEM)*NRHF(ISYMJ) 226 LWRK1 = LWORK - KEND1 227C 228 IF (LWRK1 .LT. 0) THEN 229 CALL QUIT('Insufficient work space in CC_XI') 230 ENDIF 231C 232 CALL DZERO(WORK,KEND1) 233C 234 DO 120 F = 1,NVIR(ISYMF) 235C 236C-------------------------------------- 237C Copy submatrix of T2AM out. 238C-------------------------------------- 239C 240 DO 130 L = 1,NRHF(ISYML) 241C 242 NFL = IT1AM(ISYMF,ISYML) + NVIR(ISYMF)*(L - 1) + F 243C 244 IF (ISYMEM .EQ. ISYMFL) THEN 245C 246 DO 140 NEM = 1,NT1AM(ISYMEM) 247C 248 NEMFL = IT2AM(ISYMEM,ISYMFL) + INDEX(NEM,NFL) 249 NLEM = KT2SM + NRHF(ISYML)*(NEM - 1) + L - 1 250C 251 WORK(NLEM) = T2AM(NEMFL) 252C 253 140 CONTINUE 254C 255 ELSE IF (ISYMEM .LT. ISYMFL) THEN 256C 257 DO 150 NEM = 1,NT1AM(ISYMEM) 258C 259 NEMFL = IT2AM(ISYMEM,ISYMFL) 260 & + NT1AM(ISYMEM)*(NFL - 1) + NEM 261 NLEM = KT2SM + NRHF(ISYML)*(NEM - 1) + L - 1 262C 263 WORK(NLEM) = T2AM(NEMFL) 264C 265 150 CONTINUE 266C 267 ELSE IF (ISYMEM .GT. ISYMFL) THEN 268C 269 DO 160 NEM = 1,NT1AM(ISYMEM) 270C 271 NEMFL = IT2AM(ISYMFL,ISYMEM) 272 & + NT1AM(ISYMFL)*(NEM - 1) + NFL 273 NLEM = KT2SM + NRHF(ISYML)*(NEM - 1) + L - 1 274C 275 WORK(NLEM) = T2AM(NEMFL) 276C 277 160 CONTINUE 278C 279 ENDIF 280C 281 130 CONTINUE 282C 283C----------------------------------------- 284C Copy submatrix of CLTR2 out. 285C----------------------------------------- 286C 287 DO 170 J = 1,NRHF(ISYMJ) 288C 289 NFJ = IT1AM(ISYMF,ISYMJ) + NVIR(ISYMF)*(J - 1) + F 290C 291 DO 180 NEM = 1,NT1AM(ISYMEM) 292C 293 NEMFJ = IT2SQ(ISYMEM,ISYMFJ) 294 & + NT1AM(ISYMEM)*(NFJ - 1) + NEM 295 NEMJ = KC2SM + NT1AM(ISYMEM)*(J - 1) + NEM - 1 296C 297 WORK(NEMJ) = CLTR2(NEMFJ) 298C 299 180 CONTINUE 300 170 CONTINUE 301C 302C------------------------------------------- 303C Contract the two submatrices. 304C------------------------------------------- 305C 306 KOFF1 = IMATIJ(ISYML,ISYMJ) + 1 307C 308 NTOTL = MAX(NRHF(ISYML),1) 309 NTOTEM = MAX(NT1AM(ISYMEM),1) 310C 311 CALL DGEMM('N','N',NRHF(ISYML),NRHF(ISYMJ),NT1AM(ISYMEM), 312 & ONE,WORK(KT2SM),NTOTL,WORK(KC2SM),NTOTEM,ONE, 313 & XMAT(KOFF1),NTOTL) 314C 315 120 CONTINUE 316 110 CONTINUE 317 100 CONTINUE 318C 319 IF (LOCDBG) THEN 320 WRITE(LUPRI,*) '-----CC_XI: Norm of X-intermediate:', 321 & DDOT(NMATIJ(ISYRES),XMAT,1,XMAT,1) 322 END IF 323C 324 RETURN 325 END 326C /* Deck cc_21efm */ 327 SUBROUTINE CC_21EFM(RHO1,FCKMO,ISYFCK,XMAT,YMAT,ISYMIM,WORK,LWORK) 328C 329C Written by Asger Halkier 28/7 - 1995. 330C 331C Version: 1.0 332C 333C Purpose: To calculate the Fock contributions to the E' term 334C of the 1.2-block. 335C 336C Ove Christiansen 20-6-1996 337C General symmetry of fckmo(ISYFCK) and X and Y 338C intermediates ISYMIM 339C 340#include "implicit.h" 341 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 342 DIMENSION RHO1(*),FCKMO(*),XMAT(*),YMAT(*),WORK(LWORK) 343#include "priunit.h" 344#include "ccorb.h" 345#include "ccsdsym.h" 346#include "cclr.h" 347C 348 IF (LWORK .LT. NT1AM(ISYFCK)) THEN 349 CALL QUIT('Insufficient work-space-area in CC_12EMF') 350 ENDIF 351C 352 ISYRES = MULD2H(ISYFCK,ISYMIM) 353C 354C----------------------------------- 355C Extract the fock matrix F(ie). 356C----------------------------------- 357C Note that F(kc) is stored in order c,k 358 DO 50 ISYMC = 1,NSYM 359C 360 ISYMK = MULD2H(ISYMC,ISYFCK) 361C 362 DO 60 K = 1,NRHF(ISYMK) 363C 364 DO 70 C = 1,NVIR(ISYMC) 365C 366 KOFF1 = IFCVIR(ISYMK,ISYMC) + NORB(ISYMK)*(C - 1) + K 367 KOFF2 = IT1AM(ISYMC,ISYMK) + NVIR(ISYMC)*(K - 1) + C 368C 369 WORK(KOFF2) = FCKMO(KOFF1) 370C 371 70 CONTINUE 372 60 CONTINUE 373 50 CONTINUE 374C 375C----------------------------------------------- 376C Calculate contribution from the XMAT-part. 377C----------------------------------------------- 378C 379 DO 80 ISYMA = 1,NSYM 380C 381 ISYML = MULD2H(ISYMA,ISYFCK) 382 ISYMI = MULD2H(ISYMA,ISYRES) 383C 384 KOFF1 = IT1AM(ISYMA,ISYML) + 1 385 KOFF2 = IMATIJ(ISYML,ISYMI) + 1 386 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 387C 388 NTOTL = MAX(NRHF(ISYML),1) 389 NTOTA = MAX(NVIR(ISYMA),1) 390C 391 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYML),-ONE, 392 & WORK(KOFF1),NTOTA,XMAT(KOFF2),NTOTL,ONE, 393 & RHO1(KOFF3),NTOTA) 394C 395 80 CONTINUE 396C 397C----------------------------------------------- 398C Calculate contribution from the YMAT-part. 399C----------------------------------------------- 400C 401 DO 90 ISYMA = 1,NSYM 402C 403 ISYMI = MULD2H(ISYMA,ISYRES) 404 ISYME = MULD2H(ISYMI,ISYFCK) 405C 406 KOFF1 = IMATAB(ISYME,ISYMA) + 1 407 KOFF2 = IT1AM(ISYME,ISYMI) + 1 408 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 409C 410 NTOTE = MAX(NVIR(ISYME),1) 411 NTOTA = MAX(NVIR(ISYMA),1) 412C 413 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYME),-ONE, 414 & YMAT(KOFF1),NTOTE,WORK(KOFF2),NTOTE,ONE, 415 & RHO1(KOFF3),NTOTA) 416C 417 90 CONTINUE 418C 419 RETURN 420 END 421C /* Deck cc_22ec */ 422 SUBROUTINE CC_22EC(RHO2,T2AM,XMAT,YMAT,ISYMXY,WORK,LWORK) 423C 424C Written by Asger Halkier 21/11 - 1995 425C 426C Version: 1.0 427C 428C Purpose: To calculate the coulomb part of the 22E-prime 429C contribution to RHO2. 430C 431C Ove Christiansen 17-9-1996: 432C Gen. sym. input of XMAT and YMAT: ISYMXY 433C (ia|jb) quantity is total symmetric and packed 434C as ai,bj. 435C 436C 437#include "implicit.h" 438 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 439 DIMENSION RHO2(*), T2AM(*), XMAT(*), YMAT(*), WORK(LWORK) 440#include "priunit.h" 441#include "ccorb.h" 442#include "ccsdsym.h" 443#include "cclr.h" 444C 445 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 446C 447 ISYRES = MULD2H(ISYMOP,ISYMXY) 448 ISYMAT = ISYRES 449C 450 DO 100 ISYMJ = 1,NSYM 451C 452 IF (NRHF(ISYMJ) .EQ. 0) GOTO 100 453C 454 DO 110 J = 1,NRHF(ISYMJ) 455C 456 DO 120 ISYMB = 1,NSYM 457C 458 ISYMBJ = MULD2H(ISYMB,ISYMJ) 459 ISYMAI = MULD2H(ISYMBJ,ISYRES) 460 ISYMAN = MULD2H(ISYMBJ,ISYMOP) 461 ISYMFI = ISYMAN 462C 463 IF (NVIR(ISYMB) .EQ. 0) GOTO 120 464C 465C------------------------------------ 466C Work space allocation. 467C------------------------------------ 468C 469 KSCR1 = 1 470 KSCR2 = KSCR1 + NT1AM(ISYMAI) 471 KEND1 = KSCR2 + NT1AM(ISYMAN) 472 LWRK1 = LWORK - KEND1 473C 474 IF (LWRK1 .LT. 0) THEN 475 CALL QUIT('Insufficient work space in CC_22EC') 476 ENDIF 477C 478 DO 130 B = 1,NVIR(ISYMB) 479C 480 CALL DZERO(WORK(KSCR1),NT1AM(ISYMAI)) 481C 482 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 483C 484C------------------------------------------------- 485C Copy submatrix out of integrals. 486C------------------------------------------------- 487C 488 DO 140 NAN = 1,NT1AM(ISYMAN) 489C 490 NANBJ = IT2AM(ISYMAN,ISYMBJ) + INDEX(NAN,NBJ) 491C 492 WORK(KSCR2 + NAN - 1) = T2AM(NANBJ) 493C 494 140 CONTINUE 495C 496C----------------------------------------------------------------- 497C Contraction of integrals with X- and Y-matrices. 498C----------------------------------------------------------------- 499C 500 DO 150 ISYMAX = 1,NSYM 501C 502 ISYMNX = MULD2H(ISYMAX,ISYMAN) 503 ISYMIX = MULD2H(ISYMNX,ISYMAT) 504 ISYMFY = ISYMAX 505 ISYMIY = MULD2H(ISYMFY,ISYMFI) 506 ISYMAY = MULD2H(ISYMFY,ISYMAT) 507C 508 KOFF1 = KSCR2 + IT1AM(ISYMAX,ISYMNX) 509 KOFF2 = IMATIJ(ISYMNX,ISYMIX) + 1 510 KOFF3 = KSCR1 + IT1AM(ISYMAX,ISYMIX) 511C 512 NTOTA = MAX(NVIR(ISYMAX),1) 513 NTOTN = MAX(NRHF(ISYMNX),1) 514C 515 CALL DGEMM('N','N',NVIR(ISYMAX),NRHF(ISYMIX), 516 & NRHF(ISYMNX),ONE,WORK(KOFF1),NTOTA, 517 & XMAT(KOFF2),NTOTN,ONE,WORK(KOFF3), 518 & NTOTA) 519C 520 KOFF4 = IMATAB(ISYMFY,ISYMAY) + 1 521 KOFF5 = KSCR2 + IT1AM(ISYMFY,ISYMIY) 522 KOFF6 = KSCR1 + IT1AM(ISYMAY,ISYMIY) 523C 524 NTOTF = MAX(NVIR(ISYMFY),1) 525 NTOTA = MAX(NVIR(ISYMAY),1) 526C 527 CALL DGEMM('T','N',NVIR(ISYMAY),NRHF(ISYMIY), 528 & NVIR(ISYMFY),ONE,YMAT(KOFF4),NTOTF, 529 & WORK(KOFF5),NTOTF,ONE,WORK(KOFF6), 530 & NTOTA) 531C 532 150 CONTINUE 533C 534C----------------------------------------------- 535C Storage in RHO2 result vector. 536C----------------------------------------------- 537C 538 IF (ISYMAI .EQ. ISYMBJ) THEN 539C 540 WORK(KSCR1 + NBJ - 1) = TWO*WORK(KSCR1 + NBJ - 1) 541C 542 DO 160 NAI = 1,NT1AM(ISYMAI) 543C 544 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 545 & + INDEX(NAI,NBJ) 546C 547 RHO2(NAIBJ) = RHO2(NAIBJ) 548 & - TWO*WORK(KSCR1 + NAI - 1) 549C 550 160 CONTINUE 551C 552 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 553C 554 KOFF7 = IT2AM(ISYMAI,ISYMBJ) 555 & + NT1AM(ISYMAI)*(NBJ - 1) + 1 556C 557 CALL DAXPY(NT1AM(ISYMAI),-TWO,WORK(KSCR1),1, 558 & RHO2(KOFF7),1) 559C 560 ELSE IF (ISYMAI .GT. ISYMBJ) THEN 561C 562 KOFF8 = IT2AM(ISYMBJ,ISYMAI) + NBJ 563C 564 CALL DAXPY(NT1AM(ISYMAI),-TWO,WORK(KSCR1),1, 565 & RHO2(KOFF8),NT1AM(ISYMBJ)) 566C 567 ENDIF 568C 569 130 CONTINUE 570 120 CONTINUE 571 110 CONTINUE 572 100 CONTINUE 573C 574 RETURN 575 END 576C /* Deck cc_22ee */ 577 SUBROUTINE CC_22EE(RHO2,T2AM,XMAT,YMAT,ISYMXY,WORK,LWORK) 578C 579C Written by Asger Halkier 21/11 - 1995 580C 581C Version: 1.0 582C 583C Purpose: To calculate the exchange part of the 22E-prime 584C contribution to RHO2. 585C 586C Ove Christiansen 17-9-1996: 587C Gen. sym. input of XMAT and YMAT: ISYMXY 588C (ia|jb) quantity is total symmetric and packed 589C as ai,bj. 590C 591C 592#include "implicit.h" 593 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 594 DIMENSION RHO2(*), T2AM(*), XMAT(*), YMAT(*), WORK(LWORK) 595#include "priunit.h" 596#include "ccorb.h" 597#include "ccsdsym.h" 598#include "cclr.h" 599C 600 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 601C 602 ISYRES = MULD2H(ISYMOP,ISYMXY) 603 ISYMAT = ISYRES 604C 605 DO 100 ISYMI = 1,NSYM 606C 607 IF (NRHF(ISYMI) .EQ. 0) GOTO 100 608C 609 DO 110 I = 1,NRHF(ISYMI) 610C 611 DO 120 ISYMB = 1,NSYM 612C 613 ISYMBI = MULD2H(ISYMB,ISYMI) 614 ISYMAJ = MULD2H(ISYMBI,ISYRES) 615 ISYMAN = MULD2H(ISYMBI,ISYMOP) 616 ISYMFJ = ISYMAN 617C 618 IF (NVIR(ISYMB) .EQ. 0) GOTO 120 619C 620C------------------------------------ 621C Work space allocation. 622C------------------------------------ 623C 624 KSCR1 = 1 625 KSCR2 = KSCR1 + NT1AM(ISYMAJ) 626 KEND1 = KSCR2 + NT1AM(ISYMAN) 627 LWRK1 = LWORK - KEND1 628C 629 IF (LWRK1 .LT. 0) THEN 630 CALL QUIT('Insufficient work space in CC_22EE') 631 ENDIF 632C 633 DO 130 B = 1,NVIR(ISYMB) 634C 635 CALL DZERO(WORK(KSCR1),NT1AM(ISYMAJ)) 636C 637 NBI = IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) + B 638C 639C------------------------------------------------- 640C Copy submatrix out of integrals. 641C------------------------------------------------- 642C 643 DO 140 NAN = 1,NT1AM(ISYMAN) 644C 645 NANBI = IT2AM(ISYMAN,ISYMBI) + INDEX(NAN,NBI) 646C 647 WORK(KSCR2 + NAN - 1) = T2AM(NANBI) 648C 649 140 CONTINUE 650C 651C----------------------------------------------------------------- 652C Contraction of integrals with X- and Y-matrices. 653C----------------------------------------------------------------- 654C 655 DO 150 ISYMAX = 1,NSYM 656C 657 ISYMNX = MULD2H(ISYMAX,ISYMAN) 658 ISYMJX = MULD2H(ISYMNX,ISYMAT) 659 ISYMFY = ISYMAX 660 ISYMJY = MULD2H(ISYMFY,ISYMFJ) 661 ISYMAY = MULD2H(ISYMFY,ISYMAT) 662C 663 KOFF1 = KSCR2 + IT1AM(ISYMAX,ISYMNX) 664 KOFF2 = IMATIJ(ISYMNX,ISYMJX) + 1 665 KOFF3 = KSCR1 + IT1AM(ISYMAX,ISYMJX) 666C 667 NTOTA = MAX(NVIR(ISYMAX),1) 668 NTOTN = MAX(NRHF(ISYMNX),1) 669C 670 CALL DGEMM('N','N',NVIR(ISYMAX),NRHF(ISYMJX), 671 & NRHF(ISYMNX),ONE,WORK(KOFF1),NTOTA, 672 & XMAT(KOFF2),NTOTN,ONE,WORK(KOFF3), 673 & NTOTA) 674C 675 KOFF4 = IMATAB(ISYMFY,ISYMAY) + 1 676 KOFF5 = KSCR2 + IT1AM(ISYMFY,ISYMJY) 677 KOFF6 = KSCR1 + IT1AM(ISYMAY,ISYMJY) 678C 679 NTOTF = MAX(NVIR(ISYMFY),1) 680 NTOTA = MAX(NVIR(ISYMAY),1) 681C 682 CALL DGEMM('T','N',NVIR(ISYMAY),NRHF(ISYMJY), 683 & NVIR(ISYMFY),ONE,YMAT(KOFF4),NTOTF, 684 & WORK(KOFF5),NTOTF,ONE,WORK(KOFF6), 685 & NTOTA) 686C 687 150 CONTINUE 688C 689C----------------------------------------------- 690C Storage in RHO2 result vector. 691C----------------------------------------------- 692C 693 DO 160 ISYMJ = 1,NSYM 694C 695 ISYMBJ = MULD2H(ISYMB,ISYMJ) 696 ISYMAI = MULD2H(ISYMBJ,ISYRES) 697 ISYMA = MULD2H(ISYMJ,ISYMAJ) 698C 699 DO 170 J = 1,NRHF(ISYMJ) 700C 701 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) 702 & + B 703C 704 IF (ISYMAI .EQ. ISYMBJ) THEN 705C 706 DO 180 A = 1,NVIR(ISYMA) 707C 708 NAJ = IT1AM(ISYMA,ISYMJ) 709 & + NVIR(ISYMA)*(J - 1) + A 710 NAI = IT1AM(ISYMA,ISYMI) 711 & + NVIR(ISYMA)*(I - 1) + A 712 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 713 & + INDEX(NAI,NBJ) 714C 715 RHO2(NAIBJ) = RHO2(NAIBJ) 716 & + WORK(KSCR1 + NAJ - 1) 717C 718 IF (NAI .EQ. NBJ) THEN 719C 720 RHO2(NAIBJ) = RHO2(NAIBJ) 721 & + WORK(KSCR1 + NAJ - 1) 722C 723 ENDIF 724C 725 180 CONTINUE 726C 727 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 728C 729 KOFF7 = KSCR1 + IT1AM(ISYMA,ISYMJ) 730 & + NVIR(ISYMA)*(J - 1) 731 KOFF8 = IT2AM(ISYMAI,ISYMBJ) 732 & + NT1AM(ISYMAI)*(NBJ - 1) 733 & + IT1AM(ISYMA,ISYMI) 734 & + NVIR(ISYMA)*(I - 1) + 1 735C 736 CALL DAXPY(NVIR(ISYMA),ONE,WORK(KOFF7),1, 737 & RHO2(KOFF8),1) 738C 739 ELSE IF (ISYMAI .GT. ISYMBJ) THEN 740C 741 DO 190 A = 1,NVIR(ISYMA) 742C 743 NAJ = IT1AM(ISYMA,ISYMJ) 744 & + NVIR(ISYMA)*(J - 1) + A 745 NAI = IT1AM(ISYMA,ISYMI) 746 & + NVIR(ISYMA)*(I - 1) + A 747 NAIBJ = IT2AM(ISYMBJ,ISYMAI) 748 & + NT1AM(ISYMBJ)*(NAI - 1) + NBJ 749C 750 RHO2(NAIBJ) = RHO2(NAIBJ) 751 & + WORK(KSCR1 + NAJ - 1) 752C 753 190 CONTINUE 754C 755 ENDIF 756C 757 170 CONTINUE 758 160 CONTINUE 759 130 CONTINUE 760 120 CONTINUE 761 110 CONTINUE 762 100 CONTINUE 763C 764 RETURN 765 END 766C /* Deck cc_22c */ 767 SUBROUTINE CC_22C(RHO2,CTR2,ISYVEC,XLAMDH,WORK,LWORK, 768 * ISYCIM,LUC,CFIL,IV,IOPT) 769C 770C Written by Asger Halkier 27/11 - 1995 771C 772C Version: 1.0 773C 774C Purpose: To calculate the 22C contribution to RHO2. 775C 776C Ove Christiansen 17-9-1996: 777C 778C modified to account for general 779C non. total symmetric vectors (ISYVEC) and 780C intermediates (ISYCIM). LUC and CFIL is 781C used to control from which file the 782C intermediate is obtained. 783C 784C if iopt = 1 the C intermediate is assumed 785C to be as in energy calc. 786C 787C if iopt ne. 1 we use the intermediate 788C on lud with address given according to 789C transformed vector nr iv (iv is not 1 790C if several vectors are transformed 791C simultaneously.) 792C 793C in ordinary Lambda equations: iv=1,iopt=1 794C 795#include "implicit.h" 796 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 797 DIMENSION RHO2(*), CTR2(*), XLAMDH(*), WORK(LWORK) 798 CHARACTER CFIL*(*) 799#include "priunit.h" 800#include "ccorb.h" 801#include "ccsdsym.h" 802#include "cclr.h" 803#include "maxorb.h" 804#include "ccsdio.h" 805C 806 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 807C 808 ISYRES = MULD2H(ISYVEC,ISYCIM) 809C 810 DO 100 ISYMD = 1,NSYM 811C 812 ISYMB = ISYMD 813 ISYIEM = MULD2H(ISYMD,ISYCIM) 814 ISYAJI = MULD2H(ISYMB,ISYRES) 815C 816C------------------------------ 817C Work space allocation. 818C------------------------------ 819C 820 KPINT = 1 821 KSCR1 = KPINT + NT2BCD(ISYIEM) 822 KEND1 = KSCR1 + NT2BCD(ISYAJI) 823 LWRK1 = LWORK - KEND1 824C 825 IF (LWRK1 .LT. 0) THEN 826 CALL QUIT('Insufficient work space in CC_22C') 827 ENDIF 828C 829 DO 110 IDEL = 1,NBAS(ISYMD) 830C 831C--------------------------------------------------- 832C Read P-intermediate submatrix from disc. 833C--------------------------------------------------- 834C 835 ID = IDEL + IBAS(ISYMD) 836C 837 IF ( IOPT .EQ. 1) THEN 838 IOFF = IT2DEL(ID) + 1 839 ELSE 840 IOFF = IT2DLR(ID,IV) + 1 841 ENDIF 842C 843 LEN = NT2BCD(ISYIEM) 844C 845 IF (LEN .GT. 0) THEN 846 CALL GETWA2(LUC,CFIL,WORK(KPINT),IOFF,LEN) 847 ENDIF 848C 849C-------------------------------------------------------- 850C Contraction of intermediate and trial vector. 851C-------------------------------------------------------- 852C 853 DO 120 ISYMAJ = 1,NSYM 854C 855 ISYMEM = MULD2H(ISYMAJ,ISYVEC) 856 ISYMI = MULD2H(ISYMAJ,ISYAJI) 857C 858 KOFF1 = IT2SQ(ISYMAJ,ISYMEM) + 1 859 KOFF2 = KPINT + IT2BCT(ISYMI,ISYMEM) 860 KOFF3 = KSCR1 + IT2BCD(ISYMAJ,ISYMI) 861C 862 NTOTAJ = MAX(NT1AM(ISYMAJ),1) 863 NTOTI = MAX(NRHF(ISYMI),1) 864C 865 CALL DGEMM('N','T',NT1AM(ISYMAJ),NRHF(ISYMI), 866 & NT1AM(ISYMEM),ONE,CTR2(KOFF1),NTOTAJ, 867 & WORK(KOFF2),NTOTI,ZERO,WORK(KOFF3),NTOTAJ) 868C 869 120 CONTINUE 870C 871C--------------------------------------------- 872C Resort intermediate result-matrix. 873C--------------------------------------------- 874C 875 CALL CCLT_P21I(WORK(KSCR1),ISYAJI,WORK(KEND1),LWRK1, 876 & IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR) 877C 878C----------------------------------------------------------------- 879C Final scaling with XLAMDH-element and storage in RHO2. 880C----------------------------------------------------------------- 881C 882 DO 130 B = 1,NVIR(ISYMB) 883C 884 NDB = ILMVIR(ISYMB) + NBAS(ISYMD)*(B - 1) + IDEL 885C 886 FACT = -HALF*XLAMDH(NDB) 887C 888 DO 140 ISYMJ = 1,NSYM 889C 890 ISYMAI = MULD2H(ISYMJ,ISYAJI) 891 ISYMBJ = MULD2H(ISYMAI,ISYRES) 892C 893 DO 150 J = 1,NRHF(ISYMJ) 894C 895 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 896 NJ = KSCR1 + IT2BCD(ISYMAI,ISYMJ) 897 & + NT1AM(ISYMAI)*(J - 1) 898C 899 IF (ISYMAI .EQ. ISYMBJ) THEN 900C 901 DO 160 NAI = 1,NT1AM(ISYMAI) 902C 903 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 904 & + INDEX(NAI,NBJ) 905C 906 RHO2(NAIBJ) = RHO2(NAIBJ) 907 & + FACT*WORK(NJ + NAI - 1) 908C 909 IF (NAI .EQ. NBJ) THEN 910C 911 RHO2(NAIBJ) = RHO2(NAIBJ) 912 & + FACT*WORK(NJ + NAI - 1) 913C 914 ENDIF 915C 916 160 CONTINUE 917C 918 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 919C 920 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 921 & + NT1AM(ISYMAI)*(NBJ - 1) + 1 922C 923 CALL DAXPY(NT1AM(ISYMAI),FACT,WORK(NJ),1, 924 & RHO2(NAIBJ),1) 925C 926 ELSE IF (ISYMAI .GT. ISYMBJ) THEN 927C 928 NAIBJ = IT2AM(ISYMBJ,ISYMAI) + NBJ 929C 930 CALL DAXPY(NT1AM(ISYMAI),FACT,WORK(NJ),1, 931 & RHO2(NAIBJ),NT1AM(ISYMBJ)) 932C 933 ENDIF 934C 935 150 CONTINUE 936 140 CONTINUE 937 130 CONTINUE 938 110 CONTINUE 939 100 CONTINUE 940C 941 RETURN 942 END 943C /* Deck cc_22d */ 944 SUBROUTINE CC_22D(RHO2,CTR2,ISYVEC,XLAMDH,WORK,LWORK, 945 * ISYDIM,LUD,DFIL,IV,IOPT) 946C 947C Written by Asger Halkier 27/11 - 1995 948C 949C Version: 1.0 950C 951C Purpose: To calculate the 22D contribution to RHO2! 952C 953C Ove Christiansen 17-9-1996: 954C Modified to account for general 955C non. total symmetric vectors (ISYVEC) and 956C intermediates (ISYDIM). LUD and DFIL is 957C used to control from which file the 958C intermediate is obtained. 959C 960C if iopt = 1 the D intermediate is assumed 961C to be as in energy calc. 962C 963C if iopt ne. 1 we use the intermediate 964C on lud with address given according to 965C transformed vector nr iv (iv is not 1 966C if several vectors are transformed 967C simultaneously.) 968C 969C in Lambda vector calc: iv=1,iopt=1 970C 971#include "implicit.h" 972 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 973 DIMENSION RHO2(*), CTR2(*), XLAMDH(*), WORK(LWORK) 974 CHARACTER DFIL*(*) 975#include "priunit.h" 976#include "ccorb.h" 977#include "ccsdsym.h" 978#include "cclr.h" 979#include "maxorb.h" 980#include "ccsdio.h" 981C 982 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 983C 984 ISYRES = MULD2H(ISYVEC,ISYDIM) 985C 986 DO 100 ISYMD = 1,NSYM 987C 988 ISYMB = ISYMD 989 ISYJEM = MULD2H(ISYMD,ISYDIM) 990 ISYAIJ = MULD2H(ISYMB,ISYRES) 991C 992C------------------------------ 993C Work space allocation. 994C------------------------------ 995C 996 KPINT = 1 997 KSCRN = KPINT + NT2BCD(ISYJEM) 998 KSCRT = KSCRN + NT2BCD(ISYAIJ) 999 KEND1 = KSCRT + NT2BCD(ISYAIJ) 1000 LWRK1 = LWORK - KEND1 1001C 1002 IF (LWRK1 .LT. 0) THEN 1003 CALL QUIT('Insufficient work space in CC_22D') 1004 ENDIF 1005C 1006 DO 110 IDEL = 1,NBAS(ISYMD) 1007C 1008C--------------------------------------------------- 1009C Read P-intermediate submatrix from disc. 1010C--------------------------------------------------- 1011C 1012 ID = IDEL + IBAS(ISYMD) 1013C 1014 IF (IOPT .EQ. 1) THEN 1015 IOFF = IT2DEL(ID) + 1 1016 ELSE 1017 IOFF = IT2DLR(ID,IV) + 1 1018 ENDIF 1019C 1020 LEN = NT2BCD(ISYJEM) 1021C 1022 IF (LEN .GT. 0) THEN 1023 CALL GETWA2(LUD,DFIL,WORK(KPINT),IOFF,LEN) 1024C 1025 ENDIF 1026C 1027C-------------------------------------------------------- 1028C Contraction of intermediate and trial vector. 1029C-------------------------------------------------------- 1030C 1031 DO 120 ISYMAI = 1,NSYM 1032C 1033 ISYMJ = MULD2H(ISYMAI,ISYAIJ) 1034 ISYMEM = MULD2H(ISYMAI,ISYVEC) 1035C 1036 KOFF1 = IT2SQ(ISYMAI,ISYMEM) + 1 1037 KOFF2 = KPINT + IT2BCT(ISYMJ,ISYMEM) 1038 KOFF3 = KSCRN + IT2BCD(ISYMAI,ISYMJ) 1039C 1040 NTOTAI = MAX(NT1AM(ISYMAI),1) 1041 NTOTJ = MAX(NRHF(ISYMJ),1) 1042C 1043 CALL DGEMM('N','T',NT1AM(ISYMAI),NRHF(ISYMJ), 1044 & NT1AM(ISYMEM),ONE,CTR2(KOFF1),NTOTAI, 1045 & WORK(KOFF2),NTOTJ,ZERO,WORK(KOFF3),NTOTAI) 1046C 1047 120 CONTINUE 1048C 1049C---------------------------------------- 1050C Construct resorted submatrix. 1051C---------------------------------------- 1052C 1053 CALL DCOPY(NT2BCD(ISYAIJ),WORK(KSCRN),1,WORK(KSCRT),1) 1054C 1055 CALL CCLT_P21I(WORK(KSCRT),ISYAIJ,WORK(KEND1),LWRK1, 1056 & IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR) 1057C 1058C----------------------------------------------------------------- 1059C Final scaling with XLAMDH-element and storage in RHO2. 1060C----------------------------------------------------------------- 1061C 1062 DO 130 B = 1,NVIR(ISYMB) 1063C 1064 NDB = ILMVIR(ISYMB) + NBAS(ISYMD)*(B - 1) + IDEL 1065C 1066 FACTN = XLAMDH(NDB) 1067 FACTT = -HALF*XLAMDH(NDB) 1068C 1069 DO 140 ISYMJ = 1,NSYM 1070C 1071 ISYMAI = MULD2H(ISYMJ,ISYAIJ) 1072 ISYMBJ = MULD2H(ISYMAI,ISYRES) 1073C 1074 DO 150 J = 1,NRHF(ISYMJ) 1075C 1076 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 1077 NJN = KSCRN + IT2BCD(ISYMAI,ISYMJ) 1078 & + NT1AM(ISYMAI)*(J - 1) 1079 NJT = KSCRT + IT2BCD(ISYMAI,ISYMJ) 1080 & + NT1AM(ISYMAI)*(J - 1) 1081C 1082 IF (ISYMAI .EQ. ISYMBJ) THEN 1083C 1084 DO 160 NAI = 1,NT1AM(ISYMAI) 1085C 1086 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 1087 & + INDEX(NAI,NBJ) 1088C 1089 RHO2(NAIBJ) = RHO2(NAIBJ) 1090 & + FACTN*WORK(NJN + NAI - 1) 1091 & + FACTT*WORK(NJT + NAI - 1) 1092C 1093 IF (NAI .EQ. NBJ) THEN 1094C 1095 RHO2(NAIBJ) = RHO2(NAIBJ) 1096 & + FACTN*WORK(NJN + NAI - 1) 1097 & + FACTT*WORK(NJT + NAI - 1) 1098C 1099 ENDIF 1100C 1101 160 CONTINUE 1102C 1103 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 1104C 1105 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 1106 & + NT1AM(ISYMAI)*(NBJ - 1) + 1 1107C 1108 CALL DAXPY(NT1AM(ISYMAI),FACTN,WORK(NJN),1, 1109 & RHO2(NAIBJ),1) 1110 CALL DAXPY(NT1AM(ISYMAI),FACTT,WORK(NJT),1, 1111 & RHO2(NAIBJ),1) 1112C 1113 ELSE IF (ISYMAI .GT. ISYMBJ) THEN 1114C 1115 NAIBJ = IT2AM(ISYMBJ,ISYMAI) + NBJ 1116C 1117 CALL DAXPY(NT1AM(ISYMAI),FACTN,WORK(NJN),1, 1118 & RHO2(NAIBJ),NT1AM(ISYMBJ)) 1119 CALL DAXPY(NT1AM(ISYMAI),FACTT,WORK(NJT),1, 1120 & RHO2(NAIBJ),NT1AM(ISYMBJ)) 1121C 1122 ENDIF 1123C 1124 150 CONTINUE 1125 140 CONTINUE 1126 130 CONTINUE 1127 110 CONTINUE 1128 100 CONTINUE 1129C 1130 RETURN 1131 END 1132C /* Deck cc_mi */ 1133 SUBROUTINE CC_MI(XMINT,CTR2,IC2SYM,T2AM,IT2SYM,WORK,LWORK) 1134C 1135C Written by Asger Halkier 28/10 - 1995 1136C 1137C Version: 1.0 1138C 1139C Purpose: To calculate the M-intermediat entering the G-term 1140C of both the 2.1-block. 1141C 1142C It is assumed that CTR2 is squared up, and that T2AM is not! 1143C 1144C Ove Christiansen 18-9-1996: 1145C Non-total symmetric CTR2 and T2AM amplitudes. 1146C Sym: IC2SYM IT2SYM 1147C 1148#include "implicit.h" 1149 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1150 DIMENSION XMINT(*),CTR2(*),T2AM(*),WORK(LWORK) 1151#include "priunit.h" 1152#include "ccorb.h" 1153#include "ccsdsym.h" 1154#include "cclr.h" 1155C 1156 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 1157C 1158 ISYINT = MULD2H(IC2SYM,IT2SYM) 1159C 1160C----------------------------- 1161C Initialize result array. 1162C----------------------------- 1163C 1164 CALL DZERO(XMINT,N3ORHF(ISYINT)) 1165C 1166 DO 100 ISYME = 1,NSYM 1167C 1168 DO 110 E = 1,NVIR(ISYME) 1169C 1170 DO 120 ISYML = 1,NSYM 1171C 1172 ISYMEL = MULD2H(ISYML,ISYME) 1173 ISYMDJ = MULD2H(ISYMEL,IT2SYM) 1174C 1175 DO 130 L = 1,NRHF(ISYML) 1176C 1177 NEL = IT1AM(ISYME,ISYML) + NVIR(ISYME)*(L - 1) + E 1178C 1179 DO 140 ISYMI = 1,NSYM 1180C 1181 ISYMEI = MULD2H(ISYMI,ISYME) 1182 ISYMDK = MULD2H(ISYMEI,IC2SYM) 1183 ISYJKL = MULD2H(ISYMI,ISYINT) 1184 ISYMJK = MULD2H(ISYJKL,ISYML) 1185C 1186 DO 150 ISYMJ = 1,NSYM 1187C 1188 ISYMD = MULD2H(ISYMJ,ISYMDJ) 1189 ISYMK = MULD2H(ISYMJ,ISYMJK) 1190C 1191 LET2SM = NRHF(ISYMJ)*NVIR(ISYMD) 1192C 1193 IF (LWORK .LT. LET2SM) THEN 1194 CALL QUIT('Insufficient work space in CC_MI') 1195 ENDIF 1196C 1197C----------------------------------------------------- 1198C Copy T1-amplitude out of T2AM. 1199C----------------------------------------------------- 1200C 1201 IF (ISYMDJ .EQ. ISYMEL) THEN 1202C 1203 DO 160 D = 1,NVIR(ISYMD) 1204C 1205 DO 170 J = 1,NRHF(ISYMJ) 1206C 1207 NDJ = IT1AM(ISYMD,ISYMJ) 1208 & + NVIR(ISYMD)*(J - 1) + D 1209 NDJEL = IT2AM(ISYMDJ,ISYMEL) 1210 & + INDEX(NDJ,NEL) 1211 NJD = NRHF(ISYMJ)*(D - 1) + J 1212C 1213 WORK(NJD) = T2AM(NDJEL) 1214C 1215 170 CONTINUE 1216 160 CONTINUE 1217C 1218 ELSE IF (ISYMDJ .LT. ISYMEL) THEN 1219C 1220 DO 180 J = 1,NRHF(ISYMJ) 1221C 1222 NDJ = IT1AM(ISYMD,ISYMJ) 1223 & + NVIR(ISYMD)*(J - 1) + 1 1224 NDJEL = IT2AM(ISYMDJ,ISYMEL) 1225 & + NT1AM(ISYMDJ)*(NEL - 1) + NDJ 1226C 1227 CALL DCOPY(NVIR(ISYMD),T2AM(NDJEL),1, 1228 & WORK(J),NRHF(ISYMJ)) 1229C 1230 180 CONTINUE 1231C 1232 ELSE IF (ISYMDJ .GT. ISYMEL) THEN 1233C 1234 DO 190 D = 1,NVIR(ISYMD) 1235C 1236 DO 200 J = 1,NRHF(ISYMJ) 1237C 1238 NDJ = IT1AM(ISYMD,ISYMJ) 1239 & + NVIR(ISYMD)*(J - 1) + D 1240 NELDJ = IT2AM(ISYMEL,ISYMDJ) 1241 & + NT1AM(ISYMEL)*(NDJ - 1) + NEL 1242 NJD = NRHF(ISYMJ)*(D - 1) + J 1243C 1244 WORK(NJD) = T2AM(NELDJ) 1245C 1246 200 CONTINUE 1247 190 CONTINUE 1248C 1249 ENDIF 1250C 1251C---------------------------------------------------------------------- 1252C Contraction of T2AM & CTR2 & storage in result. 1253C---------------------------------------------------------------------- 1254C 1255 DO 210 I = 1,NRHF(ISYMI) 1256C 1257 NEI = IT1AM(ISYME,ISYMI) 1258 & + NVIR(ISYME)*(I - 1) + E 1259 KOFF2 = IT2SQ(ISYMDK,ISYMEI) 1260 & + NT1AM(ISYMDK)*(NEI - 1) 1261 & + IT1AM(ISYMD,ISYMK) + 1 1262 KOFF3 = I3ORHF(ISYJKL,ISYMI) 1263 & + NMAIJK(ISYJKL)*(I - 1) 1264 & + IMAIJK(ISYMJK,ISYML) 1265 & + NMATIJ(ISYMJK)*(L - 1) 1266 & + IMATIJ(ISYMJ,ISYMK) + 1 1267C 1268 NTOTJ = MAX(NRHF(ISYMJ),1) 1269 NTOTD = MAX(NVIR(ISYMD),1) 1270C 1271 CALL DGEMM('N','N',NRHF(ISYMJ),NRHF(ISYMK), 1272 & NVIR(ISYMD),ONE,WORK,NTOTJ, 1273 & CTR2(KOFF2),NTOTD,ONE, 1274 & XMINT(KOFF3),NTOTJ) 1275C 1276 210 CONTINUE 1277 150 CONTINUE 1278 140 CONTINUE 1279 130 CONTINUE 1280 120 CONTINUE 1281 110 CONTINUE 1282 100 CONTINUE 1283C 1284 RETURN 1285 END 1286C /* Deck cc_22am */ 1287 SUBROUTINE CC_22AM(RHO2,XIAJB,XMINT,ISYMIM,WORK,LWORK) 1288C 1289C Written by Asger Halkier 16/11 - 1995 1290C 1291C Version: 1.0 1292C 1293C Purpose: To calculate the 22A-prime contribution to RHO2! 1294C 1295C Ove Christiansen 18-9-1996: 1296C ISYMIM is symmetry of intermediate. 1297C 1298#include "implicit.h" 1299 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1300 DIMENSION RHO2(*), XIAJB(*), XMINT(*), WORK(LWORK) 1301#include "priunit.h" 1302#include "ccorb.h" 1303#include "ccsdsym.h" 1304#include "cclr.h" 1305C 1306 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 1307C 1308 ISYRES = ISYMIM 1309C 1310 DO 100 ISYMJ = 1,NSYM 1311C 1312 ISYMIN = MULD2H(ISYMJ,ISYRES) 1313C 1314 IF (NRHF(ISYMJ) .EQ. 0) GOTO 100 1315C 1316 DO 110 J = 1,NRHF(ISYMJ) 1317C 1318 DO 120 ISYMB = 1,NSYM 1319C 1320 ISYMBJ = MULD2H(ISYMB,ISYMJ) 1321 ISYMAI = MULD2H(ISYMBJ,ISYRES) 1322C 1323 IF (ISYMAI .GT. ISYMBJ) GOTO 120 1324 IF (NVIR(ISYMB) .EQ. 0) GOTO 120 1325C 1326C---------------------------------------- 1327C Work space allocation one. 1328C---------------------------------------- 1329C 1330 KSCR1 = 1 1331 KEND1 = KSCR1 + NT1AM(ISYMAI) 1332 LWRK1 = LWORK - KEND1 1333C 1334 IF (LWRK1 .LT. 0) THEN 1335 CALL QUIT('Insufficient work space in CC_22AM') 1336 ENDIF 1337C 1338 DO 130 B = 1,NVIR(ISYMB) 1339C 1340 CALL DZERO(WORK(KSCR1),NT1AM(ISYMAI)) 1341C 1342 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 1343C 1344 DO 140 ISYMN = 1,NSYM 1345C 1346 ISYMBN = MULD2H(ISYMN,ISYMB) 1347 ISYMAM = MULD2H(ISYMBN,ISYMOP) 1348 ISYMMI = MULD2H(ISYMN,ISYMIN) 1349C 1350C---------------------------------------------- 1351C Work space allocation two. 1352C---------------------------------------------- 1353C 1354 KSCR2 = KEND1 1355 KEND2 = KSCR2 + NT1AM(ISYMAM) 1356 LWRK2 = LWORK - KEND2 1357C 1358 IF (LWRK2 .LT. 0) THEN 1359 CALL QUIT('Insufficient work space in CC_22AM') 1360 ENDIF 1361C 1362 DO 150 N = 1,NRHF(ISYMN) 1363C 1364 NBN = IT1AM(ISYMB,ISYMN) 1365 & + NVIR(ISYMB)*(N - 1) + B 1366C 1367C------------------------------------------------------- 1368C Copy submatrix out of integrals. 1369C------------------------------------------------------- 1370C 1371 DO 160 NAM = 1,NT1AM(ISYMAM) 1372C 1373 NAMBN = IT2AM(ISYMAM,ISYMBN) 1374 & + INDEX(NAM,NBN) 1375C 1376 WORK(KSCR2 + NAM - 1) = XIAJB(NAMBN) 1377C 1378 160 CONTINUE 1379C 1380C-------------------------------------------------------------------- 1381C Contraction of integrals with M-intermediate. 1382C-------------------------------------------------------------------- 1383C 1384 DO 170 ISYMA = 1,NSYM 1385C 1386 ISYMM = MULD2H(ISYMA,ISYMAM) 1387 ISYMI = MULD2H(ISYMA,ISYMAI) 1388C 1389 KOFF1 = KSCR2 + IT1AM(ISYMA,ISYMM) 1390 KOFF2 = I3ORHF(ISYMIN,ISYMJ) 1391 & + NMAIJK(ISYMIN)*(J - 1) 1392 & + IMAIJK(ISYMMI,ISYMN) 1393 & + NMATIJ(ISYMMI)*(N - 1) 1394 & + IMATIJ(ISYMM,ISYMI) + 1 1395 KOFF3 = KSCR1 + IT1AM(ISYMA,ISYMI) 1396C 1397 NTOTA = MAX(NVIR(ISYMA),1) 1398 NTOTM = MAX(NRHF(ISYMM),1) 1399C 1400 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI), 1401 & NRHF(ISYMM),ONE,WORK(KOFF1), 1402 & NTOTA,XMINT(KOFF2),NTOTM,ONE, 1403 & WORK(KOFF3),NTOTA) 1404C 1405 170 CONTINUE 1406C 1407 150 CONTINUE 1408 140 CONTINUE 1409C 1410C----------------------------------------------- 1411C Storage in RHO2 result vector. 1412C----------------------------------------------- 1413C 1414 IF (ISYMAI .EQ. ISYMBJ) THEN 1415C 1416 DO 180 NAI = 1,NBJ 1417C 1418 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 1419 & + INDEX(NAI,NBJ) 1420C 1421 RHO2(NAIBJ) = RHO2(NAIBJ) 1422 & + WORK(KSCR1 + NAI - 1) 1423C 1424 180 CONTINUE 1425C 1426 ELSE IF (ISYMAI .LT. ISYMBJ) THEN 1427C 1428 KOFF4 = IT2AM(ISYMAI,ISYMBJ) 1429 & + NT1AM(ISYMAI)*(NBJ - 1) + 1 1430C 1431 CALL DAXPY(NT1AM(ISYMAI),ONE,WORK(KSCR1),1, 1432 & RHO2(KOFF4),1) 1433C 1434 ENDIF 1435C 1436 130 CONTINUE 1437 120 CONTINUE 1438 110 CONTINUE 1439 100 CONTINUE 1440C 1441 RETURN 1442 END 1443