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 CC_FDD(NC1VEC,NC2VEC,LISTA,ITAMPA,LISTB,ITAMPB,
21     &                  LISTC,ITAMPC,TYAM,RESULT,WORK,LWORK)
22C
23C---------------------------------------------------------------------
24C Test routine for calculating the CC D matrix by finite difference on
25C the C matrix transformation.
26C Ch. Haettig, maj 1997, based on Oves CCLR_FDF routine
27C---------------------------------------------------------------------
28C
29#include "implicit.h"
30#include "priunit.h"
31#include "dummy.h"
32#include "maxorb.h"
33#include "iratdef.h"
34#include "ccorb.h"
35#include "aovec.h"
36#include "ccsdinp.h"
37#include "cclr.h"
38#include "ccsdsym.h"
39#include "ccsdio.h"
40#include "leinf.h"
41C
42      DIMENSION WORK(LWORK),ITADR(2),RESULT(*)
43      PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00, DELTA = 1.0D-07)
44      PARAMETER (ONE = 1.0d0, ZERO = 0.0d0, TWO = 2.0d0)
45      CHARACTER MODEL*10
46      CHARACTER FILCMA*8
47      LOGICAL L1TST,L2TST, LETST
48      INTEGER ICTRAN(4)
49C
50      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
51C
52      MODEL = 'CCSD      '
53      IF (CCS) MODEL = 'CCS       '
54      IF (CC2) MODEL = 'CC2       '
55C
56      IF (CCR12) CALL QUIT('Finite-difference D-matrix for CCR12 '//
57     &                     'not adapted')
58C
59      IF (IPRINT.GT.5) THEN
60         CALL AROUND( 'IN CC_FDD  : MAKING FINITE DIFF. CC D Matrix')
61      ENDIF
62C
63C----------------------------
64C     Work space allocations.
65C----------------------------
66C
67      ISYMTR     = 1
68      ISYMOP     = 1
69C
70      NTAMP      = NT1AM(ISYMTR) + NT2AM(ISYMTR)
71      NTAMP2     = NTAMP*(NC1VEC + NC2VEC )
72      KF         = 1
73      KRHO1      = KF       + NTAMP2
74      KRHO2      = KRHO1    + NT1AMX
75      KC1AM      = KRHO2    + MAX(NT2AMX,NT2AM(ISYMTR))
76      KC2AM      = KC1AM    + NT1AM(ISYMTR)
77      KEND1      = KC2AM
78     *           + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR),
79     *                 2*NT2ORT(ISYMTR))
80      LWRK1      = LWORK    - KEND1
81C
82      KRHO1D     = KEND1
83      KRHO2D     = KRHO1D   + NT1AMX
84      KEND2      = KRHO2D
85     *           + MAX(NT2AMX,NT2AM(ISYMTR),NT2AO(ISYMTR),
86     *                 2*NT2ORT(ISYMTR))
87      LWRK2      = LWORK      - KEND1
88C
89      IF (IPRINT .GT. 100 ) THEN
90         WRITE(LUPRI,*) ' IN CC_FDD: KF      =  ',KF
91         WRITE(LUPRI,*) ' IN CC_FDD: KRHO1   =  ',KRHO1
92         WRITE(LUPRI,*) ' IN CC_FDD: KRHO2   =  ',KRHO2
93         WRITE(LUPRI,*) ' IN CC_FDD: KC1AM   =  ',KC1AM
94         WRITE(LUPRI,*) ' IN CC_FDD: KC2AM   =  ',KC2AM
95         WRITE(LUPRI,*) ' IN CC_FDD: KRHO1D  =  ',KRHO1D
96         WRITE(LUPRI,*) ' IN CC_FDD: KRHO2D  =  ',KRHO2D
97         WRITE(LUPRI,*) ' IN CC_FDD: KEND2   =  ',KEND2
98         WRITE(LUPRI,*) ' IN CC_FDD: LWRK2   =  ',LWRK2
99      ENDIF
100      IF (LWRK2.LT.0 ) THEN
101         WRITE(LUPRI,*) 'Too little work space in CC_FDD '
102         WRITE(LUPRI,*) 'AVAILABLE: LWORK   =  ',LWORK
103         WRITE(LUPRI,*) 'NEEDED (AT LEAST)  =  ',KEND2
104         CALL QUIT('TOO LITTLE WORKSPACE IN CC_FDD ')
105      ENDIF
106      KF2   = KF      + NC1VEC*NTAMP
107C
108C---------------------
109C     Initializations.
110C---------------------
111C
112      CALL DZERO(WORK(KC1AM),NT1AMX)
113      CALL DZERO(WORK(KC2AM),NT2AMX)
114      CALL DZERO(WORK(KF),NTAMP2)
115      IF (ABS(DELTA) .GT. 1.0D-15 ) THEN
116         DELTAI = 1.0D00/DELTA
117      ELSE
118         DELTAI = 1
119      ENDIF
120      X11 = 0.0D00
121      X12 = 0.0D00
122      X21 = 0.0D00
123      X22 = 0.0D00
124      XNJ = 0.0D00
125C
126C------------------------------------------------
127C     Read the CC reference amplitudes From disk.
128C------------------------------------------------
129C
130      IOPT = 3
131      CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KC1AM),WORK(KC2AM))
132C
133C----------------------------------------------
134C     Save the CC reference amplitudes on disk.
135C----------------------------------------------
136C
137      LUTAM = -1
138      CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY,
139     *            .FALSE.)
140      REWIND(LUTAM)
141      WRITE(LUTAM) (WORK(KC1AM + I -1 ), I = 1, NT1AMX)
142      WRITE(LUTAM) (WORK(KC2AM + I -1 ), I = 1, NT2AMX)
143      CALL GPCLOSE(LUTAM,'KEEP')
144C
145      IF (IPRINT.GT.125) THEN
146         RHO1N = DDOT(NT1AMX,WORK(KC1AM),1,WORK(KC1AM),1)
147         RHO2N = DDOT(NT2AMX,WORK(KC2AM),1,WORK(KC2AM),1)
148         WRITE(LUPRI,*) 'Norm of T1AM: ',RHO1N
149         WRITE(LUPRI,*) 'Norm of T2AM: ',RHO2N
150         CALL CC_PRP(WORK(KC1AM),WORK(KC2AM),1,1,1)
151      ENDIF
152      RSPIM = .TRUE.
153C
154C------------------------------------------
155C     Calculate reference A*T vector.
156C------------------------------------------
157C
158      ICTRAN(1) = ITAMPA
159      ICTRAN(2) = ITAMPB
160      ICTRAN(3) = ITAMPC
161      ICTRAN(4) = 0
162
163      NCTRAN = 1
164      IOPT   = 1
165      FILCMA = 'CCCMAT'
166
167      CALL CC_CMAT(ICTRAN,NCTRAN,LISTA,LISTB,LISTC,IOPT,
168     &             FILCMA,IDUM,RDUM,0,WORK(KRHO1D),LWORK-KRHO1D)
169
170C
171C-------------------------
172C     Zero out components.
173C-------------------------
174C
175      IF (LCOR .OR. LSEC) THEN
176C
177         CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR)
178C
179      ENDIF
180C
181      IF (IPRINT.GT.2) THEN
182         RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1)
183         RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1)
184         WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ref'
185         WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ref'
186      ENDIF
187      IF (IPRINT.GT.125) THEN
188         CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1)
189      ENDIF
190
191      CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1),1)
192      CALL DCOPY(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2),1)
193C
194C=============================================
195C     Calculate C-matrix by finite difference.
196C=============================================
197C
198      DO 100 I = 1, NC1VEC
199         WRITE (LUPRI,*) 'singles index:',I
200C
201C----------------------------------------
202C        Add finite displadement to t and
203C        calculate new intermediates.
204C----------------------------------------
205C
206         LUTAM = -1
207         CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY,
208     *               .FALSE.)
209         READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX)
210         READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX)
211         CALL GPCLOSE(LUTAM,'KEEP')
212C
213         TI   = SECOND()
214         WORK(KC1AM +I -1) = WORK(KC1AM +I -1 ) + DELTA
215         IF (LCOR .OR. LSEC) THEN
216            CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR)
217         ENDIF
218C
219         IOPT = 3
220         CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM),
221     *                 WORK(KC2AM),WORK(KEND2),LWRK2)
222C
223         RSPIM = .TRUE.
224         CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM),
225     *               WORK(KC2AM),WORK(KEND2),LWRK2,'XXX')
226C
227C------------------
228C        Transform.
229C------------------
230C
231         ICTRAN(1) = ITAMPA
232         ICTRAN(2) = ITAMPB
233         ICTRAN(3) = ITAMPC
234         ICTRAN(4) = 0
235
236         NCTRAN = 1
237         IOPT   = 1
238         FILCMA = 'CCCMAT'
239
240         CALL CC_CMAT(ICTRAN,NCTRAN,LISTA,LISTB,LISTC,IOPT,
241     &                FILCMA,IDUM,RDUM,0,WORK(KRHO1D),LWORK-KRHO1D)
242
243
244         IF (LCOR .OR. LSEC) THEN
245            CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR)
246         ENDIF
247C
248         IF (IPRINT.GT.2) THEN
249            RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1)
250            RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1)
251            WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'ai=',I
252            WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'ai=',I
253         ENDIF
254         IF (IPRINT.GT.125) THEN
255            CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1)
256         ENDIF
257         CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1)
258         CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1)
259         CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1)
260         CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1)
261         CALL DCOPY(NT1AMX,WORK(KRHO1D),1,
262     *              WORK(KF+NTAMP*(I-1)),1)
263         CALL DCOPY(NT2AMX,WORK(KRHO2D),1,
264     *             WORK(KF+NTAMP*(I-1)+NT1AMX),1)
265         X11 = X11 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1)
266         X21 = X21 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1)
267C
268         TI   = SECOND() - TI
269         IF (IPRINT.GT.5 ) THEN
270            WRITE(LUPRI,*) '  '
271            WRITE(LUPRI,*) 'FDB ROW NR. ',I,' DONE IN ',TI,' SEC.'
272         ENDIF
273C
274 100  CONTINUE
275C
276C----------------------------------------------------------------
277C     Loop over T2 amplitudes. Take care of diagonal t2 elements
278C     is in a different convention in the energy code.
279C     Factor 1/2 from right , and factor 2 from left.
280C----------------------------------------------------------------
281C
282      DO 200 NAI = 1, NT1AMX
283        DO 300 NBJ = 1, NAI
284         I = INDEX(NAI,NBJ)
285C
286         IF (I.LE.NC2VEC) THEN
287           WRITE (LUPRI,*) 'doubles index:',I
288C
289C--------------------------------------------
290C          Add finite displacement to t and
291C          calculate new intermediates.
292C-------------------------------------------
293C
294           LUTAM = -1
295           CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',
296     *                 IDUMMY,.FALSE.)
297           READ(LUTAM) (WORK(KC1AM + J -1 ) , J = 1, NT1AMX)
298           READ(LUTAM) (WORK(KC2AM + J -1 ) , J = 1, NT2AMX)
299           CALL GPCLOSE(LUTAM,'KEEP')
300C
301           TI   = SECOND()
302           DELT = DELTA
303           IF (NAI.EQ.NBJ) DELT = 2*DELTA
304           WORK(KC2AM + I -1) = WORK(KC2AM+I -1) + DELT
305           IF (LCOR .OR. LSEC) THEN
306             CALL CC_CORE(WORK(KC1AM),WORK(KC2AM),ISYMTR)
307           ENDIF
308C
309           IOPT = 3
310           CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM),
311     *                   WORK(KC2AM),WORK(KEND2),LWRK2)
312C
313           RSPIM = .TRUE.
314           CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM),
315     *                 WORK(KC2AM),WORK(KEND2),LWRK2,'XXX')
316C
317C--------------------
318C          Transform.
319C--------------------
320C
321           ICTRAN(1) = ITAMPA
322           ICTRAN(2) = ITAMPB
323           ICTRAN(3) = ITAMPC
324           ICTRAN(4) = 0
325
326           NCTRAN = 1
327           IOPT   = 1
328           FILCMA = 'CCCMAT'
329
330           CALL CC_CMAT(ICTRAN,NCTRAN,LISTA,LISTB,LISTC,IOPT,
331     &                  FILCMA,IDUM,RDUM,0,WORK(KRHO1D),LWORK-KRHO1D)
332
333
334           IF (LCOR .OR. LSEC) THEN
335              CALL CC_CORE(WORK(KRHO1D),WORK(KRHO2D),ISYMTR)
336           ENDIF
337C
338           IF (IPRINT.GT.2) THEN
339             RHO1N = DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1)
340             RHO2N = DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1)
341             WRITE(LUPRI,*) 'Norm of RHO1: ',RHO1N,'aibj=',I
342             WRITE(LUPRI,*) 'Norm of RHO2: ',RHO2N,'aibj=',I
343           ENDIF
344           IF (IPRINT.GT.125) THEN
345            CALL CC_PRP(WORK(KRHO1D),WORK(KRHO2D),1,1,1)
346           ENDIF
347C
348           CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1)
349           CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1)
350           CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1)
351           CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1)
352           CALL DCOPY(NT1AMX,WORK(KRHO1D),1,
353     *              WORK(KF2+NTAMP*(I-1)),1)
354           CALL DCOPY(NT2AMX,WORK(KRHO2D),1,
355     *              WORK(KF2+NTAMP*(I-1)+NT1AMX),1)
356C
357           X12 = X12 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1)
358           X22 = X22 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1)
359           TI   = SECOND() - TI
360           IF (IPRINT.GT.5 ) THEN
361              WRITE(LUPRI,*) '  '
362              WRITE(LUPRI,*) 'FDB ROW NR. ',I+NT1AMX,
363     *                  ' DONE IN ',TI,' SEC.'
364           ENDIF
365C
366         ENDIF
367C
368 300    CONTINUE
369 200  CONTINUE
370C
371      WRITE(LUPRI,*) '    '
372      WRITE(LUPRI,*) '**  FINITE DIFF WITH DELTA ',DELTA, '**'
373      WRITE(LUPRI,*) '    '
374      IF ((IPRINT .GT. 4).AND.(.TRUE.)) THEN
375       CALL AROUND('FINITE DIFF. CC D*TA*TB*TC-Matrix - 11 & 21 PART' )
376       CALL OUTPUT(WORK(KF),1,NTAMP,1,NC1VEC,NTAMP,NC1VEC,1,LUPRI)
377       CALL AROUND('FINITE DIFF. CC D*TA*TB*TC-Matrix - 12 & 22 PART' )
378       CALL OUTPUT(WORK(KF+NTAMP*NC1VEC),1,NTAMP,1,NC2VEC,
379     *               NTAMP,NC2VEC,1,LUPRI)
380      ENDIF
381      IF (.TRUE.) THEN
382       XNJ = X11 + X12 + X21 + X22
383       WRITE(LUPRI,*)'  '
384       WRITE(LUPRI,*)' NORM OF FIN. DIFF. D*TA*TB*TC-Matrix.', SQRT(XNJ)
385       WRITE(LUPRI,*)'  '
386       WRITE(LUPRI,*)' NORM OF 11 PART OF FD. D*TA*TB*TC-mat.: ',
387     &      SQRT(X11)
388       WRITE(LUPRI,*)' NORM OF 21 PART OF FD. D*TA*TB*TC-mat.: ',
389     &      SQRT(X21)
390       WRITE(LUPRI,*)' NORM OF 12 PART OF FD. D*TA*TB*TC-mat.: ',
391     &      SQRT(X12)
392       WRITE(LUPRI,*)' NORM OF 22 PART OF FD. D*TA*TB*TC-mat.: ',
393     &      SQRT(X22)
394      ENDIF
395C
396C--------------------------------------
397C     Calculate Matrix times Ty vector.
398C--------------------------------------
399C
400      CALL DGEMV('N',NTAMP,NTAMP,ONE,WORK(KF),NTAMP,TYAM,1,
401     *           ZERO,RESULT,1)
402
403C--------------------------------------
404C     scale diagonal with 1/2:
405C--------------------------------------
406C     CALL CCLR_DIASCL(RESULT(NT1AM(1)+1),TWO,1)
407
408      WRITE (LUPRI,*) 'NTAMP:',NTAMP
409      WRITE (LUPRI,*) 'NORM^2 OF TXAM VECTOR:',
410     *   DDOT(NT1AM(1)+NT2AM(1),TXAM,1,TXAM,1)
411      WRITE (LUPRI,*) 'NORM^2 OF TYAM VECTOR:',
412     *   DDOT(NT1AM(1)+NT2AM(1),TYAM,1,TYAM,1)
413      WRITE (LUPRI,*) 'NORM^2 OF RESULT VECTOR:',
414     *   DDOT(NTAMP,RESULT,1,RESULT,1)
415
416C
417C-------------------------------------------------
418C     Restore the CC reference amplitudes on disk.
419C-------------------------------------------------
420C
421      LUTAM = -1
422      CALL GPOPEN(LUTAM,'TAM_SAV','UNKNOWN',' ','UNFORMATTED',IDUMMY,
423     *            .FALSE.)
424      REWIND(LUTAM)
425      READ(LUTAM) (WORK(KC1AM + I -1 ) , I = 1, NT1AMX)
426      READ(LUTAM) (WORK(KC2AM + I -1 ) , I = 1, NT2AMX)
427      CALL GPCLOSE(LUTAM,'DELETE')
428C
429      IOPT = 3
430      CALL CC_WRRSP('R0',0,1,IOPT,MODEL,DUMMY,WORK(KC1AM),
431     *              WORK(KC2AM),WORK(KEND2),LWRK2)
432C
433      RSPIM = .TRUE.
434      CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM),
435     &            WORK(KC2AM),
436     *            WORK(KEND2),LWRK2,'XXX')
437C
438      IF (IPRINT .GT. 10) THEN
439         CALL AROUND(' END OF CC_FDD ')
440      ENDIF
441C
442      RETURN
443      END
444*=====================================================================*
445