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_r12prepccsd(work,lwork) 20c-------------------------------------------------------------------- 21c purpose: prepare V intermediates for the CCSD(R12) model 22c using ansatz 1 that do not depend on cluster ampl. 23c 24c H. Fliegl, C. Haettig summer 2004 25c-------------------------------------------------------------------- 26 implicit none 27#include "priunit.h" 28#include "ccorb.h" 29#include "ccsdsym.h" 30#include "r12int.h" 31#include "ccr12int.h" 32#include "dummy.h" 33#include "ccsdinp.h" 34 35 logical locdbg,lvajkl,lvabkl,lv,lvijkl,lexist 36 parameter (locdbg = .false.) 37 38 integer lwork,ibasx(8),isym,kend0,lwrk0 39 integer klamdp,klamdh,kend1,kt1am,kend2,lwrk2 40 integer isymtr,kvabkl,kend3,lwrk3,kvajkl, 41 & isymc,isymab,isymkl,lunit,idum 42 integer ioptbas,iopt 43 44#if defined (SYS_CRAY) 45 real zero,one,work(*),ddummy,timrdao,timfr12,timintr12,ddot 46#else 47 double precision zero,one,work(*),ddummy,timrdao,timfr12,timintr12 48 double precision ddot 49#endif 50 parameter (zero = 0.0D0, one = 1.0D0) 51 52 call qenter('prepccsd') 53 if (locdbg) then 54 write(lupri,*)'entered cc_r12prepccsd' 55 end if 56c 57 kend0 = 1 58c 59c test if file fvabkl already exists, if yes exit: 60 inquire(file=fvabkl,exist=lexist) 61 if (lexist.and.ccrstr) then 62 write(lupri,*) 'Restart: Found V_(alpha beta)^(kl) on disk' 63 call qexit('prepccsd') 64 return 65 end if 66c 67c get CMO coefficients: 68c 69 klamdp = kend0 70 klamdh = klamdp + nlamdt 71 kend1 = klamdh + nlamdt 72 kt1am = kend1 73 kend2 = kt1am + nt1amx 74 lwrk2 = lwork - kend2 75 if(lwrk2.lt.0) then 76 call quit('insufficient work space in prepccsd') 77 end if 78 call dzero(work(kt1am),nt1amx) 79 call lammat(work(klamdp),work(klamdh),work(kt1am), 80 & work(kend2),lwrk2) 81 82c get V^kl_alpha_beta 83 kvabkl = kend2 84 kend2 = kvabkl + nvabkl(1) 85 lwrk2 = lwork - kend2 86 if(lwrk2.lt.0) then 87 call quit('insufficient work space in prepccsd') 88 end if 89 ! initialize V_(alpha beta)^(kl) 90 iopt = 0 91 call cc_r12mkvamkl0(work(kvabkl),nvabkl(1),iopt,ddummy,idummy, 92 & work(kend2),lwrk2) 93 94 if (locdbg) then 95 write(lupri,*)'norm^2 Sak,bl =', 96 & ddot(nvabkl(1),work(kvabkl),1,work(kvabkl),1) 97 write(lupri,*)'Sak*Sbl' 98 do isymab = 1, nsym 99 isymkl = isymab 100 call output(work(kvabkl+ivabkl(isymab,isymkl)),1, 101 & n2bst(isymab),1,nmatkl(isymkl),n2bst(isymab), 102 & nmatkl(isymkl),1,lupri) 103 end do 104 end if 105 106 if (locdbg) then 107c test if V is ok, transform all indices to occupied 108 isymc = 1 109 call cc_r12vtest(work(kvabkl),work(klamdp),isymc, 110 & work(kend2),lwrk2) 111 end if 112c 113 isymtr = 1 114 lvabkl = .true. 115 ioptbas = 2 116 call cc_mofconr12(work(klamdh),1,ddummy,ddummy,ddummy, 117 & isymtr,ddummy,ddummy,ddummy,work(kvabkl), 118 & .false.,.false.,lvabkl,ioptbas, 119 & timrdao,timfr12,timintr12, 120 & idummy,idummy,idummy,idummy,idummy,idummy, 121 & work(kend2),lwrk2) 122 123 if (locdbg) then 124 write(lupri,*)'norm^2 Vabkl: ', 125 & ddot(nvabkl(1),work(kvabkl),1,work(kvabkl),1) 126 write(lupri,*)'V^ab_kl in prepccsdr12' 127 do isymab = 1, nsym 128 isymkl = isymab 129 call output(work(kvabkl+ivabkl(isymab,isymkl)), 130 & 1,n2bst(isymab),1,nmatkl(isymkl), 131 & n2bst(isymab),nmatkl(isymkl),1,lupri) 132 end do 133 end if 134 135c symmetrize V^ kl_alpha_beta and write it on file 136 call cc_r12_symv(work(kvabkl),work(kend2),lwrk2) 137 lunit = -1 138 call gpopen(lunit,fvabkl,'unknown',' ','unformatted',idum,.false.) 139 write(lunit)(work(kvabkl-1+i),i=1,nvabkl(1)) 140 call gpclose(lunit,'KEEP') 141 142 if (locdbg) then 143 write(lupri,*)'norm^2 Vabkl sym', 144 & ddot(nvabkl(1),work(kvabkl),1,work(kvabkl),1) 145 write(lupri,*)'symmetrized V^ab_kl in prepccsdr12' 146 do isymab = 1, nsym 147 isymkl = isymab 148 call output(work(kvabkl+ivabkl(isymab,isymkl)), 149 & 1,n2bst(isymab),1,nmatkl(isymkl), 150 & n2bst(isymab),nmatkl(isymkl),1,lupri) 151 end do 152c test if V is ok, transform all indices to occupied 153 isymc = 1 154 call cc_r12vtest(work(kvabkl),work(klamdp),isymc, 155 & work(kend2),lwrk2) 156 end if 157 158c calculate V^kl_alphaj and write it on file 159 isymc = 1 160 kvajkl = kend2 161 kend3 = kvajkl + nvajkl(isymc) 162 lwrk3 = lwork - kend3 163 if(lwrk3.lt.0) then 164 call quit('insufficient work space in prepccsd') 165 end if 166 lv = .false. 167 lvajkl = .true. 168 lvijkl = .false. 169 call cc_r12mkvtf(work(kvabkl),work(kvajkl),dummy, 170 & work(klamdp),isymc,lv,lvijkl,lvajkl,fvajkl, 171 & work(kend3),lwrk3) 172 173c calculate V^ab_kl = \sum_albe L^h_ala L^h_beb V^albe_kl and 174c write it on file 175 isymc = 1 176 iopt = 1 177 call cc_r12mkvirt(work(kvabkl),work(klamdh),isymc, 178 & work(klamdh),isymc,fvcdkl,iopt, 179 & work(kend2),lwrk2) 180 181 call qexit('prepccsd') 182 return 183 end 184*====================================================================* 185 subroutine cc_r12mkvtf(vabkl,vajkl,vijkl,cmo,isymc, 186 & lv,lvijkl,lvajkl,filvajkl,work,lwork) 187c-------------------------------------------------------------------- 188c purpose: transform V^alpha_beta _kl to V^ij_kl to test V with 189c MP2-R12 190c 191c lv read V^(alpha beta)_kl from file 192c lvijkl return vijkl (result is ADDED to vijkl) 193c lvajkl write vajkl to file filvajkl 194c 195c Note: when lvijkl and lvajkl are both .FALSE. only vajkl 196c will be returned (in memory) 197c 198c 199c H. Fliegl, C. Haettig, summer 2004 200c modified C. Neiss, autumn 2005 201c-------------------------------------------------------------------- 202 implicit none 203#include "priunit.h" 204#include "ccorb.h" 205#include "ccsdsym.h" 206#include "r12int.h" 207#include "ccr12int.h" 208 209 logical lres,locdbg,lv,lvajkl,lvijkl 210 parameter (locdbg = .false.) 211 integer isymc,lwork,isymkl,isymab,isymk, 212 & isyml,isyma,isymb,isymaj,isymj, 213 & koffc,ntota,ntotb,idxkl, 214 & kvabkl,kvajkl 215 integer luvajkl,lunit,idummy 216 character*(*) filvajkl 217#if defined (SYS_CRAY) 218 real vabkl(*),vajkl(*),vijkl(*),cmo(*),work(*),one,zero,ddot 219#else 220 double precision vabkl(*),vajkl(*),vijkl(*),cmo(*),work(*),one, 221 & zero,ddot 222#endif 223 parameter (one = 1.0d0, zero = 0.0d0) 224 225 call qenter('r12mkvtf') 226 227 if (lv) then 228c read in Vabkl 229 lunit = -1 230 call gpopen(lunit,fvabkl,'unknown',' ','unformatted', 231 & idummy,.false.) 232 read (lunit)(vabkl(i),i=1,nvabkl(1)) 233 call gpclose(lunit,'KEEP') 234 end if 235 236 call dzero(vajkl,nvajkl(isymc)) 237 do isymkl = 1, nsym 238 isymab = isymkl 239 do isyma = 1, nsym 240 isymb = muld2h(isymab,isyma) 241 isymj = muld2h(isymc,isymb) 242 isymaj = muld2h(isyma,isymj) 243 do isymk = 1, nsym 244 isyml = muld2h(isymkl,isymk) 245 do k = 1, nrhfb(isymk) 246 do l = 1, nrhfb(isyml) 247 idxkl = imatkl(isymk,isyml)+nrhfb(isymk)*(l-1)+k 248 kvabkl = ivabkl(isymab,isymkl)+n2bst(isymab)*(idxkl-1)+ 249 & iaodis(isyma,isymb)+1 250 kvajkl = ivajkl(isymaj,isymkl)+nt1ao(isymaj)*(idxkl-1)+ 251 & it1ao(isyma,isymj)+1 252 253 koffc = iglmrh(isymb,isymj)+1 254 ntotb = max(mbas1(isymb),1) 255 ntota = max(mbas1(isyma),1) 256 call dgemm('N','N',mbas1(isyma),nrhf(isymj), 257 & mbas1(isymb),one,vabkl(kvabkl),ntota, 258 & cmo(koffc),ntotb,one,vajkl(kvajkl),ntota) 259 end do 260 end do 261 end do 262 end do 263 end do 264 265 if (locdbg) then 266 write(lupri,*)'norm^2 vajkl', 267 & ddot(nvajkl(isymc),vajkl,1,vajkl,1) 268c write(lupri,*)'vajkl' 269c do isymkl = 1, nsym 270c isymaj = muld2h(isymc,isymkl) 271c call output(vajkl(1+ivajkl(isymaj,isymkl)),1, 272c & nt1ao(isymaj),1,nmatkl(isymkl), 273c & nt1ao(isymaj),nmatkl(isymkl),1,lupri) 274c end do 275 end if 276c 277 if (lvajkl) then 278c ----------------------------- 279c write V(alpha j,kl) on file 280c ----------------------------- 281 luvajkl = -1 282 call gpopen(luvajkl,filvajkl,'unknown',' ','unformatted', 283 & idummy,.false.) 284 rewind(luvajkl) 285 write(luvajkl) (vajkl(i), i = 1,nvajkl(isymc)) 286 call gpclose(luvajkl,'KEEP') 287 end if 288c 289 if (lvijkl) then 290 lres = .true. 291 call cc_r12mkvijkl(vajkl,isymc,cmo,isymc,work,lwork, 292 & lres,one,vijkl) 293 294 if (locdbg) then 295 write(lupri,*)'in mkvtf: norm^2 vijkl:', 296 & ddot(ntr12sq(1),vijkl,1,vijkl,1) 297 end if 298 end if 299c 300 call qexit('r12mkvtf') 301 end 302*====================================================================* 303 subroutine cc_r12mkvabkl(vabkl,xint,idel,isymd,isydis,ibastyp, 304 & ibasx,work,lwork) 305c-------------------------------------------------------------------- 306c purpose: make V^kl_alpha_beta for CCSD(12) model with ansatz 1 307c 308c H. Fliegl, C. Haettig, summer 2004 309c 310c C. Neiss, 05.12.2005: adapted for CABS 311c-------------------------------------------------------------------- 312 implicit none 313#include "priunit.h" 314#include "maxorb.h" 315#include "ccorb.h" 316#include "ccisao.h" 317#include "ccsdsym.h" 318#include "ccsdinp.h" 319#include "r12int.h" 320#include "ccr12int.h" 321 322 logical locdbg,lauxd 323 parameter(locdbg =.false.) 324 325 integer nr1orb(8),nr1bas(8),nr1xbas(8),nr2bas,n2bst1(8) 326 integer ir1xbas(8,8),ir1orb(8,8),ir1bas(8,8),ir2bas(8,8), 327 & iaodis1(8,8),ir2xbas(8,8),irgkl(8,8),nrgkl(8), 328 & nalphaj(8),ialphaj(8,8) 329 integer nr1xorb(8),ir1xorb(8,8),nrxgkl(8),irxgkl(8,8) 330 integer idel,isymd,isydis,ibastyp,lwork,lwrk1,kend1,krgkl, 331 & ibasx(8),isym,isymg,isymab,isymkl,koffr,koffg,koffv, 332 & ntotg,ntotab,isym1,isym2,icount1,ngab(8),igab(8,8) 333 integer kscr1,kscr2,isymb,isymga,isymgab,koff1, 334 & idxab,idxga,idxgab,isyma,kend2,lwrk2,isymgkl 335 336#if defined (SYS_CRAY) 337 real work(*),xint(*),vabkl(*),one,two,factor,ddot 338#else 339 double precision work(*),xint(*),vabkl(*),one,two,factor, 340 & ddot 341#endif 342 parameter (one = 1.0d0, two = 2.0d0) 343 344 call qenter('mkvabkl') 345 346 call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 347 & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas, 348 & ir2bas,ir2xbas,irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 349c 350 do isym = 1, nsym 351 ngab(isym) = 0 352 icount1 = 0 353 do isym2 = 1, nsym 354 isym1 = muld2h(isym,isym2) 355 ngab(isym) = ngab(isym) + mbas1(isym1)*n2bst(isym2) 356 igab(isym1,isym2) = icount1 357 icount1 = icount1 + mbas1(isym1)*n2bst(isym2) 358 end do 359 end do 360c 361 krgkl = 1 362 kend1 = krgkl + nrgkl(isymd) 363 lwrk1 = lwork - kend1 364 if (lwrk1.lt.0) then 365 call quit('insufficient work space in mkvabkl') 366 end if 367 368 if (locdbg) then 369 write(lupri,*) 'ibastyp,isymd,idel:',ibastyp,isymd,idel 370 write(lupri,*) 'ibasx:',ibasx 371 end if 372 373c get r12 integrals 374 if (ibastyp.eq.1) then 375 lauxd = .false. 376 call cc_r12getrint(work(krgkl),idel,isymd,nr1bas,ir1bas, 377 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 378 & nrhfb,nmatkl,imatkl,ibasx,lauxd,.false., 379 & fnback,work(kend1),lwrk1) 380 factor = one 381 if (r12cbs) factor = -one 382 else 383 lauxd = .true. 384 call cc_r12getrint(work(krgkl),idel,isymd,nr1bas,ir1bas, 385 & nr2bas,ir2bas,nrgkl,irgkl,ir1xbas,ir2xbas, 386 & nrhfb,nmatkl,imatkl,ibasx,lauxd,.false., 387 & fnback,work(kend1),lwrk1) 388 factor = - two 389 end if 390c 391 isymgab = isydis 392c 393 kscr2 = kend1 394 kend1 = kscr2 + ngab(isymgab) 395 lwrk1 = lwork - kend1 396 if (lwrk1.lt.0) then 397 call quit('insufficient work space in mkvabkl') 398 end if 399c 400 do isymb = 1, nsym 401 isymga = muld2h(isydis,isymb) 402 403c dynamic work space allocation 404 kscr1 = kend1 405 kend2 = kscr1 + n2bst(isymga) 406 lwrk2 = lwork - kend2 407 if (lwrk2.lt.0) then 408 call quit('insufficient work space in mkvabkl') 409 end if 410c 411 do b = 1, mbas1(isymb) 412c pack triangular matrix to quadratic matrix 413 koff1 = idsaog(isymb,isydis) + 1 + nnbst(isymga)*(b-1) 414 call ccsd_symsq(xint(koff1),isymga,work(kscr1)) 415c 416 do isyma = 1, nsym 417 isymab = muld2h(isymb,isyma) 418 isymg = muld2h(isymga,isyma) 419 do a = 1, mbas1(isyma) 420 idxab = iaodis(isyma,isymb)+mbas1(isyma)*(b-1)+a 421 do g = 1, mbas1(isymg) 422 idxga = iaodis(isymg,isyma)+mbas1(isymg)*(a-1)+g 423 idxgab = igab(isymg,isymab)+mbas1(isymg)*(idxab-1)+g 424 work(kscr2-1+idxgab) = work(kscr1-1+idxga) 425 end do 426 end do 427 end do 428c 429 end do !b 430 end do ! isymb 431 432 if (locdbg) then 433 write(lupri,*)'norm^2 xint:', 434 & ddot(ndisao(isydis),xint,1,xint,1) 435 write(lupri,*)'norm^2 work(kscr2): ', 436 & ddot(ngab(isymgab),work(kscr2),1,work(kscr2),1) 437 write(lupri,*)'norm^2 work(krgkl)', 438 & ddot(nrgkl(isymd),work(krgkl),1,work(krgkl),1) 439 end if 440 441 do isymg = 1, nsym 442 isymab = muld2h(isymgab,isymg) 443 isymkl = isymab 444 isymgkl = muld2h(isymg,isymkl) 445 446 koffr = krgkl + irgkl(isymg,isymkl) 447 koffg = kscr2 + igab(isymg,isymab) 448 koffv = ivabkl(isymab,isymkl)+1 449 450 ntotab = max(n2bst(isymab),1) 451 ntotg = max(mbas1(isymg),1) 452 453 if (locdbg) then 454 write(lupri,*)'norm^2 work(koffg)', 455 & ddot(mbas1(isymg)*n2bst(isymab),work(koffg),1,work(koffg),1) 456 write(lupri,*)'norm^2 work(koffr)', 457 & ddot(mbas1(isymg)*nmatkl(isymkl),work(koffr),1,work(koffr),1) 458 end if 459 460 call dgemm('T','N',n2bst(isymab),nmatkl(isymkl),mbas1(isymg), 461 & factor,work(koffg),ntotg,work(koffr),ntotg,one, 462 & vabkl(koffv),ntotab) 463 end do 464 465 call qexit('mkvabkl') 466 end 467*====================================================================* 468 subroutine cc_r12vtest(vabkl,cmo,isymc,work,lwork) 469c-------------------------------------------------------------------- 470c purpose: test if V^ab_kl is ok 471c 472c H. Fliegl, C. Haettig, summer 2004 473c-------------------------------------------------------------------- 474 implicit none 475#include "priunit.h" 476#include "ccorb.h" 477#include "ccsdsym.h" 478#include "r12int.h" 479#include "ccr12int.h" 480 481 logical lv,lvajkl,lvijkl 482 integer kvijkl,kvajkl,isymc, 483 & kend1,lwrk1,lwork 484 485#if defined (SYS_CRAY) 486 real vabkl(*),work(*),cmo(*),ddot 487#else 488 double precision vabkl(*),work(*),cmo(*),ddot 489#endif 490 call qenter('r12vtest') 491 492 kvijkl = 1 493 kvajkl = kvijkl + nrhftb*nrhftb*nrhftb*nrhftb 494 kend1 = kvajkl + nvajkl(isymc) 495 lwrk1 = lwork - kend1 496 if (lwrk1.lt.0) then 497 call quit('insufficient work space in r12vtest') 498 end if 499 call dzero(work(kvijkl),nrhftb*nrhftb*nrhftb*nrhftb) 500 lv = .false. 501 lvajkl = .false. 502 lvijkl = .true. 503 call cc_r12mkvtf(vabkl,work(kvajkl),work(kvijkl), 504 & cmo,isymc,lv,lvijkl,lvajkl,fvajkl, 505 & work(kend1),lwrk1) 506 write(lupri,*)'Norm^2 Vijkl in r12vtest: ', 507 & ddot(nrhftb*nrhftb*nrhftb*nrhftb,work(kvijkl),1,work(kvijkl),1) 508 call output(work(kvijkl),1,nrhftb*nrhftb,1,nrhftb*nrhftb, 509 & nrhftb*nrhftb,nrhftb*nrhftb,1,lupri) 510 write(lupri,*)'Vijkl in r12vtest in triangular format:' 511 call ccr12pck2(work(kend1),isymc,.FALSE.,work(kvijkl),'T',1) 512 call cc_prpr12(work(kend1),isymc,1,.FALSE.) 513 514 call qexit('r12vtest') 515 end 516*====================================================================* 517 subroutine cc_r12_symv(vabkl,work,lwork) 518c-------------------------------------------------------------------- 519c purpose: symmetrize V^kl_ab 520c P^ab_kl * V^kl_ab = 1/2 * (V^kl_ab + V^lk_ba) 521c 522c H. Fliegl, C. Haettig, summer 2004 523c-------------------------------------------------------------------- 524 implicit none 525#include "priunit.h" 526#include "ccorb.h" 527#include "ccsdsym.h" 528#include "r12int.h" 529 530 integer lwork,kvsym,kend1,lwrk1,isyma,isymb, 531 & isymab,isymkl,isyml,isymk,idxkl,idxlk,idxab,idxba,idxabkl, 532 & idxbalk 533#if defined (SYS_CRAY) 534 real vabkl(*),work(*),half 535#else 536 double precision vabkl(*),work(*),half 537#endif 538 parameter (half = 0.5d0) 539 call qenter('r12symv') 540 541 kvsym = 1 542 kend1 = kvsym + nvabkl(1) 543 lwrk1 = lwork - kend1 544 if (lwrk1.lt.0) then 545 call quit('insufficient work space in r12symv') 546 end if 547 548 do isyma = 1, nsym 549 do isymb = 1, nsym 550 isymab = muld2h(isyma,isymb) 551 isymkl = isymab 552 do isymk = 1, nsym 553 isyml = muld2h(isymkl,isymk) 554 do k = 1, nrhfb(isymk) 555 do l = 1, nrhfb(isyml) 556 idxkl = imatkl(isymk,isyml)+nrhfb(isymk)*(l-1)+k 557 idxlk = imatkl(isyml,isymk)+nrhfb(isyml)*(k-1)+l 558 do a = 1, mbas1(isyma) 559 do b = 1, mbas1(isymb) 560 idxab = iaodis(isyma,isymb)+mbas1(isyma)*(b-1)+a 561 idxba = iaodis(isymb,isyma)+mbas1(isymb)*(a-1)+b 562 idxabkl = ivabkl(isymab,isymkl)+ 563 & n2bst(isymab)*(idxkl-1)+idxab 564 idxbalk = ivabkl(isymab,isymkl)+ 565 & n2bst(isymab)*(idxlk-1)+idxba 566 work(kvsym-1+idxabkl) = 567 & half*(vabkl(idxabkl)+vabkl(idxbalk)) 568 end do 569 end do 570 end do 571 end do 572 end do 573 end do 574 end do 575 576 call dcopy(nvabkl(1),work(kvsym),1,vabkl,1) 577 578 call qexit('r12symv') 579 return 580 end 581*====================================================================* 582 SUBROUTINE CCRHS_BP(OMEGA2,ISYMTR,IOPTOM,IAMP,FC12AM,LUFC12,IFILE, 583 & LISTR,IDLSTR,BASSCL2,WORK,LWORK) 584c-------------------------------------------------------------------- 585c purpose: calculate \sum_mn C^ij_mn * V^mn_ab 586c 587c IOPTOM = 0: dimension of OMEGA2 = 2*NT2ORT(ISYMTR) 588c = 1: dimension of OMEGA2 = NT2AO(ISYMTR) 589c 590c H. Fliegl, C. Haettig, summer 2004 591c modified C. Neiss, autumn 2005 592c-------------------------------------------------------------------- 593 implicit none 594#include "priunit.h" 595#include "ccorb.h" 596#include "ccsdsym.h" 597#include "dummy.h" 598#include "r12int.h" 599#include "ccr12int.h" 600#include "ccsdinp.h" 601 602 603 logical locdbg 604 parameter (locdbg = .false.) 605 integer lwork,kend1,lwrk1, 606 & isymv,idum,kvabkl,lunit,krpck,isymc,ntamp 607 integer isymab,iopt,iamp,lufc12,ifile,isymtr,isymij,idlstr,ioptom 608 character*10 model 609 character*8 fc12am 610 character*3 listr 611#if defined (SYS_CRAY) 612 real work(*),one,two,factor,zero,omega2(*), 613 & four,ddot,x,basscl2 614#else 615 double precision work(*),one,two, 616 & factor,zero,omega2(*),four,ddot,x,basscl2 617#endif 618 parameter (one = 1.0d0, two = 2.0d0, zero = 0.0d0, four = 4.0d0) 619 620 call qenter('ccrhs_bp') 621 622 if (locdbg.and.(ioptom.eq.0)) then 623 write(lupri,*)'Omega 2 entering ccrhs_bp' 624 do isymab = 1, nsym 625 isymij = muld2h(isymab,isymtr) 626 write(lupri,*) 'isymab, isymij: ', isymab,isymij 627 write(lupri,*) 'Omega+ =' 628 call output(omega2(it2ort(isymab,isymij)+1),1,nnbst(isymab),1, 629 & nmijp(isymij),nnbst(isymab),nmijp(isymij),1,lupri) 630 write(lupri,*)'Omega- =' 631 call output(omega2(it2ort(isymab,isymij)+1+nt2ort(isymtr)),1, 632 & nnbst(isymab),1,nmijp(isymij),nnbst(isymab), 633 & nmijp(isymij),1,lupri) 634 635 end do 636 end if 637 638C call cc_r12offset(nr1orb,nr1xorb,nr1bas,nr1xbas,nr2bas, 639C & nrgkl,nrxgkl,n2bst1,ir1orb,ir1xorb,ir1bas,ir1xbas,ir2bas,ir2xbas, 640C & irgkl,irxgkl,iaodis1,nalphaj,ialphaj) 641c 642 krpck = 1 643 kvabkl = krpck + ntr12am(isymtr) 644 kend1 = kvabkl + nvabkl(1) 645 lwrk1 = lwork - kend1 646 if (lwrk1 .lt.0) then 647 call quit('Insufficient work space in ccrhs_bp') 648 end if 649 650c read V 651 lunit = -1 652 call gpopen(lunit,fvabkl,'unknown',' ','unformatted',idum,.false.) 653 read(lunit)(work(kvabkl-1+i),i=1,nvabkl(1)) 654 call gpclose(lunit,'KEEP') 655 656c read r12 amplitudes 657 if (iamp.eq.0) then 658 iopt = 32 659 call cc_rdrsp('R0 ',0,1,iopt,model,dummy,work(krpck)) 660 else if (iamp.eq.1) then 661 call cc_rvec(lufc12,fc12am,ntr12am(isymtr),ntr12am(isymtr), 662 & ifile,work(krpck)) 663 call cclr_diasclr12(work(krpck),basscl2,isymtr) 664 else if (iamp.eq.2) then 665 iopt = 32 666 call cc_rdrsp(listr,idlstr,isymtr,iopt,model,dummy,work(krpck)) 667 call cclr_diasclr12(work(krpck),basscl2,isymtr) 668 else 669 call quit('Unknown value of IAMP in CCRHS_BP') 670 end if 671 672c calculate c*V 673 isymv = 1 674 isymc = isymtr 675 factor = one 676 call cc_r12mkcv(omega2,ioptom,work(kvabkl),work(krpck), 677 & isymv,isymc,factor,work(kend1),lwrk1) 678 679 if (locdbg) then 680 write(lupri,*)'in ccrhs_bp:' 681 write(lupri,*)'norm^2 work(kvabkl):', 682 & ddot(nvabkl(1),work(kvabkl),1,work(kvabkl),1) 683 write(lupri,*)'norm^2 work(krpck):', 684 & ddot(ntr12am(isymtr),work(krpck),1,work(krpck),1) 685 if (ioptom.eq.0) then 686 write(lupri,*)'norm^2 omega2', 687 & ddot(2*nt2ort(isymtr),omega2,1,omega2,1) 688 else if (ioptom.eq.1) then 689 write(lupri,*)'norm^2 omega2', 690 & ddot(nt2ao(isymtr),omega2,1,omega2,1) 691 end if 692 end if 693 694 if (locdbg.and.(ioptom.eq.0)) then 695 write(lupri,*)'Omega 2 leaving ccrhs_bp' 696 do isymab = 1, nsym 697 isymij = muld2h(isymab,isymtr) 698 write(lupri,*) 'isymab, isymij: ', isymab,isymij 699 write(lupri,*) 'Omega+ =' 700 call output(omega2(it2ort(isymab,isymij)+1),1,nnbst(isymab),1, 701 & nmijp(isymij),nnbst(isymab),nmijp(isymij),1,lupri) 702 write(lupri,*)'Omega- =' 703 call output(omega2(it2ort(isymab,isymij)+1+nt2ort(isymtr)),1, 704 & nnbst(isymab),1,nmijp(isymij),nnbst(isymab), 705 & nmijp(isymij),1,lupri) 706 end do 707 end if 708 709 call qexit('ccrhs_bp') 710 end 711*====================================================================* 712 subroutine cc_r12mkcv(omega2,ioptom,vabkl,c2am,isymv, 713 & isymc,factor,work,lwork) 714c-------------------------------------------------------------------- 715c purpose: calculate \sum_mn C^ij_mn * V^mn_ab 716c 717c H. Fliegl, C. Haettig, summer 2004 718c 719c modified C. Neiss, autumn 2005 720c 721c IOPTOM = 0: sort Omega2 in singlet/triplet format 722c (dim: 2*NT2ORT(ISYMOM)) 723c = 1: sort Omega2 in triangular matrix format 724c (dim: NT2AO(ISYMOM)) 725c-------------------------------------------------------------------- 726 implicit none 727#include "priunit.h" 728#include "maxorb.h" 729#include "ccorb.h" 730#include "ccsdsym.h" 731#include "symsq.h" 732#include "r12int.h" 733#include "ccr12int.h" 734 735 logical locdbg 736 parameter (locdbg = .false.) 737 738 integer lwork,isymv,isymc,isymom,isymb, 739 & isyma,isymab,isymkl,isymij,isymj,isymi,isymai,isymbj, 740 & kvabkl,kc2am,kom,kend1,lwrk1,ntoti,ntota,idxbj,idxai, 741 & koffai,maxai,naibj,index,kom2,kend2,lwrk2,idxij,idxab, 742 & idxabijp,idxabijm,idxbi,idxaj,idxaibj,idxbiaj,isymbi, 743 & isymaj,nab,kend3,lwrk3,maxj,maxb,ioptom 744 integer isymlj,idxljlj,idxlj,isymk,isyml,idxki,idxkilj 745#if defined (SYS_CRAY) 746 real work(*),vabkl(*),one,two,c2am(*),factor,zero 747 real half,omega2(*),ddot,x1,xv,xc,dummy 748#else 749 double precision work(*),vabkl(*),one,two,c2am(*), 750 & factor,zero,half,omega2(*),ddot, 751 & x1,xv,xc,dummy 752#endif 753 parameter (one = 1.0d0, two = 2.0d0, zero = 0.0d0, half = 0.5d0) 754 755 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 756 757 call qenter('r12mkcv') 758 if (.false.) then 759c small test if gamma term in ccrhs is exact V^ij_kl matrix 760c Caution!!!! d_ki, d_lj have to be put to one !!!! 761 call dzero(c2am,ntr12am(1)) 762 do isymk = 1, nsym 763 isymi = isymk 764 do isyml = 1, nsym 765 isymj = isyml 766 do i = 1, nrhf(isymi) 767 k = i 768 idxki = imatki(isymk,isymi)+nrhfb(isymk)*(i-1)+k 769 do j = 1, nrhf(isymj) 770 l = j 771 idxlj = imatki(isyml,isymj)+nrhfb(isyml)*(j-1)+l 772 idxkilj = itr12am(1,1)+index(idxki,idxlj) 773 c2am(idxkilj) = one 774 end do 775 end do 776 end do 777 end do 778 write(lupri,*)'c2am:' 779 call outpak(c2am,nmatij(1),1,lupri) 780 call dzero(omega2,2*nt2ort(1)) 781 end if 782 783 isymom = muld2h(isymv,isymc) 784 785 kend1 = 1 786 787 if (ioptom.eq.0) then 788 kom2 = 1 789 kend1 = kom2 + nt2ao(isymom) 790 lwrk1 = lwork - kend1 791 if (lwrk1.lt.0) then 792 call quit('insufficient work space in r12mkcv') 793 end if 794 call dzero(work(kom2),nt2ao(isymom)) 795 else if (ioptom.ne.1) then 796 call quit('Unknown value of "IOPTOM" in CC_R12MKCV') 797 end if 798 799 if (locdbg) then 800 x1 = 0.0d0 801 xv = 0.0d0 802 xc = 0.0d0 803 write(lupri,*)'norm^2(c2am):',ddot(ntr12am(isymc),c2am,1,c2am,1) 804 write(lupri,*)'isymom,nt2ao:',isymom,nt2ao(isymom) 805 if (ioptom.eq.0) then 806 write(lupri,*) 'norm^2 work(kom2)', 807 & ddot(nt2ao(isymom),work(kom2),1,work(kom2),1) 808 else if (ioptom.eq.1) then 809 write(lupri,*) 'norm^2 omega2', 810 & ddot(nt2ao(isymom),omega2,1,omega2,1) 811 end if 812 end if 813 814 do isymb = 1, nsym 815 do isyma = 1, nsym 816 isymab = muld2h(isyma,isymb) 817 isymkl = muld2h(isymv,isymab) 818 isymij = muld2h(isymkl,isymc) 819 820c dynamic allocation of work space 821 kvabkl = kend1 822 kend2 = kvabkl + mbas1(isyma)*nmatkl(isymkl) 823 lwrk2 = lwork - kend2 824 if (lwrk2.lt.0) then 825 call quit('insufficient work space in r12mkcv') 826 end if 827 828c 829 do b = 1, mbas1(isymb) 830 call cc_r12sortv(work(kvabkl),vabkl,b,isymb,isymv,isyma) 831 if (locdbg) then 832 xv = xv + ddot(mbas1(isyma)*nmatkl(isymkl), 833 & work(kvabkl),1,work(kvabkl),1) 834 end if 835 do isymj = 1, nsym 836 isymi = muld2h(isymij,isymj) 837 isymai = muld2h(isyma,isymi) 838 isymbj = muld2h(isymb,isymj) 839 kc2am = kend2 840 kom = kc2am + nrhf(isymi)*nmatkl(isymkl) 841 kend3 = kom + mbas1(isyma)*nrhf(isymi) 842 lwrk3 = lwork - kend3 843 if (lwrk3.lt.0) then 844 call quit('insufficient work space in r12mkcv') 845 end if 846 do j = 1, nrhf(isymj) 847 call cc_r12sortc(work(kc2am),c2am,j,isymc,isymj,isymi) 848 if (locdbg) then 849 xc = xc + ddot(nrhf(isymi)*nmatkl(isymkl), 850 & work(kc2am),1,work(kc2am),1) 851 end if 852c 853 ntoti = max(nrhf(isymi),1) 854 ntota = max(mbas1(isyma),1) 855 call dgemm('N','T',mbas1(isyma),nrhf(isymi), 856 & nmatkl(isymkl),one,work(kvabkl),ntota, 857 & work(kc2am),ntoti,zero,work(kom),ntota) 858 if (locdbg) then 859 x1 = x1 + ddot(mbas1(isyma)*nrhf(isymi), 860 & work(kom),1,work(kom),1) 861 end if 862 863c update Omega_alpha_i,beta_j which is stored as a triangular 864c matrix with bj.ge.ai 865 idxbj = it1ao(isymb,isymj)+mbas1(isymb)*(j-1)+b 866 koffai = it1ao(isyma,isymi) 867 if (isymai.eq.isymbj) then 868 maxai = min(koffai+mbas1(isyma)*nrhf(isymi),idxbj) 869 do idxai = koffai+1, maxai 870 naibj = it2ao(isymai,isymbj)+index(idxai,idxbj) 871 if (ioptom.eq.0) then 872 work(kom2-1+naibj) = work(kom2-1+naibj) + 873 & factor*work(kom-1+idxai-koffai) 874 else if (ioptom.eq.1) then 875 omega2(naibj) = omega2(naibj) + 876 & factor*work(kom-1+idxai-koffai) 877 end if 878 end do 879 else if (isymai.lt.isymbj) then 880 maxai = koffai+mbas1(isyma)*nrhf(isymi) 881 do idxai = koffai+1, maxai 882 naibj = it2ao(isymai,isymbj)+ 883 & nt1ao(isymai)*(idxbj-1) + idxai 884 if (ioptom.eq.0) then 885 work(kom2-1+naibj) = work(kom2-1+naibj) + 886 & factor*work(kom-1+idxai-koffai) 887 else if (ioptom.eq.1) then 888 omega2(naibj) = omega2(naibj) + 889 & factor*work(kom-1+idxai-koffai) 890 end if 891 end do 892 end if 893 end do 894 end do 895 896 end do 897 end do 898 end do 899 900 901 if (locdbg) then 902 write(lupri,*)'xv =', xv 903 write(lupri,*)'xc =', xc 904 write(lupri,*)'x1 =', x1 905 write(lupri,*)'in mkcv:' 906 if (ioptom.eq.0) then 907 write(lupri,*)'norm^2 work(kom2):', 908 & ddot(nt2ao(isymom),work(kom2),1,work(kom2),1) 909 call cc_prpao(dummy,work(kom2),isymom,0,1) 910 write(lupri,*)'norm^2 omega2: ', 911 & ddot(2*nt2ort(isymom),omega2,1,omega2,1) 912 else if (ioptom.eq.1) then 913 write(lupri,*)'norm^2 omega2:', 914 & ddot(nt2ao(isymom),omega2,1,omega2,1) 915 call cc_prpao(dummy,omega2,isymom,0,1) 916 end if 917 end if 918c 919 if (ioptom.eq.0) then 920c sort Omega as Omega_ab,_ij triplet and singlet ..... 921c note Omega is a triangular matrix with a.ge.b and i.ge.j 922 do isymi = 1, nsym 923 do isymj = 1, nsym 924 if (isymj.le.isymi) then 925 isymij = muld2h(isymi,isymj) 926 isymab = muld2h(isymij,isymom) 927 do isyma = 1, nsym 928 isymb = muld2h(isymab,isyma) 929 if (isymb.le.isyma) then 930 931 isymai = muld2h(isyma,isymi) 932 isymbj = muld2h(isymb,isymj) 933 isymbi = muld2h(isymb,isymi) 934 isymaj = muld2h(isyma,isymj) 935C if (isymom.ne.(muld2h(isymai,isymbj))) 936C * call quit('Symmetry error in cc_r12mkcv') 937 938 do i = 1, nrhf(isymi) 939 940 maxj = nrhf(isymj) 941 if (isymj.eq.isymi) maxj = i 942 do j = 1, maxj 943c 944 if (isymi.eq.isymj) then 945 idxij = imijp(isymi,isymj)+index(i,j) 946 else if (isymi.lt.isymj) then 947 idxij = imijp(isymi,isymj)+nrhf(isymi)*(j-1)+i 948 else 949 idxij = imijp(isymj,isymi)+nrhf(isymj)*(i-1)+j 950 end if 951 952 if (locdbg) then 953 write(lupri,'(a,11i5)') 954 & 'isymai,isymbj,ab,a,b,ij,i,j,i,j,idxij:', 955 & isymai,isymbj,isymab,isyma,isymb, 956 & isymij,isymi,isymj,i,j,idxij 957 end if 958c 959 do a = 1, mbas1(isyma) 960 maxb = mbas1(isymb) 961 if (isymb.eq.isyma) maxb = a 962 do b = 1, maxb 963c 964 if (isyma.eq.isymb) then 965 idxab = iaodpk(isyma,isymb)+index(a,b) 966 else if (isyma.lt.isymb) then 967 idxab = iaodpk(isyma,isymb)+mbas1(isyma)*(b-1)+a 968 else 969 idxab = iaodpk(isymb,isyma)+mbas1(isymb)*(a-1)+b 970 end if 971c 972 idxabijp = it2ort(isymab,isymij)+ 973 & nnbst(isymab)*(idxij-1)+idxab 974 idxabijm = nt2ort(isymom)+it2ort(isymab,isymij)+ 975 & nnbst(isymab)*(idxij-1)+idxab 976c 977 idxai = it1ao(isyma,isymi)+mbas1(isyma)*(i-1)+a 978 idxbj = it1ao(isymb,isymj)+mbas1(isymb)*(j-1)+b 979 idxbi = it1ao(isymb,isymi)+mbas1(isymb)*(i-1)+b 980 idxaj = it1ao(isyma,isymj)+mbas1(isyma)*(j-1)+a 981c 982 if (isymai.eq.isymbj) then 983 idxaibj = it2ao(isymai,isymbj)+index(idxai,idxbj) 984 idxbiaj = it2ao(isymbi,isymaj)+index(idxbi,idxaj) 985 986 omega2(idxabijp) = omega2(idxabijp) + 987 & (work(kom2-1+idxaibj)+work(kom2-1+idxbiaj)) 988 omega2(idxabijm) = omega2(idxabijm)+ 989 & (work(kom2-1+idxaibj)-work(kom2-1+idxbiaj)) 990 991 else 992 993 if (isymbj.lt.isymai) then 994 idxaibj = it2ao(isymbj,isymai)+ 995 & nt1ao(isymbj)*(idxai-1)+idxbj 996 else 997 idxaibj = it2ao(isymai,isymbj)+ 998 & nt1ao(isymai)*(idxbj-1)+idxai 999 end if 1000 1001 if (isymbi.lt.isymaj) then 1002 idxbiaj = it2ao(isymbi,isymaj)+ 1003 & nt1ao(isymbi)*(idxaj-1)+idxbi 1004 else 1005 idxbiaj = it2ao(isymaj,isymbi)+ 1006 & nt1ao(isymaj)*(idxbi-1)+idxaj 1007 end if 1008 1009 omega2(idxabijp) = omega2(idxabijp)+ 1010 & (work(kom2-1+idxaibj)+work(kom2-1+idxbiaj)) 1011 omega2(idxabijm) = omega2(idxabijm)+ 1012 & (work(kom2-1+idxaibj)-work(kom2-1+idxbiaj)) 1013 1014 end if 1015c 1016 if (.false.) then 1017 write(lupri,'(a,8i5,2g20.10)') 1018 & 'a,b,ai,bj,iaibj,ibiaj,iabijp,iabijm,O(aibj),O(biaj):', 1019 & a,b,idxai,idxbj,idxaibj,idxbiaj,idxabijp,idxabijm, 1020 & work(kom2-1+idxaibj),work(kom2-1+idxbiaj) 1021 write(lupri,*)'OM(+)', omega2(idxabijp) 1022 write(lupri,*)'OM(-)', omega2(idxabijm) 1023 end if 1024 1025 end do 1026 end do 1027c 1028 end do 1029 end do 1030 end if 1031 end do 1032 end if 1033 end do 1034 end do 1035c 1036 end if 1037 1038 if (.false.) then 1039c if (locdbg.and.(ioptom.eq.0)) then 1040 write(lupri,*)'norm^2 omega2+: ', 1041 & ddot(nt2ort(isymom),omega2,1,omega2,1) 1042 write(lupri,*)'norm^2 omega2-: ', 1043 & ddot(nt2ort(isymom),omega2(nt2ort(isymom)+1),1, 1044 & omega2(nt2ort(isymom)+1),1) 1045 end if 1046 1047 call qexit('r12mkcv') 1048 end 1049*====================================================================* 1050 subroutine cc_r12sortv(vs,vabkl,b,isymb,isymv,isyma) 1051c-------------------------------------------------------------------- 1052c purpose: sort V as V^beta_alpha,_mn with fixed beta 1053c expect V stored as a square matrix 1054c 1055c H. Fliegl, C. Haettig, summer 2004 1056c-------------------------------------------------------------------- 1057 implicit none 1058#include "priunit.h" 1059#include "ccorb.h" 1060#include "ccsdsym.h" 1061#include "r12int.h" 1062 1063 integer isymb,isymv,isyma,index,isymamn,isymmn,isymm,isymn, 1064 & isymam,isymbn, idxmn,idxamn 1065 integer idxab,idxabmn,isymab 1066#if defined (SYS_CRAY) 1067 real vs(*),vabkl(*) 1068#else 1069 double precision vs(*),vabkl(*) 1070#endif 1071C index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 1072 1073 call qenter('sortv') 1074 1075 isymab = muld2h(isyma,isymb) 1076 isymmn = muld2h(isymv,isymab) 1077 do idxmn = 1, nmatkl(isymmn) 1078 do a = 1, mbas1(isyma) 1079 idxab = iaodis(isyma,isymb)+mbas1(isyma)*(b-1)+a 1080 idxamn = mbas1(isyma)*(idxmn-1)+a 1081 idxabmn = ivabkl(isymab,isymmn)+ 1082 & n2bst(isymab)*(idxmn-1)+idxab 1083 vs(idxamn) = vabkl(idxabmn) 1084 end do 1085 end do 1086 1087 call qexit('sortv') 1088 end 1089*====================================================================* 1090 subroutine cc_r12mkvirt(vabkl,xlamd1,isymc1,xlamd2,isymc2, 1091 & filename,icon,work,lwork) 1092c-------------------------------------------------------------------- 1093c purpose: calculate V^ab_kl = 1094c \sum_{\alpha \beta} \Lambda1_{\alpha a} 1095c \Lambda2_{\beta b} V^{\alpha \beta}_kl 1096c 1097c icon = 2: calculate V^ab_kl = 1098c \sum_{\alpha \beta} [ \Lambda1_{\alpha a} 1099c \Lambda2_{\beta b} V^{\alpha \beta}_kl + 1100c Lambda2_{\alpha a} \Lambda1_{\beta b} 1101c V^{\alpha \beta}_kl] 1102c 1103c H.Fliegl, C.Haettig, summer 2004 1104c 1105c generalized for two transformation matrices: icon = 2 1106c C. Neiss, Nov. 2005 1107c-------------------------------------------------------------------- 1108 implicit none 1109#include "priunit.h" 1110#include "ccorb.h" 1111#include "ccsdsym.h" 1112#include "r12int.h" 1113#include "ccr12int.h" 1114 1115 logical locdbg 1116 parameter (locdbg = .false.) 1117 integer lwork,icount1 1118 integer kvabekl,kvcdkl,isym,idum,isym1,isym2, 1119 & kend1,lwrk1,kend2,lwrk2,kscr1,kscr2,kend3,lwrk3 1120 integer isymkl,isymalbe,isymal,isymbe,isyma,isymab,idxkl, 1121 & isymk,isyml,isymak,isymbl,idxab,idxakbl,idxak,idxbl, 1122 & kvabkl,koffc,koffr,ntotal,ntota,ntotbe,isymb,isymabe, 1123 & idxblak 1124 integer lunit,icon,index 1125 integer nck(8),ick(8,8),nvckdl(8),ivckdl(8,8),isymv,isymc1,isymc2 1126 character*(*) filename 1127#if defined (SYS_CRAY) 1128 real vabkl(*),xlamd1(*),xlamd2(*),work(*),one,zero, 1129 & ddot,x2,x3,x1 1130#else 1131 double precision vabkl(*),xlamd1(*),xlamd2(*),work(*),one,zero, 1132 & ddot,x2,x3,x1 1133#endif 1134 parameter(one = 1.0d0 , zero = 0.0d0) 1135 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 1136 1137 call qenter('mkvirt') 1138 1139c calculate some offsets and dimensions 1140 do isym = 1, nsym 1141 nck(isym) = 0 1142 icount1 = 0 1143 do isym2 = 1, nsym 1144 isym1 = muld2h(isym,isym2) 1145 nck(isym) = nck(isym) + nvir(isym1)*nrhfb(isym2) 1146 ick(isym1,isym2) = icount1 1147 icount1 = icount1 + nvir(isym1)*nrhfb(isym2) 1148 end do 1149 end do 1150 do isym = 1, nsym 1151 nvckdl(isym) = 0 1152 icount1 = 0 1153 do isym2 = 1, nsym 1154 isym1 = muld2h(isym,isym2) 1155 if (isym2.gt.isym1) then 1156 nvckdl(isym) = nvckdl(isym) + nck(isym1)*nck(isym2) 1157 ivckdl(isym1,isym2) = icount1 1158 ivckdl(isym2,isym1) = icount1 1159 icount1 = icount1 + nck(isym1)*nck(isym2) 1160 else if (isym2.eq.isym1) then 1161 nvckdl(isym) = nvckdl(isym) + nck(isym1)*(nck(isym2)+1)/2 1162 ivckdl(isym1,isym2) = icount1 1163 icount1 = icount1 + nck(isym1)*(nck(isym2)+1)/2 1164 end if 1165 end do 1166 end do 1167 1168 isymv = muld2h(isymc1,isymc2) 1169 1170 kvcdkl = 1 1171 kend1 = kvcdkl + nvckdl(isymv) 1172 lwrk1 = lwork - kend1 1173 if (lwrk1.lt.0) then 1174 call quit('insufficient work space in mkvirt') 1175 end if 1176 1177c transform to virtual 1178 call dzero(work(kvcdkl),nvckdl(isymv)) 1179 if (locdbg) then 1180 x2 = zero 1181 x3 = zero 1182 end if 1183 1184 do isymkl = 1, nsym 1185 isymalbe = isymkl 1186 do isymk = 1, nsym 1187 isyml = muld2h(isymkl,isymk) 1188 1189 do k = 1, nrhfb(isymk) 1190 do l = 1, nrhfb(isyml) 1191 idxkl = imatkl(isymk,isyml)+nrhfb(isymk)*(l-1)+k 1192 1193 do isymal = 1, nsym 1194 isyma = muld2h(isymc1,isymal) 1195 isymbe = muld2h(isymalbe,isymal) 1196 isymb = muld2h(isymc2,isymbe) 1197 isymabe = muld2h(isyma,isymbe) 1198 1199 kvabekl = kend1 1200 kscr1 = kvabekl + nvir(isyma)*nbas(isymbe) 1201 kend2 = kscr1 + nvir(isyma)*nvir(isymb) 1202 lwrk2 = lwork - kend2 1203 if (lwrk2.lt.0) then 1204 call quit('insufficient work space in mkvirt') 1205 end if 1206 1207 kvabkl = ivabkl(isymalbe,isymkl)+ 1208 & n2bst(isymalbe)*(idxkl-1)+ 1209 & iaodis(isymal,isymbe)+1 1210 koffc = iglmvi(isymal,isyma)+1 1211 1212 ntotal = max(nbas(isymal),1) 1213 ntota = max(nvir(isyma),1) 1214 1215 call dgemm('T','N',nvir(isyma),nbas(isymbe), 1216 & nbas(isymal),one,xlamd1(koffc),ntotal, 1217 & vabkl(kvabkl),ntotal,zero,work(kvabekl), 1218 & ntota) 1219 1220 koffc = iglmvi(isymbe,isymb)+1 1221 ntotbe = max(nbas(isymbe),1) 1222 ntota = max(nvir(isyma),1) 1223 1224 call dgemm('N','N',nvir(isyma),nvir(isymb),nbas(isymbe), 1225 & one,work(kvabekl),ntota,xlamd2(koffc),ntotbe, 1226 & zero,work(kscr1),ntota) 1227 1228 if (locdbg) then 1229 write(lupri,*) 'idxkl, norm^2(kscr1): ', 1230 & idxkl, ddot(nvir(isyma)*nvir(isymb), 1231 & work(kscr1),1,work(kscr1),1) 1232 x2 = x2 + ddot(nvir(isyma)*nvir(isymb), 1233 & work(kscr1),1,work(kscr1),1) 1234 end if 1235 1236C sort result as Vakbl and pack as triangular matrix 1237 isymak = muld2h(isyma,isymk) 1238 isymbl = muld2h(isymb,isyml) 1239 1240 do b = 1, nvir(isymb) 1241 idxbl = ick(isymb,isyml)+nvir(isymb)*(l-1)+b 1242 do a = 1, nvir(isyma) 1243 idxab = nvir(isyma)*(b-1)+a 1244 idxak = ick(isyma,isymk)+nvir(isyma)*(k-1)+a 1245 1246 if (isymak.eq.isymbl) then 1247 if (idxak.le.idxbl) then 1248 idxakbl = ivckdl(isymak,isymbl)+ 1249 & index(idxak,idxbl) 1250 work(kvcdkl-1+idxakbl) = work(kvcdkl-1+idxakbl)+ 1251 & work(kscr1-1+idxab) 1252 end if 1253 if ((idxbl.le.idxak).and.(icon.eq.2)) then 1254 idxblak = ivckdl(isymbl,isymak)+ 1255 & index(idxbl,idxak) 1256 work(kvcdkl-1+idxblak) = work(kvcdkl-1+idxblak)+ 1257 & work(kscr1-1+idxab) 1258 end if 1259 else if (isymak.lt.isymbl) then 1260 idxakbl = ivckdl(isymak,isymbl)+ 1261 & nck(isymak)*(idxbl-1)+idxak 1262 work(kvcdkl-1+idxakbl) = work(kvcdkl-1+idxakbl) + 1263 & work(kscr1-1+idxab) 1264 else if ((isymbl.lt.isymak).and.(icon.eq.2)) then 1265 idxblak = ivckdl(isymbl,isymak)+ 1266 & nck(isymbl)*(idxak-1)+idxbl 1267 work(kvcdkl-1+idxblak) = work(kvcdkl-1+idxblak) + 1268 & work(kscr1-1+idxab) 1269 end if 1270 1271 if (locdbg) then 1272 if (idxak.eq.idxbl) then 1273C sum diagonal elements 1274 x3 = x3 + ddot(1,work(kvcdkl-1+idxakbl),1, 1275 & work(kvcdkl-1+idxakbl),1) 1276 end if 1277 end if 1278 1279 end do 1280 end do 1281 1282 end do 1283 end do 1284 end do 1285 1286 end do 1287 end do 1288 1289 if (locdbg) then 1290 write(lupri,*)'x2 =', x2 1291 write(lupri,*)'x3 =', x3 1292 write(lupri,*)'nvckdl(1) =', nvckdl(1) 1293 1294c write(lupri,*)'norm^2 Vabkl =', 1295c & ddot(nvckdl(isymv),work(kvcdkl),1,work(kvcdkl),1) 1296 write(lupri,*)'norm^2 Vabkl =', 1297 & 2.0d0*ddot(nvckdl(isymv),work(kvcdkl),1,work(kvcdkl),1)-x3 1298 write(lupri,*)'norm^2 Vabkl triang =', 1299 & ddot(nvckdl(isymv),work(kvcdkl),1,work(kvcdkl),1) 1300 1301 if (isymv.eq.1) then 1302 call dreinorm(work(kvcdkl),nvckdl(1),nck(1),nsym,x1) 1303 write(lupri,*)'dreinorm =', x1 1304 end if 1305 do isymak = 1, nsym 1306 isymbl = isymak 1307 call outpkb(work(kvcdkl+ivckdl(isymak,isymbl)),nck(isymak), 1308 & 1,1,lupri) 1309 end do 1310 end if 1311 1312c write Vabkl on file 1313 lunit = -1 1314 call gpopen(lunit,filename,'unknown',' ','unformatted',idum, 1315 & .false.) 1316 rewind(lunit) 1317 write(lunit)(work(kvcdkl-1+i),i=1,nvckdl(isymv)) 1318 call gpclose(lunit,'KEEP') 1319 1320 call qexit('mkvirt') 1321 end 1322c====================================================================* 1323 subroutine dreinorm(a,na,ix,nsym,x) 1324c-------------------------------------------------------------------- 1325c purpose: calculate from triangular input matrix the norm 1326c of the corresponding quadratic matrix 1327c 1328c W. Klopper, summer 2004 1329c-------------------------------------------------------------------- 1330 implicit double precision (a-h,o-z) 1331 dimension a(*),ix(*) 1332 x=2d0*ddot(na,a,1,a,1) 1333 ij=0 1334 do isym=1,nsym 1335 do k=1,ix(isym) 1336 ij=ij+k 1337 x=x-a(ij)*a(ij) 1338 enddo 1339 enddo 1340 return 1341 end 1342c====================================================================* 1343 subroutine ccrhs_bpp(omega,t2am,isymt2am,lt2pcked, 1344 & filev,isymvint,work,lwork) 1345c-------------------------------------------------------------------- 1346c purpose: calculate B''-Term and update Omega_ijkl 1347c \sum_ab V^ab_kl * t^ij_ab 1348c 1349c H. Fliegl, C. Haettig, summer 2004 1350c 1351c modified C. Neiss, autumn 2005: 1352c - general symmetry of V 1353c - reduced memory requirement 1354c-------------------------------------------------------------------- 1355 implicit none 1356#include "priunit.h" 1357#include "ccorb.h" 1358#include "ccsdsym.h" 1359#include "r12int.h" 1360#include "ccr12int.h" 1361 1362 logical locdbg,lt2pcked 1363 parameter (locdbg = .false.) 1364 1365 integer isymt2am,lwork,kvpk,lwrk1,kend1,iopt,lunit, 1366 & isymvint,idum 1367 integer icount1,isym1,isym2,isym, 1368 & nck(8),ick(8,8),nvckdl(8),ivckdl(8,8) 1369 character*(*) filev 1370#if defined (SYS_CRAY) 1371 real omega(*),t2am(*),work(*),factor,ddot,one 1372#else 1373 double precision omega(*),t2am(*),work(*),factor,ddot,one 1374#endif 1375 parameter(one = 1.0d0) 1376 1377 call qenter('ccrhs_bpp') 1378 1379c calculate some offsets 1380 do isym = 1, nsym 1381 nck(isym) = 0 1382 icount1 = 0 1383 do isym2 = 1, nsym 1384 isym1 = muld2h(isym,isym2) 1385 nck(isym) = nck(isym) + nvir(isym1)*nrhfb(isym2) 1386 ick(isym1,isym2) = icount1 1387 icount1 = icount1 + nvir(isym1)*nrhfb(isym2) 1388 end do 1389 end do 1390 do isym = 1, nsym 1391 nvckdl(isym) = 0 1392 icount1 = 0 1393 do isym2 = 1, nsym 1394 isym1 = muld2h(isym,isym2) 1395 if (isym2.gt.isym1) then 1396 nvckdl(isym) = nvckdl(isym) + nck(isym1)*nck(isym2) 1397 ivckdl(isym1,isym2) = icount1 1398 icount1 = icount1 + nck(isym1)*nck(isym2) 1399 else if (isym2.eq.isym1) then 1400 nvckdl(isym) = nvckdl(isym) + nck(isym1)*(nck(isym2)+1)/2 1401 ivckdl(isym1,isym2) = icount1 1402 icount1 = icount1 + nck(isym1)*(nck(isym2)+1)/2 1403 end if 1404 end do 1405 end do 1406 1407 kvpk = 1 1408 kend1 = kvpk + max(nvckdl(isymvint),nt2am(isymt2am)) 1409 lwrk1 = lwork - kend1 1410 if (lwrk1.lt.0) then 1411 call quit('insufficient work space in ccrhs_bpp') 1412 end if 1413 1414 if (locdbg) then 1415 write(lupri,*)'Norm^2 Omega^ij_kl on entry in ccrhs_bpp', 1416 & ddot(ntr12sq(muld2h(isymvint,isymt2am)), 1417 & omega,1,omega,1) 1418 write(lupri,*) 'Omega on entry in ccrhs_bpp:' 1419 call cc_prsqr12(omega,muld2h(isymvint,isymt2am),'T',1,.FALSE.) 1420 end if 1421 1422 if (.not.lt2pcked) then 1423c pack T2 amplitudes as a triangular matrix and put on T2AM 1424 iopt = 0 1425 call cc_t2pk(work(kvpk),t2am,isymt2am,iopt) 1426 if (nt2sq(isymt2am).lt.nt2am(isymt2am)) 1427 & call quit('Internal Error in CCRHS_BPP') 1428 call dcopy (nt2am(isymt2am),work(kvpk),1,t2am,1) 1429 end if 1430 1431 if (locdbg) then 1432 write(lupri,*)'Norm^2 T2 ampl. in ccrhs_bpp', 1433 & ddot(nt2am(isymt2am),t2am,1,t2am,1) 1434 write(lupri,*) 'T2 ampl. in ccrhs_bpp:' 1435 call cc_prp(one,t2am,isymt2am,0,1) 1436 end if 1437 1438c read V^ab_kl 1439 lunit = -1 1440 call gpopen(lunit,filev,'old',' ','unformatted',idum,.false.) 1441 read(lunit)(work(kvpk-1+i),i=1,nvckdl(isymvint)) 1442 call gpclose(lunit,'KEEP') 1443 1444c calculate Omega_ijkl 1445 factor = one 1446 call cc_r12mi2(omega,work(kvpk),t2am,isymvint,isymt2am, 1447 & factor,work(kend1),lwrk1) 1448 1449 if (locdbg) then 1450 write(lupri,*)'Norm^2 Omega^ij_kl on exit in ccrhs_bpp', 1451 & ddot(ntr12sq(muld2h(isymvint,isymt2am)), 1452 & omega,1,omega,1) 1453 write(lupri,*) 'Omega on exit in ccrhs_bpp:' 1454 call cc_prsqr12(omega,muld2h(isymvint,isymt2am),'T',1,.FALSE.) 1455 end if 1456 1457 if (.not.lt2pcked) then 1458c restore T2 amplitudes as quadratic matrix 1459 call dcopy(nt2am(isymt2am),t2am,1,work(kvpk),1) 1460 call cc_t2sq(work(kvpk),t2am,isymt2am) 1461 end if 1462 1463 call qexit('ccrhs_bpp') 1464 end 1465c=====================================================================* 1466 1467c====================================================================* 1468 subroutine cc_r12mktbjpk(tbjpk,work,lwork) 1469c-------------------------------------------------------------------- 1470c purpose: calculate t(bj,p'k) = 1471c \sum_mn c_mn^jk * r_bp'^mn 1472c 1473c C. Neiss, spring 2006 1474c-------------------------------------------------------------------- 1475 implicit none 1476#include "priunit.h" 1477#include "ccorb.h" 1478#include "ccsdsym.h" 1479#include "dummy.h" 1480#include "r12int.h" 1481#include "ccr12int.h" 1482 1483 logical locdbg 1484 parameter (locdbg = .false.) 1485 1486 integer isymc,isymp,isymb,isymj,isymk,isymbmn,isymmn,isymjmn, 1487 & isympk,isymbj,isym,isym1,isym2 1488 integer nramn(8),iramn(8,8),nramnq(8),iramnq(8,8) 1489 integer kcamp,krbmn,kcjmn,idxbj,idxpk 1490 integer lwork,kend1,lwrk1,kend2,lwrk2,kend3,lwrk3 1491 integer iopt,lunit,iadr,len,koff,koffr 1492 character*10 model 1493 character*8 framnp 1494#if defined (SYS_CRAY) 1495 real work(*), tbjpk(*), ddot, one, zero 1496#else 1497 double precision work(*), tbjpk(*), ddot, one, zero 1498#endif 1499 parameter(zero=0.0D0, one=1.0D0) 1500 parameter(framnp='R12RMNAP') 1501 1502 call qenter('cc_r12mktbjpk') 1503C 1504 if (locdbg) then 1505 write(lupri,*) 'Entered cc_r12mktbjpk' 1506 end if 1507C 1508C calc. some dimensions/offsets: 1509 do isym = 1, nsym 1510 nramn(isym) = 0 1511 do isym2 = 1, nsym 1512 isym1 = muld2h(isym,isym2) 1513 iramn(isym1,isym2) = nramn(isym) 1514 nramn(isym) = nramn(isym) + nvir(isym1)*nmatkl(isym2) 1515 end do 1516 end do 1517C 1518 do isym = 1, nsym 1519 nramnq(isym) = 0 1520 do isym2 = 1, nsym 1521 isym1 = muld2h(isym,isym2) 1522 iramnq(isym1,isym2) = nramnq(isym) 1523 nramnq(isym) = nramnq(isym) + nramn(isym1)*norb2(isym2) 1524 end do 1525 end do 1526C 1527C initialisation: 1528 isymc = 1 1529 call dzero(tbjpk,ntg2sq(isymc)) 1530C 1531C allocate memory: 1532 kcamp = 1 1533 kend1 = kcamp + ntr12am(isymc) 1534 lwrk1 = lwork - kend1 1535 if (lwrk1.lt.0) then 1536 call quit('Insufficient memory in CC_R12MKTBJPK') 1537 end if 1538C 1539C read r12 amlitudes: 1540 iopt = 32 1541 call cc_rdrsp('R0 ',0,isymc,iopt,model,dummy,work(kcamp)) 1542C 1543C open file with r12-integrals r_(b,mn)^p' 1544 lunit = -1 1545 call wopen2(lunit,framnp,64,0) 1546C 1547C loop over contractions: 1548 do isymp = 1, nsym 1549 isymbmn = isymp 1550 do p = 1, norb2(isymp) 1551 len = nramn(isymbmn) 1552 krbmn = kend1 1553 kend2 = krbmn + len 1554 lwrk2 = lwork - kend2 1555 if (lwrk2.lt.0) then 1556 call quit('Insufficient memory in CC_R12MKTBJPK') 1557 end if 1558 iadr = iramnq(isymbmn,isymp) + 1559 & nramn(isymbmn)*(p-1) + 1 1560 call getwa2(lunit,framnp,work(krbmn),iadr,len) 1561C 1562C if (locdbg) then 1563C write(lupri,*) 'isymp, p, iadr= ',isymp,p,iadr 1564C write(lupri,*) 'Norm^2 (r-int):', ddot(len,work(krbmn),1, 1565C & work(krbmn),1) 1566C do isym2 = 1, nsym 1567C isym1 = muld2h(isymbmn,isym2) 1568C write(lupri,*) 'Symmetry block: ',isym1,isym2 1569C call output(work(krbmn+iramn(isym1,isym2)),1,nvir(isym1), 1570C & 1,nmatkl(isym2),nvir(isym1), 1571C & nmatkl(isym2),1,lupri) 1572C end do 1573C end if 1574C 1575 do isymk = 1, nsym 1576 isymjmn = muld2h(isymc,isymk) 1577 isympk = muld2h(isymp,isymk) 1578 do isymj = 1, nsym 1579 isymmn = muld2h(isymjmn,isymj) 1580 isymb = muld2h(isymmn,isymbmn) 1581 isymbj = muld2h(isymb,isymj) 1582 kcjmn = kend2 1583 kend3 = kcjmn + nrhfa(isymj)*nmatkl(isymmn) 1584 lwrk3 = lwork - kend3 1585 if (lwrk3.lt.0) then 1586 call quit('Insufficient memory in CC_R12MKTBJPK') 1587 end if 1588 do k = 1, nrhfa(isymk) 1589 call cc_r12sortc(work(kcjmn),work(kcamp),k,isymc,isymk, 1590 & isymj) 1591 1592 idxbj = it1am(isymb,isymj) + 1 1593 idxpk = ig1am(isymp,isymk)+norb2(isymp)*(k-1)+p 1594 koffr = krbmn + iramn(isymb,isymmn) 1595 koff = itg2sq(isymbj,isympk) + 1596 & nt1am(isymbj)*(idxpk-1) + idxbj 1597 call dgemm('N','T',nvir(isymb),nrhfa(isymj), 1598 & nmatkl(isymmn),one, 1599 & work(koffr),max(nvir(isymb),1), 1600 & work(kcjmn),max(nrhfa(isymj),1), 1601 & zero,tbjpk(koff),max(nvir(isymb),1)) 1602 1603 1604 end do 1605 end do 1606 end do 1607 end do 1608 end do 1609C 1610 if (locdbg) then 1611 write(lupri,*) 'Norm^2(c*r):',ddot(ntg2sq(isymc),tbjpk,1, 1612 & tbjpk,1) 1613 do isym2 = 1, nsym 1614 isym1 = muld2h(isymc,isym2) 1615 koff = 1 + itg2sq(isym1,isym2) 1616 write(lupri,*) 'Symmetry block:',isym1,isym2 1617 call output(tbjpk(koff),1,nt1am(isym1),1,ng1am(isym2), 1618 & nt1am(isym1),ng1am(isym2),1,lupri) 1619 end do 1620 end if 1621C 1622C close file 1623 call wclose2(lunit,framnp,'KEEP') 1624C 1625 if (locdbg) then 1626 write(lupri,*) 'Leaving cc_r12mktbjpk' 1627 end if 1628C 1629 call qexit('cc_r12mktbjpk') 1630 return 1631 end 1632c=====================================================================* 1633 1634c====================================================================* 1635 subroutine ccrhs_e1pim(e1pim,cmox,ilamdx,xlamdh,work,lwork) 1636c-------------------------------------------------------------------- 1637c purpose: add T1-independent Fock-contribution to 1638c E1'_(a delta) intermediate 1639c and transform into MO basis 1640c 1641c E1'_(a p') = Sum_(delta) [ CMOX_(delta p') * E1'_(a delta) + 1642c Sum_(gamma) XLAMDH_(gamma a) * F_(gamma delta)] 1643c (delta runs over MO- and aux-basis) 1644c 1645c C. Neiss, spring 2006 1646c-------------------------------------------------------------------- 1647 implicit none 1648#include "priunit.h" 1649#include "ccorb.h" 1650#include "ccsdsym.h" 1651#include "dummy.h" 1652#include "r12int.h" 1653#include "ccr12int.h" 1654 1655 logical locdbg,lidxtf1 1656 parameter (locdbg = .false.) 1657 1658 integer lwork,kend1,lwrk1,lunit,koff1,koff2,koff3 1659 integer kfgdp,isymfck,isymh,isyme1p,isymd,isyma,isymp 1660 integer isym,isym1,isym2,nbas2(8),ngdp(8),igdp(8,8), 1661 & nadp(8),iadp(8,8),ilamdx(8,8),napp(8),iapp(8,8) 1662#if defined (SYS_CRAY) 1663 real work(*), e1pim(*), cmox(*), xlamdh(*), one, zero 1664#else 1665 double precision work(*), e1pim(*), cmox(*), xlamdh(*), one, zero 1666#endif 1667 parameter (one=1.0D0, zero=0.0D0) 1668C 1669 call qenter('ccrhs_e1pim') 1670C 1671C calc. dimensions/offsets: 1672 do isym = 1, nsym 1673 nbas2(isym) = mbas1(isym) + mbas2(isym) 1674 end do 1675 do isym = 1, nsym 1676 ngdp(isym) = 0 1677 nadp(isym) = 0 1678 napp(isym) = 0 1679 do isym2 = 1, nsym 1680 isym1 = muld2h(isym,isym2) 1681 igdp(isym1,isym2) = ngdp(isym) 1682 ngdp(isym) = ngdp(isym) + mbas1(isym1)*nbas2(isym2) 1683 iadp(isym1,isym2) = nadp(isym) 1684 nadp(isym) = nadp(isym) + nvir(isym1)*nbas2(isym2) 1685 iapp(isym1,isym2) = napp(isym) 1686 napp(isym) = napp(isym) + nvir(isym1)*norb2(isym2) 1687 end do 1688 end do 1689C 1690C initialisations: 1691 isymh = 1 1692 isymfck = 1 1693 isyme1p = 1 1694C 1695 kfgdp = 1 1696 kend1 = kfgdp + ngdp(isymfck) 1697 lwrk1 = lwork - kend1 1698 if (lwrk1.lt.0) then 1699 call quit('Insufficient memory in CCRHS_E1PIM') 1700 end if 1701C 1702C read F_(gamma delta); delta runs over MO- and aux-basis (NBAS2) 1703 lunit = -1 1704 call gpopen(lunit,'R12FOCK','UNKNOWN',' ','UNFORMATTED',idummy, 1705 & .false.) 1706 call readt(lunit,ngdp(isymfck),work(kfgdp)) 1707 call gpclose(lunit,'KEEP') 1708C 1709C transform first index to virtual and add on E1PIM (which is 1710C half-transformed already) 1711 lidxtf1 = .TRUE. 1712 call cc_r12generaltf(work(kfgdp),e1pim,idummy,isymfck,xlamdh, 1713 & isymh,dummy,idummy,lidxtf1,iglmvi,idummy, 1714 & nbas2,nvir,nbas,igdp,idummy,idummy, 1715 & iadp,idummy,work(kend1),lwrk1) 1716 if (locdbg) then 1717 write(lupri,*)'E1PIM after addition of Fock-Matrix:' 1718 do isym2 = 1, nsym 1719 isym1 = muld2h(isyme1p,isym2) 1720 write(lupri,*)'Symmetry block:',isym1,isym2 1721 call output(e1pim(1+iadp(isym1,isym2)), 1722 & 1,nvir(isym1),1,nbas2(isym2), 1723 & nvir(isym1),nbas2(isym2),1,lupri) 1724 end do 1725 end if 1726C 1727C transform second index of e1pim to NORB2: 1728C (reusing the work space from the Fock-matrix) 1729 call dzero(work(kfgdp),ngdp(isymfck)) 1730 do isymd = 1, nsym 1731 isyma = muld2h(isyme1p,isymd) 1732 isymp = isymd !cmox is total-symmetric 1733 koff1 = 1 + iadp(isyma,isymd) 1734 koff2 = 1 + ilamdx(isymd,isymp) + nbas2(isymd)*norb1(isymp) 1735 koff3 = kfgdp + iapp(isyma,isymp) 1736 call dgemm('N','N',nvir(isyma),norb2(isymp),nbas2(isymd), 1737 & one,e1pim(koff1),max(nvir(isyma),1), 1738 & cmox(koff2),max(nbas2(isymd),1), 1739 & zero,work(koff3),max(nvir(isyma),1)) 1740 if (locdbg) then 1741 write(lupri,*) 'CMOX:',isymd,isymp 1742 call output(cmox(koff2),1,nbas2(isymd),1,norb2(isymp), 1743 & nbas2(isymd),norb2(isymp),1,lupri) 1744 end if 1745 end do 1746C 1747C copy result to E1PIM: 1748 call dcopy(napp(isyme1p),work(kfgdp),1,e1pim,1) 1749C 1750 call qexit('ccrhs_e1pim') 1751 return 1752 end 1753c=====================================================================* 1754 1755