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 fonedn */
20      SUBROUTINE FONEDN(NSIM,FCOCO,FVOCO,FCOEX,FVOEX)
21C 1)
22C SCALE FCOCO AND FVOEX TO HAVE RIGHT FACTORS
23C 2)
24C ADD CONTRIBUTIONS FROM OCCUPIED SECONDARY TO FCOCO AND FVOCO
25C (USE THAT COULOMB CONTRIBUTION IS SYMMETRIC) AND
26C 3)
27C SUM UP ONE INDEX TRANSFORMED DENSITY CONTRIBUTIONS TO FOCK
28C MATRICES ( FCOEX  = FCOCO  + FCOEX )
29C
30#include "implicit.h"
31C
32#include "maxorb.h"
33#include "maxash.h"
34#include "infinp.h"
35#include "inforb.h"
36#include "infind.h"
37#include "dftcom.h"
38#include "wrkrsp.h"
39#include "infrsp.h"
40C
41      DIMENSION FCOCO(NORBT,NORBT,*),FVOCO(NORBT,NORBT,*)
42      DIMENSION FCOEX(NORBT,NORBT,*),FVOEX(NORBT,NORBT,*)
43C
44      PARAMETER ( D2 = 2.0D0  , DP5 = 0.5D0 )
45C
46      IF (NASHT.GT.0) THEN
47         FACTOR = DP5 * HFXFAC
48         CALL DSCAL(NORBT*NORBT*NSIM,FACTOR,FVOEX,1)
49      END IF
50      IF (HFXFAC .NE. 1.0D0) THEN
51         CALL DSCAL(NORBT*NORBT*NSIM,HFXFAC,FCOEX,1)
52      END IF
53      IF (TRPLET) RETURN
54      IF (NISHT.GT.0) CALL DSCAL(NORBT*NORBT*NSIM,D2,FCOCO,1)
55      DO 100 IPSYM = 1,NSYM
56         IORBP = IORB(IPSYM) + 1
57         NOCCP = NOCC(IPSYM)
58         IQSYM = MULD2H(IPSYM,KSYMOP)
59         IORBQ = IORB(IQSYM) + 1
60         NOCCQ = NOCC(IQSYM)
61         NORBQ = NORB(IQSYM)
62         IF ((NOCCP.EQ.0).OR.((NORBQ-NOCCQ).EQ.0)) GO TO 100
63         DO 400 ISIM = 1,NSIM
64            DO 200 IP = IORBP,IORBP+NOCCP-1
65               DO 300 IQ = IORBQ+NOCCQ,IORBQ+NORBQ-1
66                  FCOCO(IP,IQ,ISIM) = FCOCO(IQ,IP,ISIM)
67                  IF (NASHT.GT.0) FVOCO(IP,IQ,ISIM) = FVOCO(IQ,IP,ISIM)
68 300           CONTINUE
69 200        CONTINUE
70 400     CONTINUE
71 100  CONTINUE
72      DO 500 IPSYM = 1,NSYM
73         IORBP = IORB(IPSYM)
74         NORBP = NORB(IPSYM)
75         IQSYM = MULD2H(IPSYM,KSYMOP)
76         IORBQ = IORB(IQSYM)
77         NORBQ = NORB(IQSYM)
78         IF ((NORBQ.EQ.0).OR.(NORBP.EQ.0)) GO TO 500
79         DO ISIM = 1,NSIM
80            DO IQ = IORBQ+1,IORBQ+NORBQ
81               DO IP = IORBP+1,IORBP+NORBP
82                  FCOEX(IP,IQ,ISIM) = FCOCO(IP,IQ,ISIM)
83     *                              + FCOEX(IP,IQ,ISIM)
84                  IF (NASHT.GT.0) FVOEX(IP,IQ,ISIM) = FVOCO(IP,IQ,ISIM)
85     *                              + FVOEX(IP,IQ,ISIM)
86               END DO
87            END DO
88         END DO
89 500  CONTINUE
90      RETURN
91      END
92C  /* Deck fonemu */
93      SUBROUTINE FONEMU(NSIM,ICI1,IDI1,H2,
94     *                  FCOEX,FVOEX,ZYMAT,DENA,DENB,WRK,LWRK)
95C
96C CALCULATE CONTRIBUTIONS TO INACTIVE AND ACTIVE FOCK MATRICES
97C WHICH ORIGINATE FROM OCCUPIED-OCCUPIED COULOMN DISTRIBUTIONS
98C (ONLY THE EXCHANGE PART OF THE ONE-INDEX TRANSFORMED
99C  DENSITY CONTRIBUTE)
100C
101C  FCOEX(P,Q) = FCOEX(P,Q) + SUM(R) (JQ/PR)*ZYMAT(R,J)  (Q OCC)
102C                          - SUM(R) (PJ/RQ)*ZYMAT(J,R)  (P OCC)
103C
104C  FVOEX(P,Q) = FVOEX(P,Q) + SUM(R) (XQ/PR)*DENA(R,X)   (Q OCC)
105C                          - SUM(R) (PY/QR)*DENB(R,Y)   (P OCC)
106C
107#include "implicit.h"
108C
109C Used from common blocks:
110C   INFORB :
111C   INFIND : ISMO,IOBTYP
112#include "maxorb.h"
113#include "maxash.h"
114#include "inforb.h"
115#include "infind.h"
116#include "wrkrsp.h"
117#include "infrsp.h"
118#include "orbtypdef.h"
119C
120      DIMENSION H2(NORBT,*)
121      DIMENSION FCOEX(NORBT,NORBT,*),FVOEX(NORBT,NORBT,*)
122      DIMENSION ZYMAT(NORBT,NORBT,*)
123      DIMENSION DENA(NORBT,NASHT,*),DENB(NORBT,NASHT,*),WRK(*)
124C
125C     Order (C,D) index such that C .ge. D
126C     in inactive-active-secondary order (using ISW)
127C
128      IF (ISW(ICI1) .GE. ISW(IDI1)) THEN
129         ICI = ICI1
130         IDI = IDI1
131      ELSE
132         ICI = IDI1
133         IDI = ICI1
134      END IF
135C
136C     Find distribution type ITYPCD =
137C     1:inactive-inactive  2:active-inactive  3:active-active
138C     4:secondary-inactive 5:secondary-active 6:secondary-secondary
139C     We only need occupied-occupied distributions, i.e. itypcd .le. 3
140C
141      ITYPC  = IOBTYP(ICI)
142      ITYPD  = IOBTYP(IDI)
143      ITYPCD = IDBTYP(ITYPC,ITYPD)
144      IF (ITYPCD .GT. 3) RETURN
145      IF (ITYPC .EQ. JTACT) NCIW = ICH(ICI)
146      IF (ITYPD .EQ. JTACT) NDIW = ICH(IDI)
147C
148      ICSYM = ISMO(ICI)
149      IDSYM = ISMO(IDI)
150      IOFFC = IORB(ICSYM) + 1
151      IOFFD = IORB(IDSYM) + 1
152      ICDSYM = MULD2H(ICSYM,IDSYM)
153      DO 100 ISIM = 1,NSIM
154         DO 200 IRSYM = 1,NSYM
155            IPSYM = MULD2H(IRSYM,ICDSYM)
156            IOFFR = IORB(IRSYM) + 1
157            IOFFP = IORB(IPSYM) + 1
158            NORBR = NORB(IRSYM)
159            NORBP = NORB(IPSYM)
160            IF ( (NORBR.EQ.0) .OR. (NORBP.EQ.0) ) GO TO 200
161C
162            ICRSYM = MULD2H(ICSYM,IRSYM)
163            IDRSYM = MULD2H(IDSYM,IRSYM)
164            IF ( KSYMOP.EQ.IDRSYM ) THEN
165C
166               IF ((ITYPCD.EQ.1) .OR. (ITYPCD.EQ.2)) THEN
167C
168C  FCOEX(P,Q) = FCOEX(P,Q) + SUM(R) (JQ/PR)*ZYMAT(R,J)  (Q OCC)
169C                                    QJ/PR
170                  CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
171     &                       H2(IOFFP,IOFFR),NORBT,
172     &                       ZYMAT(IOFFR,IDI,ISIM),NORBT,1.D0,
173     &                       FCOEX(IOFFP,ICI,ISIM),NORBT)
174C
175C                          - SUM(R) (PJ/RQ)*ZYMAT(J,R)  (P OCC)
176C
177                  CALL DGEMM('N','N',1,NORBP,NORBR,-1.D0,
178     &                       ZYMAT(IDI,IOFFR,ISIM),NORBT,
179     &                       H2(IOFFR,IOFFP),NORBT,1.D0,
180     &                       FCOEX(ICI,IOFFP,ISIM),NORBT)
181               END IF
182            END IF
183            IF (ICRSYM.EQ.KSYMOP) THEN
184               IF ((ITYPCD.EQ.1) .AND. (ICI.NE.IDI)) THEN
185C
186C  FCOEX(P,Q) = FCOEX(P,Q) + SUM(R) (JQ/PR)*ZYMAT(R,J)  (Q OCC)
187C                                    QJ/PR
188                  CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
189     &                       H2(IOFFP,IOFFR),NORBT,
190     &                       ZYMAT(IOFFR,ICI,ISIM),NORBT,1.D0,
191     &                       FCOEX(IOFFP,IDI,ISIM),NORBT)
192C
193C                          - SUM(R) (PJ/RQ)*ZYMAT(J,R)  (P OCC)
194C
195                  CALL DGEMM('N','N',1,NORBP,NORBR,-1.D0,
196     &                       ZYMAT(ICI,IOFFR,ISIM),NORBT,
197     &                       H2(IOFFR,IOFFP),NORBT,1.D0,
198     &                       FCOEX(IDI,IOFFP,ISIM),NORBT)
199               END IF
200            END IF
201            IF ( ICRSYM.EQ.KSYMOP ) THEN
202               IF ((ITYPCD.EQ.2).OR.(ITYPCD.EQ.3)) THEN
203C
204C  FVOEX(P,Q) = FVOEX(P,Q) + SUM(R) (XQ/PR)*DENA(R,X)   (Q OCC)
205C
206                  CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
207     &                       H2(IOFFP,IOFFR),NORBT,
208     &                       DENA(IOFFR,NCIW,ISIM),NORBT,1.D0,
209     &                       FVOEX(IOFFP,IDI,ISIM),NORBT)
210C
211C                          - SUM(R) (PY/QR)*DENB(R,Y)   (P OCC)
212C                                    YP/QR
213C
214                  CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
215     &                       H2(IOFFP,IOFFR),NORBT,
216     &                       DENB(IOFFR,NCIW,ISIM),NORBT,0.D0,
217     &                       WRK(1),NORBT)
218                  DO 450 IQ = IOFFP,IOFFP+NORBP-1
219                     FVOEX(IDI,IQ,ISIM) = FVOEX(IDI,IQ,ISIM)
220     *                                       - WRK(IQ-IOFFP+1)
221 450              CONTINUE
222               END IF
223            END IF
224            IF (IDRSYM.EQ.KSYMOP) THEN
225               IF ((ITYPCD.EQ.3).AND.(ICI.NE.IDI)) THEN
226C
227C  FVOEX(P,Q) = FVOEX(P,Q) + SUM(R) (XQ/PR)*DENA(R,X)   (Q OCC)
228C
229                  CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
230     &                       H2(IOFFP,IOFFR),NORBT,
231     &                       DENA(IOFFR,NDIW,ISIM),NORBT,1.D0,
232     &                       FVOEX(IOFFP,ICI,ISIM),NORBT)
233C
234C                          - SUM(R) (PY/QR)*DENB(R,Y)   (P OCC)
235C                                    YP/QR
236C
237                  CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
238     &                       H2(IOFFP,IOFFR),NORBT,
239     &                       DENB(IOFFR,NDIW,ISIM),NORBT,0.D0,
240     &                       WRK(1),NORBT)
241                  DO 451 IQ = IOFFP,IOFFP+NORBP-1
242                     FVOEX(ICI,IQ,ISIM) = FVOEX(ICI,IQ,ISIM)
243     *                                       - WRK(IQ-IOFFP+1)
244 451              CONTINUE
245               END IF
246            END IF
247 200     CONTINUE
248 100  CONTINUE
249      RETURN
250      END
251C  /* Deck fonedr */
252      SUBROUTINE FONEDR(NSIM,ICI1,IDI1,H2D,FCOCO,FVOCO,
253     *                  FCOEX,FVOEX,ZYMAT,DENA,DENB,WRK,LWRK)
254C
255C CALCULATE CONTRIBUTIONS TO INACTIVE AND ACTIVE FOCK MATRICES
256C WHICH ORIGINATE FROM OCCUPIED-OCCUPIED DIRAC DISTRIBUTIONS
257C
258C F(CV)OEX CONTAIN EXCHANCE CONTRIBUTION
259C F(CV)OMU CONTAIN COULOMN CONTRIBUTIONS
260C
261C  FCOEX(P,Q) = FCOEX(P,Q) + SUM(R) <JP/QR>*ZYMAT(R,J) P OCC Q SEC
262C                          - SUM(R) <QJ/RP>*ZYMAT(J,R) Q OCC P SEC
263C
264C  FVOEX(P,Q) = FVOEX(P,Q) + SUM(R) <XP/QR>*DENA(R,X)  P OCC Q SEC
265C                          - SUM(R) <QY/RP>*DENB(R,Y)  Q OCC P SEC
266C
267C  FCOCO(P,Q) = FCOCO(P,Q) + SUM(R) <QJ/PR>*(ZYMAT(J,R)-ZYMAT(R,J))
268C
269C  FVOCO(P,Q) = FVOCO(P,Q) + SUM(R) <QY/PR>*(DENB(R,Y)-DENA(R,Y))
270C
271#include "implicit.h"
272C
273C Used from common blocks:
274C   INFORB :
275C   INFIND : ISMO,IOBTYP
276#include "maxorb.h"
277#include "maxash.h"
278#include "inforb.h"
279#include "infind.h"
280#include "infrsp.h"
281#include "wrkrsp.h"
282#include "orbtypdef.h"
283      DIMENSION H2D(NORBT,*)
284      DIMENSION FCOCO(NORBT,NORBT,*),FVOCO(NORBT,NORBT,*)
285      DIMENSION FCOEX(NORBT,NORBT,*),FVOEX(NORBT,NORBT,*)
286      DIMENSION ZYMAT(NORBT,NORBT,*)
287      DIMENSION DENA(NORBT,NASHT,*),DENB(NORBT,NASHT,*),WRK(*)
288C
289C     Order (C,D) index such that C .ge. D
290C     in inactive-active-secondary order (using ISW)
291C
292      NDITR = 1
293      IF (ISW(ICI1) .GE. ISW(IDI1)) THEN
294         ICI = ICI1
295         IDI = IDI1
296      ELSE
297         ICI = IDI1
298         IDI = ICI1
299         CALL DGETRN(H2D,NORBT,NORBT)
300         NDITR = -NDITR
301      END IF
302C
303C     Find distribution type ITYPCD =
304C     1:inactive-inactive  2:active-inactive  3:active-active
305C     4:secondary-inactive 5:secondary-active 6:secondary-secondary
306C     We only need occupied-occupied distributions, i.e. itypcd .le. 3
307C
308      ITYPC  = IOBTYP(ICI)
309      ITYPD  = IOBTYP(IDI)
310      ITYPCD = IDBTYP(ITYPC,ITYPD)
311      IF (ITYPCD .GT. 3) RETURN
312      IF (ITYPC .EQ. JTACT) NCIW = ICH(ICI)
313      IF (ITYPD .EQ. JTACT) NDIW = ICH(IDI)
314C
315      ICSYM = ISMO(ICI)
316      IDSYM = ISMO(IDI)
317      IOFFC = IORB(ICSYM) + 1
318      IOFFD = IORB(IDSYM) + 1
319      ICDSYM = MULD2H(ICSYM,IDSYM)
320      DO 100 ISIM = 1,NSIM
321         DO 200 IRSYM = 1,NSYM
322            IPSYM = MULD2H(IRSYM,ICDSYM)
323            NORBR = NORB(IRSYM)
324            NSSHP = NSSH(IPSYM)
325            IF (NORBR.EQ.0) GO TO 200
326            IOFFR = IORB(IRSYM) + 1
327            IOFFP = IORB(IPSYM) + 1
328            NORBP = NORB(IPSYM)
329            IF (NORBP.EQ.0) GO TO 200
330            NOCCP = NOCC(IPSYM)
331            IDRSYM = MULD2H(IDSYM,IRSYM)
332            ICRSYM = MULD2H(ICSYM,IRSYM)
333            IF ( IDRSYM.EQ.KSYMOP) THEN
334               IF ((.NOT.TRPLET).AND.(ITYPCD.LE.2)) THEN
335C
336C  FCOCO(P,Q) = FCOCO(P,Q) + SUM(R) <QJ/PR>*(ZYMAT(J,R)-ZYMAT(R,J))
337C
338                  DO 1100 IR=1,NORBR
339                     WRK(IR) = ZYMAT(IDI,IOFFR-1+IR,ISIM)-
340     *                         ZYMAT(IOFFR-1+IR,IDI,ISIM)
341 1100             CONTINUE
342                  CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
343     &                       H2D(IOFFP,IOFFR),NORBT,
344     &                       WRK,NORBT,1.D0,
345     &                       FCOCO(IOFFP,ICI,ISIM),NORBT)
346               END IF
347            END IF
348            IF ( ICRSYM.EQ.KSYMOP) THEN
349               IF ((.NOT.TRPLET).AND.(ITYPCD.GE.2)) THEN
350C
351C  FVOCO(P,Q) = FVOCO(P,Q) + SUM(R) <QY/PR>*(DENB(R,Y)-DENA(R,Y))
352C                                    YQ/RP        R,Y       R,Y
353C
354                  DO 1200 IR=1,NORBR
355                     WRK(IR) = DENB(IOFFR-1+IR,NCIW,ISIM)-
356     *                         DENA(IOFFR-1+IR,NCIW,ISIM)
357 1200             CONTINUE
358                  CALL DGEMM('T','N',NORBP,1,NORBR,1.D0,
359     &                       H2D(IOFFR,IOFFP),NORBT,
360     &                       WRK,NORBT,1.D0,
361     &                       FVOCO(IOFFP,IDI,ISIM),NORBT)
362               END IF
363            END IF
364            IF (NSSHP.EQ.0) GO TO 200
365            IF ( IDRSYM.EQ.KSYMOP) THEN
366               IF ((ITYPCD.EQ.1) .OR. (ITYPCD.EQ.2)) THEN
367C
368C  FCOEX(P,Q) = FCOEX(P,Q) - SUM(R) <QJ/RP>*ZYMAT(J,R) Q OCC P SEC
369C
370
371                  DO 111 IR=1,NORBR
372                     WRK(IR) =-ZYMAT(IDI,IOFFR-1+IR,ISIM)
373 111              CONTINUE
374                  CALL DGEMM('T','N',NSSHP,1,NORBR,1.D0,
375     &                       H2D(IOFFR,IOFFP+NOCCP),NORBT,
376     &                       WRK(1),NORBR,1.D0,
377     &                       FCOEX(IOFFP+NOCCP,ICI,ISIM),NORBT)
378C
379C  FCOEX(P,Q) = FCOEX(P,Q) + SUM(R) <JP/QR>*ZYMAT(R,J) P OCC Q SEC
380C                                    PJ RQ
381C
382                  CALL DGEMM('T','N',1,NSSHP,NORBR,1.D0,
383     &                       ZYMAT(IOFFR,IDI,ISIM),NORBT,
384     &                       H2D(IOFFR,IOFFP+NOCCP),NORBT,1.D0,
385     &                       FCOEX(ICI,IOFFP+NOCCP,ISIM),NORBT)
386               END IF
387            END IF
388            IF ( ICRSYM.EQ.KSYMOP) THEN
389               IF ((ITYPCD.EQ.2).OR.(ITYPCD.EQ.3)) THEN
390C
391C  FVOEX(P,Q) = FVOEX(P,Q) - SUM(R) <QY/RP>*DENB(R,Y) Q OCC P SEC
392C                                    YQ/PR       R,Y
393                  CALL DGEMM('N','N',NSSHP,1,NORBR,-1.D0,
394     &                       H2D(IOFFP+NOCCP,IOFFR),NORBT,
395     &                       DENB(IOFFR,NCIW,ISIM),NORBT,1.D0,
396     &                       FVOEX(IOFFP+NOCCP,IDI,ISIM),NORBT)
397C
398C  FVOEX(P,Q) = FVOEX(P,Q) + SUM(R) <XP/QR>*DENA(R,X)   (P OCC)
399C
400                  CALL DGEMM('N','N',NSSHP,1,NORBR,1.D0,
401     &                       H2D(IOFFP+NOCCP,IOFFR),NORBT,
402     &                       DENA(IOFFR,NCIW,ISIM),NORBT,0.D0,
403     &                       WRK(1),NORBT)
404                  DO 350 IQ = IOFFP+NOCCP,IOFFP+NORBP-1
405                     IQ1 = IQ - IOFFP - NOCCP + 1
406                     FVOEX(IDI,IQ,ISIM) = FVOEX(IDI,IQ,ISIM) + WRK(IQ1)
407 350              CONTINUE
408               END IF
409            END IF
410 200     CONTINUE
411 100  CONTINUE
412      IF (((ITYPCD.EQ.1).OR.(ITYPCD.EQ.3)) .AND. (ICI.NE.IDI) ) THEN
413         CALL DGETRN(H2D,NORBT,NORBT)
414         NDITR = -NDITR
415         DO 110 ISIM = 1,NSIM
416            DO 210 IRSYM = 1,NSYM
417               IPSYM = MULD2H(IRSYM,ICDSYM)
418               NORBR = NORB(IRSYM)
419               NSSHP = NSSH(IPSYM)
420               IF (NORBR.EQ.0) GO TO 210
421               IOFFR = IORB(IRSYM) + 1
422               IOFFP = IORB(IPSYM) + 1
423               NORBP = NORB(IPSYM)
424               IF (NORBP.EQ.0) GO TO 210
425               NOCCP = NOCC(IPSYM)
426C
427               IDRSYM = MULD2H(IDSYM,IRSYM)
428               ICRSYM = MULD2H(ICSYM,IRSYM)
429               IF (ICRSYM.EQ.KSYMOP) THEN
430                  IF ((.NOT.TRPLET).AND.(ITYPCD.EQ.1)) THEN
431C
432C  FCOCO(P,Q) = FCOCO(P,Q) + SUM(R) <QJ/PR>*(ZYMAT(J,R)-ZYMAT(R,J))
433C
434                     DO 1101 IR=1,NORBR
435                        WRK(IR) = ZYMAT(ICI,IOFFR-1+IR,ISIM)-
436     *                         ZYMAT(IOFFR-1+IR,ICI,ISIM)
437 1101                CONTINUE
438                     CALL DGEMM('N','N',NORBP,1,NORBR,1.D0,
439     &                          H2D(IOFFP,IOFFR),NORBT,
440     &                          WRK,NORBT,1.D0,
441     &                          FCOCO(IOFFP,IDI,ISIM),NORBT)
442C
443                  END IF
444               END IF
445               IF (IDRSYM.EQ.KSYMOP) THEN
446                  IF ((.NOT.TRPLET).AND.(ITYPCD.EQ.3)) THEN
447C
448C  FVOCO(P,Q) = FVOCO(P,Q) + SUM(R) <QY/PR>*(DENB(R,Y)-DENA(R,Y))
449C                                    YQ/RP        R,Y       R,Y
450C
451                     DO 1201 IR=1,NORBR
452                        WRK(IR) = DENB(IOFFR-1+IR,NDIW,ISIM)-
453     *                            DENA(IOFFR-1+IR,NDIW,ISIM)
454 1201                CONTINUE
455                     CALL DGEMM('T','N',NORBP,1,NORBR,1.D0,
456     &                          H2D(IOFFR,IOFFP),NORBT,
457     &                          WRK,NORBT,1.D0,
458     &                          FVOCO(IOFFP,ICI,ISIM),NORBT)
459C
460                  END IF
461               END IF
462               IF (NSSHP.EQ.0) GO TO 210
463               IF (ICRSYM.EQ.KSYMOP) THEN
464                  IF (ITYPCD.EQ.1) THEN
465C
466C  FCOEX(P,Q) = FCOEX(P,Q) - SUM(R) <QJ/RP>*ZYMAT(J,R) Q OCC P SEC
467C
468
469                     DO 112 IR=1,NORBR
470                        WRK(IR) =-ZYMAT(ICI,IOFFR-1+IR,ISIM)
471 112                 CONTINUE
472                     CALL DGEMM('T','N',NSSHP,1,NORBR,1.D0,
473     &                          H2D(IOFFR,IOFFP+NOCCP),NORBT,
474     &                          WRK(1),NORBR,1.D0,
475     &                          FCOEX(IOFFP+NOCCP,IDI,ISIM),NORBT)
476C
477C  FCOEX(P,Q) = FCOEX(P,Q) + SUM(R) <JP/QR>*ZYMAT(R,J) P OCC Q SEC
478C                                    PJ RQ
479C
480                     CALL DGEMM('T','N',1,NSSHP,NORBR,1.D0,
481     &                          ZYMAT(IOFFR,ICI,ISIM),NORBT,
482     &                          H2D(IOFFR,IOFFP+NOCCP),NORBT,1.D0,
483     &                          FCOEX(IDI,IOFFP+NOCCP,ISIM),NORBT)
484C
485                  END IF
486               END IF
487               IF (IDRSYM.EQ.KSYMOP) THEN
488                  IF (ITYPCD.EQ.3) THEN
489C
490C  FVOEX(P,Q) = FVOEX(P,Q) - SUM(R) <QY/RP>*DENB(R,Y) Q OCC P SEC
491C                                    YQ/PR       R,Y
492                     CALL DGEMM('N','N',NSSHP,1,NORBR,-1.D0,
493     &                          H2D(IOFFP+NOCCP,IOFFR),NORBT,
494     &                          DENB(IOFFR,NDIW,ISIM),NORBT,1.D0,
495     &                          FVOEX(IOFFP+NOCCP,ICI,ISIM),NORBT)
496C
497C  FVOEX(P,Q) = FVOEX(P,Q) + SUM(R) <XP/QR>*DENA(R,X)   (P OCC)
498C
499                     CALL DGEMM('N','N',NSSHP,1,NORBR,1.D0,
500     &                          H2D(IOFFP+NOCCP,IOFFR),NORBT,
501     &                          DENA(IOFFR,NDIW,ISIM),NORBT,0.D0,
502     &                          WRK(1),NORBT)
503                     DO 353 IQ = IOFFP+NOCCP,IOFFP+NORBP-1
504                        IQ1 = IQ - IOFFP - NOCCP + 1
505                        FVOEX(ICI,IQ,ISIM) = FVOEX(ICI,IQ,ISIM)
506     *                                      + WRK(IQ1)
507 353                 CONTINUE
508                  END IF
509               END IF
510C
511 210        CONTINUE
512 110     CONTINUE
513      END IF
514      IF ( NDITR.LT.0) CALL DGETRN(H2D,NORBT,NORBT)
515      RETURN
516      END
517C  /* Deck ftddr */
518      SUBROUTINE FTDDR(NSIM,ICI1,IDI1,FVTD,DVT,H2D)
519C
520C CALCULATE EXCHANGE CONTRIBUTIONS TO ACTIVE FOCK MATRIX WITH A
521C TRANSITION DENSITY MATRIX
522C
523C  FV(P,Q) = FV(P,Q) -0.5*SUM(X,Y)*<XY/QP>*DVT(X,Y)
524C THE SUM X,Y ARE TAKEN BY READING IN ALL DIRAC DISTRIBUTIONS
525C
526#include "implicit.h"
527C
528C Used from common blocks:
529C   INFORB :
530C   INFIND : ISMO,IOBTYP
531#include "maxorb.h"
532#include "maxash.h"
533#include "inforb.h"
534#include "infind.h"
535#include "infrsp.h"
536#include "wrkrsp.h"
537#include "orbtypdef.h"
538C
539      DIMENSION FVTD(NORBT,NORBT,*),DVT(NASHT,NASHT,*)
540      PARAMETER ( DMP5 = -0.5D0 )
541C
542C     Order (C,D) index such that C .ge. D
543C     in inactive-active-secondary order (using ISW)
544C
545      NDITR = 1
546      IF (ISW(ICI1) .GE. ISW(IDI1)) THEN
547         ICI = ICI1
548         IDI = IDI1
549      ELSE
550         ICI = IDI1
551         IDI = ICI1
552         CALL DGETRN(H2D,NORBT,NORBT)
553         NDITR = -NDITR
554      END IF
555C
556C     Find distribution type ITYPCD =
557C     1:inactive-inactive  2:active-inactive  3:active-active
558C     4:secondary-inactive 5:secondary-active 6:secondary-secondary
559C     We only need occupied-occupied distributions, i.e. itypcd .le. 3
560C
561      ITYPC  = IOBTYP(ICI)
562      ITYPD  = IOBTYP(IDI)
563      ITYPCD = IDBTYP(ITYPC,ITYPD)
564      IF (ITYPCD .NE. 3) RETURN
565      IF (ITYPC .EQ. JTACT) NCIW = ICH(ICI)
566      IF (ITYPD .EQ. JTACT) NDIW = ICH(IDI)
567C
568      ICSYM = ISMO(ICI)
569      IDSYM = ISMO(IDI)
570      ICDSYM = MULD2H(ICSYM,IDSYM)
571      IF ( ICDSYM.EQ.KSYMOP) THEN
572         DO 100 ISIM = 1,NSIM
573            FAC = DMP5*DVT(NDIW,NCIW,ISIM)
574            CALL DAXPY(NORBT*NORBT,FAC,H2D,1,FVTD(1,1,ISIM),1)
575 100     CONTINUE
576         IF (NCIW.NE.NDIW) THEN
577            CALL DGETRN(H2D,NORBT,NORBT)
578            NDITR = -NDITR
579            DO 200 ISIM = 1,NSIM
580               FAC = DMP5*DVT(NCIW,NDIW,ISIM)
581               CALL DAXPY(NORBT*NORBT,FAC,H2D,1,FVTD(1,1,ISIM),1)
582 200        CONTINUE
583         END IF
584      END IF
585      IF ( NDITR.LT.0) CALL DGETRN(H2D,NORBT,NORBT)
586      RETURN
587      END
588C  /* Deck ftdmu */
589      SUBROUTINE FTDMU(NSIM,ICI1,IDI1,FVTD,DVT,H2)
590C
591C CALCULATE EXCHANGE CONTRIBUTIONS TO ACTIVE FOCK MATRIX WITH A
592C TRANSITION DENSITY MATRIX
593C
594C  FV(P,Q) = FV(P,Q) + SUM(X,Y)*(PQ/XY)*DVT(X,Y)
595C THE SUM X,Y ARE TAKEN BY READING IN ALL MULLIKEN DISTRIBUTIONS
596C
597#include "implicit.h"
598C
599C Used from common blocks:
600C   INFORB :
601C   INFIND : ISMO,IOBTYP
602#include "maxorb.h"
603#include "maxash.h"
604#include "inforb.h"
605#include "infind.h"
606#include "infrsp.h"
607#include "wrkrsp.h"
608#include "orbtypdef.h"
609C
610      DIMENSION FVTD(NORBT,NORBT,*),DVT(NASHT,NASHT,*)
611      DIMENSION H2(NORBT,*)
612C
613      IF (TRPLET) RETURN
614C
615C     Order (C,D) index such that C .ge. D
616C     in inactive-active-secondary order (using ISW)
617C
618      NDITR = 1
619      IF (ISW(ICI1) .GE. ISW(IDI1)) THEN
620         ICI = ICI1
621         IDI = IDI1
622      ELSE
623         ICI = IDI1
624         IDI = ICI1
625      END IF
626C
627C     Find distribution type ITYPCD =
628C     1:inactive-inactive  2:active-inactive  3:active-active
629C     4:secondary-inactive 5:secondary-active 6:secondary-secondary
630C     We only need occupied-occupied distributions, i.e. itypcd .le. 3
631C
632      ITYPC  = IOBTYP(ICI)
633      ITYPD  = IOBTYP(IDI)
634      ITYPCD = IDBTYP(ITYPC,ITYPD)
635      IF (ITYPCD .NE. 3) RETURN
636      IF (ITYPC .EQ. JTACT) NCIW = ICH(ICI)
637      IF (ITYPD .EQ. JTACT) NDIW = ICH(IDI)
638C
639      ICSYM = ISMO(ICI)
640      IDSYM = ISMO(IDI)
641      ICDSYM = MULD2H(ICSYM,IDSYM)
642      IF ( ICDSYM.EQ.KSYMOP) THEN
643         DO 100 ISIM = 1,NSIM
644            IF (NCIW.NE.NDIW) THEN
645               FAC = DVT(NCIW,NDIW,ISIM)+DVT(NDIW,NCIW,ISIM)
646            ELSE
647               FAC = DVT(NCIW,NDIW,ISIM)
648            ENDIF
649            CALL DAXPY(NORBT*NORBT,FAC,H2,1,FVTD(1,1,ISIM),1)
650 100     CONTINUE
651      END IF
652      RETURN
653      END
654