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 rpa_grad */ 20 SUBROUTINE rpa_GRAD(GRADH1,GRADFS,GRADH2,WORK,LWORK) 21C 22C Written by Sonia Coriani, May 2011, based on CC_GRAD0. 23C Purpose: To calculate the various contribution to the gradient: 24C - from derivative one-electron integrals GRADH1 25C - reorthonomalization part GRADFS 26C - from derivative two-electron integrals GRADH2 27C using the RPA (RCCD, DRCCR, SOSEX, RPAx????) 28C density matrices and the new integral program! 29C 30C Current RPA models: RCCD (singlet only), DRCCD, SOSEX 31C 32#include "implicit.h" 33#include "priunit.h" 34#include "dummy.h" 35#include "maxash.h" 36#include "maxorb.h" 37#include "mxcent.h" 38#include "maxaqn.h" 39#include "aovec.h" 40#include "iratdef.h" 41#include "nuclei.h" 42#include "symmet.h" 43#include "chrnos.h" 44 DATA CHRNOS /'0','1','2','3','4','5','6','7','8','9'/ 45#include "eridst.h" 46 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 47 PARAMETER (FOUR = 4.0D0) 48 LOGICAL LOCDBG 49 PARAMETER (LOCDBG = .FALSE.) 50C-tbp: flags for noddy code generation of two-electron density, 51C symmetrized or not. 52 LOGICAL USE_DEN2_NODDY, SYMMETRIZE_DEN2_NODDY 53 PARAMETER (USE_DEN2_NODDY=.FALSE., SYMMETRIZE_DEN2_NODDY=.FALSE.) 54C 55 DIMENSION INDEXA(MXCORB_CC), INDEXB(MXCORB_CC) 56 DIMENSION IPNTAB(MXCORB_CC,2) 57 DIMENSION IADR(MXCORB_CC,MXDIST) 58 DIMENSION WORK(LWORK) 59 DIMENSION GRADH1(*),GRADFS(*),GRADH2(*) 60 CHARACTER*8 LABEL 61 CHARACTER*10 MODEL 62 LOGICAL OLDDX, DIRSAV 63 LOGICAL REPORT 64#include "ccorb.h" 65#include "infind.h" 66#include "blocks.h" 67#include "ccfield.h" 68#include "ccsdinp.h" 69#include "ccsdsym.h" 70#include "ccsdio.h" 71#include "ccinftap.h" 72#include "inftap.h" 73#include "distcl.h" 74#include "cbieri.h" 75#include "eritap.h" 76#include "cclr.h" 77#include "ccfro.h" 78C 79 CALL QENTER('RPA_GRAD') 80C 81C------------------------------ 82C Initialization of result. 83C------------------------------ 84C 85 IF (IPRINT .GT. 9 .OR. LOCDBG .OR. USE_DEN2_NODDY) THEN 86 CALL AROUND('Entering RPA_GRAD') 87 IF (USE_DEN2_NODDY) THEN 88 IF (SYMMETRIZE_DEN2_NODDY) THEN 89 WRITE(LUPRI,'(A)') 90 & 'WARNING: using noddy code for symmetrized two-electron density' 91 ELSE 92 WRITE(LUPRI,'(A)') 93 &'WARNING: using noddy code for unsymmetrized two-electron density' 94 END IF 95 END IF 96 CALL FLSHFO(LUPRI) 97 END IF 98 IF (USE_DEN2_NODDY) THEN 99 IF (NSYM.NE.1 .OR. FROIMP .OR. FROEXP) THEN 100 CALL QUIT('RPA_GRAD: noddy not available') 101 END IF 102 END IF 103C 104 CALL DZERO(GRADH2,MXCOOR) 105 CALL DZERO(GRADH1,MXCOOR) 106 CALL DZERO(GRADFS,MXCOOR) 107C 108C----------------------------------------- 109C Initialization of timing parameters. 110C----------------------------------------- 111C 112 TIMTOT = ZERO 113 TIMTOT = SECOND() 114 TIMDEN = ZERO 115 TIMDAO = ZERO 116 TIRDAO = ZERO 117 TIMHE2 = ZERO 118 TIMONE = ZERO 119 TIMONE = SECOND() 120C 121C---------------------------------------------------- 122C Both zeta- and t-vectors are totally symmetric. 123C---------------------------------------------------- 124C 125 ISYMTR = 1 126 ISYMOP = 1 127C 128 N2BSTM = 0 129 DO ISYM = 1, NSYM 130 N2BSTM = MAX(N2BSTM,N2BST(ISYM)) 131 END DO 132C 133C----------------------------------- 134C Initial work space allocation. 135C----------------------------------- 136C 137 REPORT=LOCDBG .AND. NSYM.EQ.1 138 KEND0=1 139 IF (REPORT) THEN 140 KR1_0=KEND0 141 KR1=KR1_0+N2BST(1) 142 KEND0=KR1+N2BST(1) 143 END IF 144 KFCKEF = KEND0 145 KAODSY = KFCKEF + N2BST(1) 146 KAODEN = KAODSY + N2BSTM 147 KZ2AM = KAODEN + N2BSTM 148 KT2AM = KZ2AM + NT2SQ(1) 149 KT2AMT = KT2AM + NT2AMX 150 KLAMDP = KT2AMT + NT2AMX 151 KLAMDH = KLAMDP + NLAMDT 152 KZ2TIL = KLAMDH + NLAMDT 153 KZ2PCK = KZ2TIL + NT2SQ(1) 154 KT2SQ = KZ2PCK + NT2AMX 155 KT1AM = KT2SQ + NT2SQ(1) 156 KZ1AM = KT1AM !both are zero 157 KEND1 = KZ1AM + NT1AMX 158 LWRK1 = LWORK - KEND1 159C 160 IF (LWRK1 .LT. 0) THEN 161 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 162 CALL QUIT( 163 * 'Insufficient core for first allocation in CC_GRAD2E') 164 ENDIF 165C 166C---------------------------------------- 167C Read zero'th order zeta amplitudes. 168C---------------------------------------- 169C 170 IOPT = 2 171 CALL CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KZ1AM),WORK(KZ2AM)) 172 CALL DZERO(WORK(KZ1AM),NT1AMX) 173C-------------------------------------------------------- 174C Calculate tbar_tilde = 2C-E of Tbar in squared form 175C and save a packed copy in KZ2PCK 176C For dRCCD just 2*tbar 177C-------------------------------------------------------- 178 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KZ2PCK),1) 179 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 180 if (RCCD) then 181 ISYOPE = 1 182 IOPTTCME = 1 183 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME) 184 else !the same for DRCCD and SOSEX 185 CALL DSCAL(NT2AMX,TWO,WORK(KT2AM),1) 186 end if 187 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2TIL),1) 188 189C-------------------------------- 190C Square up zeta2 amplitudes. 191C-------------------------------- 192C 193 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 194 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 195C 196C------------------------------------------- 197C Read zero'th order cluster amplitudes. 198C------------------------------------------- 199C 200 IOPT = 2 201 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM)) 202 CALL DZERO(WORK(KT1AM),NT1AMX) 203 CALL CC_T2SQ(WORK(KT2AM),WORK(KT2SQ),1) 204C 205C----------------------------------------------------------------------- 206C Calculate the lambda matrices (redundant, t1=0, but I am lazy....) 207C FIXME Sonia 208C----------------------------------------------------------------------- 209C 210 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 211 * LWRK1) 212C 213C-------------------------------------------------- 214C Set up 2C-E of cluster amplitudes (T2 tilde). 215C for RCCD and SOSEX, otherwise just 2*ampl 216C-------------------------------------------------- 217C 218 ISYOPE = 1 219 CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1) 220 if (DRCCD) then 221 if (SOSEX) then 222 IOPTTCME = 1 223 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 224 else 225 CALL DSCAL(NT2AMX,TWO,WORK(KT2AMT),1) 226 end if 227 else !if (RCCD) then 228 IOPTTCME = 1 229 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 230 end if 231C 232C------------------------------- 233C Work space allocation one. 234C Note that D(ai) = 0 = D(ia) 235C both D(ia) and h(ia) 236C are stored transposed! 237C------------------------------- 238C 239 LENBAR = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) 240 * + 2*NCOFRO(1) 241C 242 KONEAI = KZ1AM 243 KONEAB = KONEAI + NT1AMX 244 KONEIJ = KONEAB + NMATAB(1) 245 KONEIA = KONEIJ + NMATIJ(1) 246 KXMAT = KONEIA + NT1AMX 247 KYMAT = KXMAT + NMATIJ(1) 248 KONINT = KYMAT + NMATAB(1) 249 KD1ABT = KONINT + N2BST(ISYMOP) 250 KD1IJT = KD1ABT + NMATAB(1) 251 KKABAR = KD1IJT + NMATIJ(1) 252 KDHFAO = KKABAR + LENBAR 253 KKABAO = KDHFAO + N2BST(1) 254 KCMO = KKABAO + N2BST(1) 255 KEND1 = KCMO + NLAMDS 256 LWRK1 = LWORK - KEND1 257C 258 IF (FROIMP) THEN 259 KCMOF = KEND1 260 KEND1 = KCMOF + NLAMDS 261 LWRK1 = LWORK - KEND1 262 ENDIF 263C 264 IF (LWRK1 .LT. 0) THEN 265 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 266 CALL QUIT( 267 * 'Insufficient memory for allocation 1 CC_GRAD2E') 268 ENDIF 269C 270 IF (FROIMP) THEN 271C 272C---------------------------------------------- 273C Get the FULL MO coefficient matrix. 274C---------------------------------------------- 275C 276 CALL CMO_ALL(WORK(KCMOF),WORK(KEND1),LWRK1) 277C 278 ENDIF 279C 280C------------------------------------------------------ 281C Initialize remaining one electron density arrays. 282C------------------------------------------------------ 283C 284 CALL DZERO(WORK(KONEAB),NMATAB(1)) 285 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 286 CALL DZERO(WORK(KONEIA),NT1AMX) 287 CALL DZERO(WORK(KONEAI),NT1AMX) 288C 289C-------------------------------------------------------- 290C Calculate X-intermediate of zeta- and t-amplitudes. 291C Scale by 2 before dumping on D_ij 292C-------------------------------------------------------- 293C 294 CALL DZERO(WORK(KYMAT),NMATAB(1)) 295 CALL DZERO(WORK(KXMAT),NMATIJ(1)) 296 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 297 * WORK(KEND1),LWRK1) 298 CALL DSCAL(NMATIJ(1),TWO,WORK(KXMAT),1) 299C 300C-------------------------------------------------------- 301C Calculate Y-intermediate of zeta- and t-amplitudes. 302C Scale by 2 before dumping in D_ab 303C-------------------------------------------------------- 304C 305 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 306 * WORK(KEND1),LWRK1) 307 CALL DSCAL(NMATAB(1),TWO,WORK(KYMAT),1) 308C 309C-------------------------------------------------------------- 310C Construct non-zero blocks of one electron density. 311C Note that X and Y are first 2*X and 2*Y 312C but then they are scaled back to "true" values 313C-------------------------------------------------------------- 314C 315 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 316 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1) 317 CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT)) 318! Rescale X and Y back as they should be 319 CALL DSCAL(NMATAB(1),HALF,WORK(KYMAT),1) 320 CALL DSCAL(NMATIJ(1),HALF,WORK(KXMAT),1) 321C 322C--------------------------------- 323C Set up transposed densities. 324C--------------------------------- 325C 326 CALL DCOPY(NMATAB(1),WORK(KONEAB),1,WORK(KD1ABT),1) 327 CALL DCOPY(NMATIJ(1),WORK(KONEIJ),1,WORK(KD1IJT),1) 328 CALL CC_EITR(WORK(KD1ABT),WORK(KD1IJT),WORK(KEND1),LWRK1,1) 329C 330C---------------------------------------------- 331C Read orbital relaxation vector from disc. 332C---------------------------------------------- 333C 334 CALL DZERO(WORK(KKABAR),LENBAR) 335C 336 LUBAR0 = -904 337 CALL GPOPEN(LUBAR0,'CCKABAR0','OLD',' ','UNFORMATTED', 338 * IDUMMY,.FALSE.) 339 REWIND(LUBAR0) 340 READ(LUBAR0) (WORK(KKABAR+I-1), I = 1,LENBAR) 341 CALL GPCLOSE(LUBAR0,'KEEP') 342C 343C---------------------------------------------------------- 344C Read MO-coefficients from interface file and reorder. 345C---------------------------------------------------------- 346C 347 LUSIFC = -1 348 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 349 * IDUMMY,.FALSE.) 350 REWIND LUSIFC 351 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 352 READ (LUSIFC) 353 READ (LUSIFC) 354 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 355 CALL GPCLOSE (LUSIFC,'KEEP') 356C 357 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 358C 359C-------------------------------------------------------------- 360C Calculate ao-transformed zeta-kappa-bar-0 and HF density. 361C-------------------------------------------------------------- 362C 363 KOFDIJ = KKABAR 364 KOFDAB = KOFDIJ + NMATIJ(1) 365 KOFDAI = KOFDAB + NMATAB(1) 366 KOFDIA = KOFDAI + NT1AMX 367C 368 ISDEN = 1 369 CALL DZERO(WORK(KKABAO),N2BST(1)) 370 CALL CC_DENAO(WORK(KKABAO),ISDEN,WORK(KOFDAI),WORK(KOFDAB), 371 * WORK(KOFDIJ),WORK(KOFDIA),ISDEN,WORK(KCMO),1, 372 * WORK(KCMO),1,WORK(KEND1),LWRK1) 373C 374 if (locdbg) then 375 WRITE(LUPRI,*)'AO BACKTRANSFORMED KAPPABAR MATRIX' 376 WRITE(LUPRI,*)'----------------------' 377 CALL OUTPUT(WORK(KKABAO),1,NBAST,1,NBAST, 378 & NBAST,NBAST,1,LUPRI) 379 end if 380 381 CALL CCS_D1AO(WORK(KDHFAO),WORK(KEND1),LWRK1) 382 IF (FROIMP .OR. FROEXP) THEN 383 MODEL = 'DUMMY' 384 CALL CC_FCD1AO(WORK(KDHFAO),WORK(KEND1),LWRK1,MODEL) 385 ENDIF 386C 387C------------------------------------------------------------ 388C Add orbital relaxation for effective density matrix. 389C------------------------------------------------------------ 390C 391 CALL DCOPY(N2BST(1),WORK(KKABAO),1,WORK(KAODEN),1) 392C 393C------------------------------------------------------ 394C Add frozen core contributions to AO densities. 395C------------------------------------------------------ 396C 397 IF (FROIMP) THEN 398C 399 KOFFAI = KKABAR + NMATIJ(1) + NMATAB(1) + 2*NT1AMX 400 KOFFIA = KOFFAI + NT1FRO(1) 401 KOFFIJ = KOFFIA + NT1FRO(1) 402 KOFFJI = KOFFIJ + NCOFRO(1) 403C 404 ISDEN = 1 405 ICON = 1 406 CALL CC_D1FCB(WORK(KAODEN),WORK(KOFFIJ),WORK(KOFFJI), 407 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 408 * LWRK1,ISDEN,ICON) 409C 410 ISDEN = 1 411 ICON = 2 412 CALL CC_D1FCB(WORK(KKABAO),WORK(KOFFIJ),WORK(KOFFJI), 413 * WORK(KOFFAI),WORK(KOFFIA),WORK(KEND1), 414 * LWRK1,ISDEN,ICON) 415C 416 ENDIF 417C 418C------------------------------------------------------------ 419C Backtransform the one electron density to AO-basis. 420C We thus have the entire effective one-electron density. 421C------------------------------------------------------------ 422C 423 424 if (locdbg) then 425 write(lupri,*)'The IJ density of ', MODEL(1:5) 426 CALL OUTPUT(WORK(KONEIJ),1,NRHF(1),1,NRHF(1), 427 * NRHF(1),NRHF(1),1,LUPRI) 428 write(lupri,*)'The AI density of ', MODEL(1:5) 429 CALL OUTPUT(WORK(KONEAI),1,NVIR(1),1,NRHF(1), 430 * NVIR(1),NRHF(1),1,LUPRI) 431 write(lupri,*)'The IA density of ', MODEL(1:5) 432 CALL OUTPUT(WORK(KONEIA),1,NVIR(1),1,NRHF(1), 433 * NVIR(1),NRHF(1),1,LUPRI) 434 write(lupri,*)'The AB density of ', MODEL(1:5) 435 CALL OUTPUT(WORK(KONEAB),1,NVIR(1),1,NVIR(1), 436 * NVIR(1),NVIR(1),1,LUPRI) 437 end if 438 439 ISDEN = 1 440 CALL CC_DENAO(WORK(KAODEN),ISDEN,WORK(KONEAI),WORK(KONEAB), 441 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 442 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 443 444 if (locdbg) then 445 WRITE(LUPRI,*)'AO BACKTRANSFORMED 1E DENSITY MATRIX' 446 WRITE(LUPRI,*)'----------------------' 447 CALL OUTPUT(WORK(KAODEN),1,NBAST,1,NBAST, 448 & NBAST,NBAST,1,LUPRI) 449 end if 450C 451C=========================================================== 452C Derivative one-electron integral contribution to gradient: 453C=========================================================== 454C 455C------------------------------------------------------ 456C Loop over symmetry distinct atoms and contract 457C 1-electron density with h^(1) integrals (IATOM = 458C 0 for zeroth-order effective Fock matrix/energy). 459C------------------------------------------------------ 460C 461 DO IATOM = 0, NUCIND 462C 463 464 MAXCOMP = 3 465 IF (IATOM .EQ. 0) MAXCOMP = 1 466C 467 DO ICOOR = 1, MAXCOMP 468C 469 ICORSY = 1 470C 471 IF (IATOM .GT. 0) THEN 472 ISCOOR = IPTCNT(3*(IATOM-1)+ICOOR,ICORSY-1,1) 473 ELSE 474 ISCOOR = 1 475 ENDIF 476C 477 IF (ISCOOR .GT. 0) THEN 478 K1DHAM = KEND1 479 KEND2 = K1DHAM + N2BST(ICORSY) 480 LWRK2 = LWORK - KEND2 481C 482 IF (LWRK2 .LE. 0) THEN 483 CALL QUIT('Insufficient work space in RPA_GRAD.') 484 END IF 485C 486 IF (IATOM .EQ. 0) THEN 487C 488C---------------------------------------------------- 489C Test: calculate energy contribution. 490C---------------------------------------------------- 491 492 CALL DZERO(WORK(K1DHAM),N2BST(1)) 493 CALL CCRHS_ONEAO(WORK(K1DHAM),WORK(KEND2),LWRK2) 494 ECCSD1 = DDOT(N2BST(1),WORK(K1DHAM),1,WORK(KAODEN),1) 495 IF (IPRINT .GT. 10 .OR. LOCDBG) THEN 496 WRITE(LUPRI,*) 497 & 'ECCSD1: 1-e energy contribution', ECCSD1 498 ENDIF 499C 500C------------------------------------------------------------------- 501C Calculate zeroth-order effective Fock matrix contr. 502C------------------------------------------------------------------- 503C 504 DO ISYMA = 1, NSYM 505 KOFF1 = IAODIS(ISYMA,ISYMA) 506 CALL TRSREC(NBAS(ISYMA),NBAS(ISYMA), 507 & WORK(KAODEN+KOFF1),WORK(KAODSY+KOFF1)) 508 END DO 509 CALL DAXPY(N2BST(1),ONE,WORK(KAODEN),1,WORK(KAODSY),1) 510C 511 DO ISYMA = 1, NSYM 512C 513 NBASA = MAX(NBAS(ISYMA),1) 514C 515 KOFF1 = KAODSY + IAODIS(ISYMA,ISYMA) 516 KOFF2 = K1DHAM + IAODIS(ISYMA,ISYMA) 517 KOFF3 = KFCKEF + IAODIS(ISYMA,ISYMA) 518C 519 CALL DGEMM('N','N',NBAS(ISYMA),NBAS(ISYMA), 520 & NBAS(ISYMA),HALF,WORK(KOFF1),NBASA, 521 & WORK(KOFF2),NBASA,ZERO, 522 & WORK(KOFF3),NBASA) 523C 524 END DO 525C 526C-------------------------------------------------------- 527C Test trace of the effective Fock matrix: 528C-------------------------------------------------------- 529C 530 FTRACE = ZERO 531 DO ISYM = 1, NSYM 532 KOFF1 = KFCKEF + IAODIS(ISYM,ISYM) 533 DO I = 1, NBAS(ISYM) 534 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1) 535 END DO 536 END DO 537 IF (IPRINT .GT. 10 .OR. LOCDBG) THEN 538 WRITE(LUPRI,*) 539 & 'Trace of 1el cont. to eff. Fock:',FTRACE 540 ENDIF 541C 542 ELSE 543C 544 LABEL = '1DHAM'//CHRNOS(ISCOOR/100) 545 & //CHRNOS(MOD(ISCOOR,100)/10) 546 & //CHRNOS(MOD((MOD(ISCOOR,100)),10)) 547C 548 CALL CCPRPAO(LABEL,.TRUE.,WORK(K1DHAM),IRREP,ISYM,IERR, 549 & WORK(KEND2),LWRK2) 550C 551 IF (IERR.NE.0 .OR. IRREP.NE.ICORSY) THEN 552 CALL QUIT('CC_DERIV: error while reading ' 553 & //LABEL//' integrals.') 554 END IF 555C 556C--------------------------------------------------- 557C Calculate contribution to gradient: 558C--------------------------------------------------- 559C 560 GRADH1(ISCOOR) = DDOT(N2BST(1),WORK(K1DHAM),1, 561 * WORK(KAODEN),1) 562C 563 ENDIF 564C 565 IF (IPRINT .GT. 10 .OR. LOCDBG) THEN 566 WRITE (LUPRI,*) 'IATOM :', IATOM 567 WRITE (LUPRI,*) 'ICOOR :', ICOOR 568 WRITE (LUPRI,*) 'ICORSY:', ICORSY 569 WRITE (LUPRI,*) 'ISCOOR:', ISCOOR 570 WRITE (LUPRI,*) 'GRADH1:', GRADH1(ISCOOR) 571 END IF 572 END IF 573 END DO 574C 166 CONTINUE 575 END DO 576C 577 TIMONE = SECOND() - TIMONE 578 CALL FLSHFO(LUPRI) 579 IF (REPORT) THEN 580 call dcopy(n2bst(1),work(kfckef),1,work(kr1_0),1) 581 call dzero(work(kfckef),n2bst(1)) 582 END IF 583C 584C============================================== 585C Two-electron density dependent contributions: 586C============================================== 587C---------------------------------------------------- 588C Open file for effective two electron densities: 589C---------------------------------------------------- 590C 591 LDECH = N2BSTM*2 + 1 592 LUDE = -1 593 CALL GPOPEN(LUDE,'CCTWODEN','UNKNOWN','DIRECT', 594 * 'UNFORMATTED',LDECH,OLDDX) 595C 596C------------------------------------------------ 597C Start the loop over two-electron integrals: 598C------------------------------------------------ 599C 600C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 601C Change this to have non-direct option for the 602C undifferentiated integrals!! 603C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 604C 605 DIRSAV = DIRECT 606 DIRECT = .TRUE. 607C 608 KEND1A = KEND1 609 LWRK1A = LWRK1 610 611 KCCFB1 = KEND1 612 KINDXB = KCCFB1 + MXPRIM*MXCONT 613 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 614 LWRK1 = LWORK - KEND1 615C 616 NTOSYM = 1 617 DTIME = SECOND() 618 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 619 * KODPP1,KODPP2,KRDPP1,KRDPP2, 620 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 621 * WORK(KEND1),LWRK1,IPRERI) 622 KEND1 = KFREE 623 LWRK1 = LFREE 624 DTIME = SECOND() - DTIME 625 TIMHE2 = TIMHE2 + DTIME 626C 627 KENDSV = KEND1 628 LWRKSV = LWRK1 629C 630 NTOT = MXCALL 631C 632C------------------------------------------- 633C Loop over sets of delta distributions: 634C------------------------------------------- 635C 636 IF (LOCDBG) WRITE(LUPRI,*) 'number of sets:',NTOT 637C 638 ECCSD2 = ZERO 639C 640 DO 100 ILLD = 1,NTOT 641C 642 KEND1 = KENDSV 643 LWRK1 = LWRKSV 644C 645C--------------------------------------- 646C Get delta indices for this set: 647C--------------------------------------- 648C 649 CALL ERIIDX(ILLD,INDEXA,NUMDISD,WORK(KINDXB),IPRERI) 650C 651C----------------------------------------------- 652C Compute the undifferentiated integrals: 653C----------------------------------------------- 654 655 CALL ERIDI2(ILLD,INDEXA,NUMDISD,0,0, 656 * WORK(KODCL1),WORK(KODCL2), 657 * WORK(KODBC1),WORK(KODBC2), 658 * WORK(KRDBC1),WORK(KRDBC2), 659 * WORK(KODPP1),WORK(KODPP2), 660 * WORK(KRDPP1),WORK(KRDPP2), 661 * WORK(KCCFB1),WORK(KINDXB), 662 * WORK(KEND1), LWRK1,IPRERI) 663 664 KRECNR = KEND1 665 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 666 LWRK1 = LWORK - KEND1 667 IF (LWRK1 .LT. 0) THEN 668 CALL QUIT('Insufficient core in CC_GRAD2E (ERIDI2)') 669 END IF 670 671C 672C------------------------------------------------ 673C Loop over number of delta distributions: 674C------------------------------------------------ 675C 676 DO 110 IDEL2 = 1,NUMDISD 677C 678 IDEL = INDEXA(IDEL2) 679 ISYMD = ISAO(IDEL) 680 ISYDEN = ISYMD 681C 682C------------------------------------- 683C Work space allocation two: 684C------------------------------------- 685C 686 IF (RCCD.OR.DRCCD) THEN 687 KD2IJG = KEND1 688 KD2AIG = KD2IJG + ND2IJG(ISYDEN) 689 KD2IAG = KD2AIG + ND2AIG(ISYDEN) 690 KD2ABG = KD2IAG + ND2AIG(ISYDEN) 691 KEND2 = KD2ABG + ND2ABG(ISYDEN) 692 ENDIF 693 LWRK2 = LWORK - KEND2 694C 695 IF (LWRK2 .LT. 0) THEN 696 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 697 CALL QUIT( 698 * 'Insufficient core for allocation 2 in RPA_GRAD2E') 699 ENDIF 700C 701C-------------------------------------------------- 702C Initialize two electron density arrays. 703C-------------------------------------------------- 704C 705 AUTIME = SECOND() 706C 707 CALL DZERO(WORK(KD2IJG),NF2IJG(ISYDEN)) 708 CALL DZERO(WORK(KD2IAG),ND2AIG(ISYDEN)) 709 CALL DZERO(WORK(KD2AIG),ND2AIG(ISYDEN)) 710 CALL DZERO(WORK(KD2ABG),ND2ABG(ISYDEN)) 711 CALL DZERO(WORK(KD2IJG),ND2IJG(ISYDEN)) 712C 713C---------------------------------------------------------------- 714C Calculate the two electron density d(pq,gamma;delta). 715C---------------------------------------------------------------- 716C 717! IF (RCCD.or.DRCCD) THEN 718! IF (USE_DEN2_NODDY) THEN 719! call drpa_den2_gd_noddy(WORK(KD2IJG),WORK(KD2AIG), 720! & WORK(KD2IAG),WORK(KD2ABG), 721! & WORK(KLAMDP),WORK(KEND2), 722! & LWRK2,IDEL, 723! & SYMMETRIZE_DEN2_NODDY) 724! ELSE 725 CALL CC_DEN2_RCCD(WORK(KD2IJG),WORK(KD2AIG), 726 & WORK(KD2IAG),WORK(KD2ABG), 727 & WORK(KZ2PCK),WORK(KZ2AM), 728 & WORK(KT2AM),WORK(KT2AMT),WORK(KT2SQ), 729 & WORK(KZ2TIL),WORK(KXMAT), 730 & WORK(KYMAT),WORK(KONEAB),WORK(KONEAI), 731 & WORK(KONEIA),WORK(KLAMDH),1, 732 & WORK(KLAMDP),1,WORK(KEND2),LWRK2,IDEL, 733 & ISYMD) 734! END IF 735! ELSE 736! CALL QUIT('RPA_GRAD: only RCCD or DRCCD') 737! ENDIF 738 AUTIME = SECOND() - AUTIME 739 TIMDEN = TIMDEN + AUTIME 740C 741C------------------------------------------------------------------- 742C Read in undifferentiated 2e-integrals for Eff. Fock matrix: 743C------------------------------------------------------------------- 744C 745 KXINT = KEND2 746 KEND3 = KXINT + NDISAO(ISYMD) 747 LWRK3 = LWORK - KEND3 748 749 IF (LWRK3 .LT. 0) THEN 750 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 751 CALL QUIT('Insufficient memory in RPA_GRAD') 752 END IF 753 754 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 755 * WORK(KRECNR),DIRECT) 756C 757C--------------------------------------------------- 758C Start loop over second AO-index (gamma). 759C--------------------------------------------------- 760C 761C for a 2-index approach, but it does not work, IDEL2 and ILLG are 762C interchanged!!! 763C 764C DO 120 ILLG = 1, NTOT 765C 766C CALL ERIIDX(ILLG,INDEXB,NUMDISG,WORK(KINDXB),IPRERI) 767C 768C DO 130 IGAM2 = 1, NUMDISG 769C 770C IGAM = INDEXB(IGAM2) 771C ISYMG = ISAO(IGAM) 772C ISYMPQ = MULD2H(ISYMG,ISYDEN) 773C G = IGAM - IBAS(ISYMG) 774 775 DO 120 ISYMG = 1, NSYM 776 DO 130 G = 1, NBAS(ISYMG) 777 IGAM = G + IBAS(ISYMG) 778 ISYMPQ = MULD2H(ISYMG,ISYDEN) 779C 780C-------------------------------------------------------- 781C Set addresses for 2-electron densities. 782C-------------------------------------------------------- 783C 784 AUTIME = SECOND() 785 KD2GIJ = KD2IJG + ID2IJG(ISYMPQ,ISYMG) 786 & + NMATIJ(ISYMPQ)*(G - 1) 787 KD2GAI = KD2AIG + ID2AIG(ISYMPQ,ISYMG) 788 & + NT1AM(ISYMPQ)*(G - 1) 789 KD2GAB = KD2ABG + ID2ABG(ISYMPQ,ISYMG) 790 & + NMATAB(ISYMPQ)*(G - 1) 791 KD2GIA = KD2IAG + ID2AIG(ISYMPQ,ISYMG) 792 & + NT1AM(ISYMPQ)*(G - 1) 793C 794C---------------------------------------------------------- 795C Calculate frozen core contributions to d. 796C---------------------------------------------------------- 797C 798 CALL DZERO(WORK(KAODEN),N2BST(ISYMPQ)) 799C 800! IF ((CCSD) .AND. (FROIMP)) THEN 801C 802! KFD2IJ = KEND3 803! KFD2JI = KFD2IJ + NCOFRO(ISYMPQ) 804! KFD2AI = KFD2JI + NCOFRO(ISYMPQ) 805! KFD2IA = KFD2AI + NT1FRO(ISYMPQ) 806! KFD2II = KFD2IA + NT1FRO(ISYMPQ) 807! KEND4 = KFD2II + NFROFR(ISYMPQ) 808! LWRK4 = LWORK - KEND4 809C 810! IF (LWRK4 .LT. 0) THEN 811! WRITE(LUPRI,*) 'Available:', LWORK 812! WRITE(LUPRI,*) 'Needed:', KEND4 813! CALL QUIT( 814! & 'Insufficient work space in CC_GRAD2E') 815! ENDIF 816C 817! CALL DZERO(WORK(KFD2IJ),NCOFRO(ISYMPQ)) 818! CALL DZERO(WORK(KFD2JI),NCOFRO(ISYMPQ)) 819! CALL DZERO(WORK(KFD2AI),NT1FRO(ISYMPQ)) 820! CALL DZERO(WORK(KFD2IA),NT1FRO(ISYMPQ)) 821! CALL DZERO(WORK(KFD2II),NFROFR(ISYMPQ)) 822C 823! CALL CC_FD2BL(WORK(KFD2II),WORK(KFD2IJ), 824! * WORK(KFD2JI),WORK(KFD2AI), 825! * WORK(KFD2IA),WORK(KONEIJ), 826! * WORK(KONEAB),WORK(KONEAI), 827! * WORK(KONEIA),WORK(KCMOF), 828! * WORK(KLAMDH),WORK(KLAMDP), 829! * WORK(KEND4),LWRK4,IDEL, 830! * ISYMD,G,ISYMG) 831C 832! CALL CC_FD2AO(WORK(KAODEN),WORK(KFD2II), 833! * WORK(KFD2IJ),WORK(KFD2JI), 834! * WORK(KFD2AI),WORK(KFD2IA), 835! * WORK(KCMOF),WORK(KLAMDH), 836! * WORK(KLAMDP),WORK(KEND4),LWRK4, 837! * ISYMPQ) 838C 839! CALL CC_D2GAF(WORK(KD2GIJ),WORK(KD2GAB), 840! * WORK(KD2GAI),WORK(KD2GIA), 841! * WORK(KONEIJ),WORK(KONEAB), 842! * WORK(KONEAI),WORK(KONEIA), 843! * WORK(KCMOF),IDEL,ISYMD,G,ISYMG) 844C 845! KEND5 = KEND4 846! LWRK5 = LWRK4 847C 848! ELSE 849C 850 KEND5 = KEND3 851 LWRK5 = LWRK3 852C 853! ENDIF 854 AUTIME = SECOND() - AUTIME 855 TIMDEN = TIMDEN + AUTIME 856C 857C--------------------------------------------------------- 858C Backtransform density fully to AO basis. 859C--------------------------------------------------------- 860C 861 AUTIM1 = SECOND() 862 CALL CC_DENAO(WORK(KAODEN),ISYMPQ, 863 * WORK(KD2GAI),WORK(KD2GAB), 864 * WORK(KD2GIJ),WORK(KD2GIA),ISYMPQ, 865 * WORK(KLAMDP),1,WORK(KLAMDH),1, 866 * WORK(KEND5),LWRK5) 867C 868C----------------------------------------------------- 869C Add relaxation terms to set up 870C effective density. We thus have the 871C entire effective 2-electron density. 872C----------------------------------------------------- 873C 874 IF (.NOT. CCS) THEN 875 IF (.NOT. USE_DEN2_NODDY) THEN 876 ICON = 2 877 CALL CC_D2EFF(WORK(KAODEN),G,ISYMG,IDEL,ISYMD, 878 * WORK(KKABAO),WORK(KDHFAO),ICON) 879 CALL CC_D2EFF(WORK(KAODEN),G,ISYMG,IDEL,ISYMD, 880 * WORK(KDHFAO),WORK(KKABAO),ICON) 881 END IF 882 ENDIF 883 AUTIM1 = SECOND() - AUTIM1 884 TIMDAO = TIMDAO + AUTIM1 885C---------------------------------------------------------- 886C Calculate 2e- contribution to the Energy: 887C using 2 e- density d in memory 888C---------------------------------------------------------- 889 890 KINTAO = KEND5 891 KEND6 = KINTAO + N2BST(ISYMPQ) 892 LWRK6 = LWORK - KEND6 893 894 IF (LWRK6 .LT. 0) THEN 895 WRITE(LUPRI,*) 'Available:', LWORK 896 WRITE(LUPRI,*) 'Needed:', KEND6 897 CALL QUIT('Insufficient work space in RPA_GRAD') 898 ENDIF 899 900 KOFFIN = KXINT + IDSAOG(ISYMG,ISYMD) 901 * + NNBST(ISYMPQ)*(G - 1) 902 903 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ,WORK(KINTAO)) 904 905C-tbp: 906c call header('d(alpha,beta;gamma,delta)',-1) 907c write(lupri,'(A)') '(used to compute energy)' 908c write(lupri,'(A,2I4)') 'gamma,delta=',igam,idel 909c call output(work(kaoden),1,nbas(1),1,nbas(1),nbas(1),nbas(1),1, 910c & lupri) 911 912 ECCSD2 = ECCSD2 + HALF*DDOT(N2BST(ISYMPQ), 913 * WORK(KAODEN),1,WORK(KINTAO),1) 914C-------------------------------------------------------------- 915C calculate the 2 e- density contribution to 916C the zeroth-order effective Fock matrix. 917C-------------------------------------------------------------- 918 DO ISYMP = 1, NSYM 919C 920 ISYMQ = MULD2H(ISYMP,ISYMPQ) 921C 922 KOFF1 = IAODIS(ISYMP,ISYMQ) 923 KOFF2 = IAODIS(ISYMQ,ISYMP) 924C 925 CALL TRSREC(NBAS(ISYMP),NBAS(ISYMQ), 926 & WORK(KAODEN+KOFF1),WORK(KAODSY+KOFF2)) 927C 928 END DO 929C 930 CALL DAXPY(N2BST(ISYMPQ),ONE,WORK(KAODEN),1, 931 & WORK(KAODSY),1) 932 933C-tbp: 934c call header('d(alpha,beta;gamma,delta)',-1) 935c write(lupri,'(A)') '(used to compute effective Fock matrix)' 936c write(lupri,'(A,2I4)') 'gamma,delta=',igam,idel 937c call output(work(kaodsy),1,nbas(1),1,nbas(1),nbas(1),nbas(1),1, 938c & lupri) 939C 940 DO ISYMP = 1, NSYM 941C 942 ISYMQ = MULD2H(ISYMP,ISYMPQ) 943C 944 NBASP = MAX(NBAS(ISYMP),1) 945 NBASQ = MAX(NBAS(ISYMQ),1) 946C 947 KOFF1 = KAODSY + IAODIS(ISYMP,ISYMQ) 948 KOFF2 = KINTAO + IAODIS(ISYMQ,ISYMP) 949 KOFF3 = KFCKEF + IAODIS(ISYMP,ISYMP) 950C 951 CALL DGEMM('N','N',NBAS(ISYMP), 952 & NBAS(ISYMP),NBAS(ISYMQ), 953 & HALF,WORK(KOFF1),NBASP, 954 & WORK(KOFF2),NBASQ,ONE, 955 & WORK(KOFF3),NBASP) 956C 957 END DO 958C 959C----------------------------------------------------- 960C Write effective density to disc for 961C subsequent use in integral program, 962C which performs the contraction of 963C the density with the 2 e- integrals. 964C----------------------------------------------------- 965C 966 967 AUTIME = SECOND() 968 NDAD = NBAST*(IDEL2 - 1) + IGAM 969 NDENEL = N2BST(ISYMPQ) 970 CALL DUMP2DEN(LUDE,WORK(KAODEN),NDENEL,NDAD) 971 AUTIME = SECOND() - AUTIME 972 TIRDAO = TIRDAO + AUTIME 973 974 130 CONTINUE 975 120 CONTINUE 976 110 CONTINUE 977 978C------------------------------------------------------------ 979C Derivative-two-electron-integral contribution to gradient 980C from a two-index approach. 981C------------------------------------------------------------ 982C 983C------------------------------------------------- 984C Loop over sets of gamma distributions: 985C------------------------------------------------- 986C 987 DO ILLG = 1, NTOT 988C 989C----------------------------------------------------- 990C Get sets of gammas for this set and set up 991C pointer arrays for ERI: 992C----------------------------------------------------- 993C 994 CALL ERIIDX(ILLG,INDEXB,NUMDISG,WORK(KINDXB),IPRERI) 995C 996 CALL IZERO(IPNTAB,2*MXCORB_CC) 997C 998 DO IDEL2 = 1, NUMDISD 999 IDEL = INDEXA(IDEL2) 1000 IPNTAB(IDEL,1) = IDEL2 1001 END DO 1002C 1003 DO IGAM2 = 1, NUMDISG 1004 IGAM = INDEXB(IGAM2) 1005 IPNTAB(IGAM,2) = IGAM2 1006 END DO 1007C 1008C------------------------------------------------------------- 1009C Work space allocation 1010C------------------------------------------------------------- 1011C Loop over number of delta and gamma distributions: 1012 KDENSITY = KEND1 1013 LDENSITY = NUMDISD*NUMDISG*NBAST*NBAST 1014C 1015 KEND2 = KDENSITY + LDENSITY 1016 LWRK2 = LWORK - KEND2 1017C 1018 KSCRATCH = KEND2 1019 KEND3 = KSCRATCH + NBAST*NBAST 1020 LWRK3 = LWORK - KEND3 1021C 1022 IF (LWRK3 .LT. 0) THEN 1023 WRITE(LUPRI,*) 'Available:', LWORK 1024 WRITE(LUPRI,*) 'Needed:', KEND3 1025 CALL QUIT('Insufficient work space in RPA_GRAD') 1026 ENDIF 1027C 1028 CALL DZERO(WORK(KDENSITY),LDENSITY) 1029C 1030C------------------------------------------------------------- 1031C Loop over number of delta and gamma distributions: 1032C------------------------------------------------------------- 1033C 1034 DO IDEL2 = 1,NUMDISD 1035C 1036 IDEL = INDEXA(IDEL2) 1037 ISYMD = ISAO(IDEL) 1038 ISYDEN = ISYMD 1039C 1040 DO IGAM2 = 1, NUMDISG 1041C 1042 IGAM = INDEXB(IGAM2) 1043 ISYMG = ISAO(IGAM) 1044C 1045 ISYMPQ = MULD2H(ISYMG,ISYDEN) 1046C 1047C---------------------------------------------- 1048C Read density block from disc. 1049C---------------------------------------------- 1050C 1051 AUTIME = SECOND() 1052 NDAD = NBAST*(IDEL2 - 1) + IGAM 1053 NDENEL = N2BST(ISYMPQ) 1054 CALL RETR2DEN(LUDE,WORK(KSCRATCH),NDENEL,NDAD) 1055 AUTIME = SECOND() - AUTIME 1056 TIRDAO = TIRDAO + AUTIME 1057C 1058C---------------------------------------------- 1059C expand density matrix: 1060C---------------------------------------------- 1061C 1062 IADRDEN = KDENSITY + 1063 & (NUMDISG*(IDEL2-1)+IGAM2-1)*NBAST*NBAST 1064C 1065 DO ISYMA = 1, NSYM 1066 ISYMB = MULD2H(ISYMPQ,ISYMA) 1067 1068 DO A = 1, NBAS(ISYMA) 1069 DO B = 1, NBAS(ISYMB) 1070 IALP = A + IBAS(ISYMA) 1071 IBET = B + IBAS(ISYMB) 1072 1073 KOFF1A = KSCRATCH-1 + IAODIS(ISYMA,ISYMB) + 1074 & NBAS(ISYMA)*(B-1) + A 1075 KOFF1B = KSCRATCH-1 + IAODIS(ISYMB,ISYMA) + 1076 & NBAS(ISYMB)*(A-1) + B 1077 KOFF2 = IADRDEN-1 + NBAST*(IBET-1) + IALP 1078 1079 WORK(KOFF2) = HALF*( WORK(KOFF1A) + WORK(KOFF1B) ) 1080 1081 END DO 1082 END DO 1083 END DO 1084 1085C-tbp: 1086c call header('d(alpha,beta;gamma,delta)',-1) 1087c write(lupri,'(A)') '(used to contract with derivative ERIs)' 1088c write(lupri,'(A,2I4)') 'gamma,delta=',igam,idel 1089c call output(work(iadrden),1,nbas(1),1,nbas(1),nbas(1),nbas(1),1, 1090c & lupri) 1091C 1092C---------------------------------------------- 1093C close loop over distributions: 1094C---------------------------------------------- 1095 END DO ! IGAM2 1096 END DO ! IDEL2 1097C 1098C--------------------------------------------------------------- 1099C Call ERI to compute derivative integrals and to 1100C contract them with the density 1101C--------------------------------------------------------------- 1102C 1103 DTIME = SECOND() 1104 CALL ERIDID(ILLD,ILLG,WORK(KDENSITY),IPNTAB, 1105 & NUMDISD,NUMDISG, 1106 * WORK(KODCL1),WORK(KODCL2), 1107 * WORK(KODBC1),WORK(KODBC2), 1108 * WORK(KRDBC1),WORK(KRDBC2), 1109 * WORK(KODPP1),WORK(KODPP2), 1110 * WORK(KRDPP1),WORK(KRDPP2), 1111 * WORK(KCCFB1),WORK(KINDXB), 1112 * WORK(KEND2), LWRK2,IPRERI) 1113 DTIME = SECOND() - DTIME 1114 TIMHE2 = TIMHE2 + DTIME 1115c 1116C-------------------------------------------- 1117C Close loops over sets of distributions: 1118C-------------------------------------------- 1119C 1120 END DO ! ILLG 1121C 1122 100 CONTINUE 1123C 1124 CALL GPCLOSE(LUDE,'DELETE') 1125c 1126C--------------------------------------------- 1127C Test trace of the effective Fock matrix: 1128C--------------------------------------------- 1129C 1130 FTRACE = ZERO 1131 DO ISYM = 1, NSYM 1132 KOFF1 = KFCKEF + IAODIS(ISYM,ISYM) 1133 DO I = 1, NBAS(ISYM) 1134 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1) 1135 END DO 1136 END DO 1137 1138 IF (LOCDBG) THEN 1139 FNORM = DDOT(N2BST(1),WORK(KFCKEF),1,WORK(KFCKEF),1) 1140 WRITE(LUPRI,*) 1141 & 'Norm^2 of the eff. Fock matrix before transformation:',FNORM 1142 WRITE(LUPRI,*) 1143 & 'Trace of the eff. Fock matrix before transformation:',FTRACE 1144 ! CALL OUTPUT(WORK(KFCKEF),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1145 END IF 1146C 1147C------------------------------------------------------ 1148C Transform effective Fock matrix to contravariant. 1149C 1150C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 1151C NOTE: Change this, so S^(0) is read in from disc. 1152C $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 1153C 1154C------------------------------------------------------ 1155C 1156 KEFFCO = KEND1A 1157 KCMODE = KEFFCO + N2BST(1) 1158 KEND2 = KCMODE + N2BST(1) 1159 LWRK2 = LWORK - KEND2 1160C 1161 IF (LWRK2 .LT. 0) THEN 1162 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1163 CALL QUIT( 1164 & 'Insufficient memory for initial allocation in RPA_GRAD') 1165 ENDIF 1166C 1167 CALL DZERO(WORK(KCMODE),N2BST(1)) 1168 CALL DCOPY(N2BST(1),WORK(KFCKEF),1,WORK(KEFFCO),1) 1169C 1170! IF ((CCSD) .AND. (FROIMP)) THEN 1171! CALL DCOPY(NLAMDS,WORK(KCMOF),1,WORK(KCMO),1) 1172! ENDIF 1173C 1174 DO ISYM = 1,NSYM 1175C 1176 NTOT = MAX(NBAS(ISYM),1) 1177C 1178 KOFF1 = KCMO + ILVISI(ISYM) 1179 KOFF2 = KCMODE + IAODIS(ISYM,ISYM) 1180C 1181 CALL DGEMM('N','T',NBAS(ISYM),NBAS(ISYM),NVIRS(ISYM),ONE, 1182 * WORK(KOFF1),NTOT,WORK(KOFF1),NTOT,ONE, 1183 * WORK(KOFF2),NTOT) 1184C 1185 KOFF3 = KCMO + ILRHSI(ISYM) 1186 KOFF4 = KCMODE + IAODIS(ISYM,ISYM) 1187C 1188 CALL DGEMM('N','T',NBAS(ISYM),NBAS(ISYM),NRHFS(ISYM),ONE, 1189 * WORK(KOFF3),NTOT,WORK(KOFF3),NTOT,ONE, 1190 * WORK(KOFF4),NTOT) 1191C 1192 END DO 1193C 1194 IF (LOCDBG) THEN 1195 XNORM = DDOT(N2BST(1),WORK(KCMODE),1,WORK(KCMODE),1) 1196 WRITE(LUPRI,*) 1197 & 'Norm^2 of the Q matrix:',XNORM 1198 call header('Qmat',-1) 1199 CALL OUTPUT(WORK(KCMODE),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1200 END IF 1201C 1202 DO ISYM = 1,NSYM 1203C 1204 NTOT = MAX(NBAS(ISYM),1) 1205C 1206 KOFF5 = KEFFCO + IAODIS(ISYM,ISYM) 1207 KOFF6 = KCMODE + IAODIS(ISYM,ISYM) 1208 KOFF7 = KFCKEF + IAODIS(ISYM,ISYM) 1209C 1210 CALL DGEMM('N','N',NBAS(ISYM),NBAS(ISYM),NBAS(ISYM),ONE, 1211 * WORK(KOFF5),NTOT,WORK(KOFF6),NTOT,ZERO, 1212 * WORK(KOFF7),NTOT) 1213C 1214 END DO 1215 IF (REPORT) THEN 1216 call dgemm('n','n',nbas(1),nbas(1),nbas(1),one, 1217 & work(kr1_0),nbas(1),work(kcmode),nbas(1),zero, 1218 & work(kr1),nbas(1)) 1219 call header('R1',-1) 1220 call output(work(kr1),1,nbas(1),1,nbas(1),nbas(1),nbas(1),1, 1221 & lupri) 1222 call header('R2',-1) 1223 call output(work(kfckef),1,nbas(1),1,nbas(1),nbas(1),nbas(1),1, 1224 & lupri) 1225 call daxpy(n2bst(1),1.0d0,work(kr1),1,work(kfckef),1) 1226 END IF 1227C 1228 IF (LOCDBG) THEN 1229 FNORM = DDOT(N2BST(1),WORK(KFCKEF),1,WORK(KFCKEF),1) 1230 WRITE(LUPRI,*) 1231 & 'Norm^2 of the zeroth-order eff. Fock matrix:',FNORM 1232 ! CALL OUTPUT(WORK(KFCKEF),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1233 END IF 1234C 1235C-------------------------------------------------------------- 1236C calculate reorthonormalization contributions to gradient. 1237C-------------------------------------------------------------- 1238C 1239 DO IATOM = 1, NUCIND 1240CAMT SKIP CENTRES WITH NO BASIS FUNCTIONS 1241CAMT OTHERWISE WE RUN INTO PROBLEMS WHEN USING 1242CAMT FLOATING FUNCTIONS 1243 !IF (NBASNUC(IATOM).EQ.0) THEN 1244 ! WRITE(LUPRI,*)'SKIPPING RE-ORTHO FOR ATOM',IATOM 1245 ! GOTO 176 1246 !ENDIF 1247 1248 DO ICOOR = 1, 3 1249C 1250 ICORSY = 1 1251 ISCOOR = IPTCNT(3*(IATOM-1)+ICOOR,ICORSY-1,1) 1252C 1253 IF (ISCOOR .GT. 0) THEN 1254C 1255 K1DOVL = KEND1A 1256 KEND2 = K1DOVL + N2BST(ICORSY) 1257 LWRK2 = LWORK - KEND2 1258C 1259 IF (LWRK2 .LE. 0) THEN 1260 CALL QUIT('Insufficient work space in CC_GRAD.') 1261 END IF 1262C 1263 LABEL = '1DOVL'//CHRNOS(ISCOOR/100) 1264 & //CHRNOS(MOD(ISCOOR,100)/10) 1265 & //CHRNOS(MOD((MOD(ISCOOR,100)),10)) 1266C 1267 CALL CCPRPAO(LABEL,.TRUE.,WORK(K1DOVL),IRREP,ISYM,IERR, 1268 & WORK(KEND2),LWRK2) 1269C 1270 IF (IERR.NE.0 .OR. IRREP.NE.ICORSY) THEN 1271 CALL QUIT('CC_DERIV: error while reading ' 1272 & //LABEL//' integrals.') 1273 ENDIF 1274C 1275 GRADFS(ISCOOR) = DDOT(N2BST(1),WORK(K1DOVL),1, 1276 * WORK(KFCKEF),1) 1277C 1278 IF (IPRINT .GT. 10 .OR. LOCDBG) THEN 1279 WRITE (LUPRI,*) 'IATOM :', IATOM 1280 WRITE (LUPRI,*) 'ICOOR :', ICOOR 1281 WRITE (LUPRI,*) 'ICORSY:', ICORSY 1282 WRITE (LUPRI,*) 'ISCOOR:', ISCOOR 1283 WRITE (LUPRI,*) 'GRADFS:', GRADFS(ISCOOR) 1284 WRITE(LUPRI,*) 1285 & 'Norm^2 of the derivative overlap matrix:', 1286 & DDOT(N2BST(1),WORK(K1DOVL),1,WORK(K1DOVL),1) 1287 WRITE(LUPRI,*) 1288 & 'Norm^2 of the zeroth-order eff. Fock matrix:', 1289 & DDOT(N2BST(1),WORK(KFCKEF),1,WORK(KFCKEF),1) 1290 CALL HEADER('Effective Fock matrix (R)',-1) 1291 CALL OUTPUT(WORK(KFCKEF),1,NBAST,1,NBAST,NBAST,NBAST, 1292 & 1,LUPRI) 1293c CALL HEADER('Derivative overlap matrix:',-1) 1294c CALL OUTPUT(WORK(K1DOVL),1,NBAST,1,NBAST,NBAST,NBAST, 1295c & 1,LUPRI) 1296 END IF 1297 END IF 1298 END DO 1299 176 CONTINUE 1300 END DO 1301C 1302C--------------------------------------------------- 1303C Read in potential energy and write out energy. 1304C--------------------------------------------------- 1305C 1306C CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1307C & .FALSE.) 1308C REWIND LUSIFC 1309C 1310C CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 1311C READ (LUSIFC) POTNUC 1312C CALL GPCLOSE(LUSIFC,'KEEP') 1313C 1314C ECCSD = ECCSD1 + ECCSD2 + POTNUC 1315C 1316C WRITE(LUPRI,*) ' ' 1317C IF (IPRINT .GT. 10 .OR. LOCDBG) THEN 1318C WRITE(LUPRI,*) 'Coupled Cluster energy constructed' 1319C WRITE(LUPRI,*) 'from density matrices:' 1320C IF (CCS) WRITE(LUPRI,*) 'CCS-energy:', ECCSD 1321C IF (MP2) WRITE(LUPRI,*) 'MP2-energy:', ECCSD 1322C IF (CCD) WRITE(LUPRI,*) 'CCD-energy:', ECCSD 1323C IF (CCSD) WRITE(LUPRI,*) 'CCSD-energy:', ECCSD 1324C WRITE(LUPRI,*) 'H1 energy, ECCSD1 = ',ECCSD1 1325C WRITE(LUPRI,*) 'H2 energy, ECCSD2 = ',ECCSD2 1326C WRITE(LUPRI,*) 'Nuc. Pot. energy = ',POTNUC 1327C WRITE(LUPRI,*) 'FTRACE = ',FTRACE 1328C 1329C FNORM = DDOT(N2BST(1),WORK(KFCKEF),1,WORK(KFCKEF),1) 1330C WRITE(LUPRI,*) 1331C & 'Norm^2 of the zeroth-order eff. Fock matrix:',FNORM 1332C ENDIF 1333C 1334C TIMREO = SECOND() - TIMREO 1335C 1336C 1337C=================================================================== 1338 1339 1340 IF (IPRINT.GT.3 .OR. LOCDBG) THEN 1341C 1342C -------------------------------- 1343C Write out timings and test info: 1344C -------------------------------- 1345C 1346 LUSIFC = -1 1347 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1348 & .FALSE.) 1349 REWIND LUSIFC 1350 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 1351 READ (LUSIFC) POTNUC 1352 CALL GPCLOSE(LUSIFC,'KEEP') 1353 1354 WRITE(LUPRI,*) ' ' 1355 WRITE(LUPRI,*) 'NUCLEAR :',POTNUC 1356 WRITE(LUPRI,*) '1e- ENERGY :',ECCSD1 1357 WRITE(LUPRI,*) '2e- ENERGY :',ECCSD2 1358 WRITE(LUPRI,*) '1+2e- ENERGY:',ECCSD1+ECCSD2 1359 WRITE(LUPRI,*) 'TOTAL :',POTNUC+ECCSD1+ECCSD2 1360 WRITE(LUPRI,*) 'FTRACE :',FTRACE 1361 1362 WRITE(LUPRI,*) ' ' 1363 WRITE(LUPRI,*) 'One electron and reorthonormalization'// 1364 * ' gradient calculation completed' 1365 WRITE(LUPRI,*) 'Two electron derivative gradient'// 1366 * ' calculation completed' 1367 1368 TIMTOT = SECOND() - TIMTOT 1369 WRITE(LUPRI,'(A,f10.2)') 'Total time used in CC_GRAD2E:',TIMTOT 1370 1371 IF (IPRINT .GT. 9) THEN 1372 WRITE(LUPRI,'(/6(/A,f10.2))') 1373 & 'Time used for setting up d(pq,ga,de) :',TIMDEN, 1374 & 'Time used for full AO backtransformation :',TIMDAO, 1375 & 'Time used for reading 2 e- AO-integrals :',TIRDAO, 1376 & 'Time used for calculating 2 e- AO-integrals:',TIMHE2, 1377 & 'Time used for 1 e- part & intermediates :',TIMONE 1378C & ,'Time used for reorthonormalization part :',TIMREO 1379 ENDIF 1380 CALL FLSHFO(LUPRI) 1381 END IF 1382C 1383C restore DIRECT flag 1384C 1385 DIRECT = DIRSAV 1386C 1387 CALL QEXIT('RPA_GRAD') 1388 RETURN 1389 165 CALL QUIT('Error reading CCTWODEN') 1390 END 1391