1!
2!...   Copyright (c) 2015 by the authors of Dalton (see below).
3!...   All Rights Reserved.
4!...
5!...   The source code in this file is part of
6!...   "Dalton, a molecular electronic structure program,
7!...    Release DALTON2015 (2015), see http://daltonprogram.org"
8!...
9!...   This source code is provided under a written licence and may be
10!...   used, copied, transmitted, or stored only in accord with that
11!...   written licence.
12!...
13!...   In particular, no part of the source code or compiled modules may
14!...   be distributed outside the research group of the licence holder.
15!...   This means also that persons (e.g. post-docs) leaving the research
16!...   group of the licence holder may not take any part of Dalton,
17!...   including modified files, with him/her, unless that person has
18!...   obtained his/her own licence.
19!...
20!...   For further information, including how to get a licence, see:
21!...      http://daltonprogram.org
22!
23!
24C
25c /* deck CC_EOM_XOPA */
26*=====================================================================*
27       SUBROUTINE CC_EOM_XOPA(WORK,LWORK)
28*---------------------------------------------------------------------*
29*
30*    Purpose: direct calculation of first-order transition properties
31*             (transition moments and oscillator strengths)
32*             for transitions between two excited states with the
33*             Coupled Cluster models
34*
35*                        CCS, CC2, CCSD, CC3
36*
37*             and partially with SCF and CIS
38*
39*     Written by Christof Haettig winter 2002/2003.
40*     Modified version by Sonia, to include EOM models, 2015
41*=====================================================================*
42      IMPLICIT NONE
43#include "priunit.h"
44#include "cclists.h"
45#include "ccxopainf.h"
46#include "ccsdinp.h"
47#include "dummy.h"
48#include "second.h"
49#include "ccexcinf.h"
50#include "ccorb.h"
51
52* local parameters:
53      CHARACTER*(16) MSGDBG
54      PARAMETER (MSGDBG = '[debug] CC_EOM_XOPA> ')
55
56#if defined (SYS_CRAY)
57      REAL ZERO
58#else
59      DOUBLE PRECISION ZERO
60#endif
61      PARAMETER (ZERO = 0.0d0)
62
63      CHARACTER*10 MODEL
64      INTEGER LWORK
65
66#if defined (SYS_CRAY)
67      REAL WORK(LWORK)
68      REAL TIM0, TIM1, TIMF, TIMXE1, TIMXE2
69#else
70      DOUBLE PRECISION WORK(LWORK)
71      DOUBLE PRECISION TIM0, TIM1, TIMF, TIMXE1, TIMXE2
72#endif
73
74      LOGICAL LADD
75      INTEGER NBOPA, MXFTRAN, MXATRAN, MXXTRAN, MXFVEC, MXAVEC, MXXVEC,
76     &        NFTRAN, NXE1TRAN, NXE2TRAN, NSTATES,
77     &        KRESULT, KFTRAN, KFDOTS, KFCONS, KEND0, LEND0,
78     &        KE1TRAN, KE1DOTS, KE1CONS,
79     &        KX2TRAN, KX2DOTS, KX2CONS,
80     &        IOPT, IORDER, ISYM
81      INTEGER KEOMX2TRAN, KEOMX2DOTS, KEOMX2CONS
82      INTEGER KEOML0TRAN, KEOML0DOTS, KEOML0CONS
83      INTEGER NEOMXE2TRAN, MXEOMXTRAN, MXEOMXVEC
84      INTEGER NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC
85
86* external functions: none
87
88*---------------------------------------------------------------------*
89* print header for second-order property section:
90*---------------------------------------------------------------------*
91      WRITE (LUPRI,'(7(/1X,2A),/)')
92     & '************************************',
93     &                               '******************************',
94     & '*                                   ',
95     &                               '                             *',
96     & '*<<<<<<    OUTPUT FROM COUPLED CLUST',
97     &                               'ER LINEAR RESPONSE    >>>>>>>*',
98     & '*<<<<<<  CALCULATION OF ONE-PHOTON A',
99     &                               'BSORPTION STRENGTHS  >>>>>>>*',
100     & '*<<<<<<     FOR EXCITED TO EXCITED S',
101     &                               'TATE TRANSITIONS      >>>>>>>*',
102     & '*                                   ',
103     &                               '                             *',
104     & '************************************',
105     &                               '******************************'
106
107*---------------------------------------------------------------------*
108      IF (.NOT. (CCS .OR. CC2 .OR. CCSD .OR. CC3) ) THEN
109         CALL QUIT('CC_EOM_XOPA called for unknown Coupled Cluster.')
110      END IF
111
112* print some debug/info output
113      IF(IPRINT .GT. 10) WRITE(LUPRI,*) 'CC_EOM_XOPA Workspace:',LWORK
114
115      TIM0  = SECOND()
116
117*---------------------------------------------------------------------*
118* allocate & initialize work space for property contributions:
119*---------------------------------------------------------------------*
120      ! maximum number of transition moments to compute
121      NBOPA   = 2 * NQR2OP * NXQR2ST
122
123      ! number of excited states
124      NSTATES = 0
125      DO ISYM = 1, NSYM
126        NSTATES = NSTATES + NCCEXCI(ISYM,1) + NCCEXCI(ISYM,3)
127      END DO
128
129      ! maximum number of transformations or vector calculations
130      ! NSTATES * NQR2OP   LE x Eta{X} transformations
131      ! NQR2OP             Xi{X} vectors
132      ! 2*NXQR2ST          LE x B x RE transformations
133      MXATRAN    = NSTATES * NQR2OP
134      MXXTRAN    = NQR2OP
135      MXEOMXTRAN = NQR2OP
136      MXEOML0TRAN  = 1
137      MXFTRAN    = 2*NXQR2ST
138
139      ! maximum number of vectors to dot on
140      ! NSTATES    RE vectors dotted on a LE x Eta{X} transformation
141      ! 2*NXQR2ST  N2 vectors dotted on a Xi{X} vector
142      ! NQR2OP     R1 vectors dotted on a LE x B x RE transformation
143      MXAVEC      = NSTATES
144      MXXVEC      = 2*NXQR2ST
145      MXEOMXVEC   = 2*NXQR2ST !I am not sure I understand why this number...
146      MXFVEC      = NQR2OP
147      MXEOML0VEC  = NSTATES
148
149      KRESULT  = 1
150      KEND0    = KRESULT  + NBOPA
151
152      KFTRAN   = KEND0
153      KFDOTS   = KFTRAN   + MXFTRAN * MXDIM_FTRAN
154      KFCONS   = KFDOTS   + MXFVEC  * MXFTRAN
155      KEND0    = KFCONS   + MXFVEC  * MXFTRAN
156
157      KE1TRAN  = KEND0
158      KE1DOTS  = KE1TRAN  + MXATRAN * MXDIM_XEVEC
159      KE1CONS  = KE1DOTS  + MXAVEC  * MXATRAN
160      KEND0    = KE1CONS  + MXAVEC  * MXATRAN
161
162      KX2TRAN  = KEND0
163      KX2DOTS  = KX2TRAN  + MXXTRAN * MXDIM_XEVEC
164      KX2CONS  = KX2DOTS  + MXXVEC  * MXXTRAN
165      KEND0    = KX2CONS  + MXXVEC  * MXXTRAN
166
167      KEOMX2TRAN= KEND0
168      KEOMX2DOTS= KEOMX2TRAN  + MXEOMXTRAN * MXDIM_XEVEC
169      KEOMX2CONS= KEOMX2DOTS  + MXEOMXVEC * MXEOMXTRAN
170      KEND0     = KEOMX2CONS  + MXEOMXVEC * MXEOMXTRAN
171
172      KEOML0TRAN= KEND0
173      KEOML0DOTS= KEOML0TRAN  + MXEOML0TRAN
174      KEOML0CONS= KEOML0DOTS  + MXEOML0VEC * MXEOML0TRAN
175      KEND0     = KEOML0CONS  + MXEOML0VEC * MXEOML0TRAN
176
177      LEND0 = LWORK - KEND0
178      IF (LEND0 .LT. 0) THEN
179        CALL QUIT('Insufficient memory in CC_EOM_XOPA. (1)')
180      END IF
181
182      CALL DZERO(WORK(KRESULT),NBOPA)
183
184      !sonia
185      CALL DZERO(WORK(KEOMX2CONS),MXEOMXVEC * MXEOMXTRAN)
186
187*---------------------------------------------------------------------*
188* set up lists for F transformations, ETA{O} and Xi{O} vectors:
189*---------------------------------------------------------------------*
190      LADD = .FALSE.
191
192      CALL CCXOPA_EOMSETUP(WORK(KFTRAN),WORK(KFDOTS),WORK(KFCONS),
193     &                  NFTRAN,MXFTRAN,MXFVEC,
194     &                  WORK(KE1TRAN),WORK(KE1DOTS),WORK(KE1CONS),
195     &                  NXE1TRAN,MXATRAN,MXAVEC,
196     &                  WORK(KX2TRAN),WORK(KX2DOTS),WORK(KX2CONS),
197     &                  NXE2TRAN,MXXTRAN,MXXVEC,
198     &                  WORK(KEOMX2TRAN),WORK(KEOMX2DOTS),
199     &                  WORK(KEOMX2CONS),
200     &                  NEOMXE2TRAN,MXEOMXTRAN,MXEOMXVEC,
201     &                  WORK(KEOML0TRAN),WORK(KEOML0DOTS),
202     &                  WORK(KEOML0CONS),
203     &                  NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC,
204     &                  WORK(KRESULT),NBOPA,LADD,WORK(KEND0),LEND0)
205
206*---------------------------------------------------------------------*
207* calculate F matrix contributions:
208*---------------------------------------------------------------------*
209      TIM1 = SECOND()
210
211      CALL DZERO(WORK(KFCONS),MXFVEC*NFTRAN)
212
213      IOPT = 5
214      CALL CC_FMATRIX(WORK(KFTRAN),NFTRAN,'LE ','RE ',IOPT,'R1 ',
215     &                WORK(KFDOTS),WORK(KFCONS),MXFVEC,
216     &                WORK(KEND0), LEND0)
217
218      TIMF = SECOND() - TIM1
219
220      IF (NFTRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")')
221     & '>>> Time used for',NFTRAN,' F matrix transformations:',TIMF
222      CALL FLSHFO(LUPRI)
223
224*---------------------------------------------------------------------*
225* calculate LE x A{O} x RE contributions:
226*---------------------------------------------------------------------*
227      TIM1 = SECOND()
228
229      CALL DZERO(WORK(KE1CONS),MXAVEC*NXE1TRAN)
230
231
232
233      IOPT   = 5
234      IORDER = 1
235      if (leomxopa) then
236         CALL CCEOM_XIETA( WORK(KE1TRAN),NXE1TRAN,IOPT,IORDER,'LE ',
237     &                 '---',IDUMMY,       DUMMY,
238     &                 'RE ',WORK(KE1DOTS),WORK(KE1CONS),
239     &                       MXAVEC, WORK(KEND0), LEND0 )
240      else
241       CALL CC_XIETA( WORK(KE1TRAN), NXE1TRAN, IOPT, IORDER, 'LE ',
242     &               '---',DUMMY,DUMMY,
243     &               'RE ',WORK(KE1DOTS),WORK(KE1CONS),
244     &               .FALSE.,MXAVEC, WORK(KEND0), LEND0 )
245      end if
246
247      TIMXE1 = SECOND() - TIM1
248      IF (NXE1TRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")')
249     & '>>> Time used for',NXE1TRAN,' A{X} matrix transformations:',
250     & TIMXE1
251      CALL FLSHFO(LUPRI)
252
253*---------------------------------------------------------------------*
254* calculate N2 x Xksi{O} vector contributions:
255*---------------------------------------------------------------------*
256      TIM1 = SECOND()
257
258      CALL DZERO(WORK(KX2CONS),MXXVEC*NXE2TRAN)
259
260      IOPT   = 5
261      IORDER = 1
262      CALL CC_XIETA( WORK(KX2TRAN), NXE2TRAN, IOPT, IORDER, '---',
263     &               'N2 ',WORK(KX2DOTS),WORK(KX2CONS),
264     &               '---',IDUMMY,DUMMY,
265     &               .FALSE.,MXXVEC, WORK(KEND0), LEND0 )
266
267      TIMXE2 = SECOND() - TIM1
268      IF (NXE2TRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")')
269     & '>>> Time used for',NXE2TRAN,' O1/X1 vector calculation:',TIMXE2
270      CALL FLSHFO(LUPRI)
271
272*---------------------------------------------------------------------*
273* calculate LE x Xksi{O} vector contributions:
274* calculate RE x Tbar0   vector contributions:
275* multiply them together to get the EOM contrib
276*---------------------------------------------------------------------*
277      TIM1 = SECOND()
278
279      if (leomxopa) then
280         CALL DZERO(WORK(KEOMX2CONS),MXEOMXVEC*NEOMXE2TRAN)
281         CALL DZERO(WORK(KEOML0CONS),MXEOML0VEC*NEOML0TRAN)
282
283         IOPT   = 5
284         IORDER = 1
285         CALL CC_XIETA( WORK(KEOMX2TRAN), NEOMXE2TRAN,
286     &                IOPT, IORDER, 'L0 ',
287     &               'LE ',WORK(KEOMX2DOTS),WORK(KEOMX2CONS),
288     &               '---',IDUMMY,DUMMY,
289     &               .FALSE.,MXEOMXVEC, WORK(KEND0), LEND0 )
290
291         CALL CC_DOTDRV('L0 ','RE ',NEOML0TRAN,MXEOML0VEC,
292     &                 WORK(KEOML0TRAN), WORK(KEOML0DOTS),
293     &                 WORK(KEOML0CONS),
294     &                 WORK(KEND0), LEND0 )
295
296         TIMXE2 = SECOND() - TIM1
297      IF (NEOMXE2TRAN.GT.0) WRITE (LUPRI,'(/A,I5,A,F12.2," seconds.")')
298     & '>>> Time used for',NEOMXE2TRAN,
299     &  ' O1/X1 vector calculation:',TIMXE2
300         CALL FLSHFO(LUPRI)
301      end if
302
303*---------------------------------------------------------------------*
304* collect contributions and sum them up to the final results:
305*---------------------------------------------------------------------*
306      LADD = .TRUE.
307
308      CALL CCXOPA_EOMSETUP(WORK(KFTRAN),WORK(KFDOTS),WORK(KFCONS),
309     &                  NFTRAN,MXFTRAN,MXFVEC,
310     &                  WORK(KE1TRAN),WORK(KE1DOTS),WORK(KE1CONS),
311     &                  NXE1TRAN,MXATRAN,MXAVEC,
312     &                  WORK(KX2TRAN),WORK(KX2DOTS),WORK(KX2CONS),
313     &                  NXE2TRAN,MXXTRAN,MXXVEC,
314     &                  WORK(KEOMX2TRAN),WORK(KEOMX2DOTS),
315     &                  WORK(KEOMX2CONS),
316     &                  NEOMXE2TRAN,MXEOMXTRAN,MXEOMXVEC,
317     &                  WORK(KEOML0TRAN),WORK(KEOML0DOTS),
318     &                  WORK(KEOML0CONS),
319     &                  NEOML0TRAN,MXEOML0TRAN,MXEOML0VEC,
320     &                  WORK(KRESULT),NBOPA,LADD,WORK(KEND0),LEND0)
321
322*---------------------------------------------------------------------*
323* print timing:
324*---------------------------------------------------------------------*
325      WRITE (LUPRI,'(/A,I4,A,F12.2," seconds.")') '>>> Total time for',
326     &  NBOPA,' quadratic response func.:', SECOND() - TIM0
327
328*---------------------------------------------------------------------*
329* print one-photon absorption properties and return:
330*---------------------------------------------------------------------*
331      CALL  CCOPAPRT(WORK(KRESULT),.TRUE.,NQR2OP,NXQR2ST)
332
333      CALL FLSHFO(LUPRI)
334
335      RETURN
336      END
337
338*=====================================================================*
339*              END OF SUBROUTINE CC_EOM_XOPA                              *
340*=====================================================================*
341