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
19C    /* Deck cc3_aabmat_doub */
20      SUBROUTINE CC3_AABMAT_DOUB(OMEGA2,LISTB,IDLSTB,LISTC,
21     *                     IDLSTC,WORK,LWORK)
22
23******************************************************************
24*
25* Calculate the triples correction to the doubles part of the
26* right-hand side of second-order amplitudes equation:
27*
28* omega2 = omega2 + <mu2|[X,T3Y]|HF> + <mu2|[[H0,T1X],T3Y]|HF>
29*
30* T3Y calculated with W intermmediate.
31*
32*
33*
34* F. Pawlowski, P. Jorgensen, Aarhus, Spring 2003.
35*
36* (based on CCSDT_FBMAT routine)
37*
38******************************************************************
39*
40      IMPLICIT NONE
41C
42      INTEGER LWORK
43      INTEGER ISINT1, ISINT2, ISYMT2, ISYMT3, ISYMWB
44      INTEGER LU3SRTR, LUCKJDR, LUDELDR, LUDKBCR
45      INTEGER LUCKJDR2, LUDKBCR2
46      INTEGER KT2TP, KRBJIA, KFOCKD, KEND0, LWRK0, IOPT
47      INTEGER IDLSTB, ISYMB, KT1B, KT2B
48      INTEGER ISINTR1B, ISINTR2B, KLAMDP, KLAMDH, KEND1, LWRK1
49      INTEGER KTROCY, KTROC1Y, KTROC0, KTROC02, KTROC0X
50      INTEGER KINTOC, KEND2, LWRK2
51      INTEGER IOFF, ISYMD, ISCKB1, ISCKB2, ISCKB2X
52      INTEGER KTRVIY, KTRVI1Y, KTRVI2, KTRVI3, KTRVI0, KTRVI0X
53      INTEGER KEND3, LWRK3, KEND4, LWRK4, KINTVI, LENSQ, LENSQWB
54      INTEGER ISYALJ, ISYALJ2, ISYMBD, ISCKIJ, ISCKD2, ISCKD2X
55      INTEGER ISWMATB, KSMAT, KSMAT3, KUMAT, KUMAT3, KDIAG, KDIAGWB
56      INTEGER KWMAT, KINDSQ, KINDSQWB, KINDEX, KINDEX2, KTMAT
57      INTEGER KTRVI8, KTRVI8X, KTRVI9, KTRVI10
58      INTEGER LUCKJD, LUTOC, LU3VI, LU3VI2, LUDKBC, LUDELD
59      INTEGER KSAVE, KEND00, LWRK00
60      INTEGER IDLSTC, ISYMRB, ISYMRC, ISINT1C, ISYRES
61      INTEGER KLAMPB, KLAMHB, KFCKB, LUFCK
62      INTEGER KLAMPC, KLAMHC, KT1C, ISYFCKB, ISYFCKC, KFCKC
63      INTEGER IRREP, IERR, KXIAJB, IOPTTCME, ISYMK, KOFF1, KOFF2
64      INTEGER ISYMC, IRREP2, ISCKB1C, ISYCKB
65      INTEGER ISINT1B
66      INTEGER KTROCX,KTROC1X
67      INTEGER ISCKB1B
68      INTEGER KTRVIX,KTRVI1X
69      INTEGER ISINTR2C,ISINTR1C,ISCKB2Y,ISCKD2Y
70      INTEGER KTRVI0Y,KTRVI8Y,KTROC0Y
71      INTEGER ISYMWC,ISWMATC
72      INTEGER KT2C,KFCKBB,KFCKCC,KOFF3,MAXBC,KT3MAT,KDIAGWC,KINDSQWC
73      INTEGER LENSQWC
74      INTEGER ILSTSYM
75c
76      INTEGER kx3am
77c
78#if defined (SYS_CRAY)
79      REAL OMEGA2(*)
80      REAL WORK(LWORK)
81      REAL FREQB, FREQC
82      REAL XNORMVAL, DDOT
83#else
84      DOUBLE PRECISION OMEGA2(*)
85      DOUBLE PRECISION WORK(LWORK)
86      DOUBLE PRECISION FREQB, FREQC
87      DOUBLE PRECISION XNORMVAL, DDOT
88#endif
89      CHARACTER*6 FNCKJD, FNDELD, FN3VI
90      CHARACTER*8 FNTOC,FN3VI2
91      CHARACTER*4 FNDKBC
92C
93      PARAMETER (FNCKJD ='CKJDEL' , FNTOC   = 'CCSDT_OC' )
94      PARAMETER (FNDELD = 'CKDELD', FNDKBC  = 'DKBC' )
95      PARAMETER (FN3VI  = 'CC3_VI', FN3VI2  = 'CC3_VI12')
96C
97      CHARACTER*12 FN3SRTR, FNCKJDR, FNDELDR, FNDKBCR, FNDKBCR4
98      CHARACTER*13 FNCKJDR2, FNDKBCR2
99      CHARACTER*10 MODEL
100      CHARACTER*8 LABELB, LABELC
101      CHARACTER*3 LISTB, LISTC
102      CHARACTER CDUMMY*1
103C
104#include "priunit.h"
105#include "dummy.h"
106#include "iratdef.h"
107#include "ccsdsym.h"
108#include "inftap.h"
109#include "ccinftap.h"
110#include "ccorb.h"
111#include "ccsdinp.h"
112#include "ccr1rsp.h"
113#include "ccexci.h"
114C
115      PARAMETER(FN3SRTR  = 'CCSDT_FBMAT1', FNCKJDR  = 'CCSDT_FBMAT2',
116     *          FNDELDR  = 'CCSDT_FBMAT3', FNDKBCR  = 'CCSDT_FBMAT4',
117     *          FNDKBCR4 = 'CCSDT_FBMAT5')
118C
119      PARAMETER(FNCKJDR2  = 'CCSDT_FBMAT22',FNDKBCR2  = 'CCSDT_FBMAT42')
120C
121      CALL QENTER('AABDB')
122C
123C     Lists check
124C
125      IF ( (LISTB(1:3).EQ.'R1 ') .AND. (LISTC(1:3).EQ.'R1 ') ) THEN
126         CONTINUE
127      ELSE IF ( (LISTB(1:3).EQ.'R1 ') .AND. (LISTC(1:3).EQ.'RE ') ) THEN
128         CONTINUE
129      ELSE
130         WRITE(LUPRI,*)'LISTB = ', LISTB
131         WRITE(LUPRI,*)'LISTC = ', LISTC
132         WRITE(LUPRI,*)'Only implemented when  '
133     *   //'LISTB = R1 and LISTC = R1 or RE '
134         CALL QUIT('Case not implemented in CC3_AABMAT_DOUB')
135      END IF
136
137c
138c     write(lupri,*)'OMEGA2 at the beginning of CC3_AABMAT_DOUB '
139c     call outpak(OMEGA2,NT1AM(1),1,lupri)
140
141C
142C--------------------------
143C     Save and set flags.
144C--------------------------
145C
146      IPRCC   = IPRINT
147      ISINT1  = 1
148      ISINT2  = 1
149      ISYMT2  = 1
150      ISYMT3  = MULD2H(ISINT2,ISYMT2)
151C
152C--------------------------------
153C     Open files
154C--------------------------------
155C
156      LUCKJD   = -1
157      LUTOC    = -1
158      LU3VI    = -1
159      LU3VI2   = -1
160      LUDKBC   = -1
161      LUDELD   = -1
162C
163      CALL WOPEN2(LUCKJD,FNCKJD,64,0)
164      CALL WOPEN2(LUTOC,FNTOC,64,0)
165      CALL WOPEN2(LU3VI,FN3VI,64,0)
166      CALL WOPEN2(LU3VI2,FN3VI2,64,0)
167      CALL WOPEN2(LUDKBC,FNDKBC,64,0)
168      CALL WOPEN2(LUDELD,FNDELD,64,0)
169C
170C--------------------------------
171C     Open temporary files
172C--------------------------------
173C
174      LU3SRTR   = -1
175      LUCKJDR   = -1
176      LUDELDR   = -1
177      LUDKBCR   = -1
178      LUCKJDR2  = -1
179      LUDKBCR2  = -1
180C
181      CALL WOPEN2(LU3SRTR,FN3SRTR,64,0)
182      CALL WOPEN2(LUCKJDR,FNCKJDR,64,0)
183      CALL WOPEN2(LUDELDR,FNDELDR,64,0)
184      CALL WOPEN2(LUDKBCR,FNDKBCR,64,0)
185      CALL WOPEN2(LUCKJDR2,FNCKJDR2,64,0)
186      CALL WOPEN2(LUDKBCR2,FNDKBCR2,64,0)
187C
188C----------------------------------------------
189C     Calculate the zeroth order stuff once
190C----------------------------------------------
191C
192      KT2TP  = 1
193      KFOCKD = KT2TP  + NT2SQ(ISYMT2)
194      KLAMDP = KFOCKD + NORBTS
195      KLAMDH = KLAMDP + NLAMDT
196      KXIAJB = KLAMDH + NLAMDT
197      KEND00 = KXIAJB + NT2AM(ISINT1)
198      LWRK00 = LWORK  - KEND00
199C
200      IF (LWRK00 .LT. 0) THEN
201         CALL QUIT('Out of memory in CC3_AABMAT_DOUB (zeroth allo.')
202      ENDIF
203C
204C------------------------
205C     Construct L(ia,jb).
206C------------------------
207C
208      REWIND(LUIAJB)
209      CALL READI(LUIAJB,IRAT*NT2AM(ISINT1),WORK(KXIAJB))
210      IOPTTCME = 1
211      CALL CCSD_TCMEPK(WORK(KXIAJB),1.0D0,ISINT1,IOPTTCME)
212C
213      IF ( IPRINT .GT. 55) THEN
214         XNORMVAL = DDOT(NT2AM(ISINT1),WORK(KXIAJB),1,
215     *                WORK(KXIAJB),1)
216         WRITE(LUPRI,*) 'Norm of IAJB ',XNORMVAL
217      ENDIF
218C
219C----------------------------------------------------------------
220C     Read t1 and calculate the zero'th order Lambda matrices
221C----------------------------------------------------------------
222C
223      CALL GET_LAMBDA0(WORK(KLAMDP),WORK(KLAMDH),WORK(KEND00),LWRK00)
224C
225C-------------------------------------------
226C     Read in t2 , square it and reorder
227C-------------------------------------------
228C
229      IOPT = 2
230      CALL GET_T1_T2(IOPT,.FALSE.,DUMMY,WORK(KT2TP),'R0',0,ISYMT2,
231     *               WORK(KEND00),LWRK00)
232C
233      IF (IPRINT .GT. 55) THEN
234         XNORMVAL = DDOT(NT2SQ(ISYMT2),WORK(KT2TP),1,WORK(KT2TP),1)
235         WRITE(LUPRI,*) 'Norm of T2TP ',XNORMVAL
236      ENDIF
237C
238C--------------------------------------
239C     Read in orbital energies
240C--------------------------------------
241C
242      CALL GET_ORBEN(WORK(KFOCKD),WORK(KEND00),LWRK00)
243COMMENT COMMENT COMMENT
244COMMENT COMMENT COMMENT If we want to sum the T3 amplitudes
245COMMENT COMMENT COMMENT
246      if (.false.) then
247         kx3am  = kend00
248         kend00 = kx3am + nt1amx*nt1amx*nt1amx
249         call dzero(work(kx3am),nt1amx*nt1amx*nt1amx)
250         lwrk00 = lwork - kend00
251         if (lwrk00 .lt. 0) then
252            write(lupri,*) 'Memory available : ',lwork
253            write(lupri,*) 'Memory needed    : ',kend00
254            call quit('Insufficient space (T3) in CC3_AABMAT_DOUB')
255         END IF
256      endif
257COMMENT COMMENT COMMENT
258COMMENT COMMENT COMMENT
259COMMENT COMMENT COMMENT
260C
261C-------------------
262C     Lists handling
263C-------------------
264C
265
266
267      IF (LISTB(1:3).EQ.'R1 ') THEN
268         ISYMRB = ISYLRT(IDLSTB)
269         FREQB  = FRQLRT(IDLSTB)
270         LABELB = LRTLBL(IDLSTB)
271      ELSE
272         CALL QUIT('Unknown list for LISTB in CC3_AABMAT_DOUB.')
273      END IF
274C
275      IF (LISTC(1:3).EQ.'R1 ') THEN
276         ISYMRC = ISYLRT(IDLSTC)
277         FREQC  = FRQLRT(IDLSTC)
278         LABELC = LRTLBL(IDLSTC)
279      ELSE IF (LISTC(1:3).EQ.'RE ') THEN
280         ISYMRC = ILSTSYM(LISTC,IDLSTC)
281         FREQC  = EIGVAL(IDLSTC)
282         LABELC = '- none -'
283      ELSE
284         CALL QUIT('Unknown list for LISTC in CC3_AABMAT_DOUB.')
285      END IF
286C
287      ISYMWB   = MULD2H(ISYMT3,ISYMRB)
288      ISYMWC   = MULD2H(ISYMT3,ISYMRC)
289      ISINT1C = MULD2H(ISINT1,ISYMRC)
290      ISINT1B = MULD2H(ISINT1,ISYMRB)
291      ISYRES  = MULD2H(ISYMWB,ISINT1C)
292      !Symmetry check
293      IF ( ISYRES .NE. MULD2H(ISYMWC,ISINT1B) ) THEN
294         CALL QUIT('Symmetry mismatch in CC3_AABMAT_DOUB (ISYRES)')
295      END IF
296C
297      ISYFCKB = MULD2H(ISYMOP,ISYMRB)
298      ISYFCKC = MULD2H(ISYMOP,ISYMRC)
299C
300C----------------------------------------------
301C     Read T1^B and T2^B
302C     Read T1^C and T2^C
303C     Calculate (ck|de)-tilde(B) and (ck|lm)-tilde(B)
304C     Used to construct T3^B
305C
306C     Calculate (ck|de)-tilde(C) and (ck|lm)-tilde(C)
307C     Used to construct T3^C
308C----------------------------------------------
309C
310      KT2B   = KEND00
311      KEND0  = KT2B   + NT2SQ(ISYMRB)
312      LWRK0  = LWORK  - KEND0
313C
314      KT2C   = KEND0
315      KEND0  = KT2C   + NT2SQ(ISYMRC)
316C
317      KT1B   = KEND0
318      KEND0  = KT1B   + NT1AM(ISYMRB)
319      LWRK0  = LWORK  - KEND0
320C
321      KT1C   = KEND0
322      KEND0  = KT1C   + NT1AM(ISYMRC)
323      LWRK1  = LWORK  - KEND0
324C
325      IF (LWRK0 .LT. NT2SQ(ISYMRB)) THEN
326         CALL QUIT('Out of memory in CC3_AABMAT_DOUB (TOT_T3Y) ')
327      ENDIF
328C
329C     Readin T1B and T2B
330C
331      IOPT = 3
332      CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1B),WORK(KT2B),LISTB,IDLSTB,
333     *            ISYMRB,WORK(KEND0),LWRK0)
334C
335      IF (IPRINT .GT. 55) THEN
336         XNORMVAL = DDOT(NT1AM(ISYMRB),WORK(KT1B),1,WORK(KT1B),1)
337         WRITE(LUPRI,*) 'Norm of T1B  ',XNORMVAL
338         XNORMVAL = DDOT(NT2SQ(ISYMRB),WORK(KT2B),1,WORK(KT2B),1)
339         WRITE(LUPRI,*) 'Norm of T2B  ',XNORMVAL
340      ENDIF
341C
342C     Readin T1C and T2C
343C
344      IOPT = 3
345      CALL GET_T1_T2(IOPT,.TRUE.,WORK(KT1C),WORK(KT2C),LISTC,IDLSTC,
346     *            ISYMRC,WORK(KEND0),LWRK0)
347C
348      IF (IPRINT .GT. 55) THEN
349         XNORMVAL = DDOT(NT1AM(ISYMRC),WORK(KT1C),1,WORK(KT1C),1)
350         WRITE(LUPRI,*) 'Norm of T1C  ',XNORMVAL
351         XNORMVAL = DDOT(NT2SQ(ISYMRC),WORK(KT2C),1,WORK(KT2C),1)
352         WRITE(LUPRI,*) 'Norm of T2C  ',XNORMVAL
353      ENDIF
354C
355C----------------------------------------------------
356C     Integrals (ck|de)-tilde(B) and (ck|lm)-tilde(B)
357C----------------------------------------------------
358C
359      ISINTR1B = MULD2H(ISINT1,ISYMRB)
360      ISINTR2B = MULD2H(ISINT2,ISYMRB)
361C
362      CALL CC3_BARINT(WORK(KT1B),ISYMRB,WORK(KLAMDP),
363     *                WORK(KLAMDH),WORK(KEND0),LWRK0,
364     *                LU3SRTR,FN3SRTR,LUCKJDR,FNCKJDR)
365C
366      CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISINTR1B,LU3SRTR,FN3SRTR,
367     *               LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY,
368     *               IDUMMY,CDUMMY)
369C
370      CALL CC3_SINT(WORK(KLAMDH),WORK(KEND0),LWRK0,ISINTR1B,
371     *              LUDELDR,FNDELDR,LUDKBCR,FNDKBCR)
372C
373C--------------------------------
374C    Re-use some temporary files
375C--------------------------------
376C
377      CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE')
378      CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE')
379C
380      CALL WOPEN2(LU3SRTR,FN3SRTR,64,0)
381      CALL WOPEN2(LUDELDR,FNDELDR,64,0)
382C
383C----------------------------------------------------
384C     Integrals (ck|de)-tilde(C) and (ck|lm)-tilde(C)
385C----------------------------------------------------
386C
387      ISINTR1C = MULD2H(ISINT1,ISYMRC)
388      ISINTR2C = MULD2H(ISINT2,ISYMRC)
389C
390      CALL CC3_BARINT(WORK(KT1C),ISYMRC,WORK(KLAMDP),
391     *                WORK(KLAMDH),WORK(KEND0),LWRK0,
392     *                LU3SRTR,FN3SRTR,LUCKJDR2,FNCKJDR2)
393C
394      CALL CC3_SORT1(WORK(KEND0),LWRK0,2,ISINTR1C,LU3SRTR,FN3SRTR,
395     *               LUDELDR,FNDELDR,IDUMMY,CDUMMY,IDUMMY,CDUMMY,
396     *               IDUMMY,CDUMMY)
397C
398      CALL CC3_SINT(WORK(KLAMDH),WORK(KEND0),LWRK0,ISINTR1C,
399     *              LUDELDR,FNDELDR,LUDKBCR2,FNDKBCR2)
400C
401C---------------------------------------------------
402C     Regain work space and create lambdaC and lambdaB
403C     lambda-transformed (B and C) matrices
404C     and C and B property matices.
405C---------------------------------------------------
406C
407      ! Use KEND0 again, because neither KT1B nor KT1C are needed any more
408      KFCKB  = KEND0
409      KFCKC  = KFCKB  + N2BST(ISYFCKB)
410      KEND0  = KFCKC  + N2BST(ISYFCKC)
411C
412      KFCKBB = KEND0
413      KFCKCC = KFCKBB + N2BST(ISYFCKB)
414      KEND0  = KFCKCC + N2BST(ISYFCKC)
415C
416      KLAMPC = KEND0
417      KLAMHC = KLAMPC + NLAMDT
418      KEND0  = KLAMHC + NLAMDT
419      LWRK0  = LWORK  - KEND0
420C
421      KLAMPB = KEND0
422      KLAMHB = KLAMPB + NLAMDT
423      KEND0  = KLAMHB + NLAMDT
424      LWRK0  = LWORK  - KEND0
425C
426      IF (LWRK0 .LT. 0) THEN
427         WRITE(LUPRI,*) 'Memory available : ',LWORK
428         WRITE(LUPRI,*) 'Memory needed    : ',KEND0
429         CALL QUIT('Insufficient space in CC3_AABMAT_DOUB')
430      END IF
431C
432      CALL GET_LAMBDAX(WORK(KLAMPC),WORK(KLAMHC),LISTC,IDLSTC,ISYMRC,
433     *                 WORK(KLAMDP),WORK(KLAMDH),WORK(KEND0),LWRK0)
434C
435      CALL GET_LAMBDAX(WORK(KLAMPB),WORK(KLAMHB),LISTB,IDLSTB,ISYMRB,
436     *                 WORK(KLAMDP),WORK(KLAMDH),WORK(KEND0),LWRK0)
437C
438C---------------------------------------
439C     Calculate the F^B matrix
440C---------------------------------------
441C
442      CALL DZERO(WORK(KFCKB),N2BST(ISYFCKB))
443      IF (LISTB(1:3).EQ.'R1 ') THEN
444C
445        CALL GET_FOCKX(WORK(KFCKB),LABELB,IDLSTB,ISYFCKB,
446     *                 WORK(KLAMDP),1,WORK(KLAMDH),1,WORK(KEND0),LWRK0)
447C
448        IF (IPRINT .GT. 50) THEN
449           CALL AROUND( 'In CC3_AABMAT_DOUB: Fock^B MO matrix' )
450           CALL CC_PRFCKMO(WORK(KFCKB),ISYFCKB)
451        ENDIF
452C
453      END IF
454C
455C------------------------------------------------
456C     Calculate the F^C property-operator matrix
457C------------------------------------------------
458C
459      CALL DZERO(WORK(KFCKC),N2BST(ISYFCKC))
460C
461      IF (LISTC(1:3).EQ.'R1 ') THEN
462C
463        CALL GET_FOCKX(WORK(KFCKC),LABELC,IDLSTC,ISYFCKC,
464     *                 WORK(KLAMDP),1,WORK(KLAMDH),1,WORK(KEND0),LWRK0)
465C
466        IF (IPRINT .GT. 50) THEN
467           CALL AROUND( 'In CC3_AABMAT_DOUB: F^C property matrix' )
468           CALL CC_PRFCKMO(WORK(KFCKC),ISYFCKC)
469        ENDIF
470C
471      END IF
472C
473C-------------------------------------------------------------------------
474C        Calculate the F^C Fock matrix and add C property matrix
475C        (to be able to calculate <mu2|[C,T3B]|HF> contribution also here)
476C-------------------------------------------------------------------------
477C
478      CALL CC3LR_MFOCK(WORK(KEND0),WORK(KT1C),WORK(KXIAJB),ISYFCKC)
479C
480C     Put the F_{kc} part into correct F_{pq} and add the property matrix
481C
482      DO ISYMC = 1, NSYM
483         ISYMK = MULD2H(ISYFCKC,ISYMC)
484         DO K = 1, NRHF(ISYMK)
485            DO C = 1, NVIR(ISYMC)
486               KOFF1 = KFCKC - 1
487     *               + IFCVIR(ISYMK,ISYMC)
488     *               + NORB(ISYMK)*(C - 1)
489     *               + K
490               KOFF2 = KEND0 - 1
491     *               + IT1AM(ISYMC,ISYMK)
492     *               + NVIR(ISYMC)*(K - 1)
493     *               + C
494               KOFF3 = KFCKCC - 1
495     *               + IFCVIR(ISYMK,ISYMC)
496     *               + NORB(ISYMK)*(C - 1)
497     *               + K
498C
499               WORK(KOFF3) = WORK(KOFF1) + WORK(KOFF2)
500C
501            ENDDO
502         ENDDO
503      ENDDO
504C
505C-------------------------------------------------------------------------
506C        Calculate the F^B Fock matrix and add B property matrix
507C        (to be able to calculate <mu2|[B,T3C]|HF> contribution also here)
508C-------------------------------------------------------------------------
509C
510      CALL CC3LR_MFOCK(WORK(KEND0),WORK(KT1B),WORK(KXIAJB),ISYFCKB)
511C
512C     Put the F_{kc} part into correct F_{pq} and add the property matrix
513C
514      DO ISYMC = 1, NSYM
515         ISYMK = MULD2H(ISYFCKB,ISYMC)
516         DO K = 1, NRHF(ISYMK)
517            DO C = 1, NVIR(ISYMC)
518               KOFF1 = KFCKB - 1
519     *               + IFCVIR(ISYMK,ISYMC)
520     *               + NORB(ISYMK)*(C - 1)
521     *               + K
522               KOFF2 = KEND0 - 1
523     *               + IT1AM(ISYMC,ISYMK)
524     *               + NVIR(ISYMC)*(K - 1)
525     *               + C
526               KOFF3 = KFCKBB - 1
527     *               + IFCVIR(ISYMK,ISYMC)
528     *               + NORB(ISYMK)*(C - 1)
529     *               + K
530C
531               WORK(KOFF3) = WORK(KOFF1) + WORK(KOFF2)
532C
533            ENDDO
534         ENDDO
535      ENDDO
536C
537      IF (IPRINT .GT. 50) THEN
538         CALL AROUND( 'In CC3_AABMAT_DOUB: Fock^C MO matrix' )
539         CALL CC_PRFCKMO(WORK(KFCKC),ISYFCKC)
540      ENDIF
541C
542C--------------------------
543C     Read occupied integrals.
544C--------------------------
545C
546C     Memory allocation.
547C
548C------------------------------------------------------
549C     Work allocation
550C------------------------------------------------------
551C
552      KRBJIA = KEND0
553      KTROCY  = KRBJIA + NT2SQ(ISYRES)
554      KTROC1Y = KTROCY  + NTRAOC(ISINT1C)
555      KTROC0 = KTROC1Y + NTRAOC(ISINT1C)
556      KTROC02= KTROC0 + NTRAOC(ISINT2)
557      KEND1  = KTROC02+ NTRAOC(ISINT2)
558      LWRK1  = LWORK  - KEND1
559C
560      KTROCX  = KEND1
561      KTROC1X = KTROCX  + NTRAOC(ISINT1B)
562      KEND1   = KTROC1X + NTRAOC(ISINT1B)
563C
564      KTROC0X = KEND1
565      KEND1   = KTROC0X + NTRAOC(ISINTR2B)
566C
567      KTROC0Y = KEND1
568      KEND1   = KTROC0Y + NTRAOC(ISINTR2C)
569C
570      KINTOC = KEND1
571      KEND2  = KINTOC  + NTOTOC(ISINT1)
572      LWRK2  = LWORK  - KEND2
573C
574      IF (LWRK2 .LT. 0) THEN
575         WRITE(LUPRI,*) 'Memory available : ',LWORK
576         WRITE(LUPRI,*) 'Memory needed    : ',KEND2
577         CALL QUIT('Insufficient space in CC3_XI')
578      END IF
579C
580C-------------------------------------
581C     Initialise the result vectors
582C-------------------------------------
583C
584      CALL DZERO(WORK(KRBJIA),NT2SQ(ISYRES))
585C
586C------------------------
587C     Occupied integrals.
588C
589C     Read in integrals for SMAT etc.
590C-----------------------
591C
592      CALL INTOCC_T30(LUCKJD,FNCKJD,WORK(KLAMDP),ISINT2,WORK(KTROC0),
593     *                WORK(KTROC02),WORK(KEND2),LWRK2)
594C
595C------------------------------------------
596C     B transformed Occupied integrals.
597C-----------------------------------------
598C
599      CALL INTOCC_T3X(LUCKJDR,FNCKJDR,WORK(KLAMDP),ISINTR2B,
600     *                WORK(KTROC0X),WORK(KEND2),LWRK2)
601C
602C------------------------------------------
603C     C transformed Occupied integrals.
604C-----------------------------------------
605C
606      CALL INTOCC_T3X(LUCKJDR2,FNCKJDR2,WORK(KLAMDP),ISINTR2C,
607     *                WORK(KTROC0Y),WORK(KEND2),LWRK2)
608C
609C------------------------
610C     Occupied integrals.
611C
612C     Read in integrals for WMAT etc.
613C-----------------------
614C
615      IOFF = 1
616      IF (NTOTOC(ISYMOP) .GT. 0) THEN
617         CALL GETWA2(LUTOC,FNTOC,WORK(KINTOC),IOFF,NTOTOC(ISYMOP))
618      ENDIF
619C
620C----------------------------------
621C     Write out norms of Integrals.
622C----------------------------------
623C
624      IF (IPRINT .GT. 55) THEN
625         XNORMVAL = DDOT(NTOTOC(ISYMOP),WORK(KINTOC),1,
626     *                WORK(KINTOC),1)
627         WRITE(LUPRI,*) 'Norm of CCSDT_OC-INT ',XNORMVAL
628      ENDIF
629C
630C----------------------------------------------------------------------
631C     Transform (ia|j delta) integrals to (ia|j k)-tildeC and sort as (i,j,k,a)
632C----------------------------------------------------------------------
633C
634      CALL CCLR_TROCC(WORK(KINTOC),WORK(KTROCY),
635     *                 WORK(KLAMHC),ISYMRC,
636     *                 WORK(KEND2),LWRK2)
637C
638C----------------------------------------------------------------------
639C     sort (i,j,k,a) as (a,i,j,k)
640C----------------------------------------------------------------------
641C
642         CALL CCSDT_SRTOC2(WORK(KTROCY),WORK(KTROC1Y),ISINT1C,
643     *                     WORK(KEND2),LWRK2)
644C
645C----------------------------------------------------------------------
646C     Transform (ia|j delta) integrals to (ia|j k)-tildeB and sort as (i,j,k,a)
647C----------------------------------------------------------------------
648C
649      CALL CCLR_TROCC(WORK(KINTOC),WORK(KTROCX),
650     *                 WORK(KLAMHB),ISYMRB,
651     *                 WORK(KEND2),LWRK2)
652C
653C----------------------------------------------------------------------
654C     sort (i,j,k,a) as (a,i,j,k)
655C----------------------------------------------------------------------
656C
657         CALL CCSDT_SRTOC2(WORK(KTROCX),WORK(KTROC1X),ISINT1B,
658     *                     WORK(KEND2),LWRK2)
659C
660C----------------------------
661C     General loop structure.
662C----------------------------
663C
664      DO ISYMD = 1,NSYM
665C
666         ISYCKB  = MULD2H(ISYMOP,ISYMD)
667         ISCKB2  = MULD2H(ISINT2,ISYMD)
668         ISCKB2X = MULD2H(ISINTR2B,ISYMD)
669         ISCKB2Y = MULD2H(ISINTR2C,ISYMD)
670         ISCKB1  = MULD2H(ISINT1,ISYMD)
671         ISCKB1C = MULD2H(ISINT1C,ISYMD)
672         ISCKB1B = MULD2H(ISINT1B,ISYMD)
673C
674         IF (IPRINT .GT. 55) THEN
675C
676            WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISYCKB :',ISYCKB
677            WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISCKB2 :',ISCKB2
678            WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISCKB2X:',ISCKB2X
679            WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISCKB1 :',ISCKB1
680            WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISCKB1C:',ISCKB1C
681C
682         ENDIF
683
684C
685C--------------------------
686C        Memory allocation.
687C--------------------------
688C
689         KTRVIY  = KEND1
690         KTRVI1Y = KTRVIY  + NCKATR(ISCKB1C)
691         KTRVI3 = KTRVI1Y + NCKATR(ISCKB1C)
692         KTRVI2 = KTRVI3 + NCKATR(ISCKB2)
693         KEND2  = KTRVI2 + NCKATR(ISCKB2)
694         LWRK2  = LWORK  - KEND2
695C
696         KTRVIX  = KEND2
697         KTRVI1X = KTRVIX  + NCKATR(ISCKB1B)
698         KEND2 = KTRVI1X + NCKATR(ISCKB1B)
699C
700         KTRVI0  = KEND2
701         KEND3   = KTRVI0  + NCKATR(ISCKB2)
702         LWRK3   = LWORK  - KEND3
703C
704         KTRVI0X  = KEND3
705         KEND3    = KTRVI0X  + NCKATR(ISCKB2X)
706         LWRK3    = LWORK    - KEND3
707C
708         KTRVI0Y  = KEND3
709         KEND3    = KTRVI0Y  + NCKATR(ISCKB2Y)
710         LWRK3    = LWORK    - KEND3
711C
712         KINTVI = KEND3
713         KEND4  = KINTVI + NCKA(ISYCKB)
714         LWRK4  = LWORK  - KEND4
715C
716         IF (LWRK4 .LT. 0) THEN
717            WRITE(LUPRI,*) 'Memory available : ',LWORK
718            WRITE(LUPRI,*) 'Memory needed    : ',KEND4
719            CALL QUIT('Insufficient space in CC3_XI')
720         END IF
721C
722C---------------------
723C        Sum over D
724C---------------------
725C
726         DO D = 1,NVIR(ISYMD)
727C
728C------------------------------------------------------------
729C           Read and transform integrals used in contraction
730C           with W intermediate.
731C------------------------------------------------------------
732C
733            IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1
734            IF (NCKA(ISYCKB) .GT. 0) THEN
735               CALL GETWA2(LU3VI2,FN3VI2,WORK(KINTVI),IOFF,
736     &                     NCKA(ISYCKB))
737            ENDIF
738C
739            ! Transform for C operator
740            CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVIY),
741     *                       WORK(KLAMPC),ISYMRC,
742     *                       ISYMD,D,ISYMOP,WORK(KEND4),LWRK4)
743C
744            ! Transform for B operator
745            CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVIX),
746     *                       WORK(KLAMPB),ISYMRB,
747     *                       ISYMD,D,ISYMOP,WORK(KEND4),LWRK4)
748C
749            IOFF = ICKAD(ISYCKB,ISYMD) + NCKA(ISYCKB)*(D - 1) + 1
750            IF (NCKA(ISYCKB) .GT. 0) THEN
751               CALL GETWA2(LU3VI,FN3VI,WORK(KINTVI),IOFF,
752     &                     NCKA(ISYCKB))
753            ENDIF
754C
755            ! Transform for C operator
756            CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVI1Y),
757     *                       WORK(KLAMPC),ISYMRC,
758     *                       ISYMD,D,ISYMOP,WORK(KEND4),LWRK4)
759C
760            ! Transform for B operator
761            CALL CCLR_TRVIR(WORK(KINTVI),WORK(KTRVI1X),
762     *                       WORK(KLAMPB),ISYMRB,
763     *                       ISYMD,D,ISYMOP,WORK(KEND4),LWRK4)
764
765
766C
767            CALL INTVIR_T30_D(LUDKBC,FNDKBC,LUDELD,FNDELD,ISINT2,
768     *                        WORK(KTRVI0),WORK(KTRVI2),WORK(KTRVI3),
769     *                        WORK(KLAMDH),ISYMD,D,WORK(KEND4),LWRK4)
770C
771C-----------------------------------------------------------------------
772C           Read B transformed virtual integrals used for W for TOT_T3Y
773C-----------------------------------------------------------------------
774C
775            IOFF = ICKBD(ISCKB2X,ISYMD) + NCKATR(ISCKB2X)*(D - 1) + 1
776            IF (NCKATR(ISCKB2X) .GT. 0) THEN
777               CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI0X),IOFF,
778     &                     NCKATR(ISCKB2X))
779            ENDIF
780C
781C-----------------------------------------------------------------------
782C           Read C transformed virtual integrals used for W for TOT_T3Y
783C-----------------------------------------------------------------------
784C
785            IOFF = ICKBD(ISCKB2Y,ISYMD) + NCKATR(ISCKB2Y)*(D - 1) + 1
786            IF (NCKATR(ISCKB2Y) .GT. 0) THEN
787               CALL GETWA2(LUDKBCR2,FNDKBCR2,WORK(KTRVI0Y),IOFF,
788     &                     NCKATR(ISCKB2Y))
789            ENDIF
790C
791C-----------------------------------------------------
792C           Sum over ISYMB
793C-----------------------------------------------------
794C
795            DO ISYMB = 1,NSYM
796C
797               ISYALJ  = MULD2H(ISYMB,ISYMT2)
798               ISYALJ2 = MULD2H(ISYMD,ISYMT2)
799               ISYMBD  = MULD2H(ISYMB,ISYMD)
800               ISCKIJ  = MULD2H(ISYMBD,ISYMT3)
801               ISCKD2  = MULD2H(ISINT2,ISYMB)
802               ISCKD2X = MULD2H(ISINTR2B,ISYMB)
803               ISCKD2Y = MULD2H(ISINTR2C,ISYMB)
804               ISWMATB  = MULD2H(ISYMRB,ISCKIJ)
805               ISWMATC  = MULD2H(ISYMRC,ISCKIJ)
806C
807               IF (IPRINT .GT. 55) THEN
808C
809                  WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISYMD  :',ISYMD
810                  WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISYMB  :',ISYMB
811                  WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISYALJ :',ISYALJ
812                  WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISYALJ2:',ISYALJ2
813                  WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISYMBD :',ISYMBD
814                  WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISCKIJ :',ISCKIJ
815                  WRITE(LUPRI,*) 'In CC3_AABMAT_DOUB: ISWMATB :',ISWMATB
816C
817               ENDIF
818C
819C
820               MAXBC   = MAX(NCKIJ(ISWMATB),NCKIJ(ISWMATC))
821C
822C              Can use kend3 since we do not need the integrals anymore.
823               KSMAT   = KEND3
824               KSMAT3  = KSMAT   + NCKIJ(ISCKIJ)
825               KUMAT   = KSMAT3  + NCKIJ(ISCKIJ)
826               KUMAT3  = KUMAT   + NCKIJ(ISCKIJ)
827               KDIAG   = KUMAT3  + NCKIJ(ISCKIJ)
828               KDIAGWB  = KDIAG   + NCKIJ(ISCKIJ)
829               KWMAT   = KDIAGWB  + NCKIJ(ISWMATB)
830               KINDSQWB = KWMAT   + MAXBC
831               KINDSQ  = KINDSQWB + (6*NCKIJ(ISWMATB) - 1)/IRAT + 1
832               KINDEX  = KINDSQ  + (6*NCKIJ(ISCKIJ) - 1)/IRAT + 1
833               KINDEX2 = KINDEX  + (NCKI(ISYALJ) - 1)/IRAT + 1
834               KTMAT   = KINDEX2 + (NCKI(ISYALJ2) - 1)/IRAT + 1
835               KT3MAT  = KTMAT   + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWMATB))
836               KTRVI8  = KT3MAT   + MAX(NCKIJ(ISCKIJ),NCKIJ(ISWMATC))
837               KTRVI9  = KTRVI8  + NCKATR(ISCKD2)
838               KTRVI10 = KTRVI9  + NCKATR(ISCKD2)
839               KEND4   = KTRVI10 + NCKATR(ISCKD2)
840               LWRK4   = LWORK   - KEND4
841C
842               KDIAGWC  = KEND4
843               KINDSQWC = KDIAGWC + NCKIJ(ISWMATC)
844               KEND4    = KINDSQWC + (6*NCKIJ(ISWMATC) - 1)/IRAT + 1
845C
846               KTRVI8X = KEND4
847               KEND4   = KTRVI8X + NCKATR(ISCKD2X)
848               LWRK4   = LWORK   - KEND4
849C
850               KTRVI8Y = KEND4
851               KEND4   = KTRVI8Y + NCKATR(ISCKD2Y)
852               LWRK4   = LWORK   - KEND4
853C
854               IF (LWRK4 .LT. 0) THEN
855                  WRITE(LUPRI,*) 'Memory available : ',LWORK
856                  WRITE(LUPRI,*) 'Memory needed    : ',KEND4
857                  CALL QUIT('Insufficient space in CC3_XI')
858               END IF
859C
860C---------------------------------------------
861C              Construct part of the diagonal.
862C---------------------------------------------
863C
864               CALL CC3_DIAG(WORK(KDIAG),WORK(KFOCKD),ISCKIJ)
865               CALL CC3_DIAG(WORK(KDIAGWB),WORK(KFOCKD),ISWMATB)
866               CALL CC3_DIAG(WORK(KDIAGWC),WORK(KFOCKD),ISWMATC)
867C
868               IF (IPRINT .GT. 55) THEN
869                  XNORMVAL = DDOT(NCKIJ(ISCKIJ),WORK(KDIAG),1,
870     *                    WORK(KDIAG),1)
871                  WRITE(LUPRI,*) 'Norm of DIA  ',XNORMVAL
872               ENDIF
873C
874C-------------------------------------
875C              Construct index arrays.
876C-------------------------------------
877C
878               LENSQ = NCKIJ(ISCKIJ)
879               CALL CC3_INDSQ(WORK(KINDSQ),LENSQ,ISCKIJ)
880               CALL CC3_INDEX(WORK(KINDEX),ISYALJ)
881               CALL CC3_INDEX(WORK(KINDEX2),ISYALJ2)
882               LENSQWB = NCKIJ(ISWMATB)
883               CALL CC3_INDSQ(WORK(KINDSQWB),LENSQWB,ISWMATB)
884               LENSQWC = NCKIJ(ISWMATC)
885               CALL CC3_INDSQ(WORK(KINDSQWC),LENSQWC,ISWMATC)
886C
887               DO B = 1,NVIR(ISYMB)
888C
889C------------------------------------
890C                 Integrals for t30
891C------------------------------------
892C
893                  CALL INTVIR_T30_D(LUDKBC,FNDKBC,LUDELD,FNDELD,
894     *                              ISINT2,WORK(KTRVI8),WORK(KTRVI9),
895     *                              WORK(KTRVI10),WORK(KLAMDH),
896     *                              ISYMB,B,WORK(KEND4),LWRK4)
897C
898C--------------------------------------------------------------------
899C                 Read B transformed virtual integrals used in W
900C--------------------------------------------------------------------
901C
902                  IOFF = ICKBD(ISCKD2X,ISYMB) +
903     &                   NCKATR(ISCKD2X)*(B - 1) + 1
904                  IF (NCKATR(ISCKD2X) .GT. 0) THEN
905                     CALL GETWA2(LUDKBCR,FNDKBCR,WORK(KTRVI8X),IOFF,
906     &                           NCKATR(ISCKD2X))
907                  ENDIF
908C
909C--------------------------------------------------------------------
910C                 Read C transformed virtual integrals used in W
911C--------------------------------------------------------------------
912C
913                  IOFF = ICKBD(ISCKD2Y,ISYMB) +
914     &                   NCKATR(ISCKD2Y)*(B - 1) + 1
915                  IF (NCKATR(ISCKD2Y) .GT. 0) THEN
916                     CALL GETWA2(LUDKBCR2,FNDKBCR2,WORK(KTRVI8Y),IOFF,
917     &                           NCKATR(ISCKD2Y))
918                  ENDIF
919C
920                  CALL GET_T30_BD(ISYMT3,ISINT2,WORK(KT2TP),ISYMT2,
921     *                            WORK(KTMAT),WORK(KFOCKD),
922     *                            WORK(KDIAG),WORK(KINDSQ),LENSQ,
923     *                            WORK(KSMAT),WORK(KTRVI0),
924     *                            WORK(KTRVI2),WORK(KTROC0),
925     *                            WORK(KINDEX),WORK(KSMAT3),
926     *                            WORK(KTRVI8),WORK(KTRVI9),
927     *                            WORK(KINDEX2),WORK(KUMAT),
928     *                            WORK(KTRVI3),WORK(KTROC02),
929     *                            WORK(KUMAT3),WORK(KTRVI10),
930     *                            ISYMB,B,ISYMD,D,
931     *                            ISCKIJ,WORK(KEND4),LWRK4)
932
933c                 call sum_pt3(work(KTMAT),isymb,b,isymd,d,
934c    *                        isckij,work(kx3am),4)
935
936C
937
938                  ! Copy T3^0 from KTMAT to KT3MAT
939                  ! KTMAT is used as help array and overwritten.
940                  ! But we still need T3^0  for WC (and there we use KT3MAT)
941                  CALL DCOPY(NCKIJ(ISCKIJ),WORK(KTMAT),1,
942     *                       WORK(KT3MAT),1)
943C
944
945C
946C-----------------------------------------------------------------------
947C                 Based on T30 BD amplitudes calculate WB BD intermediates
948C-----------------------------------------------------------------------
949C
950C
951                  CALL DZERO(WORK(KWMAT),MAXBC)
952C
953                  CALL GET_T3X_BD(ISYMRB,WORK(KWMAT),ISWMATB,
954     *                           WORK(KTMAT),ISCKIJ,WORK(KFCKB),
955     *                           ISYMRB,WORK(KT2TP),ISYMT2,
956     *                           WORK(KT2B),ISYMRB,
957     *                           WORK(KTRVI0),WORK(KTRVI8),
958     *                           WORK(KTROC0),ISINT2,
959     *                           WORK(KTRVI0X),WORK(KTRVI8X),
960     *                           WORK(KTROC0X),ISINTR2B,
961     *                           FREQB,WORK(KDIAGWB),WORK(KFOCKD),
962     *                           WORK(KINDSQWB),LENSQWB,
963     *                           B,ISYMB,D,ISYMD,WORK(KEND4),LWRK4)
964c
965
966c                 call sum_pt3(work(kwmat),isymb,b,isymd,d,
967c    *                        isckij,work(kx3am),4)
968c
969
970C
971C---------------------------------------------------------
972C                 Calculate the  term <mu2|[H,W^BD(3)]|HF>
973C                 ( Fock matrix cont )
974C                 added in OMEGA2
975C---------------------------------------------------------
976C
977                  CALL CC3_CY2F(OMEGA2,ISYRES,WORK(KWMAT),
978     *                          ISWMATB,WORK(KTMAT),WORK(KFCKCC),
979     *                          ISYFCKC,WORK(KINDSQWB),LENSQWB,
980     *                          WORK(KEND4),LWRK4,
981     *                          ISYMB,B,ISYMD,D)
982C
983C---------------------------------------------------------
984C                 Calculate the  term <mu2|[H,W^BD(3)]|HF>
985C                 ( Occupied  cont )
986C                 added in OMEGA2
987C---------------------------------------------------------
988C
989                  CALL CC3_CY2O(OMEGA2,ISYRES,WORK(KWMAT),
990     *                          ISWMATB,WORK(KTMAT),WORK(KTROCY),
991     *                          WORK(KTROC1Y),ISINT1C,
992     *                          WORK(KEND4),LWRK4,
993     *                          WORK(KINDSQWB),LENSQWB,
994     *                          ISYMB,B,ISYMD,D)
995C
996C---------------------------------------------------------
997C                 Calculate the  term <mu2|[H,W^BD(3)]|HF>
998C                 ( Virtual cont )
999C                 added in OMEGA2
1000C---------------------------------------------------------
1001C
1002                  CALL CC3_CY2V(OMEGA2,ISYRES,WORK(KRBJIA),
1003     *                          WORK(KWMAT),ISWMATB,WORK(KTMAT),
1004     *                          WORK(KTRVIY),WORK(KTRVI1Y),ISINT1C,
1005     *                          WORK(KEND4),LWRK4,WORK(KINDSQWB),
1006     *                          LENSQWB,ISYMB,B,ISYMD,D)
1007C
1008C-----------------------------------------------------------------------
1009C                 Based on T30 BD amplitudes calculate WC BD intermediates
1010C-----------------------------------------------------------------------
1011C
1012C
1013                  CALL DZERO(WORK(KWMAT),MAXBC)
1014C
1015                  IF (LISTC(1:3).EQ.'R1 ') THEN
1016                    CALL GET_T3X_BD(ISYMRC,WORK(KWMAT),ISWMATC,
1017     *                             WORK(KT3MAT),ISCKIJ,WORK(KFCKC),
1018     *                             ISYMRC,WORK(KT2TP),ISYMT2,
1019     *                             WORK(KT2C),ISYMRC,
1020     *                             WORK(KTRVI0),WORK(KTRVI8),
1021     *                             WORK(KTROC0),ISINT2,
1022     *                             WORK(KTRVI0Y),WORK(KTRVI8Y),
1023     *                             WORK(KTROC0Y),ISINTR2C,
1024     *                             FREQC,WORK(KDIAGWC),WORK(KFOCKD),
1025     *                             WORK(KINDSQWC),LENSQWC,
1026     *                             B,ISYMB,D,ISYMD,WORK(KEND4),LWRK4)
1027c
1028c
1029c                   call sum_pt3(work(kwmat),isymb,b,isymd,d,
1030c    *                          isckij,work(kx3am),4)
1031c
1032                  END IF
1033C
1034                  IF (LISTC(1:3).EQ.'RE ') THEN
1035                    CALL GET_E3_BD(ISYMRC,WORK(KWMAT),ISWMATC,
1036     *                             WORK(KT3MAT),
1037     *                             WORK(KT2TP),ISYMT2,
1038     *                             WORK(KT2C),ISYMRC,
1039     *                             WORK(KTRVI0),WORK(KTRVI8),
1040     *                             WORK(KTROC0),ISINT2,
1041     *                             WORK(KTRVI0Y),WORK(KTRVI8Y),
1042     *                             WORK(KTROC0Y),ISINTR2C,
1043     *                             FREQC,WORK(KDIAGWC),WORK(KFOCKD),
1044     *                             WORK(KINDSQWC),LENSQWC,
1045     *                             B,ISYMB,D,ISYMD,WORK(KEND4),LWRK4)
1046c
1047c
1048c                   call sum_pt3(work(kwmat),isymb,b,isymd,d,
1049c    *                          isckij,work(kx3am),4)
1050c
1051                  END IF
1052
1053C
1054C---------------------------------------------------------
1055C                 Calculate the  term <mu2|[H,W^BD(3)]|HF>
1056C                 ( Fock matrix cont )
1057C                 added in OMEGA2
1058C---------------------------------------------------------
1059C
1060                  CALL CC3_CY2F(OMEGA2,ISYRES,WORK(KWMAT),
1061     *                          ISWMATC,WORK(KTMAT),WORK(KFCKBB),
1062     *                          ISYFCKB,WORK(KINDSQWC),LENSQWC,
1063     *                          WORK(KEND4),LWRK4,
1064     *                          ISYMB,B,ISYMD,D)
1065c
1066c---------------------------------------------------------
1067c                 Calculate the  term <mu2|[H,W^BD(3)]|HF>
1068c                 ( Occupied  cont )
1069c                 added in OMEGA2
1070c---------------------------------------------------------
1071c
1072                  CALL CC3_CY2O(OMEGA2,ISYRES,WORK(KWMAT),
1073     *                          ISWMATC,WORK(KTMAT),WORK(KTROCX),
1074     *                          WORK(KTROC1X),ISINT1B,
1075     *                          WORK(KEND4),LWRK4,
1076     *                          WORK(KINDSQWC),LENSQWC,
1077     *                          ISYMB,B,ISYMD,D)
1078c
1079c---------------------------------------------------------
1080c                 Calculate the  term <mu2|[H,W^BD(3)]|HF>
1081c                 ( Virtual cont )
1082c                 added in OMEGA2
1083c---------------------------------------------------------
1084
1085                  CALL CC3_CY2V(OMEGA2,ISYRES,WORK(KRBJIA),
1086     *                          WORK(KWMAT),ISWMATC,WORK(KTMAT),
1087     *                          WORK(KTRVIX),WORK(KTRVI1X),ISINT1B,
1088     *                          WORK(KEND4),LWRK4,WORK(KINDSQWC),
1089     *                          LENSQWC,ISYMB,B,ISYMD,D)
1090C
1091               END DO  ! B
1092            END DO     ! ISYMB
1093C
1094         END DO        ! D
1095      END DO           ! ISYMD
1096C
1097
1098C
1099C------------------------------------------------------
1100C     Accumulate RBJIA from <mu2|[H,W^BD(3)]|HF> ( Vccupied  cont )
1101C     in OMEGA2
1102C------------------------------------------------------
1103C
1104      CALL CC3_RBJIA(OMEGA2,ISYRES,WORK(KRBJIA))
1105C
1106      IF (IPRINT .GT. 55) THEN
1107         XNORMVAL = DDOT(NT2AM(ISYRES),OMEGA2,1,OMEGA2,1)
1108         WRITE(LUPRI,*) 'Norm of final OMEGA2 ',XNORMVAL
1109      ENDIF
1110C
1111C--------------------------------
1112C     Close files
1113C--------------------------------
1114C
1115      CALL WCLOSE2(LUCKJD,FNCKJD,'KEEP')
1116      CALL WCLOSE2(LUTOC,FNTOC,'KEEP')
1117      CALL WCLOSE2(LU3VI,FN3VI,'KEEP')
1118      CALL WCLOSE2(LU3VI2,FN3VI2,'KEEP')
1119      CALL WCLOSE2(LUDKBC,FNDKBC,'KEEP')
1120      CALL WCLOSE2(LUDELD,FNDELD,'KEEP')
1121C
1122C--------------------------------
1123C     Close files for "response"
1124C--------------------------------
1125C
1126      CALL WCLOSE2(LU3SRTR,FN3SRTR,'DELETE')
1127      CALL WCLOSE2(LUCKJDR,FNCKJDR,'DELETE')
1128      CALL WCLOSE2(LUDELDR,FNDELDR,'DELETE')
1129      CALL WCLOSE2(LUDKBCR,FNDKBCR,'DELETE')
1130      CALL WCLOSE2(LUCKJDR2,FNCKJDR2,'DELETE')
1131      CALL WCLOSE2(LUDKBCR2,FNDKBCR2,'DELETE')
1132C
1133C-------------
1134C     End
1135C-------------
1136C
1137      CALL QEXIT('AABDB')
1138C
1139      RETURN
1140      END
1141