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! 17C 18C /* Deck cc_den_rccd */ 19 SUBROUTINE CC_DEN_RCCD(POTNUC,ETAAI,ZKDIA,WORK,LWORK, 20 & IOPT,IMODEL,LTSTE) 21C 22C Written by S. Coriani, based on CC_DEN_PTFC 23C Debugged version using particle-symmetrized densities 24C 25C Version: 3.0 26C 27C Current models: RCCD & DRCCD & SOSEX 28C 29C LTSTE = .true. test densities via Energy calculation 30C 31#include "implicit.h" 32#include "priunit.h" 33#include "dummy.h" 34#include "maxorb.h" 35#include "maxash.h" 36#include "mxcent.h" 37#include "aovec.h" 38#include "iratdef.h" 39 PARAMETER (ZERO = 0.0D0,HALF=0.5D0,ONE = 1.0D0,TWO = 2.0D0) 40 PARAMETER (TRE = 3.0D0, FOUR = 4.0D0) 41 DIMENSION INDEXA(MXCORB_CC) 42 DIMENSION ETAAI(*), ZKDIA(*), WORK(LWORK) 43 LOGICAL LTSTE, LETAFI, LETIFJ 44#include "ccorb.h" 45#include "ccisao.h" 46#include "r12int.h" 47#include "inftap.h" 48#include "blocks.h" 49#include "ccfield.h" 50#include "ccsdinp.h" 51#include "ccinftap.h" 52#include "ccsdsym.h" 53#include "ccsdio.h" 54#include "distcl.h" 55#include "cbieri.h" 56#include "eritap.h" 57#include "ccfro.h" 58CAMT 59C#include "dftcom.h" 60C#include "oepopt.h" 61C#include "ccandy.h" 62 63CTEST 64C#include "ccfop.h" 65 66 CHARACTER MODEL*10 67 CHARACTER NAME1*8 68 CHARACTER NAME2*8 69 70 LOGICAL LOCDBG 71 PARAMETER (LOCDBG=.FALSE.) 72C 73 CALL QENTER('CC_DEN_RCCD') 74 75C 76 IF (FROIMP) THEN 77C 78 NAME1 = 'CCFRO1IN' 79 NAME2 = 'CCFRO2IN' 80C 81 LUN1 = -1 82 LUN2 = -1 83C 84 CALL WOPEN2(LUN1,NAME1,64,0) 85 CALL WOPEN2(LUN2,NAME2,64,0) 86C 87 ENDIF 88C 89 IF (IOPT .LE. 2) THEN 90 !IF (LPRNCC) THEN 91 CALL HEADER('CC_DEN_RCCD: constructing RHS for RCCD-kapbar-0', 92 & -1) 93 call flshfo(lupri) 94 !ENDIF 95 ENDIF 96C 97C----------------------------------------- 98C Initialization of timing parameters. 99C----------------------------------------- 100C 101 TIMTOT = ZERO 102 TIMTOT = SECOND() 103 TIMDEN = ZERO 104 TIMRES = ZERO 105 TIRDAO = ZERO 106 TIMHE2 = ZERO 107 TIMONE = ZERO 108 TIMONE = SECOND() 109C 110C---------------------------------------------------- 111C Both zeta- and t-vectors are totally symmetric. 112C---------------------------------------------------- 113C 114 ISYMTR = 1 115 ISYMOP = 1 116C 117 LUNGO = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 118 * + 2*NCOFRO(1) + 2*NT1FRO(1) 119C 120C----------------------------------- 121C Initial work space allocation. 122C----------------------------------- 123C 124 IF (LTSTE) THEN 125 KD1AOB = 1 126 KSTART = KD1AOB + N2BST(1) 127 ELSE 128 KSTART = 1 129 END IF 130 131 KZ2AM = KSTART 132 KT2AM = KZ2AM + NT2SQ(1) 133 KT2AMT = KT2AM + NT2AMX 134 KLAMDP = KT2AMT + NT2AMX 135 KLAMDH = KLAMDP + NLAMDT 136 KZ2TIL = KLAMDH + NLAMDT !2C-E of multipliers 137 KZ2PCK = KZ2TIL + NT2SQ(1) 138 KT2SQ = KZ2PCK + NT2AMX 139 KEND1 = KT2SQ + NT2SQ(1) 140 LWRK1 = LWORK - KEND1 141C 142 IF (LWRK1 .LT. 0) THEN 143 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 144 CALL QUIT('Insufficient memory for initial allocation '// 145 & 'in CC_DEN_RCCD') 146 ENDIF 147C 148C---------------------------------------- 149C Read zero-th order zeta amplitudes. 150C---------------------------------------- 151C 152 IOPTRW = 2 153 CALL CC_RDRSP('L0',0,1,IOPTRW,MODEL,WORK(KEND1),WORK(KZ2AM)) 154 call flshfo(lupri) 155C 156C-------------------------------------------------------- 157C Calculate tbar_tilde = 2C-E of Tbar for RCCD and 158C for dRCCD just 2*tbar in squared form 159C and save a packed copy of Tbar in KZ2PCK 160C-------------------------------------------------------- 161C 162 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KZ2PCK),1) 163 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 164 if (RCCD) then 165 ISYOPE = 1 166 IOPTTCME = 1 167 CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME) 168 else 169 CALL DSCAL(NT2AMX,TWO,WORK(KT2AM),1) 170 end if 171 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2TIL),1) 172C 173C------------------------------------------------------------- 174C Square up zeta2 amplitudes. 175C------------------------------------------------------------- 176C 177 CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1) 178 CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1) 179C 180C------------------------------------------- 181C Read zero'th order cluster amplitudes. 182C------------------------------------------- 183C 184 IOPTRW = 2 185 CALL CC_RDRSP('R0',0,1,IOPTRW,MODEL,WORK(KEND1),WORK(KT2AM)) 186 CALL CC_T2SQ(WORK(KT2AM),WORK(KT2SQ),1) 187C 188C------------------------------------------------- 189C Set up 2C-E of cluster amplitudes (T2 tilde). 190C for RCCD and SOSEX, otherwise just 2*ampl 191C------------------------------------------------- 192C 193 ISYOPE = 1 194 CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1) 195 if (DRCCD) then 196 if (SOSEX) then 197 IOPTTCME = 1 198 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 199 else 200 CALL DSCAL(NT2AMX,TWO,WORK(KT2AMT),1) 201 end if 202 else !if (RCCD) then 203 IOPTTCME = 1 204 CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME) 205 end if 206C 207C---------------------------------------------------------------- 208C Calculate the lambda matrices. 209C Redundant, it's just CMO, but I let it go to avoid problems 210C---------------------------------------------------------------- 211C 212 KT1AM = KEND1 213 KEND1 = KT1AM + NT1AMX 214 LWRK1 = LWORK-KEND1 215 CALL DZERO(WORK(KT1AM),NT1AMX) 216 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 217 * LWRK1) 218C 219C---------------------------------------------- 220C Work space allocation one. CCSD-like part 221C Note that D(ai) = D(ia) = 0, and both 222C D(ia) and h(ia) are stored transposed! 223C---------------------------------------------- 224C 225 KONEAI = KEND1 226 KONEAB = KONEAI + NT1AMX 227 KONEIJ = KONEAB + NMATAB(1) 228 KONEIA = KONEIJ + NMATIJ(1) 229 KXMAT = KONEIA + NT1AMX 230 KYMAT = KXMAT + NMATIJ(1) 231 KONINT = KYMAT + NMATAB(1) 232C 233 KINTAI = KONINT + N2BST(ISYMOP) 234 KINTAB = KINTAI + NT1AMX 235 KINTIJ = KINTAB + NMATAB(1) 236 KINTIA = KINTIJ + NMATIJ(1) 237 KINABT = KINTIA + NT1AMX 238 KINIJT = KINABT + NMATAB(1) 239 KD1ABT = KINIJT + NMATIJ(1) 240 KD1IJT = KD1ABT + NMATAB(1) 241 KEND1 = KD1IJT + NMATIJ(1) 242 LWRK1 = LWORK - KEND1 243 244 IF (FROIMP) THEN 245 KFROII = KEND1 246 KFROIJ = KFROII + NFROFR(1) 247 KFROJI = KFROIJ + NCOFRO(1) 248 KFROAI = KFROJI + NCOFRO(1) 249 KFROIA = KFROAI + NT1FRO(1) 250 KFD1II = KFROIA + NT1FRO(1) 251 KEND1 = KFD1II + NFROFR(1) 252 LWRK1 = LWORK - KEND1 253 ENDIF 254 255 KCMO = KEND1 256 KEND1 = KCMO + NLAMDS 257 LWRK1 = LWORK - KEND1 258 259 IF (FROIMP) THEN 260 KCMOF = KEND1 261 KEND1 = KCMOF + NLAMDS 262 LWRK1 = LWORK - KEND1 263 ENDIF 264C 265 IF (LWRK1 .LT. 0) THEN 266 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 267 CALL QUIT('Insuff. memory for allocation 1 CC_DEN_RCCD') 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 Get the non-frozen MO coefficient matrix reorder. 282C------------------------------------------------- 283C 284 CALL CC_GET_CMO(WORK(KCMO)) 285 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 286C 287C------------------------------------------------------ 288C Initialize remaining one electron density arrays. 289C------------------------------------------------------ 290C 291 CALL DZERO(WORK(KONEAB),NMATAB(1)) 292 CALL DZERO(WORK(KONEIJ),NMATIJ(1)) 293 CALL DZERO(WORK(KONEIA),NT1AMX) 294 CALL DZERO(WORK(KONEAI),NT1AMX) 295C 296C-------------------------------------------------------- 297C Calculate X-intermediate of zeta- and t-amplitudes. 298C-------------------------------------------------------- 299C 300 CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 301 * WORK(KEND1),LWRK1) 302 CALL DSCAL(NMATIJ(1),TWO,WORK(KXMAT),1) 303C 304C-------------------------------------------------------- 305C Calculate Y-intermediate of zeta- and t-amplitudes. 306C-------------------------------------------------------- 307C 308 CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP, 309 * WORK(KEND1),LWRK1) 310 CALL DSCAL(NMATAB(1),TWO,WORK(KYMAT),1) 311C 312C------------------------------------------------------------------------ 313C Construct non-zero blocks of one electron density. 314C Note that X and Y are actually 2*X and 2*Y 315C------------------------------------------------------------------------ 316C 317 CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1) 318 CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1) 319 CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT)) 320 321C rescale X and Y back to "true" X and Y value 322 CALL DSCAL(NMATIJ(1),HALF,WORK(KXMAT),1) 323 CALL DSCAL(NMATAB(1),HALF,WORK(KYMAT),1) 324 325 IF (LOCDBG) THEN 326 DIJNO = DDOT(NMATIJ(1),WORK(KONEIJ),1,WORK(KONEIJ),1) 327 DAINO = DDOT(NT1AMX,WORK(KONEAI),1,WORK(KONEAI),1) 328 DIANO = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1) 329 DABNO = DDOT(NMATAB(1),WORK(KONEAB),1,WORK(KONEAB),1) 330 WRITE(LUPRI,*) 'CC_DEN_RCCD: IOPT = ', IOPT 331 WRITE(LUPRI,*) 332 & 'Norm of ',MODEL(1:5),' one electron densities in MO-basis:' 333 WRITE(LUPRI,*) DIJNO, DAINO, DIANO, DABNO 334 call flshfo(lupri) 335 336 write(lupri,*)'The IJ density of ', MODEL(1:5) 337 CALL OUTPUT(WORK(KONEIJ),1,NRHF(1),1,NRHF(1), 338 * NRHF(1),NRHF(1),1,LUPRI) 339 write(lupri,*)'The AI density of ', MODEL(1:5) 340 CALL OUTPUT(WORK(KONEAI),1,NVIR(1),1,NRHF(1), 341 * NVIR(1),NRHF(1),1,LUPRI) 342 write(lupri,*)'The IA density of ', MODEL(1:5) 343 CALL OUTPUT(WORK(KONEIA),1,NVIR(1),1,NRHF(1), 344 * NVIR(1),NRHF(1),1,LUPRI) 345 write(lupri,*)'The AB density of ', MODEL(1:5) 346 CALL OUTPUT(WORK(KONEAB),1,NVIR(1),1,NVIR(1), 347 * NVIR(1),NVIR(1),1,LUPRI) 348 349 ENDIF !locdbg 350C 351C--------------------------------- 352C Read one-electron integrals. 353C--------------------------------- 354C 355 CALL CCRHS_ONEAO(WORK(KONINT),WORK(KEND1),LWRK1) 356 357 IF (LTSTE) THEN 358 !IF (LPRNCC) write(lupri,*)'LTSTE=', LTSTE 359 !call flshfo(lupri) 360 361 CALL DZERO(WORK(KD1AOB),N2BST(1)) 362 ISDEN = 1 363 CALL CC_DENAO(WORK(KD1AOB),ISDEN,WORK(KONEAI),WORK(KONEAB), 364 * WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1, 365 * WORK(KLAMDH),1,WORK(KEND1),LWRK1) 366C 367 IF (FROIMP) THEN 368 MODEL = 'DUMMY' 369 CALL CC_FCD1AO(WORK(KD1AOB),WORK(KEND1),LWRK1,MODEL) 370 END IF 371 372 ECCSD1 = DDOT(N2BST(ISYMOP),WORK(KD1AOB),1,WORK(KONINT),1) 373 ECCSD2 = ZERO 374 375 END IF !LTSTE 376C 377C--------------------------------------------------------- 378C Ove 24-20-1996 379C Read one-electron integrals into Fock-matrix for 380C finite field. 381C--------------------------------------------------------- 382C 383 DO 13 IF = 1, NFIELD 384 FF = EFIELD(IF) 385 CALL CC_ONEP(WORK(KONINT),WORK(KEND1),LWRK1,FF,1,LFIELD(IF)) 386 13 CONTINUE 387C 388C-------------------------------------------------- 389C Transform one electron integrals to MO-basis. 390C-------------------------------------------------- 391C 392 ISYM = 1 393 CALL CCDINTMO(WORK(KINTIJ),WORK(KINTIA),WORK(KINTAB), 394 * WORK(KINTAI),WORK(KONINT),WORK(KLAMDP), 395 * WORK(KLAMDH),WORK(KEND1),LWRK1,ISYM) 396C 397C 398 IF (FROIMP) THEN 399C 400 ISYM = 1 401 !obtain integrals with frozen indices 402 ! h_Ij h_jI h_aJ h_Ia h_IJ 403 ! 404 CALL CCIFROMO(WORK(KFROIJ),WORK(KFROJI),WORK(KFROAI), 405 * WORK(KFROIA),WORK(KFROII),WORK(KONINT), 406 * WORK(KLAMDP),WORK(KLAMDH),WORK(KCMOF), 407 * WORK(KEND1),LWRK1,ISYM) 408C 409 !calculate D_II = 2 delta_IJ 410 CALL CCFD1II(WORK(KFD1II)) 411C 412 ENDIF !froimp 413C 414C-------------------------------------------------- 415C Set up transposed integrals and densities. 416C-------------------------------------------------- 417C 418 CALL DCOPY(NMATAB(1),WORK(KINTAB),1,WORK(KINABT),1) 419 CALL DCOPY(NMATIJ(1),WORK(KINTIJ),1,WORK(KINIJT),1) 420 CALL DCOPY(NMATAB(1),WORK(KONEAB),1,WORK(KD1ABT),1) 421 CALL DCOPY(NMATIJ(1),WORK(KONEIJ),1,WORK(KD1IJT),1) 422C 423 CALL CC_EITR(WORK(KINABT),WORK(KINIJT),WORK(KEND1),LWRK1,1) 424 CALL CC_EITR(WORK(KD1ABT),WORK(KD1IJT),WORK(KEND1),LWRK1,1) 425C 426C------------------------------------------------------------ 427C Calculate one electron contribution to Zeta-kappa-0. 428C------------------------------------------------------------ 429C 430 ISYM = 1 431 432 IF (IOPT.EQ.2) THEN 433 434 !I let it go thru this unaltered as T1AM is set to zero 435 !compute eta_ij 436 KOFFIJ = 1 437 !CALL DZERO(NMATIJ(1),ZKDIA(KOFFIJ)) 438 CALL RCCD_ETIJ(ZKDIA(KOFFIJ),WORK(KINTIJ),WORK(KINTAI), 439 & WORK(KINTIA),WORK(KINTAB),WORK(KONEIJ), 440 & WORK(KONEAI),WORK(KONEIA),WORK(KONEAB), 441 & WORK(KEND1),LWRK1,ISYM) 442c write(lupri,*)'Norm of eta_ij=', 443c & SQRT(ABS(DDOT(NMATIJ(1),ZKDIA(KOFFIJ),1,ZKDIA(KOFFIJ),1))) 444 445 !I let it go thru this unaltered as T1AM is set to zero 446 !compute eta_ab 447 KOFFAB = NMATIJ(1) + 1 448 !CALL DZERO(NMATAB(1),ZKDIA(KOFFAB)) 449 CALL RCCD_ETAB(ZKDIA(KOFFAB),WORK(KINTIJ),WORK(KINTAI), 450 & WORK(KINTIA),WORK(KINTAB),WORK(KONEIJ), 451 & WORK(KONEAI),WORK(KONEIA),WORK(KONEAB), 452 & WORK(KEND1),LWRK1,ISYM) 453c write(lupri,*)'Norm of eta_ab=', 454c & SQRT(ABS(DDOT(NMATAB(1),ZKDIA(KOFFAB),1,ZKDIA(KOFFAB),1))) 455 456 END IF !iopt2 457 458C------------------------------------------------------------ 459 460 CALL DZERO(ETAAI,NALLAI(1)) 461 CALL DZERO(WORK(KT1AM),NT1AMX) 462 !I let it go thru this unaltered as T1AM is set to zero 463 !eta_ai 464 CALL CCDENZK0(ETAAI,WORK(KINTIJ),WORK(KINTAI),WORK(KINTIA), 465 * WORK(KINTAB),WORK(KONEIJ),WORK(KONEAI), 466 * WORK(KONEIA),WORK(KONEAB),WORK(KINIJT), 467 * WORK(KINABT),WORK(KD1IJT),WORK(KD1ABT), 468 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM) 469C 470! write(lupri,*)'The eta_IJ one-electron of ', MODEL(1:5) 471! CALL OUTPUT(ZKDIA(KOFFIJ),1,NRHF(1),1,NRHF(1), 472! * NRHF(1),NRHF(1),1,LUPRI) 473! write(lupri,*)'The eta_AI one-electron of ', MODEL(1:5) 474! CALL OUTPUT(ETAAI,1,NVIR(1),1,NRHF(1), 475! * NVIR(1),NRHF(1),1,LUPRI) 476! write(lupri,*)'The eta_IA one-electron of RPA' 477! CALL OUTPUT(WORK(KINTIA),1,NVIR(1),1,NRHF(1), 478! * NVIR(1),NRHF(1),1,LUPRI) 479! write(lupri,*)'The eta_AB one-electron of ', MODEL(1:5) 480! CALL OUTPUT(ZKDIA(KOFFAB),1,NVIR(1),1,NVIR(1), 481! * NVIR(1),NVIR(1),1,LUPRI) 482 483 484 IF (FROIMP) THEN 485C 486C-------------------------------------------------------- 487C Calculate one-electron contribution to right- 488C hand-side of correlated-frozen multipliers. 489C-------------------------------------------------------- 490C 491 ISYM = 1 492 ICON = 1 493 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 494 ! 495 !eta_iJ (one el) 496 ! 497 CALL CC_ETIJF(ZKDIA(KOFRES),WORK(KONEIJ),WORK(KONEAB), 498 * WORK(KONEAI),WORK(KONEIA),WORK(KD1IJT), 499 * WORK(KFD1II),SK1,SK2,SK3,SK4,WORK(KINTIJ), 500 * WORK(KINTAI),WORK(KINTIA),WORK(KINIJT), 501 * WORK(KINABT),WORK(KFROIJ),WORK(KFROJI), 502 * WORK(KFROAI),WORK(KFROIA),WORK(KFROII), 503 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM,ICON) 504C 505C-------------------------------------------------------- 506C Calculate one-electron contribution to right- 507C hand-side of virtual-frozen multipliers. 508C-------------------------------------------------------- 509C 510 ISYM = 1 511 ICON = 1 512 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 513 ! 514 !eta_aI 515 ! 516 CALL CC_ETAIF(ZKDIA(KOFRES),WORK(KONEAB),WORK(KONEAI), 517 * WORK(KONEIA),WORK(KD1IJT),WORK(KD1ABT), 518 * WORK(KFD1II),SK1,SK2,SK3,SK4,WORK(KINTIJ), 519 * WORK(KINTAB),WORK(KINTAI),WORK(KINTIA), 520 * WORK(KINABT),WORK(KFROIJ),WORK(KFROJI), 521 * WORK(KFROAI),WORK(KFROIA),WORK(KFROII), 522 * WORK(KT1AM),WORK(KEND1),LWRK1,ISYM,ICON) 523 524 ENDIF !froimp 525C 526 TIMONE = SECOND() - TIMONE 527C 528C-------------------------------------------- 529C Start the loop over 2-e integrals. 530C Salva tutto quanto definito fino ad ora 531C-------------------------------------------- 532C 533 ! IF (LPRNCC) 534 !& 535 WRITE(LUPRI,*)'DONE WITH 1-E, START 2-E, 1e Energy=', ECCSD1 536 call flshfo(lupri) 537 538 KENDS2 = KEND1 539 LWRKS2 = LWRK1 540C 541 IF (DIRECT) THEN 542 IF (HERDIR) THEN 543 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 544 ELSE 545 KCCFB1 = KEND1 546 KINDXB = KCCFB1 + MXPRIM*MXCONT 547 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 548 LWRK1 = LWORK - KEND1 549 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 550 * KODPP1,KODPP2,KRDPP1,KRDPP2, 551 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 552 * WORK(KEND1),LWRK1,IPRERI) 553 KEND1 = KFREE 554 LWRK1 = LFREE 555 ENDIF 556 NTOSYM = 1 557 ELSE 558 NTOSYM = NSYM 559 ENDIF !direct 560C 561 KENDSV = KEND1 562 LWRKSV = LWRK1 563C 564 ICDEL1 = 0 565 DO 100 ISYMD1 = 1,NTOSYM 566C 567 IF (DIRECT) THEN 568 IF (HERDIR) THEN 569 NTOT = MAXSHL 570 ELSE 571 NTOT = MXCALL 572 ENDIF 573 ELSE 574 NTOT = NBAS(ISYMD1) 575 ENDIF 576C 577 DO 110 ILLL = 1,NTOT 578C 579C--------------------------------------------- 580C If direct calculate the integrals. 581C--------------------------------------------- 582C 583 IF (DIRECT) THEN 584C 585 KEND1 = KENDSV 586 LWRK1 = LWRKSV 587C 588 DTIME = SECOND() 589 IF (HERDIR) THEN 590 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 591 & IPRINT) 592 ELSE 593 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 594 * WORK(KODCL1),WORK(KODCL2), 595 * WORK(KODBC1),WORK(KODBC2), 596 * WORK(KRDBC1),WORK(KRDBC2), 597 * WORK(KODPP1),WORK(KODPP2), 598 * WORK(KRDPP1),WORK(KRDPP2), 599 * WORK(KCCFB1),WORK(KINDXB), 600 * WORK(KEND1), LWRK1,IPRERI) 601 ENDIF 602 DTIME = SECOND() - DTIME 603 TIMHE2 = TIMHE2 + DTIME 604C 605 KRECNR = KEND1 606 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 607 LWRK1 = LWORK - KEND1 608 IF (LWRK1 .LT. 0) THEN 609 CALL QUIT('Insufficient core in CC_DEN_RCCD') 610 END IF 611C 612 ELSE 613 NUMDIS = 1 614 ENDIF !direct 615C 616C----------------------------------------------------- 617C Loop over number of distributions in disk. 618C----------------------------------------------------- 619C 620 DO 120 IDEL2 = 1,NUMDIS 621C 622 IF (DIRECT) THEN 623 IDEL = INDEXA(IDEL2) 624 IF (NOAUXB) THEN 625 IDUM = 1 626 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 627 END IF 628 ISYMD = ISAO(IDEL) 629 ELSE 630 IDEL = IBAS(ISYMD1) + ILLL 631 ISYMD = ISYMD1 632 ENDIF 633C 634C---------------------------------------- 635C Work space allocation two. 636C---------------------------------------- 637C 638 ISYDEN = ISYMD 639C 640 KD2IJG = KEND1 641 KD2AIG = KD2IJG + ND2IJG(ISYDEN) 642 KD2IAG = KD2AIG + ND2AIG(ISYDEN) 643 KD2ABG = KD2IAG + ND2AIG(ISYDEN) 644 KEND2 = KD2ABG + ND2ABG(ISYDEN) 645 LWRK2 = LWORK - KEND2 646C 647 IF (LWRK2 .LT. 0) THEN 648 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 649 CALL QUIT('Insufficient space for allocation '// 650 & '2 in CC_DEN_RCCD') 651 ENDIF 652C 653C------------------------------------------------------- 654C Initialize 4 two electron density arrays. 655C------------------------------------------------------- 656C 657 CALL DZERO(WORK(KD2IJG),ND2IJG(ISYDEN)) 658 CALL DZERO(WORK(KD2AIG),ND2AIG(ISYDEN)) 659 CALL DZERO(WORK(KD2IAG),ND2AIG(ISYDEN)) 660 CALL DZERO(WORK(KD2ABG),ND2ABG(ISYDEN)) 661C 662C----------------------------------------------------------------------------- 663C Calculate the RCCD two electron density d(pq,gamma;delta). 664C----------------------------------------------------------------------------- 665C 666 AUTIME = SECOND() 667C 668 CALL CC_DEN2_RCCD(WORK(KD2IJG),WORK(KD2AIG), 669 & WORK(KD2IAG),WORK(KD2ABG), 670 & WORK(KZ2PCK),WORK(KZ2AM), 671 & WORK(KT2AM),WORK(KT2AMT),WORK(KT2SQ), 672 & WORK(KZ2TIL),WORK(KXMAT), 673 & WORK(KYMAT),WORK(KONEAB),WORK(KONEAI), 674 & WORK(KONEIA),WORK(KLAMDH),1, 675 & WORK(KLAMDP),1,WORK(KEND2),LWRK2,IDEL,ISYMD) 676C 677 AUTIME = SECOND() - AUTIME 678 TIMDEN = TIMDEN + AUTIME 679C 680C------------------------------------------ 681C Work space allocation three. 682C------------------------------------------ 683C 684 ISYDIS = MULD2H(ISYMD,ISYMOP) 685C 686 KXINT = KEND2 687 KEND3 = KXINT + NDISAO(ISYDIS) 688 LWRK3 = LWORK - KEND3 689C 690 IF (LWRK3 .LT. 0) THEN 691 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 692 CALL QUIT('Insufficient space for allocation '// 693 & '3 in CC_DEN_PTFC') 694 ENDIF 695C 696C-------------------------------------------- 697C Read AO integral distribution. 698C-------------------------------------------- 699C 700 AUTIME = SECOND() 701 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 702 * WORK(KRECNR),DIRECT) 703 AUTIME = SECOND() - AUTIME 704 TIRDAO = TIRDAO + AUTIME 705C 706C---------------------------------------------------------------------- 707C Calculate integral intermediates needed for frozen core. 708C---------------------------------------------------------------------- 709C 710 IF (FROIMP) THEN 711 712 KDSRHF = KEND3 713 KOFOIN = KDSRHF + NDSRHF(ISYMD) 714 KDSFRO = KOFOIN + NOFROO(ISYDIS) 715 KSCRAI = KDSFRO + NDSFRO(ISYDIS) 716 KSCAIF = KSCRAI + NOFROO(ISYDIS) 717 KEND3 = KSCAIF + NCOFRF(ISYDIS) 718 LWRK3 = LWORK - KEND3 719C 720 IF (LWRK3 .LT. 0) THEN 721 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 722 CALL QUIT('Insufficient space for allocation '// 723 & 'in CC_DEN_PTFC') 724 ENDIF 725C 726 CALL DZERO(WORK(KSCRAI),NOFROO(ISYDIS)) 727 CALL DZERO(WORK(KSCAIF),NCOFRF(ISYDIS)) 728C 729C------------------------------------------------------------------------- 730C Transform one index in the integral batch to correlated. 731C------------------------------------------------------------------------- 732C 733 !(alp-bet,i,delta) 734 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KLAMDP), 735 * ISYMOP,WORK(KEND3),LWRK3,ISYDIS) 736C 737C--------------------------------------------------------------------- 738C Transform one index in the integral batch to frozen. 739C--------------------------------------------------------------------- 740C 741 !(alp-bet,i,delta) 742 CALL CC_GTOFRO(WORK(KXINT),WORK(KDSFRO),WORK(KCMOF), 743 * WORK(KEND3),LWRK3,ISYDIS) 744C 745C-------------------------------------------------------------- 746C Calculate integral batch (cor fro | cor del). 747C-------------------------------------------------------------- 748C 749 CALL CC_OFROIN(WORK(KDSRHF),WORK(KOFOIN),WORK(KCMOF), 750 * WORK(KEND3),LWRK3,ISYDIS) 751 752C 753C------------------------------------------------------------------ 754C Calculate terms to ai-part from KOFOIN integrals. 755C------------------------------------------------------------------ 756C 757 CALL CC_FRCOC1(WORK(KSCRAI),WORK(KOFOIN),ISYDIS) 758 759C 760C---------------------------------------------------------------- 761C Calculate exchange parts with KDSFRO integrals. 762C---------------------------------------------------------------- 763C 764 CALL CC_FRCOMI(WORK(KSCRAI),WORK(KSCAIF), 765 * WORK(KDSFRO),WORK(KCMOF), 766 * WORK(KEND3),LWRK3,ISYMD) 767 768C 769C---------------------------------------------------- 770C Calculate coulomb part of aI block. 771C---------------------------------------------------- 772C 773 CALL CC_FRCOF1(WORK(KSCAIF),WORK(KDSFRO),WORK(KCMOF), 774 * WORK(KEND3),LWRK3,ISYMD) 775 776C 777C----------------------------------------------------- 778C Calculate exchange part of aI block. 779C----------------------------------------------------- 780C 781 CALL CC_FRCOF2(WORK(KSCAIF),WORK(KDSRHF),WORK(KCMOF), 782 * WORK(KEND3),LWRK3,ISYMD) 783 784C 785C---------------------------------------------------------- 786C Dump intermediates to disc for later use. 787C---------------------------------------------------------- 788C 789 CALL CC_FRINDU(WORK(KSCRAI),WORK(KSCAIF),IDEL,ISYMD, 790 & LUN1,LUN2) 791 792 ENDIF !froimp 793C 794C------------------------------------------------------ 795C Start loop over second AO-index (gamma). 796C------------------------------------------------------ 797C 798 AUTIME = SECOND() 799C 800 DO 130 ISYMG = 1,NSYM 801C 802 ISYMPQ = MULD2H(ISYMG,ISYDEN) 803C 804 DO 140 G = 1,NBAS(ISYMG) 805C 806 KD2GIJ = KD2IJG + ID2IJG(ISYMPQ,ISYMG) 807 * + NMATIJ(ISYMPQ)*(G - 1) 808 KD2GAI = KD2AIG + ID2AIG(ISYMPQ,ISYMG) 809 * + NT1AM(ISYMPQ)*(G - 1) 810 KD2GAB = KD2ABG + ID2ABG(ISYMPQ,ISYMG) 811 * + NMATAB(ISYMPQ)*(G - 1) 812 KD2GIA = KD2IAG + ID2AIG(ISYMPQ,ISYMG) 813 * + NT1AM(ISYMPQ)*(G - 1) 814C 815C------------------------------------------------------------------------- 816C Work space allocation four. 817C Note: d2aob is only used for the ltest case!!!!!!!!!! 818C------------------------------------------------------------------------- 819C 820 KINTAO = KEND3 821 KD2AOB = KINTAO + N2BST(ISYMPQ) 822 KEND4 = KD2AOB + N2BST(ISYMPQ) 823 LWRK4 = LWORK - KEND4 824C 825 IF (LWRK4 .LT. 0) THEN 826 WRITE(LUPRI,*) 'Available:', LWORK 827 WRITE(LUPRI,*) 'Needed:', KEND4 828 CALL QUIT('Insufficient space in CC_DEN_PTFC') 829 ENDIF 830C 831 CALL DZERO(WORK(KD2AOB),N2BST(ISYMPQ)) 832C 833C------------------------------------------------------------- 834C Calculate frozen core contributions to d. 835C------------------------------------------------------------- 836C 837 IF (FROIMP) THEN 838C 839 KFD2IJ = KEND4 840 KFD2JI = KFD2IJ + NCOFRO(ISYMPQ) 841 KFD2AI = KFD2JI + NCOFRO(ISYMPQ) 842 KFD2IA = KFD2AI + NT1FRO(ISYMPQ) 843 KFD2II = KFD2IA + NT1FRO(ISYMPQ) 844 KEND4 = KFD2II + NFROFR(ISYMPQ) 845 LWRK4 = LWORK - KEND4 846C 847 IF (LWRK4 .LT. 0) THEN 848 WRITE (LUPRI,*) 'Available:', LWORK 849 WRITE (LUPRI,*) 'Needed:', KEND4 850 CALL QUIT( 851 * 'Insufficient work space in CC_DEN_PTFC') 852 ENDIF 853C 854 CALL DZERO(WORK(KFD2IJ),NCOFRO(ISYMPQ)) 855 CALL DZERO(WORK(KFD2JI),NCOFRO(ISYMPQ)) 856 CALL DZERO(WORK(KFD2AI),NT1FRO(ISYMPQ)) 857 CALL DZERO(WORK(KFD2IA),NT1FRO(ISYMPQ)) 858 CALL DZERO(WORK(KFD2II),NFROFR(ISYMPQ)) 859C 860C To calculate the contributions to d(pq,gam;del) 861C where at least one of the two indices p & q is frozen. 862C 863 CALL CC_FD2BL(WORK(KFD2II),WORK(KFD2IJ), 864 * WORK(KFD2JI),WORK(KFD2AI), 865 * WORK(KFD2IA),WORK(KONEIJ), 866 * WORK(KONEAB),WORK(KONEAI), 867 * WORK(KONEIA),WORK(KCMOF), 868 * WORK(KLAMDH),WORK(KLAMDP), 869 * WORK(KEND4),LWRK4,IDEL, 870 * ISYMD,G,ISYMG) 871C 872C ! calculate the contributions to D2AO from d(pq,gam;del) 873C ! where at least one of the two indices p & q is frozen 874 875 CALL CC_FD2AO(WORK(KD2AOB),WORK(KFD2II), 876 * WORK(KFD2IJ),WORK(KFD2JI), 877 * WORK(KFD2AI),WORK(KFD2IA), 878 * WORK(KCMOF),WORK(KLAMDH), 879 * WORK(KLAMDP),WORK(KEND4),LWRK4, 880 * ISYMPQ) 881C 882C Purpose: To calculate the contributions to d(pq,gam;del) where 883C gamma has been backtransformed from a frozen index. 884C 885 CALL CC_D2GAF(WORK(KD2GIJ),WORK(KD2GAB), 886 * WORK(KD2GAI),WORK(KD2GIA), 887 * WORK(KONEIJ),WORK(KONEAB), 888 * WORK(KONEAI),WORK(KONEIA), 889 * WORK(KCMOF),IDEL,ISYMD,G,ISYMG) 890C 891 ENDIF !froimp 892C 893C------------------------------------------------------- 894C Square up AO-integral distribution. 895C------------------------------------------------------- 896C 897 KOFFIN = KXINT + IDSAOG(ISYMG,ISYDIS) 898 * + NNBST(ISYMPQ)*(G - 1) 899C 900 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ, 901 * WORK(KINTAO)) 902C 903C--------------------------------------------------------------------------- 904C If energy test backtransform density fully to AO basis. 905C--------------------------------------------------------------------------- 906C 907 IF (LTSTE) THEN 908C 909 CALL CC_DENAO(WORK(KD2AOB),ISYMPQ, 910 * WORK(KD2GAI),WORK(KD2GAB), 911 * WORK(KD2GIJ),WORK(KD2GIA),ISYMPQ, 912 * WORK(KLAMDP),1,WORK(KLAMDH),1, 913 * WORK(KEND4),LWRK4) 914C 915C--------------------------------------------------------------------- 916C Add relaxation terms to set up effective density. 917C--------------------------------------------------------------------- 918C 919! IF (IOPT .EQ. 3) THEN 920C 921! ICON = 1 922! CALL CC_D2EFF(WORK(KD2AOB),G,ISYMG,IDEL, 923! * ISYMD,WORK(KKABAO), 924! * WORK(KDHFAO),ICON) 925C 926! ENDIF 927C 928C---------------------------------------------------------------------- 929C Calculate the 2 e- density contribution to E_rccd 930C---------------------------------------------------------------------- 931C 932 ECCSD2 = ECCSD2 + HALF*DDOT(N2BST(ISYMPQ), 933 * WORK(KD2AOB),1,WORK(KINTAO),1) 934C 935 END IF !ltste 936C 937C----------------------------------------------- 938C Work space allocation five. 939C----------------------------------------------- 940C 941 KIJINT = KEND4 942 KAIINT = KIJINT + NMATIJ(ISYMPQ) 943 KIAINT = KAIINT + NT1AM(ISYMPQ) 944 KABINT = KIAINT + NT1AM(ISYMPQ) 945 KABTIN = KABINT + NMATAB(ISYMPQ) 946 KIJTIN = KABTIN + NMATAB(ISYMPQ) 947 KD2TAB = KIJTIN + NMATIJ(ISYMPQ) 948 KD2TIJ = KD2TAB + NMATAB(ISYMPQ) 949 KEND5 = KD2TIJ + NMATIJ(ISYMPQ) 950 LWRK5 = LWORK - KEND5 951 IF (FROIMP) THEN 952 KIIFRO = KEND5 953 KIJFRO = KIIFRO + NFROFR(ISYMPQ) 954 KJIFRO = KIJFRO + NCOFRO(ISYMPQ) 955 KAIFRO = KJIFRO + NCOFRO(ISYMPQ) 956 KIAFRO = KAIFRO + NT1FRO(ISYMPQ) 957 KEND5 = KIAFRO + NT1FRO(ISYMPQ) 958 LWRK5 = LWORK - KEND5 959 ENDIF 960C 961 IF (LWRK5 .LT. 0) THEN 962 WRITE(LUPRI,*) 'Available:', LWORK 963 WRITE(LUPRI,*) 'Needed:', KEND5 964 CALL QUIT('Insufficient work space '// 965 & 'in CC_DEN_RCCD') 966 ENDIF 967C 968C---------------------------------------------------------------- 969C Transform 2 integral indices to MO-basis. 970C---------------------------------------------------------------- 971C 972 ISYM = ISYMPQ 973 CALL CCDINTMO(WORK(KIJINT),WORK(KIAINT), 974 * WORK(KABINT),WORK(KAIINT), 975 * WORK(KINTAO),WORK(KLAMDP), 976 * WORK(KLAMDH),WORK(KEND5), 977 * LWRK5,ISYM) 978C 979 IF (FROIMP) THEN 980C 981C Prepare integrals g_pq (gam,del) where one index is frozen 982C 983 ISYM = ISYMPQ 984 CALL CCIFROMO(WORK(KIJFRO),WORK(KJIFRO), 985 * WORK(KAIFRO),WORK(KIAFRO), 986 * WORK(KIIFRO),WORK(KINTAO), 987 * WORK(KLAMDP),WORK(KLAMDH), 988 * WORK(KCMOF),WORK(KEND5), 989 * LWRK5,ISYM) 990C 991 ENDIF !froimp 992C 993C----------------------------------------------------------------- 994C Set up transposed integrals and densities. 995C----------------------------------------------------------------- 996C 997 CALL DCOPY(NMATAB(ISYMPQ),WORK(KABINT),1, 998 * WORK(KABTIN),1) 999 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KIJINT),1, 1000 * WORK(KIJTIN),1) 1001 CALL DCOPY(NMATAB(ISYMPQ),WORK(KD2GAB),1, 1002 * WORK(KD2TAB),1) 1003 CALL DCOPY(NMATIJ(ISYMPQ),WORK(KD2GIJ),1, 1004 * WORK(KD2TIJ),1) 1005C 1006 CALL CC_EITR(WORK(KABTIN),WORK(KIJTIN), 1007 * WORK(KEND5),LWRK5,ISYMPQ) 1008 CALL CC_EITR(WORK(KD2TAB),WORK(KD2TIJ), 1009 * WORK(KEND5),LWRK5,ISYMPQ) 1010C 1011C------------------------------------------------------------------- 1012C Calculate 2 e- contribution to Zeta-Kappa-0. 1013C------------------------------------------------------------------- 1014C 1015 ISYM = ISYMPQ 1016 IF (IOPT.EQ.2) THEN 1017 1018 KOFFIJ = 1 1019 CALL CC2_ETIJ(ZKDIA(KOFFIJ), 1020 & WORK(KIJINT),WORK(KAIINT), 1021 & WORK(KIAINT),WORK(KABINT), 1022 & WORK(KD2GIJ),WORK(KD2GAI), 1023 & WORK(KD2GIA),WORK(KD2GAB), 1024 & WORK(KT1AM),WORK(KEND5),LWRK5, 1025 & ISYM) 1026 1027 KOFFAB = NMATIJ(1) + 1 1028 CALL CC2_ETAB(ZKDIA(KOFFAB), 1029 & WORK(KIJINT),WORK(KAIINT), 1030 * WORK(KIAINT),WORK(KABINT), 1031 * WORK(KD2GIJ),WORK(KD2GAI), 1032 * WORK(KD2GIA),WORK(KD2GAB), 1033 * WORK(KT1AM),WORK(KEND5),LWRK5, 1034 * ISYM) 1035 END IF !iopt2 1036 1037 CALL CCDENZK0(ETAAI,WORK(KIJINT),WORK(KAIINT), 1038 * WORK(KIAINT),WORK(KABINT), 1039 * WORK(KD2GIJ),WORK(KD2GAI), 1040 * WORK(KD2GIA),WORK(KD2GAB), 1041 * WORK(KIJTIN),WORK(KABTIN), 1042 * WORK(KD2TIJ),WORK(KD2TAB), 1043 * WORK(KT1AM),WORK(KEND5),LWRK5, 1044 * ISYM) 1045C 1046 IF (FROIMP) THEN 1047C 1048 ISYM = ISYMPQ 1049 ! 1050 ! contributions to eta_ai from loop over frozen 1051 ! 1052 CALL CCFRETAI(ETAAI,WORK(KIJFRO), 1053 * WORK(KJIFRO),WORK(KAIFRO), 1054 * WORK(KIAFRO),WORK(KFD2IJ), 1055 * WORK(KFD2JI),WORK(KFD2AI), 1056 * WORK(KFD2IA),WORK(KT1AM), 1057 * WORK(KEND5),LWRK5,ISYM) 1058C 1059C----------------------------------------------------------------------- 1060C Calculate two-electron contribution to right- 1061C hand-side of correlated-frozen multipliers. 1062C----------------------------------------------------------------------- 1063C 1064 ICON = 2 1065 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 1066 * + 2*NT1FRO(1) + 1 1067 ! 1068 ! eta_iJ 1069 ! 1070 CALL CC_ETIJF(ZKDIA(KOFRES),WORK(KD2GIJ), 1071 * WORK(KD2GAB),WORK(KD2GAI), 1072 * WORK(KD2GIA),WORK(KD2TIJ), 1073 * WORK(KFD2II),WORK(KFD2IJ), 1074 * WORK(KFD2JI),WORK(KFD2AI), 1075 * WORK(KFD2IA),WORK(KIJINT), 1076 * WORK(KAIINT),WORK(KIAINT), 1077 * WORK(KIJTIN),WORK(KABTIN), 1078 * WORK(KIJFRO),WORK(KJIFRO), 1079 * WORK(KAIFRO),WORK(KIAFRO), 1080 * WORK(KIIFRO),WORK(KT1AM), 1081 * WORK(KEND5),LWRK5,ISYM,ICON) 1082C 1083C----------------------------------------------------------------------- 1084C Calculate two-electron contribution to right- 1085C hand-side of virtual-frozen multipliers. 1086C----------------------------------------------------------------------- 1087C 1088 ICON = 2 1089 KOFRE = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1090 ! 1091 ! eta_aI 1092 ! 1093 CALL CC_ETAIF(ZKDIA(KOFRE),WORK(KD2GAB), 1094 * WORK(KD2GAI),WORK(KD2GIA), 1095 * WORK(KD2TIJ),WORK(KD2TAB), 1096 * WORK(KFD2II),WORK(KFD2IJ), 1097 * WORK(KFD2JI),WORK(KFD2AI), 1098 * WORK(KFD2IA),WORK(KIJINT), 1099 * WORK(KABINT),WORK(KAIINT), 1100 * WORK(KIAINT),WORK(KABTIN), 1101 * WORK(KIJFRO),WORK(KJIFRO), 1102 * WORK(KAIFRO),WORK(KIAFRO), 1103 * WORK(KIIFRO),WORK(KT1AM), 1104 * WORK(KEND5),LWRK5,ISYM,ICON) 1105C 1106C----------------------------------------------------------------------- 1107C Calculate two-electron contribution to right- 1108C hand-side of frozen-frozen multipliers. 1109C----------------------------------------------------------------------- 1110C 1111 ICON = 2 1112 KOFRES = 2*NT1AMX + NMATIJ(1) + NMATAB(1) 1113 * + 2*NT1FRO(1) + 2*NCOFRO(1) + 1 1114c ! 1115c ! eta_ab contribs from loop over frozen I 1116c ! 1117 CALL CCFRETAB(ZKDIA,WORK(KIJFRO), 1118 * WORK(KJIFRO),WORK(KAIFRO), 1119 * WORK(KIAFRO),WORK(KFD2IJ), 1120 * WORK(KFD2JI),WORK(KFD2AI), 1121 * WORK(KFD2IA),WORK(KT1AM), 1122 * WORK(KEND5),LWRK5,ISYM) 1123 ! 1124 ! eta_ij contribs from loop over frozen L 1125 ! 1126 CALL CCFRETIJ(ZKDIA,WORK(KIJFRO), 1127 * WORK(KJIFRO),WORK(KAIFRO), 1128 * WORK(KIAFRO),WORK(KFD2IJ), 1129 * WORK(KFD2JI),WORK(KFD2AI), 1130 * WORK(KFD2IA),WORK(KT1AM), 1131 * WORK(KEND5),LWRK5,ISYM) 1132c 1133C 1134 ENDIF !froimp 1135C 1136 140 CONTINUE 1137 130 CONTINUE 1138C 1139 AUTIME = SECOND() - AUTIME 1140 TIMRES = TIMRES + AUTIME 1141C 1142 120 CONTINUE 1143 110 CONTINUE 1144 100 CONTINUE 1145 1146C 1147C----------------------- 1148C Regain work space. 1149C----------------------- 1150C 1151 KEND1 = KENDS2 1152 LWRK1 = LWRKS2 1153C 1154 if (locdbg) then 1155 KOFFIJ = 1 1156 KOFFAB = 1 + NMATIJ(1) 1157 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1158 KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1159 write(lupri,*) ' ' 1160 write(lupri,*) 'Before call to CCSD_ZKBLO' 1161 xtest=ddot(nmatij(1),zkdia(koffij),1,zkdia(koffij),1) 1162 write(lupri,*) 'Norm of ETAIJ: ', xtest 1163 xtest=ddot(2*NCOFRO(1),zkdia(kofifj),1,zkdia(kofifj),1) 1164 write(lupri,*) 'Norm of ETIFJ: ', xtest 1165 xtest=ddot(nmatab(1),zkdia(koffab),1,zkdia(koffab),1) 1166 write(lupri,*) 'Norm of ETAAB: ', xtest 1167 xtest=ddot(2*nt1amx,etaai(1),1,etaai(1),1) 1168 write(lupri,*) 'Norm of ETAAI: ', xtest 1169 xtest=ddot(2*nt1fro(1),zkdia(kofafi),1,zkdia(kofafi),1) 1170 write(lupri,*) 'Norm of ETAFI: ', xtest 1171 call flshfo(lupri) 1172 end if !locdbg 1173C 1174C------------------------------------------------------------------------ 1175C Calculate the ZK0(ij) and ZK0(ab) blocks (to be used to correct eta_ai) 1176C from eta_ij/e_i-e_j and eta_ab/e_a-e_b 1177C------------------------------------------------------------------------ 1178C 1179 IF (IOPT.EQ.2) THEN 1180 1181 CALL CCSD_ZKBLO(ZKDIA,WORK(KEND1),LWRK1) 1182 write(lupri,*) 'I came out of ccsd_zkblo' 1183 1184 if (locdbg) then 1185 KOFFIJ = 1 1186 KOFFAB = 1 + NMATIJ(1) 1187 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1188 KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1189 write(lupri,*) 'after call to CCSD_ZKBLO' 1190 xtest=ddot(nmatij(1),zkdia(koffij),1,zkdia(koffij),1) 1191 write(lupri,*) 'Norm of ETAIJ: ', xtest 1192 xtest=ddot(2*NCOFRO(1),zkdia(kofifj),1,zkdia(kofifj),1) 1193 write(lupri,*) 'Norm of ETIFJ: ', xtest 1194 xtest=ddot(nmatab(1),zkdia(koffab),1,zkdia(koffab),1) 1195 write(lupri,*) 'Norm of ETAAB: ', xtest 1196 xtest=ddot(2*nt1amx,etaai(1),1,etaai(1),1) 1197 write(lupri,*) 'Norm of ETAAI: ', xtest 1198 xtest=ddot(2*nt1fro(1),zkdia(kofafi),1,zkdia(kofafi),1) 1199 write(lupri,*) 'Norm of ETAFI: ', xtest 1200 call flshfo(lupri) 1201 end if 1202 END IF !locdbg 1203C 1204C------------------------------------------------------------------------ 1205C Add the diagonal ZK0(ii) and ZK0(aa) elements in the proper places 1206C------------------------------------------------------------------------ 1207C 1208! if (.false.) then 1209! 1210! IF (LTSTE) THEN 1211! 1212! !multiply by epsilon_p and sum over p to get the energy 1213! 1214! KFOCKDIA = KEND1 1215! KEND1 = KFOCKDIA + NORBTS 1216! LWRK1 = LWORK-KEND1 1217C 1218C------------------------------------- 1219C Read canonical orbital energies. 1220C------------------------------------- 1221C 1222! CALL DZERO(WORK(KFOCKDIA),NORBTS) 1223! CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1224! & .FALSE.) 1225! REWIND (LUSIFC) 1226C 1227! CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1228! READ (LUSIFC) 1229! READ (LUSIFC) (WORK(KFOCKDIA + I - 1), I = 1,NORBTS) 1230C 1231! CALL GPCLOSE(LUSIFC,'KEEP') 1232C 1233C---------------------------------------------------------------- 1234C Change symmetry ordering of the canonical orbital energies. 1235C---------------------------------------------------------------- 1236C 1237! IF (FROIMP) 1238! * CALL CCSD_DELFRO(WORK(KFOCKDIA),WORK(KEND1),LWRK1) 1239C 1240! CALL FOCK_REORDER(WORK(KFOCKDIA),WORK(KEND1),LWRK1) 1241C 1242C-------------------------------------------- 1243C Calculate sum_p kappabar_pp * epsilon_p 1244C Occupied block: 1245C-------------------------------------------- 1246C 1247! SKAPEP = DDOT(NORBT,WORK(KKAPII),1,WORK(KFOCKDIA),1) 1248! END IF 1249! END IF 1250c 1251! end if 1252! END IF 1253C 1254C------------------------------------------------ 1255C Correct Eta_ai with ZK0(ij) and ZK0(ab) blocks 1256C------------------------------------------------ 1257C-tbp: 1258c CALL AROUND('Eta-kappa-bar-0-ai vector before modification') 1259c do ISYM = 1,NSYM 1260c WRITE(LUPRI,*) ' ' 1261c WRITE(LUPRI,444) 'Sub-symmetry block number:', ISYM 1262c WRITE(LUPRI,555) '--------------------------' 1263c KOFF = IALLAI(ISYM,ISYM) + 1 1264c write(lupri,'(A,1P,D25.16)') 'Norm=', 1265c * sqrt(ddot(NVIR(ISYM)*NRHFS(ISYM),ETAAI(KOFF),1, 1266c * ETAAI(KOFF),1)) 1267c CALL OUTPUT(ETAAI(KOFF),1,NVIR(ISYM),1,NRHFS(ISYM), 1268c * NVIR(ISYM),NRHFS(ISYM),1,LUPRI) 1269c IF ((NVIR(ISYM) .EQ. 0) .OR. (NRHFS(ISYM) .EQ. 0)) THEN 1270c WRITE(LUPRI,*) 'This sub-symmetry is empty' 1271c ENDIF 1272c end do 1273 1274C 1275 IF (IOPT.EQ.2) THEN 1276 !SONIA: NB!! I removed the froimp stuff 1277 !IF (LPRNCC) write(lupri,*)'CC_DEN_RCCD using CCETACOR' 1278 call flshfo(lupri) 1279 CALL CCETACOR(ETAAI,ZKDIA,WORK(KEND1),LWRK1) 1280 !IF (LPRNCC) write(lupri,*)'CC_DEN_RCCD out of etacor' 1281 call flshfo(lupri) 1282 1283 KOFFIJ = 1 1284 KOFFAB = NMATIJ(1) + 1 1285 KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1286 KOFIFJ = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 2*NT1FRO(1) + 1 1287 1288 END IF 1289 1290 XZKDIA = DDOT(LUNGO,ZKDIA,1,ZKDIA,1) 1291! IF (LPRNCC) 1292! & WRITE(LUPRI,*) 'CC_DEN_RCCD: ZKDIA before CCETACOR', XZKDIA 1293! call flshfo(lupri) 1294C 1295C--------------------- 1296C Reorder results. 1297C it does nothing if it's NOT froimp 1298C--------------------- 1299C 1300! KOFAFI = 2*NT1AMX + NMATIJ(1) + NMATAB(1) + 1 1301! CALL CC_ETARE(ETAAI,ZKDIA(KOFAFI),WORK(KEND1),LWRK1) 1302C 1303C--------------------------------- 1304C Write out eta-ai and eta-aI. 1305C--------------------------------- 1306C 1307 IF ((IPRINT .GT. 20).OR.(LOCDBG)) THEN 1308C 1309 CALL AROUND('Eta-kappa-bar-0-ai vector exiting CC_DEN_RCCD') 1310C 1311 DO 20 ISYM = 1,NSYM 1312C 1313 WRITE(LUPRI,*) ' ' 1314 WRITE(LUPRI,444) 'Sub-symmetry block number:', ISYM 1315 WRITE(LUPRI,555) '--------------------------' 1316 444 FORMAT(3X,A26,2X,I1) 1317 555 FORMAT(3X,A25) 1318C 1319 KOFF = IALLAI(ISYM,ISYM) + 1 1320 write(lupri,'(A,1P,D25.16)') 'Norm=', 1321 * sqrt(ddot(NVIR(ISYM)*NRHFS(ISYM),ETAAI(KOFF),1, 1322 * ETAAI(KOFF),1)) 1323 CALL OUTPUT(ETAAI(KOFF),1,NVIR(ISYM),1,NRHFS(ISYM), 1324 * NVIR(ISYM),NRHFS(ISYM),1,LUPRI) 1325C 1326 IF ((NVIR(ISYM) .EQ. 0) .OR. (NRHFS(ISYM) .EQ. 0)) THEN 1327 WRITE(LUPRI,*) 'This sub-symmetry is empty' 1328 ENDIF 1329C 1330 20 CONTINUE 1331 ENDIF 1332C 1333 IF ((IPRINT .GT. 9).OR.(LOCDBG)) THEN 1334 ETAKAN = DDOT(NALLAI(1),ETAAI,1,ETAAI,1) 1335 WRITE(LUPRI,*) 'CC_DEN_RCCD ' 1336 WRITE(LUPRI,*) 'Norm of occupied-virtual block:', ETAKAN 1337 call flshfo(lupri) 1338 ENDIF 1339C 1340C------------------------------ 1341C Close intermediate files. 1342C------------------------------ 1343C 1344 IF (FROIMP) THEN 1345 CALL WCLOSE2(LUN1,NAME1,'KEEP') 1346 CALL WCLOSE2(LUN2,NAME2,'KEEP') 1347 ENDIF 1348C 1349C----------------------- 1350C Write out timings. 1351C----------------------- 1352C 1353 99 TIMTOT = SECOND() - TIMTOT 1354C 1355 IF (IPRINT .GT. 3) THEN 1356 WRITE (LUPRI,*) ' ' 1357 WRITE (LUPRI,*) 'Requested density dependent '// 1358 & 'quantities calculated' 1359 WRITE (LUPRI,*) 'Total time used in CC_DEN_RCCD:', TIMTOT 1360 ENDIF 1361 IF (IPRINT .GT. 9) THEN 1362 WRITE (LUPRI,*) 'Time used for setting up 2 e- density:',TIMDEN 1363 WRITE (LUPRI,*) 'Time used for contraction with integrals:', 1364 & TIMRES 1365 WRITE (LUPRI,*) 'Time used for reading 2 e- AO-integrals:', 1366 & TIRDAO 1367 WRITE (LUPRI,*) 'Time used for calculating 2 e- AO-integrals:', 1368 * TIMHE2 1369 WRITE (LUPRI,*) 'Time used for 1 e- density & intermediates:', 1370 * TIMONE 1371 ENDIF 1372C 1373C--------------------------------------------------------------- 1374C If energy test add nuclear nuclear repulsion energy and write out E-ccsd. 1375C--------------------------------------------------------------- 1376C 1377 IF (ltste) THEN 1378C 1379 ECCSD = ECCSD1 + ECCSD2 + POTNUC 1380C 1381 !IF (LPRNCC) THEN 1382 WRITE(LUPRI,*) ' ' 1383 WRITE(LUPRI,*) 'RPA energy constructed' 1384 WRITE(LUPRI,*) 'from density matrices:' 1385 !IF (CCSD) WRITE(LUPRI,*) 'CCSD-(type) energy:', ECCSD 1386 WRITE(LUPRI,'(A,f15.10)') 'H1 energy, ECCSD1(type) = ',ECCSD1 1387 WRITE(LUPRI,'(A,f15.10)') 'H2 energy, ECCSD2(type) = ',ECCSD2 1388 !WRITE(lupri,'(A,f15.10)') 'sum_p e_p kbar_pp = ',SKAPEP 1389 WRITE(LUPRI,'(A,f15.10)') 'Nuc. Pot. energy = ',POTNUC 1390 WRITE(lupri,'(A,f15.10)') 'CCSD energy ? = ', 1391 & ECCSD1+ECCSD2+ POTNUC 1392 call flshfo(lupri) 1393 WRITE(lupri,'(A,f15.10)') 'CCSD energy (2) ? = ',ECCSD 1394 !WRITE(LUPRI,*) 'OBS POTNUC is missing!!! ' 1395 !ENDIF 1396C 1397 ENDIF 1398C---------------------------------------------------------------------- 1399 CALL QEXIT('CC_DEN_RCCD') 1400 RETURN 1401 END 1402C---------------------------------------------------------------------- 1403 1404 1405C /* Deck rccd_etij */ 1406 SUBROUTINE RCCD_ETIJ(ETAIJ,XINTIJ,XINTAI,XINTIA,XINTAB,DIJ,DAI, 1407 * DIA,DAB,WORK,LWORK,ISYM) 1408C 1409C Written by S. Coriani 2010 1410C 1411C Version: 1.0 1412C 1413C Purpose: To set up the right hand side of the equation for 1414C zeta-kappa-0_ij (ETAIJ) from MO-integrals (XIN*) and RCCD 1415C densities (D*) 1416C Note that due to the symmetry in the formulas, this 1417C routine is able to handle both the one- and the two- 1418C electron contributions! 1419C ISYM is the symmetry of both the density and the 1420C integrals! 1421C 1422#include "implicit.h" 1423 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1424 DIMENSION ETAIJ(*), XINTIJ(*), XINTAI(*), XINTIA(*), XINTAB(*) 1425 DIMENSION DIJ(*), DAI(*), DIA(*), DAB(*), WORK(LWORK) 1426#include "priunit.h" 1427#include "ccorb.h" 1428#include "ccsdsym.h" 1429#include "cclr.h" 1430C 1431 CALL QENTER('RCCD_ETIJ') 1432C 1433 DO 100 ISYMI = 1,NSYM 1434C 1435C---------------------------------------------------------------- 1436C Calculate direct terms to eta_ij. 1437C---------------------------------------------------------------- 1438C 1439 ISYMJ = ISYMI 1440 ISYMK = MULD2H(ISYMI,ISYM) 1441 ISYMC = MULD2H(ISYMI,ISYM) 1442C 1443 KOFFRE = IMATIJ(ISYMI,ISYMJ) + 1 1444C 1445 NTOTRE = MAX(NRHF(ISYMI),1) 1446 NTOTI = MAX(NRHF(ISYMI),1) 1447 NTOTJ = MAX(NRHF(ISYMJ),1) 1448 NTOTK = MAX(NRHF(ISYMK),1) 1449 NTOTC = MAX(NVIR(ISYMC),1) 1450C 1451 KOFF1 = IMATIJ(ISYMK,ISYMI) + 1 1452 KOFF2 = IMATIJ(ISYMK,ISYMJ) + 1 1453 KOFF3 = IMATIJ(ISYMI,ISYMK) + 1 1454 KOFF4 = IMATIJ(ISYMJ,ISYMK) + 1 1455C 1456 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),ONE, 1457 * XINTIJ(KOFF1),NTOTK,DIJ(KOFF2),NTOTK,ONE, 1458 * ETAIJ(KOFFRE),NTOTRE) 1459C 1460 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),ONE, 1461 * XINTIJ(KOFF3),NTOTI,DIJ(KOFF4),NTOTJ,ONE, 1462 * ETAIJ(KOFFRE),NTOTRE) 1463C 1464 CALL DGEMM('N','T',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),-ONE, 1465 * DIJ(KOFF3),NTOTI,XINTIJ(KOFF4),NTOTJ,ONE, 1466 * ETAIJ(KOFFRE),NTOTRE) 1467C 1468 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NRHF(ISYMK),-ONE, 1469 * DIJ(KOFF1),NTOTK,XINTIJ(KOFF2),NTOTK,ONE, 1470 * ETAIJ(KOFFRE),NTOTRE) 1471C 1472 KOFF5 = IT1AM(ISYMC,ISYMI) + 1 1473 KOFF6 = IT1AM(ISYMC,ISYMJ) + 1 1474C 1475 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 1476 * XINTAI(KOFF5),NTOTC,DAI(KOFF6),NTOTC,ONE, 1477 * ETAIJ(KOFFRE),NTOTRE) 1478C 1479 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),ONE, 1480 * XINTIA(KOFF5),NTOTC,DIA(KOFF6),NTOTC,ONE, 1481 * ETAIJ(KOFFRE),NTOTRE) 1482C 1483 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 1484 * DIA(KOFF5),NTOTC,XINTIA(KOFF6),NTOTC,ONE, 1485 * ETAIJ(KOFFRE),NTOTRE) 1486C 1487 CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),-ONE, 1488 * DAI(KOFF5),NTOTC,XINTAI(KOFF6),NTOTC,ONE, 1489 * ETAIJ(KOFFRE),NTOTRE) 1490C 1491 100 CONTINUE 1492C 1493 CALL QEXIT('RCCD_ETIJ') 1494C 1495 RETURN 1496 END 1497C /* Deck rccd_etab */ 1498 SUBROUTINE RCCD_ETAB(ETAAB,XINTIJ,XINTAI,XINTIA,XINTAB,DIJ,DAI, 1499 * DIA,DAB,WORK,LWORK,ISYM) 1500C 1501C Written by S. Coriani 2010 1502C 1503C Version: 1.0 1504C 1505C Purpose: To set up the right hand side of the equation for 1506C zeta-kappa-0_ab (ETAAB) from MO-integrals (XIN*) & RCCD 1507C densities (D*) 1508C Note that due to the symmetry in the formulas, this 1509C routine is able to handle both the one- and the two- 1510C electron contributions! 1511C ISYM is the symmetry of both the density and the 1512C integrals! 1513C 1514#include "implicit.h" 1515 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1516 DIMENSION ETAAB(*), XINTIJ(*), XINTAI(*), XINTIA(*), XINTAB(*) 1517 DIMENSION DIJ(*), DAI(*), DIA(*), DAB(*), WORK(LWORK) 1518#include "priunit.h" 1519#include "ccorb.h" 1520#include "ccsdsym.h" 1521#include "cclr.h" 1522C 1523 CALL QENTER('RCCD_ETAB') 1524C 1525 DO 100 ISYMA = 1,NSYM 1526C 1527C---------------------------------------------------------------- 1528C Calculate direct terms to eta_ab. 1529C---------------------------------------------------------------- 1530C 1531 ISYMB = ISYMA 1532 ISYMK = MULD2H(ISYMA,ISYM) 1533 ISYMC = MULD2H(ISYMA,ISYM) 1534C 1535 KOFFRE = IMATAB(ISYMA,ISYMB) + 1 1536C 1537 NTOTRE = MAX(NVIR(ISYMA),1) 1538 NTOTA = MAX(NVIR(ISYMA),1) 1539 NTOTB = MAX(NVIR(ISYMB),1) 1540 NTOTC = MAX(NVIR(ISYMC),1) 1541 NTOTK = MAX(NRHF(ISYMK),1) 1542C 1543 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 1544 KOFF2 = IT1AM(ISYMB,ISYMK) + 1 1545C 1546 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 1547 * XINTIA(KOFF1),NTOTA,DIA(KOFF2),NTOTB,ONE, 1548 * ETAAB(KOFFRE),NTOTRE) 1549C 1550 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 1551 * DAI(KOFF1),NTOTA,XINTAI(KOFF2),NTOTB,ONE, 1552 * ETAAB(KOFFRE),NTOTRE) 1553C 1554 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),-ONE, 1555 * DIA(KOFF1),NTOTA,XINTIA(KOFF2),NTOTB,ONE, 1556 * ETAAB(KOFFRE),NTOTRE) 1557C 1558 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),ONE, 1559 * XINTAI(KOFF1),NTOTA,DAI(KOFF2),NTOTB,ONE, 1560 * ETAAB(KOFFRE),NTOTRE) 1561C 1562 KOFF3 = IMATAB(ISYMC,ISYMA) + 1 1563 KOFF4 = IMATAB(ISYMC,ISYMB) + 1 1564 KOFF5 = IMATAB(ISYMA,ISYMC) + 1 1565 KOFF6 = IMATAB(ISYMB,ISYMC) + 1 1566C 1567 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),ONE, 1568 * XINTAB(KOFF3),NTOTC,DAB(KOFF4),NTOTC,ONE, 1569 * ETAAB(KOFFRE),NTOTRE) 1570C 1571 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),-ONE, 1572 * DAB(KOFF5),NTOTA,XINTAB(KOFF6),NTOTB,ONE, 1573 * ETAAB(KOFFRE),NTOTRE) 1574C 1575 CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),-ONE, 1576 * DAB(KOFF3),NTOTC,XINTAB(KOFF4),NTOTC,ONE, 1577 * ETAAB(KOFFRE),NTOTRE) 1578C 1579 CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NVIR(ISYMC),ONE, 1580 * XINTAB(KOFF5),NTOTA,DAB(KOFF6),NTOTB,ONE, 1581 * ETAAB(KOFFRE),NTOTRE) 1582C 1583 100 CONTINUE 1584C 1585C 1586 CALL QEXIT('RCCD_ETAB') 1587C 1588 RETURN 1589 END 1590C /* Deck Rccdenzk0 */ 1591 SUBROUTINE RCCDENZK0(ETAKA,XINTIJ,XINTAI,XINTIA,XINTAB,DIJ,DAI, 1592 * DIA,DAB,XIJT,XABT,DIJT,DABT, 1593 * WORK,LWORK,ISYM) 1594C 1595C Written by Sonia Halkier 28/1 - 2011 - based on CCDENZK0 1596C 1597C Version: 1.0 1598C 1599C Purpose: To set up the right hand side of the equation for 1600C zeta-kappa-0 (ETAKA) from MO-integrals (XI*), Coupled 1601C Cluster densities (D*) and t1-amplitudes (T1AM)! 1602C Note that due to the symmetry in the formulas, this 1603C routine is able to handle both the one- and the two- 1604C electron contributions! 1605C ISYM is the symmetry of both the density and the 1606C integrals! 1607#include "implicit.h" 1608 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 1609 DIMENSION ETAKA(*), XINTIJ(*), XINTAI(*), XINTIA(*), XINTAB(*) 1610 DIMENSION DIJ(*), DAI(*), DIA(*), DAB(*), XIJT(*), XABT(*) 1611 DIMENSION DIJT(*), DABT(*), WORK(LWORK) 1612#include "priunit.h" 1613#include "ccorb.h" 1614#include "ccsdsym.h" 1615C 1616 CALL QENTER('RCCDENZK0') 1617C 1618 DO 100 ISYMA = 1,NSYM 1619C 1620C------------------------------------------------------- 1621C Calculate terms originating from [H(t1),E(ai)]. 1622C------------------------------------------------------- 1623C 1624 ISYMI = ISYMA 1625 ISYMB = MULD2H(ISYMA,ISYM) 1626 ISYMJ = MULD2H(ISYMA,ISYM) 1627C 1628 KOFFRE = IT1AM(ISYMA,ISYMI) + 1 1629C 1630 NTOTRE = MAX(NVIR(ISYMA),1) 1631 NTOTI = MAX(NRHF(ISYMI),1) 1632 NTOTB = MAX(NVIR(ISYMB),1) 1633 NTOTJ = MAX(NRHF(ISYMJ),1) 1634C 1635 KOFF1 = IMATAB(ISYMA,ISYMB) + 1 1636 KOFF2 = IT1AM(ISYMB,ISYMI) + 1 1637C 1638 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 1639 * ONE,XABT(KOFF1),NTOTRE,DAI(KOFF2),NTOTB,ONE, 1640 * ETAKA(KOFFRE),NTOTRE) 1641C 1642 KOFF3 = IMATAB(ISYMA,ISYMB) + 1 1643 KOFF4 = IT1AM(ISYMB,ISYMI) + 1 1644 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 1645 * -ONE,DAB(KOFF3),NTOTRE,XINTIA(KOFF4),NTOTB,ONE, 1646 * ETAKA(KOFFRE),NTOTRE) 1647C 1648 KOFF5 = IT1AM(ISYMA,ISYMJ) + 1 1649 KOFF6 = IMATIJ(ISYMJ,ISYMI) + 1 1650C 1651 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 1652 * ONE,XINTIA(KOFF5),NTOTRE,DIJ(KOFF6),NTOTJ,ONE, 1653 * ETAKA(KOFFRE),NTOTRE) 1654C 1655 KOFF7 = IT1AM(ISYMA,ISYMJ) + 1 1656 KOFF8 = IMATIJ(ISYMJ,ISYMI) + 1 1657C 1658 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 1659 * -ONE,DAI(KOFF7),NTOTRE,XIJT(KOFF8),NTOTJ,ONE, 1660 * ETAKA(KOFFRE),NTOTRE) 1661C 1662C------------------------------------------------------- 1663C Calculate terms originating from [H(t1),E(ia)]. 1664C------------------------------------------------------- 1665C 1666 KOFF9 = IMATAB(ISYMA,ISYMB) + 1 1667 KOFF10 = IT1AM(ISYMB,ISYMI) + 1 1668C 1669 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 1670 * -ONE,DABT(KOFF9),NTOTRE,XINTAI(KOFF10),NTOTB, 1671 * ONE,ETAKA(KOFFRE),NTOTRE) 1672C 1673 KOFF11 = IMATAB(ISYMA,ISYMB) + 1 1674 KOFF12 = IT1AM(ISYMB,ISYMI) + 1 1675C 1676 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 1677 * ONE,XINTAB(KOFF11),NTOTRE,DIA(KOFF12),NTOTB, 1678 * ONE,ETAKA(KOFFRE),NTOTRE) 1679C 1680 KOFF13 = IT1AM(ISYMA,ISYMJ) + 1 1681 KOFF14 = IMATIJ(ISYMJ,ISYMI) + 1 1682 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 1683 * -ONE,DIA(KOFF13),NTOTRE,XINTIJ(KOFF14),NTOTJ, 1684 * ONE,ETAKA(KOFFRE),NTOTRE) 1685C 1686 KOFF15 = IT1AM(ISYMA,ISYMJ) + 1 1687 KOFF16 = IMATIJ(ISYMJ,ISYMI) + 1 1688C 1689 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 1690 * ONE,XINTAI(KOFF15),NTOTRE,DIJT(KOFF16),NTOTJ, 1691 * ONE,ETAKA(KOFFRE),NTOTRE) 1692C 1693 100 CONTINUE 1694C 1695 CALL QEXIT('RCCDENZK0') 1696C 1697 RETURN 1698 END 1699 1700C /* Deck rccd_zkblo */ 1701 SUBROUTINE RCCD_ZKBLO(ZKDIA,WORK,LWORK) 1702C 1703C Written by Sonia Halkier 2011 1704C 1705C Version: 1.0 1706C 1707C Purpose: To calculate the ab & ij parts of the CCSD kappa-bar-0, 1708C from right-hand-sides (ZKDIA on input) and canonical 1709C orbital energies. 1710C Control numerical instabilities. Sonia, 2002 1711C 1712#include "implicit.h" 1713#include "dummy.h" 1714 PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 1715 PARAMETER(FOUR = 4.0D0, EIGHT=8.0D0) 1716 PARAMETER (THRDNM = 1.0D-12) 1717 DIMENSION ZKDIA(*), WORK(LWORK) 1718#include "priunit.h" 1719#include "maxorb.h" 1720#include "ccorb.h" 1721#include "iratdef.h" 1722#include "inftap.h" 1723#include "ccsdsym.h" 1724#include "ccsdio.h" 1725#include "ccsdinp.h" 1726C 1727 CALL QENTER('RCCD_ZKBLO') 1728C-tbp: 1729 write(lupri,*) 1730 write(lupri,*) 'RCCD_ZKBLO WARNING: I''m not numerically stable!' 1731 write(lupri,*) 'RCCD_ZKBLO WARNING: use CCSD_ZKBLO instead...' 1732 write(lupri,*) 1733 call flshfo(lupri) 1734C--------------------------- 1735C Work space allocation. 1736C--------------------------- 1737C 1738 KFOCKD = 1 1739 KETAIJ = KFOCKD + NORBTS 1740 KETAAB = KETAIJ + NMATIJ(1) 1741 KEND1 = KETAAB + NMATAB(1) 1742 LWRK1 = LWORK - KEND1 1743C 1744 IF (LWRK1 .LT. 0) THEN 1745 WRITE(LUPRI,*) 'Need:', KEND1, 'Available:', LWORK 1746 CALL QUIT('Insufficient memory for allocation in CCSD_ZKBLO') 1747 ENDIF 1748C 1749 CALL DZERO(WORK,KEND1) 1750 CALL DCOPY(NMATIJ(1),ZKDIA(1),1,WORK(KETAIJ),1) 1751 CALL DCOPY(NMATAB(1),ZKDIA(NMATIJ(1)+1),1,WORK(KETAAB),1) 1752 CALL DZERO(ZKDIA,NMATIJ(1)+NMATAB(1)) 1753C------------------------------------- 1754C Read canonical orbital energies. 1755C------------------------------------- 1756C 1757 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 1758 & .FALSE.) 1759 REWIND (LUSIFC) 1760C 1761 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 1762 READ (LUSIFC) 1763 READ (LUSIFC) (WORK(KFOCKD + I - 1), I = 1,NORBTS) 1764C 1765 CALL GPCLOSE(LUSIFC,'KEEP') 1766C 1767C---------------------------------------------------------------- 1768C Change symmetry ordering of the canonical orbital energies. 1769C---------------------------------------------------------------- 1770C 1771 IF (FROIMP) 1772 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK1) 1773C 1774 CALL FOCK_REORDER(WORK(KFOCKD),WORK(KEND1),LWRK1) 1775! do i=1,NORBTS 1776! write(lupri,*) 'Epsilon_',i,'=',WORK(KFOCKD+i-1) 1777! end do 1778C 1779C--------------------------- 1780C Calculate the results: 1781C Occupied block: 1782C--------------------------- 1783 DO 100 ISYMI = 1,NSYM 1784 ISYMJ = ISYMI 1785 DO 110 J = 1,NRHF(ISYMJ) 1786 KOFFJ = KFOCKD + IRHF(ISYMJ) + J - 1 1787 DO 120 I = J+1,NRHF(ISYMI) 1788 KOFFI = KFOCKD + IRHF(ISYMI) + I - 1 1789 DENOM = WORK(KOFFJ) - WORK(KOFFI) 1790 !write(lupri,*) 'Denom IJ in RCCD_ZKBLO is:', DENOM 1791 IF (ABS(DENOM) .GT. THRDNM) THEN 1792 NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 1793 NJI = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)*(I - 1) + J 1794C 1795 ZKDIA(NIJ) = HALF*WORK(KETAIJ+NIJ-1)/DENOM 1796 ZKDIA(NJI) = ZKDIA(NIJ) 1797 END IF 1798C 1799 120 CONTINUE 1800 110 CONTINUE 1801 100 CONTINUE 1802C 1803C------------------- 1804C Virtual block: 1805C------------------- 1806C 1807 DO 130 ISYMA = 1,NSYM 1808 ISYMB = ISYMA 1809 DO 140 B = 1,NVIR(ISYMB) 1810 KOFFB = KFOCKD + IVIR(ISYMB) + B - 1 1811 DO 150 A = B+1,NVIR(ISYMA) 1812 KOFFA = KFOCKD + IVIR(ISYMA) + A - 1 1813 DENOM = WORK(KOFFB) - WORK(KOFFA) 1814 !write(lupri,*) 'Denom AB in RCCD_ZKBLO is:', DENOM 1815 IF (ABS(DENOM) .GT. THRDNM) THEN 1816 NAB = IMATAB(ISYMA,ISYMB) + NVIR(ISYMA)*(B - 1) + A 1817 NBA = IMATAB(ISYMB,ISYMA) + NVIR(ISYMB)*(A - 1) + B 1818C 1819 ZKDIA(NMATIJ(1) + NAB)= HALF*WORK(KETAAB+NAB-1)/ 1820 * DENOM 1821 ZKDIA(NMATIJ(1) + NBA)= ZKDIA(NMATIJ(1) + NAB) 1822 END IF 1823C 1824 150 CONTINUE 1825 140 CONTINUE 1826 130 CONTINUE 1827C 1828 CALL QEXIT('RCCD_ZKBLO') 1829C 1830 RETURN 1831 END 1832 1833!============ 1834C /* Deck cc_den2_rccd */ 1835 SUBROUTINE CC_DEN2_RCCD(D2IJG,D2AIG,D2IAG,D2ABG, 1836 & Z2PK,Z2AM,T2AM,T2TILDE,T2SQ, 1837 & Z2TILDE, 1838 * XMAT,YMAT,DAB,DAI,DIA, 1839 * XLAMDH,ISYMLH,XLAMDP,ISYMLP, 1840 * WORK,LWORK,IDEL,ISYMD) 1841C 1842C Written by Sonia & Francesca, 18/3/2010, based on CC_DEN2 1843C Debugged by Sonia & Thomas BP, 18/6/2012 1844C Purpose: Directs the calculation of the 2 electron RCCD density 1845C d(pq,gam;del) for a given del (IDEL). The 4 blocks pq 1846C of the result is returned in D2IJG through D2ABG! 1847C Densities must be particle-symmetrized!!!!!!!!!!!!!!!!!!!!!!!!! 1848C 1849#include "implicit.h" 1850 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1851 INTEGER ISYMLP, ISYMLH 1852 DIMENSION D2IJG(*), D2AIG(*), D2IAG(*), D2ABG(*), Z2AM(*) 1853 DIMENSION T2AM(*), T2TILDE(*), Z2TILDE(*), XMAT(*), YMAT(*) 1854 DIMENSION DAB(*), DAI(*), DIA(*), Z2PK(*), XLAMDH(*) 1855 DIMENSION T2SQ(*) 1856 DIMENSION XLAMDP(*), WORK(LWORK) 1857#include "priunit.h" 1858#include "ccorb.h" 1859#include "ccsdsym.h" 1860#include "ccsdinp.h" 1861C 1862 CALL QENTER('CC_DEN2_RCCD') 1863 IF (DEBUG) WRITE(LUPRI,*) 'Entered CC_DEN2_RCCD' 1864C 1865C------------------------------- 1866C set some symmetries: 1867C------------------------------- 1868C 1869 ISYD1 = 1 ! one-particle density 1870 ISYMTR = 1 ! Z2AM 1871 ISYD2H = MULD2H(ISYMD,ISYMLH) ! lambdah backtransformed density 1872 ISYDEN = MULD2H(ISYD2H,ISYMLP) ! lambdah + lambdap transformed 1873 ISYMTI = MULD2H(ISYMLH,ISYMD) 1874C 1875C------------------------------- 1876C Work space allocation one. 1877C------------------------------- 1878C 1879 KT2DEL = 1 1880 KZ2DEL = KT2DEL + NT2BCD(ISYMTI) 1881 KEND1 = KZ2DEL + NT2BCD(ISYMTI) 1882 LWRK1 = LWORK - KEND1 1883C 1884 IF (LWRK1 .LT. 0) THEN 1885 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 1886 CALL QUIT('Insuff. space for 1st allocation, CC_DEN2_RCCD') 1887 ENDIF 1888 1889C 1890C------------------------------------------------------ 1891C Calculate the delta backtransformed amplitude 1892C t_ai,j;delta and multiplier tbar_ai,j;delta 1893C------------------------------------------------------ 1894C 1895 CALL CC_TI(WORK(KT2DEL),ISYMTI,T2AM,ISYMOP,XLAMDH,ISYMLH, 1896 * WORK(KEND1),LWRK1,IDEL,ISYMD) 1897 1898 CALL CC_TI(WORK(KZ2DEL),ISYMTI,Z2PK,ISYMOP,XLAMDH,ISYMLH, 1899 * WORK(KEND1),LWRK1,IDEL,ISYMD) 1900 1901C 1902C--------------------------------------------------- 1903C Calculate the (occ.occ,vir;del) density block. 1904C--------------------------------------------------- 1905C 1906 KD2IJA = KEND1 1907 KEND2 = KD2IJA + NMAIJA(ISYD2H) 1908 LWRK2 = LWORK - KEND2 1909C 1910 IF (LWRK2 .LT. 0) THEN 1911 WRITE (LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1912 CALL QUIT( 1913 * 'Insufficient space for second allocation in CC_DEN2') 1914 ENDIF 1915C 1916 CALL DZERO(WORK(KD2IJA),NMAIJA(ISYD2H)) 1917C 1918C------------------------------------------------ 1919C Contributions from projection against <u2|. 1920C------------------------------------------------ 1921C 1922 !Y part, note that D_ab=2*Y_ba 1923 FACTOR=2.0d0 1924 CALL RPADIJA21(WORK(KD2IJA),ISYD2H,DAB,ISYD1,XLAMDH,ISYMLH, 1925 * WORK(KEND2),LWRK2,IDEL,ISYMD,FACTOR) 1926 IF (RCCD) THEN 1927 FACTOR=2.0d0 1928 CALL RPADIJA22(WORK(KD2IJA),ISYD2H,T2SQ,WORK(KZ2DEL),ISYMTI, 1929 * WORK(KEND2),LWRK2,FACTOR) 1930 END IF 1931C 1932C----------------------------------------------------------------- 1933C Backtransform third virtual index to AO and store in result. 1934C----------------------------------------------------------------- 1935C 1936 ICON = 4 1937 CALL CCD2_PQAO(D2IJG,ISYDEN,WORK(KD2IJA),ISYD2H, 1938 & XLAMDP,ISYMLP,ICON) 1939C 1940C------------------------------------------------------------ 1941C Calculate terms of (occ.vir,occ;del) block with t2-del. 1942C Note that the Z-intermediate is overwritten by the 1943C last calculation! 1944C------------------------------------------------------------ 1945C 1946 ISYMZI = MULD2H(ISYMTI,ISYMTR) 1947C 1948 KD2IAJ = KEND1 1949 KZINT = KD2IAJ + NT2BCD(ISYD2H) 1950 KZBIN = KZINT + NT2BCD(ISYMZI) 1951 KEND2 = KZBIN + NT2BCD(ISYMZI) 1952 LWRK2 = LWORK - KEND2 1953C 1954 IF (LWRK2 .LT. 0) THEN 1955 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 1956 CALL QUIT('Insuff space for 4th allocation in CC_DEN2_RCCD') 1957 ENDIF 1958C 1959 CALL DZERO(WORK(KD2IAJ),NT2BCD(ISYD2H)) 1960C 1961C-------------------------------------------------------------- 1962C Calculate ZINT of t^delta amplitudes. 1963C-------------------------------------------------------------- 1964C 1965 ICON = 1 1966 CALL CC_ZWVI(WORK(KZINT),Z2AM,ISYMTR,WORK(KT2DEL),ISYMTI, 1967 * WORK(KEND2),LWRK2,ICON) 1968C 1969C-------------------------------------------------------------------- 1970C Calculate terms of (occ.vir,occ;del) block with Zbar (in ZINT). 1971C-------------------------------------------------------------------- 1972C 1973 CALL DSCAL(NT2BCD(ISYMZI),4.0d0,WORK(KZINT),1) !to get -4*t*Z 1974 IF (RCCD) THEN 1975 CALL DIAJ26(WORK(KD2IAJ),ISYD2H,T2AM,WORK(KZINT),ISYMZI, 1976 * WORK(KEND2),LWRK2) 1977 END IF 1978 1979 CALL DSCAL(NT2BCD(ISYMZI),2.0d0,WORK(KZINT),1) !to get 8*t*Z 1980 CALL DIAJ28(WORK(KD2IAJ),ISYD2H,T2AM,WORK(KZINT),ISYMZI, 1981 * WORK(KEND2),LWRK2) 1982C 1983C----------------------------------------------------------------- 1984C Add the backtransformed 2C-E of t ampl to (occ.vir,occ;del) 1985C (= Contribution from projection against <HF|) 1986C----------------------------------------------------------------- 1987C 1988 !Compute backtranformed 2C-E of amplitudes and add to density 1989 CALL CC_TI(WORK(KT2DEL),ISYMTI,T2TILDE,ISYMOP,XLAMDH,ISYMLH, 1990 * WORK(KEND2),LWRK2,IDEL,ISYMD) 1991 IF (ISYMTI.NE.ISYD2H) CALL QUIT('Symmetry mismatch in CC_DEN2') 1992 !scale by 2 1993 CALL DAXPY(NT2BCD(ISYMTI),TWO,WORK(KT2DEL),1,WORK(KD2IAJ),1) 1994C 1995C-------------------------------------------------------------- 1996C Calculate ZINT of tbar^delta multipliers (ZBIN) 1997C and place it in of the (vir.occ,occ;del) block. 1998C (first contribution from proj against <u2|?) 1999C-------------------------------------------------------------- 2000C 2001 ICON = 1 2002! recycle later on ZINT 2003! CALL CC_ZWVI(WORK(KZINT),T2SQ,ISYMTR,WORK(KZ2DEL),ISYMTI, 2004! * WORK(KEND2),LWRK2,ICON) 2005 call dzero(WORK(KZBIN),NT2BCD(ISYMZI)) 2006 CALL CC_ZWVI(WORK(KZBIN),T2SQ,ISYMTR,WORK(KZ2DEL),ISYMTI, 2007 * WORK(KEND2),LWRK2,ICON) 2008C 2009 !KD2AIJ = KZINT 2010 KD2AIJ = KZBIN 2011 KEND2 = KD2AIJ + NT2BCD(ISYD2H) 2012 LWRK2 = LWORK - KEND2 2013C 2014 IF (LWRK2 .LT. 0) THEN 2015 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2016 CALL QUIT( 2017 * 'Insufficient space for fourth allocation in CC_DEN2') 2018 ENDIF 2019 !4*Zbar term 2020 CALL DSCAL(NT2BCD(ISYMZI),4.0d0,WORK(KD2AIJ),1) 2021C 2022C------------------------------------------------ 2023C Second contribution from projection against <u2|. 2024C------------------------------------------------ 2025C 2026 FACTOR=1.0d0 2027 CALL RPADAIJ22(WORK(KD2AIJ),ISYD2H,DAB,ISYD1,XLAMDH,ISYMLH, 2028 * WORK(KEND2),LWRK2,IDEL,ISYMD,FACTOR) 2029C 2030C------------------------------------------------------------- 2031C Backtransform occupied index to AO and store in results. 2032C------------------------------------------------------------- 2033C 2034 ICON = 2 2035 CALL CCD2_PQAO(D2IAG,ISYDEN,WORK(KD2IAJ),ISYD2H, 2036 * XLAMDP,ISYMLP,ICON) 2037 2038 2039 CALL CCD2_PQAO(D2AIG,ISYDEN,WORK(KD2AIJ),ISYD2H, 2040 * XLAMDP,ISYMLP,ICON) 2041C 2042C-------------------------------------------------- 2043C Calculate the three blocks: (occ.vir,vir;del) 2044C (vir.occ,vir;del) & (vir.vir,occ;del). 2045C-------------------------------------------------- 2046C 2047 KD2IAB = KEND1 2048 KD2ABI = KD2IAB + NCKATR(ISYD2H) 2049 KD2AIB = KD2ABI + NCKASR(ISYD2H) 2050 KEND2 = KD2AIB + NCKATR(ISYD2H) 2051 LWRK2 = LWORK - KEND2 2052C 2053 IF (LWRK2 .LT. 0) THEN 2054 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2055 CALL QUIT('Insuff space for 5th alloc in CC_DEN2_RCCD') 2056 ENDIF 2057C 2058 CALL DZERO(WORK(KD2IAB),NCKATR(ISYD2H)) 2059 CALL DZERO(WORK(KD2ABI),NCKASR(ISYD2H)) 2060 CALL DZERO(WORK(KD2AIB),NCKATR(ISYD2H)) 2061C 2062C-------------------------------------------------------------- 2063C Backtransf 2C-E of tbar: tbartilde(ai,b;del). 2064C Note that this equals the d(ai,b;del) density block. 2065C-------------------------------------------------------------- 2066C 2067 CALL DAIB21(WORK(KD2AIB),ISYD2H,Z2TILDE,XLAMDH,ISYMLH, 2068 & WORK(KEND2),LWRK2,IDEL,ISYMD) 2069 CALL DSCAL(NCKATR(ISYD2H),1.0d0,WORK(KD2AIB),1) 2070C 2071C------------------------------------------------------- 2072C Backtransform third index of the (vir.occ,vir;del) 2073C block to AO-basis and store in result. 2074C------------------------------------------------------- 2075C 2076 ICON = 5 2077 CALL CCD2_PQAO(D2AIG,ISYDEN,WORK(KD2AIB),ISYD2H, 2078 * XLAMDP,ISYMLP,ICON) 2079C 2080C----------------------------------------------- 2081C Calculate the zeta(ai,b;del)-containing 2082C contributions to the remaining two blocks. 2083C----------------------------------------------- 2084C 2085 KZ2AO = KD2AIB !ricicla lo spazio 2086 ISYZ2AO = ISYD2H 2087 !j-backtr of tam placed on Z2AO 2088 CALL DZERO(WORK(KZ2AO),NCKATR(ISYZ2AO)) 2089 CALL DAIB21(WORK(KZ2AO),ISYD2H,T2SQ,XLAMDH,ISYMLH, 2090 * WORK(KEND2),LWRK2,IDEL,ISYMD) 2091! 2092! THIS IS THE CORRECT ONE FOR IABJ 2093! 2094 CALL DIAC22(WORK(KD2IAB),ISYD2H,Z2PK,WORK(KZ2AO),ISYZ2AO, 2095 & WORK(KEND2),LWRK2) 2096 CALL DSCAL(NCKATR(ISYD2H),4.0d0,WORK(KD2IAB),1) 2097! 2098 IF (RCCD) THEN 2099 CALL DABI22(WORK(KD2ABI),ISYD2H,Z2PK,WORK(KZ2AO),ISYZ2AO, 2100 * WORK(KEND2),LWRK2) 2101 CALL DSCAL(NCKASR(ISYD2H),2.0d0,WORK(KD2ABI),1) 2102 END IF 2103C 2104C-------------------------------------------------------------------- 2105C Calculate remaining contributions from projection against <u2|. 2106C-------------------------------------------------------------------- 2107C 2108 FACTOR=2.0d0 2109 CALL RPADABI21(WORK(KD2ABI),ISYD2H,DAB,ISYD1,XLAMDH,ISYMLH, 2110 & IDEL,ISYMD,FACTOR) 2111 2112 FACTOR=2.0d0 2113 CALL RPADIAC21(WORK(KD2IAB),ISYD2H,YMAT,ISYD1,XLAMDH,ISYMLH, 2114 & IDEL,ISYMD,FACTOR) 2115C 2116C--------------------------------------------------- 2117C Backtransform third index of the two remaining 2118C blocks to AO and store in results. 2119C--------------------------------------------------- 2120 2121 ICON = 5 2122 CALL CCD2_PQAO(D2IAG,ISYDEN,WORK(KD2IAB),ISYD2H, 2123 * XLAMDP,ISYMLP,ICON) 2124 ICON = 3 2125 CALL CCD2_PQAO(D2ABG,ISYDEN,WORK(KD2ABI),ISYD2H, 2126 * XLAMDP,ISYMLP,ICON) 2127C 2128C--------------------------------------------------- 2129C Calculate the (occ.occ,occ;del) density block. 2130C--------------------------------------------------- 2131C 2132 KD2IJK = KEND1 2133 KEND2 = KD2IJK + NMAIJK(ISYD2H) 2134 LWRK2 = LWORK - KEND2 2135C 2136 IF (LWRK2 .LT. 0) THEN 2137 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2 2138 CALL QUIT('Insufficient space for seventh allocation '// 2139 & 'in CC_DEN2') 2140 ENDIF 2141C 2142 CALL DZERO(WORK(KD2IJK),NMAIJK(ISYD2H)) 2143C 2144C------------------------------------------------ 2145C Contributions from projection against <u2|. 2146C Note that X is ''true'' X, not scaled by two!!!! 2147C Density must be defined symmetrized, so IJKL=KLIJ 2148C------------------------------------------------ 2149C 2150! Use CCSD code directly 2151 2152 CALL DIJK21(WORK(KD2IJK),ISYD2H,XMAT,XLAMDH,ISYMLH, 2153 & WORK(KEND2),LWRK2,IDEL,ISYMD) 2154 CALL DIJK23(WORK(KD2IJK),ISYD2H,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD) 2155 CALL DIJK24(WORK(KD2IJK),ISYD2H,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD) 2156 !scale by 2 to match RPA definition 2157 CALL DSCAL(NMAIJK(ISYD2H),2.0d0,WORK(KD2IJK),1) 2158C 2159C------------------------------------------------ 2160C Contributions from projection against <HF|. 2161C------------------------------------------------ 2162C 2163 CALL DIJK01(WORK(KD2IJK),ISYD2H,XLAMDH,ISYMLH,IDEL,ISYMD) 2164 CALL DIJK02(WORK(KD2IJK),ISYD2H,XLAMDH,ISYMLH,IDEL,ISYMD) 2165C 2166C------------------------------------------------------------------ 2167C Backtransform third occupied index to AO and store in result. 2168C------------------------------------------------------------------ 2169C 2170 ICON = 1 2171 CALL CCD2_PQAO(D2IJG,ISYDEN,WORK(KD2IJK),ISYD2H, 2172 * XLAMDP,ISYMLP,ICON) 2173C 2174 CALL QEXIT('CC_DEN2_RCCD') 2175C 2176 RETURN 2177 END 2178 2179 2180C /* Deck rpadija21 */ 2181 SUBROUTINE RPADIJA21(D2IJA,ISYDEN,DAB,ISYDAB,XLAMDH,ISYMLH, 2182 & WORK,LWORK,IDEL,ISYMD,FACTOR) 2183C 2184C MEMO: D_ab is 2Y already!!!! 2185C Modified version of DIJA21 2186C Purpose: To calculate the first contribution to D2IJA 2187C from projection against doubles in RPA/RCCD/DRCCD! 2188C Sonia 2189C 2190#if defined (IMPLICIT_NONE) 2191 IMPLICIT NONE 2192#else 2193# include "implicit.h" 2194#endif 2195#include "priunit.h" 2196#include "ccorb.h" 2197#include "ccsdsym.h" 2198C 2199#if defined (SYS_CRAY) 2200 REAL ZERO, ONE, TWO, FACTOR 2201#else 2202 DOUBLE PRECISION ZERO, ONE, TWO, FACTOR 2203#endif 2204 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2205 2206 INTEGER ISYDEN, ISYDAB, ISYMLH, IDEL, ISYMD, LWORK 2207 INTEGER ISYMB, ISYMA, ISYMIJ, ISYMI, ISYMJ, 2208 & KOFF1, KOFF2, NTOTA, NIJA 2209 2210#if defined (SYS_CRAY) 2211 REAL D2IJA(*), DAB(*), XLAMDH(*), WORK(LWORK) 2212#else 2213 DOUBLE PRECISION D2IJA(*), DAB(*), XLAMDH(*), WORK(LWORK) 2214#endif 2215C 2216 CALL QENTER('RPADIJA21') 2217C 2218 IF (ISYDAB.NE.1) CALL QUIT('Illegal ISYDAI in DIJA21') 2219C 2220 ISYMB = MULD2H(ISYMLH,ISYMD) 2221 ISYMA = MULD2H(ISYDAB,ISYMB) 2222 ISYMIJ = 1 2223C 2224 IF (ISYMA.NE.ISYDEN) THEN 2225 CALL QUIT('Symmetry mismatch in DIJA21. (2)') 2226 END IF 2227C 2228C------------------------------------------------------------------- 2229C Calculate contraction of transposed Y-matrix and lambda matrix. 2230C------------------------------------------------------------------- 2231C 2232 IF (LWORK .LT. NVIR(ISYMA)) THEN 2233 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NVIR(ISYMA) 2234 CALL QUIT('Insufficient space for allocation in DIJA21') 2235 ENDIF 2236C 2237 CALL DZERO(WORK,NVIR(ISYMA)) 2238C 2239 KOFF1 = IMATAB(ISYMA,ISYMB) + 1 2240 KOFF2 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD) 2241C 2242 NTOTA = MAX(NVIR(ISYMA),1) 2243C 2244 CALL DGEMV('N',NVIR(ISYMA),NVIR(ISYMB),FACTOR,DAB(KOFF1),NTOTA, 2245 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 2246C 2247 DO 100 A = 1,NVIR(ISYMA) 2248C 2249 DO 110 ISYMI = 1,NSYM 2250C 2251 ISYMJ = ISYMI 2252C 2253 DO 120 I = 1,NRHF(ISYMI) 2254C 2255 J = I 2256 NIJA = IMAIJA(ISYMIJ,ISYMA) + NMATIJ(ISYMIJ)*(A - 1) 2257 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I 2258C 2259 D2IJA(NIJA) = D2IJA(NIJA) + ONE*WORK(A) 2260C 2261 120 CONTINUE 2262 110 CONTINUE 2263 100 CONTINUE 2264C 2265 CALL QEXIT('RPADIJA21') 2266C 2267 RETURN 2268 END 2269 2270C /* Deck rpadija22 */ 2271 SUBROUTINE rpaDIJA22(D2IJA,ISYDEN,T2SQ,Z2INT,ISYMTI,WORK, 2272 & LWORK,FACTOR) 2273C 2274C Purpose: To calculate the second contribution to D2IJA 2275C from projection against doubles! 2276C in RPA, -2 sum_ck t_aj,ck tbar_ck,i;delta 2277C 2278#include "implicit.h" 2279 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, FOUR = 4.0d0) 2280 DIMENSION D2IJA(*), T2SQ(*), Z2INT(*), WORK(LWORK) 2281#include "priunit.h" 2282#include "ccorb.h" 2283#include "ccsdsym.h" 2284C 2285 CALL QENTER('RPADIJA22') 2286C 2287 IF (ISYDEN.NE.ISYMTI) CALL QUIT('Symmetry mismatch in DIJA22') 2288C 2289 ISYCKI = ISYDEN 2290 ISYIJA = ISYDEN 2291C 2292 DO 100 ISYMI = 1,NSYM 2293C 2294 ISYMCK = MULD2H(ISYMI,ISYCKI) 2295 ISYMAJ = ISYMCK 2296C 2297C------------------------------ 2298C Work space allocation. 2299C------------------------------ 2300C 2301 KSCR = 1 2302 KEND = KSCR + NT1AM(ISYMAJ)*NRHF(ISYMI) 2303 LWRK = LWORK - KEND 2304C 2305 IF (LWRK .LT. 0) THEN 2306 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND 2307 CALL QUIT('Insufficient work space for '// 2308 & 'allocation in DIJA22') 2309 ENDIF 2310C 2311 CALL DZERO(WORK(KSCR),NT1AM(ISYMAJ)*NRHF(ISYMI)) 2312C 2313C-------------------------------------------- 2314C Calculate contraction of zeta and t. 2315C-------------------------------------------- 2316C 2317 KOFF1 = IT2SQ(ISYMAJ,ISYMCK) + 1 2318 KOFF2 = IT2BCD(ISYMCK,ISYMI) + 1 2319C 2320 NTOTAJ = MAX(NT1AM(ISYMAJ),1) 2321 NTOTCK = MAX(NT1AM(ISYMCK),1) 2322C 2323 CALL DGEMM('N','N',NT1AM(ISYMAJ),NRHF(ISYMI),NT1AM(ISYMCK), 2324 * ONE,T2SQ(KOFF1),NTOTAJ,Z2INT(KOFF2),NTOTCK,ZERO, 2325 * WORK(KSCR),NTOTAJ) 2326C 2327C--------------------------------------------------------- 2328C Store properly reordered scratch-array in result. 2329C--------------------------------------------------------- 2330C 2331 DO 110 ISYMA = 1,NSYM 2332C 2333 ISYMJ = MULD2H(ISYMA,ISYMAJ) 2334 ISYMIJ = MULD2H(ISYMJ,ISYMI) 2335C 2336 DO 120 I = 1,NRHF(ISYMI) 2337C 2338 DO 130 J = 1,NRHF(ISYMJ) 2339C 2340 DO 140 A = 1,NVIR(ISYMA) 2341C 2342 NAJI = NT1AM(ISYMAJ)*(I - 1) + IT1AM(ISYMA,ISYMJ) 2343 * + NVIR(ISYMA)*(J - 1) + A 2344 NIJA = IMAIJA(ISYMIJ,ISYMA) 2345 * + NMATIJ(ISYMIJ)*(A - 1) 2346 * + IMATIJ(ISYMI,ISYMJ) 2347 * + NRHF(ISYMI)*(J - 1) + I 2348C 2349 D2IJA(NIJA) = D2IJA(NIJA) 2350 & -FACTOR*WORK(KSCR + NAJI - 1) 2351C 2352 140 CONTINUE 2353 130 CONTINUE 2354 120 CONTINUE 2355 110 CONTINUE 2356 100 CONTINUE 2357C 2358 CALL QEXIT('RPADIJA22') 2359C 2360 RETURN 2361 END 2362 2363C /* Deck rpadijk21 */ 2364 SUBROUTINE rpaDIJK21(D2IJK,ISYIJK,XMAT,XLAMDH,ISYMLH, 2365 * WORK,LWORK,IDEL,ISYMD,FACTOR) 2366C 2367C Purpose: To calculate the first contribution to D2IJK 2368C in RCCD/DRCCD, (projection against doubles)! 2369C 2370#include "implicit.h" 2371 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2372 DIMENSION D2IJK(*), XMAT(*), XLAMDH(*), WORK(LWORK) 2373#include "priunit.h" 2374#include "ccorb.h" 2375#include "ccsdsym.h" 2376C 2377 CALL QENTER('rpaDIJK21') 2378C 2379 ISYML = MULD2H(ISYMD,ISYMLH) 2380 ISYMX1 = ISYML 2381C 2382 IF (ISYML.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK21') 2383C 2384C----------------------------------------------------------------- 2385C Calculate contraction of X-intermediate & lambda hole matrix. 2386C----------------------------------------------------------------- 2387C 2388 IF (LWORK .LT. NRHF(ISYMX1)) THEN 2389 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NRHF(ISYMX1) 2390 CALL QUIT('Insufficient work space in DIJK21') 2391 ENDIF 2392C 2393 CALL DZERO(WORK,NRHF(ISYMX1)) 2394C 2395 KOFF1 = IMATIJ(ISYMX1,ISYML) + 1 2396 KOFF2 = IGLMRH(ISYMD,ISYML) + IDEL - IBAS(ISYMD) 2397C 2398 NTOTI = MAX(NRHF(ISYMX1),1) 2399C 2400 CALL DGEMV('N',NRHF(ISYMX1),NRHF(ISYML),FACTOR,XMAT(KOFF1),NTOTI, 2401 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 2402C 2403C---------------------------------- 2404C Calculate first contribution. 2405C---------------------------------- 2406C 2407 ISYMI = ISYMX1 2408C 2409 DO 100 ISYMK = 1,NSYM 2410C 2411 ISYMJ = ISYMK 2412 ISYMIJ = MULD2H(ISYMJ,ISYMI) 2413C 2414 DO 110 K = 1,NRHF(ISYMK) 2415C 2416 J = K 2417 NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) 2418 * + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + 1 2419C 2420 CALL DAXPY(NRHF(ISYMI),ONE,WORK,1,D2IJK(NIJK),1) 2421C 2422 110 CONTINUE 2423 100 CONTINUE 2424C 2425 CALL QEXIT('rpaDIJK21') 2426C 2427 RETURN 2428 END 2429 2430C /* Deck rpadaij22 */ 2431 SUBROUTINE RPADAIJ22(D2AIJ,ISYAIJ,DAB,ISYD1,XLAMDH,ISYMLH, 2432 * WORK,LWORK,IDEL,ISYMD,FACTOR) 2433C 2434C Written by Asger Halkier 10/2 - 1996 2435C FACTOR added for RPA, Sonia & Francesca, 2010 2436C Version: 1.0 2437C 2438C Purpose: To calculate the second contribution to D2AIJ 2439C from projection against doubles! 2440C 2441#include "implicit.h" 2442 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2443 DIMENSION D2AIJ(*), DAB(*), XLAMDH(*), WORK(LWORK) 2444#include "priunit.h" 2445#include "ccorb.h" 2446#include "ccsdsym.h" 2447C 2448 CALL QENTER('RPADAIJ22') 2449C 2450 ISYMB = MULD2H(ISYMD,ISYMLH) 2451 ISYMA = MULD2H(ISYMB,ISYD1) 2452 ISYMIJ = 1 2453C 2454 IF (ISYMA.NE.ISYAIJ) CALL QUIT('Symmetry mismatch in DAIJ22') 2455C 2456C----------------------------------------------------------------- 2457C Calculate contraction of Y-intermediate & lamda hole matrix. 2458C----------------------------------------------------------------- 2459C 2460 IF (LWORK .LT. NVIR(ISYMA)) THEN 2461 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NVIR(ISYMA) 2462 CALL QUIT('Insufficient work space in DAIJ22') 2463 ENDIF 2464C 2465 CALL DZERO(WORK,NVIR(ISYMA)) 2466C 2467 KOFF1 = IMATAB(ISYMA,ISYMB) + 1 2468 KOFF2 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD) 2469C 2470 NTOTA = MAX(NVIR(ISYMA),1) 2471C 2472 CALL DGEMV('N',NVIR(ISYMA),NVIR(ISYMB),ONE,DAB(KOFF1),NTOTA, 2473 * XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1) 2474C 2475C------------------------------------- 2476C Calculate contribution to D2AIJ. 2477C------------------------------------- 2478C 2479 DO 100 ISYMJ = 1,NSYM 2480C 2481 ISYMI = ISYMJ 2482 ISYMAI = MULD2H(ISYMI,ISYMA) 2483C 2484 DO 110 J = 1,NRHF(ISYMJ) 2485C 2486 I = J 2487 NAIJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) 2488 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2489C 2490 CALL DAXPY(NVIR(ISYMA),-FACTOR,WORK,1,D2AIJ(NAIJ),1) 2491C 2492 110 CONTINUE 2493 100 CONTINUE 2494C 2495 CALL QEXIT('RPADAIJ22') 2496C 2497 RETURN 2498 END 2499!-- 2500C /* Deck rpadabi21 */ 2501 SUBROUTINE rpaDABI21(D2ABI,ISYABI,DAB,ISYD1,XLAMDH,ISYMLH, 2502 & IDEL,ISYMD,FACTOR) 2503C 2504C Written by Asger Halkier 11/2 - 1996 2505C FACTOR added for RPA, Sonia e Francesca, 2010 2506C Version: 1.0 2507C 2508C Purpose: To calculate the first contribution to D2ABI 2509C from projection against doubles! 2510C 2511#include "implicit.h" 2512 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 2513 DIMENSION D2ABI(*), DAB(*), XLAMDH(*) 2514#include "ccorb.h" 2515#include "ccsdsym.h" 2516C 2517 CALL QENTER('RPADABI21') 2518C 2519 ISYMI = MULD2H(ISYMD,ISYMLH) 2520 ISYMAB = ISYD1 2521C 2522 IF (MULD2H(ISYMI,ISYMAB).NE.ISYABI) 2523 * CALL QUIT('Symmetry mismatch in DABI21') 2524C 2525C-------------------------------- 2526C Calculate the contribution. 2527C-------------------------------- 2528C 2529 DO 100 I = 1,NRHF(ISYMI) 2530C 2531 NDEI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) 2532 * + IDEL - IBAS(ISYMD) 2533 NABI = ICKASR(ISYMAB,ISYMI) + NMATAB(ISYMAB)*(I - 1) + 1 2534C 2535 FACT = FACTOR*XLAMDH(NDEI) 2536C 2537 CALL DAXPY(NMATAB(ISYMAB),FACT,DAB,1,D2ABI(NABI),1) 2538C 2539 100 CONTINUE 2540C 2541 CALL QEXIT('RPADABI21') 2542C 2543 RETURN 2544 END 2545!--- 2546C /* Deck rpadiac21 */ 2547 SUBROUTINE rpaDIAC21(D2IAC,ISYIAC,YMAT,ISYMY,XLAMDH,ISYMLH, 2548 * IDEL,ISYMD,FACTOR) 2549C 2550C Written by Asger Halkier 21/2 - 1996 2551C FACTOR introduced for RPA, Sonia & Francesca, 2010 2552C Version: 1.0 2553C 2554C Purpose: To calculate the first contribution to D2IAC 2555C from projection against doubles! 2556C 2557#include "implicit.h" 2558 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0) 2559 DIMENSION D2IAC(*), YMAT(*), XLAMDH(*) 2560#include "priunit.h" 2561#include "ccorb.h" 2562#include "ccsdsym.h" 2563C 2564 CALL QENTER('RPADIAC21') 2565C 2566 ISYMI = MULD2H(ISYMD,ISYMLH) 2567 ISYMAC = ISYMY 2568C 2569 IF (MULD2H(ISYMI,ISYMAC).NE.ISYIAC) 2570 * CALL QUIT('Symmetry mismatch in RPADIAC21') 2571C 2572C-------------------------------- 2573C Calculate the contribution. 2574C-------------------------------- 2575C 2576 DO 100 I = 1,NRHF(ISYMI) 2577C 2578 NDEI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) 2579 * + IDEL - IBAS(ISYMD) 2580C 2581 FACT = -XLAMDH(NDEI)*FACTOR 2582C 2583 DO 110 ISYMA = 1,NSYM 2584C 2585 ISYMC = MULD2H(ISYMA,ISYMAC) 2586 ISYMAI = MULD2H(ISYMA,ISYMI) 2587C 2588 DO 120 C = 1,NVIR(ISYMC) 2589C 2590 NAC = IMATAB(ISYMA,ISYMC) + NVIR(ISYMA)*(C - 1) + 1 2591 NAIC = ICKATR(ISYMAI,ISYMC) + NT1AM(ISYMAI)*(C - 1) 2592 * + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1 2593C 2594 CALL DAXPY(NVIR(ISYMA),FACT,YMAT(NAC),1, 2595 * D2IAC(NAIC),1) 2596C 2597 120 CONTINUE 2598 110 CONTINUE 2599 100 CONTINUE 2600C 2601 CALL QEXIT('RPADIAC21') 2602C 2603 RETURN 2604 END 2605