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_AAMAT */
20*=====================================================================*
21      SUBROUTINE CC_AAMAT(IATRAN,NATRAN,LISTO,LISTR,IOPTRES,FILAAM,
22     &                    IADOTS,ACONS,MXVEC,WORK,LWORK)
23*---------------------------------------------------------------------*
24*
25*    Purpose: driver for a list of A{O} matrix transformations
26*
27*           IOPTRES = 3 :  each result vector is written to its own
28*                          file by a call to CC_WRRSP, FILAAM is used
29*                          as list type and IATRAN(3,*) as list index
30*
31*           IOPTRES = 4 :  each result vector is added to a vector on
32*                          file by a call to CC_WARSP, FILAAM is used
33*                          as list type and IATRAN(3,*) as list index
34*
35*           IOPTRES = 5 : the result vectors are dotted on an array
36*                         of vectors, the type of the arrays is given
37*                         by FILAAM and the indices are taken from
38*                         IADOTS. the results are returned in ACONS
39*
40*    N.B.: this is a very first quick hack solution...
41*          should finally be incorporated in the CCCR_AA routine...
42*
43*    Written by Christof Haettig, maj 1996.
44*
45*=====================================================================*
46      IMPLICIT NONE
47#include "priunit.h"
48#include "ccorb.h"
49#include "ccsdinp.h"
50#include "ccsdsym.h"
51#include "ccroper.h"
52#include "cclists.h"
53#include "ccnoddy.h"
54#include "dummy.h"
55
56      CHARACTER*1 CDUMMY
57      PARAMETER (CDUMMY = ' ')
58
59* local parameters:
60      CHARACTER*(19) MSGDBG
61      PARAMETER (MSGDBG = '[debug] CC_AAMAT> ')
62      LOGICAL LOCDBG
63      PARAMETER (LOCDBG = .FALSE.)
64
65      INTEGER LWORK, NATRAN, MXVEC
66
67#if defined (SYS_CRAY)
68      REAL WORK(LWORK), ACONS(MXVEC,NATRAN)
69      REAL ZERO, TWO, FREQ
70#else
71      DOUBLE PRECISION WORK(LWORK), ACONS(MXVEC,NATRAN)
72      DOUBLE PRECISION ZERO, TWO, FREQ
73#endif
74      PARAMETER (ZERO = 0.0d0, TWO = 2.0D0)
75
76      CHARACTER LISTO*(*), LISTR*(*), FILAAM*(*)
77      CHARACTER*(10) MODEL
78      CHARACTER*(8)  LABEL
79      CHARACTER APROXR12*3
80      LOGICAL LTWOEL
81      INTEGER IATRAN(MXDIM_AATRAN,NATRAN), IADOTS(MXVEC,NATRAN)
82      INTEGER IOPER, IDLSTR, IFILE, IVEC, ITRAN, IOPT, IOPTRES
83      INTEGER ISYOPE, ISYTAM, ISYRES, IRELAX, IOPTE
84      INTEGER KEND1, LWRK1, KRES1, KRES2, LEN, KRES1EFF, KRES2EFF
85      INTEGER LENMOD,IOPTWR12,KRESR12
86
87
88* external functions:
89      INTEGER ILSTSYM
90      REAL*8  DDOT, FREQLST
91
92*---------------------------------------------------------------------*
93* initializations:
94*---------------------------------------------------------------------*
95      CALL QENTER('CC_AAMAT')
96      IF (LOCDBG) THEN
97        WRITE(LUPRI,*) 'Entered  cc_aamat...'
98        CALL FLSHFO(LUPRI)
99      END IF
100
101      IF ( .NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN
102        WRITE (LUPRI,'(/1x,a)') 'CC_AAMAT Called for a Coupled Cluster '
103     &          //'method not implemented in CC_AAMAT...'
104        CALL QUIT('Unknown CC method in CC_AAMAT.')
105      END IF
106
107* check list combinations:
108      IF ( (LISTO(1:2).NE.'o1' .AND. LISTO(1:2).NE.'o2')
109     &    .OR. ( LISTR(1:1).NE.'R'.AND. LISTR(1:2).NE.'ER')) THEN
110        WRITE (LUPRI,*)
111     &        'A{O} matrix driver called with a strange list set:'
112        WRITE (LUPRI,*) 'LISTO:',LISTO
113        WRITE (LUPRI,*) 'LISTR:',LISTR
114        CALL QUIT('A{O} mat. driver called with a strange '//
115     &             'set of list types.')
116      END IF
117
118      IF ( LISTO(1:2).EQ.'O2' .AND. NATRAN.NE.0) THEN
119         WRITE (LUPRI,*)
120     &        'LISTO="o2" case not yet implemented in CC_AAMAT.'
121         CALL QUIT('LISTO="o2" case not yet implemented in CC_AAMAT.')
122      END IF
123
124*---------------------------------------------------------------------*
125* start loop over all requested A{O} matrix transformations
126*---------------------------------------------------------------------*
127      DO ITRAN = 1, NATRAN
128
129        IOPER  = IATRAN(1,ITRAN)
130        IDLSTR = IATRAN(2,ITRAN)
131        IFILE  = IATRAN(3,ITRAN)
132        IRELAX = IATRAN(4,ITRAN)
133
134        LABEL  = LBLOPR(IOPER)
135        LTWOEL = LPDBSOP(IOPER)
136
137        ISYOPE = ILSTSYM(LISTO,IOPER)
138        ISYTAM = ILSTSYM(LISTR,IDLSTR)
139        ISYRES = MULD2H(ISYOPE,ISYTAM)
140
141        IF ((IOPTRES.EQ.3.OR.IOPTRES.EQ.4) .AND.
142     &       ILSTSYM(FILAAM,IFILE).NE.ISYRES) THEN
143          CALL QUIT('Symmetry mismatch in CC_AAMAT.')
144        END IF
145
146*---------------------------------------------------------------------*
147* calculate A{O} matrix times a response amplitude vector:
148*---------------------------------------------------------------------*
149
150        ! allocate memory for results vector
151        KRES1 = 1
152        KRES2 = KRES1 + NT1AM(ISYRES)
153        KEND1 = KRES2 + NT2AM(ISYRES)
154        IF (CCR12) THEN
155          KRESR12 = KRES2 + NT2AM(ISYRES)
156          KEND1   = KRESR12 + NTR12AM(ISYRES)
157        END IF
158        LWRK1 = LWORK - KEND1
159        IF (LWRK1.LT.0) CALL QUIT('Insufficient memory in CC_AAMAT.')
160
161        IF ( (IRELAX.GT.0) .OR. LTWOEL ) THEN
162
163           IF (CCSDT) CALL QUIT('relaxed with triples not available!')
164
165           CALL CC_FDAAMAT(LISTR,IDLSTR,WORK(KRES1),WORK(KRES2),
166     &                     LABEL,IRELAX,WORK(KEND1),LWRK1)
167
168        ELSE
169
170c          --------------------------------------------------
171c          do the CCS/CC2/CCSD parts:
172c          --------------------------------------------------
173           CALL CCCR_AA(LABEL,ISYOPE,LISTR,IDLSTR,DUMMY,WORK,LWORK)
174           CALL CCLR_DIASCL(WORK(KRES2),2.0d0,ISYRES)
175
176c          --------------------------------------------------
177c          do the triples parts:
178c          --------------------------------------------------
179           IF (CCSDT) THEN
180             ! find frequency associates with result vector
181             IF (IOPTRES.EQ.5) THEN
182               FREQ  = 0.0D0
183             ELSE
184               FREQ  = FREQLST(FILAAM,IFILE)
185             END IF
186
187             ! allocate memory for effective results vector
188             KRES1EFF = KEND1
189             KRES2EFF = KRES1EFF + NT1AM(ISYRES)
190             KEND1    = KRES2EFF + NT2AM(ISYRES)
191             LWRK1    = LWORK - KEND1
192             IF (LWRK1.LT.0)
193     &          CALL QUIT('Insufficient memory in CC_AAMAT. (2)')
194
195             CALL DZERO(WORK(KRES1EFF),NT1AM(ISYRES))
196             CALL DZERO(WORK(KRES2EFF),NT2AM(ISYRES))
197
198             IF (NODDY_AAMAT) THEN
199               CALL CCSDT_AAMAT_NODDY(IOPTRES,FREQ,LABEL,ISYOPE,
200     &                                LISTR,IDLSTR,.TRUE.,
201     &                                WORK(KRES1),WORK(KRES2),
202     &                                WORK(KRES1EFF),WORK(KRES2EFF),
203     &                                IADOTS,ACONS,FILAAM,ITRAN,
204     &                                NATRAN,MXVEC,WORK(KEND1),LWRK1)
205
206             ELSE
207
208               WRITE(LUPRI,*)'AAMATSD should be called from '
209               WRITE(LUPRI,*)'CC_BMAT module. Check NEW_RHS flag '
210               WRITE(LUPRI,*)'in CCRHSVEC.'
211               CALL QUIT('Real code for AAMATSD not called here...')
212
213             END IF
214
215
216           END IF
217
218        END IF
219
220*---------------------------------------------------------------------*
221* write/add to vector on file:
222*---------------------------------------------------------------------*
223        IF (CCS) THEN
224           MODEL = 'CCS       '
225           IOPT  = 1
226        ELSE IF (CC2) THEN
227           MODEL = 'CC2       '
228           IOPT  = 3
229        ELSE IF (CCSD) THEN
230           MODEL = 'CCSD      '
231           IOPT  = 3
232        ELSE IF (CC3) THEN
233           MODEL = 'CC3       '
234           IOPT  = 3
235           IOPTE = 24
236        ELSE
237           CALL QUIT('Unknown coupled cluster model in CC_AAMAT.')
238        END IF
239        IF (CCR12) THEN
240          APROXR12 = '   '
241          CALL CCSD_MODEL(MODEL,LENMOD,10,MODEL,10,APROXR12)
242          IOPTWR12 = 32
243        END IF
244
245        IF      (IOPTRES.EQ.3) THEN
246          CALL CC_WRRSP(FILAAM, IFILE, ISYRES, IOPT, MODEL,
247     &                  DUMMY,WORK(KRES1),WORK(KRES2),
248     &                  WORK(KEND1),LWRK1)
249          IF (CCR12) THEN
250            CALL CC_WRRSP(FILAAM,IFILE,ISYRES,IOPTWR12,
251     &                    MODEL,DUMMY,DUMMY,WORK(KRESR12),
252     &                    WORK(KEND1),LWRK1)
253          END IF
254          IF (CCSDT) THEN
255            CALL CC_WRRSP(FILAAM,IFILE,ISYRES,IOPTE,MODEL,DUMMY,
256     &                 WORK(KRES1EFF),WORK(KRES2EFF),WORK(KEND1),LWRK1)
257          END IF
258        ELSE IF (IOPTRES.EQ.4) THEN
259          CALL CC_WARSP(FILAAM, IFILE, ISYRES, IOPT, MODEL,
260     &                  DUMMY,WORK(KRES1),WORK(KRES2),
261     &                  WORK(KEND1),LWRK1)
262          IF (CCR12) THEN
263            CALL CC_WARSP(FILAAM,IFILE,ISYRES,IOPTWR12,
264     &                    MODEL,DUMMY,DUMMY,WORK(KRESR12),
265     &                    WORK(KEND1),LWRK1)
266          END IF
267          IF (CCSDT) THEN
268            CALL CC_WARSP(FILAAM,IFILE,ISYRES,IOPTE,MODEL,DUMMY,
269     &                 WORK(KRES1EFF),WORK(KRES2EFF),WORK(KEND1),LWRK1)
270          END IF
271        ELSE IF (IOPTRES.EQ.5) THEN
272           IF (LOCDBG) THEN
273             IVEC = 1
274             WRITE(LUPRI,*) 'ACONS TRIPLES CONTRIBUTION:'
275             DO WHILE (IADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC)
276                WRITE (LUPRI,*)
277     &                'ACONS:',IVEC,ITRAN,ACONS(IVEC,ITRAN),IOPT
278                IVEC = IVEC + 1
279             END DO
280           END IF
281
282           CALL CCDOTRSP(IADOTS,ACONS,IOPT,FILAAM,ITRAN,NATRAN,MXVEC,
283     &                   WORK(KRES1),WORK(KRES2),ISYRES,
284     &                   WORK(KEND1),LWRK1)
285           IF (CCR12) THEN
286             CALL CCDOTRSP(IADOTS,ACONS,IOPTWR12,FILAAM,ITRAN,NATRAN,
287     &                     MXVEC,DUMMY,WORK(KRESR12),ISYRES,
288     &                     WORK(KEND1),LWRK1)
289           END IF
290
291           IF (LOCDBG) THEN
292             IVEC = 1
293             DO WHILE (IADOTS(IVEC,ITRAN).NE.0 .AND. IVEC.LE.MXVEC)
294                WRITE (LUPRI,*)
295     &                'ACONS:',IVEC,ITRAN,ACONS(IVEC,ITRAN),IOPT
296                IVEC = IVEC + 1
297             END DO
298           END IF
299        ELSE
300          WRITE (LUPRI,*) 'Illegal option IOPTRES in CC_AAMAT.'
301          CALL QUIT('Illegal option IOPTRES in CC_AAMAT.')
302        END IF
303
304        IF (LOCDBG .AND. (IOPTRES.EQ.3 .OR. IOPTRES.EQ.4)) THEN
305         WRITE (LUPRI,*) MSGDBG, 'wrote ',FILAAM,':',IFILE,' to disk.'
306         LEN = NT1AM(ISYRES)
307         IF (.NOT.CCS) LEN = LEN + NT2AM(ISYRES)
308         IF (CCR12) LEN = LEN + NTR12AM(ISYRES)
309         WRITE (LUPRI,*) MSGDBG, 'NORM^2 for result vector = ',
310     &    DDOT(LEN,WORK(KRES1),1,WORK(KRES1),1)
311         WRITE (LUPRI,*) 'Listing of the result vector:'
312         CALL CC_PRP(WORK(KRES1),WORK(KRES2),ISYRES,1,1)
313         IF (CCR12) CALL CC_PRPR12(WORK(KRESR12),ISYRES,1,.TRUE.)
314         IF (CCSDT) THEN
315           WRITE (LUPRI,*) 'norm^2 of eff. result vector:',
316     &        DDOT(LEN,WORK(KRES1EFF),1,WORK(KRES1EFF),1)
317           WRITE (LUPRI,*) 'Listing of the eff. result vector:'
318           CALL CC_PRP(WORK(KRES1EFF),WORK(KRES2EFF),ISYRES,1,1)
319         END IF
320        END IF
321
322*---------------------------------------------------------------------*
323* end of loop over A{O} transformations:
324*---------------------------------------------------------------------*
325      END DO
326
327      CALL QEXIT('CC_AAMAT')
328      RETURN
329      END
330*---------------------------------------------------------------------*
331*               END OF SUBROUTINE CC_AAMAT                            *
332*---------------------------------------------------------------------*
333