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! 18 SUBROUTINE CC_2EEXP_2(WORK,LWORK,IOPREL) 19C 20C Written by S. Coriani 05/04/2003, based on CC_2EEXP 21C 22C Version: 1.0 23C 24C Purpose: To calculate the Breit-Pauli orbit-orbit 25C contribution, the two-electron Darwin 26C and the derivative two-electron DPT 27C contribution to the energy 28C 29C Current models: CCS, MP2, CCD, CCSD, CCSD(T) 30C 31C CC2 (without frozen core) by A. Halkier & S. Coriani 20/01-2000. 32C CCSD(T) by S. Coriani 05/04/2003. 33C OBS: CCSD(T) is treated as an addition to CCSD 34C CORRONLY option introduced: it yields the correlation-only 35C contribution to the property (no Hartree-Fock part) 36C for CCD, CCSD 37C 38#include "implicit.h" 39#include "priunit.h" 40#include "maxash.h" 41#include "maxorb.h" 42#include "mxcent.h" 43#include "maxaqn.h" 44#include "aovec.h" 45#include "iratdef.h" 46#include "nuclei.h" 47#include "symmet.h" 48#include "chrnos.h" 49#include "eridst.h" 50 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 51 PARAMETER (FOUR = 4.0D0) 52 LOGICAL SAVDIR, LEX, SAVHER, OLDDX 53 DIMENSION INDEXA(MXCORB_CC) 54 DIMENSION IADR(MXCORB_CC,MXDIST) 55 DIMENSION WORK(LWORK) 56 CHARACTER*8 LABEL 57 CHARACTER*10 MODEL 58! 59 INTEGER LUPTIA 60 CHARACTER*5 FNDPTIA 61 logical LCCSD_SAV, LOCDBG 62 PARAMETER (LOCDBG =.FALSE.) 63! 64#include "ccorb.h" 65#include "infind.h" 66#include "blocks.h" 67#include "ccfield.h" 68#include "ccfop.h" 69#include "ccsdinp.h" 70#include "ccsdsym.h" 71#include "ccsdio.h" 72#include "distcl.h" 73#include "cbieri.h" 74#include "eritap.h" 75#include "cclr.h" 76#include "ccfro.h" 77#include "drw2el.h" 78C 79 CALL QENTER('CC_2EEXP_2') 80C 81C------------------------------------------------------------- 82C corronly option disabled for FROZEN core as not fully tested 83C------------------------------------------------------------- 84 IF (FROIMP) THEN 85 CORRONLY = .FALSE. 86 END IF 87C 88C----------------------------- 89C Trick for CCSD(T) 90C----------------------------- 91C 92 IF (CCPT) THEN 93 LCCSD_SAV = CCSD 94 CCSD = .TRUE. 95 END IF 96C 97C------------------------------ 98C Initialization of result. 99C------------------------------ 100C 101 IF (IPRINT .GT. 9) CALL AROUND('Entering CC_2EEXP_2') 102 CALL FLSHFO(LUPRI) 103 104 RE2DAR = ZERO 105 IF (IOPREL .EQ. 1) RELCO1 = WORK(1) !the sum MV + D1E is passed inside this routine 106C 107C----------------------------------------- 108C Initialization of timing parameters. 109C----------------------------------------- 110C 111 TIMTOT = ZERO 112 TIMTOT = SECOND() 113 TIMDEN = ZERO 114 TIMDAO = ZERO 115 TIRDAO = ZERO 116 TIMHE2 = ZERO 117 TIMONE = ZERO 118 TIMONE = SECOND() 119C 120C---------------------------------------------------- 121C Both zeta- and t-vectors are totally symmetric. 122C---------------------------------------------------- 123C 124 ISYMTR = 1 125 ISYMOP = 1 126C 127 IF (CC2) THEN 128C 129C 130C----------------------------------- 131C Initial work space allocation. 132C----------------------------------- 133C 134 N2BSTM = 0 135 DO ISYM = 1, NSYM 136 N2BSTM = MAX(N2BSTM,N2BST(ISYM)) 137 END DO 138 139 KFCKEF = 1 140 KAODEN = KFCKEF + N2BST(1) 141 KCMO = KAODEN + N2BSTM 142 KT2AM = KCMO + NLAMDS 143 KZ2AM = KT2AM + NT2AMX 144 KLAMDP = KZ2AM + NT2SQ(1) 145 KLAMDH = KLAMDP + NLAMDT 146 KT1AM = KLAMDH + NLAMDT 147 KZ1AM = KT1AM + NT1AMX 148 KEND1 = KZ1AM + NT1AMX 149 LWRK1 = LWORK - KEND1 150C 151 IF (LWRK1 .LT. 0) THEN 152 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 153 CALL QUIT( 154 * 'Insufficient core for initial allocation in CC_2EEXP') 155 ENDIF 156C 157C------------------------------------------------------------- 158C Read MO-coefficients from interface file and reorder. 159C------------------------------------------------------------- 160C 161 LUSIFC = -993 162 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 163 & .FALSE.) 164 REWIND LUSIFC 165 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 166 READ (LUSIFC) 167 READ (LUSIFC) 168 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 169 CALL GPCLOSE(LUSIFC,'KEEP') 170C 171 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 172C 173C------------------------------------------- 174C Read zero'th order zeta amplitudes. 175C------------------------------------------- 176C 177 IOPT = 3 178 CALL CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KZ1AM),WORK(KZ2AM)) 179C 180C----------------------------------- 181C Square up zeta2 amplitudes. 182C----------------------------------- 183C 184 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 185 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 186C 187C---------------------------------------------- 188C Read zero'th order cluster amplitudes. 189C---------------------------------------------- 190C 191 IOPT = 3 192 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM)) 193C 194C------------------------------------- 195C Calculate the lambda matrices. 196C------------------------------------- 197C 198 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 199 * LWRK1) 200C 201C 202C----------------------------------------------- 203C Set up 2C-E of cluster amplitudes and save 204C in KT2AM, as we only need T(2c-e) below. 205C----------------------------------------------- 206C 207 ISYOPE = 1 208 IOPTTCME = 1 209 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME) 210 KT2AMT = KT2AM !for safety 211C 212C------------------------------- 213C Work space allocation one. 214C Note that D(ai) = ZETA(ai) 215C and both D(ia) and h(ia) 216C are stored transposed! 217C------------------------------- 218C 219 LENBAR = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) 220 * + 2*NCOFRO(1) 221C 222 KONEAI = KZ1AM 223 KONEAB = KONEAI + NT1AMX 224 KONEIJ = KONEAB + NMATAB(1) 225 KONEIA = KONEIJ + NMATIJ(1) 226 KONINT = KONEIA + NT1AMX 227 KKABAR = KONINT + N2BST(ISYMOP) 228 KDHFAO = KKABAR + LENBAR 229 KKABAO = KDHFAO + N2BST(1) 230 KINTIJ = KKABAO + N2BST(1) 231 KEND1 = KINTIJ + NMATIJ(1) 232 LWRK1 = LWORK - KEND1 233C 234 IF (LWRK1 .LT. 0) THEN 235 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 236 CALL QUIT('Insufficient core for allocation 1 in CC_2EEXP') 237 ENDIF 238C 239C 240C------------------------------------------------------ 241C Initialize remaining one electron density arrays. 242C------------------------------------------------------ 243C 244 CALL DZERO(WORK(KONEAB),NMATAB(1)) 245 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 246 CALL DZERO(WORK(KONEIA),NT1AMX) 247C 248C-------------------------------------------------------- 249C Construct remaining blocks of one electron density. 250C-------------------------------------------------------- 251C 252 CALL DZERO(WORK(KINTIJ),NMATIJ(1)) 253 CALL DIJGEN(WORK(KONEIJ),WORK(KINTIJ)) 254 CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI)) 255C 256C 257C-------------------------------------------------------- 258C Backtransform the one electron density to AO-basis. 259C-------------------------------------------------------- 260C 261 CALL DZERO(WORK(KAODEN),N2BST(1)) 262C 263 ISDEN = 1 264 CALL CC_DENAO(WORK(KAODEN),ISDEN,WORK(KONEAI),WORK(KONEAB), 265 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 266 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 267C 268C---------------------------------------------- 269C Read orbital relaxation vector from disc. 270C---------------------------------------------- 271C 272 CALL DZERO(WORK(KKABAR),LENBAR) 273C 274 LUCCK = -987 275 CALL GPOPEN(LUCCK,'CCKABAR0','UNKNOWN',' ','UNFORMATTED', 276 * IDUMMY,.FALSE.) 277 REWIND(LUCCK) 278 READ(LUCCK) (WORK(KKABAR+I-1), I = 1,LENBAR) 279 CALL GPCLOSE(LUCCK,'KEEP') 280C 281C-------------------------------------------------------------- 282C Calculate ao-transformed zeta-kappa-bar-0 and HF density. 283C-------------------------------------------------------------- 284C 285 KOFDIJ = KKABAR 286 KOFDAB = KOFDIJ + NMATIJ(1) 287 KOFDAI = KOFDAB + NMATAB(1) 288 KOFDIA = KOFDAI + NT1AMX 289C 290 ISDEN = 1 291 CALL DZERO(WORK(KKABAO),N2BST(1)) 292 CALL CC_DENAO(WORK(KKABAO),ISDEN,WORK(KOFDAI),WORK(KOFDAB), 293 * WORK(KOFDIJ),WORK(KOFDIA),ISDEN,WORK(KCMO),1, 294 * WORK(KCMO),1,WORK(KEND1),LWRK1) 295C 296 CALL CCS_D1AO(WORK(KDHFAO),WORK(KEND1),LWRK1) 297 IF (FROIMP .OR. FROEXP) THEN 298 MODEL = 'DUMMY' 299 CALL CC_FCD1AO(WORK(KDHFAO),WORK(KEND1),LWRK1,MODEL) 300 ENDIF 301C 302C------------------------------------------------------------ 303C Add orbital relaxation for effective density matrix. 304C------------------------------------------------------------ 305C 306 CALL DAXPY(N2BST(1),ONE,WORK(KKABAO),1,WORK(KAODEN),1) 307C 308 ELSE IF (CCSD .or. CCD .or. CCPT) THEN 309C 310C----------------------------------- 311C Initial work space allocation. 312C----------------------------------- 313C 314 if (corronly) then 315 WRITE(LUPRI,*)'WARNING: only CORRELATION PART computed' 316 end if 317 318 N2BSTM = 0 319 DO ISYM = 1, NSYM 320 N2BSTM = MAX(N2BSTM,N2BST(ISYM)) 321 END DO 322 323 KFCKEF = 1 324 KAODSY = KFCKEF + N2BST(1) 325 KAODEN = KAODSY + N2BSTM 326 KZ2AM = KAODEN + N2BSTM 327 KT2AM = KZ2AM + NT2SQ(1) 328 KT2AMT = KT2AM + NT2AMX 329 KLAMDP = KT2AMT + NT2AMX 330 KLAMDH = KLAMDP + NLAMDT 331 KT1AM = KLAMDH + NLAMDT 332 KZ1AM = KT1AM + NT1AMX 333 KEND1 = KZ1AM + NT1AMX 334 LWRK1 = LWORK - KEND1 335C 336 IF (LWRK1 .LT. 0) THEN 337 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 338 CALL QUIT( 339 * 'Insufficient core for first allocation in CC_2EEXP') 340 ENDIF 341C 342C---------------------------------------- 343C Read zero'th order zeta amplitudes. 344C---------------------------------------- 345C 346 IOPT = 3 347 CALL CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KZ1AM),WORK(KZ2AM)) 348C 349C-------------------------------- 350C Square up zeta2 amplitudes. 351C-------------------------------- 352C 353 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 354 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 355C 356C------------------------------------------- 357C Read zero'th order cluster amplitudes. 358C------------------------------------------- 359C 360 IOPT = 3 361 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM)) 362C 363C 364C------------------------------------------------ 365C Zero out single vectors in CCD-calculation. 366C------------------------------------------------ 367C 368 IF (CCD) THEN 369 CALL DZERO(WORK(KT1AM),NT1AMX) 370 CALL DZERO(WORK(KZ1AM),NT1AMX) 371 ENDIF 372C 373C---------------------------------- 374C Calculate the lambda matrices. 375C---------------------------------- 376C 377 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 378 * LWRK1) 379C 380C--------------------------------------- 381C Set up 2C-E of cluster amplitudes. 382C--------------------------------------- 383C 384 ISYOPE = 1 385C 386 CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1) 387 IOPTTCME = 1 388 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 389C 390C------------------------------- 391C Work space allocation one. 392C Note that D(ai) = ZETA(ai) 393C and both D(ia) and h(ia) 394C are stored transposed! 395C------------------------------- 396C 397 LENBAR = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) 398 * + 2*NCOFRO(1) 399C 400 KONEAI = KZ1AM 401 KONEAB = KONEAI + NT1AMX 402 KONEIJ = KONEAB + NMATAB(1) 403 KONEIA = KONEIJ + NMATIJ(1) 404 KXMAT = KONEIA + NT1AMX 405 KYMAT = KXMAT + NMATIJ(1) 406 KMINT = KYMAT + NMATAB(1) 407 KONINT = KMINT + N3ORHF(1) 408 KMIRES = KONINT + N2BST(ISYMOP) 409 KD1ABT = KMIRES + N3ORHF(1) 410 KD1IJT = KD1ABT + NMATAB(1) 411 KKABAR = KD1IJT + NMATIJ(1) 412 KDHFAO = KKABAR + LENBAR 413 KKABAO = KDHFAO + N2BST(1) 414 KCMO = KKABAO + N2BST(1) 415 KEND1 = KCMO + NLAMDS 416 LWRK1 = LWORK - KEND1 417C 418 IF (FROIMP) THEN 419 KCMOF = KEND1 420 KEND1 = KCMOF + NLAMDS 421 LWRK1 = LWORK - KEND1 422 ENDIF 423C 424 IF (LWRK1 .LT. 0) THEN 425 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 426 CALL QUIT('Insufficient memory for allocation 1 CC_2EEXP') 427 ENDIF 428C 429 IF (FROIMP) THEN 430C 431C---------------------------------------------- 432C Get the FULL MO coefficient matrix. 433C---------------------------------------------- 434C 435 CALL CMO_ALL(WORK(KCMOF),WORK(KEND1),LWRK1) 436C 437 ENDIF 438C 439C------------------------------------------------------ 440C Initialize remaining one electron density arrays. 441C------------------------------------------------------ 442C 443 CALL DZERO(WORK(KONEAB),NMATAB(1)) 444 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 445 CALL DZERO(WORK(KONEIA),NT1AMX) 446C 447C-------------------------------------------------------- 448C Calculate X-intermediate of zeta- and t-amplitudes. 449C-------------------------------------------------------- 450C 451 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 452 * WORK(KEND1),LWRK1) 453C 454C-------------------------------------------------------- 455C Calculate Y-intermediate of zeta- and t-amplitudes. 456C-------------------------------------------------------- 457C 458 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 459 * WORK(KEND1),LWRK1) 460C 461C-------------------------------------------------------------- 462C Construct three remaining blocks of one electron density. 463C-------------------------------------------------------------- 464C 465 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 466 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1) 467 CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT)) 468 CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI)) 469C 470C--------------------------------- 471C Set up transposed densities. 472C--------------------------------- 473C 474 CALL DCOPY(NMATAB(1),WORK(KONEAB),1,WORK(KD1ABT),1) 475 CALL DCOPY(NMATIJ(1),WORK(KONEIJ),1,WORK(KD1IJT),1) 476 CALL CC_EITR(WORK(KD1ABT),WORK(KD1IJT),WORK(KEND1),LWRK1,1) 477C 478C---------------------------------------------- 479C Read orbital relaxation vector from disc. 480C---------------------------------------------- 481C 482 CALL DZERO(WORK(KKABAR),LENBAR) 483C 484 LUCCK = -678 485 CALL GPOPEN(LUCCK,'CCKABAR0','UNKNOWN',' ', 486 * 'UNFORMATTED',IDUMMY,.FALSE.) 487 REWIND(LUCCK) 488 READ(LUCCK) (WORK(KKABAR+I-1), I = 1,LENBAR) 489 CALL GPCLOSE(LUCCK,'KEEP') 490C 491C---------------------------------------------------------- 492C Read MO-coefficients from interface file and reorder. 493C---------------------------------------------------------- 494C 495 LUSIFC = -1 496 CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ', 497 * 'UNFORMATTED',IDUMMY,.FALSE.) 498 REWIND LUSIFC 499 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 500 READ (LUSIFC) 501 READ (LUSIFC) 502 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 503 CALL GPCLOSE (LUSIFC,'KEEP') 504C 505 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 506C 507C-------------------------------------------------------------- 508C Calculate ao-transformed zeta-kappa-bar-0 and HF density. 509C-------------------------------------------------------------- 510C 511 KOFDIJ = KKABAR 512 KOFDAB = KOFDIJ + NMATIJ(1) 513 KOFDAI = KOFDAB + NMATAB(1) 514 KOFDIA = KOFDAI + NT1AMX 515C 516 ISDEN = 1 517 CALL DZERO(WORK(KKABAO),N2BST(1)) 518 CALL CC_DENAO(WORK(KKABAO),ISDEN,WORK(KOFDAI),WORK(KOFDAB), 519 * WORK(KOFDIJ),WORK(KOFDIA),ISDEN,WORK(KCMO),1, 520 * WORK(KCMO),1,WORK(KEND1),LWRK1) 521C 522 CALL CCS_D1AO(WORK(KDHFAO),WORK(KEND1),LWRK1) 523 IF (FROIMP .OR. FROEXP) THEN 524 !calculating fr.core contribution to D^HF 525 MODEL = 'DUMMY' 526 CALL CC_FCD1AO(WORK(KDHFAO),WORK(KEND1),LWRK1,MODEL) 527 ENDIF 528C 529C------------------------------------------------------------ 530C Add orbital relaxation for effective density matrix. 531C------------------------------------------------------------ 532C 533 CALL DCOPY(N2BST(1),WORK(KKABAO),1,WORK(KAODEN),1) 534 if (corronly) then 535! remove HF contribution from the density 536 write(lupri,*)'Warning: HF density removed from D_1e' 537 CALL DAXPY(N2BST(1),-ONE,WORK(KDHFAO),1,WORK(KAODEN),1) 538 end if 539 540C 541C------------------------------------------------------ 542C Add frozen core contributions to AO densities. 543C------------------------------------------------------ 544C 545 IF (FROIMP) THEN 546C 547 KOFFAI = KKABAR + NMATIJ(1) + NMATAB(1) + 2*NT1AMX 548 KOFFIA = KOFFAI + NT1FRO(1) 549 KOFFIJ = KOFFIA + NT1FRO(1) 550 KOFFJI = KOFFIJ + NCOFRO(1) 551C 552 ISDEN = 1 553 ICON = 1 554 CALL CC_D1FCB(WORK(KAODEN),WORK(KOFFIJ),WORK(KOFFJI), 555 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 556 * LWRK1,ISDEN,ICON) 557C 558 ISDEN = 1 559 ICON = 2 560 CALL CC_D1FCB(WORK(KKABAO),WORK(KOFFIJ),WORK(KOFFJI), 561 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 562 * LWRK1,ISDEN,ICON) 563C 564 ENDIF 565C 566C------------------------------------------------------------ 567C Backtransform the one electron density to AO-basis. 568C We thus have the entire effective one-electron density. 569C------------------------------------------------------------ 570C 571 ISDEN = 1 572 CALL CC_DENAO(WORK(KAODEN),ISDEN,WORK(KONEAI),WORK(KONEAB), 573 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 574 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 575 576C 577C-------------------------------------------------------------- 578C For CCSD(T) include contribution from the D_ia(T) density 579C We thus have the entire effective one-electron density 580C for CCSD(T) in AO basis 581C-------------------------------------------------------------- 582C 583 584 IF (CCPT) THEN 585 KONEIA_PT = KEND1 586 KEND1 = KONEIA_PT + NT1AMX 587 LWRK1 = LWORK - KEND1 588 IF (LWRK1 .LT. 0) THEN 589 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 590 CALL QUIT( 591 & 'Insufficient core for allocation in CC_2EEXP (CCSD(T))') 592 ENDIF 593C 594 LUPTIA = -1 595 FNDPTIA = 'DPTIA' 596 CALL WOPEN2(LUPTIA,FNDPTIA,64,0) 597C 598 IOFF = 1 599 CALL GETWA2(LUPTIA,FNDPTIA,WORK(KONEIA_PT),IOFF,NT1AMX) 600 CALL WCLOSE2(LUPTIA,FNDPTIA,'KEEP') 601 602 KONEAI_PT = KEND1 603 KONEIJ_PT = KONEAI_PT + NT1AMX 604 KONEAB_PT = KONEIJ_PT + NMATIJ(1) 605 KEND11 = KONEAB_PT + NMATAB(1) 606 KDTEST = KEND11 607 LWRK11 = LWORK - KEND11 608 IF (LWRK11 .LT. 0) THEN 609 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 610 CALL QUIT( 611 & 'Insufficient core for allocation in CCPT_GRAD0') 612 ENDIF 613 614 CALL DZERO(WORK(KONEAI_PT),NT1AMX) 615 CALL DZERO(WORK(KONEAB_PT),NMATAB(1)) 616 CALL DZERO(WORK(KONEIJ_PT),NMATIJ(1)) 617 618 ISDEN = 1 619 CALL CC_DENAO(WORK(KAODEN),ISDEN,WORK(KONEAI_PT), 620 * WORK(KONEAB_PT), 621 * WORK(KONEIJ_PT),WORK(KONEIA_PT),ISDEN, 622 * WORK(KCMO),1, 623 * WORK(KCMO),1,WORK(KEND11),LWRK11) 624 625 END IF 626C 627C-------------------------------------------------------- 628C Calculate M-intermediate of zeta- and t-amplitudes. 629C-------------------------------------------------------- 630C 631 CALL CC_MI(WORK(KMINT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 632 * WORK(KEND1),LWRK1) 633C 634C-------------------------------------------------------- 635C Calculate resorted M-intermediate M(imjk)->M(mkij). 636C-------------------------------------------------------- 637C 638 CALL CC_MIRS(WORK(KMIRES),WORK(KMINT)) 639C 640 ELSE IF (MP2) THEN 641C 642C--------------------------------- 643C First work space allocation. 644C--------------------------------- 645C 646 N2BSTM = 0 647 DO ISYM = 1, NSYM 648 N2BSTM = MAX(N2BSTM,N2BST(ISYM)) 649 END DO 650C 651 LENBAR = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NCOFRO(1) 652 * + 2*NT1FRO(1) 653C 654 KFCKEF = 1 655 KAODSY = KFCKEF + N2BST(1) 656 KAODEN = KAODSY + N2BSTM 657 KONEAI = KAODEN + N2BSTM 658 KONEAB = KONEAI + NT1AMX 659 KONEIJ = KONEAB + NMATAB(1) 660 KONEIA = KONEIJ + NMATIJ(1) 661 KCMO = KONEIA + NT1AMX 662 KKABAR = KCMO + NLAMDS 663 KDHFAO = KKABAR + LENBAR 664 KKABAO = KDHFAO + N2BST(1) 665 KLAMDH = KKABAO + N2BST(1) 666 KLAMDP = KLAMDH + NLAMDT 667 KEND1 = KLAMDP + NLAMDT 668 LWRK1 = LWORK - KEND1 669C 670 IF (LWRK1 .LT. 0) THEN 671 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 672 CALL QUIT( 673 * 'Insufficient memory for work allocation in CC_2EEXP') 674 ENDIF 675C 676C-------------------------- 677C Initialize arrays. 678C-------------------------- 679C 680 CALL DZERO(WORK(KONEAI),NT1AMX) 681 CALL DZERO(WORK(KONEAB),NMATAB(1)) 682 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 683 CALL DZERO(WORK(KONEIA),NT1AMX) 684 CALL DZERO(WORK(KKABAR),LENBAR) 685C 686C----------------------------------------------------------- 687C Calculate correlated part of MO coefficient matrix. 688C----------------------------------------------------------- 689C 690 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KONEAI), 691 * WORK(KEND1),LWRK1) 692 CALL DZERO(WORK(KONEAI),NT1AMX) 693C 694C------------------------------------------------- 695C Read orbital relaxation vector from disc. 696C------------------------------------------------- 697C 698 LUCCK = -6347 699 CALL GPOPEN(LUCCK,'CCKABAR0','UNKNOWN',' ', 700 * 'UNFORMATTED',IDUMMY,.FALSE.) 701 REWIND(LUCCK) 702 READ(LUCCK) (WORK(KKABAR+I-1), I = 1,LENBAR) 703 CALL GPCLOSE(LUCCK,'KEEP') 704C 705C---------------------------------------------------------------- 706C Set up the relaxation (correlation) part of the density. 707C---------------------------------------------------------------- 708C 709 CALL DCOPY(NMATIJ(1),WORK(KKABAR),1,WORK(KONEIJ),1) 710 CALL DCOPY(NMATAB(1),WORK(KKABAR+NMATIJ(1)),1,WORK(KONEAB),1) 711 CALL DCOPY(NT1AMX,WORK(KKABAR+NMATIJ(1)+NMATAB(1)),1, 712 * WORK(KONEAI),1) 713 CALL DCOPY(NT1AMX,WORK(KONEAI),1,WORK(KONEIA),1) 714C 715C------------------------------------- 716C Add the Hartree-Fock density. 717C------------------------------------- 718C 719 DO 80 ISYM = 1,NSYM 720 DO 85 I = 1,NRHF(ISYM) 721 NII = IMATIJ(ISYM,ISYM) + NRHF(ISYM)*(I - 1) + I 722 WORK(KONEIJ + NII - 1) = WORK(KONEIJ + NII - 1) + TWO 723 85 CONTINUE 724 80 CONTINUE 725C 726C-------------------------------------- 727C Transform density to AO basis. 728C-------------------------------------- 729C 730 CALL DZERO(WORK(KAODEN),N2BST(1)) 731C 732 ISDEN = 1 733 CALL CC_DENAO(WORK(KAODEN),ISDEN,WORK(KONEAI),WORK(KONEAB), 734 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 735 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 736C 737C-------------------------------------------------------------- 738C Calculate ao-transformed zeta-kappa-bar-0 and HF density. 739C-------------------------------------------------------------- 740C 741 KOFDIJ = KKABAR 742 KOFDAB = KOFDIJ + NMATIJ(1) 743 KOFDAI = KOFDAB + NMATAB(1) 744 KOFDIA = KOFDAI + NT1AMX 745C 746 ISDEN = 1 747 CALL DZERO(WORK(KKABAO),N2BST(1)) 748 CALL CC_DENAO(WORK(KKABAO),ISDEN,WORK(KOFDAI),WORK(KOFDAB), 749 * WORK(KOFDIJ),WORK(KOFDIA),ISDEN,WORK(KLAMDP),1, 750 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 751C 752 CALL CCS_D1AO(WORK(KDHFAO),WORK(KEND1),LWRK1) 753 IF (FROIMP .OR. FROEXP) THEN 754 MODEL = 'DUMMY' 755 CALL CC_FCD1AO(WORK(KDHFAO),WORK(KEND1),LWRK1,MODEL) 756 ENDIF 757C 758C------------------------------------------- 759C Get the FULL MO coefficient matrix. 760C------------------------------------------- 761C 762 CALL CMO_ALL(WORK(KCMO),WORK(KEND1),LWRK1) 763C 764C------------------------------------------------------ 765C Add frozen core contributions to AO densities. 766C------------------------------------------------------ 767C 768 IF (FROIMP) THEN 769C 770 KOFFAI = KKABAR + NMATIJ(1) + NMATAB(1) + 2*NT1AMX 771 KOFFIA = KOFFAI + NT1FRO(1) 772 KOFFIJ = KOFFIA + NT1FRO(1) 773 KOFFJI = KOFFIJ + NCOFRO(1) 774C 775 ISDEN = 1 776 ICON = 1 777 CALL CC_D1FCB(WORK(KAODEN),WORK(KOFFIJ),WORK(KOFFJI), 778 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 779 * LWRK1,ISDEN,ICON) 780C 781 ISDEN = 1 782 ICON = 2 783 CALL CC_D1FCB(WORK(KKABAO),WORK(KOFFIJ),WORK(KOFFJI), 784 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 785 * LWRK1,ISDEN,ICON) 786C 787 ENDIF 788C 789C---------------------------------- 790C Work space allocation two. 791C---------------------------------- 792C 793 KT2AM = KEND1 794 KZ2AM = KT2AM + NT2AMX 795 KSKOD = KZ2AM + NT2AMX 796 KEND1 = KSKOD + NT1AMX 797 LWRK1 = LWORK - KEND1 798C 799 IF (LWRK1 .LT. 0) THEN 800 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 801 CALL QUIT( 802 * 'Insufficient memory for work allocation in CC_2EEXP') 803 ENDIF 804C 805C---------------------------------------- 806C Read zero'th order zeta amplitudes. 807C---------------------------------------- 808C 809 IOPT = 3 810 CALL CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KSKOD),WORK(KZ2AM)) 811C 812C------------------------------------------- 813C Read zero'th order cluster amplitudes. 814C------------------------------------------- 815C 816 IOPT = 3 817 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KSKOD),WORK(KT2AM)) 818C 819C----------------------------------------------------------------------- 820C Set up special modified amplitudes needed in the integral loop. 821C (By doing it this way, we only need one packed vector in core 822C along with the integral distribution in the delta loop.) 823C----------------------------------------------------------------------- 824C 825 ISYOPE = 1 826 IOPTTCME = 1 827 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME) 828 CALL DSCAL(NT2AMX,TWO,WORK(KT2AM),1) 829 CALL DAXPY(NT2AMX,ONE,WORK(KZ2AM),1,WORK(KT2AM),1) 830C 831 KEND1 = KSKOD 832 LWRK1 = LWORK - KEND1 833C 834 ELSE IF (CCS) THEN 835C 836C--------------------------------- 837C First work space allocation. 838C--------------------------------- 839C 840 N2BSTM = 0 841 DO ISYM = 1, NSYM 842 N2BSTM = MAX(N2BSTM,N2BST(ISYM)) 843 END DO 844 845 KFCKEF = 1 846 KAODSY = KFCKEF + N2BST(1) 847 KAODEN = KAODSY + N2BSTM 848 KCMO = KAODEN + N2BSTM 849 KEND1 = KCMO + NLAMDS 850 LWRK1 = LWORK - KEND1 851C 852 IF (LWRK1 .LT. 0) THEN 853 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 854 CALL QUIT 855 * ('Insufficient memory for work allocation in CC_2EEXP') 856 ENDIF 857C 858 CALL CCS_D1AO(WORK(KAODEN),WORK(KEND1),LWRK1) 859 IF (FROIMP .OR. FROEXP) THEN 860 MODEL = 'DUMMY' 861 CALL CC_FCD1AO(WORK(KAODEN),WORK(KEND1),LWRK1,MODEL) 862 ENDIF 863C 864C------------------------------------------- 865C Get the FULL MO coefficient matrix. 866C------------------------------------------- 867C 868 CALL CMO_ALL(WORK(KCMO),WORK(KEND1),LWRK1) 869C 870 ENDIF 871C 872C----------------------------------------- 873C Test: calculate energy contribution. 874C----------------------------------------- 875C 876 IF (.FALSE.) THEN 877 KTEST1 = KEND1 878 KENDTS = KEND1 + N2BST(1) 879 LWRKTS = LWORK - KENDTS 880 CALL CCRHS_ONEAO(WORK(KTEST1),WORK(KENDTS),LWRKTS) 881 ECCSD1 = DDOT(N2BST(1),WORK(KTEST1),1,WORK(KAODEN),1) 882 ENDIF 883C 884 TIMONE = SECOND() - TIMONE 885 CALL FLSHFO(LUPRI) 886C 887C------------------------------------------------ 888C Start the loop over two-electron integrals. 889C------------------------------------------------ 890C 891 SAVDIR = DIRECT 892 SAVHER = HERDIR 893 DIRECT = .TRUE. 894 HERDIR = .TRUE. 895C 896C 897 IF ((IOPREL .EQ. 0).OR.(IOPREL .EQ. 1)) THEN 898 IF (DAR2EL) THEN 899 DO2DAR = .TRUE. 900 AD2DAR = .FALSE. 901 S4CENT = .FALSE. 902 ELSE 903 WRITE(LUPRI,*)'Inconsistency in OP. LABEL specification' 904 CALL QUIT('CC_2EEXP_2: Inconsistency LABEL specification') 905 ENDIF 906 ELSE IF (IOPREL .EQ. 2) THEN 907 DPTINT = .TRUE. 908 ELSE IF (IOPREL .EQ. 3) THEN 909 BPH2OO = .TRUE. 910 ENDIF 911C 912 KEND1A = KEND1 913 LWRK1A = LWRK1 914C 915 DTIME = SECOND() 916 917 IF (HERDIR) THEN 918 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 919 ELSE 920 KCCFB1 = KEND1 921 KINDXB = KCCFB1 + MXPRIM*MXCONT 922 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 923 LWRK1 = LWORK - KEND1 924C 925 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 926 * KODPP1,KODPP2,KRDPP1,KRDPP2, 927 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 928 * WORK(KEND1),LWRK1,IPRERI) 929 KEND1 = KFREE 930 LWRK1 = LFREE 931 ENDIF 932 DTIME = SECOND() - DTIME 933 TIMHE2 = TIMHE2 + DTIME 934 NTOSYM = 1 935C 936 KENDSV = KEND1 937 LWRKSV = LWRK1 938C 939 ICDEL1 = 0 940 IF (HERDIR) THEN 941 NTOT = MAXSHL 942 ELSE 943 NTOT = MXCALL 944 ENDIF 945C 946 DO 100 ILLL = 1,NTOT 947C 948C--------------------------------------------------------------- 949C Determine which delta's to be calculated in this round. 950C--------------------------------------------------------------- 951C 952 KEND1 = KENDSV 953 LWRK1 = LWRKSV 954C 955 DTIME = SECOND() 956 IF (HERDIR) THEN 957 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 958 & IPRERI) 959 ELSE 960 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 961 * WORK(KODCL1),WORK(KODCL2), 962 * WORK(KODBC1),WORK(KODBC2), 963 * WORK(KRDBC1),WORK(KRDBC2), 964 * WORK(KODPP1),WORK(KODPP2), 965 * WORK(KRDPP1),WORK(KRDPP2), 966 * WORK(KCCFB1),WORK(KINDXB), 967 * WORK(KEND1), LWRK1,IPRERI) 968 ENDIF 969 DTIME = SECOND() - DTIME 970 TIMHE2 = TIMHE2 + DTIME 971C 972 KRECNR = KEND1 973 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 974 LWRK1 = LWORK - KEND1 975 IF (LWRK1 .LT. 0) THEN 976 CALL QUIT('Insufficient core in CC_2EEXP') 977 END IF 978C 979C------------------------------------------------ 980C Loop over number of delta distributions. 981C------------------------------------------------ 982C 983 DO 110 IDEL2 = 1,NUMDIS 984C 985 IDEL = INDEXA(IDEL2) 986 ISYMD = ISAO(IDEL) 987 ISYDEN = ISYMD 988C 989C------------------------------------- 990C Work space allocation two. 991C------------------------------------- 992C 993C 994 IF (CCSD .OR. CCD .OR. CC2 .or. CCPT) THEN 995 KD2IJG = KEND1 996 KD2AIG = KD2IJG + ND2IJG(ISYDEN) 997 KD2IAG = KD2AIG + ND2AIG(ISYDEN) 998 KD2ABG = KD2IAG + ND2AIG(ISYDEN) 999 KEND2 = KD2ABG + ND2ABG(ISYDEN) 1000 LWRK2 = LWORK - KEND2 1001! 1002 IF (CCPT) THEN 1003 KD2IJG_PT = KEND2 1004 KD2AIG_PT = KD2IJG_PT + ND2IJG(ISYDEN) 1005 KD2IAG_PT = KD2AIG_PT + ND2AIG(ISYDEN) 1006 KD2ABG_PT = KD2IAG_PT + ND2AIG(ISYDEN) 1007 KEND2 = KD2ABG_PT + ND2ABG(ISYDEN) 1008 LWRK2 = LWORK - KEND2 1009 END IF 1010 if (corronly) then 1011 !write(lupri,*)'WRNG: allocating for the HF 2-e density' 1012 !call flshfo(lupri) 1013 !Sonia: frozen core not tested 1014 !if (SkipFCT) then 1015 !lenijkHF = ND2IJG(ISYDEN) 1016 !else 1017 lenijkHF = NF2IJG(ISYDEN) 1018 !end if 1019 1020 kd2IJgHF = kend2 1021 kend2 = kd2IJgHF + lenijkHF 1022 LWRK2 = LWORK - KEND2 1023 end if 1024 1025! 1026 ELSE IF (MP2) THEN 1027 KD2IJG = KEND1 1028 KD2IAG = KD2IJG + NF2IJG(ISYDEN) 1029 KEND2 = KD2IAG + ND2AIG(ISYDEN) 1030 LWRK2 = LWORK - KEND2 1031 ELSE IF (CCS) THEN 1032 KD2IJG = KEND1 1033 KEND2 = KD2IJG + NF2IJG(ISYDEN) 1034 LWRK2 = LWORK - KEND2 1035 ENDIF 1036C 1037 IF (LWRK2 .LT. 0) THEN 1038 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 1039 CALL QUIT( 1040 * 'Insufficient core for allocation 2 in CC_2EEXP') 1041 ENDIF 1042C 1043C-------------------------------------------------- 1044C Initialize two electron density arrays. 1045C-------------------------------------------------- 1046C 1047 AUTIME = SECOND() 1048C 1049 CALL DZERO(WORK(KD2IJG),NF2IJG(ISYDEN)) 1050 IF (.NOT. CCS) THEN 1051 CALL DZERO(WORK(KD2IAG),ND2AIG(ISYDEN)) 1052 IF (CCSD .or. CCD .OR. CC2 .or. CCPT) THEN 1053 CALL DZERO(WORK(KD2AIG),ND2AIG(ISYDEN)) 1054 CALL DZERO(WORK(KD2ABG),ND2ABG(ISYDEN)) 1055 CALL DZERO(WORK(KD2IJG),ND2IJG(ISYDEN)) 1056C 1057 IF (CCPT) THEN 1058 CALL DZERO(WORK(KD2IAG_PT),ND2AIG(ISYDEN)) 1059 CALL DZERO(WORK(KD2AIG_PT),ND2AIG(ISYDEN)) 1060 CALL DZERO(WORK(KD2ABG_PT),ND2ABG(ISYDEN)) 1061 CALL DZERO(WORK(KD2IJG_PT),ND2IJG(ISYDEN)) 1062 END IF 1063 if (corronly) then 1064 !write(lupri,*)'OBS: dzero d2ijgHF blck ', lenijkHF 1065 CALL DZERO(WORK(kd2ijgHF),lenijkHF) 1066 end if 1067C 1068 ENDIF 1069 ENDIF 1070C 1071C---------------------------------------------------------------- 1072C Calculate the two electron density d(pq,gamma;delta). 1073C---------------------------------------------------------------- 1074C 1075 IF (CCSD .or. CCD .or. CCPT) THEN 1076 if (corronly) then 1077 !warning: frozen core not correct yet, disabled 1078 !at the beginning 1079 if (froimp) then 1080 kcmoall = kcmof 1081 else 1082 kcmoall = kcmo 1083 end if 1084 !write(lupri,*)'Compute HF 2electron density' 1085 CALL CCS_DEN2(WORK(kd2ijgHF),WORK(kcmoall), 1086 & WORK(KEND2), LWRK2,IDEL,ISYMD) 1087 !write(lupri,*)'out of CCS_DEN2()' 1088 end if 1089 1090 CALL CC_DEN2(WORK(KD2IJG),WORK(KD2AIG),WORK(KD2IAG), 1091 & WORK(KD2ABG),WORK(KZ2AM),WORK(KT2AM), 1092 & WORK(KT2AMT),WORK(KMINT),WORK(KXMAT), 1093 & WORK(KYMAT),WORK(KONEAB),WORK(KONEAI), 1094 & WORK(KONEIA),WORK(KMIRES),WORK(KLAMDH),1, 1095 & WORK(KLAMDP),1,WORK(KEND2),LWRK2,IDEL, 1096 & ISYMD) 1097! 1098 IF (CCPT) THEN 1099 1100 IF (IOPREL.EQ.3) THEN 1101 !IF ((BPH2OO).and.(IOPREL.EQ.3)) THEN 1102 !IF ((BP2EOO).and.(IOPREL.EQ.3)) THEN 1103 CALL CC_DEN2_PTanti(WORK(KD2IJG_PT),WORK(KD2AIG_PT), 1104 & WORK(KD2IAG_PT),WORK(KD2ABG_PT), 1105 & WORK(KCMO),1,WORK(KEND2),LWRK2, 1106 & IDEL,ISYMD) 1107 ELSE 1108 CALL CC_DEN2_PT(WORK(KD2IJG_PT),WORK(KD2AIG_PT), 1109 & WORK(KD2IAG_PT),WORK(KD2ABG_PT), 1110 & WORK(KCMO),1,WORK(KEND2),LWRK2, 1111 & IDEL,ISYMD) 1112 END IF 1113 END IF 1114! 1115 ELSE IF (CC2) THEN 1116 CALL CC_DEN2(WORK(KD2IJG),WORK(KD2AIG),WORK(KD2IAG), 1117 & WORK(KD2ABG),WORK(KZ2AM),WORK(KT2AM), 1118 & WORK(KT2AMT),WORK(KEND2),WORK(KEND2), 1119 & WORK(KEND2),WORK(KONEAB),WORK(KONEAI), 1120 & WORK(KONEIA),WORK(KEND2),WORK(KLAMDH),1, 1121 & WORK(KLAMDP),1,WORK(KEND2),LWRK2,IDEL,ISYMD) 1122 ELSE IF (MP2) THEN 1123 CALL CCS_DEN2(WORK(KD2IJG),WORK(KCMO),WORK(KEND2), 1124 & LWRK2,IDEL,ISYMD) 1125 CALL MP2_DEN2(WORK(KD2IAG),WORK(KT2AM),WORK(KLAMDH), 1126 & WORK(KEND2),LWRK2,IDEL,ISYMD) 1127 ELSE IF (CCS) THEN 1128 CALL CCS_DEN2(WORK(KD2IJG),WORK(KCMO),WORK(KEND2), 1129 * LWRK2,IDEL,ISYMD) 1130 ENDIF 1131 AUTIME = SECOND() - AUTIME 1132 TIMDEN = TIMDEN + AUTIME 1133C 1134C---------------------------------------------- 1135C Work space allocation for integrals. 1136C---------------------------------------------- 1137C 1138 ISYDIS = MULD2H(ISYMD,ISYMOP) 1139C 1140 KXINT = KEND2 1141 KEND3 = KXINT + NDISAO(ISYDIS) 1142 LWRK3 = LWORK - KEND3 1143C 1144 IF (LWRK3 .LT. 0) THEN 1145 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 1146 CALL QUIT('Insufficient core for allocation in CC_2EEXP_1') 1147 ENDIF 1148C 1149C----------------------------------------- 1150C Read AO integral distribution. 1151C----------------------------------------- 1152C 1153 AUTIME = SECOND() 1154 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 1155 * WORK(KRECNR),DIRECT) 1156 AUTIME = SECOND() - AUTIME 1157 TIRDAO = TIRDAO + AUTIME 1158C 1159C 1160C--------------------------------------------------- 1161C Start loop over second AO-index (gamma). 1162C--------------------------------------------------- 1163C 1164 DO 120 ISYMG = 1, NSYM 1165 DO 130 G = 1, NBAS(ISYMG) 1166C 1167 IGAM = G + IBAS(ISYMG) 1168 ISYMPQ = MULD2H(ISYMG,ISYDEN) 1169C 1170C-------------------------------------------------------- 1171C Set addresses for 2-electron densities. 1172C-------------------------------------------------------- 1173C 1174 AUTIME = SECOND() 1175 IF (CCSD .or. CCD .OR. CC2 .or. CCPT) THEN 1176 KD2GIJ = KD2IJG + ID2IJG(ISYMPQ,ISYMG) 1177 * + NMATIJ(ISYMPQ)*(G - 1) 1178 KD2GAI = KD2AIG + ID2AIG(ISYMPQ,ISYMG) 1179 * + NT1AM(ISYMPQ)*(G - 1) 1180 KD2GAB = KD2ABG + ID2ABG(ISYMPQ,ISYMG) 1181 * + NMATAB(ISYMPQ)*(G - 1) 1182 KD2GIA = KD2IAG + ID2AIG(ISYMPQ,ISYMG) 1183 * + NT1AM(ISYMPQ)*(G - 1) 1184C 1185 IF (CCPT) THEN 1186 1187 KD2GIJ_PT = KD2IJG_PT + ID2IJG(ISYMPQ,ISYMG) 1188 & + NMATIJ(ISYMPQ)*(G - 1) 1189 KD2GAI_PT = KD2AIG_PT + ID2AIG(ISYMPQ,ISYMG) 1190 & + NT1AM(ISYMPQ)*(G - 1) 1191 KD2GAB_PT = KD2ABG_PT + ID2ABG(ISYMPQ,ISYMG) 1192 & + NMATAB(ISYMPQ)*(G - 1) 1193 KD2GIA_PT = KD2IAG_PT + ID2AIG(ISYMPQ,ISYMG) 1194 & + NT1AM(ISYMPQ)*(G - 1) 1195 END IF 1196 if (corronly) then 1197! KD2GIJHF = KD2IJGHF 1198! & + ID2IJG(ISYMPQ,ISYMG) 1199! & + NMATIJ(ISYMPQ)*(G - 1) 1200 KD2GIJHF = KD2IJGHF 1201 & + IF2IJG(ISYMPQ,ISYMG) 1202 & + NFROIJ(ISYMPQ)*(G - 1) 1203 end if 1204C 1205 ELSE IF (MP2) THEN 1206 KD2GIJ = KD2IJG + IF2IJG(ISYMPQ,ISYMG) 1207 * + NFROIJ(ISYMPQ)*(G - 1) 1208 KD2GIA = KD2IAG + ID2AIG(ISYMPQ,ISYMG) 1209 * + NT1AM(ISYMPQ)*(G - 1) 1210 ELSE IF (CCS) THEN 1211 KD2GIJ = KD2IJG + IF2IJG(ISYMPQ,ISYMG) 1212 * + NFROIJ(ISYMPQ)*(G - 1) 1213 ENDIF 1214C 1215C---------------------------------------------------------- 1216C Calculate frozen core contributions to d. 1217C---------------------------------------------------------- 1218C 1219 CALL DZERO(WORK(KAODEN),N2BST(ISYMPQ)) 1220C 1221 IF ((CCSD .or. CCPT) .AND. (FROIMP)) THEN 1222C 1223 KFD2IJ = KEND3 1224 KFD2JI = KFD2IJ + NCOFRO(ISYMPQ) 1225 KFD2AI = KFD2JI + NCOFRO(ISYMPQ) 1226 KFD2IA = KFD2AI + NT1FRO(ISYMPQ) 1227 KFD2II = KFD2IA + NT1FRO(ISYMPQ) 1228 KEND4 = KFD2II + NFROFR(ISYMPQ) 1229 LWRK4 = LWORK - KEND4 1230C 1231 IF (LWRK4 .LT. 0) THEN 1232 WRITE(LUPRI,*) 'Available:', LWORK 1233 WRITE(LUPRI,*) 'Needed:', KEND4 1234 CALL QUIT('Insufficient work space in CC_2EEXP_2') 1235 ENDIF 1236C 1237 CALL DZERO(WORK(KFD2IJ),NCOFRO(ISYMPQ)) 1238 CALL DZERO(WORK(KFD2JI),NCOFRO(ISYMPQ)) 1239 CALL DZERO(WORK(KFD2AI),NT1FRO(ISYMPQ)) 1240 CALL DZERO(WORK(KFD2IA),NT1FRO(ISYMPQ)) 1241 CALL DZERO(WORK(KFD2II),NFROFR(ISYMPQ)) 1242C 1243 IF (CCPT) THEN 1244 KFD2IJ_PT = KEND4 1245 KFD2JI_PT = KFD2IJ_PT + NCOFRO(ISYMPQ) 1246 KFD2AI_PT = KFD2JI_PT + NCOFRO(ISYMPQ) 1247 KFD2IA_PT = KFD2AI_PT + NT1FRO(ISYMPQ) 1248 KFD2II_PT = KFD2IA_PT + NT1FRO(ISYMPQ) 1249 KEND4 = KFD2II_PT + NFROFR(ISYMPQ) 1250 LWRK4 = LWORK - KEND4 1251 IF (LWRK4 .LT. 0) THEN 1252 WRITE(LUPRI,*) 'Available:', LWORK 1253 WRITE(LUPRI,*) 'Needed:', KEND4 1254 CALL QUIT( 1255 & 'Insufficient work space in CC_2EEXP_2') 1256 ENDIF 1257 CALL DZERO(WORK(KFD2IJ_PT),NCOFRO(ISYMPQ)) 1258 CALL DZERO(WORK(KFD2JI_PT),NCOFRO(ISYMPQ)) 1259 CALL DZERO(WORK(KFD2AI_PT),NT1FRO(ISYMPQ)) 1260 CALL DZERO(WORK(KFD2IA_PT),NT1FRO(ISYMPQ)) 1261 CALL DZERO(WORK(KFD2II_PT),NFROFR(ISYMPQ)) 1262 END IF 1263C 1264 CALL CC_FD2BL(WORK(KFD2II),WORK(KFD2IJ), 1265 & WORK(KFD2JI),WORK(KFD2AI), 1266 & WORK(KFD2IA),WORK(KONEIJ), 1267 & WORK(KONEAB),WORK(KONEAI), 1268 & WORK(KONEIA),WORK(KCMOF), 1269 & WORK(KLAMDH),WORK(KLAMDP), 1270 & WORK(KEND4),LWRK4,IDEL, 1271 & ISYMD,G,ISYMG) 1272C 1273 CALL CC_FD2AO(WORK(KAODEN),WORK(KFD2II), 1274 & WORK(KFD2IJ),WORK(KFD2JI), 1275 & WORK(KFD2AI),WORK(KFD2IA), 1276 & WORK(KCMOF),WORK(KLAMDH), 1277 & WORK(KLAMDP),WORK(KEND4),LWRK4, 1278 & ISYMPQ) 1279C 1280 CALL CC_D2GAF(WORK(KD2GIJ),WORK(KD2GAB), 1281 & WORK(KD2GAI),WORK(KD2GIA), 1282 & WORK(KONEIJ),WORK(KONEAB), 1283 & WORK(KONEAI),WORK(KONEIA), 1284 & WORK(KCMOF),IDEL,ISYMD,G,ISYMG) 1285 1286 IF (CCPT) THEN 1287 1288 CALL CCPT_FD2BL(WORK(KFD2II_PT),WORK(KFD2IJ_PT), 1289 & WORK(KFD2JI_PT),WORK(KFD2AI_PT), 1290 & WORK(KFD2IA_PT),WORK(KONEIA_PT), 1291 & WORK(KCMO),WORK(KCMOF), 1292 & WORK(KEND4),LWRK4, 1293 & IDEL,ISYMD,G,ISYMG) 1294 1295 CALL CC_FD2AO(WORK(KAODEN),WORK(KFD2II_PT), 1296 & WORK(KFD2IJ_PT),WORK(KFD2JI_PT), 1297 & WORK(KFD2AI_PT),WORK(KFD2IA_PT), 1298 & WORK(KCMOF),WORK(KCMO), 1299 & WORK(KCMO),WORK(KEND4),LWRK4, 1300 & ISYMPQ) 1301 1302 CALL CCPT_D2GAF(WORK(KD2GIA_PT), 1303 & WORK(KONEIA_PT),WORK(KCMOF), 1304 & IDEL,ISYMD,G,ISYMG) 1305 END IF 1306C 1307C 1308 KEND5 = KEND4 1309 LWRK5 = LWRK4 1310C 1311 ELSE 1312C 1313 KEND5 = KEND3 1314 LWRK5 = LWRK3 1315C 1316 ENDIF 1317 AUTIME = SECOND() - AUTIME 1318 TIMDEN = TIMDEN + AUTIME 1319C 1320C--------------------------------------------------------- 1321C Backtransform density fully to AO basis. 1322C--------------------------------------------------------- 1323C 1324 AUTIM1 = SECOND() 1325 IF (CCSD .OR. CCD .or. CC2 .or. CCPT) THEN 1326 CALL CC_DENAO(WORK(KAODEN),ISYMPQ, 1327 * WORK(KD2GAI),WORK(KD2GAB), 1328 * WORK(KD2GIJ),WORK(KD2GIA),ISYMPQ, 1329 * WORK(KLAMDP),1,WORK(KLAMDH),1, 1330 * WORK(KEND5),LWRK5) 1331C 1332 IF (CCPT) THEN 1333 CALL CC_DENAO(WORK(KAODEN),ISYMPQ, 1334 & WORK(KD2GAI_PT),WORK(KD2GAB_PT), 1335 & WORK(KD2GIJ_PT),WORK(KD2GIA_PT), 1336 & ISYMPQ, 1337 & WORK(KCMO),1,WORK(KCMO),1, 1338 & WORK(KEND5),LWRK5) 1339 END IF 1340 if (corronly) then 1341 !Removes the Hartree-Fock 2e terms from the 2e Density 1342! if (SkipFCT) then 1343! !remove the active-only contribution from HF 1344! KCMOALL = KCMO 1345! CALL CC_D2HFCAO(WORK(KAODEN),WORK(KD2GIJHF), 1346! & WORK(KCMOALL),WORK(KEND5), 1347! & LWRK5,ISYMPQ) 1348! else 1349 !remove whole HF gradient in fully correlated case 1350 !SONIA: redundant but different addressing..... 1351 if (froimp) then 1352 KCMOALL = KCMOF 1353 else 1354 KCMOALL = KCMO 1355 end if 1356 CALL CC_D2HFAO(WORK(KAODEN),WORK(KD2GIJHF), 1357 & WORK(KCMOALL),WORK(KEND5), 1358 & LWRK5,ISYMPQ) 1359! end if 1360 end if 1361 1362C 1363 ELSE 1364 CALL CCMP_DAO(WORK(KAODEN),WORK(KD2GIJ), 1365 * WORK(KD2GIA),WORK(KCMO), 1366 * WORK(KLAMDH),WORK(KEND5), 1367 * LWRK5,ISYMPQ) 1368 ENDIF 1369C 1370C----------------------------------------------------- 1371C Add relaxation terms to set up 1372C effective density. We thus have the 1373C entire effective 2-electron density. 1374C----------------------------------------------------- 1375C 1376 IF (.NOT. CCS) THEN 1377 ICON = 2 1378 CALL CC_D2EFF(WORK(KAODEN),G,ISYMG,IDEL,ISYMD, 1379 * WORK(KKABAO),WORK(KDHFAO),ICON) 1380 CALL CC_D2EFF(WORK(KAODEN),G,ISYMG,IDEL,ISYMD, 1381 * WORK(KDHFAO),WORK(KKABAO),ICON) 1382 ENDIF 1383 AUTIM1 = SECOND() - AUTIM1 1384 TIMDAO = TIMDAO + AUTIM1 1385C 1386C---------------------------------------------------------- 1387C Calculate 2e- contribution to the DPT Energy: 1388C using 2 e- density d in memory 1389C---------------------------------------------------------- 1390C 1391 KINTAO = KEND5 1392 KEND6 = KINTAO + N2BST(ISYMPQ) 1393 LWRK6 = LWORK - KEND6 1394 IF (LWRK6 .LT. 0) THEN 1395 WRITE(LUPRI,*) 'Available:', LWORK 1396 WRITE(LUPRI,*) 'Needed:', KEND6 1397 CALL QUIT('Insufficient work space in CC_2EEXP2') 1398 ENDIF 1399 1400 KOFFIN = KXINT + IDSAOG(ISYMG,ISYMD) 1401 * + NNBST(ISYMPQ)*(G - 1) 1402 1403 IF (IOPREL.EQ.3) THEN 1404 CALL CCSD_ASYMSQ(WORK(KOFFIN),ISYMPQ,WORK(KINTAO), 1405 & 0,0) 1406c & ISYMG,ISYMD) 1407 ELSE 1408 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ,WORK(KINTAO)) 1409 END IF 1410 1411 RE2DAR1 = HALF*DDOT(N2BST(ISYMPQ), 1412 * WORK(KAODEN),1,WORK(KINTAO),1) 1413 1414 RE2DAR = RE2DAR + HALF*DDOT(N2BST(ISYMPQ), 1415 * WORK(KAODEN),1,WORK(KINTAO),1) 1416 1417 if (locdbg) then 1418 xnorm = ddot(N2BST(ISYMPQ),WORK(KAODEN),1,WORK(KAODEN),1) 1419 xnorm1= ddot(N2BST(ISYMPQ),WORK(KINTAO),1,WORK(KINTAO),1) 1420 1421 write(lupri,*)'BPH2OO, IDEL, IGAM, ISYMD, ISYMG : ', 1422 & IDEL, IGAM, ISYMD, ISYMG 1423 write(lupri,*)'DENSNORM, INTNORM : ', xnorm, xnorm1 1424 write(lupri,*)'RE2DAR1 :', RE2DAR1 1425 1426 write(lupri,*)' ' 1427 write(lupri,*)'DENSITY BLOCK : ' 1428 do isymb = 1, nsym 1429 isyma = muld2h(isympq,isymb) 1430 koff = KAODEN + iaodis(isyma,isymb) 1431 if ((nbas(isyma).gt.0).and.(nbas(isymb).gt.0)) then 1432 call output(WORK(KOFF),1,nbas(ISYMA),1,nbas(ISYMB), 1433 * NBAS(ISYMA),NBAS(ISYMB),1,LUPRI) 1434 end if 1435 end do 1436 write(lupri,*)' ' 1437 write(lupri,*)'INTEGRALS BLOCK : ' 1438 do isymb = 1, nsym 1439 isyma = muld2h(isympq,isymb) 1440 koff = KINTAO + iaodis(isyma,isymb) 1441 if ((nbas(isyma).gt.0).and.(nbas(isymb).gt.0)) then 1442 write(lupri,*)'ISYMA, ISYMB ', ISYMA, ISYMB 1443 call output(WORK(KOFF),1,nbas(ISYMA),1,nbas(ISYMB), 1444 * NBAS(ISYMA),NBAS(ISYMB),1,LUPRI) 1445 end if 1446 end do 1447 end if 1448 1449 AUTIME = SECOND() 1450 TIRDAO = TIRDAO + AUTIME 1451C 1452 130 CONTINUE 1453 120 CONTINUE 1454 110 CONTINUE 1455 100 CONTINUE 1456C 1457C------------------------------------------------ 1458C Restore logical flags for integral program. 1459C------------------------------------------------ 1460C 1461 DIRECT = SAVDIR 1462 HERDIR = SAVHER 1463 IF (DAR2EL) DO2DAR = .FALSE. 1464 IF (IOPREL .EQ. 2) THEN 1465 DPTINT = .FALSE. 1466 ELSE IF (IOPREL .EQ. 3) THEN 1467 BPH2OO = .FALSE. 1468 ENDIF 1469C 1470C---------------------- 1471C Print out result. 1472C---------------------- 1473C 1474 IF ((IOPREL .EQ. 2).OR.(IOPREL .EQ. 3)) THEN 1475 1476 WORK(1) = RE2DAR 1477 1478 ELSE IF ((DAR2EL).AND.(IOPREL.LT.2)) THEN 1479C 1480 IF (IOPREL .EQ. 0) THEN 1481 CALL AROUND('Relativistic two-electron Darwin correction') 1482 ENDIF 1483C 1484 WRITE(LUPRI,*) ' ' 1485 WRITE(LUPRI,131) '2-elec. Darwin term:', RE2DAR 1486 WRITE(LUPRI,132) '------------------- ' 1487C 1488 IF (IOPREL .EQ. 1) THEN 1489 RELCO1 = RELCO1 + RE2DAR 1490 WRITE(LUPRI,*) ' ' 1491 WRITE(LUPRI,133) 'Total relativistic correction:', RELCO1 1492 WRITE(LUPRI,134) '----------------------------- ' 1493 ENDIF 1494C 1495 131 FORMAT(9X,A20,1X,F17.9) 1496 132 FORMAT(9X,A20) 1497 133 FORMAT(9X,A30,1X,F17.9) 1498 134 FORMAT(9X,A30) 1499C 1500 ENDIF 1501C 1502 IF (.FALSE.) THEN 1503C 1504 LUSIFC = -1 1505 CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ','UNFORMATTED', 1506 * IDUMMY,.FALSE.) 1507 REWIND LUSIFC 1508C 1509 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 1510 READ (LUSIFC) POTNUC 1511 CALL GPCLOSE (LUSIFC,'KEEP') 1512C 1513 ECCSD = ECCSD1 + RE2DAR + POTNUC 1514C 1515 WRITE(LUPRI,*) ' ' 1516 WRITE(LUPRI,*) 'Coupled Cluster energy constructed' 1517 WRITE(LUPRI,*) 'from density matrices:' 1518 WRITE(LUPRI,*) 'CCSD-energy:', ECCSD 1519 WRITE(LUPRI,*) 'H1 energy, ECCSD1 = ',ECCSD1 1520c WRITE(LUPRI,*) 'H2 energy, ECCSD2 = ',RE2DAR 1521 WRITE(LUPRI,*) 'Two-electron contribution to FODPT:',RE2DAR 1522 WRITE(LUPRI,*) 'Nuc. Pot. energy = ',POTNUC 1523C 1524 ENDIF 1525C 1526C----------------------- 1527C Write out timings. 1528C----------------------- 1529C 1530 99 TIMTOT = SECOND() - TIMTOT 1531C 1532 IF (IPRINT .GT. 3) THEN 1533 WRITE(LUPRI,*) ' ' 1534 WRITE(LUPRI,*) 'Two-electron first-order property'// 1535 * ' calculation completed' 1536 WRITE(LUPRI,*) 'Total time used in CC_2EEXP:', TIMTOT 1537 ENDIF 1538 IF (IPRINT .GT. 9) THEN 1539 WRITE(LUPRI,*) 1540 * 'Time used for setting up d(pq,ga,de):', TIMDEN 1541 WRITE(LUPRI,*) 1542 * 'Time used for full AO backtransformation:', TIMDAO 1543 WRITE(LUPRI,*) 1544 * 'Time used for reading and writing d and I:',TIRDAO 1545 WRITE(LUPRI,*) 1546 * 'Time used for calculating 2 e- AO-integrals:',TIMHE2 1547 WRITE(LUPRI,*) 1548 * 'Time used for 1 e- density & intermediates:',TIMONE 1549 ENDIF 1550C 1551 IF (CCPT) THEN 1552 CCSD = LCCSD_SAV 1553 END IF 1554C 1555 CALL QEXIT('CC_2EEXP_2') 1556C 1557 RETURN 1558 165 CALL QUIT('Error reading CCTWODEN') 1559 END 1560C 1561