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*=====================================================================*
20      SUBROUTINE  CCFBINT1(ITRAN,   LISTL,   IDLSTL, LISTR, IDLSTR,
21     &                     XAINT,   YAINT,   MAINT,
22     &                     ZETA1,   ZETA2,   ISYCTR,
23     &                     TA1AMP,  TA2AMP,  ISYMTA,
24     &                     XLAMDP,  XLAMDH,  ISYLAM,
25     &                     XLAMDPQ, XLAMDHQ, ISYHOP,
26     &                     BZDENA,  BZDENAB,
27     &                     FILPQAMO,LUPQAMO,  IADRPQAMO,IADRPQA,
28     &                     FILPQA0, LUPQA0,   IADRPQA0, IADRPQAI0,
29     &                     FILPQA1, LUPQA1,   IADRPQA1, IADRPQAI1,
30     &                     LRELAX,  LTWOEL,  WORK,    LWORK)
31*---------------------------------------------------------------------*
32* Purpose:
33*
34*     Precalculate some intermediates for FbTa result vector depending
35*     on ZETA, T^A and IOPER:
36*     -- the XA and YA intermediates (in memory, allocated outside)
37*     -- the effective density BZDENA for the rho^BZA intermediate (disk)
38*     -- the effective density BZDENAB for the rho^BZQA intermediate (memory)
39*     -- the PA and QA intermediates, both with 4 MO indices and
40*                                     with 1 index in AO (0 and bar)
41*
42* Note: flags LRELAX & LTWOEL not yet carried through...?
43*
44*     Sonia Coriani, February 1999 (based on CCETAINT1)
45*     Version: 08/10-1999
46*---------------------------------------------------------------------*
47      IMPLICIT NONE
48#include "priunit.h"
49#include "ccsdinp.h"
50#include "ccsdsym.h"
51#include "maxorb.h"
52#include "ccorb.h"
53
54      LOGICAL LOCDBG
55      PARAMETER (LOCDBG = .FALSE.)
56
57      LOGICAL LRELAX, LTWOEL
58      CHARACTER*(*) LISTL, LISTR
59c
60      CHARACTER*(*) FILPQAMO, FILPQA0, FILPQA1
61c
62      INTEGER ITRAN, IDLSTL, ISYCTR, IDLSTR, ISYMTA
63      INTEGER ISYLAM, ISYHOP, LWORK
64c
65      INTEGER LUPQAMO, IADRPQA,   IADRPQAMO(MXCORB_CC,*)
66      INTEGER LUPQA0,  IADRPQAI0, IADRPQA0(MXCORB_CC,*)
67      INTEGER LUPQA1,  IADRPQAI1, IADRPQA1(MXCORB_CC,*)
68
69#if defined (SYS_CRAY)
70      REAL XAINT(*), YAINT(*), MAINT(*), ZETA1(*), ZETA2(*)
71      REAL BZDENA(*), BZDENAB(*)
72      REAL XLAMDP(*), XLAMDH(*), XLAMDPQ(*), XLAMDHQ(*)
73      REAL TA1AMP(*), TA2AMP(*), WORK(*)
74      REAL DUMMY, ONE, TWO
75#else
76      DOUBLE PRECISION XAINT(*), YAINT(*), MAINT(*), ZETA1(*), ZETA2(*)
77      DOUBLE PRECISION BZDENA(*), BZDENAB(*)
78      DOUBLE PRECISION XLAMDP(*), XLAMDH(*), XLAMDPQ(*), XLAMDHQ(*)
79      DOUBLE PRECISION TA1AMP(*), TA2AMP(*), WORK(*)
80      DOUBLE PRECISION DUMMY, ONE, TWO
81#endif
82      PARAMETER (ONE = 1.0D0, TWO = 2.0D0)
83
84      CHARACTER MODEL*(10)
85      INTEGER IOPT, ISYINT, ISYINTQ, KCHI, KCHIQ, KEND1, LWRK1, IDEL
86      INTEGER IDUMMY, KZ1A, KCHIA
87      INTEGER ILLSTOLD
88      SAVE    ILLSTOLD
89      INTEGER IRLSTOLD
90      SAVE    IRLSTOLD
91
92*---------------------------------------------------------------------*
93* do some tests and set symmetries:
94*---------------------------------------------------------------------*
95* if ITRAN=1, make sure that everything will be initialzed:
96*
97      IF (ITRAN.EQ.1) THEN
98         ILLSTOLD = IDLSTL - 1
99         IRLSTOLD = IDLSTR - 1
100      END IF
101
102* set symmetries for intermediates:
103
104      ISYINT  = MULD2H(ISYCTR,ISYMTA)     !ZxTA
105      ISYINTQ = MULD2H(ISYINT,ISYHOP)     !ZxTAxOperB
106
107*---------------------------------------------------------------------*
108* if only left vector (IDLSTL) has changed,
109*    read new multipliers Z1+Z2 into memory and square Z2 up:
110*---------------------------------------------------------------------*
111      IF ((IDLSTL .NE. ILLSTOLD).AND.(IDLSTR.EQ.IRLSTOLD)) THEN
112
113         IF (LWORK .LT. NT2AM(ISYCTR)) THEN
114            CALL QUIT('Insufficient memory in CCFBINT1 (a1)')
115         END IF
116
117         IOPT = 3
118         CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL,ZETA1,WORK)
119
120         CALL CC_T2SQ(WORK,ZETA2,ISYCTR)
121
122         IF (LOCDBG) THEN
123            WRITE (LUPRI,*) 'CCFBINT1> Only NEW Zeta vector'
124            WRITE (LUPRI,*) 'CCFBINT1> the zeta2 vector:'
125            CALL CC_PRP(WORK,WORK,ISYCTR,0,1)
126            WRITE (LUPRI,*) 'CCFBINT1> the tA2amp vector:'
127            CALL CC_PRP(WORK,TA2AMP,ISYMTA,0,1)
128         END IF
129      END IF
130*---------------------------------------------------------------------*
131* if only right vector (IDLSTR) has changed,
132*    read new response amplitudes TA1+TA2 into memory:
133*---------------------------------------------------------------------*
134
135      IF ((IDLSTL .EQ. ILLSTOLD).AND.(IDLSTR.NE.IRLSTOLD)) THEN
136
137         IF (LWORK .LT. NT2AM(ISYMTA)) THEN
138            CALL QUIT('Insufficient memory in CCFBINT1 (a2)')
139         END IF
140
141         IOPT = 3
142         CALL CC_RDRSP(LISTR,IDLSTR,ISYMTA,IOPT,MODEL,TA1AMP,TA2AMP)
143
144         IF (LOCDBG) THEN
145            WRITE (LUPRI,*) 'CCFBINT1> Only NEW T^A vector'
146            WRITE (LUPRI,*) 'CCFBINT1> the zeta2 vector:'
147            CALL CC_PRP(WORK,WORK,ISYCTR,0,1)
148            WRITE (LUPRI,*) 'CCFBINT1> the tA2amp vector:'
149            CALL CC_PRP(WORK,TA2AMP,ISYMTA,0,1)
150         END IF
151*
152* scale with 2 the tA_2 part
153*
154         IF ( (IOPT.EQ.3) .AND. (.NOT. LISTR(1:2).EQ.'R0') ) THEN
155            CALL CCLR_DIASCL(TA2AMP,TWO,ISYMTA)
156         END IF
157*
158      END IF
159*---------------------------------------------------------------------*
160* if both vectors (IDLSTR,IDLSTL) have changed,
161*                 read new Zeta and response amplitudes into memory:
162*---------------------------------------------------------------------*
163
164      IF ((IDLSTL .NE. ILLSTOLD).AND.(IDLSTR.NE.IRLSTOLD)) THEN
165
166c         IF (LWORK .LT. (NT2AM(ISYCTR)+NT2AM(ISYMTA))) THEN
167         IF (LWORK .LT. NT2AM(ISYCTR)) THEN
168            CALL QUIT('Insufficient memory in CCFBINT1 (a3)')
169         END IF
170
171         IOPT = 3  !(both singles and doubles)
172         CALL CC_RDRSP(LISTR,IDLSTR,ISYMTA,IOPT,MODEL,TA1AMP,TA2AMP)
173
174c         KSCR  = 1
175c         KEND1 = KSCR + NT2AM(ISYCTR)
176
177         IOPT = 3
178         CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL,ZETA1,WORK)
179
180         CALL CC_T2SQ(WORK,ZETA2,ISYCTR)
181
182         IF (LOCDBG) THEN
183            WRITE (LUPRI,*) 'CCFBINT1> Both NEW T^A and Zeta vectors'
184            WRITE (LUPRI,*) 'CCFBINT1> the zeta2 vector:(packed)'
185            CALL CC_PRP(WORK,WORK,ISYCTR,0,1)
186            WRITE (LUPRI,*) 'CCFBINT1> the tA2amp vector:'
187            CALL CC_PRP(WORK(1),TA2AMP,ISYMTA,0,1)
188         END IF
189*
190* scale with 2 the tA_2 part
191*
192         IF ( (IOPT.EQ.3) .AND. (.NOT. LISTR(1:2).EQ.'R0') ) THEN
193            CALL CCLR_DIASCL(TA2AMP,TWO,ISYMTA)
194         END IF
195*
196       END IF
197*---------------------------------------------------------------------*
198* if either left or right vector changed, calculate
199* new XA, YA and MA intermediates:
200*---------------------------------------------------------------------*
201      IF (( IDLSTL .NE. ILLSTOLD ).OR.( IDLSTR .NE. IRLSTOLD )) THEN
202
203          CALL CC_XI(XAINT,ZETA2,ISYCTR,TA2AMP,ISYMTA,WORK,LWORK)
204
205          CALL CC_YI(YAINT,ZETA2,ISYCTR,TA2AMP,ISYMTA,WORK,LWORK)
206
207          IF (CCSD .AND. LRELAX) THEN
208             CALL CC_MI(MAINT,ZETA2,ISYCTR,TA2AMP,ISYMTA,WORK,LWORK)
209          END IF
210
211          IF (LOCDBG) THEN
212             WRITE (LUPRI,'(//A)') 'CCFBINT1> XA-intermediate:'
213             WRITE (LUPRI,'(5G15.6)') (XAINT(I),I=1,NMATIJ(ISYINT))
214             WRITE (LUPRI,'(//A)') 'CCFBINT1> YA-intermediate:'
215             WRITE (LUPRI,'(5G15.6)') (YAINT(I),I=1,NMATAB(ISYINT))
216             IF (CCSD .AND. LRELAX) THEN
217                WRITE (LUPRI,'(//A)') 'CCFBINT1> MA-intermediate:'
218                WRITE (LUPRI,'(5G15.6)') (MAINT(I),I=1,N3ORHF(ISYINT))
219             END IF
220          END IF
221
222      END IF
223
224*---------------------------------------------------------------------*
225* calculate Chi^A matrices (calN)=(breve-Zeta1 - XA) intermediate
226* transform Zeta^ci to ZetaA^ki with T^A_ck
227* (cheap, do always!)
228*---------------------------------------------------------------------*
229      IF (CCSD) THEN
230        KCHIA = 1
231        KEND1 = KCHIA + NMATIJ(ISYINT)
232        LWRK1 = LWORK - KEND1
233
234        IF (LWRK1 .LT. NT2AM(ISYCTR)) THEN
235          CALL QUIT('Insufficient memory in CCFBINT1 (a)')
236        END IF
237        CALL CCLT_Z1A(ZETA1,ISYCTR,TA1AMP,ISYMTA, ISYINT,WORK(KCHIA) )
238        IF (LOCDBG) THEN
239          WRITE (LUPRI,'(//A)') 'CCFBINT1> ChiA-intermediate 1:'
240          WRITE (LUPRI,'(5G15.6)') (WORK(KCHIA+I-1),I=1,NMATIJ(ISYINT))
241        END IF
242        CALL DAXPY(NMATIJ(ISYINT),-ONE,XAINT,1,WORK(KCHIA),1)
243        IF (LOCDBG) THEN
244          WRITE (LUPRI,'(//A)') 'CCFBINT1> ChiA-intermediate:'
245          WRITE (LUPRI,'(5G15.6)') (WORK(KCHIA+I-1),I=1,NMATIJ(ISYINT))
246        END IF
247      END IF
248
249*---------------------------------------------------------------------*
250* calculate the effective BZA density (only for new TA or ZETA)
251* D_{k\a}^{ij} (backtr. with Lambda^p) (does not depend on Lambda-bar)
252* D_{k\a}^{ij} -> BZDENA (in memory)
253*---------------------------------------------------------------------*
254      IF (CCSD .AND. LRELAX) THEN
255        IF ((IDLSTL .NE. ILLSTOLD).OR.(IDLSTR .NE. IRLSTOLD)) THEN
256
257         IF (LOCDBG) THEN
258           WRITE (LUPRI,*) 'CCFBINT1: Print LambdaP for symmetry 1'
259           CALL CC_PRLAM(XLAMDP,XLAMDH,ISYLAM)
260           WRITE (LUPRI,*) 'FBINT1: Print ZETA2 for symmetry 1'
261           CALL CC_PRSQ(ZETA1,ZETA2,ISYCTR,0,1)
262           CALL DZERO(BZDENA,NT2AOIJ(ISYINT))
263           CALL DZERO(WORK(KEND1),LWRK1)
264         END IF
265
266         CALL CC_BFDENF( ZETA2,  ISYCTR, MAINT, ISYINT,
267     *                   XLAMDP, ISYLAM, WORK(KCHIA),  ISYINT,
268     *                   TA1AMP, ISYMTA, BZDENA, WORK(KEND1), LWRK1)
269         IF (LOCDBG) THEN
270           WRITE (LUPRI,*) 'FBINT1: Print BZDENA for symmetry 1'
271           CALL OUTPUT(BZDENA,1,NT1AO(1),1,NMATIJ(1),
272     *                          NT1AO(1),NMATIJ(1),1,LUPRI)
273           WRITE (LUPRI,*) 'returned 1 from CC_BFDENF... '
274           CALL FLSHFO(LUPRI)
275         END IF
276        ENDIF
277      END IF
278
279*---------------------------------------------------------------------*
280* calculate the effective BZAB density which does depend on the
281* external perturbation for _each_ transformation (ITRAN, as it
282* depends on IOPER): (or pass IOPER?)
283* 2) Dbar_{k\a}^{ij} -> BZDENAB (send in LambdapQ zero)
284*---------------------------------------------------------------------*
285      IF (CCSD .AND. LRELAX) THEN
286         CALL CC_BFDENF( ZETA2,  ISYCTR, MAINT, ISYINT,
287     *                   XLAMDPQ, ISYHOP, WORK(KCHIA),  ISYINT,
288     *                   TA1AMP,  ISYMTA, BZDENAB, WORK(KEND1), LWRK1)
289         IF (LOCDBG) THEN
290           WRITE (LUPRI,*) 'CCFBINT1: Print BZDENQA for symmetry 1'
291           CALL OUTPUT(BZDENAB,1,NT1AO(1),1,NMATIJ(1),
292     *               NT1AO(1),NMATIJ(1),1,LUPRI)
293           WRITE (LUPRI,*) 'returned 2 from CC_BFDENF... '
294           CALL FLSHFO(LUPRI)
295         END IF
296      ENDIF
297*---------------------------------------------------------------------*
298* calculate one-index backtransformed PA and QA intermediates used
299* in the F and G terms contribution to the result vector:
300*---------------------------------------------------------------------*
301      IF (CCSD .AND. LRELAX) THEN
302
303        IF (( IDLSTL .NE. ILLSTOLD ).OR.( IDLSTR .NE. IRLSTOLD )) THEN
304
305          IOPT = 1
306          CALL CC_PQIMO(ZETA2,ISYCTR,TA2AMP,ISYMTA,YAINT,IOPT,
307     *                  FILPQAMO,LUPQAMO,IADRPQAMO(1,ITRAN),IADRPQA,
308     *                  ITRAN,WORK(KEND1),LWRK1)
309
310          CALL CC_PQIAO(FILPQAMO,LUPQAMO,IADRPQAMO(1,ITRAN),ISYINT,
311     *                  FILPQA0, LUPQA0, IADRPQA0(1,ITRAN), IADRPQAI0,
312     *                  ITRAN, XLAMDH, ISYLAM,WORK(KEND1), LWRK1)
313
314
315         ELSE
316          DO IDEL = 1, NBAST
317            !dimensioned bigger that Nvir
318            IADRPQAMO(IDEL,ITRAN) = IADRPQAMO(IDEL,ITRAN-1)
319            IADRPQA0(IDEL,ITRAN)  = IADRPQA0(IDEL,ITRAN-1)
320          END DO
321        END IF
322
323C Backtransformation with Lambda-bar always performed
324c this is not necessary. I could pass LNEWOP!!!!!!!!!!!!!
325
326        CALL CC_PQIAO(FILPQAMO,LUPQAMO,IADRPQAMO(1,ITRAN),ISYINT,
327     *                FILPQA1, LUPQA1, IADRPQA1(1,ITRAN), IADRPQAI1,
328     *                ITRAN,XLAMDHQ,ISYHOP,WORK(KEND1), LWRK1)
329
330      END IF
331*---------------------------------------------------------------------*
332* save present IDLSTL in IDLSTOLD and return:
333*---------------------------------------------------------------------*
334      ILLSTOLD = IDLSTL
335      IRLSTOLD = IDLSTR
336
337      RETURN
338
339      END
340*=====================================================================*
341