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 cc2_den */ 20 SUBROUTINE CC2_DEN(ETAAI,ZKDIA,WORK,LWORK,IOPT) 21C 22C Written by Asger Halkier & Sonia Coriani 12/1 - 2000 23C 24C Version: 1.0 25C 26C PURPOSE: 27C 28C IOPT=1: To calculate the pp, ab, & ij parts of kappa-bar-0 29C that do not need the solution of any coupled equations 30C in the case of (relaxed) CC2 calculations. 31C ZKDIA holds all the blocks pq in the following order: 32C ij, ab, ai, ia; and these are stored full blocks after 33C each other. 34C 35C IOPT=2: To calculate the complete right-hand-side for the 36C kappa-bar-0 ai block (eta_ai) - this includes both the 37C "direct" terms calculated from densities and t1-transformed 38C integrals, as well as the indirect terms coming from 39C kappa-bar-0_ab and kappa-bar-0_ij. The result is stored in 40C ETAAI. 41C 42C IMPORTANT NOTE: NO FROZEN CORE SO FAR!!! 43C 44#include "implicit.h" 45#include "priunit.h" 46#include "dummy.h" 47#include "ccinftap.h" 48#include "maxash.h" 49#include "maxorb.h" 50#include "mxcent.h" 51#include "aovec.h" 52#include "iratdef.h" 53 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 54 DIMENSION ETAAI(*), ZKDIA(*), WORK(LWORK) 55 DIMENSION INDEXA(MXCORB_CC) 56#include "ccorb.h" 57CCN !#include "infind.h" not compatible with R12! 58#include "ccisao.h" 59#include "r12int.h" 60#include "blocks.h" 61#include "ccfield.h" 62#include "ccsdinp.h" 63#include "ccsdsym.h" 64#include "ccsdio.h" 65#include "distcl.h" 66#include "cbieri.h" 67#include "eritap.h" 68#include "cclr.h" 69#include "ccfro.h" 70C 71 CHARACTER MODEL*10 72 LOGICAL LOCDBG 73 PARAMETER (LOCDBG=.FALSE.) 74C 75 CALL QENTER('CC2_DEN') 76C 77 TIMETO = SECOND() 78 TIMDEN = ZERO 79 TIMHE2 = ZERO 80 TIMRES = ZERO 81 TIRDAO = ZERO 82C 83 IF ((IPRINT .GT. 3) .AND. (IOPT .EQ. 1)) THEN 84 CALL HEADER('Calculating diagonal blocks of zeta-kappa-0',-1) 85 ENDIF 86 IF ((IPRINT .GT. 3) .AND. (IOPT .EQ. 2)) THEN 87 CALL HEADER('Calculating right-hand-side for kappa-bar_ai',-1) 88 ENDIF 89C 90C------------------------------------------------------------------ 91C Both t-vectors and tbar-vectors (zeta) are totally symmetric. 92C------------------------------------------------------------------ 93C 94 ISYMTR = 1 95 ISYMOP = 1 96C 97C------------------------------- 98C Work space allocation one. 99C------------------------------- 100C 101 KCMO = 1 102 KETAIJ = KCMO + NLAMDS 103 KETAAB = KETAIJ + NMATIJ(1) 104 KT2AM = KETAAB + NMATAB(1) 105 KZ2AM = KT2AM + NT2AMX 106 KLAMDP = KZ2AM + NT2SQ(1) 107 KLAMDH = KLAMDP + NLAMDT 108 KT1AM = KLAMDH + NLAMDT 109 KZ1AM = KT1AM + NT1AMX 110 KXMAT = KZ1AM + NT1AMX 111 KYMAT = KXMAT + NMATIJ(1) 112 KEND1 = KYMAT + NMATAB(1) 113 LWRK1 = LWORK - KEND1 114C 115 IF (LWRK1 .LT. 0) THEN 116 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 117 CALL QUIT( 118 * 'Insufficient memory for initial allocation in CC2_DEN') 119 ENDIF 120C 121C---------------------------------------------------------- 122C Read MO-coefficients from interface file and reorder. 123C---------------------------------------------------------- 124C 125 LUSIFC = -1 126 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 127 * IDUMMY,.FALSE.) 128 REWIND LUSIFC 129 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 130 READ (LUSIFC) 131 READ (LUSIFC) 132 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 133 CALL GPCLOSE (LUSIFC,'KEEP') 134C 135C 136 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 137C 138C---------------------------------------- 139C Read zero'th order zeta amplitudes. 140C---------------------------------------- 141C 142 IOPTRW = 3 143 CALL CC_RDRSP('L0',0,1,IOPTRW,MODEL,WORK(KZ1AM),WORK(KZ2AM)) 144C 145C-------------------------------- 146C Square up zeta2 amplitudes. 147C-------------------------------- 148C 149 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 150 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 151C 152C------------------------------------------- 153C Read zero'th order cluster amplitudes. 154C------------------------------------------- 155C 156 IOPTRW = 3 157 CALL CC_RDRSP('R0',0,1,IOPTRW,MODEL,WORK(KT1AM),WORK(KT2AM)) 158C 159C---------------------------------- 160C Calculate the lambda matrices. 161C---------------------------------- 162C 163 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 164 * LWRK1) 165C 166C-------------------------------------------------------- 167C Calculate X-intermediate of tbar- and t-amplitudes. 168C-------------------------------------------------------- 169C 170 IF (IOPT .EQ. 1) THEN 171 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 172 * WORK(KEND1),LWRK1) 173C 174C-------------------------------------------------------- 175C Calculate Y-intermediate of tbar- and t-amplitudes. 176C-------------------------------------------------------- 177C 178 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 179 * WORK(KEND1),LWRK1) 180C 181C-------------------------------------------------------------------------- 182C Calculate the diagonal elements ZK0(ii) = -X(ii) and ZK0(aa) = Y(aa). 183C-------------------------------------------------------------------------- 184C 185 DO 105 ISYMI = 1,NSYM 186 DO 115 I = 1,NRHF(ISYMI) 187C 188 NII = IMATIJ(ISYMI,ISYMI) + NRHF(ISYMI)*(I - 1) + I 189C 190 ZKDIA(NII) = -WORK(KXMAT + NII - 1) 191C 192 115 CONTINUE 193 105 CONTINUE 194C 195 DO 125 ISYMA = 1,NSYM 196 DO 135 A = 1,NVIR(ISYMA) 197C 198 NAA = IMATAB(ISYMA,ISYMA) + NVIR(ISYMA)*(A - 1) + A 199C 200 ZKDIA(NMATIJ(1) + NAA) = WORK(KYMAT + NAA - 1) 201C 202 135 CONTINUE 203 125 CONTINUE 204 ENDIF 205C 206C----------------------------------------------- 207C Set up 2C-E of cluster amplitudes and save 208C in KT2AM, as we only need T(2c-e) below. 209C----------------------------------------------- 210C 211 ISYOPE = 1 212 IOPTTCME = 1 213 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME) 214 KT2AMT = KT2AM !for safety 215C 216C------------------------------- 217C Work space allocation one. 218C Note that D(ai) = ZETA(ai) 219C and both D(ia) and h(ia) 220C are stored transposed! 221C------------------------------- 222C 223 KONEAI = KZ1AM 224 KONEAB = KONEAI + NT1AMX 225 KONEIJ = KONEAB + NMATAB(1) 226 KONEIA = KONEIJ + NMATIJ(1) 227 KONINT = KONEIA + NT1AMX 228 KINTAI = KONINT + N2BST(ISYMOP) 229 KINTAB = KINTAI + NT1AMX 230 KINTIJ = KINTAB + NMATAB(1) 231 KINTIA = KINTIJ + NMATIJ(1) 232 KINABT = KINTIA + NT1AMX 233 KINIJT = KINABT + NMATAB(1) 234 KD1ABT = KINIJT + NMATIJ(1) 235 KD1IJT = KD1ABT + NMATAB(1) 236 KEND1 = KD1IJT + NMATIJ(1) 237 LWRK1 = LWORK - KEND1 238C 239 IF (LWRK1 .LT. 0) THEN 240 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 241 CALL QUIT('Insufficient memory for allocation 1 CC2_DEN') 242 ENDIF 243C 244C 245C------------------------------------------------------ 246C Initialize remaining one electron density arrays. 247C------------------------------------------------------ 248C 249 CALL DZERO(WORK(KONEAB),NMATAB(1)) 250 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 251 CALL DZERO(WORK(KONEIA),NT1AMX) 252C 253C-------------------------------------------------------- 254C Construct remaining blocks of one electron density. 255C-------------------------------------------------------- 256C 257 CALL DZERO(WORK(KINTIJ),NMATIJ(1)) 258 CALL DIJGEN(WORK(KONEIJ),WORK(KINTIJ)) 259 CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI)) 260C 261 IF (LOCDBG) THEN 262 DIJNO = DDOT(NMATIJ(1),WORK(KONEIJ),1,WORK(KONEIJ),1) 263 DAINO = DDOT(NT1AMX,WORK(KONEAI),1,WORK(KONEAI),1) 264 DIANO = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 265 DABNO = DDOT(NMATAB(1),WORK(KONEAB),1,WORK(KONEAB),1) 266 WRITE(LUPRI,*) ' ' 267 WRITE(LUPRI,*) 'Norm of one electron densities in MO-basis:' 268 WRITE(LUPRI,*) DIJNO, DAINO, DIANO, DABNO 269 ENDIF 270C 271C--------------------------------- 272C Read one-electron integrals. 273C--------------------------------- 274C 275 CALL CCRHS_ONEAO(WORK(KONINT),WORK(KEND1),LWRK1) 276C 277C-------------------------------------------------- 278C Transform one electron integrals to MO-basis. 279C-------------------------------------------------- 280C 281 ISYM = 1 282 CALL CCDINTMO(WORK(KINTIJ),WORK(KINTIA),WORK(KINTAB), 283 * WORK(KINTAI),WORK(KONINT),WORK(KLAMDP), 284 * WORK(KLAMDH),WORK(KEND1),LWRK1,ISYM) 285C 286 IF (LOCDBG) THEN 287C 288 HIJNO = DDOT(NMATIJ(1),WORK(KINTIJ),1,WORK(KINTIJ),1) 289 HAINO = DDOT(NT1AMX,WORK(KINTAI),1,WORK(KINTAI),1) 290 HIANO = DDOT(NT1AMX,WORK(KINTIA),1,WORK(KINTIA),1) 291 HABNO = DDOT(NMATAB(1),WORK(KINTAB),1,WORK(KINTAB),1) 292 WRITE(LUPRI,*) ' ' 293 WRITE(LUPRI,*) 'Norm of one electron integrals in MO-basis:' 294 WRITE(LUPRI,*) HIJNO, HAINO, HIANO, HABNO 295C 296 ENDIF 297C 298C----------------------------------------------- 299C Set up transposed integrals and densities. 300C----------------------------------------------- 301C 302 CALL DCOPY(NMATAB(1),WORK(KINTAB),1,WORK(KINABT),1) 303 CALL DCOPY(NMATIJ(1),WORK(KINTIJ),1,WORK(KINIJT),1) 304 CALL DCOPY(NMATAB(1),WORK(KONEAB),1,WORK(KD1ABT),1) 305 CALL DCOPY(NMATIJ(1),WORK(KONEIJ),1,WORK(KD1IJT),1) 306C 307 CALL CC_EITR(WORK(KINABT),WORK(KINIJT),WORK(KEND1),LWRK1,1) 308 CALL CC_EITR(WORK(KD1ABT),WORK(KD1IJT),WORK(KEND1),LWRK1,1) 309C 310 IF (IOPT .EQ. 1) THEN 311C 312C--------------------------------------------------------------------- 313C Calculate direct one electron contribution to eta-ij and eta-ab. 314C--------------------------------------------------------------------- 315C 316 CALL DZERO(WORK(KETAIJ),NMATIJ(1)) 317 CALL DZERO(WORK(KETAAB),NMATAB(1)) 318C 319 ISYM = 1 320 CALL CC2_ETIJ(WORK(KETAIJ),WORK(KINTIJ),WORK(KINTAI), 321 * WORK(KINTIA),WORK(KINTAB),WORK(KONEIJ), 322 * WORK(KONEAI),WORK(KONEIA),WORK(KONEAB), 323 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 324C 325 CALL CC2_ETAB(WORK(KETAAB),WORK(KINTIJ),WORK(KINTAI), 326 * WORK(KINTIA),WORK(KINTAB),WORK(KONEIJ), 327 * WORK(KONEAI),WORK(KONEIA),WORK(KONEAB), 328 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 329C 330 ELSE IF (IOPT .EQ. 2) THEN 331C 332C---------------------------------------------------------- 333C Calculate direct one electron contribution to eta-ai. 334C---------------------------------------------------------- 335C 336 ISYM = 1 337 CALL CCDENZK0(ETAAI,WORK(KINTIJ),WORK(KINTAI),WORK(KINTIA), 338 * WORK(KINTAB),WORK(KONEIJ),WORK(KONEAI), 339 * WORK(KONEIA),WORK(KONEAB),WORK(KINIJT), 340 * WORK(KINABT),WORK(KD1IJT),WORK(KD1ABT), 341 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 342C 343 ENDIF 344C 345C TIMONE = SECOND() - TIMONE 346C 347 KEND1 = KONINT 348 LWRK1 = LWORK - KEND1 349C 350C--------------------------------------------------------------- 351C Start the loop for the 2-electron integrals and densities. 352C--------------------------------------------------------------- 353C 354 KENDS2 = KEND1 355 LWRKS2 = LWRK1 356C 357 IF (DIRECT) THEN 358 IF (HERDIR) THEN 359 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 360 ELSE 361 KCCFB1 = KEND1 362 KINDXB = KCCFB1 + MXPRIM*MXCONT 363 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 364 LWRK1 = LWORK - KEND1 365 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 366 * KODPP1,KODPP2,KRDPP1,KRDPP2, 367 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 368 * WORK(KEND1),LWRK1,IPRERI) 369 KEND1 = KFREE 370 LWRK1 = LFREE 371 ENDIF 372 NTOSYM = 1 373 ELSE 374 NTOSYM = NSYM 375 ENDIF 376C 377 KENDSV = KEND1 378 LWRKSV = LWRK1 379C 380 ICDEL1 = 0 381 DO 100 ISYMD1 = 1,NTOSYM 382C 383 IF (DIRECT) THEN 384 IF (HERDIR) THEN 385 NTOT = MAXSHL 386 ELSE 387 NTOT = MXCALL 388 ENDIF 389 ELSE 390 NTOT = NBAS(ISYMD1) 391 ENDIF 392C 393 DO 110 ILLL = 1,NTOT 394C 395C--------------------------------------------- 396C If direct calculate the integrals. 397C--------------------------------------------- 398C 399 IF (DIRECT) THEN 400C 401 KEND1 = KENDSV 402 LWRK1 = LWRKSV 403C 404 DTIME = SECOND() 405 IF (HERDIR) THEN 406 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 407 & IPRINT) 408 ELSE 409 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 410 * WORK(KODCL1),WORK(KODCL2), 411 * WORK(KODBC1),WORK(KODBC2), 412 * WORK(KRDBC1),WORK(KRDBC2), 413 * WORK(KODPP1),WORK(KODPP2), 414 * WORK(KRDPP1),WORK(KRDPP2), 415 * WORK(KCCFB1),WORK(KINDXB), 416 * WORK(KEND1), LWRK1,IPRERI) 417 ENDIF 418 DTIME = SECOND() - DTIME 419 TIMHE2 = TIMHE2 + DTIME 420C 421 KRECNR = KEND1 422 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 423 LWRK1 = LWORK - KEND1 424 IF (LWRK1 .LT. 0) THEN 425 CALL QUIT('Insufficient core in CC2_DEN') 426 END IF 427C 428 ELSE 429 NUMDIS = 1 430 ENDIF 431C 432C----------------------------------------------------- 433C Loop over number of distributions in disk. 434C----------------------------------------------------- 435C 436 DO 120 IDEL2 = 1,NUMDIS 437C 438 IF (DIRECT) THEN 439 IDEL = INDEXA(IDEL2) 440 ISYMD = ISAO(IDEL) 441 IF (NOAUXB) THEN 442 IDUM = 1 443 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 444 END IF 445 ELSE 446 IDEL = IBAS(ISYMD1) + ILLL 447 ISYMD = ISYMD1 448 ENDIF 449C 450C---------------------------------------- 451C Work space allocation two. 452C---------------------------------------- 453C 454 ISYDEN = ISYMD 455C 456 KD2IJG = KEND1 457 KD2AIG = KD2IJG + ND2IJG(ISYDEN) 458 KD2IAG = KD2AIG + ND2AIG(ISYDEN) 459 KD2ABG = KD2IAG + ND2AIG(ISYDEN) 460 KEND2 = KD2ABG + ND2ABG(ISYDEN) 461 LWRK2 = LWORK - KEND2 462C 463 IF (LWRK2 .LT. 0) THEN 464 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 465 CALL QUIT( 466 * 'Insufficient space for allocation 2 in CC2_DEN') 467 ENDIF 468C 469C------------------------------------------------------- 470C Initialize 4 two electron density arrays. 471C------------------------------------------------------- 472C 473 CALL DZERO(WORK(KD2IJG),ND2IJG(ISYDEN)) 474 CALL DZERO(WORK(KD2AIG),ND2AIG(ISYDEN)) 475 CALL DZERO(WORK(KD2IAG),ND2AIG(ISYDEN)) 476 CALL DZERO(WORK(KD2ABG),ND2ABG(ISYDEN)) 477C 478C------------------------------------------------------------------- 479C Calculate the two electron density d(pq,gamma;delta). 480C------------------------------------------------------------------- 481C 482 AUTIME = SECOND() 483C 484 CALL CC_DEN2(WORK(KD2IJG),WORK(KD2AIG),WORK(KD2IAG), 485 * WORK(KD2ABG),WORK(KZ2AM),WORK(KT2AM), 486 * WORK(KT2AMT),WORK(KEND2),WORK(KEND2), 487 * WORK(KEND2),WORK(KONEAB),WORK(KONEAI), 488 * WORK(KONEIA),WORK(KEND2),WORK(KLAMDH),1, 489 * WORK(KLAMDP),1,WORK(KEND2),LWRK2,IDEL,ISYMD) 490C 491 AUTIME = SECOND() - AUTIME 492 TIMDEN = TIMDEN + AUTIME 493C 494C------------------------------------------ 495C Work space allocation three. 496C------------------------------------------ 497C 498 ISYDIS = MULD2H(ISYMD,ISYMOP) 499C 500 KXINT = KEND2 501 KEND3 = KXINT + NDISAO(ISYDIS) 502 LWRK3 = LWORK - KEND3 503C 504 IF (LWRK3 .LT. 0) THEN 505 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 506 CALL QUIT( 507 * 'Insufficient space for allocation 3 in CC2_DEN') 508 ENDIF 509C 510C-------------------------------------------- 511C Read AO integral distribution. 512C-------------------------------------------- 513C 514 AUTIME = SECOND() 515 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 516 * WORK(KRECNR),DIRECT) 517 AUTIME = SECOND() - AUTIME 518 TIRDAO = TIRDAO + AUTIME 519C 520 IF (IOPT .EQ. 2) THEN 521C 522C-------------------------------------------------------------- 523C Calculate correction terms to eta_ai originating 524C from kbar_ij and kbar_ab. 525C-------------------------------------------------------------- 526C 527 KDSRHF = KEND3 528 K3OINT = KDSRHF + NDSRHF(ISYMD) 529 KEND3 = K3OINT + NMAIJK(ISYDIS) 530 LWRK3 = LWORK - KEND3 531C 532 IF (LWRK3 .LT. 0) THEN 533 WRITE(LUPRI,*) 'Need : ',KEND3,'Available : ',LWORK 534 CALL QUIT( 535 * 'Insufficient core for integrals in CC2_DEN') 536 ENDIF 537C 538C-------------------------------------------------------------------- 539C Transform one index in the integral batch to occupied. 540C-------------------------------------------------------------------- 541C 542 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KCMO), 543 * ISYMOP,WORK(KEND3),LWRK3,ISYDIS) 544C 545C------------------------------------------------------------------- 546C Calculate integral batch with three occupied indices. 547C Since LUDUMLOCAL < 0 the integrals are not written 548C to disk. 549C------------------------------------------------------------------- 550C 551 LUDUMLOCAL = -10 552 CALL CC_INT3O(WORK(K3OINT),WORK(KDSRHF),WORK(KCMO), 553 * ISYMOP,WORK(KCMO),WORK(KEND3),LWRK3, 554 * IDEL,ISYMD,LUDUMLOCAL,'DUMMY') 555C 556C---------------------------------------------------- 557C Calculate the actual correction terms. 558C---------------------------------------------------- 559C 560 CALL MP2_YTV(ETAAI,ZKDIA(NMATIJ(1)+1),WORK(KDSRHF), 561 * WORK(KCMO),WORK(KEND3),LWRK3,IDEL,ISYMD) 562C 563 CALL MP2_NXY(ETAAI,ZKDIA(1),ZKDIA(NMATIJ(1)+1), 564 * WORK(K3OINT),WORK(KDSRHF),WORK(KCMO), 565 * WORK(KEND3),LWRK3,IDEL,ISYMD) 566C 567 CALL MP2_XTO(ETAAI,ZKDIA(1),WORK(K3OINT), 568 * WORK(KCMO),WORK(KEND3),LWRK3,IDEL,ISYMD) 569C 570 ENDIF 571C 572C------------------------------------------------------ 573C Start loop over second AO-index (gamma). 574C------------------------------------------------------ 575C 576 AUTIME = SECOND() 577C 578 DO 130 ISYMG = 1,NSYM 579C 580 ISYMPQ = MULD2H(ISYMG,ISYDEN) 581C 582 DO 140 G = 1,NBAS(ISYMG) 583C 584 KD2GIJ = KD2IJG + ID2IJG(ISYMPQ,ISYMG) 585 * + NMATIJ(ISYMPQ)*(G - 1) 586 KD2GAI = KD2AIG + ID2AIG(ISYMPQ,ISYMG) 587 * + NT1AM(ISYMPQ)*(G - 1) 588 KD2GAB = KD2ABG + ID2ABG(ISYMPQ,ISYMG) 589 * + NMATAB(ISYMPQ)*(G - 1) 590 KD2GIA = KD2IAG + ID2AIG(ISYMPQ,ISYMG) 591 * + NT1AM(ISYMPQ)*(G - 1) 592C 593C----------------------------------------------- 594C Work space allocation four. 595C----------------------------------------------- 596C 597 KINTAO = KEND3 598 KD2AOB = KINTAO + N2BST(ISYMPQ) 599 KEND4 = KD2AOB + N2BST(ISYMPQ) 600 LWRK4 = LWORK - KEND4 601C 602 IF (LWRK4 .LT. 0) THEN 603 WRITE(LUPRI,*) 'Available:', LWORK 604 WRITE(LUPRI,*) 'Needed:', KEND4 605 CALL QUIT('Insufficient work space in CC2_DEN') 606 ENDIF 607C 608 CALL DZERO(WORK(KD2AOB),N2BST(ISYMPQ)) 609C 610C------------------------------------------------------- 611C Square up AO-integral distribution. 612C------------------------------------------------------- 613C 614 KOFFIN = KXINT + IDSAOG(ISYMG,ISYDIS) 615 * + NNBST(ISYMPQ)*(G - 1) 616C 617 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ, 618 * WORK(KINTAO)) 619C 620C----------------------------------------------- 621C Work space allocation five. 622C----------------------------------------------- 623C 624 KIJINT = KEND4 625 KAIINT = KIJINT + NMATIJ(ISYMPQ) 626 KIAINT = KAIINT + NT1AM(ISYMPQ) 627 KABINT = KIAINT + NT1AM(ISYMPQ) 628 KABTIN = KABINT + NMATAB(ISYMPQ) 629 KIJTIN = KABTIN + NMATAB(ISYMPQ) 630 KD2TAB = KIJTIN + NMATIJ(ISYMPQ) 631 KD2TIJ = KD2TAB + NMATAB(ISYMPQ) 632 KEND5 = KD2TIJ + NMATIJ(ISYMPQ) 633 LWRK5 = LWORK - KEND5 634C 635 IF (LWRK5 .LT. 0) THEN 636 WRITE(LUPRI,*) 'Available:', LWORK 637 WRITE(LUPRI,*) 'Needed:', KEND5 638 CALL QUIT('Insufficient work space in CC2_DEN') 639 ENDIF 640C 641C------------------------------------------------------------- 642C Transform 2 integral indices to MO-basis. 643C------------------------------------------------------------- 644C 645 ISYM = ISYMPQ 646 CALL CCDINTMO(WORK(KIJINT),WORK(KIAINT), 647 * WORK(KABINT),WORK(KAIINT), 648 * WORK(KINTAO),WORK(KLAMDP), 649 * WORK(KLAMDH),WORK(KEND5), 650 * LWRK5,ISYM) 651C 652C-------------------------------------------------------------- 653C Set up transposed integrals and densities. 654C-------------------------------------------------------------- 655C 656 CALL DCOPY(NMATAB(ISYMPQ),WORK(KABINT),1, 657 * WORK(KABTIN),1) 658 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KIJINT),1, 659 * WORK(KIJTIN),1) 660 CALL DCOPY(NMATAB(ISYMPQ),WORK(KD2GAB),1, 661 * WORK(KD2TAB),1) 662 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KD2GIJ),1, 663 * WORK(KD2TIJ),1) 664C 665 CALL CC_EITR(WORK(KABTIN),WORK(KIJTIN), 666 * WORK(KEND5),LWRK5,ISYMPQ) 667 CALL CC_EITR(WORK(KD2TAB),WORK(KD2TIJ), 668 * WORK(KEND5),LWRK5,ISYMPQ) 669C 670 IF (IOPT .EQ. 1) THEN 671C 672 ISYM = ISYMPQ 673 CALL CC2_ETIJ(WORK(KETAIJ),WORK(KIJINT), 674 * WORK(KAIINT),WORK(KIAINT), 675 * WORK(KABINT),WORK(KD2GIJ), 676 * WORK(KD2GAI),WORK(KD2GIA), 677 * WORK(KD2GAB),WORK(KT1AM), 678 * WORK(KEND5),LWRK5,ISYM) 679C 680 CALL CC2_ETAB(WORK(KETAAB),WORK(KIJINT), 681 * WORK(KAIINT),WORK(KIAINT), 682 * WORK(KABINT),WORK(KD2GIJ), 683 * WORK(KD2GAI),WORK(KD2GIA), 684 * WORK(KD2GAB),WORK(KT1AM), 685 * WORK(KEND5),LWRK5,ISYM) 686C 687 ELSE IF (IOPT .EQ. 2) THEN 688C 689C----------------------------------------------------------------- 690C Calculate direct 2 e- contribution to eta_ai. 691C----------------------------------------------------------------- 692C 693 ISYM = ISYMPQ 694 CALL CCDENZK0(ETAAI,WORK(KIJINT),WORK(KAIINT), 695 * WORK(KIAINT),WORK(KABINT), 696 * WORK(KD2GIJ),WORK(KD2GAI), 697 * WORK(KD2GIA),WORK(KD2GAB), 698 * WORK(KIJTIN),WORK(KABTIN), 699 * WORK(KD2TIJ),WORK(KD2TAB), 700 * WORK(KT1AM),WORK(KEND5),LWRK5, 701 * ISYM) 702C 703 ENDIF 704C 705 140 CONTINUE 706 130 CONTINUE 707C 708 AUTIME = SECOND() - AUTIME 709 TIMRES = TIMRES + AUTIME 710C 711 120 CONTINUE 712 110 CONTINUE 713 100 CONTINUE 714C 715C----------------------- 716C Regain work space. 717C----------------------- 718C 719 KEND1 = KENDS2 720 LWRK1 = LWRKS2 721C 722C-------------------------------------------------- 723C Calculate the diagonal blocks of kappa-bar-0. 724C-------------------------------------------------- 725C 726c IF ((FROIMP) .AND. (IOPT .EQ. 1)) THEN 727c CALL CCSD_ZKBLO(ZKDIA,WORK(KEND1),LWRK1) 728c ENDIF 729C 730C----------------------------------------------------------- 731C Calculate the correlated-frozen blocks of kappa-bar-0. 732C----------------------------------------------------------- 733C 734c IF (FROIMP) THEN 735c KOFFIJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 736c CALL CC_ZKFCB(ZKDIA(KOFFIJ),WORK(KEND1),LWRK1) 737C 738C---------------------------------------------------------------- 739C Calculate correction terms from correlated-frozen block. 740C---------------------------------------------------------------- 741C 742c KOFFAI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 743c CALL CC_ETACOR(ETAAI,ZKDIA(KOFFAI),ZKDIA(KOFFIJ), 744c * WORK(KCMOF),WORK(KEND1),LWRK1) 745C 746C--------------------- 747C Reorder results. 748C--------------------- 749C 750c CALL CC_ETARE(ETAAI,ZKDIA(KOFFAI),WORK(KEND1),LWRK1) 751c ENDIF 752C 753C--------------------------------- 754C Write out eta-ai and eta-aI. 755C--------------------------------- 756C 757 IF (LOCDBG) THEN 758C 759 CALL AROUND('Eta-kappa-bar-0-ai vector exiting CC2_DEN') 760C 761 DO 20 ISYM = 1,NSYM 762C 763 WRITE(LUPRI,*) ' ' 764 WRITE(LUPRI,444) 'Sub-symmetry block number:', ISYM 765 WRITE(LUPRI,555) '--------------------------' 766 444 FORMAT(3X,A26,2X,I1) 767 555 FORMAT(3X,A25) 768C 769 KOFF = IALLAI(ISYM,ISYM) + 1 770 CALL OUTPUT(ETAAI(KOFF),1,NVIR(ISYM),1,NRHFS(ISYM), 771 * NVIR(ISYM),NRHFS(ISYM),1,LUPRI) 772C 773 IF ((NVIR(ISYM) .EQ. 0) .OR. (NRHFS(ISYM) .EQ. 0)) THEN 774 WRITE(LUPRI,*) 'This sub-symmetry is empty' 775 ENDIF 776C 777 20 CONTINUE 778 ENDIF 779C 780 IF (IPRINT .GT. 5) THEN 781 ETAKAN = DDOT(NALLAI(1),ETAAI,1,ETAAI,1) 782 WRITE(LUPRI,*) ' ' 783 WRITE(LUPRI,*) 'Norm of occupied-virtual block:', ETAKAN 784 ENDIF 785C 786C------------------------------------------- 787C Write out frozen block of kappa-bar-0. 788C------------------------------------------- 789C 790c IF ((IPRINT .GT. 9) .AND. (FROIMP) .AND. (IOPT .EQ. 1)) THEN 791c CALL AROUND('Zeta-kappa-0 correlated-frozen block') 792c KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 793c ZKAPF1 = DDOT(NCOFRO(1),ZKDIA(KOFRES),1, 794c * ZKDIA(KOFRES),1) 795c ZKAPFR = ZKAPF1**0.5 796c WRITE(LUPRI,*) 'Norm of correlated-frozen block:', ZKAPFR 797c ENDIF 798C 799C----------------------- 800C Write out timings. 801C----------------------- 802C 803 99 TIMETO = SECOND() - TIMETO 804C 805 IF (IPRINT .GT. 3) THEN 806 WRITE(LUPRI,*) ' ' 807 WRITE(LUPRI,*) 808 * 'Requested density dependent quantities calculated' 809 WRITE(LUPRI,*) 'Total time used in CC2_DEN:', TIMETO 810 ENDIF 811 IF (IPRINT .GT. 5) THEN 812 WRITE(LUPRI,'(A,F15.2)') 813 & ' Time used for setting up 2 e- density :', TIMDEN 814 & ,' Time used for contraction with integrals :', TIMRES 815 & ,' Time used for reading 2 e- AO-integrals :', TIRDAO 816 & ,' Time used for calculating 2 e- AO-integrals:', TIMHE2 817C & ,' Time used for 1 e- density & intermediates :', TIMONE 818 ENDIF 819C 820C---------------------------------------------- 821C Calculate the ZK0(ab) and ZK0(ij) blocks. 822C---------------------------------------------- 823C 824 IF (IOPT .EQ. 1) THEN 825 CALL MP2_ZKBLO(ZKDIA,WORK(KETAIJ),WORK(KETAAB),WORK(KEND1), 826 * LWRK1) 827 ENDIF 828C 829C--------------------------------------------------- 830C Calculate frozen core occupied blocks ZK0(iJ). 831C--------------------------------------------------- 832C 833 IF (FROIMP) THEN 834 CALL QUIT('NO FROZEN CORE FOR RELAXED CC2 YET') 835c KOFRES = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 2*NT1FRO(1) + 1 836c CALL MP2_ZKFCB(ZKDIA(KOFRES),WORK(KZ2AM),WORK(KEND1),LWRK1) 837 ENDIF 838C 839C------------------------------------------------ 840C Write out timings and results if requested. 841C------------------------------------------------ 842C 843 IF ((IPRINT .GT. 3) .AND. (IOPT .EQ. 1)) THEN 844C 845 CALL AROUND('eta_ij and eta_ab blocks') 846 ETAIJN = DDOT(NMATIJ(1),WORK(KETAIJ),1,WORK(KETAIJ),1) 847 ETAABN = DDOT(NMATAB(1),WORK(KETAAB),1,WORK(KETAAB),1) 848 ETAIJN = ETAIJN**0.5 849 ETAABN = ETAABN**0.5 850 WRITE(LUPRI,*) ' ' 851 WRITE(LUPRI,*) 'Norm of occupied-occupied block:', ETAIJN 852 WRITE(LUPRI,*) 'Norm of virtual-virtual block :', ETAABN 853C 854 CALL AROUND('Zeta-kappa-0 diagonal blocks') 855 ZKAPI1 = DDOT(NMATIJ(1),ZKDIA(1),1,ZKDIA(1),1) 856 ZKAPA1 = DDOT(NMATAB(1),ZKDIA(NMATIJ(1)+1),1, 857 * ZKDIA(NMATIJ(1)+1),1) 858 ZKAPIJ = ZKAPI1**0.5 859 ZKAPAB = ZKAPA1**0.5 860 WRITE(LUPRI,*) ' ' 861 WRITE(LUPRI,*) 'Norm of occupied-occupied block:', ZKAPIJ 862 WRITE(LUPRI,*) 'Norm of virtual-virtual block :', ZKAPAB 863C 864c IF (FROIMP) THEN 865c ZKAPF1 = DDOT(NCOFRO(1),ZKDIA(KOFRES),1, 866c * ZKDIA(KOFRES),1) 867c ZKAPFR = ZKAPF1**0.5 868c WRITE(LUPRI,*) 'Norm of frozen-core-occupied block:', ZKAPFR 869c ENDIF 870C 871 IF (IPRINT .GT. 50) THEN 872 DO 240 ISYM = 1,NSYM 873 WRITE(LUPRI,*) ' ' 874 WRITE(LUPRI,*) 'Symmetry block:', ISYM 875 KIJ = IMATIJ(ISYM,ISYM) + 1 876 KAB = IMATAB(ISYM,ISYM) + 1 + NMATIJ(1) 877 CALL AROUND('occ-occ block') 878 CALL OUTPUT(ZKDIA(KIJ),1,NRHF(ISYM),1,NRHF(ISYM), 879 * NRHF(ISYM),NRHF(ISYM),1,LUPRI) 880 CALL AROUND('vir-vir block') 881 CALL OUTPUT(ZKDIA(KAB),1,NVIR(ISYM),1,NVIR(ISYM), 882 * NVIR(ISYM),NVIR(ISYM),1,LUPRI) 883 240 CONTINUE 884 ENDIF 885 ENDIF 886C 887 TIMETO = SECOND() - TIMETO 888C 889 IF (IPRINT .GT. 3) THEN 890 WRITE(LUPRI,*) ' ' 891 WRITE(LUPRI,*) 'Diagonal blocks of Zeta-kappa-0 calculated' 892 WRITE(LUPRI,*) 'Total time used in CC2_DEN:', TIMETO 893 ENDIF 894C 895 CALL QEXIT('CC2_DEN') 896C 897 RETURN 898 END 899C /* Deck cc2_etij */ 900 SUBROUTINE CC2_ETIJ(ETAIJ,XINTIJ,XINTAI,XINTIA,XINTAB,DIJ,DAI, 901 * DIA,DAB,T1AM,WORK,LWORK,ISYM) 902C 903C Written by A. Halkier & S. Coriani 14/1-2000 904C 905C Version: 1.0 906C 907C Purpose: To set up the right hand side of the equation for 908C zeta-kappa-0_ij (ETAIJ) from MO-integrals (XIN*), CC2 909C densities (D*) and t1-amplitudes (T1AM)! 910C Note that due to the symmetry in the formulas, this 911C routine is able to handle both the one- and the two- 912C electron contributions! 913C ISYM is the symmetry of both the density and the 914C integrals! 915C 916#include "implicit.h" 917 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 918 DIMENSION ETAIJ(*), XINTIJ(*), XINTAI(*), XINTIA(*), XINTAB(*) 919 DIMENSION DIJ(*), DAI(*), DIA(*), DAB(*), T1AM(*), WORK(LWORK) 920#include "priunit.h" 921#include "ccorb.h" 922#include "ccsdsym.h" 923#include "cclr.h" 924C 925 CALL QENTER('CC2_ETIJ') 926C 927 DO 100 ISYMI = 1,NSYM 928C 929C---------------------------------------------------------------- 930C Calculate direct terms - those without T1AM - to eta_ij. 931C---------------------------------------------------------------- 932C 933 ISYMJ = ISYMI 934 ISYMK = MULD2H(ISYMI,ISYM) 935 ISYMC = MULD2H(ISYMI,ISYM) 936C 937 KOFFRE = IMATIJ(ISYMI,ISYMJ) + 1 938C 939 NTOTRE = MAX(NRHF(ISYMI),1) 940 NTOTI = MAX(NRHF(ISYMI),1) 941 NTOTJ = MAX(NRHF(ISYMJ),1) 942 NTOTK = MAX(NRHF(ISYMK),1) 943 NTOTC = MAX(NVIR(ISYMC),1) 944C 945 KOFF1 = IMATIJ(ISYMK,ISYMI) + 1 946 KOFF2 = IMATIJ(ISYMK,ISYMJ) + 1 947 KOFF3 = IMATIJ(ISYMI,ISYMK) + 1 948 KOFF4 = IMATIJ(ISYMJ,ISYMK) + 1 949C 950 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),ONE, 951 * XINTIJ(KOFF1),NTOTK,DIJ(KOFF2),NTOTK,ONE, 952 * ETAIJ(KOFFRE),NTOTRE) 953C 954 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),ONE, 955 * XINTIJ(KOFF3),NTOTI,DIJ(KOFF4),NTOTJ,ONE, 956 * ETAIJ(KOFFRE),NTOTRE) 957C 958 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),-ONE, 959 * DIJ(KOFF3),NTOTI,XINTIJ(KOFF4),NTOTJ,ONE, 960 * ETAIJ(KOFFRE),NTOTRE) 961C 962 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),-ONE, 963 * DIJ(KOFF1),NTOTK,XINTIJ(KOFF2),NTOTK,ONE, 964 * ETAIJ(KOFFRE),NTOTRE) 965C 966 KOFF5 = IT1AM(ISYMC,ISYMI) + 1 967 KOFF6 = IT1AM(ISYMC,ISYMJ) + 1 968C 969 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 970 * XINTAI(KOFF5),NTOTC,DAI(KOFF6),NTOTC,ONE, 971 * ETAIJ(KOFFRE),NTOTRE) 972C 973 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 974 * XINTIA(KOFF5),NTOTC,DIA(KOFF6),NTOTC,ONE, 975 * ETAIJ(KOFFRE),NTOTRE) 976C 977 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 978 * DIA(KOFF5),NTOTC,XINTIA(KOFF6),NTOTC,ONE, 979 * ETAIJ(KOFFRE),NTOTRE) 980C 981 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 982 * DAI(KOFF5),NTOTC,XINTAI(KOFF6),NTOTC,ONE, 983 * ETAIJ(KOFFRE),NTOTRE) 984C 985 100 CONTINUE 986C 987 DO 110 ISYMI = 1,NSYM 988C 989C--------------------------------------------------------------------------- 990C Calculate terms involving T1AM and intermediate looping over occupied. 991C Note that the intermediate of integrals and densities is used for both 992C contributions as ISYMI=ISYMJ. 993C--------------------------------------------------------------------------- 994C 995 ISYMJ = ISYMI 996 ISYMB = ISYMJ 997 ISYMK = MULD2H(ISYMB,ISYM) 998C 999 KOFFRE = IMATIJ(ISYMI,ISYMJ) + 1 1000C 1001 NTOTRE = MAX(NRHF(ISYMI),1) 1002 NTOTI = MAX(NRHF(ISYMI),1) 1003 NTOTJ = MAX(NRHF(ISYMJ),1) 1004 NTOTK = MAX(NRHF(ISYMK),1) 1005 NTOTB = MAX(NVIR(ISYMB),1) 1006C 1007C---------------------------------- 1008C Work space allocation one. 1009C---------------------------------- 1010C 1011 KSCRHD = 1 1012 KEND1 = KSCRHD + NVIR(ISYMB)*NRHF(ISYMI) 1013 LWRK1 = LWORK - KEND1 1014C 1015 IF (LWRK1 .LT. 0) THEN 1016 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1017 CALL QUIT( 1018 * 'Insufficient memory for allocation (1) in CC2_ETIJ') 1019 ENDIF 1020C 1021 CALL DZERO(WORK(KSCRHD),NVIR(ISYMB)*NRHF(ISYMI)) 1022C 1023 KOFF1 = IT1AM(ISYMB,ISYMK) + 1 1024 KOFF2 = IMATIJ(ISYMK,ISYMI) + 1 1025 KOFF3 = IMATIJ(ISYMI,ISYMK) + 1 1026 KOFF4 = IT1AM(ISYMB,ISYMJ) + 1 1027C 1028 CALL DGEMM('N','N',NVIR(ISYMB),NRHF(ISYMI),NRHF(ISYMK),ONE, 1029 * XINTIA(KOFF1),NTOTB,DIJ(KOFF2),NTOTK,ZERO, 1030 * WORK(KSCRHD),NTOTB) 1031C 1032 CALL DGEMM('N','T',NVIR(ISYMB),NRHF(ISYMI),NRHF(ISYMK),-ONE, 1033 * DAI(KOFF1),NTOTB,XINTIJ(KOFF3),NTOTI,ONE, 1034 * WORK(KSCRHD),NTOTB) 1035C 1036 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMB),ONE, 1037 * WORK(KSCRHD),NTOTB,T1AM(KOFF4),NTOTB,ONE, 1038 * ETAIJ(KOFFRE),NTOTRE) 1039C 1040 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMB),-ONE, 1041 * T1AM(KOFF4),NTOTB,WORK(KSCRHD),NTOTB,ONE, 1042 * ETAIJ(KOFFRE),NTOTRE) 1043C 1044 110 CONTINUE 1045C 1046 DO 120 ISYMI = 1,NSYM 1047C 1048C--------------------------------------------------------------------------- 1049C Calculate terms involving T1AM and intermediate looping over virtual. 1050C Note that the intermediate of integrals and densities is used for both 1051C contributions as ISYMI=ISYMJ. 1052C--------------------------------------------------------------------------- 1053C 1054 ISYMJ = ISYMI 1055 ISYMB = ISYMJ 1056 ISYMC = MULD2H(ISYMB,ISYM) 1057C 1058 KOFFRE = IMATIJ(ISYMI,ISYMJ) + 1 1059C 1060 NTOTRE = MAX(NRHF(ISYMI),1) 1061 NTOTI = MAX(NRHF(ISYMI),1) 1062 NTOTJ = MAX(NRHF(ISYMJ),1) 1063 NTOTC = MAX(NVIR(ISYMC),1) 1064 NTOTB = MAX(NVIR(ISYMB),1) 1065C 1066C---------------------------------- 1067C Work space allocation two. 1068C---------------------------------- 1069C 1070 KSCRHD = 1 1071 KEND1 = KSCRHD + NVIR(ISYMB)*NRHF(ISYMI) 1072 LWRK1 = LWORK - KEND1 1073C 1074 IF (LWRK1 .LT. 0) THEN 1075 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1076 CALL QUIT( 1077 * 'Insufficient memory for allocation (2) in CC2_ETIJ') 1078 ENDIF 1079C 1080 CALL DZERO(WORK(KSCRHD),NVIR(ISYMB)*NRHF(ISYMI)) 1081C 1082 KOFF1 = IMATAB(ISYMC,ISYMB) + 1 1083 KOFF2 = IT1AM(ISYMC,ISYMI) + 1 1084 KOFF3 = IMATAB(ISYMB,ISYMC) + 1 1085 KOFF4 = IT1AM(ISYMB,ISYMJ) + 1 1086C 1087 CALL DGEMM('T','N',NVIR(ISYMB),NRHF(ISYMI),NVIR(ISYMC),ONE, 1088 * XINTAB(KOFF1),NTOTC,DAI(KOFF2),NTOTC,ZERO, 1089 * WORK(KSCRHD),NTOTB) 1090C 1091 CALL DGEMM('N','N',NVIR(ISYMB),NRHF(ISYMI),NVIR(ISYMC),-ONE, 1092 * DAB(KOFF3),NTOTB,XINTIA(KOFF2),NTOTC,ONE, 1093 * WORK(KSCRHD),NTOTB) 1094C 1095 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMB),ONE, 1096 * WORK(KSCRHD),NTOTB,T1AM(KOFF4),NTOTB,ONE, 1097 * ETAIJ(KOFFRE),NTOTRE) 1098C 1099 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMB),-ONE, 1100 * T1AM(KOFF4),NTOTB,WORK(KSCRHD),NTOTB,ONE, 1101 * ETAIJ(KOFFRE),NTOTRE) 1102C 1103 120 CONTINUE 1104C 1105 CALL QEXIT('CC2_ETIJ') 1106C 1107 RETURN 1108 END 1109C /* Deck cc2_etab */ 1110 SUBROUTINE CC2_ETAB(ETAAB,XINTIJ,XINTAI,XINTIA,XINTAB,DIJ,DAI, 1111 * DIA,DAB,T1AM,WORK,LWORK,ISYM) 1112C 1113C Written by A. Halkier & S. Coriani 15/1-2000 1114C 1115C Version: 1.0 1116C 1117C Purpose: To set up the right hand side of the equation for 1118C zeta-kappa-0_ab (ETAAB) from MO-integrals (XIN*), CC2 1119C densities (D*) and t1-amplitudes (T1AM)! 1120C Note that due to the symmetry in the formulas, this 1121C routine is able to handle both the one- and the two- 1122C electron contributions! 1123C ISYM is the symmetry of both the density and the 1124C integrals! 1125C 1126#include "implicit.h" 1127 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1128 DIMENSION ETAAB(*), XINTIJ(*), XINTAI(*), XINTIA(*), XINTAB(*) 1129 DIMENSION DIJ(*), DAI(*), DIA(*), DAB(*), T1AM(*), WORK(LWORK) 1130#include "priunit.h" 1131#include "ccorb.h" 1132#include "ccsdsym.h" 1133#include "cclr.h" 1134C 1135 CALL QENTER('CC2_ETAB') 1136C 1137 DO 100 ISYMA = 1,NSYM 1138C 1139C---------------------------------------------------------------- 1140C Calculate direct terms - those without T1AM - to eta_ab. 1141C---------------------------------------------------------------- 1142C 1143 ISYMB = ISYMA 1144 ISYMK = MULD2H(ISYMA,ISYM) 1145 ISYMC = MULD2H(ISYMA,ISYM) 1146C 1147 KOFFRE = IMATAB(ISYMA,ISYMB) + 1 1148C 1149 NTOTRE = MAX(NVIR(ISYMA),1) 1150 NTOTA = MAX(NVIR(ISYMA),1) 1151 NTOTB = MAX(NVIR(ISYMB),1) 1152 NTOTC = MAX(NVIR(ISYMC),1) 1153 NTOTK = MAX(NRHF(ISYMK),1) 1154C 1155 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1156 KOFF2 = IT1AM(ISYMB,ISYMK) + 1 1157C 1158 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 1159 * XINTIA(KOFF1),NTOTA,DIA(KOFF2),NTOTB,ONE, 1160 * ETAAB(KOFFRE),NTOTRE) 1161C 1162 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 1163 * DAI(KOFF1),NTOTA,XINTAI(KOFF2),NTOTB,ONE, 1164 * ETAAB(KOFFRE),NTOTRE) 1165C 1166 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 1167 * DIA(KOFF1),NTOTA,XINTIA(KOFF2),NTOTB,ONE, 1168 * ETAAB(KOFFRE),NTOTRE) 1169C 1170 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 1171 * XINTAI(KOFF1),NTOTA,DAI(KOFF2),NTOTB,ONE, 1172 * ETAAB(KOFFRE),NTOTRE) 1173C 1174 KOFF3 = IMATAB(ISYMC,ISYMA) + 1 1175 KOFF4 = IMATAB(ISYMC,ISYMB) + 1 1176 KOFF5 = IMATAB(ISYMA,ISYMC) + 1 1177 KOFF6 = IMATAB(ISYMB,ISYMC) + 1 1178C 1179 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),ONE, 1180 * XINTAB(KOFF3),NTOTC,DAB(KOFF4),NTOTC,ONE, 1181 * ETAAB(KOFFRE),NTOTRE) 1182C 1183 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),-ONE, 1184 * DAB(KOFF5),NTOTA,XINTAB(KOFF6),NTOTB,ONE, 1185 * ETAAB(KOFFRE),NTOTRE) 1186C 1187 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),-ONE, 1188 * DAB(KOFF3),NTOTC,XINTAB(KOFF4),NTOTC,ONE, 1189 * ETAAB(KOFFRE),NTOTRE) 1190C 1191 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),ONE, 1192 * XINTAB(KOFF5),NTOTA,DAB(KOFF6),NTOTB,ONE, 1193 * ETAAB(KOFFRE),NTOTRE) 1194C 1195 100 CONTINUE 1196C 1197 DO 110 ISYMA = 1,NSYM 1198C 1199C--------------------------------------------------------------------------- 1200C Calculate terms involving T1AM and intermediate looping over occupied. 1201C Note that the intermediate of integrals and densities is used for both 1202C contributions as ISYMA=ISYMB. 1203C--------------------------------------------------------------------------- 1204C 1205 ISYMB = ISYMA 1206 ISYMJ = ISYMB 1207 ISYMK = MULD2H(ISYMA,ISYM) 1208C 1209 KOFFRE = IMATAB(ISYMA,ISYMB) + 1 1210C 1211 NTOTRE = MAX(NVIR(ISYMA),1) 1212 NTOTA = MAX(NVIR(ISYMA),1) 1213 NTOTB = MAX(NVIR(ISYMB),1) 1214 NTOTK = MAX(NRHF(ISYMK),1) 1215 NTOTJ = MAX(NRHF(ISYMJ),1) 1216C 1217C---------------------------------- 1218C Work space allocation one. 1219C---------------------------------- 1220C 1221 KSCRHD = 1 1222 KEND1 = KSCRHD + NVIR(ISYMA)*NRHF(ISYMJ) 1223 LWRK1 = LWORK - KEND1 1224C 1225 IF (LWRK1 .LT. 0) THEN 1226 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1227 CALL QUIT( 1228 * 'Insufficient memory for allocation (1) in CC2_ETAB') 1229 ENDIF 1230C 1231 CALL DZERO(WORK(KSCRHD),NVIR(ISYMA)*NRHF(ISYMJ)) 1232C 1233 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1234 KOFF2 = IMATIJ(ISYMK,ISYMJ) + 1 1235 KOFF3 = IMATIJ(ISYMJ,ISYMK) + 1 1236 KOFF4 = IT1AM(ISYMB,ISYMJ) + 1 1237C 1238 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMJ),NRHF(ISYMK),ONE, 1239 * XINTIA(KOFF1),NTOTA,DIJ(KOFF2),NTOTK,ZERO, 1240 * WORK(KSCRHD),NTOTA) 1241C 1242 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMJ),NRHF(ISYMK),-ONE, 1243 * DAI(KOFF1),NTOTA,XINTIJ(KOFF3),NTOTJ,ONE, 1244 * WORK(KSCRHD),NTOTA) 1245C 1246 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMJ),ONE, 1247 * WORK(KSCRHD),NTOTA,T1AM(KOFF4),NTOTB,ONE, 1248 * ETAAB(KOFFRE),NTOTRE) 1249C 1250 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMJ),-ONE, 1251 * T1AM(KOFF4),NTOTA,WORK(KSCRHD),NTOTB,ONE, 1252 * ETAAB(KOFFRE),NTOTRE) 1253C 1254 110 CONTINUE 1255C 1256 DO 120 ISYMA = 1,NSYM 1257C 1258C--------------------------------------------------------------------------- 1259C Calculate terms involving T1AM and intermediate looping over virtual. 1260C Note that the intermediate of integrals and densities is used for both 1261C contributions as ISYMA=ISYMB. 1262C--------------------------------------------------------------------------- 1263C 1264 ISYMB = ISYMA 1265 ISYMJ = ISYMB 1266 ISYMC = MULD2H(ISYMA,ISYM) 1267C 1268 KOFFRE = IMATAB(ISYMA,ISYMB) + 1 1269C 1270 NTOTRE = MAX(NVIR(ISYMA),1) 1271 NTOTA = MAX(NVIR(ISYMA),1) 1272 NTOTB = MAX(NVIR(ISYMB),1) 1273 NTOTC = MAX(NVIR(ISYMC),1) 1274C 1275C---------------------------------- 1276C Work space allocation one. 1277C---------------------------------- 1278C 1279 KSCRHD = 1 1280 KEND1 = KSCRHD + NVIR(ISYMA)*NRHF(ISYMJ) 1281 LWRK1 = LWORK - KEND1 1282C 1283 IF (LWRK1 .LT. 0) THEN 1284 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1285 CALL QUIT( 1286 * 'Insufficient memory for allocation (2) in CC2_ETAB') 1287 ENDIF 1288C 1289 CALL DZERO(WORK(KSCRHD),NVIR(ISYMA)*NRHF(ISYMJ)) 1290C 1291 KOFF1 = IMATAB(ISYMC,ISYMA) + 1 1292 KOFF2 = IT1AM(ISYMC,ISYMJ) + 1 1293 KOFF3 = IMATAB(ISYMA,ISYMC) + 1 1294 KOFF4 = IT1AM(ISYMB,ISYMJ) + 1 1295C 1296 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMJ),NVIR(ISYMC),ONE, 1297 * XINTAB(KOFF1),NTOTC,DAI(KOFF2),NTOTC,ZERO, 1298 * WORK(KSCRHD),NTOTA) 1299C 1300 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 1301 * DAB(KOFF3),NTOTA,XINTIA(KOFF2),NTOTC,ONE, 1302 * WORK(KSCRHD),NTOTA) 1303C 1304 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMJ),ONE, 1305 * WORK(KSCRHD),NTOTA,T1AM(KOFF4),NTOTB,ONE, 1306 * ETAAB(KOFFRE),NTOTRE) 1307C 1308 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMJ),-ONE, 1309 * T1AM(KOFF4),NTOTA,WORK(KSCRHD),NTOTB,ONE, 1310 * ETAAB(KOFFRE),NTOTRE) 1311C 1312 120 CONTINUE 1313C 1314 CALL QEXIT('CC2_ETAB') 1315C 1316 RETURN 1317 END 1318