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 cc_r12mkxajkl(kvajkl,cmo,work,lwork,fel) 20c--------------------------------------------------------------------- 21c purpose: calculate X_{aj}^{kl}, X = V,T necessary for the 22c calculation of MP2-R12 contributions to the right hand side 23c of Z-Vector Equations (cc_r12grad) 24c--------------------------------------------------------------------- 25 26 implicit none 27#include "priunit.h" 28#include "dummy.h" 29#include "ccorb.h" 30#include "ccsdsym.h" 31#include "r12int.h" 32#if defined (SYS_CRAY) 33 real zero,one 34#else 35 double precision zero,one 36#endif 37 logical ldum,test,fel 38 parameter (test = .false.,zero = 0.0d0, one = 1.0d0) 39 40 integer lwork,luvajkl,lwrk1,lwrk2,isymai 41 integer isymkl,isymij,isymmu,isyml,isymk,isymak,isymaj,icount7, 42 & ntota,isymv,isymmuj,koffl,koffv,kvmujkl, 43 & isymj,isyma,idxkl,xajkl,iglmrv(8,8) 44 integer kend1,kend2,koffc,isymi,ntklaj,yajkl 45 integer icou2,isymmua,ntotmu,nr1bas(8),ir1bas(8,8),icount2 46 integer noffro,ntota1,noffset,isym,nvir0(8) 47 integer idxij,idxji,kcmoa,idum,ntotgu,ntotj,lu43 48 integer nklaj(8),nr1vir(8),iglmro(8,8),icou3 49 character*10 fnelena 50 data fnelena/'DDR12xxxxx'/ 51 real*8 cmo(*),work(*),kvajkl(*) 52c 53 call qenter('cc_r12mkxajkl') 54 fnelena = fnelena(1:5)//fnvajkl(6:10) 55 56 if (fnvajkl .eq. 'CCR12QAJKL' .and. fel .eqv. .true.) then 57 lu43 = -43 58 call dzero(cmo,nlamds) 59 if (R12HYB) then 60 call gpopen(lu43,'GUMAT.1','UNKNOWN',' ','UNFORMATTED', 61 & IDUM,LDUM) 62 else 63 call quit('Sorry, calculation of R12 corrections to'// 64 & ' first order properties not implemented'// 65 & ' with .NO HYB') 66 endif 67 REWIND(LU43) 68 READ(LU43) NTOTGU 69 READ(LU43) (cmo(I), I = 1, NTOTGU) 70 CALL GPCLOSE(LU43,'KEEP') 71 endif 72c 73 if (fnvajkl .eq. 'CCR12QAJKL' .and. fel .eqv. .false.) then 74 kend1 = 1 75 lwrk1 = lwork - kend1 76 lu43 = -43 77 call dzero(cmo,nlamds) 78 call gpopen(lu43,'SIRIFC','OLD',' ','UNFORMATTED', 79 & IDUMMY,.FALSE.) 80 REWIND(LU43) 81 CALL MOLLAB(LABEL,LU43,LUPRI) 82 READ(LU43) 83 READ(LU43) 84 READ(LU43) (cmo(I), I = 1, NLAMDS) 85 CALL GPCLOSE(LU43,'KEEP') 86 CALL CMO_REORDER(CMO,WORK(KEND1),LWRK1) 87 88 endif 89c 90 91 do isymak = 1, nsym 92 nr1bas(isymak) = 0 93 nr1vir(isymak) = 0 94 do isymk = 1, nsym 95 isyma = muld2h(isymak,isymk) 96 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 97 nr1bas(isymak) = nr1bas(isymak) +mbas1(isyma)*nrhf(isymk) 98 nr1vir(isymak) = nr1vir(isymak) +nvir0(isyma)*nrhf(isymk) 99 end do 100 end do 101 102 do isymak = 1, nsym 103 nvajkl(isymak) = 0 104 icount7 = 0 105 icount2 = 0 106 do isymk = 1, nsym 107 isyma = muld2h(isymak,isymk) 108 ivajkl(isyma,isymk) = icount7 109 nvajkl(isymak) = nvajkl(isymak)+ nr1bas(isyma)*nmatij(isymk) 110 icount7 = icount7 + nr1bas(isyma)*nmatij(isymk) 111 ir1bas(isyma,isymk) = icount2 112 icount2 = icount2 + nrhf(isymk)*mbas1(isyma) 113 end do 114 end do 115 116 ntklaj = 0 117 do isymkl = 1,nsym 118 isymaj = isymkl 119 nklaj(isymaj) = 0 120 do isymj = 1,nsym 121 isyma = muld2h(isymaj,isymj) 122 nklaj(isymaj) = nklaj(isymaj) + nmatij(isymj)*nr1vir(isyma) 123 ntklaj = ntklaj + nmatij(isymkl)*nr1vir(isymj) 124 enddo 125 enddo 126 127 do isymmua = 1,nsym 128 icou2 = 0 129 icou3 = 0 130 do isyma = 1,nsym 131 isymmu = muld2h(isymmua,isyma) 132 iglmrv(isymmu,isyma) = icou2 133 iglmro(isymmu,isyma) = icou3 134 icou2 = icou2 + nbas(isymmu)*(nvir(isyma)-nrhffr(isyma)) 135 icou3 = icou3 + nbas(isymmu)*(norb1(isyma)-nrhffr(isyma)) 136 enddo 137 enddo 138 139 140 noffro = 0 141 do isym = 1,nsym 142 noffro = noffro + nbas(isym)*nrhf(isym) 143 enddo 144 145 146 kend1 = 1 147 yajkl = kend1 148 kcmoa = yajkl + ntklaj 149 kend2 = kcmoa + nlamds 150 lwrk2 = lwork - kend2 151 152 if (lwrk2 .lt.0) then 153 call quit('Insufficient work space in cc_r12mkxajkl') 154 end if 155c 156c 157 call dzero(work(yajkl),ntklaj) 158 xajkl = yajkl 159 do isymkl = 1, nsym 160 isymaj = isymkl 161 isymmuj = isymkl 162 do isyml =1, nsym 163 isymk = muld2h(isymkl,isyml) 164 do k = 1, nrhf(isymk) 165 do l = 1, nrhf(isyml) 166 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 167 noffset = 0 168 do isymmu =1, nsym 169 isyma = muld2h(1,isymmu) 170 isymj = muld2h(isymmuj,isymmu) 171 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 172 if (fel) then 173 koffc = 1 + iglmro(isymmu,isyma) 174 & +nrhf(isymmu)*nbas(isymmu) 175 else 176 noffset = noffset+ nbas(isymmu)*nrhffr(isymmu) 177 koffc = 1 + nrhf(isymmu)*nbas(isymmu) 178 & +iglmrv(isymmu,isyma)+noffset 179 & +noffro 180 endif 181 kvmujkl = 1+ivajkl(isymmuj,isymkl)+ 182 & nr1bas(isymmuj)*(idxkl-1) 183 & +ir1bas(isymmu,isymj) 184 ntotmu = max(1,mbas1(isymmu)) 185 ntota1 = max(1,nbas(isymmu)) 186 ntota = max(1,nvir0(isyma)) 187 ntotj = max(1,nrhf(isymj)) 188c 189 call dgemm('T','N',nvir0(isyma), 190 & nrhf(isymj),mbas1(isymmu),one,cmo(koffc), 191 & ntota1,kvajkl(kvmujkl),ntotmu,zero, 192 & work(xajkl),ntota) 193 194 xajkl = xajkl + nvir0(isyma)*nrhf(isymj) 195 end do 196 end do 197 end do 198 end do 199 end do 200 201 if (r12cbs) call dscal(nvajkl(1),1.00d0,work(yajkl),1) 202 203 204 luvajkl = -1 205 call gpopen(luvajkl,fnelena,'unknown',' ','unformatted', 206 & idummy,.false.) 207 rewind(luvajkl) 208 write(luvajkl) (work(yajkl+i-1), i = 1,ntklaj) 209 call gpclose(luvajkl,'KEEP') 210 211 212 call qexit('cc_r12mkxajkl') 213 end 214C=====================================================================* 215C /* Deck cc_r12grad */ 216 subroutine cc_r12grad(ipdd,orb,v2am,work,lwork) 217C Calculate MP2-R12 contributions to the right hand side of Z-Vector 218C Equations 219 implicit none 220#include "ccsdsym.h" 221#include "ccorb.h" 222#include "dummy.h" 223#include "cbirea.h" 224#include "r12int.h" 225#include "priunit.h" 226#include "ccfro.h" 227 logical locdbg 228 logical lauxd,ldum 229 parameter (locdbg = .false.) 230 character*10 CCR12YAJIJ 231 integer kend1,kccr12,kend2,lwrk1,lwork,nrhftria,isymij 232 integer lu43,idum,n2,lwrk2,ipdd 233 integer krsing,krtrip,kend3,lwrk3,nr1bas(8) 234 integer isymkl,isyml,isymk,idxij,idxji,idxklij,idxklji 235 integer isymai,isymv,isymi,isymj,idxkl,kctil,ntklaj 236 integer kbajkl,kbijal,kvajkl,kvijal,isymaj,isyma,luvajkl 237 integer lcvai,koffc,nvir0(8),isymmu,ntota,ntotmu,kcvai 238 integer ir1bas(8,8),isymmuj 239 integer dklij,ntotk,lcvia,kvmujkl,ntoti 240 integer koffd,bdajij,lccr12 241 integer dklijz,kajij,kcvia,ntotj,imataj(8,8),icou3 242 integer nr1orb(8),nr1xbas(8),nrgkl,n2bst1(8),ir1orb(8,8) 243 integer ir2bas(8,8),ir2xbas(8,8),irgkl,iaodis1(8,8),nr2bas(8) 244 integer ir1xbas(8,8),nijkl(8) 245 integer koffe,kctil2,idxijkl,lctil,ncvai 246 integer res,kscr1 247 248#if defined (SYS_CRAY) 249 real work(lwork),one,two,zero 250 real orb(*),v2am(*) 251#else 252 double precision work(lwork),one,zero,two,orb(*),v2am(*) 253#endif 254 parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0) 255 256 257 call qenter('cc_r12grad') 258 259 260 nrhftria = nrhft*(nrhft+1)/2 261 n2 = nrhftria * nrhftria 262 263 kccr12 = 1 264 krsing = kccr12 + ngamsq(1) 265 krtrip = krsing + n2 266 kctil = krtrip + n2 267 kctil2 = kctil + ngamsq(1) 268 kend1 = kctil2 + ngamsq(1) 269 lwrk1 = lwork - kend1 270 271 call dzero(work(kccr12),ngamsq(1)) 272 call dzero(work(krsing),n2) 273 call dzero(work(krtrip),n2) 274 call dzero(work(kctil),ngamsq(1)) 275 call dzero(work(kctil2),ngamsq(1)) 276 277 isymv = 1 278 279 if (lwrk1 .lt. 0) then 280 call quit('Insufficient work space in cc_r12grad') 281 END IF 282 283c ---------------------------------------- 284c read R12 amplitudes 285c ---------------------------------------- 286 287 288 lu43 = -43 289 call gpopen(lu43,'CCR12_C','unknown',' ','formatted', 290 & idummy,.false.) 291 read(lu43,'(4e30.20)') (work(krsing+i), i = 0, n2-1 ) 292 read(lu43,'(4e30.20)') (work(krtrip+i), i = 0, n2-1 ) 293 call gpclose(lu43,'KEEP') 294 295 call dscal(nrhftria*nrhftria,0.5d0,work(krsing),1) 296 call dscal(nrhftria*nrhftria,0.5d0,work(krtrip),1) 297 298 call cc_r12vpack(work(kccr12),isymv,work(krsing),work(krtrip), 299 & nrhf,nrhfa,nijkl) 300 301c --------------------------------- 302C get c_tilde = 2*c_ij^kl - c_ji^kl 303c --------------------------------- 304 305 do isymkl = 1, nsym 306 isymij = muld2h(isymkl,isymv) 307 do isyml =1, nsym 308 isymk = muld2h(isymkl,isyml) 309 do k = 1, nrhf(isymk) 310 do l = 1, nrhf(isyml) 311 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 312 do isymi = 1, nsym 313 isymj = muld2h(isymij,isymi) 314 do j = 1, nrhf(isymj) 315 do i = 1, nrhf(isymi) 316 idxij = imatij(isymi,isymj)+ 317 & nrhf(isymi)*(j-1) + i 318 idxji = imatij(isymj,isymi)+ 319 & nrhf(isymj)*(i-1) + j 320 idxklij = igamsq(isymkl,isymij)+ 321 & nmatij(isymkl)*(idxij-1)+idxkl 322 idxklji =igamsq(isymkl,isymij)+ 323 & nmatij(isymkl)*(idxji-1)+idxkl 324 work(kctil-1+idxklij) = 325 & two*work(kccr12-1+idxklij) - 326 & work(kccr12-1+idxklji) 327 end do 328 end do 329 330 end do 331 end do 332 end do 333 end do 334 end do 335 336 if (locdbg) then 337 write(lupri,*) 'resorted R12 amplitudes (ctilde):' 338 do isymkl = 1, nsym 339 isymij = muld2h(isymkl,1) 340 write(lupri,*) 'symmetry block ',isymkl,isymij 341 call output(work(kctil+igamsq(isymkl,isymij)), 342 & 1,nmatij(isymkl),1,nmatij(isymij), 343 & nmatij(isymkl),nmatij(isymij),1,lupri) 344 end do 345 endif 346 347c --------------------------------- 348C Transpone c_tilde 349c --------------------------------- 350 351 do isymkl = 1, nsym 352 isymij = muld2h(isymkl,isymv) 353 do isyml =1, nsym 354 isymk = muld2h(isymkl,isyml) 355 do k = 1, nrhf(isymk) 356 do l = 1, nrhf(isyml) 357 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 358 do isymi = 1, nsym 359 isymj = muld2h(isymij,isymi) 360 do j = 1, nrhf(isymj) 361 do i = 1, nrhf(isymi) 362 idxij = imatij(isymi,isymj)+ 363 & nrhf(isymi)*(j-1) + i 364 idxji = imatij(isymj,isymi)+ 365 & nrhf(isymj)*(i-1) + j 366 idxklij = igamsq(isymij,isymkl)+ 367 & nmatij(isymkl)*(idxij-1)+idxkl 368 idxklji =igamsq(isymij,isymkl)+ 369 & nmatij(isymkl)*(idxji-1)+idxkl 370 idxijkl = igamsq(isymkl,isymij)+ 371 & nmatij(isymij)*(idxkl-1) +idxij 372 work(kctil2-1+idxklij) = 373 & work(kctil-1+idxijkl) 374 end do 375 end do 376 377 end do 378 end do 379 end do 380 end do 381 end do 382 383 384c 385C--------------------------------------------- 386C Calculate some offsets and dimensions 387C--------------------------------------------- 388 ntklaj = 0 389 do isymkl = 1, nsym 390 isymaj = isymkl 391 do isyml =1, nsym 392 isymk = muld2h(isymkl,isyml) 393 do k = 1, nrhf(isymk) 394 do l = 1, nrhf(isyml) 395 do isyma =1, nsym 396 isymj = muld2h(isymaj,isyma) 397 isymi = isyma 398 ntklaj = ntklaj + nrhf(isymj) * 399 & (norb1(isyma) - nrhfs(isyma)) 400 enddo 401 enddo 402 enddo 403 enddo 404 enddo 405 406 407 icou3 = 0 408 do isymai = 1,nsym 409 icou3 = 0 410 do isyma = 1,nsym 411 isymi = muld2h(isymai,isyma) 412 imataj(isyma,isymi)= icou3 413 icou3 = icou3+nrhf(isymi)*(norb1(isyma) - nrhfs(isyma)) 414 enddo 415 enddo 416 417 do isymaj = 1,nsym 418 isymij = isymaj 419 ncvai = 0 420 do isyma = 1,nsym 421 isymj = muld2h(isymaj,isyma) 422 isymi = muld2h(isymij,isymj) 423 ncvai = ncvai + (norb1(isyma) - nrhfs(isyma))*nrhf(isymi) 424 enddo 425 enddo 426 427 428 kbajkl = kend1 429 kbijal = kbajkl + ntklaj 430 kvajkl = kbijal + ntklaj 431 kvijal = kvajkl + ntklaj 432 kcvai = kvijal + ntklaj 433 kcvia = kcvai + ncvai 434 kajij = kcvia + ncvai 435 dklij = kajij + ncvai 436 res = dklij + ngamsq(1) 437 kend2 = res + ncvai 438 lwrk2 = lwork - kend2 439 440c ---------------------------------------------------- 441c Read V_{aj}^{kl},V_{kl}^{aj},B_{aj}^{kl},B_{kl}^{aj} 442c and transform it with Ctilde 443c ----------------------------------------------------- 444 445 if (lwrk2 .lt.0) then 446 call quit('Insufficient work space in cc_r12grad') 447 end if 448 449 call dzero(work(kbajkl),ntklaj) 450 call dzero(work(kbijal),ntklaj) 451 luvajkl = -1 452 call gpopen(luvajkl,'DDR12BAJKL','unknown',' ','unformatted', 453 & idummy,.false.) 454 rewind(luvajkl) 455 read(luvajkl) (work(kbajkl+i-1), i = 1,ntklaj) 456 if (ipdd .eq. 5) then 457 call gpclose(luvajkl,'DELETE') 458 else 459 call gpclose(luvajkl,'KEEP') 460 endif 461 462 463 call gpopen(luvajkl,'DDR12BIJAL','unknown',' ','unformatted', 464 & idummy,.false.) 465 rewind(luvajkl) 466 read(luvajkl) (work(kbijal+i-1), i = 1,ntklaj) 467 if (ipdd .eq. 5) then 468 call gpclose(luvajkl,'DELETE') 469 else 470 call gpclose(luvajkl,'KEEP') 471 endif 472c 473 call daxpy(ntklaj,one,work(kbajkl),1,work(kbijal),1) 474 475 call dzero(work(kvajkl),ntklaj) 476 call dzero(work(kvijal),ntklaj) 477 call gpopen(luvajkl,'DDR12VAJKL','unknown',' ','unformatted', 478 & idummy,.false.) 479 rewind(luvajkl) 480 read(luvajkl) (work(kvajkl+i-1), i = 1,ntklaj) 481 if (ipdd .eq. 5) then 482 call gpclose(luvajkl,'DELETE') 483 else 484 call gpclose(luvajkl,'KEEP') 485 endif 486 487 call gpopen(luvajkl,'DDR12VIJAL','unknown',' ','unformatted', 488 & idummy,.false.) 489 rewind(luvajkl) 490 read(luvajkl) (work(kvijal+i-1), i = 1,ntklaj) 491 if (ipdd .eq. 5) then 492 call gpclose(luvajkl,'DELETE') 493 else 494 call gpclose(luvajkl,'KEEP') 495 endif 496 497 498 499c Calculate d_{kl}^{ij} = c_{kl}^{mn}*ctilde_{ij}^{mn} 500 501 call dzero(work(dklij),ngamsq(1)) 502 do isymkl = 1,nsym 503 isymij = isymkl 504 lccr12 = kccr12 + igamsq(isymij,isymkl) 505 lctil = kctil + igamsq(isymij,isymkl) 506 dklijz = dklij + igamsq(isymij,isymkl) 507 ntoti = max(1,nmatij(isymij)) 508 ntotk = max(1,nmatij(isymkl)) 509 call dgemm('T','N',nmatij(isymij),nmatij(isymkl), 510 & nmatij(isymkl),one,work(lccr12), 511 & ntotk,work(lctil),ntoti,zero, 512 & work(dklijz),ntotk) 513 end do 514 515 516 call dzero(work(kcvai),ncvai) 517 lcvai = kcvai 518 call dzero(work(kcvia),ncvai) 519 lcvia = kcvia 520 call dzero(work(kajij),ncvai) 521 bdajij = kajij 522 523c Calculate V_{aj}^{kl}*ctilde_{kl}^{ij} 524 525 526 do isymkl = 1, nsym 527 isymaj = isymkl 528 isymij = isymkl 529 do isyml =1, nsym 530 isymk = muld2h(isymkl,isyml) 531 do k = 1, nrhf(isymk) 532 do l = 1, nrhf(isyml) 533 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 534 do isyma =1, nsym 535 isymj = muld2h(isymaj,isyma) 536 isymi = muld2h(isymij,isymj) 537 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 538 koffc = kctil + igamsq(isymij,isymkl) 539 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 540 ntotmu = max(1,nvir0(isyma)) 541 ntotj = max(1,nrhf(isymj)) 542 ntoti = max(1,nrhf(isymi)) 543 lcvai = kcvai + imataj(isyma,isymi) 544 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 545 & nrhf(isymj),one,work(kvajkl), 546 & ntotmu,work(koffc),ntoti,one, 547 & work(lcvai),ntotmu) 548 549 kvajkl = kvajkl + nvir0(isyma)*nrhf(isymj) 550 551c Calculate V_{kl}^{aj}*ctilde_{ij}^{kl} 552 553 lcvia = kcvia + imataj(isyma,isymi) 554 koffe = kctil2+ igamsq(isymkl,isymij) 555 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 556 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 557 & nrhf(isymj),one,work(kvijal), 558 & ntotmu,work(koffe),ntoti,one, 559 & work(lcvia),ntotmu) 560 kvijal = kvijal + nvir0(isyma)*nrhf(isymj) 561c 562c Calculate B_{aj}^{kl}*d_{kl}^{ij} 563c 564 koffd = dklij + igamsq(isymij,isymkl) 565 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 566 bdajij = kajij + imataj(isyma,isymi) 567 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 568 & nrhf(isymj),one,work(kbijal), 569 & ntotmu,work(koffd),ntoti,one, 570 & work(bdajij),ntotmu) 571 kbijal = kbijal + nvir0(isyma)*nrhf(isymj) 572 end do 573 end do 574 end do 575 end do 576 end do 577 call dscal(ncvai,two,work(kcvai),1) 578 call dscal(ncvai,two,work(kcvia),1) 579 call daxpy(ncvai,one,work(kcvai),1,work(kcvia),1) 580 call daxpy(ncvai,one,work(kcvia),1,work(kajij),1) 581 call dscal(ncvai,two,work(kajij),1) 582 if (.not. lmulbs) call dscal(ncvai,-one,work(kajij),1) 583 584 if (locdbg) then 585 write(lupri,*) 'Matrixelemente MP2-R12 Gradient',ipdd 586 do isymaj = 1,nsym 587 do isyma = 1,nsym 588 isymi =isyma 589 write(lupri,*) 'symmetry block',imataj(isyma,isymi) 590 call output(work(kajij+imataj(isyma,isymi)), 591 & 1,nvir0(isyma),1,nrhf(isymi), 592 & nvir0(isyma),nrhf(isymi),1,lupri) 593 594 enddo 595 enddo 596 end if 597 598 if (ipdd .eq. 2) then 599 luvajkl = -1 600 call gpopen(luvajkl,'CCR12YAJIJ','unknown',' ','unformatted', 601 & idummy,.false.) 602 write(luvajkl) (work(kajij+i-1), i = 1,ncvai) 603 call gpclose(luvajkl,'KEEP') 604 else if (ipdd .eq. 3) then 605 call dzero(work(res),ncvai) 606 call r12gradap(ipdd,work(res),orb,work(kend2),v2am,lwrk2) 607 call dscal(ncvai,2.0D0,work(res),1) 608 call daxpy(ncvai,one,work(kajij),1,work(res),1) 609 luvajkl = -1 610 call gpopen(luvajkl,'CCR12ZAJIJ','unknown',' ','unformatted', 611 & idummy,.false.) 612 write(luvajkl) (work(res+i-1), i = 1,ncvai) 613 call gpclose(luvajkl,'KEEP') 614 else if (ipdd .eq. 5) then 615 call dzero(work(res),ncvai) 616 call r12gradap(ipdd,work(res),orb,work(kend2),v2am,lwrk2) 617 call dscal(ncvai,2.00D0,work(res),1) 618 call daxpy(ncvai,one,work(kajij),1,work(res),1) 619 luvajkl = -1 620 call gpopen(luvajkl,'CCR12XAJIJ','unknown',' ','unformatted', 621 & idummy,.false.) 622 write(luvajkl) (work(res+i-1), i = 1,ncvai) 623 call gpclose(luvajkl,'KEEP') 624 625 endif 626c 627c 628 call qexit('cc_r12grad') 629 end 630c------------------------------------------------------------------------ 631C=====================================================================* 632C /* Deck r12gradap */ 633 subroutine r12gradap(ipdd,res,orb,work,v2am,lwork) 634 implicit none 635#include "ccsdsym.h" 636#include "ccorb.h" 637#include "dummy.h" 638#include "r12int.h" 639#include "priunit.h" 640#include "ccfro.h" 641#include "ccsdinp.h" 642 integer kccr12,krsing,krtrip,kctil,kend1,lwrk1,nrhftria,n2 643 integer lwork,lu43,idxkl,isymv,isymi,isymj,isymk,idxij 644 integer isyml,isymkl,isymij,idxji,idxklij,idxijkl,idxklji 645 integer kxajkl,lunit,luvajkl,ntklaj,isyma,isymaj,kend2 646 integer nmatan(8),isyman,isymn,ntotij,ntotan,nvir0(8) 647 integer lctil,klijan,ntota,ntoti,ntotn,lcvai,isymai 648 integer kijan,dklij,dklijz,ntotk,koffc,lccr12,koffd 649 integer idxijlk,idxlk,isymin,kbajkl,ncvai,imatan(8,8),icou3 650 integer lwrk2,lctil1,klij,kcotil,kdotil,koccr,koctil 651 integer idxlkji,idxlkij,idxjilk,krklaj,krajkl,ntaj 652 integer kxsing,kxtrip,kxxr12,klijmn,lxxr12,ij,al,idx 653 integer isym,nvir0t,blajkl,klijz,krajklz,klijanz,kxajklz 654 integer ipdd,kxtil,koffk,koffl,koffj,koffi,isymmu,isymmuj 655 integer icou4,imataj(8,8),blajklf,kl,aj,ajkl,nrhfst 656 integer imatmn(8,8),icou6,nmatkl1(8),imatkl1(8,8),icou5 657 integer imatka(8,8),icou7,nmatak(8),nijkl(8) 658#if defined (SYS_CRAY) 659 real res(*),work(*),orb(*),v2am(*),one,two,zero 660#else 661 double precision work(*),orb(*),v2am(*),one,zero,two 662 double precision res(*) 663#endif 664 parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0) 665 666 call qenter('r12gradap') 667 668 isymv = 1 669 670 nrhftria = nrhft*(nrhft+1)/2 671 n2 = nrhftria * nrhftria 672 673 nrhfst = 0 674 do isymi = 1,nsym 675 nrhfst = nrhfst + nrhfs(isymi) 676 enddo 677 678 kccr12 = 1 679 krsing = kccr12 + ngamsq(1) 680 krtrip = krsing + n2 681 kctil = krtrip + n2 682 kend1 = kctil + ngamsq(1) 683 lwrk1 = lwork - kend1 684 685 call dzero(work(krsing),2*n2) 686 call dzero(work(kccr12),ngamsq(1)) 687 call dzero(work(kctil),ngamsq(1)) 688 689 if (lwrk1 .lt. 0) then 690 call quit('Insufficient work space in r12gradap') 691 end if 692 693c ---------------------------------------- 694c read R12 amplitudes 695c ---------------------------------------- 696 lu43 = -43 697 call gpopen(lu43,'CCR12_C','unknown',' ','formatted', 698 & idummy,.false.) 699 read(lu43,'(4e30.20)') (work(krsing+i), i = 0, n2-1 ) 700 read(lu43,'(4e30.20)') (work(krtrip+i), i = 0, n2-1 ) 701 call gpclose(lu43,'KEEP') 702 703 704 call dscal(nrhftria*nrhftria,0.5d0,work(krsing),1) 705 call dscal(nrhftria*nrhftria,0.5d0,work(krtrip),1) 706 707 call cc_r12vpack(work(kccr12),isymv,work(krsing),work(krtrip), 708 & nrhf,nrhfa,nijkl) 709 710c --------------------------------- 711C get c_tilde = 2*c_ij^kl - c_ji^kl 712c --------------------------------- 713 714 do isymkl = 1, nsym 715 isymij = muld2h(isymkl,isymv) 716 do isyml =1, nsym 717 isymk = muld2h(isymkl,isyml) 718 do k = 1, nrhf(isymk) 719 do l = 1, nrhf(isyml) 720 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 721 do isymi = 1, nsym 722 isymj = muld2h(isymij,isymi) 723 do j = 1, nrhf(isymj) 724 do i = 1, nrhf(isymi) 725 idxij = imatij(isymi,isymj)+ 726 & nrhf(isymi)*(j-1) + i 727 idxji = imatij(isymj,isymi)+ 728 & nrhf(isymj)*(i-1) + j 729 idxklij = igamsq(isymij,isymkl)+ 730 & nmatij(isymkl)*(idxij-1)+idxkl 731 idxklji =igamsq(isymkl,isymij)+ 732 & nmatij(isymkl)*(idxji-1)+idxkl 733 idxijkl = igamsq(isymkl,isymij)+ 734 & nmatij(isymij)*(idxkl-1) +idxij 735 work(kctil-1+idxklij) = 736 & two*work(kccr12-1+idxklij) - 737 & work(kccr12-1+idxklji) 738 end do 739 end do 740 741 end do 742 end do 743 end do 744 end do 745 end do 746 747 748 749 ntklaj = 0 750 do isymkl = 1, nsym 751 isymaj = isymkl 752 do isyml =1, nsym 753 isymk = muld2h(isymkl,isyml) 754 do k = 1, nrhf(isymk) 755 do l = 1, nrhf(isyml) 756 do isyma =1, nsym 757 isymj = muld2h(isymaj,isyma) 758 isymi = isyma 759 ntklaj = ntklaj + nrhf(isymj) * 760 & (norb1(isyma) - nrhfs(isyma)) 761 enddo 762 enddo 763 enddo 764 enddo 765 enddo 766 767 do isyman = 1,nsym 768 icou3 = 0 769 icou4 = 0 770 icou5 = 0 771 do isyma = 1,nsym 772 isymn = muld2h(isyman,isyma) 773 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 774 icou3 = icou3 + nvir0(isyma)*nrhfs(isymn) 775 icou4 = icou4 + nrhfs(isyma)*nrhfs(isymn) 776 icou5 = icou5 + nvir0(isyma)*nrhf(isymn) 777 enddo 778 nmatan(isyman) = icou3 779 nmatkl1(isyman) = icou4 780 nmatak(isyman) = icou5 781 enddo 782 783 icou3 = 0 784 do isyman = 1,nsym 785 icou3 = 0 786 icou4 = 0 787 icou5 = 0 788 icou6 = 0 789 icou7 = 0 790 do isyma = 1,nsym 791 isymn = muld2h(isyman,isyma) 792 imatan(isyma,isymn)= icou3 793 imataj(isymn,isyma)= icou4 794 imatkl1(isyma,isymn)= icou5 795 imatmn(isymn,isyma)= icou6 796 imatka(isymn,isyma)= icou7 797 icou3 = icou3+nrhf(isymn)*(norb1(isyma) - nrhfs(isyma)) 798 icou4 = icou4+nrhfs(isymn)*(norb1(isyma) - nrhfs(isyma)) 799 icou5 = icou5+nmatan(isymn)* nmatkl1(isyma) 800 icou6 = icou6+nrhfs(isyma)*nrhfs(isymn) 801 icou7 = icou7+nmatak(isyma)*nmatij(isymn) 802 enddo 803 enddo 804 do isyman = 1,nsym 805 do isyma = 1,nsym 806 isymn = muld2h(isyma,isyman) 807 enddo 808 enddo 809 810 do isymaj = 1,nsym 811 isymij = isymaj 812 ncvai = 0 813 do isyma = 1,nsym 814 isymj = muld2h(isymaj,isyma) 815 isymi = muld2h(isymij,isymj) 816 ncvai = ncvai + (norb1(isyma) - nrhfs(isyma))*nrhf(isymi) 817 enddo 818 enddo 819 820 ntaj = 0 821 do isymkl = 1,nsym 822 isymij = isymkl 823 do isymj = 1,nsym 824 isyma = muld2h(isymij,isymj) 825 ntaj = ntaj + nrhf(isymj) 826 & *(norb1(isyma) - nrhfs(isyma)) 827 enddo 828 enddo 829 830 nvir0t = 0 831 do isym = 1, nsym 832 nvir0(isym) = norb1(isym) - nrhfs(isym) 833 nvir0t = nvir0t + nvir0(isym) 834 end do 835 836 837 kxajkl = kend1 838 klijan = kxajkl + ntklaj 839 dklij = klijan + ncvai 840 koctil = dklij + ngamsq(1) 841 kdotil = koctil + ngamsq(1) 842 klij = kdotil + ngamsq(1) 843 krajkl = klij + ngamsq(1) 844 kxxr12 = krajkl + nvir0t*nrhfst**3 845 kxsing = kxxr12 + ngamsq(1) 846 kxtrip = kxsing + n2 847 klijmn = kxtrip + n2 848 blajkl = klijmn + ngamsq(1) 849 kxtil = blajkl + nvir0t*nrhfst**3 850 blajklf = kxtil + ngamsq(1) 851 kend2 = blajklf + ntklaj 852 lwrk2 = lwork - kend2 853 854 855 856 857 call dzero(work(koctil),ngamsq(1)) 858 do isymkl = 1, nsym 859 isymij = muld2h(isymkl,isymv) 860 do isyml =1, nsym 861 isymk = muld2h(isymkl,isyml) 862 do k = 1, nrhf(isymk) 863 do l = 1, nrhf(isyml) 864 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 865 idxlk=imatij(isyml,isymk)+nrhf(isyml)*(k-1)+l 866 do isymi = 1, nsym 867 isymj = muld2h(isymij,isymi) 868 do j = 1, nrhf(isymj) 869 koffj = irhf(isymj)+j 870 do i = 1, nrhf(isymi) 871 ij = (i-1)*nrhf(isymj)+j 872 koffi = irhf(isymi)+i 873 idxij = imatij(isymi,isymj)+ 874 & nrhf(isymi)*(j-1) + i 875 idxji = imatij(isymj,isymi)+ 876 & nrhf(isymj)*(i-1) + j 877 idxklij = igamsq(isymij,isymkl)+ 878 & nmatij(isymkl)*(idxij-1)+idxkl 879 idxklji =igamsq(isymkl,isymij)+ 880 & nmatij(isymkl)*(idxji-1)+idxkl 881 idxijkl = igamsq(isymkl,isymij)+ 882 & nmatij(isymij)*(idxkl-1) +idxij 883 idxijlk = igamsq(isymij,isymkl)+ 884 & nmatij(isymij)*(idxlk-1) +idxij 885 work(koctil-1+idxijlk) = 886 & work(koctil-1+idxijlk) -two* 887 & work(kctil-1+idxijlk) * orb(koffi) 888 work(koctil-1+idxijlk) = 889 & work(koctil-1+idxijlk) -two* 890 & work(kctil-1+idxijlk) * orb(koffj) 891 end do 892 end do 893 894 end do 895 end do 896 end do 897 end do 898 end do 899 900c Calculate d_{kl}^{ij} = c_{kl}^{mn}*ctilde_{ij}^{mn} 901 902 call dzero(work(dklij),ngamsq(1)) 903 call dzero(work(klij),ngamsq(1)) 904 do isymkl = 1,nsym 905 isymij = isymkl 906 lccr12 = kccr12 + igamsq(isymij,isymkl) 907 lctil = koctil + igamsq(isymij,isymkl) 908 lctil1 = kctil + igamsq(isymij,isymkl) 909 dklijz = dklij + igamsq(isymij,isymkl) 910 klijz = klij + igamsq(isymij,isymkl) 911 ntoti = max(1,nmatij(isymij)) 912 ntotk = max(1,nmatij(isymkl)) 913 call dgemm('T','N',nmatij(isymij),nmatij(isymkl), 914 & nmatij(isymkl),one,work(lccr12), 915 & ntotk,work(lctil),ntoti,zero, 916 & work(dklijz),ntotk) 917 call dgemm('T','N',nmatij(isymij),nmatij(isymkl), 918 & nmatij(isymkl),one,work(lccr12), 919 & ntotk,work(lctil1),ntoti,zero, 920 & work(klijz),ntotk) 921 end do 922 923 924 925 call dzero(work(kdotil),ngamsq(1)) 926 do isymkl = 1, nsym 927 isymij = muld2h(isymkl,isymv) 928 do isyml =1, nsym 929 isymk = muld2h(isymkl,isyml) 930 do k = 1, nrhf(isymk) 931 koffk = irhf(isymk)+k 932 do l = 1, nrhf(isyml) 933 koffl = irhf(isyml)+l 934 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 935 idxlk=imatij(isyml,isymk)+nrhf(isyml)*(k-1)+l 936 do isymi = 1, nsym 937 isymj = muld2h(isymij,isymi) 938 do j = 1, nrhf(isymj) 939 koffj = irhf(isymj)+j 940 do i = 1, nrhf(isymi) 941 koffi = irhf(isymi)+i 942 idxij = imatij(isymi,isymj)+ 943 & nrhf(isymi)*(j-1) + i 944 idxji = imatij(isymj,isymi)+ 945 & nrhf(isymj)*(i-1) + j 946 idxklij = igamsq(isymij,isymkl)+ 947 & nmatij(isymkl)*(idxij-1)+idxkl 948 idxklji =igamsq(isymkl,isymij)+ 949 & nmatij(isymkl)*(idxji-1)+idxkl 950 idxijkl = igamsq(isymkl,isymij)+ 951 & nmatij(isymij)*(idxkl-1) +idxij 952 idxijlk = igamsq(isymkl,isymij)+ 953 & nmatij(isymij)*(idxlk-1) +idxij 954 work(kdotil-1+idxijlk) = 955 & work(kdotil-1+idxijlk) + 956 & work(klij-1+idxijlk) * orb(koffi) 957 work(kdotil-1+idxijlk) = 958 & work(kdotil-1+idxijlk) + 959 & work(klij-1+idxijlk) * orb(koffj) 960 work(kdotil-1+idxijlk) = 961 & work(kdotil-1+idxijlk) + 962 & work(klij-1+idxijlk) * orb(koffk) 963 work(kdotil-1+idxijlk) = 964 & work(kdotil-1+idxijlk) + 965 & work(klij-1+idxijlk) * orb(koffl) 966 end do 967 end do 968 969 end do 970 end do 971 end do 972 end do 973 end do 974 975 call daxpy(ngamsq(1),one,work(dklij),1,work(kdotil),1) 976 977 978 luvajkl = -1 979 call dzero(work(kxajkl),ntklaj) 980 call gpopen(luvajkl,'DDR12XAJKL','unknown',' ','unformatted', 981 & idummy,.false.) 982 rewind(luvajkl) 983 read(luvajkl) (work(kxajkl+i-1), i = 1,ntklaj) 984 call gpclose(luvajkl,'KEEP') 985 986 987 988 call dzero(work(krajkl),nvir0t*nrhfst**3) 989 luvajkl = -1 990 call gpopen(luvajkl,'AUXQA12','unknown',' ','formatted', 991 & idummy,.false.) 992 rewind(luvajkl) 993 read(luvajkl,'(4E30.20)') (work(krajkl+i-1), i=1,nvir0t*nrhfst**3) 994 call gpclose(luvajkl,'KEEP') 995 996 call dzero(work(blajkl),ntklaj) 997 call cc_r12xsort(work(krajkl),work(blajkl),nvir0) 998 999 if (froimp) then 1000 call dzero(work(blajklf),ntklaj) 1001 idx = 0 1002 do isymkl = 1, nsym 1003 isymaj = isymkl 1004 isymij = isymkl 1005 do isyml =1, nsym 1006 isymk = muld2h(isymkl,isyml) 1007 do k = nrhffr(isymk)+1, nrhfs(isymk) 1008 do l = nrhffr(isyml)+1, nrhfs(isyml) 1009 do isyma =1, nsym 1010 isymj = muld2h(isymaj,isyma) 1011 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 1012 do j = nrhffr(isymj)+1, nrhfs(isymj) 1013 do a = 1, nvir0(isyma) 1014 idx = idx + 1 1015 kl =imatmn(isymk,isyml)+ 1016 & (k-1)*nrhfs(isyml) + l 1017 aj = imataj(isymj,isyma)+ 1018 & (j-1)*nvir0(isyma) + a 1019 ajkl = imatkl1(isymaj,isymkl)+ 1020 & aj +(kl-1)*nmatan(isymaj) 1021 work(blajklf+idx-1) = work(blajkl+ajkl-1) 1022 enddo 1023 enddo 1024 enddo 1025 enddo 1026 enddo 1027 enddo 1028 enddo 1029 endif 1030 1031 1032 if (froimp) then 1033 call daxpy(ntklaj,one,work(blajklf),1,work(kxajkl),1) 1034 else 1035 call daxpy(ntklaj,one,work(blajkl),1,work(kxajkl),1) 1036 endif 1037 1038 call dzero(work(klijan),ncvai) 1039 do isymkl = 1, nsym 1040 isyman = isymkl 1041 isymij = isymkl 1042 do isyml =1, nsym 1043 isymk = muld2h(isymkl,isyml) 1044 do k = 1, nrhf(isymk) 1045 do l = 1, nrhf(isyml) 1046 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 1047 do isyma =1, nsym 1048 isymn = muld2h(isyman,isyma) 1049 isymi = muld2h(isymij,isymn) 1050 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 1051 ntota = max(1,nvir0(isyma)) 1052 ntotn = max(1,nrhf(isymn)) 1053 ntoti = max(1,nrhf(isymi)) 1054 koffd = kdotil+ igamsq(isymij,isymkl) 1055 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymn) 1056 klijanz = klijan + imatan(isyma,isymi) 1057 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 1058 & nrhf(isymn),one,work(kxajkl), 1059 & ntota,work(koffd),ntoti,one, 1060 & work(klijanz),ntota) 1061 kxajkl = kxajkl + nvir0(isyma)*nrhf(isymn) 1062 blajklf = blajklf + nvir0(isyma)*nrhf(isymn) 1063 1064 enddo 1065 enddo 1066 enddo 1067 enddo 1068 enddo 1069 1070 1071 1072 call r12gradaxd(ipdd,work(kccr12),work(kctil),work(klij), 1073 & work(kend2), 1074 & v2am,res,lwrk2) 1075 call daxpy(ncvai,one,work(klijan),1,res(1),1) 1076 1077 1078 call qexit('r12gradap') 1079 end 1080c-------------------------------------------------------------------------- 1081C=====================================================================* 1082C /* Deck r12gradaxd */ 1083 subroutine r12gradaxd(ipdd,kccr12,kctil,dklij,work,v2am,res, 1084 & lwork) 1085 implicit none 1086#include "ccsdsym.h" 1087#include "ccorb.h" 1088#include "dummy.h" 1089#include "r12int.h" 1090#include "priunit.h" 1091#include "ccfro.h" 1092 character*8 CCR12YIJ 1093 integer krsing,krtrip,kend1,lwrk1,nrhftria,n2 1094 integer lwork,lu43,idxkl,isymv,isymi,isymj,isymk,idxij 1095 integer isyml,isymkl,isymij,idxji,idxklij,idxijkl,idxklji 1096 integer isymia,luvajkl,ntklaj,isymaj,kend2 1097 integer nmatan(8),isyman,isyma,isymn,ntotij,ntotan,nvir0(8) 1098 integer lctil,klijan,ntota,ntoti,ntotn,lcvai,isymai 1099 integer kijan,ntotk,koffc,lccr12,koffd 1100 integer idxijlk,idxlk,isymin,kbajkl,ncvai,imatan(8,8),icou3 1101 integer lwrk2,lctil1,kcotil,kdotil,koccr,koctil 1102 1103 integer kxijkl,kaklij,nak,isymik,nai,naiki,nkj,naikj,isymjk 1104 integer ntotmu,dklijz,eklij,eklijz,lxxr12,kxxr12,kxsing,kxtrip 1105 integer kiakj,nia,njk,niajk,idx,ntotj,kgiakj,kgaijk,kaxdia 1106 integer fklij,fklijz,kxcijkl,kxcijklz,ntklij,ntaj 1107 integer naj,nki,najki,nji,naijk,nakji,kxcijklt,nka 1108 integer rklij,koffh,dijkl,hklij,kaxcia,kaxdai,kiajk,kgaikj 1109 integer dlkmnz,dlkmn,kctil1,kcxijklt,kcxijkl,krcc12,lunit 1110 integer kcxklmnz,kcxmnklz,kcxmnkl,kcxklmn,kxxr12n,idxjikl,idxjilk 1111 integer aijkl,hijkl,hklijz,lxxr12n,dijklz,eijkl,eijklz 1112 integer njaik,nik,nja,najik,nikja,aklij,kctil2,aijklz,lrcc12,bijkl 1113 integer isymja,isymka,ntoto,isymoj,isymo,idxlkij,kres,nikaj 1114 integer kiakjz,aj,ajkl,kl,isymkj,kiajkz,rklijz,index 1115 integer ipdd,resb,kxtil,kresz,icoun4,nmataj(8),ncvij,nijkl(8) 1116#if defined (SYS_CRAY) 1117 real work(*),dklij(*),v2am(*),kctil(*),kccr12(*),one,two,zero 1118 real res(*) 1119#else 1120 double precision work(*),dklij(*),v2am(*),kctil(*),one,zero,two 1121 double precision kccr12(*),res(*) 1122#endif 1123 parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0) 1124 1125 index(i,j) = max(i,j)*(max(i,j)-3)/2 +i+j 1126 1127 call qenter('r12gradaxd') 1128 1129C Calculate some offsets 1130 1131 isymv = 1 1132 1133 nrhftria = nrhft*(nrhft+1)/2 1134 n2 = nrhftria * nrhftria 1135 1136 1137 1138 do isymaj = 1,nsym 1139 isymij = isymaj 1140 ncvij = 0 1141 do isyma = 1,nsym 1142 isymj = muld2h(isymaj,isyma) 1143 isymi = muld2h(isymij,isymj) 1144 ncvij = ncvij + nrhf(isyma)*nrhf(isymi) 1145 enddo 1146 enddo 1147 1148 1149 1150 1151 do isymaj = 1,nsym 1152 isymij = isymaj 1153 ntaj = 0 1154 do isyma = 1,nsym 1155 isymj = muld2h(isymaj,isyma) 1156 isymi = muld2h(isymij,isymj) 1157 ntaj = ntaj + (norb1(isyma) - nrhfs(isyma))*nrhf(isymi) 1158 enddo 1159 enddo 1160 1161 1162 ntklaj = 0 1163 do isymkl = 1,nsym 1164 isymaj = isymkl 1165 do isymj = 1,nsym 1166 isyma = muld2h(isymaj,isymj) 1167 ntklaj = ntklaj + nmatij(isymkl)*nrhf(isymj) 1168 & *(norb1(isyma) - nrhfs(isyma)) 1169 enddo 1170 enddo 1171 1172 1173 1174 kxxr12 = 1 1175 kxsing = kxxr12 + ngamsq(1) 1176 kxtrip = kxsing + n2 1177 eklij = kxtrip + n2 1178 kgiakj = eklij + ncvij 1179 kgaijk = kgiakj + ntklaj 1180 kxcijkl = kgaijk + ntklaj 1181 fklij = kxcijkl + ngamsq(1) 1182 kiajk = fklij + ncvij 1183 kgaikj = kiajk + ntklaj 1184 rklij = kgaikj + ntklaj 1185 kcxijkl = rklij + ncvij 1186 eijkl = kcxijkl + ngamsq(1) 1187 kxtil = eijkl + ncvij 1188 resb = kxtil + ngamsq(1) 1189 kend2 = resb + ntaj 1190 lwrk2 = lwork - kend2 1191 1192C Read matrix X 1193 1194 lu43 = -43 1195 call gpopen(lu43,'CCR12_XP','unknown',' ','formatted', 1196 & idummy,.false.) 1197 read(lu43,*) i 1198 read(lu43,'(4e30.20)') (work(kxsing+i), i = 0, n2-1 ) 1199 read(lu43,'(4e30.20)') (work(kxtrip+i), i = 0, n2-1 ) 1200 if (ipdd.eq.3) then 1201 call gpclose(lu43,'KEEP') 1202 else 1203 call gpclose(lu43,'DELETE') 1204 endif 1205 call dscal(n2,2d0,work(kxsing),1) 1206 call dscal(n2,2d0/3d0,work(kxtrip),1) 1207 1208 call cc_r12vpack(work(kxxr12),isymv,work(kxsing),work(kxtrip), 1209 & nrhf,nrhfa,nijkl) 1210 1211 1212C Calculate XC_{ij}^{kl} = X_{ij}^{mn} * Ctilde_{mn}^{kl} 1213 1214 call dzero(work(kxcijkl),ngamsq(1)) 1215 do isymkl = 1,nsym 1216 isymij = isymkl 1217 lxxr12 = kxxr12 + igamsq(isymij,isymkl) 1218 lctil = 1 + igamsq(isymij,isymkl) 1219 kxcijklz = kxcijkl + igamsq(isymij,isymkl) 1220 ntoti = max(1,nmatij(isymij)) 1221 ntotk = max(1,nmatij(isymkl)) 1222 call dgemm('N','T',nmatij(isymij),nmatij(isymkl), 1223 & nmatij(isymkl),one,work(lxxr12), 1224 & ntotk,kctil(lctil),ntoti,zero, 1225 & work(kxcijklz),ntotk) 1226 end do 1227c --------------------------------- 1228C Transpone CX_{ij}^{kl} 1229c --------------------------------- 1230 1231 call dzero(work(kcxijkl),ngamsq(1)) 1232 do isymkl = 1, nsym 1233 isymij = muld2h(isymkl,isymv) 1234 do isyml =1, nsym 1235 isymk = muld2h(isymkl,isyml) 1236 do k = 1, nrhf(isymk) 1237 do l = 1, nrhf(isyml) 1238 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 1239 do isymi = 1, nsym 1240 isymj = muld2h(isymij,isymi) 1241 do j = 1, nrhf(isymj) 1242 do i = 1, nrhf(isymi) 1243 idxij = imatij(isymi,isymj)+ 1244 & nrhf(isymi)*(j-1) + i 1245 idxklij = igamsq(isymkl,isymij)+ 1246 & nmatij(isymkl)*(idxij-1)+idxkl 1247 idxijkl = igamsq(isymij,isymkl)+ 1248 & nmatij(isymij)*(idxkl-1) +idxij 1249 work(kcxijkl-1+idxijkl) = 1250 & work(kxcijkl-1+idxklij) 1251 end do 1252 end do 1253 1254 end do 1255 end do 1256 end do 1257 end do 1258 end do 1259 1260 1261 1262 call dzero(work(eklij),ncvij) 1263 call dzero(work(eijkl),ncvij) 1264 call dzero(work(fklij),ncvij) 1265 do isymkl = 1, nsym 1266 isymoj = isymkl 1267 isymij = isymkl 1268 do isyml =1, nsym 1269 isymk = muld2h(isymkl,isyml) 1270 do k = 1, nrhf(isymk) 1271 do l = 1, nrhf(isyml) 1272 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 1273 do isymo =1, nsym 1274 isymj = muld2h(isymoj,isymo) 1275 isymi = muld2h(1,isymo) 1276 1277C Calculate X_{ol}^{ij} * D_{ij}^{kl} 1278 1279 koffd = 1 + igamsq(isymij,isymkl) 1280 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 1281 lxxr12 = kxxr12 + igamsq(isymij,isymkl) 1282 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 1283 ntotj = max(1,nrhf(isymj)) 1284 ntoti = max(1,nrhf(isymi)) 1285 ntoto = max(1,nrhf(isymo)) 1286 eklijz = eklij + imatij(isymo,isymi) 1287 call dgemm('N','T',nrhf(isymo),nrhf(isymi), 1288 & nrhf(isymj),one,work(lxxr12), 1289 & ntoto,dklij(koffd),ntoti,one, 1290 & work(eklijz),ntoto) 1291 eijklz = eijkl + imatij(isymi,isymo) 1292 call dgemm('N','T',nrhf(isymi),nrhf(isymo), 1293 & nrhf(isymj),one,dklij(koffd), 1294 & ntoti,work(lxxr12),ntoto,one, 1295 & work(eijklz),ntoti) 1296 1297C Calculate F_{}^{} = C_{ij}^{kl} * XC_{ij}^{ol} 1298 1299 lccr12 = 1+ igamsq(isymij,isymkl) 1300 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 1301 kcxijklt = kcxijkl + igamsq(isymij,isymkl) 1302 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 1303 fklijz = fklij + imatij(isymo,isymi) 1304 call dgemm('N','T',nrhf(isymi),nrhf(isymo), 1305 & nrhf(isymj),one,kccr12(lccr12), 1306 & ntoti,work(kcxijklt),ntoti,one, 1307 & work(fklijz),ntoti) 1308 1309 enddo 1310 enddo 1311 enddo 1312 enddo 1313 end do 1314 1315 call dzero(work(rklij),ncvij) 1316 call daxpy(ncvij,one,work(eklij),1,work(rklij),1) 1317 call daxpy(ncvij,one,work(eijkl),1,work(rklij),1) 1318 call daxpy(ncvij,-two,work(fklij),1,work(rklij),1) 1319 call dscal(ncvij,two,work(rklij),1) 1320 lunit = -1 1321 if (ipdd .eq. 3) then 1322 call gpopen(lunit,'CCR12YIJ','unknown',' ','unformatted', 1323 & idummy,.false.) 1324 else if (ipdd .eq. 5) then 1325 call gpopen(lunit,'CCR12YIJB','unknown',' ','unformatted', 1326 & idummy,.false.) 1327 endif 1328 write(lunit)(work(rklij+i-1),i=1,ncvij) 1329 call gpclose(lunit,'KEEP') 1330 1331 1332 call dzero(work(kgiakj),ntklaj) 1333 call dzero(work(kgaijk),ntklaj) 1334 call dzero(work(kgaikj),ntklaj) 1335 call dzero(work(kiajk),ntklaj) 1336 1337 idx = 0 1338 do isymka = 1,nsym 1339 isymij = isymka 1340 do isyma = 1, nsym 1341 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 1342 isymk = muld2h(isymka,isyma) 1343 do isymj = 1,nsym 1344 isymi = muld2h(isymij,isymj) 1345 do i = 1,nrhf(isymi) 1346 do j = 1,nrhf(isymj) 1347 isymia = muld2h(isymi,isyma) 1348 isymkj = muld2h(isymk,isymj) 1349 isymja = muld2h(isymj,isyma) 1350 isymik = muld2h(isymi,isymk) 1351 do k = 1, nrhf(isymk) 1352 do a = 1,nvir0(isyma) 1353 idx = idx + 1 1354 nai = it1am(isyma,isymi) 1355 * + nvir(isyma)*(i-1) 1356 & + a + nrhfs(isyma) 1357 njk = it1am(isymj,isymk) 1358 * + nvir(isymj)*(k-1) + j+nrhffr(isymj) 1359 naijk = it2am(isymia,isymkj) 1360 * + index(nai,njk) 1361 nak = it1am(isyma,isymk) 1362 * + nvir(isyma)*(k-1) 1363 & + a + nrhfs(isyma) 1364 nji = it1am(isymj,isymi) 1365 * + nvir(isymj)*(i-1) + j+nrhffr(isymj) 1366 nakji = it2am(isymka,isymij) 1367 * + index(nak,nji) 1368 naj = it1am(isyma,isymj) 1369 * + nvir(isyma)*(j-1) 1370 & + a + nrhfs(isyma) 1371 nki = it1am(isymk,isymi) 1372 * + nvir(isymk)*(i-1) + k +nrhffr(isymk) 1373 nik = it1am(isymi,isymk) 1374 * + nvir(isymi)*(k-1) + i 1375 najki = it2am(isymja,isymik) 1376 * + index(naj,nki) 1377 work(kgiakj+idx-1) = v2am(naijk) 1378 work(kgaijk+idx-1) = v2am(nakji) 1379 work(kgaikj+idx-1) = v2am(najki) 1380 work(kiajk+idx-1) = 4.0d0*work(kgaijk+idx-1) 1381 & - work(kgaikj+idx-1) - work(kgiakj+idx-1) 1382 enddo 1383 enddo 1384 enddo 1385 enddo 1386 enddo 1387 enddo 1388 enddo 1389 1390 1391 1392 do isymi = 1,nsym 1393 nvir0(isymi) = norb1(isymi)-nrhfs(isymi) 1394 nmataj(isymi) = nvir0(isymi)*nrhf(isymi) 1395 enddo 1396 1397 1398 1399 call dzero(res(1),ntaj) 1400 kres = 1 1401 do isymij = 1,nsym 1402 isymkl = isymij 1403 isymaj = isymij 1404 do isyma =1, nsym 1405 isymj = muld2h(isymaj,isyma) 1406 isymi = muld2h(1,isyma) 1407 nvir0(isyma) = norb1(isyma)-nrhfs(isyma) 1408 ntotk = max(1,nvir0(isyma)*nrhf(isymj)) 1409 ntoti = max(1,nmatij(isymkl)) 1410 call dgemm('N','N',nrhf(isymj)*nvir0(isyma),1, 1411 & nmatij(isymkl),one,work(kiajk),ntotk, 1412 & work(rklij),ntoti,zero,res(kres),ntotk) 1413 kres = kres + nvir0(isyma)*nrhf(isymj) 1414 kiajk = kiajk + nrhf(isymj)*nvir0(isyma) 1415 & * nmatij(isymkl) 1416 enddo 1417 enddo 1418 1419 call dscal(ntaj,0.5d0,res(1),1) 1420 1421 1422 if (ipdd .eq. 5) then 1423 call dscal(ntaj,1.0d0,res(1),1) 1424 call dzero(work(resb),ntaj) 1425 call r12gradb(work(resb),dklij,work(kend2),ntklaj,ntaj,lwrk2) 1426 1427 call daxpy(ntaj,one,work(resb),1,res(1),1) 1428 endif 1429 call qexit('r12gradaxd') 1430 end 1431c-------------------------------------------------------------------------- 1432C=====================================================================* 1433C /* Deck r12gradb */ 1434 subroutine r12gradb(res,dklij,work,ntklaj,ncvai,lwork) 1435 implicit none 1436#include "ccsdsym.h" 1437#include "ccorb.h" 1438#include "dummy.h" 1439#include "r12int.h" 1440#include "priunit.h" 1441#include "ccfro.h" 1442 logical locdbg 1443 parameter (locdbg = .false.) 1444 integer kqajkl,kqijal,kuijal,krajkl,kend1 1445 integer isym,luvajkl,nvir0t,lwork,lwrk1,ntklaj 1446 integer icou3,imataj(8,8),koffd,ksajkl,nvir0(8),ncvai 1447 integer kcqia,kcqai,lcqia,lcqai,isymai,isyma,isymi 1448 integer isymaj,isymkl,isymij,isyml,isymk,isymj,idxkl,ntotj 1449 integer isymv,idxlk,idxklij,idxlkij,idxij,idxji 1450 integer idxjikl,idxijkl,idxijlk,koff,ntota,ntoti 1451 integer kuajkl,ktajkl,ktijal,kcuia,kcuai,lcuai,lcuia 1452 integer kssijal,ksrajkl,kstijal,idxklji,koff1,kdjikl 1453 integer kuijals,kqijals,kssajkl 1454 1455 integer kccr12,krsing,krtrip,kctil,lccr12,lctil 1456 integer nrhftria,n2,lu43,ntotk 1457#if defined (SYS_CRAY) 1458 real dklij(*),res(*),work(*),one,two,zero 1459#else 1460 double precision dklij(*),res(*),work(*),one,zero,two 1461#endif 1462 parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0) 1463 1464 call qenter('r12gradb') 1465 1466 isymv = 1 1467 1468 nrhftria = nrhft*(nrhft+1)/2 1469 n2 = nrhftria * nrhftria 1470 1471 1472 nvir0t = 0 1473 do isym = 1, nsym 1474 nvir0(isym) = norb1(isym) - nrhfs(isym) 1475 nvir0t = nvir0t + nvir0(isym) 1476 end do 1477 1478 icou3 = 0 1479 do isymai = 1,nsym 1480 icou3 = 0 1481 do isyma = 1,nsym 1482 isymi = muld2h(isymai,isyma) 1483 imataj(isyma,isymi)= icou3 1484 icou3 = icou3+nrhf(isymi)*(norb1(isyma) - nrhfs(isyma)) 1485 enddo 1486 enddo 1487 1488 1489 kqajkl = 1 1490 kqijal = kqajkl + ntklaj 1491 kuijal = kqijal + ntklaj 1492 krajkl = kuijal + ntklaj 1493 ksajkl = krajkl + nvir0t*nrhft**3 1494 kcqai = ksajkl + nvir0t*nrhft**3 1495 kcqia = kcqai + ncvai 1496 ktijal = kcqia + ncvai 1497 kcuai = ktijal + nvir0t*nrhft**3 1498 kuajkl = kcuai + ncvai 1499 kcuia = kuajkl + ntklaj 1500 kstijal= kcuia + ncvai 1501 ksrajkl= kstijal+ ntklaj 1502 kssijal= ksrajkl+ ntklaj 1503 kssajkl= kssijal+ ntklaj 1504 kuijals= kssajkl+ ntklaj 1505 kqijals= kuijals+ ntklaj 1506 kend1 = kqijals+ ntklaj 1507 lwrk1 = lwork - kend1 1508 1509 1510 luvajkl = -1 1511 call dzero(work(kqajkl),ntklaj) 1512 call gpopen(luvajkl,'DDR12QAJKL','unknown',' ','unformatted', 1513 & idummy,.false.) 1514 rewind(luvajkl) 1515 read(luvajkl) (work(kqajkl+i-1), i = 1,ntklaj) 1516 call gpclose(luvajkl,'DELETE') 1517 1518 1519 luvajkl = -1 1520 call dzero(work(kqijal),ntklaj) 1521 call gpopen(luvajkl,'DDR12QIJAL','unknown',' ','unformatted', 1522 & idummy,.false.) 1523 rewind(luvajkl) 1524 read(luvajkl) (work(kqijal+i-1), i = 1,ntklaj) 1525 call gpclose(luvajkl,'DELETE') 1526 1527 call dzero(work(kqijals),ntklaj) 1528 1529 luvajkl = -1 1530 call dzero(work(kuijal),ntklaj) 1531 call gpopen(luvajkl,'DDR12UIJAL','unknown',' ','unformatted', 1532 & idummy,.false.) 1533 rewind(luvajkl) 1534 read(luvajkl) (work(kuijal+i-1), i = 1,ntklaj) 1535 call gpclose(luvajkl,'DELETE') 1536 1537 call dzero(work(kuijals),ntklaj) 1538 1539 1540 luvajkl = -1 1541 call dzero(work(kuajkl),ntklaj) 1542 call gpopen(luvajkl,'DDR12UAJKL','unknown',' ','unformatted', 1543 & idummy,.false.) 1544 rewind(luvajkl) 1545 read(luvajkl) (work(kuajkl+i-1), i = 1,ntklaj) 1546 call gpclose(luvajkl,'DELETE') 1547 1548 call dzero(work(krajkl),nvir0t*nrhft**3) 1549 luvajkl = -1 1550 call gpopen(luvajkl,'AUXQSA12','unknown',' ','formatted', 1551 & idummy,.false.) 1552 rewind(luvajkl) 1553 read(luvajkl,'(4E30.20)') (work(krajkl+i-1), i=1,nvir0t*nrhft**3) 1554 call gpclose(luvajkl,'DELETE') 1555 1556 call dzero(work(ksrajkl),ntklaj) 1557 call cc_r12xsort2(work(krajkl),work(ksrajkl),nvir0) 1558 1559 1560 call dzero(work(ksajkl),nvir0t*nrhft**3) 1561 luvajkl = -1 1562 call gpopen(luvajkl,'AUXQSAJ12','unknown',' ','formatted', 1563 & idummy,.false.) 1564 rewind(luvajkl) 1565 read(luvajkl,'(4E30.20)') (work(ksajkl+i-1), i=1,nvir0t*nrhft**3) 1566 call gpclose(luvajkl,'DELETE') 1567 1568 call dzero(work(kssijal),ntklaj) 1569 call cc_r12xsort2(work(ksajkl),work(kssijal),nvir0) 1570 1571 call dzero(work(kssajkl),ntklaj) 1572 call cc_r12xsort3(work(ksajkl),work(kssajkl),nvir0) 1573 1574 1575 1576 call dzero(work(ktijal),nvir0t*nrhft**3) 1577 luvajkl = -1 1578 call gpopen(luvajkl,'AUXQSAI12','unknown',' ','formatted', 1579 & idummy,.false.) 1580 rewind(luvajkl) 1581 read(luvajkl,'(4E30.20)') (work(ktijal+i-1), i=1,nvir0t*nrhft**3) 1582 call gpclose(luvajkl,'DELETE') 1583 1584 call dzero(work(kstijal),ntklaj) 1585 call cc_r12xsort3(work(ktijal),work(kstijal),nvir0) 1586 1587 1588 call daxpy(ntklaj,one,work(kssijal),1,work(kqijal),1) 1589 call daxpy(ntklaj,one,work(kssajkl),1,work(kuijal),1) 1590 call daxpy(ntklaj,one,work(kstijal),1,work(kuajkl),1) 1591 call daxpy(ntklaj,one,work(ksrajkl),1,work(kqajkl),1) 1592 1593 1594 1595 call dzero(work(kcqai),ncvai) 1596 lcqai = kcqai 1597 call dzero(work(kcqia),ncvai) 1598 lcqia = kcqia 1599 call dzero(work(kcuai),ncvai) 1600 lcuai = kcuai 1601 call dzero(work(kcuia),ncvai) 1602 lcuia = kcuia 1603 1604c Calculate Q_{aj}^{kl}*D_{kl}^{ij} 1605 1606 1607 do isymkl = 1, nsym 1608 isymaj = isymkl 1609 isymij = isymkl 1610 do isyml =1, nsym 1611 isymk = muld2h(isymkl,isyml) 1612 do k = 1, nrhf(isymk) 1613 do l = 1, nrhf(isyml) 1614 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 1615 do isyma =1, nsym 1616 isymj = muld2h(isymaj,isyma) 1617 isymi = muld2h(isymij,isymj) 1618 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 1619 koffd = 1 + igamsq(isymij,isymkl) 1620 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 1621 ntota = max(1,nvir0(isyma)) 1622 ntotj = max(1,nrhf(isymj)) 1623 ntoti = max(1,nrhf(isymi)) 1624 lcqai = kcqai + imataj(isyma,isymi) 1625 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 1626 & nrhf(isymj),one,work(kqajkl), 1627 & ntota,dklij(koffd),ntoti,one, 1628 & work(lcqai),ntota) 1629 1630 kqajkl = kqajkl + nvir0(isyma)*nrhf(isymj) 1631 ksrajkl = ksrajkl + nvir0(isyma)*nrhf(isymj) 1632 1633 ntota = max(1,nvir0(isyma)) 1634 ntotj = max(1,nrhf(isymj)) 1635 ntoti = max(1,nrhf(isymi)) 1636 lcqia = kcqia + imataj(isyma,isymi) 1637 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 1638 & nrhf(isymj),one,work(kqijal), 1639 & ntota,dklij(koffd),ntoti,one, 1640 & work(lcqia),ntota) 1641 1642 kqijal = kqijal + nvir0(isyma)*nrhf(isymj) 1643 1644 lcuai = kcuai + imataj(isyma,isymi) 1645 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 1646 & nrhf(isymj),one,work(kuijal), 1647 & ntota,dklij(koffd),ntoti,one, 1648 & work(lcuai),ntota) 1649 1650 kuijal = kuijal + nvir0(isyma)*nrhf(isymj) 1651 kssijal = kssijal + nvir0(isyma)*nrhf(isymj) 1652 kssajkl = kssajkl + nvir0(isyma)*nrhf(isymj) 1653 1654 lcuia = kcuia + imataj(isyma,isymi) 1655 call dgemm('N','T',nvir0(isyma),nrhf(isymi), 1656 & nrhf(isymj),one,work(kuajkl), 1657 & ntota,dklij(koffd),ntoti,one, 1658 & work(lcuia),ntota) 1659 1660 kuajkl = kuajkl + nvir0(isyma)*nrhf(isymj) 1661 kstijal = kstijal + nvir0(isyma)*nrhf(isymj) 1662 ktijal = ktijal + nvir0(isyma)*nrhf(isymj) 1663 1664 enddo 1665 enddo 1666 enddo 1667 enddo 1668 enddo 1669 1670 1671 call r12gradbd(res,dklij,work(kend1),lwrk1) 1672 1673 call daxpy(ncvai,one,work(kcqai),1,res,1) 1674 call daxpy(ncvai,one,work(kcuai),1,res,1) 1675 call daxpy(ncvai,one,work(kcqia),1,res,1) 1676 call daxpy(ncvai,one,work(kcuia),1,res,1) 1677 1678 if (locdbg) then 1679 write(lupri,*) 'Ergebnisse r12gradb' 1680 do isymaj = 1,nsym 1681 do isyma = 1,nsym 1682 isymi =isyma 1683 write(lupri,*) 'symmetry block',isyma 1684 call output(res(1+imataj(isyma,isymi)), 1685 & 1,nvir0(isyma),1,nrhf(isymi), 1686 & nvir0(isyma),nrhf(isymi),1,lupri) 1687 enddo 1688 enddo 1689 endif 1690 1691 1692 call qexit('r12gradb') 1693 end 1694c------------------------------------------------------------------------ 1695C=====================================================================* 1696C /* Deck r12gradbd */ 1697 subroutine r12gradbd(res,dklij,work,lwork) 1698 implicit none 1699#include "ccsdsym.h" 1700#include "ccorb.h" 1701#include "dummy.h" 1702#include "r12int.h" 1703#include "priunit.h" 1704#include "ccfro.h" 1705#include "ccsdinp.h" 1706 integer lxdmmu,nr1bast,isymak,isymk,isyma,nr1bas(8) 1707 integer lwork,kxmujkl,kxdmmu,kend1,nmum,luvajkl 1708 integer isymkl,isymij,isymi,isymj,idxkl,koffd,isymaj 1709 integer isymai,icou3,imatmuj(8,8),ntota,ntoti,ntotj,isyml 1710 integer smujkl,koffc,lusifc,lwrk1,kcmo,ntota1,kxdnumu 1711 integer kxpjkl,ntklpj,isymmua,iglmrv(8,8),isymmu 1712 integer kxajkl,ntotmu,isymmuj,kxpjklz,icount2,icou2,icount7 1713 integer nr1xbas(8),koff,idxji,idxij,kzmunu,kxdpmu,ir1bas(8,8) 1714 integer nrhftria,kxsing,kxtrip,n2,isymv,lu43,nxajkl(8),kcmoa 1715 integer isym,noffset,kspjkl,nvir0t,nvir0(8),ntklaj,krajkl 1716 integer imbas1(8),kdens,lxdnumu,icou4,imatmunu(8,8),lxdpmu 1717 integer nrhfst,ksqijklz,ksqjkl,ksajkl,sajkl,icoun1,norb1t 1718 integer imataj(8,8),idx,kl,aj,ajkl,kspjklf,ksajklf 1719 integer imatmn(8,8),icou6,nmatkl1(8),imatkl1(8,8),icou5 1720 integer nmatan(8),icou7 1721 integer imatoi(8,8),icou8,isyman,isymn 1722 integer iglmrv1(8,8),noffset1,icou9,nrhffrt 1723 integer ksijklf,ij,ijkl,ksijkl,sijklf 1724#if defined (SYS_CRAY) 1725 real res(*),work(*),dklij(*),zero,one,two 1726#else 1727 double precision res(*),work(*),dklij(*),zero,one,two 1728#endif 1729 parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0) 1730 1731 call qenter('r12gradbd') 1732 1733 isymv = 1 1734 1735 nrhftria = nrhft*(nrhft+1)/2 1736 n2 = nrhftria * nrhftria 1737 1738 1739 nvir0t = 0 1740 norb1t = 0 1741 nrhfst = 0 1742 nrhffrt = 0 1743 do isymi = 1, nsym 1744 nvir0(isymi) = norb1(isymi) - nrhfs(isymi) 1745 nvir0t = nvir0t + nvir0(isymi) 1746 norb1t = norb1t + norb1(isymi) 1747 nrhfst = nrhfst + nrhfs(isymi) 1748 nrhffrt = nrhffrt+ nrhffr(isymi) 1749 end do 1750 1751 1752 do isymmua = 1,nsym 1753 icou2 = 0 1754 icou3 = 0 1755 do isyma = 1,nsym 1756 isymmu = muld2h(isymmua,isyma) 1757 iglmrv(isymmu,isyma) = icou2 1758 iglmrv1(isymmu,isyma) = icou3 1759 icou2 = icou2 + nbas(isymmu)*(nvir(isyma)) 1760 icou3 = icou3 + nbas(isymmu)*(nvir(isyma)-nrhffr(isyma)) 1761 enddo 1762 enddo 1763 1764 do isyman = 1,nsym 1765 icou3 = 0 1766 icou4 = 0 1767 icou7 = 0 1768 nmatkl1(isyman) = 0 1769 do isyma = 1,nsym 1770 isymn = muld2h(isyman,isyma) 1771 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 1772 icou3 = icou3 + nvir0(isyma)*nrhfs(isymn) 1773 icou4 = icou4 + nrhfs(isyma)*nrhfs(isymn) 1774 enddo 1775 nmatan(isyman) = icou3 1776 nmatkl1(isyman) = icou4 1777 enddo 1778 1779 do isyman = 1,nsym 1780 icou4 = 0 1781 icou5 = 0 1782 icou6 = 0 1783 icou8 = 0 1784 do isyma = 1,nsym 1785 isymn = muld2h(isyman,isyma) 1786 imataj(isyma,isymn)= icou4 1787 imatkl1(isyma,isymn)= icou5 1788 imatmn(isyma,isymn)= icou6 1789 imatoi(isyma,isymn)= icou8 1790 icou4 = icou4+nrhfs(isymn)*(norb1(isyma) - nrhfs(isyma)) 1791 icou5 = icou5+nmatan(isymn)* nmatkl1(isyma) 1792 icou6 = icou6+nrhfs(isyma)*nrhfs(isymn) 1793 icou8 = icou8+nmatkl1(isyma)*nmatkl1(isymn) 1794 enddo 1795 enddo 1796 1797 1798 1799 do isymak = 1, nsym 1800 nr1bas(isymak) = 0 1801 nr1xbas(isymak) = 0 1802 do isymk = 1, nsym 1803 isyma = muld2h(isymak,isymk) 1804 nr1bas(isymak) = nr1bas(isymak) +mbas1(isyma)*nrhf(isymk) 1805 nr1xbas(isymak) = nr1xbas(isymak)+mbas2(isyma)*nrhf(isymk) 1806 end do 1807 end do 1808 1809 nr1bast = 0 1810 do isyma = 1,nsym 1811 nr1bast = nr1bast + nr1bas(isyma) 1812 enddo 1813 1814 ntklpj = 0 1815 do isymkl = 1,nsym 1816 isymaj = isymkl 1817 do isymj = 1,nsym 1818 isyma = muld2h(isymaj,isymj) 1819 ntklpj = ntklpj + nmatij(isymkl)*nrhf(isymj)*norb1(isyma) 1820 enddo 1821 enddo 1822 1823 1824 1825 do isymak = 1, nsym 1826 nvajkl(isymak) = 0 1827 icount7 = 0 1828 icount2 = 0 1829 do isymk = 1, nsym 1830 isyma = muld2h(isymak,isymk) 1831 ivajkl(isyma,isymk) = icount7 1832 nvajkl(isymak) = nvajkl(isymak)+ nr1bas(isyma)*nmatij(isymk) 1833 icount7 = icount7 + nr1bas(isyma)*nmatij(isymk) 1834 ir1bas(isyma,isymk) = icount2 1835 icount2 = icount2 + nrhf(isymk)*mbas1(isyma) 1836 end do 1837 end do 1838 1839 icou3 = 0 1840 do isymai = 1,nsym 1841 icou3 = 0 1842 icou4 = 0 1843 do isyma = 1,nsym 1844 isymi = muld2h(isyma,isymai) 1845 imatmuj(isyma,isymi)= icou3 1846 icou3 = icou3+nrhf(isymi)*norb1(isyma) 1847 imatmunu(isymi,isyma) = icou4 1848 icou4 = icou4+mbas1(isymi)*mbas1(isyma) 1849 enddo 1850 enddo 1851 1852 ntklaj = 0 1853 do isymkl = 1, nsym 1854 isymaj = isymkl 1855 do isyml =1, nsym 1856 isymk = muld2h(isymkl,isyml) 1857 do k = 1, nrhf(isymk) 1858 do l = 1, nrhf(isyml) 1859 do isyma =1, nsym 1860 isymj = muld2h(isymaj,isyma) 1861 isymi = isyma 1862 ntklaj = ntklaj + nrhf(isymj) * 1863 & (norb1(isyma) - nrhfs(isyma)) 1864 enddo 1865 enddo 1866 enddo 1867 enddo 1868 enddo 1869 1870 1871 1872 kxdmmu = 1 1873 kxdnumu = kxdmmu + nr1bast 1874 kcmo = kxdnumu + 10*mbas1t*mbas1t 1875 smujkl = kcmo + nlamds 1876 kxpjkl = smujkl + nrhfst**4 1877 kxdpmu = kxpjkl + ntklpj 1878 kzmunu = kxdpmu + mbas1t*nrhft 1879 kcmoa = kzmunu + n2basx 1880 kxmujkl = kcmoa + nlamds 1881 kspjkl = kxmujkl + nvajkl(1) 1882 kdens = kspjkl + ntklpj 1883 sajkl = kdens + n2basx 1884 ksajkl = sajkl + nvir0t*nrhfst**3 1885 ksqjkl = ksajkl + nvir0t*nrhfst**3 1886 ksajklf = ksqjkl + ntklpj 1887 kspjklf = ksajklf + nvir0t*nrhfst**3 1888 ksijklf = kspjklf + nrhfst**4 1889 ksijkl = ksijklf + nrhfst*nrhfst**3 1890 kend1 = ksijkl + nrhfst*nrhfst**3 1891 lwrk1 = lwork - kend1 1892 1893 1894 call dzero(work(kxmujkl),nvajkl(1)) 1895 luvajkl = -1 1896 call gpopen(luvajkl,'CCR12QAJKL','unknown',' ','unformatted', 1897 & idummy,.false.) 1898 rewind(luvajkl) 1899 read(luvajkl) (work(kxmujkl+i-1), i = 1,nvajkl(1)) 1900 call gpclose(luvajkl,'KEEP') 1901c 1902 call dzero(work(smujkl),nrhfst**4) 1903 luvajkl = -1 1904 1905 call gpopen(luvajkl,'AUXQ12','unknown',' ','formatted', 1906 & idummy,.false.) 1907 rewind(luvajkl) 1908 read(luvajkl,'(4E30.20)') (work(smujkl+i-1), i = 1, 1909 & nrhfst**4) 1910 call gpclose(luvajkl,'KEEP') 1911c 1912 luvajkl = -1 1913 call gpopen(luvajkl,'AUXQA12','unknown',' ','formatted', 1914 & idummy,.false.) 1915 rewind(luvajkl) 1916 read(luvajkl,'(4E30.20)') (work(sajkl+i-1), i = 1, 1917 & nvir0t*nrhfst**3) 1918 call gpclose(luvajkl,'DELETE') 1919 1920 1921 1922 call dzero(work(kspjkl),ntklpj) 1923 call dzero(work(ksijkl),nrhfst**4) 1924 call cc_r12xsort1(work(smujkl),work(kspjkl)) 1925 call cc_r12xsort(work(sajkl),work(ksajkl),nvir0) 1926 call cc_r12xsort4(work(smujkl),work(ksijkl)) 1927 1928 if (froimp) then 1929 call dzero(work(ksijklf),nrhfst**4) 1930 idx = 0 1931 do isymkl = 1, nsym 1932 isymij = isymkl 1933 do isyml =1, nsym 1934 isymk = muld2h(isymkl,isyml) 1935 do k = nrhffr(isymk)+1, nrhfs(isymk) 1936 do l = nrhffr(isyml)+1, nrhfs(isyml) 1937 idxkl=imatij(isymk,isyml)+nrhfs(isymk)*(l-1)+k 1938 do isymi =1, nsym 1939 isymj = muld2h(isymij,isymi) 1940 do j = nrhffr(isymj)+1, nrhfs(isymj) 1941 do i = 1, nrhffr(isymi) 1942 idx = idx + 1 1943 kl = imatmn(isymk,isyml)+ 1944 & (l-1)*nrhfs(isymk) + k 1945 ij = imatmn(isymi,isymj)+ 1946 & (j-1)*nrhfs(isymi) + i 1947 ijkl = imatoi(isymij,isymkl)+ 1948 & (kl-1)*nmatkl1(isymij) + ij 1949 work(ksijklf+idx-1) = work(ksijkl+ijkl-1) 1950 enddo 1951 enddo 1952 enddo 1953 enddo 1954 enddo 1955 enddo 1956 enddo 1957 1958 1959 1960 call dzero(work(ksajklf),ntklaj) 1961 idx = 0 1962 do isymkl = 1, nsym 1963 isymaj = isymkl 1964 isymij = isymkl 1965 do isyml =1, nsym 1966 isymk = muld2h(isymkl,isyml) 1967 do k = nrhffr(isymk)+1, nrhfs(isymk) 1968 do l = nrhffr(isyml)+1, nrhfs(isyml) 1969 idxkl=imatij(isymk,isyml)+nrhfs(isymk)*(l-1)+k 1970 do isyma =1, nsym 1971 isymj = muld2h(isymaj,isyma) 1972 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 1973 do j = nrhffr(isymj)+1, nrhfs(isymj) 1974 do a = 1, nvir0(isyma) 1975 idx = idx + 1 1976 kl =imatmn(isyml,isymk)+ 1977 & (k-1)*nrhfs(isyml) + l 1978 aj = imataj(isyma,isymj)+ 1979 & (j-1)*nvir0(isyma) + a 1980 ajkl = imatkl1(isymaj,isymkl)+ 1981 & aj +(kl-1)*nmatan(isymaj) 1982 work(ksajklf+idx-1) = work(ksajkl+ajkl-1) 1983 enddo 1984 enddo 1985 enddo 1986 enddo 1987 enddo 1988 enddo 1989 enddo 1990 1991 call dzero(work(kspjklf),nrhfst**4) 1992 idx = 0 1993 do isymkl = 1, nsym 1994 isymaj = isymkl 1995 isymij = isymkl 1996 do isyml =1, nsym 1997 isymk = muld2h(isymkl,isyml) 1998 do k = nrhffr(isymk)+1, nrhfs(isymk) 1999 do l = nrhffr(isyml)+1, nrhfs(isyml) 2000 do isyma =1, nsym 2001 isymj = muld2h(isymaj,isyma) 2002 nvir0(isyma) = norb1(isyma) - nrhfs(isyma) 2003 do j = nrhffr(isymj)+1, nrhfs(isymj) 2004 do a = nrhffr(isyma)+1, nrhfs(isyma) 2005 idx = idx + 1 2006 kl =imatmn(isyml,isymk)+ 2007 & (k-1)*nrhfs(isyml) + l 2008 aj = imatmn(isyma,isymj)+ 2009 & (j-1)*nrhfs(isyma) + a 2010 ajkl = imatoi(isymaj,isymkl)+ 2011 & aj +(kl-1)*nmatkl1(isymaj) 2012 work(kspjklf+idx-1) = work(kspjkl+ajkl-1) 2013 enddo 2014 enddo 2015 enddo 2016 enddo 2017 enddo 2018 enddo 2019 enddo 2020 2021 endif 2022 2023c Read MO-coefficients 2024 2025 call dzero(work(kcmo),nlamds) 2026 lusifc = -1 2027 call gpopen(lusifc,'SIRIFC','OLD',' ','UNFORMATTED', 2028 & idummy,.false.) 2029 rewind(lusifc) 2030 call mollab(label,lusifc,lupri) 2031 read(lusifc) 2032 read(lusifc) 2033 read(lusifc) (work(kcmo+i-1),i=1,nlamds) 2034 call gpclose(lusifc,'KEEP') 2035 2036! reorder to occupied, symmetry; virtual, symmetry 2037! from symmetry, (occupied,virtual) 2038 2039 call cmo_reorder(work(kcmo),work(kend1),lwrk1) 2040 2041 noffset = 0 2042 noffset1 = 0 2043 do isym = 1,nsym 2044 noffset = noffset + nbas(isym)*nrhf(isym) 2045 noffset1 = noffset1 + nbas(isym)*nrhfs(isym) 2046 enddo 2047 2048c Transform into orbital basis 2049 2050 call dzero(work(kxpjkl),ntklpj) 2051 kxpjklz = kxpjkl 2052 do isymkl = 1, nsym 2053 isymaj = isymkl 2054 isymmuj = isymkl 2055 do isyml =1, nsym 2056 isymk = muld2h(isymkl,isyml) 2057 do k = 1, nrhf(isymk) 2058 do l = 1, nrhf(isyml) 2059 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 2060 do isymmu =1, nsym 2061 isyma = muld2h(1,isymmu) 2062 isymj = muld2h(isymmuj,isymmu) 2063 koffc = kcmo + noffset+ 2064 & iglmrv(isymmu,isyma) 2065 kxajkl = kxmujkl+ivajkl(isymmuj,isymkl)+ 2066 & nr1bas(isymmuj)*(idxkl-1) 2067 & +ir1bas(isymmu,isymj) 2068 ntotmu = max(1,mbas1(isymmu)) 2069 ntota1 = max(1,nbas(isymmu)) 2070 ntota = max(1,norb1(isyma)) 2071 ntotj = max(1,nrhf(isymj)) 2072c 2073 call dgemm('T','N',norb1(isyma), 2074 & nrhf(isymj),mbas1(isymmu),one,work(koffc), 2075 & ntota1,work(kxajkl),ntotmu,zero, 2076 & work(kxpjkl),ntota) 2077 kxpjkl = kxpjkl + norb1(isyma)*nrhf(isymj) 2078 end do 2079 end do 2080 end do 2081 end do 2082 end do 2083 2084 2085 if (froimp) then 2086 ksajkl = ksajklf 2087 kspjkl = kspjklf 2088 ksijkl = ksijklf 2089 endif 2090 ksqijklz = ksqjkl 2091 do isymkl = 1, nsym 2092 isymaj = isymkl 2093 do isyml =1, nsym 2094 isymk = muld2h(isymkl,isyml) 2095 do k = 1, nrhf(isymk) 2096 do l = 1, nrhf(isyml) 2097 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 2098 do isym =1, nsym 2099 isyma = muld2h(1,isym) 2100 isymj = muld2h(isymkl,isym) 2101 do i = 1,nrhf(isymj) 2102 call dcopy(nrhffr(isym),work(ksijklf),1, 2103 & work(ksqjkl),1) 2104 ksqjkl = ksqjkl + nrhffr(isym) 2105 ksijklf = ksijklf+ nrhffr(isym) 2106 call dcopy(nrhf(isym),work(kspjkl),1, 2107 & work(ksqjkl),1) 2108 ksqjkl = ksqjkl + nrhf(isym) 2109 kspjkl = kspjkl + nrhf(isym) 2110 call dcopy(nvir0(isym),work(ksajkl),1, 2111 & work(ksqjkl),1) 2112 ksqjkl = ksqjkl + nvir0(isym) 2113 ksajkl = ksajkl + nvir0(isym) 2114 enddo 2115 enddo 2116 enddo 2117 enddo 2118 enddo 2119 enddo 2120 2121 2122 call daxpy(ntklpj,one,work(ksqijklz),1,work(kxpjklz),1) 2123c 2124 2125c-------------------------------------------------------------------- 2126c Calculate X_{pj}^{kl}*D_{kl}^{ij} 2127c--------------------------------------------------------------------- 2128 2129 call dzero(work(kxdmmu),nr1bast) 2130 do isymkl = 1, nsym 2131 isymaj = isymkl 2132 isymij = isymkl 2133 do isyml =1, nsym 2134 isymk = muld2h(isymkl,isyml) 2135 do k = 1, nrhf(isymk) 2136 do l = 1, nrhf(isyml) 2137 idxkl=imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 2138 do isyma =1, nsym 2139 isymj = muld2h(isymaj,isyma) 2140 isymi = muld2h(isymij,isymj) 2141 koffd = 1 + igamsq(isymij,isymkl) 2142 & + nmatij(isymij)*(idxkl-1)+imatij(isymi,isymj) 2143 ntota = max(1,norb1(isyma)) 2144 ntotj = max(1,nrhf(isymj)) 2145 ntoti = max(1,nrhf(isymi)) 2146 ntota1= max(1,nbas(isyma)) 2147 lxdmmu= kxdmmu + imatmuj(isyma,isymi) 2148 call dgemm('N','T',norb1(isyma),nrhf(isymi), 2149 & nrhf(isymj),one,work(kxpjklz), 2150 & ntota,dklij(koffd),ntoti,one, 2151 & work(lxdmmu),ntota) 2152 2153 kxpjklz = kxpjklz + norb1(isyma)*nrhf(isymj) 2154 ksqijklz= ksqijklz+ norb1(isyma)*nrhf(isymj) 2155 ksijkl = ksijkl+ nrhffr(isyma)*nrhf(isymj) 2156 enddo 2157 enddo 2158 enddo 2159 enddo 2160 enddo 2161 2162 2163 2164 2165 call dzero(work(kxdpmu),mbas1t*nrhft) 2166 call dzero(work(kxdnumu),10*mbas1t**2) 2167 do isymaj = 1, nsym 2168 isymkl = isymaj 2169 isymij = isymaj 2170 noffset1 = 0 2171 do isyma =1, nsym 2172 isymmu = muld2h(1,isyma) 2173 isymj = muld2h(isymaj,isyma) 2174 isymi = muld2h(isymij,isymj) 2175 ntota = max(1,mbas1(isymi)) 2176 ntotj = max(1,norb1(isyma)) 2177 ntoti = max(1,mbas1(isymi)) 2178 ntota1= max(1,nbas(isymi)) 2179 koffc = kcmo + noffset + 2180 & iglmrv(isyma,isymi) 2181 lxdmmu = kxdmmu + imatmuj(isyma,isymi) 2182 lxdpmu = kxdpmu + imatmuj(isyma,isymi) 2183 lxdnumu = kxdnumu + imatmunu(isymi,isymi) 2184 2185 call dgemm('N','N',mbas1(isymi), 2186 & nrhf(isyma),norb1(isyma),one, 2187 & work(koffc), 2188 & ntota1,work(lxdmmu),ntotj,zero, 2189 & work(lxdpmu),ntota) 2190 2191 noffset1 = noffset1 + nbas(isymi)*nrhffr(isymi) 2192 koffd = kcmo + noffset + noffset1 + 2193 & iglmrv1(isyma,isymi) 2194 2195 call dgemm('N','T',mbas1(isymi), 2196 & mbas1(isymi),nrhf(isyma),one,work(koffd), 2197 & ntota1,work(lxdpmu),ntoti,zero, 2198 & work(lxdnumu),ntota) 2199 2200 enddo 2201 enddo 2202 2203 2204 2205 call dzero(work(kzmunu),n2basx) 2206 call dzero(work(kdens),n2basx) 2207 koff = 0 2208 do isymij = 1,nsym 2209 isymj = isymij 2210 isymi = isymij 2211 koff = imatmunu(isymi,isymj) 2212 do j = 1,mbas1(isymj) 2213 do i = 1,mbas1(isymi) 2214 koff = koff+1 2215 idxij = iaodis(isymi,isymj)+nbas(isymi)* 2216 & (j-1)+i 2217 work(kzmunu+idxij-1) = work(kxdnumu+koff-1) 2218 enddo 2219 enddo 2220 enddo 2221 2222 2223 call cc_r12fsym(imatmunu,work(kzmunu),work(kdens)) 2224 2225 2226 call r12exdens(res,work(kcmo),work(kdens),work(kend1),lwrk1) 2227 2228 2229 call qexit('r12gradbd') 2230 end 2231c------------------------------------------------------------------------ 2232C=====================================================================* 2233C /* Deck r12exdens*/ 2234 subroutine r12exdens(res,cmo,dens,work,lwork) 2235 implicit none 2236#include "ccsdsym.h" 2237#include "ccorb.h" 2238#include "dummy.h" 2239#include "r12int.h" 2240#include "priunit.h" 2241#include "ccfro.h" 2242#include "inftap.h" 2243#include "ccsdinp.h" 2244ccc#include "inforb.h" 2245 integer kcmo,lwork,kfrsav,kfree,lfree,kfvao,kfcao,kdcao,kdvao,kdv 2246 integer pq,kcao,isym,ioff,ioff1,kexia,isymi,NNORB1 2247 integer kctao,koff,index,ISYMAI,isyma,nnbast 2248 integer idxkl,idxlk,idxij,idxji,idxklij,idxijkl,idxijlk 2249 integer idxjikl,isymj,isymij,isymkl,isymk,isyml,idxklji 2250 integer koffc,kexiaz,noffset,nvir0(8),kextia,kfctao 2251 integer icoun1,imatpq(8,8),icoun2,imatsq(8,8),ij,ji 2252 integer kajfr,isymaj,luvajkl,ncfai,kexfa,kexfaz 2253 logical onlyfc 2254c 2255 integer icoun3,imatfq(8,8),kxij,kxijij,ncvij,nvirf(8) 2256#include "maxorb.h" 2257#include "infinp.h" 2258 real*8 res(*),cmo(*),emcmy,zero,one,two,work(*),dens(*),densx 2259 parameter (zero = 0.0d0,one = 1.0d0, two = 2.0d0) 2260 2261 index(i,j) = max(i,j)*(max(i,j)-3)/2 +i+j 2262 2263 call qenter('r12exdens') 2264 2265 2266 nrhft = 0 2267 do isymi = 1, nsym 2268 nvir0(isymi) = norb1(isymi) - nrhfs(isymi) 2269 nvirf(isymi) = nvir(isymi) - nrhffr(isymi) 2270 nrhft = nrhft + nrhf(isymi) 2271 end do 2272 2273 do isymaj = 1,nsym 2274 isymij = isymaj 2275 ncfai = 0 2276 do isyma = 1,nsym 2277 isymj = muld2h(isymaj,isyma) 2278 isymi = muld2h(isymij,isymj) 2279 ncfai = ncfai + (norb1(isyma) - nrhfs(isyma))*nrhffr(isymi) 2280 enddo 2281 enddo 2282 2283 2284 do isymi = 1,nsym 2285 icoun1 = 0 2286 icoun2 = 0 2287 icoun3 = 0 2288 do isymj = 1,nsym 2289 isymk = muld2h(isymj,isymi) 2290 imatpq(isymk,isymj) = icoun1 2291 imatsq(isymk,isymj) = icoun2 2292 imatfq(isymk,isymj) = icoun3 2293 icoun1 = icoun1 + nvir(isymk)*(nvir(isymj)+1)/2 2294 icoun2 = icoun2 + nrhf(isymk)*nvir0(isymj) 2295 icoun3 = icoun3 + nrhffr(isymk)*nvir0(isymj) 2296 enddo 2297 enddo 2298 2299 do isymaj = 1,nsym 2300 isymij = isymaj 2301 ncvij = 0 2302 do isyma = 1,nsym 2303 isymj = muld2h(isymaj,isyma) 2304 isymi = muld2h(isymij,isymj) 2305 ncvij = ncvij + nrhffr(isyma)*nrhf(isymi) 2306 enddo 2307 enddo 2308 2309 2310 2311 2312 NBAST = 0 2313 NNBAST = 0 2314 DO I = 1,NSYM 2315 NBAST = NBAST + NBAS(I) 2316 NNBAST = NNBAST + NBAS(I)*(NBAS(I)+1)/2 2317 ENDDO 2318 2319 noffset = 0 2320 do isym = 1,nsym 2321 noffset = noffset+nbas(isym)*nrhf(isym) 2322 enddo 2323 2324 kfrsav = 1 2325 KFREE = 1 2326 LFREE = lwork 2327 CALL MEMGET('REAL',kfcao,N2BASX,WORK,KFREE,LFREE) 2328 CALL MEMGET('REAL',kfvao, 0,WORK,KFREE,LFREE) 2329 CALL MEMGET('REAL',kdcao,N2BASX,WORK,KFREE,LFREE) 2330 CALL MEMGET('REAL',kdvao, 0,WORK,KFREE,LFREE) 2331 CALL MEMGET('REAL',kdv , 0,WORK,KFREE,LFREE) 2332 CALL MEMGET('REAL',kcmo, nlamds,WORK,KFREE,LFREE) 2333 CALL MEMGET('REAL',kexia,nnbast,WORK,KFREE,LFREE) 2334 CALL MEMGET('REAL',kexfa,nnbast,WORK,KFREE,LFREE) 2335 CALL MEMGET('REAL',kajfr,ncfai ,WORK,KFREE,LFREE) 2336 CALL MEMGET('REAL',kxijij,ncvij,WORK,KFREE,LFREE) 2337 CALL MEMGET('REAL',kxij ,ncvij,WORK,KFREE,LFREE) 2338 2339 call erisfk(dens,nbast,1) 2340 call dscal(n2basx,0.50d0,dens,1) 2341 2342 2343 DIRFCK = .TRUE. 2344 R12TRA = .TRUE. 2345 R12INT = .FALSE. 2346 V12INT = .TRUE. 2347 NORXR = NORXR .OR. R12HYB 2348 mbsmax = 4 2349 2350 LUSUPM = -1 2351 call dzero(work(kfcao),n2basx) 2352 onlyfc = .TRUE. 2353 call fckmao(onlyfc,emcmy,work(kfcao),work(kfvao), 2354 & dens,work(kdvao),work(kdv),cmo(1+noffset),work(kfree),lfree) 2355 2356! Note: dens is destroyed in fckmao because of dirfck, but that doesn't 2357! matter as it is not used any more /Jan 2011 2358 2359c 2360 WRITE(LUPRI,'(/A,I1,A)') 2361 & ' Computation of exchange matrix done [',MBSMAX,']' 2362 2363 2364 R12INT = .TRUE. 2365 V12INT = .FALSE. 2366 call dzero(WORK(KDCAO),N2BASX) 2367 call dzero(WORK(kexia),NNBAST) 2368 CALL DCOPY(N2BASX,WORK(KFCAO),1,WORK(KDCAO),1) 2369c N2BASX = NBAST*NBAST! 2370 CALL DGETSP(NBAST,WORK(KDCAO),WORK(KFCAO)) 2371 IF (NSYM .GT. 1) CALL PKSYM1(WORK(KFCAO),WORK(KFCAO),NBAS,NSYM,2) 2372 CALL DZERO(WORK(kexia),NNBAST) 2373 kexiaz = kexia 2374 koffc= 1 2375 CALL UTHUB(WORK(KFCAO),WORK(kexia), 2376 * CMO(1+noffset),WORK(kfree),NSYM,nbas,nvir) 2377 2378 2379 koff = 0 2380 do isym = 1,nsym 2381 do p = 1+nrhfs(isym),norb1(isym) 2382 do q = 1+nrhffr(isym),nrhfs(isym) 2383 koff = imatsq(isym,isym)+nvir0(isym)*(q-nrhffr(isym)-1) 2384 & +p-nrhfs(isym) 2385 pq = imatpq(isym,isym)+index(p,q) 2386 res(koff) = 4.00d0*work(kexia+pq-1) 2387 enddo 2388 enddo 2389 enddo 2390 2391 if (froimp) then 2392 2393 koff = 0 2394 do isym = 1,nsym 2395 do p = 1+nrhfs(isym),norb1(isym) 2396 do q = 1,nrhffr(isym) 2397 koff = imatfq(isym,isym)+nvir0(isym) 2398 & *(q-1)+p - nrhfs(isym) 2399 pq = imatpq(isym,isym)+index(p,q) 2400 work(kajfr+koff-1) = 8.00d0*work(kexia+pq-1) 2401 enddo 2402 enddo 2403 enddo 2404 2405 luvajkl = -1 2406 call gpopen(luvajkl,'CCR12YAIFR','unknown',' ','unformatted', 2407 & idummy,.false.) 2408 write(luvajkl) (work(kajfr+i-1), i = 1,ncfai) 2409 call gpclose(luvajkl,'KEEP') 2410 2411 2412 endif 2413 2414 2415 CALL MEMREL('r12exdens',WORK,KFRSAV,KFRSAV,KFREE,LFREE) 2416 call qexit('r12exdens') 2417 end 2418c-------------------------------------------------------------------------- 2419