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! 18*======================================================================* 19 subroutine ccsdr12ao(ccsdr12, 20 & t2am,xlambdap,xlambdah, 21 & fniadj,luiadj,fnijda,luijda, 22 & cpfil,lucp,dpfil,ludp,e1pim, 23 & timintr12,timrdao,timtrbt, 24 & timc,timd,timt2tr,timt2bt, 25 & work,lwork) 26*----------------------------------------------------------------------* 27* Purpose: compute contributions to CCSDR12 requiring the calculation 28* of integrals with auxiliary basis function indices 29* 30* C. Haettig, C. Neiss, spring 2006 31*----------------------------------------------------------------------* 32 implicit none 33#include "priunit.h" 34#include "dummy.h" 35#include "maxorb.h" 36#include "mxcent.h" 37#include "aovec.h" 38#include "iratdef.h" 39#include "ccorb.h" 40#include "ccisao.h" 41#include "blocks.h" 42#include "ccsdsym.h" 43#include "ccsdinp.h" 44#include "cbieri.h" 45#include "distcl.h" 46#include "eritap.h" 47#include "r12int.h" 48#include "ccsdio.h" 49#include "second.h" 50 logical locdbg 51 parameter (locdbg = .false.) 52 53 integer isym0 54 parameter (isym0 = 1) 55 56 real*8 zero,half,one,two,xmhalf,xmone 57 parameter (zero = 0.0d0, half = 0.5d0, one = 1.0d0, two = 2.0d0) 58 parameter (xmhalf = -0.5d0, xmone= -1.0d0 ) 59 60* input: 61 logical ccsdr12 62 integer lwork, lucp, ludp, luiadj, luijda 63 character*(*) cpfil, dpfil, fniadj, fnijda 64 real*8 work(*), xlambdap(*), xlambdah(*), t2am(*), e1pim(*), 65 & timintr12, timrdao, timtrbt, timc, timd, timt2tr, timt2bt 66 67* local 68 logical lauxg, temp_direct 69 integer indexa(MXCORB_CC) 70 integer ioff, ibasx(8), kend1, lwrk1, kendsv, lwrksv, 71 & icdel1, isymd1, ntot, illl, numdis, idel2, idel, isymd, 72 & isydis, irecord, leniaj, isygam, isyalbe, igam, iadr, 73 & icon, iv, isym, nviraop(8), iviraop(8,8), isym1, isym2, 74 & kxint, kxiadj, kxijda, kend3, lwrk3, koffg, kdsrhf, 75 & kt2amt, kend2, lwrk2, kfckvaop, isyfck, lunit, ioptr12 76 integer kccfb1,kindxb,kfree,lfree,ntosym,icount,ibasd, iopte, 77 & kodcl1,kodcl2,kodbc1,kodbc2,krdbc1,krdbc2, 78 & kodpp1,kodpp2,krdpp1,krdpp2,krecnr 79 real*8 dtime, factc, factd, ddot 80 81 call qenter('ccsdr12ao') 82 if (.not.dumpcd) call quit('CCSDR12AO requires DUMPCD=.TRUE.') 83 84 ! switch to integrals with delta index from auxiliary basis: 85 mbsmax = 5 86 loopdp = .true. 87 88 ! integrals with auxiliary basis function are not on file 89 ! -> switch locally to direct mode to calculate them 90 temp_direct = direct 91 direct = .true. 92 93 ! substract aux. functions when calculating index of g index 94 ! within irrep from index running over both basis sets and 95 ! all symmetries 96 ! (needed since isao is overwritten) 97 lauxg = .true. 98 99 ioff = 0 100 ibasx(1) = 0 101 do isym = 1,nsym 102 if (isym.gt.1) ibasx(isym) = ibasx(isym-1)+mbas2(isym-1) 103 if (ioff+mbas1(isym)+mbas2(isym).gt.MXCORB_CC) 104 & call quit('CCSDR12AO') 105 do i = 1,mbas1(isym)+mbas2(isym) 106 ioff = ioff + 1 107 isao(ioff) = isym 108 end do 109 end do 110 111 kend1 = 1 112 lwrk1 = lwork 113 114*----------------------------------------------------------------------* 115* initialization of symmetry arrays and allocation of work space 116* for Fhat(del',a) 117*----------------------------------------------------------------------* 118 do isym = 1, nsym 119 icount = 0 120 do isym2 = 1, nsym 121 isym1 = muld2h(isym2,isym) 122 iviraop(isym1,isym2) = icount 123 icount = icount + nvir(isym1)*mbas2(isym2) 124 end do 125 nviraop(isym) = icount 126 end do 127 128 isyfck = 1 ! symmetry of the Fhat matrix that will be computed 129 130 kfckvaop = kend1 131 kend1 = kfckvaop + nviraop(isyfck) 132 lwrk1 = lwork - kend1 133 if (lwrk1 .lt. 0) call quit('Insufficient core in CCSDR12AO') 134 135 call dzero(work(kfckvaop),nviraop(isyfck)) 136 137*----------------------------------------------------------------------* 138* prepare cluster amplutides with transposed occupied indices: 139*----------------------------------------------------------------------* 140 if ((.not. direct) .and. t2tcor) then 141 kt2amt = kend1 142 kend1 = kt2amt + nt2sq(1) 143 lwrk1 = lwork - kend1 144 if (lwrk1 .lt. 0) call quit('Insufficient core in CCSDR12AO') 145 call dcopy(nt2sq(1),t2am,1,work(kt2amt),1) 146 call ccsd_t2tp(work(kt2amt),work(kend1),lwrk1,isym0) 147 end if 148 149*----------------------------------------------------------------------* 150* intialize integral program 151*----------------------------------------------------------------------* 152 if (direct) then 153 dtime = second() 154 if (herdir) then 155 call herdi1(work(kend1),lwrk1,ipreri) 156 else 157 kccfb1 = kend1 158 kindxb = kccfb1 + mxprim*mxcont 159 kend1 = kindxb + (8*mxshel*mxcont + 1)/irat 160 lwrk1 = lwork - kend1 161 if (lwrk1 .lt.0) 162 & call quit('Insufficient work space in CC_MOFCONR12') 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 timintr12 = timintr12 + ( second() - dtime ) 171 ntosym = 1 172 else 173 ntosym = nsym 174 endif 175 176*----------------------------------------------------------------------* 177* start the loop over integral distributions: 178*----------------------------------------------------------------------* 179 180 kendsv = kend1 181 lwrksv = lwrk1 182 183 icdel1 = 0 184 do isym = 1, nsym 185 icdel1 = icdel1 + nt2bcd(isym)*mbas1(isym) 186 end do 187 188 do isymd1 = 1,ntosym 189 190 if (direct) then 191 if (herdir) then 192 ntot = maxshl 193 else 194 ntot = mxcall 195 endif 196 else 197 ntot = nbas(isymd1) 198 endif 199 200 do illl = 1,ntot 201 202 if (direct) then 203 dtime = second() 204 kend1 = kendsv 205 lwrk1 = lwrksv 206 if (herdir) then 207 call herdi2(work(kend1),lwrk1,indexa,illl,numdis, 208 & ipreri) 209 else 210 call eridi2(illl,indexa,numdis,0,0, 211 & work(kodcl1),work(kodcl2),work(kodbc1), 212 & work(kodbc2),work(krdbc1),work(krdbc2), 213 & work(kodpp1),work(kodpp2),work(krdpp1), 214 & work(krdpp2),work(kccfb1),work(kindxb), 215 & work(kend1), lwrk1,ipreri) 216 endif 217 218 krecnr = kend1 219 kend1 = krecnr + (nbufx(0) - 1)/irat + 1 220 lwrk1 = lwork - kend1 221 if (lwrk1 .lt.0) 222 & call quit('Insufficient work space in CCSDR12AO') 223 timintr12 = timintr12 + ( second() - dtime ) 224 225 if (t2tcor) then 226 kt2amt = kend1 227 kend1 = kt2amt + nt2sq(1) 228 lwrk1 = lwork - kend1 229 if (lwrk1.lt.0) call quit('Insuff. core in CCSDR12AO') 230 call dcopy(nt2sq(1),t2am,1,work(kt2amt),1) 231 call ccsd_t2tp(work(kt2amt),work(kend1),lwrk1,isym0) 232 end if 233 else 234 numdis = 1 235 endif 236 237c ------------------------------------------------------ 238c loop over AOs delta included in the distribution 239c calculated in the above call to the integrals program: 240c ------------------------------------------------------ 241 do idel2 = 1,numdis 242 243 if (direct) then 244 idel = indexa(idel2) 245 isymd = isao(idel) 246 else 247 idel = ibas(isymd1) + ibasx(isymd1) + illl 248 isymd = isymd1 249 endif 250 251 isydis = isymd 252 253C ------------------------------------------------------------ 254C record number and start addresses for intermediates on file: 255C append the records for auxiliary basis functions after those 256C where delta is primary basis function 257C ------------------------------------------------------------ 258 irecord = mbas1t + idel - mbas1(isymd) - ibas(isymd) 259 if (locdbg) 260 & write(lupri,*) 'in CCSDR12AO: irecord = ',irecord 261 it2del(irecord) = icdel1 262 icdel1 = icdel1 + nt2bcd(isydis) 263 264C -------------------------------------------------- 265C allocate memory for 3-index batch of AO integrals 266C and read batch into memory: 267C -------------------------------------------------- 268 kxint = kend1 269 kend2 = kxint + ndisao(isydis) 270 lwrk2 = lwork - kend2 271 if (lwrk2 .lt. 0) 272 & call quit('Insufficient work space in CCSDR12AO') 273 274 dtime = second() 275 call ccrdao(work(kxint),idel,idel2,work(kend2),lwrk2, 276 & work(krecnr),direct) 277 dtime = second() - dtime 278 timrdao = timrdao + dtime 279 280*----------------------------------------------------------------------* 281* compute integrals (ia|delta' j) and (ij|delta' a), where i is obtained 282* by transformation with lambda-particle and j and a by transformation 283* with lambda-hole. the integtrals are saved on disk. 284*----------------------------------------------------------------------* 285 if (ccsdr12 .and. (ianr12.eq.2.or.ianr12.eq.3)) then 286 leniaj = nt2bcd(isydis) 287 288 kxiadj = kend2 289 kxijda = kxiadj + leniaj 290 kend3 = kxijda + leniaj 291 lwrk3 = lwork - kend3 292 if (lwrk3 .lt. 0) 293 & call quit('Insufficient space in CCSDR12AO') 294 295 call dzero(work(kxiadj),leniaj) 296 call dzero(work(kxijda),leniaj) 297CCN 298C write(lupri,*) 'in CCSDR12AO:' 299C write(lupri,*) 'Norm^2 of kxint = ', 300C & ddot(ndisao(isydis),work(kxint),1,work(kxint),1) 301CCN 302 303 do isygam = 1, nsym 304 isyalbe = muld2h(isydis,isygam) 305 do g = 1, mbas1(isygam) 306 igam = g + ibas(isygam) + ibasx(isygam) 307 308 koffg = kxint + idsaog(isygam,isydis) 309 & + nnbst(isyalbe)*(g-1) 310 311 call cc_iajb( work(koffg), isyalbe, dummy, isym0, 312 & idel, igam, lauxg, ibasx, 313 & dummy, work(kxiadj), work(kxijda), 314 & dummy, dummy, dummy, 315 & xlambdap, xlambdah, isym0, 316 & dummy, dummy, isym0, 317 & xlambdap, xlambdah, isym0, 318 & dummy, dummy, isym0, 319 & work(kend3), lwrk3, 3, 320 & .false., .false., .true., 321 & .false., .false., 0 ) 322 end do 323 end do 324 325c ------------------------------------ 326c update Fhat_{del a}: 327c ------------------------------------ 328 ibasd = idel - mbas1(isymd) - ibas(isymd) - ibasx(isymd) 329 if (locdbg) then 330 write(lupri,*) 'in CCSDR12AO: isymd, ibasd = ', 331 & isymd, ibasd 332 end if 333 334 call cc_fckdela(ibasd,isymd,work(kfckvaop),isyfck, 335 & work(kxijda),work(kxiadj),iviraop) 336 337c ------------------------------------ 338c transform (ia|del j) to L(ia|del j): 339c ------------------------------------ 340 call dscal(leniaj, two,work(kxiadj),1) 341 call daxpy(leniaj,-one,work(kxijda),1,work(kxiadj),1) 342 343c -------------------------------------------- 344c write 3-index transformed integrals to disk: 345c -------------------------------------------- 346 iadr = it2del(irecord) + 1 347 call putwa2(luiadj,fniadj,work(kxiadj),iadr,leniaj) 348 call putwa2(luijda,fnijda,work(kxijda),iadr,leniaj) 349 350 if (locdbg) then 351 write(lupri,*) 'in CCSDR12AO: iadr = ', iadr 352 write(lupri,*) 'in CCSDR12AO: Norm^2 of XIADJ: ', 353 & ddot(leniaj,work(kxiadj),1,work(kxiadj),1) 354 write(lupri,*) 'in CCSDR12AO: Norm^2 of XIJDA: ', 355 & ddot(leniaj,work(kxijda),1,work(kxijda),1) 356 end if 357 358 end if 359 360*----------------------------------------------------------------------* 361* Compute the C' and D' intermediates with delta' an auxiliary 362* basis function using CCRHS_C and CCRHS_D. These routines require 363* the amplitudes as full square matrix on t2am 364*----------------------------------------------------------------------* 365 if (ccsdr12 .and. (ianr12.eq.2.or.ianr12.eq.3)) then 366 367C ------------------------------------------------------- 368C transform gamma index in the integral batch to occupied 369C using the lambda-partical = CMO coefficients: 370C ------------------------------------------------------- 371 kdsrhf = kend2 372 kend3 = kdsrhf + ndsrhf(isymd) 373 lwrk3 = lwork - kend3 374 if (lwrk3 .lt. 0) 375 & call quit('Insufficient space in CCSDR12AO') 376 377 dtime = second() 378 call cctrbt(work(kxint),work(kdsrhf),xlambdap, 379 * isym0,work(kend3),lwrk3,isydis) 380 dtime = second() - dtime 381 timtrbt = timtrbt + dtime 382 383 384C ----------------------------- 385C calculate the C intermediate: 386C ----------------------------- 387 dtime = second() 388 factc = xmone 389 icon = 2 390 ioptr12 = 0 !only calculate ONE C-intermediate 391 iopte = 1 392 iv = 1 393 if (.not. t2tcor) then 394 call ccrhs_c(work(kxint),work(kdsrhf),dummy, 395 * t2am,isym0,xlambdap,dummy, 396 * xlambdah,xlambdap,isym0, 397 * xlambdap,isym0, 398 * dummy,e1pim,work(kend3),lwrk3, 399 * irecord,isymd,factc,icon,ioptr12,iopte, 400 * lucp,cpfil,idummy,dummy,iv) 401 else 402 call ccrhs_c(work(kxint),work(kdsrhf),dummy, 403 * work(kt2amt),isym0, 404 * xlambdap,dummy, 405 * xlambdah,xlambdap,isym0, 406 * xlambdap,isym0, 407 * dummy,e1pim,work(kend3),lwrk3, 408 * irecord,isymd,factc,icon,ioptr12,iopte, 409 * lucp,cpfil,idummy,dummy,iv) 410 end if 411 dtime = second() - dtime 412 timc = timc + dtime 413 414c ------------------------- 415c transform T2 to 2T2 - T2. 416c ------------------------- 417 dtime = second() 418 if (t2tcor) then 419 call dscal(nt2sq(1),two,t2am,1) 420 call daxpy(nt2sq(1),-one,work(kt2amt),1,t2am,1) 421 else 422 call ccrhs_t2tr(t2am,work(kend3),lwrk3,isym0) 423 end if 424 dtime = second() - dtime 425 timt2tr = timt2tr + dtime 426 427c ----------------------------- 428c calculate the D intermediate: 429c ----------------------------- 430 dtime = second() 431 factd = one 432 icon = 2 433 ioptr12 = 0 434 iopte = 1 435 iv = 1 436 call ccrhs_d(work(kxint),work(kdsrhf),dummy, 437 * t2am,isym0, 438 * xlambdap,dummy, 439 * xlambdah,xlambdap,isym0, 440 * xlambdah,isym0, 441 * dummy,e1pim,work(kend3),lwrk3, 442 * irecord,isymd,factd,icon,ioptr12,iopte, 443 * ludp,dpfil,idummy,dummy,iv) 444 dtime = second() - dtime 445 timd = timd + dtime 446 447 448c ------------------------- 449c restore T2 from 2T2 - T2: 450c ------------------------- 451 dtime = second() 452 if (t2tcor) then 453 call daxpy(nt2sq(1),one,work(kt2amt),1,t2am,1) 454 call dscal(nt2sq(1),half,t2am,1) 455 else 456 call ccrhs_t2bt(t2am,work(kend3),lwrk3,isym0) 457 end if 458 dtime = second() - dtime 459 timt2bt = timt2bt + dtime 460 461 end if 462 463 end do ! idel2 464 end do ! illl 465 end do ! isymd1 466 467*----------------------------------------------------------------------* 468* end of the loop over integral distributions... 469* restore default settings for the integral evaluation: 470*----------------------------------------------------------------------* 471 472 mbsmax = 4 473 loopdp = .false. 474 direct = temp_direct 475 476 ioff = 0 477 do isym = 1,nsym 478 do i = 1,nbas(isym) 479 ioff = ioff + 1 480 isao(ioff) = isym 481 end do 482 end do 483 484*----------------------------------------------------------------------* 485* save Fhat(a,del') on file: 486*----------------------------------------------------------------------* 487 lunit = -1 488 call gpopen(lunit,'CCFHATADEL','UNKNOWN',' ','UNFORMATTED', 489 & idummy,.false.) 490 read(lunit) 491 write(lunit) (work(kfckvaop-1+i),i=1,nviraop(isyfck)) 492 call gpclose(lunit,'KEEP') 493 494*----------------------------------------------------------------------* 495* return 496*----------------------------------------------------------------------* 497 call qexit('ccsdr12ao') 498 return 499 end 500*----------------------------------------------------------------------* 501* END OF SUBROUTINE CCSDR12AO * 502*======================================================================* 503