1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18*=====================================================================* 19 subroutine cc_r12whtf(xint,idel,isymd,j,isymj,isymbi,lauxbeta, 20 & lunitr12,filer12,lunitr12_2,filer12_2, 21 & work,lwork) 22c---------------------------------------------------------------------- 23c purpose: get half transformed r12 integrals and write them 24c on file first delta then delta' in ccsd_iabj2 25c 26c xint array of half transformed 12 integrals 27c (i beta|r12|j delta) with j, delta fixed 28c idel number of delta (running over all symmetries; in 29c contrast, idelta starts with 1 for every symmetry!) 30c isymd symmetry of delta 31c j number of j (starts with 1 for every symmetry) 32c isymj symmetry of j 33c isymbi symmetry of (i beta) 34c lauxbeta flag to write out also integrals with beta being 35c auxiliary basis function (needed for R12 response) 36c lunitr12, filer12 unit number, file name to be written to 37c lunitr12_2, filer12_2 unit number, file name for integrals with 38c two auxiliary basis functions 39c 40c H. Fliegl, C. Haettig, winter 2003 41c modified C. Neiss, fall 2004: write also integrals with beta 42c auxiliary basis functions if desired 43c---------------------------------------------------------------------- 44 implicit none 45#include "priunit.h" 46#include "ccorb.h" 47#include "ccsdsym.h" 48#include "r12int.h" 49#include "dummy.h" 50 51 logical locdbg,lauxbeta 52 parameter (locdbg = .false.) 53 integer lunitr12, lunitr12_2, lwork, idel, isymd, isymi, isymb 54 integer isymbi, isymdj, isymj, ndelj, iadr, ilen 55 integer isym, isym1, isym2, icount2, icount3, icount4, icount8 56 integer idxbi, idxdj, idxbi2, koff1 57 integer nr1bas(8),nr1xbas(8),nr2bas,nr2bas2 58 integer ir1bas(8,8), ir1xbas(8,8), ir2bas(8,8) 59 integer ir2xbas(8,8),idelta 60 integer ir2bas2(8,8),ir2xbas2(8,8) 61 integer nrhf1(8) 62 63 character*(*) filer12, filer12_2 64 65#if defined (SYS_CRAY) 66 real work(*), xint(*) 67#else 68 double precision work(*), xint(*) 69#endif 70 71 call qenter('r12whtf') 72 73 if (locdbg) then 74 write(lupri,*) 'Entered CC_R12WHTF' 75 write(lupri,*) 'lauxbeta,r12int = ',lauxbeta,r12int 76 write(lupri,'(A,I5,2X,A)') 77 & ' LUNITR12, FILER12 :',lunitr12,filer12, 78 & ' LUNITR12_2, FILER12_2:',lunitr12_2,filer12_2 79 write(lupri,*) 'V12INT = ',V12INT 80 end if 81 82 if (V12INT) then 83 ! use active occupied orbitals 84 call icopy(8,nrhfa,1,nrhf1,1) 85 else 86 ! use r12 orbitals 87 call icopy(8,nrhfb,1,nrhf1,1) 88 end if 89 90 if (j.gt.nrhf1(isymj)) then 91 call qexit('r12whtf') 92 return 93 end if 94 95 do isym = 1, nsym 96 nr1bas(isym) = 0 97 nr1xbas(isym) = 0 98 do isym2 = 1, nsym 99 isym1 = muld2h(isym,isym2) 100 nr1bas(isym) = nr1bas(isym) + mbas1(isym1)*nrhf1(isym2) 101 nr1xbas(isym)= nr1xbas(isym)+ mbas2(isym1)*nrhf1(isym2) 102 end do 103 end do 104 nr2bas = 0 105 do isym = 1, nsym 106 nr2bas = nr2bas + nr1bas(isym)*nr1bas(isym) 107 end do 108 do isym = 1, nsym 109 icount2 = 0 110 icount3 = 0 111 icount4 = 0 112 icount8 = 0 113 do isym2 = 1, nsym 114 isym1 = muld2h(isym,isym2) 115 ir1bas(isym1,isym2) = icount2 116 ir1xbas(isym1,isym2)= icount8 117 ir2bas(isym1,isym2) = icount3 118 ir2xbas(isym1,isym2)= icount4 119 icount2 = icount2 + nrhf1(isym2)*mbas1(isym1) 120 icount8 = icount8 + nrhf1(isym2)*mbas2(isym1) 121 icount3 = icount3 + nr1bas(isym1)*nr1bas(isym2) 122 icount4 = icount4 + nr1bas(isym1)*nr1xbas(isym2) 123 end do 124 end do 125 126 if (lwork.lt.nr1bas(isymbi)) then 127 write(lupri,*) 'lwork, nr1bas(isymbi):',lwork, nr1bas(isymbi) 128 call quit('Insufficient work space in cc_r12whtf') 129 end if 130 131 idelta = idel - ibas(isymd) 132 133 koff1 = 1 134 isymdj = muld2h(isymd,isymj) 135 do isymi = 1, nsym 136 isymb = muld2h(isymbi,isymi) 137 do i = 1, nrhf1(isymi) 138 do b = 1, mbas1(isymb) 139 idxbi = ir1bas(isymb,isymi) + mbas1(isymb)*(i-1)+b 140 idxbi2 = koff1 -1 + nbas(isymb)*(i-1) +b 141 work(idxbi) = xint(idxbi2) 142 end do 143 end do 144 if (locdbg) then 145 write(lupri,*)'idel, j, isymdj =',idel,j,isymdj 146 write(lupri,*) 'symmetry block (beta,i):',isymb,isymi 147 call around('integrals before reordering') 148 call output(xint(koff1),1,nbas(isymb),1,nrhf(isymi), 149 & nbas(isymb),nrhf(isymi),1,lupri) 150 call around('integrals after reordering') 151 call output(work(ir1bas(isymb,isymi)+1),1,mbas1(isymb), 152 & 1,nrhf1(isymi),mbas1(isymb),nrhf1(isymi),1,lupri) 153 write(lupri,*) 154 end if 155 koff1 = koff1 + nbas(isymb)*nrhf(isymi) 156 end do 157 158 159 if (idelta.le.mbas1(isymd)) then 160 ndelj = ir1bas(isymd,isymj) + mbas1(isymd)*(j-1) + idelta 161 iadr = ir2bas(isymbi,isymdj) + nr1bas(isymbi)*(ndelj-1)+1 162 else 163 ndelj = ir1xbas(isymd,isymj) + 164 & mbas2(isymd)*(j-1) + idelta - mbas1(isymd) 165 iadr = nr2bas + ir2xbas(isymbi,isymdj) + 166 & nr1bas(isymbi)*(ndelj-1)+1 167 end if 168 ilen = nr1bas(isymbi) 169 call putwa2(lunitr12,filer12,work,iadr,ilen) 170 171C write also integrals with beta being auxiliary function on extra 172C file if desired: 173 if (lauxbeta) then 174 if (mbsmax.lt.6) then 175 write (lupri,*) 'MBSMAX = ',MBSMAX 176 call quit("MBSMAX must at least be 6 for "// 177 & "(k beta'|r12|l delta')") 178 end if 179 180 ! nr2bas2 = lenghts for nr1bas x nr1xbas over all symmetries 181 nr2bas2 = 0 182 do isym = 1, nsym 183 nr2bas2 = nr2bas2 + nr1bas(isym)*nr1xbas(isym) 184 end do 185 ! ir2bas2 = symmetry offsets for nr1xbas x nr1bas 186 ! ir2xbas2 = symmetry offsets for nr1xbas x nr1xbas 187 do isym = 1, nsym 188 icount3 = 0 189 icount4 = 0 190 do isym2 = 1, nsym 191 isym1 = muld2h(isym,isym2) 192 ir2bas2(isym1,isym2) = icount3 193 ir2xbas2(isym1,isym2) = icount4 194 icount3 = icount3 + nr1xbas(isym1)*nr1bas(isym2) 195 icount4 = icount4 + nr1xbas(isym1)*nr1xbas(isym2) 196 end do 197 end do 198 199 ! reallocate work space: 200 if (lwork.lt.nr1xbas(isymbi)) then 201 call quit('Insufficient work space in cc_r12whtf') 202 end if 203 204 koff1 = 1 205 isymdj = muld2h(isymd,isymj) 206 do isymi = 1, nsym 207 isymb = muld2h(isymbi,isymi) 208 do i = 1, nrhf1(isymi) 209 do b = mbas1(isymb) + 1, mbas1(isymb)+mbas2(isymb) 210 idxbi = ir1xbas(isymb,isymi) + mbas2(isymb)*(i-1) + b - 211 & mbas1(isymb) 212 idxbi2 = koff1 -1 + nbas(isymb)*(i-1) + b 213 work(idxbi) = xint(idxbi2) 214 end do 215 end do 216 if (locdbg) then 217 call around('integrals after reordering') 218 call output(work(ir1xbas(isymb,isymi)+1),1,mbas2(isymb), 219 & 1,nrhf1(isymi),mbas2(isymb),nrhf1(isymi),1,lupri) 220 write(lupri,*) 221 end if 222 koff1 = koff1 + nbas(isymb)*nrhf(isymi) 223 end do 224 225 if (idelta.le.mbas1(isymd)) then 226 ndelj = ir1bas(isymd,isymj) + mbas1(isymd)*(j-1) + idelta 227 iadr = ir2bas2(isymbi,isymdj) + nr1xbas(isymbi)*(ndelj-1)+1 228 else 229 ndelj = ir1xbas(isymd,isymj) + 230 & mbas2(isymd)*(j-1) + idelta - mbas1(isymd) 231 iadr = nr2bas2 + ir2xbas2(isymbi,isymdj) + 232 & nr1xbas(isymbi)*(ndelj -1) + 1 233 end if 234 ilen = nr1xbas(isymbi) 235 call putwa2(lunitr12_2,filer12_2,work,iadr,ilen) 236 end if 237 238 if (locdbg) then 239 write(lupri,*) 'Leaving CC_R12WHTF' 240 end if 241 242 call qexit('r12whtf') 243 end 244*=====================================================================* 245 subroutine cc_r12intf2(vitjtkl,lamdh,isymh,lamdhs,isymv, 246 & lamdhcs,isymhcs,work,lwork) 247c---------------------------------------------------------------------- 248c purpose: calculate V^itjt_kl for ansatz 2 249c 250c H. Fliegl, C. Haettig winter 2003 251c 252c \sum_MN r^MtNt_kl * g^alpha jt_MN 253c---------------------------------------------------------------------- 254 implicit none 255#include "priunit.h" 256#include "ccorb.h" 257#include "ccsdsym.h" 258#include "r12int.h" 259#include "ccr12int.h" 260#include "dummy.h" 261 262 logical ldum,lauxd,locdbg 263 parameter (locdbg = .false.) 264 265 integer isymv,krmtntkl,kgmtntkl,nt2aox,isymhs,isymhcs 266 integer krhtf, kend1, lwrk1, lwork, idum, idelta, isymd 267 integer isymk, isymdk, isymm, isymij, isymkl, isyml, koffv, 268 & kend2, koffr, koffg, isymnm, ntotmn, ntotij 269 integer nr1orb(8), nr1bas(8), nr1xbas(8), nr2bas, n2bst1(8) 270 integer ir1xbas(8,8),ibasx(8),isymdm,idxkl,idxlk,iadr1,iadr2 271 integer ir1orb(8,8), ir1bas(8,8), ir2bas(8,8), iaodis1(8,8) 272 integer ir2xbas(8,8), irgkl(8,8), irxgkl(8,8) 273 integer nrgkl(8), nrxgkl(8), lwrk2 274 integer idxd,idxdp,isym,icount1,isym1,isym2,kgtf,isyma,koffc, 275 & koff1,koff2,ntota,ntotm,isymn,idxmn,idxdn 276 integer iglmrhs(8,8),iglmvis(8,8),nglmds(8),igamsm2(8,8), 277 & icmo(8,8),ncmo(8),imaklm(8,8),nmaklm(8),ngamsm2(8),kgmkl, 278 & imatijm(8,8),nmatijm(8),icount5,ngamsm(8),igamsm(8,8), 279 & irgijs(8,8),nrgijs(8),ir1basm(8,8),nr1basm(8), 280 & ir2basm(8,8),nr2basm,ir1xbasm(8,8),nr1xbasm(8),nvajkls(8), 281 & ir2xbasm(8,8),imatf(8,8),nmatf(8),icount2,ivajkls(8,8) 282 integer nr1xorb(8),ir1xorb(8,8), 283 & nalphaj(8),ialphaj(8,8),nmaijm(8),imaijm(8,8) 284 integer isymaj,ntotaj,isymh,lunit,krsym,isymg,isymr 285 286 double precision one, zero,two 287 parameter (one = 1.0d0, zero = 0.0d0, two = 2.0d0) 288#if defined (SYS_CRAY) 289 real lamdhs(*), vitjtkl(*), work(*),ddot,lamdhcs(*), 290 & factor, lamdh(*),x1,x2,x3 291#else 292 double precision lamdhs(*), vitjtkl(*), work(*), ddot,lamdhcs(*), 293 & factor, lamdh(*),x1,x2,x3 294#endif 295 296 call qenter('r12intf2') 297 298 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 299 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas, 300 & ir2bas,ir2xbas,irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 301 call cc_r12offs23(iglmrhs,iglmvis,nglmds,icmo,ncmo, 302 & imaijm,nmaijm,imaklm,nmaklm, 303 & imatijm,nmatijm,igamsm,ngamsm,irgijs,nrgijs, 304 & ir1basm,nr1basm,ir2basm,nr2basm,ir1xbasm, 305 & nr1xbasm,ir2xbasm,imatf,nmatf) 306 307 if (locdbg) then 308 write(lupri,*)'V before tf2' 309 do isymij = 1, nsym 310 isymkl = muld2h(isymv,isymij) 311 write(lupri,*)'symmetry block (ij,kl)',isymij,isymkl 312 call output(vitjtkl(1+itr12sqt(isymij,isymkl)), 313 & 1,nmatij(isymij),1,nmatkl(isymkl), 314 & nmatij(isymij),nmatkl(isymkl),1,lupri) 315 end do 316 write(lupri,*)'norm^2 V before tf2', 317 & ddot(ntr12sq(isymv),vitjtkl,1,vitjtkl,1) 318 end if 319 320 ibasx(1) = 0 321 do isym = 2, nsym 322 ibasx(isym) = ibasx(isym-1) + mbas2(isym-1) 323 end do 324 325 do isym = 1, nsym 326 icount1 = 0 327 icount2 = 0 328 do isym2 = 1, nsym 329 isym1 = muld2h(isym,isym2) 330 igamsm2(isym1,isym2) = icount1 331 ivajkls(isym1,isym2) = icount2 332 icount1 = icount1 + nmatij(isym1)*nmatijm(isym2) 333 icount2 = icount2 + nalphaj(isym1)*nmatijm(isym2) 334 end do 335 ngamsm2(isym) = icount1 336 nvajkls(isym) = icount2 337 end do 338 339 isymhs = isymh 340 isymg = muld2h(isymh,isymh) 341 isymr = muld2h(isymhcs,isymhs) 342 343 krmtntkl = 1 344 kgmtntkl = krmtntkl + ngamsm(isymr) 345 kend1 = kgmtntkl + ngamsm2(isymg) 346 lwrk1 = lwork - kend1 347 if (lwrk1.lt.0) then 348 write(lupri,*)'lwork,lwrk1,kend1',lwork,lwrk1,kend1 349 call quit('Insufficient work space in cc_r12intf2') 350 end if 351 352 call dzero(work(krmtntkl),ngamsm(isymr)) 353 call dzero(work(kgmtntkl),ngamsm2(isymg)) 354 355 do isymd = 1, nsym 356 krhtf = kend1 357 kend2 = krhtf + max(nrgkl(isymd),nrgijs(isymd)) 358 lwrk2 = lwork - kend2 359 if (lwrk2.lt.0) then 360 write(lupri,*)'lwrk2, lwork, kend2', lwrk2, lwork, kend2 361 call quit('Insufficient work space in cc_r12intf2') 362 end if 363 364 do idelta = 1, mbas1(isymd) 365c read (ka|r12|lb) from file and sort as R^d_(ak,l) 366 lauxd = .false. 367 idxd = idelta+ibas(isymd) 368 call cc_r12getrint(work(krhtf),idxd,isymd,nr1bas,ir1bas, 369 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 370 & nrhfb,nmatkl,imatkl, 371 & ibasx,lauxd,.false.,frhtf,work(kend2),lwrk2) 372 if (locdbg) then 373 write(lupri,*)'Rhtf in tf2' 374 write(lupri,*)' norm^2 Rhtf in tf2', 375 & ddot(nrgkl(isymd),work(krhtf),1,work(krhtf),1) 376 end if 377 378c case response calculation: lres = .true. 379c get rmtntkl 380 call cc_r12generaltf(work(krhtf),work(krmtntkl),idelta, 381 & isymd,lamdhcs,isymhcs,lamdhs,isymhs, 382 & .false.,iglmrhs,iglmrhs,nmatkl,nrhfsa, 383 & nbas,irgkl,imatijm,nmatijm,imaklm, 384 & igamsm,work(kend2),lwrk2) 385 386c read (Ka|Lb) from file and sort as I^d_(aK,L) 387 lauxd = .false. 388 idxd = idelta+ibas(isymd) 389 call cc_r12getrints(work(krhtf),idxd,isymd,nr1basm,ir1basm, 390 & nr2basm,ir2basm,nrgijs,irgijs,ir1xbasm,ir2xbasm, 391 & ibasx,imatijm,nmatijm,lauxd,fccgmnab,work(kend2),lwrk2) 392 393c get gmtntKL 394 call cc_r12generaltf(work(krhtf),work(kgmtntkl),idelta, 395 & isymd,lamdh,isymh,lamdh,isymh, 396 & .false.,iglmrh,iglmrh,nmatijm,nrhf, 397 & nbas,irgijs,imatij,nmatij,imatf, 398 & igamsm2,work(kend2),lwrk2) 399 400 end do 401 end do 402c 403 if (locdbg) then 404c print result 405 write(lupri,*)'Rmtntkl', 406 & ddot(ngamsm(isymr),work(krmtntkl),1,work(krmtntkl),1) 407 write(lupri,*)'Imtntkl', 408 & ddot(ngamsm2(isymg),work(kgmtntkl),1,work(kgmtntkl),1) 409 end if 410 411 x1 = 0.0d0 412 x2 = 0.0d0 413 x3 = 0.0d0 414c 415c make Vitjtkl 416c 417 factor = one 418 do isymnm = 1, nsym 419 isymij = isymnm 420 isymkl = muld2h(isymv,isymij) 421 422 koffv = 1 + itr12sqt(isymij,isymkl) 423 koffr = krmtntkl + igamsm(isymnm,isymkl) 424 koffg = kgmtntkl + igamsm2(isymij,isymnm) 425 426 ntotmn = max(nmatijm(isymnm),1) 427 ntotij = max(nmatij(isymij),1) 428 429 call dgemm('N','N',nmatij(isymij),nmatkl(isymkl), 430 & nmatijm(isymnm),factor,work(koffg),ntotij, 431 & work(koffr),ntotmn,one,vitjtkl(koffv),ntotij) 432 433 if (locdbg) then 434 x1 = x1 + 435 & ddot(nmatijm(isymnm)*nmatij(isymij),work(koffg),1, 436 & work(koffg),1) 437 x2 = x2 + 438 & ddot(nmatijm(isymnm)*nmatkl(isymkl),work(koffr),1, 439 & work(koffr),1) 440 x3 = x3 + 441 & ddot(nmatij(isymij)*nmatkl(isymkl),vitjtkl(koffv),1, 442 & vitjtkl(koffv),1) 443 end if 444 445 end do 446c 447 if (locdbg) then 448 write(lupri,*)'V after tf2' 449 do isymij = 1, nsym 450 isymkl = muld2h(isymv,isymij) 451 write(lupri,*)'symmetry block (ij,kl)',isymij,isymkl 452 call output(vitjtkl(1+itr12sqt(isymij,isymkl)), 453 & 1,nmatij(isymij),1,nmatkl(isymkl), 454 & nmatij(isymij),nmatkl(isymkl),1,lupri) 455 end do 456 write(lupri,*)'norm^2 V after tf2', 457 & ddot(ntr12sq(isymv),vitjtkl,1,vitjtkl,1) 458 write(lupri,*)'x1 = ', x1 459 write(lupri,*)'x2 = ', x2 460 write(lupri,*)'x3 = ', x3 461 end if 462c 463 call qexit('r12intf2') 464 end 465*======================================================================* 466 subroutine ccr12prep2(work,lwork) 467*----------------------------------------------------------------------* 468c Purpose: make some intermediates used for ansatz 2 & 3 469c 470c H. Fliegl, C. Haettig, winter 2003 471c adapted for ansatz 3, Christof Haettig, spring 2005 472c adapted for CABS, Christian Neiss, winter 2006 473*----------------------------------------------------------------------* 474 implicit none 475#include "priunit.h" 476#include "ccorb.h" 477#include "ccsdsym.h" 478#include "r12int.h" 479#include "ccr12int.h" 480#include "dummy.h" 481 482#if defined (SYS_CRAY) 483 real zero,one,half 484#else 485 double precision zero,one,half 486#endif 487 parameter (zero = 0.0D0, one = 1.0D0, half = 0.5d0) 488 489 logical noauxbsave, locdbg, lauxd, do_virt 490 parameter (locdbg = .false.) 491 492 integer nr1orb(8),nr1bas(8),nr1xbas(8),nr2bas,n2bst1(8) 493 integer ir1xbas(8,8),basoff(8) 494 integer ir1orb(8,8),ir1bas(8,8),ir2bas(8,8),iaodis1(8,8) 495 integer ir2xbas(8,8),irgkl(8,8),nrgkl(8),ibasx(8) 496 integer nbasx(8),nnbastx,isym,isymdk,n2bst2(8),isymd,isymk, 497 & icount,kend0,lusifc,norbsx(8),nbastx,nlamdsx, 498 & nrhfsx(8),kcmox,lwrk0,lwork,ksfull 499 integer norbtsx,isympq,isymp,isymq,ksjbp,ksbbp,ksmupbp 500 integer isymbp,krhtf,ib,isymdp,idelp,idxbpdp,krdpakl 501 integer isymakl,kend1,lwrk1,koff1,klamdp,klamdh,kend2,lwrk2, 502 & isymj,koffs,kcmo,ntotc,ntots,ntots2,kgdakl, 503 & kghtf,idxddp,koff2,idel,kgdpakl,idxd,idxdp 504 integer kt1am,isymmu,kscr1,kstep,luvajkl,icoun1 505 integer kiajb, idxb,krtot, icount7, kvajkl, kend3,lwrk3 506 integer iglmrhs(8,8),iglmvis(8,8),nglmds(8),icmo(8,8),ncmo(8), 507 & imaklm(8,8),nmaklm(8),imatijm(8,8),nmatijm(8),igamsm(8,8), 508 & ngamsm(8),klamdhs,klamdps,irgijs(8,8),nrgijs(8), 509 & ir1basm(8,8),nr1basm(8),ir2basm(8,8),nr2basm, 510 & ir1xbasm(8,8),nr1xbasm(8),ir2xbasm(8,8),imatf(8,8), 511 & nmatf(8),isbbp(8,8),nsbbp(8),isjbp(8,8),nsjbp(8),nbas2(8) 512 integer nr1xorb(8),ir1xorb(8,8),nrxgkl(8),irxgkl(8,8), 513 & nalphaj(8),ialphaj(8,8),nmaijm(8),imaijm(8,8), 514 & nrgij(8),irgij(8,8),nt1xao(8),it1xao(8,8),it2xaos(8,8) 515#if defined (SYS_CRAY) 516 real work(*),ddot,xs 517#else 518 double precision work(*),ddot,xs 519#endif 520 521c---------------------------------------------------------------------- 522c some initializations: 523c---------------------------------------------------------------------- 524 call qenter('r12prep2') 525 if(locdbg) write(lupri,*) 'entered ccr12prep2' 526 527 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 528 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas,ir2bas, 529 & ir2xbas,irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 530 531 call cc_r12offs23(iglmrhs,iglmvis,nglmds,icmo,ncmo, 532 & imaijm,nmaijm,imaklm,nmaklm, 533 & imatijm,nmatijm,igamsm,ngamsm,irgijs,nrgijs, 534 & ir1basm,nr1basm,ir2basm,nr2basm,ir1xbasm, 535 & nr1xbasm,ir2xbasm,imatf,nmatf) 536c 537 if (r12cbs) then 538 do isym = 1, nsym 539 nbas2(isym) = mbas1(isym) + mbas2(isym) 540 basoff(isym) = 0 541 end do 542 else 543 call icopy(8,mbas2,1,nbas2,1) 544 call icopy(8,mbas1,1,basoff,1) 545 end if 546c 547c calculate number of auxbasis functions in previous symmetries 548 ibasx(1) = 0 549 do isym = 2,nsym 550 ibasx(isym) = ibasx(isym-1)+mbas2(isym-1) 551 end do 552c 553 nnbastx = 0 554 do isym = 1, nsym 555 nbasx(isym) = mbas1(isym) + mbas2(isym) 556 nnbastx = nnbastx + nbasx(isym)*(nbasx(isym)+1)/2 557 end do 558c 559c do isymdk = 1, nsym 560c n2bst2(isymdk) = 0 561c do isymk = 1, nsym 562c isymd = muld2h(isymdk,isymk) 563c n2bst2(isymdk) = n2bst2(isymdk) + mbas2(isymd)*mbas2(isymk) 564c end do 565c end do 566c 567c do isymdk = 1, nsym 568c icount = 0 569c do isymk = 1, nsym 570c isymd = muld2h(isymdk,isymk) 571c iaodis2(isymd,isymk) = icount 572c icount = icount + mbas2(isymd)*mbas2(isymk) 573c end do 574c end do 575c 576 do isympq = 1, nsym 577 nsjbp(isympq)= 0 578 nsbbp(isympq)= 0 579 do isymp = 1, nsym 580 isymq = muld2h(isympq,isymp) 581 isjbp(isymp,isymq) = nsjbp(isympq) 582 nsjbp(isympq) = nsjbp(isympq)+nrhfsa(isymp)*nbas2(isymq) 583 isbbp(isymp,isymq) = nsbbp(isympq) 584 nsbbp(isympq) = nsbbp(isympq)+nvir(isymp)*nbas2(isymq) 585 end do 586 end do 587c 588 kend0 = 1 589 590c---------------------------------------------------------------------- 591c get MO coefficients defined for combined orbital + aux. basis: 592c---------------------------------------------------------------------- 593c lusifc = -1 594c call gpopen(lusifc,'SIRIFC','OLD',' ','UNFORMATTED', 595c & idummy,.false.) 596c 597c read dimensions for CMO coefficients for full basis (nlamdsx) 598c rewind(lusifc) 599c call mollab('FULLBAS ',lusifc,lupri) 600c read(lusifc) nsymx,norbtsx,nbastx,nlamdsx,(nrhfsx(i),I=1,nsym), 601c & (norbsx(i),i=1,nsym) 602c 603c 604c allocate work space for CMO coefficients for full basis 605c kcmox = kend0 606c kend0 = kcmox + nlamdsx 607c lwrk0 = lwork - kend0 608c if (lwrk0 .lt.0) then 609c call quit('Insufficient work space in ccr12prep2') 610c end if 611c 612c read CMO coefficients and close file 613c read(lusifc) 614c read(lusifc) (work(kcmox+i-1),i=1,nlamdsx) 615c call gpclose(lusifc,'KEEP') 616c 617c test if CMO is ok 618c if (locdbg) then 619c write(lupri,*)'CMO out of ccr12prep2' 620c kscr1 = 1 621c do isymp = 1, nsym 622c isymmu = isymp 623c write(lupri,*) 'symmetry block,nbasx,norbsx:',isymp, 624c & nbasx(isymmu), norbsx(isymp) 625c call output(work(kcmox+kscr1-1),1,nbasx(isymmu), 626c & 1,norbsx(isymp),nbasx(isymmu),norbsx(isymp), 627c & 1,lupri) 628c kscr1 = kscr1 + nbasx(isymmu)*norbsx(isymp) 629c end do 630c end if 631 632c order CMO coefficients 633 klamdp = kend0 634 klamdh = klamdp + nlamdt 635 kend1 = klamdh + nlamdt 636 kt1am = kend1 637 kend2 = kt1am + nt1amx 638 lwrk2 = lwork - kend2 639 if(lwrk2.lt.0) then 640 call quit('insufficient work space in cc_r12prep2') 641 end if 642 call dzero(work(kt1am),nt1amx) 643 call lammat(work(klamdp),work(klamdh),work(kt1am), 644 & work(kend2),lwrk2) 645 646c----------------------------------------------------------------------- 647c read overlap matrix for full basis from file: 648c----------------------------------------------------------------------- 649 650 ksfull = kend1 651 kend1 = ksfull + nnbastx 652 lwrk1 = lwork - kend1 653 if (lwrk1 .lt.0) then 654 call quit('Insufficient work space in ccr12prep2') 655 end if 656 657 noauxbsave = noauxb 658 noauxb = .false. ! read integrals for the full basis 659 call rdonel('OVERLAP ',.true.,work(ksfull),nnbastx) 660 noauxb = noauxbsave 661 662 if (locdbg) then 663 !test if S is ok compare with sirius/r12aux.F ! 664 write(lupri,*)'in r12prep nnbastx = ', nnbastx 665 write(lupri,*)'in r12prep S' 666 call outpkb(work(ksfull),nbasx,1,1,lupri) 667 write(lupri,*)'norm^2 ksfull:', 668 & ddot(nnbastx,work(ksfull),1,work(ksfull),1) 669 end if 670 671c get S_Jb' and for ansatz 3 also S_bb' 672 ksjbp = kend1 673 kend1 = ksjbp + nsjbp(1) 674 675 ksbbp = kend1 676 if (ianr12.eq.3) then 677 kend1 = ksbbp + nsbbp(1) 678 do_virt = .true. 679 else 680 do_virt = .false. 681 end if 682 683 klamdhs = kend1 684 klamdps = klamdhs + nglmds(1) 685 kend1 = klamdps + nglmds(1) 686 lwrk1 = lwork - kend1 687 if (lwrk1 .lt.0) then 688 call quit('Insufficient work space in ccr12prep2') 689 end if 690 691 call lammats(work(klamdps),work(klamdhs),dummy, 692 & 1,.true.,.false., 693 & nglmds,iglmrhs,iglmvis,icmo,work(kend1),lwrk1) 694 695 call cc_r12mksjbp(work(ksfull),work(ksjbp),work(ksbbp),do_virt, 696 & work(klamdhs),nbasx,nbas2,iglmrhs,iglmvis, 697 & isjbp,isbbp,work(kend1),lwrk1) 698 699 if (locdbg) then 700 write(lupri,*)'norm^2 S_Jbp', 701 & ddot(nsjbp(1),work(ksjbp),1,work(ksjbp),1) 702 end if 703 704c---------------------------------------------------------------------- 705c make I^dp_a,kl 706c---------------------------------------------------------------------- 707c calculate new offsets: 708 do isymdk = 1, nsym 709 icount = 0 710 icoun1 = 0 711 do isymk = 1, nsym 712 isymd = muld2h(isymdk,isymk) 713 irgij(isymd,isymk) = icount 714 it1xao(isymd,isymk) = icoun1 715 icount = icount + mbas1(isymd)*nmatij(isymk) 716 icoun1 = icoun1 + mbas2(isymd)*nrhf(isymk) 717 end do 718 nrgij(isymdk) = icount 719 nt1xao(isymdk) = icoun1 720 end do 721 do isymdk = 1, nsym 722 icount = 0 723 do isymk = 1, nsym 724 isymd = muld2h(isymdk,isymk) 725 it2xaos(isymd,isymk) = icount 726 icount = icount + nt1ao(isymd)*nt1xao(isymk) 727 end do 728 end do 729c 730 do isymdp = 1, nsym 731 isymj = isymdp 732 isymd = isymdp 733 isymakl = isymd 734c 735 koffs = kend1 736 kghtf = koffs + mbas1(isymd)*nbas2(isymdp) 737 kgdpakl = kghtf + nrgij(isymakl) 738 kend2 = kgdpakl + nrgij(isymakl)*nbas2(isymd) 739 lwrk2 = lwork - kend2 740 if (lwrk2.lt.0) then 741 call quit('Insufficient work space in ccr12prep') 742 end if 743 744c get S_ddp 745 kcmo = klamdhs + iglmrhs(isymd,isymj) 746 koff1 = ksjbp + isjbp(isymdp,isymj) 747 ntotc = max(mbas1(isymd),1) 748 ntots = max(nrhfsa(isymj),1) 749 750 call dgemm('N','N',mbas1(isymd),nbas2(isymdp), 751 & nrhfsa(isymj),one,work(kcmo),ntotc,work(koff1), 752 & ntots,zero,work(koffs),ntotc) 753 754 if (ianr12.eq.3) then 755 kcmo = klamdhs + iglmvis(isymd,isymj) 756 koff1 = ksbbp + isbbp(isymj,isymdp) 757 ntots = max(nvir(isymj),1) 758 call dgemm('N','N',mbas1(isymd),nbas2(isymdp),nvir(isymj), 759 & one,work(kcmo),ntotc,work(koff1),ntots, 760 & one,work(koffs),ntotc) 761 end if 762c 763c get I^dp_akl 764 do idelp = 1, nbas2(isymdp) 765 lauxd = .false. 766 idxdp = idelp+ibas(isymdp)+basoff(isymdp) 767 koff2 = kgdpakl-1 + nrgij(isymakl)*(idelp-1)+1 768 call cc_r12getrint(work(koff2),idxdp,isymdp, 769 & nt1ao,it1ao,nt2aos,it2aos,nrgij,irgij,it1xao,it2xaos, 770 & nrhf,nmatij,imatij, 771 & ibasx,lauxd,.false.,fghtf,work(kend2),lwrk2) 772 if (locdbg) then 773 write(lupri,*)'norm^2 Idpakl:', 774 & ddot(nrgij(isymd),work(koff2),1,work(koff2),1) 775 end if 776 end do 777c add -I*S to I^dp_akl 778 do idel = 1, mbas1(isymd) 779 lauxd = .false. 780 idxd = idel+ibas(isymd) 781 call cc_r12getrint(work(kghtf),idxd,isymd,nt1ao,it1ao, 782 & nt2aos,it2aos,nrgij,irgij,it1xao,it2xaos, 783 & nrhf,nmatij,imatij, 784 & ibasx,lauxd,.false.,fghtf,work(kend2),lwrk2) 785 if (locdbg) then 786 write(lupri,*)'norm^2 Idakl:', 787 & ddot(nrgij(isymd),work(kghtf),1,work(kghtf),1) 788 end if 789 790 do idelp = 1, nbas2(isymdp) 791 idxddp = mbas1(isymd)*(idelp-1)+idel 792 koff2 = kgdpakl-1 + nrgij(isymakl)*(idelp-1)+1 793 call daxpy(nrgij(isymakl),-work(koffs-1+idxddp), 794 & work(kghtf),1,work(koff2),1) 795 end do 796 end do 797c save I^d'_akl on file 798 do idelp = 1, nbas2(isymdp) 799 koff2 = kgdpakl-1 + nrgij(isymakl)*(idelp-1)+1 800 lauxd = .false. 801 idxdp = idelp + ibas(isymdp) + basoff(isymdp) 802c 803 if (locdbg) then 804 write(lupri,*)'mbas1(isymd)', mbas1(isymd) 805 write(lupri,*) "save I^d'akl... isymdp,idelp,idxdp :", 806 & isymdp,idelp,idxdp 807 write(lupri,*)'norm^2 before write Idpakl:', 808 & ddot(nrgij(isymd),work(koff2),1,work(koff2),1) 809 write(lupri,*)"I^d'akl" 810 call output(work(koff2),1,nrgij(isymd),1,nrgij(isymd), 811 & nrgij(isymd),nrgij(isymd),1,lupri) 812 end if 813c 814 call cc_r12putrint(work(koff2),idxdp,isymdp,nt1ao,it1ao, 815 & nt2aos,it2aos,nrgij,irgij,it1xao,it2xaos, 816 & nrhf,nmatij,imatij, 817 & ibasx,lauxd,fgdpakl,work(kend2),lwrk2) 818 end do 819 end do 820c 821c---------------------------------------------------------------------- 822c calculate \sum_mn part for V_alphaj,kl 823c---------------------------------------------------------------------- 824c 825 kvajkl = kend1 826 kend2 = kvajkl + nvajkl(1) 827 lwrk2 = lwork - kend2 828 if (lwrk2.lt.0) then 829 call quit('Insufficient work space in prep2') 830 end if 831c 832c read V_alphajkl from file 833 luvajkl = -1 834 call gpopen(luvajkl,fvajkl,'unknown',' ','unformatted', 835 & idummy,.false.) 836 rewind(luvajkl) 837 read(luvajkl)(work(kvajkl+i-1),i=1,nvajkl(1)) 838 call gpclose(luvajkl,'KEEP') 839 840 if (locdbg) then 841 write(lupri,*)'Valphajkl after read in prep2:', 842 & (work(kvajkl+i-1),i=1,nvajkl(1)) 843 end if 844c 845 call cc_r12mkvaj2(work(kvajkl),1,work(klamdh),1, 846 & work(klamdhs),1,work(kend2),lwrk2) 847 if (locdbg) then 848 write(lupri,*)'Valphajkl out of prep2:', 849 & (work(kvajkl-1+i),i=1,nvajkl(1)) 850 end if 851c write V_alphajkl on file 852 luvajkl = -1 853 call gpopen(luvajkL,fvajkl,'unknown',' ','unformatted', 854 & idummy,.false.) 855 rewind(luvajkl) 856 write(luvajkl)(work(kvajkl+i-1),i=1,nvajkl(1)) 857 call gpclose(luvajkl,'KEEP') 858 859c test if V_alphajkl is ok 860 if (locdbg) then 861 write(lupri,*)'norm^2 Valphajkl in prep2:', 862 & ddot(nvajkl(1),work(kvajkl),1,work(kvajkl),1) 863 call cc_r12mkvijkl(work(kvajkl),1,work(klamdh),1, 864 & work(kend2),lwrk2,.false.,idummy,idummy) 865 write(lupri,*)'Vajkl for ansatz 2 finished in prep2' 866 end if 867 868c---------------------------------------------------------------------- 869c close file, clean trace and return: 870c---------------------------------------------------------------------- 871 if(locdbg) write(lupri,*)'leaving ccr12prep2' 872 call qexit('r12prep2') 873 874 end 875*======================================================================* 876 subroutine cc_r12mksjbp(smat,sjbp,sbbp,do_virt,cmo, 877 & nbasx,nbas2,iglmrhs,iglmvis, 878 & isjbp,isbbp,work,lwork) 879*----------------------------------------------------------------------* 880c purpose: get S_abp symmetry blocked out of S packed as a 881c triangular matrix and transform to S_Jbp and Sbbp 882c 883c H. Fliegl, C. Haettig, winter 2003 884c added transformation to Sbbp, Christof Haettig, spring 2005 885c adapted for CABS, Christian Neiss, winter 2006 886*----------------------------------------------------------------------* 887 implicit none 888#include "priunit.h" 889#include "ccorb.h" 890#include "ccsdsym.h" 891#include "r12int.h" 892 893 logical locdbg 894 parameter (locdbg = .false.) 895 896 logical do_virt 897 integer lwork,isjbp(8,8),iglmrhs(8,8),iglmvis(8,8),isbbp(8,8) 898 integer koff1, kcmo, ks, koffs, ntotb1, ntotbx, ntotrhs, ntotvir 899 integer nbasx(8), nbas2(8), ndim, isym 900 901#if defined (SYS_CRAY) 902 real work(*), smat(*), sjbp(*), sbbp(*), cmo(*) 903 real one, zero 904#else 905 double precision work(*), smat(*), sjbp(*), sbbp(*), cmo(*) 906 double precision one, zero 907#endif 908 parameter(zero = 0.0d0, one = 1.0d0) 909 910 call qenter('r12mksjbp') 911 912 koff1 = 1 913 do isym = 1, nsym 914 915 ndim = nbasx(isym) 916 if (lwork.lt.ndim*ndim) then 917 call quit('Insufficient work space in cc_r12mksjbp') 918 end if 919 920c -------------------------- 921c unpack triangular S matrix 922c -------------------------- 923 call sqmatr(ndim,smat(koff1),work) 924 koff1 = koff1 + ndim*(ndim+1)/2 925 926c -------------------------------------- 927c transform S(alpha,beta') to S(J,beta') 928c -------------------------------------- 929 kcmo = 1 + iglmrhs(isym,isym) 930 ks = 1 + isjbp(isym,isym) 931c koffs = 1 + nbasx(isym)*mbas1(isym) 932 koffs = 1 + nbasx(isym)*(nbasx(isym)-nbas2(isym)) 933 934 ntotb1 = max(mbas1(isym),1) 935 ntotbx = max(nbasx(isym),1) 936 ntotrhs = max(nrhfsa(isym),1) 937 938 call dgemm('T','N',nrhfsa(isym),nbas2(isym),mbas1(isym), 939 & one,cmo(kcmo),ntotb1,work(koffs),ntotbx, 940 & zero,sjbp(ks),ntotrhs) 941 942c -------------------------------------- 943c transform S(alpha,beta') to S(b,beta') 944c -------------------------------------- 945 if (do_virt) then 946 kcmo = 1 + iglmvis(isym,isym) 947 ks = 1 + isbbp(isym,isym) 948 949 ntotvir = max(nvir(isym),1) 950 951 call dgemm('T','N',nvir(isym),nbas2(isym),mbas1(isym), 952 & one,cmo(kcmo),ntotb1,work(koffs),ntotbx, 953 & zero,sbbp(ks),ntotvir) 954 end if 955 956 end do 957 958 call qexit('r12mksjbp') 959 end 960*=====================================================================* 961 subroutine cc_r12mksmupbp(cmo,smupbp,iaodis2,work,lwork) 962c--------------------------------------------------------------------- 963c purpose: make D_mu'b' = \sum_p' C_b'p' * C_mu'p' for ansatz 2 964c for CABS: sum over MO and aux. basis 965c 966c H. Fliegl, C. Haettig, winter 2003 967c C. Neiss, winter 2006: adapted for CABS 968c--------------------------------------------------------------------- 969 implicit none 970#include "priunit.h" 971#include "ccorb.h" 972#include "ccsdsym.h" 973#include "r12int.h" 974 logical locdbg 975 parameter (locdbg = .false.) 976 977 integer lwork, isymp, isymq, isympq, isymbp, isym, 978 & isymmup, kcmo, ntotc, ntots, 979 & iaodis2(8,8), ks, isympbp 980 integer idxbpp, idxbpp2, ip, ib, koff0, idxp, idxb 981 integer nbasx(8),norbx(8),basoff(8),orboff(8) 982 983#if defined (SYS_CRAY) 984 real cmo(*), smupbp(*), work(*) 985#else 986 double precision cmo(*), smupbp(*), work(*) 987#endif 988 double precision one, zero 989 parameter(one = 1.0d0, zero = 0.0d0) 990 991 call qenter('r12mksmupbp') 992c 993 if (r12cbs) then 994 do isym = 1, nsym 995 nbasx(isym) = mbas1(isym)+mbas2(isym) 996 norbx(isym) = norb1(isym)+norb2(isym) 997 basoff(isym) = 0 998 orboff(isym) = 0 999 end do 1000 else 1001 call icopy(8,mbas2,1,nbasx,1) 1002 call icopy(8,norb2,1,norbx,1) 1003 call icopy(8,mbas1,1,basoff,1) 1004 call icopy(8,norb1,1,orboff,1) 1005 end if 1006c 1007 koff0 = 1 1008 do isymp = 1, nsym 1009 isymbp = isymp 1010 isymmup = isymp 1011c isympbp = muld2h(isymp,isymbp) ! eq 1 1012c get C_p'b' out of CMO 1013 if (locdbg) then 1014 write(lupri,*)'mbas1(isymbp),mbas2(isymbp)', 1015 & mbas1(isymbp),mbas2(isymbp) 1016 write(lupri,*)'norb1(isymp),norb2(isymp)', 1017 & norb1(isymp),norb2(isymp) 1018 write(lupri,*)'norb2(isymp)*mbas2(isymbp)', 1019 & norb2(isymp)*mbas2(isymbp) 1020 write(lupri,*)'nrhfsb(isymp)',nrhfsb(isymp) 1021 write(lupri,*)'isymbp, isymp', isymbp, isymp 1022 end if 1023 do ip=orboff(isymp)+1+nrhfsb(isymp), 1024 & (norb1(isymp)+norb2(isymp)+nrhfsb(isymp)) 1025 idxp = ip - nrhfsb(isymp) - orboff(isymp) 1026 do ib = basoff(isymbp)+1,(mbas1(isymbp)+mbas2(isymbp)) 1027 idxb = ib - basoff(isymp) 1028 idxbpp = koff0-1 + (mbas1(isymbp)+mbas2(isymbp))*(ip-1)+ib 1029 idxbpp2 = nbasx(isymbp)*(idxp-1)+idxb 1030 work(idxbpp2) = cmo(idxbpp) 1031 if (locdbg) then 1032 write(lupri,*)'wrk,idxbpp,idxbpp2',idxbpp, idxbpp2 1033 write(lupri,*)work(idxbpp2) 1034 end if 1035 end do 1036 end do 1037 koff0 = koff0 + (mbas1(isymbp)+mbas2(isymbp)) * 1038 & (nrhfsb(isymp)+norb1(isymp)+norb2(isymp)) 1039 1040c transform to D 1041 1042 ks = 1+iaodis2(isymbp,isymmup) 1043 ntotc = max(nbasx(isymbp),1) 1044 1045 call dgemm('N','T',nbasx(isymbp),nbasx(isymmup),norbx(isymp), 1046 & one,work,ntotc,work,ntotc,zero, 1047 & smupbp(ks),ntotc) 1048 if (locdbg) then 1049 write(lupri,*)'resorted CMO in cc_r12mksmupbp' 1050 call output(work,1,nbasx(isymbp),1, 1051 & norbx(isymp),nbasx(isymbp),norbx(isymp),1,lupri) 1052 write(lupri,*)'S_bp,mup in cc_r12mksmupbp' 1053 call output(smupbp(ks),1,nbasx(isymbp),1, 1054 & nbasx(isymmup),nbasx(isymbp),nbasx(isymmup),1,lupri) 1055 end if 1056 end do 1057 1058 call qexit('r12mksmupbp') 1059 end 1060*=====================================================================* 1061 subroutine cc_r12getkiajb(xkipjq,xkiajb,oneaux,nt2amx, 1062 & nh2amx,ih2am,nh1am,ih1am,norb1, 1063 & nrhf1,nrhfs1) 1064c--------------------------------------------------------------------- 1065c purpose: get K^ij_ab out of K^ij_pq stored as a lower 1066c symmetry packed triangular matrix 1067c 1068c H. Fliegl, C. Haettig, spring 2004 1069c--------------------------------------------------------------------- 1070 implicit none 1071#include "priunit.h" 1072#include "ccorb.h" 1073 logical locdbg,oneaux 1074 parameter (locdbg = .false.) 1075 integer nh2amx, ih2am(8,8), nh1am(8), ih1am(8,8), nvircc(8), 1076 & nt2amx, nt1am(8), it1am(8,8), it2am(8,8), isympi, 1077 & isymi, isym, icount1, icount2, isymj, isymk, ldr12(8), 1078 & isymqj, isymq, isymp, nqj, npi, nbj, nai, idxpiqj, 1079 & idxaibj, index, i,j,p,a,b,q, norb1(8),nrhf1(8),nrhfs1(8) 1080#if defined (SYS_CRAY) 1081 real xkiajb(*), xkipjq(*) 1082#else 1083 double precision xkiajb(*), xkipjq(*), ddot 1084#endif 1085 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 1086 1087 call qenter('getkiajb') 1088 1089c calculate some offsets and dimensions 1090c calculate number of active virtual orbitals in CC 1091c nvirfr = frozen virtual orbitals 1092 do isym = 1, nsym 1093 nvircc(isym) = norb1(isym) - nrhfsa(isym) - nvirfr(isym) 1094 end do 1095 1096 if (locdbg) then 1097 write(lupri,* ) 'isym nrhfs1 nrhf1 nvircc norb1' 1098 do isym = 1, nsym 1099 write(lupri,*) isym,nrhfs1(isym),nrhf1(isym), 1100 & nvircc(isym),norb1(isym) 1101 end do 1102 end if 1103 1104 nt2amx = 0 1105 do isympi = 1, nsym 1106 nt1am(isympi) = 0 1107 do isymi = 1, nsym 1108 isymp = muld2h(isympi,isymi) 1109 nt1am(isympi) = nt1am(isympi) + nvircc(isymp)*nrhf1(isymi) 1110 end do 1111 nt2amx = nt2amx + nt1am(isympi)*(nt1am(isympi)+1)/2 1112 end do 1113 1114 do isymk = 1, nsym 1115 icount1 = 0 1116 icount2 = 0 1117 do isymj = 1, nsym 1118 isymi = muld2h(isymj,isymk) 1119 it1am(isymi,isymj) = icount1 1120 icount1 = icount1 + nrhf1(isymj)*nvircc(isymi) 1121 if (isymj .gt. isymi) then 1122 it2am(isymi,isymj) = icount2 1123 it2am(isymj,isymi) = icount2 1124 icount2 = icount2 + nt1am(isymi)*nt1am(isymj) 1125 else if (isymk .eq. 1) then 1126 it2am(isymi,isymj) = icount2 1127 icount2 = icount2 + nt1am(isymi)*(nt1am(isymi)+1)/2 1128 end if 1129 end do 1130 end do 1131 1132 do isym = 1, nsym 1133 if (oneaux) then 1134 ldr12(isym) = norb1(isym) 1135 else 1136 ldr12(isym) = nvir(isym) 1137 end if 1138 end do 1139c 1140 if (locdbg) then 1141 write(lupri,*) 'norm^2(kpiqj):',ddot(nh2amx,xkipjq,1,xkipjq,1) 1142 end if 1143c 1144c --------------------------------------------------------------- 1145c loop over symmetries and indices pick out p,q as virtual MO's 1146c and reorder integrals 1147c 1148c reorder integrals 1149c npi : pair p,i with i occupied (r12 orbital) and p MO 1150c virtual dimensions as needed for 1151c MP2-R12 code 1152c nai : pair a,i dimensions as needed for 1153c CC code 1154c idxpiqj: quadrupel pi,qj dimensions as 1155c needed for MP2-R12 code 1156c idxaibj: quadrupel ai,bj dimensions as 1157c needed for CC code 1158c --------------------------------------------------------------- 1159 do isympi = 1, nsym 1160 isymqj = isympi 1161 do isymp = 1, nsym 1162 isymi = muld2h(isympi,isymp) 1163 do isymq =1, nsym 1164 isymj = muld2h(isymqj,isymq) 1165c 1166 do i = 1, nrhf1(isymi) 1167 do p = nrhfs1(isymp)+1, nrhfs1(isymp)+nvircc(isymp) 1168 a = p - nrhfs1(isymp) 1169 do j = 1, nrhf1(isymj) 1170 do q = nrhfs1(isymq)+1, nrhfs1(isymq)+ 1171 & nvircc(isymq) 1172 b = q - nrhfs1(isymq) 1173 nqj = ih1am(isymq,isymj)+ ldr12(isymq)*(j-1)+q 1174 npi = ih1am(isymp,isymi)+ ldr12(isymp)*(i-1)+p 1175 nbj = it1am(isymq,isymj)+nvircc(isymq)*(j-1)+b 1176 nai = it1am(isymp,isymi)+nvircc(isymp)*(i-1)+a 1177 idxpiqj = ih2am(isympi,isymqj)+ index(npi,nqj) 1178 idxaibj = it2am(isympi,isymqj)+ index(nai,nbj) 1179 xkiajb(idxaibj) = xkipjq(idxpiqj) 1180 end do 1181 end do 1182 end do 1183 end do 1184c 1185 end do 1186 end do 1187 end do 1188 1189 if (locdbg) then 1190 write(lupri,*)'xkiajb' 1191 do isymi =1, nsym 1192 call outpak(xkiajb(1+it2am(isymi,isymi)),nt1am(isymi),1,lupri) 1193 end do 1194 write(lupri,*) 'norm^2(kiajb):',ddot(nt2amx,xkiajb,1,xkiajb,1) 1195 end if 1196c 1197 call qexit('getkiajb') 1198 end 1199*======================================================================* 1200 subroutine ccrhs_epp(omega,t2am,isymt2am,work,lwork,aproxr12,lres, 1201 & lufc2,fc2am,ifile) 1202*----------------------------------------------------------------------* 1203c purpose: add E'' term to Omega^ij_kl (stored in omega) for the 1204c ansaetze 2 and 3 1205c 1206c omega^ij_kl <-- omega^ij_kl 1207c + sum_ab (K^ab_kl - t^ab_kl) * T^ab_ij 1208c 1209c the contrib. of K is skipped for approximation A 1210c 1211c for the approximations ??? is in ansatz 2 added: 1212c 1213c omega^ij_kl <-- omega^ij_kl 1214c + sum_ab (e_k+e_l-e_i-e_j) r^ab_kl * T^ab_ij 1215c 1216c while in ansatz 3 the following term is added: 1217c 1218c omega^ij_kl <-- omega^ij_kl 1219c + sum_ab (e_k+e_l-e_a-e_b) r^ab_ij * T^ab_ij 1220c 1221c 1222c H. Fliegl, C. Haettig, spring 2004 1223c adapted for ansatz 3, Christof Haettig, spring 2005 1224*----------------------------------------------------------------------* 1225 implicit none 1226#include "priunit.h" 1227#include "ccorb.h" 1228#include "ccsdsym.h" 1229#include "r12int.h" 1230#include "ccr12int.h" 1231 logical locdbg,lres 1232 parameter (locdbg = .false.) 1233 integer lwork,lunit,isymt2am,isymkint,kt2pck,kpck,lwrk1,idum, 1234 & kend1,iopt,idummy,kap,lufc2,ifile,kend2,lwrk2 1235 character*(*) fc2am 1236#if defined (SYS_CRAY) 1237 real work(*), t2am(*), omega(*),factor,ddot,two,one 1238#else 1239 double precision work(*), t2am(*), omega(*),factor,ddot,two,one 1240#endif 1241 parameter(two=2.0D0, one=1.0D0) 1242 character*(*) aproxr12 1243 1244 call qenter('rhs_epp') 1245 1246 kt2pck = 1 1247 kpck = kt2pck + max(nt2am(isymt2am),nt2r12(isymt2am),nt2r12(1)) ! if lres isymt2am = isymtr! 1248 kend1 = kpck + nt2r12(1) 1249 lwrk1 = lwork - kend1 1250 1251 kend2 = kend1 1252 lwrk2 = lwork - kend2 1253 1254 if (lwrk1.lt.0 .or. lwrk2.lt.0) then 1255 call quit('Insufficient work space in ccrhs_epp') 1256 end if 1257 1258c --------------------------------------------- 1259c read t^ab_kl integrals from file 1260c --------------------------------------------- 1261 lunit = -1 1262 call gpopen(lunit,fkr12,'unknown',' ','unformatted',idum,.false.) 1263 read(lunit)(work(kpck-1+i),i=1,nt2r12(1)) 1264 call gpclose(lunit,'KEEP') 1265 1266c ---------------------------------------------------- 1267c for approximation B substract K^ab_kl from t^ab_kl 1268c ---------------------------------------------------- 1269 if (aproxr12(1:1).eq.'B') then 1270 lunit = -1 1271 call gpopen(lunit,ftr12,'unknown',' ','unformatted', 1272 & idum,.false.) 1273 read(lunit)(work(kt2pck+i-1),i=1,nt2r12(1)) 1274 call gpclose(lunit,'KEEP') 1275 1276 ! substract (t - K) => work(kpck) 1277 call daxpy(nt2r12(1),-1.0d0,work(kt2pck),1,work(kpck),1) 1278 end if 1279 1280c ------------------------------------------------------------ 1281c in ansatz 3 add for approximations A' and B the term 1282c (eps_k + eps_l - eps_a - eps_b) * r^ab_kl 1283c ------------------------------------------------------------ 1284 if (ianr12.eq.3 .and. 1285 & (aproxr12(1:2).eq."A'".or.aproxr12(1:1).eq."B") ) then 1286 lunit = -1 1287 call gpopen(lunit,fr12r12,'unknown',' ','unformatted', 1288 & idum,.false.) 1289 read(lunit)(work(kt2pck+i-1),i=1,nt2r12(1)) 1290 call gpclose(lunit,'KEEP') 1291 call cc_t2xe(work(kt2pck),work(kend2),lwrk2) 1292 call daxpy(nt2r12(1),+1.0d0,work(kt2pck),1,work(kpck),1) 1293 end if 1294 1295c ---------------------------------------------------------- 1296c get doubles amplitudes R^mn_ab stored as triangular matrix 1297c ---------------------------------------------------------- 1298 if (lres) then 1299 call cc_rvec(lufc2,fc2am,nt2am(isymt2am),nt2am(isymt2am), 1300 & ifile,work(kt2pck)) 1301 call cclr_diascl(work(kt2pck),two,isymt2am) 1302 else 1303 iopt = 0 1304 call cc_t2pk(work(kt2pck),t2am,isymt2am,iopt) 1305 end if 1306 1307c ------------------------------------------------------ 1308c calculate Omega^ij_kl = \sum_ab k^ab_kl * t^ij_ab 1309c ------------------------------------------------------ 1310 if (locdbg) then 1311 write(lupri,*) 'in ccrhs_epp before cc_r12mi2:' 1312 write(lupri,*) 'work(kt2pck):' 1313 call outpak(work(kt2pck),nt1am(1),1,lupri) 1314 end if 1315 isymkint = 1 1316 factor = 1.0d0 1317 call cc_r12mi2(omega,work(kpck),work(kt2pck),isymkint,isymt2am, 1318 & factor,work(kend1),lwrk1) 1319 1320 if (locdbg) then 1321 write(lupri,*)'Norm^2 Omega^ij_kl in ccrhs_epp', 1322 & ddot(ntr12sq(isymt2am),omega,1,omega,1) 1323 end if 1324 1325c ------------------------------------------------------------ 1326c in ansatz 2 add for approximation A' and B the term 1327c omega^ij_kl += sum_ab (e_k+e_l-e_i-e_j) r^ab_kl * c^kl_ij 1328c ------------------------------------------------------------ 1329 if (ianr12.eq.2 .and. 1330 & (aproxr12(1:2).eq."A'".or.aproxr12(1:1).eq."B") ) then 1331 kap = kend1 1332 kend1 = kap + ntr12sq(isymt2am) 1333 lwrk1 = lwork - kend1 1334 if (lwrk1.lt.0) then 1335 call quit('Insufficient work space in ccrhs_epp') 1336 end if 1337 1338 lunit = -1 1339 call gpopen(lunit,fr12r12,'unknown',' ','unformatted', 1340 & idum,.false.) 1341 read(lunit)(work(kpck-1+i),i=1,nt2r12(1)) 1342 call gpclose(lunit,'KEEP') 1343 1344 ! calculate \sum_ab r^ab_kl*t^ij_ab 1345 factor = one 1346 isymkint = 1 1347 call dzero(work(kap),ntr12sq(isymt2am)) 1348 call cc_r12mi2(work(kap),work(kpck),work(kt2pck),isymkint, 1349 & isymt2am,factor,work(kend1),lwrk1) 1350 1351 call cc_r12mkediff(work(kap),isymt2am,work(kend1),lwrk1) 1352 1353 call daxpy(ntr12sq(isymt2am),1.0d0,work(kap),1,omega,1) 1354 end if 1355 1356c ----------------------------------------------------------------- 1357c that's it; restore T2 amplitudues, print debug output and return: 1358c ----------------------------------------------------------------- 1359 if (.not.lres) then 1360c pack T2 amplitudes back to a quadratic matrix 1361 call cc_t2sq(work(kt2pck),t2am,isymt2am) 1362 end if 1363 1364 if (locdbg) then 1365 write(lupri,*) 'omega^ij_kl after ccrhs_epp contribution:' 1366 write(lupri,*) 'norm^2 omega^ij,kl', 1367 & ddot(ntr12sq(isymt2am),omega,1,omega,1) 1368 end if 1369 1370 call qexit('rhs_epp') 1371 end 1372*=====================================================================* 1373 subroutine cc_r12mkediff(r12int,isymr,work,lwork) 1374c--------------------------------------------------------------------- 1375c purpose: calculate r12int*(e_k+e_l-e_i-e_j) 1376c 1377c H. Fliegl, C. Haettig, spring 2004 1378c modified C. Neiss, september 2005 1379c--------------------------------------------------------------------- 1380 implicit none 1381#include "priunit.h" 1382#include "ccorb.h" 1383#include "ccsdsym.h" 1384#include "ccsdinp.h" 1385#include "inftap.h" 1386 integer lwork, kscr1,kend1,lwrk1,isymkl,isymij,keij,kekl, 1387 & idxijkl,isymr,ij,kl,idummy,norbtsx 1388 integer nij,nkl,inmatij(8),inmatkl(8),isym 1389#if defined (SYS_CRAY) 1390 real work(*), r12int(*) 1391#else 1392 double precision work(*), r12int(*) 1393#endif 1394c 1395 call qenter('r12mkediff') 1396c 1397c------------------------------------- 1398c Read canonical orbital energies. 1399c------------------------------------- 1400c 1401 call gpopen(lusifc,'SIRIFC','UNKNOWN',' ','UNFORMATTED',idummy, 1402 & .false.) 1403 rewind lusifc 1404c 1405 call mollab('FULLBAS ',lusifc,lupri) 1406 read (lusifc) idummy,norbtsx 1407 kscr1 = 1 1408 kend1 = kscr1 + norbtsx 1409 lwrk1 = lwork - kend1 1410 if (lwrk1.lt.0) then 1411 call quit('Insufficient work space in cc_r12mkediff') 1412 end if 1413 read (lusifc) (work(kscr1+i-1), i=1,norbtsx) 1414 call gpclose(lusifc,'KEEP') 1415c 1416 nkl = 0 1417 nij = 0 1418 do isym = 1, nsym 1419 inmatkl(isym) = nkl 1420 inmatij(isym) = nij 1421 nkl = nkl + nmatkl(isym) 1422 nij = nij + nmatij(isym) 1423 end do 1424c 1425 keij = kend1 1426 kekl = keij + nij 1427 kend1 = kekl + nkl 1428 lwrk1 = lwork - kend1 1429 if (lwrk1.lt.0) then 1430 call quit('Insufficient work space in cc_r12mkediff') 1431 end if 1432c 1433 call cc_r12epair(work(kekl),nkl,work(kscr1),inmatkl,imatkl,nrhfb) 1434 call cc_r12epair(work(keij),nij,work(kscr1),inmatij,imatij,nrhf) 1435c 1436 do isymkl = 1, nsym 1437 isymij = muld2h(isymr,isymkl) 1438 do kl = 1, nmatkl(isymkl) 1439 do ij = 1, nmatij(isymij) 1440 idxijkl = itr12sqt(isymij,isymkl)+nmatij(isymij)*(kl-1)+ij 1441 r12int(idxijkl) = r12int(idxijkl)* 1442 & (work(kekl-1+kl+inmatkl(isymkl))- 1443 & work(keij-1+ij+inmatij(isymij))) 1444 end do 1445 end do 1446 end do 1447 1448 call qexit('r12mkediff') 1449 end 1450*=====================================================================* 1451 subroutine cc_r12mkediff2(c2am,isymc,work,lwork) 1452c--------------------------------------------------------------------- 1453c purpose: calculate c12*(e_k+e_l-e_i-e_j) 1454c 1455c H. Fliegl, C. Haettig, spring 2004 1456c--------------------------------------------------------------------- 1457 implicit none 1458#include "priunit.h" 1459#include "ccorb.h" 1460#include "ccsdsym.h" 1461#include "ccsdinp.h" 1462#include "inftap.h" 1463 integer lwork, kscr1,kend1,lwrk1, isymc,index,idummy, 1464 & kend2,lwrk2,isymk,isymi,isyml,isymj, 1465 & mk,mi,mj,ml 1466 integer isymki,isymlj,keki,kelj,lj,ki,idxkilj 1467#if defined (SYS_CRAY) 1468 real work(*), c2am(*),el,ei,ek,ej 1469#else 1470 double precision work(*), c2am(*),el,ei,ek,ej 1471#endif 1472 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 1473 1474 call qenter('r12mkediff2') 1475 1476 kscr1 = 1 1477 kend1 = kscr1 + norbts 1478 lwrk1 = lwork - kend1 1479 if (lwrk1.lt.0) then 1480 call quit('Insufficient work space in cc_r12mkediff2') 1481 end if 1482c------------------------------------- 1483c Read canonical orbital energies. 1484c------------------------------------- 1485c 1486 call gpopen(lusifc,'SIRIFC','UNKNOWN',' ','UNFORMATTED',idummy, 1487 & .false.) 1488 rewind lusifc 1489c 1490 call mollab('TRCCINT ',lusifc,lupri) 1491 read (lusifc) 1492 read (lusifc) (work(kscr1+i-1), i=1,norbts) 1493 call gpclose(lusifc,'KEEP') 1494c 1495 if (froimp .or. froexp) 1496 * call ccsd_delfro(work(kscr1),work(kend1),lwrk1) 1497c 1498 do isymki = 1, nsym 1499 isymlj = muld2h(isymc,isymki) 1500 if (isymki.le.isymlj) then 1501 keki = kend1 1502 kelj = keki + nmatki(isymki) 1503 kend2 = kelj + nmatki(isymlj) 1504 lwrk2 = lwork - kend2 1505 if (lwrk2.lt.0) then 1506 call quit('Insufficient work space in cc_r12mkediff2') 1507 end if 1508 do isymk = 1, nsym 1509 isymi = muld2h(isymki,isymk) 1510 do isyml = 1, nsym 1511 isymj = muld2h(isymlj,isyml) 1512 do k = 1, nrhfb(isymk) 1513 mk = iorb(isymk)+k 1514 ek = work(kscr1+mk-1) 1515 do i = 1, nrhfa(isymi) 1516 ki = imatki(isymk,isymi)+nrhfb(isymk)*(i-1)+k 1517 mi = iorb(isymi)+i 1518 ei = work(kscr1+mi-1) 1519 work(keki+ki-1) = ek - ei 1520 end do 1521 end do 1522 do l = 1, nrhfb(isyml) 1523 ml = iorb(isyml) +l 1524 el = work(kscr1+ml-1) 1525 do j = 1, nrhfa(isymj) 1526 lj = imatki(isyml,isymj)+nrhfb(isyml)*(j-1)+l 1527 mj = iorb(isymj)+j 1528 ej = work(kscr1+mj-1) 1529 work(kelj+lj-1) = el - ej 1530 end do 1531 end do 1532 end do 1533 end do 1534c 1535 if (isymki.eq.isymlj) then 1536 do lj = 1, nmatki(isymlj) 1537 do ki = 1, lj 1538 idxkilj = itr12am(isymki,isymlj)+index(ki,lj) 1539 c2am(idxkilj) = c2am(idxkilj) * 1540 & ( work(keki-1+ki) + work(kelj-1+lj) ) 1541 end do 1542 end do 1543 else 1544 do lj = 1, nmatki(isymlj) 1545 do ki = 1, nmatki(isymki) 1546 idxkilj = itr12am(isymki,isymlj)+nmatki(isymki)*(lj-1)+ 1547 & ki 1548 c2am(idxkilj) = c2am(idxkilj)* 1549 & ( work(keki-1+ki) + work(kelj-1+lj) ) 1550 end do 1551 end do 1552 end if 1553 end if 1554 end do 1555 1556 call qexit('r12mkediff2') 1557 return 1558 end 1559*======================================================================* 1560 subroutine cc_t2xe(t2am,work,lwork) 1561*----------------------------------------------------------------------* 1562* Purpose: multiply vector on t2am in place with orbital 1563* energy differences 1564* 1565* t2(ak,bl) <-- t2(ak,bl) * (e_k + e_l - e_a - e_b) 1566* 1567* amplitudes are addressed as triangle 1568* 1569* C. Haettig, spring 2005 1570* modified for general r12 orbitals: C. Neiss, september 2005 1571*----------------------------------------------------------------------* 1572 implicit none 1573#include "priunit.h" 1574#include "ccsdinp.h" 1575#include "ccorb.h" 1576#include "ccsdsym.h" 1577 1578 logical locdbg, delfro 1579 parameter (locdbg = .false.) 1580 parameter (delfro = .true.) 1581 1582 integer isymbj, isymai, isymj, isymi, isymb, isyma, index, 1583 & koffa, koffb, nbj, nai, naibj, lwork, idum 1584 integer nkl,inmatkl(8),iak(8,8),icount,kekl,kend1,lwrk1,lunit, 1585 & koffij,isymij,keps,norbtsx 1586#if defined (SYS_CRAY) 1587 real t2am(*), work(lwork) 1588#else 1589 double precision t2am(*), work(lwork) 1590#endif 1591 1592 index(i,j) = max(i,j)*(max(i,j) - 3)/2 + i + j 1593 1594 call qenter('cc_t2xe') 1595 1596 nkl = 0 1597 do isymij = 1, nsym 1598 inmatkl(isymij) = nkl 1599 nkl = nkl + nmatkl(isymij) 1600 end do 1601c 1602 do isymai = 1, nsym 1603 icount = 0 1604 do isymi = 1, nsym 1605 isyma = muld2h(isymai,isymi) 1606 iak(isyma,isymi) = icount 1607 icount = icount + nvir(isyma)*nrhfb(isymi) 1608 end do 1609 end do 1610c 1611 kekl = 1 1612 kend1 = kekl + nkl 1613 lwrk1 = lwork - kend1 1614 if (lwrk1.lt.0) call quit('Insufficient memory in CC_T2XE') 1615c 1616 lunit = -1 1617 call gpopen(lunit,'SIRIFC','OLD',' ','UNFORMATTED',idum, 1618 & .false.) 1619 rewind lunit 1620 call mollab('FULLBAS ',lunit,lupri) 1621 read (lunit) idum,norbtsx 1622 keps = kend1 1623 kend1 = keps + norbtsx 1624 lwrk1 = lwork - kend1 1625 if (lwrk1.lt.0) call quit('Insufficient memory in CC_T2XE') 1626 read (lunit) (work(keps+i-1), i=1,norbtsx) 1627c call gpclose(lunit,'KEEP') 1628c 1629 call cc_r12epair(work(kekl),nkl,work(keps),inmatkl,imatkl,nrhfb) 1630c 1631 rewind lunit 1632 call mollab('TRCCINT ',lunit,lupri) 1633 read (lunit) 1634 read (lunit) (work(keps+i-1), i=1,norbts) 1635 call gpclose(lunit,'KEEP') 1636c 1637 if (froimp .or. froexp) 1638 * call ccsd_delfro(work(keps),work(kend1),lwrk1) 1639c 1640 call fock_reorder(work(keps),work(kend1),lwrk1) 1641c 1642 do isymbj = 1,nsym 1643 isymai = isymbj 1644 do isymj = 1,nsym 1645 isymb = muld2h(isymj,isymbj) 1646 do isymi = 1,nsym 1647 isyma = muld2h(isymi,isymai) 1648 do j = 1,nrhfb(isymj) 1649 do b = 1,nvir(isymb) 1650 koffb = ivir(isymb) + b 1651 nbj = iak(isymb,isymj)+nvir(isymb)*(j-1)+b 1652 do i = 1,nrhfb(isymi) 1653 do a = 1,nvir(isyma) 1654 koffa = ivir(isyma) + a 1655 nai = iak(isyma,isymi)+nvir(isyma)*(i-1)+a 1656 1657 if (nai .le. nbj) then 1658 isymij = muld2h(isymi,isymj) 1659 koffij = inmatkl(isymij)+imatkl(isymi,isymj)+ 1660 & nrhfb(isymi)*(j-1)+i-1 1661 naibj = it2r12(isymai,isymbj) + index(nai,nbj) 1662 t2am(naibj) = t2am(naibj) * 1663 & (work(kekl+koffij)-work(keps+koffa-1)- 1664 & work(keps+koffb-1)) 1665 end if 1666 1667 end do 1668 end do 1669 end do 1670 end do 1671 end do 1672 end do 1673 end do 1674 1675 call qexit('cc_t2xe') 1676 return 1677 end 1678*=====================================================================* 1679 subroutine cc_r12mi2(omega,xkint,t2am,isymkint,isymt2am, 1680 & factor,work,lwork) 1681c--------------------------------------------------------------------- 1682c purpose: calculate Omega^ij_kl = \sum_ab k^ab_kl * t^ij_ab 1683c input: xkint and t2am packed as a triangular matrix 1684c 1685c H. Fliegl, C. Haettig, spring 2004 1686c--------------------------------------------------------------------- 1687 implicit none 1688#include "priunit.h" 1689#include "ccorb.h" 1690#include "ccsdsym.h" 1691 logical locdbg 1692 parameter (locdbg = .false.) 1693 integer isymkint, isymt2am, lwork, ntoto, ntotk, 1694 & isyma, isymij, isymkl, kxint, kt2am, kend1, lwrk1, 1695 & isymab, isymint, isymb, kom 1696#if defined (SYS_CRAY) 1697 real work(*), xkint(*),t2am(*), omega(*),factor 1698 real one,ddot,xk,xt 1699#else 1700 double precision work(*), xkint(*),t2am(*), omega(*),factor 1701 double precision one,ddot,xk,xt 1702#endif 1703 parameter (one = 1.0d0) 1704 1705 call qenter('r12mi2') 1706c 1707 isymint = muld2h(isymkint,isymt2am) !isymint = isymomg 1708 1709 if (locdbg) then 1710 xk = 0.0d0 1711 xt = 0.0d0 1712 end if 1713 do isymb = 1, nsym 1714 do isyma = 1, nsym 1715 isymab = muld2h(isyma,isymb) 1716 isymkl = muld2h(isymkint,isymab) 1717 isymij = muld2h(isymkl,isymint) 1718c 1719 kxint = 1 1720 kt2am = kxint + nvir(isyma)*nmatkl(isymkl) 1721 kend1 = kt2am + nvir(isyma)*nmatij(isymij) 1722 lwrk1 = lwork - kend1 1723 if (lwrk1.lt.0) then 1724 call quit('Insufficient work space in cc_r12mi2') 1725 end if 1726c 1727 do b = 1, nvir(isymb) 1728 1729C call cc_r12sortk(work(kxint),xkint,b,isymkint,isymb, 1730C & isyma,nrhfb,.false.) 1731 call cc_r12sort(xkint,isymkint,work(kxint),b,isymb, 1732 & 1,nvir(isyma),isyma,nvir,nrhfb) 1733 if (locdbg) then 1734 xk = xk + 1735 & ddot(nvir(isyma)*nmatkl(isymkl),work(kxint),1, 1736 & work(kxint),1) 1737 end if 1738 1739C call cc_r12sortk(work(kt2am),t2am,b,isymt2am,isymb, 1740C & isyma,nrhf,.false.) 1741 call cc_r12sort(t2am,isymt2am,work(kt2am),b,isymb, 1742 & 1,nvir(isyma),isyma,nvir,nrhf) 1743 if (locdbg) then 1744 xt = xt + 1745 & ddot(nvir(isyma)*nmatij(isymij),work(kt2am),1, 1746 & work(kt2am),1) 1747 end if 1748c 1749 ntotk = max(nvir(isyma),1) 1750 ntoto = max(nmatij(isymij),1) 1751c 1752 kom = itr12sqt(isymij,isymkl) +1 1753 call dgemm('T','N',nmatij(isymij),nmatkl(isymkl), 1754 & nvir(isyma),factor,work(kt2am),ntotk,work(kxint), 1755 & ntotk,one,omega(kom),ntoto) 1756 end do 1757 end do 1758 end do 1759 1760 if (locdbg) then 1761 write(lupri,*)'xk =', xk 1762 write(lupri,*)'xt =', xt 1763 end if 1764 1765 call qexit('r12mi2') 1766 end 1767*=====================================================================* 1768 subroutine cc_r12sortk(xk,xkint,b,isymkint,isymb,isyma,nr12orb, 1769 & lopt) 1770c--------------------------------------------------------------------- 1771c purpose: sort k^ab_kl as three index array with fixed b 1772c expect that k is stored as a lower triangular 1773c matrix 1774c lopt = .false. : xkint stored with it2am/it1am 1775c lopt = .true. : xkint stored with ih2am/ih1am 1776c nr12orb = # of k (nrhf or nrhfb) 1777c 1778c H. Fliegl, C. Haettig, spring 2004 1779c modified C. Neiss august 2005 1780c--------------------------------------------------------------------- 1781 implicit none 1782#include "priunit.h" 1783#include "ccorb.h" 1784#include "ccsdsym.h" 1785#include "r12int.h" 1786 integer isymkint, index, isymakl, isyma, isyml, 1787 & isymk, isymak, isymbl, isymkl, idxbl, idxkl, idxakl, 1788 & idxabkl, idxak, isymb 1789 integer isym,isym1,isym2,icoun1,icoun2,iak(8,8),nak(8),iakbl(8,8), 1790 & ipk(8,8),npk(8),ipkql(8,8),nr12orb(8),nkl(8),ikl(8,8), 1791 & icoun3,nakbl(8) 1792#if defined (SYS_CRAY) 1793 real xkint(*), xk(*), ddot 1794#else 1795 double precision xkint(*), xk(*), ddot 1796#endif 1797 logical lopt,locdbg 1798 parameter (locdbg = .false.) 1799 1800 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 1801 call qenter('r12sortk') 1802c 1803 do isym = 1, nsym 1804 icoun1 = 0 1805 icoun2 = 0 1806 icoun3 = 0 1807 do isym2 = 1, nsym 1808 isym1 = muld2h(isym,isym2) 1809 iak(isym1,isym2) = icoun1 1810 ipk(isym1,isym2) = icoun2 1811 ikl(isym1,isym2) = icoun3 1812 icoun1 = icoun1 + nvir(isym1)*nr12orb(isym2) 1813 icoun2 = icoun2 + norb1(isym1)*nr12orb(isym2) 1814 icoun3 = icoun3 + nr12orb(isym1)*nr12orb(isym2) 1815 end do 1816 nak(isym) = icoun1 1817 npk(isym) = icoun2 1818 nkl(isym) = icoun3 1819 end do 1820 do isym = 1, nsym 1821 icoun1 = 0 1822 icoun2 = 0 1823 do isym2 = 1, nsym 1824 isym1 = muld2h(isym,isym2) 1825 if (isym2.gt.isym1) then 1826 iakbl(isym1,isym2) = icoun1 1827 iakbl(isym2,isym1) = icoun1 1828 icoun1 = icoun1 + nak(isym1)*nak(isym2) 1829 ipkql(isym1,isym2) = icoun2 1830 ipkql(isym2,isym1) = icoun2 1831 icoun2 = icoun2 + npk(isym1)*npk(isym2) 1832 else if (isym2.eq.isym1) then 1833 iakbl(isym1,isym2) = icoun1 1834 icoun1 = icoun1 + nak(isym1)*(nak(isym2)+1)/2 1835 ipkql(isym1,isym2) = icoun2 1836 icoun2 = icoun2 + npk(isym1)*(npk(isym2)+1)/2 1837 end if 1838 end do 1839 nakbl(isym) = icoun1 1840 end do 1841c 1842 if (locdbg) then 1843 write(lupri,*) 'Entered CC_R12SORTK' 1844 write(lupri,*)'norm^2(xkint): ', 1845 & ddot(nakbl(isymkint),xkint,1,xkint,1) 1846C if (nsym.eq.1) 1847C & call outpak(xkint,nak(1),1,lupri) 1848 end if 1849c 1850 isymakl = muld2h(isymb,isymkint) 1851 isymkl = muld2h(isyma,isymakl) 1852 do isymk = 1, nsym 1853 isyml = muld2h(isymkl,isymk) 1854 isymak = muld2h(isyma,isymk) 1855 isymbl = muld2h(isymb,isyml) 1856 if (.not.lopt) then 1857 do l = 1, nr12orb(isyml) 1858 idxbl = iak(isymb,isyml)+nvir(isymb)*(l-1)+b 1859 do k = 1, nr12orb(isymk) 1860 idxkl = ikl(isymk,isyml)+nr12orb(isymk)*(l-1)+k 1861 if (isymak.eq.isymbl) then 1862 do a = 1, nvir(isyma) 1863 idxakl = nvir(isyma)*(idxkl-1)+a 1864 idxak = iak(isyma,isymk)+nvir(isyma)*(k-1)+a 1865 idxabkl = iakbl(isymak,isymbl)+index(idxak,idxbl) 1866 xk(idxakl) = xkint(idxabkl) 1867 end do 1868 else if (isymak.lt.isymbl) then 1869 idxakl = nvir(isyma)*(idxkl-1)+1 1870 idxak = iak(isyma,isymk)+nvir(isyma)*(k-1)+1 1871 idxabkl = iakbl(isymak,isymbl)+ 1872 & nak(isymak)*(idxbl-1)+idxak 1873 call dcopy(nvir(isyma),xkint(idxabkl),1, 1874 & xk(idxakl),1) 1875 else if (isymak.gt.isymbl) then 1876 do a = 1, nvir(isyma) 1877 idxakl = nvir(isyma)*(idxkl-1)+a 1878 idxak = iak(isyma,isymk)+nvir(isyma)*(k-1)+a 1879 idxabkl = iakbl(isymbl,isymak)+ 1880 & nak(isymbl)*(idxak-1)+idxbl 1881 xk(idxakl) = xkint(idxabkl) 1882 end do 1883 end if 1884 end do 1885 end do 1886 else 1887 do l = 1, nr12orb(isyml) 1888 idxbl=ipk(isymb,isyml)+norb1(isymb)*(l-1)+nrhfsa(isymb)+b 1889 do k = 1, nr12orb(isymk) 1890 idxkl = ikl(isymk,isyml)+nr12orb(isymk)*(l-1)+k 1891 if (isymak.eq.isymbl) then 1892 do a = 1, nvir(isyma) 1893 idxakl = nvir(isyma)*(idxkl-1)+a 1894 idxak = ipk(isyma,isymk) + 1895 & norb1(isyma)*(k-1)+nrhfsa(isyma)+a 1896 idxabkl = ipkql(isymak,isymbl)+index(idxak,idxbl) 1897 xk(idxakl) = xkint(idxabkl) 1898 end do 1899 else if (isymak.lt.isymbl) then 1900 idxakl = nvir(isyma)*(idxkl-1)+1 1901 idxak = ipk(isyma,isymk) + 1902 & norb1(isyma)*(k-1)+nrhfsa(isyma)+1 1903 idxabkl = ipkql(isymak,isymbl)+ 1904 & npk(isymak)*(idxbl-1)+idxak 1905 call dcopy(nvir(isyma),xkint(idxabkl),1, 1906 & xk(idxakl),1) 1907 else if (isymak.gt.isymbl) then 1908 do a = 1, nvir(isyma) 1909 idxakl = nvir(isyma)*(idxkl-1)+a 1910 idxak = ipk(isyma,isymk) + 1911 & norb1(isyma)*(k-1)+nrhfsa(isyma)+a 1912 idxabkl = ipkql(isymbl,isymak)+ 1913 & npk(isymbl)*(idxak-1)+idxbl 1914 xk(idxakl) = xkint(idxabkl) 1915 end do 1916 end if 1917 end do 1918 end do 1919 end if 1920 end do 1921c 1922 if (locdbg) then 1923 write(lupri,*)'norm^2 sorted integrals: ', 1924 & ddot(nvir(isyma)*nkl(isymkl),xk,1,xk,1) 1925C if (nsym.eq.1) then 1926C write(lupri,'(a,3i5)') 1927C & 'integrals for fixed isymb,b,isyma:',isymb,b,isyma 1928C call output(xk,1,nvir(isyma),1,nkl(1), 1929C & nvir(isyma),nkl(1),1,lupri) 1930C end if 1931 end if 1932c 1933 call qexit('r12sortk') 1934 end 1935*======================================================================* 1936 subroutine ccrhs_eppp(omega2,work,lwork,aproxr12,lres, 1937 & lufc12,fc12am,ifile,isymtr) 1938*----------------------------------------------------------------------* 1939c Purpose: add E''' term to vector function in omega2 for the 1940c ansaetze 2 and 3 1941c 1942c omega^ab_ij <-- omega^ab_ij 1943c + sum_kl (K^ab_kl - t^ab_kl) * c^kl_ij 1944c 1945c the contrib. of K is skipped for the approximations A, A' 1946c 1947c for the approximations ??? is in ansatz 2 added: 1948c 1949c omega^ab_ij <-- omega^ab_ij 1950c + sum_kl (e_k+e_l-e_i-e_j) r^ab_kl * c^kl_ij 1951c 1952c while in ansatz 3 the following term is added: 1953c 1954c omega^ab_ij <-- omega^ab_ij 1955c + sum_kl (e_k+e_l-e_a-e_b) r^ab_ij * c^kl_ij 1956c 1957c H. Fliegl, C. Haettig, spring 2004 1958c adapted for ansatz 3, Christof Haettig, spring 2005 1959*----------------------------------------------------------------------* 1960 implicit none 1961#include "priunit.h" 1962#include "ccorb.h" 1963#include "ccsdsym.h" 1964#include "dummy.h" 1965#include "r12int.h" 1966#include "ccr12int.h" 1967 1968 logical locdbg, lres, needkrint 1969 parameter (locdbg = .false.) 1970 character*(*) fc12am 1971 character*(*) aproxr12 1972 character*10 model 1973 integer lwork, lwrk1, kend1, ktint, kcamp, krint 1974 integer lunit, idum, isymkint, iopt, kend2, lwrk2 1975 integer ifile, lufc12, isymtr 1976 1977#if defined (SYS_CRAY) 1978 real omega2(*), work(*), factor, ddot 1979#else 1980 double precision omega2(*), work(*), factor, ddot 1981#endif 1982 1983 call qenter('rhs_eppp') 1984 1985 if (locdbg) then 1986 write(lupri,*) 'entered ccrhs_eppp... isymtr =',isymtr 1987 end if 1988 1989 ! is a second integral array of length nt2r12(1) needed? 1990 needkrint = (aproxr12(1:1).eq."B") 1991 if (ianr12.eq.3 .and. 1992 & (aproxr12(1:2).eq."A'".or.aproxr12(1:1).eq."B") ) then 1993 needkrint=.true. 1994 end if 1995 1996 ktint = 1 1997 kcamp = ktint + nt2r12(1) 1998 kend1 = kcamp + ntr12am(isymtr) 1999 lwrk1 = lwork - kend1 2000 2001 kend2 = kend1 2002 if (needkrint) then 2003 krint = kend1 2004 kend2 = krint + nt2r12(1) 2005 end if 2006 lwrk2 = lwork - kend2 2007 2008 if (lwrk1.lt.0 .or. lwrk2.lt.0) then 2009 call quit('Insufficient work space in ccrhs_eppp') 2010 end if 2011 2012c --------------------------------------------- 2013c read t^ab_kl integrals from file 2014c --------------------------------------------- 2015 lunit = -1 2016 call gpopen(lunit,fkr12,'unknown',' ','unformatted',idum,.false.) 2017 read(lunit)(work(ktint+i-1),i=1,nt2r12(1)) 2018 call gpclose(lunit,'KEEP') 2019 2020c ---------------------------------------------------- 2021c for approximation B substract K^ab_kl from t^ab_kl 2022c ---------------------------------------------------- 2023 if (aproxr12(1:1).eq."B") then 2024 ! read K^ab_kl from file 2025 lunit = -1 2026 call gpopen(lunit,ftr12,'unknown',' ','unformatted', 2027 & idum,.false.) 2028 read(lunit)(work(krint+i-1),i=1,nt2r12(1)) 2029 call gpclose(lunit,'KEEP') 2030 2031 ! substract (t - K) => work(ktint) 2032 call daxpy(nt2r12(1),-1.0d0,work(krint),1,work(ktint),1) 2033 end if 2034 2035c ------------------------------------------------------------ 2036c in ansatz 3 add for approximation A' and B the term 2037c (eps_k + eps_l - eps_a - eps_b) * r^ab_kl 2038c ------------------------------------------------------------ 2039 if (ianr12.eq.3 .and. 2040 & (aproxr12(1:2).eq."A'".or.aproxr12(1:1).eq."B") ) then 2041 lunit = -1 2042 call gpopen(lunit,fr12r12,'unknown',' ','unformatted', 2043 & idum,.false.) 2044 read(lunit)(work(krint+i-1),i=1,nt2r12(1)) 2045 call gpclose(lunit,'KEEP') 2046 2047 call cc_t2xe(work(krint),work(kend2),lwrk2) 2048 2049 call daxpy(nt2r12(1),+1.0d0,work(krint),1,work(ktint),1) 2050 end if 2051 2052c ------------------------------------------------------ 2053c get R12 amplitudes c^mn_kl stored as triangular matrix 2054c ------------------------------------------------------ 2055 if (lres) then 2056 call cc_rvec(lufc12,fc12am,ntr12am(isymtr),ntr12am(isymtr), 2057 & ifile,work(kcamp)) 2058 call cclr_diasclr12(work(kcamp),ketscl,isymtr) 2059 else 2060c read r12 amplitudes from file 2061 iopt = 32 2062 call cc_rdrsp('R0 ',0,1,iopt,model,dummy,work(kcamp)) 2063 end if 2064 2065 if (locdbg) then 2066 write(lupri,*) 'in ccrhs_eppp norm^2 r12 amplitudes', 2067 & ddot(ntr12am(isymtr),work(kcamp),1,work(kcamp),1) 2068 write(lupri,*)'Omega_aibj before contribution of ccrhs_eppp:', 2069 & ddot(nt2am(isymtr),omega2,1,omega2,1) 2070 call flshfo(lupri) 2071 end if 2072 2073c ------------------------------------------------------ 2074c calculate Omega^ab_ij += - sum_kl k^ab_kl * c^kl_ij 2075c ------------------------------------------------------ 2076 isymkint = 1 2077 factor = 1.0d0 2078 call cc_r12mi22(omega2,work(ktint),work(kcamp),isymkint,isymtr, 2079 & factor,work(kend1),lwrk1) 2080 2081 if (locdbg) then 2082 write(lupri,*)'norm^2 omega_aibj after mi22:', 2083 & ddot(nt2am(isymtr),omega2,1,omega2,1) 2084 call flshfo(lupri) 2085 end if 2086 2087c ------------------------------------------------------------ 2088c in ansatz 2 add for approximation A' and B the term 2089c omega^ab_ij += sum_kl (e_k+e_l-e_i-e_j) r^ab_kl * c^kl_ij 2090c ------------------------------------------------------------ 2091 if (ianr12.eq.2 .and. 2092 & (aproxr12(1:2).eq."A'".or.aproxr12(1:1).eq."B") ) then 2093 2094 if (locdbg) then 2095 write(lupri,*) 'R12 amplitudes before cc_r12mkediff2:' 2096 call cc_prpr12(work(kcamp),isymtr,1,.false.) 2097 end if 2098 ! multiply R12 amplitudes with orbital energy differences: 2099 ! c^kl_ij <-- (e_k + el - e_i - e_j) * c^kl_ij 2100 call cc_r12mkediff2(work(kcamp),isymtr,work(kend1),lwrk1) 2101 if (locdbg) then 2102 write(lupri,*) 'R12 amplitudes after cc_r12mkediff2:' 2103 call cc_prpr12(work(kcamp),isymtr,1,.false.) 2104 end if 2105 2106 lunit = -1 2107 call gpopen(lunit,fr12r12,'unknown',' ','unformatted', 2108 & idum,.false.) 2109 read(lunit)(work(ktint+i-1),i=1,nt2r12(1)) 2110 call gpclose(lunit,'KEEP') 2111 2112 isymkint = 1 2113 factor = 1.0d0 2114 call cc_r12mi22(omega2,work(ktint),work(kcamp),isymkint,isymtr, 2115 & factor,work(kend1),lwrk1) 2116 end if 2117 2118c ---------------------------------------------- 2119c that's it; print some debug output and return: 2120c ---------------------------------------------- 2121 if (locdbg) then 2122 write(lupri,*)'Omega_aibj after ccrhs_eppp contribution:', 2123 & ddot(nt2am(isymtr),omega2,1,omega2,1) 2124 end if 2125 2126 call qexit('rhs_eppp') 2127 end 2128*=====================================================================* 2129 subroutine cc_r12mi22(omega2,xkint,c2am,isymkint,isymc, 2130 & factor,work,lwork) 2131c-------------------------------------------------------------------- 2132c purpose: make Omega_ab_ij 2133c 2134c H. Fliegl, C. Haettig, spring 2004 2135c-------------------------------------------------------------------- 2136 implicit none 2137#include "priunit.h" 2138#include "ccorb.h" 2139#include "ccsdsym.h" 2140 logical locdbg 2141 parameter (locdbg = .false.) 2142 integer lwork, isymkint, isymc, isymb, isyma, isymab, isymkl, 2143 & isymij, isymj, isymi, kxint, kc2am, kend1, isymint 2144 integer ntoti, ntota, kom, isymai, isymbj, maxai, idxai, idxbj, 2145 & naibj, koffai, lwrk1 2146#if defined (SYS_CRAY) 2147 real work(*), xkint(*),c2am(*), omega2(*),factor 2148 real one, zero,xk 2149#else 2150 double precision work(*), xkint(*),c2am(*), omega2(*),factor 2151 double precision one, zero,ddot,x1,x2,x3,xk 2152#endif 2153 parameter (zero = 0.0d0, one = 1.0d0) 2154 2155 integer index 2156 index(i,j) = max(i,j)*(max(i,j) - 3)/2 + i + j 2157 2158 call qenter('r12mi22') 2159 isymint = muld2h(isymkint,isymc) 2160 2161 if (locdbg) then 2162 write(lupri,*) 'entered cc_r12mi22' 2163 x1 = 0.0d0 2164 xk = 0.0d0 2165 x2 = 0.0d0 2166 x3 = 0.0d0 2167 end if 2168 2169 do isymb = 1, nsym 2170 do isyma = 1, nsym 2171 isymab = muld2h(isyma,isymb) 2172 isymkl = muld2h(isymkint,isymab) 2173 isymij = muld2h(isymkl,isymc) 2174 do isymj = 1, nsym 2175 isymi = muld2h(isymij,isymj) 2176 isymai = muld2h(isyma,isymi) 2177 isymbj = muld2h(isymb,isymj) 2178 2179 kxint = 1 2180 kc2am = kxint + nvir(isyma)*nmatkl(isymkl) 2181 kom = kc2am + nrhf(isymi)*nmatkl(isymkl) 2182 kend1 = kom + nvir(isyma)*nrhf(isymi) 2183 lwrk1 = lwork - kend1 2184c 2185 if (lwrk1.lt.0) then 2186 call quit('Insufficient memory in cc_r12mi22') 2187 end if 2188c 2189 do b = 1, nvir(isymb) 2190C call cc_r12sortk(work(kxint),xkint,b,isymkint,isymb, 2191C & isyma,nrhfb,.false.) 2192 call cc_r12sort(xkint,isymkint,work(kxint),b,isymb, 2193 & 1,nvir(isyma),isyma,nvir,nrhfb) 2194 xk = xk + ddot(nvir(isyma)*nmatkl(isymkl), 2195 & work(kxint),1,work(kxint),1) 2196 do j = 1, nrhf(isymj) 2197 call cc_r12sortc(work(kc2am),c2am,j,isymc,isymj,isymi) 2198 2199 ntoti = max(nrhf(isymi),1) 2200 ntota = max(nvir(isyma),1) 2201 call dgemm('N','T',nvir(isyma),nrhf(isymi), 2202 & nmatkl(isymkl),one,work(kxint),ntota, 2203 & work(kc2am),ntoti,zero,work(kom),ntota) 2204c 2205 if (locdbg) then 2206 x1 = x1 + ddot(nvir(isyma)*nmatkl(isymkl), 2207 & work(kxint),1,work(kxint),1) 2208 x2 = x2 + ddot(nrhf(isymi)*nmatkl(isymkl), 2209 & work(kc2am),1,work(kc2am),1) 2210 x3 = x3 + ddot(nvir(isyma)*nrhf(isymi), 2211 & work(kom),1,work(kom),1) 2212 end if 2213 2214 ! subtract the result from Omega(ai,bj) which is 2215 ! stored as a triangular matrix 2216 if (isymbj.eq.isymai) then 2217 idxbj = it1am(isymb,isymj) + nvir(isymb)*(j-1) + b 2218 koffai = it1am(isyma,isymi) 2219 maxai = min(koffai+nvir(isyma)*nrhf(isymi),idxbj) 2220 2221 do idxai = koffai+1, maxai 2222 naibj = it2am(isymai,isymbj) + index(idxai,idxbj) 2223 omega2(naibj)= omega2(naibj)+ 2224 & factor*work(kom-1+idxai-koffai) 2225 end do 2226c 2227 else if (isymbj.gt.isymai) then 2228c ! omega not totally symmetric ... only needed for 2229c ! response 2230 idxbj = it1am(isymb,isymj) + nvir(isymb)*(j-1) + b 2231 koffai = it1am(isyma,isymi) 2232 maxai = koffai+nvir(isyma)*nrhf(isymi) 2233 2234 do idxai = koffai+1, maxai 2235 naibj = it2am(isymai,isymbj) + 2236 & nt1am(isymai)*(idxbj-1) + idxai 2237 omega2(naibj)= omega2(naibj)+ 2238 & factor*work(kom-1+idxai-koffai) 2239 end do 2240 end if 2241 end do 2242 end do 2243c 2244 end do 2245 end do 2246 end do 2247 2248 if (locdbg) then 2249 write(lupri,*)'xk =', xk 2250 write(lupri,*)'x1 =', x1 2251 write(lupri,*)'x2 =', x2 2252 write(lupri,*)'x3 =', x3 2253 end if 2254 2255 call qexit('r12mi22') 2256 end 2257*=====================================================================* 2258 subroutine cc_r12sortc(xc,c2am,j,isymc,isymj,isymi) 2259c-------------------------------------------------------------------- 2260c purpose: sort c^kl_ij as c^j_i,kl 2261c 2262c H. Fliegl, C. Haettig, spring 2004 2263c-------------------------------------------------------------------- 2264 implicit none 2265#include "priunit.h" 2266#include "ccorb.h" 2267#include "ccsdsym.h" 2268 logical locdbg 2269 parameter ( locdbg = .false. ) 2270 integer isymc, index, isymj, isymi, isymikl, isymkl, isymk, 2271 & isyml, isymki, isymlj, idxlj, idxkl, idxikl, idxki, 2272 & idxkilj 2273#if defined (SYS_CRAY) 2274 real c2am(*), xc(*) 2275 real ddot 2276#else 2277 double precision c2am(*), xc(*) 2278 double precision ddot 2279#endif 2280 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 2281 call qenter('r12sortc') 2282 2283 if (locdbg) then 2284 write(lupri,*)'entered cc_r12sortc..' 2285 write(lupri,*)'norm^2(c2am):',ddot(ntr12am(isymc),c2am,1,c2am,1) 2286 call cc_prpr12(c2am,isymc,1,.false.) 2287 end if 2288 2289 isymikl = muld2h(isymj,isymc) 2290 isymkl = muld2h(isymi,isymikl) 2291 do isymk = 1, nsym 2292 isyml = muld2h(isymkl,isymk) 2293 isymki = muld2h(isymi,isymk) 2294 isymlj = muld2h(isymj,isyml) 2295 2296 do l = 1, nrhfb(isyml) 2297 idxlj = imatki(isyml,isymj) + nrhfb(isyml)*(j-1)+l 2298 do k = 1, nrhfb(isymk) 2299 idxkl = imatkl(isymk,isyml) + nrhfb(isymk)*(l-1)+k 2300 if (isymki.eq.isymlj) then 2301 do i =1, nrhf(isymi) 2302 idxikl = nrhf(isymi)*(idxkl-1)+i 2303 idxki = imatki(isymk,isymi)+nrhfb(isymk)*(i-1)+k 2304 idxkilj = itr12am(isymki,isymlj)+index(idxki,idxlj) 2305 xc(idxikl) = c2am(idxkilj) 2306 end do 2307 else if (isymki.lt.isymlj) then 2308 idxikl = nrhf(isymi)*(idxkl-1)+1 2309 idxki = imatki(isymk,isymi)+k 2310 idxkilj = itr12am(isymki,isymlj)+ 2311 & nmatki(isymki)*(idxlj-1)+idxki 2312 call dcopy(nrhf(isymi),c2am(idxkilj),nrhfb(isymk), 2313 & xc(idxikl),1) 2314 else if (isymki.gt.isymlj) then 2315 do i =1, nrhf(isymi) 2316 idxikl = nrhf(isymi)*(idxkl-1)+i 2317 idxki = imatki(isymk,isymi)+nrhfb(isymk)*(i-1)+k 2318 idxkilj = itr12am(isymlj,isymki)+ 2319 & nmatki(isymlj)*(idxki-1)+idxlj 2320 xc(idxikl) = c2am(idxkilj) 2321 end do 2322 end if 2323 end do 2324 end do 2325 end do 2326 2327 if (locdbg) then 2328 write(lupri,'(a,3i5)') 2329 & 'R12 amplitudes for fixed isymj,j,isymi:',isymj,j,isymi 2330 call output(xc,1,nrhf(isymi),1,nmatkl(isymkl), 2331 & nrhf(isymi),nmatkl(isymkl),1,lupri) 2332 end if 2333 2334 call qexit('r12sortc') 2335 end 2336*=====================================================================* 2337 subroutine cc_r12batf2(work,lwork) 2338c-------------------------------------------------------------------- 2339c purpose: prepare R^d_akl for ansatz 2 to obtain (V^kl_aj)^+ 2340c 2341c H. Fliegl, C. Haettig, sping 2004 2342c 2343c C. Neiss, winter 2006: adapted for CABS 2344c-------------------------------------------------------------------- 2345 implicit none 2346#include "priunit.h" 2347#include "ccorb.h" 2348#include "ccsdsym.h" 2349#include "r12int.h" 2350#include "ccr12int.h" 2351 2352#if defined (SYS_CRAY) 2353 real zero,one, work(*),ddot 2354#else 2355 double precision zero,one, work(*),ddot 2356#endif 2357 parameter (zero = 0.0D0, one = 1.0D0) 2358 2359 logical locdbg,lauxd 2360 parameter (locdbg = .false.) 2361 2362 integer nr1orb(8),nr1bas(8),nr1xbas(8),nr2bas,n2bst1(8) 2363 integer ir1xbas(8,8),nxx2(8),lwork,isymmu,kscr1 2364 integer ir1orb(8,8),ir1bas(8,8),ir2bas(8,8),iaodis1(8,8) 2365 integer ir2xbas(8,8),irgkl(8,8),nrgkl(8) 2366 integer nnbastx,n2bst2(8),isymd,isymk,isym, 2367 & kend0,lusifc,nsymx,norbsx(8),nbastx,nlamdsx, 2368 & nrhfsx(8),kcmox,lwrk0,ksfull,iaodis2(8,8) 2369 integer idummy,norbtsx,isympq,isymp,isymq,ksjbp,ksmupbp 2370 integer kend1,lwrk1,isymbp,isymdp,isymakl,krdpakl,krtot,kend2, 2371 & lwrk2,ib,id,idxbp,idelp,idxbpdp,koff1,idxd, 2372 & nalphaj(8),ialphaj(8,8),nbasx(8),basoff(8) 2373 integer nr1xorb(8),ir1xorb(8,8),nrxgkl(8),irxgkl(8,8) 2374 2375 call qenter('batf2') 2376 if(locdbg) write(lupri,*) 'entered batf2' 2377 2378 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 2379 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas,ir2bas, 2380 & ir2xbas,irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 2381c 2382 if (r12cbs) then 2383 do isym = 1, nsym 2384 nbasx(isym) = mbas1(isym)+mbas2(isym) 2385 basoff(isym) = 0 2386 end do 2387 else 2388 call icopy(8,mbas2,1,nbasx,1) 2389 call icopy(8,mbas1,1,basoff,1) 2390 end if 2391c 2392 do isympq = 1, nsym 2393 nxx2(isympq) = 0 2394 do isymq = 1, nsym 2395 isymp = muld2h(isympq,isymq) 2396 iaodis2(isymp,isymq) = nxx2(isympq) 2397 nxx2(isympq) = nxx2(isympq) + nbasx(isymp)*nbasx(isymq) 2398 end do 2399 end do 2400c 2401 kend0 = 1 2402 2403c---------------------------------------------------------------------- 2404c get MO coefficients defined for combined orbital + aux. basis: 2405c---------------------------------------------------------------------- 2406 lusifc = -1 2407 call gpopen(lusifc,'SIRIFC','OLD',' ','UNFORMATTED', 2408 & idummy,.false.) 2409 2410c read dimensions for CMO coefficients for full basis (nlamdsx) 2411 rewind(lusifc) 2412 call mollab('FULLBAS ',lusifc,lupri) 2413 read(lusifc) nsymx,norbtsx,nbastx,nlamdsx,(nrhfsx(i),i=1,nsym), 2414 & (norbsx(i),i=1,nsym) 2415 2416 2417c allocate work space for CMO coefficients for full basis 2418 kcmox = kend0 2419 kend0 = kcmox + nlamdsx 2420 lwrk0 = lwork - kend0 2421 if (lwrk0 .lt.0) then 2422 call quit('Insufficient work space in batf2') 2423 end if 2424 2425c read CMO coefficients and close file 2426 read(lusifc) 2427 read(lusifc) (work(kcmox+i-1),i=1,nlamdsx) 2428 call gpclose(lusifc,'KEEP') 2429 2430c test if CMO is ok 2431 if (locdbg) then 2432 write(lupri,*)'CMO out batf2' 2433 kscr1 = 1 2434 do isymp = 1, nsym 2435 isymmu = isymp 2436 write(lupri,*) 'symmetry block,nbas,norbs:',isymp, 2437 & nbas(isymmu), norbs(isymp) 2438 call output(work(kcmox+kscr1-1),1,nbas(isymmu), 2439 & 1,norbs(isymp),nbas(isymmu),norbs(isymp), 2440 & 1,lupri) 2441 2442 kscr1 = kscr1 + nbas(isymmu)*norbs(isymp) 2443 end do 2444 end if 2445 2446c get S_mu'b' 2447 ksmupbp = kend0 2448 kend1 = ksmupbp + nxx2(1) 2449 lwrk1 = lwork - kend1 2450 if (lwrk1 .lt.0) then 2451 call quit('Insufficient work space in batf2') 2452 end if 2453 call cc_r12mksmupbp(work(kcmox),work(ksmupbp),iaodis2, 2454 & work(kend1),lwrk1) 2455 2456 if (locdbg) then 2457 write(lupri,*) "norm^2(S_mu'b'):",ddot(nxx2(1),work(ksmupbp),1, 2458 & work(ksmupbp),1) 2459 do isymp = 1, nsym 2460 write(lupri,*) 'symmetry block:',isymp 2461 call output(work(ksmupbp+iaodis2(isymp,isymp)), 2462 & 1,nbasx(isymp),1,nbasx(isymp), 2463 & nbasx(isymp),nbasx(isymp),1,lupri) 2464 end do 2465 end if 2466 2467c----------------------------------------------------------------------- 2468c make R^d'_a,kl 2469c----------------------------------------------------------------------- 2470 do isymdp = 1, nsym 2471 isymbp = isymdp 2472 isymakl = isymdp 2473 2474 krdpakl = kend1 2475 krtot = krdpakl + nrgkl(isymakl) 2476 kend2 = krtot + nrgkl(isymakl) 2477 lwrk2 = lwork - kend2 2478 if(lwrk2.lt.0) then 2479 call quit('Insufficient works space in batf2 ') 2480 end if 2481 2482 do id = 1, mbas1(isymdp)+mbas2(isymdp) 2483 lauxd = .false. 2484 idxd = id + ibas(isymdp) 2485 if (id .gt. basoff(isymdp)) then 2486 !make backtransformation... 2487 call dzero(work(krdpakl),nrgkl(isymakl)) 2488 idelp = id - basoff(isymbp) 2489 do ib = 1, nbasx(isymbp) 2490 lauxd = .false. 2491 idxbp = ib + ibas(isymbp) + basoff(isymbp) 2492 call cc_r12getrint(work(krtot),idxbp,isymbp,nr1bas,ir1bas, 2493 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 2494 & nrhfb,nmatkl,imatkl, 2495 & idummy,lauxd,.false.,frhtf,work(kend2),lwrk2) 2496 idxbpdp = iaodis2(isymbp,isymdp) + 2497 & nbasx(isymbp)*(idelp-1) + ib 2498 call daxpy(nrgkl(isymakl),work(ksmupbp-1+idxbpdp), 2499 & work(krtot),1,work(krdpakl),1) 2500 end do 2501 else 2502 !only copy... 2503 call cc_r12getrint(work(krdpakl),idxd,isymdp,nr1bas,ir1bas, 2504 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 2505 & nrhfb,nmatkl,imatkl, 2506 & idummy,lauxd,.false.,frhtf,work(kend2),lwrk2) 2507 end if 2508 call cc_r12putrint(work(krdpakl),idxd,isymdp,nr1bas,ir1bas, 2509 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 2510 & nrhfb,nmatkl,imatkl, 2511 & idummy,lauxd,fnback2,work(kend2),lwrk2) 2512 end do 2513 end do 2514c 2515 call qexit('batf2') 2516 end 2517*=====================================================================* 2518 subroutine cc_r12mkvaj2(vajkl,isymv,cmo,isymc,cmos,isymcs, 2519 & work,lwork) 2520c--------------------------------------------------------------------- 2521c purpose: make \sum_mn r_mn^kl*g_aj^mn 2522c 2523c H. Fliegl, C. Haettig, spring 2004 2524c--------------------------------------------------------------------- 2525 implicit none 2526#include "priunit.h" 2527#include "ccorb.h" 2528#include "ccsdsym.h" 2529#include "r12int.h" 2530#include "ccr12int.h" 2531#include "dummy.h" 2532 logical locdbg, lauxd 2533 parameter (locdbg = .false.) 2534 integer isymv, isymc, lwork, lwrk1, lwrk2, kend1, kend2 2535 integer krmnkl, kgmakl, krhtf, koffv, koffg, koffr, idxkl, idxd, 2536 & ntotaj, ntotmn, idelta, isymd, isymnm, isymaj, isymkl 2537 integer nr1orb(8), nr1bas(8), nr1xbas(8), nr2bas, n2bst1(8) 2538 integer ir1xbas(8,8), ibasx(8),kmkl,iadr1,iadr2,isymdm,isymm,ndim 2539 integer iglmrhs(8,8),iglmvis(8,8),nglmds(8),icmo(8,8),ncmo(8), 2540 & imaklm(8,8),nmaklm(8),imatijm(8,8),nmatijm(8), 2541 & igamsm(8,8),ngamsm(8),ivajkls(8,8),nvajkls(8) 2542 integer ir1orb(8,8), ir1bas(8,8), ir2bas(8,8), iaodis1(8,8) 2543 integer ir2xbas(8,8), irgkl(8,8), nrgkl(8), isyma,isymj 2544 integer irgijs(8,8),nrgijs(8),ir1basm(8,8),nr1basm(8), 2545 & ir2basm(8,8),nr2basm,ir1xbasm(8,8),nr1xbasm(8), 2546 & ir2xbasm(8,8),isym,icount1,isym2,isym1,isymg,koffc, 2547 & koff1,ntotg,ntotm,imatf(8,8),nmatf(8) 2548 integer ioffg(8,8),isymn,isymmn,idxmn,kgtf,idxdn,idxlk,isymk, 2549 & isyml,isymcs 2550 integer nr1xorb(8),ir1xorb(8,8),nrxgkl(8),irxgkl(8,8), 2551 & nalphaj(8),ialphaj(8,8),nmaijm(8),imaijm(8,8) 2552 double precision one, zero 2553 parameter (one = 1.0d0, zero = 0.0d0) 2554#if defined (SYS_CRAY) 2555 real work(*), vajkl(*), cmo(*),ddot,ddummy,cmos(*),factor 2556#else 2557 double precision work(*),vajkl(*),cmo(*),ddot,ddummy,cmos(*), 2558 & factor 2559#endif 2560 2561 call qenter('r12mkvaj2') 2562 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 2563 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas, 2564 & ir2bas,ir2xbas,irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 2565 call cc_r12offs23(iglmrhs,iglmvis,nglmds,icmo,ncmo, 2566 & imaijm,nmaijm,imaklm,nmaklm, 2567 & imatijm,nmatijm,igamsm,ngamsm,irgijs,nrgijs, 2568 & ir1basm,nr1basm,ir2basm,nr2basm,ir1xbasm, 2569 & nr1xbasm,ir2xbasm,imatf,nmatf) 2570c 2571 ibasx(1) = 0 2572 do isym = 2, nsym 2573 ibasx(isym) = ibasx(isym-1) + mbas2(isym-1) 2574 end do 2575c 2576 do isym = 1, nsym 2577 icount1 = 0 2578 do isym2 = 1, nsym 2579 isym1 = muld2h(isym,isym2) 2580 ivajkls(isym1,isym2) = icount1 2581 icount1 = icount1 + nalphaj(isym1)*nmatijm(isym2) 2582 end do 2583 nvajkls(isym) = icount1 2584 end do 2585c 2586 if (locdbg) then 2587 do isym = 1, nsym 2588 icount1 = 0 2589 do isym2 = 1, nsym 2590 isym1 = muld2h(isym,isym2) 2591 ioffg(isym1,isym2) = icount1 2592 icount1 = icount1 + nmatij(isym1)*nmatijm(isym2) 2593 end do 2594 end do 2595 end if 2596 2597 krmnkl = 1 2598 kgmakl = krmnkl + ngamsm(1) 2599 kend1 = kgmakl + nvajkls(1) 2600 lwrk1 = lwork - kend1 2601 if (lwrk1.lt.0) then 2602 call quit('Insufficient work space in mkvaj2') 2603 end if 2604 2605 call dzero(work(krmnkl),ngamsm(1)) 2606 call dzero(work(kgmakl),nvajkls(1)) 2607c 2608 do isymd = 1, nsym 2609 krhtf = kend1 2610 kmkl = krhtf + max(nrgkl(isymd),nrgijs(isymd)) 2611 kend2 = kmkl + nmatf(isymd) 2612 lwrk2 = lwork - kend2 2613 if (lwrk2.lt.0) then 2614 call quit('Insufficient work space in mkvaj2') 2615 end if 2616 do idelta = 1, mbas1(isymd) 2617 lauxd = .false. 2618 idxd = idelta+ibas(isymd) 2619 call cc_r12getrint(work(krhtf),idxd,isymd,nr1bas,ir1bas, 2620 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 2621 & nrhfb,nmatkl,imatkl, 2622 & ibasx,lauxd,.false.,frhtf,work(kend2),lwrk2) 2623 2624c get rMNkl 2625 call cc_r12generaltf(work(krhtf),work(krmnkl),idelta,isymd, 2626 & cmos,isymcs,cmos,isymcs,.false.,iglmrhs, 2627 & iglmrhs,nmatkl,nrhfsa,nbas,irgkl, 2628 & imatijm,nmatijm,imaklm,igamsm, 2629 & work(kend2),lwrk2) 2630c 2631 lauxd = .false. 2632 idxd = idelta+ibas(isymd) 2633c get g^delta_gamma,KL 2634 call cc_r12getrints(work(krhtf),idxd,isymd,nr1basm,ir1basm, 2635 & nr2basm,ir2basm,nrgijs,irgijs,ir1xbasm,ir2xbasm, 2636 & ibasx,imatijm,nmatijm,lauxd,fccgmnab,work(kend2),lwrk2) 2637c 2638 if (locdbg) then 2639 write(lupri,*) 'norm^2(Gmbnd) after read:',idxd, 2640 & ddot(nrgijs(isymd),work(krhtf),1,work(krhtf),1) 2641 end if 2642c 2643c transform g^delta_gamma,KL to g^delta_m,KL 2644 do isymg = 1, nsym 2645 isymkl = muld2h(isymg,isymd) 2646 isymm = muld2h(isymc,isymg) 2647 2648c koffc = iglmrhs(isymg,isymm)+1+nrhffr(isymm)*nbas(isymg) 2649 koffc = iglmrh(isymg,isymm)+1 2650 koffg = krhtf + irgijs(isymg,isymkl) 2651 koff1 = kmkl + imatf(isymkl,isymm) 2652 2653 ntotg = max(nbas(isymg),1) 2654 ntotm = max(nrhf(isymm),1) 2655 2656 call dgemm('T','N',nrhf(isymm),nmatijm(isymkl),nbas(isymg), 2657 & one,cmo(koffc),ntotg,work(koffg),ntotg,zero, 2658 & work(koff1),ntotm) 2659 end do 2660c 2661 do isymkl = 1, nsym 2662 isymm = muld2h(isymkl,isymd) 2663 isymdm = muld2h(isymd,isymm) 2664 do isymk = 1, nsym 2665 isyml = muld2h(isymkl,isymk) 2666 do k = 1, nrhfsa(isymk) 2667 do l = 1, nrhfsa(isyml) 2668 idxkl = imatijm(isymk,isyml)+nrhfsa(isymk)*(l-1)+k 2669 idxlk = imatijm(isyml,isymk)+nrhfsa(isyml)*(k-1)+l 2670 iadr1 = kgmakl-1+ivajkls(isymdm,isymkl)+ 2671 & nalphaj(isymdm)*(idxlk-1)+ 2672 & ialphaj(isymd,isymm)+idelta 2673 iadr2 = kmkl+imatf(isymkl,isymm)+nrhf(isymm)*(idxkl-1) 2674 call dcopy(nrhf(isymm),work(iadr2),1, 2675 & work(iadr1),mbas1(isymd)) 2676 end do 2677 end do 2678 end do 2679 end do 2680c 2681 end do 2682 end do 2683 2684 if (locdbg) then 2685 write(lupri,*)'work(krmnkl)' 2686c write(lupri,*)(work(krmnkl-1+i),i=1,ngamsm(1)) 2687 write(lupri,*)'norm^2 work(krmnkl):', 2688 & ddot(ngamsm(1),work(krmnkl),1,work(krmnkl),1) 2689 write(lupri,*)'work(kgmakl)' 2690c write(lupri,*)(work(kgmakl-1+i),i=1,nvajkls(1)) 2691 write(lupri,*)'norm^2 work(kgmakl):', 2692 & ddot(nvajkls(1),work(kgmakl),1,work(kgmakl),1) 2693 end if 2694c 2695c make Vajkl 2696c 2697 factor = one 2698 do isymnm = 1, nsym 2699 isymkl = isymnm 2700 isymaj = muld2h(isymv,isymkl) 2701 2702 koffv = 1 + ivajkl(isymaj,isymkl) 2703 koffr = krmnkl + igamsm(isymnm,isymkl) 2704 koffg = kgmakl + ivajkls(isymaj,isymnm) 2705 2706 ntotmn = max(1,nmatijm(isymnm)) 2707 ntotaj = max(1,nalphaj(isymaj)) 2708 2709 call dgemm('N','N',nalphaj(isymaj),nmatkl(isymkl), 2710 & nmatijm(isymnm),factor,work(koffg),ntotaj, 2711 & work(koffr),ntotmn,one,vajkl(koffv),ntotaj) 2712 if (locdbg) then 2713 write(lupri,*)'DEBUG: Vajkl in mkvaj2' 2714 call output(vajkl(koffv),1,nalphaj(isymaj),1, 2715 & nmatkl(isymkl),nalphaj(isymaj), 2716 & nmatkl(isymkl),1,lupri) 2717 end if 2718 end do 2719c 2720 call qexit('r12mkvaj2') 2721 end 2722*=====================================================================* 2723 subroutine ccrhs_hp(omega1,lamdh,isymh,lamdh0,isymh0,work, 2724 & lwork,iamp,isymc,fc12am,lufc12,ifile,iopte) 2725c--------------------------------------------------------------------- 2726c purpose: update Omega_ai - = 2727c \sum_alpha Lambda^h_alphaa * Omega_alphai 2728c 2729c H. Fliegl, C. Haettig, spring 2004 2730c adapted for CABS, C. Neiss, winter 2006 2731c 2732c iopte = 1 make E_ab^(R12) intermediate instead of Omega_ai^H' 2733c C. Neiss, spring 2006 2734c--------------------------------------------------------------------- 2735 implicit none 2736#include "priunit.h" 2737#include "ccorb.h" 2738#include "ccsdsym.h" 2739#include "r12int.h" 2740#include "ccr12int.h" 2741#include "dummy.h" 2742 integer lwork,isymc,isymh,ifile,lufc12,komai, isymom, 2743 & isymh0,kctil,kend1,lwrk1,isymd,isymilk,isymimn, 2744 & idelp,kint,kgtf,kmint,kend2,lwrk2,isymkl, 2745 & isymi,isymmn,isyma,isymb,koffr,koffm,koffom, 2746 & koffom1,ntota,ntotb,koffc,koffg,ntotkl,ntotmn, 2747 & koffl,ibasx(8),isym,nbas2(8),basoff(8) 2748 integer nr1orb(8),nr1bas(8),nr1xbas(8),nr2bas,n2bst1(8) 2749 integer ir1xbas(8,8) 2750 integer ir1orb(8,8),ir1bas(8,8),ir2bas(8,8),iaodis1(8,8) 2751 integer ir2xbas(8,8),irgkl(8,8),nrgkl(8),nrxgkl(8),irxgkl(8,8) 2752 integer nr1xorb(8),ir1xorb(8,8) 2753 integer idxdp,ntoti,nalphaj(8),ialphaj(8,8),nr1(8) 2754 integer icount,isym1,isym2,nmaijx(8),imaijx(8,8) 2755 integer nmaklx(8),imaklx(8,8),imatax(8,8),nmatax(8) 2756 integer nrgij(8),irgij(8,8),it1xao(8,8),nt1xao(8),it2xaos(8,8) 2757 integer iamp,iopte 2758#if defined (SYS_CRAY) 2759 real work(*), omega1(*), lamdh(*), lamdh0(*), one 2760 real zero,ddot 2761#else 2762 double precision work(*), omega1(*), lamdh(*),lamdh0(*), one 2763 double precision zero,ddot 2764#endif 2765 parameter(one = 1.0d0, zero = 0.0d0) 2766 character*(*) fc12am 2767 character cdummy*(8) 2768 logical locdbg, lauxd 2769 parameter(locdbg = .false.) 2770 2771 call qenter('ccrhs_hp') 2772 2773 isymom = muld2h(isymh,isymc) 2774 2775 if (locdbg) then 2776 if (iopte.eq.1) then 2777 write(lupri,*)'in ccrhs_hp: E1IM on entry ' 2778 call cc_prei(omega1,dummy,isymom,1) 2779 else 2780 write(lupri,*)'in ccrhs_hp: omega1 on entry ' 2781 call cc_prp(omega1,dummy,isymom,1,0) 2782 end if 2783 end if 2784 2785 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 2786 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas,ir2bas, 2787 & ir2xbas,irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 2788 2789 if (r12cbs) then 2790 do isym = 1, nsym 2791 nbas2(isym) = mbas1(isym) + mbas2(isym) 2792 basoff(isym) = 0 2793 end do 2794 else 2795 call icopy(8,mbas2,1,nbas2,1) 2796 call icopy(8,mbas1,1,basoff,1) 2797 end if 2798 2799 if (iopte.eq.1) then 2800 call icopy(8,nvir,1,nr1,1) 2801 else 2802 call icopy(8,nrhf,1,nr1,1) 2803 end if 2804 2805 do isym = 1, nsym 2806 nmaijx(isym) = 0 2807 nmaklx(isym) = 0 2808 nmatax(isym) = 0 2809 do isym2 = 1,nsym 2810 isym1 = muld2h(isym,isym2) 2811 imaijx(isym1,isym2) = nmaijx(isym) 2812 imaklx(isym1,isym2) = nmaklx(isym) 2813 imatax(isym1,isym2) = nmatax(isym) 2814 nmaijx(isym) = nmaijx(isym) + nmatij(isym1)*nr1(isym2) 2815 nmaklx(isym) = nmaklx(isym) + nmatkl(isym1)*nr1(isym2) 2816 nmatax(isym) = nmatax(isym) + mbas1(isym1)*nr1(isym2) 2817 end do 2818 end do 2819 do isym = 1, nsym 2820 nrgij(isym) = 0 2821 nt1xao(isym)= 0 2822 do isym2 = 1, nsym 2823 isym1 = muld2h(isym,isym2) 2824 irgij(isym1,isym2) = nrgij(isym) 2825 it1xao(isym1,isym2) = nt1xao(isym) 2826 nrgij(isym) = nrgij(isym) + mbas1(isym1)*nmatij(isym2) 2827 nt1xao(isym)= nt1xao(isym) + mbas2(isym1)*nrhf(isym2) 2828 end do 2829 end do 2830 do isym = 1, nsym 2831 icount = 0 2832 do isym2 = 1, nsym 2833 isym1 = muld2h(isym,isym2) 2834 it2xaos(isym1,isym2) = icount 2835 icount = icount + nt1ao(isym1)*nt1xao(isym2) 2836 end do 2837 end do 2838 2839 komai = 1 2840 kctil = komai + nmatax(isymom) 2841 kend1 = kctil + ntr12sq(isymc) 2842 lwrk1 = lwork - kend1 2843 if (lwrk1 .lt.0) then 2844 call quit('Insufficient work space in ccrhs_hp') 2845 end if 2846 2847 ibasx(1) = 0 2848 do isym = 2,nsym 2849 ibasx(isym) = ibasx(isym-1)+mbas2(isym-1) 2850 end do 2851c 2852 call dzero(work(komai),nmatax(isymom)) 2853 2854c read r12 amplitudes or trial vector 2855 call cc_r12getct(work(kctil),isymc,iamp,ketscl,.true.,'T', 2856 & lufc12,fc12am,ifile,cdummy,idummy,work(kend1), 2857 & lwrk1) 2858 2859 do isymd = 1, nsym 2860 isymilk = muld2h(isymh,isymd) 2861 isymimn = muld2h(isymc,isymilk) 2862 kint = kend1 2863 kgtf = kint + max(nrgkl(isymd),nrgij(isymd)) 2864 kmint = kgtf + nmaijx(isymilk) 2865 kend2 = kmint + nmaklx(isymimn) 2866 lwrk2 = lwork - kend2 2867 if (lwrk2.lt.0) then 2868 call quit('Insufficient work space in ccrhs_hp') 2869 end if 2870c 2871 do idelp = 1, nbas2(isymd) 2872 call dzero(work(kgtf),nmaijx(isymilk)) 2873c get I^dp_akl 2874 lauxd = .false. 2875 idxdp = idelp+ibas(isymd)+basoff(isymd) 2876 call cc_r12getrint(work(kint),idxdp,isymd,nt1ao,it1ao, 2877 & nt2aos,it2aos,nrgij,irgij,it1xao,it2xaos, 2878 & nrhf,nmatij,imatij, 2879 & ibasx,lauxd,.false.,fgdpakl,work(kend2),lwrk2) 2880c if IOPTE.EQ.1 transform I^dp_akl to I^dp_bkl 2881 if (iopte.eq.1) then 2882 call cc_r12generaltf(work(kint),work(kgtf),idelp,isymd, 2883 & lamdh,isymh,lamdh,isymh,.true.,iglmvi, 2884 & iglmvi,nmatij,nr1,nbas,irgij,imatij, 2885 & nmatij,imaijx,idummy,work(kend2),lwrk2) 2886 else 2887c transform I^dp_akl to I^dp_ikl 2888 call cc_r12generaltf(work(kint),work(kgtf),idelp,isymd, 2889 & lamdh,isymh,lamdh,isymh,.true.,iglmrh, 2890 & iglmrh,nmatij,nrhf,nbas,irgij,imatij, 2891 & nmatij,imaijx,idummy,work(kend2),lwrk2) 2892 end if 2893c 2894 if (locdbg) then 2895 write(lupri,*)'norm^2 coulomb integrals:', 2896 & ddot(nrgij(isymd),work(kint),1,work(kint),1) 2897 write(lupri,*)'norm^2 kgtf:', 2898 & ddot(nmaijx(isymilk),work(kgtf),1,work(kgtf),1) 2899 end if 2900 2901c transform M^dpi_mn = /sum_kl C^tilde * I^dp_ikl 2902 do isymkl = 1, nsym 2903 isymi = muld2h(isymilk,isymkl) 2904 isymmn = muld2h(isymimn,isymi) 2905 2906 koffc = kctil + itr12sqt(isymkl,isymmn) 2907 koffg = kgtf + imaijx(isymkl,isymi) 2908 koffm = kmint + imaklx(isymmn,isymi) 2909 2910 ntotkl = max(nmatij(isymkl),1) 2911 ntotmn = max(nmatkl(isymmn),1) 2912 ntoti = max(nr1(isymi),1) 2913 2914 call dgemm('T','T',nmatkl(isymmn),nr1(isymi), 2915 & nmatij(isymkl),one,work(koffc),ntotkl, 2916 & work(koffg),ntoti,zero,work(koffm),ntotmn) 2917 end do 2918 if (locdbg) then 2919 write(lupri,*)'norm^2 kmint:', 2920 & ddot(nmaklx(isymimn),work(kmint),1,work(kmint),1) 2921 end if 2922 2923c get R^dp_amn 2924 lauxd = .false. 2925 idxdp = idelp+ibas(isymd)+basoff(isymd) 2926 call cc_r12getrint(work(kint),idxdp,isymd,nr1bas,ir1bas, 2927 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 2928 & nrhfb,nmatkl,imatkl, 2929 & ibasx,lauxd,.false.,fnback2,work(kend2),lwrk2) 2930 2931 if (locdbg) then 2932 write(lupri,*)'norm^2 kint:', 2933 & ddot(nrgkl(isymd),work(kint),1,work(kint),1) 2934 end if 2935 2936 do isymmn = 1, nsym 2937 isyma = muld2h(isymmn,isymd) 2938 isymi = muld2h(isymom,isyma) 2939 2940 koffr = kint + irgkl(isyma,isymmn) 2941 koffm = kmint + imaklx(isymmn,isymi) 2942 koffom = komai + imatax(isyma,isymi) 2943 2944 ntota = max(mbas1(isyma),1) 2945 ntotmn = max(nmatkl(isymmn),1) 2946 2947 call dgemm('N','N',mbas1(isyma),nr1(isymi), 2948 & nmatkl(isymmn),-one,work(koffr),ntota, 2949 & work(koffm),ntotmn,one,work(koffom),ntota) 2950 end do 2951 end do 2952 end do 2953 if (locdbg) then 2954 write(lupri,*)'norm^2 komai:', 2955 & ddot(nmatax(isymom),work(komai),1,work(komai),1) 2956 end if 2957c 2958c update Omega_ai 2959 do isyma = 1, nsym 2960 isymb = muld2h(isymh0,isyma) 2961 isymi = muld2h(isymom,isyma) 2962 koffl = 1 + iglmvi(isyma,isymb) 2963 koffom = komai + imatax(isyma,isymi) 2964 if (iopte.eq.1) then 2965 koffom1 = 1 + imatab(isymb,isymi) 2966 else 2967 koffom1 = 1 + it1am(isymb,isymi) 2968 end if 2969 2970 ntota = max(mbas1(isyma),1) 2971 ntotb = max(nvir(isymb),1) 2972 2973 call dgemm('T','N',nvir(isymb),nr1(isymi),mbas1(isyma),one, 2974 & lamdh0(koffl),ntota,work(koffom),ntota,one, 2975 & omega1(koffom1),ntotb) 2976 end do 2977c 2978 if (locdbg) then 2979 if (iopte.eq.1) then 2980 write(lupri,*)'in ccrhs_hp: E1IM on exit ' 2981 call cc_prei(omega1,dummy,isymom,1) 2982 else 2983 write(lupri,*)'in ccrhs_hp: omega1 on exit ' 2984 call cc_prp(omega1,dummy,isymom,1,0) 2985 end if 2986 end if 2987 2988 call qexit('ccrhs_hp') 2989 end 2990*=====================================================================* 2991 subroutine ccrhs_ip(omega1,t1am,isymt,lamdh,isymh,iamp,isymc, 2992 & fc12am,lufc12,ifile,work,lwork) 2993c--------------------------------------------------------------------- 2994c purpose: update Omega_ai = /sum_klm ctilde * M^ak_mn 2995c 2996c H. Fliegl, C. Haettig, spring 2004 2997c adapted for CABS, C. Neiss, winter 2006 2998c--------------------------------------------------------------------- 2999 implicit none 3000#include "priunit.h" 3001#include "ccorb.h" 3002#include "ccsdsym.h" 3003#include "r12int.h" 3004#include "ccr12int.h" 3005#include "dummy.h" 3006 logical lauxd, locdbg 3007 parameter (locdbg = .false.) 3008 integer lwork,ifile,lufc12,isymt,isymh,isymc,isymt1ao, 3009 & isymak,icount7,isymk,isyma,kt1ao,kmint,kend1, 3010 & lwrk1,isymb,isyml,koffl,kofft,kofftam,ntota, 3011 & ntotb,isymd,idelp,kint,kfint,kend2,lwrk2, 3012 & ibasx(8),isymkl,idxkl,idxlk,kofft1am,koffgc, 3013 & koffge,isymmn,idxmn,koffr,koffm,kctil,isymomg,isym 3014 integer nr1orb(8),nr1bas(8),nr1xbas(8),nr2bas,n2bst1(8) 3015 integer ir1xbas(8,8),nbas2(8),basoff(8) 3016 integer ir1orb(8,8),ir1bas(8,8),ir2bas(8,8),iaodis1(8,8) 3017 integer ir2xbas(8,8),irgkl(8,8),nrgkl(8),nrxgkl(8),irxgkl(8,8) 3018 integer nr1xorb(8),ir1xorb(8,8),nalphaj(8),ialphaj(8,8) 3019 integer idxdp,nrgij(8),irgij(8,8), 3020 & nt1xao(8),it1xao(8,8),it2xaos(8,8),isymdk,icount,icoun1 3021 integer iamp 3022#if defined (SYS_CRAY) 3023 real omega1(*), work(*), lamdh(*), t1am(*), 3024 & ddot, er12, one, zero 3025#else 3026 double precision omega1(*), work(*), lamdh(*), t1am(*), 3027 & ddot, er12, one, zero 3028#endif 3029 parameter (one = 1.0d0, zero = 0.0d0) 3030 character*(*) fc12am 3031 character cdummy*(8) 3032 3033 call qenter('ccrhs_ip') 3034 3035 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 3036 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas, 3037 & ir2bas,ir2xbas,irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 3038 3039 isymt1ao = muld2h(isymt,isymh) 3040 3041 if (locdbg) then 3042 write(lupri,*)'in ccrhs_ip: omega1 on entry ' 3043 call cc_prp(omega1,dummy,isymt1ao,1,0) 3044 end if 3045c 3046 if (r12cbs) then 3047 do isym = 1, nsym 3048 nbas2(isym) = mbas1(isym) + mbas2(isym) 3049 basoff(isym) = 0 3050 end do 3051 else 3052 call icopy(8,mbas2,1,nbas2,1) 3053 call icopy(8,mbas1,1,basoff,1) 3054 end if 3055c 3056 ibasx(1) = 0 3057 do isym = 2,nsym 3058 ibasx(isym) = ibasx(isym-1)+mbas2(isym-1) 3059 end do 3060 do isymdk = 1, nsym 3061 icount = 0 3062 icoun1 = 0 3063 do isymk = 1, nsym 3064 isymd = muld2h(isymdk,isymk) 3065 irgij(isymd,isymk) = icount 3066 it1xao(isymd,isymk) = icoun1 3067 icount = icount + mbas1(isymd)*nmatij(isymk) 3068 icoun1 = icoun1 + mbas2(isymd)*nrhf(isymk) 3069 end do 3070 nrgij(isymdk) = icount 3071 nt1xao(isymdk) = icoun1 3072 end do 3073 do isymdk = 1, nsym 3074 icount = 0 3075 do isymk = 1, nsym 3076 isymd = muld2h(isymdk,isymk) 3077 it2xaos(isymd,isymk) = icount 3078 icount = icount + nt1ao(isymd)*nt1xao(isymk) 3079 end do 3080 end do 3081c 3082 kt1ao = 1 3083 kmint = kt1ao + nt1ao(isymt1ao) 3084 kend1 = kmint + nvajkl(isymt1ao) 3085 lwrk1 = lwork - kend1 3086 if (lwrk1 .lt.0) then 3087 call quit('Insufficient work space in ccrhs_ip') 3088 end if 3089c 3090 call dzero(work(kmint),nvajkl(isymt1ao)) 3091 call dzero(work(kt1ao),nt1ao(isymt1ao)) 3092c 3093c calculate t_al = /sum_b t_bl * Lamda^h_ab 3094 do isymb = 1, nsym 3095 isyma = muld2h(isymh,isymb) 3096 isyml = muld2h(isymt,isymb) 3097 3098 koffl = 1 + iglmvi(isyma,isymb) 3099 kofft = 1 + it1am(isymb,isyml) 3100 kofftam = kt1ao + it1ao(isyma,isyml) 3101 3102 ntota = max(mbas1(isyma),1) 3103 ntotb = max(nvir(isymb),1) 3104 3105 call dgemm('N','N',mbas1(isyma),nrhf(isyml),nvir(isymb),one, 3106 & lamdh(koffl),ntota,t1am(kofft),ntotb,zero, 3107 & work(kofftam),ntota) 3108 end do 3109 3110 do isymd = 1, nsym 3111 isymk = muld2h(isymt1ao,isymd) 3112 do idelp = 1, nbas2(isymd) 3113 kint = kend1 3114 kfint = kint + nrgkl(isymd) 3115 kend2 = kfint + nrhf(isymk) 3116 lwrk2 = lwork - kend2 3117 if (lwrk2.lt.0) then 3118 call quit('Insufficient work space in ccrhs_hp') 3119 end if 3120 3121c get I^dp_akl 3122 lauxd = .false. 3123 idxdp = idelp+ibas(isymd)+basoff(isymd) 3124 call cc_r12getrint(work(kint),idxdp,isymd,nt1ao,it1ao, 3125 & nt2aos,it2aos,nrgij,irgij,it1xao,it2xaos, 3126 & nrhf,nmatij,imatij, 3127 & ibasx,lauxd,.false.,fgdpakl,work(kend2),lwrk2) 3128 if (locdbg) then 3129 write(lupri,*)'norm^2 kint:', 3130 & ddot(nrgij(isymd),work(kint),1,work(kint),1) 3131 end if 3132 3133c get F intermediate 3134 call dzero(work(kfint),nrhf(isymk)) 3135 do isyml = 1, nsym 3136 isyma = muld2h(isymt1ao,isyml) 3137 isymkl = muld2h(isymk,isyml) 3138c 3139 do k = 1, nrhf(isymk) 3140 do l = 1, nrhf(isyml) 3141 idxkl = imatij(isymk,isyml)+nrhf(isymk)*(l-1)+k 3142 idxlk = imatij(isyml,isymk)+nrhf(isyml)*(k-1)+l 3143 3144 kofft1am = kt1ao + it1ao(isyma,isyml)+ 3145 & mbas1(isyma)*(l-1) 3146 koffgc = kint + irgij(isyma,isymkl)+ 3147 & mbas1(isyma)*(idxlk-1) 3148 koffge = kint + irgij(isyma,isymkl)+ 3149 & mbas1(isyma)*(idxkl-1) 3150 3151 work(kfint-1+k) = work(kfint-1+k) + 3152 & 2.0d0*ddot(mbas1(isyma), 3153 & work(kofft1am),1,work(koffgc),1) - 3154 & ddot(mbas1(isyma),work(kofft1am),1, 3155 & work(koffge),1) 3156 end do 3157 end do 3158 end do 3159 if (locdbg) then 3160 write(lupri,*)'norm^2: kfint', 3161 & ddot(nrhf(isymk),work(kfint),1,work(kfint),1) 3162 end if 3163 3164c get R^dp_amn 3165 lauxd = .false. 3166 idxdp = idelp+ibas(isymd)+basoff(isymd) 3167 call cc_r12getrint(work(kint),idxdp,isymd,nr1bas,ir1bas, 3168 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 3169 & nrhfb,nmatkl,imatkl, 3170 & ibasx,lauxd,.false.,fnback2,work(kend2),lwrk2) 3171 3172c get M intermediate 3173 do isymmn = 1, nsym 3174 isyma = muld2h(isymd,isymmn) 3175 isymak = muld2h(isyma,isymk) 3176 do k = 1, nrhf(isymk) 3177 do idxmn = 1, nmatkl(isymmn) 3178 koffr = kint+irgkl(isyma,isymmn)+ 3179 & mbas1(isyma)*(idxmn-1) 3180 koffm = kmint+ivajkl(isymak,isymmn)+ 3181 & nt1ao(isymak)*(idxmn-1)+it1ao(isyma,isymk)+ 3182 & mbas1(isyma)*(k-1) 3183 3184 call daxpy(mbas1(isyma),work(kfint-1+k),work(koffr), 3185 & 1,work(koffm),1) 3186 end do 3187 end do 3188 end do 3189c 3190 end do 3191 end do 3192 if (locdbg) then 3193 write(lupri,*) 'norm^2(M):',ddot(nvajkl(isymt1ao), 3194 & work(kmint),1,work(kmint),1) 3195 end if 3196 3197 kctil = kend1 3198 kend2 = kctil + ntr12sq(isymc) 3199 lwrk2 = lwork - kend2 3200 if (lwrk2.lt.0) then 3201 call quit('Insufficient work space in ccrhs_ip') 3202 end if 3203 3204c read r12 amplitudes or trial vector 3205 call cc_r12getct(work(kctil),isymc,iamp,ketscl,.true.,'T', 3206 & lufc12,fc12am,ifile,cdummy,idummy,work(kend2), 3207 & lwrk2) 3208 3209c isymomg = muld2h(isymc,muld2h(isymt1ao,isymh)) 3210 isymomg = muld2h(isymh,muld2h(isymt1ao,isymc)) 3211 call ccrhs_gp0(omega1,isymomg,work(kmint),isymt1ao,lamdh,isymh, 3212 & work(kctil),isymc,.false.,er12,locdbg,work(kend2),lwrk2) 3213 3214 if (locdbg) then 3215 write(lupri,*)'in ccrhs_ip: omega1 on exit ' 3216 call cc_prp(omega1,dummy,isymt1ao,1,0) 3217 end if 3218 3219 call qexit('ccrhs_ip') 3220 end 3221*======================================================================* 3222 subroutine cc_r12offs23(iglmrhs,iglmvis,nglmds,icmo,ncmo, 3223 & imaijm,nmaijm, 3224 & imaklm,nmaklm,imatijm,nmatijm, 3225 & igamsm,ngamsm,irgijs,nrgijs, 3226 & ir1basm,nr1basm,ir2basm,nr2basm,ir1xbasm, 3227 & nr1xbasm,ir2xbasm,imatf,nmatf) 3228*----------------------------------------------------------------------* 3229c purpose: calculate some offsets needet for ansatz 2 & 3 to include 3230c all active and inactive orbitals 3231c 3232c H. Fliegl, C. Haettig, spring 2004 3233*----------------------------------------------------------------------* 3234 implicit none 3235#include "priunit.h" 3236#include "ccorb.h" 3237#include "ccsdsym.h" 3238#include "r12int.h" 3239 integer iglmrhs(8,8),iglmvis(8,8),nglmds(8),isym, 3240 & icmo(8,8),ncmo(8),imaklm(8,8),nmaklm(8), 3241 & imatijm(8,8),nmatijm(8),ngamsm(8),igamsm(8,8), 3242 & isym2,isym1,irgijs(8,8),nrgijs(8),ir1basm(8,8),nr1basm(8), 3243 & ir2basm(8,8),nr2basm,ir1xbasm(8,8),nr1xbasm(8), 3244 & ir2xbasm(8,8),imatf(8,8),nmatf(8),imaijm(8,8),nmaijm(8), 3245 & icount1,icount2,icount2a,icount3,icount4,icount5, 3246 & icount6,icount7 3247 3248 call qenter('ccoffs23') 3249 3250 do isym = 1, nsym 3251 3252 ! precalculate dimensions of transformation matrices 3253 ! and the offset for the virtual block 3254 icount1 = 0 3255 icount2a= 0 3256 do isym2 = 1, nsym 3257 isym1 = muld2h(isym,isym2) 3258 3259 icmo(isym1,isym2) = icount1 3260 3261 icount1 = icount1 + nbas(isym1)*norbs(isym2) 3262 icount2a= icount2a+ nbas(isym1)*nrhfsa(isym2) 3263 end do 3264 ncmo(isym) = icount1 3265 nglmds(isym) = icount1 3266 3267 icount2 = 0 3268 icount3 = 0 3269 icount4 = 0 3270 icount5 = 0 3271 icount6 = 0 3272 icount7 = 0 3273 do isym2 = 1, nsym 3274 isym1 = muld2h(isym,isym2) 3275 iglmrhs(isym1,isym2) = icount2 3276 iglmvis(isym1,isym2) = icount2a 3277 imaklm(isym1,isym2) = icount3 3278 imaijm(isym1,isym2) = icount7 3279 imatijm(isym1,isym2) = icount4 3280 ir1basm(isym1,isym2) = icount5 3281 ir1xbasm(isym1,isym2) = icount6 3282 3283 icount2 = icount2 + nbas(isym1)*nrhfsa(isym2) 3284 icount2a= icount2a+ nbas(isym1)*nvirs(isym2) 3285 icount3 = icount3 + nmatkl(isym1)*nrhfsa(isym2) 3286 icount7 = icount7 + nmatij(isym1)*nrhfsa(isym2) 3287 icount4 = icount4 + nrhfs(isym1)*nrhfsa(isym2) 3288 icount5 = icount5 + mbas1(isym1)*nrhfsa(isym2) 3289 icount6 = icount6 + mbas2(isym1)*nrhfsa(isym2) 3290 end do 3291 nmaklm(isym) = icount3 3292 nmaijm(isym) = icount7 3293 nmatijm(isym) = icount4 3294 nr1basm(isym) = icount5 3295 nr1xbasm(isym) = icount6 3296 end do 3297 3298 do isym = 1, nsym 3299 icount7 = 0 3300 do isym2 = 1, nsym 3301 isym1 = muld2h(isym,isym2) 3302 imatf(isym1,isym2) = icount7 3303 icount7 = icount7 + nmatijm(isym1)*nrhf(isym2) 3304 end do 3305 nmatf(isym) = icount7 3306 end do 3307 3308 do isym = 1, nsym 3309 icount1 = 0 3310 icount2 = 0 3311 icount3 = 0 3312 icount4 = 0 3313 do isym2 = 1, nsym 3314 isym1 = muld2h(isym,isym2) 3315 igamsm(isym1,isym2) = icount1 3316 irgijs(isym1,isym2) = icount2 3317 ir2basm(isym1,isym2) = icount3 3318 ir2xbasm(isym1,isym2) = icount4 3319 3320 icount1 = icount1 + nmatijm(isym1)*nmatkl(isym2) 3321 icount2 = icount2 + mbas1(isym1)*nmatijm(isym2) 3322 icount3 = icount3 + nr1basm(isym1)*nr1basm(isym2) 3323 icount4 = icount4 + nr1basm(isym1)*nr1xbasm(isym2) 3324 end do 3325 ngamsm(isym) = icount1 3326 nrgijs(isym) = icount2 3327c nr2basm(isym) = icount3 3328 end do 3329c 3330 nr2basm = 0 3331 do isym = 1, nsym 3332 nr2basm = nr2basm + nr1basm(isym)*nr1basm(isym) 3333 end do 3334 3335 call qexit('ccoffs23') 3336 3337 end 3338*=====================================================================* 3339 subroutine cc_r12mkgmnab(xint,cmos,isymc,idel,isymd, 3340 & iglmrhs,nglmds,lunitr12,filer12, 3341 & work,lwork) 3342c--------------------------------------------------------------------- 3343c purpose: make g^MN_alpha_beta intermediate 3344c M and N run over all active and inactive occupied 3345c molecular orbitals 3346c 3347c H. Fliegl, C. Haettig, spring 2004 3348c--------------------------------------------------------------------- 3349 implicit none 3350#include "priunit.h" 3351#include "ccsdsym.h" 3352#include "ccorb.h" 3353#include "r12int.h" 3354 logical locdbg 3355 parameter (locdbg = .false.) 3356 integer icount1,isym,isym1,isym2 3357 integer idelta,idel,isymd,iglmrhs(8,8),nglmds(8),lunitr12,lwork, 3358 & isymn,isymdn,isymab,isymbm,koff1,kscr1,kend1,lwrk1,koffg, 3359 & koffc,ntotab,ntotg,idxabn,isymm,isymb,isyma,koff5,kscr2, 3360 & koff3,koff4,ntota,ntotb,isymc,isymg 3361 integer nrhftsorb 3362 3363 character*(*) filer12 3364 3365#if defined (SYS_CRAY) 3366 real work(*), cmos(*), xint(*), one, zero, xnrm, ddot 3367 real xnrm1, xnrm2 3368#else 3369 double precision work(*), cmos(*), xint(*), one, zero, xnrm, ddot 3370 double precision xnrm1, xnrm2 3371#endif 3372 parameter (one = 1.0d0, zero = 0.0d0) 3373 3374 idelta = idel - ibas(isymd) 3375 3376c skip if idelta not in orbital basis 3377 if (idelta.gt.mbas1(isymd)) return 3378 3379 call qenter('r12mkgmnab') 3380 3381 nrhftsorb = 0 3382 do isym = 1, nsym 3383 nrhftsorb = nrhftsorb + nrhfsa(isym) 3384 end do 3385c 3386c transform first index 3387c 3388 if (locdbg) then 3389 xnrm = 0.0d0 3390 xnrm1 = 0.0d0 3391 xnrm2 = 0.0d0 3392 end if 3393c 3394 do isymg = 1, nsym 3395 isymn = muld2h(isymc,isymg) ! isymn = isymg 3396 isymdn = muld2h(isymd,isymn) 3397 isymab = isymdn 3398 isymbm = isymdn 3399 3400 koff1 = 1 3401 kscr1 = koff1 + nnbst(isymab)*nrhfsa(isymn) 3402 kscr2 = kscr1 + nbast*nbast 3403 kend1 = kscr2 + nbast*nrhftsorb 3404 lwrk1 = lwork - kend1 3405 if (lwrk1.lt.0) then 3406 call quit('Insufficient work space in cc_r12mkgmnab') 3407 end if 3408 3409 koffg = idsaog(isymg,isymd) + 1 3410 koffc = iglmrhs(isymg,isymn) + 1 3411 3412 ntotab = max(nnbst(isymab),1) 3413 ntotg = max(nbas(isymg),1) 3414 3415 call dgemm('N','N',nnbst(isymab),nrhfsa(isymn),mbas1(isymg), 3416 & one,xint(koffg),ntotab,cmos(koffc),ntotg,zero, 3417 & work(koff1),ntotab) 3418c 3419 if (locdbg) then 3420 xnrm1 = xnrm1+ddot(nnbst(isymab)*nrhfsa(isymn), 3421 & work(koff1),1,work(koff1),1) 3422 end if 3423c 3424 3425 do n = 1, nrhfsa(isymn) 3426 idxabn = koff1 + nnbst(isymab)*(n-1) 3427c unpack integrals as a quadratic matrix 3428 call ccsd_symsq(work(idxabn),isymab,work(kscr1)) 3429c 3430 if (locdbg) then 3431 xnrm2=xnrm2+ddot(n2bst(isymab),work(kscr1),1,work(kscr1),1) 3432 end if 3433 3434c transform second index to M 3435 koff5 = kscr2 3436 do isyma = 1, nsym 3437 isymm = muld2h(isymc,isyma) !isymm = isyma 3438 isymb = muld2h(isymab,isyma) 3439 3440 koff3 = kscr1 + iaodis(isyma,isymb) 3441 koff4 = iglmrhs(isyma,isymm) + 1 3442 3443 ntota = max(nbas(isyma),1) 3444 ntotb = max(nbas(isymb),1) 3445 3446 call dgemm('T','N',nbas(isymb),nrhfsa(isymm),mbas1(isyma), 3447 & one,work(koff3),ntota,cmos(koff4),ntota,zero, 3448 & work(koff5),ntotb) 3449c 3450 if (locdbg) then 3451 xnrm = xnrm+ddot(nbas(isymb)*nrhfsa(isymm),work(koff5),1, 3452 & work(koff5),1) 3453 end if 3454c 3455 3456 koff5 = koff5 + nbas(isymb)*nrhfsa(isymm) 3457 end do 3458 3459c write result on file 3460 call cc_r12whtfs(work(kscr2),idel,isymd,n,isymn,isymbm, 3461 & lunitr12,filer12,work(kend1),lwrk1) 3462 end do 3463 end do 3464c 3465 if (locdbg) then 3466 write(lupri,*) 'norm^2(cmo)',ddot(nglmds,cmos,1,cmos,1) 3467 write(lupri,*) 'Gmbnd for idel=',idel,xnrm 3468 write(lupri,*) 'Gabnd for idel=',idel,xnrm1,xnrm2 3469 end if 3470c 3471 call qexit('r12mkgmnab') 3472 end 3473*=====================================================================* 3474 subroutine cc_r12whtfs(xint,idel,isymd,j,isymj,isymai, 3475 & lunitr12,filer12,work,lwork) 3476c---------------------------------------------------------------------- 3477c purpose: get half transformed coulomb integrals and write them 3478c on file first delta in ccsd_aibj 3479c note: M,N run over all active and inactive occupied 3480c molecular orbitals 3481c 3482c H. Fliegl, C. Haettig, spring 2004 3483c---------------------------------------------------------------------- 3484 implicit none 3485#include "priunit.h" 3486#include "ccorb.h" 3487#include "ccsdsym.h" 3488#include "r12int.h" 3489#include "dummy.h" 3490 3491 logical locdbg 3492 parameter (locdbg = .false.) 3493 integer lunitr12, lwork, idel, isymd, isymi, isyma 3494 integer isymai, isymdj, isymj, ndelj, iadr, ilen 3495 integer idxai, idxdj, idxai2, koff1, idelta 3496 integer isym,isym1,isym2,icount1,ng1bas(8),ig2bas(8,8),ig1bas(8,8) 3497 3498 character*(*) filer12 3499 3500#if defined (SYS_CRAY) 3501 real work(*), xint(*) 3502#else 3503 double precision work(*), xint(*) 3504#endif 3505 3506 call qenter('whtfs') 3507 3508c calculate some offsets and dimensions 3509 do isym = 1, nsym 3510 ng1bas(isym) = 0 3511 icount1 = 0 3512 do isym2 = 1, nsym 3513 isym1 = muld2h(isym,isym2) 3514 ng1bas(isym) = ng1bas(isym)+mbas1(isym1)*nrhfsa(isym2) 3515 ig1bas(isym1,isym2) = icount1 3516 icount1 = icount1 + mbas1(isym1)*nrhfsa(isym2) 3517 end do 3518 end do 3519c 3520 do isym = 1, nsym 3521 icount1 = 0 3522 do isym2 = 1, nsym 3523 isym1 = muld2h(isym,isym2) 3524 ig2bas(isym1,isym2) = icount1 3525 icount1 = icount1 + ng1bas(isym1)*ng1bas(isym2) 3526 end do 3527 end do 3528c 3529 if (lwork.lt.ng1bas(isymai)) then 3530 call quit('Insufficient work space in cc_r12whtfs') 3531 end if 3532 3533 idelta = idel - ibas(isymd) 3534 3535 koff1 = 1 3536 isymdj = muld2h(isymd,isymj) 3537 do isymi = 1, nsym 3538 isyma = muld2h(isymai,isymi) 3539 do i = 1, nrhfsa(isymi) 3540 do a = 1, mbas1(isyma) 3541 idxai = ig1bas(isyma,isymi) + mbas1(isyma)*(i-1)+a 3542 idxai2 = koff1 -1 + nbas(isyma)*(i-1) +a 3543 work(idxai) = xint(idxai2) 3544 end do 3545 end do 3546 koff1 = koff1 + nbas(isyma)*nrhfsa(isymi) 3547 end do 3548 3549 if (locdbg) then 3550 write(lupri,*)'idel=',idel 3551 write(lupri,*)'in whtfs' 3552 do isymi = 1, nsym 3553 isyma = muld2h(isymai,isymi) 3554 write(lupri,*) 'symmetry block (a,i):',isyma,isymi 3555 call output(work(ig1bas(isyma,isymi)+1),1,mbas1(isyma), 3556 & 1,nrhfsa(isymi),mbas1(isyma),nrhfsa(isymi),1,lupri) 3557 end do 3558 end if 3559 3560 if (idelta.le.mbas1(isymd)) then 3561 ndelj = ig1bas(isymd,isymj) + mbas1(isymd)*(j-1) + idelta 3562 iadr = ig2bas(isymai,isymdj) + ng1bas(isymai)*(ndelj-1)+1 3563 else 3564 call quit('There is something wrong in cc_r12whtfs'// 3565 & ' with the auxiliary basis') 3566 end if 3567 ilen = ng1bas(isymai) 3568 3569 call putwa2(lunitr12,filer12,work,iadr,ilen) 3570 3571 call qexit('whtfs') 3572 end 3573*=====================================================================* 3574 subroutine cc_r12getrints(r12bkl,idel,isymd,nr1basm,ir1basm, 3575 & nr2basm,ir2basm,nrgkls,irgkls,ir1xbasm,ir2xbasm, 3576 & ibasx,imatijm,nmatijm,lauxd,fconst,work,lwork) 3577c--------------------------------------------------------------------- 3578c purpose: get g^MN_gamma_delta integrals, needed for ansatz 2 3579c 3580c H. Fliegl, C. Haettig spring 2004 3581c--------------------------------------------------------------------- 3582 implicit none 3583#include "priunit.h" 3584#include "ccorb.h" 3585#include "ccsdsym.h" 3586#include "r12int.h" 3587 3588 logical lauxd,locdbg 3589 parameter(locdbg=.false.) 3590 3591 character*(*) fconst 3592 3593#if defined (SYS_CRAY) 3594 real zero,one 3595#else 3596 double precision zero,one 3597#endif 3598 parameter (zero = 0.0d0, one = 1.0d0) 3599 3600 integer nr1basm(8),ir1basm(8,8),nr2basm,ir2basm(8,8) 3601 integer ir2xbasm(8,8),nrgkls(8),imatijm(8,8) 3602 integer irgkls(8,8),ir1xbasm(8,8),nmatijm(8) 3603 integer idelta,idel,lwork,lwrk1,ibasx(8) 3604 integer lur2back,kscr2,kend2,lwrk2,iadr,ndell,ilen 3605 integer isymb,isymk,isymkl,isyml,idxbk,idxkl,idxbkl 3606 integer isymd,isymdl,isymbk,kend1 3607 3608 3609#if defined (SYS_CRAY) 3610 real r12bkl(*),work(*),ddot 3611#else 3612 double precision r12bkl(*),work(*),ddot 3613#endif 3614 3615 kend1 = 1 3616 3617 idelta = idel - ibas(isymd) 3618 if (lauxd) idelta = idelta - ibasx(isymd) 3619 3620 lur2back = -1 3621 call wopen2(lur2back,fconst,64,0) 3622 3623 do isyml = 1, nsym 3624 isymdl = muld2h(isymd,isyml) 3625 isymbk = isymdl 3626 3627 kscr2 = kend1 3628 kend2 = kscr2 + nr1basm(isymbk) 3629 lwrk2 = lwork - kend2 3630 if (lwrk2 .lt.0) then 3631 write(lupri,*)'lwork',lwork 3632 write(lupri,*)'lwrk2,kend2,kscr2', lwrk2,kend2,kscr2 3633 call quit('Insufficient work space in cc_r12getrints') 3634 end if 3635c ------------------------------------------------------------ 3636c read r^{beta,delta}_{k,l} from file and pack as 3637c R^{delta}(beta,kl) or R^{delta'}(beta,kl) 3638c ------------------------------------------------------------ 3639 do l = 1, nrhfsa(isyml) 3640 if (idelta.le.mbas1(isymd)) then 3641 ndell = ir1basm(isymd,isyml) + 3642 & mbas1(isymd)*(l-1)+idelta 3643 iadr = ir2basm(isymbk,isymdl)+ 3644 & nr1basm(isymbk)*(ndell-1)+1 3645 else 3646 ndell = ir1xbasm(isymd,isyml)+ 3647 & mbas2(isymd)*(l-1)+idelta-mbas1(isymd) 3648 iadr = nr2basm + ir2xbasm(isymbk,isymdl)+ 3649 & nr1basm(isymbk)*(ndell-1)+1 3650 3651 if(locdbg) then 3652 write(lupri,*)'ndell,iadr,l,idelta',ndell,iadr, 3653 & l,idelta 3654 end if 3655 3656 end if 3657 ilen = nr1basm(isymbk) 3658 call getwa2(lur2back,fconst,work(kscr2),iadr,ilen) 3659c 3660 do isymb =1, nsym 3661 isymk = muld2h(isymbk,isymb) 3662 isymkl = muld2h(isymk,isyml) 3663 do k=1, nrhfsa(isymk) 3664 do b=1, mbas1(isymb) 3665 idxbk = ir1basm(isymb,isymk) + mbas1(isymb)*(k-1)+b 3666 idxkl = imatijm(isymk,isyml) +nrhfsa(isymk)*(l-1)+k 3667 idxbkl = irgkls(isymb,isymkl) + 3668 & mbas1(isymb)*(idxkl-1)+b 3669 3670 if(locdbg) then 3671 write(lupri,*)'idxbk,idxkl,idxbkl', 3672 & idxbk,idxkl,idxbkl 3673 end if 3674 3675 r12bkl(idxbkl) = work(kscr2 -1 + idxbk) 3676 3677 end do 3678 end do 3679 3680 if (locdbg) then 3681 write(lupri,'(2a,4i5)')' R^{delta}(beta,kl) for ', 3682 & 'delta,l,isymb,isymkl',idelta,l,isymb,isymkl 3683 idxbkl = irgkls(isymb,isymkl) + 1 3684 call output(r12bkl(idxbkl),1,mbas1(isymb),1, 3685 & nmatijm(isymkl),mbas1(isymb),nmatijm(isymkl),1,lupri) 3686 end if 3687 3688 end do 3689 3690 end do 3691 end do 3692 call wclose2(lur2back,fconst,'KEEP') 3693 3694 if (locdbg) then 3695 write(lupri,*) 'leaving cc_r12getrints... norm^2(r12bkl):', 3696 & ddot(nrgkls(isymd),r12bkl,1,r12bkl,1) 3697 end if 3698 3699 end 3700*=====================================================================* 3701 subroutine cc_r12fdvint(c1am,c2am,work,lwork,aproxr12, 3702 & fc12am,lufc12,ivec) 3703c----------------------------------------------------------------------- 3704c purpose: Test routine to calculate V bar numerically used for 3705c CC2-R12 ansatz 2 3706c 3707c H. Fliegl, C. Haettig, winter 2004 3708c----------------------------------------------------------------------- 3709 implicit none 3710#include "maxorb.h" 3711#include "ccorb.h" 3712#include "ccsdsym.h" 3713#include "priunit.h" 3714#include "dummy.h" 3715#include "second.h" 3716#include "ccsdinp.h" 3717#include "cclr.h" 3718#include "r12int.h" 3719#include "ccr12int.h" 3720 logical omegorsv, testvajtkl, testvbar, testrbar 3721 parameter (testvajtkl = .false., testvbar = .false., 3722 & testrbar = .false.) 3723 3724 integer lwork,kvintp,kvintm,kend1,krho1,krho2,kt1am, 3725 & kt2am,lwrk1,iopt,luvijkl,isymij,isymkl,kvajkl,luvajkl 3726 integer kvsym,kvsing,kvtrip,nrhftria,krp,krm,ngamsm(8), 3727 & imatijm(8,8),nmatijm(8),igamsm(8,8),isym,isym1,isym2, 3728 & icount1,icount4,krsym,lufc12,krho1p,krho2p,nrho2,ivec 3729 integer khintm,khintp 3730#if defined (SYS_CRAY) 3731 real work(*),c1am(*),c2am(*) 3732 real xhalf, xmtwo, delta 3733#else 3734 double precision work(*),c1am(*),ddot,c2am(*) 3735 double precision xhalf, xmtwo, delta, deltai 3736#endif 3737 parameter (xhalf = 0.5d00,xmtwo = -2.0d00, 3738 & delta = 1.0d-7 , deltai = 1.0d00/delta) 3739 character*10 model 3740 character*(*) aproxr12,fc12am 3741 3742C integer index 3743C index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 3744c 3745 call qenter('r12fdvint') 3746 do isym = 1, nsym 3747 icount4 = 0 3748 do isym2 = 1, nsym 3749 isym1 = muld2h(isym,isym2) 3750 imatijm(isym1,isym2) = icount4 3751 icount4 = icount4 + nrhfsa(isym1)*nrhfsa(isym2) 3752 end do 3753 nmatijm(isym) = icount4 3754 end do 3755 do isym = 1, nsym 3756 icount1 = 0 3757 do isym2 = 1, nsym 3758 isym1 = muld2h(isym,isym2) 3759 igamsm(isym1,isym2) = icount1 3760 icount1 = icount1 + nmatijm(isym1)*nmatkl(isym2) 3761 end do 3762 ngamsm(isym) = icount1 3763 end do 3764 3765 if (omegsq) call quit('cclr_fdvint INCOMPATIBLE WITH omegsq.') 3766 if (lcor.or.lsec) 3767 & call quit('cclr_fdvint INCOMPATIBLE WITH fcore AND fsecon.') 3768c 3769 call around( 'in CCLR_FDVINT: making finite diff. v') 3770 omegorsv = omegor 3771 omegor =.true. 3772c 3773c isymtr = 1 ! in cclr.h 3774 3775 krho1 = 1 3776 krho2 = krho1 + nt1amx 3777 kt1am = krho2 + max(nt2amx,nt2ao(isymop),2*nt2ort(isymop)) 3778 kt2am = kt1am + nt1amx 3779 kvintm = kt2am + max(nt2sq(isymop),nt2r12(1)) 3780 kvintp = kvintm + ngamsq(isymtr) 3781 kvajkl = kvintp + ngamsq(isymtr) 3782 kvsym = kvajkl + nvajkl(isymtr) 3783 krm = kvsym + ngamsq(isymtr) 3784 krp = krm + ngamsm(isymtr) 3785 krsym = krp + ngamsm(isymtr) 3786 krho1p = krsym + ngamsm(isymtr) 3787 krho2p = krho1p + nt1amx 3788 khintm = krho2p + max(nt2amx,nt2ao(isymop),2*nt2ort(isymop)) 3789 khintp = khintm + nt1amx 3790 kend1 = khintp + nt1amx 3791 lwrk1 = lwork - kend1 3792 if (lwrk1.lt.0) then 3793 write (lupri,*) 'available: lwork = ',lwork 3794 write (lupri,*) 'needed (at least) = ',kend1 3795 call quit('insufficient work space in cclr_fdvint ') 3796 endif 3797 nrho2 = max(nt2amx,nt2ao(isymop),2*nt2ort(isymop)) 3798 3799 if (testvajtkl) then 3800c calculate V^alpha_j_ kl 3801 call dzero(work(kt1am),nt1am(isymtr)) 3802 rspim = .true. 3803 call ccrhsn(work(krho1),work(krho2),work(kt1am), 3804 & work(kt2am), 3805 * work(kend1),lwrk1,aproxr12) 3806 3807c read Valpajtildekl 3808 luvajkl = -1 3809 call gpopen(luvajkl,fvajtkl,'unknown',' ','unformatted', 3810 & idummy,.false.) 3811 read(luvajkl)(work(kvajkl-1+i),i=1,nvajkl(isymtr)) 3812 call gpclose(luvajkl,'KEEP') 3813 3814 write(lupri,*)'Vajtkl out of fdvint:', 3815 & (work(kvajkl-1+i),i=1,nvajkl(isymtr)) 3816 end if 3817 3818c read the CC reference amplitudes From disk. 3819 iopt = 3 3820 call cc_rdrsp('R0',0,1,iopt,model,work(kt1am),work(kt2am)) 3821 3822 rspim = .false. 3823 3824c--------------------------------------------- 3825c Calculate V bar by finite difference. 3826c--------------------------------------------- 3827 3828c Compute V(t - 0.5 delta * C) 3829 call daxpy(nt1amx,-0.5d0*delta,c1am,1,work(kt1am),1) 3830 call ccrhsn(work(krho1),work(krho2),work(kt1am), 3831 & work(kt2am), 3832 * work(kend1),lwrk1,aproxr12) 3833 3834c Read V(t - 0.5 delta * C) into WORK(KVINTM) 3835 luvijkl = -1 3836 call gpopen(luvijkl,fvijkl,'unknown',' ','unformatted', 3837 & idummy,.false.) 3838 read(luvijkl)(work(kvintm-1+i),i=1,ngamsq(isymtr)) 3839 call gpclose(luvijkl,'KEEP') 3840 3841c Read HP or IP TERM 'CCR12HTST' or 'CCR12ITST' 3842c luvijkl = -1 3843c call gpopen(luvijkl,'CCR12ITST','unknown',' ','unformatted', 3844c & idummy,.false.) 3845c read(luvijkl)(work(khintm-1+i),i=1,nt1amx) 3846c call gpclose(luvijkl,'KEEP') 3847 3848 if (testrbar) then 3849c Read R(t - 0.5 delta * C) into WORK(KRM) 3850 luvijkl = -1 3851 call gpopen(luvijkl,'RMTNTKL','unknown',' ','unformatted', 3852 & idummy,.false.) 3853 read(luvijkl)(work(krm-1+i),i=1,ngamsm(isymtr)) 3854 call gpclose(luvijkl,'KEEP') 3855 end if 3856c Compute V(t + 0.5 delta * C) 3857 call daxpy(nt1amx,delta,c1am,1,work(kt1am),1) 3858 call ccrhsn(work(krho1p),work(krho2p),work(kt1am), 3859 & work(kt2am), 3860 * work(kend1),lwrk1,aproxr12) 3861 3862c Read V(t + 0.5 delta * C) into WORK(KVINTP) 3863 luvijkl = -1 3864 call gpopen(luvijkl,fvijkl,'unknown',' ','unformatted', 3865 & idummy,.false.) 3866 read(luvijkl)(work(kvintp-1+i),i=1,ngamsq(isymtr)) 3867 call gpclose(luvijkl,'KEEP') 3868 3869c Read HP or IP TERM 3870c luvijkl = -1 3871c call gpopen(luvijkl,'CCR12ITST','unknown',' ','unformatted', 3872c & idummy,.false.) 3873c read(luvijkl)(work(khintp-1+i),i=1,nt1amx) 3874c call gpclose(luvijkl,'KEEP') 3875 3876 if (testrbar) then 3877c Read R(t - 0.5 delta * C) into WORK(KRP) 3878 luvijkl = -1 3879 call gpopen(luvijkl,'RMTNTKL','unknown',' ','unformatted', 3880 & idummy,.false.) 3881 read(luvijkl)(work(krp-1+i),i=1,ngamsm(isymtr)) 3882 call gpclose(luvijkl,'KEEP') 3883 end if 3884 3885c Restore t amplitudes 3886 call daxpy(nt1amx,-0.5d0*delta,c1am,1,work(kt1am),1) 3887 3888c Calculate [V(t+0.5*delta*C)-V(t-0.5*delta*C)]/delta 3889 call daxpy(ngamsq(1),-1.0d00,work(kvintm),1,work(kvintp),1) 3890 call dscal(ngamsq(1),deltai,work(kvintp),1) 3891 3892c Calculate [RHO1(t+0.5*delta*C)-RHO1(t-0.5*delta*C)]/delta 3893 call daxpy(nt1amx,-1.0d00,work(krho1),1,work(krho1p),1) 3894 call dscal(nt1amx,deltai,work(krho1p),1) 3895 3896c Calculate [RHO2(t+0.5*delta*C)-RHO2(t-0.5*delta*C)]/delta 3897 call daxpy(nrho2,-1.0d00,work(krho2),1,work(krho2p),1) 3898 call dscal(nrho2,deltai,work(krho2p),1) 3899 3900c Calculate finite diff HP BAR 3901c call daxpy(nt1amx,-1.0d00,work(khintm),1,work(khintp),1) 3902c call dscal(nt1amx,deltai,work(khintp),1) 3903 3904c write(lupri,*)'IBAR numerically' 3905c do isymij = 1, nsym 3906c isymkl = muld2h(isymij,isymtr) 3907c call output(work(khintp+it1am(isymij,isymkl)),1, 3908c & nvir(isymij),1,nrhf(isymkl),nvir(isymij), 3909c & nrhf(isymkl),1,lupri) 3910c end do 3911 3912 3913c Print dOmega/dt_1 Result: 3914 write(lupri,*)'numerical: RHO1 and RHO2:' 3915 call cc_prp(work(krho1p),work(krho2p),1,1,1) 3916 write(lupri,*)'numerical: RHOR12:' 3917 do isymij = 1, nsym 3918 isymkl = muld2h(isymij,isymtr) 3919 call output(work(kvintp+igamsq(isymij,isymkl)),1, 3920 & nmatij(isymij),1,nmatij(isymkl),nmatij(isymij), 3921 & nmatij(isymkl),1,lupri) 3922 end do 3923 3924 if (testrbar) then 3925c Calculate [R(t+0.5*delta*C)-R(t-0.5*delta*C)]/delta 3926 call daxpy(ngamsm(1),-1.0d00,work(krm),1,work(krp),1) 3927 call dscal(ngamsm(1),deltai,work(krp),1) 3928 end if 3929 3930 if (testvbar) then 3931c symmetrize V bar 3932 call symV(work(kvintp),isymtr,work(kvsym), 3933 & nrhf,imatij,igamsq,nmatij,work(kend1),lwrk1) 3934 3935c Print the result: 3936 write(lupri,*)'norm^2 of V bar numerically:', 3937 & ddot(ngamsq(isymtr),work(kvsym),1,work(kvsym),1) 3938 do isymij = 1, nsym 3939 isymkl = muld2h(isymij,isymtr) 3940 write(lupri,*) 'isymij,isymkl:',isymij,isymkl 3941 write(lupri,*)'CC2-R12 <V bar numerically>' 3942 call output(work(kvsym+igamsq(isymij,isymkl)),1, 3943 & nmatij(isymij),1,nmatij(isymkl),nmatij(isymij), 3944 & nmatij(isymkl),1,lupri) 3945 end do 3946c split into triplett and singulett 3947 nrhftria = nrhftb*(nrhftb+1)/2 3948 kvsing = kend1 3949 kvtrip = kvsing + nrhftria*nrhftria 3950 kend1 = kvtrip + nrhftria*nrhftria 3951 write(lupri,*)'Vsing and Vtrip out of fdvint' 3952 call cc_r12vunpack(work(kvsym),isymtr,work(kvsing), 3953 & work(kvtrip),.true.,nrhfb,nrhf) 3954 end if 3955 3956 if (testrbar) then 3957c Print result for R: 3958c symmetrize R: 3959 call symV(work(krp),isymtr,work(krsym), 3960 & nrhfsa,imatijm,igamsm,nmatijm,work(kend1),lwrk1) 3961 write(lupri,*)'norm^2 of R bar numerically:', 3962 & ddot(ngamsm(isymtr),work(krsym),1,work(krsym),1) 3963 do isymij = 1, nsym 3964 isymkl = muld2h(isymij,isymtr) 3965 write(lupri,*) 'isymij,isymkl:',isymij,isymkl 3966 write(lupri,*)'CC2-R12 <R bar numerically>' 3967 call output(work(krsym+igamsm(isymij,isymkl)),1, 3968 & nmatijm(isymij),1,nmatij(isymkl),nmatijm(isymij), 3969 & nmatij(isymkl),1,lupri) 3970 end do 3971 end if 3972 3973 omegor = omegorsv 3974 3975 write(lupri,*)'leaving cc_r12fdvint' 3976 call qexit('r12fdvint') 3977 end 3978*=====================================================================* 3979 subroutine symV(vijkl,isymv,vsym,ndim1,ioff1,ioff2,ndim2, 3980 & work,lwork) 3981c--------------------------------------------------------------------- 3982c purpose: symmetrize Vijkl = 1/2 (Vijkl + Vjilk) 3983c 3984c H. Fliegl, C. Haettig, winter 2004 3985c 3986c ndim1 --> nrhf / nrhfs 3987c ioff1 --> imatij / imatijm 3988c ioff2 --> igamsq / igamsm 3989c ndim2 --> nmatij / nmatijm 3990c 3991c--------------------------------------------------------------------- 3992 implicit none 3993#include "priunit.h" 3994#include "ccsdsym.h" 3995#include "ccorb.h" 3996 3997 integer lwork,isymv,lwrk1,kend1,isymkl,isymij,isymk, 3998 & isyml,isymi,isymj,idxkl,idxlk,idxij,idxji,idxijkl,idxjilk 3999 integer ndim1(8),ioff1(8,8),ioff2(8,8),ndim2(8) 4000 4001#if defined (SYS_CRAY) 4002 real vijkl(*),work(*),half,vsym(*) 4003#else 4004 double precision vijkl(*),work(*),half,vsym(*) 4005#endif 4006 parameter ( half = 0.5d0 ) 4007 4008 call qenter('symV') 4009 4010 do isymkl = 1, nsym 4011 isymij = muld2h(isymv,isymkl) 4012 do isymk = 1, nsym 4013 isyml = muld2h(isymkl,isymk) 4014 do k = 1, nrhfb(isymk) 4015 do l = 1, nrhfb(isyml) 4016 idxkl = imatkl(isymk,isyml)+nrhfb(isymk)*(l-1)+k 4017 idxlk = imatkl(isyml,isymk)+nrhfb(isyml)*(k-1)+l 4018 do isymi = 1, nsym 4019 isymj = muld2h(isymij,isymi) 4020 do i = 1, ndim1(isymi) 4021 do j = 1, ndim1(isymj) 4022 idxij = ioff1(isymi,isymj)+ndim1(isymi)*(j-1)+i 4023 idxji = ioff1(isymj,isymi)+ndim1(isymj)*(i-1)+j 4024 idxijkl = ioff2(isymij,isymkl)+ 4025 & ndim2(isymij)*(idxkl-1)+idxij 4026 idxjilk = ioff2(isymij,isymkl)+ 4027 & ndim2(isymij)*(idxlk-1)+idxji 4028 vsym(idxijkl)=half*(vijkl(idxijkl)+vijkl(idxjilk)) 4029 end do 4030 end do 4031 end do 4032 end do 4033 end do 4034 end do 4035 end do 4036 4037 call qexit('symV') 4038 end 4039*=====================================================================* 4040 subroutine cc_r12stabmat(amat,maxred,len1,len2,len3,work,lwork) 4041c---------------------------------------------------------------------- 4042c purpose: get stability matrix and calculate the corresponding 4043c eigenvalues 4044c 4045c e C^T 4046c C B 4047c 4048c H. Fliegl, C. Haettig 4049c---------------------------------------------------------------------- 4050 implicit none 4051#include "priunit.h" 4052 4053 integer maxred,len1,len2,len3,lwork,i,j,idxi,idxj,ierr 4054 integer kevr,kevi,kevec,kiv1,kfv1,kend1,lwrk1,matz,nm 4055#if defined (SYS_CRAY) 4056 real amat(maxred,maxred),work(*),stabmat(len2+len3,len2+len3) 4057 real zero,stabsave(len2+len3,len2+len3) 4058#else 4059 double precision amat(maxred,maxred),work(*), 4060 & stabmat(len2+len3,len2+len3),zero 4061 double precision stabsave(len2+len3,len2+len3) 4062#endif 4063 parameter(zero = 0.0d0) 4064 4065 call qenter('stabmat') 4066c 4067 nm = len2+len3 4068 4069 kevr = 1 4070 kevi = kevr + nm 4071 kevec = kevi + nm 4072 kiv1 = kevec + nm*nm 4073 kfv1 = kiv1 + nm*nm 4074 kend1 = kfv1 + nm*nm 4075 lwrk1 = lwork - kend1 4076 if (lwrk1.lt.0) then 4077 call quit('Insufficient work space in stabmat') 4078 end if 4079 4080c get stability matrix out of jacobian 4081 do i = len1+1, len1+len2+len3 4082 idxi = i-len1 4083 do j = len1+1, len1+len2+len3 4084 idxj = j-len1 4085 stabmat(idxi,idxj) = amat(i,j) 4086 end do 4087 end do 4088 4089 write(lupri,*)'STABILITY-MATRIX' 4090 call output(stabmat,1,len2+len3,1,len2+len3,len2+len3, 4091 & len2+len3,1,lupri) 4092c save matrix 4093 call dcopy(nm*nm,stabmat,1,stabsave,1) 4094 4095c get eigenvalues of stability matrix 4096 matz = 1 ! get eigenvalues and eigenvectors 4097 call rg(nm,nm,stabmat,work(kevr),work(kevi),matz,work(kevec), 4098 & work(kiv1),work(kfv1),ierr) 4099c 4100 do i = 1, nm 4101 if (work(kevr-1+i).lt.zero) then 4102 write(lupri,*) 4103 & 'WARNING: NEGATIVE EIGENVALUES OF STABILITY-MATRIX' 4104 end if 4105 end do 4106c 4107 write(lupri,*)'Eigenvalues of STABILITY-MATRIX' 4108 call output(work(kevr),1,len2+len3,1,1,len2+len3, 4109 & len2+len3,1,lupri) 4110c 4111 call dcopy(nm*nm,stabsave,1,stabmat,1) 4112c set nondiagonal blocks to zero and recalculate the eigenvalues 4113 do i = len2+1, nm 4114 do j = 1, len2 4115 stabmat(i,j) = zero 4116 stabmat(j,i) = zero 4117 end do 4118 end do 4119 write(lupri,*)'STABILITY-MATRIX without C-MAT' 4120 call output(stabmat,1,len2+len3,1,len2+len3,len2+len3, 4121 & len2+len3,1,lupri) 4122 4123 call rg(nm,nm,stabmat,work(kevr),work(kevi),matz,work(kevec), 4124 & work(kiv1),work(kfv1),ierr) 4125c 4126 do i = 1, nm 4127 if (work(kevr-1+i).lt.zero) then 4128 write(lupri,*) 4129 & 'WARNING: NEGATIVE EIGENVALUES OF STABILITY-MATRIX' 4130 end if 4131 end do 4132c 4133 write(lupri,*)'Eigenvalues of STABILITY-MATRIX without C-MAT' 4134 call output(work(kevr),1,len2+len3,1,1,len2+len3, 4135 & len2+len3,1,lupri) 4136c 4137 4138 call qexit('stabmat') 4139 4140 end 4141*=====================================================================* 4142*=====================================================================* 4143 subroutine cc_r12getramnp(xpmqn,framnp,oneaux,norb1,norb2,nrhf1, 4144 & nrhfs1,ih1am,ih2am,work,lwork) 4145c--------------------------------------------------------------------- 4146c purpose: get r_a,mn^q' out of r^mn_pq 4147c and put on file 4148c 4149c C. Neiss, spring 2006 4150c--------------------------------------------------------------------- 4151 implicit none 4152#include "priunit.h" 4153#include "ccorb.h" 4154 logical locdbg,oneaux 4155 parameter (locdbg = .false.) 4156 integer ih2am(8,8), ih1am(8,8), nvircc(8),ldr12(8), 4157 & norb1(8),norb2(8),nrhf1(8),nrhfs1(8), 4158 & nmatkl(8),imatkl(8,8),nramn(8),iramn(8,8), 4159 & nramnq(8),iramnq(8,8) 4160 integer isym,isym1,isym2,isymp,isymq,isympmn,isymmn,isymm,isymn, 4161 & isympm,isymqn 4162 integer lwork,kend1,lwrk1,kramn,lunit,len,iadr,i,j,a,p,q,m,n 4163 integer nmn,npm,nqn,idxpmqn,idxamn,index 4164 character*(*) framnp 4165#if defined (SYS_CRAY) 4166 real xpmqn(*),work(*),ddot 4167#else 4168 double precision xpmqn(*),work(*),ddot 4169#endif 4170 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 4171 4172 call qenter('cc_r12getramnp') 4173 4174c calculate some offsets and dimensions 4175c calculate number of active virtual orbitals in CC 4176c nvirfr = frozen virtual orbitals 4177 do isym = 1, nsym 4178 nvircc(isym) = norb1(isym) - nrhfsa(isym) - nvirfr(isym) 4179 end do 4180 4181 if (locdbg) then 4182 write(lupri,*) 'Entered cc_r12getramnp' 4183 write(lupri,*) 'isym nrhfs1 nrhf1 nvircc norb1' 4184 do isym = 1, nsym 4185 write(lupri,*) isym,nrhfs1(isym),nrhf1(isym), 4186 & nvircc(isym),norb1(isym) 4187 end do 4188 end if 4189 4190 do isym = 1, nsym 4191 nmatkl(isym) = 0 4192 do isym2 = 1, nsym 4193 isym1 = muld2h(isym,isym2) 4194 imatkl(isym1,isym2) = nmatkl(isym) 4195 nmatkl(isym) = nmatkl(isym) + nrhf1(isym1)*nrhf1(isym2) 4196 end do 4197 end do 4198 4199 do isym = 1, nsym 4200 nramn(isym) = 0 4201 do isym2 = 1, nsym 4202 isym1 = muld2h(isym,isym2) 4203 iramn(isym1,isym2) = nramn(isym) 4204 nramn(isym) = nramn(isym) + nvircc(isym1)*nmatkl(isym2) 4205 end do 4206 end do 4207 4208 do isym = 1, nsym 4209 nramnq(isym) = 0 4210 do isym2 = 1, nsym 4211 isym1 = muld2h(isym,isym2) 4212 iramnq(isym1,isym2) = nramnq(isym) 4213 nramnq(isym) = nramnq(isym) + nramn(isym1)*norb2(isym2) 4214 end do 4215 end do 4216 4217 do isym = 1, nsym 4218 if (oneaux) then 4219 ldr12(isym) = norb1(isym) 4220 else 4221 ldr12(isym) = nvir(isym) 4222 end if 4223 end do 4224 4225c 4226c --------------------------------------------------------------- 4227c open output file 4228c --------------------------------------------------------------- 4229 lunit = -1 4230 call wopen2(lunit,framnp,64,0) 4231c 4232c --------------------------------------------------------------- 4233c reorder integrals 4234c --------------------------------------------------------------- 4235 do isymq = 1, nsym 4236 isympmn = isymq 4237 do q = norb1(isymq)+1, norb1(isymq)+norb2(isymq) 4238 len = nramn(isympmn) 4239 kramn = 1 4240 kend1 = kramn + len 4241 lwrk1 = lwork - kend1 4242 if (lwrk1.lt.0) then 4243 call quit('Insufficient memory in CC_R12GETRAMNP') 4244 end if 4245 do isymp = 1, nsym 4246 isymmn = muld2h(isympmn,isymp) 4247c 4248 do isymm = 1, nsym 4249 isymn = muld2h(isymmn,isymm) 4250 isympm = muld2h(isymp,isymm) 4251 isymqn = muld2h(isymq,isymn) 4252 do n = 1, nrhf1(isymn) 4253 do m = 1, nrhf1(isymm) 4254 do p = nrhfs1(isymp)+1, nrhfs1(isymp)+nvircc(isymp) 4255 a = p - nrhfs1(isymp) 4256 npm = ih1am(isymp,isymm)+ldr12(isymp)*(m-1)+p 4257 nqn = ih1am(isymq,isymn)+ldr12(isymq)*(n-1)+q 4258 nmn = imatkl(isymm,isymn)+nrhf1(isymm)*(n-1)+m 4259 idxpmqn = ih2am(isympm,isymqn) + index(npm,nqn) 4260 idxamn = iramn(isymp,isymmn) + 4261 & nvircc(isymp)*(nmn-1)+a 4262 work(kramn+idxamn-1) = xpmqn(idxpmqn) 4263 end do 4264 end do 4265 end do 4266 end do 4267 end do 4268 iadr = iramnq(isympmn,isymq) + 4269 & nramn(isympmn)*(q-norb1(isymq)-1) + 1 4270 call putwa2(lunit,framnp,work(kramn),iadr,len) 4271 if (locdbg) then 4272 write(lupri,*) 'isymq, q, iadr= ',isymq,q,iadr 4273 write(lupri,*) 'Norm^2:', ddot(len,work(kramn),1, 4274 & work(kramn),1) 4275 do isym2 = 1, nsym 4276 isym1 = muld2h(isympmn,isym2) 4277 write(lupri,*) 'Symmetry block: ',isym1,isym2 4278 call output(work(kramn+iramn(isym1,isym2)), 4279 & 1,nvircc(isym1),1,nmatkl(isym2), 4280 & nvircc(isym1),nmatkl(isym2),1,lupri) 4281 end do 4282 end if 4283c 4284 end do 4285 end do 4286c 4287c --------------------------------------------------------------- 4288c close file 4289c --------------------------------------------------------------- 4290 call wclose2(lunit,framnp,'KEEP') 4291c 4292 if (locdbg) then 4293 write(lupri,*) 'Leaving cc_r12getramnp' 4294 end if 4295c 4296 call qexit('cc_r12getramnp') 4297 return 4298 end 4299*=====================================================================* 4300