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
19C  /* Deck cc_21i3 */
20*--------------------------------------------------------------
21      SUBROUTINE CC_21I3(RHO1, XINT0, ISYDIS0, XINT1, ISYDIS1,
22     *                   PINT0,QINT0,ISYPQ0,PINT1,QINT1,ISYPQ1,
23     *                   PAINT0,QAINT0,ISYPQA0,
24     *                   PAINT1,QAINT1,ISYPQA1,
25     *                   XLAMP0,XLAMH0, ISYML0,XLAMPQ,ISYMLQ,
26     *                   XLAMPA,ISYMLA,XLAMPQA,ISYMLQA,
27     *                   WORK,LWORK,IOPT,LRHOAO,LDERIV,LXI1SQ)
28*--------------------------------------------------------------
29*
30*     Generalization of the CC_21I and CC_21I2 routines
31*     for the F contribution (21I contr) to the F^bT^A result vector
32*     (Calculation of the FQA intermediate)
33*
34*     Sonia Coriani February 1999
35*
36*     Version: 2.0
37*
38*     XINT0, ISYDIS0  -- (* *|* delta) batch of integrals, its symmetry
39*     XINT1, ISYDIS1 -- the same for derivative integrals
40*     PINT0, QINT0   -- P_ci,j (del), Q_ci,j (del)      (ISYPQ0)
41*     PINT1, QINT1   -- Pbar_ci,j (del), Qbar_ci,j(del) (ISYPQ1)
42*     PAINT0,QAINT0  -- PA_ci,j (del) QA_ci,j(del)      (ISYPQA0)
43*     PAINT1,QAINT1  -- PAbar_ci,j QAbar_ci,j (del)     (ISYPQA1)
44*
45*     LRHOAO = .TRUE.    return result with a index in AO
46*     LRHOAO = .FALSE.   return result with a index in MO
47*     LDERIV = .FALSE.   omit contribution from derivative integrals
48*     LXI1SQ = .TRUE.    (a b|* *) of XINT1 is a full matrix (fx LAO)
49*
50*     IOPT = 1   The FA intermediate (test case)
51*     IOPT = 2   The FQA intermediate
52*------------------------------------------------------------------
53*
54#if defined (IMPLICIT_NONE)
55      IMPLICIT NONE
56#else
57#   include "implicit.h"
58#endif
59#include "priunit.h"
60#include "ccorb.h"
61#include "ccsdsym.h"
62C
63      LOGICAL LRHOAO, LDERIV, LXI1SQ
64      INTEGER LWORK,ISYDIS0,ISYDIS1,IOPT
65      INTEGER ISYPQ0,ISYPQ1,ISYPQA0,ISYPQA1, ip
66      INTEGER ISYML0, ISYMLQ, ISYMLA, ISYMLQA, ISYM
67C
68#if defined (SYS_CRAY)
69      REAL RHO1(*), PINT0(*), PINT1(*), QINT0(*), QINT1(*)
70      REAL PAINT0(*), PAINT1(*), QAINT0(*), QAINT1(*)
71      REAL XINT0(*), XINT1(*), WORK(LWORK)
72      REAL XLAMP0(*), XLAMPA(*), XLAMPQ(*), XLAMPQA(*)
73      REAL XLAMH0(*)
74      REAL DUMMY, ZERO, ONE, TWO, DDOT, XNORM
75#else
76      DOUBLE PRECISION RHO1(*), PINT0(*), PINT1(*), QINT0(*), QINT1(*)
77      DOUBLE PRECISION PAINT0(*), PAINT1(*), QAINT0(*), QAINT1(*)
78      DOUBLE PRECISION XINT0(*), XINT1(*), WORK(LWORK)
79      DOUBLE PRECISION XLAMP0(*), XLAMPA(*), XLAMPQ(*), XLAMPQA(*)
80      DOUBLE PRECISION XLAMH0(*)
81      DOUBLE PRECISION DUMMY, ZERO, ONE, TWO, DDOT, XNORM
82#endif
83C
84      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
85C
86      INTEGER NT3BAO(8), IT3BAO(8,8)
87      INTEGER ISYGIJ1, ISYGIJ2, ISYDENBA, ISYDENA
88      INTEGER ISYRES, ISYMJ, ISYMGI, ISYMAL, ISY3AO, ISYMG
89      INTEGER ISYALI, ISYMI, ISALBE, ISYMBE, ISYMA, ISYGAM
90      INTEGER KRHOAO, KPQMO, KPQONE, KPQTWO, KDENBA, KDENA
91      INTEGER KAOINT,NTOTA,NPQMO, NPQMO01, NPQMOA1, NTOTGA
92      INTEGER KEND1,LWRK1,KOFF1,KOFF2,KOFF3,KOFF4,KOFF5
93      INTEGER ICOUNT, NTOTAL, NTOTGI, NGI, NAGI, NAIG, NTOTBE, NTOTAJ
94      INTEGER IOPTDEN, IOPTCON
95C
96C     set some symmetries used globally in this routine:
97C
98      ISYGIJ1   = MULD2H(ISYPQ0,ISYMLA)
99*      IF (IOPT.EQ.2) THEN
100*      ISYGIJ2   = MULD2H(ISYPQ0,ISYMLQA)
101*      END IF
102      ISYGIJ2   = MULD2H(ISYPQA0,ISYMLQ)
103      ISYDENBA  = MULD2H(ISYGIJ1,ISYMLQ)
104      ISYDENA   = MULD2H(ISYGIJ1,ISYML0)
105      ISYRES    = MULD2H(ISYDENBA,ISYDIS0)
106C
107*      IF ( IOPT.GE.2 .AND. MULD2H(ISYVWZ2,ISYXL0).NE.ISYGIJ ) THEN
108*            CALL QUIT('Symmetry mismatch in CC_21I3NEW.')
109*      END IF
110C
111C     calculate index arrays for intermediates with 3 indices in AO:
112C
113      DO ISY3AO = 1, NSYM
114        ICOUNT = 0
115        DO ISYMG = 1, NSYM
116           ISYALI = MULD2H(ISY3AO,ISYMG)
117           IT3BAO(ISYALI,ISYMG) = ICOUNT
118           ICOUNT = ICOUNT + NBAS(ISYMG) * NT1AO(ISYALI)
119        END DO
120        NT3BAO(ISY3AO) = ICOUNT
121      END DO
122C
123C----------------------------------------------------------------
124C     allocate work space for
125C          --  intermediate with one index backtransformed
126C          --  intermediate with one more index backtransformed
127C          --  intermediate with two more indices backtransformed
128C----------------------------------------------------------------
129C
130      NPQMO = 1
131      DO ISYM = 1, NSYM
132        NPQMO = MAX(NPQMO,NT2BCD(ISYM))
133      END DO
134
135C
136C      IF (IOPT.GE.2) NPQMO = MAX(NPQMO,NT2BCD(ISYVWZ2))
137C
138      KRHOAO = 1
139      KPQMO  = KRHOAO + NT1AO(ISYRES)
140      KPQONE = KPQMO  + NPQMO
141      IF (IOPT.EQ.2) THEN
142        KPQTWO = KPQONE + NT2BGD(ISYGIJ1)
143      ELSE
144        KPQTWO = KPQONE
145*        ISYGIJ2 = ISYGIJ1
146      END IF
147
148      KDENBA = KPQTWO + NT2BGD(ISYGIJ2)
149      KEND1  = KDENBA + NT3BAO(ISYDENBA)
150
151      LWRK1  = LWORK  - KEND1
152
153      IF ((IOPT.EQ.2).AND.LDERIV) THEN
154         KDENA = KEND1
155         KEND1 = KDENA + NT3BAO(ISYDENA)
156      END IF
157      LWRK1   = LWORK - KEND1
158C
159      IF ( LWRK1 .LT. 0) THEN
160         CALL QUIT('Insufficient work space in CC_21I3.')
161      ENDIF
162C
163C--------------------------------------------------------------------
164C     1. Exchange like contribution: add P and Q intermediates
165C     together and transform virtual index "c" back to AO
166C     with the LambdaA^p matrix
167C     add PA and QA intermediates together and backtransform
168C     "c" with the Lambda^p matrix. Add to previous
169C     Repeat for Pbar (P1) and Qbar (Q1)
170C     Repeat for PAbar (PA1) and QAbar (QA1)
171C--------------------------------------------------------------------
172C
173      CALL DCOPY(NT2BCD(ISYPQ0),PINT0,1,WORK(KPQMO),1)
174      CALL DAXPY(NT2BCD(ISYPQ0),ONE,QINT0,1,WORK(KPQMO),1)
175C
176      !LambdaA*PQ -> KPQONE
177      CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQ0,XLAMPA,ISYMLA,
178     &             1,1,DUMMY,1,DUMMY,1,ZERO,0)
179c
180c      XNORM = DDOT(NT2BGD(ISYGIJ1),WORK(KPQONE),1,WORK(KPQONE),1)
181c      WRITE (LUPRI,*) 'CC_21I3> Norm(PQ-AO)1=',XNORM
182c
183      IF (IOPT.EQ.2) THEN
184        !LambdaQA*PQ -> KPQTWO
185        CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ0,XLAMPQA,ISYMLQA,
186     &                    1,1,DUMMY,1,DUMMY,1,ZERO,0)
187      END IF
188
189      CALL DCOPY(NT2BCD(ISYPQA0),PAINT0,1,WORK(KPQMO),1)
190c      WRITE(LUPRI,*) 'CC_21I3> test PQ'
191c      do ip = 1, 5
192c       WRITE(LUPRI,*) PAINT0(ip)
193c      end do
194      CALL DAXPY(NT2BCD(ISYPQA0),ONE,QAINT0,1,WORK(KPQMO),1)
195
196      !Lambda0*calPQ -> + KPQONE
197      CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQA0,XLAMP0,ISYML0,
198     &             1,1,DUMMY,1,DUMMY,1,ONE,0)
199c
200c      XNORM = DDOT(NT2BGD(ISYGIJ1),WORK(KPQONE),1,WORK(KPQONE),1)
201c      WRITE (LUPRI,*) 'CC_21I3> Norm(PQA-AO)1=',XNORM
202c
203
204      IF (IOPT.EQ.2) THEN
205        !LambdaQ*calPQ -> + KPQTWO
206        CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA0,XLAMPQ,ISYMLQ,
207     &               1,1,DUMMY,1,DUMMY,1,ONE,0)
208C--------------------------------------------------------------------
209C     Backtransform also barred PQ  and calPQ and add to KPQTWO
210C--------------------------------------------------------------------
211
212        CALL DCOPY(NT2BCD(ISYPQ1),PINT1,1,WORK(KPQMO),1)
213        CALL DAXPY(NT2BCD(ISYPQ1),ONE,QINT1,1,WORK(KPQMO),1)
214C
215        !LambdaA*barPQ -> + KPQTWO
216        CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ1,XLAMPA,ISYMLA,
217     &                1,1,DUMMY,1,DUMMY,1,ONE,0)
218
219        CALL DCOPY(NT2BCD(ISYPQA1),PAINT1,1,WORK(KPQMO),1)
220        !LambdaA*bar-calPQ -> + KPQTWO
221        CALL DAXPY(NT2BCD(ISYPQA1),ONE,QAINT1,1,WORK(KPQMO),1)
222
223        CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA1,XLAMP0,ISYML0,
224     &                1,1,DUMMY,1,DUMMY,1,ONE,0)
225      END IF
226C
227C--------------------------------------------------------------------
228C     transform occupied index j back to AO (alpha)
229C--------------------------------------------------------------------
230C
231      IOPTDEN = 0   !initialize density with actual contribution
232      !backt j with XLAMPQ
233      CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMPQ,ISYMLQ,WORK(KDENBA),
234     *               ISYDENBA,IT3BAO,WORK(KEND1),LWRK1,'EXC',IOPTDEN)
235
236      IF (IOPT.EQ.2) THEN
237        IOPTDEN = 1   !add to density
238        !backt j with XLAMP0
239        CALL CC_21IDEN(WORK(KPQTWO),ISYGIJ2,XLAMP0,ISYML0,
240     *                  WORK(KDENBA),ISYDENBA,IT3BAO,
241     *                  WORK(KEND1),LWRK1,'EXC',IOPTDEN)
242
243        IF (LDERIV) THEN
244          IOPTDEN = 0
245          CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMP0,ISYML0,
246     *                     WORK(KDENA),ISYDENA,IT3BAO,
247     *                     WORK(KEND1),LWRK1,'EXC',IOPTDEN)
248        END IF
249      END IF
250C
251c     XNORM = DDOT(NT3BAO(ISYGIJ1),WORK(KDENBA),1,WORK(KPQDEN),1)
252c     WRITE (LUPRI,*) 'Norm(PQDEN)1=',XNORM
253C
254C--------------------------------------------------------------------
255C     2. Coulomb like contribution: scale P intermediate with -2
256C     and transform virtual "c" back to AO
257C--------------------------------------------------------------------
258C
259      CALL DCOPY(NT2BCD(ISYPQ0),PINT0,1,WORK(KPQMO),1)
260      CALL DSCAL(NT2BCD(ISYPQ0),-TWO,WORK(KPQMO),1)
261
262      CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQ0,XLAMPA,ISYMLA,
263     &             1,1,DUMMY,1,DUMMY,1,ZERO,0)
264C
265c     XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1)
266c     WRITE (LUPRI,*) 'Norm(PQAO)3=',XNORM
267C
268C
269      IF (IOPT.EQ.2) THEN
270        CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ0,XLAMPQA,
271     &               ISYMLQA,1,1,DUMMY,1,DUMMY,1,ZERO,0)
272      END IF
273C
274
275      CALL DCOPY(NT2BCD(ISYPQA0),PAINT0,1,WORK(KPQMO),1)
276      CALL DSCAL(NT2BCD(ISYPQA0),-TWO,WORK(KPQMO),1)
277
278      CALL CC_BFAO(WORK(KPQONE),WORK(KPQMO),ISYPQA0,XLAMP0,ISYML0,
279     &                                  1,1,DUMMY,1,DUMMY,1,ONE,0)
280
281      IF (IOPT.EQ.2) THEN
282         CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA0,XLAMPQ,ISYMLQ,
283     &                                     1,1,DUMMY,1,DUMMY,1,ONE,0)
284
285C
286c     XNORM = DDOT(NT2BGD(ISYGIJ),WORK(KPQAO),1,WORK(KPQAO),1)
287c     WRITE (LUPRI,*) 'Norm(PQAO)4=',XNORM
288C
289
290         CALL DCOPY(NT2BCD(ISYPQ1),PINT1,1,WORK(KPQMO),1)
291         CALL DSCAL(NT2BCD(ISYPQ1),-TWO,WORK(KPQMO),1)
292
293         CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQ1,XLAMPA,ISYMLA,
294     &                                    1,1,DUMMY,1,DUMMY,1,ONE,0)
295
296         CALL DCOPY(NT2BCD(ISYPQA1),PAINT1,1,WORK(KPQMO),1)
297         CALL DSCAL(NT2BCD(ISYPQA1),-TWO,WORK(KPQMO),1)
298
299         CALL CC_BFAO(WORK(KPQTWO),WORK(KPQMO),ISYPQA1,XLAMP0,ISYML0,
300     &                                     1,1,DUMMY,1,DUMMY,1,ONE,0)
301
302      END IF
303C--------------------------------------------------------------------
304C     transform occupied index j back to AO (gamma), add result to
305C     the effective density stored in KPQDEN
306C--------------------------------------------------------------------
307C
308
309      IOPTDEN = 1
310      CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMPQ,ISYMLQ,WORK(KDENBA),
311     *             ISYDENBA,IT3BAO,WORK(KEND1),LWRK1,'COU',IOPTDEN)
312
313      IF (IOPT.EQ.2) THEN
314      IOPTDEN = 1
315      CALL CC_21IDEN(WORK(KPQTWO),ISYGIJ2,XLAMP0,ISYML0,
316     *                  WORK(KDENBA),ISYDENBA,IT3BAO,
317     *                  WORK(KEND1),LWRK1,'COU',IOPTDEN)
318
319      IF (LDERIV) THEN
320         IOPTDEN = 1
321         CALL CC_21IDEN(WORK(KPQONE),ISYGIJ1,XLAMP0,ISYML0,
322     *                     WORK(KDENA),ISYDENA,IT3BAO,
323     *                     WORK(KEND1),LWRK1,'COU',IOPTDEN)
324      END IF
325      END IF
326C
327C--------------------------------------------------------------------
328C     contract the effective density with the 2-el integrals:
329C--------------------------------------------------------------------
330C
331c     XNORM = DDOT(NT3BAO(ISYGIJ),WORK(KPQDEN),1,WORK(KPQDEN),1)
332c     WRITE (LUPRI,*) 'Norm(PQDEN)2=',XNORM
333C
334      CALL DZERO(WORK(KRHOAO),NT1AO(ISYRES))
335
336      CALL CC_21ICON(WORK(KRHOAO),ISYRES,WORK(KDENBA),ISYDENBA,
337     *               XINT0,ISYDIS0,IT3BAO,WORK(KEND1),LWRK1,1)
338
339      IF ((IOPT.EQ.2).AND.LDERIV) THEN
340
341         IF (LXI1SQ) THEN
342            IOPTCON = 2             !X2INT has full (a b| already
343         ELSE
344            IOPTCON = 1             !X2INT is squared inside CC_21ICON1
345         END IF
346
347         CALL CC_21ICON(WORK(KRHOAO),ISYRES,WORK(KDENA),ISYDENA,
348     *            XINT1,ISYDIS1,IT3BAO,WORK(KEND1),LWRK1,IOPTCON)
349      END IF
350C
351c     XNORM = DDOT(NT1AO(ISYRES),RHO1,1,RHO1,1)
352c     WRITE (LUPRI,*) 'XNORM=',XNORM
353c
354C---------------------------------------------
355C     transformation and/or storage in result:
356C---------------------------------------------
357C
358      IF (LRHOAO) THEN
359        CALL DAXPY(NT1AO(ISYRES),ONE,WORK(KRHOAO),1,RHO1,1)
360      ELSE
361        CALL CC_T1AM(RHO1,ISYRES,WORK(KRHOAO),ISYRES,XLAMH0,ISYML0,ONE)
362      END IF
363C
364      RETURN
365      END
366C
367*=====================================================================*
368