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
19*----------------------------------------------------------------*
20      SUBROUTINE CC_INT3O2(X3OINT,DSRHF,ISYDRHF,
21     *                     XLAMDH,ISYMH,XLAMDP,ISYMP,
22     *                     WORK,LWORK,IDEL,ISYMD,
23     *                     FAC,LWRDSK,LUO3,FILO3,ITRAN,LSQRAB)
24*----------------------------------------------------------------*
25C   Purpose: To calculate (and write to disc) an integral batch with
26C            three occupied indices for a given delta.
27C
28C    Written by Henrik Koch * Asger Halkier 27/7 - 1995
29C
30C    Modified by Asger Halkier to return integrals in X3OINT as
31C    well as writing them to disc 28/10 - 1995.
32C
33C    Generalized to also calculated integrals (j del i k) where k is barred
34C    thus lamda is not total symmetric.
35C    Ove Christiansen 13-6-1996
36C
37C    Generalized to handle LambdaP and LambdaH matrices of different
38C    symmetry. The symmetry of the distribution is passed from
39C    outside.
40C    LWRDSK = .FALSE. the result is NOT written on file
41C    FAC = 1.0  the integrals are ADDED to previous result of call
42C    (second DGEMM).
43C    ALL deltas are dumped on file for the given ITRAN
44C    to allow later sequencial reading of the whole ITRAN batch
45C    of deltas and transformation to occupied using DGEM
46C    Generalized to handle a full (al bet| space of DSRHF (LSQRAB = true)
47C    Sonia Coriani, 10-03-1999
48C
49#include "implicit.h"
50      PARAMETER(ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0)
51      DIMENSION X3OINT(*),DSRHF(*),WORK(LWORK)
52      DIMENSION XLAMDH(*),XLAMDP(*)
53      CHARACTER FILO3*(*)
54      LOGICAL LSQRAB,LWRDSK
55#include "ccorb.h"
56#include "ccsdinp.h"
57#include "ccsdsym.h"
58#include "cclr.h"
59C
60C-------------------------------------------
61C     Work space allocation 1 * outer loops.
62C-------------------------------------------
63C
64      N3ODMX = 0
65      DO ISYM  = 1, NSYM
66        N3ODMX = MAX(N3ODMX,N3ODEL(ISYM))
67      END DO
68
69      ISALBEJ = ISYDRHF
70      ISYIKJ = MULD2H(ISALBEJ,MULD2H(ISYMP,ISYMH))
71C
72      DO 100 ISYMJ = 1,NSYM
73C
74         ISALBE = MULD2H(ISALBEJ,ISYMJ)
75C
76         DO 110 J = 1,NRHF(ISYMJ)
77C
78C------------------------------------------------------------
79C           Work space allocation 1 * unpacking of integrals
80C           if LSQRAB = .FALSE.
81C------------------------------------------------------------
82C
83            IF (.NOT.LSQRAB) THEN
84              KUNINT = 1
85              KEND1  = KUNINT + N2BST(ISALBE)
86              LWRK1  = LWORK  - KEND1
87C
88              IF (LWRK1 .LT. 0) THEN
89                 CALL QUIT(
90     &              '1-Insufficient work space area in CC_INT3O2')
91              ENDIF
92C
93              KOFF1 = IDSRHF(ISALBE,ISYMJ) + NNBST(ISALBE)*(J-1)+1
94C
95              CALL CCSD_SYMSQ(DSRHF(KOFF1),ISALBE,WORK(KUNINT))
96            ELSE
97              KUNINT = IDSRHFSQ(ISALBE,ISYMJ) + N2BST(ISALBE)*(J-1)+1
98              KEND1 = 1
99            END IF
100C
101            DO 120 ISYMK = 1,NSYM
102C
103C-----------------------------------------------------------------------
104C              Transform remaining AO-indices of integrals to occ. space
105C-----------------------------------------------------------------------
106C
107               ISYMBE = MULD2H(ISYMK,ISYMH)
108               ISYMAL = MULD2H(ISALBE,ISYMBE)
109               ISYMI  = MULD2H(ISYMAL,ISYMP)
110               ISYMIK = MULD2H(ISYMI,ISYMK)
111C
112               KINMA1 = KEND1
113               KEND2  = KINMA1 + NBAS(ISYMAL)*NRHF(ISYMK)
114               LWRK2  = LWORK - KEND2
115C
116               IF (LWRK2 .LT. 0) THEN
117                  CALL QUIT(
118     &              '2-Insufficient work space area in CCINT3O2')
119               ENDIF
120C
121               KOFF1 = KUNINT + IAODIS(ISYMAL,ISYMBE)
122               KOFF2 = IGLMRH(ISYMBE,ISYMK) + 1
123               KOFF3 = KINMA1
124C
125               NTOTA = MAX(NBAS(ISYMAL),1)
126               NTOTB = MAX(NBAS(ISYMBE),1)
127C
128               IF (LSQRAB) THEN
129                  CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMK),
130     &                       NBAS(ISYMBE),ONE,DSRHF(KOFF1),NTOTA,
131     &                       XLAMDH(KOFF2),NTOTB,ZERO,WORK(KOFF3),NTOTA)
132               ELSE
133                  CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMK),
134     &                       NBAS(ISYMBE),ONE,WORK(KOFF1),NTOTA,
135     &                       XLAMDH(KOFF2),NTOTB,ZERO,WORK(KOFF3),NTOTA)
136               END IF
137C
138               KOFF1 = IGLMRH(ISYMAL,ISYMI) + 1
139               KOFF2 = KINMA1
140               KOFF3 = IMAIJK(ISYMIK,ISYMJ) + NMATIJ(ISYMIK)*(J - 1)
141     *               + IMATIJ(ISYMI,ISYMK) + 1
142C
143               NTOTA = MAX(NBAS(ISYMAL),1)
144               NTOTI = MAX(NRHF(ISYMI),1)
145C
146               CALL DGEMM('T','N',NRHF(ISYMI),NRHF(ISYMK),NBAS(ISYMAL),
147     *                    ONE,XLAMDP(KOFF1),NTOTA,WORK(KOFF2),NTOTA,
148     *                    FAC,X3OINT(KOFF3),NTOTI)
149C
150  120       CONTINUE
151C
152  110    CONTINUE
153C
154  100 CONTINUE
155C
156      IF ((.NOT. MP2).AND.(LWRDSK)) THEN
157C
158C--------------------------------
159C        Write integrals on disc.
160C--------------------------------
161C
162         D    = IDEL - IBAS(ISYMD)
163C
164         NTOT = NMAIJK(ISYIKJ)
165         IOFF = N3ODMX*(ITRAN - 1) +        !skip previous ITRAN distrib.
166     &          I3ODEL(ISYIKJ,ISYMD) + NTOT*(D - 1) + 1
167C
168         IF (NTOT .GT. 0) THEN
169            CALL PUTWA2(LUO3,FILO3,X3OINT,IOFF,NTOT)
170         ENDIF
171      ENDIF
172C
173      RETURN
174      END
175