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 cc_d1ao */
20      SUBROUTINE CC_D1AO(IPDD,R12PRP,AODEN,ZKAM,WORK,LWORK,MODEL,LLIST,
21     &                   ILLNR,NO,FNDPTIA, FNDPTIA2, FNDPTAB, FNDPTIJ)
22C
23C     Written by Asger Halkier January 1996.
24C
25C     Version: 1.0
26C
27C     Purpose: To calculate the one electron Coupled Cluster
28C              density matrix and returning it backtransformed
29C              to AO-basis in AODEN!
30C
31C     Current models:
32C     MP2, CCD, CCSD, CCSD(T), CC3, DRCCD (DRPA), RCCD (singlet RPA)
33C
34C     NO = .true.  compute natural occupation numbers
35C
36C     OC: 15-4-1997: general llist.
37C
38C     May 1998: MP2 Frozen core density by Asger Halkier.
39C
40C     Spring 2002 : CC3 section by K. Hald.
41C     Summer 2010 : CCD + RCCD/DRCCD (RPA) session by Sonia Coriani
42C
43#include "implicit.h"
44#include "priunit.h"
45#include "dummy.h"
46#include "maxash.h"
47#include "maxorb.h"
48#include "mxcent.h"
49#include "aovec.h"
50#include "iratdef.h"
51      logical r12prp
52      PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0)
53      DIMENSION INDEXA(MXCORB_CC)
54      DIMENSION AODEN(*), ZKAM(*), WORK(LWORK)
55#include "ccorb.h"
56#include "infind.h"
57#include "inftap.h"
58#include "blocks.h"
59#include "ccsdinp.h"
60#include "ccsdsym.h"
61#include "ccsdio.h"
62#include "distcl.h"
63#include "cbieri.h"
64#include "eribuf.h"
65#include "ccfop.h"
66#include "ccfro.h"
67C
68      CHARACTER MODEL*10,MODDUM*10
69      CHARACTER LLIST*(*)
70      CHARACTER*(*) FNDPTIA, FNDPTAB, FNDPTIJ, FNDPTIA2
71      LOGICAL NO
72C
73C----------------------------------------
74C     Initialization of timing parameter.
75C----------------------------------------
76C
77      CALL QENTER('CC_D1AO')
78      TIMTOT = ZERO
79      TIMTOT = SECOND()
80C
81      IF ((CC2) .AND. (RELORB)) THEN
82C
83         CALL CC2_D1AO(AODEN,ZKAM,WORK,LWORK)
84         CALL QEXIT('CC_D1AO')
85         RETURN
86C
87      ELSE IF (MP2) THEN
88C
89         KONEAI = 1
90         KONEAB = KONEAI + NT1AMX
91         KONEIJ = KONEAB + NMATAB(1)
92         KONEIA = KONEIJ + NMATIJ(1)
93         KT1AM  = KONEIA + NT1AMX
94         KLAMDH = KT1AM  + NT1AMX
95         KLAMDP = KLAMDH + NLAMDT
96         KRMAT  = KLAMDP + NLAMDT
97         KEND1  = KRMAT  + NMATIJ(1)
98         LWRK1  = LWORK  - KEND1
99C
100         IF (LWRK1 .LT. 0) THEN
101            WRITE(LUPRI,*) 'MP2 - Available:', LWORK, 'Needed:', KEND1
102            CALL QUIT('Insufficient memory for work '//
103     &           'allocation in CC_D1AO for MP2')
104         ENDIF
105C
106C-----------------------------------------------------------------
107C        Initialize arrays (note that the t1-amplitudes are zero).
108C-----------------------------------------------------------------
109C
110         CALL DZERO(WORK(KONEAI),NT1AMX)
111         CALL DZERO(WORK(KONEAB),NMATAB(1))
112         CALL DZERO(WORK(KONEIJ),NMATIJ(1))
113         CALL DZERO(WORK(KONEIA),NT1AMX)
114         CALL DZERO(WORK(KT1AM),NT1AMX)
115         CALL DZERO(WORK(KRMAT),NMATIJ(1))
116C
117C--------------------------------------------------
118C        Set up the relaxation part of the density.
119C--------------------------------------------------
120C
121         CALL DCOPY(NMATIJ(1),ZKAM(1),1,WORK(KONEIJ),1)
122         CALL DCOPY(NMATAB(1),ZKAM(NMATIJ(1)+1),1,WORK(KONEAB),1)
123         CALL DCOPY(NT1AMX,ZKAM(NMATIJ(1)+NMATAB(1)+1),1,WORK(KONEAI),1)
124         CALL DCOPY(NT1AMX,WORK(KONEAI),1,WORK(KONEIA),1)
125C
126C-------------------------------------
127C        Add the Hartree-Fock density.
128C        Add R12/A' contribution (Elena)
129C-------------------------------------
130C
131celena
132         IF (R12PRP .AND. (IPDD .EQ. 3 .OR. IPDD .EQ. 5)) THEN
133            lunit = -1
134            IF (IPDD .EQ. 3) THEN
135               call gpopen(lunit,'CCR12YIJ','unknown',' ','unformatted',
136     &                     idummy,.false.)
137             ELSEIF (IPDD .EQ. 5) THEN
138               call gpopen(lunit,'CCR12YIJB','unknown',' ',
139     &                     'unformatted',idummy,.false.)
140             ENDIF
141
142           do isymaj = 1,nsym
143              isymij = isymaj
144              ncvij = 0
145              do isyma = 1,nsym
146                 isymj = muld2h(isymaj,isyma)
147                 isymi = muld2h(isymij,isymj)
148                 ncvij = ncvij +  nrhf(isyma)*nrhf(isymi)
149              enddo
150            enddo
151
152            read(lunit)(work(krmat+i-1),i=1,ncvij)
153            call gpclose(lunit,'KEEP')
154
155            DO ISYMI = 1,NSYM
156               ISYMJ = ISYMI
157               DO  J = 1,NRHF(ISYMJ)
158                   DO  I = J+1,NRHF(ISYMI)
159                       NIJ   = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)
160     *                          *(J - 1) + I
161                       NJI   = IMATIJ(ISYMJ,ISYMI) + NRHF(ISYMJ)
162     *                          *(I - 1) + J
163                       WORK(KONEIJ + NIJ - 1) = WORK(KONEIJ + NIJ - 1)
164     *                          + work(krmat + nij - 1)
165                       WORK(KONEIJ + NJI - 1) = WORK(KONEIJ + NJI - 1)
166     *                          + work(krmat + nji - 1)
167C
168                   ENDDO
169               ENDDO
170            ENDDO
171
172         ENDIF
173
174celena
175
176         DO 80 ISYM = 1,NSYM
177            DO 90 I = 1,NRHF(ISYM)
178               NII = IMATIJ(ISYM,ISYM) + NRHF(ISYM)*(I - 1) + I
179               WORK(KONEIJ + NII - 1) = WORK(KONEIJ + NII - 1) + TWO
180               if (r12prp .and. (ipdd .eq. 3 .or. ipdd .eq. 5)) then
181                  WORK(KONEIJ + NII - 1) = WORK(KONEIJ + NII - 1) +
182     &                   1.0d0*work(krmat + nii - 1)
183               endif
184   90       CONTINUE
185   80    CONTINUE
186C
187C----------------------------------------
188C        Calculate MO coefficient matrix.
189C----------------------------------------
190C
191         CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),
192     *               WORK(KEND1),LWRK1)
193C
194C--------------------------------------
195C        Transform density to AO basis.
196C--------------------------------------
197C
198         CALL DZERO(AODEN,N2BST(1))
199C
200         ISDEN = 1
201         CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB),
202     *                 WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1,
203     *                 WORK(KLAMDH),1,WORK(KEND1),LWRK1)
204C
205C----------------------------------------------------
206C        Add frozen core contributions to AO density.
207C----------------------------------------------------
208C
209         IF (FROIMP) THEN
210C
211            KOFFAI = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 1
212            KOFFIA = KOFFAI    + NT1FRO(1)
213            KOFFIJ = KOFFIA    + NT1FRO(1)
214            KOFFJI = KOFFIJ    + NCOFRO(1)
215C
216            ISDEN = 1
217            ICON  = 1
218            CALL CC_D1FCB(AODEN,ZKAM(KOFFIJ),ZKAM(KOFFJI),ZKAM(KOFFAI),
219     *                    ZKAM(KOFFIA),WORK(KEND1),LWRK1,ISDEN,ICON)
220C
221         ENDIF
222C
223      ELSE IF ((CCSD .OR. CCPT .OR. CC3 .OR. CCD)
224     *          .OR. ((CC2) .AND. (.NOT. RELORB))) THEN
225C
226C-------------------------------------------------------
227C        Both zeta- and t-vectors are totally symmetric.
228C-------------------------------------------------------
229C
230         ISYMTR = 1
231         ISYMOP = 1
232C
233C--------------------------------------
234C        Initial work space allocation.
235C--------------------------------------
236C
237         NHELP = NT2AMX
238C
239         KZ2AM  = 1
240         KT2AM  = KZ2AM  + NT2SQ(1)
241         KT2AMT = KT2AM  + NT2AMX
242         KLAMDP = KT2AMT + NHELP
243         KLAMDH = KLAMDP + NLAMDT
244         KT1AM  = KLAMDH + NLAMDT
245         KZ1AM  = KT1AM  + NT1AMX
246         KEND1  = KZ1AM  + NT1AMX
247         LWRK1  = LWORK  - KEND1
248C
249         IF (LWRK1 .LT. 0) THEN
250            WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1
251            CALL QUIT('Insufficient core for initial '//
252     &           'allocation in CC_D1AO')
253         ENDIF
254C
255C-------------------------------------------
256C        Read zero'th order zeta amplitudes.
257C-------------------------------------------
258C
259         ISYML = 1
260         IOPT = 3
261         CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODDUM,
262     &                 WORK(KZ1AM),WORK(KT2AM))
263         IF ( IPRINT .GT. 10 ) THEN
264            RHO1N = DDOT(NT1AM(ISYML),WORK(KZ1AM),1,WORK(KZ1AM),1)
265            RHO2N = DDOT(NT2AM(ISYML),WORK(KT2AM),1,WORK(KT2AM),1)
266            WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO1N
267            WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO2N
268         ENDIF
269
270C
271C-----------------------------------
272C        Square up zeta2 amplitudes.
273C-----------------------------------
274C
275         CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1)
276C
277C----------------------------------------------
278C        Read zero'th order cluster amplitudes.
279C----------------------------------------------
280C
281         IOPT = 3
282         CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM))
283C
284C---------------------------------------------------
285C        Zero out single vectors in CCD-calculation.
286C---------------------------------------------------
287C
288         IF (CCD) THEN
289            CALL DZERO(WORK(KT1AM),NT1AMX)
290            CALL DZERO(WORK(KZ1AM),NT1AMX)
291         ENDIF
292C
293C-------------------------------------
294C        Calculate the lamda matrices.
295C-------------------------------------
296C
297         CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1),
298     *               LWRK1)
299C
300C------------------------------------------
301C        Set up 2C-E of cluster amplitudes.
302C------------------------------------------
303C
304         CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1)
305         ISYOPE = 1
306         IOPTTCME = 1
307         CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME)
308C
309C----------------------------------
310C        Work space allocation one.
311C        Note that D(ai)=ZETA(ai) &
312C        D(ia) is stored transposed
313C----------------------------------
314C
315         KONEAI = KZ1AM
316         KONEAB = KONEAI + NT1AMX
317         KONEIJ = KONEAB + NMATAB(1)
318         KONEIA = KONEIJ + NMATIJ(1)
319         KXMAT  = KONEIA + NT1AMX
320         KYMAT  = KXMAT  + NMATIJ(1)
321         KCMO   = KYMAT  + NMATAB(1)
322         KEND1  = KCMO   + NLAMDS
323         LWRK1  = LWORK  - KEND1
324C
325         IF (LWRK1 .LT. 0) THEN
326            WRITE(LUPRI,*) 'CC_D1AO; Available:',LWORK, 'Needed:',KEND1
327            CALL QUIT(
328     &           'Insufficient memory for first allocation in CC_D1AO')
329         ENDIF
330C
331C---------------------------------------------------------
332C        Initialize remaining one electron density arrays.
333C---------------------------------------------------------
334C
335         CALL DZERO(WORK(KONEAB),NMATAB(1))
336         CALL DZERO(WORK(KONEIJ),NMATIJ(1))
337         CALL DZERO(WORK(KONEIA),NT1AMX)
338C
339C-----------------------------------------------------------
340C        Calculate X-intermediate of zeta- and t-amplitudes.
341C-----------------------------------------------------------
342C
343         CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP,
344     *              WORK(KEND1),LWRK1)
345C
346C-----------------------------------------------------------
347C        Calculate Y-intermediate of zeta- and t-amplitudes.
348C-----------------------------------------------------------
349C
350         CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP,
351     *              WORK(KEND1),LWRK1)
352C
353C-----------------------------------------------------------------
354C        Construct three remaining blocks of one electron density.
355C-----------------------------------------------------------------
356C
357         CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1)
358         CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1)
359         CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT))
360         CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI))
361C
362C------------------------------------------------------------
363C        Get the contribution from CC3.
364C------------------------------------------------------------
365C
366         IF (CC3) THEN
367C
368            IF (LWRK1 .LT. MAX(NT1AMX,NMATIJ(1),NMATAB(1))) THEN
369               CALL QUIT('Out of memory in CC_D1AO (CC3)')
370            ENDIF
371C
372            LUPTIA = -1
373            CALL WOPEN2(LUPTIA,FNDPTIA,64,0)
374            IOFF = 1 + NT1AMX*(ILLNR-1)
375            CALL GETWA2(LUPTIA,FNDPTIA,WORK(KEND1),IOFF,NT1AMX)
376            CALL WCLOSE2(LUPTIA,FNDPTIA,'KEEP')
377            CALL DAXPY(NT1AMX,1.0D0,WORK(KEND1),1,WORK(KONEIA),1)
378C
379            LUPTAB = -1
380            CALL WOPEN2(LUPTAB,FNDPTAB,64,0)
381            IOFF = 1 + NMATAB(1)*(ILLNR-1)
382            CALL GETWA2(LUPTAB,FNDPTAB,WORK(KEND1),IOFF,NMATAB(1))
383            CALL WCLOSE2(LUPTAB,FNDPTAB,'KEEP')
384            CALL DAXPY(NMATAB(1),1.0D0,WORK(KEND1),1,WORK(KONEAB),1)
385C
386            LUPTIJ = -1
387            CALL WOPEN2(LUPTIJ,FNDPTIJ,64,0)
388            IOFF = 1 + NMATIJ(1)*(ILLNR-1)
389            CALL GETWA2(LUPTIJ,FNDPTIJ,WORK(KEND1),IOFF,NMATIJ(1))
390            CALL WCLOSE2(LUPTIJ,FNDPTIJ,'KEEP')
391            CALL DAXPY(NMATIJ(1),1.0D0,WORK(KEND1),1,WORK(KONEIJ),1)
392C
393            LUPTIA2 = -1
394            CALL WOPEN2(LUPTIA2,FNDPTIA2,64,0)
395            IOFF = 1 + NT1AMX*(ILLNR-1)
396            CALL GETWA2(LUPTIA2,FNDPTIA2,WORK(KEND1),IOFF,NT1AMX)
397            CALL WCLOSE2(LUPTIA2,FNDPTIA2,'KEEP')
398            CALL DAXPY(NT1AMX,1.0D0,WORK(KEND1),1,WORK(KONEIA),1)
399         ENDIF
400C
401C-----------------------------------------------------------
402C        Backtransform the one electron density to AO-basis.
403C-----------------------------------------------------------
404C
405         CALL DZERO(AODEN,N2BST(1))
406C
407         ISDEN = 1
408         CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB),
409     *                 WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1,
410     *                 WORK(KLAMDH),1,WORK(KEND1),LWRK1)
411C
412         IF (RELORB) THEN
413C
414C----------------------------------------------------------------
415C           Read MO-coefficients from interface file and reorder.
416C----------------------------------------------------------------
417C
418            LUSIFC = -1
419            CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
420     &                  .FALSE.)
421            REWIND LUSIFC
422            CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
423            READ (LUSIFC)
424            READ (LUSIFC)
425            READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS)
426            CALL GPCLOSE(LUSIFC,'KEEP')
427C
428            CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1)
429C
430C-------------------------------------
431C        Add orbital relaxation terms.
432C-------------------------------------
433C
434c           CALL CC_D1ORRE(AODEN,ZKAM(NMATIJ(1)+NMATAB(1)+1),
435c    *                     WORK(KEND1),LWRK1)
436C
437            KOFDIJ = 1
438            KOFDAB = KOFDIJ + NMATIJ(1)
439            KOFDAI = KOFDAB + NMATAB(1)
440            KOFDIA = KOFDAI + NT1AMX
441C
442            ISDEN = 1
443            CALL CC_DENAO(AODEN,1,ZKAM(KOFDAI),ZKAM(KOFDAB),
444     *                    ZKAM(KOFDIJ),ZKAM(KOFDIA),ISDEN,WORK(KCMO),1,
445     *                    WORK(KCMO),1,WORK(KEND1),LWRK1)
446C
447C
448C-------------------------------------------------------
449C           Add frozen core contributions to AO density.
450C-------------------------------------------------------
451C
452            IF (FROIMP) THEN
453C
454               KOFFAI = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 1
455               KOFFIA = KOFFAI    + NT1FRO(1)
456               KOFFIJ = KOFFIA    + NT1FRO(1)
457               KOFFJI = KOFFIJ    + NCOFRO(1)
458C
459               ISDEN = 1
460               ICON  = 2
461               CALL CC_D1FCB(AODEN,ZKAM(KOFFIJ),ZKAM(KOFFJI),
462     *                       ZKAM(KOFFAI),ZKAM(KOFFIA),WORK(KEND1),
463     *                       LWRK1,ISDEN,ICON)
464C
465            ENDIF
466         ENDIF
467C
468C
469      ELSE IF ((RCCD).or.(DRCCD)) THEN
470C
471C-------------------------------------------------------
472C        Both zeta- and t-vectors are totally symmetric.
473C-------------------------------------------------------
474C
475         ISYMTR = 1
476         ISYMOP = 1
477C
478C--------------------------------------
479C        Initial work space allocation.
480C--------------------------------------
481C
482         NHELP = NT2AMX
483C
484         KZ2AM  = 1
485         KT2AM  = KZ2AM  + NT2SQ(1)
486         KT2AMT = KT2AM  + NT2AMX
487         KLAMDP = KT2AMT + NHELP
488         KLAMDH = KLAMDP + NLAMDT
489         KT1AM  = KLAMDH + NLAMDT
490         KZ1AM  = KT1AM  + NT1AMX
491         KEND1  = KZ1AM  + NT1AMX
492         LWRK1  = LWORK  - KEND1
493C
494         IF (LWRK1 .LT. 0) THEN
495            WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1
496            CALL QUIT('Insufficient core for initial '//
497     &           'allocation in CC_D1AO')
498         ENDIF
499C
500C-------------------------------------------
501C        Read zero'th order zeta amplitudes.
502C-------------------------------------------
503C
504         ISYML = 1
505         IOPT = 2
506         CALL CC_RDRSP(LLIST,ILLNR,ISYML,IOPT,MODDUM,
507     &                 WORK(KZ1AM),WORK(KT2AM))
508         IF ( IPRINT .GT. 10 ) THEN
509            RHO1N = DDOT(NT1AM(ISYML),WORK(KZ1AM),1,WORK(KZ1AM),1)
510            RHO2N = DDOT(NT2AM(ISYML),WORK(KT2AM),1,WORK(KT2AM),1)
511            WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO1N
512            WRITE(LUPRI,*) 'CC_D1AO: Norm of L. vector :',RHO2N
513         ENDIF
514
515C
516C-----------------------------------
517C        Square up zeta2 amplitudes.
518C-----------------------------------
519C
520         CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1)
521C
522C----------------------------------------------
523C        Read zero'th order cluster amplitudes.
524C----------------------------------------------
525C
526         IOPT = 2
527         CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM))
528C
529C---------------------------------------------------
530C        Zero out single vectors (do no exist)
531C---------------------------------------------------
532C
533         CALL DZERO(WORK(KT1AM),NT1AMX)
534         CALL DZERO(WORK(KZ1AM),NT1AMX)
535C
536C-------------------------------------------------
537C        Calculate the lambda matrices.
538C        FIXME Useless (since T1=0) but I am lazy. Sonia
539C-------------------------------------------------
540C
541         CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1),
542     *               LWRK1)
543C
544C------------------------------------------
545C        Set up 2C-E of cluster amplitudes.
546C------------------------------------------
547C
548         CALL DCOPY(NT2AMX,WORK(KT2AM),1,WORK(KT2AMT),1)
549         ISYOPE = 1
550         IOPTTCME = 1
551         CALL CCSD_TCMEPK(WORK(KT2AMT),1.0D0,ISYOPE,IOPTTCME)
552C
553C----------------------------------
554C        Work space allocation one.
555C        Note that D(ai)=ZETA(ai) &
556C        D(ia) is stored transposed
557C----------------------------------
558C
559         KONEAI = KZ1AM
560         KONEAB = KONEAI + NT1AMX
561         KONEIJ = KONEAB + NMATAB(1)
562         KONEIA = KONEIJ + NMATIJ(1)
563         KXMAT  = KONEIA + NT1AMX
564         KYMAT  = KXMAT  + NMATIJ(1)
565         KCMO   = KYMAT  + NMATAB(1)
566         KEND1  = KCMO   + NLAMDS
567         LWRK1  = LWORK  - KEND1
568C
569         IF (LWRK1 .LT. 0) THEN
570            WRITE(LUPRI,*) 'CC_D1AO; Available:',LWORK, 'Needed:',KEND1
571            CALL QUIT(
572     &           'Insufficient memory for first allocation in CC_D1AO')
573         ENDIF
574C
575C---------------------------------------------------------
576C        Initialize remaining one electron density arrays.
577C---------------------------------------------------------
578C
579         CALL DZERO(WORK(KONEAB),NMATAB(1))
580         CALL DZERO(WORK(KONEIJ),NMATIJ(1))
581         CALL DZERO(WORK(KONEIA),NT1AMX)
582         CALL DZERO(WORK(KONEAI),NT1AMX)
583C
584C-----------------------------------------------------------
585C        Calculate X-intermediate of zeta- and t-amplitudes.
586C-----------------------------------------------------------
587C
588         CALL CC_XI(WORK(KXMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP,
589     *              WORK(KEND1),LWRK1)
590         CALL DSCAL(NMATIJ(1),TWO,WORK(KXMAT),1)
591C
592C-----------------------------------------------------------
593C        Calculate Y-intermediate of zeta- and t-amplitudes.
594C-----------------------------------------------------------
595C
596         CALL CC_YI(WORK(KYMAT),WORK(KZ2AM),ISYMTR,WORK(KT2AM),ISYMOP,
597     *              WORK(KEND1),LWRK1)
598         CALL DSCAL(NMATAB(1),TWO,WORK(KYMAT),1)
599C
600C-----------------------------------------------------------------
601C        Construct nonzero blocks of one electron density.
602C-----------------------------------------------------------------
603C
604         CALL DCOPY(NMATAB(1),WORK(KYMAT),1,WORK(KONEAB),1)
605         CALL CC_EITR(WORK(KONEAB),WORK(KONEIJ),WORK(KEND1),LWRK1,1)
606         !pure Hartree-Fock contribution to the density is added inside DIJGEN
607         CALL DIJGEN(WORK(KONEIJ),WORK(KXMAT))
608C
609C-----------------------------------------------------------
610C        Backtransform the one electron density to AO-basis.
611C-----------------------------------------------------------
612C
613         CALL DZERO(AODEN,N2BST(1))
614C
615         ISDEN = 1
616         CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB),
617     *                 WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1,
618     *                 WORK(KLAMDH),1,WORK(KEND1),LWRK1)
619C
620         IF (RELORB) THEN
621C
622C----------------------------------------------------------------
623C           Read MO-coefficients from interface file and reorder.
624C----------------------------------------------------------------
625C
626            LUSIFC = -1
627            CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
628     &                  .FALSE.)
629            REWIND LUSIFC
630            CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
631            READ (LUSIFC)
632            READ (LUSIFC)
633            READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS)
634            CALL GPCLOSE(LUSIFC,'KEEP')
635C
636            CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1)
637C
638C-------------------------------------
639C        Add orbital relaxation terms.
640C        (backtransform kappabar_pq passed from outside)
641C-------------------------------------
642C
643c           CALL CC_D1ORRE(AODEN,ZKAM(NMATIJ(1)+NMATAB(1)+1),
644c    *                     WORK(KEND1),LWRK1)
645            KOFDIJ = 1
646            KOFDAB = KOFDIJ + NMATIJ(1)
647            KOFDAI = KOFDAB + NMATAB(1)
648            KOFDIA = KOFDAI + NT1AMX
649C
650!            write(lupri,*)'The IJ kappabar of ', MODEL(1:5)
651!         CALL OUTPUT(ZKAM(KOFDIJ),1,NRHF(1),1,NRHF(1),
652!     *                  NRHF(1),NRHF(1),1,LUPRI)
653!         write(lupri,*)'The AI kappabar of ', MODEL(1:5)
654!         CALL OUTPUT(ZKAM(KOFDAI),1,NVIR(1),1,NRHF(1),
655!     *                  NVIR(1),NRHF(1),1,LUPRI)
656!
657!!         CALL DSCAL(NT1AMX,-ONE,ZKAM(KOFDIA),1)
658!         write(lupri,*)'The IA kappabar of ', MODEL(1:5)
659!         CALL OUTPUT(ZKAM(KOFDIA),1,NVIR(1),1,NRHF(1),
660!     *                  NVIR(1),NRHF(1),1,LUPRI)
661!
662!         write(lupri,*)'The AB kappabar of ', MODEL(1:5)
663!         CALL OUTPUT(ZKAM(KOFDAB),1,NVIR(1),1,NVIR(1),
664!     *                  NVIR(1),NVIR(1),1,LUPRI)
665
666            ISDEN = 1
667            CALL CC_DENAO(AODEN,1,ZKAM(KOFDAI),ZKAM(KOFDAB),
668     *                    ZKAM(KOFDIJ),ZKAM(KOFDIA),ISDEN,WORK(KCMO),1,
669     *                    WORK(KCMO),1,WORK(KEND1),LWRK1)
670C
671C-------------------------------------------------------
672C           Add frozen core contributions to AO density.
673C-------------------------------------------------------
674C
675!            IF (FROIMP) THEN
676C
677!               KOFFAI = NMATIJ(1) + NMATAB(1) + 2*NT1AMX + 1
678!               KOFFIA = KOFFAI    + NT1FRO(1)
679!               KOFFIJ = KOFFIA    + NT1FRO(1)
680!               KOFFJI = KOFFIJ    + NCOFRO(1)
681C
682!               ISDEN = 1
683!               ICON  = 2
684!               CALL CC_D1FCB(AODEN,ZKAM(KOFFIJ),ZKAM(KOFFJI),
685!     *                       ZKAM(KOFFAI),ZKAM(KOFFIA),WORK(KEND1),
686!     *                       LWRK1,ISDEN,ICON)
687C
688!            ENDIF
689         ENDIF
690C
691      ENDIF
692C
693      IF (FROIMP) THEN
694         CALL QEXIT('CC_D1AO')
695         RETURN
696      END IF
697C
698C-----------------------------------------------------------------
699C     Calculate additional terms to the one electron density
700C     needed for the evaluation of the natural occupation numbers.
701C     Only if NO, else skip.
702C-----------------------------------------------------------------
703C
704      IF (.NOT. NO) THEN
705         CALL QEXIT('CC_D1AO')
706         RETURN
707      END IF
708C
709      IF (.NOT. MP2) THEN
710C
711         CALL CCD1T1CO(WORK(KONEAI),WORK(KONEAB),WORK(KONEIJ),
712     *                 WORK(KONEIA),WORK(KT1AM),WORK(KYMAT),
713     *                 WORK(KEND1),LWRK1)
714C
715      ENDIF
716C
717C--------------------------------------------------
718C     Reorder and diagonalize one electron density.
719C--------------------------------------------------
720C
721      KDENFO = KEND1
722      KEND2  = KDENFO + N2BST(1)
723      LWRK2  = LWORK  - KEND2
724C
725      IF (LWRK2 .LT. 0) THEN
726         WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2
727         CALL QUIT('Insufficient memory for second '//
728     &        'allocation in CC_D1AO')
729      ENDIF
730C
731      CALL DZERO(WORK(KDENFO),N2BST(1))
732C
733      CALL CCD1REOR(WORK(KDENFO),WORK(KONEAI),WORK(KONEAB),
734     *              WORK(KONEIJ),WORK(KONEIA))
735C
736      IF (IPRINT .GT. 10) THEN
737C
738         CALL AROUND('Nondiagonalised CC one electron density')
739C
740         DO 100 ISYM = 1,NSYM
741C
742            WRITE(LUPRI,333) 'Symmetry block number:', ISYM
743  333       FORMAT(3X,A22,2X,I1)
744C
745            KOFFAI = KONEAI + IT1AM(ISYM,ISYM)
746            KOFFAB = KONEAB + IMATAB(ISYM,ISYM)
747            KOFFIJ = KONEIJ + IMATIJ(ISYM,ISYM)
748            KOFFIA = KONEIA + IT1AM(ISYM,ISYM)
749C
750            CALL AROUND('DAI')
751            CALL OUTPUT(WORK(KOFFAI),1,NVIR(ISYM),1,NRHF(ISYM),
752     *                  NVIR(ISYM),NRHF(ISYM),1,LUPRI)
753C
754            CALL AROUND('DAB')
755            CALL OUTPUT(WORK(KOFFAB),1,NVIR(ISYM),1,NVIR(ISYM),
756     *                  NVIR(ISYM),NVIR(ISYM),1,LUPRI)
757C
758            CALL AROUND('DIJ')
759            CALL OUTPUT(WORK(KOFFIJ),1,NRHF(ISYM),1,NRHF(ISYM),
760     *                  NRHF(ISYM),NRHF(ISYM),1,LUPRI)
761C
762            CALL AROUND('DIA')
763            CALL OUTPUT(WORK(KOFFIA),1,NVIR(ISYM),1,NRHF(ISYM),
764     *                  NVIR(ISYM),NRHF(ISYM),1,LUPRI)
765C
766  100    CONTINUE
767C
768      ENDIF
769C
770      CALL CC_DEDIAN(WORK(KDENFO),MODEL,WORK(KEND2),LWRK2)
771C
772C-----------------------
773C     write out timings.
774C-----------------------
775C
776      TIMTOT = SECOND() -TIMTOT
777C
778      IF (IPRINT .GT. 9) THEN
779         WRITE(LUPRI,*) ' '
780         WRITE(LUPRI,*) ' '
781         WRITE(LUPRI,*) 'CC one electron density '//
782     &        'calculated & diagonalized'
783         WRITE(LUPRI,*) 'Total time used in CC_D1AO:', TIMTOT
784      ENDIF
785C
786      CALL QEXIT('CC_D1AO')
787      RETURN
788      END
789C  /* Deck cc_denao */
790      SUBROUTINE CC_DENAO(AODEN,ISYDAO,DENAI,DENAB,DENIJ,DENIA,ISYDMO,
791     *                    XLAMDP,ISYMLP,XLAMDH,ISYMLH,WORK,LWORK)
792C
793C     Written by Asger Halkier 26/1 - 1996
794C
795C     Version: 1.0
796C
797C     Symmetries:   ISYDAO  --  AODEN
798C                   ISYDMO  --  DENAI,DENAB,DENIJ,DENIA
799C                   ISYMLP  --  XLAMDP
800C                   ISYMLH  --  XLAMDH
801C
802C     Purpose: To transform all four blocks of the Coupled
803C              Cluster one electron density to AO-basis and add
804C              the results!
805C
806#include "implicit.h"
807      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
808      DIMENSION AODEN(*), DENAI(*), DENAB(*), DENIJ(*), DENIA(*)
809      DIMENSION XLAMDP(*), XLAMDH(*), WORK(LWORK)
810#include "priunit.h"
811#include "ccorb.h"
812#include "ccsdsym.h"
813C
814      CALL QENTER('CC_DENAO')
815C
816      DO ISYM1 = 1,NSYM
817C
818         ISYM2 = MULD2H(ISYDMO,ISYM1)
819         ISYMA = MULD2H(ISYMLP,ISYM1)
820         ISYMB = MULD2H(ISYMLH,ISYM2)
821C
822         IF (ISYDAO.NE.MULD2H(ISYMA,ISYMB)) THEN
823           CALL QUIT('Symmetry mismatch in CC_DENAO')
824         END IF
825C
826         LVIR  = MAX(NBAS(ISYMA)*NVIR(ISYM2),NBAS(ISYMB)*NVIR(ISYM1))
827         LRHF  = MAX(NBAS(ISYMA)*NRHF(ISYM2),NBAS(ISYMB)*NRHF(ISYM2))
828         LNEED = MAX(LVIR,LRHF)
829C
830         IF (LWORK .LT. LNEED) THEN
831            WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', LNEED
832            CALL QUIT('Insufficient work space in CC_DENAO')
833         ENDIF
834C
835         CALL DZERO(WORK,LNEED)
836C
837         NBASA  = MAX(NBAS(ISYMA),1)
838         NBASB  = MAX(NBAS(ISYMB),1)
839         NTOVI1 = MAX(NVIR(ISYM1),1)
840         NTOVI2 = MAX(NVIR(ISYM2),1)
841         NTORH1 = MAX(NRHF(ISYM1),1)
842C
843C-----------------------------------------
844C        Transform virtual-occupied block.
845C-----------------------------------------
846C
847         KOFF1  = IGLMVI(ISYMA,ISYM1) + 1
848         KOFF2  =  IT1AM(ISYM1,ISYM2) + 1
849C
850         CALL DGEMM('N','N',NBAS(ISYMA),NRHF(ISYM2),NVIR(ISYM1),
851     *              ONE,XLAMDP(KOFF1),NBASA,DENAI(KOFF2),NTOVI1,
852     *              ZERO,WORK,NBASA)
853C
854         KOFF3  = IGLMRH(ISYMB,ISYM2) + 1
855         KOFF4  = IAODIS(ISYMA,ISYMB) + 1
856C
857         CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYM2),
858     *              ONE,WORK,NBASA,XLAMDH(KOFF3),NBASB,ONE,
859     *              AODEN(KOFF4),NBASA)
860C
861C----------------------------------------
862C        Transform virtual-virtual block.
863C----------------------------------------
864C
865         KOFF1  = IGLMVI(ISYMA,ISYM1) + 1
866         KOFF2  = IMATAB(ISYM1,ISYM2) + 1
867C
868         CALL DGEMM('N','N',NBAS(ISYMA),NVIR(ISYM2),NVIR(ISYM1),
869     *              ONE,XLAMDP(KOFF1),NBASA,DENAB(KOFF2),NTOVI1,
870     *              ZERO,WORK,NBASA)
871C
872         KOFF3  = IGLMVI(ISYMB,ISYM2) + 1
873         KOFF4  = IAODIS(ISYMA,ISYMB) + 1
874C
875         CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NVIR(ISYM2),
876     *              ONE,WORK,NBASA,XLAMDH(KOFF3),NBASB,ONE,
877     *              AODEN(KOFF4),NBASA)
878C
879C------------------------------------------
880C        Transform occupied-occupied block.
881C------------------------------------------
882C
883         KOFF1  = IGLMRH(ISYMA,ISYM1) + 1
884         KOFF2  = IMATIJ(ISYM1,ISYM2) + 1
885C
886C
887         CALL DGEMM('N','N',NBAS(ISYMA),NRHF(ISYM2),NRHF(ISYM1),
888     *              ONE,XLAMDP(KOFF1),NBASA,DENIJ(KOFF2),NTORH1,
889     *              ZERO,WORK,NBASA)
890C
891         KOFF3  = IGLMRH(ISYMB,ISYM2) + 1
892         KOFF4  = IAODIS(ISYMA,ISYMB) + 1
893C
894         CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYM2),
895     *              ONE,WORK,NBASA,XLAMDH(KOFF3),NBASB,ONE,
896     *              AODEN(KOFF4),NBASA)
897C
898C-------------------------------------------------------------
899C        Transform occupied-virtual block (stored transposed).
900C-------------------------------------------------------------
901C
902         KOFF1  = IGLMVI(ISYMB,ISYM2) + 1
903         KOFF2  =  IT1AM(ISYM2,ISYM1) + 1
904C
905         CALL DGEMM('N','N',NBAS(ISYMB),NRHF(ISYM1),NVIR(ISYM2),
906     *              ONE,XLAMDH(KOFF1),NBASB,DENIA(KOFF2),NTOVI2,
907     *              ZERO,WORK,NBASB)
908C
909         KOFF3  = IGLMRH(ISYMA,ISYM1) + 1
910         KOFF4  = IAODIS(ISYMA,ISYMB) + 1
911C
912         CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYM1),
913     *              ONE,XLAMDP(KOFF3),NBASA,WORK,NBASB,ONE,
914     *              AODEN(KOFF4),NBASA)
915C
916      END DO
917C
918      CALL QEXIT('CC_DENAO')
919C
920      RETURN
921      END
922C  /* Deck dijgen */
923      SUBROUTINE DIJGEN(DIJ,XMAT)
924C
925C     Written by Asger Halkier 26/1 - 1996
926C
927C     Version: 1.0
928C
929C     Purpose: To set up the (occ,occ) part of the Coupled
930C              Cluster one electron density.
931C
932#include "implicit.h"
933      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
934      DIMENSION DIJ(*), XMAT(*)
935#include "priunit.h"
936#include "ccorb.h"
937#include "ccsdsym.h"
938C
939C--------------------------
940C     Set up density block.
941C--------------------------
942C
943      CALL QENTER('DIJGEN')
944C
945      CALL DAXPY(NMATIJ(1),-ONE,XMAT,1,DIJ,1)
946C
947      DO 100 ISYMI = 1,NSYM
948C
949         ISYMJ = ISYMI
950C
951         DO 110 I = 1,NRHF(ISYMI)
952C
953            NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(I - 1) + I
954C
955            DIJ(NIJ) = DIJ(NIJ) + TWO
956C
957  110    CONTINUE
958C
959  100 CONTINUE
960C
961      CALL QEXIT('DIJGEN')
962C
963      RETURN
964      END
965C  /* Deck diagen */
966      SUBROUTINE DIAGEN(DIA,T2AMT,Z1AM)
967C
968C     Written by Asger Halkier 26/1 - 1996
969C
970C     Version: 1.0
971C
972C     Purpose: To set up the (occ,vir) part of the Coupled
973C              Cluster one electron density.
974C
975C     NOTE: This block is stored transposed, i.e. like a t1-amplitude!
976C
977#include "implicit.h"
978      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
979      DIMENSION DIA(*), T2AMT(*), Z1AM(*)
980#include "priunit.h"
981#include "ccorb.h"
982#include "ccsdsym.h"
983C
984      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
985C
986      CALL QENTER('DIAGEN')
987C
988      ISYMDK = 1
989      ISYMAI = ISYMDK
990C
991C--------------------------
992C     Set up density block.
993C--------------------------
994C
995      DO 100 NAI = 1,NT1AM(ISYMAI)
996C
997         DO 110 NDK = 1,NT1AM(ISYMDK)
998C
999            NDKAI = IT2AM(ISYMDK,ISYMAI) + INDEX(NDK,NAI)
1000C
1001            DIA(NAI) = DIA(NAI) + T2AMT(NDKAI)*Z1AM(NDK)
1002C
1003  110    CONTINUE
1004  100 CONTINUE
1005C
1006      CALL QEXIT('DIAGEN')
1007C
1008      RETURN
1009      END
1010C  /* Deck dijk01 */
1011      SUBROUTINE DIJK01(D2IJK,ISYIJK,XLAMDH,ISYMLH,IDEL,ISYMD)
1012C
1013C     Written by Asger Halkier 29/1 - 1996
1014C
1015C     Version: 1.0
1016C
1017C     Purpose: To calculate the first contribution to D2IJK!
1018C
1019#include "implicit.h"
1020      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, FOUR = 4.0D0)
1021      DIMENSION D2IJK(*), XLAMDH(*)
1022#include "priunit.h"
1023#include "ccorb.h"
1024#include "ccsdsym.h"
1025C
1026      CALL QENTER('DIJK01')
1027C
1028      ISYMK  = MULD2H(ISYMD,ISYMLH)
1029      ISYMIJ = MULD2H(ISYIJK,ISYMK)
1030C
1031      IF (ISYMIJ.NE.1) CALL QUIT('Symmetry mismatch in DIJK01')
1032C
1033C-----------------------------------------------
1034C     Calculate special adresses and add result.
1035C-----------------------------------------------
1036C
1037      DO 100 K = 1,NRHF(ISYMK)
1038C
1039         DO 110 ISYMJ = 1,NSYM
1040C
1041            DO 120 J = 1,NRHF(ISYMJ)
1042C
1043               NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
1044     *              + IMATIJ(ISYMJ,ISYMJ) + NRHF(ISYMJ)*(J - 1) + J
1045               NDEL = IGLMRH(ISYMD,ISYMK) + NBAS(ISYMD)*(K - 1)
1046     *              + IDEL - IBAS(ISYMD)
1047C
1048               D2IJK(NIJK) = D2IJK(NIJK) + FOUR*XLAMDH(NDEL)
1049C
1050  120       CONTINUE
1051  110    CONTINUE
1052  100 CONTINUE
1053C
1054      CALL QEXIT('DIJK01')
1055C
1056      RETURN
1057      END
1058C  /* Deck dijk02 */
1059      SUBROUTINE DIJK02(D2IJK,ISYIJK,XLAMDH,ISYMLH,IDEL,ISYMD)
1060C
1061C     Written by Asger Halkier 29/1 - 1996
1062C
1063C     Version: 1.0
1064C
1065C     Purpose: To calculate the second contribution to D2IJK!
1066C
1067#include "implicit.h"
1068      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
1069      DIMENSION D2IJK(*), XLAMDH(*)
1070#include "priunit.h"
1071#include "ccorb.h"
1072#include "ccsdsym.h"
1073C
1074      CALL QENTER('DIJK02')
1075C
1076      ISYMI  = MULD2H(ISYMD,ISYMLH)
1077C
1078      IF (ISYMI.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK02')
1079C
1080C-----------------------------------------------
1081C     Calculate special adresses and add result.
1082C-----------------------------------------------
1083C
1084      DO 100 ISYMK = 1,NSYM
1085C
1086         ISYMIJ = MULD2H(ISYMI,ISYMK)
1087C
1088         DO 110 K = 1,NRHF(ISYMK)
1089C
1090            DO 120 I = 1,NRHF(ISYMI)
1091C
1092               NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
1093     *              + IMATIJ(ISYMI,ISYMK) + NRHF(ISYMI)*(K - 1) + I
1094               NDEL = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1)
1095     *              + IDEL - IBAS(ISYMD)
1096C
1097               D2IJK(NIJK) = D2IJK(NIJK) - TWO*XLAMDH(NDEL)
1098C
1099  120       CONTINUE
1100  110    CONTINUE
1101  100 CONTINUE
1102C
1103      CALL QEXIT('DIJK02')
1104C
1105      RETURN
1106      END
1107C  /* Deck cc_d2pao */
1108      SUBROUTINE CC_D2PAO(D2AOIJ,D2AOAI,D2AOAB,D2AOIA,XLAMDP,D2IJK,
1109     *                    D2AIJ,D2IAJ,D2ABI,G,ISYMG,ISYDEN)
1110C
1111C     Written by Asger Halkier 30/1 - 1996
1112C
1113C     Version: 1.0
1114C
1115C     Purpose: To backtransform the four 2 electron densities
1116C              d(pq,k;del) to AO basis for a given gamma d(pq;gam;del)!
1117C
1118#include "implicit.h"
1119      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
1120      DIMENSION D2AOIJ(*), D2AOAI(*), D2AOAB(*), D2AOIA(*), XLAMDP(*)
1121      DIMENSION D2IJK(*), D2AIJ(*), D2IAJ(*), D2ABI(*)
1122#include "priunit.h"
1123#include "ccorb.h"
1124#include "ccsdsym.h"
1125C
1126      CALL QENTER('CC_D2PAO')
1127C
1128      ISYMK  = ISYMG
1129      ISYMPQ = MULD2H(ISYMG,ISYDEN)
1130C
1131      KOFFGK = IGLMRH(ISYMG,ISYMK) + G
1132C
1133C-------------------------------
1134C     Transform (occ,occ) block.
1135C-------------------------------
1136C
1137      KOFF1  = IMAIJK(ISYMPQ,ISYMK) + 1
1138C
1139      NTOTIJ = MAX(NMATIJ(ISYMPQ),1)
1140C
1141      CALL DGEMV('N',NMATIJ(ISYMPQ),NRHF(ISYMK),ONE,D2IJK(KOFF1),
1142     *           NTOTIJ,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOIJ,1)
1143C
1144C----------------------------------------------
1145C     Transform (vir,occ) and (occ,vir) blocks.
1146C----------------------------------------------
1147C
1148      KOFF2  = IT2BCD(ISYMPQ,ISYMK) + 1
1149C
1150      NTOTAI = MAX(NT1AM(ISYMPQ),1)
1151C
1152      CALL DGEMV('N',NT1AM(ISYMPQ),NRHF(ISYMK),ONE,D2AIJ(KOFF2),
1153     *           NTOTAI,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOAI,1)
1154C
1155      CALL DGEMV('N',NT1AM(ISYMPQ),NRHF(ISYMK),ONE,D2IAJ(KOFF2),
1156     *           NTOTAI,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOIA,1)
1157C
1158C-------------------------------
1159C     Transform (vir,vir) block.
1160C-------------------------------
1161C
1162      KOFF3  = ICKASR(ISYMPQ,ISYMK) + 1
1163C
1164      NTOTAB = MAX(NMATAB(ISYMPQ),1)
1165C
1166      CALL DGEMV('N',NMATAB(ISYMPQ),NRHF(ISYMK),ONE,D2ABI(KOFF3),
1167     *           NTOTAB,XLAMDP(KOFFGK),NBAS(ISYMG),ZERO,D2AOAB,1)
1168C
1169      CALL QEXIT('CC_D2PAO')
1170C
1171      RETURN
1172      END
1173C  /* Deck ccsd_d2en */
1174      SUBROUTINE CCSD_D2EN(ECCSD,XINT,XLAMDH,XLAMDP,D2AOIJ,D2AOAI,
1175     *                     D2AOAB,D2AOIA,G,ISYMG,ISYMD,WORK,LWORK)
1176C
1177C     Written by Asger Halkier 30/1 - 1996
1178C
1179C     Version: 1.0
1180C
1181C     Purpose: To calculate the contribution from the 2 electron
1182C              Coupled Cluster density to the ccsd energy!
1183C
1184#include "implicit.h"
1185      PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0)
1186      DIMENSION ECCSD(*), XINT(*), XLAMDH(*), XLAMDP(*)
1187      DIMENSION D2AOIJ(*), D2AOAI(*), D2AOAB(*), D2AOIA(*), WORK(LWORK)
1188#include "priunit.h"
1189#include "ccorb.h"
1190#include "ccsdsym.h"
1191C
1192      CALL QENTER('CCSD_D2EN')
1193C
1194      ISYINT = ISYMD
1195      ISYDEN = ISYMD
1196      ISYMPQ = MULD2H(ISYMG,ISYDEN)
1197      ISALBE = ISYMPQ
1198C
1199C-------------------------------
1200C     Work space allocation one.
1201C-------------------------------
1202C
1203      KAOINT = 1
1204      KEND1  = KAOINT + N2BST(ISALBE)
1205      LWRK1  = LWORK  - KEND1
1206C
1207      IF (LWRK1. LT. 0) THEN
1208         WRITE(LUPRI,*) 'Needed:', KEND1, 'Available:', LWORK
1209         CALL QUIT('Insufficient memory for allocation in CCSD_D2EN')
1210      ENDIF
1211C
1212C-------------------------------------
1213C     Square up integral distribution.
1214C-------------------------------------
1215C
1216      KOFF1 = IDSAOG(ISYMG,ISYINT) + NNBST(ISALBE)*(G - 1) + 1
1217C
1218      CALL CCSD_SYMSQ(XINT(KOFF1),ISALBE,WORK(KAOINT))
1219C
1220C-------------------------------
1221C     Work space allocation two.
1222C-------------------------------
1223C
1224      LINTMO = MAX(NT1AM(ISYMPQ),NMATIJ(ISYMPQ),NMATAB(ISYMPQ))
1225C
1226      KINTMO = KEND1
1227      KEND2  = KINTMO + LINTMO
1228      LWRK2  = LWORK  - KEND2
1229C
1230      IF (LWRK2. LT. 0) THEN
1231         WRITE(LUPRI,*) 'Needed:', KEND2, 'Available:', LWORK
1232         CALL QUIT('Insufficient memory for allocation in CCSD_D2EN')
1233      ENDIF
1234C
1235C--------------------------------------
1236C     Calculate (occ,occ) contribution.
1237C--------------------------------------
1238C
1239      CALL DZERO(WORK(KINTMO),NMATIJ(ISYMPQ))
1240C
1241      DO 100 ISYMAL = 1,NSYM
1242C
1243         ISYMBE = MULD2H(ISYMAL,ISALBE)
1244         ISYMJ  = ISYMBE
1245         ISYMI  = ISYMAL
1246C
1247         LNEED  = LWORK - KEND2 - NBAS(ISYMAL)*NRHF(ISYMJ)
1248C
1249         IF (LNEED .LT. 0) THEN
1250            WRITE(LUPRI,*) 'Work exceeded by:', -LNEED
1251            CALL QUIT('Insufficient work space in CCSD_D2EN')
1252         ENDIF
1253C
1254         CALL DZERO(WORK(KEND2),NBAS(ISYMAL)*NRHF(ISYMJ))
1255C
1256         KOFF2  = KAOINT + IAODIS(ISYMAL,ISYMBE)
1257         KOFF3  = IGLMRH(ISYMBE,ISYMJ) + 1
1258C
1259         NTOTAL = MAX(NBAS(ISYMAL),1)
1260         NTOTBE = MAX(NBAS(ISYMBE),1)
1261C
1262         CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMJ),NBAS(ISYMBE),ONE,
1263     *              WORK(KOFF2),NTOTAL,XLAMDH(KOFF3),NTOTBE,ZERO,
1264     *              WORK(KEND2),NTOTAL)
1265C
1266         KOFF4  = IGLMRH(ISYMAL,ISYMI) + 1
1267         KOFF5  = KINTMO + IMATIJ(ISYMI,ISYMJ)
1268C
1269         NTOTAL = MAX(NBAS(ISYMAL),1)
1270         NTOTI  = MAX(NRHF(ISYMI),1)
1271C
1272         CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NBAS(ISYMAL),HALF,
1273     *              XLAMDP(KOFF4),NTOTAL,WORK(KEND2),NTOTAL,ZERO,
1274     *              WORK(KOFF5),NTOTI)
1275C
1276  100 CONTINUE
1277C
1278      ECCSD(1) = ECCSD(1) + DDOT(NMATIJ(ISYMPQ),WORK(KINTMO),1,
1279     *                           D2AOIJ,1)
1280C
1281C--------------------------------------
1282C     Calculate (vir,vir) contribution.
1283C--------------------------------------
1284C
1285      CALL DZERO(WORK(KINTMO),NMATAB(ISYMPQ))
1286C
1287      DO 110 ISYMAL = 1,NSYM
1288C
1289         ISYMBE = MULD2H(ISYMAL,ISALBE)
1290         ISYMA  = ISYMAL
1291         ISYMB  = ISYMBE
1292C
1293         LNEED  = LWORK - KEND2 - NBAS(ISYMAL)*NVIR(ISYMB)
1294C
1295         IF (LNEED .LT. 0) THEN
1296            WRITE(LUPRI,*) 'Work exceeded by:', -LNEED
1297            CALL QUIT('Insufficient work space in CCSD_D2EN')
1298         ENDIF
1299C
1300         CALL DZERO(WORK(KEND2),NBAS(ISYMAL)*NVIR(ISYMB))
1301C
1302         KOFF6  = KAOINT + IAODIS(ISYMAL,ISYMBE)
1303         KOFF7  = IGLMVI(ISYMBE,ISYMB) + 1
1304C
1305         NTOTAL = MAX(NBAS(ISYMAL),1)
1306         NTOTBE = MAX(NBAS(ISYMBE),1)
1307C
1308         CALL DGEMM('N','N',NBAS(ISYMAL),NVIR(ISYMB),NBAS(ISYMBE),ONE,
1309     *              WORK(KOFF6),NTOTAL,XLAMDH(KOFF7),NTOTBE,ZERO,
1310     *              WORK(KEND2),NTOTAL)
1311C
1312         KOFF8  = IGLMVI(ISYMAL,ISYMA) + 1
1313         KOFF9  = KINTMO + IMATAB(ISYMA,ISYMB)
1314C
1315         NTOTAL = MAX(NBAS(ISYMAL),1)
1316         NTOTA  = MAX(NVIR(ISYMA),1)
1317C
1318         CALL DGEMM('T','N',NVIR(ISYMA),NVIR(ISYMB),NBAS(ISYMAL),HALF,
1319     *              XLAMDP(KOFF8),NTOTAL,WORK(KEND2),NTOTAL,ZERO,
1320     *              WORK(KOFF9),NTOTA)
1321C
1322  110 CONTINUE
1323C
1324      ECCSD(1) = ECCSD(1) + DDOT(NMATAB(ISYMPQ),WORK(KINTMO),1,
1325     *                           D2AOAB,1)
1326C
1327C------------------------------------------
1328C     Calculate the (vir,occ) contribution.
1329C------------------------------------------
1330C
1331      CALL DZERO(WORK(KINTMO),NT1AM(ISYMPQ))
1332C
1333      DO 120 ISYMAL = 1,NSYM
1334C
1335         ISYMBE = MULD2H(ISYMAL,ISALBE)
1336         ISYMI  = ISYMBE
1337         ISYMA  = ISYMAL
1338C
1339         LNEED  = LWORK - KEND2 - NBAS(ISYMAL)*NRHF(ISYMI)
1340C
1341         IF (LNEED .LT. 0) THEN
1342            WRITE(LUPRI,*) 'Work exceeded by:', -LNEED
1343            CALL QUIT('Insufficient work space in CCSD_D2EN')
1344         ENDIF
1345C
1346         CALL DZERO(WORK(KEND2),NBAS(ISYMAL)*NRHF(ISYMI))
1347C
1348         KOFF10 = KAOINT + IAODIS(ISYMAL,ISYMBE)
1349         KOFF11 = IGLMRH(ISYMBE,ISYMI) + 1
1350C
1351         NTOTAL = MAX(NBAS(ISYMAL),1)
1352         NTOTBE = MAX(NBAS(ISYMBE),1)
1353C
1354         CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMI),NBAS(ISYMBE),ONE,
1355     *              WORK(KOFF10),NTOTAL,XLAMDH(KOFF11),NTOTBE,ZERO,
1356     *              WORK(KEND2),NTOTAL)
1357C
1358         KOFF12 = IGLMVI(ISYMAL,ISYMA) + 1
1359         KOFF13 = KINTMO + IT1AM(ISYMA,ISYMI)
1360C
1361         NTOTAL = MAX(NBAS(ISYMAL),1)
1362         NTOTA  = MAX(NVIR(ISYMA),1)
1363C
1364         CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMAL),HALF,
1365     *              XLAMDP(KOFF12),NTOTAL,WORK(KEND2),NTOTAL,ZERO,
1366     *              WORK(KOFF13),NTOTA)
1367C
1368  120 CONTINUE
1369C
1370      ECCSD(1) = ECCSD(1) + DDOT(NT1AM(ISYMPQ),WORK(KINTMO),1,
1371     *                           D2AOAI,1)
1372C
1373C------------------------------------------
1374C     Calculate the (occ,vir) contribution.
1375C------------------------------------------
1376C
1377      CALL DZERO(WORK(KINTMO),NT1AM(ISYMPQ))
1378C
1379      DO 130 ISYMAL = 1,NSYM
1380C
1381         ISYMBE = MULD2H(ISYMAL,ISALBE)
1382         ISYMI  = ISYMAL
1383         ISYMA  = ISYMBE
1384C
1385         LNEED = LWORK - KEND2 - NBAS(ISYMBE)*NRHF(ISYMI)
1386C
1387         IF (LNEED .LT. 0) THEN
1388            WRITE(LUPRI,*) 'Work exceeded by:', -LNEED
1389            CALL QUIT('Insufficient work space in CCSD_D2EN')
1390         ENDIF
1391C
1392         CALL DZERO(WORK(KEND2),NBAS(ISYMBE)*NRHF(ISYMI))
1393C
1394         KOFF14 = KAOINT + IAODIS(ISYMAL,ISYMBE)
1395         KOFF15 = IGLMRH(ISYMAL,ISYMI) + 1
1396C
1397         NTOTAL = MAX(NBAS(ISYMAL),1)
1398         NTOTBE = MAX(NBAS(ISYMBE),1)
1399C
1400         CALL DGEMM('T','N',NBAS(ISYMBE),NRHF(ISYMI),NBAS(ISYMAL),ONE,
1401     *              WORK(KOFF14),NTOTAL,XLAMDP(KOFF15),NTOTAL,ZERO,
1402     *              WORK(KEND2),NTOTBE)
1403C
1404         KOFF16 = IGLMVI(ISYMBE,ISYMA) + 1
1405         KOFF17 = KINTMO + IT1AM(ISYMA,ISYMI)
1406C
1407         NTOTBE = MAX(NBAS(ISYMBE),1)
1408         NTOTA  = MAX(NVIR(ISYMA),1)
1409C
1410         CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMBE),HALF,
1411     *              XLAMDH(KOFF16),NTOTBE,WORK(KEND2),NTOTBE,ZERO,
1412     *              WORK(KOFF17),NTOTA)
1413C
1414  130 CONTINUE
1415C
1416      ECCSD(1) = ECCSD(1) + DDOT(NT1AM(ISYMPQ),WORK(KINTMO),1,
1417     *                           D2AOIA,1)
1418C
1419      CALL QEXIT('CCSD_D2EN')
1420C
1421      RETURN
1422      END
1423C  /* Deck dijk11 */
1424      SUBROUTINE DIJK11(D2IJK,ISYIJK,XLAMDH,ISYMLH,DIA,ISYD1,
1425     *                  WORK,LWORK,IDEL,ISYMD)
1426C
1427C     Written by Asger Halkier 2/2 - 1996
1428C
1429C     Version: 1.0
1430C
1431C     Purpose: To calculate the first two contribution to D2IJK
1432C              from projection against singles!
1433C
1434#include "implicit.h"
1435      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
1436      DIMENSION D2IJK(*), XLAMDH(*), DIA(*), WORK(LWORK)
1437#include "priunit.h"
1438#include "ccorb.h"
1439#include "ccsdsym.h"
1440C
1441      CALL QENTER('DIJK11')
1442C
1443      ISYMA  = MULD2H(ISYMD,ISYMLH)
1444      ISYMK  = MULD2H(ISYMA,ISYD1)
1445      ISYMIJ = MULD2H(ISYMK,ISYIJK)
1446C
1447      IF (ISYMK.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK11 (1)')
1448C
1449      IF (LWORK .LT. NRHF(ISYMK)) THEN
1450         WRITE(LUPRI,*) 'Available:', LWORK, 'NEEDED:', NRHF(ISYMK)
1451         CALL QUIT('Insufficient space for allocation in DIJK11')
1452      ENDIF
1453C
1454      CALL DZERO(WORK,NRHF(ISYMK))
1455C
1456C-------------------------------------------
1457C     Contract D(ia) with lamda hole matrix.
1458C-------------------------------------------
1459C
1460      KOFF1 = IT1AM(ISYMA,ISYMK) + 1
1461      KOFF2 = IGLMVI(ISYMD,ISYMA) + IDEL - IBAS(ISYMD)
1462C
1463      NTOTA = MAX(NVIR(ISYMA),1)
1464C
1465      CALL DGEMV('T',NVIR(ISYMA),NRHF(ISYMK),ONE,DIA(KOFF1),NTOTA,
1466     *           XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1)
1467C
1468C---------------------------------------
1469C     Calculate the DIJK11 contribution.
1470C---------------------------------------
1471C
1472      DO 100 K = 1,NRHF(ISYMK)
1473C
1474         DO 110 ISYMI = 1,NSYM
1475C
1476            ISYMJ = ISYMI
1477C
1478            DO 120 I = 1,NRHF(ISYMI)
1479C
1480               J    = I
1481               NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
1482     *              + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I
1483C
1484               D2IJK(NIJK) = D2IJK(NIJK) + TWO*WORK(K)
1485C
1486  120       CONTINUE
1487  110    CONTINUE
1488  100 CONTINUE
1489C
1490C---------------------------------------
1491C     Calculate the DIJK12 contribution.
1492C---------------------------------------
1493C
1494      ISYMI  = MULD2H(ISYMA,ISYD1)
1495      ISYMJK = MULD2H(ISYMI,ISYIJK)
1496
1497      IF (ISYMI.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK11 (2)')
1498C
1499      DO 130 ISYMK = 1,NSYM
1500C
1501         ISYMJ  = ISYMK
1502         ISYMIJ = MULD2H(ISYMI,ISYMJ)
1503C
1504         DO 140 K = 1,NRHF(ISYMK)
1505C
1506            J = K
1507C
1508            DO 150 I = 1,NRHF(ISYMI)
1509C
1510               NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
1511     *              + IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I
1512C
1513               D2IJK(NIJK) = D2IJK(NIJK) - WORK(I)
1514C
1515  150       CONTINUE
1516  140    CONTINUE
1517  130 CONTINUE
1518C
1519      CALL QEXIT('DIJK11')
1520C
1521      RETURN
1522      END
1523C  /* Deck dijk13 */
1524      SUBROUTINE DIJK13(D2IJK,ISYIJK,DENAI,ISYD1,T2TIN,ISYTIN)
1525C
1526C     Written by Asger Halkier 2/2 - 1996
1527C
1528C     Version: 1.0
1529C
1530C     Purpose: To calculate the third contribution to D2IJK
1531C              from projection against singles!
1532C
1533#include "implicit.h"
1534      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
1535      DIMENSION D2IJK(*), DENAI(*), T2TIN(*)
1536#include "ccorb.h"
1537#include "ccsdsym.h"
1538C
1539      CALL QENTER('DIJK13')
1540C
1541      IF (MULD2H(ISYTIN,ISYD1).NE.ISYIJK)
1542     &  CALL QUIT('Symmetry mismatch in DIJK13')
1543C
1544C----------------------------
1545C     Calculate contribution.
1546C----------------------------
1547C
1548      DO 100 ISYMK = 1,NSYM
1549C
1550         ISYMEI = MULD2H(ISYMK,ISYTIN)
1551         ISYMIJ = MULD2H(ISYMK,ISYIJK)
1552C
1553         DO 110 K = 1,NRHF(ISYMK)
1554C
1555            DO 120 ISYMI = 1,NSYM
1556C
1557               ISYME = MULD2H(ISYMI,ISYMEI)
1558               ISYMJ = MULD2H(ISYME,ISYD1)
1559C
1560               KOFF1 = IT2BCD(ISYMEI,ISYMK) + NT1AM(ISYMEI)*(K - 1)
1561     *               + IT1AM(ISYME,ISYMI) + 1
1562               KOFF2 = IT1AM(ISYME,ISYMJ) + 1
1563               KOFF3 = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
1564     *               + IMATIJ(ISYMI,ISYMJ) + 1
1565C
1566               NTOTE = MAX(NVIR(ISYME),1)
1567               NTOTI = MAX(NRHF(ISYMI),1)
1568C
1569               CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYME),
1570     *                    -ONE,T2TIN(KOFF1),NTOTE,DENAI(KOFF2),NTOTE,
1571     *                    ONE,D2IJK(KOFF3),NTOTI)
1572C
1573  120       CONTINUE
1574  110    CONTINUE
1575  100 CONTINUE
1576C
1577      CALL QEXIT('DIJK13')
1578C
1579      RETURN
1580      END
1581C  /* Deck daij11 */
1582      SUBROUTINE DAIJ11(D2AIJ,ISYAIJ,DAI,ISYD1,XLAMDH,ISYMLH,IDEL,ISYMD)
1583C
1584C     Written by Asger Halkier 2/2 - 1996
1585C
1586C     Version: 1.0
1587C
1588C     Purpose: To calculate the first contribution to D2AIJ
1589C              from projection against singles!
1590C
1591#include "implicit.h"
1592      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
1593      DIMENSION D2AIJ(*), XLAMDH(*), DAI(*)
1594#include "ccorb.h"
1595#include "ccsdsym.h"
1596C
1597      CALL QENTER('DAIJ11')
1598C
1599      ISYMJ  = MULD2H(ISYMD,ISYMLH)
1600      ISYMAI = ISYD1
1601      IF (MULD2H(ISYMJ,ISYMAI).NE.ISYAIJ) THEN
1602        CALL QUIT('Symmetry mismatch in DAIJ11')
1603      END IF
1604C
1605C----------------------------
1606C     Calculate contribution.
1607C----------------------------
1608C
1609      DO 100 J = 1,NRHF(ISYMJ)
1610C
1611         NDJ = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J - 1) +
1612     *         IDEL - IBAS(ISYMD)
1613C
1614         FACT = TWO*XLAMDH(NDJ)
1615C
1616         KOFF1 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) + 1
1617C
1618         CALL DAXPY(NT1AM(ISYMAI),FACT,DAI,1,D2AIJ(KOFF1),1)
1619C
1620  100 CONTINUE
1621C
1622      CALL QEXIT('DAIJ11')
1623C
1624      RETURN
1625      END
1626C  /* Deck daij12 */
1627      SUBROUTINE DAIJ12(D2AIJ,ISYAIJ,DAI,ISYD1,XLAMDH,ISYMLH,
1628     *                  WORK,LWORK,IDEL,ISYMD)
1629C
1630C     Written by Asger Halkier 2/2 - 1996
1631C
1632C     Version: 1.0
1633C
1634C     Purpose: To calculate the second contribution to D2AIJ
1635C              from projection against singles!
1636C
1637#include "implicit.h"
1638      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
1639      DIMENSION D2AIJ(*), XLAMDH(*), DAI(*), WORK(LWORK)
1640#include "priunit.h"
1641#include "ccorb.h"
1642#include "ccsdsym.h"
1643C
1644      CALL QENTER('DAIJ12')
1645C
1646      ISYMK  = MULD2H(ISYMD,ISYMLH)
1647      ISYMA  = MULD2H(ISYMK,ISYD1)
1648C
1649      IF (ISYMA.NE.ISYAIJ) THEN
1650        CALL QUIT('Symmetry mismatch in DAIJ12')
1651      END IF
1652C
1653      IF (LWORK .LT. NVIR(ISYMA)) THEN
1654         WRITE(LUPRI,*) 'Available:', LWORK, 'NEEDED:', NVIR(ISYMA)
1655         CALL QUIT('Insufficient space for allocation in DAIJ12')
1656      ENDIF
1657C
1658      CALL DZERO(WORK,NVIR(ISYMA))
1659C
1660C----------------------------------------------------------
1661C     Calculate contraction of D(ai) and lamda hole matrix.
1662C----------------------------------------------------------
1663C
1664      KOFF1 = IT1AM(ISYMA,ISYMK) + 1
1665      KOFF2 = IGLMRH(ISYMD,ISYMK) + IDEl - IBAS(ISYMD)
1666C
1667      NTOTA = MAX(NVIR(ISYMA),1)
1668C
1669      CALL DGEMV('N',NVIR(ISYMA),NRHF(ISYMK),ONE,DAI(KOFF1),NTOTA,
1670     *           XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1)
1671C
1672C----------------------------
1673C     Calculate contribution.
1674C----------------------------
1675C
1676      DO 100 ISYMI = 1,NSYM
1677C
1678         ISYMJ  = ISYMI
1679         ISYMAI = MULD2H(ISYMI,ISYMA)
1680C
1681         DO 110 I = 1,NRHF(ISYMI)
1682C
1683            J    = I
1684            NAIJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1)
1685     *           + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1
1686C
1687            CALL DAXPY(NVIR(ISYMA),-ONE,WORK,1,D2AIJ(NAIJ),1)
1688C
1689  110    CONTINUE
1690  100 CONTINUE
1691C
1692      CALL QEXIT('DAIJ12')
1693C
1694      RETURN
1695      END
1696C  /* Deck diaj11 */
1697      SUBROUTINE DIAJ11(D2IAJ,ISYIAJ,DIA,ISYD1,XLAMDH,ISYMLH,IDEL,ISYMD)
1698C
1699C     Written by Asger Halkier 4/2 - 1996
1700C
1701C     Version: 1.0
1702C
1703C     Purpose: To calculate the first and second contribution to D2IAJ
1704C              from projection against singles!
1705C
1706#include "implicit.h"
1707      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
1708      DIMENSION D2IAJ(*), XLAMDH(*), DIA(*)
1709#include "ccorb.h"
1710#include "ccsdsym.h"
1711C
1712      CALL QENTER('DIAJ11')
1713C
1714      ISYMJ  = MULD2H(ISYMD,ISYMLH)
1715      ISYMAI = ISYD1
1716
1717      IF (MULD2H(ISYMJ,ISYMAI).NE.ISYIAJ)
1718     *   CALL QUIT('Symmetry mismatch in DIAJ11. (1)')
1719C
1720C--------------------------------------
1721C     Calculate the first contribution.
1722C--------------------------------------
1723C
1724      DO 100 J = 1,NRHF(ISYMJ)
1725C
1726         NDJ = IGLMRH(ISYMD,ISYMJ) + NBAS(ISYMD)*(J - 1) +
1727     *         IDEL - IBAS(ISYMD)
1728C
1729         FACT = TWO*XLAMDH(NDJ)
1730C
1731         KOFF = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) + 1
1732C
1733         CALL DAXPY(NT1AM(ISYMAI),FACT,DIA,1,D2IAJ(KOFF),1)
1734C
1735  100 CONTINUE
1736C
1737      ISYMI  = MULD2H(ISYMD,ISYMLH)
1738      ISYMAJ = ISYD1
1739
1740      IF (MULD2H(ISYMI,ISYMAJ).NE.ISYIAJ)
1741     *   CALL QUIT('Symmetry mismatch in DIAJ11. (2)')
1742C
1743C---------------------------------------
1744C     Calculate the second contribution.
1745C---------------------------------------
1746C
1747      DO 110 ISYMJ = 1,NSYM
1748C
1749         ISYMA  = MULD2H(ISYMJ,ISYMAJ)
1750         ISYMAI = MULD2H(ISYMA,ISYMI)
1751C
1752         DO 120 J = 1,NRHF(ISYMJ)
1753C
1754            DO 130 I = 1,NRHF(ISYMI)
1755C
1756               NDI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I- 1) + IDEL
1757     *             - IBAS(ISYMD)
1758               NAJ  = IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J - 1) + 1
1759               NIAJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1)
1760     *              + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + 1
1761C
1762               CALL DAXPY(NVIR(ISYMA),-XLAMDH(NDI),DIA(NAJ),1,
1763     *                    D2IAJ(NIAJ),1)
1764C
1765  130       CONTINUE
1766  120    CONTINUE
1767  110 CONTINUE
1768C
1769      CALL QEXIT('DIAJ11')
1770C
1771      RETURN
1772      END
1773C  /* Deck dabi11 */
1774      SUBROUTINE DABI11(D2ABI,ISYABI,DAI,ISYD1,T2TDEL,ISYMTI)
1775C
1776C     Written by Asger Halkier 4/2 - 1996
1777C
1778C     Version: 1.0
1779C
1780C     Purpose: To calculate the one and only contribution to D2ABI
1781C              from projection against singles!
1782C
1783#include "implicit.h"
1784      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
1785      DIMENSION D2ABI(*), DAI(*), T2TDEL(*)
1786#include "ccorb.h"
1787#include "ccsdsym.h"
1788C
1789      CALL QENTER('DABI11')
1790C
1791      IF (MULD2H(ISYMTI,ISYD1).NE.ISYABI)
1792     &  CALL QUIT('Symmetry in DABI11')
1793C
1794C----------------------------
1795C     Calculate contribution.
1796C----------------------------
1797C
1798      DO 100 ISYMI = 1,NSYM
1799C
1800         ISYMBM = MULD2H(ISYMI,ISYMTI)
1801         ISYMAB = MULD2H(ISYMI,ISYABI)
1802C
1803         DO 110 I = 1,NRHF(ISYMI)
1804C
1805            DO 120 ISYMA = 1,NSYM
1806C
1807               ISYMM = MULD2H(ISYMA,ISYD1)
1808               ISYMB = MULD2H(ISYMM,ISYMBM)
1809C
1810               KOFF1 = IT1AM(ISYMA,ISYMM)   + 1
1811               KOFF2 = IT2BCD(ISYMBM,ISYMI) + NT1AM(ISYMBM)*(I - 1)
1812     *               + IT1AM(ISYMB,ISYMM)   + 1
1813               KOFF3 = ICKASR(ISYMAB,ISYMI) + NMATAB(ISYMAB)*(I - 1)
1814     *               + IMATAB(ISYMA,ISYMB)  + 1
1815C
1816               NTOTA = MAX(NVIR(ISYMA),1)
1817               NTOTB = MAX(NVIR(ISYMB),1)
1818C
1819               CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMM),
1820     *                    ONE,DAI(KOFF1),NTOTA,T2TDEL(KOFF2),NTOTB,ONE,
1821     *                    D2ABI(KOFF3),NTOTA)
1822C
1823  120       CONTINUE
1824  110    CONTINUE
1825  100 CONTINUE
1826C
1827      CALL QEXIT('DABI11')
1828C
1829      RETURN
1830      END
1831C  /* Deck dijg11 */
1832      SUBROUTINE DIJG11(D2AOIJ,XLAMDH,Z1AO,IDEL,ISYMD,G,ISYMG)
1833C
1834C     Written by Asger Halkier 4/2 - 1996
1835C
1836C     Version: 1.0
1837C
1838C     Purpose: To calculate the first and second contribution to D2AOIJ
1839C              from projection against singles!
1840C
1841#include "implicit.h"
1842      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
1843      DIMENSION D2AOIJ(*), XLAMDH(*), Z1AO(*)
1844#include "ccorb.h"
1845#include "ccsdsym.h"
1846C
1847      CALL QENTER('DIJG11')
1848C
1849C--------------------------------------
1850C     Calculate the first contribution.
1851C--------------------------------------
1852C
1853      IF (ISYMG .EQ. ISYMD) THEN
1854C
1855         ISYMK = ISYMG
1856C
1857         NDK = IGLMRH(ISYMD,ISYMK) + IDEL - IBAS(ISYMD)
1858         NGK = IGLMRH(ISYMG,ISYMK) + G
1859C
1860         FACT = DDOT(NRHF(ISYMK),Z1AO(NGK),NBAS(ISYMG),XLAMDH(NDK),
1861     *               NBAS(ISYMD))
1862C
1863         DO 100 ISYMI = 1,NSYM
1864C
1865            ISYMJ = ISYMI
1866C
1867            DO 110 I = 1,NRHF(ISYMI)
1868C
1869               J   = I
1870               NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I
1871C
1872               D2AOIJ(NIJ) = D2AOIJ(NIJ) + TWO*FACT
1873C
1874  110       CONTINUE
1875  100    CONTINUE
1876C
1877      ENDIF
1878C
1879C-----------------------------------
1880C     Calculate second contribution.
1881C-----------------------------------
1882C
1883      ISYMI = ISYMD
1884      ISYMJ = ISYMG
1885C
1886      DO 120 J = 1,NRHF(ISYMJ)
1887C
1888         NGJ = IGLMRH(ISYMG,ISYMJ) + NBAS(ISYMG)*(J - 1) + G
1889C
1890         DO 130 I = 1,NRHF(ISYMI)
1891C
1892            NDI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1) + IDEL
1893     *          - IBAS(ISYMD)
1894            NIJ = IMATIJ(ISYMI,ISYMJ) + NRHF(ISYMI)*(J - 1) + I
1895C
1896            D2AOIJ(NIJ) = D2AOIJ(NIJ) - XLAMDH(NDI)*Z1AO(NGJ)
1897C
1898  130    CONTINUE
1899  120 CONTINUE
1900C
1901      CALL QEXIT('DIJG11')
1902C
1903      RETURN
1904      END
1905C  /* Deck diaj13 */
1906      SUBROUTINE DIAJ13(D2IAJ,ISYIAJ,DAI,ISYD1,T2AMT,XLAMDH,ISYMLH,
1907     *                  WORK,LWORK,IDEL,ISYMD)
1908C
1909C     Written by Asger Halkier 5/2 - 1996
1910C
1911C     Version: 1.0
1912C
1913C     Purpose: To calculate the third contribution to D2IAJ from
1914C              projection against singles!
1915C
1916#include "implicit.h"
1917      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
1918      DIMENSION D2IAJ(*), DAI(*), T2AMT(*), XLAMDH(*), WORK(LWORK)
1919#include "priunit.h"
1920#include "ccorb.h"
1921#include "ccsdsym.h"
1922C
1923      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
1924C
1925      CALL QENTER('DIAJ13')
1926C
1927      IF (MULD2H(ISYMLH,ISYMD).NE.MULD2H(ISYIAJ,ISYD1)) THEN
1928        CALL QUIT('Symmetry mismatch in DIAJ13.')
1929      END IF
1930C
1931      ISYMK  = MULD2H(ISYMD,ISYMLH)
1932      ISYME  = MULD2H(ISYMK,ISYD1)
1933C
1934C-----------------------------------------
1935C     Spaghetti goto's if no contribution.
1936C-----------------------------------------
1937C
1938      IF (NRHF(ISYMK) .EQ. 0) GOTO 101
1939      IF (NVIR(ISYME) .EQ. 0) GOTO 101
1940      IF (NBAS(ISYMD) .EQ. 0) GOTO 101
1941C
1942C-------------------------------
1943C     Work space allocation one.
1944C-------------------------------
1945C
1946      KVECE = 1
1947      KEND1 = KVECE + NVIR(ISYME)
1948      LWRK1 = LWORK - KEND1
1949C
1950      IF (LWRK1 .LT. 0) THEN
1951         WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1
1952         CALL QUIT('Insufficient work space for allocation in DIAJ13')
1953      ENDIF
1954C
1955      CALL DZERO(WORK(KVECE),NVIR(ISYME))
1956C
1957C-----------------------------------------------
1958C     Calculate contraction of zeta1 and xlamdh.
1959C-----------------------------------------------
1960C
1961      KOFF1 = IT1AM(ISYME,ISYMK) + 1
1962      KOFF2 = IGLMRH(ISYMD,ISYMK) + IDEL - IBAS(ISYMD)
1963C
1964      NTOTE = MAX(NVIR(ISYME),1)
1965C
1966      CALL DGEMV('N',NVIR(ISYME),NRHF(ISYMK),ONE,DAI(KOFF1),NTOTE,
1967     *           XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK(KVECE),1)
1968C
1969      DO 100 ISYMJ = 1,NSYM
1970C
1971         ISYMAI = MULD2H(ISYMJ,ISYIAJ)
1972         ISYMEJ = ISYMAI
1973C
1974C----------------------------------
1975C        Work space allocation two.
1976C----------------------------------
1977C
1978         KT2SUB = KEND1
1979         KEND2  = KT2SUB + NT1AM(ISYMAI)
1980         LWRK2  = LWORK  - KEND2
1981C
1982         IF (LWRK2 .LT. 0) THEN
1983            WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND2
1984            CALL QUIT('Insufficient work space for '//
1985     &           'allocation in DIAJ13')
1986         ENDIF
1987C
1988         DO 110 J = 1,NRHF(ISYMJ)
1989C
1990            DO 120 E = 1,NVIR(ISYME)
1991C
1992               NEJ = IT1AM(ISYME,ISYMJ) + NVIR(ISYME)*(J - 1) + E
1993C
1994C--------------------------------------------------
1995C              Copy t1-submatrix (ai) out of t2amt.
1996C--------------------------------------------------
1997C
1998               DO 130 NAI = 1,NT1AM(ISYMAI)
1999C
2000                  NAIEJ = IT2AM(ISYMAI,ISYMEJ) + INDEX(NAI,NEJ)
2001C
2002                  WORK(KT2SUB + NAI - 1) = T2AMT(NAIEJ)
2003C
2004  130          CONTINUE
2005C
2006C------------------------------------------------------
2007C              Final contraction and storage in result.
2008C------------------------------------------------------
2009C
2010               FACT  = - WORK(KVECE + E - 1)
2011C
2012               KOFF3 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1) + 1
2013C
2014               CALL DAXPY(NT1AM(ISYMAI),FACT,WORK(KT2SUB),1,
2015     *                    D2IAJ(KOFF3),1)
2016C
2017  120       CONTINUE
2018  110    CONTINUE
2019  100 CONTINUE
2020C
2021  101 CONTINUE
2022C
2023      CALL QEXIT('DIAJ13')
2024C
2025      RETURN
2026      END
2027C  /* Deck diag11 */
2028      SUBROUTINE DIAG11(D2AOIA,Z1AO,T2TIN,ISYMD,G,ISYMG)
2029C
2030C     Written by Asger Halkier 7/2 - 1996
2031C
2032C     Version: 1.0
2033C
2034C     Purpose: To calculate the one and only contribution to D2AOIA
2035C              from projection against singles!
2036C
2037#include "implicit.h"
2038      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2039      DIMENSION D2AOIA(*), Z1AO(*), T2TIN(*)
2040#include "ccorb.h"
2041#include "ccsdsym.h"
2042C
2043      CALL QENTER('DIAG11')
2044C
2045      ISYAIM = ISYMD
2046      ISYMM  = ISYMG
2047      ISYMAI = MULD2H(ISYMM,ISYAIM)
2048C
2049C----------------------------
2050C     Calculate contribution.
2051C----------------------------
2052C
2053      KOFF1  = IT2BCD(ISYMAI,ISYMM) + 1
2054      KOFF2  = IGLMRH(ISYMG,ISYMM) + G
2055C
2056      NTOTAI = MAX(NT1AM(ISYMAI),1)
2057C
2058      CALL DGEMV('N',NT1AM(ISYMAI),NRHF(ISYMM),ONE,T2TIN(KOFF1),
2059     *           NTOTAI,Z1AO(KOFF2),NBAS(ISYMG),ONE,D2AOIA,1)
2060C
2061      CALL QEXIT('DIAG11')
2062C
2063      RETURN
2064      END
2065C  /* Deck dijk25 */
2066      SUBROUTINE DIJK25(D2IJK,ISYIJK,XMINT,XLAMDH,ISYMLH,IDEL,ISYMD)
2067C
2068C     Written by Asger Halkier 10/2 - 1996
2069C
2070C     Version: 1.0
2071C
2072C     Purpose: To calculate the fifth contribution to D2IJK
2073C              from projection against doubles!
2074C
2075#include "implicit.h"
2076      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2077      DIMENSION D2IJK(*), XMINT(*), XLAMDH(*)
2078#include "ccorb.h"
2079#include "ccsdsym.h"
2080C
2081      CALL QENTER('DIJK25')
2082C
2083      ISYML  = MULD2H(ISYMD,ISYMLH)
2084      IF (ISYIJK.NE.ISYML) CALL QUIT('Symmetry mismatch in DIJK25')
2085C
2086C--------------------------------
2087C     Calculate the contribution.
2088C--------------------------------
2089C
2090      KOFF1  = I3ORHF(ISYIJK,ISYML) + 1
2091      KOFF2  = IGLMRH(ISYMD,ISYML) + IDEL - IBAS(ISYMD)
2092C
2093      NTOIJK = MAX(NMAIJK(ISYIJK),1)
2094C
2095      CALL DGEMV('N',NMAIJK(ISYIJK),NRHF(ISYML),ONE,XMINT(KOFF1),
2096     *           NTOIJK,XLAMDH(KOFF2),NBAS(ISYMD),ONE,D2IJK,1)
2097C
2098      CALL QEXIT('DIJK25')
2099C
2100      RETURN
2101      END
2102C  /* Deck dijk24 */
2103      SUBROUTINE DIJK24(D2IJK,ISYIJK,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD)
2104C
2105C     Written by Asger Halkier 10/2 - 1996
2106C
2107C     Version: 1.0
2108C
2109C     Purpose: To calculate the fourth contribution to D2IJK
2110C              from projection against doubles!
2111C
2112#include "implicit.h"
2113      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2114      DIMENSION D2IJK(*), XMAT(*), XLAMDH(*)
2115#include "ccorb.h"
2116#include "ccsdsym.h"
2117C
2118      CALL QENTER('DIJK24')
2119C
2120      ISYMI  = MULD2H(ISYMD,ISYMLH)
2121      ISYMKJ = 1
2122C
2123      IF (ISYMI.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK24')
2124C
2125C--------------------------------
2126C     Calculate the contribution.
2127C--------------------------------
2128C
2129      DO 100 ISYMK = 1,NSYM
2130C
2131         ISYMJ  = MULD2H(ISYMK,ISYMKJ)
2132         ISYMIJ = MULD2H(ISYMI,ISYMJ)
2133C
2134         DO 110 K = 1,NRHF(ISYMK)
2135C
2136            DO 120 J = 1,NRHF(ISYMJ)
2137C
2138               NKJ  = IMATIJ(ISYMK,ISYMJ)  + NRHF(ISYMK)*(J - 1) + K
2139               NDEI = IGLMRH(ISYMD,ISYMI) + IDEL - IBAS(ISYMD)
2140               NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
2141     *              + IMATIJ(ISYMI,ISYMJ)  + NRHF(ISYMI)*(J - 1) + 1
2142C
2143               CALL DAXPY(NRHF(ISYMI),XMAT(NKJ),XLAMDH(NDEI),
2144     *                    NBAS(ISYMD),D2IJK(NIJK),1)
2145C
2146  120       CONTINUE
2147  110    CONTINUE
2148  100 CONTINUE
2149C
2150      CALL QEXIT('DIJK24')
2151C
2152      RETURN
2153      END
2154C  /* Deck dijk23 */
2155      SUBROUTINE DIJK23(D2IJK,ISYIJK,XMAT,XLAMDH,ISYMLH,IDEL,ISYMD)
2156C
2157C     Written by Asger Halkier 10/2 - 1996
2158C
2159C     Version: 1.0
2160C
2161C     Purpose: To calculate the third contribution to D2IJK
2162C              from projection against doubles!
2163C
2164#include "implicit.h"
2165      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
2166      DIMENSION D2IJK(*), XMAT(*), XLAMDH(*)
2167#include "ccorb.h"
2168#include "ccsdsym.h"
2169C
2170      CALL QENTER('DIJK23')
2171C
2172      ISYMK  = MULD2H(ISYMD,ISYMLH)
2173      ISYMIJ = 1
2174C
2175      IF (ISYMK.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK23')
2176C
2177C--------------------------------
2178C     Calculate the contribution.
2179C--------------------------------
2180C
2181      DO 100 K = 1,NRHF(ISYMK)
2182C
2183         NDEK = IGLMRH(ISYMD,ISYMK) + NBAS(ISYMD)*(K - 1)
2184     *        + IDEL - IBAS(ISYMD)
2185         NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1) + 1
2186C
2187         FACT = -TWO*XLAMDH(NDEK)
2188C
2189         CALL DAXPY(NMATIJ(ISYMIJ),FACT,XMAT,1,D2IJK(NIJK),1)
2190C
2191  100 CONTINUE
2192C
2193      CALL QEXIT('DIJK23')
2194C
2195      RETURN
2196      END
2197C  /* Deck dijk21 */
2198      SUBROUTINE DIJK21(D2IJK,ISYIJK,XMAT,XLAMDH,ISYMLH,
2199     *                  WORK,LWORK,IDEL,ISYMD)
2200C
2201C     Written by Asger Halkier 10/2 - 1996
2202C
2203C     Version: 1.0
2204C
2205C     Purpose: To calculate the first and second contribution to D2IJK
2206C              from projection against doubles!
2207C
2208#include "implicit.h"
2209      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
2210      DIMENSION D2IJK(*), XMAT(*), XLAMDH(*), WORK(LWORK)
2211#include "priunit.h"
2212#include "ccorb.h"
2213#include "ccsdsym.h"
2214C
2215      CALL QENTER('DIJK21')
2216C
2217      ISYML  = MULD2H(ISYMD,ISYMLH)
2218      ISYMX1 = ISYML
2219C
2220      IF (ISYML.NE.ISYIJK) CALL QUIT('Symmetry mismatch in DIJK21')
2221C
2222C-----------------------------------------------------------------
2223C     Calculate contraction of X-intermediate & lamda hole matrix.
2224C-----------------------------------------------------------------
2225C
2226      IF (LWORK .LT. NRHF(ISYMX1)) THEN
2227         WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NRHF(ISYMX1)
2228         CALL QUIT('Insufficient work space in DIJK21')
2229      ENDIF
2230C
2231      CALL DZERO(WORK,NRHF(ISYMX1))
2232C
2233      KOFF1 = IMATIJ(ISYMX1,ISYML) + 1
2234      KOFF2 = IGLMRH(ISYMD,ISYML) + IDEL - IBAS(ISYMD)
2235C
2236      NTOTI = MAX(NRHF(ISYMX1),1)
2237C
2238      CALL DGEMV('N',NRHF(ISYMX1),NRHF(ISYML),ONE,XMAT(KOFF1),NTOTI,
2239     *           XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1)
2240C
2241C----------------------------------
2242C     Calculate first contribution.
2243C----------------------------------
2244C
2245      ISYMI = ISYMX1
2246C
2247      DO 100 ISYMK = 1,NSYM
2248C
2249         ISYMJ  = ISYMK
2250         ISYMIJ = MULD2H(ISYMJ,ISYMI)
2251C
2252         DO 110 K = 1,NRHF(ISYMK)
2253C
2254            J    = K
2255            NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
2256     *           + IMATIJ(ISYMI,ISYMJ)  + NRHF(ISYMI)*(J - 1) + 1
2257C
2258            CALL DAXPY(NRHF(ISYMI),ONE,WORK,1,D2IJK(NIJK),1)
2259C
2260  110    CONTINUE
2261  100 CONTINUE
2262C
2263C-----------------------------------
2264C     Calculate second contribution.
2265C-----------------------------------
2266C
2267      ISYMK  = ISYMX1
2268      ISYMIJ = 1
2269C
2270      DO 120 K = 1,NRHF(ISYMK)
2271C
2272         DO 130 ISYMJ = 1,NSYM
2273C
2274            ISYMI = ISYMJ
2275C
2276            DO 140 J = 1,NRHF(ISYMJ)
2277C
2278               I    = J
2279               NIJK = IMAIJK(ISYMIJ,ISYMK) + NMATIJ(ISYMIJ)*(K - 1)
2280     *              + IMATIJ(ISYMI,ISYMJ)  + NRHF(ISYMI)*(J - 1) + I
2281C
2282               D2IJK(NIJK) = D2IJK(NIJK) - TWO*WORK(K)
2283C
2284  140       CONTINUE
2285  130    CONTINUE
2286  120 CONTINUE
2287C
2288      CALL QEXIT('DIJK21')
2289C
2290      RETURN
2291      END
2292C  /* Deck daij22 */
2293      SUBROUTINE DAIJ22(D2AIJ,ISYAIJ,DAB,ISYD1,XLAMDH,ISYMLH,
2294     *                  WORK,LWORK,IDEL,ISYMD)
2295C
2296C     Written by Asger Halkier 10/2 - 1996
2297C
2298C     Version: 1.0
2299C
2300C     Purpose: To calculate the second contribution to D2AIJ
2301C              from projection against doubles!
2302C
2303#include "implicit.h"
2304      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2305      DIMENSION D2AIJ(*), DAB(*), XLAMDH(*), WORK(LWORK)
2306#include "priunit.h"
2307#include "ccorb.h"
2308#include "ccsdsym.h"
2309C
2310      CALL QENTER('DAIJ22')
2311C
2312      ISYMB  = MULD2H(ISYMD,ISYMLH)
2313      ISYMA  = MULD2H(ISYMB,ISYD1)
2314      ISYMIJ = 1
2315C
2316      IF (ISYMA.NE.ISYAIJ) CALL QUIT('Symmetry mismatch in DAIJ22')
2317C
2318C-----------------------------------------------------------------
2319C     Calculate contraction of Y-intermediate & lamda hole matrix.
2320C-----------------------------------------------------------------
2321C
2322      IF (LWORK .LT. NVIR(ISYMA)) THEN
2323         WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', NVIR(ISYMA)
2324         CALL QUIT('Insufficient work space in DAIJ22')
2325      ENDIF
2326C
2327      CALL DZERO(WORK,NVIR(ISYMA))
2328C
2329      KOFF1 = IMATAB(ISYMA,ISYMB) + 1
2330      KOFF2 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD)
2331C
2332      NTOTA = MAX(NVIR(ISYMA),1)
2333C
2334      CALL DGEMV('N',NVIR(ISYMA),NVIR(ISYMB),ONE,DAB(KOFF1),NTOTA,
2335     *           XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK,1)
2336C
2337C-------------------------------------
2338C     Calculate contribution to D2AIJ.
2339C-------------------------------------
2340C
2341      DO 100 ISYMJ = 1,NSYM
2342C
2343         ISYMI  = ISYMJ
2344         ISYMAI = MULD2H(ISYMI,ISYMA)
2345C
2346         DO 110 J = 1,NRHF(ISYMJ)
2347C
2348            I    = J
2349            NAIJ = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1)
2350     *           + IT1AM(ISYMA,ISYMI)   + NVIR(ISYMA)*(I - 1) + 1
2351C
2352            CALL DAXPY(NVIR(ISYMA),-ONE,WORK,1,D2AIJ(NAIJ),1)
2353C
2354  110    CONTINUE
2355  100 CONTINUE
2356C
2357      CALL QEXIT('DAIJ22')
2358C
2359      RETURN
2360      END
2361C  /* Deck dabi21 */
2362      SUBROUTINE DABI21(D2ABI,ISYABI,DAB,ISYD1,XLAMDH,ISYMLH,IDEL,ISYMD)
2363C
2364C     Written by Asger Halkier 11/2 - 1996
2365C
2366C     Version: 1.0
2367C
2368C     Purpose: To calculate the first contribution to D2ABI
2369C              from projection against doubles!
2370C
2371#include "implicit.h"
2372      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
2373      DIMENSION D2ABI(*), DAB(*), XLAMDH(*)
2374#include "ccorb.h"
2375#include "ccsdsym.h"
2376C
2377      CALL QENTER('DABI21')
2378C
2379      ISYMI  = MULD2H(ISYMD,ISYMLH)
2380      ISYMAB = ISYD1
2381C
2382      IF (MULD2H(ISYMI,ISYMAB).NE.ISYABI)
2383     *  CALL QUIT('Symmetry mismatch in DABI21')
2384C
2385C--------------------------------
2386C     Calculate the contribution.
2387C--------------------------------
2388C
2389      DO 100 I = 1,NRHF(ISYMI)
2390C
2391         NDEI = IGLMRH(ISYMD,ISYMI) + NBAS(ISYMD)*(I - 1)
2392     *        + IDEL - IBAS(ISYMD)
2393         NABI = ICKASR(ISYMAB,ISYMI) + NMATAB(ISYMAB)*(I - 1) + 1
2394C
2395         FACT = TWO*XLAMDH(NDEI)
2396C
2397         CALL DAXPY(NMATAB(ISYMAB),FACT,DAB,1,D2ABI(NABI),1)
2398C
2399  100 CONTINUE
2400C
2401      CALL QEXIT('DABI21')
2402C
2403      RETURN
2404      END
2405C  /* Deck diaj21 */
2406      SUBROUTINE DIAJ21(D2IAJ,ISYAIJ,YMAT,T2TINT,ISYMTI)
2407C
2408C     Written by Asger Halkier 11/2 - 1996
2409C
2410C     Version: 1.0
2411C
2412C     Purpose: To calculate the first contribution to D2IAJ
2413C              from projection against doubles!
2414C
2415#include "implicit.h"
2416      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2417      DIMENSION D2IAJ(*), YMAT(*), T2TINT(*)
2418#include "ccorb.h"
2419#include "ccsdsym.h"
2420C
2421      CALL QENTER('DIAJ21')
2422C
2423      ISYMAE = 1
2424C
2425      IF (ISYMTI.NE.ISYAIJ) CALL QUIT('Symmetry mismatch in DIAJ21.')
2426C
2427C--------------------------------
2428C     Calculate the contribution.
2429C--------------------------------
2430C
2431      DO 100 ISYMJ = 1,NSYM
2432C
2433         ISYMEI = MULD2H(ISYMJ,ISYMTI)
2434         ISYMAI = MULD2H(ISYMJ,ISYAIJ)
2435C
2436         DO 110 J = 1,NRHF(ISYMJ)
2437C
2438            DO 120 ISYMA = 1,NSYM
2439C
2440               ISYME = MULD2H(ISYMA,ISYMAE)
2441               ISYMI = MULD2H(ISYMA,ISYMAI)
2442C
2443               KOFF1 = IMATAB(ISYMA,ISYME)  + 1
2444               KOFF2 = IT2BCD(ISYMEI,ISYMJ) + NT1AM(ISYMEI)*(J - 1)
2445     *               + IT1AM(ISYME,ISYMI)   + 1
2446               KOFF3 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1)
2447     *               + IT1AM(ISYMA,ISYMI)   + 1
2448C
2449               NTOTA = MAX(NVIR(ISYMA),1)
2450               NTOTE = MAX(NVIR(ISYME),1)
2451C
2452               CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYME),
2453     *                    -ONE,YMAT(KOFF1),NTOTA,T2TINT(KOFF2),NTOTE,
2454     *                    ONE,D2IAJ(KOFF3),NTOTA)
2455C
2456  120       CONTINUE
2457  110    CONTINUE
2458  100 CONTINUE
2459C
2460      CALL QEXIT('DIAJ21')
2461C
2462      RETURN
2463      END
2464C  /* Deck diaj22 */
2465      SUBROUTINE DIAJ22(D2IAJ,ISYAIJ,XMAT,T2TINT,ISYMTI)
2466C
2467C     Written by Asger Halkier 11/2 - 1996
2468C
2469C     Version: 1.0
2470C
2471C     Purpose: To calculate the second contribution to D2IAJ
2472C              from projection against doubles!
2473C
2474#include "implicit.h"
2475      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2476      DIMENSION D2IAJ(*), XMAT(*), T2TINT(*)
2477#include "ccorb.h"
2478#include "ccsdsym.h"
2479C
2480      CALL QENTER('DIAJ22')
2481C
2482      IF (ISYAIJ.NE.ISYMTI) CALL QUIT('Symmetry mismatch in DIAJ22.')
2483C
2484C--------------------------------
2485C     Calculate the contribution.
2486C--------------------------------
2487C
2488      DO 100 ISYMJ = 1,NSYM
2489C
2490         ISYMAM = MULD2H(ISYMJ,ISYMTI)
2491         ISYMAI = MULD2H(ISYMJ,ISYAIJ)
2492C
2493         DO 110 J = 1,NRHF(ISYMJ)
2494C
2495            DO 120 ISYMA = 1,NSYM
2496C
2497               ISYMM = MULD2H(ISYMA,ISYMAM)
2498               ISYMI = MULD2H(ISYMA,ISYMAI)
2499C
2500               KOFF1 = IT2BCD(ISYMAM,ISYMJ) + NT1AM(ISYMAM)*(J - 1)
2501     *               + IT1AM(ISYMA,ISYMM)   + 1
2502               KOFF2 = IMATIJ(ISYMI,ISYMM)  + 1
2503               KOFF3 = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1)
2504     *               + IT1AM(ISYMA,ISYMI)   + 1
2505C
2506               NTOTA = MAX(NVIR(ISYMA),1)
2507               NTOTI = MAX(NRHF(ISYMI),1)
2508C
2509               CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMM),
2510     *                    -ONE,T2TINT(KOFF1),NTOTA,XMAT(KOFF2),NTOTI,
2511     *                    ONE,D2IAJ(KOFF3),NTOTA)
2512C
2513  120       CONTINUE
2514  110    CONTINUE
2515  100 CONTINUE
2516C
2517      CALL QEXIT('DIAJ22')
2518C
2519      RETURN
2520      END
2521C  /* Deck diaj23 */
2522      SUBROUTINE DIAJ23(D2IAJ,ISYAIJ,XMAT,T2TINT,ISYMTI)
2523C
2524C     Written by Asger Halkier 11/2 - 1996
2525C
2526C     Version: 1.0
2527C
2528C     Purpose: To calculate the third contribution to D2IAJ
2529C              from projection against doubles!
2530C
2531#include "implicit.h"
2532      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2533      DIMENSION D2IAJ(*), XMAT(*), T2TINT(*)
2534#include "ccorb.h"
2535#include "ccsdsym.h"
2536C
2537      CALL QENTER('DIAJ23')
2538C
2539      ISYMJM = 1
2540C
2541      IF (ISYMTI.NE.ISYAIJ) CALL QUIT('Symmetry mismatch in DIAJ23')
2542C
2543C--------------------------------
2544C     Calculate the contribution.
2545C--------------------------------
2546C
2547      DO 100 ISYMAI = 1,NSYM
2548C
2549         ISYMM = MULD2H(ISYMAI,ISYMTI)
2550         ISYMJ = MULD2H(ISYMAI,ISYAIJ)
2551C
2552         KOFF1 = IT2BCD(ISYMAI,ISYMM) + 1
2553         KOFF2 = IMATIJ(ISYMJ,ISYMM)  + 1
2554         KOFF3 = IT2BCD(ISYMAI,ISYMJ) + 1
2555C
2556         NTOAI = MAX(NT1AM(ISYMAI),1)
2557         NTOTJ = MAX(NRHF(ISYMJ),1)
2558C
2559         CALL DGEMM('N','T',NT1AM(ISYMAI),NRHF(ISYMJ),NRHF(ISYMM),-ONE,
2560     *              T2TINT(KOFF1),NTOAI,XMAT(KOFF2),NTOTJ,ONE,
2561     *              D2IAJ(KOFF3),NTOAI)
2562C
2563  100 CONTINUE
2564C
2565      CALL QEXIT('DIAJ23')
2566C
2567      RETURN
2568      END
2569C  /* Deck diaj24 */
2570      SUBROUTINE DIAJ24(D2IAJ,ISYIAJ,T2AMT,DAB,ISYD1,XLAMDH,ISYMLH,
2571     *                  WORK,LWORK,IDEL,ISYMD)
2572C
2573C     Written by Asger Halkier 13/2 - 1996
2574C
2575C     Version: 1.0
2576C
2577C     Purpose: To calculate the fourth contribution to D2IAJ
2578C              from projection against doubles!
2579C
2580#include "implicit.h"
2581      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2582      DIMENSION D2IAJ(*), T2AMT(*), DAB(*), XLAMDH(*), WORK(LWORK)
2583#include "priunit.h"
2584#include "ccorb.h"
2585#include "ccsdsym.h"
2586C
2587      INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J
2588C
2589      CALL QENTER('DIAJ24')
2590C
2591      ISYMB  = MULD2H(ISYMD,ISYMLH)
2592      ISYMC  = MULD2H(ISYMB,ISYD1)
2593C
2594      IF (MULD2H(ISYMLH,ISYMD).NE.MULD2H(ISYD1,ISYIAJ)) THEN
2595        CALL QUIT('Symmetry mismatch in DIAJ24.')
2596      END IF
2597C
2598      IF (NVIR(ISYMB) .EQ. 0) THEN
2599        CALL QEXIT('DIAJ24')
2600        RETURN
2601      END IF
2602C
2603C-------------------------------
2604C     Work space allocation one.
2605C-------------------------------
2606C
2607      KVECC = 1
2608      KEND1 = KVECC + NVIR(ISYMC)
2609      LWRK1 = LWORK - KEND1
2610C
2611      IF (LWRK1 .LT. 0) THEN
2612         WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1
2613         CALL QUIT('Insufficient work space for allocation in DIAJ24')
2614      ENDIF
2615C
2616      CALL DZERO(WORK(KVECC),NVIR(ISYMC))
2617C
2618C-------------------------------------------------------------------
2619C     Calculate contraction of transposed Y-matrix (DAB) and XLAMDH.
2620C-------------------------------------------------------------------
2621C
2622      KOFF1 = IMATAB(ISYMC,ISYMB) + 1
2623      KOFF2 = IGLMVI(ISYMD,ISYMB) + IDEL - IBAS(ISYMD)
2624C
2625      NTOTC = MAX(NVIR(ISYMC),1)
2626C
2627      CALL DGEMV('N',NVIR(ISYMC),NVIR(ISYMB),ONE,DAB(KOFF1),NTOTC,
2628     *           XLAMDH(KOFF2),NBAS(ISYMD),ZERO,WORK(KVECC),1)
2629C
2630C--------------------------------
2631C     Calculate the contribution.
2632C--------------------------------
2633C
2634      DO 100 ISYMJ = 1,NSYM
2635C
2636         ISYMAI = MULD2H(ISYMJ,ISYIAJ)
2637         ISYMCJ = ISYMAI
2638C
2639         DO 110 J = 1,NRHF(ISYMJ)
2640C
2641            DO 120 C = 1,NVIR(ISYMC)
2642C
2643               NCJ = IT1AM(ISYMC,ISYMJ) + NVIR(ISYMC)*(J - 1) + C
2644C
2645               DO 130 NAI = 1,NT1AM(ISYMAI)
2646C
2647                  NAICJ = IT2AM(ISYMAI,ISYMCJ) + INDEX(NAI,NCJ)
2648                  NAIJ  = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J - 1)
2649     *                  + NAI
2650C
2651                  D2IAJ(NAIJ) = D2IAJ(NAIJ)
2652     *                        - T2AMT(NAICJ)*WORK(KVECC + C - 1)
2653C
2654  130          CONTINUE
2655  120       CONTINUE
2656  110    CONTINUE
2657  100 CONTINUE
2658C
2659      CALL QEXIT('DIAJ24')
2660C
2661      RETURN
2662      END
2663C  /* Deck ccd1reor */
2664      SUBROUTINE CCD1REOR(DFOCK,DAI,DAB,DIJ,DIA)
2665C
2666C     Written by Asger Halkier 28/2 - 1996
2667C
2668C     Version: 1.0
2669C
2670C     Purpose: To reorder the four separate blocks DAI through DIA
2671C              of the CC one electron density to one matrix stored
2672C              as a directly diagonalisable Fock-matrix (DFOCK)!
2673C
2674#include "implicit.h"
2675      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2676      DIMENSION DFOCK(*), DAI(*), DAB(*), DIJ(*), DIA(*)
2677#include "ccorb.h"
2678#include "ccsdsym.h"
2679C
2680      CALL QENTER('CCD1REOR')
2681C
2682      DO 100 ISYM = 1,NSYM
2683C
2684C---------------------------------------
2685C        Do the occupied occupied block.
2686C---------------------------------------
2687C
2688         DO 110 J = 1,NRHF(ISYM)
2689C
2690            NIJOR = IMATIJ(ISYM,ISYM) + NRHF(ISYM)*(J - 1) + 1
2691            NIJNE = IFCKDO(ISYM) + NORB(ISYM)*(J - 1) + 1
2692C
2693            CALL DCOPY(NRHF(ISYM),DIJ(NIJOR),1,DFOCK(NIJNE),1)
2694C
2695  110    CONTINUE
2696C
2697C--------------------------------------
2698C        Do the virtual occupied block.
2699C--------------------------------------
2700C
2701         DO 120 I = 1,NRHF(ISYM)
2702C
2703            NAIOR = IT1AM(ISYM,ISYM)  + NVIR(ISYM)*(I - 1) + 1
2704            NAINE = IFCKDO(ISYM) + NORB(ISYM)*(I - 1)
2705     *            + NRHF(ISYM) + 1
2706C
2707            CALL DCOPY(NVIR(ISYM),DAI(NAIOR),1,DFOCK(NAINE),1)
2708C
2709  120    CONTINUE
2710C
2711C--------------------------------------
2712C        Do the occupied virtual block.
2713C--------------------------------------
2714C
2715         DO 130 A = 1,NVIR(ISYM)
2716C
2717            NIAOR = IT1AM(ISYM,ISYM)  + A
2718            NIANE = IFCKDV(ISYM) + NORB(ISYM)*(A - 1) + 1
2719C
2720            CALL DCOPY(NRHF(ISYM),DIA(NIAOR),NVIR(ISYM),
2721     *                 DFOCK(NIANE),1)
2722C
2723  130    CONTINUE
2724C
2725C-------------------------------------
2726C        Do the virtual virtual block.
2727C-------------------------------------
2728C
2729         DO 140 B = 1,NVIR(ISYM)
2730C
2731            NABOR = IMATAB(ISYM,ISYM) + NVIR(ISYM)*(B - 1) + 1
2732            NABNE = IFCKDV(ISYM) + NORB(ISYM)*(B - 1)
2733     *            + NRHF(ISYM) + 1
2734C
2735            CALL DCOPY(NVIR(ISYM),DAB(NABOR),1,DFOCK(NABNE),1)
2736C
2737  140    CONTINUE
2738  100 CONTINUE
2739C
2740      CALL QEXIT('CCD1REOR')
2741C
2742      RETURN
2743      END
2744C  /* Deck ccd1t1co */
2745      SUBROUTINE CCD1T1CO(DAI,DAB,DIJ,DIA,T1AM,YMAT,WORK,LWORK)
2746C
2747C     Written by Asger Halkier 29/2 - 1996
2748C
2749C     Version: 1.0
2750C
2751C     Purpose: To add contributions from the t1-amplitudes to
2752C              the CC one electron density in MO-basis.
2753C                                             --
2754C
2755#include "implicit.h"
2756      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2757      DIMENSION DAI(*), DAB(*), DIJ(*), DIA(*)
2758      DIMENSION T1AM(*), YMAT(*), WORK(LWORK)
2759#include "priunit.h"
2760#include "ccorb.h"
2761#include "ccsdsym.h"
2762C
2763      CALL QENTER('CCD1T1CO')
2764C
2765C----------------------------------------------
2766C     Add contributions to the (occ,vir) block.
2767C----------------------------------------------
2768C
2769      DO 100 ISYMA = 1,NSYM
2770C
2771         ISYMI = ISYMA
2772         ISYMK = ISYMA
2773         ISYMC = ISYMA
2774C
2775         KOFF1 = IT1AM(ISYMA,ISYMK)  + 1
2776         KOFF2 = IMATIJ(ISYMI,ISYMK) + 1
2777         KOFF3 = IT1AM(ISYMA,ISYMI)  + 1
2778C
2779         NTOTA = MAX(NVIR(ISYMA),1)
2780         NTOTI = MAX(NRHF(ISYMI),1)
2781C
2782         CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK),
2783     *              ONE,T1AM(KOFF1),NTOTA,DIJ(KOFF2),NTOTI,ONE,
2784     *              DIA(KOFF3),NTOTA)
2785C
2786         KOFF4 = IMATAB(ISYMA,ISYMC) + 1
2787         KOFF5 = IT1AM(ISYMC,ISYMI)  + 1
2788         KOFF6 = IT1AM(ISYMA,ISYMI)  + 1
2789C
2790         NTOTA = MAX(NVIR(ISYMA),1)
2791         NTOTC = MAX(NVIR(ISYMC),1)
2792C
2793         CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMC),
2794     *              -ONE,YMAT(KOFF4),NTOTA,T1AM(KOFF5),NTOTC,ONE,
2795     *              DIA(KOFF6),NTOTA)
2796C
2797C------------------------------
2798C        Work space allocation.
2799C------------------------------
2800C
2801         KSCR  = 1
2802         KEND1 = KSCR  + NRHF(ISYMK)*NRHF(ISYMI)
2803         LWRK1 = LWORK - KEND1
2804C
2805         IF (LWRK1 .LT. 0) THEN
2806            WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1
2807            CALL QUIT('Insufficient work space for '//
2808     &           'allocation in CCD1T1CO')
2809         ENDIF
2810C
2811         CALL DZERO(WORK(KSCR),NRHF(ISYMK)*NRHF(ISYMI))
2812C
2813         KOFF7 = IT1AM(ISYMC,ISYMK) + 1
2814         KOFF8 = IT1AM(ISYMC,ISYMI) + 1
2815C
2816         NTOTC = MAX(NVIR(ISYMC),1)
2817         NTOTK = MAX(NRHF(ISYMK),1)
2818C
2819         CALL DGEMM('T','N',NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC),
2820     *              ONE,DAI(KOFF7),NTOTC,T1AM(KOFF8),NTOTC,ZERO,
2821     *              WORK(KSCR),NTOTK)
2822C
2823         KOFF9  = IT1AM(ISYMA,ISYMK) + 1
2824         KOFF10 = IT1AM(ISYMA,ISYMI) + 1
2825C
2826         NTOTA  = MAX(NVIR(ISYMA),1)
2827         NTOTK  = MAX(NRHF(ISYMK),1)
2828C
2829         CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK),
2830     *              -ONE,T1AM(KOFF9),NTOTA,WORK(KSCR),NTOTK,ONE,
2831     *              DIA(KOFF10),NTOTA)
2832C
2833  100 CONTINUE
2834C
2835C----------------------------------------------
2836C     Add contributions to the (vir,vir) block.
2837C----------------------------------------------
2838C
2839      DO 110 ISYMA = 1,NSYM
2840C
2841         ISYMB  = ISYMA
2842         ISYMK  = ISYMA
2843C
2844         KOFF11 = IT1AM(ISYMA,ISYMK)  + 1
2845         KOFF12 = IT1AM(ISYMB,ISYMK)  + 1
2846         KOFF13 = IMATAB(ISYMA,ISYMB) + 1
2847C
2848         NTOTA  = MAX(NVIR(ISYMA),1)
2849         NTOTB  = MAX(NVIR(ISYMB),1)
2850C
2851         CALL DGEMM('N','T',NVIR(ISYMA),NVIR(ISYMB),NRHF(ISYMK),
2852     *              ONE,DAI(KOFF11),NTOTA,T1AM(KOFF12),NTOTB,ONE,
2853     *              DAB(KOFF13),NTOTA)
2854C
2855  110 CONTINUE
2856C
2857C----------------------------------------------
2858C     Add contributions to the (occ,occ) block.
2859C----------------------------------------------
2860C
2861      DO 120 ISYMI = 1,NSYM
2862C
2863         ISYMJ  = ISYMI
2864         ISYMC  = ISYMI
2865C
2866         KOFF14 = IT1AM(ISYMC,ISYMI)  + 1
2867         KOFF15 = IT1AM(ISYMC,ISYMJ)  + 1
2868         KOFF16 = IMATIJ(ISYMI,ISYMJ) + 1
2869C
2870         NTOTC  = MAX(NVIR(ISYMC),1)
2871         NTOTI  = MAX(NRHF(ISYMI),1)
2872C
2873         CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMJ),NVIR(ISYMC),
2874     *              -ONE,T1AM(KOFF14),NTOTC,DAI(KOFF15),NTOTC,ONE,
2875     *              DIJ(KOFF16),NTOTI)
2876C
2877  120 CONTINUE
2878C
2879      CALL QEXIT('CCD1T1CO')
2880C
2881      RETURN
2882      END
2883C  /* Deck sortash */
2884      SUBROUTINE SORTASH(XARR,YARR,N)
2885C
2886C     Written by Asger Halkier 4/2 - 1996
2887C
2888C     Version: 1.0
2889C
2890C     Purpose: To reorder the array xarr after size, and
2891C              do the same reordering of the concomitant
2892C              array yarr!
2893C              (This routine is based on Bubble-sort from NUM.REC.)
2894C
2895#include "implicit.h"
2896      DIMENSION XARR(*), YARR(*)
2897#include "ccorb.h"
2898#include "ccsdsym.h"
2899C
2900      CALL QENTER('SORTASH')
2901C
2902C----------------------
2903C     Do the resorting.
2904C----------------------
2905C
2906      DO 100 J = 2,N
2907C
2908         XA = XARR(J)
2909         XB = YARR(J)
2910C
2911         DO 110 I = J-1,1,-1
2912C
2913            IF (XARR(I) .LE. XA) GOTO 120
2914C
2915            XARR(I + 1) = XARR(I)
2916            YARR(I + 1) = YARR(I)
2917C
2918  110    CONTINUE
2919C
2920         I=0
2921C
2922  120    XARR(I + 1) = XA
2923         YARR(I + 1) = XB
2924C
2925  100 CONTINUE
2926C
2927      CALL QEXIT('SORTASH')
2928C
2929      RETURN
2930      END
2931C  /* Deck ccnaocan */
2932      SUBROUTINE CCNAOCAN(XARR,YARR)
2933C
2934C     Written by Asger Halkier 4/2 - 1996
2935C
2936C     Version: 1.0
2937C
2938C     Purpose: To analyze and check the natural occupation
2939C              numbers in xarr and yarr (real and imag)!
2940C
2941#include "implicit.h"
2942      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
2943      DIMENSION XARR(*), YARR(*)
2944#include "priunit.h"
2945#include "ccorb.h"
2946#include "ccsdsym.h"
2947#include "ccsdinp.h"
2948C
2949      CALL QENTER('CCNAOCAN')
2950C
2951      SUMRHF = ZERO
2952      SUMALL = ZERO
2953      SUMIMA = ZERO
2954C
2955      DO 100 I = 1,NORBT
2956C
2957         SUMALL = SUMALL + XARR(I)
2958         SUMIMA = SUMIMA + YARR(I)
2959C
2960  100 CONTINUE
2961C
2962      CHECK = 2*NRHFT - SUMALL
2963      CHEPR = NRHFT*2.0D0
2964      THR   = 1.0D-06
2965C
2966      IF (CHECK .LT. THR) THEN
2967C
2968         WRITE(LUPRI,*) ' '
2969         WRITE(LUPRI,*) ' '
2970         WRITE(LUPRI,111) 'Total Sum of natural occupation numbers:',
2971     &        CHEPR
2972C
2973      ELSE
2974C
2975         WRITE(LUPRI,*) ' '
2976         WRITE(LUPRI,*) ' '
2977         WRITE(LUPRI,111) 'Total Sum of natural occupation numbers:',
2978     &        SUMALL
2979C
2980      ENDIF
2981C
2982      IF (IPRINT .GT. 20) THEN
2983C
2984         WRITE(LUPRI,*) ' '
2985         WRITE(LUPRI,111) 'Total Sum of imaginary parts of nat.occ:',
2986     &        SUMIMA
2987C
2988      ENDIF
2989C
2990  111 FORMAT(3X,A40,2X,F9.6)
2991C
2992      DO 110 I = 1,NRHFT
2993C
2994         SUMRHF = SUMRHF + XARR(NORBT - I + 1)
2995C
2996  110 CONTINUE
2997C
2998      TOMOVE = 2*NRHFT - SUMRHF
2999C
3000      WRITE(LUPRI,*) ' '
3001      WRITE (LUPRI,222) 'Dynamical correlation move:',
3002     *                   TOMOVE,'electrons'
3003C
3004  222 FORMAT(3X,A27,F9.6,1X,A9)
3005C
3006      CALL QEXIT('CCNAOCAN')
3007C
3008      RETURN
3009      END
3010C  /* Deck cc_fcd1ao */
3011      SUBROUTINE CC_FCD1AO(AODEN,WORK,LWORK,MODEL)
3012C
3013C     Written by Ove Christiansen 2-may 1996.
3014C
3015C     Version: 1.0
3016C
3017C     Purpose: To calculate the Frozen core contribution to the one electron Coupled Cluster
3018C              density matrix and add it to AODEN.
3019C
3020#include "implicit.h"
3021#include "priunit.h"
3022#include "dummy.h"
3023#include "maxash.h"
3024#include "maxorb.h"
3025#include "mxcent.h"
3026#include "aovec.h"
3027#include "iratdef.h"
3028      PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0)
3029      DIMENSION AODEN(*), WORK(LWORK)
3030#include "ccorb.h"
3031#include "infind.h"
3032#include "inftap.h"
3033#include "blocks.h"
3034#include "ccsdinp.h"
3035#include "ccsdsym.h"
3036#include "ccsdio.h"
3037#include "distcl.h"
3038#include "cbieri.h"
3039#include "eribuf.h"
3040C
3041      CHARACTER MODEL*10
3042C
3043      CALL QENTER('CC_FCD1AO')
3044C
3045      IF (IPRINT .GT. 10) THEN
3046         CALL AROUND( 'Calculate core contribution to density')
3047      ENDIF
3048C
3049      KCMO  = 1
3050      KEND  = KCMO  + NLAMDS
3051      LWRK1 = LWORK - KEND
3052C
3053      IF (LWRK1 .LT. 0) THEN
3054         CALL QUIT('Insufficient spaces in CC_FCD1AO')
3055      ENDIF
3056C
3057C----------------------------------------------
3058C     Read MO-coefficients from interface file.
3059C----------------------------------------------
3060C
3061      LUSIFC = -1
3062      CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,
3063     &            .FALSE.)
3064      REWIND LUSIFC
3065C
3066      CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
3067      READ (LUSIFC)
3068C
3069      READ (LUSIFC)
3070      READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS)
3071C
3072      CALL GPCLOSE(LUSIFC,'KEEP')
3073C
3074      ICMOI  = KCMO - 1
3075      DO ISYMI = 1, NSYM
3076        ISYMA = ISYMI
3077        ISYMB = ISYMI
3078        DO IFR = 1, NRHFFR(ISYMI)
3079          I = KFRRHF(IFR,ISYMI)
3080          DO IAL = 1, NBAS(ISYMA)
3081            ICMOIA = ICMOI + NBAS(ISYMA)*(I-1) + IAL
3082            DO IBE = 1, NBAS(ISYMA)
3083              IDAB = IAODIS(ISYMA,ISYMB)+NBAS(ISYMA)*(IBE-1)+IAL
3084              ICMOIB = ICMOI + NBAS(ISYMB)*(I-1) + IBE
3085              AODEN(IDAB) = AODEN(IDAB) + TWO*WORK(ICMOIA)*WORK(ICMOIB)
3086            ENDDO
3087          ENDDO
3088        ENDDO
3089        ICMOI  = ICMOI  + NBAS(ISYMI)*NRHFS(ISYMI)
3090        ICMOI  = ICMOI  + NBAS(ISYMI)*NVIRS(ISYMI)
3091      ENDDO
3092C
3093      CALL QEXIT('CC_FCD1AO')
3094C
3095      END
3096C  /* Deck cc2_d1ao */
3097      SUBROUTINE CC2_D1AO(AODEN,ZKAM,WORK,LWORK)
3098C
3099C     Written by A. Halkier & S. Coriani 13/01-2000.
3100C
3101C     Version: 1.0
3102C
3103C     Purpose: To calculate the relaxed one electron CC2
3104C              density matrix and return it backtransformed
3105C              to AO-basis in AODEN! ZKAM is kappa-bar-0
3106C              (input), which contributes to the relaxed density.
3107C
3108C     Important note: NO FROZEN CORE SO FAR!!!
3109C
3110#include "implicit.h"
3111#include "priunit.h"
3112#include "dummy.h"
3113#include "ccinftap.h"
3114#include "maxash.h"
3115#include "maxorb.h"
3116#include "mxcent.h"
3117#include "aovec.h"
3118#include "iratdef.h"
3119      PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0)
3120      DIMENSION INDEXA(MXCORB_CC)
3121      DIMENSION AODEN(*), ZKAM(*), WORK(LWORK)
3122      CHARACTER*10 MODEL
3123#include "ccorb.h"
3124#include "infind.h"
3125#include "blocks.h"
3126#include "ccsdinp.h"
3127#include "ccsdsym.h"
3128#include "ccsdio.h"
3129#include "distcl.h"
3130#include "cbieri.h"
3131#include "eribuf.h"
3132#include "ccfop.h"
3133#include "ccfro.h"
3134C
3135      LOGICAL LOCDBG
3136      PARAMETER ( LOCDBG=.FALSE. )
3137C
3138      CALL QENTER('CC2_D1AO')
3139C
3140C------------------------------------------------------------------
3141C     Both t-vectors and tbar-vectors (zeta) are totally symmetric.
3142C------------------------------------------------------------------
3143C
3144      ISYMTR = 1
3145      ISYMOP = 1
3146C
3147C-------------------------------
3148C     Work space allocation one.
3149C-------------------------------
3150C
3151      KT2AM  = 1
3152      KZ2AM  = KT2AM  + NT2AMX
3153      KLAMDP = KZ2AM  + NT2SQ(1)
3154      KLAMDH = KLAMDP + NLAMDT
3155      KT1AM  = KLAMDH + NLAMDT
3156      KZ1AM  = KT1AM  + NT1AMX
3157      KEND1  = KZ1AM  + NT1AMX
3158      LWRK1  = LWORK  - KEND1
3159C
3160      IF (LWRK1 .LT. 0) THEN
3161         WRITE(LUPRI,*) 'Available:',
3162     *                   LWORK, 'Needed:', KEND1
3163         CALL QUIT(
3164     *    'Insufficient memory for initial allocation in CC2_D1AO')
3165      ENDIF
3166C
3167C----------------------------------------
3168C     Read zero'th order zeta amplitudes.
3169C----------------------------------------
3170C
3171      IOPT   = 3
3172      CALL CC_RDRSP('L0',0,1,IOPT,MODEL,WORK(KZ1AM),WORK(KZ2AM))
3173C
3174C--------------------------------
3175C     Square up zeta2 amplitudes.
3176C--------------------------------
3177C
3178      CALL DCOPY(NT2AMX,WORK(KZ2AM),1,WORK(KT2AM),1)
3179      CALL CC_T2SQ(WORK(KT2AM),WORK(KZ2AM),1)
3180C
3181C
3182C-------------------------------------------
3183C     Read zero'th order cluster amplitudes.
3184C-------------------------------------------
3185C
3186      IOPT = 3
3187      CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KT1AM),WORK(KT2AM))
3188C
3189C----------------------------------
3190C     Calculate the lamda matrices.
3191C----------------------------------
3192C
3193      CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1),
3194     *            LWRK1)
3195C
3196C-----------------------------------------------
3197C     Set up 2C-E of cluster amplitudes and save
3198C     in KT2AM, as we only need T(2c-e) below.
3199C-----------------------------------------------
3200C
3201      ISYOPE = 1
3202      IOPTTCME = 1
3203      CALL CCSD_TCMEPK(WORK(KT2AM),1.0D0,ISYOPE,IOPTTCME)
3204      KT2AMT = KT2AM                  !for safety
3205C
3206c     write(6,*) 'after initial stuff'
3207c     CALL FLSHFO(LUPRI)
3208C
3209C-------------------------------
3210C     Work space allocation one.
3211C     Note that D(ai) = ZETA(ai)
3212C     and both D(ia) and h(ia)
3213C     are stored transposed!
3214C-------------------------------
3215C
3216      KONEAI = KZ1AM
3217      KONEAB = KONEAI + NT1AMX
3218      KONEIJ = KONEAB + NMATAB(1)
3219      KONEIA = KONEIJ + NMATIJ(1)
3220      KCMO   = KONEIA + NT1AMX
3221      KDUM   = KCMO   + NLAMDS
3222      KEND1  = KDUM   + NMATIJ(1)
3223      LWRK1  = LWORK  - KEND1
3224C
3225      IF (LWRK1 .LT. 0) THEN
3226         WRITE(LUPRI,*) 'Available:', LWORK, 'Needed:', KEND1
3227         CALL QUIT('Insufficient memory for allocation 1 CC2_D1AO')
3228      ENDIF
3229C
3230C
3231C------------------------------------------------------
3232C     Initialize remaining one electron density arrays.
3233C------------------------------------------------------
3234C
3235      CALL DZERO(WORK(KONEAB),NMATAB(1))
3236      CALL DZERO(WORK(KONEIJ),NMATIJ(1))
3237      CALL DZERO(WORK(KONEIA),NT1AMX)
3238C
3239C--------------------------------------------------------
3240C     Construct remaining blocks of one electron density.
3241C--------------------------------------------------------
3242C
3243      CALL DZERO(WORK(KDUM),NMATIJ(1))
3244      CALL DIJGEN(WORK(KONEIJ),WORK(KDUM))
3245      CALL DIAGEN(WORK(KONEIA),WORK(KT2AMT),WORK(KONEAI))
3246C
3247c     write(6,*) 'after D-lambda'
3248c     CALL FLSHFO(LUPRI)
3249C
3250      IF (LOCDBG) THEN
3251         DIJNO = DDOT(NMATIJ(1),WORK(KONEIJ),1,WORK(KONEIJ),1)
3252         DAINO = DDOT(NT1AMX,WORK(KONEAI),1,WORK(KONEAI),1)
3253         DIANO = DDOT(NT1AMX,WORK(KONEIA),1,WORK(KONEIA),1)
3254         DABNO = DDOT(NMATAB(1),WORK(KONEAB),1,WORK(KONEAB),1)
3255         WRITE(LUPRI,*) ' '
3256         WRITE(LUPRI,*) 'Norm of one electron densities in MO-basis:'
3257         WRITE(LUPRI,*) DIJNO, DAINO, DIANO, DABNO
3258      ENDIF
3259C
3260C--------------------------------------------------------
3261C     Backtransform the one electron density to AO-basis.
3262C--------------------------------------------------------
3263C
3264      CALL DZERO(AODEN,N2BST(1))
3265C
3266      ISDEN = 1
3267      CALL CC_DENAO(AODEN,1,WORK(KONEAI),WORK(KONEAB),
3268     *              WORK(KONEIJ),WORK(KONEIA),ISDEN,WORK(KLAMDP),1,
3269     *              WORK(KLAMDH),1,WORK(KEND1),LWRK1)
3270C
3271C----------------------------------------------------------
3272C     Read MO-coefficients from interface file and reorder.
3273C----------------------------------------------------------
3274C
3275      LUSIFC = -1
3276      CALL GPOPEN(LUSIFC,'SIRIFC','UNKNOWN',' ','UNFORMATTED',
3277     *            IDUMMY,.FALSE.)
3278      REWIND LUSIFC
3279      CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
3280      READ (LUSIFC)
3281      READ (LUSIFC)
3282      READ (LUSIFC) (WORK(KCMO+I-1), I=1,NLAMDS)
3283      CALL GPCLOSE (LUSIFC,'KEEP')
3284C
3285      CALL CMO_REORDER(WORK(KCMO),WORK(KEND1),LWRK1)
3286C
3287C----------------------------------
3288C     Add orbital relaxation terms.
3289C----------------------------------
3290C
3291      KOFDIJ = 1
3292      KOFDAB = KOFDIJ + NMATIJ(1)
3293      KOFDAI = KOFDAB + NMATAB(1)
3294      KOFDIA = KOFDAI + NT1AMX
3295C
3296      ISDEN = 1
3297      CALL CC_DENAO(AODEN,1,ZKAM(KOFDAI),ZKAM(KOFDAB),
3298     *              ZKAM(KOFDIJ),ZKAM(KOFDIA),ISDEN,WORK(KCMO),1,
3299     *              WORK(KCMO),1,WORK(KEND1),LWRK1)
3300C
3301      CALL QEXIT('CC2_D1AO')
3302      RETURN
3303      END
3304