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 r12batf(r2am,fback,fback2,lauxbeta,projvir,work,lwork) 20c---------------------------------------------------------------------- 21c Purpose: back transforms two indices of the integrals 22c (kp|r_12|lq) integrals to the contravariant AO basis 23c and store them on file 24c 25c (k gamma|r_12|l delta) = sum_pq C_gamma,p C_delta,q (kp|r_12|lq) 26c 27c where: k,l active occupied / r12 MOs 28c p,q general MOs (orbital + auxiliary basis) 29c gamma AO in orbital basis 30c delta AO in orbital or auxiliary basis 31c 32c r2am contains (kp|r_12|lq) stored as symmetry 33c packed triangular matrix using it2am,nt1am,it1am 34c fback file name for back transformed integrals 35c fback2 file name for back transformed integrals with 36c two indices from the auxiliary basis 37c lauxbeta flag whether to calculate also integrals with 38c p transformed to the auxiliary basis 39c projvir include only virtual orbitals in transformation 40c 41c Note: this routine is used in the R12 environment with a 42c number of dimension defined differently form the 43c definitions used in the Coupled Cluster program: 44c 45c mbas1 AOs in orbital basis 46c mbas2 AOs in auxiliary basis 47c norb1 orbital basis functions per irrep 48c norb2 auxiliary basis functions per irrep 49c nvir general MOs per irrep ( = norb1 + norb2 ) 50c nrhf active occupied / r12 MOs 51c nrhfa active occupied MOs 52c nvircc active virtual MOs 53c nbas mbas1 + mbas2 54c 55c Heike Fliegl, Christof Haettig spring 2003 56c adapted for max. two auxiliary functions: Christian Neiss, fall 2004 57c introduced projvir option: Christof Haettig, spring 2005 58c adapted for CABS, Ansatz 3: Christian Neiss, winter 2006 59c---------------------------------------------------------------------- 60 implicit none 61#include "priunit.h" 62#include "ccorb.h" 63#include "ccsdsym.h" 64#include "r12int.h" 65 66 integer isymr2 67 parameter (isymr2 = 1) 68 69 logical locdbg 70 parameter (locdbg = .false.) 71 72 character*(*) fback,fback2 73 logical lauxbeta, projvir, modifyinttyp1 74 75#if defined (SYS_CRAY) 76 real zero,one 77#else 78 double precision zero,one 79#endif 80 parameter (zero = 0.0D0, one = 1.0D0) 81 82 integer nr1orb(8),nrgkl(8),irgkl(8,8) 83 integer ir1orb(8,8),ir1bas(8,8),nr1bas(8),ir1xbas(8,8),nr2xbas(8) 84 integer nr1xbas(8),ir2bas(8,8),ir2xbas(8,8),n2bst1(8),iaodis1(8,8) 85 integer nr1xorb(8),ir1xorb(8,8),nrxgkl(8),irxgkl(8,8) 86 integer nr2bas2,ir2bas2(8,8),ir2xbas2(8,8),nvircc(8) 87 integer isymql,lur2back,isymp,isympk,isymk,isyml 88 integer lwork,isymq,isymd,idelta,idel,isym,isym1,isym2 89 integer icount3,icount4 90 integer kcmo,lwrk0,lusifc,idummy 91 integer nr2bas,kcmobas,kcmoaux,lwrk1,koff1,ntotg1 92 integer kr2sm,kend2,kr2gk,kr2pk,lwrk2,istartq,iendq,ntotdk,isymg 93 integer koffl,kofft,ntotp,ntotg,ndell,iadr,len 94 integer kend0,kend1,kend3,lwrk3 95 integer iendtyp,inttyp,istartdelta,ienddelta,istartgamma 96 integer istartp,iendp,norbp,norbq,nbasg,lur2back2 97 integer nalphaj(8),ialphaj(8,8) 98celena 99 integer isymdk,kr2sm2,kr2gk2,kr2pk2,icount1 100 integer koffs 101celena 102#if defined (SYS_CRAY) 103 real r2am(*), work(*), dnrm2, ddot 104#else 105 double precision r2am(*), work(*), dnrm2, ddot 106#endif 107 108c----------------------------------------------------------------------- 109c some initializations: 110c set label for trace routines 111c calculate a number of symmetry offsets and dimensions 112c----------------------------------------------------------------------- 113 call qenter('r12batf') 114 if(locdbg) write(lupri,*)'entered r12batf' 115 116 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 117 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas,ir2bas,ir2xbas, 118 & irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 119 120 ! nr2bas2 = lenghts for nr1bas x nr1xbas over all symmetries 121 ! nvircc = number of active virtual orbitals in a symmetry 122 nr2bas2 = 0 123 do isym = 1, nsym 124 nr2bas2 = nr2bas2 + nr1bas(isym)*nr1xbas(isym) 125 nvircc(isym) = norb1(isym) - nrhfsa(isym) - nvirfr(isym) 126 end do 127 ! ir2bas2 = symmetry offsets for nr1xbas x nr1bas 128 ! ir2xbas2 = symmetry offsets for nr1xbas x nr1xbas 129 do isym = 1, nsym 130 icount3 = 0 131 icount4 = 0 132 do isym2 = 1, nsym 133 isym1 = muld2h(isym,isym2) 134 ir2bas2(isym1,isym2) = icount3 135 ir2xbas2(isym1,isym2) = icount4 136 icount3 = icount3 + nr1xbas(isym1)*nr1bas(isym2) 137 icount4 = icount4 + nr1xbas(isym1)*nr1xbas(isym2) 138 end do 139 end do 140 141celena 142c do isymdk = 1, nsym 143c nr2orb(isymdk) = 0 144c icount1 = 0 145c do isymk = 1, nsym 146c isymd = muld2h(isymdk,isymk) 147c ir2orb(isymd,isymk) = icount1 148c icount1 = icount1 + nrhf(isymk)*norb2(isymd) 149c nr2orb(isymdk) = nr2orb(isymdk) + norb2(isymd)*nrhf(isymk) 150c end do 151c end do 152celena 153c nr2orb:=nr1xorb 154c ir2orb:=ir1xorb 155 156 157 if (locdbg) then 158 do i = 1, nsym 159 write(lupri,*) 'nrhf(',i,') = ',nrhf(i) 160 write(lupri,*) 'norb1(',i,') = ',norb1(i) 161 write(lupri,*) 'norb2(',i,') = ',norb2(i) 162 write(lupri,*) 'nvir(',i,') = ',nvir(i) 163 write(lupri,*) 'mbas1(',i,') = ',mbas1(i) 164 write(lupri,*) 'mbas2(',i,') = ',mbas2(i) 165 write(lupri,*) 'nbas(',i,') = ',nbas(i) 166 end do 167 end if 168 169c----------------------------------------------------------------------- 170c get MO coefficients (defined for combined orbital + aux. basis): 171c----------------------------------------------------------------------- 172 kcmo = 1 173 kend0 = kcmo + nlamds 174 lwrk0 = lwork - kend0 175 if (lwrk0 .lt.0) then 176 call quit('Insufficient work space in R12BATF') 177 end if 178 179 lusifc = -1 180 call gpopen(lusifc,'SIRIFC','OLD',' ','UNFORMATTED', 181 & idummy,.false.) 182 rewind(lusifc) 183 call mollab('FULLBAS ',lusifc,lupri) 184 read(lusifc) 185 read(lusifc) 186 read(lusifc) (work(kcmo+i-1),i=1,nlamds) 187 call gpclose(lusifc,'KEEP') 188 !Reorder C_MO coefficients such that first all occ. orbitals 189 !in all symmetries, then all virtual(here in MP2-R12: = ALL orbitals!!) 190 !orbitals in all symmetries; by this, the redundant occ. 191 !orbitals are grouped in one single block! 192 call cmo_reorder(work(kcmo),work(kend0),lwrk0) 193 194C----------------------------------------------------------------------- 195C Loop over types of integrals inttyp 196C----------------------------------------------------------------------- 197 if (projvir) then 198 iendtyp = 1 ! integrals of type (kp|r12|lq) 199 else if (.not. lauxbeta) then 200 iendtyp = 2 ! integrals of type (kp|r12|lq) and (kp|r12|lq') 201 if (mbsmax.lt.5) then 202 write (lupri,*) 'MBSMAX = ',MBSMAX 203 call quit("MBSMAX must at least be 5 for (kp|r12|lq')") 204 end if 205 else if (lauxbeta) then 206 iendtyp = 4 ! additionally (kp'|r12|lq) and (kp'|r12|lq') 207 if (mbsmax.lt.6) then 208 write (lupri,*) 'MBSMAX = ',MBSMAX 209 call quit("MBSMAX must at least be 6 for (kp'|r12|lq')") 210 end if 211 end if 212 213 modifyinttyp1 = ianr12.eq.1 .and. r12cbs 214 215 if (locdbg) then 216 write (lupri,*) 'MBSMAX = ', MBSMAX 217 end if 218 219 do inttyp = 1, iendtyp 220 221c----------------------------------------------------------------------- 222c open file for back transformed r12 integrals 223c----------------------------------------------------------------------- 224 if ((inttyp .eq. 1) .or. (inttyp .eq. 2)) then 225 lur2back = -1 226 call wopen2(lur2back,fback,64,0) 227 else if ((inttyp .eq. 3) .or. (inttyp .eq. 4)) then 228 lur2back2 = -1 229 call wopen2(lur2back2,fback2,64,0) 230 end if 231 232C----------------------------------------------------------------------- 233C Loop over symmetries of delta 234C----------------------------------------------------------------------- 235 do isymd = 1, nsym 236 isymq = isymd 237 238C----------------------------------------------------------------------- 239C Set start and end values for orbital spaces (part1) 240C----------------------------------------------------------------------- 241 if (modifyinttyp1 .and. (inttyp .eq. 1)) then 242 istartdelta = 1 243 ienddelta = mbas1(isymd) 244 istartq = 1 245 iendq = norb1(isymq) + norb2(isymq) 246 norbq = norb1(isymq) + norb2(isymq) 247 elseif ((inttyp .eq. 1) .or. (inttyp .eq. 3)) then 248 istartdelta = 1 249 ienddelta = mbas1(isymd) 250 if (projvir) then 251 istartq = nrhfsa(isymq) + 1 252 iendq = nrhfsa(isymq) + nvircc(isymq) 253 norbq = nvircc(isymq) 254 else 255 istartq = 1 256 iendq = norb1(isymq) 257 norbq = norb1(isymq) 258 end if 259 else if ((inttyp .eq. 2) .or. (inttyp .eq. 4)) then 260 istartdelta = mbas1(isymd) + 1 261 ienddelta = mbas1(isymd) + mbas2(isymd) 262 istartq = norb1(isymq) + 1 263 iendq = norb1(isymq) + norb2(isymq) 264 norbq = norb2(isymq) 265 end if 266 267C----------------------------------------------------------------------- 268C Start loop over delta 269C----------------------------------------------------------------------- 270 do idelta = istartdelta, ienddelta 271 idel = ibas(isymd) + idelta 272 273 if (locdbg) write(lupri,*) 'isymd, idelta, idel: ', 274 & isymd, idelta, idel 275 276c ------------------------------------------------------------ 277c get row vector out of transformation matrix C_delta,q 278c for fixed delta: 279c ------------------------------------------------------------ 280 kcmobas = kend0 281 if (modifyinttyp1 .and. inttyp .eq. 1) then 282 kend1 = kcmobas + norb1(isymd) + norb2(isymd) 283 koff1 = kcmo - 1 + iglmvi(isymd,isymq) + idelta 284 elseif ((inttyp.eq.1) .or. (inttyp.eq.3)) then 285 if (projvir) then 286 kend1 = kcmobas + norbq 287 koff1 = kcmo - 1 + iglmvi(isymd,isymq) + 288 & nbas(isymd)*(istartq-1) + idelta 289 else 290 kend1 = kcmobas + norb1(isymd) 291 koff1 = kcmo - 1 + iglmvi(isymd,isymq) + idelta 292 end if 293 else 294 kend1 = kcmobas + norb2(isymd) 295 koff1 = kcmo - 1 + iglmvi(isymd,isymq) + 296 & nbas(isymd)*norb1(isymq) + idelta 297 end if 298 lwrk1 = lwork - kend1 299 if (lwrk1 .lt. 0) then 300 call quit('Insufficient work space in R12BATF') 301 end if 302 call dcopy(norbq,work(koff1),nbas(isymd), 303 & work(kcmobas),1) 304 if (locdbg) then 305 write(lupri,*) 'C_MO vector for fixed delta:' 306 write(lupri,*) 'idelta, C_MO: ',idelta 307 call output(work(kcmobas),1,1,1,kend1-kcmobas, 308 & 1,kend1-kcmobas,1,lupri) 309 end if 310 311c ------------------------------------------------------------ 312c loop over symmetry blocks of the pair index (pk) -> 313c determines also symmetry of pair (ql) and of index l 314c ------------------------------------------------------------ 315 do isympk = 1, nsym 316 isymql = muld2h(isympk,isymr2) 317 isyml = muld2h(isymql,isymq) 318 319 if (nrhf(isyml) .ne. 0) then 320 321C -------------------------------------------------------- 322C allocate memory for: 323C kr2gk : R(gamma,k) matrix AO * occupied 324C -------------------------------------------------------- 325 kr2gk = kend1 326 if ((inttyp .eq. 1) .or. (inttyp .eq. 2)) then 327 if (modifyinttyp1 .and. (inttyp .eq. 1)) then 328 kr2gk2= kr2gk + nr1bas(isympk) 329 kend2 = kr2gk2+ nr1bas(isympk) 330 else 331 kend2 = kr2gk + nr1bas(isympk) 332 endif 333 else if ((inttyp .eq. 3) .or. (inttyp .eq. 4)) then 334 kend2 = kr2gk + nr1xbas(isympk) 335 end if 336 lwrk2 = lwork - kend2 337 if (lwrk2 .lt. 0) then 338 call quit('Insufficient work in R12BATF') 339 end if 340 341C -------------------------------------------------------- 342C Loop over occupied orbitals l 343C -------------------------------------------------------- 344 do l = 1, nrhf(isyml) 345 346C ------------------------------------------------------ 347C Loop over symmetries of gamma 348C ------------------------------------------------------ 349 do isymg = 1, nsym 350 isymp = isymg 351 isymk = muld2h(isympk,isymp) 352 353C ---------------------------------------------------- 354C Set start and end values for orbital spaces (part 2) 355C ---------------------------------------------------- 356 if ((inttyp .eq. 1) .or. (inttyp .eq. 2)) then 357 istartgamma = 1 358 nbasg = mbas1(isymg) 359 if (projvir) then 360 istartp = nrhfsa(isymp) + 1 361 iendp = nrhfsa(isymp) + nvircc(isymp) 362 norbp = nvircc(isymp) 363 else 364 istartp = 1 365 iendp = norb1(isymp) 366 norbp = norb1(isymp) 367 end if 368 else if ((inttyp .eq. 3) .or. (inttyp .eq. 4)) then 369 istartgamma = mbas1(isymg) + 1 370 nbasg = mbas2(isymg) 371 istartp = norb1(isymp) + 1 372 iendp = norb1(isymp) + norb2(isymp) 373 norbp = norb2(isymp) 374 end if 375 376c ---------------------------------------------------- 377c allocate work space for: 378c kr2sm : r(pk,q) three-index array batch of r2am 379c kr2pk : R(pk) matrix MO * occupied 380c ---------------------------------------------------- 381 kr2sm = kend2 382 if (modifyinttyp1 .and. inttyp .eq. 1) then 383 kr2sm2= kr2sm + norb1(isymp)*nrhf(isymk)* 384 & (norb1(isymq)+norb2(isymq)) 385 kr2pk = kr2sm2+ norb2(isymp)*nrhf(isymk)* 386 & norb1(isymq) 387 kr2pk2= kr2pk + norb1(isymp)*nrhf(isymk) 388 kend3 = kr2pk2+ norb2(isymp)*nrhf(isymk) 389 else 390 kr2pk = kr2sm + norbp * nrhf(isymk) * norbq 391 kend3 = kr2pk + norbp * nrhf(isymk) 392 endif 393 lwrk3 = lwork - kend3 394 if (lwrk3 .lt. 0) then 395 call quit('Insufficient work space in R12BATF') 396 end if 397 398 call cc_r12getrpkq(fback,work(kr2sm),l,isyml, 399 & istartq,iendq,isymq,istartp, 400 & iendp,isymp,r2am, 401 & isympk,isymql,nr1orb,ir1orb) 402 if (modifyinttyp1 .and. (inttyp .eq.1)) then 403 call cc_r12getrpkq(fback,work(kr2sm2),l,isyml, 404 & 1,norb1(isymq),isymq, 405 & norb1(isymp)+1, 406 & norb1(isymp) +norb2(isymp), 407 & isymp,r2am, 408 & isympk,isymql,nr1orb,ir1orb) 409 endif 410 if (locdbg) then 411 write(lupri,*) 'r^(l)_(pkq):' 412 call output(work(kr2sm),1,norbp*nrhf(isymk),1, 413 & norbq,norbp*nrhf(isymk),norbq, 414 & 1,lupri) 415 if (modifyinttyp1 .and. (inttyp .eq.1)) then 416 write(lupri,*) 'r^(l)_(pkq) 2:' 417 call output(work(kr2sm2),1,norb2(isymp)* 418 & nrhf(isymk),1,norb1(isymq), 419 & norb2(isymp)*nrhf(isymk), 420 & norb1(isymq),1,lupri) 421 end if 422 end if 423 424c ---------------------------------------------------- 425c transform q to delta (contravariant AO basis) 426c ---------------------------------------------------- 427 ntotdk = max(norbp*nrhf(isymk),1) 428 !zero out result (dgemv may give wrong result)! 429 call dzero(work(kr2pk),norbp*nrhf(isymk)) 430 call dgemv('N',norbp*nrhf(isymk), 431 & norbq,one, 432 & work(kr2sm),ntotdk, 433 & work(kcmobas),1, 434 & zero,work(kr2pk),1) 435 if (modifyinttyp1 .and. (inttyp .eq.1)) then 436 ntotdk = max(norb2(isymp)*nrhf(isymk),1) 437 !zero out result (dgemv may give wrong result)! 438 call dzero(work(kr2pk2),norb2(isymp)*nrhf(isymk)) 439 call dgemv('N',norb2(isymp)*nrhf(isymk), 440 & norb1(isymq),one, 441 & work(kr2sm2),ntotdk, 442 & work(kcmobas),1, 443 & zero,work(kr2pk2),1) 444 endif 445 446c ---------------------------------------------------- 447c transform p to gamma (contravariant AO basis) 448c ---------------------------------------------------- 449 koffl = kcmo + iglmvi(isymg,isymp) 450 & + nbas(isymg)*(istartp-1) + istartgamma-1 451 if ((inttyp .eq. 1) .or. (inttyp .eq. 2)) then 452 kofft = kr2gk + ir1bas(isymg,isymk) 453 else if ((inttyp .eq. 3) .or. (inttyp .eq. 4)) then 454 kofft = kr2gk + ir1xbas(isymg,isymk) 455 end if 456 457 ntotp = max(norbp,1) 458 ntotg = max(nbas(isymg),1) 459 ntotg1 = max(nbasg,1) 460 call dgemm('N','N',nbasg,nrhf(isymk),norbp, 461 & one,work(koffl),ntotg, 462 & work(kr2pk),ntotp, 463 & zero,work(kofft),ntotg1) 464 if (modifyinttyp1 .and. (inttyp .eq. 1)) then 465 koffl=kcmo + norb1(isymp)*nbas(isymg) 466 & + iglmvi(isymg,isymp) 467 koffs=kr2gk2 + ir1bas(isymg,isymk) 468 ntotp=max(norb2(isymp),1) 469 call dgemm('N','N',mbas1(isymg),nrhf(isymk), 470 & norb2(isymp),one,work(koffl),ntotg, 471 & work(kr2pk2),ntotp,zero,work(koffs),ntotg1) 472 call daxpy(mbas1(isymg)*nrhf(isymk),one, 473 & work(koffs),1,work(kofft),1) 474 endif 475 476 end do ! loop over symmetries of gamma 477 478c ---------------------------------------------------- 479c store R(gamma k, l, delta) as R_(gamma k),(delta l) 480c on file: 481c (note: isympk = isymgk, isymql = isymdl) 482c ---------------------------------------------------- 483 if (inttyp .eq. 1) then 484 ndell = ir1bas(isymd,isyml) + 485 & mbas1(isymd)*(l-1) + idelta 486 iadr = ir2bas(isympk,isymql) + 487 & nr1bas(isympk)*(ndell -1) + 1 488 len = nr1bas(isympk) 489 call putwa2(lur2back,fback,work(kr2gk), 490 & iadr,len) 491 else if (inttyp .eq. 2) then 492 ndell = ir1xbas(isymd,isyml) + 493 & mbas2(isymd)*(l-1) + idelta - mbas1(isymd) 494 iadr = nr2bas + ir2xbas(isympk,isymql) + 495 & nr1bas(isympk)*(ndell -1) + 1 496 len = nr1bas(isympk) 497 call putwa2(lur2back,fback,work(kr2gk), 498 & iadr,len) 499 else if (inttyp .eq. 3) then 500 ndell = ir1bas(isymd,isyml) + 501 & mbas1(isymd)*(l-1) + idelta 502 iadr = ir2bas2(isympk,isymql) + 503 & nr1xbas(isympk)*(ndell -1) + 1 504 len = nr1xbas(isympk) 505 call putwa2(lur2back2,fback2,work(kr2gk), 506 & iadr,len) 507 else if (inttyp .eq. 4) then 508 ndell = ir1xbas(isymd,isyml) + 509 & mbas2(isymd)*(l-1) + idelta - mbas1(isymd) 510 iadr = nr2bas2 + ir2xbas2(isympk,isymql) + 511 & nr1xbas(isympk)*(ndell -1) + 1 512 len = nr1xbas(isympk) 513 call putwa2(lur2back2,fback2,work(kr2gk), 514 & iadr,len) 515 end if 516 if (locdbg) then 517 write(lupri,'(a,2i5,g20.10)') 518 & 'iadr,len,norm^2(R)=',iadr,len, 519 & ddot(len,work(kr2gk),1,work(kr2gk),1) 520 write(lupri,*) 'R^(l delta)_(gamma k):' 521 call output(work(kr2gk),1,len,1,1, 522 & len,1,1,lupri) 523 write(lupri,*)'nr2bas,ir2xbas',nr2bas, 524 & ir2xbas(isympk,isymql) 525 write(lupri,*)'ndell',ndell 526 write(lupri,*)'ir1xbas',ir1xbas(isymd,isyml) 527 write(lupri,*)'mbas2',mbas2(isymd) 528 write(lupri,*) 529 end if 530 531 end do ! loop over l 532 end if 533 end do ! loop over symmetry blocks (pk) 534 end do ! loop over delta 535 end do ! loop over symmetries of delta 536 537c----------------------------------------------------------------------- 538c close file, clean trace and return: 539c----------------------------------------------------------------------- 540 if ((inttyp .eq. 1) .or. (inttyp .eq. 2)) then 541 call wclose2(lur2back,fback,'KEEP') 542 else 543 call wclose2(lur2back2,fback2,'KEEP') 544 end if 545 546 end do ! loop over integral types 547 548 if(locdbg) write(lupri,*)'leaving cc_r12batf' 549 call qexit('r12batf') 550 551 return 552 end 553*=====================================================================* 554 subroutine cc_r12getrpkq(fback,rpkq,l,isyml,istartq,iendq,isymq, 555 & istartp,iendp,isymp, 556 & r2am,isympk,isymql,nr1orb,ir1orb) 557c---------------------------------------------------------------------- 558c Purpose: collect a three index batch r^l_{pk,q} with fixed l 559c and q and symmetry of p, ranging from istartq to iendq 560c and from istartp to iendp out of the 561c symmetry packed lower triangular matrix r_{pk,ql} 562c stored in r2am 563c 564c rpkq = r^l_{pk,q} with : p MOs ranging from istartp ... iendp 565c k occupied / r12 MOs 566c q MOs ranging from istartq ... iendq 567c 568c r2am = r_{pk,ql} with : p,q general MOs (orbital + aux. basis) 569c k,l occupied / r12 MOs 570c 571c H. Fliegl, C. Haettig spring 2003 572c modified C. Neiss autumn 2004 573c---------------------------------------------------------------------- 574 implicit none 575#include "priunit.h" 576#include "ccorb.h" 577#include "ccsdsym.h" 578#include "r12int.h" 579!for first order properties, app. B (Elena): 580#include "ccr12int.h" 581 582 logical locdbg 583 character*(*) fback !for first order properties, app. B (Elena) 584 parameter (locdbg=.false.) 585 integer ir1orb(8,8),nr1orb(8),npk,npk2,npkq,isymk,isymp 586 integer npkql,istartq,iendq,nql,istartp,iendp,norbp 587 integer isymq,isyml,isympk,isymql, iqloff 588 integer index 589#if defined (SYS_CRAY) 590 real rpkq(*), r2am(*) 591#else 592 double precision rpkq(*), r2am(*) 593#endif 594 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 595c 596 if (locdbg) then 597 write(lupri,*) 'ONEAUX = ',ONEAUX 598 end if 599 600c 601c number of p orbitals 602 norbp = iendp - istartp + 1 603c 604c symmetry of k 605 isymk = muld2h(isymp,isympk) 606 607 iqloff = nh1am(isymql) * (nh1am(isymql)+1) / 2 608 609 if (locdbg) then 610 write(lupri,*) 'istartq,iendq,istartp,iendp,l: ', 611 & istartq,iendq,istartp,iendp,l 612 write(lupri,*) 613 write(lupri,*) 'Input: r_{pk,ql}' 614 call outpak(r2am,nbast*nrhf(isymk),1,lupri) 615 end if 616 617c ------------------------------------------------ 618c loop over q and index k 619c ------------------------------------------------ 620 do q = istartq, iendq 621 622 if (.not.oneaux) then 623 nql = it1am(isymq,isyml) + nvir(isymq)*(l-1) + q 624 else if (q. le. norb1(isymq)) then 625 nql = ih1am(isymq,isyml) + norb1(isymq)*(l-1) + q 626 else 627 nql = ig1am(isymq,isyml) + norb2(isymq)*(l-1) + q-norb1(isymq) 628 end if 629 630 do k = 1, nrhf(isymk) 631 632 if (locdbg) then 633 write(lupri,*)'nql,isymk,k',nql,isymk,k 634 end if 635 636c --------------------------------------------------------- 637c reorder integrals, the following indeces are used: 638c npk : pair p,k with k occupied, p gen. MO (orb.+aux.) 639c npk2 : pair p,k with k occupied, p MO (orbital only) 640c npkql : quadrupel pk,ql as needed for r2am 641c npkq : tripel pk,q as needed for rpkq 642c --------------------------------------------------------- 643 if (isympk.eq.isymql) then 644 do p = istartp, iendp 645 if (oneaux) then 646 npk = ih1am(isymp,isymk)+norb1(isymp)*(k-1)+p 647 if (q. le. norb1(isymq)) then 648 npkql = ih2am(isymql,isympk)+index(nql,npk) 649 else 650 npkql = ih2am(isymql,isympk)+iqloff 651 * + nh1am(isympk)*(nql-1)+npk 652 end if 653 else 654 npk = it1am(isymp,isymk) + nvir(isymp)*(k-1)+p 655 if (fback.eq. fq12back) then 656 npkql = it2sq(isymql,isympk)+ nt1am(isymql) 657 & *(npk-1)+nql 658 elseif (fback.eq. fu12back) then 659 npkql = it2sq(isympk,isymql)+ nt1am(isympk) 660 & *(nql-1)+npk 661 else 662 npkql = it2am(isympk,isymql)+index(npk,nql) 663 endif 664 end if 665 npk2 = norbp*(k-1)+(p-istartp+1) 666 npkq = norbp*nrhf(isymk)*(q-istartq)+npk2 667 668 if (locdbg) then 669 write(lupri,'(a,5i5)')'p,npk,npk2,npkql,npkq', 670 & p,npk,npk2,npkql,npkq 671 call flshfo(lupri) 672 end if 673 674 rpkq(npkq) = r2am(npkql) 675 end do 676 else if (isympk.lt.isymql) then 677 if (oneaux) 678 & call quit('isympk.ne.isymql has not been implemented') 679 do p = istartp, iendp 680 npk = it1am(isymp,isymk) + nvir(isymp)*(k-1)+p 681 npk2 = norbp*(k-1)+(p-istartp+1) 682 npkql = it2am(isympk,isymql)+ nt1am(isympk)*(nql-1)+npk 683 npkq = norbp*nrhf(isymk)*(q-istartq)+npk2 684 rpkq(npkq) = r2am(npkql) 685 end do 686 else if (isympk.gt.isymql) then 687 if (oneaux) 688 & call quit('isympk.ne.isymql has not been implemented') 689 do p = istartp, iendp 690 npk = it1am(isymp,isymk) + nvir(isymp)*(k-1)+p 691 npk2 = norbp*(k-1)+(p-istartp+1) 692 npkql = it2am(isymql,isympk)+nt1am(isymql)*(npk-1)+nql 693 npkq = norbp*nrhf(isymk)*(q-istartq)+npk2 694 rpkq(npkq) = r2am(npkql) 695 end do 696 end if 697 end do 698 end do 699 700 end 701*=====================================================================* 702 subroutine cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 703 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas,ir2bas,ir2xbas, 704 & irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 705c---------------------------------------------------------------------- 706c set some symmetry offsets and dimensions needed only when 707c working in the R12 environment and the standard CC offset 708c cannot be used because some dimensions were overwritten 709c 710c H. Fliegl, C. Haettig spring 2003 711c 712c adapted for more general R12-MOs: 713c C. Neiss, summer 2005 714c---------------------------------------------------------------------- 715 implicit none 716#include "priunit.h" 717#include "ccorb.h" 718#include "ccsdsym.h" 719#include "r12int.h" 720 integer nr1orb(8),nr1xorb(8),nr1bas(8),nr1xbas(8),nr2bas,n2bst1(8) 721 integer ir1xbas(8,8),ir1xorb(8,8) 722 integer ir1orb(8,8),ir1bas(8,8),ir2bas(8,8),iaodis1(8,8) 723 integer ir2xbas(8,8),irgkl(8,8),nrgkl(8) 724 integer nrxgkl(8),irxgkl(8,8) 725 integer isymdk,isymk,isymd,isym,isym1,isym2 726 integer icount1,icount2,icount3,icount4,icount5,icount6 727 integer icount7,icount8,icount9,icount10,icount11 728 integer nalphaj(8),ialphaj(8,8) 729 logical locdbg 730 parameter (locdbg = .FALSE.) 731 732 call qenter('cc_r12offset') 733 if (locdbg) then 734 write(LUPRI,*) 'Entered CC_R12OFFSET' 735 call flshfo(LUPRI) 736 end if 737 738 do isymdk = 1, nsym 739 nr1orb(isymdk) = 0 740 nr1xorb(isymdk) = 0 741 nr1bas(isymdk) = 0 742 nr1xbas(isymdk) = 0 743 n2bst1(isymdk) = 0 744 nalphaj(isymdk) = 0 745 do isymk = 1, nsym 746 isymd = muld2h(isymdk,isymk) 747 nr1orb(isymdk) = nr1orb(isymdk) + norb1(isymd)*nrhfb(isymk) 748 nr1xorb(isymdk)= nr1xorb(isymdk)+ norb2(isymd)*nrhfb(isymk) 749 nr1bas(isymdk) = nr1bas(isymdk) + mbas1(isymd)*nrhfb(isymk) 750 nr1xbas(isymdk)= nr1xbas(isymdk)+ mbas2(isymd)*nrhfb(isymk) 751 n2bst1(isymdk) = n2bst1(isymdk) + mbas1(isymd)*mbas1(isymk) 752 nalphaj(isymdk)= nalphaj(isymdk)+ mbas1(isymd)*nrhfa(isymk) 753 end do 754 if (locdbg) then 755 write (lupri,*) 'nrhf, nrhfa, nrhfb:', 756 & nrhf(isymdk), nrhfa(isymdk), nrhfb(isymdk) 757 write (lupri,*) 'nr1orb, nr1xorb, nr1bas, nr1xbas: ', 758 & nr1orb(isymdk),nr1xorb(isymdk), 759 & nr1bas(isymdk),nr1xbas(isymdk),n2bst1(isymdk) 760 end if 761 end do 762 763c ------------------------------------ 764c three index array for r12-integrals 765c ------------------------------------ 766 do isymdk = 1, nsym 767 nrgkl(isymdk) = 0 768 nrxgkl(isymdk) = 0 769 nvajkl(isymdk) = 0 770 do isymk = 1, nsym 771 isymd = muld2h(isymdk,isymk) 772 nrgkl(isymdk) = nrgkl(isymdk) + mbas1(isymd)*nmatkl(isymk) 773 nrxgkl(isymdk) = nrxgkl(isymdk) + mbas2(isymd)*nmatkl(isymk) 774 nvajkl(isymdk) = nvajkl(isymdk) + nalphaj(isymd)*nmatkl(isymk) 775 end do 776 end do 777 778c 779c do isym = 1, nsym 780c nvabkl(isym) = 0 781c icount1 = 0 782c do isym2 = 1, nsym 783c isym1 = muld2h(isym,isym2) 784c nvabkl(isym) = nvabkl(isym) + n2bst(isym1)*nmatkl(isym2) 785c ivabkl(isym1,isym2) = icount1 786c icount1 = icount1 + n2bst(isym1)*nmatkl(isym2) 787c end do 788c end do 789 790c --------------------------------------- 791c lenghts for nr1bas over all symmetries 792c --------------------------------------- 793 nr2bas = 0 794 do isym = 1, nsym 795 nr2bas = nr2bas + nr1bas(isym)*nr1bas(isym) 796 end do 797 798c -------------------------------------------- 799c calculate now the offsets for r12 integrals 800c -------------------------------------------- 801 do isymdk = 1, nsym 802 icount1 = 0 803 icount2 = 0 804 icount3 = 0 805 icount4 = 0 806 icount5 = 0 807 icount6 = 0 808 icount7 = 0 809 icount8 = 0 810 icount9 = 0 811 icount10 = 0 812 icount11 = 0 813 do isymk = 1, nsym 814 isymd = muld2h(isymdk,isymk) 815 ir1orb(isymd,isymk) = icount1 !occ+vir 816 ir1xorb(isymd,isymk)= icount9 817 ir1bas(isymd,isymk) = icount2 !AO 818 ir1xbas(isymd,isymk)= icount8 819 ir2bas(isymd,isymk) = icount3 !AO+aux 820 ir2xbas(isymd,isymk)= icount4 !aux 821 irgkl(isymd,isymk) = icount5 822 irxgkl(isymd,isymk) = icount10 823 iaodis1(isymd,isymk)= icount6 824 ivajkl(isymd,isymk) = icount7 825 ialphaj(isymd,isymk)= icount11 826 icount1 = icount1 + nrhfb(isymk)*norb1(isymd) 827 icount9 = icount9 + nrhfb(isymk)*norb2(isymd) 828 icount2 = icount2 + nrhfb(isymk)*mbas1(isymd) 829 icount8 = icount8 + nrhfb(isymk)*mbas2(isymd) 830 icount3 = icount3 + nr1bas(isymd)*nr1bas(isymk) 831 icount4 = icount4 + nr1bas(isymd)*nr1xbas(isymk) 832 icount5 = icount5 + mbas1(isymd)*nmatkl(isymk) 833 icount6 = icount6 + mbas1(isymd)*mbas1(isymk) 834 icount7 = icount7 + nalphaj(isymd)*nmatkl(isymk) 835 icount10 = icount10 + mbas2(isymd)*nmatkl(isymk) 836 icount11 = icount11 + nrhfa(isymk)*mbas1(isymd) 837 end do 838 end do 839 if (locdbg) then 840 write(LUPRI,*) 'Leaving CC_R12OFFSET' 841 call flshfo(LUPRI) 842 end if 843 call qexit('cc_r12offset') 844 return 845 end 846*======================================================================* 847