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!
18C
19      SUBROUTINE CCETACOR(ETAAI,ZKDIA,WORK,LWORK)
20C
21C     Purpose: to calculate the correction to the RHS (eta_ai and eta_aI)
22C     of the z-kappa-bar0_ai equations originating from the
23C     "diagonal" z-kappa-bar0_ij and z-kappa-bar0_ab blocks
24C     over a loop over delta
25C
26C     Written by Sonia Coriani, Fall 2001. Based on CC2_DEN
27C
28C     Current models: CCD, CCSD, CCSD(T)
29C
30#include "implicit.h"
31#include "priunit.h"
32#include "dummy.h"
33#include "maxorb.h"
34#include "maxash.h"
35#include "mxcent.h"
36#include "aovec.h"
37#include "iratdef.h"
38      PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0)
39      PARAMETER (FOUR = 4.0D0)
40      DIMENSION INDEXA(MXCORB_CC)
41      DIMENSION ETAAI(*), ZKDIA(*), WORK(LWORK)
42#include "ccorb.h"
43#include "ccisao.h"
44#include "r12int.h"
45#include "inftap.h"
46#include "blocks.h"
47#include "ccfield.h"
48#include "ccsdinp.h"
49#include "ccinftap.h"
50#include "ccsdsym.h"
51#include "ccsdio.h"
52#include "distcl.h"
53#include "cbieri.h"
54#include "eritap.h"
55#include "ccfro.h"
56C
57      CHARACTER MODEL*10
58      CHARACTER NAME1*8
59      CHARACTER NAME2*8
60
61      LOGICAL LOCDBG
62      PARAMETER (LOCDBG=.FALSE.)
63C
64      CALL QENTER('CCETACOR')
65
66C--------------------------------------------------------------
67C Initialize work space
68C--------------------------------------------------------------
69
70      KCMO = 1
71      KEND1 = KCMO + NLAMDS
72      LWRK1 = LWORK - KEND1
73C
74C----------------------------------------------------------
75C     Read MO-coefficients from interface file and reorder.
76C----------------------------------------------------------
77C
78      LUSIFC = -1
79      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',
80     *            IDUMMY,.FALSE.)
81      REWIND LUSIFC
82      CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
83      READ (LUSIFC)
84      READ (LUSIFC)
85      READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS)
86      CALL GPCLOSE (LUSIFC,'KEEP')
87C
88      CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1)
89
90      IF (FROIMP) THEN
91         KCMOF = KEND1
92         KEND1 = KCMOF + NLAMDS
93         LWRK1 = LWORK - KEND1
94         CALL CMO_ALL(WORK(KCMOF),WORK(KEND1),LWRK1)
95      END IF
96C
97C---------------------------------------------------------------
98C     Start the loop for the 2-electron integrals and densities.
99C---------------------------------------------------------------
100C
101      KENDS2 = KEND1
102      LWRKS2 = LWRK1
103C
104      IF (DIRECT) THEN
105         IF (HERDIR) THEN
106            CALL HERDI1(WORK(KEND1),LWRK1,IPRERI)
107         ELSE
108            KCCFB1 = KEND1
109            KINDXB = KCCFB1 + MXPRIM*MXCONT
110            KEND1  = KINDXB + (8*MXSHEL*MXCONT + 1)/IRAT
111            LWRK1  = LWORK  - KEND1
112            CALL ERIDI1(KODCL1,KODCL2,KODBC1,KODBC2,KRDBC1,KRDBC2,
113     *                  KODPP1,KODPP2,KRDPP1,KRDPP2,
114     *                  KFREE,LFREE,KEND1,WORK(KCCFB1),WORK(KINDXB),
115     *                  WORK(KEND1),LWRK1,IPRERI)
116            KEND1 = KFREE
117            LWRK1 = LFREE
118         ENDIF
119         NTOSYM = 1
120      ELSE
121         NTOSYM = NSYM
122      ENDIF
123C
124      KENDSV = KEND1
125      LWRKSV = LWRK1
126C
127      ICDEL1 = 0
128      DO 100 ISYMD1 = 1,NTOSYM
129C
130         IF (DIRECT) THEN
131            IF (HERDIR) THEN
132               NTOT = MAXSHL
133            ELSE
134               NTOT = MXCALL
135            ENDIF
136         ELSE
137            NTOT = NBAS(ISYMD1)
138         ENDIF
139C
140         DO 110 ILLL = 1,NTOT
141C
142C---------------------------------------------
143C           If direct calculate the integrals.
144C---------------------------------------------
145C
146            IF (DIRECT) THEN
147C
148               KEND1 = KENDSV
149               LWRK1 = LWRKSV
150C
151c              DTIME  = SECOND()
152               IF (HERDIR) THEN
153                  CALL HERDI2(WORK(KEND1),LWRK1,INDEXA,ILLL,NUMDIS,
154     &                        IPRINT)
155               ELSE
156                  CALL ERIDI2(ILLL,INDEXA,NUMDIS,0,0,
157     *                        WORK(KODCL1),WORK(KODCL2),
158     *                        WORK(KODBC1),WORK(KODBC2),
159     *                        WORK(KRDBC1),WORK(KRDBC2),
160     *                        WORK(KODPP1),WORK(KODPP2),
161     *                        WORK(KRDPP1),WORK(KRDPP2),
162     *                        WORK(KCCFB1),WORK(KINDXB),
163     *                        WORK(KEND1), LWRK1,IPRERI)
164               ENDIF
165c              DTIME  = SECOND() - DTIME
166c              TIMHE2 = TIMHE2   + DTIME
167C
168               KRECNR = KEND1
169               KEND1  = KRECNR + (NBUFX(0) - 1)/IRAT + 1
170               LWRK1  = LWORK  - KEND1
171               IF (LWRK1 .LT. 0) THEN
172                  CALL QUIT('Insufficient core in CC2_DEN')
173               END IF
174C
175            ELSE
176               NUMDIS = 1
177            ENDIF
178C
179C-----------------------------------------------------
180C           Loop over number of distributions in disk.
181C-----------------------------------------------------
182C
183            DO 120 IDEL2 = 1,NUMDIS
184C
185               IF (DIRECT) THEN
186                  IDEL  = INDEXA(IDEL2)
187                  IF (NOAUXB) THEN
188                     IDUM = 1
189                     CALL IJKAUX(IDEL,IDUM,IDUM,IDUM)
190                  END IF
191                  ISYMD = ISAO(IDEL)
192               ELSE
193                  IDEL  = IBAS(ISYMD1) + ILLL
194                  ISYMD = ISYMD1
195               ENDIF
196C
197C----------------------------------------
198C              Work space allocation two.
199C----------------------------------------
200C
201               ISYDEN = ISYMD
202               ISYDIS = MULD2H(ISYMD,ISYMOP)
203C
204               KXINT  = KEND1
205               KEND2  = KXINT  + NDISAO(ISYDIS)
206               LWRK2  = LWORK  - KEND2
207C
208               IF (LWRK2 .LT. 0) THEN
209                  WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:',KEND2
210                  CALL QUIT(
211     *         'Insufficient space for allocation 2 in CCETAAICORR')
212               ENDIF
213C
214C--------------------------------------------
215C              Read AO integral distribution.
216C--------------------------------------------
217C
218               AUTIME = SECOND()
219               CALL CCRDAO(WORK(KXINT),IDEL,IDEL2,WORK(KEND2),LWRK2,
220     *                     WORK(KRECNR),DIRECT)
221C
222C--------------------------------------------------------------
223C              Calculate correction terms to eta_ai originating
224C              from kbar_ij and kbar_ab.
225C--------------------------------------------------------------
226C
227               KDSRHF = KEND2
228               K3OINT = KDSRHF + NDSRHF(ISYMD)
229               KEND3  = K3OINT + NMAIJK(ISYDIS)
230               IF (FROIMP) THEN
231                  KDSFRO = KEND3
232                  KOFOIN = KDSFRO + NDSFRO(ISYDIS)
233                  KEND3  = KOFOIN + NOFROO(ISYDIS)
234               ENDIF
235               LWRK3  = LWORK  - KEND3
236C
237               IF (LWRK3 .LT. 0) THEN
238                  WRITE(LUPRI,*) 'Need : ',KEND3,'Available : ',LWORK
239                  CALL QUIT(
240     *              'Insufficient core for integrals in CCETAC')
241               ENDIF
242C
243C--------------------------------------------------------------------
244C              Transform one index in the integral batch to occupied.
245C--------------------------------------------------------------------
246C
247               CALL CCTRBT(WORK(KXINT),WORK(KDSRHF),WORK(KCMO),
248     *                        ISYMOP,WORK(KEND3),LWRK3,ISYDIS)
249csonia
250C
251C------------------------------------------------------------------
252C              Transform one index in the integral batch to frozen.
253C------------------------------------------------------------------
254C
255               IF (FROIMP) THEN
256C
257                  CALL CC_GTOFRO(WORK(KXINT),WORK(KDSFRO),WORK(KCMOF),
258     *                           WORK(KEND3),LWRK3,ISYDIS)
259C
260C--------------------------------------------------------------
261C                 Calculate integral batch (cor fro | cor del).
262C--------------------------------------------------------------
263C
264                  CALL CC_OFROIN(WORK(KDSRHF),WORK(KOFOIN),WORK(KCMOF),
265     *                           WORK(KEND3),LWRK3,ISYDIS)
266C
267C-------------------------------------------------------------------------
268C                 Calculate correction to eta_aI from kappabar_ab
269C-------------------------------------------------------------------------
270C
271                  if (.true.) then
272                  KAFROI = 1 + NMATIJ(1) + NMATAB(1) + 2*NT1AMX
273                  CALL MP2_EIDV1(ZKDIA(KAFROI),WORK(KDSFRO),
274     *                           ZKDIA(NMATIJ(1)+1),WORK(KCMOF),
275     *                           WORK(KEND3),LWRK3,IDEL,ISYMD)
276C
277                  CALL MP2_EIDV2(ZKDIA(KAFROI),WORK(KDSFRO),
278     *                           ZKDIA(NMATIJ(1)+1),WORK(KCMOF),
279     *                           WORK(KEND3),LWRK3,IDEL,ISYMD)
280C
281C-------------------------------------------------------------------------
282C                 Calculate correction to eta_aI from kappabar_ij
283C-------------------------------------------------------------------------
284C
285                  CALL MP2_EIDC1(ZKDIA(KAFROI),WORK(KDSFRO),
286     *                           ZKDIA(1),WORK(KCMOF),WORK(KEND3),
287     *                           LWRK3,IDEL,ISYMD)
288C
289                  CALL MP2_EIDC2(ZKDIA(KAFROI),WORK(KOFOIN),
290     *                           ZKDIA(1),WORK(KCMOF),WORK(KEND3),
291     *                           LWRK3,IDEL,ISYMD)
292
293                  end if
294               END IF
295csonia
296C
297C-------------------------------------------------------------------
298C              Calculate integral batch with three occupied indices.
299C              Since LUDUMLOCAL < 0 the integrals are not written
300C              to disk.
301C-------------------------------------------------------------------
302C
303               LUDUMLOCAL = -10
304               CALL CC_INT3O(WORK(K3OINT),WORK(KDSRHF),WORK(KCMO),
305     *                         ISYMOP,WORK(KCMO),WORK(KEND3),LWRK3,
306     *                         IDEL,ISYMD,LUDUMLOCAL,'DUMMY')
307C
308C-------------------------------------------------------
309C              Calculate the correction terms to eta_ai
310C              from kappabar_ab and kappabar_ij
311C-------------------------------------------------------
312C
313               CALL MP2_YTV(ETAAI,ZKDIA(NMATIJ(1)+1),WORK(KDSRHF),
314     *                      WORK(KCMO),WORK(KEND3),LWRK3,IDEL,ISYMD)
315C
316               CALL MP2_NXY(ETAAI,ZKDIA(1),ZKDIA(NMATIJ(1)+1),
317     *                      WORK(K3OINT),WORK(KDSRHF),WORK(KCMO),
318     *                      WORK(KEND3),LWRK3,IDEL,ISYMD)
319C
320               CALL MP2_XTO(ETAAI,ZKDIA(1),WORK(K3OINT),
321     *                      WORK(KCMO),WORK(KEND3),LWRK3,IDEL,ISYMD)
322C
323               AUTIME = SECOND() - AUTIME
324C
325  120       CONTINUE
326  110    CONTINUE
327  100 CONTINUE
328
329C-----------------------------------------------------
330C     Some print out
331C-----------------------------------------------------
332      IF (LOCDBG) THEN
333         ETAKAN = DDOT(NALLAI(1),ETAAI,1,ETAAI,1)
334         WRITE(LUPRI,*) 'CCETACOR: eta_ai '
335         WRITE(LUPRI,*) 'Norm of corrected occ-vir block:', ETAKAN
336      ENDIF
337C
338C-----------------------
339C     Regain work space.
340C-----------------------
341C
342      KEND1 = KENDS2
343      LWRK1 = LWRKS2
344C
345C---------------------------------------------------------------
346      CALL QEXIT('CCETACOR')
347      RETURN
348      END
349