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!
18      SUBROUTINE CCSDT_AAMAT_NODDY(IOPTRES,FREQRES,LABELB,ISYMB,
19     &                             LISTC,IDLSTC,OEFF_INIT,
20     &                             OMEGA1,OMEGA2,
21     &                             OMEGA1EFF,OMEGA2EFF,
22     &                             IDOTS,DOTPROD,LISTDP,ITRAN,
23     &                             NXTRAN,MXVEC,WORK,LWORK)
24*---------------------------------------------------------------------*
25*
26*    Purpose: compute triples contribution to A{B} transformed vector
27*
28*    (A{B} T^C)^eff_1,2 = (A{B} T^C)_1,2(CCSD) + (A{B} T^C)_1,2(T3)
29*                            - A_1,2;3 (w_3 - w)^1 (A{B} T^C)_3
30*
31*
32*   Written by Christof Haettig, April 2002, based on CCSDT_XI_NODDY
33*
34*=====================================================================*
35      IMPLICIT NONE
36#include "priunit.h"
37#include "ccsdinp.h"
38#include "maxorb.h"
39#include "ccsdsym.h"
40#include "ccfield.h"
41#include "ccorb.h"
42
43      LOGICAL LOCDBG, PRINT_T3
44      PARAMETER (LOCDBG=.FALSE., PRINT_T3=.FALSE.)
45      INTEGER ISYM0
46      PARAMETER (ISYM0 = 1)
47
48      CHARACTER LISTDP*3, LABELB*8, LISTC*3
49      LOGICAL OEFF_INIT
50      INTEGER LWORK, IOPTRES, ITRAN, MXVEC, NXTRAN, ISYMB, IDLSTC
51      INTEGER IDOTS(MXVEC,NXTRAN)
52
53#if defined (SYS_CRAY)
54      REAL DOTPROD(MXVEC,NXTRAN), FREQC, FREQRES
55      REAL WORK(LWORK), ONE, TWO, DUMMY, FF, SIGN, DDOT
56      REAL OMEGA1(*),OMEGA2(*)
57      REAL OMEGA1EFF(*),OMEGA2EFF(*)
58#else
59      DOUBLE PRECISION DOTPROD(MXVEC,NXTRAN), FREQC, FREQRES
60      DOUBLE PRECISION WORK(LWORK), ONE, TWO, DUMMY, FF, SIGN, DDOT
61      DOUBLE PRECISION OMEGA1(*),OMEGA2(*)
62      DOUBLE PRECISION OMEGA1EFF(*),OMEGA2EFF(*)
63#endif
64      PARAMETER( ONE = 1.0D0, TWO = 2.0D0 )
65
66      CHARACTER*10 MODEL
67      INTEGER KFOCKB, KT3AM, KT2AM0, KT2AMC, KEND2, KEND1, KOMEGA1,
68     &        KOMEGA2, KOMEGA3, KSCR1, KFOCKD, KFOCK0, KT1AMC,
69     &        KFOCKB_AO, KFOCKBC, KLAMPC, KLAMHC, KFOCKC, KINT1SC,
70     &        KINT2SC, KINT1S, KINT2S, KXIAJB, KYIAJB, KINT1T, KINT2T,
71     &        LWRK1, KDUM, LWRK2, KT1AMP0, KLAMP0, KLAMH0
72      INTEGER IJ, NIJ, LUSIFC, INDEX, IDUMMY, ILSTSYM, ISYMC, LUFOCK,
73     &        IRREP, IERR, ILLL, IDEL, ISYDIS, IOPT, ISYMD, IVEC,
74     &        IDLSTD, KFCKBUF, KFIELD, KFIELDAO
75
76      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
77
78*---------------------------------------------------------------------*
79*     Begin:
80*---------------------------------------------------------------------*
81      CALL QENTER('CCSDT_AAMAT_NODDY')
82
83      IF (LOCDBG) WRITE(LUPRI,*) 'Entered CCSDT_AAMAT_NODDY...'
84
85      IF(DIRECT)CALL QUIT('DIRECT NOT IMPLEMENTED IN CCSDT_AAMAT_NODDY')
86
87*---------------------------------------------------------------------*
88*     Memory allocation:
89*---------------------------------------------------------------------*
90      KEND1   = 1
91
92      KOMEGA1 = KEND1
93      KOMEGA2 = KOMEGA1 + NT1AMX
94      KOMEGA3 = KOMEGA2 + NT1AMX*NT1AMX
95      KEND1   = KOMEGA3 + NT1AMX*NT1AMX*NT1AMX
96
97      KSCR1   = KEND1
98      KFOCKD  = KSCR1  + NT1AMX
99      KLAMP0  = KFOCKD + NORBT
100      KLAMH0  = KLAMP0 + NLAMDT
101      KFOCK0  = KLAMH0 + NLAMDT
102      KT1AMP0 = KFOCK0 + NORBT*NORBT
103      KEND1   = KT1AMP0+ NT1AMX
104
105
106      KFOCKB    = KEND1
107      KFOCKB_AO = KFOCKB    + NORBT*NORBT
108      KFOCKBC   = KFOCKB_AO + NORBT*NORBT
109      KFCKBUF   = KFOCKBC   + NORBT*NORBT
110      KEND1     = KFCKBUF   + NORBT*NORBT
111
112      IF (NONHF) THEN
113        KFIELD   = KEND1
114        KFIELDAO = KFIELD   + NORBT*NORBT
115        KEND1    = KFIELDAO + NORBT*NORBT
116      END IF
117
118      KLAMPC  = KEND1
119      KLAMHC  = KLAMPC + NLAMDT
120      KFOCKC  = KLAMHC + NLAMDT
121      KEND1   = KFOCKC + NORBT*NORBT
122
123      KINT1S  = KEND1
124      KINT2S  = KINT1S  + NT1AMX*NVIRT*NVIRT
125      KEND1   = KINT2S  + NRHFT*NRHFT*NT1AMX
126
127      KINT1SC = KEND1
128      KINT2SC = KINT1SC + NT1AMX*NVIRT*NVIRT
129      KEND1   = KINT2SC + NRHFT*NRHFT*NT1AMX
130
131      KXIAJB  = KEND1
132      KYIAJB  = KXIAJB  + NT1AMX*NT1AMX
133      KEND1   = KYIAJB  + NT1AMX*NT1AMX
134
135      KINT1T = KEND1
136      KINT2T = KINT1T + NT1AMX*NVIRT*NVIRT
137      KEND1  = KINT2T + NRHFT*NRHFT*NT1AMX
138
139      LWRK1  = LWORK  - KEND1
140      IF (LWRK1 .LT. 0) THEN
141         CALL QUIT('Insufficient space in CCSDT_AAMAT_NODDY')
142      ENDIF
143
144      CALL DZERO(WORK(KOMEGA1),NT1AMX)
145
146*---------------------------------------------------------------------*
147*     Get zeroth-order Lambda matrices:
148*---------------------------------------------------------------------*
149      IOPT   = 1
150      KDUM = KEND1
151      Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT1AMP0),WORK(KDUM))
152
153      Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT1AMP0),
154     &            WORK(KEND1),LWRK1)
155
156*---------------------------------------------------------------------*
157*     Read SCF orbital energies from file:
158*---------------------------------------------------------------------*
159      LUSIFC = -1
160      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
161     &            .FALSE.)
162      REWIND LUSIFC
163      CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
164      READ (LUSIFC)
165      READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBT)
166      CALL GPCLOSE(LUSIFC,'KEEP')
167
168*---------------------------------------------------------------------*
169*     read zeroth-order AO Fock matrix from file:
170*---------------------------------------------------------------------*
171      LUFOCK = -1
172      CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY,
173     &            .FALSE.)
174      REWIND(LUFOCK)
175      READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0))
176      CALL GPCLOSE(LUFOCK,'KEEP')
177
178      CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0),
179     &              WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0)
180
181*---------------------------------------------------------------------*
182*     If needed get external field:
183*---------------------------------------------------------------------*
184      IF ((NONHF) .AND. NFIELD .GT. 0) THEN
185         CALL DZERO(WORK(KFIELDAO),NORBT*NORBT)
186         DO I = 1, NFIELD
187          FF = EFIELD(I)
188          CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND1),LWRK1,FF,1,LFIELD(I))
189         ENDDO
190
191         CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFIELD),1)
192         CALL CC_FCKMO(WORK(KFIELD),WORK(KLAMP0),WORK(KLAMH0),
193     *                 WORK(KEND1),LWRK1,1,1,1)
194      ENDIF
195
196*---------------------------------------------------------------------*
197*     Matrix with property integrals in MO basis:
198*---------------------------------------------------------------------*
199      ! read property integrals from file:
200      CALL CCPRPAO(LABELB,.TRUE.,WORK(KFOCKB_AO),IRREP,ISYMB,IERR,
201     &             WORK(KEND1),LWRK1)
202      CALL DCOPY(NORBT*NORBT,WORK(KFOCKB_AO),1,WORK(KFOCKB),1)
203      IF ((IERR.GT.0) .OR. (IERR.EQ.0 .AND. IRREP.NE.ISYMB)) THEN
204        CALL QUIT('CCSDT_AAMAT_NODDY: error reading operator '//LABELB)
205      ELSE IF (IERR.LT.0) THEN
206        CALL DZERO(WORK(KFOCKB),N2BST(ISYMB))
207      END IF
208
209      ! transform property integrals to Lambda-MO basis
210      CALL CC_FCKMO(WORK(KFOCKB),WORK(KLAMP0),WORK(KLAMH0),
211     &              WORK(KEND1),LWRK1,ISYMB,1,1)
212
213*---------------------------------------------------------------------*
214*     Compute some integrals:
215*           XINT1S  =  (CK|BD)
216*           XINT2S  =  (CK|LJ)
217*           XINT1T  =  (KC|BD)
218*           XINT2T  =  (KC|LJ)
219*           XIAJB   = 2(IA|JB) - (IB|JA)
220*           YIAJB   =  (IA|JB)
221*---------------------------------------------------------------------*
222      CALL CCSDT_INTS0_NODDY(.TRUE.,WORK(KXIAJB),WORK(KYIAJB),
223     &                       .TRUE.,WORK(KINT1S),WORK(KINT2S),
224     &                       .TRUE.,WORK(KINT1T),WORK(KINT2T),
225     &                       WORK(KLAMP0),WORK(KLAMH0),
226     &                       WORK(KEND1),LWRK1)
227
228*---------------------------------------------------------------------*
229*     allocate work space for a triples and two doubles vectors and
230*     read zeroth-order and response doubles amplitudes from file,
231*     compute response Lambda matrices:
232*---------------------------------------------------------------------*
233
234      KT3AM  = KEND1
235      KT2AM0 = KT3AM  + NT1AMX*NT1AMX*NT1AMX
236      KT2AMC = KT2AM0 + NT1AMX*NT1AMX
237      KT1AMC = KT2AMC + NT1AMX*NT1AMX
238      KEND2  = KT1AMC + NT1AMX
239      LWRK2  = LWORK  - KEND2
240      IF (LWRK2 .LT. 0) THEN
241         CALL QUIT('Insufficient space in CCSDT_AAMAT_NODDY')
242      ENDIF
243
244      IOPT = 2
245      CALL CC_RDRSP('R0',0,ISYMOP,IOPT,MODEL,DUMMY,WORK(KT3AM))
246      CALL CC_T2SQ(WORK(KT3AM),WORK(KT2AM0),ISYMOP)
247
248      ISYMC = ILSTSYM(LISTC,IDLSTC)
249      IOPT = 3
250      CALL CC_RDRSP(LISTC,IDLSTC,ISYMC,IOPT,MODEL,
251     &              WORK(KT1AMC),WORK(KT3AM))
252      CALL CCLR_DIASCL(WORK(KT3AM),TWO,ISYMC)
253      CALL CC_T2SQ(WORK(KT3AM),WORK(KT2AMC),ISYMC)
254
255      CALL CCLR_LAMTRA(WORK(KLAMP0),WORK(KLAMPC),
256     &                 WORK(KLAMH0),WORK(KLAMHC),WORK(KT1AMC),ISYMC)
257
258*---------------------------------------------------------------------*
259*     compute zeroth-order triples response amplitudes and its
260*     contributions to the result vector
261*---------------------------------------------------------------------*
262
263C     --------------------------------------------------------------
264C     compute zeroth-order triples; note: CCSDT_T3AM uses in non-hf
265C     case WORK(KOMEGA3) as scratch vector for the iterative
266C     solution of the triples equations!
267C     --------------------------------------------------------------
268      CALL DZERO(WORK(KT3AM),NT1AMX*NT1AMX*NT1AMX)
269      CALL CCSDT_T03AM(WORK(KT3AM),WORK(KINT1S),WORK(KINT2S),
270     &                 WORK(KT2AM0),WORK(KSCR1),WORK(KFOCKD),
271     &                 WORK(KFIELD),WORK(KOMEGA3))
272
273C     --------------------------------------------------------------
274C     calculate one-index transf. property integrals: B^C = [B,T1^C]
275C     --------------------------------------------------------------
276      CALL DCOPY(NORBT*NORBT,WORK(KFOCKB_AO),1,WORK(KFCKBUF),1)
277      CALL DCOPY(NORBT*NORBT,WORK(KFOCKB_AO),1,WORK(KFOCKBC),1)
278      CALL CC_FCKMO(WORK(KFCKBUF),WORK(KLAMPC),WORK(KLAMH0),
279     &              WORK(KEND2),LWRK2,ISYMB,ISYMC,ISYM0)
280      CALL CC_FCKMO(WORK(KFOCKBC),WORK(KLAMP0),WORK(KLAMHC),
281     &              WORK(KEND2),LWRK2,ISYMB,ISYM0,ISYMC)
282      CALL DAXPY(NORBT*NORBT,ONE,WORK(KFCKBUF),1,WORK(KFOCKBC),1)
283
284C     ---------------------------------
285C     Initialize triples result vector:
286C     ---------------------------------
287      CALL DZERO(WORK(KOMEGA3),NT1AMX*NT1AMX*NT1AMX)
288
289C     ----------------------------------------------
290C     add triples contribution: <mu_3|[B^C,T3^0]|HF>
291C     ----------------------------------------------
292      CALL CCSDT_XKSI3_2(WORK(KOMEGA3),WORK(KFOCKBC),WORK(KT3AM))
293
294*---------------------------------------------------------------------*
295*     compute contribution from doubles amplitudes:
296*---------------------------------------------------------------------*
297
298C     ---------------------------------------------------
299C     add triples contribution: <mu_3|[[B,T2^C],T2^0]|HF>
300C     ---------------------------------------------------
301      CALL CCSDT_XKSI3_1(WORK(KOMEGA3),WORK(KFOCKB),
302     &                   WORK(KT2AM0),WORK(KT2AMC),ONE)
303      CALL CCSDT_XKSI3_1(WORK(KOMEGA3),WORK(KFOCKB),
304     &                   WORK(KT2AMC),WORK(KT2AM0),ONE)
305
306*---------------------------------------------------------------------*
307*     compute triples response amplitudes and its contributions
308*     to the result vector
309*---------------------------------------------------------------------*
310      IF      (LISTC(1:3).EQ.'R1 ' .OR. LISTC(1:3).EQ.'RE ' .OR.
311     &         LISTC(1:3).EQ.'RC '                              ) THEN
312
313C       -----------------------------------------
314C       calculate first-order triples amplitudes:
315C       -----------------------------------------
316        KDUM = KEND2
317        CALL CCSDT_T31_NODDY(WORK(KT3AM),LISTC,IDLSTC,FREQC,.FALSE.,
318     &                       .FALSE.,WORK(KINT1S),WORK(KINT2S),
319     &                       .FALSE.,WORK(KDUM),WORK(KDUM),
320     &                       .FALSE.,WORK(KDUM),WORK(KDUM),
321     &                               WORK(KINT1SC),WORK(KINT2SC),
322     &                       WORK(KLAMPC),WORK(KLAMHC),WORK(KFOCKC),
323     &                       WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0),
324     &                       WORK(KDUM),WORK(KFOCKD),
325     &                       WORK(KEND2),LWRK2)
326
327      ELSE IF (LISTC(1:3).EQ.'R2 ' .OR. LISTC(1:3).EQ.'ER1') THEN
328
329C       ------------------------------------------
330C       calculate second-order triples amplitudes:
331C       ------------------------------------------
332        CALL CCSDT_T32_NODDY(WORK(KT3AM),LISTC,IDLSTC,FREQC,
333     &                       WORK(KINT1S),WORK(KINT2S),
334     &                       WORK(KLAMP0),WORK(KLAMH0),WORK(KFOCK0),
335     &                       WORK(KFOCKD),WORK(KFIELDAO),WORK(KFIELD),
336     &                       WORK(KSCR1),WORK(KEND2),LWRK2)
337
338      ELSE
339
340        CALL QUIT('Unknown list '//LISTC//' in CCSDT_AA_NODDY.')
341
342      END IF
343
344      CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KT3AM),1)
345
346      IF (LOCDBG) THEN
347         WRITE(LUPRI,*) 'CCSDT_AAMAT_NODDY> T^C:'
348         WRITE(LUPRI,*) 'CCSDT_AAMAT_NODDY> LISTC,IDLSTC:',LISTC,IDLSTC
349         CALL PRINT_PT3_NODDY(WORK(KT3AM))
350       END IF
351
352C     -----------------------------------------------
353C     add contribution to doubles: <mu_2|[B,T3^C]|HF>
354C     -----------------------------------------------
355      CALL DZERO(WORK(KOMEGA2),NT1AMX*NT1AMX)
356      CALL CCSDT_XKSI2_2(WORK(KOMEGA2),WORK(KFOCKB),WORK(KT3AM))
357
358      DO I = 1,NT1AMX
359         DO J = 1,I
360            IJ = NT1AMX*(I-1) + J
361            NIJ = INDEX(I,J)
362            OMEGA2(NIJ) = OMEGA2(NIJ) + WORK(KOMEGA2+IJ-1)
363            IF (.NOT. OEFF_INIT)
364     &        OMEGA2EFF(NIJ) = OMEGA2EFF(NIJ) + WORK(KOMEGA2+IJ-1)
365         END DO
366      END DO
367
368C     -------------------------------------------
369C     contribution to triples: <mu_3|[B,T3^C]|HF>
370C     -------------------------------------------
371      CALL CCSDT_XKSI3_2(WORK(KOMEGA3),WORK(KFOCKB),WORK(KT3AM))
372
373      if (print_t3) then
374        write(lupri,*)'CCSDT_AA_NODDY> vector type:',listc
375        write(lupri,*)'CCSDT_AA_NODDY> idlst:',idlstc
376        write(lupri,*)'CCSDT_AA_NODDY> freq:',freqres
377        call ccsdt_clean_t3(work(komega3),nt1amx,nvirt,nrhft)
378        call print_pt3_noddy(work(komega3))
379      end if
380
381*---------------------------------------------------------------------*
382*     Now we split:
383*       for IOPTRES < 5 we compute the effective result vector
384*       for IOPTRES = 5 we compute the contractions Tbar^D A{B} T^C
385*---------------------------------------------------------------------*
386      IF (IOPTRES.GE.1 .AND. IOPTRES.LE.4) THEN
387
388         IF (OEFF_INIT) THEN
389           CALL DCOPY(NT1AMX,OMEGA1,1,OMEGA1EFF,1)
390           CALL DCOPY(NT2AMX,OMEGA2,1,OMEGA2EFF,1)
391         END IF
392
393         CALL CC_RHPART_NODDY(OMEGA1EFF,OMEGA2EFF,WORK(KOMEGA3),FREQRES,
394     &                        WORK(KFOCKD),WORK(KFOCK0),WORK(KFIELD),
395     &                        WORK(KXIAJB),WORK(KINT1T),WORK(KINT2T),
396     &                        WORK(KEND1),LWRK1)
397
398      ELSE IF (IOPTRES.EQ.5) THEN
399
400        SIGN = -1.0D0
401        CALL CCDOTRSP_NODDY(WORK(KOMEGA1),WORK(KOMEGA2),
402     &                      WORK(KOMEGA3),SIGN,
403     &                      ITRAN,LISTDP,IDOTS,DOTPROD,MXVEC,
404     &                      WORK(KLAMP0),WORK(KLAMH0),
405     &                      WORK(KFOCK0),WORK(KFOCKD),
406     &                      WORK(KXIAJB),WORK(KYIAJB),
407     &                      WORK(KINT1T),WORK(KINT2T),
408     &                      WORK(KINT1S),WORK(KINT2S),
409     &                      'CCSDT_AAMAT_NODDY',LOCDBG,LOCDBG,.FALSE.,
410     &                      WORK(KEND1),LWRK1)
411
412      ELSE
413        CALL QUIT('Illegal value for IOPTRES IN CCSDT_AAMAT_NODDY')
414      END IF
415
416      CALL QEXIT('CCSDT_AAMAT_NODDY')
417      RETURN
418      END
419
420*---------------------------------------------------------------------*
421*              END OF SUBROUTINE CCSDT_AAMAT_NODDY                    *
422*---------------------------------------------------------------------*
423