1C
2C  /* Deck dc_calc */
3      SUBROUTINE DC_CALC(MODEL,ISYMTR,NEXCI,EXVAL,LEXVAL,
4     &                   DENSIJ,LDENSIJ,DENSAB,LDENSAB,
5     &                   DENSAI,LDENSAI,
6     &                   T2MP,LT2MP,FOCKD,LFOCKD,
7     &                   WORK,LWORK)
8C
9C     This routine is part of the atomic integral direct SOPPA program.
10C
11C     Keld Bak, June 1997
12C     Stephan P. A. Sauer: 10.11.2003: merge with Dalton 2.0
13C
14C     PURPOSE: Determine "double corrected RPA" excitation energies.
15C
16#ifdef VAR_MPI
17      use so_parutils, only: parsoppa_do_eres, my_mpi_integer,
18     &                       soppa_nint, sop_master, one_mpi
19#endif
20      use so_info, only: sop_excita, so_full_name, so_model_number
21C
22#include "implicit.h"
23#ifdef VAR_MPI
24#include "mpif.h"
25C  IRAT in order to assign space for load-balancing
26#include "iratdef.h"
27#endif
28#include "priunit.h"
29C
30#include "soppinf.h"
31#include "ccsdsym.h"
32C
33      PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0)
34C
35      LOGICAL   NONEWT
36      CHARACTER*8 PDENS_LABEL
37CPi
38      CHARACTER*5 MODEL
39      CHARACTER*11 FULL_NAME
40C
41      DIMENSION EXVAL(LEXVAL)
42      DIMENSION DENSIJ(LDENSIJ), DENSAB(LDENSAB), DENSAI(LDENSAI)
43      DIMENSION T2MP(LT2MP),     FOCKD(LFOCKD),   WORK(LWORK)
44C
45      INTEGER NIT, NOLDTR, NNEWTR, IDTYPE, IMODEL
46#ifdef VAR_MPI
47      INTEGER   CP_ISYMTR
48      INTEGER   MAXNUMJOBS
49C     This array is only there to ensure that the four above variables
50C     are allocated consecutively, so that it can be send together. Only
51C     use it for this purpose.
52C     The definition must match that in soppa_nodedriver
53      INTEGER   INFO_ARRAY(6)
54      EQUIVALENCE (info_array(1), cp_isymtr), (info_array(2),nit),
55     &            (info_array(3), nnewtr),    (info_array(4),noldtr),
56     &            (info_array(5), idtype), (info_array(6),imodel)
57      INTEGER(MPI_INTEGER_KIND) :: ierr_mpi, numprocs_mpi
58C      IF (numprocs_mpi .GT. 1) CALL PARQUIT('RPA(D)')
59      IDTYPE = 0
60#endif
61C
62C------------------
63C     Add to trace.
64C------------------
65C
66      CALL QENTER('DC_CALC')
67C
68      FULL_NAME = SO_FULL_NAME(MODEL)
69      IMODEL = SO_MODEL_NUMBER(MODEL)
70C
71      KEND1 = 1
72      LWORK1 = LWORK
73C
74C--------------------------------------------------------------
75C     Initialize iteration counter, number of old trialvectors,
76C     and number of new trialvectors.
77C--------------------------------------------------------------
78C
79      NIT    = 1
80C
81      NOLDTR = 0
82C
83      NNEWTR = NEXCI
84C
85#ifdef VAR_MPI
86C------------------------------------------------------------------
87C     For MPI, we need some space in which to store the indices each
88C     process is to work with in so_eres.
89C------------------------------------------------------------------
90C
91      call mpi_comm_size( mpi_comm_world, numprocs_mpi, ierr_mpi)
92      maxnumjobs = soppa_nint - min(soppa_nint, numprocs_mpi) + 1
93      if ( numprocs_mpi .eq. 1 ) then
94C Not a real parallel job, don't bother
95         lAssignedIndices = 1
96         kAssignedIndices = 0
97      else
98         lAssignedIndices = (maxnumjobs + irat-1) /IRAT
99         kAssignedIndices = KEND1
100         KEND1 = kAssignedIndices + lAssignedIndices
101         LWORK1 = LWORK - KEND1
102         CALL SO_MEMMAX ('DC_CALC.1A',LWORK1)
103         IF (LWORK1 .LT. 0) CALL STOPIT('DC_CALC.1A',' ',KEND1,LWORK)
104      endif
105#endif
106C------------------
107C     Write banner.
108C------------------
109CPi
110      WRITE(LUPRI,'(/,2X,A)') '*********************************'//
111     &                        '*********************************'
112         WRITE(LUPRI,'(11X,A,A,I1)') ADJUSTR(FULL_NAME),
113     &        ' calculation, Excitation symmetry ',ISYMTR
114      WRITE(LUPRI,'(2X,A)') '*********************************'//
115     &                        '*********************************'
116C
117C---------------------------------------------------------
118C     Make E[2] linear transformation of RPA eigenvectors.
119C---------------------------------------------------------
120C
121#ifdef VAR_MPI
122C In parallel, send slaves to so_eres
123C
124         call mpi_bcast( parsoppa_do_eres, one_mpi,
125     &                   my_mpi_integer, sop_master,
126     &                   mpi_comm_world, ierr_mpi )
127C ISYMTR is a non-local parameter, we need to copy it to the info-array
128         CP_ISYMTR = ISYMTR
129         CALL MPI_BCAST( INFO_ARRAY, 6_mpi_integer_kind, MY_MPI_INTEGER,
130     &                   sop_master, MPI_COMM_WORLD, ierr_mpi)
131#endif
132CPi 18.08.16: Copy from so_lrsolv
133      CALL GETTIM (DUMMY,WTIMES)
134      DTIME      = SECOND()
135      CALL SO_ERES(MODEL,0,NEXCI,DENSIJ,LDENSIJ,DENSAB,LDENSAB,
136     &             T2MP,LT2MP,FOCKD,LFOCKD,DENSAI,LDENSAI,1,ISYMTR,
137     &             0,
138#ifdef VAR_MPI
139     &             WORK(kAssignedIndices), maxnumjobs,
140#endif
141     &             WORK(KEND1),LWORK1)
142      DTIME      = SECOND()   - DTIME
143      SOTIME(35) = SOTIME(35) + DTIME
144      CALL GETTIM (DUMMY,WTIMEE)
145      SOWTIM(1)  = SOWTIM(1)  + WTIMEE - WTIMES
146C
147C-----------------------------------------------------------
148C     Make S[2] linear transformation of trialvectors giving
149C     resultvectors.
150C-----------------------------------------------------------
151C
152      DTIME      = SECOND()
153      CALL DC_SRES(0,NEXCI,DENSIJ,LDENSIJ,DENSAB,LDENSAB,
154     &             ISYMTR,WORK,LWORK)
155      DTIME      = SECOND()   - DTIME
156      SOTIME(40) = SOTIME(40) + DTIME
157C
158C-----------------------------------------------------------------
159C     Calculate diagonal parts of E[2] and S[2] and write to disk.
160C-----------------------------------------------------------------
161C
162      DTIME      = SECOND()
163      CALL SO_DIAG(MODEL,FOCKD,LFOCKD,DENSIJ,LDENSIJ,DENSAB,LDENSAB,
164     &             ISYMTR,WORK,LWORK)
165      DTIME      = SECOND()   - DTIME
166      SOTIME(31) = SOTIME(31) + DTIME
167C
168C------------------------------
169C     Allocation of work space.
170C------------------------------
171C
172      LTR1E = NT1AM(ISYMTR)
173      LTR1D = NT1AM(ISYMTR)
174      LTR2E = N2P2HOP(ISYMTR)
175      LTR2D = N2P2HOP(ISYMTR)
176C
177      KTR1   = 1
178      KRES1  = KTR1   + LTR1E
179      KRESO1 = KRES1  + LTR1E
180      KRES2  = KRESO1 + LTR1E
181      KEND1  = KRES2  + LTR2E
182      LWORK1 = LWORK  - KEND1
183C
184      CALL SO_MEMMAX ('DC_CALC',LWORK1)
185      IF (LWORK1 .LT. 0) CALL STOPIT('DC_CALC',' ',KEND1,LWORK)
186C
187C----------------
188C     Open files.
189C----------------
190C
191      CALL SO_OPEN(LUTR1E,FNTR1E,LTR1E)
192      CALL SO_OPEN(LUTR1D,FNTR1D,LTR1D)
193      CALL SO_OPEN(LUTR2E,FNTR2E,LTR2E)
194      CALL SO_OPEN(LUTR2D,FNTR2D,LTR2D)
195      CALL SO_OPEN(LURS1E,FNRS1E,LTR1E)
196      CALL SO_OPEN(LURS1D,FNRS1D,LTR1D)
197      CALL SO_OPEN(LURO1E,FNRO1E,LTR1E)
198      CALL SO_OPEN(LURO1D,FNRO1D,LTR1D)
199      CALL SO_OPEN(LURS2E,FNRS2E,LTR2E)
200      CALL SO_OPEN(LURS2D,FNRS2D,LTR2D)
201C
202C---------------------------
203C     Loop over excitations.
204C---------------------------
205      WRITE(LUPRI,9010)
206      WRITE(LUPRI,9011)
207      WRITE(LUPRI,9010)
208C
209      DO 100 IEXCI = 1,NEXCI
210C
211C-----------------------------------------------------------------------
212C        Calculate 1p-1h part of double corrected RPA excitation energy.
213C-----------------------------------------------------------------------
214C
215         CALL SO_READ(WORK(KTR1),  LTR1E,LUTR1E,FNTR1E,IEXCI)
216         CALL SO_READ(WORK(KRES1), LTR1E,LURS1E,FNRS1E,IEXCI)
217         CALL SO_READ(WORK(KRESO1),LTR1E,LURO1E,FNRO1E,IEXCI)
218C
219         IF (MODEL .EQ. 'DCRPA' .OR. MODEL .EQ. 'SDCHR') THEN
220C
221            CALL DAXPY(LTR1E,-EXVAL(IEXCI),WORK(KRESO1),1,WORK(KRES1),1)
222C
223         END IF
224C
225         EX1P1H = DDOT(LTR1E,WORK(KRES1),1,WORK(KTR1),1)
226C
227C-----------------------------------------------------------------------
228C        Calculate 1h-1p part of double corrected RPA excitation energy.
229C-----------------------------------------------------------------------
230C
231         CALL SO_READ(WORK(KTR1),  LTR1D,LUTR1D,FNTR1D,IEXCI)
232         CALL SO_READ(WORK(KRES1), LTR1D,LURS1D,FNRS1D,IEXCI)
233         CALL SO_READ(WORK(KRESO1),LTR1D,LURO1D,FNRO1D,IEXCI)
234C
235         IF (MODEL .EQ. 'DCRPA' .OR. MODEL .EQ. 'SDCHR') THEN
236C
237            CALL DAXPY(LTR1D,-EXVAL(IEXCI),WORK(KRESO1),1,WORK(KRES1),1)
238C
239         END IF
240C
241         EX1H1P = DDOT(LTR1D,WORK(KRES1),1,WORK(KTR1),1)
242C
243C-----------------------------------------------------------------------
244C        Calculate 2p-2h part of double corrected RPA excitation energy.
245C-----------------------------------------------------------------------
246C
247         CALL SO_READ(WORK(KRES2), LTR2E,LURS2E,FNRS2E,IEXCI)
248C
249         IF (TRIPLET) THEN
250C
251            CALL DC_OMEC(EX2P2H,WORK(KRES2),LTR2E,EXVAL(IEXCI),
252     &                ISYMTR,WORK(KEND1),LWORK1)
253C
254         ELSE
255C
256            CALL CC_OMEC(EX2P2H,WORK(KRES2),EXVAL(IEXCI),
257     &                WORK(KEND1),LWORK1,ISYMTR)
258C
259C--------------------------------------------------------------------
260C        Calculate the first order 2p-2h eigenvector in RPA(D) theory
261C        and write to output.
262C--------------------------------------------------------------------
263C
264            CALL DC_R1VEC(WORK(KRES2), LTR2E,EXVAL(IEXCI),ISYMTR,
265     &                    WORK(KEND1),LWORK1)
266C
267         END IF
268C
269         CALL SO_WRITE(WORK(KRES2), LTR2E, LUTR2E,FNTR2E, IEXCI)
270C
271C-----------------------------------------------------------------------
272C        Calculate 2h-2p part of double corrected RPA excitation energy.
273C-----------------------------------------------------------------------
274C
275         CALL SO_READ(WORK(KRES2), LTR2D,LURS2D,FNRS2D,IEXCI)
276C
277         IF (TRIPLET) THEN
278C
279            CALL DC_OMEC(EX2H2P,WORK(KRES2),LTR2D,EXVAL(IEXCI),
280     &                ISYMTR,WORK(KEND1),LWORK1)
281C
282         ELSE
283C
284            CALL CC_OMEC(EX2H2P,WORK(KRES2),EXVAL(IEXCI),
285     &                WORK(KEND1),LWORK1,ISYMTR)
286C
287C--------------------------------------------------------------------
288C        Calculate the first order 2h-2p eigenvector in RPA(D) theory
289C        and write to output.
290C--------------------------------------------------------------------
291C
292            CALL DC_R1VEC(WORK(KRES2), LTR2D,EXVAL(IEXCI),ISYMTR,
293     &                    WORK(KEND1),LWORK1)
294C
295         END IF
296C
297         CALL SO_WRITE(WORK(KRES2), LTR2D, LUTR2D,FNTR2D, IEXCI)
298C
299C-----------------------------------------------------------
300C        Add contributions to give RPA(D) excitation energy.
301C-----------------------------------------------------------
302C
303         IF (TRIPLET) THEN
304C
305            WRITE(LUPRI,9012) IEXCI,
306     &                        EX1P1H + EX1H1P + EX2P2H + EX2H2P,
307     &                        EXVAL(IEXCI),
308     &                        (EX1P1H + EX1H1P) - EXVAL(IEXCI),
309     &                        EX2P2H + EX2H2P
310C
311            EXVAL(IEXCI) = EX1P1H + EX1H1P + EX2P2H + EX2H2P
312C
313         ELSE
314C
315            WRITE(LUPRI,9012) IEXCI,
316     &                        EX1P1H + EX1H1P + (EX2P2H + EX2H2P)/TWO,
317     &                        EXVAL(IEXCI),
318     &                        (EX1P1H + EX1H1P)- EXVAL(IEXCI),
319     &                        (EX2P2H + EX2H2P)/TWO
320C
321            EXVAL(IEXCI) = EX1P1H + EX1H1P + EX2P2H/TWO + EX2H2P/TWO
322C
323         END IF
324C
325  100 CONTINUE
326C
327      WRITE(LUPRI,9010)
328C
329C-----------------
330C     Close files.
331C-----------------
332C
333      CALL SO_CLOSE(LUTR1E,FNTR1E,'KEEP')
334      CALL SO_CLOSE(LUTR1D,FNTR1D,'KEEP')
335      CALL SO_CLOSE(LUTR2E,FNTR2E,'KEEP')
336      CALL SO_CLOSE(LUTR2D,FNTR2D,'KEEP')
337      CALL SO_CLOSE(LURS1E,FNRS1E,'KEEP')
338      CALL SO_CLOSE(LURS1D,FNRS1D,'KEEP')
339      CALL SO_CLOSE(LURO1E,FNRO1E,'KEEP')
340      CALL SO_CLOSE(LURO1D,FNRO1D,'KEEP')
341      CALL SO_CLOSE(LURS2E,FNRS2E,'KEEP')
342      CALL SO_CLOSE(LURS2D,FNRS2D,'KEEP')
343C
344C--------------------------------------------------------------
345C     Calculate and save the pertubed density matrix
346C-------------------------------------------------------------
347CPi Solution vectors are not normalized as in so_rpslex based on
348CPi excitation weights !!!
349C
350      WRITE(PDENS_LABEL,'(A7,I1)') 'EXCITA ',ISYMTR
351      CALL SO_PERTDENS('DCRPA',SOP_EXCITA,NEXCI,
352     &                 EXVAL,NEXCI,PDENS_LABEL,
353     &                 ISYMTR,.FALSE.,1.0D0,
354     &                 T2MP,LT2MP,DENSIJ,LDENSIJ,
355     &                 DENSAB,LDENSAB,DENSAI,LDENSAI,
356     &                 WORK(KEND1),LWORK1)
357C
358C
359C------------------------------------
360C     Flush the standard output file.
361C------------------------------------
362C
363      CALL FLSHFO(LUPRI)
364C
365C-----------------------
366C     Remove from trace.
367C-----------------------
368C
369      CALL QEXIT('DC_CALC')
370C
371      RETURN
372C
373 9010 FORMAT(1X,'----------------------------------------',
374     &       '-------------------------')
375 9011 FORMAT(1X,'Excitation   Energy (au)      ph(0)     ',
376     &       '   ph(2)       2p2h(2)')
377 9012 FORMAT(2X,I5,6X,F12.8,1X,F12.8,1X,F12.8,1X,F12.8)
378C
379      END
380