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 19 SUBROUTINE CCETACOR(ETAAI,ZKDIA,WORK,LWORK) 20C 21C Purpose: to calculate the correction to the RHS (eta_ai and eta_aI) 22C of the z-kappa-bar0_ai equations originating from the 23C "diagonal" z-kappa-bar0_ij and z-kappa-bar0_ab blocks 24C over a loop over delta 25C 26C Written by Sonia Coriani, Fall 2001. Based on CC2_DEN 27C 28C Current models: CCD, CCSD, CCSD(T) 29C 30#include "implicit.h" 31#include "priunit.h" 32#include "dummy.h" 33#include "maxorb.h" 34#include "maxash.h" 35#include "mxcent.h" 36#include "aovec.h" 37#include "iratdef.h" 38 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 39 PARAMETER (FOUR = 4.0D0) 40 DIMENSION INDEXA(MXCORB_CC) 41 DIMENSION ETAAI(*), ZKDIA(*), WORK(LWORK) 42#include "ccorb.h" 43#include "ccisao.h" 44#include "r12int.h" 45#include "inftap.h" 46#include "blocks.h" 47#include "ccfield.h" 48#include "ccsdinp.h" 49#include "ccinftap.h" 50#include "ccsdsym.h" 51#include "ccsdio.h" 52#include "distcl.h" 53#include "cbieri.h" 54#include "eritap.h" 55#include "ccfro.h" 56C 57 CHARACTER MODEL*10 58 CHARACTER NAME1*8 59 CHARACTER NAME2*8 60 61 LOGICAL LOCDBG 62 PARAMETER (LOCDBG=.FALSE.) 63C 64 CALL QENTER('CCETACOR') 65 66C-------------------------------------------------------------- 67C Initialize work space 68C-------------------------------------------------------------- 69 70 KCMO = 1 71 KEND1 = KCMO + NLAMDS 72 LWRK1 = LWORK - KEND1 73C 74C---------------------------------------------------------- 75C Read MO-coefficients from interface file and reorder. 76C---------------------------------------------------------- 77C 78 LUSIFC = -1 79 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 80 * IDUMMY,.FALSE.) 81 REWIND LUSIFC 82 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 83 READ (LUSIFC) 84 READ (LUSIFC) 85 READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS) 86 CALL GPCLOSE (LUSIFC,'KEEP') 87C 88 CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1) 89 90 IF (FROIMP) THEN 91 KCMOF = KEND1 92 KEND1 = KCMOF + NLAMDS 93 LWRK1 = LWORK - KEND1 94 CALL CMO_ALL(WORK(KCMOF),WORK(KEND1),LWRK1) 95 END IF 96C 97C--------------------------------------------------------------- 98C Start the loop for the 2-electron integrals and densities. 99C--------------------------------------------------------------- 100C 101 KENDS2 = KEND1 102 LWRKS2 = LWRK1 103C 104 IF (DIRECT) THEN 105 IF (HERDIR) THEN 106 CALL HERDI1(WORK(KEND1),LWRK1,IPRERI) 107 ELSE 108 KCCFB1 = KEND1 109 KINDXB = KCCFB1 + MXPRIM*MXCONT 110 KEND1 = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT 111 LWRK1 = LWORK - KEND1 112 CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2, 113 * KODPP1,KODPP2,KRDPP1,KRDPP2, 114 * KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB), 115 * WORK(KEND1),LWRK1,IPRERI) 116 KEND1 = KFREE 117 LWRK1 = LFREE 118 ENDIF 119 NTOSYM = 1 120 ELSE 121 NTOSYM = NSYM 122 ENDIF 123C 124 KENDSV = KEND1 125 LWRKSV = LWRK1 126C 127 ICDEL1 = 0 128 DO 100 ISYMD1 = 1,NTOSYM 129C 130 IF (DIRECT) THEN 131 IF (HERDIR) THEN 132 NTOT = MAXSHL 133 ELSE 134 NTOT = MXCALL 135 ENDIF 136 ELSE 137 NTOT = NBAS(ISYMD1) 138 ENDIF 139C 140 DO 110 ILLL = 1,NTOT 141C 142C--------------------------------------------- 143C If direct calculate the integrals. 144C--------------------------------------------- 145C 146 IF (DIRECT) THEN 147C 148 KEND1 = KENDSV 149 LWRK1 = LWRKSV 150C 151c DTIME = SECOND() 152 IF (HERDIR) THEN 153 CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS, 154 & IPRINT) 155 ELSE 156 CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0, 157 * WORK(KODCL1),WORK(KODCL2), 158 * WORK(KODBC1),WORK(KODBC2), 159 * WORK(KRDBC1),WORK(KRDBC2), 160 * WORK(KODPP1),WORK(KODPP2), 161 * WORK(KRDPP1),WORK(KRDPP2), 162 * WORK(KCCFB1),WORK(KINDXB), 163 * WORK(KEND1), LWRK1,IPRERI) 164 ENDIF 165c DTIME = SECOND() - DTIME 166c TIMHE2 = TIMHE2 + DTIME 167C 168 KRECNR = KEND1 169 KEND1 = KRECNR + (NBUFX(0) - 1)/IRAT + 1 170 LWRK1 = LWORK - KEND1 171 IF (LWRK1 .LT. 0) THEN 172 CALL QUIT('Insufficient core in CC2_DEN') 173 END IF 174C 175 ELSE 176 NUMDIS = 1 177 ENDIF 178C 179C----------------------------------------------------- 180C Loop over number of distributions in disk. 181C----------------------------------------------------- 182C 183 DO 120 IDEL2 = 1,NUMDIS 184C 185 IF (DIRECT) THEN 186 IDEL = INDEXA(IDEL2) 187 IF (NOAUXB) THEN 188 IDUM = 1 189 CALL IJKAUX(IDEL,IDUM,IDUM,IDUM) 190 END IF 191 ISYMD = ISAO(IDEL) 192 ELSE 193 IDEL = IBAS(ISYMD1) + ILLL 194 ISYMD = ISYMD1 195 ENDIF 196C 197C---------------------------------------- 198C Work space allocation two. 199C---------------------------------------- 200C 201 ISYDEN = ISYMD 202 ISYDIS = MULD2H(ISYMD,ISYMOP) 203C 204 KXINT = KEND1 205 KEND2 = KXINT + NDISAO(ISYDIS) 206 LWRK2 = LWORK - KEND2 207C 208 IF (LWRK2 .LT. 0) THEN 209 WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2 210 CALL QUIT( 211 * 'Insufficient space for allocation 2 in CCETAAICORR') 212 ENDIF 213C 214C-------------------------------------------- 215C Read AO integral distribution. 216C-------------------------------------------- 217C 218 AUTIME = SECOND() 219 CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2, 220 * WORK(KRECNR),DIRECT) 221C 222C-------------------------------------------------------------- 223C Calculate correction terms to eta_ai originating 224C from kbar_ij and kbar_ab. 225C-------------------------------------------------------------- 226C 227 KDSRHF = KEND2 228 K3OINT = KDSRHF + NDSRHF(ISYMD) 229 KEND3 = K3OINT + NMAIJK(ISYDIS) 230 IF (FROIMP) THEN 231 KDSFRO = KEND3 232 KOFOIN = KDSFRO + NDSFRO(ISYDIS) 233 KEND3 = KOFOIN + NOFROO(ISYDIS) 234 ENDIF 235 LWRK3 = LWORK - KEND3 236C 237 IF (LWRK3 .LT. 0) THEN 238 WRITE(LUPRI,*) 'Need : ',KEND3,'Available : ',LWORK 239 CALL QUIT( 240 * 'Insufficient core for integrals in CCETAC') 241 ENDIF 242C 243C-------------------------------------------------------------------- 244C Transform one index in the integral batch to occupied. 245C-------------------------------------------------------------------- 246C 247 CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KCMO), 248 * ISYMOP,WORK(KEND3),LWRK3,ISYDIS) 249csonia 250C 251C------------------------------------------------------------------ 252C Transform one index in the integral batch to frozen. 253C------------------------------------------------------------------ 254C 255 IF (FROIMP) THEN 256C 257 CALL CC_GTOFRO(WORK(KXINT),WORK(KDSFRO),WORK(KCMOF), 258 * WORK(KEND3),LWRK3,ISYDIS) 259C 260C-------------------------------------------------------------- 261C Calculate integral batch (cor fro | cor del). 262C-------------------------------------------------------------- 263C 264 CALL CC_OFROIN(WORK(KDSRHF),WORK(KOFOIN),WORK(KCMOF), 265 * WORK(KEND3),LWRK3,ISYDIS) 266C 267C------------------------------------------------------------------------- 268C Calculate correction to eta_aI from kappabar_ab 269C------------------------------------------------------------------------- 270C 271 if (.true.) then 272 KAFROI = 1 + NMATIJ(1) + NMATAB(1) + 2*NT1AMX 273 CALL MP2_EIDV1(ZKDIA(KAFROI),WORK(KDSFRO), 274 * ZKDIA(NMATIJ(1)+1),WORK(KCMOF), 275 * WORK(KEND3),LWRK3,IDEL,ISYMD) 276C 277 CALL MP2_EIDV2(ZKDIA(KAFROI),WORK(KDSFRO), 278 * ZKDIA(NMATIJ(1)+1),WORK(KCMOF), 279 * WORK(KEND3),LWRK3,IDEL,ISYMD) 280C 281C------------------------------------------------------------------------- 282C Calculate correction to eta_aI from kappabar_ij 283C------------------------------------------------------------------------- 284C 285 CALL MP2_EIDC1(ZKDIA(KAFROI),WORK(KDSFRO), 286 * ZKDIA(1),WORK(KCMOF),WORK(KEND3), 287 * LWRK3,IDEL,ISYMD) 288C 289 CALL MP2_EIDC2(ZKDIA(KAFROI),WORK(KOFOIN), 290 * ZKDIA(1),WORK(KCMOF),WORK(KEND3), 291 * LWRK3,IDEL,ISYMD) 292 293 end if 294 END IF 295csonia 296C 297C------------------------------------------------------------------- 298C Calculate integral batch with three occupied indices. 299C Since LUDUMLOCAL < 0 the integrals are not written 300C to disk. 301C------------------------------------------------------------------- 302C 303 LUDUMLOCAL = -10 304 CALL CC_INT3O(WORK(K3OINT),WORK(KDSRHF),WORK(KCMO), 305 * ISYMOP,WORK(KCMO),WORK(KEND3),LWRK3, 306 * IDEL,ISYMD,LUDUMLOCAL,'DUMMY') 307C 308C------------------------------------------------------- 309C Calculate the correction terms to eta_ai 310C from kappabar_ab and kappabar_ij 311C------------------------------------------------------- 312C 313 CALL MP2_YTV(ETAAI,ZKDIA(NMATIJ(1)+1),WORK(KDSRHF), 314 * WORK(KCMO),WORK(KEND3),LWRK3,IDEL,ISYMD) 315C 316 CALL MP2_NXY(ETAAI,ZKDIA(1),ZKDIA(NMATIJ(1)+1), 317 * WORK(K3OINT),WORK(KDSRHF),WORK(KCMO), 318 * WORK(KEND3),LWRK3,IDEL,ISYMD) 319C 320 CALL MP2_XTO(ETAAI,ZKDIA(1),WORK(K3OINT), 321 * WORK(KCMO),WORK(KEND3),LWRK3,IDEL,ISYMD) 322C 323 AUTIME = SECOND() - AUTIME 324C 325 120 CONTINUE 326 110 CONTINUE 327 100 CONTINUE 328 329C----------------------------------------------------- 330C Some print out 331C----------------------------------------------------- 332 IF (LOCDBG) THEN 333 ETAKAN = DDOT(NALLAI(1),ETAAI,1,ETAAI,1) 334 WRITE(LUPRI,*) 'CCETACOR: eta_ai ' 335 WRITE(LUPRI,*) 'Norm of corrected occ-vir block:', ETAKAN 336 ENDIF 337C 338C----------------------- 339C Regain work space. 340C----------------------- 341C 342 KEND1 = KENDS2 343 LWRK1 = LWRKS2 344C 345C--------------------------------------------------------------- 346 CALL QEXIT('CCETACOR') 347 RETURN 348 END 349