1*=====================================================================*
2c /* deck ccxopa_eomsetup */
3*=====================================================================*
4      SUBROUTINE CCXOPA_EOMSETUP(IFTRAN, IFDOTS, FCONS,
5     &                       NFTRAN, MXFTRAN, MXFVEC,
6     &                       IEATRAN, IEADOTS, EACONS,
7     &                       NXE1TRAN,MXATRAN, MXAVEC,
8     &                       IXE2TRAN,IX2DOTS, X2CONS,
9     &                       NXE2TRAN,MXXTRAN, MXXVEC,
10     &                       IEOMXE2TRAN,IEOMX2DOTS, EOMX2CONS,
11     &                       NEOMXE2TRAN,MXEOMXTRAN, MXEOMXVEC,
12     &                       IEOML0TRAN,IEOML0DOTS,EOML0CONS,
13     &                       NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC,
14     &                       RESULT, MXOPA, LADD, WORK, LWORK )
15*---------------------------------------------------------------------*
16*
17*    Purpose: set up for CC first-order transition moments
18*         - list of B matrix transformations with eigenvectors
19*         - list of A{X} matrix transformations with eigenvectors
20*         - list of XKSI vector contractions with Nbar multipliers
21*
22*     Written by Christof Haettig, Oct 2003
23*
24*     For EOM XOPA
25*         - list of XKSI vector contractions with LE vectors
26*           and of dot products tbar*RE)
27*     Sonia Coriani, Nov 2015
28*
29*=====================================================================*
30      IMPLICIT NONE
31#include "priunit.h"
32#include "cclists.h"
33#include "ccxopainf.h"
34#include "ccroper.h"
35#include "ccexci.h"
36#include "ccsdinp.h"
37#include "ccorb.h"
38
39* local parameters:
40      CHARACTER*(22) MSGDBG
41      PARAMETER (MSGDBG = '[debug] CCXOPA_EOMSETUP> ')
42      LOGICAL LOCDBG
43      PARAMETER (LOCDBG = .false.)
44
45      LOGICAL LADD
46      INTEGER MXOPA,MXFTRAN,MXFVEC,MXATRAN,MXAVEC,MXXTRAN,MXXVEC
47      INTEGER MXEOMXTRAN,MXEOMXVEC
48      INTEGER MXEOML0TRAN,MXEOML0VEC
49
50      INTEGER IFTRAN(MXDIM_FTRAN,MXFTRAN)
51      INTEGER IFDOTS(MXFVEC,MXFTRAN)
52      INTEGER IEATRAN(MXDIM_XEVEC,MXATRAN)
53      INTEGER IEADOTS(MXAVEC,MXATRAN)
54      INTEGER IXE2TRAN(MXDIM_XEVEC,MXXTRAN)
55      INTEGER IX2DOTS(MXXVEC,MXXTRAN)
56!
57! EOM is one of these
58!
59      INTEGER IEOMXE2TRAN(MXDIM_XEVEC,MXEOMXTRAN)
60      INTEGER IEOMX2DOTS(MXEOMXVEC,MXEOMXTRAN)
61      INTEGER IEOML0TRAN(MXEOML0TRAN)
62      INTEGER IEOML0DOTS(MXEOML0VEC,MXEOML0VEC)
63
64      INTEGER NFTRAN, NXE1TRAN, NXE2TRAN, LWORK
65      INTEGER NEOMXE2TRAN, NEOML0TRAN
66
67#if defined (SYS_CRAY)
68      REAL RESULT(MXOPA)
69      REAL FCONS(MXFVEC,MXFTRAN)
70      REAL EACONS(MXAVEC,MXATRAN)
71      REAL X2CONS(MXXVEC,MXXTRAN)
72      REAL EOMX2CONS(MXEOMXVEC,MXEOMXTRAN)     !(csi*LE)
73      REAL EOML0CONS(MXEOML0VEC,MXEOML0VEC) !(Tbar*RE)
74      !REAL EOMCONS(MXXVEC,MXEOMXTRAN)     !(csi*LE).(Tbar*RE)
75      REAL WORK(LWORK)
76      REAL ZERO, SIGN, EIGVI, EIGVF
77      REAL WIAF, WXINIF, WIBF, WLiXi, WL0Rf, WLiXiRf
78#else
79      DOUBLE PRECISION RESULT(MXOPA)
80      DOUBLE PRECISION FCONS(MXFVEC,MXFTRAN)
81      DOUBLE PRECISION EACONS(MXAVEC,MXATRAN)
82      DOUBLE PRECISION X2CONS(MXXVEC,MXXTRAN)
83      DOUBLE PRECISION EOMX2CONS(MXEOMXVEC,MXEOMXTRAN)!(csi*LE)
84      DOUBLE PRECISION EOML0CONS(MXEOML0VEC,MXEOML0TRAN) !(Tbar*RE)
85      !DOUBLE PRECISION EOMCONS(MXEOMXVEC,MXEOMXTRAN)  !(csi*LE).(Tbar*RE)
86      DOUBLE PRECISION WORK(LWORK)
87      DOUBLE PRECISION ZERO, SIGN, EIGVI, EIGVF
88      DOUBLE PRECISION WIAF, WXINIF, WIBF, WLiXi, WL0Rf, WLiXiRf
89#endif
90      PARAMETER (ZERO = 0.0D0)
91
92      CHARACTER LABEL*(8)
93      LOGICAL LORX, LPDBS
94      INTEGER ITRAN, I, IRSD, IRSDX, ISTATEI, ISTATEF, ISYMI, ISYMF,
95     &        ISTISY, ISTFSY, IOP, IOPER, ISYMO, ISYME, ITURN,
96     &        IKAP, MXEAVEC, MXE2VEC, IN2VEC, IR1VEC, MFVEC,
97     &        ITMIF, IVEC, NBOPA, IDUM,
98     &        IMULI, IMULF
99
100      INTEGER MXEOME2VEC, MXEOM2
101      INTEGER NL0TRAN
102
103* external functions:
104      INTEGER IR1TAMP
105      INTEGER IN2AMP
106
107*---------------------------------------------------------------------*
108* initializations:
109* initialize for EOM as well ....
110*---------------------------------------------------------------------*
111      DO ITRAN = 1, MXATRAN
112       IEATRAN(1,ITRAN)  = 0
113       IEATRAN(2,ITRAN)  = 0
114       IEATRAN(3,ITRAN)  = -1
115       IEATRAN(4,ITRAN)  = -1
116       IEATRAN(5,ITRAN)  = 0
117       DO IVEC = 1, MXAVEC
118        IEADOTS(IVEC,ITRAN) = 0
119       END DO
120      END DO
121
122      DO ITRAN = 1, MXXTRAN
123       IXE2TRAN(1,ITRAN)  = 0
124       IXE2TRAN(2,ITRAN)  = 0
125       IXE2TRAN(3,ITRAN)  = -1
126       IXE2TRAN(4,ITRAN)  = -1
127       IXE2TRAN(5,ITRAN)  = 0
128       DO IVEC = 1, MXXVEC
129        IX2DOTS(IVEC,ITRAN) = 0
130       END DO
131      END DO
132
133      DO ITRAN = 1, MXEOMXTRAN
134       IEOMXE2TRAN(1,ITRAN)  = 0
135       IEOMXE2TRAN(2,ITRAN)  = 0
136       IEOMXE2TRAN(3,ITRAN)  = -1
137       IEOMXE2TRAN(4,ITRAN)  = -1
138       IEOMXE2TRAN(5,ITRAN)  = 0
139       DO IVEC = 1, MXEOMXVEC
140        IEOMX2DOTS(IVEC,ITRAN) = 0
141       END DO
142      END DO
143
144      !megaredundant but I am too tired to make it smarter...
145      DO ITRAN = 1, MXEOML0TRAN
146         IEOML0TRAN(ITRAN) = 0
147         DO IVEC = 1, MXEOML0VEC
148            IEOML0DOTS(IVEC,ITRAN) = 0
149         END DO
150      END DO
151
152      DO ITRAN = 1, MXFTRAN
153       DO I = 1, 3
154        IFTRAN(I,ITRAN)  = 0
155       END DO
156       DO IVEC = 1, MXFVEC
157        IFDOTS(IVEC,ITRAN)  = 0
158       END DO
159      END DO
160
161      NFTRAN   = 0
162      NXE1TRAN = 0
163      NXE2TRAN = 0
164      NEOMXE2TRAN = 0
165      NEOML0TRAN = 0
166
167      NBOPA   = 0
168      MFVEC   = 0
169      MXE2VEC = 0
170      MXEAVEC = 0
171
172      MXEOME2VEC = 0
173      MXEOM2 = 0
174      NL0TRAN = 0
175
176! mi manca qualcosa qua sopra...
177
178*---------------------------------------------------------------------*
179* start loop over all requested transition moments:
180*---------------------------------------------------------------------*
181      DO IRSDX  = 1, 2*NXQR2ST
182       ITURN = 1 + (IRSDX-1)/NXQR2ST
183       IRSD  = IRSDX - (ITURN-1)*NXQR2ST
184
185       IF (ITURN.EQ.1) THEN
186         ISTATEI = IQR2ST(IRSD,1)
187         ISTATEF = IQR2ST(IRSD,2)
188       ELSE IF (ITURN.EQ.2) THEN
189         ! switch state indices (and thereby also the sign of the freqs)
190         ! to get the conjugated transition moments
191         ISTATEI = IQR2ST(IRSD,2)
192         ISTATEF = IQR2ST(IRSD,1)
193       ELSE
194         CALL QUIT('Error in CCXOPA_EOMSETUP')
195       END IF
196
197       ISYMI   = ISYEXC(ISTATEI)
198       ISYMF   = ISYEXC(ISTATEF)
199       ISYME   = MULD2H(ISYMI,ISYMF)
200       ISTISY  = ISTATEI - ISYOFE(ISYMI)
201       ISTFSY  = ISTATEF - ISYOFE(ISYMF)
202       EIGVI   = EIGVAL(ISTATEI)
203       EIGVF   = EIGVAL(ISTATEF)
204CRF    We need to know multiplicities as well.
205       IMULI   = IMULTE(ISTATEI)
206       IMULF   = IMULTE(ISTATEF)
207
208       IF (LOCDBG) THEN
209         WRITE(LUPRI,*) 'CCXOPA_EOMSETUP:'
210         WRITE(LUPRI,*) 'ITURN,IRSD:',ITURN,IRSD
211         WRITE(LUPRI,*) 'ISTATEI,ISTATEF:',ISTATEI,ISTATEF
212         WRITE(LUPRI,*) 'ISYMI,ISYMF:',ISYMI,ISYMF
213         WRITE(LUPRI,*) 'ISTISY,ISTFSY:',ISTISY,ISTFSY
214         WRITE(LUPRI,*) 'EIGVI,EIGVF:',EIGVI,EIGVF
215         WRITE(LUPRI,*) 'IMULI,IMULF:',IMULI,IMULF
216       END IF
217
218       IF (IMULI.NE.IMULF) THEN
219          WRITE(LUPRI,*) 'Only singlet operators are currently'
220     &                   //' implemented!'
221          CYCLE
222       END IF
223
224       DO IOP = 1, NQR2OP
225        IOPER = IQR2OP(IOP)
226        LORX  = .FALSE.
227        ISYMO = ISYOPR(IOPER)
228        LABEL = LBLOPR(IOPER)
229        LPDBS = LPDBSOP(IOPER)
230        IKAP  = 0
231
232        IF (LPDBS) CALL QUIT('perturbation-dependent basis sets not '//
233     &              'implemented in CCXOPA_EOMSETUP.')
234
235        IF (ISYMO.EQ.ISYME) THEN
236
237          NBOPA = NBOPA + 1
238
239          IF (NBOPA.GT.MXOPA) THEN
240             CALL QUIT('NBOPA out of range in CCXOPA_EOMSETUP.')
241          END IF
242
243*---------------------------------------------------------------------*
244*         in all cases we need LE x A{X} x RE
245*---------------------------------------------------------------------*
246          !write(lupri,*) "Call CC_SETXE('Eta')"
247          !call flshfo(lupri)
248
249          CALL CC_SETXE('Eta',IEATRAN,IEADOTS,MXATRAN,MXAVEC,
250     &                  ISTATEI,IOPER,IKAP,0,0,0,ISTATEF,ITRAN,IVEC)
251          NXE1TRAN = MAX(NXE1TRAN,ITRAN)
252          MXEAVEC  = MAX(MXEAVEC, IVEC)
253          WIAF     = EACONS(IVEC,ITRAN)
254
255*---------------------------------------------------------------------*
256*         add N2 * Xksi{X} or LE * B * RE * R1, depending on QR22N1
257*---------------------------------------------------------------------*
258          WXINIF  = ZERO
259          WIBF    = ZERO
260          WLiXi   = ZERO
261          WL0Rf   = ZERO
262          WLiXiRf = ZERO
263
264          IF (.NOT.CIS) THEN
265            IF (QR22N1) THEN
266              IN2VEC=IN2AMP(ISTATEI,-EIGVI,ISYMI,ISTATEF,+EIGVF,ISYMF)
267              CALL CC_SETXE('Xi ',IXE2TRAN,IX2DOTS,MXXTRAN,MXXVEC,
268     &                      0,IOPER,IKAP,0,0,0,IN2VEC,ITRAN,IVEC)
269              NXE2TRAN = MAX(NXE2TRAN,ITRAN)
270              MXE2VEC  = MAX(MXE2VEC, IVEC)
271              WXINIF   = X2CONS(IVEC,ITRAN)
272            ELSE IF (LEOMXOPA.AND.(IMULF.EQ.1).AND.(ISYMF.EQ.1)) then
273              ! Only do this term if the F state is totally symmetric in
274              ! both spin and space!
275              !
276              ! recover index of LE and do dot prod with Xi(X)
277              !
278              CALL CC_SETXE('Xi ',IEOMXE2TRAN,IEOMX2DOTS,
279     &                            MXEOMXTRAN,MXEOMXVEC,
280     &             0,IOPER,IKAP,0,0,0,ISTATEI,ITRAN,IVEC)
281              NEOMXE2TRAN = MAX(NEOMXE2TRAN,ITRAN)
282              MXEOME2VEC  = MAX(MXEOME2VEC, IVEC)
283              WLiXi  = EOMX2CONS(IVEC,ITRAN)
284              !write(lupri,*) "WLiXi:", WLiXi
285              !
286              !now we fill in the list required for the dot products
287              !
288              CALL CC_SETDOT(IEOML0TRAN,IEOML0DOTS,
289     &                       MXEOML0TRAN,MXEOML0VEC,
290     &                       0,ISTATEF,ITRAN,IVEC)
291              NEOML0TRAN  = MAX(NEOML0TRAN,ITRAN)
292              MXEOM2      = MAX(MXEOM2,IVEC)
293
294              WL0Rf    = EOML0CONS(IVEC,ITRAN)
295              !write(lupri,*) "WL0Rf:", WL0Rf
296              WLiXiRf = WLiXi*WL0Rf
297            ELSE IF(.NOT.LEOMXOPA) THEN
298              IR1VEC = IR1TAMP(LABEL,LORX,EIGVI-EIGVF,IDUM)
299              CALL CC_SETF12(IFTRAN,IFDOTS,MXFTRAN,MXFVEC,
300     &                       ISTATEI,ISTATEF,IR1VEC,ITRAN,IVEC)
301              NFTRAN = MAX(NFTRAN,ITRAN)
302              MFVEC  = MAX(MFVEC, IVEC)
303              WIBF   = FCONS(IVEC,ITRAN)
304            END IF
305            !end if
306          END IF
307
308*---------------------------------------------------------------------*
309*          add contributions together:
310*---------------------------------------------------------------------*
311           IF (LADD) THEN
312
313              ITMIF = (NQR2OP*(IRSD-1) + IOP-1)*2 + ITURN
314
315              RESULT(ITMIF) = WIAF + WXINIF + WIBF - WLIXiRF
316
317              IF (LOCDBG) THEN
318                WRITE (LUPRI,*) '----- Summary after add ------'
319                WRITE (LUPRI,*) 'IDX of result = ',ITMIF
320                WRITE (LUPRI,*) 'ISTATEI, EIGVI:',ISTATEI,EIGVI
321                WRITE (LUPRI,*) 'ISTATEF, EIGVF:',ISTATEF,EIGVF
322                WRITE (LUPRI,*) 'OPERATOR:',LABEL
323                WRITE (LUPRI,*) 'L^i A{X} x R^f :',WIAF
324                WRITE (LUPRI,*) 'N^if x Xksi{X}:',WXINIF
325                WRITE (LUPRI,*) 'L^i x B x R^f x R^X:',WIBF
326                WRITE (LUPRI,*) 'L^i x Xksi{X}:',WLiXi
327                WRITE (LUPRI,*) 'R^f x L0:',WL0Rf
328                WRITE (LUPRI,*) '(L^i x Xksi{X})(R^f x L0):',WLiXiRF
329                WRITE (LUPRI,*) 'Total result:',RESULT(ITMIF)
330                WRITE (LUPRI,*) '------------------------------'
331              END IF
332
333           END IF
334
335*---------------------------------------------------------------------*
336*       end loop over transition moments
337*---------------------------------------------------------------------*
338
339        END IF
340       END DO
341      END DO
342
343      IF      (MFVEC.GT.MXFVEC) THEN
344         CALL QUIT('MFVEC has been out of bounds CCXOPA_EOMSETUP.')
345      ELSE IF (MXEAVEC.GT.MXAVEC) THEN
346         CALL QUIT('MXEAVEC has been out of bounds CCXOPA_EOMSETUP.')
347      ELSE IF (MXE2VEC.GT.MXXVEC) THEN
348         CALL QUIT('MXE2VEC has been out of bounds CCXOPA_EOMSETUP.')
349      ELSE IF (MXEOME2VEC.GT.MXEOMXVEC) THEN
350         CALL QUIT('MXEOME2VEC has been out of bounds CCXOPA_EOMSETUP.')
351      ELSE IF (MXEOM2.GT.MXEOML0VEC) THEN
352         CALL QUIT('MXEOM2 has been out of bounds CCXOPA_EOMSETUP.')
353      ELSE IF (NFTRAN.GT.MXFTRAN) THEN
354         CALL QUIT('NFTRAN has been out of bounds CCXOPA_EOMSETUP.')
355      ELSE IF (NXE1TRAN.GT.MXATRAN) THEN
356         CALL QUIT('NXE1TRAN has been out of bounds CCXOPA_EOMSETUP.')
357      ELSE IF (NXE2TRAN.GT.MXXTRAN) THEN
358         CALL QUIT('NXE2TRAN has been out of bounds CCXOPA_EOMSETUP.')
359      ELSE IF (NEOMXE2TRAN.GT.MXEOMXTRAN) THEN
360         CALL QUIT('NEOMXE2TRAN has been out bounds CCXOPA_EOMSETUP.')
361      ELSE IF (NEOML0TRAN.GT.MXEOML0TRAN) THEN
362         CALL QUIT('NEOML0TRAN has been out bounds CCXOPA_EOMSETUP.')
363      END IF
364
365*---------------------------------------------------------------------*
366* print the lists:
367*---------------------------------------------------------------------*
368* general statistics:
369      IF ((.NOT.LADD) .OR. LOCDBG) THEN
370       WRITE(LUPRI,'(/,/3X,A,I3,A)') 'For the requested',NBOPA,
371     &      ' transition moments'
372       WRITE(LUPRI,'((8X,A,I3,A))')
373     & ' - ',NFTRAN,  ' F matrix transformations with RE vectors',
374     & ' - ',NXE1TRAN,' A{X} matrix transformations with LE vectors',
375     & ' - ',NXE2TRAN,' extra XKSI vector calculations ',
376     & ' - ',NEOMXE2TRAN,' extra XKSI vector calculations (EOM)',
377     & ' - ',NEOML0TRAN,' extra L0 vector (EOM)'
378       WRITE(LUPRI,'(3X,A,/,/)') 'will be performed.'
379      END IF
380
381      IF (LOCDBG) THEN
382
383         ! F matrix transformations:
384         WRITE(LUPRI,*)'List of F matrix transformations:'
385         DO ITRAN = 1, NFTRAN
386           WRITE(LUPRI,'(A,2I5,5X,(25I3,20X))') MSGDBG,
387     &      (IFTRAN(I,ITRAN),I=1,2),(IFDOTS(I,ITRAN),I=1,MFVEC)
388         END DO
389         WRITE(LUPRI,*)
390
391         ! LE x A{X} vector calculations:
392         WRITE(LUPRI,*) 'List of A{O} matrix transformations:'
393         DO ITRAN = 1, NXE1TRAN
394           WRITE(LUPRI,'(A,5I5,5X,(25I3,20X))') MSGDBG,
395     &      (IEATRAN(I,ITRAN),I=1,5),(IEADOTS(I,ITRAN),I=1,MXEAVEC)
396         END DO
397         WRITE(LUPRI,*)
398
399         ! extra Xi{O} vector calculations:
400         WRITE(LUPRI,*) 'List of extra Xi{O} vector calculations:'
401         DO ITRAN = 1, NXE2TRAN
402           WRITE(LUPRI,'(A,5I5,5X,(25I3,20X))') MSGDBG,
403     &      (IXE2TRAN(I,ITRAN),I=1,5),(IX2DOTS(I,ITRAN),I=1,MXE2VEC)
404         END DO
405         WRITE(LUPRI,*)
406
407         ! extra Xi{O} vector calculations for EOM:
408         WRITE(LUPRI,*) 'List of EOM Xi{O} vector calculations:'
409         DO ITRAN = 1, NEOMXE2TRAN
410           WRITE(LUPRI,'(A,5I5,5X,(25I3,20X))') MSGDBG,
411     &      (IEOMXE2TRAN(I,ITRAN),I=1,5),
412     &             (IEOMX2DOTS(I,ITRAN),I=1,MXEOME2VEC)
413         END DO
414         WRITE(LUPRI,*)
415
416
417         !extra L0 x RE vector dot products:
418         WRITE (LUPRI,*) 'List of L0 x RE dot products:'
419         DO ITRAN = 1, NEOML0TRAN
420           WRITE(LUPRI,'(A,I5,5X,(12I5,20X))') MSGDBG,
421     &     IEOML0TRAN(ITRAN),(IEOML0DOTS(I,ITRAN),I=1,MXEOML0VEC)
422         END DO
423         WRITE (LUPRI,*)
424
425      END IF
426
427      RETURN
428      END
429
430*---------------------------------------------------------------------*
431*              END OF SUBROUTINE CCXOPA_EOMSETUP                      *
432*---------------------------------------------------------------------*
433