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 eriave */
20      SUBROUTINE ERIAVE(SO,PMAT,DMAT,D2MAT,ID2MAT,IPNTCR,IODDCC,IPNTUV,
21     &                  INDXBT,WORK,LWORK,IPRINT)
22C
23C     tuh march 99 - april 01
24C
25#include "implicit.h"
26#include "priunit.h"
27#include "iratdef.h"
28C
29      DIMENSION SO(*), PMAT(*), DMAT(*), D2MAT(*), ID2MAT(*),
30     &          IPNTCR(*), IPNTUV(*), IODDCC(*),
31     &          INDXBT(*), WORK(*)
32#include "ericom.h"
33C
34C     Allocations
35C
36      KPOINT = 1
37      KLAST  = KPOINT + (5*NCCS - 1)/IRAT + 1
38      IF (KLAST .GT. LWORK) CALL STOPIT('ERIAVE',' ',KLAST,LWORK)
39      CALL ERIAV1(SO,DMAT,D2MAT,ID2MAT,IPNTCR,IODDCC,IPNTUV,
40     &            INDXBT,WORK(KPOINT),IPRINT)
41      RETURN
42      END
43C  /* Deck eriav1 */
44      SUBROUTINE ERIAV1(SO,DMAT,D2MAT,ID2MAT,IPNTCR,IODDCC,IPNTUV,
45     &                  INDXBT,IPOINT,IPRINT)
46#include "implicit.h"
47#include "priunit.h"
48#include "iratdef.h"
49#include "maxaqn.h"
50#include "mxcent.h"
51#include "maxorb.h"
52#include "aovec.h"
53      PARAMETER (D1=1.0D0, DP5=0.5D0)
54      INTEGER A, B, C, D, R, S, T
55      LOGICAL DOREP(0:7,4), CTRIAB, CTRICD, DAB, DCD
56      DIMENSION IPOINT(NCCS,5),
57     &          IPNTCR(MAXBCH,4), INDXBT(MXSHEL*MXCONT,0:7),
58     &          IODDCC(NRTOP), IPNTUV(KC2MAX,0:NRDER,2),
59     &          IPNRST(0:7,3), IODXYZ(3),
60     &          IADCMP(MXAQN,MXAQN,2),
61     &          SO(NCCS,MLTPR,MLTPS,MLTPT,KHKTAB,KHKTCD,12),
62     &          DMAT(NBASE,NBASE),
63     &          D2MAT(NBASE,NBASE,NPDIMB,NPDIMA),
64     &          ID2MAT(MXCORB,2)
65#include "cbieri.h"
66#include "ericom.h"
67#include "erithr.h"
68#include "eribuf.h"
69#include "aobtch.h"
70#include "hertop.h"
71#include "symmet.h"
72#include "nuclei.h"
73#include "energy.h"
74C
75
76      IBTEST(I,J,K,L) = IAND(I,IEOR(J,ISYMAO(K,L)))
77C
78      IF (IPRINT .GT. 6) CALL HEADER('Subroutine ERIAV1',-1)
79C
80      IF (GRDZER) CALL DZERO(GRADEE,MXCOOR)
81C
82      CALL PRPREP(DOREP(0,1),NHKTA,KHKTA,ISTBLA)
83      CALL PRPREP(DOREP(0,2),NHKTB,KHKTB,ISTBLB)
84      CALL PRPREP(DOREP(0,3),NHKTC,KHKTC,ISTBLC)
85      CALL PRPREP(DOREP(0,4),NHKTD,KHKTD,ISTBLD)
86C
87      IF (IPRINT .GT. 10) THEN
88         WRITE (LUPRI,'(/,2X,A,8L2)')'DOREP A  ',(DOREP(I,1),I=0,MAXREP)
89         WRITE (LUPRI,'(2X,A,8L2)')  'DOREP B  ',(DOREP(I,2),I=0,MAXREP)
90         WRITE (LUPRI,'(2X,A,8L2)')  'DOREP C  ',(DOREP(I,3),I=0,MAXREP)
91         WRITE (LUPRI,'(2X,A,8L2)')  'DOREP D  ',(DOREP(I,4),I=0,MAXREP)
92      END IF
93C
94      CALL CMPADR(IADCMP(1,1,1),KHKTA,KHKTB,TKMPAB)
95      CALL CMPADR(IADCMP(1,1,2),KHKTC,KHKTD,TKMPCD)
96C
97      CALL GETRST(IPNRST(0,1),ISTBLR)
98      CALL GETRST(IPNRST(0,2),ISTBLS)
99      CALL GETRST(IPNRST(0,3),ISTBLT)
100C
101      DO A = 0, MAXREP
102      IF (DOREP(A,1)) THEN
103      DO B = 0, MAXREP
104      IF (DOREP(B,2)) THEN
105      DO C = 0, MAXREP
106      IF (DOREP(C,3) .AND. DOREP(IEOR(IEOR(A,B),C),4)) THEN
107         D = IEOR(IEOR(A,B),C)
108         IF (DIAGAB .AND. B.GT.A) GO TO 100
109         IF (DIAGCD .AND. D.GT.C) GO TO 100
110C
111         CTRIAB = DIAGAB .AND. A.EQ.B
112         CTRICD = DIAGCD .AND. C.EQ.D
113C
114         R = IPNRST(B,1)
115         S = IPNRST(D,2)
116         T = IPNRST(IEOR(C,D),3)
117C
118         CALL ERIPNT(IPOINT,A,B,C,D,IPNTCR,INDXBT,0)
119C
120         IA   = -1
121         MAXB = KHKTB
122         MAXD = KHKTD
123         DO ICMPA = 1, KHKTA
124         IVARA = IBTEST(ISTBLA,A,NHKTA,ICMPA)
125         IF (IVARA.EQ.0) THEN
126            IA = IA + 1
127            IB = -1
128            IF (CTRIAB) MAXB = ICMPA
129            DO ICMPB = 1, MAXB
130            IVARB = IBTEST(ISTBLB,B,NHKTB,ICMPB)
131            IF (IVARB.EQ.0) THEN
132               IB = IB + 1
133               IC = -1
134               ICMPAB = IADCMP(ICMPA,ICMPB,1)
135               IODXYZ(1) = IODDCC(IPNTUV(ICMPAB,1,1))
136               IODXYZ(2) = IODDCC(IPNTUV(ICMPAB,2,1))
137               IODXYZ(3) = IODDCC(IPNTUV(ICMPAB,3,1))
138               DO ICMPC = 1, KHKTC
139               IVARC = IBTEST(ISTBLC,C,NHKTC,ICMPC)
140               IF (IVARC.EQ.0) THEN
141                  IC = IC + 1
142                  ID = -1
143                  IF (CTRICD) MAXD = ICMPC
144                  DO ICMPD = 1, MAXD
145                  IVARD = IBTEST(ISTBLD,D,NHKTD,ICMPD)
146                  IF (IVARD.EQ.0) THEN
147                     ID = ID + 1
148                     ICMPCD = IADCMP(ICMPC,ICMPD,2)
149                     IODDCD = IODDCC(IPNTUV(ICMPCD,0,2))
150                     IF(IODDCD.EQ.IODXYZ(1).OR.IODDCD.EQ.IODXYZ(2)
151     &                                     .OR.IODDCD.EQ.IODXYZ(3))THEN
152                        DAB = GCONAB .AND. CTRIAB.AND.ICMPA.EQ.ICMPB
153                        DCD = GCONCD .AND. CTRICD.AND.ICMPC.EQ.ICMPD
154                        FAC = D1
155                        IF (DIAGPQ) FAC = DP5*FAC
156                        CALL ERIFRC(IA,IB,IC,ID,
157     &                              SO(1,R,S,T,ICMPAB,ICMPCD,1),
158     &                              IPOINT,DMAT,D2MAT,ID2MAT,FAC,
159     &                              IODDCD,IODXYZ,DAB,DCD)
160                     END IF
161                  END IF
162                  END DO
163               END IF
164               END DO
165            END IF
166            END DO
167         END IF
168         END DO
169  100    CONTINUE
170      END IF
171      END DO
172      END IF
173      END DO
174      END IF
175      END DO
176C
177      IF (IPRINT .GE. 10) THEN
178         CALL HEADER('GRADEE in ERIAV1',-1)
179         CALL OUTPUT(GRADEE,1,3,1,NUCDEP,3,NUCDEP,1,LUPRI)
180      END IF
181C
182      RETURN
183      END
184C  /* Deck erifrc */
185      SUBROUTINE ERIFRC(IA,IB,IC,ID,SO,IPOINT,DMAT,D2MAT,ID2MAT,
186     &                  DFAC,IODDCD,IODXYZ,DAB,DCD)
187#include "implicit.h"
188#include "priunit.h"
189#include "maxaqn.h"
190#include "mxcent.h"
191#include "maxorb.h"
192#include "aovec.h"
193      PARAMETER (D4=4.0D0, D2=2.0D0, D1=1.0D0, DP5=0.5D0, DP25=0.25D0)
194      LOGICAL DAB, DCD, NOTRAS
195      DIMENSION SO(NAOINT,3,4), IPOINT(NCCS,5), DMAT(NBASE,NBASE),
196     &          D2MAT(NBASE,NBASE,NPDIMB,NPDIMA), ID2MAT(MXCORB_CC,2),
197     &          IODXYZ(3)
198#include "cbieri.h"
199#include "ericom.h"
200#include "eribuf.h"
201#include "aobtch.h"
202#include "hertop.h"
203#include "nuclei.h"
204#include "symmet.h"
205#include "energy.h"
206C
207      NOTRAS = .TRUE.
208      DO I = 1, NCCS
209         KA = IPOINT(I,1) + IA
210         KB = IPOINT(I,2) + IB
211         KC = IPOINT(I,3) + IC
212         KD = IPOINT(I,4) + ID
213         IF (CCRUN) THEN
214            PMAT = DFAC*D2MAT(KD,KC,ID2MAT(KB,2),ID2MAT(KA,1))
215            IF (DCD .OR. KC.EQ.KD) PMAT = DP5*PMAT
216         ELSE
217            PMAT = DFAC*D4*(DMAT(KA,KB)*DMAT(KC,KD)
218     &               - DP25*DMAT(KA,KC)*DMAT(KB,KD)
219     &               - DP25*DMAT(KA,KD)*DMAT(KB,KC))
220            IF (DAB .OR. KA .EQ. KB) PMAT = DP5*PMAT
221            IF (DCD .OR. KC .EQ. KD) PMAT = DP5*PMAT
222         END IF
223C
224         JA = 3*ICNTAO(KA) - 3
225         JB = 3*ICNTAO(KB) - 3
226         JC = 3*ICNTAO(KC) - 3
227         JD = 3*ICNTAO(KD) - 3
228C
229         IF (NOTRAS) THEN
230            J = IPOINT(I,5)
231            DO ICOOR = 1, 3
232            IF (IODDCD .EQ. IODXYZ(ICOOR)) THEN
233               ISCORA = IPTCNT(JA+ICOOR,0,1)
234               ISCORB = IPTCNT(JB+ICOOR,0,1)
235               ISCORC = IPTCNT(JC+ICOOR,0,1)
236               ISCORD = IPTCNT(JD+ICOOR,0,1)
237               IF (ISCORA.GT.0) GRADEE(ISCORA)
238     &                        = GRADEE(ISCORA) + PMAT*SO(J,ICOOR,1)
239               IF (ISCORB.GT.0) GRADEE(ISCORB)
240     &                        = GRADEE(ISCORB) + PMAT*SO(J,ICOOR,2)
241               IF (ISCORC.GT.0) GRADEE(ISCORC)
242     &                        = GRADEE(ISCORC) + PMAT*SO(I,ICOOR,3)
243               IF (ISCORD.GT.0) GRADEE(ISCORD)
244     &                        = GRADEE(ISCORD) + PMAT*SO(I,ICOOR,4)
245            END IF
246            END DO
247         ELSE
248            DO ICOOR = 1, 3
249            IF (IODDCD .EQ. IODXYZ(ICOOR)) THEN
250               FA = PMAT*SO(J,ICOOR,1)
251               FB = PMAT*SO(I,ICOOR,2)
252               FC = PMAT*SO(I,ICOOR,3)
253               FD = - FA - FB - FC
254               ISCORA = IPTCNT(JA+ICOOR,0,1)
255               ISCORB = IPTCNT(JB+ICOOR,0,1)
256               ISCORC = IPTCNT(JC+ICOOR,0,1)
257               ISCORD = IPTCNT(JD+ICOOR,0,1)
258               IF (ISCORA.GT.0) GRADEE(ISCORA) = GRADEE(ISCORA) + FA
259               IF (ISCORB.GT.0) GRADEE(ISCORB) = GRADEE(ISCORB) + FB
260               IF (ISCORC.GT.0) GRADEE(ISCORC) = GRADEE(ISCORC) + FC
261               IF (ISCORD.GT.0) GRADEE(ISCORD) = GRADEE(ISCORD) + FD
262            END IF
263            END DO
264         END IF
265      END DO
266      RETURN
267      END
268