1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19 20*======================================================================= 21 subroutine cc_r12vxint(work,lwork,ltstxint) 22*----------------------------------------------------------------------- 23* Purpose: Calculate VX-intermediate for CC-R12 response and store 24* on file 25* 26* filer12 filename for half transformed r12 integrals with 27* no or one auxliary basis function 28* filer12_2 filename for half transformed r12 integrals with 29* two auxiliary basis functions 30* ltstxint only for testing: 31* if true: ignore V (i.e. set V to unit matrix), 32* and compute ONLY that contribution, which involves 33* only sums over the AO-basis (NOT the aux-basis). 34* The result should be the same as the corresponding 35* part of the R12 overlap matrix. 36* 37* Note: This routine does NOT work using CABS! 38* 39* by C. Neiss winter 2004/2005 40*----------------------------------------------------------------------- 41 42 implicit none 43#include "priunit.h" 44#include "ccorb.h" 45#include "ccsdsym.h" 46#include "ccsdinp.h" 47#include "r12int.h" 48#include "ccr12int.h" 49#include "ccrspprp.h" 50#include "second.h" 51 52C----------------------------------------------------------------------- 53C variables/parameter: 54C----------------------------------------------------------------------- 55 logical locdbg,lauxdelta,lauxd,mkvamkl,lhtf,lauxbeta,ltstxint 56 parameter (locdbg = .false.) 57 character*8 label1,lab123(3) 58 character*8 filer12,filer12_2 59 integer nr1orb(8),nr1xorb(8),nr1bas(8),nr1xbas(8),nr2bas,n2bst1(8) 60 integer ir1xbas(8,8),ir1xorb(8,8),ibasx(8) 61 integer ir1orb(8,8),ir1bas(8,8),ir2bas(8,8),iaodis1(8,8) 62 integer ir2xbas(8,8),irgkl(8,8),nrgkl(8) 63 integer nrxgkl(8),irxgkl(8,8),nalphaj(8),ialphaj(8,8) 64 integer nr2bas2,ir2bas2(8,8),ir2xbas2(8,8) 65 integer iaodis2(8,8),nallbas,nbasxtd(8),norbxtd(8),nlamdsxtd 66 integer iaomo(8,8),naomo(8),n2bstxtd(8),iaomox(8,8) 67 integer icount1,icount2,icount3,icount4 68 integer nsymx,norbtsx,nbastx,nlamdsx,nrhfsx(8),norbsx(8) 69 integer isymb,isymd,isymbd,isymp,isymkappa,isymbkl,isymkl,isymmn, 70 & isymmk,isymnl 71 integer lwork,kend1,kend2,kend3,kend4,kend5, 72 & lwrk1,lwrk2,lwrk3,lwrk4,lwrk5 73 integer kcmo,kcmox,k12back,kr12bkl,kr12bkl2,krvbkl,krvbkl2 74 integer kprpao,kprpmo,kvxintsq,kvxint,koff1,koff2 75 integer koffv,koffr,koffrv,koffrv2,koffc,koffres,koffx 76 integer isympt,imatrix,ierr,iprp,lusifc,idel,idelta,idum,iopt 77 integer luvxint 78 integer kappa,beta,ienddelta,istartdelta 79 80#if defined (SYS_CRAY) 81 real zero,one,work(*),fac,ddot,time,dummy 82#else 83 double precision zero,one,work(*),fac,ddot,time,dummy 84#endif 85 parameter (zero = 0.0D0, one = 1.0D0, dummy = -1.23D45) 86 data lab123/'********','********','********'/ 87 88C----------------------------------------------------------------------- 89C initialisation 90C----------------------------------------------------------------------- 91 call qenter ('cc_r12vxint') 92 if (locdbg) then 93 write(lupri,*) 'Entered CC_R12VXINT' 94 write(lupri,*) 'memory available: ',lwork,' double words' 95 call flshfo(lupri) 96 end if 97 time = second() 98C 99 if (R12CBS) call quit('CC_R12VXINT does not work with CABS!') 100C 101 filer12 = frhtf 102 filer12_2 = frhtf2 103C 104C additional dimensions/offsets: 105 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 106 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas, 107 & ir1xbas,ir2bas,ir2xbas,irgkl,irxgkl,iaodis1, 108 & nalphaj,ialphaj) 109 110 nallbas = 0 111 nlamdsxtd = 0 112 do i = 1, nsym 113 nbasxtd(i) = mbas1(i) + mbas2(i) 114 norbxtd(i) = norb1(i) + norb2(i) 115 nallbas = nallbas + nbasxtd(i) 116 nlamdsxtd = nlamdsxtd + nbasxtd(i)*norbxtd(i) 117 end do 118 119 ibasx(1) = 0 120 do i = 2, nsym 121 ibasx(i) = ibasx(i-1) + mbas2(i-1) 122 end do 123 124 do isymbd = 1, nsym 125 icount1 = 0 126 icount2 = 0 127 icount3 = 0 128 icount4 = 0 129 do isymd = 1, nsym 130 isymb = muld2h(isymbd,isymd) 131 iaodis2(isymb,isymd) = icount1 132 iaomo(isymb,isymd) = icount2 133 ir2bas2(isymb,isymd) = icount3 134 ir2xbas2(isymb,isymd) = icount4 135 icount1 = icount1 + nbasxtd(isymb)*nbasxtd(isymd) 136 icount2 = icount2 + nbasxtd(isymb)*norbxtd(isymd) 137 icount3 = icount3 + nr1xbas(isymb)*nr1bas(isymd) 138 icount4 = icount4 + nr1xbas(isymb)*nr1xbas(isymd) 139 end do 140 n2bstxtd(isymbd) = icount1 141 naomo(isymbd) = icount2 142 end do 143 144 nr2bas2 = 0 145 do i = 1, nsym 146 nr2bas2 = nr2bas2 + nr1bas(i)*nr1xbas(i) 147 end do 148 149C if locdbg do some printout of variables: 150 if (locdbg) then 151 do i = 1, nsym 152 write(lupri,*) 'nbas(',i,') = ',nbas(i) 153 write(lupri,*) 'nbasxtd(',i,') = ',nbasxtd(i) 154 write(lupri,*) 'mbas1(',i,') = ',mbas1(i) 155 write(lupri,*) 'ibas(',i,') = ',ibas(i) 156 write(lupri,*) 'mbas2(',i,') = ',mbas2(i) 157 write(lupri,*) 'ibasx(',i,') = ',ibasx(i) 158 write(lupri,*) 'nr1bas(',i,') = ',nr1bas(i) 159 write(lupri,*) 'nr1xbas(',i,') = ',nr1xbas(i) 160 write(lupri,*) 'nrhf(',i,') = ',nrhf(i) 161 write(lupri,*) 'nvir(',i,') = ',nvir(i) 162 write(lupri,*) 'norb(',i,') = ',norb(i) 163 write(lupri,*) 'norbxtd(',i,') = ',norbxtd(i) 164 write(lupri,*) 'norb1(',i,') = ',norb1(i) 165 write(lupri,*) 'norb2(',i,') = ',norb2(i) 166 write(lupri,*) 'nmatij(',i,') = ',nmatij(i) 167 write(lupri,*) 'nr1orb(',i,') = ',nr1orb(i) 168 write(lupri,*) 'nr1xorb(',i,') = ',nr1xorb(i) 169 write(lupri,*) 'nrgkl(',i,') = ',nrgkl(i) 170 write(lupri,*) 'nrxgkl(',i,') = ',nrxgkl(i) 171 write(lupri,*) 'ntr12am(',i,') = ',ntr12am(i) 172 write(lupri,*) 'ntr12sq(',i,') = ',ntr12sq(i) 173 end do 174 do i = 1, nsym 175 do j = 1, nsym 176 write(lupri,*) 'irxgkl(',i,',',j,') = ',irxgkl(i,j) 177 write(lupri,*) 'imatav(',i,',',j,') = ',imatav(i,j) 178 write(lupri,*) 'iaodis(',i,',',j,') = ',iaodis(i,j) 179 write(lupri,*) 'iaodis2(',i,',',j,') = ',iaodis2(i,j) 180 end do 181 end do 182 write(lupri,*) 'n2basx = ',n2basx 183 write(lupri,*) 'nlamdt = ',nlamdt 184 write(lupri,*) 'nlamdsxtd = ',nlamdsxtd 185 call FLSHFO(LUPRI) 186 end if 187 188C----------------------------------------------------------------------- 189C open output file 190C----------------------------------------------------------------------- 191 luvxint = -1 192 call gpopen(luvxint,'CCR12VXINT','UNKNOWN',' ','UNFORMATTED', 193 & idum,.false.) 194 195C----------------------------------------------------------------------- 196C read in MO coefficients 197C----------------------------------------------------------------------- 198 lusifc = -1 199 call gpopen(lusifc,'SIRIFC','OLD',' ','UNFORMATTED',idum,.false.) 200 ! read dimensions for CMO coefficients for full basis (nlamdsx) 201 rewind(lusifc) 202 call mollab('FULLBAS ',lusifc,lupri) 203 read(lusifc) nsymx,norbtsx,nbastx,nlamdsx,(nrhfsx(i),i=1,nsym), 204 & (norbsx(i),i=1,nsym) 205C allocate memory for MO coefficients: 206 kcmo = 1 207 kend1 = kcmo + nlamdsxtd 208 kcmox = kend1 209 kend2 = kcmox + nlamdsx 210 lwrk2 = lwork - kend2 211 if (lwrk2 .lt. 0) then 212 call quit ('Insufficient work memory in CC_R12VXINT at #1') 213 end if 214 call dzero(work(kcmo),nlamdsxtd) 215 call dzero(work(kcmox),nlamdsx) 216 read(lusifc) 217 read(lusifc) (work(kcmox+i-1),i=1,nlamdsx) 218 call gpclose(lusifc,'KEEP') 219 if (locdbg) then 220 write(lupri,*) 'nsymx, norbtsx, nbastx, nlamdsx: ', 221 & nsymx, norbtsx, nbastx, nlamdsx 222 do i = 1, nsym 223 write(lupri,*) 'nrhfsx(',i,') = ', nrhfsx(i) 224 write(lupri,*) 'norbsx(',i,') = ', norbsx(i) 225 end do 226 do isymbd = 1, nsym 227 icount2 = 0 228 do isymd = 1, nsym 229 isymb = muld2h(isymbd,isymd) 230 iaomox(isymb,isymd) = icount2 231 icount2 = icount2 + nbasxtd(isymb)*norbsx(isymd) 232 end do 233 end do 234 235 call around('MO-coefficient matrix incl. redundant orbitals') 236 do i = 1, nsym 237 koffc = iaomox(i,i) + kcmox 238 write(lupri,*) 'Symmetry number:', i 239 if (norbsx(i) .eq. 0) then 240 write(lupri,*) 'This symmetry is empty' 241 else 242 call output(work(koffc),1,nbasxtd(i),1,norbsx(i), 243 & nbasxtd(i),norbsx(i),1,lupri) 244 end if 245 end do 246 end if 247 248C delete redundant occupied orbitals from CMO matrix: 249 koff1 = kcmox 250 koff2 = kcmo 251 do i = 1, nsym 252 koff1 = koff1 + nbasxtd(i)*nrhfsx(i) 253 call dcopy(nbasxtd(i)*norbxtd(i),work(koff1),1,work(koff2),1) 254 koff1 = koff1 + nbasxtd(i)*norbxtd(i) 255 koff2 = koff2 + nbasxtd(i)*norbxtd(i) 256 end do 257 if (locdbg) then 258 call around('MO-coefficient matrix') 259 do i = 1, nsym 260 koffc = iaomo(i,i) + kcmo 261 write(lupri,*) 'Symmetry number:', i 262 if (norbxtd(i) .eq. 0) then 263 write(lupri,*) 'This symmetry is empty' 264 else 265 call output(work(koffc),1,nbasxtd(i),1,norbxtd(i), 266 & nbasxtd(i),norbxtd(i),1,lupri) 267 end if 268 end do 269 end if 270C 271C----------------------------------------------------------------------- 272C start loop over operators V 273C----------------------------------------------------------------------- 274 IF (LOCDBG) THEN 275 CALL AROUND('LABELS on Operator-List:') 276 DO IPRP = 1, NPRLBL_CC 277 WRITE(LUPRI,*) PRPLBL_CC(IPRP) 278 END DO 279 END IF 280 281 DO IPRP = 1, NPRLBL_CC 282 LABEL1 = PRPLBL_CC(IPRP) 283 IF (LABEL1(1:5).NE.'HAM0 ') THEN 284C 285C----------------------------------------------------------------------- 286C read in V integrals 287C----------------------------------------------------------------------- 288C allocate memory for V integrals: 289 !since we do not know the symmetry of V in advance, it is more 290 !memory allocated than needed in general! 291 kprpao = kend1 292 kend2 = kprpao + nallbas*nallbas 293 lwrk2 = lwork - kend2 294 if (locdbg) write(lupri,*) 'lwrk2 at #2 = ',lwrk2 295 if (lwrk2 .lt. 0) then 296 call quit ('Insufficient work memory in CC_R12VXINT at #2') 297 end if 298 !just to be on the safe side: 299 call dzero(work(kprpao),nallbas*nallbas) 300 301 !symmetry of V is determined here 302 CALL CCPRPAO(LABEL1,.FALSE.,WORK(KPRPAO),ISYMPT,IMATRIX,IERR, 303 & WORK(KEND2),LWRK2) 304 305 if (ierr .gt. 0 ) then 306 call quit('property not found on file') 307 end if 308 309 if (isympt .eq. 0) then 310 write(lupri,*) 'WARNING: Symmetry for operator could not be'// 311 & ' determined!' 312 write(lupri,*) 'WARNING: Will use IRREP = 1' 313 isympt = 1 314 end if 315 316 if (locdbg) then 317 call around('Original V integrals') 318 write(lupri,*) 'LABEL = ',label1, 'ISYMPT = ',isympt 319 do isymkappa = 1, nsym 320 isymb = muld2h(isympt,isymkappa) 321 koffv = iaodis2(isymb,isymkappa) + kprpao 322 write(lupri,*) 'Symmetry block:', isymb, isymkappa 323 if (nbasxtd(isymkappa) .eq. 0) then 324 write(lupri,*) 'This symmetry is empty' 325 else 326 call output(work(koffv),1,nbasxtd(isymb),1, 327 & nbasxtd(isymkappa),nbasxtd(isymb), 328 & nbasxtd(isymkappa),1,lupri) 329 end if 330 end do 331 end if 332 333C----------------------------------------------------------------------- 334C transform one index to MO 335C----------------------------------------------------------------------- 336 kprpmo = kend2 337 kend3 = kprpmo + n2bstxtd(isympt) 338 lwrk3 = lwork - kend3 339 if (locdbg) write(lupri,*) 'lwrk3 at #3 = ',lwrk3 340 if (lwrk3 .lt. 0) then 341 call quit ('Insufficient work memory in CC_R12VXINT at #3') 342 end if 343 344 do isymkappa = 1, nsym 345 isymb = muld2h(isympt,isymkappa) 346 isymp = isymkappa 347 koffv = iaodis2(isymb,isymkappa) + kprpao 348 koffc = iaomo(isymkappa,isymp) + kcmo 349 koffres = iaomo(isymb,isymp) + kprpmo 350 call dgemm('N','N',nbasxtd(isymb),norbxtd(isymp), 351 & nbasxtd(isymkappa),one,work(koffv), 352 & max(1,nbasxtd(isymb)),work(koffc), 353 & max(1,nbasxtd(isymkappa)), 354 & zero,work(koffres),max(1,nbasxtd(isymb))) 355c if (locdbg) then 356c call around('Half transformed V integrals') 357c write(lupri,*) 'symmetry block: ',isymb, isymkappa 358c call output(work(koffres),1,nbasxtd(isymb),1, 359c & norbxtd(isymp),nbasxtd(isymb),norbxtd(isymp), 360c & 1,lupri) 361c call around('MO coefficient sub-matrix') 362c call output(work(koffc),1,nbasxtd(isymkappa),1, 363c & norbxtd(isymp),nbasxtd(isymkappa), 364c & norbxtd(isymp),1,lupri) 365c end if 366 367C----------------------------------------------------------------------- 368C backtransform to AO => one index in contravariant basis 369C----------------------------------------------------------------------- 370 call dgemm('N','T',nbasxtd(isymb),nbasxtd(isymkappa), 371 & norbxtd(isymp), 372 & one,work(koffres),max(1,nbasxtd(isymb)), 373 & work(koffc),max(1,nbasxtd(isymkappa)), 374 & zero,work(koffv),max(1,nbasxtd(isymb))) 375 376 end do 377 378 if (ltstxint) then 379 !V should be the (-1)-unit matrix 380 isympt = 1 381 call dzero(work(kprpao),nallbas*nallbas) 382 do isymkappa = 1, nsym 383 isymb = isymkappa 384 do kappa = 1, nbasxtd(isymkappa) 385 beta = kappa 386 koffv = iaodis2(isymb,isymkappa) + kprpao-1 + 387 & nbasxtd(isymb)*(kappa-1) + beta 388 work(koffv) = -one 389 end do 390 end do 391 end if 392 393 if (locdbg) then 394 call around('Half back transformed V integrals') 395 write(lupri,*) 'LABEL = ',label1 396 do isymkappa = 1, nsym 397 isymb = muld2h(isympt,isymkappa) 398 koffv = iaodis2(isymb,isymkappa) + kprpao 399 write(lupri,*) 'Symmetry block:', isymb, isymkappa 400 if (nbasxtd(isymkappa) .eq. 0) then 401 write(lupri,*) 'This symmetry is empty' 402 else 403 call output(work(koffv),1,nbasxtd(isymb),1, 404 & nbasxtd(isymkappa),nbasxtd(isymb), 405 & nbasxtd(isymkappa),1,lupri) 406 end if 407 end do 408 end if 409 410C----------------------------------------------------------------------- 411C allocate memory for result arrays (square matrix and lower triangle) 412C----------------------------------------------------------------------- 413 kvxintsq = kend2 414 kvxint = kvxintsq + nr12r12sq(isympt) 415 kend3 = kvxint + nr12r12p(isympt) 416 lwrk3 = lwork - kend3 417 if (locdbg) write(lupri,*) 'lwrk3 at #4 = ',lwrk3 418 if (lwrk3 .lt. 0) then 419 call quit('Insufficient memory in CC_R12VXINT at #4') 420 end if 421C zero out result array: 422 call dzero(work(kvxintsq),nr12r12sq(isympt)) 423 call dzero(work(kvxint),nr12r12p(isympt)) 424 425C----------------------------------------------------------------------- 426C start loop over delta/delta': 427C----------------------------------------------------------------------- 428 do isymd = 1, nsym 429 430C allocate memory for "(r*V)" integrals: 431 krvbkl = kend3 432 krvbkl2 = krvbkl + nrgkl(muld2h(isympt,isymd)) 433 kend4 = krvbkl2 + nrxgkl(muld2h(isympt,isymd)) 434 lwrk4 = lwork - kend4 435 if (locdbg) write(lupri,*) 'lwrk4 at #5 = ',lwrk4 436 if (lwrk4 .lt. 0) then 437 call quit('Insufficient work memory in CC_R12VXINT at #5') 438 end if 439c call dzero(work(krvbkl),nrgkl(isymd)) 440c call dzero(work(krvbkl2),nrxgkl(isymd)) 441 442C allocate memory for r integrals: 443 kr12bkl = kend4 444 kr12bkl2 = kr12bkl + nrgkl(isymd) 445 kend5 = kr12bkl2 + nrxgkl(isymd) 446 lwrk5 = lwork - kend5 447 if (locdbg) write(lupri,*) 'lwrk5 at #6 = ',lwrk5 448 if (lwrk5 .lt. 0) then 449 call quit('Insufficient work memory in CC_R12VXINT at #6') 450 end if 451c call dzero(work(kr12bkl),nrgkl(isymd)) 452c call dzero(work(kr12bkl2),nrxgkl(isymd)) 453 454 istartdelta = 1 455 ienddelta = nbasxtd(isymd) 456 if (ltstxint) ienddelta = mbas1(isymd) 457c 458 do idelta = istartdelta, ienddelta 459 !is delta auxiliary function? 460 if (idelta .gt. mbas1(isymd)) then 461 lauxdelta = .true. 462 else 463 lauxdelta = .false. 464 end if 465C determine idel: 466 idel = idelta + ibas(isymd) + ibasx(isymd) 467 468C----------------------------------------------------------------------- 469C read in r integrals 470C----------------------------------------------------------------------- 471 lauxd = .TRUE. 472 473 lauxbeta = .FALSE. 474 ! read r^(delta)_(kappa,mn) or r^(delta')_(kappa,mn) 475 call cc_r12getrint(work(kr12bkl),idel,isymd,nr1bas,ir1bas, 476 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 477 & nrhfb,nmatkl,imatkl,ibasx,lauxd,lauxbeta, 478 & filer12,work(kend5),lwrk5) 479 if (.not. ltstxint) then 480 ! read r^(delta)_(kappa',mn) or r^(delta')_(kappa',mn) 481 lauxbeta = .TRUE. 482 call cc_r12getrint(work(kr12bkl2),idel,isymd,nr1xbas, 483 & ir1bas,nr2bas2,ir2bas2,nrxgkl,irxgkl, 484 & ir1xbas,ir2xbas2, 485 & nrhfb,nmatkl,imatkl,ibasx,lauxd, 486 & lauxbeta,filer12_2,work(kend5),lwrk5) 487 end if 488 489 if (locdbg) then 490 call around('Half transformed r integrals') 491 write(lupri,*) 'idelta, idel, isymd: ',idelta,idel,isymd 492 do isymkappa = 1, nsym 493 isymmn = muld2h(isymd,isymkappa) 494 write(lupri,*) 'Symmetry kappa, mn:',isymkappa,isymmn 495 write(lupri,*) 'kappa is AO basis function:' 496 koffr = irgkl(isymkappa,isymmn) + kr12bkl 497 call output(work(koffr),1,mbas1(isymkappa),1, 498 & nmatkl(isymmn),mbas1(isymkappa), 499 & nmatkl(isymmn),1,lupri) 500 if (min(mbas1(isymkappa),nmatkl(isymmn)).eq.0) then 501 write(lupri,*) 'This symmetry block is empty' 502 end if 503 if (.not. ltstxint) then 504 write(lupri,*) 'kappa is auxiliary basis function:' 505 koffr = irxgkl(isymkappa,isymmn) + kr12bkl2 506 call output(work(koffr),1,mbas2(isymkappa),1, 507 & nmatkl(isymmn),mbas2(isymkappa), 508 & nmatkl(isymmn),1,lupri) 509 if (min(mbas2(isymkappa),nmatkl(isymmn)).eq.0) then 510 write(lupri,*) 'This symmetry block is empty' 511 end if 512 end if 513 end do 514 end if 515 516 517C----------------------------------------------------------------------- 518C contract r with V 519C----------------------------------------------------------------------- 520 if (locdbg) then 521 call around('"V*r" integrals') 522 write(lupri,*) 'idelta, idel: ',idelta,idel 523 end if 524 do isymb = 1, nsym 525 isymkappa = muld2h(isympt,isymb) 526 isymmn = muld2h(isymd,isymkappa) 527 if (lauxdelta) then 528 fac = one 529 else 530 fac = -one 531 end if 532C 533 koffv = iaodis2(isymb,isymkappa) + kprpao 534 koffr = irgkl(isymkappa,isymmn) + kr12bkl 535 koffrv = irgkl(isymb,isymmn) + krvbkl 536 call dgemm('N','N', mbas1(isymb),nmatkl(isymmn), 537 & mbas1(isymkappa),fac,work(koffv), 538 & max(1,nbasxtd(isymb)), 539 & work(koffr),max(1,mbas1(isymkappa)),zero, 540 & work(koffrv),max(1,mbas1(isymb))) 541c if (locdbg) then 542c call output(work(koffrv),1,mbas1(isymb),1, 543c & nmatkl(isymmn),mbas1(isymb), 544c & nmatkl(isymmn),1,lupri) 545c end if 546 koffv = iaodis2(isymb,isymkappa) + 547 & nbasxtd(isymb)*mbas1(isymkappa) 548 & + kprpao 549 koffr = irxgkl(isymkappa,isymmn) + kr12bkl2 550 if (.not. ltstxint) then 551 call dgemm('N','N',mbas1(isymb),nmatkl(isymmn), 552 & mbas2(isymkappa),-fac,work(koffv), 553 & max(1,nbasxtd(isymb)), 554 & work(koffr),max(1,mbas2(isymkappa)),one, 555 & work(koffrv),max(1,mbas1(isymb))) 556 koffv = iaodis2(isymb,isymkappa) + mbas1(isymb) + kprpao 557 koffr = irgkl(isymkappa,isymmn) + kr12bkl 558 koffrv2 = irxgkl(isymb,isymmn) + krvbkl2 559 call dgemm('N','N',mbas2(isymb),nmatkl(isymmn), 560 & mbas1(isymkappa),-fac,work(koffv), 561 & max(1,nbasxtd(isymb)), 562 & work(koffr),max(1,mbas1(isymkappa)),zero, 563 & work(koffrv2),max(1,mbas2(isymb))) 564c if (locdbg) then 565c call output(work(koffrv2),1,mbas2(isymb),1, 566c & nmatkl(isymmn),mbas2(isymb), 567c & nmatkl(isymmn),1,lupri) 568c end if 569 koffv = iaodis2(isymb,isymkappa) + 570 & nbasxtd(isymb)*mbas1(isymkappa) + mbas1(isymb) + 571 & kprpao 572 koffr = irxgkl(isymkappa,isymmn) + kr12bkl2 573 call dgemm('N','N',mbas2(isymb),nmatkl(isymmn), 574 & mbas2(isymkappa),fac,work(koffv), 575 & max(1,nbasxtd(isymb)), 576 & work(koffr),max(1,mbas2(isymkappa)),one, 577 & work(koffrv2),max(1,mbas2(isymb))) 578 end if !ltstxint 579 580 if (locdbg) then 581 write(lupri,*) 'Symmetry beta, mn:',isymb,isymmn 582 write(lupri,*) 'beta is AO basis function:' 583 call output(work(koffrv),1,mbas1(isymb),1, 584 & nmatkl(isymmn),mbas1(isymb), 585 & nmatkl(isymmn),1,lupri) 586 if (min(mbas1(isymb),nmatkl(isymmn)).eq.0) then 587 write(lupri,*) 'This symmetry block is empty' 588 end if 589 if (.not. ltstxint) then 590 write(lupri,*) 'beta is auxiliary basis function:' 591 call output(work(koffrv2),1,mbas2(isymb),1, 592 & nmatkl(isymmn),mbas2(isymb), 593 & nmatkl(isymmn),1,lupri) 594 if (min(mbas2(isymb),nmatkl(isymmn)).eq.0) then 595 write(lupri,*) 'This symmetry block is empty' 596 end if 597 end if 598 end if 599 end do 600 601C----------------------------------------------------------------------- 602C read in R integrals (we use the memory where the r integrals were 603C stored, since these are not longer needed) 604C----------------------------------------------------------------------- 605 lauxd = .TRUE. 606 607 lauxbeta = .FALSE. 608 ! read R^(delta)_(beta,kl) or R^(delta')_(beta,kl) 609 call cc_r12getrint(work(kr12bkl),idel,isymd,nr1bas,ir1bas, 610 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 611 & nrhfb,nmatkl,imatkl,ibasx,lauxd,lauxbeta, 612 & fnback,work(kend5),lwrk5) 613 if (.not. ltstxint) then 614 lauxbeta = .TRUE. 615 ! read R^(delta)_(beta',kl) or R^(delta')_(beta',kl) 616 call cc_r12getrint(work(kr12bkl2),idel,isymd,nr1xbas, 617 & ir1bas,nr2bas2,ir2bas2,nrxgkl,irxgkl, 618 & ir1xbas,ir2xbas2, 619 & nrhfb,nmatkl,imatkl,ibasx,lauxd,lauxbeta, 620 & fnback2,work(kend5),lwrk5) 621 end if 622 623 if (locdbg) then 624 call around('Half back transformed R integrals') 625 write(lupri,*) 'idelta, idel: ',idelta,idel 626 do isymb = 1, nsym 627 isymkl = muld2h(isymd,isymb) 628 write(lupri,*) 'Symmetry beta, kl:',isymb,isymkl 629 write(lupri,*) 'beta is AO basis function:' 630 koffr = irgkl(isymb,isymkl) + kr12bkl 631 call output(work(koffr),1,mbas1(isymb),1, 632 & nmatkl(isymkl),mbas1(isymb), 633 & nmatkl(isymkl),1,lupri) 634 if (min(mbas1(isymb),nmatkl(isymkl)).eq.0) then 635 write(lupri,*) 'This symmetry block is empty' 636 end if 637 if (.not. ltstxint) then 638 write(lupri,*) 'beta is auxiliary basis function:' 639 koffr = irxgkl(isymb,isymkl) + kr12bkl2 640 call output(work(koffr),1,mbas2(isymb),1, 641 & nmatkl(isymkl),mbas2(isymb), 642 & nmatkl(isymkl),1,lupri) 643 if (min(mbas2(isymb),nmatkl(isymkl)).eq.0) then 644 write(lupri,*) 'This symmetry block is empty' 645 end if 646 end if 647 end do 648 end if 649 650C---------------------------------------------------------------------- 651C contract and add up over all deltas 652C---------------------------------------------------------------------- 653 if (locdbg) then 654 call around('VXINT before packing') 655 end if 656 do isymb = 1, nsym 657 isymbd = muld2h(isymd,isymb) 658 isymmn = muld2h(isymbd,isympt) 659 isymkl = muld2h(isymbd,1) 660 koffr = irgkl(isymb,isymkl) + kr12bkl 661 koffrv = irgkl(isymb,isymmn) + krvbkl 662 koffx = ir12r12sq(isymkl,isymmn) + kvxintsq 663c call output(work(koffr),1,mbas1(isymb),1,nmatkl(isymkl), 664c & mbas1(isymb),nmatkl(isymkl),1,lupri) 665c call output(work(koffrv),1,mbas1(isymb),1,nmatkl(isymmn), 666c & mbas1(isymb),nmatkl(isymmn),1,lupri) 667 call dgemm('T','N',nmatkl(isymkl),nmatkl(isymmn), 668 & mbas1(isymb),one,work(koffr),max(1,mbas1(isymb)), 669 & work(koffrv),max(1,mbas1(isymb)),one,work(koffx), 670 & max(1,nmatkl(isymkl))) 671 koffr = irxgkl(isymb,isymkl) + kr12bkl2 672 koffrv2 = irxgkl(isymb,isymmn) + krvbkl2 673c call output(work(koffr),1,mbas2(isymb),1,nmatkl(isymkl), 674c & mbas2(isymb),nmatkl(isymkl),1,lupri) 675c call output(work(koffrv2),1,mbas2(isymb),1,nmatkl(isymmn), 676c & mbas2(isymb),nmatkl(isymmn),1,lupri) 677 if (.not. ltstxint) then 678 call dgemm('T','N',nmatkl(isymkl),nmatkl(isymmn), 679 & mbas2(isymb),one,work(koffr),max(1,mbas2(isymb)), 680 & work(koffrv2),max(1,mbas2(isymb)),one, 681 & work(koffx),max(1,nmatkl(isymkl))) 682 end if 683 if (locdbg) then 684 write(lupri,*) 'idelta, idel, isymd: ',idelta,idel,isymd 685 write(lupri,*) 'isymkl, isymmn :',isymkl,isymmn 686 write(lupri,*) 'nr12r12sq(isympt), 687 & ir12r12sq(isymkl,isymmn): ', 688 & nr12r12sq(isympt), 689 & ir12r12sq(isymkl,isymmn) 690 call output(work(koffx),1,nmatkl(isymkl),1, 691 & nmatkl(isymmn),nmatkl(isymkl),nmatkl(isymmn), 692 & 1,lupri) 693 end if 694 end do 695 696 697C---------------------------------------------------------------------- 698C end of loops 699C---------------------------------------------------------------------- 700 end do !idelta 701 end do !isymd 702 703C----------------------------------------------------------------------- 704C resort result into a symmetry packed triangular matrix and 705C apply the projection operator P_(mn)^(kl): 706C----------------------------------------------------------------------- 707 iopt = 2 708 call ccr12pck2(work(kvxint),isympt,.TRUE.,work(kvxintsq), 709 & 'N',iopt) 710 711 if (locdbg) then 712 call around('final VXINT before packing') 713 do isymmn = 1, nsym 714 isymkl = muld2h(isympt,isymmn) 715 koffx = ir12r12sq(isymkl,isymmn) + kvxintsq 716 write(lupri,*) 'Symmetry block:', isymkl, isymmn 717 if (nmatkl(isymkl).eq.0 .or. nmatkl(isymmn).eq.0) then 718 write(lupri,*) 'This symmetry is empty' 719 else 720 call output(work(koffx),1,nmatkl(isymkl),1, 721 & nmatkl(isymmn),nmatkl(isymkl),nmatkl(isymmn), 722 & 1,lupri) 723 end if 724 end do 725 write(lupri,*) 726 write(lupri,*) 'Norm^2 over symmetries: ', 727 & ddot(nr12r12sq(isympt),work(kvxintsq),1,work(kvxintsq),1) 728 call around('final VXINT after packing') 729 do isymnl = 1, nsym 730 isymmk = muld2h(isympt,isymnl) 731 koffx = ir12r12p(isymmk,isymnl) + kvxint 732 if (isymmk .lt. isymnl) then 733 write(lupri,*) 'Symmetry block: ', isymmk, isymnl 734 if (nmatkl(isymmk).eq.0 .or. nmatkl(isymnl).eq.0) then 735 write(lupri,*) 'This symmetry is empty' 736 else 737 call output(work(koffx),1,nmatkl(isymmk),1, 738 & nmatkl(isymnl),nmatkl(isymmk), 739 & nmatkl(isymnl),1,lupri) 740 end if 741 else if (isymmk .eq. isymnl) then 742 write(lupri,*) 'Symmetry block: ', isymmk, isymnl 743 if (nmatkl(isymmk) .eq. 0) then 744 write(lupri,*) 'This symmetry is empty' 745 else 746 call outpak(work(koffx),nmatkl(isymmk),1,lupri) 747 end if 748 end if 749 end do 750 end if 751 752C----------------------------------------------------------------------- 753C write result on file for future use 754C----------------------------------------------------------------------- 755 !place date and time in lab123 756 call getdat(lab123(2),lab123(3)) 757 758 !need to write always at least 4 ("long") items, otherwise 759 !problems with subroutine MOLLAB 760 write(luvxint) lab123, label1 761 write(luvxint) isympt,dummy,dummy,dummy,dummy 762 write(luvxint) (work(kvxint-1+i), i=1, max(4,nr12r12p(isympt))) 763 764 if (locdbg) then 765 write(lupri,*) 766 write(lupri,*) 'nr12r12p(isympt) = ',nr12r12p(isympt) 767 write(lupri,*) label1, isympt 768 write(lupri,*) 'norm:',ddot(nr12r12p(isympt),work(kvxint),1, 769 & work(kvxint),1) 770 write(lupri,*) (work(kvxint-1+i), i=1, nr12r12p(isympt)) 771 end if 772 773C----------------------------------------------------------------------- 774C end loop over operators V 775C----------------------------------------------------------------------- 776 end if 777 end do !operator V 778 779C----------------------------------------------------------------------- 780C close output file 781C----------------------------------------------------------------------- 782 call gpclose(luvxint,'KEEP') 783 784C----------------------------------------------------------------------- 785C end subroutine 786C----------------------------------------------------------------------- 787 time = second() - time 788 WRITE(LUPRI,'(1X,A)') 789 & 'Computation of CC-R12 VXINT response intermediates done' 790 WRITE(LUPRI,'(/1X,A,F7.2,A)') 791 & ' Time used for VXINT is ',time,' seconds' 792 793 if (locdbg) then 794 write(lupri,*) 'Leaving CC_R12VXINT' 795 call flshfo(lupri) 796 end if 797 798 call qexit('cc_r12vxint') 799 800 return 801 end 802 803*======================================================================= 804*======================================================================= 805 subroutine cc_r12vxint2(rpkql,work,lwork) 806*----------------------------------------------------------------------- 807* Purpose: Calculate VX-intermediate for CC-R12 response and store 808* on file 809* 810* !!!for use in MP2-R12!!! 811* 812* by C. Neiss jan. 2006 813*----------------------------------------------------------------------- 814 815 implicit none 816#include "priunit.h" 817#include "ccorb.h" 818#include "ccsdsym.h" 819#include "ccsdinp.h" 820#include "r12int.h" 821#include "ccr12int.h" 822#include "ccrspprp.h" 823#include "second.h" 824 825C----------------------------------------------------------------------- 826C variables/parameter: 827C----------------------------------------------------------------------- 828 logical locdbg,lauxq 829 parameter (locdbg = .false.) 830 character*8 label1,lab123(3) 831 integer iaodis2(8,8),nallbas,nbasxtd(8),norbxtd(8),nlamdsxtd 832 integer iaomo(8,8),naomo(8),n2bstxtd(8),iaomox(8,8) 833 integer icount1,icount2,idum 834 integer lwork,kend1,kend2,kend3,kend4, 835 & lwrk1,lwrk2,lwrk3,lwrk4 836 integer kcmo,kprpao,kprpmo,kvxintsq,kvxint 837 integer koff1,koff2,koffc,koffv,koffres,kr12pkl,kr12pkl2, 838 & krvpmn,krvpmn2,koffx 839 integer isym,isym1,isym2,isympt,isymb,isymkappa,isymp,isymq,isyms, 840 & isymkl,isymmn,isympq,idxq,istartq 841 integer ierr,iprp,imatrix,iopt 842 integer luvxint 843 844#if defined (SYS_CRAY) 845 real zero,one,work(*),rpkql(*),fac,ddot,time,dummy 846#else 847 double precision zero,one,work(*),rpkql(*),fac,ddot, 848 & time,dummy 849#endif 850 parameter (zero = 0.0D0, one = 1.0D0, dummy = -1.23D45) 851 data lab123/'********','********','********'/ 852 853C----------------------------------------------------------------------- 854C initialisation 855C----------------------------------------------------------------------- 856 call qenter ('cc_r12vxint') 857 if (locdbg) then 858 write(lupri,*) 'Entered CC_R12VXINT' 859 write(lupri,*) 'memory available: ',lwork,' double words' 860 call flshfo(lupri) 861 end if 862 time = second() 863C 864C additional dimensions/offsets: 865C 866 nallbas = 0 867 nlamdsxtd = 0 868 do isym = 1, nsym 869 nbasxtd(isym) = mbas1(isym) + mbas2(isym) 870 norbxtd(isym) = norb1(isym) + norb2(isym) 871 nallbas = nallbas + nbasxtd(isym) 872 nlamdsxtd = nlamdsxtd + nbasxtd(isym)*norbxtd(isym) 873 end do 874 875 do isym = 1, nsym 876 icount1 = 0 877 icount2 = 0 878 do isym2 = 1, nsym 879 isym1 = muld2h(isym,isym2) 880 iaodis2(isym1,isym2) = icount1 881 iaomo(isym1,isym2) = icount2 882 icount1 = icount1 + nbasxtd(isym1)*nbasxtd(isym2) 883 icount2 = icount2 + nbasxtd(isym1)*norbxtd(isym2) 884 end do 885 n2bstxtd(isym) = icount1 886 naomo(isym) = icount2 887 end do 888 889 if (locdbg) then 890 call around('r_12 integrals:') 891 call cc_prp(dummy,rpkql,1,0,1) 892 end if 893 894C----------------------------------------------------------------------- 895C open output file 896C----------------------------------------------------------------------- 897 luvxint = -1 898 call gpopen(luvxint,'CCR12VXINT','UNKNOWN',' ','UNFORMATTED', 899 & idum,.false.) 900 901C----------------------------------------------------------------------- 902C read in MO coefficients 903C----------------------------------------------------------------------- 904C allocate memory for MO coefficients: 905 kcmo = 1 906 kend1 = kcmo + nlamdsxtd 907 lwrk1 = lwork - kend1 908 if (lwrk1 .lt. 0) then 909 call quit ('Insufficient work memory in CC_R12VXINT at #1') 910 end if 911 912 call cc_r12cmo(work(kcmo),work(kend1),lwrk1) 913C 914C----------------------------------------------------------------------- 915C start loop over operators V 916C----------------------------------------------------------------------- 917 IF (LOCDBG) THEN 918 CALL AROUND('LABELS on Operator-List:') 919 DO IPRP = 1, NPRLBL_CC 920 WRITE(LUPRI,*) PRPLBL_CC(IPRP) 921 END DO 922 END IF 923 924 DO IPRP = 1, NPRLBL_CC 925 LABEL1 = PRPLBL_CC(IPRP) 926 IF (LABEL1(1:5).NE.'HAM0 ') THEN 927C 928C----------------------------------------------------------------------- 929C read in V integrals 930C----------------------------------------------------------------------- 931C allocate memory for V integrals: 932 !since we do not know the symmetry of V in advance, it is more 933 !memory allocated than needed in general! 934 kprpao = kend1 935 kend2 = kprpao + nallbas*nallbas 936 lwrk2 = lwork - kend2 937 if (locdbg) write(lupri,*) 'lwrk2 at #2 = ',lwrk2 938 if (lwrk2 .lt. 0) then 939 call quit ('Insufficient work memory in CC_R12VXINT at #2') 940 end if 941 !just to be on the safe side: 942 call dzero(work(kprpao),nallbas*nallbas) 943 944 !symmetry of V is determined here 945 CALL CCPRPAO(LABEL1,.FALSE.,WORK(KPRPAO),ISYMPT,IMATRIX,IERR, 946 & WORK(KEND2),LWRK2) 947 948 if (ierr .gt. 0 ) then 949 call quit('property not found on file') 950 end if 951 952 if (isympt .eq. 0) then 953 write(lupri,*) 'WARNING: Symmetry for operator could not be'// 954 & ' determined!' 955 write(lupri,*) 'WARNING: Will use IRREP = 1' 956 isympt = 1 957 end if 958 959 if (locdbg) then 960 call around('Original V integrals') 961 write(lupri,*) 'LABEL = ',label1, 'ISYMPT = ',isympt 962 do isymkappa = 1, nsym 963 isymb = muld2h(isympt,isymkappa) 964 koffv = iaodis2(isymb,isymkappa) + kprpao 965 write(lupri,*) 'Symmetry block:', isymb, isymkappa 966 if (nbasxtd(isymkappa) .eq. 0) then 967 write(lupri,*) 'This symmetry is empty' 968 else 969 call output(work(koffv),1,nbasxtd(isymb),1, 970 & nbasxtd(isymkappa),nbasxtd(isymb), 971 & nbasxtd(isymkappa),1,lupri) 972 end if 973 end do 974 end if 975 976C----------------------------------------------------------------------- 977C transform first index to MO 978C----------------------------------------------------------------------- 979 kprpmo = kend2 980 kend3 = kprpmo + n2bstxtd(isympt) 981 lwrk3 = lwork - kend3 982 if (locdbg) write(lupri,*) 'lwrk3 at #3 = ',lwrk3 983 if (lwrk3 .lt. 0) then 984 call quit ('Insufficient work memory in CC_R12VXINT at #3') 985 end if 986 987 do isymkappa = 1, nsym 988 isymb = muld2h(isympt,isymkappa) 989 isyms = isymkappa 990 koffv = iaodis2(isymb,isymkappa) + kprpao 991 koffc = iaomo(isymkappa,isyms) + kcmo 992 koffres = iaomo(isymb,isyms) + kprpmo 993 call dgemm('N','N',nbasxtd(isymb),norbxtd(isyms), 994 & nbasxtd(isymkappa),one,work(koffv), 995 & max(1,nbasxtd(isymb)),work(koffc), 996 & max(1,nbasxtd(isymkappa)), 997 & zero,work(koffres),max(1,nbasxtd(isymb))) 998 999C----------------------------------------------------------------------- 1000C transform second index to MO 1001C----------------------------------------------------------------------- 1002 isymp = isymb 1003 !overwrite original values 1004 call dzero(work(koffv),nbasxtd(isymp)*nbasxtd(isyms)) 1005 if (norbxtd(isymp)*norbxtd(isyms).gt. 1006 & nbasxtd(isymb)*nbasxtd(isymkappa)) then 1007 call quit('Dimension error in CC_R12VXINT2') 1008 end if 1009 koffc = iaomo(isymb,isymp) + kcmo 1010 call dgemm('T','N',norbxtd(isymp),norbxtd(isyms), 1011 & nbasxtd(isymb), 1012 & one,work(koffc),max(1,nbasxtd(isymb)), 1013 & work(koffres),max(1,nbasxtd(isymb)), 1014 & zero,work(koffv),max(1,norbxtd(isymp))) 1015 1016 end do 1017 1018 if (locdbg) then 1019 call around('V integrals in MO basis') 1020 write(lupri,*) 'LABEL = ',label1 1021 do isym2 = 1, nsym 1022 isym1 = muld2h(isympt,isym2) 1023 koffv = iaodis2(isym1,isym2) + kprpao 1024 write(lupri,*) 'Symmetry block:', isym1, isym2 1025 if (norbxtd(isym2) .eq. 0) then 1026 write(lupri,*) 'This symmetry is empty' 1027 else 1028 call output(work(koffv),1,norbxtd(isym1),1, 1029 & norbxtd(isym2),norbxtd(isym1), 1030 & norbxtd(isym2),1,lupri) 1031 end if 1032 end do 1033 end if 1034 1035C----------------------------------------------------------------------- 1036C allocate memory for result arrays (square matrix and lower triangle) 1037C----------------------------------------------------------------------- 1038 kvxintsq = kend2 1039 kvxint = kvxintsq + nr12r12sq(isympt) 1040 kend3 = kvxint + nr12r12p(isympt) 1041 lwrk3 = lwork - kend3 1042 if (locdbg) write(lupri,*) 'lwrk3 at #4 = ',lwrk3 1043 if (lwrk3 .lt. 0) then 1044 call quit('Insufficient memory in CC_R12VXINT at #4') 1045 end if 1046C zero out result array: 1047 call dzero(work(kvxintsq),nr12r12sq(isympt)) 1048 call dzero(work(kvxint),nr12r12p(isympt)) 1049 1050C----------------------------------------------------------------------- 1051C start loop over q/q': 1052C----------------------------------------------------------------------- 1053 do isymq = 1, nsym 1054 do isymp = 1, nsym 1055 isympq = muld2h(isymp,isymq) 1056 isymkl = isympq 1057 isymmn = muld2h(isympt,isympq) 1058 isyms = muld2h(isympt,isymp) 1059 1060C allocate memory for "(r*V)", r_12 integrals: 1061 if (r12cbs) then 1062 krvpmn2 = kend3 1063 kr12pkl2 = krvpmn2 + norb2(isymp)*nmatkl(isymmn) 1064 else 1065 krvpmn = kend3 1066 krvpmn2 = krvpmn + norb1(isymp)*nmatkl(isymmn) 1067 kr12pkl = krvpmn2 + norb2(isymp)*nmatkl(isymmn) 1068 kr12pkl2 = kr12pkl + max(norb1(isymp)*nmatkl(isymkl), 1069 & norb1(isyms)*nmatkl(isymmn)) 1070 end if 1071 kend4 = kr12pkl2+ max(norb2(isymp)*nmatkl(isymkl), 1072 & norb2(isyms)*nmatkl(isymmn)) 1073 lwrk4 = lwork - kend4 1074 if (locdbg) write(lupri,*) 'lwrk4 at #5 = ',lwrk4 1075 if (lwrk4 .lt. 0) then 1076 call quit('Insufficient work memory in CC_R12VXINT at #5') 1077 end if 1078 1079 if (r12cbs) then 1080 istartq = norb1(isymq) + 1 1081 else 1082 istartq = 1 1083 end if 1084 1085 do idxq = istartq, norbxtd(isymq) 1086 !is q auxiliary function? 1087 if (idxq .gt. norb1(isymq)) then 1088 lauxq = .true. 1089 else 1090 lauxq = .false. 1091 end if 1092 1093C----------------------------------------------------------------------- 1094C contract r_12 with V 1095C----------------------------------------------------------------------- 1096 if (lauxq) then 1097 fac = one 1098 else 1099 fac = -one 1100 end if 1101C 1102 if (.not.r12cbs) 1103 & call cc_r12sort(rpkql,1,work(kr12pkl),idxq,isymq, 1104 & 1,norb1(isyms),isyms,norbxtd,nrhfb) 1105 call cc_r12sort(rpkql,1,work(kr12pkl2),idxq,isymq, 1106 & norb1(isyms)+1,norb1(isyms)+norb2(isyms), 1107 & isyms,norbxtd,nrhfb) 1108 if (locdbg) then 1109 call around('r_(s,mn)^q integrals') 1110 write(lupri,*) 'isyms, isymmn: ',isyms, isymmn 1111 if (.not.r12cbs) then 1112 write(lupri,*) 's is MO basis function:' 1113 call output(work(kr12pkl),1,norb1(isyms), 1114 & 1,nmatkl(isymmn),norb1(isyms), 1115 & nmatkl(isymmn),1,lupri) 1116 write(lupri,*) 1117 end if 1118 write(lupri,*) 's is aux. function:' 1119 call output(work(kr12pkl2),1,norb2(isyms), 1120 & 1,nmatkl(isymmn),norb2(isyms), 1121 & nmatkl(isymmn),1,lupri) 1122 end if 1123C 1124 if (.not.r12cbs) then 1125 koffv = iaodis2(isymp,isyms) + kprpao 1126 call dgemm('N','N', norb1(isymp),nmatkl(isymmn), 1127 & norb1(isyms),fac,work(koffv), 1128 & max(1,norbxtd(isymp)), 1129 & work(kr12pkl),max(1,norb1(isyms)),zero, 1130 & work(krvpmn),max(1,norb1(isymp))) 1131 koffv = iaodis2(isymp,isyms) + 1132 & norbxtd(isymp)*norb1(isyms) + kprpao 1133 call dgemm('N','N',norb1(isymp),nmatkl(isymmn), 1134 & norb2(isyms),-fac,work(koffv), 1135 & max(1,norbxtd(isymp)), 1136 & work(kr12pkl2),max(1,norb2(isyms)),one, 1137 & work(krvpmn),max(1,norb1(isymp))) 1138 koffv = iaodis2(isymp,isyms) + norb1(isymp) + kprpao 1139 call dgemm('N','N',norb2(isymp),nmatkl(isymmn), 1140 & norb1(isyms),-fac,work(koffv), 1141 & max(1,norbxtd(isymp)), 1142 & work(kr12pkl),max(1,norb1(isyms)),zero, 1143 & work(krvpmn2),max(1,norb2(isymp))) 1144 else 1145 call dzero(work(krvpmn2),norb2(isymp)*nmatkl(isymmn)) 1146 end if 1147 koffv = iaodis2(isymp,isyms) + 1148 & norbxtd(isymp)*norb1(isyms) + norb1(isymp) + 1149 & kprpao 1150 call dgemm('N','N',norb2(isymp),nmatkl(isymmn), 1151 & norb2(isyms),fac,work(koffv), 1152 & max(1,norbxtd(isymp)), 1153 & work(kr12pkl2),max(1,norb2(isyms)),one, 1154 & work(krvpmn2),max(1,norb2(isymp))) 1155 1156 if (locdbg) then 1157 call around('r*V integrals') 1158 write(lupri,*) 'Symmetry p, mn:',isymp,isymmn 1159 if (.not.r12cbs) then 1160 write(lupri,*) 'p is MO basis function:' 1161 call output(work(krvpmn),1,norb1(isymp),1, 1162 & nmatkl(isymmn),norb1(isymp), 1163 & nmatkl(isymmn),1,lupri) 1164 if (min(norb1(isymp),nmatkl(isymmn)).eq.0) then 1165 write(lupri,*) 'This symmetry block is empty' 1166 end if 1167 end if 1168 write(lupri,*) 'p is aux. function:' 1169 call output(work(krvpmn2),1,norb2(isymp),1, 1170 & nmatkl(isymmn),norb2(isymp), 1171 & nmatkl(isymmn),1,lupri) 1172 if (min(norb2(isymp),nmatkl(isymmn)).eq.0) then 1173 write(lupri,*) 'This symmetry block is empty' 1174 end if 1175 end if 1176 1177C----------------------------------------------------------------------- 1178C get r_12 integrals for second time: 1179C----------------------------------------------------------------------- 1180 1181 if (.not.r12cbs) 1182 & call cc_r12sort(rpkql,1,work(kr12pkl),idxq,isymq, 1183 & 1,norb1(isymp),isymp,norbxtd,nrhfb) 1184 call cc_r12sort(rpkql,1,work(kr12pkl2),idxq,isymq, 1185 & norb1(isymp)+1,norb1(isymp)+norb2(isymp), 1186 & isymp,norbxtd,nrhfb) 1187 1188 if (locdbg) then 1189 call around('r_(p,kl)^q integrals') 1190 write(lupri,*) 'isymp, isymkl: ',isymp, isymkl 1191 if (.not.r12cbs) then 1192 write(lupri,*) 'p is MO basis function:' 1193 call output(work(kr12pkl),1,norb1(isymp), 1194 & 1,nmatkl(isymkl),norb1(isymp), 1195 & nmatkl(isymkl),1,lupri) 1196 write(lupri,*) 1197 end if 1198 write(lupri,*) 'p is aux. function:' 1199 call output(work(kr12pkl2),1,norb2(isymp), 1200 & 1,nmatkl(isymkl),norb2(isymp), 1201 & nmatkl(isymkl),1,lupri) 1202 end if 1203 1204C---------------------------------------------------------------------- 1205C contract and add up over all q: 1206C---------------------------------------------------------------------- 1207 koffx = ir12r12sq(isymkl,isymmn) + kvxintsq 1208 if (.not.r12cbs) 1209 & call dgemm('T','N',nmatkl(isymkl),nmatkl(isymmn), 1210 & norb1(isymp),one,work(kr12pkl),max(1,norb1(isymp)), 1211 & work(krvpmn),max(1,norb1(isymp)),one,work(koffx), 1212 & max(1,nmatkl(isymkl))) 1213 call dgemm('T','N',nmatkl(isymkl),nmatkl(isymmn), 1214 & norb2(isymp),one,work(kr12pkl2),max(1,norb2(isymp)), 1215 & work(krvpmn2),max(1,norb2(isymp)),one, 1216 & work(koffx),max(1,nmatkl(isymkl))) 1217 if (locdbg) then 1218 call around ('intermediate VX:') 1219 write(lupri,*) 'idxq, isymq: ',idxq,isymq 1220 write(lupri,*) 'isymkl, isymmn :',isymkl,isymmn 1221 write(lupri,*) 'nr12r12sq(isympt),'// 1222 & 'ir12r12sq(isymkl,isymmn): ', 1223 & nr12r12sq(isympt), 1224 & ir12r12sq(isymkl,isymmn) 1225 call output(work(koffx),1,nmatkl(isymkl),1, 1226 & nmatkl(isymmn),nmatkl(isymkl),nmatkl(isymmn), 1227 & 1,lupri) 1228 end if 1229 end do 1230 1231C---------------------------------------------------------------------- 1232C end of loops 1233C---------------------------------------------------------------------- 1234 end do !isymp 1235 end do !isymq 1236 1237C----------------------------------------------------------------------- 1238C resort result into a symmetry packed triangular matrix and 1239C apply the projection operator P_(mn)^(kl): 1240C----------------------------------------------------------------------- 1241 iopt = 2 1242 call ccr12pck2(work(kvxint),isympt,.TRUE.,work(kvxintsq), 1243 & 'N',iopt) 1244 1245 if (locdbg) then 1246 call around('final VXINT before packing') 1247 do isymmn = 1, nsym 1248 isymkl = muld2h(isympt,isymmn) 1249 koffx = ir12r12sq(isymkl,isymmn) + kvxintsq 1250 write(lupri,*) 'Symmetry block:', isymkl, isymmn 1251 if (nmatkl(isymkl).eq.0 .or. nmatkl(isymmn).eq.0) then 1252 write(lupri,*) 'This symmetry is empty' 1253 else 1254 call output(work(koffx),1,nmatkl(isymkl),1, 1255 & nmatkl(isymmn),nmatkl(isymkl),nmatkl(isymmn), 1256 & 1,lupri) 1257 end if 1258 end do 1259 write(lupri,*) 1260 write(lupri,*) 'Norm^2 over symmetries: ', 1261 & ddot(nr12r12sq(isympt),work(kvxintsq),1,work(kvxintsq),1) 1262 call around('final VXINT after packing') 1263 do isym2 = 1, nsym 1264 isym1 = muld2h(isympt,isym2) 1265 koffx = ir12r12p(isym1,isym2) + kvxint 1266 if (isym1 .lt. isym2) then 1267 write(lupri,*) 'Symmetry block: ', isym1, isym2 1268 if (nmatkl(isym1).eq.0 .or. nmatkl(isym2).eq.0) then 1269 write(lupri,*) 'This symmetry is empty' 1270 else 1271 call output(work(koffx),1,nmatkl(isym1),1, 1272 & nmatkl(isym2),nmatkl(isym1), 1273 & nmatkl(isym2),1,lupri) 1274 end if 1275 else if (isym1 .eq. isym2) then 1276 write(lupri,*) 'Symmetry block: ', isym1, isym2 1277 if (nmatkl(isym1) .eq. 0) then 1278 write(lupri,*) 'This symmetry is empty' 1279 else 1280 call outpak(work(koffx),nmatkl(isym1),1,lupri) 1281 end if 1282 end if 1283 end do 1284 end if 1285 1286C----------------------------------------------------------------------- 1287C write result on file for future use 1288C----------------------------------------------------------------------- 1289 !place date and time in lab123 1290 call getdat(lab123(2),lab123(3)) 1291 1292 !need to write always at least 4 ("long") items, otherwise 1293 !problems with subroutine MOLLAB 1294 write(luvxint) lab123, label1 1295 write(luvxint) isympt,dummy,dummy,dummy,dummy 1296 write(luvxint) (work(kvxint-1+i), i=1, max(4,nr12r12p(isympt))) 1297 1298 if (locdbg) then 1299 write(lupri,*) 1300 write(lupri,*) 'nr12r12p(isympt) = ',nr12r12p(isympt) 1301 write(lupri,*) label1, isympt 1302 write(lupri,*) 'norm:',ddot(nr12r12p(isympt),work(kvxint),1, 1303 & work(kvxint),1) 1304 write(lupri,*) (work(kvxint-1+i), i=1, nr12r12p(isympt)) 1305 end if 1306 1307C----------------------------------------------------------------------- 1308C end loop over operators V 1309C----------------------------------------------------------------------- 1310 end if 1311 end do !operator V 1312 1313C----------------------------------------------------------------------- 1314C close output file 1315C----------------------------------------------------------------------- 1316 call gpclose(luvxint,'KEEP') 1317 1318C----------------------------------------------------------------------- 1319C end subroutine 1320C----------------------------------------------------------------------- 1321 time = second() - time 1322 WRITE(LUPRI,'(1X,A)') 1323 & 'Computation of CC-R12 VXINT response intermediates done' 1324 WRITE(LUPRI,'(/1X,A,F7.2,A)') 1325 & ' Time used for VXINT is',time,' seconds' 1326 WRITE(LUPRI,*) 1327 1328 if (locdbg) then 1329 write(lupri,*) 'Leaving CC_R12VXINT' 1330 call flshfo(lupri) 1331 end if 1332 1333 call qexit('cc_r12vxint') 1334 1335 return 1336 end 1337 1338*======================================================================= 1339 1340*=====================================================================* 1341 subroutine cc_r12xi(XIR12,ISYMXI,TRXI,TR12,ISYMC,XINT,VXINT,ISYMV, 1342 & PRPAO,CMO,LAMH,TRV,WORK,LWORK) 1343C---------------------------------------------------------------------- 1344C purpose: calculate Xi vector for CCR12 response 1345C 1346C XIR12 result array of dimension ntr12sq(isymxi) 1347C ISYMXI symmetry of Xi 1348C TRXI r12 index pair leading in Xi: 'N' 1349C occ. index pair leading in Xi: 'T' 1350C TR12 R12 amplitudes of dimension ntr12sq(isymc) 1351C ISYMC symmetry of amplitudes 1352C XINT X-intermediate (R12 overlap integrals) of dim. nr12r12sq(1) 1353C VXINT VX-intermediate of dimension nr12r12sq(isymv) 1354C ISYMV symmetry of perturbation 1355C PRPAO one-electron perturbation operator integrals in AO-basis 1356C LAMH Lambda^h matrix 1357C CMO MO coefficient matrix 1358C TRV use normal ('N') or transposed ('T') matrix V_m^jtilde 1359C when constructing vc-intermediate 1360C 1361C Christian Neiss 2005 1362C---------------------------------------------------------------------- 1363 implicit none 1364#include "priunit.h" 1365#include "ccorb.h" 1366#include "ccsdsym.h" 1367#include "r12int.h" 1368#include "dummy.h" 1369 1370 logical locdbg 1371 parameter(locdbg = .false.) 1372 1373#if defined (SYS_CRAY) 1374 real zero,one 1375#else 1376 double precision zero,one 1377#endif 1378 parameter (one = 1.0d0) 1379 1380 integer isymxi,isymv,isymc 1381 integer kvcint 1382 integer lwork,kend1,lwrk1 1383 character*1 trv,trxi !allowed are 'N' and 'T' 1384 1385#if defined (SYS_CRAY) 1386 real xir12(*),tr12(*),xint(*),vxint(*),prpao(*),lamh(*), 1387 & cmo(*),work(*) 1388#else 1389 double precision xir12(*),tr12(*),xint(*),vxint(*),prpao(*), 1390 & lamh(*),cmo(*),work(*) 1391#endif 1392 1393 call qenter('cc_r12xi') 1394 1395 if (locdbg) then 1396 write(lupri,*) 'Entered CC_R12XI' 1397 write(lupri,*) 'memory available: ',lwork,' double words' 1398 end if 1399 1400 !check symmetry: 1401 if (isymxi.ne.muld2h(isymc,isymv)) then 1402 call quit('Symmetry error in CC_R12XI') 1403 end if 1404 1405C----------------------------------------------------------------------- 1406C allocate memory for intermediates needed until end 1407C----------------------------------------------------------------------- 1408 kvcint = 1 !VC-intermediate 1409 kend1 = kvcint + ntr12sq(isymxi) 1410 lwrk1 = lwork - kend1 1411 if (lwrk1 .lt. 0) then 1412 call quit('Insufficient memory in CC_R12XI') 1413 end if 1414 1415C----------------------------------------------------------------------- 1416C calculate Sum_m (c_{kl}^{im}*V_m^jtilde + c_{kl}^{mj}*V_m^itilde) 1417C or Sum_m (c_{kl}^{im}*V_j^mtilde + c_{kl}^{mj}*V_j^mtilde) 1418C and apply P_kl^ij: 1419C----------------------------------------------------------------------- 1420 call cc_r12xi2a(work(kvcint),isymxi,TR12,ISYMC,PRPAO,ISYMV, 1421 & CMO,LAMH,1,TRV,WORK(KEND1),LWRK1) 1422 1423C----------------------------------------------------------------------- 1424C finally calculate Xi 1425C----------------------------------------------------------------------- 1426 !zero out Xi: 1427 call dzero(xir12,ntr12sq(isymxi)) 1428 1429 call cc_r12xi2b(xir12,trxi,vxint,isymv,'N', 1430 & tr12,isymc,'N',one) 1431 call cc_r12xi2b(xir12,trxi,xint,1,'N', 1432 & work(kvcint),isymxi,'N',-one) 1433 1434 if (locdbg) then 1435 call around('Result in CC_R12XI') 1436 call cc_prsqr12(xir12,isymxi,'N',1,.false.) 1437 end if 1438 1439 if (locdbg) write(lupri,*) 'Leaving CC_R12XI' 1440 call qexit('cc_r12xi') 1441 return 1442 end 1443 1444*======================================================================= 1445 1446*=====================================================================* 1447 subroutine cc_r12xi2a(VCINT,ISYMVC,TR12,ISYMC,PRPAO,ISYMV, 1448 & CMO,LAMH,ISYMH,TRV,WORK,LWORK) 1449C---------------------------------------------------------------------- 1450C calculate Sum_m (c_{kl}^{im}*V_m^jtilde + c_{kl}^{mj}*V_m^itilde) 1451C or Sum_m (c_{kl}^{im}*V_j^mtilde + c_{kl}^{mj}*V_j^mtilde) 1452C and apply P_kl^ij 1453C 1454C VCINT result array of dimension ntr12sq(isymvc) 1455C ISYMVC symmetry of result 1456C TR12 R12 amplitudes of dimension ntr12sq(isymc) 1457C ISYMC symmetry of amplitudes 1458C PRPAO one-electron perturbation operator integrals in AO-basis 1459C ISYMV symmetry of perturbation 1460C LAMH Lambda^h matrix 1461C ISYMH Symmetry of Lambda^h matrix 1462C CMO MO coefficient matrix 1463C TRV use normal ('N') or transposed ('T') matrix V_m^jtilde 1464C when constructing vc-intermediate 1465C 1466C Christian Neiss 2005 1467C---------------------------------------------------------------------- 1468 implicit none 1469#include "priunit.h" 1470#include "ccorb.h" 1471#include "ccsdsym.h" 1472#include "r12int.h" 1473#include "dummy.h" 1474 1475 logical locdbg 1476 parameter(locdbg = .false.) 1477 1478#if defined (SYS_CRAY) 1479 real zero,one 1480#else 1481 double precision zero,one 1482#endif 1483 parameter (zero = 0.0d0, one = 1.0d0) 1484 1485 integer isymvc,isymv,isymij,isymkl,isymc,isymh,isym1,isym2, 1486 & isymi,isymj,isymm,isymim,isymmj 1487 integer kscr1,kvint,koffc,koffv, 1488 & koffvc,kofftr12 1489 integer lwork,kend1,kend2,lwrk1,lwrk2 1490 character*1 trv !allowed are 'N' and 'T' 1491 1492#if defined (SYS_CRAY) 1493 real vcint(*),tr12(*),prpao(*),lamh(*),cmo(*),work(*),ddot 1494#else 1495 double precision vcint(*),tr12(*),prpao(*),lamh(*),cmo(*),work(*), 1496 & ddot 1497#endif 1498 1499 call qenter('cc_r12xi2a') 1500 1501 if (locdbg) then 1502 write(lupri,*) 'Entered CC_R12XI2A' 1503 write(lupri,*) 'memory available: ',lwork,' double words' 1504 end if 1505 1506 !check variable trv: 1507 if (.not.((trv.eq.'T').or.(trv.eq.'N'))) then 1508 write(lupri,*) 'TRV = ',trv 1509 call quit('Forbidden value of TRV in CC_R12XI2A') 1510 end if 1511 1512 !check symmetry: 1513 if (isymvc.ne.muld2h(isymc,muld2h(isymh,isymv))) then 1514 call quit('Symmetry error in CC_R12XI2A') 1515 end if 1516 1517C----------------------------------------------------------------------- 1518C print out R12 amplitudes (debug) 1519C----------------------------------------------------------------------- 1520 if (locdbg) then 1521 call around('R12 amplitudes in CC_R12XI after unpacking') 1522 call cc_prsqr12(tr12,isymc,'N',1,.false.) 1523 write(lupri,*) 'Norm^2 of R12 amplitudes: ', 1524 & ddot(ntr12sq(isymc),tr12,1,tr12,1) 1525 end if 1526 1527C----------------------------------------------------------------------- 1528C transform perturbation integrals in MO-basis (occ. orbitals only) 1529C----------------------------------------------------------------------- 1530 kvint = 1 1531 kend1 = kvint + nmatij(muld2h(isymh,isymv)) 1532 1533 kscr1 = kend1 1534 kend2 = kscr1 + nt1ao(isymv) 1535 lwrk2 = lwork - kend2 1536 if (lwrk2 .lt. 0) then 1537 call quit('Insufficient memory in CC_R12XI2A') 1538 end if 1539 1540 do isym2 = 1, nsym 1541 isym1 = muld2h(isymv,isym2) 1542 koffv = 1 + iaodis(isym1,isym2) 1543 koffc = 1 + iglmrh(isym1,isym1) 1544 !first transformation 1545 call dgemm('T','N',nrhf(isym1),nbas(isym2),nbas(isym1), 1546 & one,cmo(koffc),max(1,nbas(isym1)),prpao(koffv), 1547 & max(1,nbas(isym1)),zero,work(kscr1), 1548 & max(1,nrhf(isym1))) 1549C call output(work(kscr1),1,nrhf(isym1),1,nbas(isym2),nrhf(isym1), 1550C & nbas(isym2),1,lupri) 1551 !second transformation 1552 koffv = kvint + imatij(isym1,muld2h(isymh,isym2)) 1553 koffc = 1 + iglmrh(isym2,muld2h(isymh,isym2)) 1554 call dgemm('N','N',nrhf(isym1),nrhf(muld2h(isymh,isym2)), 1555 & nbas(isym2),one,work(kscr1),max(1,nrhf(isym1)), 1556 & lamh(koffc),max(1,nbas(isym2)),zero,work(koffv), 1557 & max(1,nrhf(isym1))) 1558 end do 1559 1560 1561 if (locdbg) then 1562 call around('Original V integrals:') 1563 do isym2 = 1, nsym 1564 isym1 = muld2h(isymv,isym2) 1565 koffv = 1 + iaodis(isym1,isym2) 1566 write(lupri,*) 'Symmetry block: ',isym1, isym2 1567 if (nbas(isym1).eq.0 .or. nbas(isym2).eq.0) then 1568 write(lupri,*) 'This symmetry is empty' 1569 else 1570 call output(prpao(koffv),1,nbas(isym1),1, 1571 & nbas(isym2),nbas(isym1),nbas(isym2),1,lupri) 1572 end if 1573 end do 1574 call around('CMO matrix (occ. part):') 1575 do isym1 = 1, nsym 1576 koffc = 1 + iglmrh(isym1,isym1) 1577 write(lupri,*) 'Symmetry block: ',isym1, isym1 1578 if (nbas(isym1).eq.0 .or. nrhf(isym1).eq.0) then 1579 write(lupri,*) 'This symmetry is empty' 1580 else 1581 call output(cmo(koffc),1,nbas(isym1),1,nrhf(isym1), 1582 & nbas(isym1),nrhf(isym1),1,lupri) 1583 end if 1584 end do 1585 call around('Lambda^h matrix (occ. part):') 1586 do isym2 = 1, nsym 1587 isym1 = muld2h(isymh,isym2) 1588 koffc = 1 + iglmrh(isym1,isym2) 1589 write(lupri,*) 'Symmetry block: ',isym1, isym2 1590 if (nbas(isym1).eq.0 .or. nrhf(isym2).eq.0) then 1591 write(lupri,*) 'This symmetry is empty' 1592 else 1593 call output(lamh(koffc),1,nbas(isym1),1,nrhf(isym2), 1594 & nbas(isym1),nrhf(isym2),1,lupri) 1595 end if 1596 end do 1597 call around('Transformed V integrals:') 1598 do isym1 = 1, nsym 1599 isym2 = muld2h(muld2h(isymh,isymv),isym1) 1600 koffv = kvint + imatij(isym1,isym2) 1601 write(lupri,*) 'Symmetry block: ',isym1, isym2 1602 if (nrhf(isym1).eq.0 .or. nrhf(isym2).eq.0) then 1603 write(lupri,*) 'This symmetry is empty' 1604 else 1605 call output(work(koffv),1,nrhf(isym1),1,nrhf(isym2), 1606 & nrhf(isym1),nrhf(isym2),1,lupri) 1607 end if 1608 end do 1609 write(lupri,*) 'Norm^2 of Transformed V integrals: ', 1610 & ddot(nmatij(muld2h(isymh,isymv)),work(kvint),1,work(kvint),1) 1611 call flshfo(lupri) 1612 end if 1613 1614C----------------------------------------------------------------------- 1615C calculate Sum_m (c_{kl}^{im}*V_m^jtilde + c_{kl}^{mj}*V_m^itilde) 1616C or Sum_m (c_{kl}^{im}*V_j^mtilde + c_{kl}^{mj}*V_j^mtilde) 1617C----------------------------------------------------------------------- 1618 do isymkl = 1, nsym 1619 isymim = muld2h(isymkl,isymc) 1620 isymmj = muld2h(isymh,isymv) 1621 isymij = muld2h(isymim,isymmj) 1622 do isymm = 1, nsym 1623 isymi = muld2h(isymim,isymm) 1624 isymj = muld2h(isymmj,isymm) 1625 kofftr12 = 1 + itr12sq(isymkl,isymim) + 1626 & nmatkl(isymkl)*imatij(isymi,isymm) 1627c koffv = kvint + imatij(isymm,isymj) 1628 koffvc = 1 + itr12sq(isymkl,isymij) + 1629 & nmatkl(isymkl)*imatij(isymi,isymj) 1630c write(lupri,*) 'symmetry isymkl, isymi, isymm, isymj: ', 1631c & isymkl,isymi,isymm,isymj 1632c write(lupri,*) 'itr12sq(isymkl,isymij): ', 1633c & itr12sq(isymkl,isymij) 1634c write(lupri,*) 'nmatkl(isymkl), nrhf(isymi), nrhf(isymm)', 1635c & nmatkl(isymkl), nrhf(isymi), nrhf(isymm) 1636c call output(tr12(kofftr12),1,nmatkl(isymkl)*nrhf(isymi),1, 1637c & nrhf(isymm),nmatkl(isymkl)*nrhf(isymi), 1638c & nrhf(isymm),1,lupri) 1639 if (trv.eq.'N') then 1640 koffv = kvint + imatij(isymm,isymj) 1641 call dgemm('N','N',nmatkl(isymkl)*nrhf(isymi),nrhf(isymj), 1642 & nrhf(isymm),one,tr12(kofftr12), 1643 & max(1,nmatkl(isymkl)*nrhf(isymi)), 1644 & work(koffv),max(1,nrhf(isymm)),zero, 1645 & vcint(koffvc),max(1,nmatkl(isymkl)*nrhf(isymi))) 1646 else if (trv.eq.'T') then 1647 koffv = kvint + imatij(isymj,isymm) 1648 call dgemm('N','T',nmatkl(isymkl)*nrhf(isymi),nrhf(isymj), 1649 & nrhf(isymm),one,tr12(kofftr12), 1650 & max(1,nmatkl(isymkl)*nrhf(isymi)), 1651 & work(koffv),max(1,nrhf(isymj)),zero, 1652 & vcint(koffvc),max(1,nmatkl(isymkl)*nrhf(isymi))) 1653 else 1654 call quit('Forbidden value of TRV in CC_R12XI2A') 1655 end if 1656 end do 1657 end do 1658 1659 if (locdbg) then 1660 call around('"c*V" before P_{kl}^{ij}') 1661 do isymkl = 1, nsym 1662 isymij = muld2h(isymkl,isymvc) 1663 koffvc = 1 + itr12sq(isymkl,isymij) 1664 write(lupri,*) 'Symmetry block isymkl,isymij: ',isymkl,isymij 1665 if (nmatkl(isymkl).eq.0 .or. nmatij(isymij).eq.0) then 1666 write(lupri,*) 'This symmetry is empty' 1667 else 1668 call output(vcint(koffvc),1,nmatkl(isymkl),1,nmatij(isymij), 1669 & nmatkl(isymkl),nmatij(isymij),1,lupri) 1670 end if 1671 end do 1672C write(lupri,*) (vcint(i), i=1, ntr12sq(isymvc)) 1673 write(lupri,*) 'Norm^2 of c*V before P_{kl}^{ij}: ', 1674 & ddot(ntr12sq(isymvc),vcint,1,vcint,1) 1675 end if 1676 1677 !apply P_{kl}^{ij} 1678 call cc_r12pklij(vcint,isymvc,'N',work,lwork) 1679 1680 if (locdbg) then 1681 call around('VC-intermediate') 1682 do isymkl = 1, nsym 1683 isymij = muld2h(isymkl,isymvc) 1684 koffvc = 1 + itr12sq(isymkl,isymij) 1685 write(lupri,*) 'Symmetry block isymkl,isymij: ',isymkl,isymij 1686 if (nmatkl(isymkl).eq.0 .or. nmatij(isymij).eq.0) then 1687 write(lupri,*) 'This symmetry is empty' 1688 else 1689 call output(vcint(koffvc),1,nmatkl(isymkl),1,nmatij(isymij), 1690 & nmatkl(isymkl),nmatij(isymij),1,lupri) 1691 end if 1692 end do 1693 write(lupri,*) 'Norm^2 of VC-intermediate ', 1694 & ddot(ntr12sq(isymvc),vcint,1,vcint,1) 1695 end if 1696 1697 if (locdbg) write(lupri,*) 'Leaving CC_R12XI2A' 1698 call qexit('cc_r12xi2a') 1699 return 1700 end 1701 1702*======================================================================= 1703 1704*=====================================================================* 1705 subroutine cc_r12xi2b(RESULT,TRR,MAT1,ISYM1,TR1,MAT2,ISYM2,TR2, 1706 & FAC) 1707C---------------------------------------------------------------------- 1708C purpose: calculate FAC*Sum_kl X_{mn}^{kl}*c_{kl}^{ij} and add it 1709C to a given array (RESULT) 1710C 1711C RESULT result array of dimension ntr12sq(muld2h(isym1,isym2)) 1712C TRR R12-index pair (mn) leading: 'N' 1713C occ. index pair (ij) leading: 'T' 1714C MAT1 Matrix 1 [X] dim: nr12r12sq(isym1) 1715C ISYM1 Symmetry of MAT1 1716C MAT2 Matrix 2 [c] dim: ntr12sq(isym2) 1717C ISYM2 Symmetry of MAT2 1718C FAC Factor 1719C TR1 index pair (mn) leading: 'N' 1720C index pair (kl) leading: 'T' 1721C TR2 index pair (kl) leading: 'N' 1722C index pair (ij) leading: 'T' 1723C 1724C Christian Neiss 2005 1725C---------------------------------------------------------------------- 1726 implicit none 1727#include "priunit.h" 1728#include "ccorb.h" 1729#include "ccsdsym.h" 1730 1731 logical locdbg 1732 parameter(locdbg = .false.) 1733 1734#if defined (SYS_CRAY) 1735 real zero,one 1736#else 1737 double precision zero,one 1738#endif 1739 parameter (zero = 0.0d0, one = 1.0d0) 1740 1741 integer isym1,isym2,isymmn,isymij,isymkl 1742 integer koff1,koff2,koffr 1743 character*1 tr1,tr2,trr 1744 1745#if defined (SYS_CRAY) 1746 real result(*),mat1(*),mat2(*),fac,ddot 1747#else 1748 double precision result(*),mat1(*),mat2(*),fac,ddot 1749#endif 1750 1751 call qenter('cc_r12xi2b') 1752 1753 if (locdbg) then 1754 write(lupri,*) 'Entered CC_R12XI2B' 1755 end if 1756 1757 !check TR1, TR2, TRR 1758 IF (((TR1.NE.'N').AND.(TR1.NE.'T')).OR. 1759 & ((TR2.NE.'N').AND.(TR2.NE.'T')).OR. 1760 & ((TRR.NE.'N').AND.(TRR.NE.'T'))) THEN 1761 WRITE(LUPRI,*) 'TR1, TR2: ',TR1,TR2 1762 CALL QUIT('Illegal value for TR1 or TR2 in CC_R12XI2B') 1763 END IF 1764 1765 if (locdbg) then 1766 call around('MAT1') 1767 write(lupri,*) 'Norm^2: ',DDOT(nr12r12sq(isym1),mat1,1,mat1,1) 1768C do isymmn = 1, nsym 1769C isymkl = muld2h(isymmn,isym1) 1770C koff1 = 1 + ir12r12sq(isymmn,isymkl) 1771C write(lupri,*) 'Symmetry block ',isymmn,isymkl 1772C if (nmatkl(isymmn).eq.0 .or. nmatkl(isymkl).eq.0) then 1773C write(lupri,*) 'This symmetry is empty' 1774C else 1775C call output(mat1(koff1),1,nmatkl(isymmn),1,nmatkl(isymkl), 1776C & nmatkl(isymmn),nmatkl(isymkl),1,lupri) 1777C end if 1778C end do 1779 call around('MAT2') 1780 write(lupri,*) 'Norm^2: ',DDOT(ntr12sq(isym2),mat2,1,mat2,1) 1781C if (tr2.eq.'N') then 1782C do isymkl = 1, nsym 1783C isymij = muld2h(isymkl,isym2) 1784C koff2 = 1 + itr12sq(isymkl,isymij) 1785C write(lupri,*) 'Symmetry block ',isymkl,isymij 1786C if (nmatkl(isymkl).eq.0 .or. nmatij(isymij).eq.0) then 1787C write(lupri,*) 'This symmetry is empty' 1788C else 1789C call output(mat2(koff2),1,nmatkl(isymkl),1,nmatij(isymij), 1790C & nmatkl(isymkl),nmatij(isymij),1,lupri) 1791C end if 1792C end do 1793C else if (tr2.eq.'T') then 1794C do isymij = 1, nsym 1795C isymkl = muld2h(isymij,isym2) 1796C koff2 = 1 + itr12sqt(isymij,isymkl) 1797C write(lupri,*) 'Symmetry block ',isymij,isymkl 1798C if (nmatkl(isymkl).eq.0 .or. nmatij(isymij).eq.0) then 1799C write(lupri,*) 'This symmetry is empty' 1800C else 1801C call output(mat2(koff2),1,nmatij(isymij),1,nmatkl(isymkl), 1802C & nmatij(isymij),nmatkl(isymkl),1,lupri) 1803C end if 1804C end do 1805C end if 1806 write(lupri,*) 'TR1, TR2: ',TR1,TR2 1807 call flshfo(lupri) 1808 end if 1809 1810 if (TRR.eq.'N') then 1811 do isymmn = 1, nsym 1812 isymkl = muld2h(isymmn,isym1) 1813 isymij = muld2h(isymkl,isym2) 1814 if ((tr1.eq.'N').and.(tr2.eq.'N')) then 1815 koff1 = 1 + ir12r12sq(isymmn,isymkl) 1816 koff2 = 1 + itr12sq(isymkl,isymij) 1817 koffr = 1 + itr12sq(isymmn,isymij) 1818 call dgemm(tr1,tr2,nmatkl(isymmn),nmatij(isymij), 1819 & nmatkl(isymkl),fac,mat1(koff1),max(1,nmatkl(isymmn)), 1820 & mat2(koff2),max(1,nmatkl(isymkl)),one, 1821 & result(koffr),max(1,nmatkl(isymmn))) 1822 else if ((tr1.eq.'T').and.(tr2.eq.'T')) then 1823 koff1 = 1 + ir12r12sq(isymkl,isymmn) 1824 koff2 = 1 + itr12sqt(isymij,isymkl) 1825 koffr = 1 + itr12sq(isymmn,isymij) 1826 call dgemm(tr1,tr2,nmatkl(isymmn),nmatij(isymij), 1827 & nmatkl(isymkl),fac,mat1(koff1),max(1,nmatkl(isymkl)), 1828 & mat2(koff2),max(1,nmatij(isymij)),one, 1829 & result(koffr),max(1,nmatkl(isymmn))) 1830 else if ((tr1.eq.'T').and.(tr2.eq.'N')) then 1831 koff1 = 1 + ir12r12sq(isymkl,isymmn) 1832 koff2 = 1 + itr12sq(isymkl,isymij) 1833 koffr = 1 + itr12sq(isymmn,isymij) 1834 call dgemm(tr1,tr2,nmatkl(isymmn),nmatij(isymij), 1835 & nmatkl(isymkl),fac,mat1(koff1),max(1,nmatkl(isymkl)), 1836 & mat2(koff2),max(1,nmatkl(isymkl)),one, 1837 & result(koffr),max(1,nmatkl(isymmn))) 1838 else if ((tr1.eq.'N').and.(tr2.eq.'T')) then 1839 koff1 = 1 + ir12r12sq(isymmn,isymkl) 1840 koff2 = 1 + itr12sqt(isymij,isymkl) 1841 koffr = 1 + itr12sq(isymmn,isymij) 1842 call dgemm(tr1,tr2,nmatkl(isymmn),nmatij(isymij), 1843 & nmatkl(isymkl),fac,mat1(koff1),max(1,nmatkl(isymmn)), 1844 & mat2(koff2),max(1,nmatij(isymij)),one, 1845 & result(koffr),max(1,nmatkl(isymmn))) 1846 end if 1847 end do 1848 else if (TRR.eq.'T') then 1849 do isymmn = 1, nsym 1850 isymkl = muld2h(isymmn,isym1) 1851 isymij = muld2h(isymkl,isym2) 1852 if ((tr1.eq.'N').and.(tr2.eq.'N')) then 1853 koff1 = 1 + ir12r12sq(isymmn,isymkl) 1854 koff2 = 1 + itr12sq(isymkl,isymij) 1855 koffr = 1 + itr12sqt(isymij,isymmn) 1856 call dgemm('T','T',nmatij(isymij),nmatkl(isymmn), 1857 & nmatkl(isymkl),fac,mat2(koff2),max(1,nmatkl(isymkl)), 1858 & mat1(koff1),max(1,nmatkl(isymmn)),one, 1859 & result(koffr),max(1,nmatij(isymij))) 1860 else if ((tr1.eq.'T').and.(tr2.eq.'T')) then 1861 koff1 = 1 + ir12r12sq(isymkl,isymmn) 1862 koff2 = 1 + itr12sqt(isymij,isymkl) 1863 koffr = 1 + itr12sqt(isymij,isymmn) 1864 call dgemm('N','N',nmatij(isymij),nmatkl(isymmn), 1865 & nmatkl(isymkl),fac,mat2(koff2),max(1,nmatij(isymij)), 1866 & mat1(koff1),max(1,nmatkl(isymkl)),one, 1867 & result(koffr),max(1,nmatij(isymij))) 1868 else if ((tr1.eq.'T').and.(tr2.eq.'N')) then 1869 koff1 = 1 + ir12r12sq(isymkl,isymmn) 1870 koff2 = 1 + itr12sq(isymkl,isymij) 1871 koffr = 1 + itr12sqt(isymij,isymmn) 1872 call dgemm('T','N',nmatij(isymij),nmatkl(isymmn), 1873 & nmatkl(isymkl),fac,mat2(koff2),max(1,nmatkl(isymkl)), 1874 & mat1(koff1),max(1,nmatkl(isymkl)),one, 1875 & result(koffr),max(1,nmatij(isymij))) 1876 else if ((tr1.eq.'N').and.(tr2.eq.'T')) then 1877 koff1 = 1 + ir12r12sq(isymmn,isymkl) 1878 koff2 = 1 + itr12sqt(isymij,isymkl) 1879 koffr = 1 + itr12sqt(isymij,isymmn) 1880 call dgemm('N','T',nmatij(isymij),nmatkl(isymmn), 1881 & nmatkl(isymkl),fac,mat2(koff2),max(1,nmatij(isymij)), 1882 & mat1(koff1),max(1,nmatkl(isymmn)),one, 1883 & result(koffr),max(1,nmatij(isymij))) 1884 end if 1885 end do 1886 else 1887 call quit('Unknown value for TRR in CC_R12XI2B') 1888 end if 1889 1890 if (locdbg) then 1891 call around('Result of CC_R12XI2B') 1892 write(lupri,*) 'Norm^2: ',DDOT(ntr12sq(muld2h(isym1,isym2)), 1893 & result,1,result,1) 1894C if (TRR.eq.'N') then 1895C do isymmn = 1, nsym 1896C isymij = muld2h(isymmn,muld2h(isym1,isym2)) 1897C koffr = 1 + itr12sq(isymmn,isymij) 1898C write(lupri,*) 'Symmetry block ',isymmn,isymij 1899C if (nmatkl(isymmn).eq.0 .or. nmatij(isymij).eq.0) then 1900C write(lupri,*) 'This symmetry is empty' 1901C else 1902C call output(result(koffr),1,nmatkl(isymmn),1,nmatij(isymij), 1903C & nmatkl(isymmn),nmatij(isymij),1,lupri) 1904C end if 1905C end do 1906C else if (TRR.eq.'T') then 1907C do isymij = 1, nsym 1908C isymmn = muld2h(isymij,muld2h(isym1,isym2)) 1909C koffr = 1 + itr12sqt(isymij,isymmn) 1910C write(lupri,*) 'Symmetry block ',isymij,isymmn 1911C if (nmatkl(isymmn).eq.0 .or. nmatij(isymij).eq.0) then 1912C write(lupri,*) 'This symmetry is empty' 1913C else 1914C call output(result(koffr),1,nmatij(isymij),1,nmatkl(isymmn), 1915C & nmatij(isymij),nmatkl(isymmn),1,lupri) 1916C end if 1917C end do 1918C end if 1919 end if 1920 1921 if (locdbg) write(lupri,*) 'Leaving CC_R12XI2B' 1922 call qexit('cc_r12xi2b') 1923 return 1924 end 1925 1926*======================================================================= 1927 1928*=====================================================================* 1929 subroutine cc_r12rdvxint(matrix,work,lwork,ff,isympt,labpt) 1930C---------------------------------------------------------------------- 1931C purpose: read in VXINT, calculate FF*VXINT(isympt) and 1932C add it to a given array (MATRIX) 1933C 1934C MATRIX array of dimension nr12r12sq(isympt) 1935C ISYMPT symmetry of perturbation 1936C LABPT label of perturbation 1937C FF Factor (Fieldstrength) 1938C 1939C Christian Neiss, Feb. 2005, based on CC_ONEP 1940C---------------------------------------------------------------------- 1941 implicit none 1942#include "priunit.h" 1943#include "ccorb.h" 1944#include "ccsdsym.h" 1945#include "dummy.h" 1946 1947 logical locdbg 1948 parameter(locdbg = .false.) 1949 1950#if defined (SYS_CRAY) 1951 real zero,one 1952#else 1953 double precision zero,one 1954#endif 1955 parameter (zero = 0.0d0, one = 1.0d0) 1956 1957 character*8 labpt 1958 integer isympt,isymv 1959 integer kvxint,kvxintsq,luvxint 1960 integer lwork,kend1,lwrk1,idum,iopt 1961 1962#if defined (SYS_CRAY) 1963 real matrix(*),work(*),ff,ddot 1964#else 1965 double precision matrix(*),work(*),ff,ddot 1966#endif 1967 1968 call qenter('cc_r12rdvxint') 1969 1970 if (locdbg) then 1971 write(lupri,*) 'Entered CC_R12RDVXINT' 1972 end if 1973 1974 kvxint = 1 1975 kvxintsq = kvxint + nr12r12p(isympt) 1976 kend1 = kvxintsq + nr12r12sq(isympt) 1977 lwrk1 = lwork - kend1 1978 if (lwrk1 .lt. 0) then 1979 call quit('Insufficient memory in CC_R12RDVXINT') 1980 end if 1981 1982C----------------------------------------------------------------------- 1983C read in VXINT 1984C----------------------------------------------------------------------- 1985 luvxint = -1 1986 call gpopen(luvxint,'CCR12VXINT','OLD',' ','UNFORMATTED', 1987 & idum,.false.) 1988 rewind(luvxint) 1989 call mollab(labpt,luvxint,lupri) 1990 read(luvxint) isymv 1991 if (isymv .ne. isympt) then 1992 write(lupri,*) 'LABEL, ISYMV, ISYMPT: ',LABPT, ISYMV, ISYMPT 1993 call quit('Symmetry mismatch when reading CCR12VXINT in '// 1994 & 'CC_R12RDVXINT') 1995 end if 1996 read(luvxint) (work(kvxint+i-1),i=1,nr12r12p(isympt)) 1997 call gpclose(luvxint,'KEEP') 1998 !unpack to square: 1999 iopt = 2 2000 call ccr12unpck2(work(kvxint),isympt,work(kvxintsq),'N',iopt) 2001 2002 if (locdbg) then 2003 write(lupri,*) 'LABEL = ',labpt 2004 call around('VXINT in CC_R12RDVXINT') 2005 write(lupri,*) 'Norm^2 (packed): ',ddot(ntr12sq(isympt), 2006 & work(kvxint),1,work(kvxint),1) 2007 write(lupri,*) 'Norm^2 (squared): ',ddot(nr12r12sq(isympt), 2008 & work(kvxintsq),1,work(kvxintsq),1) 2009C call cc_prsqr12(work(kvxintsq),isympt,1,.false.) 2010 end if 2011 2012C----------------------------------------------------------------------- 2013C multiply with FF and add to MATRIX 2014C----------------------------------------------------------------------- 2015 call daxpy(nr12r12sq(isympt),ff,work(kvxintsq),1,matrix,1) 2016 2017 if (locdbg) write(lupri,*) 'Leaving CC_R12RDVXINT' 2018 call qexit('cc_r12rdvxint') 2019 return 2020 end 2021 2022*======================================================================= 2023 2024*=====================================================================* 2025 subroutine cc_r12pklij(MATRIX,ISYM,TRANS,WORK,LWORK) 2026C---------------------------------------------------------------------- 2027C purpose: calculate P_{kl}^{ij}X_{kl}^{ij} 2028C = X_{kl}^{ij} + X_{lk}^{ji} 2029C 2030C MATRIX Matrix of dimension ntr12sq(ISYM) with {kl} and {ij} as 2031C packed pair indices; MATRIX is input and output! 2032C 2033C Christian Neiss, Feb. 2005 2034C---------------------------------------------------------------------- 2035 implicit none 2036#include "priunit.h" 2037#include "ccorb.h" 2038#include "ccsdsym.h" 2039 2040 logical locdbg 2041 parameter(locdbg = .false.) 2042 2043#if defined (SYS_CRAY) 2044 real zero,one 2045#else 2046 double precision zero,one 2047#endif 2048 parameter (zero = 0.0d0, one = 1.0d0) 2049 2050 character*1 trans 2051 integer isym,isymkl,isymij,isymk,isyml,isymi,isymj 2052 integer idxkl,idxlk,idxij,idxji,idxklij,idxlkji 2053 integer kscr1,kscr2,koff 2054 integer lwork,kend1,lwrk1,iopt 2055 2056#if defined (SYS_CRAY) 2057 real work(*),matrix(*),dnrm2 2058#else 2059 double precision work(*),matrix(*),dnrm2 2060#endif 2061 2062 2063 call qenter('cc_r12pklij') 2064 2065 if (locdbg) then 2066 write(lupri,*) 'Entered CC_R12PKLIJ' 2067 write(lupri,*) 'memory available: ',lwork,' double words' 2068 write(lupri,*) 'Input in CC_R12PKLIJ:' 2069 call cc_prsqr12(matrix,isym,trans,1,.false.) 2070 end if 2071 2072 !only for testing needed, see below! 2073 if (locdbg) then 2074 kscr1 = 1 2075 kscr2 = kscr1 + ntr12am(isym) 2076 kend1 = kscr2 + ntr12sq(isym) 2077 lwrk1 = lwork - kend1 2078 if (lwrk1 .lt. 0) then 2079 call quit('Insufficient memory in CC_R12PKLIJ') 2080 end if 2081 2082 call dcopy(ntr12sq(isym),matrix,1,work(kscr2),1) 2083 end if 2084 2085C calculate X_{lk}^{ji} 2086 do isymkl = 1, nsym 2087 isymij = muld2h(isym,isymkl) 2088 do isymk = 1, nsym 2089 isyml = muld2h(isymkl,isymk) 2090 do isymi = 1, nsym 2091 isymj = muld2h(isymij,isymi) 2092 do k = 1, nrhfb(isymk) 2093 do l = 1, nrhfb(isyml) 2094 idxkl = imatkl(isymk,isyml) + nrhfb(isymk)*(l-1) + k 2095 idxlk = imatkl(isyml,isymk) + nrhfb(isyml)*(k-1) + l 2096 do i = 1, nrhf(isymi) 2097 do j = 1, nrhf(isymj) 2098 idxij = imatij(isymi,isymj) + nrhf(isymi)*(j-1)+i 2099 idxji = imatij(isymj,isymi) + nrhf(isymj)*(i-1)+j 2100 if (trans.eq.'T') then 2101 idxklij = itr12sqt(isymij,isymkl) + 2102 & nmatij(isymij)*(idxkl-1) + idxij 2103 idxlkji = itr12sqt(isymij,isymkl) + 2104 & nmatij(isymij)*(idxlk-1) + idxji 2105 else if (trans.eq.'N') then 2106 idxklij = itr12sq(isymkl,isymij) + 2107 & nmatkl(isymkl)*(idxij-1) + idxkl 2108 idxlkji = itr12sq(isymkl,isymij) + 2109 & nmatkl(isymkl)*(idxji-1) + idxlk 2110 else 2111 call quit('Illeagl value for "TRANS" in '// 2112 & 'CC_R12PKLIJ') 2113 end if 2114 if (idxklij.le.idxlkji) then 2115 matrix(idxklij)=matrix(idxklij)+matrix(idxlkji) 2116 matrix(idxlkji)=matrix(idxklij) 2117 end if 2118 end do 2119 end do 2120 end do 2121 end do 2122 end do 2123 end do 2124 end do 2125 2126 if (locdbg) then 2127 call around('Result in CC_R12PKLIJ') 2128 call cc_prsqr12(matrix,isym,trans,1,.false.) 2129 end if 2130 2131 !Test: calling these two routines should give the same as before! 2132 if (locdbg) then 2133 iopt = 1 2134 call ccr12pck2(work(kscr1),isym,.TRUE.,work(kscr2),trans, 2135 & iopt) 2136 call ccr12unpck2(work(kscr1),isym,work(kscr2),trans,iopt) 2137 call around('alternatively calculated result in CC_R12PKLIJ') 2138 call cc_prsqr12(work(kscr2),isym,trans,1,.false.) 2139 2140 call daxpy(ntr12sq(isym),-one,matrix,1,work(kscr2),1) 2141 write(lupri,*) 'Norm^2 of difference: ', 2142 & dnrm2(ntr12sq(isym),work(kscr2),1) 2143 call flshfo(lupri) 2144 end if 2145 2146 if (locdbg) write(lupri,*) 'Leaving CC_R12PKLIJ' 2147 call qexit('cc_r12pklij') 2148 return 2149 end 2150 2151*======================================================================= 2152 2153*=====================================================================* 2154 subroutine cc_r12tstform(WORK,LWORK) 2155C---------------------------------------------------------------------- 2156C purpose: test the routines for reordering of amplitudes etc. 2157C pure test routine! 2158C 2159C Christian Neiss, Feb. 2005 2160C---------------------------------------------------------------------- 2161 implicit none 2162#include "priunit.h" 2163#include "ccorb.h" 2164#include "ccsdsym.h" 2165#include "dummy.h" 2166#include "r12int.h" 2167#include "ccr12int.h" 2168 2169#if defined (SYS_CRAY) 2170 real zero,one,work(*) 2171#else 2172 double precision zero,one,work(*) 2173#endif 2174 parameter (zero = 0.0d0, one = 1.0d0) 2175 2176 integer lwork, lwrk0, kend0 2177 integer KTR12SQ,KTR12PK1,KTR12PK2,KTAMP12S,KTAMP12T 2178 integer NRHFTRIA 2179 integer lunit,idum,ij,n2,nkilj(8) 2180 logical ldum 2181 2182 2183 call qenter('cc_r12tstform') 2184 2185 2186 nrhftria = nrhftb*(nrhftb+1)/2 2187 n2 = nrhftria*nrhftria 2188 2189 KTR12SQ = 1 2190 KEND0 = KTR12SQ + NTR12SQ(1) 2191 2192 KTR12PK1 = KEND0 2193 KTR12PK2 = KTR12PK1 + NTR12AM(1) 2194 KEND0 = KTR12PK2 + NTR12AM(1) 2195 2196 KTAMP12S = KEND0 2197 KTAMP12T = KTAMP12S + NRHFTRIA*NRHFTRIA 2198 KEND0 = KTAMP12T + NRHFTRIA*NRHFTRIA 2199 2200 LWRK0 = LWORK - KEND0 2201 if (lwrk0.lt.0) then 2202 call quit('Not enough memory in CC_R12TSTFORM') 2203 end if 2204 2205 2206C method 1: 2207 CALL CC_R12GETCT(WORK(KTR12SQ),1,0,ketscl,.FALSE.,'T', 2208 & DUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 2209 & WORK(KEND0),LWRK0) 2210 call cc_prsqr12(WORK(KTR12SQ),1,'T',1,.true.) 2211 CALL CCR12PCK2(WORK(KTR12PK1),1,.FALSE.,WORK(KTR12SQ),'T',1) 2212 call cc_prpr12(WORK(KTR12PK1),1,1,.true.) 2213 2214 2215C method 2: 2216C CALL GPOPEN(lunit,fccr12c,'old',' ','formatted',idum,ldum) 2217C REWIND(LUNIT) 2218C READ(LUNIT,'(4E30.20)') (WORK(KTAMP12S+IJ), IJ = 0, N2-1) 2219C READ(LUNIT,'(4E30.20)') (WORK(KTAMP12T+IJ), IJ = 0, N2-1) 2220C CALL GPCLOSE(LUNIT,'KEEP') 2221C call around('R12 singlet part') 2222C call output(work(KTAMP12S),1,nrhftria,1,nrhftria, 2223C & nrhftria,nrhftria,1,lupri) 2224C call around('R12 triplet part') 2225C call output(work(KTAMP12T),1,nrhftria,1,nrhftria, 2226C & nrhftria,nrhftria,1,lupri) 2227C CALL CCR12PCK(WORK(KTR12PK2),1,WORK(KTAMP12S),WORK(KTAMP12T), 2228C & nrhfb,nrhf,nkilj) 2229C call cc_prpr12(WORK(KTR12PK2),1,1,.true.) 2230 2231C reconstruct singlet/triplet format from method 1: 2232 call ccr12unpck(WORK(KTR12PK1),1,WORK(KTAMP12S),WORK(KTAMP12T), 2233 & nrhfb,nrhf) 2234 call around('R12 singlet part') 2235 call output(work(KTAMP12S),1,nrhftria,1,nrhftria, 2236 & nrhftria,nrhftria,1,lupri) 2237 call around('R12 triplet part') 2238 call output(work(KTAMP12T),1,nrhftria,1,nrhftria, 2239 & nrhftria,nrhftria,1,lupri) 2240 2241C resonstruct square matrix format from method 2: 2242 call ccr12unpck2(WORK(KTR12PK2),1,WORK(KTR12SQ),'T',1) 2243 call cc_prsqr12(WORK(KTR12SQ),1,'T',1,.true.) 2244 2245C test P_kl^ij: 2246 ! Method 1: 2247 call cc_r12pklij(WORK(KTR12SQ),1,'T',work(kend0),lwrk0) 2248 call cc_prsqr12(WORK(KTR12SQ),1,'T',1,.true.) 2249 ! Method 2: 2250 call cc_r12vunpack(WORK(KTR12SQ),1,work(KTAMP12S),work(KTAMP12T), 2251 & .TRUE.,nrhfb,nrhf) 2252 call cc_prsqr12(WORK(KTR12SQ),1,'T',1,.true.) 2253 2254 2255 call quit('TEST DONE') 2256 call qexit('cc_r12tstform') 2257 return 2258 end 2259 2260*======================================================================= 2261 2262*======================================================================= 2263 subroutine cc_r12eta0(etar12sq,cmo,isyres,work,lwork) 2264C----------------------------------------------------------------------- 2265C purpose: r12 contribution to eta^(0) 2266C 2267C etar12sq r12 part of eta^(0) vector of dim. ntr12sq(isyres) 2268C (occ. index pair leading) 2269C cmo MO coefficient or lambda matrix 2270C isyres symmetry of result (for eta = 1) 2271C 2272C Note that the result is ADDED to the input! 2273C 2274C Christian Neiss Mar. 2005 2275C----------------------------------------------------------------------- 2276 2277 implicit none 2278#include "priunit.h" 2279#include "ccorb.h" 2280#include "ccsdsym.h" 2281#include "r12int.h" 2282#include "ccr12int.h" 2283#include "dummy.h" 2284 2285 logical locdbg 2286 parameter (locdbg = .false.) 2287 2288 integer isyres,isymak,isymk,isyma 2289 integer kvajkl,kvijkl,kend1,lwrk1 2290 integer lwork,luvajkl 2291 2292#if defined (SYS_CRAY) 2293 real cmo(*),etar12sq(*),WORK(LWORK),one,two,ddot 2294#else 2295 double precision cmo(*),etar12sq(*),WORK(LWORK),one,two,ddot 2296#endif 2297 2298 parameter (one = 1.0D0, two = 2.0D0) 2299 2300 call qenter('cc_r12eta0') 2301 if (locdbg) then 2302 write(lupri,*) 'Entered CC_R12ETA0' 2303 write(lupri,*) 'memory available: ',lwork,' double words' 2304 end if 2305 2306 kvajkl = 1 2307 kvijkl = kvajkl + nvajkl(1) 2308 kend1 = kvijkl + ntr12sq(isyres) 2309 lwrk1 = lwork - kend1 2310 if (lwrk1 .lt. 0) then 2311 call quit('Insufficient work space in cc_r12eta0') 2312 end if 2313 2314C----------------------------------------------------------------------- 2315C Read in V(alpha j,kl) 2316C----------------------------------------------------------------------- 2317 luvajkl = -1 2318 call gpopen(luvajkl,fvajkl,'old',' ','unformatted', 2319 & idummy,.false.) 2320 rewind(luvajkl) 2321 read(luvajkl) (work(kvajkl+i-1), i = 1,nvajkl(1)) 2322 call gpclose(luvajkl,'KEEP') 2323 if (locdbg) then 2324 write(lupri,*) 'fvajkl: ', fvajkl 2325 write(lupri,*) 'norm^2(V(alpha j,kl)):', 2326 & ddot(nvajkl(1),work(kvajkl),1,work(kvajkl),1) 2327 end if 2328 2329C----------------------------------------------------------------------- 2330C Calculate CMO*V 2331C----------------------------------------------------------------------- 2332 call dzero(work(kvijkl),ntr12sq(isyres)) 2333 call cc_r12mkvijkl(work(kvajkl),1,cmo,isyres,work(kend1), 2334 & lwrk1,.true.,one,work(kvijkl)) 2335 if (locdbg) then 2336 write(lupri,*) 'norm^2(vijkl) after cc_r12mkvijkl:', 2337 & ddot(ntr12sq(isyres),work(kvijkl),1,work(kvijkl),1) 2338 end if 2339 2340C----------------------------------------------------------------------- 2341C Make 2*Coulomb-Exchange: 2V(ij,kl)-V(ji,kl) 2342C----------------------------------------------------------------------- 2343 call cc_r12tcmesq(work(kvijkl),isyres,'T',.false.) 2344 2345 call daxpy(ntr12sq(isyres),two,work(kvijkl),1,etar12sq,1) 2346 2347 if (locdbg) write(lupri,*) 'Leaving CC_R12ETA0' 2348 call qexit('cc_r12eta0') 2349 return 2350 end 2351*======================================================================= 2352 2353*======================================================================= 2354 subroutine cc_r12prop(propr12,labelh,aproxr12,work,lwork) 2355C----------------------------------------------------------------------- 2356C purpose: r12 contribution to expectation value of property <label> 2357C 2358C propr12 r12 part of property 2359C labelh label of property 2360C aproxr12 r12 approximation 2361C 2362C Christian Neiss Mar. 2005 2363C----------------------------------------------------------------------- 2364 2365 implicit none 2366#include "priunit.h" 2367#include "ccorb.h" 2368#include "ccsdsym.h" 2369#include "ccsdinp.h" 2370#include "r12int.h" 2371#include "ccr12int.h" 2372#include "dummy.h" 2373 2374 logical locdbg 2375 parameter (locdbg = .false.) 2376 2377 character*8 labelh 2378 CHARACTER*3 APROXR12, LISTL, filxi 2379 CHARACTER*10 MODEL 2380 integer idxxi,idlstl,lenmod,iopt,isymxi,isymv,isym,idx 2381 integer kxir12,kzetar12,kzeta1,kzeta2,kzetar12sq,kxir12sq 2382 integer ian,luxint,ierr 2383 integer ktr12,ktr12sq,kxint,kxintsq,kvxint,kvxintsq, 2384 & kprpao,klamp0,klamh0,kt1amp0 2385 integer lwork,lwrk1,kend1,lwrk2,kend2 2386 2387 ! external function: 2388 integer IRHSR1,ILSTSYM 2389 2390#if defined (SYS_CRAY) 2391 real work(*),propr12,ddot,one,diacont 2392#else 2393 double precision work(*),propr12,ddot,one,diacont 2394#endif 2395 2396 parameter(one = 1.0D0) 2397 2398 call qenter('cc_r12prop') 2399 if (locdbg) then 2400 write(lupri,*) 'Entered CC_R12PROP' 2401 write(lupri,*) 'memory available: ',lwork,' double words' 2402 end if 2403 2404 ! only total symmetric Xi 2405 isymxi = 1 2406 2407 kxir12 = 1 2408 kzetar12 = kxir12 + ntr12am(isymxi) 2409 kzetar12sq = kzetar12 + ntr12am(1) 2410 kend1 = kzetar12sq + ntr12sq(1) 2411 lwrk1 = lwork - kend1 2412 if (locdbg) then 2413 kzeta1 = kend1 2414 kzeta2 = kzeta1 + nt1amx 2415 kend1 = kzeta2 + nt2amx 2416 end if 2417 if (lwrk1 .lt. 0) then 2418 call quit ('Insufficient work memory in CC_R12PROP') 2419 end if 2420 2421C----------------------------------------------------------------------- 2422C Read Xi vector from file if it exists... 2423C----------------------------------------------------------------------- 2424 ! check if this is correct ! 2425c idxxi = indprpcc(labelh) 2426c idxxi = IRHSR1(labelh,.FALSE.,0.0D0,isym) 2427c if (idxxi .ne. -1) then 2428c if (isym .ne. isymxi) then 2429c call quit('Symmetry error in CC_R12PROP') 2430c end if 2431c iopt = 32 2432c filxi = 'R1 ' 2433c call cc_rdrsp(filxi,idxxi,isymxi,iopt,model,dummy,work(kxir12)) 2434c else 2435C----------------------------------------------------------------------- 2436C ...otherwise calculate Xi vector 2437C----------------------------------------------------------------------- 2438C write(lupri,*) 'Will calculate R12 part for ',labelh 2439 call flshfo(lupri) 2440 ! allocate memory 2441 ktr12 = kend1 2442 ktr12sq = ktr12 + ntr12am(1) 2443 kxir12sq = ktr12sq + ntr12sq(1) 2444 kxint = kxir12sq + nr12r12sq(1) 2445 kxintsq = kxint + nr12r12p(1) 2446 kvxintsq = kxintsq + nr12r12sq(1) 2447 kprpao = kvxintsq + nr12r12sq(isymxi) 2448 kt1amp0 = kprpao + n2bst(isymxi) 2449 klamp0 = kt1amp0 + NT1AMX 2450 klamh0 = klamp0 + NLAMDT 2451 kend2 = klamh0 + NLAMDT 2452 LWRK2 = LWORK - KEND2 2453 IF (LWRK2 .LT. 0) THEN 2454 CALL QUIT('Insufficient work space in CC_R12PROP') 2455 END IF 2456 2457 ! read R12 amplitudes from disk 2458 iopt = 32 2459 call cc_rdrsp('R0 ',0,1,iopt,model,dummy,work(ktr12)) 2460 iopt = 1 2461 call ccr12unpck2(work(ktr12),1,work(ktr12sq),'N',iopt) 2462 2463 ! read X-integrals (R12 overlap matrix) 2464 luxint = -1 2465 call gpopen(luxint,fccr12x,'old',' ','unformatted',idummy, 2466 & .false.) 2467 rewind(luxint) 2468 9999 read(luxint) ian 2469 read(luxint) (work(kxint+i), i=0, nr12r12p(1)-1 ) 2470 if (ian.ne.ianr12) goto 9999 2471 call gpclose(luxint,'KEEP') 2472 iopt = 2 2473 call ccr12unpck2(work(kxint),1,work(kxintsq),'N',iopt) 2474 2475 ! read in VXINT 2476 call dzero(work(kvxintsq),nr12r12sq(isymxi)) 2477 call cc_r12rdvxint(work(kvxintsq),work(kend2),lwrk2,one, 2478 & isymxi,labelh) 2479 2480 ! read in V (perturbation operator) matrix in AO-basis 2481 call ccprpao(labelh,.TRUE.,work(kprpao),isymv,isym, 2482 & ierr,work(kend2),lwrk2) 2483 IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. isymv.NE.isymxi)) THEN 2484 CALL QUIT('CC_XIETA1: error while reading operator '//LABELH) 2485 ELSE IF (IERR.LT.0) THEN 2486 CALL DZERO(work(kprpao),N2BST(isymxi)) 2487 END IF 2488 2489 ! generate Lambda matrices 2490 iopt = 1 2491 call cc_rdrsp('R0 ',0,1,iopt,model,work(kt1amp0),dummy) 2492 call lammat(work(klamp0),work(klamh0),work(kt1amp0), 2493 & work(kend2),lwrk2) 2494 2495 ! calulate R12 part of Xi 2496 call cc_r12xi(work(kxir12sq),isymxi,'N',work(ktr12sq),1, 2497 & work(kxintsq),work(kvxintsq),isymxi,work(kprpao), 2498 & work(klamp0),work(klamh0),'N',work(kend2),lwrk2) 2499 2500 ! pack Xi to triangular format 2501 iopt = 1 2502 call ccr12pck2(work(kxir12),isymxi,.false.,work(kxir12sq), 2503 & 'N',iopt) 2504 call cclr_diasclr12(work(kxir12),brascl,isymxi) 2505 2506c end if 2507 2508 if (locdbg) then 2509 call around('Xi R12 part in CC_R12PROP') 2510 call cc_prpr12(work(kxir12),isymxi,1,.FALSE.) 2511 call FLSHFO(LUPRI) 2512 end if 2513 2514C----------------------------------------------------------------------- 2515C Read Lagrangian multipliers from file 2516C----------------------------------------------------------------------- 2517 listl = 'L0 ' 2518 idlstl = 0 2519 iopt = 32 2520 call cc_rdrsp(listl,idlstl,1,iopt,model,dummy,work(kzetar12)) 2521 2522 if (locdbg) then 2523 !Read conventional multipliers (for comparison) 2524 IOPT = 3 2525 CALL CC_RDRSP(listl,idlstl,1,IOPT,MODEL,WORK(kzeta1), 2526 & WORK(kzeta2)) 2527 2528 call around('Lagrangian multipliers') 2529 call cc_prp(work(kzeta1),work(kzeta2),1,1,1) 2530 call cc_prpr12(work(kzetar12),1,1,.TRUE.) 2531 end if 2532 2533C----------------------------------------------------------------------- 2534C Calculate dot product 2535C----------------------------------------------------------------------- 2536 propr12 = ddot(ntr12am(1),work(kzetar12),1,work(kxir12),1) 2537 2538 if (locdbg) then 2539 diacont = 0.0D0 2540 do isym = 1, nsym 2541 do i = 1, nmatki(isym) 2542 idx = itr12am(isym,isym) + i*(i+1)/2 - 1 2543 diacont = diacont + work(kzetar12+idx)*work(kxir12+idx) 2544 end do 2545 end do 2546 write(lupri,*) 'R12 contribution analysis for operator ', 2547 & labelh 2548 write(lupri,*) 'propr12 = ', propr12 2549 write(lupri,*) 'Diagonal contribution = ',diacont 2550 write(lupri,*) 'Off-diagonal cont. = ',propr12-diacont 2551 end if 2552 2553 if (locdbg) write(lupri,*) 'Leaving CC_R12PROP' 2554 call qexit('cc_r12prop') 2555 return 2556 end 2557*======================================================================= 2558 2559*======================================================================= 2560 subroutine cc_r12etaa(ETAA,ISYRES,CTR12,ISYCTR,TR12,ISYT12,XINT, 2561 & PRPAO,ISYMV,LAMDP,LAMDH,LAO,WORK,LWORK) 2562C----------------------------------------------------------------------- 2563C purpose: r12 contribution to Eta{A} singles part 2564C 2565C ETAA singles part of Eta{A} 2566C ISYRES symmetry of Eta{A} 2567C CTR12 R12 doubles Lagr. multipliers (ntr12sq(isyctr)) 2568C ISYCTR symmetry of Lagr. multipliers 2569C TR12 R12 doubles amplitudes / trial vector (ntr12sq(isyt12)) 2570C ISYT12 R12 amplitudes symmetry 2571C XINT R12 overlap matrix (nr12r12sq(1)) 2572C PRPAO Perturbation operator integrals in AO basis 2573C ISYMV Symmtetry of perturbation operator 2574C LAMDP CMO matrix used for occ. index (assumed symm. 1) 2575C LAMDH CMO matrix used for vir. index (assumed symm. 1) 2576C LAO flag: compute Eta{A}_(alpha i), i.e. "vitual" index in AO 2577C basis (skip second transformation of PRPAO) 2578C 2579C Christian Neiss April 2005 2580C----------------------------------------------------------------------- 2581 2582 implicit none 2583#include "priunit.h" 2584#include "ccorb.h" 2585#include "ccsdsym.h" 2586#include "ccsdinp.h" 2587#include "r12int.h" 2588#include "ccr12int.h" 2589#include "dummy.h" 2590 2591 logical locdbg, lao 2592 parameter (locdbg = .false.) 2593 2594 integer isyres,isyctr,isym1,isym2,isymkl,isymmj,isymmi,isymm, 2595 & isymi,isymj,isyma,isyt12,isymv 2596 integer lwork,kend1,kend2,lwrk1,lwrk2,kvint,kscr1,kscr2 2597 integer koffv,koffc,koffscr1,koffscr2,koffres,koffctr 2598 2599#if defined (SYS_CRAY) 2600 real work(*),etaa(*),ctr12(*),tr12(*),xint(*),prpao(*), 2601 & lamdp(*),lamdh(*),zero,one 2602#else 2603 double precision work(*),etaa(*),ctr12(*),tr12(*),xint(*), 2604 & prpao(*),lamdp(*),lamdh(*),zero,one 2605#endif 2606 2607 parameter(zero = 0.0D0, one = 1.0D0) 2608 2609 call qenter('cc_r12etaa') 2610 if (locdbg) then 2611 write(lupri,*) 'Entered CC_R12ETAA' 2612 write(lupri,*) 'memory available: ',lwork,' double words' 2613 end if 2614 2615 !check symmetry: 2616 if (isyres .ne. muld2h(isyctr,muld2h(isyt12,isymv))) then 2617 call quit('Symmetry error in CC_R12ETAA') 2618 end if 2619 2620 if (locdbg) then 2621 call around('Eta{O} singles before R12 contribution') 2622 call cc_prp(etaa,dummy,isyres,1,0) 2623 end if 2624 2625C-------------------- 2626C allocate memory 2627C-------------------- 2628 kvint = 1 2629 if (lao) then 2630 kend1 = kvint + nt1ao(isymv) 2631 else 2632 kend1 = kvint + nt1am(isymv) 2633 end if 2634 lwrk1 = lwork - kend1 2635 if (lwrk1 .lt. 0) then 2636 call quit('Insufficient work space in CC_R12ETAA') 2637 end if 2638 2639C----------------------------------------------------- 2640C transform perturbation integrals in MO-basis: V(a,j) 2641C----------------------------------------------------- 2642 kscr1 = kend1 2643 kend2 = kscr1 + nt1ao(isymv) 2644 lwrk2 = lwork - kend2 2645 if (lwrk2 .lt. 0) then 2646 call quit('Insufficient work space in CC_R12ETAA') 2647 end if 2648 2649 do isym2 = 1, nsym 2650 isym1 = muld2h(isymv,isym2) 2651 koffv = 1 + iaodis(isym1,isym2) 2652 koffscr1 = kscr1 + it1ao(isym1,isym2) 2653 koffc = 1 + iglmrh(isym2,isym2) 2654 !first transformation 2655 call dgemm('N','N',nbas(isym1),nrhf(isym2),nbas(isym2), 2656 & one,prpao(koffv),max(1,nbas(isym1)),lamdp(koffc), 2657 & max(1,nbas(isym2)),zero,work(koffscr1), 2658 & max(1,nbas(isym1))) 2659 !second transformation 2660 if (lao) then 2661 call dcopy(nt1ao(isymv),work(koffscr1),1,work(kvint+ 2662 & it1ao(isym1,isym2)),1) 2663 else 2664 koffv = kvint + it1am(isym1,isym2) 2665 koffc = 1 + iglmvi(isym1,isym1) 2666 !Note that the virtual index is leading in V(a,j)! 2667 call dgemm('T','N',nvir(isym1),nrhf(isym2),nbas(isym1), 2668 & one,lamdh(koffc),max(1,nbas(isym1)),work(koffscr1), 2669 & max(1,nbas(isym1)),zero,work(koffv), 2670 & max(1,nvir(isym1))) 2671 end if 2672 end do 2673 2674 if (locdbg ) then 2675 call around('Lambda^p matrix, all active orbitals:') 2676 do isym2 = 1, nsym 2677 isym1 = muld2h(isymv,isym2) 2678 write(lupri,*) 'Symmetry block: ',isym1, isym2 2679 call output(lamdp(1+iglmrh(isym1,isym2)),1,nbas(isym1),1, 2680 & norb(isym2),nbas(isym1),norb(isym2),1,lupri) 2681 end do 2682 call around('Lambda^p matrix, occ. part:') 2683 do isym2 = 1, nsym 2684 isym1 = muld2h(isymv,isym2) 2685 write(lupri,*) 'Symmetry block: ',isym1, isym2 2686 call output(lamdp(1+iglmrh(isym1,isym2)),1,nbas(isym1),1, 2687 & nrhf(isym2),nbas(isym1),nrhf(isym2),1,lupri) 2688 end do 2689 call around('Lambda^p matrix, virt. part:') 2690 do isym2 = 1, nsym 2691 isym1 = muld2h(isymv,isym2) 2692 write(lupri,*) 'Symmetry block: ',isym1, isym2 2693 call output(lamdp(1+iglmvi(isym1,isym2)),1,nbas(isym1),1, 2694 & nvir(isym2),nbas(isym1),nvir(isym2),1,lupri) 2695 end do 2696 call around('Original V integrals:') 2697 do isym2 = 1, nsym 2698 isym1 = muld2h(isymv,isym2) 2699 write(lupri,*) 'Symmetry block: ',isym1, isym2 2700 call output(prpao(1+iaodis(isym1,isym2)),1,nbas(isym1),1, 2701 & nbas(isym2),nbas(isym1),nbas(isym2),1,lupri) 2702 end do 2703 call around('Half transformed V integrals:') 2704 do isym2 = 1, nsym 2705 isym1 = muld2h(isymv,isym2) 2706 write(lupri,*) 'Symmetry block: ',isym1, isym2 2707 call output(work(kscr1+it1ao(isym1,isym2)),1,nbas(isym1), 2708 & 1,nrhf(isym2),nbas(isym1),nrhf(isym2),1,lupri) 2709 end do 2710 end if 2711 2712 if (locdbg.and.(.not.lao)) then 2713 call around('Perturbation operator integrals V(a,j)') 2714 call cc_prp(work(kvint),dummy,isymv,1,0) 2715 end if 2716 2717 if (locdbg) then 2718 call around('R12 ground state amplitudes') 2719 call cc_prsqr12(tr12,isyt12,'N',1,.false.) 2720 call around('R12 Lagrangian multipliers') 2721 call cc_prsqr12(ctr12,isyctr,'N',1,.false.) 2722 end if 2723 2724C---------------------- 2725C make contractions 2726C---------------------- 2727 kscr1 = kend1 2728 kscr2 = kscr1 + ntr12sq(isyt12) 2729 kend2 = kscr2 + nmatij(muld2h(isyctr,isyt12)) 2730 lwrk2 = lwork - kend2 2731 if (lwrk2 .lt. 0) then 2732 call quit('Insufficient work space in CC_R12ETAA') 2733 end if 2734 2735 ! first contraction over k'l': 2736 call dzero(work(kscr1),ntr12sq(isyt12)) 2737 call cc_r12xi2b(work(kscr1),'N',xint,1,'N',tr12,isyt12,'N',one) 2738 2739 if (locdbg) then 2740 call around('c*X') 2741 call cc_prsqr12(work(kscr1),isyt12,'N',1,.false.) 2742 end if 2743 2744 ! second contraction over klm: 2745 call dzero(work(kscr2),nmatij(muld2h(isyctr,isyt12))) 2746 do isymkl = 1, nsym 2747 isymmj = muld2h(isymkl,isyt12) 2748 isymmi = muld2h(isymkl,isyctr) 2749 do isymm = 1, nsym 2750 isymj = muld2h(isymmj,isymm) 2751 isymi = muld2h(isymmi,isymm) 2752 koffscr1 = kscr1 + itr12sq(isymkl,isymmj) + 2753 & nmatkl(isymkl)*imatij(isymm,isymj) 2754 koffctr = 1 + itr12sq(isymkl,isymmi) + 2755 & nmatkl(isymkl)*imatij(isymm,isymi) 2756 koffscr2 = kscr2 + imatij(isymj,isymi) 2757C write(lupri,*) 'symmetry isymkl, isymm: ', 2758C & isymkl,isymm 2759C write(lupri,*) 'SCR1:' 2760C call output(work(koffscr1),1,nmatkl(isymkl)*nrhf(isymm),1, 2761C & nrhf(isymj),nmatkl(isymkl)*nrhf(isymm), 2762C & nrhf(isymj),1,lupri) 2763C write(lupri,*) 'CTR:' 2764C call output(ctr12(koffctr),1,nmatkl(isymkl)*nrhf(isymm),1, 2765C & nrhf(isymi),nmatkl(isymkl)*nrhf(isymm), 2766C & nrhf(isymi),1,lupri) 2767 call dgemm('T','N',nrhf(isymj),nrhf(isymi), 2768 & nmatkl(isymkl)*nrhf(isymm),one,work(koffscr1), 2769 & max(1,nmatkl(isymkl)*nrhf(isymm)),ctr12(koffctr), 2770 & max(1,nmatkl(isymkl)*nrhf(isymm)),one, 2771 & work(koffscr2),max(1,nrhf(isymj))) 2772C write(lupri,*)'SCR2: imatij(isymj,isymi)=',imatij(isymj,isymi) 2773C call output(work(koffscr2),1,nrhf(isymj),1, 2774C & nrhf(isymi),nrhf(isymj),nrhf(isymi),1,lupri) 2775 end do 2776 end do 2777 2778 if (locdbg) then 2779 call around('ctr*c*X') 2780 do isymi = 1, nsym 2781 isymj = muld2h(isymi,muld2h(isyctr,isyt12)) 2782 koffscr2 = kscr2 + imatij(isymj,isymi) 2783 write(lupri,*) 'Symmetry block: ',isymj, isymi 2784 if (nrhf(isymj).eq.0 .or. nrhf(isymi).eq.0) then 2785 write(lupri,*) 'This symmetry is empty' 2786 else 2787 call output(work(koffscr2),1,nrhf(isymj),1,nrhf(isymi), 2788 & nrhf(isymj),nrhf(isymi),1,lupri) 2789 end if 2790 end do 2791 end if 2792 2793 ! last contraction over j: 2794 if (lao) then 2795 do isymj = 1, nsym 2796 isyma = muld2h(isymv,isymj) 2797 isymi = muld2h(isymj,muld2h(isyctr,isyt12)) 2798 koffscr2 = kscr2 + imatij(isymj,isymi) 2799 koffv = kvint + it1ao(isyma,isymj) 2800 koffres = 1 + it1ao(isyma,isymi) 2801c call output(work(koffv),1,nbas(isyma),1,nrhf(isymj), 2802c & nbas(isyma),nrhf(isymj),1,lupri) 2803c call output(work(koffscr2),1,nrhf(isymj),1,nrhf(isymi), 2804c & nrhf(isymj),nrhf(isymi),1,lupri) 2805 call dgemm('N','N',nbas(isyma),nrhf(isymi),nrhf(isymj),-one, 2806 & work(koffv),max(1,nbas(isyma)),work(koffscr2), 2807 & max(1,nrhf(isymj)),one,etaa(koffres), 2808 & max(1,nbas(isyma))) 2809 end do 2810 else 2811 do isymj = 1, nsym 2812 isyma = muld2h(isymv,isymj) 2813 isymi = muld2h(isymj,muld2h(isyctr,isyt12)) 2814 koffscr2 = kscr2 + imatij(isymj,isymi) 2815 koffv = kvint + it1am(isyma,isymj) 2816 koffres = 1 + it1am(isyma,isymi) 2817 call dgemm('N','N',nvir(isyma),nrhf(isymi),nrhf(isymj),-one, 2818 & work(koffv),max(1,nvir(isyma)),work(koffscr2), 2819 & max(1,nrhf(isymj)),one,etaa(koffres), 2820 & max(1,nvir(isyma))) 2821 end do 2822 end if 2823 2824 if (locdbg) then 2825 if (lao) then 2826 call around('Eta{O}_(alpha j)') 2827 do isym2 = 1, nsym 2828 isym1 = muld2h(isyres,isym2) 2829 write(lupri,*) 'Symmetry block: ',isym1, isym2 2830 call output(etaa(1+it1ao(isym1,isym2)),1,nbas(isym1),1, 2831 & nrhf(isym2),nbas(isym1),nrhf(isym2),1,lupri) 2832 end do 2833 else 2834 call around('Eta{O} singles after R12 contribution') 2835 call cc_prp(etaa,dummy,isyres,1,0) 2836 end if 2837 end if 2838 2839 if (locdbg) write(lupri,*) 'Leaving CC_R12ETAA' 2840 call qexit('cc_r12etaa') 2841 return 2842 end 2843*======================================================================= 2844 2845*======================================================================* 2846 subroutine cc_r12sort(rxkyl,isymr,rxkly,y,isymy,istartx,iendx, 2847 & isymx,dimx,dimk) 2848c---------------------------------------------------------------------- 2849c purpose: sort R_xy^kl as R_x,kl^y with fixed y 2850c expect R_xy^kl stored as a lower triangular matrix 2851c substitute routine cc_r12sortk, more flexible 2852c 2853c it is assumed that dimx=dimy, dimk=diml 2854c 2855c rxkyl R_xy^kl 2856c rxkly R_x,kl^y with fixed y 2857c isymr sym. of R_xy^kl 2858c y index y within isymy 2859c isymy sym. of y 2860c istartx startvalue for x in rxkly within isymx 2861c iendx endvalue for x in rxkly within isymx 2862c dimx dimension of x in rxkyl (=dimy) 2863c dimk dimension of k in rxkyl (=diml) 2864c 2865c C. Neiss, jan. 2006 2866c---------------------------------------------------------------------- 2867 implicit none 2868#include "priunit.h" 2869#include "ccorb.h" 2870 2871 integer isymr,isymy,istartx,iendx,isymx,dimx(8),dimk(8),x,y 2872 integer nxk(8),ixk(8,8),nkl(8),ikl(8,8),nxkyl(8),ixkyl(8,8) 2873 integer isym,isym1,isym2,isymkl,isymxkl,isymk,isyml, 2874 & isymxk,isymyl,xdimout 2875 integer idxkl,idxxk,idxyl,idxxkl,idxxkyl 2876 integer i,j,k,l,index 2877#if defined (SYS_CRAY) 2878 real rxkyl(*),rxkly(*) 2879#else 2880 double precision rxkyl(*),rxkly(*) 2881#endif 2882 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 2883 2884 call qenter('cc_r12sort') 2885 2886C calculate dimensions and offsets: 2887 do isym = 1, nsym 2888 nxk(isym) = 0 2889 nkl(isym) = 0 2890 do isym2 = 1, nsym 2891 isym1 = muld2h(isym,isym2) 2892 ixk(isym1,isym2) = nxk(isym) 2893 ikl(isym1,isym2) = nkl(isym) 2894 nxk(isym) = nxk(isym) + dimx(isym1)*dimk(isym2) 2895 nkl(isym) = nkl(isym) + dimk(isym1)*dimk(isym2) 2896 end do 2897 end do 2898C 2899 do isym = 1, nsym 2900 nxkyl(isym) = 0 2901 do isym2 = 1, nsym 2902 isym1 = muld2h(isym,isym2) 2903 if (isym2.gt.isym1) then 2904 ixkyl(isym1,isym2) = nxkyl(isym) 2905 ixkyl(isym2,isym1) = nxkyl(isym) 2906 nxkyl(isym) = nxkyl(isym) + nxk(isym1)*nxk(isym2) 2907 else if (isym2.eq.isym1) then 2908 ixkyl(isym1,isym2) = nxkyl(isym) 2909 nxkyl(isym) = nxkyl(isym) + nxk(isym1)*(nxk(isym1)+1)/2 2910 end if 2911 end do 2912 end do 2913C 2914 xdimout = iendx - istartx + 1 2915C 2916 isymxkl = muld2h(isymr,isymy) 2917 isymkl = muld2h(isymxkl,isymx) 2918 call dzero(rxkly,xdimout*nkl(isymkl)) 2919C 2920C start resort: 2921 do isymk = 1, nsym 2922 isyml = muld2h(isymkl,isymk) 2923 isymxk = muld2h(isymx,isymk) 2924 isymyl = muld2h(isymy,isyml) 2925 do l = 1, dimk(isyml) 2926 idxyl = ixk(isymy,isyml) + dimx(isymy)*(l-1) + y 2927 do k = 1, dimk(isymk) 2928 idxkl = ikl(isymk,isyml) + dimk(isymk)*(l-1) + k 2929 if (isymxk.eq.isymyl) then 2930 do x = istartx, iendx 2931 idxxkl = xdimout*(idxkl-1) + x - istartx + 1 2932 idxxk = ixk(isymx,isymk) + dimx(isymx)*(k-1) + x 2933 idxxkyl = ixkyl(isymxk,isymyl) + index(idxxk,idxyl) 2934 rxkly(idxxkl) = rxkyl(idxxkyl) 2935 end do 2936 else if (isymxk.lt.isymyl) then 2937 idxxkl = xdimout*(idxkl-1) + 1 2938 idxxk = ixk(isymx,isymk) + dimx(isymx)*(k-1) + istartx 2939 idxxkyl = ixkyl(isymxk,isymyl) + nxk(isymxk)*(idxyl-1) + 2940 & idxxk 2941 call dcopy(xdimout,rxkyl(idxxkyl),1,rxkly(idxxkl),1) 2942 else if (isymxk.gt.isymyl) then 2943 do x = istartx, iendx 2944 idxxkl = xdimout*(idxkl-1) + x - istartx + 1 2945 idxxk = ixk(isymx,isymk) + dimx(isymx)*(k-1) + x 2946 idxxkyl = ixkyl(isymyl,isymxk) + nxk(isymyl)*(idxxk-1) + 2947 & idxyl 2948 rxkly(idxxkl) = rxkyl(idxxkyl) 2949 end do 2950 end if 2951 end do 2952 end do 2953 end do 2954 2955 call qexit('cc_r12sort') 2956 end 2957*======================================================================* 2958 2959*======================================================================* 2960 subroutine cc_r12cmo(cmo,work,lwork) 2961c---------------------------------------------------------------------- 2962c purpose: read CMO-Matrix for LABEL=FULLBAS and delete redundant 2963c orbitals 2964c 2965c C. Neiss, march 2006 2966c---------------------------------------------------------------------- 2967 implicit none 2968#include "priunit.h" 2969#include "ccsdinp.h" 2970#include "ccorb.h" 2971#include "ccsdsym.h" 2972#include "r12int.h" 2973 2974 integer lusifc,idum,isym,nsymx,norbtsx,nbastx,nlamdsx,nlamdsxtd, 2975 & nrhfsx(8),norbsx(8),nbasxtd(8),norbxtd(8) 2976 integer kcmox,kend1,lwrk1,koffc,koff1,koff2,lwork 2977 logical locdbg 2978 parameter (locdbg=.FALSE.) 2979#if defined (SYS_CRAY) 2980 real cmo(*),work(*) 2981#else 2982 double precision cmo(*),work(*) 2983#endif 2984 2985 call qenter('cc_r12cmo') 2986C----------------------------------------------------------------------- 2987C define dimensions/offsets 2988C----------------------------------------------------------------------- 2989 if (.not. CCR12) THEN 2990 do isym = 1, nsym 2991 mbas1(isym) = nbas(isym) 2992 mbas2(isym) = 0 2993 norb1(isym) = norb(isym) 2994 norb2(isym) = 0 2995 end do 2996 end if 2997 2998 nlamdsxtd = 0 2999 do isym = 1, nsym 3000 nbasxtd(isym) = mbas1(isym) + mbas2(isym) 3001 norbxtd(isym) = norb1(isym) + norb2(isym) 3002 nlamdsxtd = nlamdsxtd + nbasxtd(isym)*norbxtd(isym) 3003 end do 3004 3005C----------------------------------------------------------------------- 3006C read in MO coefficients 3007C----------------------------------------------------------------------- 3008 lusifc = -1 3009 call gpopen(lusifc,'SIRIFC','OLD',' ','UNFORMATTED',idum,.false.) 3010 ! read dimensions for CMO coefficients for full basis (nlamdsx) 3011 rewind(lusifc) 3012 call mollab('FULLBAS ',lusifc,lupri) 3013 read(lusifc) nsymx,norbtsx,nbastx,nlamdsx,(nrhfsx(i),i=1,nsym), 3014 & (norbsx(i),i=1,nsym) 3015C allocate memory for MO coefficients: 3016 kcmox = 1 3017 kend1 = kcmox + nlamdsx 3018 lwrk1 = lwork - kend1 3019 if (lwrk1 .lt. 0) then 3020 call quit ('Insufficient work memory in CC_R12CMO') 3021 end if 3022 call dzero(cmo,nlamdsxtd) 3023 call dzero(work(kcmox),nlamdsx) 3024 read(lusifc) 3025 read(lusifc) (work(kcmox+i-1),i=1,nlamdsx) 3026 call gpclose(lusifc,'KEEP') 3027 if (locdbg) then 3028 write(lupri,*) 'nsymx, norbtsx, nbastx, nlamdsx: ', 3029 & nsymx, norbtsx, nbastx, nlamdsx 3030 do isym = 1, nsym 3031 write(lupri,*) 'nrhfsx(',isym,') = ', nrhfsx(isym) 3032 write(lupri,*) 'norbsx(',isym,') = ', norbsx(isym) 3033 end do 3034 3035 call around('MO-coefficient matrix incl. redundant orbitals') 3036 do isym = 1, nsym 3037 koffc = kcmox 3038 write(lupri,*) 'Symmetry number:', isym 3039 if (norbsx(isym) .eq. 0) then 3040 write(lupri,*) 'This symmetry is empty' 3041 else 3042 call output(work(koffc),1,nbasxtd(isym),1,norbsx(isym), 3043 & nbasxtd(isym),norbsx(isym),1,lupri) 3044 end if 3045 koffc = koffc + nbasxtd(isym)*norbsx(isym) 3046 end do 3047 end if 3048 3049C delete redundant occupied orbitals from CMO matrix: 3050 koff1 = kcmox 3051 koff2 = 1 3052 do isym= 1, nsym 3053 koff1 = koff1 + nbasxtd(isym)*nrhfsx(isym) 3054 call dcopy(nbasxtd(isym)*norbxtd(isym),work(koff1),1, 3055 & cmo(koff2),1) 3056 koff1 = koff1 + nbasxtd(isym)*norbxtd(isym) 3057 koff2 = koff2 + nbasxtd(isym)*norbxtd(isym) 3058 end do 3059 if (locdbg) then 3060 call around('MO-coefficient matrix') 3061 do isym = 1, nsym 3062 koffc = 1 3063 write(lupri,*) 'Symmetry number:', isym 3064 if (norbxtd(isym) .eq. 0) then 3065 write(lupri,*) 'This symmetry is empty' 3066 else 3067 call output(cmo(koffc),1,nbasxtd(isym),1,norbxtd(isym), 3068 & nbasxtd(isym),norbxtd(isym),1,lupri) 3069 end if 3070 koffc = koffc + nbasxtd(isym)*norbxtd(isym) 3071 end do 3072 end if 3073 call qexit('cc_r12cmo') 3074 end 3075*======================================================================* 3076 3077