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