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 ccsdr12cd(ccsdr12, 20 & omegar12,isymom,t2am,isymt2,isymim, 21 & fniadj,luiadj,fnijda,luijda,it2del, 22 & xlamdah,isymlam, 23 & cmox,ilamdx, 24 & work,lwork) 25*----------------------------------------------------------------------* 26* Purpose: compute the C and D contributions to the R12 result vector 27* 28* C. Haettig, C. Neiss, spring 2006 29*----------------------------------------------------------------------* 30 implicit none 31#include "priunit.h" 32#include "dummy.h" 33#include "ccorb.h" 34#include "ccsdsym.h" 35#include "r12int.h" 36#include "ccr12int.h" 37 38 logical locdbg 39 parameter (locdbg = .FALSE.) 40 integer isym0 41 parameter ( isym0 = 1) 42 real*8 zero,one 43 parameter ( zero=0.0d0, one=1.0d0 ) 44 45* input: 46 logical ccsdr12 47 character fnijda*8, fniadj*8 48 integer isymlam, ilamdx(8,8), luijda, luiadj, isymom, isymt2, 49 & it2del(*), lwork, isymim 50 real*8 xlamdah(*),omegar12(*),t2am(*),cmox(*),work(*) 51 52* local: 53 logical lcbs, ldum 54 character cdummy*8 55 integer isym, nbas2(8), iviraop(8,8), 56 & imatpb(8,8), nviraop(8), nmatpb(8), 57 & kend1, lwrk1, kemat1, kfckvao, kfckvaop, kfgdp, kend2, 58 & lwrk2, lunit, isymp, isydel, isygam, isymc, kscr1, kscr2, 59 & kend3, lwrk3, koffv, len1, len2, koffl, koffgdp, nbasg, 60 & nvirc, koffe1, koffcx, nbasd, norbp, kcdint, iopt, 61 & komegpk, isym1, isym2, igdp(8,8), ngdp(8), iopttcme, 62 & ioptr12, idum, mtotbas 63 real*8 fac, ddot 64 65*----------------------------------------------------------------------* 66* precompute non-standard symmetry offsets and dimensions: 67*----------------------------------------------------------------------* 68 mtotbas = 0 69 do isym = 1, nsym 70 ! total number of basis functions in primary + aux. basis 71 nbas2(isym) = mbas1(isym) + mbas2(isym) 72 mtotbas = mtotbas + nbas2(isym) 73 end do 74 75 do isym = 1, nsym 76 nviraop(isym) = 0 77 ngdp(isym) = 0 78 nmatpb(isym) = 0 79 do isym2 = 1, nsym 80 isym1 = muld2h(isym2,isym) 81 iviraop(isym1,isym2) = nviraop(isym) 82 nviraop(isym) = nviraop(isym) + nvir(isym1)*mbas2(isym2) 83 imatpb(isym1,isym2) = nmatpb(isym) 84 nmatpb(isym) = nmatpb(isym) + norb2(isym1)*nvir(isym2) 85 igdp(isym1,isym2) = ngdp(isym) 86 ngdp(isym) = ngdp(isym) + mbas1(isym1)*nbas2(isym2) 87 end do 88 end do 89 90 kend1 = 1 91 lwrk1 = lwork 92 93 if (isymom.ne.muld2h(isymim,isymt2)) then 94 call quit('Symmetry mismatch in ccsdr12cd!') 95 end if 96 97 ! shift it2del by one (offset vs. start address) 98 do i = 1, mtotbas 99 it2del(i) = it2del(i) + 1 100c write(lupri,*) 'i,it2del(i):',i,it2del(i) 101 end do 102 103*----------------------------------------------------------------------* 104* transform the leading index of the E1 intermediate (for CCSD(R12) 105* identical with the correlation contribution to the Fock-hat matrix) 106* from the AO to the orthogonal complementary basis: 107*----------------------------------------------------------------------* 108 ! allocate work space for E(p',c) intermediate 109 kemat1 = kend1 110 kend1 = kemat1 + nmatpb(isymim) 111 lwrk1 = lwork - kend1 112 if (lwrk1 .lt. 0) call quit('Insufficient memory in CCSDR12CD 1') 113 114 ! allocate work space for the half-transformed and AO 115 ! intermediates needed in the following section: 116 kfckvao = kend1 117 kfckvaop = kfckvao + nemat1(isymim) 118 kfgdp = kfckvaop + nviraop(isymim) 119 kend2 = kfgdp + ngdp(isym0) 120 lwrk2 = lwork - kend2 121 if (lwrk2 .lt. 0) call quit('Insufficient memory in CCSDR12CD 2') 122 123 ! read the Fhat(c,delta) and Fhat(c,delta-p) matrices 124 lunit = -1 125 call gpopen(lunit,'CCFHATADEL','UNKNOWN',' ','UNFORMATTED', 126 & idummy,.false.) 127 read(lunit) (work(kfckvao-1+i), i=1,nemat1(isymim)) 128 read(lunit) (work(kfckvaop-1+i),i=1,nviraop(isymim)) 129 call gpclose(lunit,'KEEP') 130 131C if (locdbg) then 132C write(lupri,*) 'Fhat(c,delta) and Fhat(c,delta-p):' 133C do isym2 = 1, nsym 134C isym1 = muld2h(isym2,isymim) 135C write(lupri,*) 'Symmetry block ',isym1,isym2 136C call output(work(kfckvao+iemat1(isym1,isym2)), 137C & 1,nvir(isym1),1,mbas1(isym2), 138C & nvir(isym1),mbas1(isym2),1,lupri) 139C call output(work(kfckvaop+iviraop(isym1,isym2)), 140C & 1,nvir(isym1),1,mbas2(isym2), 141C & nvir(isym1),mbas2(isym2),1,lupri) 142C end do 143C end if 144 145 ! read the F^val(gamma,delta-p) matrix, which contains the 146 ! two-electron contribution of a Fock matrix computed from 147 ! the SCF valence electron density (i.e. w/o core orbitals) 148 lunit = -1 149 call gpopen(lunit,'R12FOCK','UNKNOWN',' ','UNFORMATTED',idummy, 150 & .false.) 151 call readt(lunit,ngdp(1),work(kfgdp)) 152 call gpclose(lunit,'KEEP') 153 154 155c ------------------------------------------------------------- 156c transform leading index to virtual and substract it from 157c the corresponding matrix computed from the CC Lambda density, 158c and then transform the remaining AO index to the orthogonal 159c complementary basis to get the matrix E1(p',c): 160c ------------------------------------------------------------- 161 if (isymlam.ne.isymim) call quit('symmetry mismatch in ccsdr12cd') 162 do isymp = 1, nsym 163 isydel = isymp 164 isygam = isydel 165 isymc = muld2h(isymlam,isygam) 166 167 len1 = nvir(isymc) * mbas1(isydel) 168 len2 = nvir(isymc) * mbas2(isydel) 169 170 kscr1 = kend2 171 kscr2 = kscr1 + len1 172 kend3 = kscr2 + len2 173 lwrk3 = lwork - kend3 174 if (lwrk3 .lt. 0) call quit('Insufficient core in CCSDR12CD 3') 175 176 ! ......................................................... 177 ! initialize scratch array with the half-transformed E1 178 ! or fock matrix computed from the Lambda density, with 179 ! the primary AOs followed immediately by the aux. AOs 180 ! ......................................................... 181 koffv = kfckvao + iemat1(isymc,isydel) 182 call dcopy(len1,work(koffv),1,work(kscr1),1) 183 koffv = kfckvaop + iviraop(isymc,isydel) 184 call dcopy(len2,work(koffv),1,work(kscr2),1) 185 186 ! ......................................................... 187 ! substract the contribution from the SCF density matrix 188 ! with the leading index transformed to the virtual basis 189 ! ......................................................... 190 koffl = iglmvi(isygam,isymc) + 1 191 koffgdp = kfgdp + igdp(isygam,isydel) 192 193 nbasg = max(mbas1(isygam),1) 194 nvirc = max(nvir(isymc),1) 195 196 call dgemm('T','N',nvir(isymc),nbas2(isydel),mbas1(isygam), 197 & -one,xlamdah(koffl),nbasg,work(koffgdp),nbasg, 198 & one,work(kscr1),nvirc) 199 200C if (locdbg) then 201C write(lupri,*) 'in CCSDR12CD: Norm^2(FGDP) = ', 202C & ddot(mbas1(isygam)*nbas2(isydel),work(koffgdp),1, 203C & work(koffgdp),1) 204C call output(work(koffgdp),1,mbas1(isygam),1,nbas2(isydel), 205C & mbas1(isygam),nbas2(isydel),1,lupri) 206C write(lupri,*) 'in CCSDR12CD: Norm^2(SCR) = ', 207C & ddot(nvir(isymc)*nbas2(isydel),work(kscr1),1,work(kscr1),1) 208C call output(work(kscr1),1,nvir(isymc),1,nbas2(isydel), 209C & nvir(isymc),nbas2(isydel),1,lupri) 210C end if 211 212 ! ......................................................... 213 ! finally transform the outer index to the orthogonal 214 ! complementary basis and store at work(kemat1) 215 ! ......................................................... 216 koffe1 = kemat1 + imatpb(isymp,isymc) 217 koffcx = 1 + ilamdx(isydel,isymp) + nbas2(isydel)*norb1(isymp) 218 219 nbasd = max(nbas2(isydel),1) 220 norbp = max(norb2(isymp),1) 221 222 call dgemm('T','T',norb2(isymp),nvir(isymc),nbas2(isydel), 223 & one,cmox(koffcx),nbasd,work(kscr1),nvirc, 224 & zero,work(koffe1),norbp) 225 226 end do 227 228 if (locdbg) then 229 write(lupri,*) 'in CCSDR12CD: Final Norm^2(EMAT1) = ', 230 & ddot(nmatpb(isymim),work(kemat1),1,work(kemat1),1) 231 end if 232 233*----------------------------------------------------------------------* 234* allocate work space for C & D intermediates: 235*----------------------------------------------------------------------* 236 kcdint = kend1 237 kend1 = kcdint + ntg2sq(isymim) 238 lwrk1 = lwork - kend1 239 if (lwrk1 .lt. 0) call quit('Insufficient memory in CCSDR12CD 4') 240 241*----------------------------------------------------------------------* 242* initialize result vector: 243*----------------------------------------------------------------------* 244 call dzero(omegar12,ntg2sq(isymom)) 245 246*----------------------------------------------------------------------* 247* transform C intermediate to orthonormal complementary basis 248*----------------------------------------------------------------------* 249 iopt = 1 250 lcbs = .true. 251 call dzero(work(kcdint),ntg2sq(isymim)) 252 call cc_iajb2(work(kcdint), isymim, iopt, .false., .false., lcbs, 253 & luijda, fnijda, it2del, cmox, isym0, 254 & idummy, cdummy, idummy, dummy, idummy, 255 & work(kend1), lwrk1 ) 256 257 if (locdbg) then 258 write(lupri,*) 'Norm^2 of C-Int. after CC_IAJB2: ', 259 & ddot(ntg2sq(isymim),work(kcdint),1,work(kcdint),1) 260 end if 261 262*----------------------------------------------------------------------* 263* add E1(p',c) intermediate to the k=1 diagonal of the C intermediate: 264*----------------------------------------------------------------------* 265 fac = -0.5d0 266 call cc_cdbar2(work(kcdint),work(kemat1),fac,.true.,isymim) 267 if (locdbg) then 268 write(lupri,*) 'Norm^2 of C-Int. after CC_CDBAR2: ', 269 & ddot(ntg2sq(isymim),work(kcdint),1,work(kcdint),1) 270 end if 271 272*----------------------------------------------------------------------* 273* calculate the C term and add to result vector: 274*----------------------------------------------------------------------* 275 ioptr12 = 1 276 call cc_cd('C', +1, ioptr12, 277 & omegar12, isymom, t2am, isymt2, 278 & work(kcdint), isymim, work(kend1), lwrk1 ) 279 280 if (locdbg) then 281 write(lupri,*) 'Norm^2 of OmegaR12-Int. after CC_CD(C): ', 282 & ddot(ntg2sq(isymom),omegar12,1,omegar12,1) 283 end if 284 285*----------------------------------------------------------------------* 286* transform D intermediate to orthonormal complementary basis 287*----------------------------------------------------------------------* 288 iopt = 1 289 lcbs = .true. 290 call dzero(work(kcdint),ntg2sq(isymim)) 291 call cc_iajb2(work(kcdint), isymim, iopt, .false., .false., lcbs, 292 & luiadj, fniadj, it2del, cmox, isym0, 293 & idummy, cdummy, idummy, dummy, idummy, 294 & work(kend1), lwrk1 ) 295 296*----------------------------------------------------------------------* 297* add E1(p',c) intermediate to the k=1 diagonal of the D intermediate: 298*----------------------------------------------------------------------* 299 fac = 0.5d0 300 call cc_cdbar2(work(kcdint),work(kemat1),fac,.true.,isymim) 301 302*----------------------------------------------------------------------* 303* calculate the D term and add to result vector: 304*----------------------------------------------------------------------* 305 ! form in place 2C-E combination of T2 306 iopttcme = 1 307 call ccsd_tcmepk(t2am,one,isymt2,iopttcme) 308 309 ioptr12 = 1 310 call cc_cd('D', +1, ioptr12, 311 & omegar12, isymom, t2am, isymt2, 312 & work(kcdint), isymim, work(kend1), lwrk1 ) 313 314 if (locdbg) then 315 write(lupri,*) 'Norm^2 of OmegaR12-Int. after CC_CD(D): ', 316 & ddot(ntg2sq(isymom),omegar12,1,omegar12,1) 317 end if 318 319*----------------------------------------------------------------------* 320* contract Omega(ai,p'j) with r12 integrals to Omega(xi,yj): 321*----------------------------------------------------------------------* 322 komegpk = 1 323 kend1 = komegpk + ntr12am(isymom) 324 lwrk1 = lwork - kend1 325 if (lwrk1 .lt. 0) call quit('Insufficient memory in CCSDR12CD 5') 326 327 ! read present result vector from file: 328 lunit = -1 329 call gpopen(lunit,'CC_OMEGAR12','unknown',' ','unformatted', 330 & idum,ldum) 331 read(lunit) (work(komegpk+i), i=0, ntr12am(isymom)-1 ) 332 call gpclose(lunit,'KEEP') 333 334 if (locdbg) then 335 write(lupri,*) 'in CCSDR12CD: Norm^2 of OMEGAR12 before '// 336 & 'CD contributions: ', 337 & ddot(ntr12am(isymom),work(komegpk),1,work(komegpk),1) 338 end if 339 340 call ccsdr12oxr(work(komegpk),omegar12,isymom,work(kend1),lwrk1) 341 342 if (locdbg) then 343 write(lupri,*) 'in CCSDR12CD: Norm^2 of OMEGAR12 after '// 344 & 'CD contributions: ', 345 & ddot(ntr12am(isymom),work(komegpk),1,work(komegpk),1) 346 end if 347 348 ! write updated result vector back to file: 349 lunit = -1 350 call gpopen(lunit,'CC_OMEGAR12','unknown',' ','unformatted', 351 & idum,ldum) 352 write(lunit) (work(komegpk+i), i=0, ntr12am(isymom)-1 ) 353 call gpclose(lunit,'KEEP') 354 355*----------------------------------------------------------------------* 356* restore it2del 357*----------------------------------------------------------------------* 358 ! shift it2del by one (offset vs. start address) 359 do i = 1, mtotbas 360 it2del(i) = it2del(i) - 1 361 end do 362 363 return 364 end 365*----------------------------------------------------------------------* 366* END OF SUBROUTINE CCSDR12CD * 367*======================================================================* 368