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_aofock2 */
20      SUBROUTINE CC_AOFOCK2(XINT,DENSIT,DENSPK,FOCK,WORK,LWORK,IDEL,
21     *                      ISYDIS,ISYMD,ISYDEN,SQRINT)
22C
23C     Purpose: Calculate the two electron contribution to the
24C              AO-fock matrix using matrix vector routines.
25C
26C     Written by Asger Halkier and Henrik Koch 27-4-95.
27C
28C     Debugged Ove Christiansen august 1995
29C
30C     Christof Haettig summer 1998:
31C         arbitrary point group symmetry of the integrals
32C
33C     Sonia Coriani autumn 1999:
34C         Distinguish between packed and squared
35C         integral distributions (for London integrals)
36C         SQRINT = .FALSE., packed integral distribution  and
37C                           Coulomb part calculated with packed density
38C         SQRINT = .TRUE.,  squared integral distribution and
39C                           Coulomb part calculated with squared density
40C
41C     Obs: It can be done as F(g>=d) = G(a>=b) I(a>=b,g,d) where
42C          G(a>=b) = D(a,b) + D(b,a), the diagonal properly scaled
43C
44#include "implicit.h"
45#include "maxorb.h"
46      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, HALF = 0.5D0)
47      DIMENSION XINT(*),DENSIT(*), DENSPK(*)
48      DIMENSION FOCK(*),WORK(LWORK)
49      LOGICAL SQRINT
50#include "ccorb.h"
51#include "ccsdsym.h"
52#include "symsq.h"
53C
54      ISYHOP = MULD2H(ISYMD,ISYDIS)
55C
56C--------------------------------------------
57C     start loop over symmetry blocks:
58C--------------------------------------------
59C
60      DO 100 ISYMG = 1,NSYM
61C
62         IF (NBAS(ISYMG) .EQ. 0) GOTO 100
63C
64         D   = IDEL - IBAS(ISYMD)
65C
66         ISYMAB = MULD2H(ISYMG,ISYDIS)
67
68C
69C---------------------------------------------------------
70C        calculate coulomb contribution:
71C        For SQRINT=.FALSE. use packed form of density and integrals
72C        For SQRINT=.TRUE.  use squared form of density and integrals
73C---------------------------------------------------------
74C
75         IF (ISYMAB .EQ. ISYDEN) THEN
76C
77            KGD = IAODIS(ISYMG,ISYMD) + NBAS(ISYMG)*(D - 1) + 1
78C
79            IF (.NOT.SQRINT) THEN
80               KOFF1 = IDSAOG(ISYMG,ISYDIS) + 1
81               NTOBST = MAX(NNBST(ISYMAB),1)
82               CALL DGEMV('T',NNBST(ISYMAB),NBAS(ISYMG),
83     *                           TWO,XINT(KOFF1),NTOBST,
84     *                           DENSPK,1,ONE,FOCK(KGD),1)
85            ELSE IF (SQRINT) THEN
86               KOFF1  = IDSAOGSQ(ISYMG,ISYDIS) + 1
87               NTOBST = MAX(N2BST(ISYMAB),1)
88               CALL DGEMV('T',N2BST(ISYMAB),NBAS(ISYMG),
89     *                        TWO,XINT(KOFF1),NTOBST,
90     *                        DENSIT,1,ONE,FOCK(KGD),1)
91            END IF
92C
93         END IF
94C
95C-----------------------------------------------------------------
96C        calculate exchange contribution in a loop over g:
97C        For SQRINT=.FALSE. use packed form of density and integrals
98C        For SQRINT=.TRUE.  use squared form of density and integrals
99C-----------------------------------------------------------------
100C
101         IF (LWORK .LT. N2BST(ISYMAB)) THEN
102            CALL QUIT('Insufficient work space in CC_AOFOCK2')
103         ENDIF
104C
105         ISYMA = MULD2H(ISYMD,ISYDEN)
106         ISYMB = MULD2H(ISYMA,ISYMAB)
107C
108         KAD = IAODIS(ISYMA,ISYMD) + NBAS(ISYMA)*(D - 1) + 1
109         NTOTA = MAX(NBAS(ISYMA),1)
110         NTOTG = MAX(NBAS(ISYMG),1)
111C
112         DO G = 1, NBAS(ISYMG)
113C
114            KGB = IAODIS(ISYMG,ISYMB) + G
115
116            IF (.NOT.SQRINT) THEN
117
118               KOFF1 = IDSAOG(ISYMG,ISYDIS) + NNBST(ISYMAB)*(G - 1) + 1
119
120               CALL CCSD_SYMSQ(XINT(KOFF1),ISYMAB,WORK)
121C
122               KAB = IAODIS(ISYMA,ISYMB) + 1
123               CALL DGEMV('T',NBAS(ISYMA),NBAS(ISYMB),-ONE,WORK(KAB),
124     *                 NTOTA,DENSIT(KAD),1,ONE,FOCK(KGB),NTOTG)
125
126            ELSE IF (SQRINT) THEN
127
128               KOFF1 = IDSAOGSQ(ISYMG,ISYDIS) + N2BST(ISYMAB)*(G - 1) +
129     &                 IAODIS(ISYMA,ISYMB) + 1
130
131               CALL DGEMV('T',NBAS(ISYMA),NBAS(ISYMB),-ONE,XINT(KOFF1),
132     *                 NTOTA,DENSIT(KAD),1,ONE,FOCK(KGB),NTOTG)
133
134            END IF
135
136C
137         END DO
138C
139  100 CONTINUE
140C
141      RETURN
142      END
143C  /* Deck cc_dnspk */
144      SUBROUTINE CC_DNSPK(DENSQ,DENSPK,ISYDEN)
145C
146C Purpose: construct lower triangular packed density matrix from
147C          a full square density matrix
148C
149C written by Christof Haettig, July 1998
150C
151
152      use dyn_iadrpk
153
154#include "implicit.h"
155#include "maxorb.h"
156#include "ccorb.h"
157#include "symsq.h"
158#include "ccsdsym.h"
159C
160      DIMENSION DENSQ(*), DENSPK(*)
161C
162C--------------------------------------------
163C        construct triangular density matrix:
164C--------------------------------------------
165C
166      CALL DZERO(DENSPK,NNBST(ISYDEN))
167      DO ISYMB = 1, NSYM
168        ISYMA = MULD2H(ISYDEN,ISYMB)
169        DO B = 1, NBAS(ISYMB)
170          DO A = 1, NBAS(ISYMA)
171            KABSQ = IAODIS(ISYMA,ISYMB) + NBAS(ISYMA)*(B-1) + A
172            KABPK = IADRPK( I2BST(ISYDEN) + KABSQ )
173            DENSPK(KABPK) = DENSPK(KABPK) + DENSQ(KABSQ)
174          END DO
175        END DO
176      END DO
177
178      RETURN
179      END
180