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 CCREO2CON(IREO,NREO,IRDOTS,RCONS,
21     &                     MXVEC,WORK,LWORK)
22*---------------------------------------------------------------------*
23*
24*     Purpose: calculate contributions from second-order orbital
25*              reorthogonalization and the coupling of reorthog.
26*              and relaxation to response functions involving
27*              perturbation-dependent basis sets.
28*
29*              IREO    --  array with the perturbation indeces
30*              NREO    --  length of IREO
31*              IRDOTS  --  matrix containing the second indeces
32*              RCONS   --  matrix with the contributions
33*              MXVEC   --  leading dimension of IRDOTS and RCONS
34*
35*     Christof Haettig 11-6-1999
36*
37*---------------------------------------------------------------------*
38      IMPLICIT NONE
39#include "priunit.h"
40#include "dummy.h"
41#include "ccorb.h"
42#include "ccsdsym.h"
43#include "ccexpfck.h"
44#include "cc1dxfck.h"
45#include "ccr1rsp.h"
46#include "ccfro.h"
47#include "ccroper.h"
48#include "iratdef.h"
49#include "inftap.h"
50
51      LOGICAL LOCDBG
52      PARAMETER (LOCDBG = .FALSE.)
53
54      INTEGER ISYM0, LUFCK
55      PARAMETER( ISYM0 = 1)
56      CHARACTER LABEL0*(8)
57      PARAMETER( LABEL0 = 'HAM0    ' )
58
59      LOGICAL LORX, LFOCK0
60      INTEGER LWORK, NREO, MXVEC
61
62      INTEGER IREO(NREO), IRDOTS(MXVEC,NREO)
63
64#if defined (SYS_CRAY)
65      REAL TWO, FREQ, RCONS(MXVEC,NREO), WORK(LWORK)
66#else
67      DOUBLE PRECISION TWO, FREQ, RCONS(MXVEC,NREO), WORK(LWORK)
68#endif
69      PARAMETER (TWO = 2.0D0)
70
71      CHARACTER*(10) MODEL
72      CHARACTER*(8)  LABELA, LABELB
73      LOGICAL LORXA, LORXB, LPDBSA, LPDBSB, NOKAPPA
74      INTEGER NCMOT, NASHT, N2ASHX, LCINDX, IDXREO, IKAPPA
75      INTEGER KCMO, KFOCK0, KOVERLP, KEND1, LWRK1, IFOCK, IADRF, ISYM
76      INTEGER ISYMA, ISYMB, IOPERA, IOPERB, KR2EFF, KRMATA, KAPASQ
77      INTEGER KAPPAA, KAPBSQ, KAPPAB, KSCR1, KEND2, LWRK2, KRMATB
78      INTEGER KQMATPA, KQMATPB
79      INTEGER IOPT, IVEC, IKAPPB, ISAMA, ISAMB
80
81* external functions:
82      INTEGER IEFFFOCK
83      INTEGER ILSTSYM, ILSTSYMRLX
84      INTEGER IROPER
85#if defined (SYS_CRAY)
86      REAL DDOT
87#else
88      DOUBLE PRECISION DDOT
89#endif
90
91*---------------------------------------------------------------------*
92*     check, if there is anything at all to do:
93*---------------------------------------------------------------------*
94
95      IF (NREO.LE.0) RETURN
96
97*---------------------------------------------------------------------*
98*     get some constants from sirius common block:
99*---------------------------------------------------------------------*
100
101      CALL CC_SIRINF(NCMOT,NASHT,N2ASHX,LCINDX)
102
103*---------------------------------------------------------------------*
104*     allocate memory for perturbation-independent stuff needed:
105*---------------------------------------------------------------------*
106      KCMO    = 1
107      KFOCK0  = KCMO    + NCMOT
108      KOVERLP = KFOCK0  + N2BST(ISYM0)
109      KEND1   = KOVERLP + N2BST(ISYM0)
110      LWRK1   = LWORK   - KEND1
111
112      IF (LWRK1.LT.0) THEN
113         CALL QUIT('Insufficient work space in CCREO2CON.')
114      END IF
115
116*---------------------------------------------------------------------*
117*     read MO coefficients from file:
118*---------------------------------------------------------------------*
119
120      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
121     &            .FALSE.)
122      REWIND LUSIFC
123      CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI)
124      READ (LUSIFC)
125      READ (LUSIFC)
126      CALL READI(LUSIFC,IRAT*NCMOT,WORK(KCMO))
127      CALL GPCLOSE(LUSIFC,'KEEP')
128
129*---------------------------------------------------------------------*
130*     read overlap matrix from file:
131*---------------------------------------------------------------------*
132
133      IF (LWRK1.LT.NBAST) THEN
134         CALL QUIT('Insufficient work space in CCREO2CON.')
135      END IF
136
137      CALL RDONEL('OVERLAP ',.TRUE.,WORK(KEND1),NBAST)
138      CALL CCSD_SYMSQ(WORK(KEND1),ISYM0,WORK(KOVERLP))
139
140*---------------------------------------------------------------------*
141*     read zeroth-order effective Fock matrix if available:
142*---------------------------------------------------------------------*
143      IFOCK = IEFFFOCK(LABEL0,ISYM,1)
144
145      IF (LEXPFCK(2,IFOCK)) THEN
146         IADRF = IADRFCK(1,IFOCK)
147
148         LUFCK = -1
149         CALL WOPEN2(LUFCK,FILFCKEFF,64,0)
150         CALL GETWA2(LUFCK,FILFCKEFF,WORK(KFOCK0),IADRF,N2BST(ISYM0))
151         CALL WCLOSE2(LUFCK,FILFCKEFF,'KEEP')
152
153         CALL CC_EFFCKMO(WORK(KFOCK0),ISYM0,WORK(KCMO),WORK(KOVERLP),
154     &                   WORK(KEND1),LWRK1)
155
156         LFOCK0 = .TRUE.
157      ELSE
158         LFOCK0 = .FALSE.
159      END IF
160
161*---------------------------------------------------------------------*
162*     start loop over first index:
163*---------------------------------------------------------------------*
164      DO IDXREO = 1, NREO
165
166         IKAPPA = IREO(IDXREO)
167         ISYMA  = ILSTSYMRLX('R1',IKAPPA)
168         LABELA = LRTHFLBL(IKAPPA)
169         IOPERA = IROPER(LABELA,ISYMA)
170         ISAMA  = ISYMAT(IOPERA)
171         LORXA  = .TRUE.
172         LPDBSA = LPDBSOP(IOPERA)
173
174         IF (LOCDBG) THEN
175            WRITE (LUPRI,*) 'CCREO2CON> IDXREO:',IDXREO
176            WRITE (LUPRI,*) 'CCREO2CON> IKAPPA:',IKAPPA
177            WRITE (LUPRI,*) 'CCREO2CON> LABELA:',LABELA
178            WRITE (LUPRI,*) 'CCREO2CON> LORXA :',LORXA
179            WRITE (LUPRI,*) 'CCREO2CON> LPDBSA:',LPDBSA
180         END IF
181
182         KR2EFF  = KEND1
183         KQMATPA = KR2EFF  + N2BST(ISYM0)
184         KRMATA  = KQMATPA + N2BST(ISYMA)
185         KAPASQ  = KRMATA  + N2BST(ISYMA)
186         KAPPAA  = KAPASQ  + N2BST(ISYMA)
187         KQMATPB = KAPPAA  + 2*NALLAI(ISYMA)
188         KRMATB  = KQMATPB + N2BST(ISYMA)
189         KAPBSQ  = KRMATB  + N2BST(ISYMA)
190         KAPPAB  = KAPBSQ  + N2BST(ISYMA)
191         KSCR1   = KAPPAB  + NALLAI(ISYMA)
192         KEND2   = KSCR1   + N2BST(ISYMA)
193
194         LWRK2  = LWORK  - KEND2
195
196         IF (LWRK2 .LT. 0) THEN
197            CALL QUIT('Insufficient memory in CCREO2CON.')
198         END IF
199
200         IF (LORXA) THEN
201           CALL CC_RDHFRSP('R1 ',IKAPPA,ISYMA,WORK(KAPPAA))
202           CALL CCKAPPASQ(WORK(KAPASQ),WORK(KAPPAA),ISYMA,'N')
203         END IF
204
205         IF (LPDBSA) THEN
206           CALL CC_GET_RMAT(WORK(KSCR1),IOPERA,1,ISYMA,
207     &                      WORK(KEND2),LWRK2)
208
209           NOKAPPA = .TRUE.
210           CALL CC_QMAT(WORK(KQMATPA),WORK(KRMATA),WORK(KSCR1),DUMMY,
211     &                ISAMA,ISYMA,NOKAPPA,WORK(KCMO),WORK(KEND2),LWRK2)
212         END IF
213
214*---------------------------------------------------------------------*
215*        loop over second index:
216*---------------------------------------------------------------------*
217         IVEC = 1
218
219         DO WHILE (IRDOTS(IVEC,IDXREO).NE.0 .AND. IVEC.LE.MXVEC)
220
221C          WRITE (LUPRI,*) 'CCREO2CON> IVEC = ',IVEC
222
223           IKAPPB = IRDOTS(IVEC,IDXREO)
224           ISYMB  = ILSTSYMRLX('R1',IKAPPB)
225           LABELB = LRTHFLBL(IKAPPB)
226           IOPERB = IROPER(LABELB,ISYMB)
227           ISAMB  = ISYMAT(IOPERB)
228           LORXB  = .TRUE.
229           LPDBSB = LPDBSOP(IOPERB)
230
231           IF (LOCDBG) THEN
232              WRITE (LUPRI,*) 'CCREO2CON> IVEC  :',IVEC
233              WRITE (LUPRI,*) 'CCREO2CON> IKAPPB:',IKAPPB
234              WRITE (LUPRI,*) 'CCREO2CON> LABELB:',LABELB
235              WRITE (LUPRI,*) 'CCREO2CON> LORXB :',LORXB
236              WRITE (LUPRI,*) 'CCREO2CON> LPDBSB:',LPDBSB
237           END IF
238
239           IF (ISYMB.NE.ISYMA) THEN
240              WRITE (LUPRI,*) IDXREO, ISYMA, IVEC, ISYMB
241              CALL QUIT('symmetry mismatch in CCREO2CON.')
242           END IF
243
244           IF (LORXB) THEN
245             CALL CC_RDHFRSP('R1 ',IKAPPB,ISYMB,WORK(KAPPAB))
246             CALL CCKAPPASQ(WORK(KAPBSQ),WORK(KAPPAB),ISYMB,'N')
247           END IF
248
249           IF (LPDBSB) THEN
250             CALL CC_GET_RMAT(WORK(KSCR1),IOPERB,1,ISYMB,
251     &                        WORK(KEND2),LWRK2)
252
253             NOKAPPA = .TRUE.
254             CALL CC_QMAT(WORK(KQMATPB),WORK(KRMATB),WORK(KSCR1),DUMMY,
255     &                ISAMB,ISYMB,NOKAPPA,WORK(KCMO),WORK(KEND2),LWRK2)
256           END IF
257
258           CALL DZERO(WORK(KR2EFF),N2BST(ISYM0))
259
260           IF (LORXB.AND.LPDBSA) THEN
261              CALL CC_MMOMMO('N','N',+1.0D0,WORK(KAPBSQ),ISYMB,
262     &                       WORK(KRMATA),ISYMA,1.0D0,WORK(KR2EFF),1)
263              CALL CC_MMOMMO('N','N',-1.0D0,WORK(KRMATA),ISYMA,
264     &                       WORK(KAPBSQ),ISYMB,1.0D0,WORK(KR2EFF),1)
265           END IF
266
267           IF (LORXA.AND.LPDBSB) THEN
268              CALL CC_MMOMMO('N','N',+1.0D0,WORK(KAPASQ),ISYMA,
269     &                       WORK(KRMATB),ISYMB,1.0D0,WORK(KR2EFF),1)
270              CALL CC_MMOMMO('N','N',-1.0D0,WORK(KRMATB),ISYMB,
271     &                       WORK(KAPASQ),ISYMA,1.0D0,WORK(KR2EFF),1)
272           END IF
273
274           IF (LPDBSA.AND.LPDBSB) THEN
275              CALL QUIT('CCREO2CON NOT YET IMPLEMENTED '//
276     &             'FOR 2. DERIVATIVES.')
277           END IF
278
279           RCONS(IVEC,IDXREO) = TWO*DDOT(N2BST(ISYM0),WORK(KR2EFF),1,
280     &                                                WORK(KFOCK0),1)
281C          WRITE (LUPRI,*) 'CCREO2CON>',RCONS(IVEC,IDXREO)
282
283           IVEC = IVEC + 1
284         END DO
285
286      END DO
287
288*---------------------------------------------------------------------*
289*     print the results:
290*---------------------------------------------------------------------*
291      IF (LOCDBG) THEN
292         WRITE(LUPRI,*)
293     &        'CCREO2CON> results for X intermediate contribs.:'
294         IF (MXVEC.NE.0) THEN
295            DO IDXREO = 1, NREO
296               WRITE (LUPRI,*) 'IDXREO = ',IDXREO
297               IVEC = 1
298               DO WHILE(IRDOTS(IVEC,IDXREO).NE.0 .AND. IVEC.LE.MXVEC)
299                  WRITE(LUPRI,'(A,2I5,2X,E19.12)') 'CCREO2CON> ',
300     &              IREO(IDXREO),IRDOTS(IVEC,IDXREO),RCONS(IVEC,IDXREO)
301                  IVEC = IVEC + 1
302               END DO
303            END DO
304         ELSE
305            WRITE (LUPRI,*) 'MXVEC.EQ.0 --> nothing calculated.'
306         END IF
307      END IF
308
309      RETURN
310      END
311*======================================================================*
312