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_gam */ 20 SUBROUTINE CC_GAM(DSRHF,GAMMA,XLAMDP,XLAMDH, 21 * XLAMPC,XLAMHC,ISYMPC,SCRM,ISYMCM, 22 * WORK,LWORK,IDEL,ISYMD,IOPT) 23C 24C Written by Henrik Koch 3-Jan-1994 25C Symmetry by Henrik Koch and Alfredo Sanchez. 21-July-1994 26C 27C Purpose: Calculate the gamma intermediate. 28C 29C Ove Christiansen 18-9-1996: 30C General symmetry for F-matrix construction. 31C 32#include "implicit.h" 33 DIMENSION DSRHF(*),GAMMA(*),SCRM(*) 34 DIMENSION WORK(LWORK) 35 DIMENSION XLAMDP(*),XLAMDH(*),XLAMPC(*),XLAMHC(*) 36#include "priunit.h" 37#include "ccorb.h" 38#include "ccsdsym.h" 39C 40C------------------------ 41C Dynamic allocation. 42C------------------------ 43C 44 ISYMJB = MULD2H(ISYMD,ISYMPC) 45 KLAMDA = 1 46 KLAMDC = KLAMDA + NRHF(ISYMD) 47 KEND1 = KLAMDC + NRHF(ISYMJB) 48 LWRK1 = LWORK - KEND1 49C 50 IF (LWRK1 .LT. 0) THEN 51 WRITE(LUPRI,*) 'Need : ',KEND1,'Available : ',LWORK 52 CALL QUIT('Insufficient space in CC_GAM') 53 ENDIF 54C 55C--------------------------------------- 56C Copy XLAMDH vector for given IDEL. 57C--------------------------------------- 58C 59 KOFF1 = ILMRHF(ISYMD) + IDEL - IBAS(ISYMD) 60 CALL DCOPY(NRHF(ISYMD),XLAMDH(KOFF1),NBAS(ISYMD),WORK(KLAMDA),1) 61C 62C-------------------------------- 63C Calculate the contribution. 64C-------------------------------- 65C 66 ISYDIS = MULD2H(ISYMD,ISYMOP) 67C 68 DO 100 ISYML = 1,NSYM 69C 70 ISYMAG = MULD2H(ISYML,ISYDIS) 71 ISYMKI = MULD2H(ISYMAG,ISYMPC) 72C 73C--------------------------- 74C Dynamic allocation. 75C--------------------------- 76C 77 KSCR1 = KEND1 78 KSCR2 = KSCR1 + N2BST(ISYMAG) 79 KSCR3 = KSCR2 + NT1AO(ISYMAG) 80 KSCR4 = KSCR3 + NT1AM(ISYMAG) 81 KSCR5 = KSCR4 + NMATIJ(ISYMAG) 82C 83 IF (IOPT .EQ. 1) THEN 84 KEND2 = KSCR5 85 ELSE 86 KEND2 = KSCR5 + NMATIJ(ISYMKI) 87 ENDIF 88C 89 LWRK2 = LWORK - KEND2 90C 91 IF (LWRK2 .LT. 0) THEN 92 WRITE(LUPRI,*) 'Need : ',KEND2,'Available : ',LWORK 93 CALL QUIT('Insufficient space in CC_GAM') 94 ENDIF 95C 96 CALL CC_GAM1(DSRHF,GAMMA,SCRM,ISYMCM,WORK(KLAMDA),WORK(KLAMDC), 97 * XLAMDP,XLAMDH,XLAMPC,XLAMHC,ISYMPC, 98 * WORK(KSCR1),WORK(KSCR2),WORK(KSCR3),WORK(KSCR4), 99 * WORK(KSCR5),WORK(KEND2),LWRK2,ISYMD,IDEL,ISYML, 100 * ISYMAG,IOPT) 101C 102 100 CONTINUE 103C 104 RETURN 105 END 106 SUBROUTINE CC_GAM1(DSRHF,GAMMA,SCRM,ISYMCM,XLAM,XLAMC, 107 * XLAMDP,XLAMDH,XLAMPC,XLAMHC,ISYMPC, 108 * SCR1,SCR2,SCR3,SCR4,SCR5,WORK, 109 * LWORK,ISYMD,IDEL,ISYML,ISYMAG,IOPT) 110C 111C Written by Henrik Koch 3-Jan-1994 112C 113C Generalized by Ove Christiansen 18-9-1996 for 114C calculation of Gamma-intermediate for F matrix. 115C F-matrix: IOPT = 2 and ISYMCM is symmetry of SCRM 116C and ISYMPC is symmetry of XLAMPC and XLAMHC. 117C (They should be the same) 118C 119C Purpose: Calculate the gamma intermediate. 120C 121#include "implicit.h" 122 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 123 DIMENSION DSRHF(*),GAMMA(*),SCRM(*),XLAM(*),XLAMC(*) 124 DIMENSION SCR1(*),SCR2(*),SCR3(*),SCR4(*),SCR5(*),WORK(*) 125 DIMENSION XLAMDP(*),XLAMDH(*),XLAMPC(*),XLAMHC(*) 126#include "priunit.h" 127#include "ccorb.h" 128#include "ccsdsym.h" 129C 130C INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 131C 132 ISYMJB = MULD2H(ISYMD,ISYMPC) 133 ISKILJ = MULD2H(ISYMCM,ISYMOP) 134 IF (ISYMPC .NE. ISYMCM) CALL QUIT('Symmetry mismatch in CC_GAM1') 135C 136 ISYMKC = ISYMAG 137C 138 DO 100 L = 1,NRHF(ISYML) 139C 140 KOFF1 = IDSRHF(ISYMAG,ISYML) + NNBST(ISYMAG)*(L - 1) + 1 141C 142 CALL CCSD_SYMSQ(DSRHF(KOFF1),ISYMAG,SCR1) 143C 144 DO 110 ISYMG = 1,NSYM 145C 146 ISYMA = MULD2H(ISYMG,ISYMAG) 147 ISYMK = ISYMA 148 ISYMC = ISYMG 149 ISYMI = ISYMG 150C 151 NBASA = MAX(NBAS(ISYMA),1) 152 NBASG = MAX(NBAS(ISYMG),1) 153 NRHFK = MAX(NRHF(ISYMK),1) 154 NVIRC = MAX(NVIR(ISYMC),1) 155C 156 KOFF2 = ILMRHF(ISYMK) + 1 157 KOFF3 = IAODIS(ISYMA,ISYMG) + 1 158 KOFF4 = IT1AOT(ISYMK,ISYMG) + 1 159C 160 CALL DGEMM('T','N',NRHF(ISYMK),NBAS(ISYMG),NBAS(ISYMA), 161 * ONE,XLAMDP(KOFF2),NBASA,SCR1(KOFF3),NBASA, 162 * ZERO,SCR2(KOFF4),NRHFK) 163C 164 KOFF5 = ILMVIR(ISYMC) + 1 165 KOFF6 = IT1AMT(ISYMK,ISYMC) + 1 166C 167 CALL DGEMM('N','N',NRHF(ISYMK),NVIR(ISYMC),NBAS(ISYMG), 168 * ONE,SCR2(KOFF4),NRHFK,XLAMDH(KOFF5),NBASG, 169 * ZERO,SCR3(KOFF6),NRHFK) 170C 171 KOFF7 = ILMRHF(ISYMI) + 1 172 KOFF8 = IMATIJ(ISYMK,ISYMI) + 1 173C 174 CALL DGEMM('N','N',NRHF(ISYMK),NRHF(ISYMI),NBAS(ISYMG), 175 * ONE,SCR2(KOFF4),NRHFK,XLAMDH(KOFF7),NBASG, 176 * ZERO,SCR4(KOFF8),NRHFK) 177C 178 IF (IOPT .EQ. 2) THEN 179C 180C------------------------------------------------------------- 181C Only for calculating Ik,i-bar,l delta integral. 182C------------------------------------------------------------- 183C 184 ISYMA = MULD2H(ISYMG,ISYMAG) 185 ISYMK = ISYMA 186 ISYMI = MULD2H(ISYMG,ISYMPC) 187C 188 KOFF7 = IGLMRH(ISYMG,ISYMI) + 1 189 KOFF8 = IMATIJ(ISYMK,ISYMI) + 1 190C 191 CALL DGEMM('N','N',NRHF(ISYMK),NRHF(ISYMI),NBAS(ISYMG), 192 * ONE,SCR2(KOFF4),NRHFK,XLAMHC(KOFF7),NBASG, 193 * ZERO,SCR5(KOFF8),NRHFK) 194C 195 ENDIF 196C 197 110 CONTINUE 198C 199 DO 120 ISYMJ = 1,NSYM 200C 201 ISYMLJ = MULD2H(ISYML,ISYMJ) 202 ISYMKI = MULD2H(ISYMLJ,ISKILJ) 203 ISYDVI = ISYMD 204 ISYDVJ = MULD2H(ISYDVI,ISYMJ) 205 ISYMCI = MULD2H(ISYMCM,ISYDVJ) 206C 207 IF (ISYMKI .GT. ISYMLJ) GOTO 120 208C 209 KLC = IGLMRH(ISYMD,ISYMJ) + IDEL - IBAS(ISYMD) 210 CALL DCOPY(NRHF(ISYMJB),XLAMHC(KLC),NBAS(ISYMD), 211 * XLAMC,1) 212C 213 KSCR5 = 1 214 KEND1 = KSCR5 + NMATIJ(ISYMKI) 215C 216 DO 130 J = 1,NRHF(ISYMJ) 217C 218 DO 140 ISYMI = 1,NSYM 219C 220 ISYMC = MULD2H(ISYMI,ISYMCI) 221 ISYMK = MULD2H(ISYMI,ISYMKI) 222C 223 NVIRC = MAX(NVIR(ISYMC),1) 224 NRHFK = MAX(NRHF(ISYMK),1) 225C 226 KOFF2 = IT1AMT(ISYMK,ISYMC) + 1 227 KOFF3 = IT2BCD(ISYMCI,ISYMJ) 228 * + NT1AM(ISYMCI)*(J - 1) 229 * + IT1AM(ISYMC,ISYMI) + 1 230 KOFF4 = KSCR5 + IMATIJ(ISYMK,ISYMI) 231C 232 CALL DGEMM('N','N',NRHF(ISYMK),NRHF(ISYMI), 233 * NVIR(ISYMC),ONE,SCR3(KOFF2),NRHFK, 234 * SCRM(KOFF3),NVIRC,ZERO,WORK(KOFF4),NRHFK) 235C 236 140 CONTINUE 237C 238 IF ( IOPT .EQ. 2 ) THEN 239 IF (ISYMJ .EQ. ISYMD) THEN 240 CALL DAXPY(NMATIJ(ISYMKI),XLAM(J),SCR5,1, 241 * WORK(KSCR5),1) 242 ENDIF 243 IF (MULD2H(ISYMJ,ISYMD).EQ.ISYMPC) THEN 244 CALL DAXPY(NMATIJ(ISYMKI),XLAMC(J),SCR4,1, 245 * WORK(KSCR5),1) 246 ENDIF 247 ELSE 248 IF (ISYMJ .EQ. ISYMD) THEN 249 CALL DAXPY(NMATIJ(ISYMKI),XLAM(J),SCR4,1, 250 * WORK(KSCR5),1) 251 ENDIF 252 ENDIF 253C 254 NLJ = IMATIJ(ISYML,ISYMJ) + NRHF(ISYML)*(J - 1) + L 255C 256 IF (ISKILJ .EQ. 1) THEN 257 KKILJ = IGAMMA(ISYMKI,ISYMLJ) + NLJ*(NLJ-1)/2 258 DO 150 NKI = 1,NLJ 259C 260 KOFF = KSCR5 + NKI - 1 261 NKILJ = KKILJ + NKI 262 GAMMA(NKILJ) = GAMMA(NKILJ) + WORK(KOFF) 263C 264 150 CONTINUE 265 ELSE 266 KKILJ = IGAMMA(ISYMKI,ISYMLJ) 267 * + NMATIJ(ISYMKI)*(NLJ - 1) 268 DO 160 NKI = 1,NMATIJ(ISYMKI) 269C 270 KOFF = KSCR5 + NKI - 1 271 NKILJ = KKILJ + NKI 272 GAMMA(NKILJ) = GAMMA(NKILJ) + WORK(KOFF) 273C 274 160 CONTINUE 275 END IF 276C 277 130 CONTINUE 278 120 CONTINUE 279C 280 IF (IOPT .EQ. 2) THEN 281C 282C------------------------------------ 283C Extra F-matrix term section. 284C------------------------------------ 285C 286 ENDIF 287C 288 100 CONTINUE 289C 290 RETURN 291 END 292C /* Deck cc_zwvi */ 293 SUBROUTINE CC_ZWVI(ZINT,CTR2,ISYMC2 ,TINT,ISYTIN, 294 * WORK,LWORK,IOPT) 295C 296C Written by Asger Halkier 26/10 - 1995 297C 298C Version: 1.0 299C 300C Purpose: To calculate the intermediates entering some of the 301C terms in the 2.1-block. 302C 303C IOPT equals 2 if transposition of the occupied indices of 304C TINT intermediate is needed, and 1 if not! 305C 306C Ove Christiansen 1-10-1996: 307C general symmetry of ctr2 (isymc2) 308C 309#include "implicit.h" 310 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 311 DIMENSION ZINT(*), CTR2(*), TINT(*), WORK(LWORK) 312#include "priunit.h" 313#include "ccorb.h" 314#include "ccsdsym.h" 315#include "cclr.h" 316C 317C----------------------------- 318C Initialize result array. 319C----------------------------- 320C 321 ISYZIN = MULD2H(ISYMC2,ISYTIN) 322C 323 CALL DZERO(ZINT,NT2BCD(ISYZIN)) 324C 325C----------------------------------------------------- 326C Transpose occupied indices of TINT if requested. 327C----------------------------------------------------- 328C 329 IF (IOPT .EQ. 2) THEN 330C 331 CALL CCLT_P21I(TINT,ISYTIN,WORK,LWORK, 332 & IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR) 333C 334 ENDIF 335C 336C------------------------ 337C Do the calculation. 338C------------------------ 339C 340 DO 100 ISYMDK = 1,NSYM 341C 342 ISYMEI = MULD2H(ISYMDK,ISYMC2) 343 ISYMJ = MULD2H(ISYMDK,ISYTIN) 344C 345 KOFF1 = IT2SQ(ISYMEI,ISYMDK) + 1 346 KOFF2 = IT2BCD(ISYMDK,ISYMJ) + 1 347 KOFF3 = IT2BCD(ISYMEI,ISYMJ) + 1 348C 349 NTOTEI = MAX(NT1AM(ISYMEI),1) 350 NTOTDK = MAX(NT1AM(ISYMDK),1) 351C 352 CALL DGEMM('N','N',NT1AM(ISYMEI),NRHF(ISYMJ),NT1AM(ISYMDK), 353 * ONE,CTR2(KOFF1),NTOTEI,TINT(KOFF2),NTOTDK,ZERO, 354 * ZINT(KOFF3),NTOTEI) 355C 356 100 CONTINUE 357C 358C------------------------------------------- 359C Restore TINT intermediate if necessary 360C------------------------------------------- 361C 362 IF (IOPT .EQ. 2) THEN 363C 364 CALL CCLT_P21I(TINT,ISYTIN,WORK,LWORK, 365 & IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR) 366C 367 ENDIF 368C 369 RETURN 370 END 371C /* Deck cc_ti */ 372 SUBROUTINE CC_TI(TINT,ISYTIN,T2AM,ISYMT2,XLAMDH,ISYMLH, 373 * WORK,LWORK,IDEL,ISYMD) 374C 375C Written by Asger Halkier 26/10 - 1995 376C 377C Version: 1.0 378C 379C Purpose: To calculate the T-intermediate entering some of the 380C terms in the 2.1-block. 381C 382C Ove Christiansen 30-9-1996: Generalised symmetry of T2AM: ISYMT2 383C Still it is assumed that XLAMDH is total symmetric. 384C 385#include "implicit.h" 386 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 387 INTEGER ISYTIN, ISYMT2, ISYMLH, IDEL, ISYMD 388 DIMENSION TINT(*), T2AM(*), XLAMDH(*), WORK(LWORK) 389#include "priunit.h" 390#include "ccorb.h" 391#include "ccsdsym.h" 392#include "cclr.h" 393C 394 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 395C 396 CALL QENTER('CC_TI') 397C 398C----------------------------- 399C Initialize result array. 400C----------------------------- 401C 402 ISYMF = MULD2H(ISYMD,ISYMLH) 403 IF (ISYTIN.NE.MULD2H(ISYMT2,ISYMF)) THEN 404 WRITE(LUPRI,*) 'ISYTIN,ISYMT2,ISYMLH,ISYMD:', 405 & ISYTIN,ISYMT2,ISYMLH,ISYMD 406 CALL QUIT('SYMMETRY MISMATCH IN CC_TI') 407 END IF 408C 409 CALL DZERO (TINT,NT2BCD(ISYTIN)) 410C 411C---------------------------------- 412C Work space allocation one. 413C---------------------------------- 414C 415 KLAHF = 1 416 KEND1 = KLAHF + NVIR(ISYMF) 417 LWRK1 = LWORK - KEND1 418C 419 IF (LWRK1 .LT. 0) THEN 420 CALL QUIT('Insufficient work space in CC_TI') 421 ENDIF 422C 423C------------------------------------------ 424C Copy vector out of lamda hole matrix. 425C------------------------------------------ 426C 427 KOFF1 = IGLMVI(ISYMD,ISYMF) + IDEL - IBAS(ISYMD) 428C 429 CALL DCOPY(NVIR(ISYMF),XLAMDH(KOFF1),NBAS(ISYMD),WORK(KLAHF),1) 430C 431 DO 100 ISYMDK = 1,NSYM 432C 433 ISYMFJ = MULD2H(ISYMDK,ISYMT2) 434 ISYMJ = MULD2H(ISYMFJ,ISYMF) 435C 436 IF (NRHF(ISYMJ) .EQ. 0) GOTO 100 437C 438C---------------------------------- 439C Work space allocation two. 440C---------------------------------- 441C 442 KT2SM = KEND1 443 KEND2 = KT2SM + NT1AM(ISYMDK)*NVIR(ISYMF) 444 LWRK2 = LWORK - KEND2 445C 446 IF (LWRK2 .LT. 0) THEN 447 CALL QUIT('Insufficient work space in CC_TI') 448 ENDIF 449C 450 DO 110 J = 1,NRHF(ISYMJ) 451C 452C--------------------------------------------- 453C Copy submatrix out of packed T2AM. 454C--------------------------------------------- 455C 456 DO 120 F = 1,NVIR(ISYMF) 457C 458 NFJ = IT1AM(ISYMF,ISYMJ) + NVIR(ISYMF)*(J - 1) + F 459C 460 IF (ISYMDK .EQ. ISYMFJ) THEN 461C 462 DO 130 NDK = 1,NT1AM(ISYMDK) 463C 464 NDKFJ = IT2AM(ISYMDK,ISYMFJ) + INDEX(NDK,NFJ) 465 NDKF = KT2SM + NT1AM(ISYMDK)*(F - 1) + NDK - 1 466C 467 WORK(NDKF) = T2AM(NDKFJ) 468C 469 130 CONTINUE 470C 471 ELSE IF (ISYMDK .LT. ISYMFJ) THEN 472C 473 NDKFJ = IT2AM(ISYMDK,ISYMFJ) 474 * + NT1AM(ISYMDK)*(NFJ - 1) + 1 475 NDKF = KT2SM + NT1AM(ISYMDK)*(F - 1) 476C 477 CALL DCOPY(NT1AM(ISYMDK),T2AM(NDKFJ),1,WORK(NDKF),1) 478C 479 ELSE IF (ISYMDK .GT. ISYMFJ) THEN 480C 481 NDKFJ = IT2AM(ISYMFJ,ISYMDK) + NFJ 482 NDKF = KT2SM + NT1AM(ISYMDK)*(F - 1) 483C 484 CALL DCOPY(NT1AM(ISYMDK),T2AM(NDKFJ),NT1AM(ISYMFJ), 485 * WORK(NDKF),1) 486C 487 ENDIF 488C 489 120 CONTINUE 490C 491C----------------------------------------------------- 492C Contraction of T2AM-submatrix with XLAMDH. 493C----------------------------------------------------- 494C 495 KOFF2 = IT2BCD(ISYMDK,ISYMJ) + NT1AM(ISYMDK)*(J - 1) + 1 496C 497 NTOTDK = MAX(NT1AM(ISYMDK),1) 498C 499 CALL DGEMV('N',NT1AM(ISYMDK),NVIR(ISYMF),ONE,WORK(KT2SM), 500 * NTOTDK,WORK(KLAHF),1,ZERO,TINT(KOFF2),1) 501C 502 110 CONTINUE 503 100 CONTINUE 504C 505 CALL QEXIT('CC_TI') 506C 507 RETURN 508 END 509C /* Deck cc_21a */ 510 SUBROUTINE CC_21A(RHO1,DSRHF,YMAT,ISYMY,YDEN,ISYMYD, 511 * XLAMDH,XLAMDP,ISYMLP,WORK, 512 * LWORK,IDEL,ISYMD) 513C 514C Written by Asger Halkier 4/10 - 1995. 515C 516C Version: 1.0 517C 518C Purpose: To calculate the 21A contribution to rho1! 519C 520C Ove Christiansen 1-10-1996: 521C Generalization to general symmetry YMAT and YDEN 522C for F-matrix. 523C 524#include "implicit.h" 525 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 526 DIMENSION RHO1(*), DSRHF(*), YMAT(*), YDEN(*), XLAMDH(*), 527 * XLAMDP(*), WORK(LWORK) 528#include "priunit.h" 529#include "ccorb.h" 530#include "ccsdsym.h" 531#include "cclr.h" 532C 533 ISYRES = MULD2H(ISYMYD,ISYMOP) 534 IF (MULD2H(ISYMY,ISYMLP).NE.ISYRES) THEN 535 WRITE(LUPRI,*) ' Symmetry mismatch in CC_21A' 536 CALL QUIT(' Symmetry mismatch in CC_21A') 537 ENDIF 538 ISYINT = MULD2H(ISYMD,ISYMOP) 539C 540C================================ 541C Calculate the coulomb part. 542C================================ 543C 544 ISYMA = ISYMD 545 ISYMI = MULD2H(ISYMA,ISYRES) 546 ISALBE = MULD2H(ISYMI,ISYINT) 547C 548C---------------------------------------------------------------------- 549C Work space allocation 1 & options for contraction with integrals. 550C---------------------------------------------------------------------- 551C 552 KINT = 1 553 KVECIS = KINT + N2BST(ISALBE) 554 KENSMA = KVECIS + NRHF(ISYMI) 555 KVECIB = KINT + NRHF(ISYMI)*N2BST(ISALBE) 556 KENBIG = KVECIB + NRHF(ISYMI) 557 LWRSMA = LWORK - KENSMA 558 LWRBIG = LWORK - KENBIG 559C 560 IF (LWRSMA .LT. 0) THEN 561 CALL QUIT('Insufficient work space in CC_21A') 562 ENDIF 563C 564 IF (LWRBIG .LE. 0) THEN 565C 566 DO 110 I = 1,NRHF(ISYMI) 567C 568 KOFF1 = IDSRHF(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 569C 570 CALL CCSD_SYMSQ(DSRHF(KOFF1),ISALBE,WORK(KINT)) 571C 572 KOFF2 = KVECIS + I - 1 573C 574 WORK(KOFF2) = DDOT(N2BST(ISALBE),YDEN,1, 575 * WORK(KINT),1) 576C 577 110 CONTINUE 578C 579 ELSE 580C 581 DO 120 I = 1,NRHF(ISYMI) 582C 583 KOFF1 = IDSRHF(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 584 KOFF2 = KINT + N2BST(ISALBE)*(I - 1) 585C 586 CALL CCSD_SYMSQ(DSRHF(KOFF1),ISALBE,WORK(KOFF2)) 587C 588 120 CONTINUE 589C 590 NTALBE = MAX(N2BST(ISALBE),1) 591C 592 CALL DGEMV('T',N2BST(ISALBE),NRHF(ISYMI),TWO,WORK(KINT),NTALBE, 593 * YDEN,1,ZERO,WORK(KVECIB),1) 594C 595 CALL DCOPY(NRHF(ISYMI),WORK(KVECIB),1,WORK(KVECIS),1) 596C 597 ENDIF 598C 599C------------------------------------------------- 600C Scale with XLAMDH-element and add to result. 601C------------------------------------------------- 602C 603 DO 130 A = 1,NVIR(ISYMA) 604C 605 KOFF1 = ILMVIR(ISYMA) + NBAS(ISYMD)*(A - 1) 606 * + IDEL - IBAS(ISYMD) 607 KOFF2 = IT1AM(ISYMA,ISYMI) + A 608C 609 CALL DAXPY(NRHF(ISYMI),XLAMDH(KOFF1),WORK(KVECIS),1, 610 * RHO1(KOFF2),NVIR(ISYMA)) 611C 612 130 CONTINUE 613C 614C================================= 615C Calculate the exchange part. 616C================================= 617C 618 ISYME = ISYMD 619 ISYMF = MULD2H(ISYME,ISYMY) 620 ISYMAL = MULD2H(ISYMF,ISYMLP) 621 ISYBEI = MULD2H(ISYMAL,ISYINT) 622C 623 DO 140 ISYMI = 1,NSYM 624C 625 ISYMBE = MULD2H(ISYMI,ISYBEI) 626 ISALBE = MULD2H(ISYMBE,ISYMAL) 627 ISYMA = MULD2H(ISYMI,ISYRES) 628C 629C-------------------------------- 630C Work space allocation 2. 631C-------------------------------- 632C 633 KINT = 1 634 KSCR1 = KINT + N2BST(ISALBE) 635 KSCR2 = KSCR1 + NBAS(ISYMAL)*NVIR(ISYMA) 636 KSCR3 = KSCR2 + NVIR(ISYMA)*NVIR(ISYMF) 637 KEND1 = KSCR3 + NVIR(ISYMA)*NVIR(ISYME) 638 LWRK1 = LWORK - KEND1 639C 640 IF (LWRK1 .LT. 0) THEN 641 CALL QUIT('Insufficient work space in CC_21A') 642 ENDIF 643C 644 DO 150 I = 1,NRHF(ISYMI) 645C 646 KOFFPI = IDSRHF(ISALBE,ISYMI) + NNBST(ISALBE)*(I - 1) + 1 647C 648 CALL CCSD_SYMSQ(DSRHF(KOFFPI),ISALBE,WORK(KINT)) 649C 650C------------------------------------------- 651C Calculate intermediate matrices. 652C------------------------------------------- 653C 654 KOFFUI = KINT + IAODIS(ISYMAL,ISYMBE) 655 KOFF1 = ILMVIR(ISYMA) + 1 656C 657 NTOTAL = MAX(NBAS(ISYMAL),1) 658 NTOTBE = MAX(NBAS(ISYMBE),1) 659C 660 CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMA),NBAS(ISYMBE), 661 * ONE,WORK(KOFFUI),NTOTAL,XLAMDH(KOFF1),NTOTBE, 662 * ZERO,WORK(KSCR1),NTOTAL) 663C 664 KOFF2 = IGLMVI(ISYMAL,ISYMF) + 1 665C 666 NTOTAL = MAX(NBAS(ISYMAL),1) 667 NTOTF = MAX(NVIR(ISYMF),1) 668C 669 CALL DGEMM('T','N',NVIR(ISYMF),NVIR(ISYMA),NBAS(ISYMAL), 670 * ONE,XLAMDP(KOFF2),NTOTAL,WORK(KSCR1),NTOTAL, 671 * ZERO,WORK(KSCR2),NTOTF) 672C 673 KOFF3 = IMATAB(ISYME,ISYMF) + 1 674C 675 NTOTE = MAX(NVIR(ISYME),1) 676 NTOTF = MAX(NVIR(ISYMF),1) 677C 678 CALL DGEMM('N','N',NVIR(ISYME),NVIR(ISYMA),NVIR(ISYMF), 679 * ONE,YMAT(KOFF3),NTOTE,WORK(KSCR2),NTOTF, 680 * ZERO,WORK(KSCR3),NTOTE) 681C 682C------------------------------------------- 683C Add contribution to RHO1 vector. 684C------------------------------------------- 685C 686 KOFF4 = ILMVIR(ISYME) + IDEL - IBAS(ISYMD) 687 KOFF5 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 688C 689 NTOTE = MAX(NVIR(ISYME),1) 690C 691 CALL DGEMV('T',NVIR(ISYME),NVIR(ISYMA),-ONE,WORK(KSCR3), 692 * NTOTE,XLAMDH(KOFF4),NBAS(ISYMD),ONE, 693 * RHO1(KOFF5),1) 694C 695 150 CONTINUE 696C 697 140 CONTINUE 698C 699 RETURN 700 END 701C /* Deck cc_yd */ 702 SUBROUTINE CC_YD(YDEN,YMAT,ISYMY,XLAMDH,XLAMDP,ISYMLP, 703 * WORK,LWORK) 704C 705C Written by Asger Halkier 8/12 - 1995. 706C 707C Version: 1.0 708C 709C Purpose: To transform the Y-matrix to AO-basis! 710C 711C Ove Christiansen 1-10-1996: 712C General symmetry of YMAT (ISYMY) and 713C XLAMDP (ISYMLP) 714C XLAMDH is assumed to have the symmetry ISYMOP 715C (it is a hole-virtuel index.) 716C 717#include "implicit.h" 718 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 719 DIMENSION YMAT(*), YDEN(*), XLAMDH(*), XLAMDP(*), WORK(LWORK) 720#include "priunit.h" 721#include "ccorb.h" 722#include "ccsdsym.h" 723#include "cclr.h" 724C 725 ISYMYD = MULD2H(ISYMLP,ISYMY) 726C 727C------------------------------------ 728C Transform Y-matrix to AO-basis. 729C------------------------------------ 730C 731 DO 100 ISYMAL = 1,NSYM 732C 733 ISYMF = MULD2H(ISYMAL,ISYMLP) 734 ISYMBE = MULD2H(ISYMAL,ISYMYD) 735 ISYME = MULD2H(ISYMBE,ISYMOP) 736C 737 LWRKCH = LWORK - NBAS(ISYMBE)*NVIR(ISYMF) 738C 739 IF (LWRKCH .LT. 0) THEN 740 CALL QUIT('Insufficient work space in CC_21A') 741 ENDIF 742C 743 KOFF1 = ILMVIR(ISYME) + 1 744 KOFF2 = IMATAB(ISYME,ISYMF) + 1 745C 746 NTOTBE = MAX(NBAS(ISYMBE),1) 747 NTOTE = MAX(NVIR(ISYME),1) 748C 749 CALL DGEMM('N','N',NBAS(ISYMBE),NVIR(ISYMF),NVIR(ISYME), 750 * ONE,XLAMDH(KOFF1),NTOTBE,YMAT(KOFF2),NTOTE, 751 * ZERO,WORK,NTOTBE) 752C 753 KOFF1 = IGLMVI(ISYMAL,ISYMF) + 1 754 KOFF2 = IAODIS(ISYMAL,ISYMBE) + 1 755C 756 NTOTAL = MAX(NBAS(ISYMAL),1) 757 NTOTBE = MAX(NBAS(ISYMBE),1) 758C 759 CALL DGEMM('N','T',NBAS(ISYMAL),NBAS(ISYMBE),NVIR(ISYMF), 760 * ONE,XLAMDP(KOFF1),NTOTAL,WORK,NTOTBE, 761 * ZERO,YDEN(KOFF2),NTOTAL) 762C 763 100 CONTINUE 764C 765 RETURN 766 END 767C /* Deck cc_21h */ 768 SUBROUTINE CC_21H(RHO1,ISYRHO,VINT,WINT,ZINT,ISYVWZ,X3OINT, 769 * ISYINT,WORK,LWORK,ISYMD) 770C 771C Written by Asger Halkier 30/10 - 1995 772C 773C Version: 1.0 774C 775C Purpose: To calculate the 21H contribution to RHO1! 776C 777C Ove Christiansen 2-10-1996: 778C Generalisation to general symmetry of intermediates 779C as well as general symmetry of integrals. 780C ISYVWZ is symmetry of V,W, and Z intermediates. 781C ISYINT is symmetry of integrals. 782C 783#include "implicit.h" 784 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 785 DIMENSION RHO1(*), VINT(*), WINT(*), ZINT(*), X3OINT(*), 786 & WORK(LWORK) 787#include "priunit.h" 788#include "ccorb.h" 789#include "ccsdsym.h" 790#include "cclr.h" 791C 792 ISYJLI = MULD2H(ISYINT,ISYMD) 793 ISYRES = MULD2H(ISYJLI,ISYVWZ) 794 IF (ISYRES .NE. ISYRHO) THEN 795 CALL QUIT('Symmetry mismatch in CC_21H') 796 ENDIF 797C 798C------------------------------- 799C Work space allocation one. 800C------------------------------- 801C 802 KVZINT = 1 803 KEND1 = KVZINT + NT2BCD(ISYVWZ) 804 LWRK1 = LWORK - KEND1 805C 806 IF (LWRK1 .LT. 0) THEN 807 CALL QUIT('Insufficient work space in CC_21H') 808 ENDIF 809C 810C--------------------------------------------------- 811C Resort and add together V- and Z-intermediate. 812C--------------------------------------------------- 813C 814 CALL CCLT_S21I(WORK(KVZINT),VINT,ZINT,ISYVWZ,ONE) 815C 816 DO 100 ISYMA = 1,NSYM 817C 818 ISYMI = MULD2H(ISYMA,ISYRES) 819 ISYMJL = MULD2H(ISYMI,ISYJLI) 820C 821C-------------------------------------- 822C Calculate the VZ-contribution. 823C-------------------------------------- 824C 825 KOFF1 = KVZINT + IT2AIJ(ISYMA,ISYMJL) 826 KOFF2 = IMAIJK(ISYMJL,ISYMI) + 1 827 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 828C 829 NTOTA = MAX(NVIR(ISYMA),1) 830 NTOTJL = MAX(NMATIJ(ISYMJL),1) 831C 832 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMJL), 833 & ONE,WORK(KOFF1),NTOTA,X3OINT(KOFF2),NTOTJL,ONE, 834 & RHO1(KOFF3),NTOTA) 835C 836 100 CONTINUE 837C 838C------------------------------- 839C Work space allocation two. 840C------------------------------- 841C 842 KWZINT = 1 843 KI3OTR = KWZINT + NT2BCD(ISYVWZ) 844 KEND2 = KI3OTR + NMAIJK(ISYJLI) 845 LWRK2 = LWORK - KEND2 846C 847 IF (LWRK2 .LT. 0) THEN 848 CALL QUIT('Insufficient work space in CC_21H') 849 ENDIF 850C 851C--------------------------------------------------------- 852C Prepare intermediates and integrals for contraction. 853C--------------------------------------------------------- 854C 855 CALL CCLT_S21I(WORK(KWZINT),WINT,ZINT,ISYVWZ,-TWO) 856C 857 CALL CCLT_PI3O(WORK(KI3OTR),X3OINT,ISYJLI) 858C 859 DO 110 ISYMA = 1,NSYM 860C 861 ISYMI = MULD2H(ISYMA,ISYRES) 862 ISYMJL = MULD2H(ISYMI,ISYJLI) 863C 864C---------------------------------- 865C Calculate WZ-contribution. 866C---------------------------------- 867C 868 KOFF1 = KWZINT + IT2AIJ(ISYMA,ISYMJL) 869 KOFF2 = KI3OTR + IMAIJK(ISYMJL,ISYMI) 870 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 871C 872 NTOTA = MAX(NVIR(ISYMA),1) 873 NTOTJL = MAX(NMATIJ(ISYMJL),1) 874C 875 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NMATIJ(ISYMJL), 876 & ONE,WORK(KOFF1),NTOTA,WORK(KOFF2),NTOTJL,ONE, 877 & RHO1(KOFF3),NTOTA) 878C 879 110 CONTINUE 880C 881 RETURN 882 END 883C /* Deck cc_21g */ 884 SUBROUTINE CC_21G(RHO1,XMINT,ISYMM,XLAMDH,WORK,LWORK, 885 * ISYINT,LUO3,O3FIL) 886C 887C Written by Asger Halkier 30/10 - 1995 888C 889C Version: 1.0 890C 891C Purpose: To calculate the 21G contribution to RHO1. 892C 893C Ove Christiansen 3-10-1996: 894C General symmetries for F-matrix 895C isyint is symmetry of integrals. 896C 897C 898#include "implicit.h" 899 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 900 DIMENSION RHO1(*), XMINT(*), XLAMDH(*), WORK(LWORK) 901 CHARACTER O3FIL*(*) 902#include "priunit.h" 903#include "ccorb.h" 904#include "ccsdsym.h" 905#include "cclr.h" 906C 907 ISYRES = MULD2H(ISYMM,ISYINT) 908C 909 DO 100 ISYMA = 1,NSYM 910C 911 ISYMD = ISYMA 912 ISYMI = MULD2H(ISYMA,ISYRES) 913 ISYJKL = MULD2H(ISYMD,ISYINT) 914C 915 IF (NBAS(ISYMD) .EQ. 0) GOTO 100 916C 917C---------------------------------- 918C Work space allocation one. 919C---------------------------------- 920C 921 KMOINT = 1 922 KAOINT = KMOINT + NVIR(ISYMA)*NMAIJK(ISYJKL) 923 KEND1 = KAOINT + NBAS(ISYMD)*NMAIJK(ISYJKL) 924 LWRK1 = LWORK - KEND1 925C 926 IF (LWRK1 .LT. 0) THEN 927 CALL QUIT('Insufficient work space in CC_21G') 928 ENDIF 929C 930C------------------------------------------- 931C Read integrals (jk|ldel) from disc. 932C------------------------------------------- 933C 934 NTOT = NBAS(ISYMD)*NMAIJK(ISYJKL) 935 IOFF = I3ODEL(ISYJKL,ISYMD) + 1 936C 937 CALL GETWA2(LUO3,O3FIL,WORK(KAOINT),IOFF,NTOT) 938C 939C----------------------------------------------------- 940C Transform AO integral index to virtual space. 941C----------------------------------------------------- 942C 943 KOFF1 = ILMVIR(ISYMA) + 1 944C 945 NTOJKL = MAX(NMAIJK(ISYJKL),1) 946 NTOTD = MAX(NBAS(ISYMD),1) 947C 948 CALL DGEMM('N','N',NMAIJK(ISYJKL),NVIR(ISYMA),NBAS(ISYMD), 949 & ONE,WORK(KAOINT),NTOJKL,XLAMDH(KOFF1),NTOTD,ZERO, 950 & WORK(KMOINT),NTOJKL) 951C 952C------------------------------------------------------------ 953C Contraction with M-intermediate & storage in result. 954C------------------------------------------------------------ 955C 956 KOFF2 = I3ORHF(ISYJKL,ISYMI) + 1 957 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 958C 959 NTOJKL = MAX(NMAIJK(ISYJKL),1) 960 NTOTA = MAX(NVIR(ISYMA),1) 961C 962 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NMAIJK(ISYJKL), 963 & ONE,WORK(KMOINT),NTOJKL,XMINT(KOFF2),NTOJKL,ONE, 964 & RHO1(KOFF3),NTOTA) 965C 966 100 CONTINUE 967C 968 RETURN 969 END 970