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_r12vunpack(vijkl,isymv,vsing,vtrip,lprojv, 20 & nr12orb,nrhforb) 21*----------------------------------------------------------------------* 22c purpose: make V^{it jt }_{kl} = P^{kl}_{ij} * V^{it jt }_{kl} 23c and set it equal to the vector function omega for s = 0,1 24c (no symmetry is used!) 25c 26c H. Fliegl, C. Haettig spring 2003 27c 28c modified by C. Neiss, summer 2005 29c---------------------------------------------------------------------- 30 implicit none 31#include "priunit.h" 32#include "ccsdsym.h" 33#include "ccorb.h" 34 35 logical locdbg,lprojv 36 parameter(locdbg = .false.) 37#if defined (SYS_CRAY) 38 real half 39#else 40 double precision half 41#endif 42 parameter ( half = 0.5d0 ) 43 44 integer idxkl,idxlk,idxij,idxijkl,idxijlk,nrhftria,idxjikl,idxjilk 45 integer isymkl,isymij,isyml,isymk,isymi,isymj,idxji 46 integer it,jt,kt,lt,kl,ij,klij,isymv,icoun1,icoun2,isym 47 integer nrhforb(8),irhforb(8),nij(8),iij(8,8),iijkl(8,8) 48 integer nr12orb(8),ir12orb(8),nkl(8),ikl(8,8),nr12t 49 integer index 50#if defined (SYS_CRAY) 51 real vijkl(*),vsing(*) 52 real vtrip(*),ff 53#else 54 double precision vijkl(*),vsing(*) 55 double precision vtrip(*),ff 56#endif 57 58 index(i,j) = max(i,j)*(max(i,j) - 3)/2 + i + j 59 60 call qenter('cc_r12vunpack') 61C 62 ff = 1.0d0/sqrt(2.0d0) 63C 64 nr12t = 0 65 icoun1 = 0 66 icoun2 = 0 67 do isym = 1, nsym 68 irhforb(isym) = icoun1 69 ir12orb(isym) = icoun2 70 icoun1 = icoun1 + nrhforb(isym) 71 icoun2 = icoun2 + nr12orb(isym) 72 nr12t = nr12t + nr12orb(isym) 73 end do 74 75 do isym = 1, nsym 76 icoun1 = 0 77 icoun2 = 0 78 do isymj = 1, nsym 79 isymi = muld2h(isymj,isym) 80 iij(isymi,isymj) = icoun1 81 ikl(isymi,isymj) = icoun2 82 icoun1 = icoun1 + nrhforb(isymi)*nrhforb(isymj) 83 icoun2 = icoun2 + nr12orb(isymi)*nr12orb(isymj) 84 end do 85 nij(isym) = icoun1 86 nkl(isym) = icoun2 87 end do 88 89 do isym = 1, nsym 90 icoun1 = 0 91 do isymij = 1, nsym 92 isymkl = muld2h(isymij,isym) 93 iijkl(isymij,isymkl) = icoun1 94 icoun1 = icoun1 + nij(isymij)*nkl(isymkl) 95 end do 96 end do 97C 98 nrhftria = nr12t*(nr12t+1)/2 99 call dzero(vsing,nrhftria*nrhftria) 100 call dzero(vtrip,nrhftria*nrhftria) 101 102 do isymkl = 1, nsym 103 isymij = muld2h(isymkl,isymv) 104 do isyml =1, nsym 105 isymk = muld2h(isymkl,isyml) 106 107 do k = 1, nr12orb(isymk) 108 kt = ir12orb(isymk) + k 109 do l = 1, nr12orb(isyml) 110 lt = ir12orb(isyml) + l 111 if (lt.le.kt) then 112 kl = index(kt,lt) 113 idxkl=ikl(isymk,isyml)+nr12orb(isymk)*(l-1)+k 114 idxlk=ikl(isyml,isymk)+nr12orb(isyml)*(k-1)+l 115 do isymi =1, nsym 116 isymj = muld2h(isymij,isymi) 117 118 do j = 1, nrhforb(isymj) 119 jt = irhforb(isymj)+j 120 do i = 1, nrhforb(isymi) 121 it = irhforb(isymi)+i 122 if (jt.le.it) then 123 ij = index(it,jt) 124 klij = nrhftria*(ij-1) + kl 125 idxij = iij(isymi,isymj)+ 126 & nrhforb(isymi)*(j-1) + i 127 idxijkl = iijkl(isymij,isymkl)+ 128 & nij(isymij)*(idxkl-1)+idxij 129 idxijlk = iijkl(isymij,isymkl)+ 130 & nij(isymij)*(idxlk-1)+idxij 131 132 vsing(klij)=vijkl(idxijkl)+ 133 & vijkl(idxijlk) 134 vtrip(klij)=vijkl(idxijkl)- 135 & vijkl(idxijlk) 136c ----------------------------------------- 137c V is not symmetric in ik <-> jl, thus 138c use 0.5 (V(ijkl) + V(jilk)) 139c ----------------------------------------- 140 if (lprojv) then 141 idxji = iij(isymj,isymi)+ 142 & nrhforb(isymj)*(i-1) + j 143 idxjilk = iijkl(isymij,isymkl)+ 144 & nij(isymij)*(idxlk-1)+idxji 145 idxjikl = iijkl(isymij,isymkl)+ 146 & nij(isymij)*(idxkl-1)+idxji 147 vsing(klij) = half*(vsing(klij) 148 & + vijkl(idxjikl)+vijkl(idxjilk)) 149 vtrip(klij) = half*(vtrip(klij) 150 & + vijkl(idxjilk)-vijkl(idxjikl)) 151 end if 152 153 if (it.eq.jt) 154 & vsing(klij)=ff*vsing(klij) 155 if (kt.eq.lt) 156 & vsing(klij)=ff*vsing(klij) 157 if (it.eq.jt) 158 & vtrip(klij)=ff*vtrip(klij) 159 if (kt.eq.lt) 160 & vtrip(klij)=ff*vtrip(klij) 161 end if 162 end do 163 end do 164 165 end do 166 end if 167 168 end do 169 end do 170 end do 171 end do 172 173 if (locdbg) then 174 write(lupri,*)'in cc_r12vunpack:' 175 write(lupri,*)'Vsing' 176 call output(vsing,1,nrhftria,1,nrhftria,nrhftria, 177 & nrhftria,1,lupri) 178 write(lupri,*)'Vtrip' 179 call output(vtrip,1,nrhftria,1,nrhftria,nrhftria, 180 & nrhftria,1,lupri) 181 end if 182 183 call qexit('cc_r12vunpack') 184 end 185 186*=====================================================================* 187 subroutine cc_r12vpack(vijkl,isymv,vsing,vtrip,nr12orb,nrhforb, 188 & nijkl) 189*---------------------------------------------------------------------* 190c purpose: transform omega for s = 0,1 back to omega(ki,lj) 191c 192c modified by C. Neiss, summer 2005 193c see CCR12PCK for further comments 194c--------------------------------------------------------------------- 195 implicit none 196#include "ccsdsym.h" 197#include "priunit.h" 198#include "ccorb.h" 199 integer idxkl,idxlk,idxij,idxijkl,idxijlk,nrhftria,idxjikl,idxjilk 200 integer isymkl,isymij,isyml,isymk,isymi,isymj,idxji 201 integer it,jt,kt,lt,index,kl,ij,klij,isymv,isym,icoun1,icoun2 202 integer nr12orb(8),nrhforb(8),irhforb(8),ir12orb(8), 203 & nij(8),nkl(8),iij(8,8),ikl(8,8),nijkl(8),iijkl(8,8),nr12t 204#if defined (SYS_CRAY) 205 real vijkl(*),vsing(*) 206 real vtrip(*),ff 207 real ffkl,half 208#else 209 double precision vijkl(*),vsing(*) 210 double precision vtrip(*),ff 211 double precision ffkl,half 212#endif 213 logical locdbg 214 parameter(locdbg = .false.) 215 parameter (half = 0.5D0) 216 217 index(i,j) = max(i,j)*(max(i,j) - 3)/2 + i + j 218 219 call qenter('cc_r12vpack') 220C 221 nr12t = 0 222 icoun1 = 0 223 icoun2 = 0 224 do isym = 1, nsym 225 irhforb(isym) = icoun1 226 ir12orb(isym) = icoun2 227 icoun1 = icoun1 + nrhforb(isym) 228 icoun2 = icoun2 + nr12orb(isym) 229 nr12t = nr12t + nr12orb(isym) 230 end do 231 232 do isym = 1, nsym 233 icoun1 = 0 234 icoun2 = 0 235 do isymj = 1, nsym 236 isymi = muld2h(isymj,isym) 237 iij(isymi,isymj) = icoun1 238 ikl(isymi,isymj) = icoun2 239 icoun1 = icoun1 + nrhforb(isymi)*nrhforb(isymj) 240 icoun2 = icoun2 + nr12orb(isymi)*nr12orb(isymj) 241 end do 242 nij(isym) = icoun1 243 nkl(isym) = icoun2 244 end do 245 246 do isym = 1, nsym 247 icoun1 = 0 248 do isymij = 1, nsym 249 isymkl = muld2h(isymij,isym) 250 iijkl(isymij,isymkl) = icoun1 251 icoun1 = icoun1 + nij(isymij)*nkl(isymkl) 252 end do 253 nijkl(isym) = icoun1 254 end do 255C 256 nrhftria = nr12t*(nr12t+1)/2 257 call dzero(vijkl,nijkl(isymv)) 258 259 do isymkl = 1, nsym 260 isymij = muld2h(isymkl,isymv) 261 do isyml =1, nsym 262 isymk = muld2h(isymkl,isyml) 263 264 do k = 1, nr12orb(isymk) 265 kt = ir12orb(isymk) + k 266 do l = 1, nr12orb(isyml) 267 lt = ir12orb(isyml) + l 268 269 if (lt.eq.kt) then 270 ffkl = sqrt(2d0) 271 else 272 ffkl = 1d0 273 end if 274 275 if (lt.le.kt) then 276 kl = index(kt,lt) 277 idxkl=ikl(isymk,isyml)+nr12orb(isymk)*(l-1)+k 278 idxlk=ikl(isyml,isymk)+nr12orb(isyml)*(k-1)+l 279 280 do isymi =1, nsym 281 isymj = muld2h(isymij,isymi) 282 283 do j = 1, nrhforb(isymj) 284 jt = irhforb(isymj)+j 285 do i = 1, nrhforb(isymi) 286 it = irhforb(isymi)+i 287 288 if (jt.eq.it) then 289 ff = ffkl * sqrt(2d0) 290 else 291 ff = ffkl 292 end if 293 294 if (jt.le.it) then 295 ij = index(it,jt) 296 klij = nrhftria*(ij-1)+kl 297 idxij = iij(isymi,isymj)+ 298 & nrhforb(isymi)*(j-1) + i 299 idxijkl = iijkl(isymij,isymkl)+ 300 & nij(isymij)*(idxkl-1)+idxij 301 idxijlk = iijkl(isymij,isymkl)+ 302 & nij(isymij)*(idxlk-1)+idxij 303 idxji = iij(isymj,isymi)+ 304 & nrhforb(isymj)*(i-1) + j 305 idxjilk = iijkl(isymij,isymkl)+ 306 & nij(isymij)*(idxlk-1)+idxji 307 idxjikl = iijkl(isymij,isymkl)+ 308 & nij(isymij)*(idxkl-1)+idxji 309 vijkl(idxijkl) = ff*half*(vsing(klij) 310 & + vtrip(klij)) 311 vijkl(idxijlk) = ff*half*(vsing(klij) 312 & - vtrip(klij)) 313 vijkl(idxjilk) = vijkl(idxijkl) 314 vijkl(idxjikl) = vijkl(idxijlk) 315 end if 316 317 end do 318 end do 319 320 end do 321 322 end if 323 324 end do 325 end do 326 327 end do 328 end do 329 330 if (locdbg) then 331 write(lupri,*) 'Result in CC_R12VPACK:' 332 do isymij = 1, nsym 333 isymkl = muld2h(isymij,isymv) 334 write(lupri,*) 'Symmetry block (ij,kl):', 335 & isymij,isymkl 336 if (nkl(isymkl).eq.0 .or. nij(isymij).eq.0) then 337 write(lupri,*) 'This symmetry is empty' 338 else 339 call output(vijkl(iijkl(isymij,isymkl)+1),1,nij(isymij), 340 & 1,nkl(isymkl),nij(isymij),nkl(isymkl),1,lupri) 341 end if 342 end do 343 end if 344 345 call qexit('cc_r12vpack') 346 end 347 348*=====================================================================* 349 subroutine ccrhs_ep_old(vijkl,isymv,lcont,cont,work,lwork, 350 & iamp,fc12am,lufc12,frho12,lufr12, 351 & ifile,aproxr12,basscl1,basscl2) 352*---------------------------------------------------------------------* 353c purpose: make the CC2-R12 vector function omega(ki,lj) 354c 355c H. Fliegl, (WK/UniKA/02-05-2003) 356c 357c C. Neiss April 2005: 358c adapted for left hand transformations 359c basscl1 = brascl for right hand transf. 360c = 0.5*ketscl for left hand transf. 361c basscl2 = ketscl for right hand transf. 362c = 2*brascl for left hand transf. 363c 364c NOTE: Does NOT work with .R12ORB 365c--------------------------------------------------------------------- 366 implicit none 367#include "ccsdsym.h" 368#include "priunit.h" 369#include "ccsdinp.h" 370#include "ccorb.h" 371#include "dummy.h" 372#include "r12int.h" 373#include "ccr12int.h" 374 375 logical lprojv,ldum,lcont 376 integer lwork,isymv,nrhftria,ksing,ktrip,kend1,lwrk1 377 integer kbsing,kbtrip,kcsing,kctrip,kevl,kxsing,kxtrip,ij 378 integer idum,ijcsvc,ijctvc,ijvsvc,ijvtvc,n2,lunit,ifile,kpck 379 integer lufc12,lufr12,ian,iap,kpckv,kpckc,kbmat,kxint 380 integer iamp,iopt 381 integer nkilj(8),nijkl(8) 382 character*(*) frho12, fc12am 383 character*3 aprox 384 character*10 model 385 character*(*) aproxr12 386#if defined (SYS_CRAY) 387 real vijkl(*),work(*), cont(2), ddot, basscl1, basscl2 388#else 389 double precision vijkl(*),work(*), cont(2), ddot, basscl1, 390 & basscl2 391#endif 392 logical locdbg 393 parameter (locdbg = .false.) 394 395 call qenter('ccrhs_ep') 396 397 if (locdbg) write(lupri,*) 'Entered CCRHS_EP' 398 399 nrhftria = nrhftb*(nrhftb+1)/2 400 n2 = nrhftria * nrhftria 401 402 ksing = 1 403 ktrip = ksing + n2 404 kbsing = ktrip + n2 405 kbtrip = kbsing + n2 406 kxsing = kbtrip + n2 407 kxtrip = kxsing + n2 408 kcsing = kxtrip + n2 409 kctrip = kcsing + n2 410 kevl = kctrip + n2 411 kpckc = kevl + nrhftria 412 kbmat = kpckc + ntr12am(isymv) 413 kxint = kbmat + nr12r12p(1) 414 kend1 = kxint + nr12r12p(1) 415 416 if (iamp .EQ. 1) then 417 kpck = kend1 418 kend1 = kpck + ntr12am(isymv) 419 end if 420 if (lcont) then 421 kpckv = kend1 422 kend1 = kpckv + ntr12am(isymv) 423 end if 424 425 lwrk1 = lwork - kend1 426 if (lwrk1 .lt. 0) then 427 call quit('Insufficient work space in ccrhs_ep') 428 end if 429 430 call dzero(work(kcsing),n2) 431 call dzero(work(kctrip),n2) 432 433c ---------------------------- 434c get V and omega for s = 0,1 435c ---------------------------- 436 lprojv = .true. 437 call cc_r12vunpack(vijkl,isymv,work(ksing),work(ktrip),lprojv, 438 & nrhfb,nrhf) 439 if (lcont) then 440 call ccr12pck(work(kpckv),isymv,work(ksing),work(ktrip),nrhfb, 441 & nrhf,nkilj) 442 end if 443 if (locdbg) then 444c call around('Vtilde integrals (singlet):') 445c call output(work(ksing),1,nrhftria,1,nrhftria, 446c & nrhftria,nrhftria,1,lupri) 447c call around('Vtilde integrals (triplet):') 448c call output(work(ktrip),1,nrhftria,1,nrhftria, 449c & nrhftria,nrhftria,1,lupri) 450 call cc_r12vpack(vijkl,isymv,work(ksing),work(ktrip),nrhfb, 451 & nrhf,nijkl) 452 call around('Vtilde on entry (after 0.5*P_kl^ij):') 453 call cc_prsqr12(vijkl,isymv,'T',1,.false.) 454 end if 455 456c ---------------------- 457c read orbital energies 458c --------------------- 459 lunit = -1 460 ldum = .false. 461 call gpopen(lunit,fccr12e,'old',' ','formatted', 462 & idum,ldum) 463 read(lunit,'(4e30.20)') (work(kevl+ij), ij = 0, nrhftria-1) 464 call gpclose(lunit,'KEEP') 465 466c ---------------- 467c read B matrices 468c ---------------- 469 lunit = -1 470 call gpopen(lunit,fccr12b,'old',' ','unformatted', 471 & idum,ldum) 472 8888 read(lunit) ian, iap, aprox 473 read(lunit) (work(kbmat+i), i=0, nr12r12p(1)-1 ) 474 if ((ian.ne.ianr12).or.(iap.ne.iapr12)) goto 8888 475 call gpclose(lunit,'KEEP') 476 if (aproxr12(1:3).ne.aprox) then 477 write(lupri,*)'aproxr12 =', aproxr12(1:3) 478 write(lupri,*)'aprox =', aprox(1:3) 479 call quit('--- Warning: inconsistent R12 approximation') 480 end if 481 call ccr12unpck(work(kbmat),1,work(kbsing),work(kbtrip),nrhfb, 482 & nrhfb) 483c ------------------------- 484c read R12 overlap matrices 485c ------------------------- 486 lunit = -1 487 call gpopen(lunit,fccr12x,'old',' ','unformatted', 488 & idum,ldum) 489 9999 read(lunit) ian 490 read(lunit) (work(kxint+i), i=0, nr12r12p(1)-1 ) 491 if (ian.ne.ianr12) goto 9999 492 call gpclose(lunit,'KEEP') 493 call ccr12unpck(work(kxint),1,work(kxsing),work(kxtrip),nrhfb, 494 & nrhfb) 495 496 if (iamp .EQ. 1) then 497c -------------------- 498c read trial vector R2 499c -------------------- 500 call cc_rvec(lufc12,fc12am,ntr12am(isymv),ntr12am(isymv), 501 * ifile,work(kpck)) 502 if (lcont) then 503 call dcopy(ntr12am(isymv),work(kpck),1,work(kpckc),1) 504 end if 505c if (locdbg) then 506c call cc_prpr12(work(kpck),isymv,1,.TRUE.) 507c end if 508 call cclr_diasclr12(work(kpck),basscl2,isymv) 509 call ccr12unpck(work(kpck),isymv,work(kcsing),work(kctrip), 510 & nrhfb,nrhf) 511 else 512c -------------------- 513c read R12 amplitudes 514c -------------------- 515c lunit = -1 516c call gpopen(lunit,fccr12c,'unknown',' ','unformatted', 517c & idum,ldum) 518c read(lunit) (work(kpckc+i), i = 0, ntr12am(1)-1 ) 519c call gpclose(lunit,'KEEP') 520 if (isymv.ne.1) call quit('Symmetry error in CCRHS_EP') 521 iopt = 32 522 call cc_rdrsp('R0 ',0,1,iopt,model,dummy,work(kpckc)) 523 call ccr12unpck(work(kpckc),isymv,work(kcsing),work(kctrip), 524 & nrhfb,nrhf) 525 end if 526 527 if (lcont) then 528 cont(1) = ddot(ntr12am(isymv),work(kpckc),1,work(kpckv),1) 529 end if 530 531c if (locdbg) then 532c call around('r12 vector function (singlet)') 533c call output(work(ksing),1,nrhftria,1,nrhftria, 534c & nrhftria,nrhftria,1,lupri) 535c call around('r12 vector function (triplet)') 536c call output(work(ktrip),1,nrhftria,1,nrhftria, 537c & nrhftria,nrhftria,1,lupri) 538c call around('r12 amplitudes (singlet)') 539c call output(work(kcsing),1,nrhftria,1,nrhftria, 540c & nrhftria,nrhftria,1,lupri) 541c call around('r12 amplitudes (triplet)') 542c call output(work(kctrip),1,nrhftria,1,nrhftria, 543c & nrhftria,nrhftria,1,lupri) 544c call around('b matrix (singlet)') 545c call output(work(kbsing),1,nrhftria,1,nrhftria, 546c & nrhftria,nrhftria,1,lupri) 547c call around('b matrix (triplet)') 548c call output(work(kbtrip),1,nrhftria,1,nrhftria, 549c & nrhftria,nrhftria,1,lupri) 550c call around('x matrix (singlet)') 551c call output(work(kxsing),1,nrhftria,1,nrhftria, 552c & nrhftria,nrhftria,1,lupri) 553c call around('x matrix (triplet)') 554c call output(work(kxtrip),1,nrhftria,1,nrhftria, 555c & nrhftria,nrhftria,1,lupri) 556c call around('orbital energies') 557c call outpak(work(kevl),nrhft,1,lupri) 558c end if 559 if ( debug ) then 560 write(lupri,*) 'ccrhs_ep> norm^2(omega_12 singlet):', 561 & ddot(n2,work(ksing),1,work(ksing),1) 562 write(lupri,*) 'ccrhs_ep> norm^2(omega_12 triplet):', 563 & ddot(n2,work(ktrip),1,work(ktrip),1) 564 if (iamp.EQ.1) write(lupri,*) 'ccrhs_ep> norm^2(C12 packed):', 565 & ddot(ntr12am(isymv),work(kpck),1,work(kpck),1) 566 write(lupri,*) 'ccrhs_ep> norm^2(C12 singlet):', 567 & ddot(n2,work(kcsing),1,work(kcsing),1) 568 write(lupri,*) 'ccrhs_ep> norm^2(C12 triplet):', 569 & ddot(n2,work(kctrip),1,work(kctrip),1) 570 end if 571 572c ------------------------------------------------- 573c loop over pairs and add B*c to omega for s = 0,1 574c ------------------------------------------------- 575 ijcsvc = kcsing 576 ijctvc = kctrip 577 ijvsvc = ksing 578 ijvtvc = ktrip 579 do ij = 1, nrhftria 580 call ccrhs_ep0(ij,nrhftria,work(kevl), 581 & work(ijvsvc),work(ijvtvc), 582 & work(ijcsvc),work(ijctvc), 583 & work(kbsing),work(kbtrip), 584 & work(kxsing),work(kxtrip)) 585 ijcsvc = ijcsvc + nrhftria 586 ijctvc = ijctvc + nrhftria 587 ijvsvc = ijvsvc + nrhftria 588 ijvtvc = ijvtvc + nrhftria 589 end do 590 591c -------------- 592c print results: 593c -------------- 594 if (locdbg) then 595c call around('r12 vector function (singlet)') 596c call output(work(ksing),1,nrhftria,1,nrhftria, 597c * nrhftria,nrhftria,1,lupri) 598c call around('r12 vector function (triplet)') 599c call output(work(ktrip),1,nrhftria,1,nrhftria, 600c * nrhftria,nrhftria,1,lupri) 601 call cc_r12vpack(vijkl,isymv,work(ksing),work(ktrip),nrhfb, 602 & nrhf,nijkl) 603 call around('OmegaR12 on exit') 604 call cc_prsqr12(vijkl,isymv,'T',1,.false.) 605 write(lupri,*) 'ccrhs_ep> norm^2(OmegaR12):', 606 & ddot(ntr12sq(isymv),vijkl,1,vijkl,1) 607 end if 608 if (locdbg) then 609 write(lupri,*) 'ccrhs_ep> norm^2(omega_12 singlet):', 610 & ddot(n2,work(ksing),1,work(ksing),1) 611 write(lupri,*) 'ccrhs_ep> norm^2(omega_12 triplet):', 612 & ddot(n2,work(ktrip),1,work(ktrip),1) 613 end if 614 615 if (iamp .EQ. 1) then 616c ------------------------------------------------------------ 617c for jacobian right transformations scale diagonal to 618c transform to pseudo-orthonormal basis 619c ------------------------------------------------------------ 620 call ccr12pck(work(kpck),isymv,work(ksing),work(ktrip),nrhfb, 621 & nrhf,nkilj) 622 call cclr_diasclr12(work(kpck),basscl1,isymv) 623 call cc_wvec(lufr12,frho12,ntr12am(isymv),ntr12am(isymv), 624 * ifile,work(kpck)) 625 if (locdbg) then 626 call around('Jacobian transformed vector in CCRHS_EP') 627 call cc_prpr12(work(kpck),isymv,1,.false.) 628 end if 629 else 630c ------------------------------------------------------------ 631c write result vector as ground state vector function to file: 632c ------------------------------------------------------------ 633 call ccr12pck(work(kpckc),isymv,work(ksing),work(ktrip),nrhfb, 634 & nrhf,nkilj) 635 lunit = -1 636 call gpopen(lunit,'CC_OMEGAR12','unknown',' ','unformatted', 637 & idum,ldum) 638 write(lunit) (work(kpckc+i), i=0, ntr12am(isymv)-1 ) 639 call gpclose(lunit,'KEEP') 640 end if 641c 642 if (lcont) then 643 call daxpy(ntr12am(isymv),-1.0d0,work(kpckv),1,work(kpck),1) 644 cont(2) = ddot(ntr12am(isymv),work(kpckc),1,work(kpck),1) 645 end if 646 647 if (locdbg) write(lupri,*) 'Leaving CCRHS_EP' 648 649 call qexit('ccrhs_ep') 650 end 651*=====================================================================* 652 subroutine ccrhs_ep0(ij,n,e,vsing,vtrip,csing,ctrip, 653 & bsing,btrip,xsing,xtrip) 654*---------------------------------------------------------------------* 655c purpose: make the CC2-R12 vector function 656c omega = sum_{m.le.n} B^{ij}(kl,mn)*c(mn,ij) for s = 0,1 657c 658c H. Fliegl, (WK/UniKA/02-05-2003) 659c--------------------------------------------------------------------- 660 implicit none 661#include "priunit.h" 662#include "r12int.h" 663 integer n,ij,kl,mn 664#if defined (SYS_CRAY) 665 real e(n),vsing(n),vtrip(n),csing(n),ctrip(n), 666 & bsing(n,n),btrip(n,n),xsing(n,n),xtrip(n,n), 667 & b,half,two,fact 668#else 669 double precision e(n),vsing(n),vtrip(n),csing(n),ctrip(n), 670 & bsing(n,n),btrip(n,n),xsing(n,n),xtrip(n,n), 671 & b,half,two,fact 672#endif 673 parameter (half = 0.5d0, two = 2d0) 674 675 call qenter('ccrhs_ep0') 676 677 if (ianr12.eq.1) then 678c if (iapr12.le.2) then 679 if (iapr12.le.2.or.iapr12.eq.4.or.iapr12.eq.6) then 680 fact = 0.0D0 681 else 682 fact = half 683 end if 684 else if (ianr12.eq.2 .or. ianr12.eq.3) then 685 if (iapr12.eq.2.or.iapr12.eq.5) then 686 fact = 0.0D0 687 else 688 fact = half 689 end if 690 end if 691c 692 do kl = 1, n 693 do mn = 1, n 694 b = - bsing(kl,mn) 695 & + 696 & fact * (e(kl) + e(mn) - two * e(ij)) * xsing(kl,mn) 697 vsing(kl) = vsing(kl) + b * csing(mn) 698 b = - btrip(kl,mn) 699 & + 700 & fact * (e(kl) + e(mn) - two * e(ij)) * xtrip(kl,mn) 701 vtrip(kl) = vtrip(kl) + b * ctrip(mn) 702 end do 703 end do 704 705 call qexit('ccrhs_ep0') 706 end 707 708*=====================================================================* 709 subroutine ccrhs_ep(vijkl,isymv,lcont,cont,work,lwork, 710 & iamp,fc12am,lufc12,frho12,lufr12, 711 & ifile,aproxr12,basscl1,basscl2) 712*---------------------------------------------------------------------* 713c purpose: make the CC2-R12 vector function omega(ki,lj): 714c take OmegaR12 = Vijkl and add OmegaR12_E' contribution, 715c write final OmegaR12 on disk. 716c alternative subroutine for old ccrhs_ep without 717c using singlet/triplet matrices in matrix 718c multiplication 719c 720c vijkl = OmegaR12 721c isymv = symmetry of result 722c basscl1 = brascl for right hand transf. 723c = 0.5*ketscl for left hand transf. 724c basscl2 = ketscl for right hand transf. 725c = 2*brascl for left hand transf. 726c 727c Christian Neiss, April. 2005, based on old ccrhs_ep 728c--------------------------------------------------------------------- 729 implicit none 730#include "ccsdsym.h" 731#include "priunit.h" 732#include "ccsdinp.h" 733#include "ccorb.h" 734#include "dummy.h" 735#include "r12int.h" 736#include "ccr12int.h" 737 738 logical ldum,lcont,lbmat 739 integer lwork,isymv,kend1,lwrk1,kend2,lwrk2,kevl,kekl,keij 740 integer iamp,idum,lunit,ifile,ian,iap,isym,iopt 741 integer lufc12,lufr12,kpck,kpckv,ktr12,kxintsq 742 integer ktr12sq,kbmatsq,koffv,kofftr12,koffb,kbscr 743 integer isymi,isymj,isymij,isymmn,isymkl,idxij 744 integer nkl,nij,inmatkl(8),inmatij(8),norbtsx 745 character*(*) frho12, fc12am 746 character*3 cdum 747 character*10 model 748 character*3 aproxr12,aprox 749#if defined (SYS_CRAY) 750 real vijkl(*),work(*), cont(2), ddot, basscl1, basscl2, one, 751 & factor, half 752#else 753 double precision vijkl(*),work(*), cont(2), ddot, basscl1, 754 & basscl2, one, factor, half 755#endif 756 logical locdbg 757 parameter (locdbg = .false.) 758 parameter (one = 1.0D0, half = 0.5D0) 759 760 call qenter('ccrhs_ep') 761 762 if (locdbg) write(lupri,*) 'Entered CCRHS_EP' 763 764 ktr12 = 1 765 ktr12sq = ktr12 + ntr12am(isymv) 766 kbmatsq = ktr12sq + ntr12sq(isymv) 767 kpck = kbmatsq + nr12r12sq(1) 768 kend1 = kpck + max(nr12r12p(1),ntr12am(isymv)) 769 770 if (lcont) then 771 kpckv = kend1 772 kend1 = kpckv + ntr12am(isymv) 773 end if 774 775 lwrk1 = lwork - kend1 776 if (lwrk1 .lt. 0) then 777 call quit('Insufficient work space in ccrhs_ep') 778 end if 779 780c --------------- 781c get V (= Omega) 782c --------------- 783 !apply 0.5*P_kl^ij: 784 call cc_r12pklij(vijkl,isymv,'T',work(kend1),lwrk1) 785 call dscal(ntr12sq(isymv),0.5D0,vijkl,1) 786 if (lcont) then 787 iopt = 1 788 call ccr12pck2(work(kpckv),isymv,.false.,vijkl,'T',iopt) 789 end if 790 if (locdbg) then 791 call around('OmegaR12 on entry (after 0.5*P_kl^ij):') 792 call cc_prsqr12(vijkl,isymv,'T',1,.false.) 793 write(lupri,*) 'ccrhs_ep> norm^2(OmegaR12):', 794 & ddot(ntr12sq(isymv),vijkl,1,vijkl,1) 795 end if 796 797c ----------------------------------- 798c read R12 amplitudes or trial vector 799c ----------------------------------- 800 call cc_r12getct(work(ktr12sq),isymv,iamp,basscl2,.false.,'T', 801 & lufc12,fc12am,ifile,cdum,idum,work(kend1),lwrk1) 802 if (locdbg) then 803 call around('R12 amplitudes') 804 call cc_prsqr12(work(ktr12sq),isymv,'T',1,.false.) 805 end if 806 807 if (lcont) then 808 iopt = 1 809 call ccr12pck2(work(ktr12),isymv,.false.,work(ktr12sq),'T', 810 & iopt) 811 cont(1) = ddot(ntr12am(isymv),work(ktr12),1,work(kpckv),1) 812 end if 813 814c ---------------- 815c read B matrices 816c ---------------- 817 lunit = -1 818 call gpopen(lunit,fccr12b,'old',' ','unformatted', 819 & idum,ldum) 820 8888 read(lunit) ian, iap, aprox 821 read(lunit) (work(kpck+i), i=0, nr12r12p(1)-1 ) 822 if ((ian.ne.ianr12).or.(iap.ne.iapr12)) goto 8888 823 call gpclose(lunit,'KEEP') 824 if (aproxr12(1:3).ne.aprox) then 825 write(lupri,*)'aproxr12 =', aproxr12(1:3) 826 write(lupri,*)'aprox =', aprox(1:3) 827 call quit('--- Warning: inconsistent R12 approximation') 828 end if 829 iopt = 2 830 call ccr12unpck2(work(kpck),1,work(kbmatsq),'N',iopt) 831 call dscal(nr12r12sq(1),-one,work(kbmatsq),1) 832 833 if (locdbg) then 834 call around('B-Matrix read from file') 835 do isymkl = 1, nsym 836 isymmn = isymkl 837 call output(work(kbmatsq+ir12r12sq(isymkl,isymmn)),1, 838 & nmatkl(isymkl),1,nmatkl(isymmn),nmatkl(isymkl), 839 & nmatkl(isymmn),1,lupri) 840 end do 841 end if 842 843c ---------------------------------------------------- 844c for some cases we are nearly done, since B does not 845c depend on (ij) 846c ---------------------------------------------------- 847 if (ianr12.eq.1) then 848 if (iapr12.le.2 .or. iapr12.eq.4 .or. iapr12.eq.6) then 849 lbmat = .false. 850 else 851 lbmat = .true. 852 factor = half 853 end if 854 else 855 if (iapr12.eq.2 .or. iapr12.eq.5) then 856 lbmat = .false. 857 else 858 lbmat = .true. 859 factor = half 860 end if 861 end if 862 863 if (.not. lbmat) then 864 !matrix multiplication: 865 call cc_r12xi2b(vijkl,'T',work(kbmatsq),1,'N',work(ktr12sq), 866 & isymv,'T',one) 867 else 868c ----------------------------------------------- 869c in other cases: 870c - read R12 overlap matrix, orbital energies 871c - loop over pairs (ij) and add B*c to OmegaR12 872c ----------------------------------------------- 873 !calculate # of indices (kl) / (ij) and offset over all symmetries: 874 nkl = 0 875 nij = 0 876 do isym = 1, nsym 877 inmatkl(isym) = nkl 878 inmatij(isym) = nij 879 nkl = nkl + nmatkl(isym) 880 nij = nij + nmatij(isym) 881 end do 882c 883 kbscr = kend1 884 kxintsq = kbscr + nr12r12sq(1) 885c kevl = kxintsq + nr12r12sq(1) 886c kekl = kevl + norbts 887 kekl = kxintsq + nr12r12sq(1) 888 keij = kekl + nkl 889 keij = kekl + nkl 890 kend2 = keij + nij 891 lwrk2 = lwork - kend2 892 if (lwrk2 .lt. 0) then 893 call quit('Insufficient work space in ccrhs_ep') 894 end if 895c 896 lunit = -1 897 call gpopen(lunit,fccr12x,'old',' ','unformatted', 898 & idum,ldum) 899 9999 read(lunit) ian 900 read(lunit) (work(kpck+i), i=0, nr12r12p(1)-1 ) 901 if (ian.ne.ianr12) goto 9999 902 call gpclose(lunit,'KEEP') 903 iopt = 2 904 call ccr12unpck2(work(kpck),1,work(kxintsq),'N',iopt) 905 if (locdbg) then 906 call around('R12 Overlap-Matrix read from file') 907 do isymkl = 1, nsym 908 isymmn = isymkl 909 call output(work(kxintsq+ir12r12sq(isymkl,isymmn)),1, 910 & nmatkl(isymkl),1,nmatkl(isymmn),nmatkl(isymkl), 911 & nmatkl(isymmn),1,lupri) 912 end do 913 end if 914c 915c ------------------------------------ 916c Read orbital energies. 917c ------------------------------------- 918c 919 lunit = -1 920 call gpopen(lunit,'SIRIFC','OLD',' ','UNFORMATTED',idummy, 921 & .false.) 922 rewind lunit 923 call mollab('FULLBAS ',lunit,lupri) 924 read (lunit) idum,norbtsx 925 kevl = kend2 926 kend2 = kevl + norbtsx 927 lwrk2 = lwork - kend2 928 if (lwrk2 .lt. 0) then 929 call quit('Insufficient work space in ccrhs_ep') 930 end if 931 read (lunit) (work(kevl+i-1), i=1,norbtsx) 932 call gpclose(lunit,'KEEP') 933c 934 call cc_r12epair(work(kekl),nkl,work(kevl),inmatkl,imatkl,nrhfb) 935 call cc_r12epair(work(keij),nij,work(kevl),inmatij,imatij,nrhf) 936c 937 do isymij = 1, nsym 938 isymmn = muld2h(isymv,isymij) 939 isymkl = isymmn 940c if (locdbg) then 941c write(lupri,*) 'in loop over idxij: isymij = ',isymij 942c end if 943 do idxij = 1, nmatij(isymij) 944 koffb = ir12r12sq(isymkl,isymmn) 945 call cc_r12buildbmat(work(kbscr+koffb),isymkl,isymmn, 946 & isymij,idxij,work(kbmatsq+koffb), 947 & work(kxintsq+koffb),work(kekl),inmatkl, 948 & work(keij),inmatij,factor) 949 kofftr12 = ktr12sq + itr12sqt(isymij,isymmn) + idxij-1 950 koffv = itr12sqt(isymij,isymkl) + idxij 951 if (locdbg) then 952 write(lupri,*) 'idxij = ',idxij 953 write(lupri,*) 'R12 amplitudes part used:' 954 do i = 1, nmatkl(isymmn) 955 write(lupri,*) kofftr12+(i-1)*nmatij(isymij), 956 & work(kofftr12+(i-1)*nmatij(isymij)) 957 end do 958 write(lupri,*) 'Omega part used:' 959 do i = 1, nmatkl(isymkl) 960 write(lupri,*) koffv+(i-1)*nmatij(isymij), 961 & vijkl(koffv+(i-1)*nmatij(isymij)) 962 end do 963 end if 964 call dgemv('N',nmatkl(isymkl),nmatkl(isymmn), 965 & one,work(kbscr+koffb),max(nmatkl(isymkl),1), 966 & work(kofftr12),max(nmatij(isymij),1), 967 & one,vijkl(koffv),max(nmatij(isymij),1)) 968 end do 969 end do 970 end if 971 972c -------------- 973c print results: 974c -------------- 975 if (locdbg) then 976 call around('OmegaR12 on exit') 977 call cc_prsqr12(vijkl,isymv,'T',1,.false.) 978 write(lupri,*) 'ccrhs_ep> norm^2(OmegaR12):', 979 & ddot(ntr12sq(isymv),vijkl,1,vijkl,1) 980 end if 981 982 if (iamp .EQ. 1) then 983c ------------------------------------------------------------ 984c for jacobian right transformations scale diagonal to 985c transform to pseudo-orthonormal basis 986c ------------------------------------------------------------ 987 iopt = 1 988 call ccr12pck2(work(kpck),isymv,.false.,vijkl,'T',iopt) 989 call cclr_diasclr12(work(kpck),basscl1,isymv) 990 call cc_wvec(lufr12,frho12,ntr12am(isymv),ntr12am(isymv), 991 * ifile,work(kpck)) 992 if (locdbg) then 993 call around('Jacobian transformed vector in CCRHS_EP') 994 call cc_prpr12(work(kpck),isymv,1,.false.) 995 end if 996 else 997c ------------------------------------------------------------ 998c write result vector as ground state vector function to file: 999c ------------------------------------------------------------ 1000 iopt = 1 1001 call ccr12pck2(work(ktr12),isymv,.false.,vijkl,'T',iopt) 1002 lunit = -1 1003 call gpopen(lunit,'CC_OMEGAR12','unknown',' ','unformatted', 1004 & idum,ldum) 1005 write(lunit) (work(ktr12+i), i=0, ntr12am(isymv)-1 ) 1006 call gpclose(lunit,'KEEP') 1007 end if 1008c 1009 if (lcont) then 1010 call daxpy(ntr12am(isymv),-one,work(kpckv),1,work(kpck),1) 1011 cont(2) = ddot(ntr12am(isymv),work(ktr12),1,work(kpck),1) 1012 end if 1013 1014 if (locdbg) write(lupri,*) 'Leaving CCRHS_EP' 1015 1016 call qexit('ccrhs_ep') 1017 end 1018*=====================================================================* 1019 1020*=====================================================================* 1021 subroutine cc_r12nxtam_old(omeg12,isym,tamp12,lcceq, 1022 & er12,work,lwork) 1023c---------------------------------------------------------------------- 1024c purpose: get new start R12 amplitudes for each iteration 1025c 1026c H. Fliegl, C. Haettig, (WK/UniKA/08-05-2003) 1027c modified by C. Neiss April 2005, September 2005 1028c---------------------------------------------------------------------- 1029 implicit none 1030#include "priunit.h" 1031#include "ccorb.h" 1032#include "ccsdsym.h" 1033#include "dummy.h" 1034#include "r12int.h" 1035#include "ccr12int.h" 1036#include "ccsdinp.h" 1037#include "second.h" 1038 1039 logical locdbg,ldum,lcceq 1040 parameter (locdbg = .false.) 1041 1042 integer kbsing,kbtrip,lwork,nrhftria,n2,kxsing,kxtrip,kevl 1043 integer kend1,lwrk1,lunit,idum,ijv,ij,kvsing,kvtrip 1044 integer komgsing,komgtrip,kcsing,kctrip 1045 integer kscr,kvint,lu43,kepsij,kepsoff 1046 integer kbb,kuu,kpp,kqq,krr,ian,iap,iopt,isym 1047 integer nkilj(8),nij,isymij,inmatij(8),it,jt,isymi,isymj,norbtsx 1048#if defined (SYS_CRAY) 1049 real tamp12(*),omeg12(*),er12,work(*) 1050 real ensing,entrip,delta,xf,fs,ft,ws,wt,ddot,timnxtam 1051#else 1052 double precision tamp12(*),omeg12(*),er12,work(*) 1053 double precision ensing,entrip,delta,xf,fs,ft,ws,wt,ddot,timnxtam 1054#endif 1055 character*3 aprox 1056 character*10 model 1057 integer index 1058c 1059 index(i,j) = max(i,j)*(max(i,j)-3)/2 + i + j 1060c 1061 call qenter('r12nxtam') 1062 if (locdbg) then 1063 write(lupri,*) 'Entered CC_R12NXTAM' 1064 call flshfo(lupri) 1065 end if 1066c 1067 timnxtam = second() 1068c 1069 nrhftria = nrhftb*(nrhftb+1)/2 1070 n2 = nrhftria * nrhftria 1071c 1072 nij = 0 1073 do isymij = 1, nsym 1074 inmatij(isymij) = nij 1075 nij = nij + nmatij(isymij) 1076 end do 1077c 1078 kend1 = 1 1079c 1080 komgsing = kend1 1081 komgtrip = komgsing + n2 1082 kcsing = komgtrip + n2 1083 kctrip = kcsing + n2 1084 kbsing = kctrip + n2 1085 kbtrip = kbsing + n2 1086 kxsing = kbtrip + n2 1087 kxtrip = kxsing + n2 1088 kend1 = kxtrip + n2 1089 1090 if (lcceq) then 1091 kvsing = kend1 1092 kvtrip = kvsing + n2 1093 kend1 = kvtrip + n2 1094 end if 1095 1096 kevl = kend1 1097 kbb = kevl + nrhftria 1098 kuu = kbb + n2 1099 kpp = kuu + n2 1100 kqq = kpp + n2 1101 krr = kqq + n2 1102 kend1 = krr + n2 1103c 1104 kepsij = kend1 1105 kend1 = kepsij + nij 1106c 1107 kscr = kend1 1108 kend1 = kscr + nr12r12p(1) 1109c 1110 lwrk1 = lwork - kend1 1111 1112 if (lwrk1 .lt. 0) then 1113 call quit('Insufficient work space in cc_r12nxtam') 1114 end if 1115 1116c ------------- 1117c test symmetry 1118c ------------- 1119 if (lcceq .and. (isym.ne.1)) then 1120 call quit('Symmetry error in CC_R12NXTAM') 1121 end if 1122 1123c ---------------------- 1124c read orbital energies 1125c --------------------- 1126 lunit = -1 1127 ldum = .false. 1128 call gpopen(lunit,fccr12e,'old',' ','formatted', 1129 & idum,ldum) 1130 read(lunit,'(4e30.20)') (work(kevl+ij), ij = 0, nrhftria-1) 1131 call gpclose(lunit,'KEEP') 1132 if (locdbg) then 1133 write(lupri,*) 'Read orbital energies done:' 1134 call outpak(work(kevl),nrhftria,1,lupri) 1135 call flshfo(lupri) 1136 end if 1137 1138 if (lcceq) then 1139c --------------- 1140c read V matrices 1141c --------------- 1142 lunit = -1 1143 call gpopen(lunit,fccr12v,'old',' ','unformatted', 1144 & idum,ldum) 1145 7777 read(lunit) ian 1146 read(lunit) (work(kscr+i), i=0, ntr12am(1)-1) 1147 if (ian.ne.ianr12) goto 7777 1148 call gpclose(lunit,'KEEP') 1149 if (locdbg) then 1150 write(lupri,*) 'Read V-intermediate done' 1151 call flshfo(lupri) 1152 end if 1153 call ccr12unpck(work(kscr),1,work(kvsing),work(kvtrip),nrhfb, 1154 & nrhf) 1155 if (locdbg) then 1156 write(lupri,*) 'Unpack V-intermediate done' 1157 write(lupri,*) 'Singlet V-intermediate:' 1158 call output(work(kvsing),1,nrhftria,1,nrhftria,nrhftria, 1159 & nrhftria,1,lupri) 1160 write(lupri,*) 'Triplet V-intermediate:' 1161 call output(work(kvtrip),1,nrhftria,1,nrhftria,nrhftria, 1162 & nrhftria,1,lupri) 1163 call flshfo(lupri) 1164 end if 1165CCN 1166C write(lupri,*) 'V-Matrix in CC_R12NXTAM:' 1167C call ccr12unpck2(work(kscr),1,work(kend1),'T',1) 1168CCN 1169 end if 1170 1171c --------------- 1172c read B matrices 1173c --------------- 1174 lunit = -1 1175 call gpopen(lunit,fccr12b,'old',' ','unformatted', 1176 & idum,ldum) 1177 8888 read(lunit) ian, iap, aprox 1178 read(lunit) (work(kscr+i), i=0, nr12r12p(1)-1) 1179 if ((ian.ne.ianr12).or.(iap.ne.iapr12)) goto 8888 1180 call gpclose(lunit,'KEEP') 1181 call ccr12unpck(work(kscr),1,work(kbsing),work(kbtrip),nrhfb, 1182 & nrhfb) 1183 1184c ---------------------------- 1185c read R12 overlap matrices X 1186c ---------------------------- 1187 lunit = -1 1188 call gpopen(lunit,fccr12x,'old',' ','unformatted', 1189 & idum,ldum) 1190 9999 read(lunit) ian 1191 read(lunit) (work(kscr+i), i=0, nr12r12p(1)-1) 1192 if (ian.ne.ianr12) goto 9999 1193 call gpclose(lunit,'KEEP') 1194 call ccr12unpck(work(kscr),1,work(kxsing),work(kxtrip),nrhfb, 1195 & nrhfb) 1196 1197 if (lcceq) then 1198c -------------------- 1199c read R12 amplitudes 1200c -------------------- 1201 iopt =32 1202 call cc_rdrsp('R0 ',0,1,iopt,model,dummy,tamp12) 1203 1204 if (iprint .ge. 5) then 1205 write(lupri,529) 'Norm^2 of r12am in NXTAM:', 1206 * ddot(ntr12am(1),tamp12,1,tamp12,1) 1207 end if 1208 1209c ------------------------------- 1210c read vector function omega_R12 1211c ------------------------------- 1212 lunit = -1 1213 call gpopen(lunit,'CC_OMEGAR12','old',' ','unformatted', 1214 & idum,ldum) 1215 read(lunit) (omeg12(i), i=1, ntr12am(1)) 1216 call gpclose(lunit,'KEEP') 1217chf 1218 write(lupri,529)'Norm^2 of Omegar12 in NXTAM:', 1219 & ddot(ntr12am(1),omeg12,1,omeg12,1) 1220chf 1221 end if 1222 !repack omega_R12 1223 call ccr12unpck(omeg12,isym,work(komgsing),work(komgtrip),nrhfb, 1224 & nrhf) 1225 1226 if (locdbg) then 1227 write(lupri,*) 'Singlet Omega-R12:' 1228 call output(work(komgsing),1,nrhftria,1,nrhftria,nrhftria, 1229 & nrhftria,1,lupri) 1230 write(lupri,*) 'Triplet Omega-R12:' 1231 call output(work(komgtrip),1,nrhftria,1,nrhftria,nrhftria, 1232 & nrhftria,1,lupri) 1233 call around('b in nxtam matrix (singlet)') 1234 call output(work(kbsing),1,nrhftria,1,nrhftria, 1235 & nrhftria,nrhftria,1,lupri) 1236 call around('b in nxtam matrix (triplet)') 1237 call output(work(kbtrip),1,nrhftria,1,nrhftria, 1238 & nrhftria,nrhftria,1,lupri) 1239 call around('x in nxtam matrix (singlet)') 1240 call output(work(kxsing),1,nrhftria,1,nrhftria, 1241 & nrhftria,nrhftria,1,lupri) 1242 call around('x in nxtam matrix (triplet)') 1243 call output(work(kxtrip),1,nrhftria,1,nrhftria, 1244 & nrhftria,nrhftria,1,lupri) 1245 call flshfo(lupri) 1246 end if 1247 1248 529 format(7x,a,d24.10) 1249 1250 if (locdbg) then 1251 write(lupri,529) 'norm^2 of omega r12 singlet:', 1252 * ddot(n2,work(komgsing),1,work(komgsing),1) 1253 write(lupri,529) 'norm^2 of omega r12 triplet:', 1254 * ddot(n2,work(komgtrip),1,work(komgtrip),1) 1255 end if 1256 1257 call dzero(work(kcsing),n2) 1258 call dzero(work(kctrip),n2) 1259 1260 if (locdbg) write(lupri,*) 'R12SVD,R12DIA:',R12SVD,R12DIA 1261c ------------------------------------------------------ 1262c get CC2-R12 contributions to the ground state energy 1263c 1264c xf =0.5d0 for approximation A' 1265c xf = 0 for for approximation A 1266c ------------------------------------------------------ 1267 xf = 0.0d0 1268 delta = 0d0 1269 if (ianr12.eq.1) then 1270 if (iapr12.le.2.or.iapr12.eq.4.or.iapr12.eq.6) then 1271 xf = 0 1272 else 1273 xf = 0.5d0 1274 end if 1275 else if (ianr12.eq.2 .or. ianr12.eq.3) then 1276c if (iapr12.le.2) then 1277 if (iapr12.eq.2.or.iapr12.eq.5) then 1278 xf = 0 1279 else 1280 xf = 0.5d0 1281 end if 1282 endif 1283c ----------------------------------------------- 1284c Get orbital energy pairs for OCCUPIED orbitals: 1285c ----------------------------------------------- 1286 lunit = -1 1287 call gpopen(lunit,'SIRIFC','OLD',' ','UNFORMATTED',idum, 1288 & .false.) 1289 rewind lunit 1290 call mollab('FULLBAS ',lunit,lupri) 1291 read (lunit) idum,norbtsx 1292 kend1 = kscr + norbtsx 1293 lwrk1 = lwork - kend1 1294 if (lwrk1.lt.0) then 1295 call quit('Memory allocation error in CC_R12NXTAM') 1296 end if 1297 read (lunit) (work(kscr+i-1), i=1, norbtsx) 1298 call gpclose(lunit,'KEEP') 1299c 1300 call cc_r12epair(work(kepsij),nij,work(kscr),inmatij,imatij,nrhf) 1301c 1302 do isymj = 1, nsym 1303 do isymi = 1, nsym 1304 do j = 1, nrhf(isymj) 1305 jt = irhf(isymj) + j 1306 do i = 1, nrhf(isymi) 1307 it = irhf(isymi) + i 1308 if (it.le.jt) then 1309 ij = index(it,jt) 1310 ijv = nrhftria*(ij-1) 1311 isymij = muld2h(isymi,isymj) 1312 kepsoff = inmatij(isymij)+imatij(isymi,isymj)+ 1313 & nrhf(isymi)*(j-1)+i-1 1314 if (R12SVD) then 1315 CALL R12MP2(work(komgsing+ijv),work(komgtrip+ijv), 1316 & work(komgsing+ijv),work(komgtrip+ijv), 1317 & work(kxsing),work(kxtrip), 1318 & work(kepsij+kepsoff),nrhftb,nrhftria, 1319 & fs,ft,work(kbb),work(kuu),work(kpp), 1320 & work(kqq), 1321 & work(kbsing),work(kbtrip), 1322 & work(kevl),xf,work(krr),ij,ws,wt, 1323 & delta,work(kcsing+ijv),work(kctrip+ijv)) 1324 else if (R12DIA) then 1325 CALL MP2R12(work(komgsing+ijv),work(komgtrip+ijv), 1326 & work(komgsing+ijv),work(komgtrip+ijv), 1327 & work(kxsing),work(kxtrip), 1328 & work(kepsij+kepsoff),nrhftb,nrhftria, 1329 & fs,ft,work(kbb),work(kuu),work(kpp), 1330 & work(kqq), 1331 & work(kbsing),work(kbtrip), 1332 & work(kevl),xf,work(krr),ij,ws,wt, 1333 & delta,work(kcsing+ijv),work(kctrip+ijv)) 1334 end if 1335 end if 1336 end do 1337 end do 1338 end do 1339 end do 1340 1341c ------------------------------------------------ 1342c copy update for amplitudes into omgsing/omgtrip 1343c ------------------------------------------------ 1344 call dcopy(n2,work(kcsing),1,work(komgsing),1) 1345 call dcopy(n2,work(kctrip),1,work(komgtrip),1) 1346 if (locdbg) then 1347 write(lupri,*) 'Update for Singlet R12 amplitudes:' 1348 call output(work(kcsing),1,nrhftria,1,nrhftria,nrhftria, 1349 & nrhftria,1,lupri) 1350 write(lupri,*) 'Update for Triplet R12 amplitudes:' 1351 call output(work(kctrip),1,nrhftria,1,nrhftria,nrhftria, 1352 & nrhftria,1,lupri) 1353 end if 1354 1355 if (lcceq) then 1356 call ccr12unpck(tamp12,1,work(kcsing),work(kctrip),nrhfb, 1357 & nrhf) 1358 1359 if (locdbg) then 1360 write(lupri,*) 'Singlet R12 amplitudes:' 1361 call output(work(kcsing),1,nrhftria,1,nrhftria,nrhftria, 1362 & nrhftria,1,lupri) 1363 write(lupri,*) 'Triplet R12 amplitudes:' 1364 call output(work(kctrip),1,nrhftria,1,nrhftria,nrhftria, 1365 & nrhftria,1,lupri) 1366 end if 1367 1368c ---------------------------------------------------- 1369c add old amplitudes to update to get new amplitudes 1370c ---------------------------------------------------- 1371 call daxpy(n2,1.0d0,work(kcsing),1,work(komgsing),1) 1372 call daxpy(n2,1.0d0,work(kctrip),1,work(komgtrip),1) 1373 1374 ensing = ddot(n2,work(kvsing),1,work(komgsing),1) 1375 entrip = 3.0D0*ddot(n2,work(kvtrip),1,work(komgtrip),1) 1376 er12 = ensing + entrip 1377 1378 if (locdbg) then 1379 write(lupri,*) 'Singlet contribution to R12 energy:',ensing 1380 write(lupri,*) 'Triplet contribution to R12 energy:',entrip 1381 end if 1382 1383 end if ! lcceq 1384 1385 if (locdbg) then 1386 write(lupri,*) 'Updated Singlet R12 amplitudes:' 1387 call output(work(komgsing),1,nrhftria,1,nrhftria,nrhftria, 1388 & nrhftria,1,lupri) 1389 write(lupri,*) 'Updated Triplet R12 amplitudes:' 1390 call output(work(komgtrip),1,nrhftria,1,nrhftria,nrhftria, 1391 & nrhftria,1,lupri) 1392 end if 1393 1394 call ccr12pck(omeg12,isym,work(komgsing),work(komgtrip),nrhfb, 1395 & nrhf,nkilj) 1396 1397 if (locdbg) write(lupri,*) 'Leaving CC_R12NXTAM' 1398c 1399 timnxtam = second() - timnxtam 1400c write(lupri,111) 'CC_R12NXTAM', timnxtam 1401 111 FORMAT(7x,'Time used in',2x,A12,2x,': ',f10.2,' seconds') 1402c 1403 call qexit('r12nxtam') 1404 end 1405*=====================================================================* 1406 subroutine cc_r12metric_old(isymr,basscl1,basscl2,work,lwork, 1407 & fc2am,lufc2,fc12am,lufc12,fs12am,lufs12, 1408 & fs2am,lufs2,ifile) 1409*---------------------------------------------------------------------* 1410c purpose: make the CC2-R12 S*R contibution for excitation 1411c energies 1412c 1413c H. Fliegl, C. Haettig summer 2003 1414c 1415c nondiagonal elements added for ansatz 2, spring 2005 1416c --> S*R written on fs2am 1417c 1418c C. Neiss, 05.04.2005: adapted for left hand transformation 1419c basscl1 = brascl for right hand transf. 1420c = 0.5*ketscl for left hand transf. 1421c basscl2 = ketscl for right hand transf. 1422c = 2*brascl for left hand transf. 1423c--------------------------------------------------------------------- 1424 implicit none 1425#include "ccsdsym.h" 1426#include "priunit.h" 1427#include "ccsdinp.h" 1428#include "ccorb.h" 1429#include "r12int.h" 1430#include "ccr12int.h" 1431 1432 logical ldum,locdbg,lprojv 1433 parameter (locdbg = .false.) 1434 1435 character*(*) fc12am,fs12am,fc2am,fs2am 1436 1437 integer lwork,isymr,nrhftria,kend1,lwrk1,lufc12,lufs12,ifile 1438 integer lufc2,lufs2 1439 integer krsing,krtrip,kxsing,kxtrip,ij,ksing,ktrip,kpck,kxint 1440 integer idum,n2,lunit,ijrsvc,ijrtvc,ijsvc,ijtvc 1441 integer ian, iap,ks2,krabkl,krmnab,ksr 1442 integer nkilj(8) 1443#if defined (SYS_CRAY) 1444 real work(*), ddot,two,one,half,basscl1,basscl2 1445#else 1446 double precision work(*), ddot,two,one,half,basscl1,basscl2 1447#endif 1448 parameter (one = 1.0d0, two = 2.0d0, half = 0.5d0) 1449 1450 call qenter('cc_r12metric') 1451 if (locdbg) then 1452 write(lupri,*) 'Entered CC_R12METRIC' 1453 end if 1454 1455 nrhftria = nrhftb*(nrhftb+1)/2 1456 n2 = nrhftria * nrhftria 1457 1458 ksing = 1 1459 ktrip = ksing + n2 1460 kxsing = ktrip + n2 1461 kxtrip = kxsing + n2 1462 krsing = kxtrip + n2 1463 krtrip = krsing + n2 1464 kxint = krtrip + n2 1465 kend1 = kxint + nr12r12p(1) 1466 1467 kpck = kend1 1468 kend1 = kpck + ntr12am(isymr) 1469 if (ianr12.eq.2) then 1470 krabkl = kend1 1471 krmnab = krabkl + nt2r12(1) 1472 ksr = krmnab + nt2am(isymr) 1473 kend1 = ksr + ntr12sq(isymr) 1474 end if 1475 1476 lwrk1 = lwork - kend1 1477 1478 if (lwrk1 .lt. 0) then 1479 write(lupri,*) 'lwork, lwrk1: ',lwork, lwrk1 1480 call quit('Insufficient work space in cc_r12metric') 1481 end if 1482 1483 call dzero(work(krsing),n2) 1484 call dzero(work(krtrip),n2) 1485 1486 call dzero(work(krsing),n2) 1487 call dzero(work(krtrip),n2) 1488 1489c ------------------------- 1490c read R12 overlap matrices 1491c ------------------------- 1492 lunit = -1 1493 call gpopen(lunit,fccr12x,'unknown',' ','unformatted', 1494 & idum,ldum) 1495 9999 read(lunit) ian 1496 read(lunit) (work(kxint+i), i=0, nr12r12p(1)-1 ) 1497 if (ian.ne.ianr12) goto 9999 1498 call gpclose(lunit,'KEEP') 1499 call ccr12unpck(work(kxint),1,work(kxsing),work(kxtrip),nrhfb, 1500 & nrhfb) 1501 1502c ------------------------ 1503c read try vector R^ij_kl 1504c ------------------------ 1505 lunit = -1 1506 call cc_rvec(lufc12,fc12am,ntr12am(isymr),ntr12am(isymr), 1507 * ifile,work(kpck)) 1508 call cclr_diasclr12(work(kpck),basscl2,isymr) 1509 call ccr12unpck(work(kpck),isymr,work(krsing),work(krtrip), 1510 & nrhfb,nrhf) 1511 1512 if (locdbg) then 1513 write(lupri,*) 'fc12am,lufc12,ifile:',fc12am,lufc12,ifile 1514 write(lupri,*) 'cc_r12metric> norm^2(R12 packed):', 1515 & ddot(ntr12am(isymr),work(kpck),1,work(kpck),1) 1516 write(lupri,*) 'cc_r12metric> norm^2(R12 singlet):', 1517 & ddot(n2,work(krsing),1,work(krsing),1) 1518 write(lupri,*) 'cc_r12metric> norm^2(R12 triplet):', 1519 & ddot(n2,work(krtrip),1,work(krtrip),1) 1520 end if 1521 if (locdbg) then 1522 call around('r12 try amplitudes (singlet)') 1523 call output(work(krsing),1,nrhftria,1,nrhftria, 1524 & nrhftria,nrhftria,1,lupri) 1525 call around('r12 try amplitudes (triplet)') 1526 call output(work(krtrip),1,nrhftria,1,nrhftria, 1527 & nrhftria,nrhftria,1,lupri) 1528 call around('x matrix (singlet)') 1529 call output(work(kxsing),1,nrhftria,1,nrhftria, 1530 & nrhftria,nrhftria,1,lupri) 1531 call around('x matrix (triplet)') 1532 call output(work(kxtrip),1,nrhftria,1,nrhftria, 1533 & nrhftria,nrhftria,1,lupri) 1534 end if 1535 1536c ------------------------------------------------- 1537c loop over pairs and make S*R for s = 0,1 1538c ------------------------------------------------- 1539 call dzero(work(ksing),n2) 1540 call dzero(work(ktrip),n2) 1541 1542 if (ianr12.eq.2) then 1543c read r^ab_kl integrals 1544 lunit = -1 1545 call gpopen(lunit,fr12r12,'unknown',' ','unformatted', 1546 & idum,.false.) 1547 read(lunit)(work(krabkl-1+i),i=1,nt2r12(1)) 1548 call gpclose(lunit,'KEEP') 1549c read R^mn_ab trial vector (stored as triangular matrix) 1550 call cc_rvec(lufc2,fc2am,nt2am(isymr),nt2am(isymr), 1551 & ifile,work(krmnab)) 1552 call cclr_diascl(work(krmnab),two,isymr) 1553c calculate \sum_ab r^ab_kl * R^mn_ab 1554 call dzero(work(ksr),ntr12sq(isymr)) 1555 call cc_r12mi2(work(ksr),work(krabkl),work(krmnab),1,isymr,one, 1556 & work(kend1),lwrk1) 1557c get singlet and triplet sr contribution 1558c and write it on work(ksing) and work(ktrip) 1559 lprojv = .true. 1560 call cc_r12vunpack(work(ksr),isymr,work(ksing),work(ktrip), 1561 & lprojv,nrhfb,nrhf) 1562 end if 1563 1564 ijrsvc = krsing 1565 ijrtvc = krtrip 1566 ijsvc = ksing 1567 ijtvc = ktrip 1568 do ij = 1, nrhftria 1569 call cc_r12metric0(ij,nrhftria, 1570 & work(ijsvc),work(ijtvc), 1571 & work(ijrsvc),work(ijrtvc), 1572 & work(kxsing),work(kxtrip)) 1573 ijrsvc = ijrsvc + nrhftria 1574 ijrtvc = ijrtvc + nrhftria 1575 ijsvc = ijsvc + nrhftria 1576 ijtvc = ijtvc + nrhftria 1577 end do 1578c ----------------------- 1579c print and pack results 1580c ----------------------- 1581 if (locdbg) then 1582 call around('r12 S*R (singlet)') 1583 call output(work(ksing),1,nrhftria,1,nrhftria, 1584 * nrhftria,nrhftria,1,lupri) 1585 call around('r12 S*R (triplet)') 1586 call output(work(ktrip),1,nrhftria,1,nrhftria, 1587 * nrhftria,nrhftria,1,lupri) 1588 end if 1589 if ( locdbg ) then 1590 write(lupri,*) 'cc_r12metric> norm^2(S * R12 singlet):', 1591 & ddot(n2,work(ksing),1,work(ksing),1) 1592 write(lupri,*) 'cc_r12metric> norm^2(S * R12 triplet):', 1593 & ddot(n2,work(ktrip),1,work(ktrip),1) 1594 end if 1595 1596c ------------------ 1597c write S*R on file 1598c ------------------ 1599 call ccr12pck(work(kpck),isymr,work(ksing),work(ktrip),nrhfb, 1600 & nrhf,nkilj) 1601 call cclr_diasclr12(work(kpck),basscl1,isymr) 1602 if (locdbg) then 1603 call around('R12 S*R after packing') 1604 call cc_prpr12(work(kpck),isymr,1,.false.) 1605 end if 1606 call cc_wvec(lufs12,fs12am,ntr12am(isymr),ntr12am(isymr), 1607 * ifile,work(kpck)) 1608c 1609 if (ianr12.eq.2) then 1610c get \sum_kl r^kl_ab * R^ij_kl + R^ij_ab and put it on file 1611 call cc_r12getsr22(fc2am,lufc2,fc12am,lufc12,fs2am, 1612 & lufs2,isymr,ifile,work(kend1),lwrk1) 1613 end if 1614c 1615 if (locdbg.and.ianr12.eq.2) then 1616 call cc_rvec(lufs12,fs12am,ntr12am(isymr),ntr12am(isymr),ifile, 1617 & work(kpck)) 1618 write(lupri,*)'FS12AM' 1619 call cc_prpr12(work(kpck),isymr,1,.false.) 1620 end if 1621 1622 if (locdbg) write(lupri,*) 'Leaving CC_R12METRIC' 1623 call qexit('cc_r12metric') 1624 end 1625*=====================================================================* 1626 subroutine cc_r12metric0(ij,n,srsing,srtrip, 1627 & rsing,rtrip,xsing,xtrip) 1628*---------------------------------------------------------------------* 1629c purpose: make the CC2-R12 S*R part for s = 0,1 1630c 1631c H. Fliegl, C. Haettig summer 2003 1632c--------------------------------------------------------------------- 1633 implicit none 1634#include "priunit.h" 1635 integer n,ij,kl,mn 1636#if defined (SYS_CRAY) 1637 real rsing(n),rtrip(n),xsing(n,n),xtrip(n,n),srsing(n),srtrip(n) 1638 & 1639#else 1640 double precision rsing(n),rtrip(n),xsing(n,n), 1641 & xtrip(n,n),srsing(n),srtrip(n) 1642#endif 1643 call qenter('cc_r12metric0') 1644 1645 do kl = 1, n 1646 do mn = 1, n 1647 srsing(kl) = srsing(kl)+ xsing(kl,mn)*rsing(mn) 1648 srtrip(kl) = srtrip(kl)+ xtrip(kl,mn)*rtrip(mn) 1649 end do 1650 end do 1651 1652 call qexit('cc_r12metric0') 1653 end 1654*=====================================================================* 1655 subroutine cc_r12metric(isymr,basscl1,basscl2,work,lwork, 1656 & fc2am,lufc2,fc12am,lufc12,fs12am,lufs12, 1657 & fs2am,lufs2,ifile,lnoread,vec) 1658*---------------------------------------------------------------------* 1659c purpose: make the CC-R12 S*R contibution for excitation 1660c energies 1661c alternative subroutine for cc_r12metric_old without 1662c using singlet/triplet matrices in matrix 1663c multiplication 1664c 1665c Christian Neiss, Feb. 2005, based on old cc_r12metric 1666c 1667c H. Fliegl, spring 2005: nondiagonal elements added for ansatz 2 1668c --> S*R written on fs2am 1669c 1670c C. Neiss, 05.04.2005: adapted for left hand transformation 1671c basscl1 = brascl for right hand transf. 1672c = 0.5*ketscl for left hand transf. 1673c basscl2 = ketscl for right hand transf. 1674c = 2*brascl for left hand transf. 1675c 1676c C. Neiss, 01.07.2004: possibility to get input vector via 1677c call-list, not via file: 1678c vec input vector 1679c dimension(vec): nt1am(isymr)+nt2am(isymr)+ntr12am(isymr) 1680c--------------------------------------------------------------------- 1681 implicit none 1682#include "ccsdsym.h" 1683#include "priunit.h" 1684#include "ccsdinp.h" 1685#include "ccorb.h" 1686#include "r12int.h" 1687#include "ccr12int.h" 1688 1689 logical locdbg 1690 parameter (locdbg=.FALSE.) 1691 1692 logical ldum,lnoread 1693 character*(*) fc12am,fs12am,fc2am,fs2am 1694 integer lwork,isymr,kend1,lwrk1 1695 integer lufc12,lufs12,lufc2,lufs2 1696 integer kxint,ij,kpck,kxintsq,kr12sq,kres,krabkl,krmnab 1697 integer idum,lunit,ifile,ian,iopt 1698#if defined (SYS_CRAY) 1699 real work(*), ddot, one, half, two, basscl1, basscl2, vec(*) 1700#else 1701 double precision work(*), ddot, one, half, two, basscl1, basscl2, 1702 & vec(*) 1703#endif 1704 parameter (one = 1.0D0, half = 0.5D0, two = 2.0D0) 1705 1706 call qenter('cc_r12metric') 1707 if (locdbg) then 1708 write(lupri,*) 'Entered CC_R12METRIC' 1709 end if 1710 1711 kres = 1 1712 kxint = kres + ntr12sq(isymr) 1713 kxintsq = kxint + nr12r12p(1) 1714 kr12sq = kxintsq + nr12r12sq(1) 1715 kpck = kr12sq + ntr12sq(isymr) 1716 kend1 = kpck + ntr12am(isymr) 1717 1718 if (ianr12.eq.2) then 1719 krabkl = kend1 1720 krmnab = krabkl + nt2r12(1) 1721 kend1 = krmnab + nt2am(isymr) 1722 end if 1723 1724 lwrk1 = lwork - kend1 1725 1726 if (lwrk1 .lt. 0) then 1727 write(lupri,*) 'lwork, lwrk1: ',lwork, lwrk1 1728 call quit('Insufficient work space in CC_R12METRIC') 1729 end if 1730 1731c -------------------------------------------------------- 1732c read R12 overlap matrices and pack to symmetrized square 1733c -------------------------------------------------------- 1734 lunit = -1 1735 call gpopen(lunit,fccr12x,'unknown',' ','unformatted', 1736 & idum,ldum) 1737 9999 read(lunit) ian 1738 read(lunit) (work(kxint+i), i=0, nr12r12p(1)-1 ) 1739 if (ian.ne.ianr12) goto 9999 1740 call gpclose(lunit,'KEEP') 1741 iopt = 2 1742 call ccr12unpck2(work(kxint),1,work(kxintsq),'T',iopt) 1743 1744c ------------------------------------------------- 1745c read trial vector and pack to symmetrized square 1746c ------------------------------------------------- 1747 if (.not. lnoread) then 1748 call cc_rvec(lufc12,fc12am,ntr12am(isymr),ntr12am(isymr), 1749 * ifile,work(kpck)) 1750 else 1751 call dcopy(ntr12am(isymr),vec(nt1am(isymr)+nt2am(isymr)+1), 1752 & 1,work(kpck),1) 1753 end if 1754 call cclr_diasclr12(work(kpck),basscl2,isymr) 1755 iopt = 1 1756 call ccr12unpck2(work(kpck),isymr,work(kr12sq),'T',iopt) 1757 1758 if (locdbg) then 1759 write(lupri,*) 'fc12am,lufc12,ifile:',fc12am,lufc12,ifile 1760 write(lupri,*) 'cc_r12metric> norm^2(R12 packed triangle):', 1761 & ddot(ntr12am(isymr),work(kpck),1,work(kpck),1) 1762 write(lupri,*) 'cc_r12metric> norm^2(R12 packed square):', 1763 & ddot(ntr12sq(isymr),work(kr12sq),1,work(kr12sq),1) 1764 end if 1765 if (locdbg) then 1766 call around('R12 trial amplitudes') 1767 call cc_prsqr12(work(kr12sq),isymr,'T',1,.false.) 1768 end if 1769 1770c ------------- 1771c calculate S*R 1772c ------------- 1773 call dzero(work(kres),ntr12sq(isymr)) 1774 1775 ! conv. doubles contributions for Ansatz 2: 1776 if (ianr12.eq.2) then 1777c read r^ab_kl integrals 1778 lunit = -1 1779 call gpopen(lunit,fr12r12,'unknown',' ','unformatted', 1780 & idum,.false.) 1781 read(lunit)(work(krabkl-1+i),i=1,nt2r12(1)) 1782 call gpclose(lunit,'KEEP') 1783c read R^mn_ab trial vector (stored as triangular matrix) 1784 if (.not. lnoread) then 1785 call cc_rvec(lufc2,fc2am,nt2am(isymr),nt2am(isymr), 1786 & ifile,work(krmnab)) 1787 else 1788 call dcopy(nt2am(isymr),vec(nt1am(isymr)+1),1, 1789 & work(krmnab),1) 1790 end if 1791 call cclr_diascl(work(krmnab),two,isymr) 1792c calculate \sum_ab r^ab_kl * R^mn_ab 1793 call cc_r12mi2(work(kres),work(krabkl),work(krmnab),1,isymr,one, 1794 & work(kend1),lwrk1) 1795c apply projection operator 0.5*P_kl^ij: 1796 call cc_r12pklij(work(kres),isymr,'T',work(kend1),lwrk1) 1797 call dscal(ntr12sq(isymr),half,work(kres),1) 1798 end if 1799 1800 ! add the r12 doubles contribution: 1801 call cc_r12xi2b(work(kres),'T',work(kxintsq),1,'T', 1802 & work(kr12sq),isymr,'T',one) 1803 1804c ----------------------- 1805c print and pack results 1806c ----------------------- 1807 if (locdbg) then 1808 call around('R12 S*R before packing') 1809 call cc_prsqr12(work(kres),isymr,'T',1,.false.) 1810 end if 1811 if ( debug ) then 1812 write(lupri,*) 'cc_r12metric> norm^2(S * R12):', 1813 & ddot(ntr12am(isymr),work(kres),1,work(kres),1) 1814 end if 1815 iopt = 1 1816 call ccr12pck2(work(kpck),isymr,.false.,work(kres),'T',iopt) 1817 call cclr_diasclr12(work(kpck),basscl1,isymr) 1818 if (locdbg) then 1819 call around('R12 S*R after packing') 1820 call cc_prpr12(work(kpck),isymr,1,.false.) 1821 end if 1822 1823c ------------------ 1824c write S*R on file 1825c ------------------ 1826 call cc_wvec(lufs12,fs12am,ntr12am(isymr),ntr12am(isymr), 1827 * ifile,work(kpck)) 1828c 1829 if (ianr12.eq.2) then 1830c get \sum_kl r^kl_ab * R^ij_kl + R^ij_ab and put it on file 1831 call cc_r12getsr22(fc2am,lufc2,fc12am,lufc12,fs2am, 1832 & lufs2,isymr,ifile,lnoread,vec, 1833 & work(kend1),lwrk1) 1834 end if 1835c 1836 if (locdbg) write(lupri,*) 'Leaving CC_R12METRIC' 1837 call qexit('cc_r12metric') 1838 end 1839*=====================================================================* 1840 1841*=====================================================================* 1842 subroutine cc_r12getsr22(fc2am,lufc2,fc12am,lufc12, 1843 & fs2am,lufs2,isymt,ifile,lnoread, 1844 & vec,work,lwork) 1845*---------------------------------------------------------------------* 1846c purpose: get S_mu2,nu2' * R^ij_kl = \sum_kl r^kl_ab * R^ij_kl 1847c needed for ansatz 2 1848c 1849c H. Fliegl, C. Haettig spring 2005 1850c C. Neiss, 07.07.2005: modified, see cc_r12metric 1851*---------------------------------------------------------------------* 1852 implicit none 1853#include "priunit.h" 1854#include "ccsdsym.h" 1855c 1856#include "ccsdinp.h" 1857#include "ccorb.h" 1858c 1859#include "r12int.h" 1860#include "ccr12int.h" 1861 logical locdbg,lnoread 1862 parameter (locdbg = .false.) 1863 1864 character*(*) fc2am,fc12am,fs2am 1865 1866 integer lufc2,lufc12,lufs2,ifile,isymt,lwork,kend1,lwrk1 1867 integer krabkl,krmnab,krijkl,lunit,idummy 1868#if defined (SYS_CRAY) 1869 real work(*),two,one,half,factor,vec(*) 1870#else 1871 double precision work(*),two,one,half,factor,vec(*) 1872#endif 1873 parameter (one = 1.0d0, two = 2.0d0, half = 0.5d0) 1874 1875 call qenter('getsr22') 1876 1877 krabkl = 1 1878 krmnab = krabkl + nt2r12(1) 1879 krijkl = krmnab + nt2am(isymt) 1880 kend1 = krijkl + ntr12am(isymt) 1881 lwrk1 = lwork - kend1 1882 if (lwrk1.lt.0) then 1883 call quit('isufficient work space in getrs22') 1884 end if 1885 1886c read r^ab_kl integrals 1887 lunit = -1 1888 call gpopen(lunit,fr12r12,'unknown',' ','unformatted', 1889 & idummy,.false.) 1890 read(lunit)(work(krabkl-1+i),i=1,nt2r12(1)) 1891 call gpclose(lunit,'KEEP') 1892 1893 if (.not. lnoread) then 1894c read R^mn_ab trial vector (stored as triangular matrix) 1895 call cc_rvec(lufc2,fc2am,nt2am(isymt),nt2am(isymt), 1896 & ifile,work(krmnab)) 1897c read R^ij_kl 1898 call cc_rvec(lufc12,fc12am,ntr12am(isymt),ntr12am(isymt), 1899 * ifile,work(krijkl)) 1900 else 1901 call dcopy(nt2am(isymt),vec(nt1am(isymt)+1),1, 1902 & work(krmnab),1) 1903 call dcopy(ntr12am(isymt),vec(nt1am(isymt)+nt2am(isymt)+1), 1904 & 1,work(krijkl),1) 1905 end if 1906 call cclr_diascl(work(krmnab),two,isymt) 1907 call cclr_diasclr12(work(krijkl),ketscl,isymt) 1908 1909c get \sum_kl r^kl_ab * R^ij_kl and add it to work(krmnab) 1910 factor = one 1911 call cc_r12mi22(work(krmnab),work(krabkl),work(krijkl),1,isymt, 1912 & factor,work(kend1),lwrk1) 1913c call dscal(nt2am(isymt),half,work(krmnab),isymt) 1914 call cclr_diascl(work(krmnab),half,isymt) 1915c write result on file fs2am 1916 call cc_wvec(lufs2,fs2am,nt2am(isymt),nt2am(isymt),ifile, 1917 & work(krmnab)) 1918 if (locdbg) then 1919 call cc_rvec(lufs2,fs2am,nt2am(isymt),nt2am(isymt),ifile, 1920 & work(krmnab)) 1921 write(lupri,*)'FS2AM' 1922 write(lupri,*)'nt2am(isymt):', nt2am(isymt) 1923 write(lupri,*)'ntr12am(isymt):',ntr12am(isymt) 1924 call outpak(work(krmnab),nt2am(isymt),1,lupri) 1925 end if 1926 1927 call qexit('getsr22') 1928 end 1929*=====================================================================* 1930 1931*=====================================================================* 1932 subroutine cc_r12buildbmat(bmatsq,isymkl,isymmn, 1933 & isymij,idxij,bmatsq0, 1934 & xintsq,epskl,inmatkl,epsij,inmatij, 1935 & factor) 1936*---------------------------------------------------------------------* 1937c purpose: Build up R12 B-matrix B_(kl,mn)^(ij) for 1938c fixed index pair (ij) 1939c 1940c bmatsq B-matrix block of dimension nmatkl(isymkl)*nmatkl(isymmn) 1941c analog: bmatsq0, xintsq 1942c isymij symmetry of pair (ij) 1943c idxij sorted pair index (ij) 1944c 1945c Christian Neiss, April 2005, restructured September 2005 1946c--------------------------------------------------------------------- 1947 implicit none 1948#include "ccsdsym.h" 1949#include "priunit.h" 1950#include "ccsdinp.h" 1951#include "ccorb.h" 1952#include "r12int.h" 1953#include "ccr12int.h" 1954 1955 logical locdbg 1956 parameter (locdbg=.false.) 1957 1958 integer idxij,idxkl,idxmn,idxklmn 1959 integer isymij,isymkl,isymmn 1960 integer inmatkl(8),inmatij(8) 1961 1962#if defined (SYS_CRAY) 1963 real bmatsq(*),bmatsq0(*),xintsq(*),epskl(*),epsij(*),factor,two 1964#else 1965 double precision bmatsq(*),bmatsq0(*),xintsq(*),epskl(*), 1966 & epsij(*),factor,two 1967#endif 1968 parameter (two = 2.0D0) 1969 1970 call qenter('cc_r12buildbmat') 1971 if (locdbg) then 1972 write(lupri,*)'Entered CC_R12BUILDBMAT' 1973 end if 1974 1975c -------------- 1976c build B-matrix 1977c -------------- 1978 do idxmn = 1, nmatkl(isymmn) 1979 do idxkl = 1, nmatkl(isymkl) 1980 idxklmn = nmatkl(isymkl)*(idxmn-1) + idxkl 1981c write(lupri,*) 'idxklmn: ',idxklmn 1982 bmatsq(idxklmn) = bmatsq0(idxklmn) + 1983 & factor*(epskl(inmatkl(isymkl)+idxkl)+ 1984 & epskl(inmatkl(isymmn)+idxmn)- 1985 & two*epsij(inmatij(isymij)+idxij))* 1986 & xintsq(idxklmn) 1987 end do 1988 end do 1989 1990 if (locdbg) then 1991 call around('B-Matrix in CC_R12BUILDBMAT') 1992 write(lupri,*) 'idxij, epsilon(ij): ',idxij, 1993 & epsij(inmatkl(isymij)+idxij) 1994 call output(bmatsq,1,nmatkl(isymkl),1,nmatkl(isymmn), 1995 & nmatkl(isymkl),nmatkl(isymmn),1,lupri) 1996 end if 1997 1998 if (locdbg) then 1999 write(lupri,*)'Leaving CC_R12BUILDBMAT' 2000 end if 2001 2002 call qexit('cc_r12buildbmat') 2003 return 2004 end 2005*=====================================================================* 2006 2007*=====================================================================* 2008 subroutine cc_r12nxtam(omeg12,isymom,tamp12,lcceq, 2009 & er12,work,lwork) 2010c---------------------------------------------------------------------- 2011c purpose: get new start R12 amplitudes for each iteration 2012c by inversion of the B matrix. 2013c substitute for old cc_r12nxtam: working without 2014c singlet/triplet format, similiar to conventional 2015c part of vector function 2016c 2017c C. Neiss, december 2005 2018c---------------------------------------------------------------------- 2019 implicit none 2020#include "priunit.h" 2021#include "ccorb.h" 2022#include "ccsdsym.h" 2023#include "r12int.h" 2024#include "ccr12int.h" 2025#include "ccsdinp.h" 2026#include "iratdef.h" 2027#include "second.h" 2028 2029 logical locdbg,ldum,lcceq,etest 2030 parameter (locdbg = .false., etest = .false.) 2031 2032 integer kbmatsq,komegsq,kend1,lwrk1,kend2,lwrk2 2033 integer kxintsq,kscr,keij,kekl,kevl,kpvt,kzvec,koffom,koffb 2034 integer lwork,ian,iap,iopt,isym,lunit,idum 2035 integer isymij,isymmn,isymkl,idxij,isymom 2036 integer nij,nkl,inmatij(8),inmatkl(8),norbtsx 2037#if defined (SYS_CRAY) 2038 real tamp12(*),omeg12(*),er12,work(*),ddot,one,half,zero,factor, 2039 & rcond,thrzero,timnxtam 2040#else 2041 double precision tamp12(*),omeg12(*),er12,work(*),ddot,one,half, 2042 & zero,factor,rcond,thrzero,timnxtam 2043#endif 2044 character*3 aprox 2045 character*10 model 2046 parameter (one = 1.0D0, half = 0.5D0, zero = 0.0D0) 2047 parameter (thrzero = 1.0D-12) 2048 2049 call qenter('r12nxtam') 2050 if (locdbg) then 2051 write(lupri,*) 'Entered CC_R12NXTAM' 2052 call flshfo(lupri) 2053 end if 2054 2055 timnxtam = second() 2056 2057 kbmatsq = 1 2058 komegsq = kbmatsq+ nr12r12sq(1) 2059 kend1 = komegsq+ ntr12sq(isymom) 2060 2061 kscr = kend1 2062 kend1 = kscr + max(nr12r12sq(1),ntr12am(1)) 2063 2064 lwrk1 = lwork - kend1 2065 2066 if (lwrk1 .lt. 0) then 2067 call quit('Insufficient work space in cc_r12nxtam') 2068 end if 2069 2070c ------------- 2071c test symmetry 2072c ------------- 2073 if (lcceq .and. (isymom.ne.1)) then 2074 call quit('Symmetry error in CC_R12NXTAM') 2075 end if 2076c 2077c ---------------- 2078c read B matrices 2079c ---------------- 2080 lunit = -1 2081 call gpopen(lunit,fccr12b,'old',' ','unformatted', 2082 & idum,ldum) 2083 8881 read(lunit) ian, iap, aprox 2084 read(lunit) (work(kscr+i), i=0, nr12r12p(1)-1 ) 2085 if ((ian.ne.ianr12).or.(iap.ne.iapr12)) goto 8881 2086 call gpclose(lunit,'KEEP') 2087 iopt = 2 2088 call ccr12unpck2(work(kscr),1,work(kbmatsq),'N',iopt) 2089c call dscal(nr12r12sq(1),-one,work(kbmatsq),1) 2090c 2091 if (lcceq) then 2092c -------------------- 2093c read R12 amplitudes 2094c -------------------- 2095 iopt =32 2096 call cc_rdrsp('R0 ',0,1,iopt,model,idum,tamp12) 2097c 2098 if (iprint .ge. 5) then 2099 write(lupri,529) 'Norm^2 of r12am in NXTAM:', 2100 * ddot(ntr12am(1),tamp12,1,tamp12,1) 2101 end if 2102c 2103c ------------------------------- 2104c read vector function omega_R12 2105c ------------------------------- 2106 lunit = -1 2107 call gpopen(lunit,'CC_OMEGAR12','old',' ','unformatted', 2108 & idum,ldum) 2109 read(lunit) (omeg12(i), i=1, ntr12am(1)) 2110 call gpclose(lunit,'KEEP') 2111c 2112 write(lupri,529)'Norm^2 of Omegar12 in NXTAM:', 2113 & ddot(ntr12am(1),omeg12,1,omeg12,1) 2114c 2115 end if 2116c 2117 !repack Omega as Omega_(kl,ij) 2118 iopt = 1 2119 call ccr12unpck2(omeg12,isymom,work(komegsq),'N',iopt) 2120 2121 529 format(7x,a,d24.10) 2122c 2123 if (ianr12.eq.1) then 2124 if (iapr12.le.2 .or. iapr12.eq.4 .or. iapr12.eq.6) then 2125 factor = zero 2126 else 2127 factor = -half 2128 end if 2129 else 2130 if (iapr12.eq.2 .or. iapr12.eq.5) then 2131 factor = zero 2132 else 2133 factor = -half 2134 end if 2135 end if 2136c 2137c ----------------------------------------------- 2138c - read R12 overlap matrix, orbital energies 2139c - loop over pairs (ij) 2140c ----------------------------------------------- 2141c 2142 !calculate # of indices (kl) / (ij) and offset over all symmetries: 2143 nkl = 0 2144 nij = 0 2145 do isym = 1, nsym 2146 inmatkl(isym) = nkl 2147 inmatij(isym) = nij 2148 nkl = nkl + nmatkl(isym) 2149 nij = nij + nmatij(isym) 2150 end do 2151c 2152 kxintsq = kend1 2153 kekl = kxintsq + nr12r12sq(1) 2154 keij = kekl + nkl 2155 keij = kekl + nkl 2156 kend1 = keij + nij 2157 lwrk1 = lwork - kend1 2158 if (lwrk1 .lt. 0) then 2159 call quit('Insufficient work space in cc_r12nxtam') 2160 end if 2161c 2162 lunit = -1 2163 call gpopen(lunit,fccr12x,'old',' ','unformatted', 2164 & idum,ldum) 2165 9991 read(lunit) ian 2166 read(lunit) (work(kscr+i), i=0, nr12r12p(1)-1 ) 2167 if (ian.ne.ianr12) goto 9991 2168 call gpclose(lunit,'KEEP') 2169 iopt = 2 2170 call ccr12unpck2(work(kscr),1,work(kxintsq),'N',iopt) 2171c 2172 !read orbital energies 2173 lunit = -1 2174 call gpopen(lunit,'SIRIFC','OLD',' ','UNFORMATTED',idum, 2175 & .false.) 2176 rewind lunit 2177 call mollab('FULLBAS ',lunit,lupri) 2178 read (lunit) idum,norbtsx 2179 kevl = kend1 2180 kend2 = kevl + norbtsx 2181 lwrk2 = lwork - kend2 2182 if (lwrk2 .lt. 0) then 2183 call quit('Insufficient work space in cc_r12nxtam') 2184 end if 2185 read (lunit) (work(kevl+i-1), i=1,norbtsx) 2186 call gpclose(lunit,'KEEP') 2187c 2188 !calculate pair orbital energies: 2189 call cc_r12epair(work(kekl),nkl,work(kevl),inmatkl,imatkl,nrhfb) 2190 call cc_r12epair(work(keij),nij,work(kevl),inmatij,imatij,nrhf) 2191c 2192c ------------------------------- 2193c calculate update for amplitudes 2194c ------------------------------- 2195c 2196 do isymij = 1, nsym 2197 isymmn = muld2h(isymom,isymij) 2198 isymkl = isymmn 2199 if (nmatkl(isymkl).gt.0) then 2200 do idxij = 1, nmatij(isymij) 2201c call dzero(work(kscr),nmatkl(isymkl)*nmatkl(isymmn)) 2202 koffb = ir12r12sq(isymkl,isymmn) 2203 call cc_r12buildbmat(work(kscr),isymkl,isymmn, 2204 & isymij,idxij,work(kbmatsq+koffb), 2205 & work(kxintsq+koffb),work(kekl),inmatkl, 2206 & work(keij),inmatij,factor) 2207c 2208c invert B-matrix and contract with omega: 2209c 2210 kpvt = kend1 2211 kzvec = kpvt + nmatkl(isymkl)/irat + 1 2212 kend2 = kzvec + nmatkl(isymkl) 2213 lwrk2 = lwork - kend2 2214 if (lwrk2 .lt. 0) then 2215 call quit('Insufficient work space in cc_r12nxtam') 2216 end if 2217c 2218 call dgeco(work(kscr),nmatkl(isymkl),nmatkl(isymkl), 2219 & work(kpvt),rcond,work(kzvec)) 2220 if (rcond.lt.thrzero) then 2221 call quit('B-Matrix is singular in CC_R12NXTAM') 2222 end if 2223c 2224 koffom = komegsq + itr12sq(isymkl,isymij) + 2225 & nmatkl(isymkl)*(idxij-1) 2226 2227c write(lupri,*) 'Omega before dgesl:' 2228c call output(work(koffom),1,nmatkl(isymkl),1,1, 2229c & nmatkl(isymkl),1,1,lupri) 2230 2231 call dgesl(work(kscr),nmatkl(isymkl),nmatkl(isymkl), 2232 & work(kpvt),work(koffom),0) 2233 2234c write(lupri,*) 'Omega after dgesl:' 2235c call output(work(koffom),1,nmatkl(isymkl),1,1, 2236c & nmatkl(isymkl),1,1,lupri) 2237c 2238 end do 2239 end if 2240 end do 2241c 2242 iopt = 1 2243 call ccr12pck2(omeg12,isymom,.false.,work(komegsq),'N',iopt) 2244c 2245 if (lcceq) then 2246c ---------------------------------------------------- 2247c add old amplitudes to update to get new amplitudes 2248c ---------------------------------------------------- 2249 call daxpy(ntr12am(isymom),one,tamp12,1,omeg12,1) 2250c 2251 if (locdbg) then 2252 write(lupri,*) 'R12 amplitudes:' 2253 call cc_prpr12(tamp12,isymom,1,.false.) 2254 write(lupri,*) 'Updated R12 amplitudes:' 2255 call cc_prpr12(omeg12,isymom,1,.false.) 2256 end if 2257c 2258 if (etest) then 2259c ------------- 2260c read V matrix 2261c ------------- 2262 lunit = -1 2263 call gpopen(lunit,fccr12v,'old',' ','unformatted', 2264 & idum,ldum) 2265 7771 read(lunit) ian 2266 read(lunit) (work(kscr+i), i=0, ntr12am(1)-1) 2267 if (ian.ne.ianr12) goto 7771 2268 call gpclose(lunit,'KEEP') 2269c 2270 call cc_r12tcmepk(work(kscr),1,.false.) 2271 call cclr_diasclr12(work(kscr),half,isymom) 2272 er12 = 2.0D0*ddot(ntr12am(isymom),omeg12,1,work(kscr),1) 2273c 2274 write(lupri,*) 'R12 contribution to energy: ',er12 2275 end if 2276c 2277 end if ! lcceq 2278c 2279 if (locdbg) write(lupri,*) 'Leaving CC_R12NXTAM' 2280c 2281 timnxtam = second() - timnxtam 2282c write(lupri,111) 'CC_R12NXTAM', timnxtam 2283 111 FORMAT(7x,'Time used in',2x,A12,2x,': ',f10.2,' seconds') 2284c 2285 call qexit('r12nxtam') 2286 end 2287*=====================================================================* 2288 2289*=====================================================================* 2290 subroutine cc_r12epair(evlout,length,evlin,inmat,imat,nrhf1) 2291c---------------------------------------------------------------------- 2292c purpose: get pair orbital energies 2293c 2294c it is assumed that the R12 orbital energies are leading 2295c (incl. frozen orbitals, as for "FULLBAS" in SIRIFC) 2296c 2297c C. Neiss september 2005 2298c---------------------------------------------------------------------- 2299 implicit none 2300#include "priunit.h" 2301#include "ccorb.h" 2302#include "ccsdsym.h" 2303#include "r12int.h" 2304#include "ccsdinp.h" 2305 2306 logical locdbg 2307 parameter (locdbg = .false.) 2308 2309 integer kend1,lwrk1,length,kscr1 2310 integer lwork,idummy,icount,isym 2311 integer isymij,isymi,isymj,mi,mj,idxij 2312 integer nrhf1(8),inmat(8),imat(8,8),ioff(8) 2313#if defined (SYS_CRAY) 2314 real evlout(length),evlin(*),ei,ej 2315#else 2316 double precision evlout(length),evlin(*),ei,ej 2317#endif 2318c 2319 call qenter('cc_r12epair') 2320 if (locdbg) then 2321 write(lupri,*) 'Entered CC_R12EPAIR' 2322 call flshfo(lupri) 2323 end if 2324c 2325 icount = 0 2326 do isym = 1, nsym 2327 icount = icount + nrhffr(isym) 2328 ioff(isym) = icount 2329 icount = icount + nrhfb(isym) + norb1(isym) + norb2(isym) 2330 end do 2331c 2332 call dzero(evlout,length) 2333c 2334 do isymij = 1, nsym 2335 do isymi = 1, nsym 2336 isymj = muld2h(isymi,isymij) 2337 do i = 1, nrhf1(isymi) 2338c mi = iorb(isymi)+i 2339 mi = ioff(isymi)+i 2340 ei = evlin(mi) 2341 do j = 1, nrhf1(isymj) 2342 idxij = inmat(isymij) + imat(isymi,isymj) + 2343 & nrhf1(isymi)*(j-1)+i 2344c mj = iorb(isymj)+j 2345 mj = ioff(isymj)+j 2346 ej = evlin(mj) 2347 evlout(idxij) = ei + ej 2348 end do 2349 end do 2350 end do 2351 end do 2352c 2353 if (locdbg) then 2354 call around('CC_R12EPAIR: pair orbital energies') 2355 do isymij = 1, nsym 2356 do isymi = 1, nsym 2357 isymj = muld2h(isymij,isymi) 2358 write(lupri,*) 'Symmetry block ',isymi, isymj 2359 call output(evlout(1+inmat(isymij)+imat(isymi,isymj)), 2360 & 1,nrhf1(isymi),1,nrhf1(isymj), 2361 & nrhf1(isymi),nrhf1(isymj),1,lupri) 2362 end do 2363 end do 2364 end if 2365c 2366 if (locdbg) write(lupri,*) 'Leaving CC_R12EPAIR' 2367 2368 call qexit('cc_r12epair') 2369 end 2370*=====================================================================* 2371 2372