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_pt2 */ 20 SUBROUTINE CC_DEN_PT2(ETAAI,ETAIJ,ETAAB,WORK,LWORK, 21 & IOPT,LTSTEN,en2pt) 22C 23C Written by S. Coriani, based on CC_DEN_PT 24C January 2002 25C 26C Version: 1.0 27C 28C Purpose: 29C drive the calculation of the "pure d_pqrs(T)" contributions to the 30C ^kappabar-eta_pq RHS of the orbital multipliers 31C LTSTEN = true, test densities via energy calculation 32C 33#include "implicit.h" 34#include "priunit.h" 35#include "dummy.h" 36#include "maxorb.h" 37#include "maxash.h" 38#include "mxcent.h" 39#include "aovec.h" 40#include "iratdef.h" 41 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 42 PARAMETER (FOUR = 4.0D0) 43 DIMENSION INDEXA(MXCORB_CC) 44 DIMENSION ETAAI(*), ETAIJ(*), ETAAB(*), WORK(LWORK) 45 LOGICAL LTSTEN 46#include "ccorb.h" 47#include "ccisao.h" 48#include "r12int.h" 49#include "inftap.h" 50#include "blocks.h" 51#include "ccfield.h" 52#include "ccsdinp.h" 53#include "ccinftap.h" 54#include "ccsdsym.h" 55#include "ccsdio.h" 56#include "distcl.h" 57#include "cbieri.h" 58#include "eritap.h" 59#include "ccfro.h" 60C 61 CHARACTER MODEL*10 62 CHARACTER NAME1*8 63 CHARACTER NAME2*8 64 65 LOGICAL LOCDBG 66 PARAMETER (LOCDBG=.FALSE.) 67C 68 CALL QENTER('CC_DEN_PT2') 69C 70 CALL HEADER('Construct part of rhs for CCSD(T)-kappa-0',-1) 71C 72C----------------------------------------- 73C Initialization of timing parameters. 74C----------------------------------------- 75C 76 TIMTOT = ZERO 77 TIMTOT = SECOND() 78 TIMDEN = ZERO 79 TIMRES = ZERO 80 TIRDAO = ZERO 81 TIMHE2 = ZERO 82 TIMONE = ZERO 83 TIMONE = SECOND() 84 85 86 IF (LTSTEN) EN2PT=ZERO 87C 88C---------------------------------------------------- 89C Both zeta- and t-vectors are totally symmetric. 90C---------------------------------------------------- 91C 92 ISYMTR = 1 93 ISYMOP = 1 94C 95C---------------------------------------- 96C Get CMO coefficients 97C---------------------------------------- 98C 99 KT1AM_PT = 1 100 KEND0 = KT1AM_PT + NT1AM(ISYMOP) 101 CALL DZERO(WORK(KT1AM_PT),NT1AM(ISYMOP)) 102 103 !KEND0 = 1 104 105 KCMO = KEND0 106 KEND1 = KCMO + NLAMDS 107 LWRK1 = LWORK - KEND1 108 IF (LWRK1 .LT. 0) THEN 109 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1 110 CALL QUIT('Insufficient memory for allocation 1 CC_DEN_PT2') 111 ENDIF 112C 113! IF (FROIMP) THEN 114C 115C------------------------------------------- 116C Get the FULL MO coefficient matrix. 117C------------------------------------------- 118C 119! CALL CMO_ALL(WORK(KCMO),WORK(KEND1),LWRK1) 120C 121! ELSE 122! !Sonia: get coefficients for (T) part and reorder 123! CALL CC_GET_CMO(WORK(KCMO)) 124! CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 125! ENDIF 126C---------------------------------------------------------- 127C Read MO-coefficients from interface file and reorder. 128C---------------------------------------------------------- 129 130 LUSIFC = -1 131 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 132 * IDUMMY,.FALSE.) 133 REWIND LUSIFC 134 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 135 READ (LUSIFC) 136 READ (LUSIFC) 137 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 138 CALL GPCLOSE (LUSIFC,'KEEP') 139C 140 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 141C 142C------------------------------------- 143C Two-electron part starts here..... 144C------------------------------------- 145C 146 TIMONE = SECOND() - TIMONE 147C 148C----------------------------------- 149C Start the loop over integrals. 150C----------------------------------- 151C 152 KENDS2 = KEND1 153 LWRKS2 = LWRK1 154C 155 IF (DIRECT) THEN 156 IF (HERDIR) THEN 157 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 158 ELSE 159 KCCFB1 = KEND1 160 KINDXB = KCCFB1 + MXPRIM*MXCONT 161 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 162 LWRK1 = LWORK - KEND1 163 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 164 * KODPP1,KODPP2,KRDPP1,KRDPP2, 165 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 166 * WORK(KEND1),LWRK1,IPRERI) 167 KEND1 = KFREE 168 LWRK1 = LFREE 169 ENDIF 170 NTOSYM = 1 171 ELSE 172 NTOSYM = NSYM 173 ENDIF 174C 175 KENDSV = KEND1 176 LWRKSV = LWRK1 177C 178 ICDEL1 = 0 179 180 181 xnaigd = zero 182 xnijgd = zero 183 184 DO 100 ISYMD1 = 1,NTOSYM 185C 186 IF (DIRECT) THEN 187 IF (HERDIR) THEN 188 NTOT = MAXSHL 189 ELSE 190 NTOT = MXCALL 191 ENDIF 192 ELSE 193 NTOT = NBAS(ISYMD1) 194 ENDIF 195C 196 DO 110 ILLL = 1,NTOT 197C 198C--------------------------------------------- 199C If direct calculate the integrals. 200C--------------------------------------------- 201C 202 IF (DIRECT) THEN 203C 204 KEND1 = KENDSV 205 LWRK1 = LWRKSV 206C 207 DTIME = SECOND() 208 IF (HERDIR) THEN 209 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 210 & IPRINT) 211 ELSE 212 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 213 * WORK(KODCL1),WORK(KODCL2), 214 * WORK(KODBC1),WORK(KODBC2), 215 * WORK(KRDBC1),WORK(KRDBC2), 216 * WORK(KODPP1),WORK(KODPP2), 217 * WORK(KRDPP1),WORK(KRDPP2), 218 * WORK(KCCFB1),WORK(KINDXB), 219 * WORK(KEND1), LWRK1,IPRERI) 220 ENDIF 221 DTIME = SECOND() - DTIME 222 TIMHE2 = TIMHE2 + DTIME 223C 224 KRECNR = KEND1 225 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 226 LWRK1 = LWORK - KEND1 227 IF (LWRK1 .LT. 0) THEN 228 CALL QUIT('Insufficient core in CC_DEN_PT2') 229 END IF 230C 231 ELSE 232 NUMDIS = 1 233 ENDIF 234C 235C----------------------------------------------------- 236C Loop over number of distributions in disk. 237C----------------------------------------------------- 238C 239 DO 120 IDEL2 = 1,NUMDIS 240C 241 IF (DIRECT) THEN 242 IDEL = INDEXA(IDEL2) 243 IF (NOAUXB) THEN 244 IDUM = 1 245 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 246 END IF 247 ISYMD = ISAO(IDEL) 248 ELSE 249 IDEL = IBAS(ISYMD1) + ILLL 250 ISYMD = ISYMD1 251 ENDIF 252C 253C--------------------------------------------------------- 254C Sonia 255C Work space allocation for the (T) densities 256C with third index backtransformed to gamma 257C All gammas together 258C--------------------------------------------------------- 259C 260 ISYDEN = ISYMD 261C 262 KD2IJG_PT = KEND1 263 KD2AIG_PT = KD2IJG_PT + ND2IJG(ISYDEN) 264 KD2IAG_PT = KD2AIG_PT + ND2AIG(ISYDEN) 265 KD2ABG_PT = KD2IAG_PT + ND2AIG(ISYDEN) 266 KEND2 = KD2ABG_PT + ND2ABG(ISYDEN) 267 LWRK2 = LWORK - KEND2 268C 269 IF (LWRK2 .LT. 0) THEN 270 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 271 CALL QUIT('Insufficient space for allocation '// 272 & '2and1/2 in CC_DEN_PT2') 273 ENDIF 274C 275C------------------------------------------------------- 276C Initialize 4 two electron density arrays. 277C------------------------------------------------------- 278C 279 CALL DZERO(WORK(KD2IJG_PT),ND2IJG(ISYDEN)) 280 CALL DZERO(WORK(KD2AIG_PT),ND2AIG(ISYDEN)) 281 CALL DZERO(WORK(KD2IAG_PT),ND2AIG(ISYDEN)) 282 CALL DZERO(WORK(KD2ABG_PT),ND2ABG(ISYDEN)) 283C 284C------------------------------------------------------------------- 285C Calculate the two electron density d(pq,gamma;delta). 286C------------------------------------------------------------------- 287C 288 AUTIME = SECOND() 289C 290 if (.true.) then 291 CALL CC_DEN2_PT(WORK(KD2IJG_PT),WORK(KD2AIG_PT), 292 & WORK(KD2IAG_PT),WORK(KD2ABG_PT), 293 & WORK(KCMO),1, 294 & WORK(KEND2),LWRK2, 295 & IDEL,ISYMD) 296 end if 297C 298 299 if (.false.) then 300 xtest = ddot(ND2IJG(ISYDEN),WORK(KD2IJG_PT),1, 301 & WORK(KD2IJG_PT),1) 302 write(lupri,*)'norm of ij,gamma;delta', 303 & xtest, isymd, idel 304 xnijgd = xnijgd + xtest 305 xtest = ddot(ND2AIG(ISYDEN),WORK(KD2AIG_PT),1, 306 & WORK(KD2AIG_PT),1) 307 write(lupri,*)'norm of ai,gamma;delta', 308 & xtest, isymd, idel 309 xtest = ddot(ND2AIG(ISYDEN),WORK(KD2IAG_PT),1, 310 & WORK(KD2IAG_PT),1) 311 write(lupri,*)'norm of ia,gamma;delta', 312 & xtest, isymd, idel 313 314 xnaigd = xnaigd + xtest 315 316 xtest = ddot(ND2ABG(ISYDEN),WORK(KD2ABG_PT),1, 317 & WORK(KD2ABG_PT),1) 318 write(lupri,*)'norm of ab;gamma;delta', 319 & xtest, isymd, idel 320 end if 321 322 AUTIME = SECOND() - AUTIME 323 TIMDEN = TIMDEN + AUTIME 324C 325C------------------------------------------ 326C Work space allocation three. 327C------------------------------------------ 328C 329 ISYDIS = MULD2H(ISYMD,ISYMOP) 330C 331 KXINT = KEND2 332 KEND3 = KXINT + NDISAO(ISYDIS) 333 LWRK3 = LWORK - KEND3 334C 335 IF (LWRK3 .LT. 0) THEN 336 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND3 337 CALL QUIT('Insufficient space for allocation '// 338 & '3 in CC_DEN_PT2') 339 ENDIF 340C 341C-------------------------------------------- 342C Read AO integral distribution. 343C-------------------------------------------- 344C 345 AUTIME = SECOND() 346 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND3),LWRK3, 347 * WORK(KRECNR),DIRECT) 348 AUTIME = SECOND() - AUTIME 349 TIRDAO = TIRDAO + AUTIME 350C 351C------------------------------------------------------ 352C Start loop over second AO-index (gamma). 353C------------------------------------------------------ 354C 355 AUTIME = SECOND() 356C 357 DO 130 ISYMG = 1,NSYM 358C 359 ISYMPQ = MULD2H(ISYMG,ISYDEN) 360C 361 DO 140 G = 1,NBAS(ISYMG) 362C 363 KD2GIJ = KD2IJG_PT + ID2IJG(ISYMPQ,ISYMG) 364 * + NMATIJ(ISYMPQ)*(G - 1) 365 KD2GAI = KD2AIG_PT + ID2AIG(ISYMPQ,ISYMG) 366 * + NT1AM(ISYMPQ)*(G - 1) 367 KD2GIA = KD2IAG_PT + ID2AIG(ISYMPQ,ISYMG) 368 * + NT1AM(ISYMPQ)*(G - 1) 369 KD2GAB = KD2ABG_PT + ID2ABG(ISYMPQ,ISYMG) 370 * + NMATAB(ISYMPQ)*(G - 1) 371C 372C----------------------------------------------- 373C Work space allocation four. 374C----------------------------------------------- 375C 376 KINTAO = KEND3 377 KD2AOB = KINTAO + N2BST(ISYMPQ) 378 KEND4 = KD2AOB + N2BST(ISYMPQ) 379 LWRK4 = LWORK - KEND4 380C 381 IF (LWRK4 .LT. 0) THEN 382 WRITE(LUPRI,*) 'Available:', LWORK 383 WRITE(LUPRI,*) 'Needed:', KEND4 384 CALL QUIT('Insufficient space in CC_DEN_PT2') 385 ENDIF 386C 387 CALL DZERO(WORK(KD2AOB),N2BST(ISYMPQ)) 388C 389C------------------------------------------------------- 390C Square up AO-integral distribution. 391C------------------------------------------------------- 392C 393 KOFFIN = KXINT + IDSAOG(ISYMG,ISYDIS) 394 * + NNBST(ISYMPQ)*(G - 1) 395C 396 CALL CCSD_SYMSQ(WORK(KOFFIN),ISYMPQ, 397 * WORK(KINTAO)) 398 399C 400C--------------------------------------------------------- 401C 402C--------------------------------------------------------- 403C 404 if (ltsten) then 405 406 CALL CC_DENAO(WORK(KD2AOB),ISYMPQ, 407 * WORK(KD2GAI),WORK(KD2GAB), 408 * WORK(KD2GIJ),WORK(KD2GIA),ISYMPQ, 409 * WORK(KCMO),1,WORK(KCMO),1, 410 * WORK(KEND4),LWRK4) 411C 412C--------------------------------------------------------------------- 413C Add relaxation terms to set up effective density. 414C--------------------------------------------------------------------- 415C 416! IF (IOPT .EQ. 3) THEN 417C 418! ICON = 1 419! CALL CC_D2EFF(WORK(KD2AOB),G,ISYMG,IDEL, 420! * ISYMD,WORK(KKABAO), 421! * WORK(KDHFAO),ICON) 422C 423! ENDIF 424C 425C---------------------------------------------------------------------- 426C Calculate the 2 e- density contribution to E-ccsd. 427C---------------------------------------------------------------------- 428C 429 EN2PT = EN2PT + HALF*DDOT(N2BST(ISYMPQ), 430 * WORK(KD2AOB),1,WORK(KINTAO),1) 431C 432 end if 433C 434C----------------------------------------------- 435C Work space allocation five. 436C----------------------------------------------- 437C 438 KIJINT = KEND4 439 KAIINT = KIJINT + NMATIJ(ISYMPQ) 440 KIAINT = KAIINT + NT1AM(ISYMPQ) 441 KABINT = KIAINT + NT1AM(ISYMPQ) 442 KEND5 = KABINT + NMATAB(ISYMPQ) 443 LWRK5 = LWORK - KEND5 444C 445 IF (LWRK5 .LT. 0) THEN 446 WRITE(LUPRI,*) 'Available:', LWORK 447 WRITE(LUPRI,*) 'Needed:', KEND5 448 CALL QUIT('Insufficient work space '// 449 & 'in CC_DEN_PT2') 450 ENDIF 451C 452C---------------------------------------------------------------- 453C Transform 2 integral indices to MO-basis. 454C---------------------------------------------------------------- 455C 456 ISYM = ISYMPQ 457 CALL CCDINTMO(WORK(KIJINT),WORK(KIAINT), 458 * WORK(KABINT),WORK(KAIINT), 459 * WORK(KINTAO),WORK(KCMO), 460 * WORK(KCMO),WORK(KEND5), 461 * LWRK5,ISYM) 462C 463C------------------------------------------------------------------- 464C Calculate 2 e- contribution to Zeta-Kappa-0. 465C------------------------------------------------------------------- 466C 467 ISYM = ISYMPQ 468 469 IF (IOPT.EQ.2) THEN 470 471 CALL CCPT_ETARS_2E(ETAIJ,ETAAB, 472 & WORK(KIJINT),WORK(KAIINT), 473 & WORK(KIAINT),WORK(KABINT), 474 & WORK(KD2GIJ),WORK(KD2GAI), 475 & WORK(KD2GIA),WORK(KD2GAB), 476 & WORK(KEND5),LWRK5,ISYM) 477 if (.false.) then 478 XETAIJ = DDOT(NMATIJ(1),ETAIJ(1),1, 479 & ETAIJ(1),1) 480 WRITE(LUPRI,*) 'Norm of eta_ij (2e) after (T)', XETAIJ 481 XETAAB = DDOT(NMATAB(1),ETAAB(1),1, 482 & ETAAB(1),1) 483 WRITE(LUPRI,*) 'Norm of eta_ab (2e) after (T)', XETAAB 484 end if 485 486 END IF 487 488 CALL CCPT_ETAAI_2E(ETAAI, 489 & WORK(KIJINT),WORK(KAIINT), 490 & WORK(KIAINT),WORK(KABINT), 491 & WORK(KD2GIJ),WORK(KD2GAI), 492 & WORK(KD2GIA),WORK(KD2GAB), 493 & WORK(KEND5),LWRK5, 494 & ISYM) 495 if (.false.) then 496 XETAAI = DDOT(NALLAI(1),ETAAI(1),1, 497 & ETAAI(1),1) 498 WRITE(LUPRI,*) 'Norm of eta_ai (2e) after (T)', XETAAI 499 end if 500C 501 502 140 CONTINUE 503 130 CONTINUE 504C 505 AUTIME = SECOND() - AUTIME 506 TIMRES = TIMRES + AUTIME 507C 508 120 CONTINUE 509 110 CONTINUE 510 100 CONTINUE 511 512C 513C----------------------- 514C Regain work space. 515C----------------------- 516C 517 KEND1 = KENDS2 518 LWRK1 = LWRKS2 519C 520C------------------------ 521C 522C------------------------ 523C 524 if (ltsten) then 525 write(lupri,*)'CC_DEN_PT2--> EN2PT: ', EN2PT 526 end if 527C 528C----------------------- 529C Write out timings. 530C----------------------- 531C 532 99 TIMTOT = SECOND() - TIMTOT 533C 534 IF (IPRINT .GT. 3) THEN 535 WRITE (LUPRI,*) ' ' 536 WRITE (LUPRI,*) 'Requested density dependent '// 537 & 'quantities calculated' 538 WRITE (LUPRI,*) 'Total time used in CC_DEN_PT2:', TIMTOT 539 ENDIF 540 IF (IPRINT .GT. 9) THEN 541 WRITE (LUPRI,*) 'Time used for setting up 2 e- density:',TIMDEN 542 WRITE (LUPRI,*) 'Time used for contraction with integrals:', 543 & TIMRES 544 WRITE (LUPRI,*) 'Time used for reading 2 e- AO-integrals:', 545 & TIRDAO 546 WRITE (LUPRI,*) 'Time used for calculating 2 e- AO-integrals:', 547 * TIMHE2 548 WRITE (LUPRI,*) 'Time used for 1 e- density & intermediates:', 549 * TIMONE 550 ENDIF 551C 552 CALL QEXIT('CC_DEN_PT2') 553 RETURN 554 END 555