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 CCXIMCON(IXIM,NXIM,LISTQ,IXDOTS,XCONS,
21     &                    MXVEC,WORK,LWORK)
22*---------------------------------------------------------------------*
23*
24*     Purpose: calculate contributions from the generalized relaxation
25*              & reorthogonalization X intermediate to response
26*              functions involving orbital relaxation or perturbation-
27*              dependent basis sets.
28*
29*              IXIM    --  array with the perturbation indeces
30*              NXIM    --  length of IXIM
31*              LISTQ   --  type of vectors the X intermediates are to
32*                          be dotted on
33*              IXDOTS  --  matrix of indeces for the vectors with which
34*                          the dot products are to be calculated
35*              XCONS   --  matrix with the dot product results
36*              MXVEC   --  leading dimension of IXDOTS and XCONS
37*
38*     Christof Haettig 1-6-1999
39*
40*     N.B.: this routine is not yet adapted for QMATP diff. from QMATH
41*
42*---------------------------------------------------------------------*
43      IMPLICIT NONE
44#include "priunit.h"
45#include "dummy.h"
46#include "ccorb.h"
47#include "ccsdsym.h"
48#include "ccexpfck.h"
49#include "cc1dxfck.h"
50#include "ccr1rsp.h"
51#include "ccfro.h"
52#include "ccroper.h"
53#include "inftap.h"
54#include "iratdef.h"
55
56      LOGICAL LOCDBG
57      PARAMETER (LOCDBG = .FALSE.)
58
59      INTEGER ISYM0, LUFCK
60      PARAMETER( ISYM0 = 1)
61      CHARACTER LABEL0*(8)
62      PARAMETER( LABEL0 = 'HAM0    ' )
63
64      LOGICAL LORX, LFOCK0, LORXQ, LPDBS
65      CHARACTER*(3) LISTQ
66      INTEGER ISYXIM, LWORK, NXIM, MXVEC
67
68      INTEGER IXIM(NXIM), IXDOTS(MXVEC,NXIM)
69
70#if defined (SYS_CRAY)
71      REAL FREQ, XCONS(MXVEC,NXIM), WORK(LWORK)
72      REAL HALF, ONE, TWO, ZERO, FTRACE
73#else
74      DOUBLE PRECISION FREQ, XCONS(MXVEC,NXIM), WORK(LWORK)
75      DOUBLE PRECISION HALF, ONE, TWO, ZERO, FTRACE
76#endif
77      PARAMETER( ONE = 1.0D0, TWO = 2.0D0, ZERO = 0.0D0, HALF = 0.5D0 )
78
79      CHARACTER*(10) MODEL
80      CHARACTER*(8)  LABEL, LABELQ
81      LOGICAL NOKAPPA
82      INTEGER NCMOT, NASHT, N2ASHX, LCINDX
83      INTEGER IFOCK, IEXPV, IADRF, KOVERLP, KFOCK0, KFOCK1, KCMO, KAPPA
84      INTEGER IKAPPA, IOPT, IOPER, ISYM, IDXXIM, KOFF1, IFCK1
85      INTEGER KEND1,LWRK1,KEND2,LWRK2,KXIM,KSCR1,KRMAT,KQMATP,KQMATH
86      INTEGER IVEC, IFILE, ISYMQ, IOPERQ, KEND3, LWRK3
87      INTEGER ISYM1, ISYM2, KOFF2, KQTRP, IREAL
88
89* external functions:
90      INTEGER IEFFFOCK
91      INTEGER IEXPECT
92      INTEGER IROPER
93      INTEGER I1DXFCK
94      INTEGER ILSTSYMRLX
95#if defined (SYS_CRAY)
96      REAL DDOT
97#else
98      DOUBLE PRECISION DDOT
99#endif
100
101*---------------------------------------------------------------------*
102*     check, if there is anything at all to do:
103*---------------------------------------------------------------------*
104
105      IF (NXIM.LE.0) RETURN
106
107*---------------------------------------------------------------------*
108*     get some constants from sirius common block:
109*---------------------------------------------------------------------*
110
111      CALL CC_SIRINF(NCMOT,NASHT,N2ASHX,LCINDX)
112
113*---------------------------------------------------------------------*
114*     allocate memory for perturbation-independent stuff needed:
115*---------------------------------------------------------------------*
116      KCMO    = 1
117      KFOCK0  = KCMO    + NCMOT
118      KOVERLP = KFOCK0  + N2BST(ISYM0)
119      KEND1   = KOVERLP + N2BST(ISYM0)
120      LWRK1   = LWORK   - KEND1
121
122      IF (LWRK1.LT.0) THEN
123         CALL QUIT('Insufficient work space in CCXIMCON.')
124      END IF
125
126*---------------------------------------------------------------------*
127*     read MO coefficients from file:
128*---------------------------------------------------------------------*
129
130      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
131     &            .FALSE.)
132      REWIND LUSIFC
133      CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI)
134      READ (LUSIFC)
135      READ (LUSIFC)
136      CALL READI(LUSIFC,IRAT*NCMOT,WORK(KCMO))
137      CALL GPCLOSE(LUSIFC,'KEEP')
138
139*---------------------------------------------------------------------*
140*     read overlap matrix from file:
141*---------------------------------------------------------------------*
142
143      IF (LWRK1.LT.NBAST) THEN
144         CALL QUIT('Insufficient work space in CCXIMCON.')
145      END IF
146
147      CALL RDONEL('OVERLAP ',.TRUE.,WORK(KEND1),NBAST)
148      CALL CCSD_SYMSQ(WORK(KEND1),ISYM0,WORK(KOVERLP))
149
150*---------------------------------------------------------------------*
151*     read zeroth-order effective Fock matrix if available:
152*---------------------------------------------------------------------*
153      call quit('hjaaj: code error, ISYM not defined, please fix!')
154      IFOCK = IEFFFOCK(LABEL0,ISYM,1)
155
156      IF (LEXPFCK(2,IFOCK)) THEN
157         IADRF = IADRFCK(1,IFOCK)
158
159         LUFCK = -1
160         CALL WOPEN2(LUFCK,FILFCKEFF,64,0)
161         CALL GETWA2(LUFCK,FILFCKEFF,WORK(KFOCK0),IADRF,N2BST(ISYM0))
162         CALL WCLOSE2(LUFCK,FILFCKEFF,'KEEP')
163
164         CALL CC_EFFCKMO(WORK(KFOCK0),ISYM0,WORK(KCMO),WORK(KOVERLP),
165     &                   WORK(KEND1),LWRK1)
166
167         LFOCK0 = .TRUE.
168
169         IF (LOCDBG) THEN
170            WRITE (LUPRI,*)
171     &           'CCXIMCON> zeroth-order effective Fock in ao/AO:'
172            CALL CC_PRONELAO(WORK(KFOCK0),ISYM0)
173         END IF
174
175      ELSE
176         LFOCK0 = .FALSE.
177      END IF
178
179*---------------------------------------------------------------------*
180*     start loop over X intermediates:
181*---------------------------------------------------------------------*
182      DO IDXXIM = 1, NXIM
183
184         IKAPPA = IXIM(IDXXIM)
185         LABEL  = LRTHFLBL(IKAPPA)
186         IOPER  = IROPER(LABEL,ISYXIM)
187         LORX   = .TRUE.
188         LPDBS  = LPDBSOP(IOPER)
189         ISYXIM = ILSTSYMRLX('R1',IKAPPA)
190
191         IF (LOCDBG) THEN
192            WRITE (LUPRI,*) 'CCXIMCON> IDXXIM:',IDXXIM
193            WRITE (LUPRI,*) 'CCXIMCON> IKAPPA:',IKAPPA
194            WRITE (LUPRI,*) 'CCXIMCON> LABEL :',LABEL
195            WRITE (LUPRI,*) 'CCXIMCON> LORX  :',LORX
196         END IF
197
198         KXIM   = KEND1
199         KRMAT  = KXIM   + N2BST(ISYXIM)
200         KQMATP = KRMAT  + N2BST(ISYXIM)
201         KQMATH = KQMATP + N2BST(ISYXIM)
202         KQTRP  = KQMATH + N2BST(ISYXIM)
203         KAPPA  = KQTRP  + N2BST(ISYXIM)
204         KEND2  = KAPPA  + 2*NALLAI(ISYXIM)
205
206         LWRK2  = LWORK  - KEND2
207
208         IF (LWRK2 .LT. 0) THEN
209            CALL QUIT('Insufficient memory in CCXIMCON.')
210         END IF
211
212         ! read first-order effective Fock matrix from file
213         call quit('hjaaj: code error, ISYM not defined, please fix!')
214         IFOCK = IEFFFOCK(LABEL,ISYM,1)
215         IEXPV = IEXPECT(LABEL,ISYM)
216         IADRF = IADRFCK(1,IFOCK)
217
218         CALL WOPEN2(LUFCK,FILFCKEFF,64,0)
219         CALL GETWA2(LUFCK,FILFCKEFF,WORK(KXIM),IADRF,N2BST(ISYXIM))
220         CALL WCLOSE2(LUFCK,FILFCKEFF,'KEEP')
221
222         IF (LOCDBG) THEN
223            FTRACE = ZERO
224            IF (ISYXIM.EQ.1) THEN
225               DO ISYM = 1, NSYM
226                  KOFF1 = KXIM + IAODIS(ISYM,ISYM)
227                  DO I = 1, NBAS(ISYM)
228                    FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1)
229                  END DO
230               END DO
231            END IF
232            WRITE (LUPRI,*) 'LABEL:',LABEL
233            WRITE (LUPRI,*) 'ISYXIM,IFOCK,IEXPV:',ISYXIM,IFOCK,IEXPV
234            WRITE (LUPRI,*) 'FTRACE of read matrix:',FTRACE
235            WRITE (LUPRI,*) 'one-electron expect:',EXPVALUE(1,IEXPV)
236            WRITE (LUPRI,*) 'two-electron expect:',EXPVALUE(2,IEXPV)
237         END IF
238
239         ! transform effective Fock matrix to MO basis
240         CALL CC_EFFCKMO(WORK(KXIM),ISYXIM,WORK(KCMO),WORK(KOVERLP),
241     &                   WORK(KEND2),LWRK2)
242
243         IF (LOCDBG) THEN
244            FTRACE = ZERO
245            IF (ISYXIM.EQ.1) THEN
246               DO ISYM = 1, NSYM
247                  KOFF1 = KXIM + IAODIS(ISYM,ISYM)
248                  DO I = 1, NBAS(ISYM)
249                    FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1)
250                  END DO
251               END DO
252            END IF
253            WRITE (LUPRI,*) 'LABEL:',LABEL
254            WRITE (LUPRI,*) 'ISYXIM:',ISYXIM
255            WRITE (LUPRI,*)
256     &           'FTRACE of matrix generated in CC_EFCKMO:',FTRACE
257         END IF
258
259         ! scale with 2 to get contribution to X intermediate
260         CALL DSCAL(N2BST(ISYXIM),TWO,WORK(KXIM),1)
261
262         IF (LOCDBG) THEN
263            WRITE (LUPRI,*)
264     &           'CCXIMCON> direct contribution to X intermediate:'
265            CALL CC_PRONELAO(WORK(KXIM),ISYXIM)
266         END IF
267
268         IF (LORX .OR. LPDBS) THEN
269
270           KSCR1  = KEND2
271           KFOCK1 = KSCR1  + N2BST(ISYXIM)
272           KEND3  = KFOCK1 + N2BST(ISYXIM)
273           LWRK3  = LWORK  - KEND3
274
275           IF (LWRK3 .LT. 0) THEN
276              CALL QUIT('Insufficient memory in CCXIMCON.')
277           END IF
278
279           IF (LORX) THEN
280             CALL CC_RDHFRSP('R1 ',IKAPPA,ISYXIM,WORK(KAPPA))
281           ELSE
282             CALL DZERO(WORK(KAPPA),2*NALLAI(ISYXIM))
283           END IF
284
285           IOPER = IROPER(LABEL,ISYXIM)
286           CALL CC_GET_RMAT(WORK(KRMAT),IOPER,1,ISYXIM,
287     &                      WORK(KEND3),LWRK3)
288
289           NOKAPPA = .FALSE.
290           IREAL   = ISYMAT(IOPER)
291           CALL CC_QMAT(WORK(KQMATP),WORK(KQMATH),
292     &                  WORK(KRMAT),WORK(KAPPA),
293     &                  IREAL,ISYXIM,NOKAPPA,
294     &                  WORK(KCMO),WORK(KEND3),LWRK3)
295
296           DO ISYM1 = 1, NSYM
297              ISYM2 = MULD2H(ISYM1,ISYXIM)
298              KOFF1 = KQMATP + IAODIS(ISYM1,ISYM2)
299              KOFF2 = KQTRP  + IAODIS(ISYM2,ISYM1)
300              CALL TRSREC(NBAS(ISYM1),NBAS(ISYM2),
301     &                    WORK(KOFF1),WORK(KOFF2))
302           END DO
303           CALL DSCAL(N2BST(ISYXIM),-ONE,WORK(KQTRP ),1)
304
305           CALL CC_MMOMMO('N','N',ONE,WORK(KFOCK0),ISYM0,
306     &                    WORK(KQTRP ),ISYXIM,ZERO,WORK(KSCR1),ISYXIM)
307
308           IF (LOCDBG) THEN
309              WRITE (LUPRI,*)
310     &             'CCXIMCON> relax. contrib. 1 to X intermediate:'
311              CALL CC_PRONELAO(WORK(KSCR1),ISYXIM)
312              FTRACE = ZERO
313              DO ISYM = 1, NSYM
314                 KOFF1 = KSCR1 + IAODIS(ISYM,ISYM)
315                 DO I = 1, NBAS(ISYM)
316                   FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1)
317                 END DO
318              END DO
319              WRITE (LUPRI,*) 'trace:',FTRACE
320           END IF
321
322           ! add to result matrix:
323           CALL DAXPY(N2BST(ISYXIM),ONE,WORK(KSCR1),1,WORK(KXIM),1)
324
325           ! read contribution from 'one-index' transformed density
326           ! from file
327           FREQ  = FRQLRTHF(IKAPPA)
328           IFCK1 = I1DXFCK('HAM0    ','R1 ',LABEL,FREQ,ISYXIM)
329           IADRF = IADR1DXF(1,IFCK1)
330           CALL WOPEN2(LUFCK,FIL1DXFCK,64,0)
331           CALL GETWA2(LUFCK,FIL1DXFCK,WORK(KFOCK1),IADRF,N2BST(ISYXIM))
332           CALL WCLOSE2(LUFCK,FIL1DXFCK,'KEEP')
333
334           CALL CC_EFFCKMO(WORK(KFOCK1),ISYXIM,WORK(KCMO),WORK(KOVERLP),
335     &                     WORK(KEND3),LWRK3)
336
337           IF (LOCDBG) THEN
338              WRITE (LUPRI,*)
339     &             'CCXIMCON> relax. contrib. 2 to X intermediate:'
340              CALL CC_PRONELAO(WORK(KFOCK1),ISYXIM)
341              FTRACE = ZERO
342              DO ISYM = 1, NSYM
343                 KOFF1 = KFOCK1 + IAODIS(ISYM,ISYM)
344                 DO I = 1, NBAS(ISYM)
345                   FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I-1)
346                 END DO
347              END DO
348              WRITE (LUPRI,*) 'trace:',FTRACE
349           END IF
350
351           ! add to result matrix:
352           CALL DAXPY(N2BST(ISYXIM),ONE,WORK(KFOCK1),1,WORK(KXIM),1)
353
354         END IF
355
356         IF (LOCDBG) THEN
357            WRITE (LUPRI,*) 'CCXIMCON> final result for X intermediate:'
358            CALL CC_PRONELAO(WORK(KXIM),ISYXIM)
359            FTRACE = ZERO
360            DO ISYM = 1, NSYM
361               KOFF1 = KXIM-1 + IAODIS(ISYM,ISYM)
362               DO I = 1, NBAS(ISYM)
363                 FTRACE = FTRACE + WORK(KOFF1+NBAS(ISYM)*(I-1)+I)
364               END DO
365            END DO
366            WRITE (LUPRI,*) 'trace of X intermediate:',FTRACE
367         END IF
368
369*---------------------------------------------------------------------*
370*        calculate all required dot products with this X intermediate:
371*---------------------------------------------------------------------*
372         IVEC = 1
373
374         DO WHILE (IXDOTS(IVEC,IDXXIM).NE.0 .AND. IVEC.LE.MXVEC)
375
376
377           IFILE = IXDOTS(IVEC,IDXXIM)
378           ISYMQ = ILSTSYMRLX(LISTQ,IFILE)
379
380           IF (ISYMQ.NE.ISYXIM) THEN
381              WRITE (LUPRI,*) IDXXIM, ISYXIM, LISTQ, IFILE, ISYMQ
382              CALL QUIT('symmetry mismatch in CCXIMCON.')
383           END IF
384
385           IF (LISTQ(1:3).EQ.'R1 ') THEN
386              LORXQ  = .TRUE.
387              LABELQ = LRTHFLBL(IFILE)
388           ELSE
389              CALL QUIT('Unknown LISTQ in CCXIMCON.')
390           END IF
391
392           IOPERQ = IROPER(LABELQ,ISYMQ)
393           IREAL  = ISYMAT(IOPERQ)
394
395           IF (LORXQ) THEN
396             CALL CC_RDHFRSP(LISTQ,IFILE,ISYMQ,WORK(KAPPA))
397           ELSE
398             CALL DZERO(WORK(KAPPA),2*NALLAI(ISYMQ))
399           END IF
400
401           CALL CC_GET_RMAT(WORK(KRMAT),IOPERQ,1,ISYMQ,
402     &                      WORK(KEND2),LWRK2)
403
404           NOKAPPA = .FALSE.
405           CALL CC_QMAT(WORK(KQMATP),WORK(KQMATH),
406     &                  WORK(KRMAT),WORK(KAPPA),
407     &                  IREAL,ISYMQ,NOKAPPA,
408     &                  WORK(KCMO),WORK(KEND2),LWRK2)
409
410           DO ISYM1 = 1, NSYM
411              ISYM2 = MULD2H(ISYM1,ISYMQ)
412              KOFF1 = KQMATH + IAODIS(ISYM1,ISYM2)
413              KOFF2 = KQTRP  + IAODIS(ISYM2,ISYM1)
414              CALL TRSREC(NBAS(ISYM1),NBAS(ISYM2),
415     &                    WORK(KOFF1),WORK(KOFF2))
416           END DO
417           CALL DCOPY(N2BST(ISYMQ),WORK(KQTRP),1,WORK(KQMATH),1)
418           CALL DSCAL(N2BST(ISYMQ),HALF,WORK(KQMATH),1)
419
420           DO ISYM1 = 1, NSYM
421              ISYM2 = MULD2H(ISYM1,ISYMQ)
422              KOFF1 = KQMATP + IAODIS(ISYM1,ISYM2)
423              KOFF2 = KQTRP  + IAODIS(ISYM2,ISYM1)
424              CALL TRSREC(NBAS(ISYM1),NBAS(ISYM2),
425     &                    WORK(KOFF1),WORK(KOFF2))
426           END DO
427           CALL DCOPY(N2BST(ISYMQ),WORK(KQTRP),1,WORK(KQMATP),1)
428           CALL DSCAL(N2BST(ISYMQ),HALF,WORK(KQMATP),1)
429
430           XCONS(IVEC,IDXXIM) =
431     &                  DDOT(N2BST(ISYMQ),WORK(KQMATH),1,WORK(KXIM),1)+
432     &      DBLE(IREAL)*DDOT(N2BST(ISYMQ),WORK(KQMATP),1,WORK(KXIM),1)
433
434
435           IVEC = IVEC + 1
436         END DO
437
438      END DO
439
440*---------------------------------------------------------------------*
441*     print the results:
442*---------------------------------------------------------------------*
443      IF (LOCDBG) THEN
444         WRITE(LUPRI,*) 'CCXIMCON> results for '//
445     &        'X intermediate contribs.:'
446         IF (MXVEC.NE.0) THEN
447            DO IDXXIM = 1, NXIM
448               WRITE (LUPRI,*) 'IDXXIM = ',IDXXIM
449               IVEC = 1
450               DO WHILE(IXDOTS(IVEC,IDXXIM).NE.0 .AND. IVEC.LE.MXVEC)
451                  WRITE(LUPRI,'(A,2I5,2X,E19.12)') 'CCXIMCON> ',
452     &              IXIM(IDXXIM),IXDOTS(IVEC,IDXXIM),XCONS(IVEC,IDXXIM)
453                  IVEC = IVEC + 1
454               END DO
455            END DO
456         ELSE
457            WRITE (LUPRI,*) 'MXVEC.EQ.0 --> nothing calculated.'
458         END IF
459      END IF
460
461      RETURN
462      END
463*======================================================================*
464