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 ccsdt_tran1_r */
20      SUBROUTINE CCSDT_TRAN1_R(XINT1,XINT2,XLAMDP0,XLAMDH0,
21     &                         XLAMDP1,XLAMDH1,AOINT,IDEL)
22C
23C     XINT1 = XINT1 + (K^p0 C^h0|B^p1 D^h0)
24C     XINT2 = XINT2 + (K^p0 C^h0|L^p0 J^h1)
25C
26C
27#include "implicit.h"
28#include "priunit.h"
29#include "inforb.h"
30#include "ccsdinp.h"
31#include "ccsdsym.h"
32      DIMENSION XINT1(NT1AMX,NVIRT,NVIRT), XINT2(NT1AMX,NRHFT,NRHFT)
33      DIMENSION XLAMDP0(NBAST,NORBT), XLAMDH0(NBAST,NORBT)
34      DIMENSION XLAMDP1(NBAST,NORBT), XLAMDH1(NBAST,NORBT)
35      DIMENSION AOINT(NNBAST,NBAST)
36C
37      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
38C
39      DO 100 IGAM = 1,NBAST
40         DO 110 IBET = 1,NBAST
41            DO 120 IALP = 1,NBAST
42               NAB = INDEX(IALP,IBET)
43C
44               if (aoint(nab,IGAM) .eq. 0.0d0) goto 120
45               DO 200 D = 1,NVIRT
46                  DO 210 B = 1,NVIRT
47                     DO 220 K = 1,NRHFT
48                        DO 230 C = 1,NVIRT
49C
50                           NCK = NVIRT*(K-1) + C
51C
52                           XINT1(NCK,B,D)=XINT1(NCK,B,D)+AOINT(NAB,IGAM)
53     &                     *XLAMDP0(IBET,K)      *XLAMDH0(IALP,NRHFT+C)
54     &                     *XLAMDP1(IGAM,NRHFT+B)*XLAMDH0(IDEL,NRHFT+D)
55C
56  230                   CONTINUE
57  220                CONTINUE
58  210             CONTINUE
59  200          CONTINUE
60C
61               DO 300 J = 1,NRHFT
62                  DO 310 L = 1,NRHFT
63                     DO 320 K = 1,NRHFT
64                        DO 330 C = 1,NVIRT
65C
66                           NCK = NVIRT*(K-1) + C
67C
68                           XINT2(NCK,L,J)=XINT2(NCK,L,J)+AOINT(NAB,IGAM)
69     &                      * XLAMDP0(IBET,K) * XLAMDH0(IALP,NRHFT+C)
70     &                      * XLAMDP0(IGAM,L) * XLAMDH1(IDEL,J)
71C
72  330                   CONTINUE
73  320                CONTINUE
74  310             CONTINUE
75  300          CONTINUE
76C
77  120       CONTINUE
78  110    CONTINUE
79  100 CONTINUE
80C
81      RETURN
82      END
83C  /* Deck ccsdt_tran3_r */
84      SUBROUTINE CCSDT_TRAN3_R(XINT1,XINT2,XLAMDP0,XLAMDH0,
85     &                         XLAMDP1,XLAMDH1,XLAMDP2,XLAMDH2,
86     &                         AOINT,IDEL)
87C
88C     XINT1 = XINT1 + (C^p1 K^h1|B^p2 D^h0)
89C     XINT2 = XINT2 + (C^p1 K^h1|L^p0 J^h2)
90C
91#include "implicit.h"
92#include "priunit.h"
93#include "inforb.h"
94#include "ccsdinp.h"
95#include "ccsdsym.h"
96      DIMENSION XINT1(NT1AMX,NVIRT,NVIRT), XINT2(NT1AMX,NRHFT,NRHFT)
97      DIMENSION XLAMDP0(NBAST,NORBT), XLAMDH0(NBAST,NORBT)
98      DIMENSION XLAMDP1(NBAST,NORBT), XLAMDH1(NBAST,NORBT)
99      DIMENSION XLAMDP2(NBAST,NORBT), XLAMDH2(NBAST,NORBT)
100      DIMENSION AOINT(NNBAST,NBAST)
101C
102      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
103C
104      DO 100 G = 1,NBAST
105         DO 110 IB = 1,NBAST
106            DO 120 A = 1,NBAST
107               NAB = INDEX(A,IB)
108C
109               if (aoint(nab,g) .eq. 0.0d0) goto 120
110               DO 200 D = 1,NVIRT
111                  DO 210 B = 1,NVIRT
112                     DO 220 K = 1,NRHFT
113                        DO 230 C = 1,NVIRT
114C
115                           NCK = NVIRT*(K-1) + C
116C
117                           XINT1(NCK,B,D) = XINT1(NCK,B,D)+AOINT(NAB,G)
118     &                      * XLAMDP1(A,NRHFT+C) * XLAMDH1(IB,K)
119     &                      * XLAMDP2(G,NRHFT+B) * XLAMDH0(IDEL,NRHFT+D)
120C
121  230                   CONTINUE
122  220                CONTINUE
123  210             CONTINUE
124  200          CONTINUE
125C
126               DO 300 J = 1,NRHFT
127                  DO 310 L = 1,NRHFT
128                     DO 320 K = 1,NRHFT
129                        DO 330 C = 1,NVIRT
130C
131                           NCK = NVIRT*(K-1) + C
132C
133                           XINT2(NCK,L,J) = XINT2(NCK,L,J)+AOINT(NAB,G)
134     &                      * XLAMDP1(A,NRHFT+C) * XLAMDH1(IB,K)
135     &                      * XLAMDP0(G,L)       * XLAMDH2(IDEL,J)
136C
137  330                   CONTINUE
138  320                CONTINUE
139  310             CONTINUE
140  300          CONTINUE
141C
142  120       CONTINUE
143  110    CONTINUE
144  100 CONTINUE
145C
146      RETURN
147      END
148C  /* Deck ccfop_tran1 */
149      SUBROUTINE CCFOP_TRAN1_R(XINT1,XINT2,XINT3,XINT4,
150     &                         XLAMDP0,XLAMDH0,
151     &                         XLAMDP1,XLAMDH1,
152     &                         XLAMDP2,XLAMDH2,
153     &                         AOINT,IDEL)
154C
155C     XINT1 = (K^p0 C^h0|D^p1 L^h2)    "(O-0 V-0|V-1 O-2)"
156C     XINT2 = (K^p0 L^h1|C^p2 D^h0)    "(O-0 O-1|V-2 V-0)"
157C     XINT3 = (K^p0 L^h1|M^p0 N^h2)    "(O-0 O-1|O-0 O-2)"
158C     XINT4 = (C^p1 D^h0|E^p2 F^h0)    "(V-1 V-0|V-2 V-0)"
159C
160#include "implicit.h"
161#include "priunit.h"
162#include "inforb.h"
163#include "ccsdinp.h"
164#include "ccsdsym.h"
165      DIMENSION XINT1(NRHFT,NVIRT,NVIRT,NRHFT)
166      DIMENSION XINT2(NRHFT,NRHFT,NVIRT,NVIRT)
167      DIMENSION XINT3(NRHFT,NRHFT,NRHFT,NRHFT)
168      DIMENSION XINT4(NVIRT,NVIRT,NVIRT,NVIRT)
169      DIMENSION XLAMDP0(NBAST,NORBT), XLAMDH0(NBAST,NORBT)
170      DIMENSION XLAMDP1(NBAST,NORBT), XLAMDH1(NBAST,NORBT)
171      DIMENSION XLAMDP2(NBAST,NORBT), XLAMDH2(NBAST,NORBT)
172      DIMENSION AOINT(NNBAST,NBAST)
173C
174      LOGICAL LDEBUG
175C
176      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
177C
178      LDEBUG = .TRUE.
179C
180C----------------------------------------
181C     Calculate integrals :
182C----------------------------------------
183C
184      DO 100 G = 1,NBAST
185         DO 110 IB = 1,NBAST
186            DO 120 A = 1,NBAST
187               NAB = INDEX(A,IB)
188C
189               if (aoint(nab,g) .eq. 0.0d0) goto 120
190C
191               DO NC = 1,NVIRT
192                  DO ND = 1,NVIRT
193                     DO NK = 1,NRHFT
194                        DO NL = 1,NRHFT
195C
196                           XINT1(NK,NC,ND,NL) = XINT1(NK,NC,ND,NL)
197     &                       + AOINT(NAB,G) *
198     &                          XLAMDP0(A,NK) *
199     &                          XLAMDH0(IB,NRHFT+NC) *
200     &                          XLAMDP1(G,NRHFT+ND) *
201     &                          XLAMDH2(IDEL,NL)
202C
203                        ENDDO
204                     ENDDO
205                  ENDDO
206               ENDDO
207C
208               DO NC = 1,NVIRT
209                  DO ND = 1,NVIRT
210                     DO NK = 1,NRHFT
211                        DO NL = 1,NRHFT
212C
213                           XINT2(NK,NL,NC,ND) = XINT2(NK,NL,NC,ND)
214     &                       + AOINT(NAB,G) *
215     &                          XLAMDP0(A,NK) *
216     &                          XLAMDH1(IB,NL) *
217     &                          XLAMDP2(G,NRHFT+NC) *
218     &                          XLAMDH0(IDEL,NRHFT+ND)
219C
220                        ENDDO
221                     ENDDO
222                  ENDDO
223               ENDDO
224C
225               DO NK = 1,NRHFT
226                  DO NL = 1,NRHFT
227                     DO NM = 1,NRHFT
228                        DO NN = 1,NRHFT
229C
230                           XINT3(NK,NL,NM,NN) = XINT3(NK,NL,NM,NN)
231     &                       + AOINT(NAB,G) *
232     &                          XLAMDP0(A,NK) *
233     &                          XLAMDH1(IB,NL) *
234     &                          XLAMDP0(G,NM) *
235     &                          XLAMDH2(IDEL,NN)
236C
237                        ENDDO
238                     ENDDO
239                  ENDDO
240               ENDDO
241C
242               DO NC = 1,NVIRT
243                  DO ND = 1,NVIRT
244                     DO NE = 1,NVIRT
245                        DO NF = 1,NVIRT
246C
247                           XINT4(NC,ND,NE,NF) = XINT4(NC,ND,NE,NF)
248     &                       + AOINT(NAB,G) *
249     &                          XLAMDP1(A,NRHFT+NC) *
250     &                          XLAMDH0(IB,NRHFT+ND) *
251     &                          XLAMDP2(G,NRHFT+NE) *
252     &                          XLAMDH0(IDEL,NRHFT+NF)
253C
254                        ENDDO
255                     ENDDO
256                  ENDDO
257               ENDDO
258C
259  120       CONTINUE
260  110    CONTINUE
261  100 CONTINUE
262C
263      RETURN
264      END
265*=====================================================================*
266      SUBROUTINE CCSDT_INTS0_NODDY(DO_INTXY,XIAJB,YIAJB,
267     &                             DO_INTS0,XINT1S,XINT2S,
268     &                             DO_INTT0,XINT1T,XINT2T,
269     &                             XLAMDP,XLAMDH,WORK,LWORK)
270*---------------------------------------------------------------------*
271*
272*    Purpose: compute some standard integrals for CC3
273*
274*             if (do_intxy) compute
275*                 XIAJB  = 2(kc|ld) - (kd|lc)
276*                 YIAJB  =  (kc|ld)
277*
278*             if (do_ints0) compute
279*                 XINT1S = (ck|bd) and
280*                 XINT2S = (ck|lj)
281*
282*             if (do_intt0) compute
283*                 XINT1T = (kc|bd) and
284*                 XINT2T = (kc|lj)
285*
286* Written by Christof Haettig, November 2002, based on CCSDT_XI3_NODDY.
287*
288*=====================================================================*
289#if defined (IMPLICIT_NONE)
290      IMPLICIT NONE
291#else
292#  include "implicit.h"
293#endif
294#include "priunit.h"
295#include "ccsdinp.h"
296#include "maxorb.h"
297#include "ccsdsym.h"
298#include "ccfield.h"
299#include "ccorb.h"
300#include "dummy.h"
301
302      LOGICAL LOCDBG
303      PARAMETER (LOCDBG=.FALSE.)
304      INTEGER ISYM0
305      PARAMETER (ISYM0=1)
306
307      LOGICAL DO_INTXY, DO_INTS0, DO_INTT0
308      INTEGER LWORK
309      REAL*8  WORK(LWORK)
310      REAL*8  XINT1T(NT1AMX*NVIRT*NVIRT)
311      REAL*8  XINT2T(NRHFT*NRHFT*NT1AMX)
312      REAL*8  XINT1S(NT1AMX*NVIRT*NVIRT)
313      REAL*8  XINT2S(NRHFT*NRHFT*NT1AMX)
314      REAL*8  XIAJB(NT1AMX*NT1AMX)
315      REAL*8  YIAJB(NT1AMX*NT1AMX)
316      REAL*8  XLAMDP(*), XLAMDH(*)
317
318      INTEGER ISYMD, ILLL, IDEL, ISYDIS, KXINT, KEND1, LWRK1
319
320      CALL QENTER('CCSDT_INTS0_NODDY')
321
322      IF (DIRECT)
323     &  CALL QUIT('DIRECT NOT IMPLEMENTED IN CCSDT_INTS0_NODDY')
324
325*---------------------------------------------------------------------*
326*     Loop over distributions of integrals:
327*---------------------------------------------------------------------*
328      IF (DO_INTS0) THEN
329        CALL DZERO(XINT1S,NT1AMX*NVIRT*NVIRT)
330        CALL DZERO(XINT2S,NT1AMX*NRHFT*NRHFT)
331      END IF
332      IF (DO_INTT0) THEN
333         CALL DZERO(XINT1T,NT1AMX*NVIRT*NVIRT)
334         CALL DZERO(XINT2T,NT1AMX*NRHFT*NRHFT)
335      END IF
336      IF (DO_INTXY) THEN
337         CALL DZERO(XIAJB,NT1AMX*NT1AMX)
338         CALL DZERO(YIAJB,NT1AMX*NT1AMX)
339      END IF
340
341      DO ISYMD = 1, NSYM
342         DO ILLL = 1,NBAS(ISYMD)
343            IDEL   = IBAS(ISYMD) + ILLL
344            ISYDIS = MULD2H(ISYMD,ISYM0)
345
346C           ----------------------------
347C           Work space allocation no. 2.
348C           ----------------------------
349            KXINT  = 1
350            KEND1  = KXINT + NDISAO(ISYDIS)
351            LWRK1  = LWORK - KEND1
352            IF (LWRK1 .LT. 0) THEN
353               WRITE(LUPRI,*) 'Need : ',KEND1,'Available : ',LWORK
354               CALL QUIT('Insufficient space in CCSDT_INTS0_NODDY')
355            ENDIF
356
357C           ---------------------------
358C           Read in batch of integrals.
359C           ---------------------------
360            CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND1),LWRK1,
361     *                  IDUMMY,DIRECT)
362
363C           ----------------------------------
364C           Calculate integrals needed in CC3:
365C           ----------------------------------
366            IF (DO_INTT0) THEN
367               CALL CCSDT_TRAN1(XINT1T,XINT2T,XLAMDP,XLAMDH,
368     *                          WORK(KXINT),IDEL)
369            END IF
370
371            IF (DO_INTXY) THEN
372               CALL CC3_TRAN2(XIAJB,YIAJB,XLAMDP,XLAMDH,
373     *                        WORK(KXINT),IDEL)
374            END IF
375
376            IF (DO_INTS0) THEN
377              CALL CCSDT_TRAN3(XINT1S,XINT2S,XLAMDP,XLAMDH,
378     *                         WORK(KXINT),IDEL)
379            END IF
380
381         END DO
382      END DO
383*---------------------------------------------------------------------*
384*     End Loop over distributions of integrals.
385*---------------------------------------------------------------------*
386
387      CALL QEXIT('CCSDT_INTS0_NODDY')
388      RETURN
389      END
390
391*---------------------------------------------------------------------*
392*              END OF SUBROUTINE CCSDT_INTS0_NODDY                    *
393*---------------------------------------------------------------------*
394      SUBROUTINE CCSDT_INIT_NODDY(WORK,LWORK,DO_L03AM)
395*---------------------------------------------------------------------*
396*
397* Purpose: precompute some intermediates for CC3
398*
399*          for DO_L03AM=.true. the triples part of the zeroth-order
400*          langrangian multipliers will be computed and save on file.
401*          for this case it will be assumed that all integrals
402*          are already present of file
403*
404* Written by Christof Haettig, Mai 2003
405*
406*=====================================================================*
407      IMPLICIT NONE
408#include "priunit.h"
409#include "ccsdinp.h"
410#include "maxorb.h"
411#include "ccsdsym.h"
412#include "ccfield.h"
413#include "ccorb.h"
414#include "ccnoddy.h"
415#include "dummy.h"
416
417      LOGICAL LOCDBG
418      PARAMETER (LOCDBG=.FALSE.)
419      INTEGER ISYM0
420      PARAMETER (ISYM0=1)
421
422      INTEGER LWORK
423      REAL*8  WORK(LWORK), FF, DDOT
424
425      LOGICAL DO_L03AM
426
427      CHARACTER*10 MODEL
428      INTEGER KSCR1, KFOCKD, KEND1, LWRK1, KFOCK0, KFIELD, KFIELDAO,
429     &        KLAMP0, KLAMH0, KINT1T0, KINT2T0, KINT1S0, KINT2S0,
430     &        KXIAJB, KYIAJB, K0IOVVO, K0IOOVV, K0IOOOO, K0IVVVV,
431     &        KT01AM, KT02AM, KT03AM, KT3SCR, LUSIFC, IOPT, LUFOCK,
432     &        LUTEMP, ILLL, IDEL, ISYMD, ISYDIS, KEND2, LWRK2,
433     &        KXINT
434
435      CALL QENTER('CCTINI')
436
437      IF (DIRECT)
438     &  CALL QUIT('DIRECT NOT IMPLEMENTED IN CCSDT_INIT_NODDY')
439
440*---------------------------------------------------------------------*
441*     Memory allocation:
442*---------------------------------------------------------------------*
443      KSCR1   = 1
444      KFOCKD  = KSCR1  + NT1AMX
445      KEND1   = KFOCKD + NORBT
446
447      KFOCK0  = KEND1
448      KEND1   = KFOCK0  + NORBT*NORBT
449
450      IF (NONHF) THEN
451        KFIELD   = KEND1
452        KFIELDAO = KFIELD   + NORBT*NORBT
453        KEND1    = KFIELDAO + NORBT*NORBT
454      END IF
455
456      KLAMP0  = KEND1
457      KLAMH0  = KLAMP0  + NLAMDT
458      KEND1   = KLAMH0  + NLAMDT
459
460      KINT1T0 = KEND1
461      KINT2T0 = KINT1T0 + NT1AMX*NVIRT*NVIRT
462      KEND1   = KINT2T0 + NRHFT*NRHFT*NT1AMX
463
464      KINT1S0 = KEND1
465      KINT2S0 = KINT1S0 + NT1AMX*NVIRT*NVIRT
466      KEND1   = KINT2S0 + NRHFT*NRHFT*NT1AMX
467
468      KXIAJB  = KEND1
469      KYIAJB  = KXIAJB  + NT1AMX*NT1AMX
470      KEND1   = KYIAJB  + NT1AMX*NT1AMX
471
472      K0IOVVO = KEND1
473      K0IOOVV = K0IOVVO + NRHFT*NVIRT*NVIRT*NRHFT
474      K0IOOOO = K0IOOVV + NRHFT*NVIRT*NVIRT*NRHFT
475      K0IVVVV = K0IOOOO + NRHFT*NRHFT*NRHFT*NRHFT
476      KEND1   = K0IVVVV + NVIRT*NVIRT*NVIRT*NVIRT
477
478      KT01AM  = KEND1
479      KT02AM  = KT01AM + NT1AMX
480      KT03AM  = KT02AM + NT1AMX*NT1AMX
481      KT3SCR  = KT03AM + NT1AMX*NT1AMX*NT1AMX
482      KEND1   = KT3SCR + NT1AMX*NT1AMX*NT1AMX
483
484      LWRK1   = LWORK  - KEND1
485      IF (LWRK1 .LT. 0) THEN
486         CALL QUIT('Insufficient space in CCSDT_INIT_NODDY')
487      ENDIF
488
489*---------------------------------------------------------------------*
490*     Read SCF orbital energies from file:
491*---------------------------------------------------------------------*
492      LUSIFC = -1
493      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
494     &            .FALSE.)
495      REWIND LUSIFC
496      CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
497      READ (LUSIFC)
498      READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBT)
499      CALL GPCLOSE(LUSIFC,'KEEP')
500
501*---------------------------------------------------------------------*
502*     Get zeroth-order Lambda matrices:
503*---------------------------------------------------------------------*
504      IOPT   = 1
505      Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,WORK(KT01AM),DUMMY)
506
507      Call LAMMAT(WORK(KLAMP0),WORK(KLAMH0),WORK(KT01AM),
508     &            WORK(KEND1),LWRK1)
509
510*---------------------------------------------------------------------*
511*     read zeroth-order AO Fock matrix from file:
512*---------------------------------------------------------------------*
513      LUFOCK = -1
514      CALL GPOPEN(LUFOCK,'CC_FCKH','OLD',' ','UNFORMATTED',IDUMMY,
515     &            .FALSE.)
516      REWIND(LUFOCK)
517      READ(LUFOCK) (WORK(KFOCK0-1+I),I=1,N2BST(ISYM0))
518      CALL GPCLOSE(LUFOCK,'KEEP')
519
520      CALL CC_FCKMO(WORK(KFOCK0),WORK(KLAMP0),WORK(KLAMH0),
521     &              WORK(KEND1),LWRK1,ISYM0,ISYM0,ISYM0)
522
523*---------------------------------------------------------------------*
524*     If needed get external field:
525*---------------------------------------------------------------------*
526      IF (NONHF) THEN
527        CALL DZERO(WORK(KFIELDAO),NORBT*NORBT)
528        DO I = 1, NFIELD
529          FF = EFIELD(I)
530          CALL CC_ONEP(WORK(KFIELDAO),WORK(KEND1),LWRK1,FF,1,LFIELD(I))
531        ENDDO
532        CALL DCOPY(NORBT*NORBT,WORK(KFIELDAO),1,WORK(KFIELD),1)
533
534        ! calculate external field in zero-order lambda basis
535        CALL CC_FCKMO(WORK(KFIELD),WORK(KLAMP0),WORK(KLAMH0),
536     *                WORK(KEND1),LWRK1,1,1,1)
537
538        IF (LOCDBG) WRITE(LUPRI,*) 'NORM^2(FIELD):',
539     &     DDOT(NORBT*NORBT,WORK(KFIELD),1,WORK(KFIELD),1)
540      ENDIF
541
542*---------------------------------------------------------------------*
543*     Compute some integrals:
544*           XINT1S0 =  (CK|BD)
545*           XINT2S0 =  (CK|LJ)
546*           XINT1T0 =  (KC|BD)
547*           XINT2T0 =  (KC|LJ)
548*           XIAJB   = 2(IA|JB) - (IB|JA)
549*           YIAJB   =  (IA|JB)
550*---------------------------------------------------------------------*
551      IF (DO_L03AM) THEN
552        CALL CCSDT_READ_NODDY(.FALSE.,WORK(KFOCKD),WORK(KFOCK0),
553     &                                WORK(KFIELD),WORK(KFIELDAO),
554     &                        .TRUE., WORK(KXIAJB),WORK(KYIAJB),
555     &                        .TRUE., WORK(KINT1S0),WORK(KINT2S0),
556     &                        .TRUE., WORK(KINT1T0),WORK(KINT2T0),
557     &                        .FALSE.,DUMMY,DUMMY,DUMMY,DUMMY,
558     &                        NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX)
559
560      ELSE
561        CALL CCSDT_INTS0_NODDY(.TRUE.,WORK(KXIAJB), WORK(KYIAJB),
562     &                         .TRUE.,WORK(KINT1S0),WORK(KINT2S0),
563     &                         .TRUE.,WORK(KINT1T0),WORK(KINT2T0),
564     &                         WORK(KLAMP0),WORK(KLAMH0),
565     &                         WORK(KEND1),LWRK1)
566      END IF
567
568*---------------------------------------------------------------------*
569*     compute and save zeroth-order triples amplitudes:
570*---------------------------------------------------------------------*
571      IF (.NOT. DO_L03AM) THEN
572        ! read T^0 doubles amplitudes from file and square up
573        IOPT   = 2
574        Call CC_RDRSP('R0',0,ISYM0,IOPT,MODEL,DUMMY,WORK(KT03AM))
575        CALL CC_T2SQ(WORK(KT03AM),WORK(KT02AM),ISYM0)
576
577        ! compute triples amplitudes
578        CALL CCSDT_T03AM(WORK(KT03AM),WORK(KINT1S0),WORK(KINT2S0),
579     *                   WORK(KT02AM),WORK(KSCR1),WORK(KFOCKD),
580     *                   WORK(KFIELD),WORK(KT3SCR))
581
582        CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KT03AM),1)
583
584        LUTEMP = -1
585        CALL GPOPEN(LUTEMP,FILNODT30,'UNKNOWN',' ','UNFORMATTED',
586     &              IDUMMY,.FALSE.)
587        WRITE(LUTEMP) (WORK(KT03AM+I-1), I=1,NT1AMX*NT1AMX*NT1AMX)
588        CALL GPCLOSE(LUTEMP,'KEEP')
589      END IF
590
591*---------------------------------------------------------------------*
592*     compute and save zeroth-order triples multipliers:
593*     (note that this section overwrite the cluster amplitudes on
594*      KT01AM, KT02AM and KT03AM...)
595*---------------------------------------------------------------------*
596      IF (DO_L03AM) THEN
597        ! read L^0 singles and doubles multipliers from file
598        IOPT   = 3
599        Call CC_RDRSP('L0',0,ISYM0,IOPT,MODEL,
600     *                 WORK(KT01AM),WORK(KT03AM))
601        CALL CC_T2SQ(WORK(KT03AM),WORK(KT02AM),ISYM0)
602
603        ! compute triples multipliers
604        CALL CCSDT_L03AM(WORK(KT03AM),WORK(KINT1T0),WORK(KINT2T0),
605     *                   WORK(KXIAJB),WORK(KFOCK0),WORK(KT01AM),
606     *                   WORK(KT02AM),WORK(KSCR1),WORK(KFOCKD),
607     *                   WORK(KFIELD),WORK(KT3SCR))
608
609        CALL DSCAL(NT1AMX*NT1AMX*NT1AMX,-1.0D0,WORK(KT03AM),1)
610
611        LUTEMP = -1
612        CALL GPOPEN(LUTEMP,FILNODL30,'UNKNOWN',' ','UNFORMATTED',
613     &              IDUMMY,.FALSE.)
614        WRITE(LUTEMP) (WORK(KT03AM+I-1), I=1,NT1AMX*NT1AMX*NT1AMX)
615        CALL GPCLOSE(LUTEMP,'KEEP')
616
617        CALL QEXIT('CCTINI')
618        RETURN
619      END IF
620
621*---------------------------------------------------------------------*
622*     Compute integrals needed for the following contributions:
623*---------------------------------------------------------------------*
624      CALL DZERO(WORK(K0IOVVO),NRHFT*NVIRT*NVIRT*NRHFT)
625      CALL DZERO(WORK(K0IOOVV),NRHFT*NVIRT*NVIRT*NRHFT)
626      CALL DZERO(WORK(K0IOOOO),NRHFT*NRHFT*NRHFT*NRHFT)
627      CALL DZERO(WORK(K0IVVVV),NVIRT*NVIRT*NVIRT*NVIRT)
628
629      DO ISYMD = 1, NSYM
630         DO ILLL = 1,NBAS(ISYMD)
631            IDEL   = IBAS(ISYMD) + ILLL
632            ISYDIS = MULD2H(ISYMD,ISYMOP)
633
634C           ----------------------------
635C           Work space allocation no. 2.
636C           ----------------------------
637            KXINT  = KEND1
638            KEND2  = KXINT + NDISAO(ISYDIS)
639            LWRK2  = LWORK - KEND2
640            IF (LWRK2 .LT. 0) THEN
641               CALL QUIT('Insufficient space in CCSDT_INIT_NODDY')
642            ENDIF
643
644C           ---------------------------
645C           Read in batch of integrals.
646C           ---------------------------
647            CALL CCRDAO(WORK(KXINT),IDEL,1,WORK(KEND2),LWRK2,
648     *                  IDUMMY,DIRECT) ! NB! DIRECT always .false. here
649
650            CALL CCFOP_TRAN1_R(WORK(K0IOVVO),WORK(K0IOOVV),
651     &                         WORK(K0IOOOO),WORK(K0IVVVV),
652     &                         WORK(KLAMP0),WORK(KLAMH0),
653     &                         WORK(KLAMP0),WORK(KLAMH0),
654     &                         WORK(KLAMP0),WORK(KLAMH0),
655     &                         WORK(KXINT),IDEL)
656
657         END DO
658      END DO
659
660*---------------------------------------------------------------------*
661*     Save fock-like intermediates and Lambda matrices
662*---------------------------------------------------------------------*
663      LUTEMP = -1
664      CALL GPOPEN(LUTEMP,FILNODFOCK,'UNKNOWN',' ','UNFORMATTED',
665     &            IDUMMY,.FALSE.)
666      WRITE(LUTEMP) (WORK(KFOCKD+I-1),  I=1,NORBT)
667      WRITE(LUTEMP) (WORK(KFOCK0+I-1),  I=1,NORBT*NORBT)
668      IF (NONHF) THEN
669        WRITE(LUTEMP) (WORK(KFIELD+I-1),   I=1,NORBT*NORBT)
670        WRITE(LUTEMP) (WORK(KFIELDAO+I-1), I=1,NORBT*NORBT)
671      END IF
672      CALL GPCLOSE(LUTEMP,'KEEP')
673
674*---------------------------------------------------------------------*
675*     Save integrals XINT1T0 and XINT2T0
676*---------------------------------------------------------------------*
677      LUTEMP = -1
678      CALL GPOPEN(LUTEMP,FILNODINTT,'UNKNOWN',' ','UNFORMATTED',
679     &            IDUMMY,.FALSE.)
680      WRITE(LUTEMP) (WORK(KINT1T0+I-1), I=1,NT1AMX*NVIRT*NVIRT)
681      WRITE(LUTEMP) (WORK(KINT2T0+I-1), I=1,NT1AMX*NRHFT*NRHFT)
682      CALL GPCLOSE(LUTEMP,'KEEP')
683
684*---------------------------------------------------------------------*
685*     Save integrals XINT1S0 and XINT2S0
686*---------------------------------------------------------------------*
687      LUTEMP = -1
688      CALL GPOPEN(LUTEMP,FILNODINTS,'UNKNOWN',' ','UNFORMATTED',
689     &            IDUMMY,.FALSE.)
690      WRITE(LUTEMP) (WORK(KINT1S0+I-1), I=1,NT1AMX*NVIRT*NVIRT)
691      WRITE(LUTEMP) (WORK(KINT2S0+I-1), I=1,NT1AMX*NRHFT*NRHFT)
692      CALL GPCLOSE(LUTEMP,'KEEP')
693
694*---------------------------------------------------------------------*
695*     Save integrals XIAJB and YIAJB
696*---------------------------------------------------------------------*
697      LUTEMP = -1
698      CALL GPOPEN(LUTEMP,FILNODINTX,'UNKNOWN',' ','UNFORMATTED',
699     &            IDUMMY,.FALSE.)
700      WRITE(LUTEMP) (WORK(KXIAJB+I-1),  I=1,NT1AMX*NT1AMX)
701      WRITE(LUTEMP) (WORK(KYIAJB+I-1),  I=1,NT1AMX*NT1AMX)
702      CALL GPCLOSE(LUTEMP,'KEEP')
703
704*---------------------------------------------------------------------*
705*     Save integrals XOOOO, XOVVO, XOOVV and XVVVV
706*---------------------------------------------------------------------*
707      LUTEMP = -1
708      CALL GPOPEN(LUTEMP,FILNODINTV,'UNKNOWN',' ','UNFORMATTED',
709     &            IDUMMY,.FALSE.)
710      WRITE(LUTEMP) (WORK(K0IOVVO+I-1), I=1,NRHFT*NVIRT*NVIRT*NRHFT)
711      WRITE(LUTEMP) (WORK(K0IOOVV+I-1), I=1,NRHFT*NVIRT*NVIRT*NRHFT)
712      WRITE(LUTEMP) (WORK(K0IOOOO+I-1), I=1,NRHFT*NRHFT*NRHFT*NRHFT)
713      WRITE(LUTEMP) (WORK(K0IVVVV+I-1), I=1,NVIRT*NVIRT*NVIRT*NVIRT)
714      CALL GPCLOSE(LUTEMP,'KEEP')
715
716      CALL QEXIT('CCTINI')
717      RETURN
718      END
719
720*---------------------------------------------------------------------*
721*              END OF SUBROUTINE CCSDT_INIT_NODDY                     *
722*---------------------------------------------------------------------*
723      SUBROUTINE CCSDT_READ_NODDY(RD_FOCK0,FOCKD,FOCK0,FIELD,FIELDAO,
724     &                            RD_INTXY,XIAJB,YIAJB,
725     &                            RD_INTS0,XINT1S0,XINT2S0,
726     &                            RD_INTT0,XINT1T0,XINT2T0,
727     &                            RD_INTV0,XOVVO,XOOVV,XOOOO,XVVVV,
728     &                            NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX)
729*---------------------------------------------------------------------*
730*
731* Purpose: read precomputed intermediates for CC3
732*
733* Written by Christof Haettig, Mai 2003
734*
735*=====================================================================*
736      IMPLICIT NONE
737#include "priunit.h"
738#include "ccfield.h"
739#include "ccnoddy.h"
740#include "ccsdinp.h"
741#include "dummy.h"
742
743      LOGICAL LOCDBG
744      PARAMETER (LOCDBG=.FALSE.)
745
746      INTEGER NORBT,NLAMDT,NRHFT,NVIRT,NT1AMX, LUTEMP, I
747
748      LOGICAL RD_FOCK0, RD_INTXY, RD_INTS0, RD_INTT0, RD_INTV0, RD_XLAM0
749
750      REAL*8  FOCKD(*), FOCK0(*)
751      REAL*8  FIELD(*), FIELDAO(*), XIAJB(*), YIAJB(*)
752      REAL*8  XINT1S0(*), XINT2S0(*), XINT1T0(*), XINT2T0(*)
753      REAL*8  XOVVO(*), XOOVV(*), XOOOO(*), XVVVV(*)
754
755      CALL QENTER('CCTRDI')
756
757      IF (DIRECT)
758     &  CALL QUIT('DIRECT NOT IMPLEMENTED IN CCSDT_READ_NODDY')
759
760*---------------------------------------------------------------------*
761*     Read integrals and triples cluster amplitudes on file
762*---------------------------------------------------------------------*
763      IF (RD_FOCK0) THEN
764        LUTEMP = -1
765        CALL GPOPEN(LUTEMP,FILNODFOCK,'UNKNOWN',' ','UNFORMATTED',
766     &              IDUMMY,.FALSE.)
767        READ(LUTEMP) (FOCKD(I),I=1,NORBT)
768        READ(LUTEMP) (FOCK0(I),I=1,NORBT*NORBT)
769        IF (NONHF) THEN
770          READ(LUTEMP) (FIELD(I),  I=1,NORBT*NORBT)
771          READ(LUTEMP) (FIELDAO(I),I=1,NORBT*NORBT)
772        END IF
773        CALL GPCLOSE(LUTEMP,'KEEP')
774      END IF
775
776      IF (RD_INTT0) THEN
777        LUTEMP = -1
778        CALL GPOPEN(LUTEMP,FILNODINTT,'UNKNOWN',' ','UNFORMATTED',
779     &              IDUMMY,.FALSE.)
780        READ(LUTEMP) (XINT1T0(I),I=1,NT1AMX*NVIRT*NVIRT)
781        READ(LUTEMP) (XINT2T0(I),I=1,NT1AMX*NRHFT*NRHFT)
782        CALL GPCLOSE(LUTEMP,'KEEP')
783      ENDIF
784
785      IF (RD_INTS0) THEN
786        LUTEMP = -1
787        CALL GPOPEN(LUTEMP,FILNODINTS,'UNKNOWN',' ','UNFORMATTED',
788     &              IDUMMY,.FALSE.)
789        READ(LUTEMP) (XINT1S0(I),I=1,NT1AMX*NVIRT*NVIRT)
790        READ(LUTEMP) (XINT2S0(I),I=1,NT1AMX*NRHFT*NRHFT)
791        CALL GPCLOSE(LUTEMP,'KEEP')
792      END IF
793
794      IF (RD_INTXY) THEN
795        LUTEMP = -1
796        CALL GPOPEN(LUTEMP,FILNODINTX,'UNKNOWN',' ','UNFORMATTED',
797     &              IDUMMY,.FALSE.)
798        READ(LUTEMP) (XIAJB(I),I=1,NT1AMX*NT1AMX)
799        READ(LUTEMP) (YIAJB(I),I=1,NT1AMX*NT1AMX)
800        CALL GPCLOSE(LUTEMP,'KEEP')
801      END IF
802
803      IF (RD_INTV0) THEN
804        LUTEMP = -1
805        CALL GPOPEN(LUTEMP,FILNODINTV,'UNKNOWN',' ','UNFORMATTED',
806     &              IDUMMY,.FALSE.)
807        READ(LUTEMP) (XOVVO(I),I=1,NRHFT*NVIRT*NVIRT*NRHFT)
808        READ(LUTEMP) (XOOVV(I),I=1,NRHFT*NVIRT*NVIRT*NRHFT)
809        READ(LUTEMP) (XOOOO(I),I=1,NRHFT*NRHFT*NRHFT*NRHFT)
810        READ(LUTEMP) (XVVVV(I),I=1,NVIRT*NVIRT*NVIRT*NVIRT)
811        CALL GPCLOSE(LUTEMP,'KEEP')
812      END IF
813
814      CALL QEXIT('CCTRDI')
815      RETURN
816      END
817
818*---------------------------------------------------------------------*
819*              END OF SUBROUTINE CCSDT_READ_NODDY                     *
820*---------------------------------------------------------------------*
821