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