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_den_ptfc */ 20 SUBROUTINE CC_DEN_PTFC(POTNUC,ETAAI,ZKDIA,WORK,LWORK, 21 & IOPT,IMODEL,LTSTE) 22C 23C Written by S. Coriani, based on CC_DEN 24C Fall 2001 25C 26C Version: 2.0 27C 28C Purpose: 29C 30C For IOPT = 1 calculate the (non corrected) right hand side of the 31C equation for the zeroth order orbital rotation lagrange multiplier 32C response (Zeta-Kappa-0)_ai and return it in the array ETAAI 33C For test purposes ... 34C The (P) densities are read in from file 35C 36C For IOPT = 2 calculate ALSO the "diagonal" ZKAPPA_pp (returned in 37C ZKDIA), as well as the "diagonal" RHS eta_ij, eta_ab (local arrays) 38C and the corresponding ZKAPPA_ij, ZKAPPA_ab (also returned in ZKDIA) 39C and use them to correct ETAAI (I moved inside the call to ETACORR!!!) 40C 41C ZKAPPA_ij, ZKAPPA_ab are simply equal to eta_ij eta_ab scaled by 42C orbital energies..... 43C 44C Current models: CCSD (IMODEL=0, test purposes), CCSD(T) (IMODEL=1) 45C 46C LTSTE = .true. test densities via Energy calculation 47C 48#include "implicit.h" 49#include "priunit.h" 50#include "dummy.h" 51#include "maxorb.h" 52#include "maxash.h" 53#include "mxcent.h" 54#include "aovec.h" 55#include "iratdef.h" 56 PARAMETER (ZERO = 0.0D0,HALF=0.5D0,ONE = 1.0D0,TWO = 2.0D0) 57 PARAMETER (TRE = 3.0D0, FOUR = 4.0D0) 58 DIMENSION INDEXA(MXCORB_CC) 59 DIMENSION ETAAI(*), ZKDIA(*), WORK(LWORK) 60 LOGICAL LTSTE, LETAFI, LETIFJ 61#include "ccorb.h" 62#include "ccisao.h" 63#include "r12int.h" 64#include "inftap.h" 65#include "blocks.h" 66#include "ccfield.h" 67#include "ccsdinp.h" 68#include "ccinftap.h" 69#include "ccsdsym.h" 70#include "ccsdio.h" 71#include "distcl.h" 72#include "cbieri.h" 73#include "eritap.h" 74#include "ccfro.h" 75CTEST 76C#include "ccfop.h" 77 78 CHARACTER MODEL*10 79 CHARACTER NAME1*8 80 CHARACTER NAME2*8 81 CHARACTER*5 FNKAPAB, FNKAPIJ, FNDPTIA 82 83 LOGICAL LOCDBG 84 PARAMETER (LOCDBG=.FALSE.) 85C 86 CALL QENTER('CC_DEN_PTFC') 87 88C 89 IF (FROIMP) THEN 90C 91 NAME1 = 'CCFRO1IN' 92 NAME2 = 'CCFRO2IN' 93C 94 LUN1 = -1 95 LUN2 = -1 96C 97 CALL WOPEN2(LUN1,NAME1,64,0) 98 CALL WOPEN2(LUN2,NAME2,64,0) 99C 100 ENDIF 101C 102 IF (IOPT .LE. 2) THEN 103 CALL HEADER('Constructing RHS for CCSD(T)-kapbar-0',-1) 104 CALL HEADER('Inside ccpt_den()',-1) 105 ENDIF 106C 107C----------------------------------------- 108C Initialization of timing parameters. 109C----------------------------------------- 110C 111 TIMTOT = ZERO 112 TIMTOT = SECOND() 113 TIMDEN = ZERO 114 TIMRES = 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 LUNGO = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 128 * + 2*NCOFRO(1) + 2*NT1FRO(1) 129C 130C----------------------------------- 131C Initial work space allocation. 132C----------------------------------- 133C 134 IF (LTSTE) THEN 135 KD1AOB = 1 136 KSTART = KD1AOB + N2BST(1) 137 ELSE 138 KSTART = 1 139 END IF 140 141 KZ2AM = KSTART 142 KT2AM = KZ2AM + NT2SQ(1) 143 KT2AMT = KT2AM + NT2AMX 144 KLAMDP = KT2AMT + NT2AMX 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('Insufficient memory for initial allocation '// 154 & 'in CC_DEN_PTFC') 155 ENDIF 156C 157C---------------------------------------- 158C Read zero-th order zeta amplitudes. 159C---------------------------------------- 160C 161 IOPTRW = 3 162 CALL CC_RDRSP('L0',0,1,IOPTRW,MODEL,WORK(KZ1AM),WORK(KZ2AM)) 163C 164C-------------------------------- 165C Square up zeta2 amplitudes. 166C-------------------------------- 167C 168 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 169 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 170C 171C------------------------------------------- 172C Read zero-th order cluster amplitudes. 173C------------------------------------------- 174C 175 IOPTRW = 3 176 CALL CC_RDRSP('R0',0,1,IOPTRW,MODEL,WORK(KT1AM),WORK(KT2AM)) 177C 178C---------------------------------- 179C Calculate the lambda matrices. 180C Used in CCSD-like part 181C---------------------------------- 182C 183 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 184 * LWRK1) 185C 186C--------------------------------------- 187C Set up 2C-E of cluster amplitudes. 188C--------------------------------------- 189C 190 ISYOPE = 1 191C 192 CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1) 193 IOPTTCME = 1 194 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 195C 196C---------------------------------------------- 197C Work space allocation one. CCSD-like part 198C Note that D(ai) = ZETA(ai), and both 199C D(ia) and h(ia) are stored transposed! 200C---------------------------------------------- 201C 202 KONEAI = KZ1AM 203 KONEAB = KONEAI + NT1AMX 204 KONEIJ = KONEAB + NMATAB(1) 205 KONEIA = KONEIJ + NMATIJ(1) 206 KXMAT = KONEIA + NT1AMX 207 KYMAT = KXMAT + NMATIJ(1) 208 KMINT = KYMAT + NMATAB(1) 209 KONINT = KMINT + N3ORHF(1) 210 KMIRES = KONINT + N2BST(ISYMOP) 211C 212 KINTAI = KMIRES + N3ORHF(1) 213 KINTAB = KINTAI + NT1AMX 214 KINTIJ = KINTAB + NMATAB(1) 215 KINTIA = KINTIJ + NMATIJ(1) 216 KINABT = KINTIA + NT1AMX 217 KINIJT = KINABT + NMATAB(1) 218 KD1ABT = KINIJT + NMATIJ(1) 219 KD1IJT = KD1ABT + NMATAB(1) 220 KEND1 = KD1IJT + NMATIJ(1) 221 LWRK1 = LWORK - KEND1 222 223 IF (FROIMP) THEN 224 KFROII = KEND1 225 KFROIJ = KFROII + NFROFR(1) 226 KFROJI = KFROIJ + NCOFRO(1) 227 KFROAI = KFROJI + NCOFRO(1) 228 KFROIA = KFROAI + NT1FRO(1) 229 KFD1II = KFROIA + NT1FRO(1) 230 KEND1 = KFD1II + NFROFR(1) 231 LWRK1 = LWORK - KEND1 232 ENDIF 233 234 KCMO = KEND1 235 KEND1 = KCMO + NLAMDS 236 LWRK1 = LWORK - KEND1 237 238 IF (FROIMP) THEN 239 KCMOF = KEND1 240 KEND1 = KCMOF + NLAMDS 241 LWRK1 = LWORK - KEND1 242 ENDIF 243C 244 IF (LWRK1 .LT. 0) THEN 245 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 246 CALL QUIT('Insuff. memory for allocation 1 CC_DEN_PTFC') 247 ENDIF 248C 249 IF (FROIMP) THEN 250C 251C------------------------------------------- 252C Get the FULL MO coefficient matrix. 253C------------------------------------------- 254C 255 CALL CMO_ALL(WORK(KCMOF),WORK(KEND1),LWRK1) 256C 257 ENDIF 258 259C 260C------------------------------------------------- 261C Get the non-frozen MO coefficient matrix 262C for (T) part and reorder. 263C------------------------------------------------- 264C 265 CALL CC_GET_CMO(WORK(KCMO)) 266 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 267C 268C------------------------------------------------------ 269C Initialize remaining one electron density arrays. 270C------------------------------------------------------ 271C 272 CALL DZERO(WORK(KONEAB),NMATAB(1)) 273 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 274 CALL DZERO(WORK(KONEIA),NT1AMX) 275C 276C-------------------------------------------------------- 277C Calculate X-intermediate of zeta- and t-amplitudes. 278C-------------------------------------------------------- 279C 280 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 281 * WORK(KEND1),LWRK1) 282C 283C-------------------------------------------------------- 284C Calculate Y-intermediate of zeta- and t-amplitudes. 285C-------------------------------------------------------- 286C 287 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 288 * WORK(KEND1),LWRK1) 289C 290C------------------------------------------------------------------------ 291C Construct three remaining blocks of CCSD-like one electron density. 292C------------------------------------------------------------------------ 293C 294 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 295 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1) 296 CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT)) 297 CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI)) 298 299 IF (LOCDBG) THEN 300 DIJNO = DDOT(NMATIJ(1),WORK(KONEIJ),1,WORK(KONEIJ),1) 301 DAINO = DDOT(NT1AMX,WORK(KONEAI),1,WORK(KONEAI),1) 302 DIANO = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 303 DABNO = DDOT(NMATAB(1),WORK(KONEAB),1,WORK(KONEAB),1) 304 WRITE(LUPRI,*) 'CC_DEN_PTFC: IOPT = ', IOPT 305 WRITE(LUPRI,*) 306 & 'Norm of CCSD-like one electron densities in MO-basis:' 307 WRITE(LUPRI,*) DIJNO, DAINO, DIANO, DABNO 308 ENDIF 309C 310C--------------------------------- 311C Read one-electron integrals. 312C--------------------------------- 313C 314 CALL CCRHS_ONEAO(WORK(KONINT),WORK(KEND1),LWRK1) 315 316 IF (LTSTE) THEN 317 318 CALL DZERO(WORK(KD1AOB),N2BST(1)) 319 ISDEN = 1 320 CALL CC_DENAO(WORK(KD1AOB),ISDEN,WORK(KONEAI),WORK(KONEAB), 321 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 322 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 323C 324 IF (FROIMP) THEN 325 MODEL = 'DUMMY' 326 CALL CC_FCD1AO(WORK(KD1AOB),WORK(KEND1),LWRK1,MODEL) 327 END IF 328 329 ECCSD1 = DDOT(N2BST(ISYMOP),WORK(KD1AOB),1,WORK(KONINT),1) 330 ECCSD2 = ZERO 331 332 END IF 333C 334C--------------------------------------------------------- 335C Ove 24-20-1996 336C Read one-electron integrals into Fock-matrix for 337C finite field. 338C--------------------------------------------------------- 339C 340 DO 13 IF = 1, NFIELD 341c DTIME = SECOND() 342 FF = EFIELD(IF) 343 CALL CC_ONEP(WORK(KONINT),WORK(KEND1),LWRK1,FF,1,LFIELD(IF)) 344c DTIME = DTIME - SECOND() 345c TIMFCK = TIMFCK + DTIME 346 13 CONTINUE 347C 348C-------------------------------------------------- 349C Transform one electron integrals to MO-basis. 350C-------------------------------------------------- 351C 352 ISYM = 1 353 CALL CCDINTMO(WORK(KINTIJ),WORK(KINTIA),WORK(KINTAB), 354 * WORK(KINTAI),WORK(KONINT),WORK(KLAMDP), 355 * WORK(KLAMDH),WORK(KEND1),LWRK1,ISYM) 356C 357 IF (LOCDBG) THEN 358 359 HIJNO = DDOT(NMATIJ(1),WORK(KINTIJ),1,WORK(KINTIJ),1) 360 HAINO = DDOT(NT1AMX,WORK(KINTAI),1,WORK(KINTAI),1) 361 HIANO = DDOT(NT1AMX,WORK(KINTIA),1,WORK(KINTIA),1) 362 HABNO = DDOT(NMATAB(1),WORK(KINTAB),1,WORK(KINTAB),1) 363 WRITE(LUPRI,*) 'CC_DEN_PTFC: IOPT = ', IOPT 364 WRITE(LUPRI,*) 365 & 'UFFA: Norm of Lambda one electron integrals in MO-basis:' 366 WRITE(LUPRI,*) HIJNO, HAINO, HIANO, HABNO 367 368 TAMNOR = DDOT(NT1AM(1),WORK(KT1AM),1,WORK(KT1AM),1) 369 WRITE(LUPRI,*) 'Norm of T1AM amplitudes', TAMNOR 370 371 ENDIF 372C 373 IF (FROIMP) THEN 374C 375 ISYM = 1 376 !obtain integrals with frozen indices 377 ! h_Ij h_jI h_aJ h_Ia h_IJ 378 ! 379 CALL CCIFROMO(WORK(KFROIJ),WORK(KFROJI),WORK(KFROAI), 380 * WORK(KFROIA),WORK(KFROII),WORK(KONINT), 381 * WORK(KLAMDP),WORK(KLAMDH),WORK(KCMOF), 382 * WORK(KEND1),LWRK1,ISYM) 383C 384 !calculare D_II = 2 delta_IJ 385 CALL CCFD1II(WORK(KFD1II)) 386C 387 ENDIF 388C 389C-------------------------------------------------- 390C Set up transposed integrals and densities. 391C-------------------------------------------------- 392C 393 CALL DCOPY(NMATAB(1),WORK(KINTAB),1,WORK(KINABT),1) 394 CALL DCOPY(NMATIJ(1),WORK(KINTIJ),1,WORK(KINIJT),1) 395 CALL DCOPY(NMATAB(1),WORK(KONEAB),1,WORK(KD1ABT),1) 396 CALL DCOPY(NMATIJ(1),WORK(KONEIJ),1,WORK(KD1IJT),1) 397C 398 CALL CC_EITR(WORK(KINABT),WORK(KINIJT),WORK(KEND1),LWRK1,1) 399 CALL CC_EITR(WORK(KD1ABT),WORK(KD1IJT),WORK(KEND1),LWRK1,1) 400C 401C------------------------------------------------------------ 402C Calculate one electron contribution to Zeta-kappa-0. 403C------------------------------------------------------------ 404C 405 ISYM = 1 406 407 IF (IOPT.EQ.2) THEN 408 409 KOFFIJ = 1 410 CALL CC2_ETIJ(ZKDIA(KOFFIJ),WORK(KINTIJ),WORK(KINTAI), 411 & WORK(KINTIA),WORK(KINTAB),WORK(KONEIJ), 412 & WORK(KONEAI),WORK(KONEIA),WORK(KONEAB), 413 & WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 414 415 416 KOFFAB = NMATIJ(1) + 1 417 CALL CC2_ETAB(ZKDIA(KOFFAB),WORK(KINTIJ),WORK(KINTAI), 418 & WORK(KINTIA),WORK(KINTAB),WORK(KONEIJ), 419 & WORK(KONEAI),WORK(KONEIA),WORK(KONEAB), 420 & WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 421 422 423 END IF 424 425C------------------------------------------------------------ 426 427 CALL DZERO(ETAAI,NALLAI(1)) 428 CALL CCDENZK0(ETAAI,WORK(KINTIJ),WORK(KINTAI),WORK(KINTIA), 429 * WORK(KINTAB),WORK(KONEIJ),WORK(KONEAI), 430 * WORK(KONEIA),WORK(KONEAB),WORK(KINIJT), 431 * WORK(KINABT),WORK(KD1IJT),WORK(KD1ABT), 432 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 433C 434 435 IF (FROIMP) THEN 436C 437C-------------------------------------------------------- 438C Calculate one-electron contribution to right- 439C hand-side of correlated-frozen multipliers. 440C-------------------------------------------------------- 441C 442 ISYM = 1 443 ICON = 1 444 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 445 ! 446 !eta_iJ (one el) 447 ! 448 CALL CC_ETIJF(ZKDIA(KOFRES),WORK(KONEIJ),WORK(KONEAB), 449 * WORK(KONEAI),WORK(KONEIA),WORK(KD1IJT), 450 * WORK(KFD1II),SK1,SK2,SK3,SK4,WORK(KINTIJ), 451 * WORK(KINTAI),WORK(KINTIA),WORK(KINIJT), 452 * WORK(KINABT),WORK(KFROIJ),WORK(KFROJI), 453 * WORK(KFROAI),WORK(KFROIA),WORK(KFROII), 454 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM,ICON) 455C 456C-------------------------------------------------------- 457C Calculate one-electron contribution to right- 458C hand-side of virtual-frozen multipliers. 459C-------------------------------------------------------- 460C 461 ISYM = 1 462 ICON = 1 463 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 464 ! 465 !eta_aI 466 ! 467 CALL CC_ETAIF(ZKDIA(KOFRES),WORK(KONEAB),WORK(KONEAI), 468 * WORK(KONEIA),WORK(KD1IJT),WORK(KD1ABT), 469 * WORK(KFD1II),SK1,SK2,SK3,SK4,WORK(KINTIJ), 470 * WORK(KINTAB),WORK(KINTAI),WORK(KINTIA), 471 * WORK(KINABT),WORK(KFROIJ),WORK(KFROJI), 472 * WORK(KFROAI),WORK(KFROIA),WORK(KFROII), 473 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM,ICON) 474 475 ENDIF 476 477C 478C----------------------------------------------------------------------- 479C Transform now one-electron integrals with regular CMO coefficients 480C----------------------------------------------------------------------- 481C 482 IF (IMODEL.EQ.1) THEN 483 ISYM = 1 484 CALL CCDINTMO(WORK(KINTIJ),WORK(KINTIA),WORK(KINTAB), 485 * WORK(KINTAI),WORK(KONINT),WORK(KCMO), 486 * WORK(KCMO),WORK(KEND1),LWRK1,ISYM) 487C 488 IF (LOCDBG) THEN 489 490 HIJNO = DDOT(NMATIJ(1),WORK(KINTIJ),1,WORK(KINTIJ),1) 491 HAINO = DDOT(NT1AMX,WORK(KINTAI),1,WORK(KINTAI),1) 492 HIANO = DDOT(NT1AMX,WORK(KINTIA),1,WORK(KINTIA),1) 493 HABNO = DDOT(NMATAB(1),WORK(KINTAB),1,WORK(KINTAB),1) 494 WRITE(LUPRI,*) 'CC_DEN_PTFC: IOPT = ', IOPT 495 WRITE(LUPRI,*) 496 & 'Norm of CMO one electron integrals in MO-basis:' 497 WRITE(LUPRI,*) HIJNO, HAINO, HIANO, HABNO 498 499 ENDIF 500C 501C-------------------------------------------------- 502C Read in D^(T)_ia from file, contract with h_pq 503C and add to result 504C-------------------------------------------------- 505 506 KONEIA_PT = KEND1 507 KEND1 = KONEIA_PT + NT1AMX 508 LWRK1 = LWORK - KEND1 509C 510 LUPTIA = -1 511 FNDPTIA = 'DPTIA' 512 CALL WOPEN2(LUPTIA,FNDPTIA,64,0) 513C 514 IOFF = 1 515 CALL GETWA2(LUPTIA,FNDPTIA,WORK(KONEIA_PT),IOFF,NT1AMX) 516 CALL WCLOSE2(LUPTIA,FNDPTIA,'KEEP') 517 518 if (locdbg) then 519 write(lupri,*)'Norm of (T) D_ia:', 520 & ddot(NT1AMX,WORK(KONEIA_PT),1,WORK(KONEIA_PT),1) 521 end if 522 523 IF (LTSTE) THEN 524 ! 525 !D(T)_ia contribution to the energy: sum_ia D_ia h_ia 526 ! 527 EN1PT = DDOT(NT1AMX,WORK(KONEIA_PT),1,WORK(KINTIA),1) 528 529 END IF 530 531C 532C------------------------------------------------------------ 533C Add one electron contribution to Zeta-kappa-0. 534C from the (T) one electron density... 535C------------------------------------------------------------ 536C 537 ISYM = 1 538 539 IF (IOPT.EQ.2) THEN 540 541 KOFFIJ = 1 542 KOFFAB = 1 + NMATIJ(1) 543 CALL CCPT_ETARS_1E(ZKDIA(KOFFIJ),ZKDIA(KOFFAB), 544 & WORK(KINTIJ),WORK(KINTAI), 545 & WORK(KINTIA),WORK(KINTAB),WORK(KONEIA_PT), 546 & WORK(KEND1),LWRK1,ISYM) 547 548 549 END IF 550C 551C------------------------------------------------------------ 552C Calculate (T) contribution to ETA_ai and transfer outside 553C------------------------------------------------------------ 554C 555 CALL CCPT_ETAAI_1E(ETAAI,WORK(KINTIJ),WORK(KINTAI), 556 & WORK(KINTIA),WORK(KINTAB), 557 & WORK(KONEIA_PT),WORK(KEND1), 558 & LWRK1,ISYM) 559 560 LUNGO = NMATIJ(1) + NMATAB(1) + 2*NT1AMX 561 * + 2*NCOFRO(1) + 2*NT1FRO(1) 562 XZKDIA = DDOT(LUNGO,ZKDIA,1,ZKDIA,1) 563 WRITE(LUPRI,*) 564 & 'CC_DEN_PTFC: ZKDIA before TWO electron part', XZKDIA 565 566 IF (FROIMP) THEN 567C 568 ISYM = 1 569 !obtain integrals with frozen indices on regular MO basis 570 ! h_Ij h_jI h_aJ h_Ia h_IJ 571 ! 572 CALL CCIFROMO(WORK(KFROIJ),WORK(KFROJI),WORK(KFROAI), 573 * WORK(KFROIA),WORK(KFROII),WORK(KONINT), 574 * WORK(KCMO),WORK(KCMO),WORK(KCMOF), 575 * WORK(KEND1),LWRK1,ISYM) 576 IF (LOCDBG) THEN 577 578 HIJFN = DDOT(NCOFRO(1),WORK(KFROIJ),1,WORK(KFROIJ),1) 579 HJIFN = DDOT(NCOFRO(1),WORK(KFROJI),1,WORK(KFROJI),1) 580 HAIFN = DDOT(NT1FRO(1),WORK(KFROAI),1,WORK(KFROAI),1) 581 HIAFN = DDOT(NT1FRO(1),WORK(KFROIA),1,WORK(KFROIA),1) 582 WRITE(LUPRI,*) 'CC_DEN_PTFC: IOPT = ', IOPT 583 WRITE(LUPRI,*) 584 & 'Norm of CMO one ele fro integrals in MO-basis:' 585 WRITE(LUPRI,*) HIJFN, HJIFN,HAIFN, HIAFN 586 587 ENDIF 588 589 KOFRE1 = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 590 KOFRE2 = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 591 592 CALL CC_ETFRO1E(ZKDIA(KOFRE1),ZKDIA(KOFRE2), 593 & WORK(KONEIA_PT), 594 & WORK(KFROIJ),WORK(KFROJI), 595 & WORK(KFROAI),WORK(KFROIA), 596 & WORK(KEND1),LWRK1,ISYM) 597C 598 END IF 599 END IF 600 601C 602C-------------------------------------------------------- 603C Calculate M-intermediate of zeta- and t-amplitudes. 604C-------------------------------------------------------- 605C 606 CALL CC_MI(WORK(KMINT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 607 * WORK(KEND1),LWRK1) 608C 609C-------------------------------------------------------- 610C Calculate resorted M-intermediate M(imjk)->M(mkij). 611C-------------------------------------------------------- 612C 613 CALL CC_MIRS(WORK(KMIRES),WORK(KMINT)) 614C 615 TIMONE = SECOND() - TIMONE 616C 617C-------------------------------------------- 618C Start the loop over integrals. 619C Salva tutto quanto definito fino ad ora 620C-------------------------------------------- 621C 622 KENDS2 = KEND1 623 LWRKS2 = LWRK1 624C 625 IF (DIRECT) THEN 626 IF (HERDIR) THEN 627 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 628 ELSE 629 KCCFB1 = KEND1 630 KINDXB = KCCFB1 + MXPRIM*MXCONT 631 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 632 LWRK1 = LWORK - KEND1 633 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 634 * KODPP1,KODPP2,KRDPP1,KRDPP2, 635 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 636 * WORK(KEND1),LWRK1,IPRERI) 637 KEND1 = KFREE 638 LWRK1 = LFREE 639 ENDIF 640 NTOSYM = 1 641 ELSE 642 NTOSYM = NSYM 643 ENDIF 644C 645 KENDSV = KEND1 646 LWRKSV = LWRK1 647C 648 ICDEL1 = 0 649 DO 100 ISYMD1 = 1,NTOSYM 650C 651 IF (DIRECT) THEN 652 IF (HERDIR) THEN 653 NTOT = MAXSHL 654 ELSE 655 NTOT = MXCALL 656 ENDIF 657 ELSE 658 NTOT = NBAS(ISYMD1) 659 ENDIF 660C 661 DO 110 ILLL = 1,NTOT 662C 663C--------------------------------------------- 664C If direct calculate the integrals. 665C--------------------------------------------- 666C 667 IF (DIRECT) THEN 668C 669 KEND1 = KENDSV 670 LWRK1 = LWRKSV 671C 672 DTIME = SECOND() 673 IF (HERDIR) THEN 674 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 675 & IPRINT) 676 ELSE 677 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 678 * WORK(KODCL1),WORK(KODCL2), 679 * WORK(KODBC1),WORK(KODBC2), 680 * WORK(KRDBC1),WORK(KRDBC2), 681 * WORK(KODPP1),WORK(KODPP2), 682 * WORK(KRDPP1),WORK(KRDPP2), 683 * WORK(KCCFB1),WORK(KINDXB), 684 * WORK(KEND1), LWRK1,IPRERI) 685 ENDIF 686 DTIME = SECOND() - DTIME 687 TIMHE2 = TIMHE2 + DTIME 688C 689 KRECNR = KEND1 690 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 691 LWRK1 = LWORK - KEND1 692 IF (LWRK1 .LT. 0) THEN 693 CALL QUIT('Insufficient core in CC_DEN_PTFC') 694 END IF 695C 696 ELSE 697 NUMDIS = 1 698 ENDIF 699C 700C----------------------------------------------------- 701C Loop over number of distributions in disk. 702C----------------------------------------------------- 703C 704 DO 120 IDEL2 = 1,NUMDIS 705C 706 IF (DIRECT) THEN 707 IDEL = INDEXA(IDEL2) 708 IF (NOAUXB) THEN 709 IDUM = 1 710 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 711 END IF 712 ISYMD = ISAO(IDEL) 713 ELSE 714 IDEL = IBAS(ISYMD1) + ILLL 715 ISYMD = ISYMD1 716 ENDIF 717C 718C---------------------------------------- 719C Work space allocation two. 720C---------------------------------------- 721C 722 ISYDEN = ISYMD 723C 724 KD2IJG = KEND1 725 KD2AIG = KD2IJG + ND2IJG(ISYDEN) 726 KD2IAG = KD2AIG + ND2AIG(ISYDEN) 727 KD2ABG = KD2IAG + ND2AIG(ISYDEN) 728 KEND2 = KD2ABG + ND2ABG(ISYDEN) 729 LWRK2 = LWORK - KEND2 730C 731 IF (LWRK2 .LT. 0) THEN 732 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 733 CALL QUIT('Insufficient space for allocation '// 734 & '2 in CC_DEN_PTFC') 735 ENDIF 736C 737C------------------------------------------------------- 738C Initialize 4 two electron density arrays. 739C------------------------------------------------------- 740C 741 CALL DZERO(WORK(KD2IJG),ND2IJG(ISYDEN)) 742 CALL DZERO(WORK(KD2AIG),ND2AIG(ISYDEN)) 743 CALL DZERO(WORK(KD2IAG),ND2AIG(ISYDEN)) 744 CALL DZERO(WORK(KD2ABG),ND2ABG(ISYDEN)) 745C 746C----------------------------------------------------------------------------- 747C Calculate the CCSD-like two electron density d(pq,gamma;delta). 748C----------------------------------------------------------------------------- 749C 750 AUTIME = SECOND() 751C 752 CALL CC_DEN2(WORK(KD2IJG),WORK(KD2AIG),WORK(KD2IAG), 753 * WORK(KD2ABG),WORK(KZ2AM),WORK(KT2AM), 754 * WORK(KT2AMT),WORK(KMINT),WORK(KXMAT), 755 * WORK(KYMAT),WORK(KONEAB),WORK(KONEAI), 756 * WORK(KONEIA),WORK(KMIRES),WORK(KLAMDH),1, 757 * WORK(KLAMDP),1,WORK(KEND2),LWRK2,IDEL,ISYMD) 758C 759 AUTIME = SECOND() - AUTIME 760 TIMDEN = TIMDEN + AUTIME 761C 762C------------------------------------------ 763C Work space allocation three. 764C------------------------------------------ 765C 766 ISYDIS = MULD2H(ISYMD,ISYMOP) 767C 768 KXINT = KEND2 769 KEND3 = KXINT + NDISAO(ISYDIS) 770 LWRK3 = LWORK - KEND3 771C 772 IF (LWRK3 .LT. 0) THEN 773 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 774 CALL QUIT('Insufficient space for allocation '// 775 & '3 in CC_DEN_PTFC') 776 ENDIF 777 778C 779C-------------------------------------------- 780C Read AO integral distribution. 781C-------------------------------------------- 782C 783 AUTIME = SECOND() 784 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 785 * WORK(KRECNR),DIRECT) 786 AUTIME = SECOND() - AUTIME 787 TIRDAO = TIRDAO + AUTIME 788C 789C---------------------------------------------------------------------- 790C Calculate integral intermediates needed for frozen core. 791C---------------------------------------------------------------------- 792C 793 IF (FROIMP) THEN 794 795 KDSRHF = KEND3 796 KOFOIN = KDSRHF + NDSRHF(ISYMD) 797 KDSFRO = KOFOIN + NOFROO(ISYDIS) 798 KSCRAI = KDSFRO + NDSFRO(ISYDIS) 799 KSCAIF = KSCRAI + NOFROO(ISYDIS) 800 KEND3 = KSCAIF + NCOFRF(ISYDIS) 801 LWRK3 = LWORK - KEND3 802C 803 IF (LWRK3 .LT. 0) THEN 804 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 805 CALL QUIT('Insufficient space for allocation '// 806 & 'in CC_DEN_PTFC') 807 ENDIF 808C 809 CALL DZERO(WORK(KSCRAI),NOFROO(ISYDIS)) 810 CALL DZERO(WORK(KSCAIF),NCOFRF(ISYDIS)) 811C 812C------------------------------------------------------------------------- 813C Transform one index in the integral batch to correlated. 814C------------------------------------------------------------------------- 815C 816 !(alp-bet,i,delta) 817 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMDP), 818 * ISYMOP,WORK(KEND3),LWRK3,ISYDIS) 819C 820C--------------------------------------------------------------------- 821C Transform one index in the integral batch to frozen. 822C--------------------------------------------------------------------- 823C 824 !(alp-bet,i,delta) 825 CALL CC_GTOFRO(WORK(KXINT),WORK(KDSFRO),WORK(KCMOF), 826 * WORK(KEND3),LWRK3,ISYDIS) 827C 828C-------------------------------------------------------------- 829C Calculate integral batch (cor fro | cor del). 830C-------------------------------------------------------------- 831C 832 CALL CC_OFROIN(WORK(KDSRHF),WORK(KOFOIN),WORK(KCMOF), 833 * WORK(KEND3),LWRK3,ISYDIS) 834 835C 836C------------------------------------------------------------------ 837C Calculate terms to ai-part from KOFOIN integrals. 838C------------------------------------------------------------------ 839C 840 CALL CC_FRCOC1(WORK(KSCRAI),WORK(KOFOIN),ISYDIS) 841 842C 843C---------------------------------------------------------------- 844C Calculate exchange parts with KDSFRO integrals. 845C---------------------------------------------------------------- 846C 847 CALL CC_FRCOMI(WORK(KSCRAI),WORK(KSCAIF), 848 * WORK(KDSFRO),WORK(KCMOF), 849 * WORK(KEND3),LWRK3,ISYMD) 850 851C 852C---------------------------------------------------- 853C Calculate coulomb part of aI block. 854C---------------------------------------------------- 855C 856 CALL CC_FRCOF1(WORK(KSCAIF),WORK(KDSFRO),WORK(KCMOF), 857 * WORK(KEND3),LWRK3,ISYMD) 858 859C 860C----------------------------------------------------- 861C Calculate exchange part of aI block. 862C----------------------------------------------------- 863C 864 CALL CC_FRCOF2(WORK(KSCAIF),WORK(KDSRHF),WORK(KCMOF), 865 * WORK(KEND3),LWRK3,ISYMD) 866 867C 868C---------------------------------------------------------- 869C Dump intermediates to disc for later use. 870C---------------------------------------------------------- 871C 872 CALL CC_FRINDU(WORK(KSCRAI),WORK(KSCAIF),IDEL,ISYMD, 873 & LUN1,LUN2) 874 875 ENDIF 876C 877C------------------------------------------------------ 878C Start loop over second AO-index (gamma). 879C------------------------------------------------------ 880C 881 AUTIME = SECOND() 882C 883 DO 130 ISYMG = 1,NSYM 884C 885 ISYMPQ = MULD2H(ISYMG,ISYDEN) 886C 887 DO 140 G = 1,NBAS(ISYMG) 888C 889 KD2GIJ = KD2IJG + ID2IJG(ISYMPQ,ISYMG) 890 * + NMATIJ(ISYMPQ)*(G - 1) 891 KD2GAI = KD2AIG + ID2AIG(ISYMPQ,ISYMG) 892 * + NT1AM(ISYMPQ)*(G - 1) 893 KD2GAB = KD2ABG + ID2ABG(ISYMPQ,ISYMG) 894 * + NMATAB(ISYMPQ)*(G - 1) 895 KD2GIA = KD2IAG + ID2AIG(ISYMPQ,ISYMG) 896 * + NT1AM(ISYMPQ)*(G - 1) 897C 898C------------------------------------------------------------------------- 899C Work space allocation four. 900C Note: d2aob is only used for the ltest case!!!!!!!!!! 901C------------------------------------------------------------------------- 902C 903 KINTAO = KEND3 904 KD2AOB = KINTAO + N2BST(ISYMPQ) 905 KEND4 = KD2AOB + N2BST(ISYMPQ) 906 LWRK4 = LWORK - KEND4 907C 908 IF (LWRK4 .LT. 0) THEN 909 WRITE(LUPRI,*) 'Available:', LWORK 910 WRITE(LUPRI,*) 'Needed:', KEND4 911 CALL QUIT('Insufficient space in CC_DEN_PTFC') 912 ENDIF 913C 914 CALL DZERO(WORK(KD2AOB),N2BST(ISYMPQ)) 915C 916C------------------------------------------------------------- 917C Calculate frozen core contributions to d. 918C------------------------------------------------------------- 919C 920 IF (FROIMP) THEN 921C 922 KFD2IJ = KEND4 923 KFD2JI = KFD2IJ + NCOFRO(ISYMPQ) 924 KFD2AI = KFD2JI + NCOFRO(ISYMPQ) 925 KFD2IA = KFD2AI + NT1FRO(ISYMPQ) 926 KFD2II = KFD2IA + NT1FRO(ISYMPQ) 927 KEND4 = KFD2II + NFROFR(ISYMPQ) 928 LWRK4 = LWORK - KEND4 929C 930 IF (LWRK4 .LT. 0) THEN 931 WRITE (LUPRI,*) 'Available:', LWORK 932 WRITE (LUPRI,*) 'Needed:', KEND4 933 CALL QUIT( 934 * 'Insufficient work space in CC_DEN_PTFC') 935 ENDIF 936C 937 CALL DZERO(WORK(KFD2IJ),NCOFRO(ISYMPQ)) 938 CALL DZERO(WORK(KFD2JI),NCOFRO(ISYMPQ)) 939 CALL DZERO(WORK(KFD2AI),NT1FRO(ISYMPQ)) 940 CALL DZERO(WORK(KFD2IA),NT1FRO(ISYMPQ)) 941 CALL DZERO(WORK(KFD2II),NFROFR(ISYMPQ)) 942C 943C To calculate the contributions to d(pq,gam;del) 944C where at least one of the two indices p & q is frozen. 945C 946 CALL CC_FD2BL(WORK(KFD2II),WORK(KFD2IJ), 947 * WORK(KFD2JI),WORK(KFD2AI), 948 * WORK(KFD2IA),WORK(KONEIJ), 949 * WORK(KONEAB),WORK(KONEAI), 950 * WORK(KONEIA),WORK(KCMOF), 951 * WORK(KLAMDH),WORK(KLAMDP), 952 * WORK(KEND4),LWRK4,IDEL, 953 * ISYMD,G,ISYMG) 954C 955C ! calculate the contributions to D2AO from d(pq,gam;del) 956C ! where at least one of the two indices p & q is frozen 957 958 CALL CC_FD2AO(WORK(KD2AOB),WORK(KFD2II), 959 * WORK(KFD2IJ),WORK(KFD2JI), 960 * WORK(KFD2AI),WORK(KFD2IA), 961 * WORK(KCMOF),WORK(KLAMDH), 962 * WORK(KLAMDP),WORK(KEND4),LWRK4, 963 * ISYMPQ) 964C 965C Purpose: To calculate the contributions to d(pq,gam;del) where 966C gamma has been backtransformed from a frozen index. 967C 968 CALL CC_D2GAF(WORK(KD2GIJ),WORK(KD2GAB), 969 * WORK(KD2GAI),WORK(KD2GIA), 970 * WORK(KONEIJ),WORK(KONEAB), 971 * WORK(KONEAI),WORK(KONEIA), 972 * WORK(KCMOF),IDEL,ISYMD,G,ISYMG) 973C 974 ENDIF 975C 976C------------------------------------------------------- 977C Square up AO-integral distribution. 978C------------------------------------------------------- 979C 980 KOFFIN = KXINT + IDSAOG(ISYMG,ISYDIS) 981 * + NNBST(ISYMPQ)*(G - 1) 982C 983 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ, 984 * WORK(KINTAO)) 985C 986C--------------------------------------------------------------------------- 987C If energy test backtransform density fully to AO basis. 988C--------------------------------------------------------------------------- 989C 990 IF (LTSTE) THEN 991C 992 CALL CC_DENAO(WORK(KD2AOB),ISYMPQ, 993 * WORK(KD2GAI),WORK(KD2GAB), 994 * WORK(KD2GIJ),WORK(KD2GIA),ISYMPQ, 995 * WORK(KLAMDP),1,WORK(KLAMDH),1, 996 * WORK(KEND4),LWRK4) 997C 998C--------------------------------------------------------------------- 999C Add relaxation terms to set up effective density. 1000C--------------------------------------------------------------------- 1001C 1002! IF (IOPT .EQ. 3) THEN 1003C 1004! ICON = 1 1005! CALL CC_D2EFF(WORK(KD2AOB),G,ISYMG,IDEL, 1006! * ISYMD,WORK(KKABAO), 1007! * WORK(KDHFAO),ICON) 1008C 1009! ENDIF 1010C 1011C---------------------------------------------------------------------- 1012C Calculate the 2 e- density contribution to E-ccsd. 1013C---------------------------------------------------------------------- 1014C 1015 ECCSD2 = ECCSD2 + HALF*DDOT(N2BST(ISYMPQ), 1016 * WORK(KD2AOB),1,WORK(KINTAO),1) 1017C 1018 end if 1019C 1020C----------------------------------------------- 1021C Work space allocation five. 1022C----------------------------------------------- 1023C 1024 KIJINT = KEND4 1025 KAIINT = KIJINT + NMATIJ(ISYMPQ) 1026 KIAINT = KAIINT + NT1AM(ISYMPQ) 1027 KABINT = KIAINT + NT1AM(ISYMPQ) 1028 KABTIN = KABINT + NMATAB(ISYMPQ) 1029 KIJTIN = KABTIN + NMATAB(ISYMPQ) 1030 KD2TAB = KIJTIN + NMATIJ(ISYMPQ) 1031 KD2TIJ = KD2TAB + NMATAB(ISYMPQ) 1032 KEND5 = KD2TIJ + NMATIJ(ISYMPQ) 1033 LWRK5 = LWORK - KEND5 1034 IF (FROIMP) THEN 1035 KIIFRO = KEND5 1036 KIJFRO = KIIFRO + NFROFR(ISYMPQ) 1037 KJIFRO = KIJFRO + NCOFRO(ISYMPQ) 1038 KAIFRO = KJIFRO + NCOFRO(ISYMPQ) 1039 KIAFRO = KAIFRO + NT1FRO(ISYMPQ) 1040 KEND5 = KIAFRO + NT1FRO(ISYMPQ) 1041 LWRK5 = LWORK - KEND5 1042 ENDIF 1043C 1044 IF (LWRK5 .LT. 0) THEN 1045 WRITE(LUPRI,*) 'Available:', LWORK 1046 WRITE(LUPRI,*) 'Needed:', KEND5 1047 CALL QUIT('Insufficient work space '// 1048 & 'in CC_DEN_PTFC') 1049 ENDIF 1050C 1051C---------------------------------------------------------------- 1052C Transform 2 integral indices to MO-basis. 1053C---------------------------------------------------------------- 1054C 1055 ISYM = ISYMPQ 1056 CALL CCDINTMO(WORK(KIJINT),WORK(KIAINT), 1057 * WORK(KABINT),WORK(KAIINT), 1058 * WORK(KINTAO),WORK(KLAMDP), 1059 * WORK(KLAMDH),WORK(KEND5), 1060 * LWRK5,ISYM) 1061C 1062 IF (FROIMP) THEN 1063C 1064C Prepare integrals g_pq (gam,del) where one index is frozen 1065C 1066 ISYM = ISYMPQ 1067 CALL CCIFROMO(WORK(KIJFRO),WORK(KJIFRO), 1068 * WORK(KAIFRO),WORK(KIAFRO), 1069 * WORK(KIIFRO),WORK(KINTAO), 1070 * WORK(KLAMDP),WORK(KLAMDH), 1071 * WORK(KCMOF),WORK(KEND5), 1072 * LWRK5,ISYM) 1073C 1074 ENDIF 1075C 1076C----------------------------------------------------------------- 1077C Set up transposed integrals and densities. 1078C----------------------------------------------------------------- 1079C 1080 CALL DCOPY(NMATAB(ISYMPQ),WORK(KABINT),1, 1081 * WORK(KABTIN),1) 1082 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KIJINT),1, 1083 * WORK(KIJTIN),1) 1084 CALL DCOPY(NMATAB(ISYMPQ),WORK(KD2GAB),1, 1085 * WORK(KD2TAB),1) 1086 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KD2GIJ),1, 1087 * WORK(KD2TIJ),1) 1088C 1089 CALL CC_EITR(WORK(KABTIN),WORK(KIJTIN), 1090 * WORK(KEND5),LWRK5,ISYMPQ) 1091 CALL CC_EITR(WORK(KD2TAB),WORK(KD2TIJ), 1092 * WORK(KEND5),LWRK5,ISYMPQ) 1093C 1094C------------------------------------------------------------------- 1095C Calculate 2 e- contribution to Zeta-Kappa-0. 1096C------------------------------------------------------------------- 1097C 1098 ISYM = ISYMPQ 1099 IF (IOPT.EQ.2) THEN 1100 1101 KOFFIJ = 1 1102 CALL CC2_ETIJ(ZKDIA(KOFFIJ), 1103 & WORK(KIJINT),WORK(KAIINT), 1104 & WORK(KIAINT),WORK(KABINT), 1105 & WORK(KD2GIJ),WORK(KD2GAI), 1106 & WORK(KD2GIA),WORK(KD2GAB), 1107 & WORK(KT1AM),WORK(KEND5),LWRK5, 1108 & ISYM) 1109 1110 KOFFAB = NMATIJ(1) + 1 1111 CALL CC2_ETAB(ZKDIA(KOFFAB), 1112 & WORK(KIJINT),WORK(KAIINT), 1113 * WORK(KIAINT),WORK(KABINT), 1114 * WORK(KD2GIJ),WORK(KD2GAI), 1115 * WORK(KD2GIA),WORK(KD2GAB), 1116 * WORK(KT1AM),WORK(KEND5),LWRK5, 1117 * ISYM) 1118 END IF 1119 1120 CALL CCDENZK0(ETAAI,WORK(KIJINT),WORK(KAIINT), 1121 * WORK(KIAINT),WORK(KABINT), 1122 * WORK(KD2GIJ),WORK(KD2GAI), 1123 * WORK(KD2GIA),WORK(KD2GAB), 1124 * WORK(KIJTIN),WORK(KABTIN), 1125 * WORK(KD2TIJ),WORK(KD2TAB), 1126 * WORK(KT1AM),WORK(KEND5),LWRK5, 1127 * ISYM) 1128C 1129 IF (FROIMP) THEN 1130C 1131 ISYM = ISYMPQ 1132 ! 1133 ! contributions to eta_ai from loop over frozen 1134 ! 1135 CALL CCFRETAI(ETAAI,WORK(KIJFRO), 1136 * WORK(KJIFRO),WORK(KAIFRO), 1137 * WORK(KIAFRO),WORK(KFD2IJ), 1138 * WORK(KFD2JI),WORK(KFD2AI), 1139 * WORK(KFD2IA),WORK(KT1AM), 1140 * WORK(KEND5),LWRK5,ISYM) 1141C 1142C----------------------------------------------------------------------- 1143C Calculate two-electron contribution to right- 1144C hand-side of correlated-frozen multipliers. 1145C----------------------------------------------------------------------- 1146C 1147 ICON = 2 1148 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 1149 * + 2*NT1FRO(1) + 1 1150 ! 1151 ! eta_iJ 1152 ! 1153 CALL CC_ETIJF(ZKDIA(KOFRES),WORK(KD2GIJ), 1154 * WORK(KD2GAB),WORK(KD2GAI), 1155 * WORK(KD2GIA),WORK(KD2TIJ), 1156 * WORK(KFD2II),WORK(KFD2IJ), 1157 * WORK(KFD2JI),WORK(KFD2AI), 1158 * WORK(KFD2IA),WORK(KIJINT), 1159 * WORK(KAIINT),WORK(KIAINT), 1160 * WORK(KIJTIN),WORK(KABTIN), 1161 * WORK(KIJFRO),WORK(KJIFRO), 1162 * WORK(KAIFRO),WORK(KIAFRO), 1163 * WORK(KIIFRO),WORK(KT1AM), 1164 * WORK(KEND5),LWRK5,ISYM,ICON) 1165C 1166C----------------------------------------------------------------------- 1167C Calculate two-electron contribution to right- 1168C hand-side of virtual-frozen multipliers. 1169C----------------------------------------------------------------------- 1170C 1171 ICON = 2 1172 KOFRE = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1173 ! 1174 ! eta_aI 1175 ! 1176 CALL CC_ETAIF(ZKDIA(KOFRE),WORK(KD2GAB), 1177 * WORK(KD2GAI),WORK(KD2GIA), 1178 * WORK(KD2TIJ),WORK(KD2TAB), 1179 * WORK(KFD2II),WORK(KFD2IJ), 1180 * WORK(KFD2JI),WORK(KFD2AI), 1181 * WORK(KFD2IA),WORK(KIJINT), 1182 * WORK(KABINT),WORK(KAIINT), 1183 * WORK(KIAINT),WORK(KABTIN), 1184 * WORK(KIJFRO),WORK(KJIFRO), 1185 * WORK(KAIFRO),WORK(KIAFRO), 1186 * WORK(KIIFRO),WORK(KT1AM), 1187 * WORK(KEND5),LWRK5,ISYM,ICON) 1188C 1189C----------------------------------------------------------------------- 1190C Calculate two-electron contribution to right- 1191C hand-side of frozen-frozen multipliers. 1192C----------------------------------------------------------------------- 1193C 1194 ICON = 2 1195 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 1196 * + 2*NT1FRO(1) + 2*NCOFRO(1) + 1 1197c ! 1198c ! eta_ab contribs from loop over frozen I 1199c ! 1200 CALL CCFRETAB(ZKDIA,WORK(KIJFRO), 1201 * WORK(KJIFRO),WORK(KAIFRO), 1202 * WORK(KIAFRO),WORK(KFD2IJ), 1203 * WORK(KFD2JI),WORK(KFD2AI), 1204 * WORK(KFD2IA),WORK(KT1AM), 1205 * WORK(KEND5),LWRK5,ISYM) 1206 ! 1207 ! eta_ij contribs from loop over frozen L 1208 ! 1209 CALL CCFRETIJ(ZKDIA,WORK(KIJFRO), 1210 * WORK(KJIFRO),WORK(KAIFRO), 1211 * WORK(KIAFRO),WORK(KFD2IJ), 1212 * WORK(KFD2JI),WORK(KFD2AI), 1213 * WORK(KFD2IA),WORK(KT1AM), 1214 * WORK(KEND5),LWRK5,ISYM) 1215c 1216C 1217 ENDIF !froimp 1218! END IF !ltste 1219C 1220 140 CONTINUE 1221 130 CONTINUE 1222C 1223 AUTIME = SECOND() - AUTIME 1224 TIMRES = TIMRES + AUTIME 1225C 1226 120 CONTINUE 1227 110 CONTINUE 1228 100 CONTINUE 1229 1230C 1231C----------------------- 1232C Regain work space. 1233C----------------------- 1234C 1235 KEND1 = KENDS2 1236 LWRK1 = LWRKS2 1237C 1238C-------------------------------------------- 1239C Add contribution from 2-e (T) densities 1240C-------------------------------------------- 1241C 1242 if (locdbg) then 1243 KOFFIJ = 1 1244 KOFFAB = 1 + NMATIJ(1) 1245 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1246 KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1247 write(lupri,*) ' ' 1248 write(lupri,*) 'Before call to CC_DEN_PT2' 1249 xtest=ddot(nmatij(1),zkdia(koffij),1,zkdia(koffij),1) 1250 write(lupri,*) 'Norm of ETAIJ: ', xtest 1251 xtest=ddot(2*NCOFRO(1),zkdia(kofifj),1,zkdia(kofifj),1) 1252 write(lupri,*) 'Norm of ETIFJ: ', xtest 1253 xtest=ddot(nmatab(1),zkdia(koffab),1,zkdia(koffab),1) 1254 write(lupri,*) 'Norm of ETAAB: ', xtest 1255 xtest=ddot(2*nt1amx,etaai(1),1,etaai(1),1) 1256 write(lupri,*) 'Norm of ETAAI: ', xtest 1257 xtest=ddot(2*nt1fro(1),zkdia(kofafi),1,zkdia(kofafi),1) 1258 write(lupri,*) 'Norm of ETAFI: ', xtest 1259 end if 1260 1261 IF (IMODEL.EQ.1) THEN 1262 1263 IOPT1 = 2 1264 en2pt = zero 1265 1266 if (.false.) then 1267 KOFFIJ = 1 1268 KOFFAB = 1 + NMATIJ(1) 1269 CALL CC_DEN_PT2(ETAAI,ZKDIA(KOFFIJ), 1270 & ZKDIA(KOFFAB),WORK(KEND1),LWRK1,IOPT1, 1271 & ltste,en2pt) 1272 else 1273 1274 CALL CCPT_DEN2FC(ETAAI,ZKDIA, 1275 & WORK(KONEIA_PT), 1276 & WORK(KEND1),LWRK1, 1277 & IOPT1,LTSTE,EN2PT) 1278 end if 1279 1280 END IF 1281 if (locdbg) then 1282 KOFFIJ = 1 1283 KOFFAB = 1 + NMATIJ(1) 1284 KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1285 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1286 write(lupri,*) ' ' 1287 write(lupri,*) 'After call to CC_DEN_PT2' 1288 xtest=ddot(nmatij(1),zkdia(koffij),1,zkdia(koffij),1) 1289 write(lupri,*) 'Norm of ETAIJ: ', xtest 1290 xtest=ddot(NCOFRO(1),zkdia(kofifj),1,zkdia(kofifj),1) 1291 write(lupri,*) 'Norm of ETIFJ: ', xtest 1292 xtest=ddot(nmatab(1),zkdia(koffab),1,zkdia(koffab),1) 1293 write(lupri,*) 'Norm of ETAAB: ', xtest 1294 xtest=ddot(nt1amx,etaai(1),1,etaai(1),1) 1295 write(lupri,*) 'Norm of ETAAI: ', xtest 1296 xtest=ddot(nt1fro(1),zkdia(kofafi),1,zkdia(kofafi),1) 1297 write(lupri,*) 'Norm of ETAFI: ', xtest 1298 end if 1299C 1300C------------------------------------------------------------------------ 1301C Calculate the ZK0(ij) and ZK0(ab) blocks (to be used to correct eta_ai) 1302C from eta_ij/e_i-e_j and eta_ab/e_a-e_b 1303C------------------------------------------------------------------------ 1304C 1305 IF (IOPT.EQ.2) THEN 1306 1307 CALL CCSD_ZKBLO(ZKDIA,WORK(KEND1),LWRK1) 1308 1309C 1310C------------------------------------------------------------------------ 1311C Add the diagonal ZK0(ii) and ZK0(aa) elements in the proper places 1312C------------------------------------------------------------------------ 1313C 1314 if (.true.) then 1315 IF (IMODEL.EQ.1) THEN 1316 1317 LUKAPAB = -1 1318 LUKAPIJ = -1 1319 FNKAPAB = 'KAPAB' 1320 FNKAPIJ = 'KAPIJ' 1321 1322 1323 CALL WOPEN2(LUKAPAB,FNKAPAB,64,0) 1324 CALL WOPEN2(LUKAPIJ,FNKAPIJ,64,0) 1325C 1326 KKAPII = KEND1 1327 KEND1 = KKAPII + NRHFT 1328 LWRK1 = LWORK - KEND1 1329 1330 IF (NRHFT .GT. 0) THEN 1331 IOFF = 1 1332 CALL GETWA2(LUKAPIJ,FNKAPIJ,WORK(KKAPII),IOFF,NRHFT) 1333 ENDIF 1334 1335 DO ISYMI = 1, NSYM 1336 DO I = 1, NRHF(ISYMI) 1337 NI = IRHF(ISYMI) + I 1338 NII = IMATIJ(ISYMI,ISYMI) + NRHF(ISYMI)*(I-1) + I 1339 ZKDIA(NII) = WORK(KKAPII+NI-1) 1340 END DO 1341 END DO 1342 1343 KKAPAA = KEND1 1344 KEND1 = KKAPAA + NVIRT 1345 LWRK1 = LWORK - KEND1 1346C 1347 IF (NVIRT .GT. 0) THEN 1348 IOFF = 1 1349 CALL GETWA2(LUKAPAB,FNKAPAB,WORK(KKAPAA),IOFF,NVIRT) 1350 ENDIF 1351 1352! Bemaerk, IVIR(ISYM) initialized to NRHFT 1353 1354 DO ISYMA = 1, NSYM 1355 DO A = 1, NVIR(ISYMA) 1356 NA = IVIR(ISYMA) + A - NRHFT 1357 NAA = IMATAB(ISYMA,ISYMA) + NVIR(ISYMA)*(A-1) + A 1358 ZKDIA(NMATIJ(1) + NAA) = WORK(KKAPAA+NA-1) 1359 END DO 1360 END DO 1361 1362 if (locdbg) then 1363 XZKDIA = DDOT(LUNGO,ZKDIA,1,ZKDIA,1) 1364 WRITE(LUPRI,*) 'CC_DEN_PTFC: ZKDIA after ii, aa ', XZKDIA 1365 1366 XZKDIA = DDOT(LUNGO,ZKDIA,1,ZKDIA,1) 1367 WRITE(LUPRI,*) 'CC_DEN_PTFC: ZKDIA after ii, aa ', XZKDIA 1368 end if 1369 1370 CALL WCLOSE2(LUKAPAB,FNKAPAB,'KEEP') 1371 CALL WCLOSE2(LUKAPIJ,FNKAPIJ,'KEEP') 1372 1373 1374 IF (LTSTE) THEN 1375 1376 !multiply by epsilon_p and sum over p to get the energy 1377 1378 KFOCKDIA = KEND1 1379 KEND1 = KFOCKDIA + NORBTS 1380 LWRK1 = LWORK-KEND1 1381C 1382C------------------------------------- 1383C Read canonical orbital energies. 1384C------------------------------------- 1385C 1386 CALL DZERO(WORK(KFOCKDIA),NORBTS) 1387 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1388 & .FALSE.) 1389 REWIND (LUSIFC) 1390C 1391 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1392 READ (LUSIFC) 1393 READ (LUSIFC) (WORK(KFOCKDIA + I - 1), I = 1,NORBTS) 1394C 1395 CALL GPCLOSE(LUSIFC,'KEEP') 1396C 1397C---------------------------------------------------------------- 1398C Change symmetry ordering of the canonical orbital energies. 1399C---------------------------------------------------------------- 1400C 1401 IF (FROIMP) 1402 * CALL CCSD_DELFRO(WORK(KFOCKDIA),WORK(KEND1),LWRK1) 1403C 1404 CALL FOCK_REORDER(WORK(KFOCKDIA),WORK(KEND1),LWRK1) 1405C 1406C-------------------------------------------- 1407C Calculate sum_p kappabar_pp * epsilon_p 1408C Occupied block: 1409C-------------------------------------------- 1410C 1411 SKAPEP = DDOT(NORBT,WORK(KKAPII),1,WORK(KFOCKDIA),1) 1412 END IF 1413 END IF 1414c 1415 END IF 1416 end if 1417C 1418C------------------------------------------------ 1419C Correct Eta_ai with ZK0(ij) and ZK0(ab) blocks 1420C------------------------------------------------ 1421C 1422 IF ((IOPT.EQ.2).OR.(FROIMP)) THEN 1423 IF (FROIMP) THEN 1424 1425 KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1426 ! 1427 !calculate kappabar_iJ = eta_iJ/e_i-e_J 1428 ! 1429 CALL CC_ZKFCB(ZKDIA(KOFIFJ),WORK(KEND1),LWRK1) 1430 !xnorm = ddot(ncofro(1),zkdia(kofifj),1,zkdia(kofifj),1) 1431 !write(lupri,*) ' Norm of ZKbar_iJ : ', xnorm 1432 END IF 1433 1434 !if (.true.) then 1435 write(lupri,*)'CC_DEN_PTFC using CCETACOR' 1436 CALL CCETACOR(ETAAI,ZKDIA,WORK(KEND1),LWRK1) 1437 !else 1438 ! write(lupri,*)'CC_DEN_PTFC using CCETACORNEW' 1439 ! CALL CCETACORNEW(ETAAI,ZKDIA,WORK(KEND1),LWRK1) 1440 !end if 1441 1442 KOFFIJ = 1 1443 KOFFAB = NMATIJ(1) + 1 1444 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1445 KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1446 1447 !write(lupri,*) 'After call to CCETACORNEW' 1448 !xtest=ddot(nmatij(1),zkdia(koffij),1,zkdia(koffij),1) 1449 !write(lupri,*) 'Norm of ETAIJ: ', xtest 1450 !xtest=ddot(NCOFRO(1),zkdia(kofifj),1,zkdia(kofifj),1) 1451 !write(lupri,*) 'Norm of ETIFJ: ', xtest 1452 !xtest=ddot(nmatab(1),zkdia(koffab),1,zkdia(koffab),1) 1453 !write(lupri,*) 'Norm of ETAAB: ', xtest 1454 !xtest=ddot(nt1amx,etaai(1),1,etaai(1),1) 1455 !write(lupri,*) 'Norm of ETAAI: ', xtest 1456 !xtest=ddot(nt1fro(1),zkdia(kofafi),1,zkdia(kofafi),1) 1457 !write(lupri,*) 'Norm of ETAFI: ', xtest 1458 1459 END IF 1460 1461 XZKDIA = DDOT(LUNGO,ZKDIA,1,ZKDIA,1) 1462 WRITE(LUPRI,*) 'CC_DEN_PTFC: ZKDIA before CCETACOR', XZKDIA 1463C 1464C------------------------------------------------ 1465C Write out frozen block of eta-kappa-bar-0. 1466C eta_iJ 1467C------------------------------------------------ 1468C 1469 IF (((IPRINT .GT. 9).OR.(LOCDBG)) .AND. (FROIMP)) THEN 1470 CALL AROUND('Eta-kappa-bar-0 correlated-frozen block') 1471 KOFRES = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 2*NT1FRO(1) + 1 1472 ZKAPF1 = DDOT(NCOFRO(1),ZKDIA(KOFRES),1, 1473 * ZKDIA(KOFRES),1) 1474 ZKAPFR = ZKAPF1**0.5 1475 WRITE(LUPRI,*) 'Norm of correlated-frozen block:', ZKAPFR 1476 ENDIF 1477C 1478C----------------------------------------------------------- 1479C Calculate the correlated-frozen blocks of kappa-bar-0. 1480C----------------------------------------------------------- 1481C 1482 IF (FROIMP) THEN 1483 1484 !KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1485 ! 1486 !calculate kappabar_iJ = eta_iJ/e_i-e_J 1487 ! 1488 !CALL CC_ZKFCB(ZKDIA(KOFIFJ),WORK(KEND1),LWRK1) 1489 !xnorm = ddot(ncofro(1),zkdia(kofifj),1,zkdia(kofifj),1) 1490 !write(lupri,*) ' Norm of ZKbar_iJ : ', xnorm 1491C 1492C---------------------------------------------------------------- 1493C Calculate correction terms from correlated-frozen block 1494C eta_iJ to both eta_ai and eta_aI 1495C---------------------------------------------------------------- 1496C 1497 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1498 CALL CC_ETACOR(ETAAI,ZKDIA(KOFAFI),ZKDIA(KOFIFJ), 1499 * WORK(KCMOF),LUN1,LUN2,WORK(KEND1),LWRK1) 1500C 1501 ENDIF 1502C 1503C--------------------- 1504C Reorder results. 1505C--------------------- 1506C 1507 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1508 CALL CC_ETARE(ETAAI,ZKDIA(KOFAFI),WORK(KEND1),LWRK1) 1509C 1510C--------------------------------- 1511C Write out eta-ai and eta-aI. 1512C--------------------------------- 1513C 1514 IF ((IPRINT .GT. 20).OR.(LOCDBG)) THEN 1515C 1516 CALL AROUND('Eta-kappa-bar-0-ai vector exiting CC_DEN_PT') 1517C 1518 DO 20 ISYM = 1,NSYM 1519C 1520 WRITE(LUPRI,*) ' ' 1521 WRITE(LUPRI,444) 'Sub-symmetry block number:', ISYM 1522 WRITE(LUPRI,555) '--------------------------' 1523 444 FORMAT(3X,A26,2X,I1) 1524 555 FORMAT(3X,A25) 1525C 1526 KOFF = IALLAI(ISYM,ISYM) + 1 1527 CALL OUTPUT(ETAAI(KOFF),1,NVIR(ISYM),1,NRHFS(ISYM), 1528 * NVIR(ISYM),NRHFS(ISYM),1,LUPRI) 1529C 1530 IF ((NVIR(ISYM) .EQ. 0) .OR. (NRHFS(ISYM) .EQ. 0)) THEN 1531 WRITE(LUPRI,*) 'This sub-symmetry is empty' 1532 ENDIF 1533C 1534 20 CONTINUE 1535 ENDIF 1536C 1537 IF ((IPRINT .GT. 9).OR.(LOCDBG)) THEN 1538 ETAKAN = DDOT(NALLAI(1),ETAAI,1,ETAAI,1) 1539 WRITE(LUPRI,*) 'CC_DEN_PTFC ' 1540 WRITE(LUPRI,*) 'Norm of occupied-virtual block:', ETAKAN 1541 ENDIF 1542C 1543C------------------------------------------- 1544C Write out frozen block of kappa-bar-0. 1545C------------------------------------------- 1546C 1547 IF (((IPRINT .GT. 9).OR.(LOCDBG)) .AND. (FROIMP)) THEN 1548 CALL AROUND('Zeta-kappa-0 correlated-frozen block') 1549 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1550 ZKAPF1 = DDOT(NCOFRO(1),ZKDIA(KOFRES),1, 1551 * ZKDIA(KOFRES),1) 1552 ZKAPFR = ZKAPF1**0.5 1553 WRITE(LUPRI,*) 'Norm of correlated-frozen block:', ZKAPFR 1554 ENDIF 1555C 1556C------------------------------ 1557C Close intermediate files. 1558C------------------------------ 1559C 1560 IF (FROIMP) THEN 1561 CALL WCLOSE2(LUN1,NAME1,'KEEP') 1562 CALL WCLOSE2(LUN2,NAME2,'KEEP') 1563 ENDIF 1564C 1565C----------------------- 1566C Write out timings. 1567C----------------------- 1568C 1569 99 TIMTOT = SECOND() - TIMTOT 1570C 1571 IF (IPRINT .GT. 3) THEN 1572 WRITE (LUPRI,*) ' ' 1573 WRITE (LUPRI,*) 'Requested density dependent '// 1574 & 'quantities calculated' 1575 WRITE (LUPRI,*) 'Total time used in CC_DEN_PTFC:', TIMTOT 1576 ENDIF 1577 IF (IPRINT .GT. 9) THEN 1578 WRITE (LUPRI,*) 'Time used for setting up 2 e- density:',TIMDEN 1579 WRITE (LUPRI,*) 'Time used for contraction with integrals:', 1580 & TIMRES 1581 WRITE (LUPRI,*) 'Time used for reading 2 e- AO-integrals:', 1582 & TIRDAO 1583 WRITE (LUPRI,*) 'Time used for calculating 2 e- AO-integrals:', 1584 * TIMHE2 1585 WRITE (LUPRI,*) 'Time used for 1 e- density & intermediates:', 1586 * TIMONE 1587 ENDIF 1588C 1589C--------------------------------------------------------------- 1590C If energy test add nuclear nuclear repulsion energy and write out E-ccsd. 1591C--------------------------------------------------------------- 1592C 1593 IF (ltste) THEN 1594C 1595 ECCSD = ECCSD1 + ECCSD2 + POTNUC 1596C 1597 WRITE(LUPRI,*) ' ' 1598 WRITE(LUPRI,*) 'Coupled Cluster energy constructed' 1599 WRITE(LUPRI,*) 'from density matrices:' 1600 !IF (CCSD) WRITE(LUPRI,*) 'CCSD-(type) energy:', ECCSD 1601 WRITE(LUPRI,'(A,f15.10)') 'H1 energy, ECCSD1(type) = ',ECCSD1 1602 WRITE(LUPRI,'(A,f15.10)') 'H2 energy, ECCSD2(type) = ',ECCSD2 1603 WRITE(lupri,'(A,f15.10)') 'sum_p e_p kbar_pp = ',SKAPEP 1604 WRITE(lupri,'(A,f15.10)') 'D(T)_ia contrib inc = ',EN1PT 1605 WRITE(lupri,'(A,f15.10)') 'd(T)_pqrs contribution = ',EN2PT 1606 WRITE(LUPRI,'(A,f15.10)') 'Nuc. Pot. energy = ',POTNUC 1607 WRITE(lupri,'(A,f15.10)') 'CCSD(T) energy ? = ', 1608 & ECCSD1+EN1PT+ECCSD2+EN2PT+SKAPEP + POTNUC 1609 !WRITE(LUPRI,*) 'OBS POTNUC is missing!!! ' 1610 1611C 1612 ENDIF 1613C---------------------------------------------------------------------- 1614 CALL QEXIT('CC_DEN_PTFC') 1615 RETURN 1616 END 1617C---------------------------------------------------------------------- 1618C /* Deck ccpt_den2fc */ 1619 SUBROUTINE CCPT_DEN2FC(ETAAI,ZKDIA, 1620 & D1PTIA, 1621 & WORK,LWORK, 1622 & IOPT,LTSTEN,en2pt) 1623C 1624C Written by S. Coriani, based on CC_DEN_PT 1625C January 2002 1626C 1627C Version: 1.0 1628C 1629C Purpose: 1630C drive the calculation of the "pure d_pqrs(T)" contributions to the 1631C ^kappabar-eta_pq RHS of the orbital multipliers 1632C LTSTEN = true, test densities via energy calculation 1633C 1634#include "implicit.h" 1635#include "priunit.h" 1636#include "dummy.h" 1637#include "maxorb.h" 1638#include "maxash.h" 1639#include "mxcent.h" 1640#include "aovec.h" 1641#include "iratdef.h" 1642 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 1643 PARAMETER (FOUR = 4.0D0) 1644 DIMENSION INDEXA(MXCORB_CC) 1645 DIMENSION ETAAI(*), ZKDIA(*), WORK(LWORK) 1646 DIMENSION D1PTIA(*) 1647 LOGICAL LTSTEN 1648#include "ccorb.h" 1649CCN #include "infind.h" not compatible with R12! 1650#include "ccisao.h" 1651#include "r12int.h" 1652#include "inftap.h" 1653#include "blocks.h" 1654#include "ccfield.h" 1655#include "ccsdinp.h" 1656#include "ccinftap.h" 1657#include "ccsdsym.h" 1658#include "ccsdio.h" 1659#include "distcl.h" 1660#include "cbieri.h" 1661#include "eritap.h" 1662#include "ccfro.h" 1663C 1664 CHARACTER MODEL*10 1665 CHARACTER NAME1*8 1666 CHARACTER NAME2*8 1667 1668 LOGICAL LETAB, LETIJ, LETAI 1669 LOGICAL LOCDBG 1670 PARAMETER (LOCDBG=.FALSE.) 1671C 1672 CALL QENTER('CCPT_DEN2FC') 1673C 1674 CALL HEADER('Construct part of rhs for CCSD(T)-kappa-0',-1) 1675C 1676C----------------------------------------- 1677C Initialization of timing parameters. 1678C----------------------------------------- 1679C 1680 TIMTOT = ZERO 1681 TIMTOT = SECOND() 1682 TIMDEN = ZERO 1683 TIMRES = ZERO 1684 TIRDAO = ZERO 1685 TIMHE2 = ZERO 1686 TIMONE = ZERO 1687 TIMONE = SECOND() 1688 1689 1690 IF (LTSTEN) EN2PT=ZERO 1691C 1692C---------------------------------------------------- 1693C Both zeta- and t-vectors are totally symmetric. 1694C---------------------------------------------------- 1695C 1696 ISYMTR = 1 1697 ISYMOP = 1 1698C 1699C---------------------------------------- 1700C Get CMO coefficients 1701C---------------------------------------- 1702C 1703 KT1AM = 1 1704 KD1PTAI = KT1AM + NT1AMX 1705 KD1PTAB = KD1PTAI + NT1AM(1) 1706 KD1PTIJ = KD1PTAB + NMATAB(1) 1707 KEND0 = KD1PTIJ + NMATIJ(1) 1708 1709 call dzero(work(kd1ptai),nt1am(1)) 1710 call dzero(work(kd1ptij),nmatij(1)) 1711 call dzero(work(kd1ptab),nmatab(1)) 1712 call dzero(work(kt1am),nt1amx) 1713 1714 KCMO = KEND0 1715 KEND1 = KCMO + NLAMDS 1716 LWRK1 = LWORK - KEND1 1717 IF (LWRK1 .LT. 0) THEN 1718 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1719 CALL QUIT('Insufficient memory for allocation 1 CCPT_DEN2') 1720 ENDIF 1721C 1722 IF (FROIMP) THEN 1723 KCMOF = KEND1 1724 KEND1 = KCMOF + NLAMDS 1725 LWRK1 = LWORK - KEND1 1726 IF (LWRK1 .LT. 0) THEN 1727 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1728 CALL QUIT('Insufficient memory for allocation 2 CCPT_DEN2') 1729 ENDIF 1730C 1731C------------------------------------------- 1732C Get the FULL MO coefficient matrix. 1733C------------------------------------------- 1734C 1735 CALL CMO_ALL(WORK(KCMOF),WORK(KEND1),LWRK1) 1736C 1737 END IF 1738C 1739C---------------------------------------------------------- 1740C Read MO-coefficients from interface file and reorder. 1741C---------------------------------------------------------- 1742C 1743 LUSIFC = -1 1744 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 1745 * IDUMMY,.FALSE.) 1746 REWIND LUSIFC 1747 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1748 READ (LUSIFC) 1749 READ (LUSIFC) 1750 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 1751 CALL GPCLOSE (LUSIFC,'KEEP') 1752C 1753 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 1754 1755C 1756C------------------------------------- 1757C Two-electron part starts here..... 1758C------------------------------------- 1759C 1760 TIMONE = SECOND() - TIMONE 1761C 1762C----------------------------------- 1763C Start the loop over integrals. 1764C----------------------------------- 1765C 1766 KENDS2 = KEND1 1767 LWRKS2 = LWRK1 1768C 1769 IF (DIRECT) THEN 1770 IF (HERDIR) THEN 1771 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 1772 ELSE 1773 KCCFB1 = KEND1 1774 KINDXB = KCCFB1 + MXPRIM*MXCONT 1775 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 1776 LWRK1 = LWORK - KEND1 1777 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 1778 * KODPP1,KODPP2,KRDPP1,KRDPP2, 1779 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 1780 * WORK(KEND1),LWRK1,IPRERI) 1781 KEND1 = KFREE 1782 LWRK1 = LFREE 1783 ENDIF 1784 NTOSYM = 1 1785 ELSE 1786 NTOSYM = NSYM 1787 ENDIF 1788C 1789 KENDSV = KEND1 1790 LWRKSV = LWRK1 1791C 1792 ICDEL1 = 0 1793 1794 DO 100 ISYMD1 = 1,NTOSYM 1795C 1796 IF (DIRECT) THEN 1797 IF (HERDIR) THEN 1798 NTOT = MAXSHL 1799 ELSE 1800 NTOT = MXCALL 1801 ENDIF 1802 ELSE 1803 NTOT = NBAS(ISYMD1) 1804 ENDIF 1805C 1806 DO 110 ILLL = 1,NTOT 1807C 1808C--------------------------------------------- 1809C If direct calculate the integrals. 1810C--------------------------------------------- 1811C 1812 IF (DIRECT) THEN 1813C 1814 KEND1 = KENDSV 1815 LWRK1 = LWRKSV 1816C 1817 DTIME = SECOND() 1818 IF (HERDIR) THEN 1819 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 1820 & IPRINT) 1821 ELSE 1822 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 1823 * WORK(KODCL1),WORK(KODCL2), 1824 * WORK(KODBC1),WORK(KODBC2), 1825 * WORK(KRDBC1),WORK(KRDBC2), 1826 * WORK(KODPP1),WORK(KODPP2), 1827 * WORK(KRDPP1),WORK(KRDPP2), 1828 * WORK(KCCFB1),WORK(KINDXB), 1829 * WORK(KEND1), LWRK1,IPRERI) 1830 ENDIF 1831 DTIME = SECOND() - DTIME 1832 TIMHE2 = TIMHE2 + DTIME 1833C 1834 KRECNR = KEND1 1835 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 1836 LWRK1 = LWORK - KEND1 1837 IF (LWRK1 .LT. 0) THEN 1838 CALL QUIT('Insufficient core in CC_DEN_PT2') 1839 END IF 1840C 1841 ELSE 1842 NUMDIS = 1 1843 ENDIF 1844C 1845C----------------------------------------------------- 1846C Loop over number of distributions in disk. 1847C----------------------------------------------------- 1848C 1849 DO 120 IDEL2 = 1,NUMDIS 1850C 1851 IF (DIRECT) THEN 1852 IDEL = INDEXA(IDEL2) 1853 IF (NOAUXB) THEN 1854 IDUM = 1 1855 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 1856 END IF 1857 ISYMD = ISAO(IDEL) 1858 ELSE 1859 IDEL = IBAS(ISYMD1) + ILLL 1860 ISYMD = ISYMD1 1861 ENDIF 1862C 1863C--------------------------------------------------------- 1864C Sonia 1865C Work space allocation for the (T) densities 1866C with third index backtransformed to gamma 1867C All gammas together 1868C--------------------------------------------------------- 1869C 1870 ISYDEN = ISYMD 1871C 1872 KD2IJG_PT = KEND1 1873 KD2AIG_PT = KD2IJG_PT + ND2IJG(ISYDEN) 1874 KD2IAG_PT = KD2AIG_PT + ND2AIG(ISYDEN) 1875 KD2ABG_PT = KD2IAG_PT + ND2AIG(ISYDEN) 1876 KEND2 = KD2ABG_PT + ND2ABG(ISYDEN) 1877 LWRK2 = LWORK - KEND2 1878C 1879 IF (LWRK2 .LT. 0) THEN 1880 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 1881 CALL QUIT('Insufficient space for allocation '// 1882 & '2and1/2 in CCPT_DEN2') 1883 ENDIF 1884C 1885C------------------------------------------------------- 1886C Initialize 4 two electron density arrays. 1887C------------------------------------------------------- 1888C 1889 CALL DZERO(WORK(KD2IJG_PT),ND2IJG(ISYDEN)) 1890 CALL DZERO(WORK(KD2AIG_PT),ND2AIG(ISYDEN)) 1891 CALL DZERO(WORK(KD2IAG_PT),ND2AIG(ISYDEN)) 1892 CALL DZERO(WORK(KD2ABG_PT),ND2ABG(ISYDEN)) 1893C 1894C------------------------------------------------------------------- 1895C Calculate the two electron density d(pq,gamma;delta). 1896C------------------------------------------------------------------- 1897C 1898 AUTIME = SECOND() 1899C 1900 CALL CC_DEN2_PT(WORK(KD2IJG_PT),WORK(KD2AIG_PT), 1901 & WORK(KD2IAG_PT),WORK(KD2ABG_PT), 1902 & WORK(KCMO),1, 1903 & WORK(KEND2),LWRK2, 1904 & IDEL,ISYMD) 1905C 1906 AUTIME = SECOND() - AUTIME 1907 TIMDEN = TIMDEN + AUTIME 1908C 1909C------------------------------------------ 1910C Work space allocation three. 1911C------------------------------------------ 1912C 1913 ISYDIS = MULD2H(ISYMD,ISYMOP) 1914C 1915 KXINT = KEND2 1916 KEND3 = KXINT + NDISAO(ISYDIS) 1917 LWRK3 = LWORK - KEND3 1918C 1919 IF (LWRK3 .LT. 0) THEN 1920 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 1921 CALL QUIT('Insufficient space for allocation '// 1922 & '3 in CC_DEN_PT2') 1923 ENDIF 1924C 1925C-------------------------------------------- 1926C Read AO integral distribution. 1927C-------------------------------------------- 1928C 1929 AUTIME = SECOND() 1930 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 1931 * WORK(KRECNR),DIRECT) 1932 AUTIME = SECOND() - AUTIME 1933 TIRDAO = TIRDAO + AUTIME 1934 1935C 1936C------------------------------------------------------ 1937C Start loop over second AO-index (gamma). 1938C------------------------------------------------------ 1939C 1940 AUTIME = SECOND() 1941C 1942 DO 130 ISYMG = 1,NSYM 1943C 1944 ISYMPQ = MULD2H(ISYMG,ISYDEN) 1945C 1946 DO 140 G = 1,NBAS(ISYMG) 1947C 1948 KD2GIJ = KD2IJG_PT + ID2IJG(ISYMPQ,ISYMG) 1949 * + NMATIJ(ISYMPQ)*(G - 1) 1950 KD2GAI = KD2AIG_PT + ID2AIG(ISYMPQ,ISYMG) 1951 * + NT1AM(ISYMPQ)*(G - 1) 1952 KD2GIA = KD2IAG_PT + ID2AIG(ISYMPQ,ISYMG) 1953 * + NT1AM(ISYMPQ)*(G - 1) 1954 KD2GAB = KD2ABG_PT + ID2ABG(ISYMPQ,ISYMG) 1955 * + NMATAB(ISYMPQ)*(G - 1) 1956C 1957C----------------------------------------------- 1958C Work space allocation four. 1959C----------------------------------------------- 1960C 1961 KINTAO = KEND3 1962 KD2AOB = KINTAO + N2BST(ISYMPQ) 1963 KEND4 = KD2AOB + N2BST(ISYMPQ) 1964 LWRK4 = LWORK - KEND4 1965C 1966 IF (LWRK4 .LT. 0) THEN 1967 WRITE(LUPRI,*) 'Available:', LWORK 1968 WRITE(LUPRI,*) 'Needed:', KEND4 1969 CALL QUIT('Insufficient space in CC_DEN_PT2') 1970 ENDIF 1971C 1972 CALL DZERO(WORK(KD2AOB),N2BST(ISYMPQ)) 1973c 1974cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 1975C 1976C------------------------------------------------------------- 1977C Calculate frozen core contributions to d. 1978C------------------------------------------------------------- 1979C 1980 IF (FROIMP) THEN 1981C 1982 KFD2IJ = KEND4 1983 KFD2JI = KFD2IJ + NCOFRO(ISYMPQ) 1984 KFD2AI = KFD2JI + NCOFRO(ISYMPQ) 1985 KFD2IA = KFD2AI + NT1FRO(ISYMPQ) 1986 KFD2II = KFD2IA + NT1FRO(ISYMPQ) 1987 KEND4 = KFD2II + NFROFR(ISYMPQ) 1988 LWRK4 = LWORK - KEND4 1989C 1990 IF (LWRK4 .LT. 0) THEN 1991 WRITE (LUPRI,*) 'Available:', LWORK 1992 WRITE (LUPRI,*) 'Needed:', KEND4 1993 CALL QUIT( 1994 * 'Insufficient work space in CC_DEN') 1995 ENDIF 1996C 1997 CALL DZERO(WORK(KFD2IJ),NCOFRO(ISYMPQ)) 1998 CALL DZERO(WORK(KFD2JI),NCOFRO(ISYMPQ)) 1999 CALL DZERO(WORK(KFD2AI),NT1FRO(ISYMPQ)) 2000 CALL DZERO(WORK(KFD2IA),NT1FRO(ISYMPQ)) 2001 CALL DZERO(WORK(KFD2II),NFROFR(ISYMPQ)) 2002C 2003 CALL CCPT_FD2BL(WORK(KFD2II),WORK(KFD2IJ), 2004 * WORK(KFD2JI),WORK(KFD2AI), 2005 * WORK(KFD2IA),D1PTIA, 2006 & WORK(KCMO),WORK(KCMOF), 2007 * WORK(KEND4),LWRK4, 2008 & IDEL,ISYMD,G,ISYMG) 2009 2010 CALL CC_FD2AO(WORK(KD2AOB),WORK(KFD2II), 2011 * WORK(KFD2IJ),WORK(KFD2JI), 2012 * WORK(KFD2AI),WORK(KFD2IA), 2013 * WORK(KCMOF),WORK(KCMO), 2014 * WORK(KCMO),WORK(KEND4),LWRK4, 2015 * ISYMPQ) 2016 2017 CALL CCPT_D2GAF(WORK(KD2GIA), 2018 * D1PTIA,WORK(KCMOF), 2019 & IDEL,ISYMD,G,ISYMG) 2020 END IF 2021 2022cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 2023C 2024C------------------------------------------------------- 2025C Square up AO-integral distribution. 2026C------------------------------------------------------- 2027C 2028 KOFFIN = KXINT + IDSAOG(ISYMG,ISYDIS) 2029 * + NNBST(ISYMPQ)*(G - 1) 2030C 2031 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ, 2032 * WORK(KINTAO)) 2033 2034C 2035C--------------------------------------------------------- 2036C 2037C--------------------------------------------------------- 2038C 2039 if (LTSTEN) then 2040 2041 CALL CC_DENAO(WORK(KD2AOB),ISYMPQ, 2042 * WORK(KD2GAI),WORK(KD2GAB), 2043 * WORK(KD2GIJ),WORK(KD2GIA),ISYMPQ, 2044 * WORK(KCMO),1,WORK(KCMO),1, 2045 * WORK(KEND4),LWRK4) 2046C 2047C---------------------------------------------------------------------- 2048C Calculate the 2 e- density contribution to E-ccsd. 2049C---------------------------------------------------------------------- 2050C 2051 EN2PT = EN2PT + HALF*DDOT(N2BST(ISYMPQ), 2052 * WORK(KD2AOB),1,WORK(KINTAO),1) 2053C 2054 end if 2055C 2056C----------------------------------------------- 2057C Work space allocation five. 2058C----------------------------------------------- 2059C 2060 KIJINT = KEND4 2061 KAIINT = KIJINT + NMATIJ(ISYMPQ) 2062 KIAINT = KAIINT + NT1AM(ISYMPQ) 2063 KABINT = KIAINT + NT1AM(ISYMPQ) 2064 KEND5 = KABINT + NMATAB(ISYMPQ) 2065 LWRK5 = LWORK - KEND5 2066!Sonia: skod allocations 2067 KABTIN = KEND5 2068 KIJTIN = KABTIN + NMATAB(ISYMPQ) 2069 KD2TAB = KIJTIN + NMATIJ(ISYMPQ) 2070 KD2TIJ = KD2TAB + NMATAB(ISYMPQ) 2071 KEND5 = KD2TIJ + NMATIJ(ISYMPQ) 2072 LWRK5 = LWORK - KEND5 2073c 2074 IF (FROIMP) THEN 2075 KIIFRO = KEND5 2076 KIJFRO = KIIFRO + NFROFR(ISYMPQ) 2077 KJIFRO = KIJFRO + NCOFRO(ISYMPQ) 2078 KAIFRO = KJIFRO + NCOFRO(ISYMPQ) 2079 KIAFRO = KAIFRO + NT1FRO(ISYMPQ) 2080 KEND5 = KIAFRO + NT1FRO(ISYMPQ) 2081 LWRK5 = LWORK - KEND5 2082 END IF 2083c 2084 IF (LWRK5 .LT. 0) THEN 2085 WRITE(LUPRI,*) 'Available:', LWORK 2086 WRITE(LUPRI,*) 'Needed:', KEND5 2087 CALL QUIT('Insufficient work space '// 2088 & 'in CC_DEN_PT2') 2089 ENDIF 2090C 2091C---------------------------------------------------------------- 2092C Transform 2 integral indices to MO-basis. 2093C---------------------------------------------------------------- 2094C 2095 ISYM = ISYMPQ 2096 CALL CCDINTMO(WORK(KIJINT),WORK(KIAINT), 2097 * WORK(KABINT),WORK(KAIINT), 2098 * WORK(KINTAO),WORK(KCMO), 2099 * WORK(KCMO),WORK(KEND5), 2100 * LWRK5,ISYM) 2101 2102 IF (FROIMP) THEN 2103C 2104C Prepare integrals g_pq (gam,del) where one index is frozen 2105C 2106 ISYM = ISYMPQ 2107 CALL CCIFROMO(WORK(KIJFRO),WORK(KJIFRO), 2108 * WORK(KAIFRO),WORK(KIAFRO), 2109 * WORK(KIIFRO),WORK(KINTAO), 2110 * WORK(KCMO),WORK(KCMO), 2111 * WORK(KCMOF),WORK(KEND5), 2112 * LWRK5,ISYM) 2113C 2114 ENDIF 2115C 2116C----------------------------------------------------------------- 2117C A lot of unnecessary allocs!!!! 2118C Set up transposed integrals and densities. 2119C----------------------------------------------------------------- 2120C 2121 CALL DCOPY(NMATAB(ISYMPQ),WORK(KABINT),1, 2122 * WORK(KABTIN),1) 2123 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KIJINT),1, 2124 * WORK(KIJTIN),1) 2125 CALL DCOPY(NMATAB(ISYMPQ),WORK(KD2GAB),1, 2126 * WORK(KD2TAB),1) 2127 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KD2GIJ),1, 2128 * WORK(KD2TIJ),1) 2129C 2130 CALL CC_EITR(WORK(KABTIN),WORK(KIJTIN), 2131 * WORK(KEND5),LWRK5,ISYMPQ) 2132 CALL CC_EITR(WORK(KD2TAB),WORK(KD2TIJ), 2133 * WORK(KEND5),LWRK5,ISYMPQ) 2134C 2135C------------------------------------------------------------------- 2136C Calculate 2 e- contribution to Zeta-Kappa-0. 2137C------------------------------------------------------------------- 2138C 2139 ISYM = ISYMPQ 2140 2141 IF (IOPT.EQ.2) THEN 2142 2143 KETAIJ = 1 2144 KETAAB = 1 + NMATIJ(1) 2145 2146 CALL CCPT_ETARS_2E( 2147 & ZKDIA(KETAIJ),ZKDIA(KETAAB), 2148 & WORK(KIJINT),WORK(KAIINT), 2149 & WORK(KIAINT),WORK(KABINT), 2150 & WORK(KD2GIJ),WORK(KD2GAI), 2151 & WORK(KD2GIA),WORK(KD2GAB), 2152 & WORK(KEND5),LWRK5,ISYM) 2153 END IF 2154 2155 CALL CCPT_ETAAI_2E(ETAAI, 2156 & WORK(KIJINT),WORK(KAIINT), 2157 & WORK(KIAINT),WORK(KABINT), 2158 & WORK(KD2GIJ),WORK(KD2GAI), 2159 & WORK(KD2GIA),WORK(KD2GAB), 2160 & WORK(KEND5),LWRK5, 2161 & ISYM) 2162 2163 IF (FROIMP) THEN 2164 2165 ISYM = ISYMPQ 2166! 2167!Contributi a eta_ai dai loop sui frozen 2168! 2169 CALL CCFRETAI(ETAAI,WORK(KIJFRO), 2170 * WORK(KJIFRO),WORK(KAIFRO), 2171 * WORK(KIAFRO),WORK(KFD2IJ), 2172 * WORK(KFD2JI),WORK(KFD2AI), 2173 * WORK(KFD2IA),WORK(KT1AM), 2174 * WORK(KEND5),LWRK5,ISYM) 2175C 2176C----------------------------------------------------------------------- 2177C Calculate two-electron contribution to right- 2178C hand-side of correlated-frozen multipliers. 2179C----------------------------------------------------------------------- 2180C 2181 ICON = 2 2182 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 2183 * + 2*NT1FRO(1) + 1 2184C 2185 CALL CC_ETIJF(ZKDIA(KOFRES),WORK(KD2GIJ), 2186 * WORK(KD2GAB),WORK(KD2GAI), 2187 * WORK(KD2GIA),WORK(KD2TIJ), 2188 * WORK(KFD2II),WORK(KFD2IJ), 2189 * WORK(KFD2JI),WORK(KFD2AI), 2190 * WORK(KFD2IA),WORK(KIJINT), 2191 * WORK(KAIINT),WORK(KIAINT), 2192 * WORK(KIJTIN),WORK(KABTIN), 2193 * WORK(KIJFRO),WORK(KJIFRO), 2194 * WORK(KAIFRO),WORK(KIAFRO), 2195 * WORK(KIIFRO),WORK(KT1AM), 2196 * WORK(KEND5),LWRK5,ISYM,ICON) 2197C 2198C----------------------------------------------------------------------- 2199C Calculate two-electron contribution to right- 2200C hand-side of virtual-frozen multipliers. et_aI 2201C----------------------------------------------------------------------- 2202C 2203 ICON = 2 2204 KOFRE = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 2205C 2206 CALL CC_ETAIF(ZKDIA(KOFRE),WORK(KD2GAB), 2207 * WORK(KD2GAI),WORK(KD2GIA), 2208 * WORK(KD2TIJ),WORK(KD2TAB), 2209 * WORK(KFD2II),WORK(KFD2IJ), 2210 * WORK(KFD2JI),WORK(KFD2AI), 2211 * WORK(KFD2IA),WORK(KIJINT), 2212 * WORK(KABINT),WORK(KAIINT), 2213 * WORK(KIAINT),WORK(KABTIN), 2214 * WORK(KIJFRO),WORK(KJIFRO), 2215 * WORK(KAIFRO),WORK(KIAFRO), 2216 * WORK(KIIFRO),WORK(KT1AM), 2217 * WORK(KEND5),LWRK5,ISYM,ICON) 2218 2219! CALL CCDIAZK0(ZKDIA,WORK(KIJINT), 2220! * WORK(KAIINT), 2221! * WORK(KIAINT),WORK(KABINT), 2222! * WORK(KD2GIJ),WORK(KD2GAI), 2223! * WORK(KD2GIA),WORK(KD2GAB), 2224! * WORK(KIJTIN),WORK(KABTIN), 2225! * WORK(KD2TIJ),WORK(KD2TAB), 2226! * WORK(KT1AM),WORK(KEND5), 2227! * LWRK5,ISYM) 2228 2229 CALL CCFRETAB(ZKDIA,WORK(KIJFRO), 2230 * WORK(KJIFRO),WORK(KAIFRO), 2231 * WORK(KIAFRO),WORK(KFD2IJ), 2232 * WORK(KFD2JI),WORK(KFD2AI), 2233 * WORK(KFD2IA),WORK(KT1AM), 2234 * WORK(KEND5),LWRK5,ISYM) 2235 2236 CALL CCFRETIJ(ZKDIA,WORK(KIJFRO), 2237 * WORK(KJIFRO),WORK(KAIFRO), 2238 * WORK(KIAFRO),WORK(KFD2IJ), 2239 * WORK(KFD2JI),WORK(KFD2AI), 2240 * WORK(KFD2IA),WORK(KT1AM), 2241 * WORK(KEND5),LWRK5,ISYM) 2242C 2243 END IF 2244 140 CONTINUE 2245 130 CONTINUE 2246C 2247 AUTIME = SECOND() - AUTIME 2248 TIMRES = TIMRES + AUTIME 2249 2250 120 CONTINUE 2251 110 CONTINUE 2252 100 CONTINUE 2253C 2254C----------------------- 2255C Regain work space. 2256C----------------------- 2257C 2258 KEND1 = KENDS2 2259 LWRK1 = LWRKS2 2260C 2261C------------------------ 2262C 2263C------------------------ 2264C 2265 if (ltsten) then 2266 write(lupri,*)'CC_DEN_PT2--> EN2PT: ', EN2PT 2267 end if 2268C 2269C----------------------- 2270C Write out timings. 2271C----------------------- 2272C 2273 99 TIMTOT = SECOND() - TIMTOT 2274C 2275 IF (IPRINT .GT. 3) THEN 2276 WRITE (LUPRI,*) ' ' 2277 WRITE (LUPRI,*) 'Requested density dependent '// 2278 & 'quantities calculated' 2279 WRITE (LUPRI,*) 'Total time used in CC_DEN_PT2:', TIMTOT 2280 ENDIF 2281 IF (IPRINT .GT. 9) THEN 2282 WRITE (LUPRI,*) 'Time used for setting up 2 e- density:',TIMDEN 2283 WRITE (LUPRI,*) 'Time used for contraction with integrals:', 2284 & TIMRES 2285 WRITE (LUPRI,*) 'Time used for reading 2 e- AO-integrals:', 2286 & TIRDAO 2287 WRITE (LUPRI,*) 'Time used for calculating 2 e- AO-integrals:', 2288 * TIMHE2 2289 WRITE (LUPRI,*) 'Time used for 1 e- density & intermediates:', 2290 * TIMONE 2291 ENDIF 2292C 2293 CALL QEXIT('CCPT_DEN2FC') 2294 RETURN 2295 END 2296