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 ccsdr12ao(ccsdr12,
20     &                     t2am,xlambdap,xlambdah,
21     &                     fniadj,luiadj,fnijda,luijda,
22     &                     cpfil,lucp,dpfil,ludp,e1pim,
23     &                     timintr12,timrdao,timtrbt,
24     &                     timc,timd,timt2tr,timt2bt,
25     &                     work,lwork)
26*----------------------------------------------------------------------*
27*  Purpose: compute contributions to CCSDR12 requiring the calculation
28*           of integrals with auxiliary basis function indices
29*
30* C. Haettig, C. Neiss, spring 2006
31*----------------------------------------------------------------------*
32      implicit none
33#include "priunit.h"
34#include "dummy.h"
35#include "maxorb.h"
36#include "mxcent.h"
37#include "aovec.h"
38#include "iratdef.h"
39#include "ccorb.h"
40#include "ccisao.h"
41#include "blocks.h"
42#include "ccsdsym.h"
43#include "ccsdinp.h"
44#include "cbieri.h"
45#include "distcl.h"
46#include "eritap.h"
47#include "r12int.h"
48#include "ccsdio.h"
49#include "second.h"
50      logical locdbg
51      parameter (locdbg = .false.)
52
53      integer isym0
54      parameter (isym0 = 1)
55
56      real*8  zero,half,one,two,xmhalf,xmone
57      parameter (zero = 0.0d0, half = 0.5d0, one = 1.0d0, two = 2.0d0)
58      parameter (xmhalf = -0.5d0, xmone= -1.0d0 )
59
60* input:
61      logical ccsdr12
62      integer lwork, lucp, ludp, luiadj, luijda
63      character*(*) cpfil, dpfil, fniadj, fnijda
64      real*8  work(*), xlambdap(*), xlambdah(*), t2am(*), e1pim(*),
65     &     timintr12, timrdao, timtrbt, timc, timd, timt2tr, timt2bt
66
67* local
68      logical lauxg, temp_direct
69      integer indexa(MXCORB_CC)
70      integer ioff, ibasx(8), kend1, lwrk1, kendsv, lwrksv,
71     &        icdel1, isymd1, ntot, illl, numdis, idel2, idel, isymd,
72     &        isydis, irecord, leniaj, isygam, isyalbe, igam, iadr,
73     &        icon, iv, isym, nviraop(8), iviraop(8,8), isym1, isym2,
74     &        kxint, kxiadj, kxijda, kend3, lwrk3, koffg, kdsrhf,
75     &        kt2amt, kend2, lwrk2, kfckvaop, isyfck, lunit, ioptr12
76      integer kccfb1,kindxb,kfree,lfree,ntosym,icount,ibasd, iopte,
77     &        kodcl1,kodcl2,kodbc1,kodbc2,krdbc1,krdbc2,
78     &        kodpp1,kodpp2,krdpp1,krdpp2,krecnr
79      real*8  dtime, factc, factd, ddot
80
81      call qenter('ccsdr12ao')
82      if (.not.dumpcd) call quit('CCSDR12AO requires DUMPCD=.TRUE.')
83
84      ! switch to integrals with delta index from auxiliary basis:
85      mbsmax = 5
86      loopdp = .true.
87
88      ! integrals with auxiliary basis function are not on file
89      ! -> switch locally to direct mode to calculate them
90      temp_direct = direct
91      direct = .true.
92
93      ! substract aux. functions when calculating index of g index
94      ! within irrep from index running over both basis sets and
95      ! all symmetries
96      ! (needed since isao is overwritten)
97      lauxg = .true.
98
99      ioff    = 0
100      ibasx(1) = 0
101      do isym = 1,nsym
102        if (isym.gt.1) ibasx(isym) = ibasx(isym-1)+mbas2(isym-1)
103        if (ioff+mbas1(isym)+mbas2(isym).gt.MXCORB_CC)
104     &     call quit('CCSDR12AO')
105        do i = 1,mbas1(isym)+mbas2(isym)
106          ioff = ioff + 1
107          isao(ioff) = isym
108        end do
109      end do
110
111      kend1 = 1
112      lwrk1 = lwork
113
114*----------------------------------------------------------------------*
115*     initialization of symmetry arrays and allocation of work space
116*     for Fhat(del',a)
117*----------------------------------------------------------------------*
118      do isym = 1, nsym
119        icount = 0
120        do isym2 = 1, nsym
121          isym1 = muld2h(isym2,isym)
122          iviraop(isym1,isym2) = icount
123          icount = icount + nvir(isym1)*mbas2(isym2)
124        end do
125        nviraop(isym) = icount
126      end do
127
128      isyfck = 1 ! symmetry of the Fhat matrix that will be computed
129
130      kfckvaop = kend1
131      kend1   = kfckvaop + nviraop(isyfck)
132      lwrk1   = lwork   - kend1
133      if (lwrk1 .lt. 0) call quit('Insufficient core in CCSDR12AO')
134
135      call dzero(work(kfckvaop),nviraop(isyfck))
136
137*----------------------------------------------------------------------*
138*     prepare cluster amplutides with transposed occupied indices:
139*----------------------------------------------------------------------*
140      if ((.not. direct) .and. t2tcor) then
141         kt2amt = kend1
142         kend1  = kt2amt + nt2sq(1)
143         lwrk1  = lwork  - kend1
144         if (lwrk1 .lt. 0) call quit('Insufficient core in CCSDR12AO')
145         call dcopy(nt2sq(1),t2am,1,work(kt2amt),1)
146         call ccsd_t2tp(work(kt2amt),work(kend1),lwrk1,isym0)
147      end if
148
149*----------------------------------------------------------------------*
150*     intialize integral program
151*----------------------------------------------------------------------*
152      if (direct) then
153         dtime  = second()
154         if (herdir) then
155            call herdi1(work(kend1),lwrk1,ipreri)
156         else
157            kccfb1 = kend1
158            kindxb = kccfb1 + mxprim*mxcont
159            kend1  = kindxb + (8*mxshel*mxcont + 1)/irat
160            lwrk1  = lwork  - kend1
161            if (lwrk1 .lt.0)
162     &        call quit('Insufficient work space in CC_MOFCONR12')
163            call eridi1(kodcl1,kodcl2,kodbc1,kodbc2,krdbc1,krdbc2,
164     &                kodpp1,kodpp2,krdpp1,krdpp2,
165     &                kfree,lfree,kend1,work(kccfb1),work(kindxb),
166     &                work(kend1),lwrk1,ipreri)
167            kend1 = kfree
168            lwrk1 = lfree
169         endif
170         timintr12 = timintr12 + ( second() - dtime )
171         ntosym = 1
172      else
173         ntosym = nsym
174      endif
175
176*----------------------------------------------------------------------*
177*     start the loop over integral distributions:
178*----------------------------------------------------------------------*
179
180      kendsv = kend1
181      lwrksv = lwrk1
182
183      icdel1 = 0
184      do isym = 1, nsym
185        icdel1 = icdel1  + nt2bcd(isym)*mbas1(isym)
186      end do
187
188      do isymd1 = 1,ntosym
189
190        if (direct) then
191           if (herdir) then
192              ntot = maxshl
193           else
194              ntot = mxcall
195           endif
196        else
197           ntot = nbas(isymd1)
198        endif
199
200        do illl = 1,ntot
201
202          if (direct) then
203             dtime = second()
204             kend1 = kendsv
205             lwrk1 = lwrksv
206             if (herdir) then
207                call herdi2(work(kend1),lwrk1,indexa,illl,numdis,
208     &                      ipreri)
209             else
210                call eridi2(illl,indexa,numdis,0,0,
211     &                    work(kodcl1),work(kodcl2),work(kodbc1),
212     &                    work(kodbc2),work(krdbc1),work(krdbc2),
213     &                    work(kodpp1),work(kodpp2),work(krdpp1),
214     &                    work(krdpp2),work(kccfb1),work(kindxb),
215     &                    work(kend1), lwrk1,ipreri)
216             endif
217
218             krecnr = kend1
219             kend1  = krecnr + (nbufx(0) - 1)/irat + 1
220             lwrk1  = lwork  - kend1
221             if (lwrk1 .lt.0)
222     &         call quit('Insufficient work space in CCSDR12AO')
223             timintr12 = timintr12 + ( second() - dtime )
224
225             if (t2tcor) then
226                kt2amt = kend1
227                kend1  = kt2amt + nt2sq(1)
228                lwrk1  = lwork  - kend1
229                if (lwrk1.lt.0) call quit('Insuff. core in CCSDR12AO')
230                call dcopy(nt2sq(1),t2am,1,work(kt2amt),1)
231                call ccsd_t2tp(work(kt2amt),work(kend1),lwrk1,isym0)
232             end if
233          else
234             numdis = 1
235          endif
236
237c         ------------------------------------------------------
238c         loop over AOs delta included in the distribution
239c         calculated in the above call to the integrals program:
240c         ------------------------------------------------------
241          do idel2 = 1,numdis
242
243            if (direct) then
244               idel  = indexa(idel2)
245               isymd = isao(idel)
246            else
247               idel  = ibas(isymd1) + ibasx(isymd1) + illl
248               isymd = isymd1
249            endif
250
251            isydis = isymd
252
253C           ------------------------------------------------------------
254C           record number and start addresses for intermediates on file:
255C           append the records for auxiliary basis functions after those
256C           where delta is primary basis function
257C           ------------------------------------------------------------
258            irecord = mbas1t + idel - mbas1(isymd) - ibas(isymd)
259            if (locdbg)
260     &        write(lupri,*) 'in CCSDR12AO: irecord = ',irecord
261            it2del(irecord) = icdel1
262            icdel1 = icdel1 + nt2bcd(isydis)
263
264C           --------------------------------------------------
265C           allocate memory for 3-index batch of AO integrals
266C           and read batch into memory:
267C           --------------------------------------------------
268            kxint  = kend1
269            kend2  = kxint + ndisao(isydis)
270            lwrk2  = lwork  - kend2
271            if (lwrk2 .lt. 0)
272     &        call quit('Insufficient work space in CCSDR12AO')
273
274            dtime   = second()
275            call ccrdao(work(kxint),idel,idel2,work(kend2),lwrk2,
276     &                  work(krecnr),direct)
277            dtime   = second() - dtime
278            timrdao = timrdao  + dtime
279
280*----------------------------------------------------------------------*
281* compute integrals (ia|delta' j) and (ij|delta' a), where i is obtained
282* by transformation with lambda-particle and j and a by transformation
283* with lambda-hole. the integtrals are saved on disk.
284*----------------------------------------------------------------------*
285            if (ccsdr12 .and. (ianr12.eq.2.or.ianr12.eq.3)) then
286               leniaj = nt2bcd(isydis)
287
288               kxiadj = kend2
289               kxijda = kxiadj + leniaj
290               kend3  = kxijda + leniaj
291               lwrk3  = lwork  - kend3
292               if (lwrk3 .lt. 0)
293     &           call quit('Insufficient space in CCSDR12AO')
294
295               call dzero(work(kxiadj),leniaj)
296               call dzero(work(kxijda),leniaj)
297CCN
298C              write(lupri,*) 'in CCSDR12AO:'
299C              write(lupri,*) 'Norm^2 of kxint = ',
300C    &           ddot(ndisao(isydis),work(kxint),1,work(kxint),1)
301CCN
302
303               do isygam = 1, nsym
304                  isyalbe = muld2h(isydis,isygam)
305               do g = 1, mbas1(isygam)
306                  igam = g + ibas(isygam) + ibasx(isygam)
307
308                  koffg = kxint + idsaog(isygam,isydis)
309     &                    + nnbst(isyalbe)*(g-1)
310
311                  call cc_iajb( work(koffg), isyalbe, dummy, isym0,
312     &                          idel, igam, lauxg, ibasx,
313     &                          dummy, work(kxiadj), work(kxijda),
314     &                          dummy, dummy, dummy,
315     &                          xlambdap, xlambdah, isym0,
316     &                          dummy, dummy, isym0,
317     &                          xlambdap, xlambdah, isym0,
318     &                          dummy, dummy, isym0,
319     &                          work(kend3), lwrk3,   3,
320     &                          .false., .false.,  .true.,
321     &                          .false., .false.,  0      )
322               end do
323               end do
324
325c              ------------------------------------
326c              update Fhat_{del a}:
327c              ------------------------------------
328               ibasd = idel - mbas1(isymd) - ibas(isymd) - ibasx(isymd)
329               if (locdbg) then
330                 write(lupri,*) 'in CCSDR12AO: isymd, ibasd = ',
331     &                                         isymd, ibasd
332               end if
333
334               call cc_fckdela(ibasd,isymd,work(kfckvaop),isyfck,
335     &                         work(kxijda),work(kxiadj),iviraop)
336
337c              ------------------------------------
338c              transform (ia|del j) to L(ia|del j):
339c              ------------------------------------
340               call dscal(leniaj, two,work(kxiadj),1)
341               call daxpy(leniaj,-one,work(kxijda),1,work(kxiadj),1)
342
343c              --------------------------------------------
344c              write 3-index transformed integrals to disk:
345c              --------------------------------------------
346               iadr = it2del(irecord) + 1
347               call putwa2(luiadj,fniadj,work(kxiadj),iadr,leniaj)
348               call putwa2(luijda,fnijda,work(kxijda),iadr,leniaj)
349
350               if (locdbg) then
351                 write(lupri,*) 'in CCSDR12AO: iadr = ', iadr
352                 write(lupri,*) 'in CCSDR12AO: Norm^2 of XIADJ: ',
353     &              ddot(leniaj,work(kxiadj),1,work(kxiadj),1)
354                 write(lupri,*) 'in CCSDR12AO: Norm^2 of XIJDA: ',
355     &              ddot(leniaj,work(kxijda),1,work(kxijda),1)
356               end if
357
358            end if
359
360*----------------------------------------------------------------------*
361* Compute the C' and D' intermediates with delta' an auxiliary
362* basis function using CCRHS_C and CCRHS_D. These routines require
363* the amplitudes as full square matrix on t2am
364*----------------------------------------------------------------------*
365            if (ccsdr12 .and. (ianr12.eq.2.or.ianr12.eq.3)) then
366
367C              -------------------------------------------------------
368C              transform gamma index in the integral batch to occupied
369C              using the lambda-partical = CMO coefficients:
370C              -------------------------------------------------------
371               kdsrhf = kend2
372               kend3  = kdsrhf + ndsrhf(isymd)
373               lwrk3  = lwork  - kend3
374               if (lwrk3 .lt. 0)
375     &           call quit('Insufficient space in CCSDR12AO')
376
377               dtime   = second()
378               call cctrbt(work(kxint),work(kdsrhf),xlambdap,
379     *                     isym0,work(kend3),lwrk3,isydis)
380               dtime   = second() - dtime
381               timtrbt = timtrbt + dtime
382
383
384C              -----------------------------
385C              calculate the C intermediate:
386C              -----------------------------
387               dtime = second()
388               factc = xmone
389               icon  = 2
390               ioptr12 = 0 !only calculate ONE C-intermediate
391               iopte = 1
392               iv    = 1
393               if (.not. t2tcor) then
394                  call ccrhs_c(work(kxint),work(kdsrhf),dummy,
395     *                         t2am,isym0,xlambdap,dummy,
396     *                         xlambdah,xlambdap,isym0,
397     *                         xlambdap,isym0,
398     *                         dummy,e1pim,work(kend3),lwrk3,
399     *                         irecord,isymd,factc,icon,ioptr12,iopte,
400     *                         lucp,cpfil,idummy,dummy,iv)
401               else
402                  call ccrhs_c(work(kxint),work(kdsrhf),dummy,
403     *                         work(kt2amt),isym0,
404     *                         xlambdap,dummy,
405     *                         xlambdah,xlambdap,isym0,
406     *                         xlambdap,isym0,
407     *                         dummy,e1pim,work(kend3),lwrk3,
408     *                         irecord,isymd,factc,icon,ioptr12,iopte,
409     *                         lucp,cpfil,idummy,dummy,iv)
410               end if
411               dtime   = second() - dtime
412               timc    = timc     + dtime
413
414c              -------------------------
415c              transform T2 to 2T2 - T2.
416c              -------------------------
417               dtime   = second()
418               if (t2tcor) then
419                  call dscal(nt2sq(1),two,t2am,1)
420                  call daxpy(nt2sq(1),-one,work(kt2amt),1,t2am,1)
421               else
422                  call ccrhs_t2tr(t2am,work(kend3),lwrk3,isym0)
423               end if
424               dtime   = second() - dtime
425               timt2tr = timt2tr  + dtime
426
427c              -----------------------------
428c              calculate the D intermediate:
429c              -----------------------------
430               dtime = second()
431               factd = one
432               icon  = 2
433               ioptr12 = 0
434               iopte = 1
435               iv    = 1
436               call ccrhs_d(work(kxint),work(kdsrhf),dummy,
437     *                      t2am,isym0,
438     *                      xlambdap,dummy,
439     *                      xlambdah,xlambdap,isym0,
440     *                      xlambdah,isym0,
441     *                      dummy,e1pim,work(kend3),lwrk3,
442     *                      irecord,isymd,factd,icon,ioptr12,iopte,
443     *                      ludp,dpfil,idummy,dummy,iv)
444               dtime   = second() - dtime
445               timd    = timd     + dtime
446
447
448c              -------------------------
449c              restore T2 from 2T2 - T2:
450c              -------------------------
451               dtime   = second()
452               if (t2tcor) then
453                  call daxpy(nt2sq(1),one,work(kt2amt),1,t2am,1)
454                  call dscal(nt2sq(1),half,t2am,1)
455               else
456                  call ccrhs_t2bt(t2am,work(kend3),lwrk3,isym0)
457               end if
458               dtime   = second() - dtime
459               timt2bt = timt2bt  + dtime
460
461            end if
462
463          end do ! idel2
464        end do ! illl
465      end do ! isymd1
466
467*----------------------------------------------------------------------*
468*     end of the loop over integral distributions...
469*     restore default settings for the integral evaluation:
470*----------------------------------------------------------------------*
471
472      mbsmax = 4
473      loopdp = .false.
474      direct = temp_direct
475
476      ioff   = 0
477      do isym = 1,nsym
478        do i = 1,nbas(isym)
479          ioff = ioff + 1
480          isao(ioff) = isym
481        end do
482      end do
483
484*----------------------------------------------------------------------*
485*     save Fhat(a,del') on file:
486*----------------------------------------------------------------------*
487      lunit = -1
488      call gpopen(lunit,'CCFHATADEL','UNKNOWN',' ','UNFORMATTED',
489     &            idummy,.false.)
490      read(lunit)
491      write(lunit) (work(kfckvaop-1+i),i=1,nviraop(isyfck))
492      call gpclose(lunit,'KEEP')
493
494*----------------------------------------------------------------------*
495*     return
496*----------------------------------------------------------------------*
497      call qexit('ccsdr12ao')
498      return
499      end
500*----------------------------------------------------------------------*
501*                     END OF SUBROUTINE CCSDR12AO                      *
502*======================================================================*
503