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_aodens3 */
20      SUBROUTINE CC_AODENS3(XLAMDP1,ISYMP1,XLAMDH2,ISYMH2,
21     &                      XLAMDP3,ISYMP3,XLAMDH4,ISYMH4,
22     &                      DENS,ISYDEN,ICORE,IOPT,WORK,LWORK)
23C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
24C
25C     Calculate (special) AO-density matrix used in contructing
26C     the (special) AO Fock matrix.
27C
28C     IOPT = 1 : XLAMDP1 x XLAMDH2
29C     IOPT = 2 : XLAMDP1 x XLAMDH2 + XLAMDP3 x XLAMDH4
30C
31C     Sonia Coriani  10-Feb-1999, based on AODENS2
32C
33C     Careful: No special treatment of ICORE.
34C     ORDER IS IMPORTANT: always P1*H2 (+ P3*H4)
35C     Debug 12.8.99 OK
36C     Reinsert in AODENS2 at some point
37C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
38#include "implicit.h"
39      PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
40      DIMENSION XLAMDP1(*), XLAMDH2(*), XLAMDP3(*), XLAMDH4(*)
41      DIMENSION DENS(*), WORK(LWORK)
42#include "ccorb.h"
43#include "ccsdinp.h"
44#include "ccsdsym.h"
45#include "priunit.h"
46#include "dummy.h"
47C
48*
49* Symmetry test
50*
51      ISYTOT1 = MULD2H(ISYMP1,ISYMH2)
52      IF (ISYTOT1.NE.ISYDEN)
53     *    CALL QUIT('Symmetry mismatch 1 in AODENS3' )
54      IF (IOPT.EQ.2) THEN
55        ISYTOT2 = MULD2H(ISYMP3,ISYMH4)
56        IF (ISYTOT1.NE.ISYTOT2)
57     *      CALL QUIT('Symmetry mismatch 1 in AODENS3' )
58      END IF
59*
60      IF (IOPT.GE.1) THEN
61         CALL DZERO(DENS,N2BST(ISYDEN))
62      END IF
63*
64      DO ISYMK = 1,NSYM
65C
66        IF (IOPT.GE.1) THEN
67C
68            ISYMA = MULD2H(ISYMP1,ISYMK)      ! XLAMDP1
69            ISYMB = MULD2H(ISYMH2,ISYMK)      ! XLAMDH2
70C
71            KOFF1 = 1 + IGLMRH(ISYMA,ISYMK)   ! XLAMDP1
72            KOFF2 = 1 + IGLMRH(ISYMB,ISYMK)   ! XLAMDH2
73            KOFF3 = 1 + IAODIS(ISYMA,ISYMB)
74            NBASA = MAX(NBAS(ISYMA),1)
75            NBASB = MAX(NBAS(ISYMB),1)
76C
77            CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYMK),ONE,
78     *                 XLAMDP1(KOFF1),NBASA,XLAMDH2(KOFF2),NBASB,ONE,
79     *                 DENS(KOFF3),NBASA)
80C
81         END IF
82C
83         IF (IOPT.EQ.2) THEN                     !2nd contribution
84C
85            ISYMA = MULD2H(ISYMP3,ISYMK)         ! XLAMDP3
86            ISYMB = MULD2H(ISYMH4,ISYMK)         ! XLAMDH4
87C
88            KOFF1 = 1 + IGLMRH(ISYMA,ISYMK)      ! XLAMDP3
89            KOFF2 = 1 + IGLMRH(ISYMB,ISYMK)      ! XLAMDH4
90            KOFF3 = 1 + IAODIS(ISYMA,ISYMB)
91            NBASA = MAX(NBAS(ISYMA),1)
92            NBASB = MAX(NBAS(ISYMB),1)
93C
94            CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYMK),ONE,
95     *                 XLAMDP3(KOFF1),NBASA,XLAMDH4(KOFF2),NBASB,ONE,
96     *                 DENS(KOFF3),NBASA)
97C
98         END IF !result has been added to previous (beta=1)
99C
100      END DO
101C
102C
103C-----------------------------
104C     Include frozen orbitals.
105C-----------------------------
106C
107      IF ( (FROIMP.OR.FROEXP) .AND. (ICORE .EQ. 1) ) THEN
108C
109         IF (IOPT.NE.0) THEN
110           WRITE (LUPRI,*)
111     *              'CC_AODENS3: ICORE=1 not yet available for IOPT>0.'
112           CALL QUIT(
113     *              'CC_AODENS3: ICORE=1 not yet available for IOPT>0.')
114         END IF
115C
116         IF (LWORK .LT. NLAMDS) THEN
117            CALL QUIT('Insufficient space in CCSD_AODENS')
118         ENDIF
119C
120C-------------------------------------------------
121C        Read MO-coefficients from interface file.
122C-------------------------------------------------
123C
124         LUSIFC = -1
125         CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',
126     *               IDUMMY,.FALSE.)
127         REWIND LUSIFC
128C
129         CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI)
130         READ (LUSIFC)
131C
132         READ (LUSIFC)
133         READ (LUSIFC) (WORK(I), I=1,NLAMDS)
134C
135         CALL GPCLOSE(LUSIFC,'KEEP')
136C
137C-------------------------------------------------------
138C        Add contribution from frozen occupied orbitals.
139C-------------------------------------------------------
140C
141         KOFF1 = 0
142         KOFF2 = 0
143         DO ISYMK = 1,NSYM
144C
145CCH the following has been changed because ISYMP0 is not defined...?
146CCH         ISYMA = MULD2H(ISYMP0,ISYMK)
147CCH         ISYMB = MULD2H(ISYMP0,ISYMK)
148CCH
149CCH Sonia is this correct?
150            ISYMA = ISYMK
151            ISYMB = ISYMK
152CCH
153C
154            DO II = 1,NRHFFR(ISYMK)
155C
156               K = KFRRHF(II,ISYMK)
157C
158               DO B = 1,NBAS(ISYMB)
159                  DO A = 1,NBAS(ISYMA)
160C
161                     NAK = KOFF1 + NBAS(ISYMA)*(K - 1) + A
162                     NBK = KOFF1 + NBAS(ISYMB)*(K - 1) + B
163                     NAB = KOFF2 + NBAS(ISYMA)*(B - 1) + A
164C
165                     DENS(NAB) = DENS(NAB) + WORK(NAK)*WORK(NBK)
166C
167                  END DO
168               END DO
169C
170            END DO
171C
172            KOFF1 = KOFF1 + NBAS(ISYMK)*NORBS(ISYMK)
173            KOFF2 = KOFF2 + NBAS(ISYMA)*NBAS(ISYMB)
174C
175         END DO
176C
177      ENDIF
178C
179      END
180