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      SUBROUTINE C3FOCK(IBEQC,FDTI,VECB,VECC,
20     *                 ZYMB,ZYMC,KZYVB,KZYVC,
21     *                 ISYMA,ISYMB,ISYMC,
22     *                 CMO,FC,MJWOP,WRK,LWRK)
23C
24C PURPOSE:
25C To calculate Fock matrix (FDTI) with doubly one-index transformed
26C integrals to be used in direct cubic RPA.
27C March 1997: tsaue Just one call to SIRFCK, big speedup !
28C
29#include "implicit.h"
30C
31      DIMENSION FDTI(N2ORBX)
32      DIMENSION VECB(*),VECC(*)
33      DIMENSION ZYMB(NORBT,NORBT),ZYMC(NORBT,NORBT)
34      DIMENSION CMO(*),FC(*),WRK(*),ISYMDM(3)
35C
36C  INFDIM : NASHDI
37C  INFINP : DIRFCK
38C  WRKRSP : KSYMOP
39C
40#include "maxorb.h"
41#include "maxash.h"
42#include "priunit.h"
43#include "inforb.h"
44#include "infinp.h"
45#include "infvar.h"
46      DIMENSION MJWOP(2,MAXWOP,8)
47#include "wrkrsp.h"
48#include "infrsp.h"
49C
50      CALL QENTER('C3FOCK')
51      IF (8*N2BASX .GT. LWRK) THEN
52         WRITE (LUPRI,'(A)') ' Due to lack of available memory, I '//
53     &     'use a slower C3FOCK routine'
54         CALL C3FCKM(IBEQC,FDTI,VECB,VECC,
55     *               ZYMB,ZYMC,KZYVB,KZYVC,
56     *               ISYMA,ISYMB,ISYMC,
57     *               CMO,FC,MJWOP,WRK,LWRK)
58      ELSE
59         CALL C3FCKO(IBEQC,FDTI,VECB,VECC,
60     *               ZYMB,ZYMC,KZYVB,KZYVC,
61     *               ISYMA,ISYMB,ISYMC,
62     *               CMO,FC,MJWOP,WRK(1),LWRK)
63      END IF
64      CALL QEXIT('C3FOCK')
65      RETURN
66      END
67C
68      SUBROUTINE C3FCKO(IBEQC,FDTI,VECB,VECC,
69     *                 ZYMB,ZYMC,KZYVB,KZYVC,
70     *                 ISYMA,ISYMB,ISYMC,
71     *                 CMO,FC,MJWOP,WRK,LWRK)
72C
73C PURPOSE:
74C To calculate Fock matrix (FDTI) with doubly one-index transformed
75C integrals to be used in direct cubic RPA.
76C March 1997: tsaue Just one call to SIRFCK, big speedup !
77C
78#include "implicit.h"
79#include "dummy.h"
80#include "thrzer.h"
81C
82      DIMENSION FDTI(N2ORBX)
83      DIMENSION VECB(*),VECC(*)
84      DIMENSION ZYMB(NORBT,NORBT),ZYMC(NORBT,NORBT)
85      DIMENSION CMO(*),FC(*),WRK(*),ISYMDM(3),IFCTYP(3)
86C
87C  INFDIM : NASHDI
88C  INFINP : DIRFCK
89C  WRKRSP : KSYMOP
90C
91#include "maxorb.h"
92#include "maxash.h"
93#include "priunit.h"
94#include "inforb.h"
95#include "infinp.h"
96#include "infvar.h"
97      DIMENSION MJWOP(2,MAXWOP,8)
98#include "wrkrsp.h"
99#include "infrsp.h"
100#include "infpri.h"
101C
102      PARAMETER ( D1 = 1.0D0, D2 = 2.0D0, D4 = 4.0D0 )
103C
104C Allocate work space for Fock matrices
105C Corresponding Fock matrices in equations
106C
107C   F-1 <--> F12 + F21
108C   F-2 <--> F1
109C   F-3 <--> F2
110C
111      CALL QENTER('C3FCKO')
112      KF1AO = 1
113      KF2AO = KF1AO + N2BASX
114      KF3AO = KF2AO + N2BASX
115      KD1AO = KF3AO + N2BASX
116      KD2AO = KD1AO + N2BASX
117      KD3AO = KD2AO + N2BASX
118      KFREE = KD3AO + N2BASX
119      LFREE = LWRK  - KFREE
120      IF (LFREE.LT.0) CALL ERRWRK('C3FCKO',LFREE,LWRK)
121C
122C     Unpack the response vectors
123C
124      CALL GTZYMT(1,VECB,KZYVB,ISYMB,ZYMB,MJWOP)
125      CALL GTZYMT(1,VECC,KZYVC,ISYMC,ZYMC,MJWOP)
126C
127C
128C Construct density matrices in AO basis
129C
130C
131      CALL DZERO(WRK(KF1AO),6*N2BASX)
132C
133C Alternative construction need to assign KTMP, KBCZYM, ISYMBC
134C Further scale D1AO with factor 2
135C
136      ISYMDM(1) = MULD2H(ISYMB,ISYMC)
137      ISYMDM(2) = ISYMB
138      ISYMDM(3) = ISYMC
139C
140      KBCZYM = KF1AO
141      KTMP   = KF2AO
142C
143      CALL ZYTRA1(ISYMB,ZYMB,CMO,WRK(KTMP),WRK(KD2AO))
144      CALL ZYTRA1(ISYMC,ZYMC,CMO,WRK(KTMP),WRK(KD3AO))
145C
146      CALL DZERO(WRK(KTMP),N2BASX)
147      CALL ZYMUL2(ISYMB,ISYMC,ZYMB,ZYMC,WRK(KBCZYM))
148      CALL ZYMUL2(ISYMC,ISYMB,ZYMC,ZYMB,WRK(KBCZYM))
149      CALL ZYTRA2(ISYMDM(1),WRK(KBCZYM),CMO,WRK(KTMP),WRK(KD1AO))
150C
151C
152C Transpose and scale by 2 to agree with convention used by SIRFCK
153C Further scale by 2 due to permutations. D2AO and D3AO scaled in one call.
154C
155      CALL DSCAL(3*N2BASX,D4,WRK(KD1AO),1)
156C
157      IF (IPRRSP.GT.200) THEN
158         WRITE(LUPRI,'(//A)')'D12 in C3FOCK'
159         CALL OUTPUT(WRK(KD1AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
160         WRITE(LUPRI,'(//A)')'D1 in C3FOCK'
161         CALL OUTPUT(WRK(KD2AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
162         WRITE(LUPRI,'(//A)')'D2 in C3FOCK'
163         CALL OUTPUT(WRK(KD3AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
164      END IF
165C
166C Construct Fock matrices in AO basis
167C
168      NFMATX = 1
169C
170C     IFCTYP = XY
171C       X indicates symmetry about diagonal
172C         X = 0 No symmetry
173C         X = 1 Symmetric
174C         X = 2 Anti-symmetric
175C       Y indicates contributions
176C         Y = 0 no contribution !
177C         Y = 1 Coulomb
178C         Y = 2 Exchange
179C         Y = 3 Coulomb + Exchange
180C
181C     Check if density matrix is unsymmetric (IX=0),
182C     symmetric (IX=10), antisymmetric (IX=20), or zero matrix (IX=30)
183C     to threshold THRZER
184C
185      JD1AO = KD1AO
186      DO I = 1,3
187         IX = 10 * MATSYM(NBAST,NBAST,WRK(JD1AO),THRZER)
188C        INTEGER FUNCTION MATSYM(N,NDIM,AMAT,THRZER)
189C
190         IF (IX .EQ. 30) THEN
191C           zero density matrix, do nothing !
192            IFCTYP(I) = 0
193CCCC     ELSE IF ("triplet Fock matrix") THEN
194C           only exchange ! hjaaj 990505: triplet not implemented!
195CCCC        IFCTYP(I) = IX + 2
196         ELSE IF (IX .EQ. 20) THEN
197C           Only exchange if antisymmetric density matrix
198            IFCTYP(I) = IX + 2
199         ELSE
200C           Coulomb+exchange
201            IFCTYP(I) = IX + 3
202         END IF
203         JD1AO = JD1AO + N2BASX
204      END DO
205C
206      CALL DZERO(WRK(KF1AO),3*N2BASX)
207      IF (IBEQC.EQ.1) THEN
208        NFMATX = 2
209        CALL SIRFCK(WRK(KF1AO),WRK(KD1AO),NFMATX,ISYMDM,IFCTYP,
210     *            DIRFCK,WRK(KFREE),LFREE)
211        CALL DCOPY(N2BASX,WRK(KF2AO),1,WRK(KF3AO),1)
212      ELSE
213        NFMATX = 3
214        CALL SIRFCK(WRK(KF1AO),WRK(KD1AO),NFMATX,ISYMDM,IFCTYP,
215     *            DIRFCK,WRK(KFREE),LFREE)
216      ENDIF
217C
218      IF (IPRRSP.GT.200) THEN
219         WRITE(LUPRI,'(//A)')'F12 in AO basis in C3FOCK'
220         CALL OUTPUT(WRK(KF1AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
221         WRITE(LUPRI,'(//A)')'F1 in AO basis in C3FOCK'
222         CALL OUTPUT(WRK(KF2AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
223         WRITE(LUPRI,'(//A)')'F2 in AO basis in C3FOCK'
224         CALL OUTPUT(WRK(KF3AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
225      END IF
226C
227C Transform Fock matrices to MO basis
228C Reuse memory used for density matrices
229C
230      KF2MO = KD1AO
231      KF3MO = KF2MO + N2ORBX
232C
233      CALL FAOMO(ISYMDM(1),WRK(KF1AO),FDTI,CMO,WRK(KFREE),LFREE)
234      CALL FAOMO(ISYMDM(2),WRK(KF2AO),WRK(KF2MO),CMO,WRK(KFREE),LFREE)
235      CALL FAOMO(ISYMDM(3),WRK(KF3AO),WRK(KF3MO),CMO,WRK(KFREE),LFREE)
236C
237      IF (IPRRSP.GT.200) THEN
238         WRITE(LUPRI,'(//A)')'F12 in MO basis in C3FOCK'
239         CALL OUTPUT(FDTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
240         WRITE(LUPRI,'(//A)')'F1 in MO basis in C3FOCK'
241         CALL OUTPUT(WRK(KF2MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
242         WRITE(LUPRI,'(//A)')'F2 in MO basis in C3FOCK'
243         CALL OUTPUT(WRK(KF3MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
244      END IF
245C
246C Evaluate FDTI see formula
247C
248C 2*F2 + [k2,Fc] --> F-3
249C 2*F1 + [k1,Fc] --> F-2
250C
251      KSYMOP = ISYMC
252      CALL FCKOIN(1,FC,DUMMY,ZYMC,WRK(KF3MO),DUMMY)
253      KSYMOP = ISYMB
254      CALL FCKOIN(1,FC,DUMMY,ZYMB,WRK(KF2MO),DUMMY)
255C
256C F12  + [k2,F-3] --> FDTI
257C FDTI + [k1,F-2] --> FDTI
258C
259      CALL OITH1(ISYMC,ZYMC,WRK(KF2MO),FDTI,ISYMDM(2))
260      CALL OITH1(ISYMB,ZYMB,WRK(KF3MO),FDTI,ISYMDM(3))
261C
262C Scale with 1/2 ( from definition of E3 )
263C
264      CALL DSCAL(NORBT*NORBT,1/D2,FDTI,1)
265C
266      IF (IPRRSP.GT.100) THEN
267         WRITE(LUPRI,'(//A)')'Final result in C3FCKO'
268         CALL OUTPUT(FDTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
269      END IF
270C
271      CALL QEXIT('C3FCKO')
272      RETURN
273      END
274C
275      SUBROUTINE C3FCKM(IBEQC,FDTI,VECB,VECC,
276     *                 ZYMB,ZYMC,KZYVB,KZYVC,
277     *                 ISYMA,ISYMB,ISYMC,
278     *                 CMO,FC,MJWOP,WRK,LWRK)
279C
280C PURPOSE:
281C To calculate Fock matrix (FDTI) with doubly one-index transformed
282C integrals to be used in direct cubic RPA.
283C
284#include "implicit.h"
285#include "dummy.h"
286#include "thrzer.h"
287C
288      DIMENSION FDTI(N2ORBX)
289      DIMENSION VECB(*),VECC(*)
290      DIMENSION ZYMB(NORBT,NORBT),ZYMC(NORBT,NORBT)
291      DIMENSION CMO(*),FC(*),WRK(*)
292C
293C  INFDIM : NASHDI
294C  INFINP : DIRFCK
295C  WRKRSP : KSYMOP
296C
297#include "maxorb.h"
298#include "maxash.h"
299#include "priunit.h"
300#include "inforb.h"
301#include "infinp.h"
302#include "infvar.h"
303      DIMENSION MJWOP(2,MAXWOP,8)
304#include "wrkrsp.h"
305#include "infrsp.h"
306#include "infpri.h"
307C
308      PARAMETER ( D1 = 1.0D0, D2 = 2.0D0, D4 = 4.0D0 )
309C
310      NFMATX = 1
311C
312C Allocate work space for Fock matrices
313C Corresponding Fock matrices in equations
314C
315C   F-1 <--> F12 + F21
316C   F-2 <--> F1
317C   F-3 <--> F2
318C
319      CALL QENTER('C3FCKM')
320      KFAO  = 1
321      KDAO  = KFAO  + N2BASX
322      KFRE1 = KDAO  + N2BASX
323      KFREE = KFRE1 + N2BASX
324      LFREE = LWRK  - KFREE
325      IF (LFREE.LT.0) CALL ERRWRK('C3FCKM',LFREE,LWRK)
326C
327      KFREE = KFRE1
328      LFREE = LWRK - KFREE
329C
330C     Unpack the response vectors
331C
332      CALL GTZYMT(1,VECB,KZYVB,ISYMB,ZYMB,MJWOP)
333      CALL GTZYMT(1,VECC,KZYVC,ISYMC,ZYMC,MJWOP)
334C
335C
336C Construct density matrices in AO basis
337C
338C
339      CALL DZERO(WRK(KFAO),3*N2BASX)
340      CALL DZERO(FDTI,N2ORBX)
341      CALL CDENS2(ISYMB,ISYMC,CMO,ZYMB,ZYMC,
342     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),FDTI)
343      CALL CDENS2(ISYMC,ISYMB,CMO,ZYMC,ZYMB,
344     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),FDTI)
345      CALL DSCAL(N2BASX,D2,WRK(KDAO),1)
346      IF (IPRRSP.GT.200) THEN
347         WRITE(LUPRI,'(//A)')'D12 in C3FOCK'
348         CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
349      END IF
350C
351      ISYMDM = MULD2H(ISYMB,ISYMC)
352      CALL DZERO(WRK(KFAO),N2BASX)
353Chj
354C     Check if density matrix is unsymmetric (IX=0),
355C     symmetric (IX=10), antisymmetric (IX=20), or zero matrix (IX=30)
356C     to threshold THRZER
357C
358      IX = 10 * MATSYM(NBAST,NBAST,WRK(KDAO),THRZER)
359C
360      IF (IX .EQ. 30) THEN
361C        zero density matrix, do nothing !
362         IFCTYP = 0
363      ELSE IF (IX .EQ. 20) THEN
364C        Only exchange if antisymmetric density matrix
365         IFCTYP = IX + 2
366      ELSE
367C        Coulomb+exchange
368         IFCTYP = IX + 3
369      END IF
370      IF (IFCTYP .NE. 0) THEN
371         CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
372     *               DIRFCK,WRK(KFREE),LFREE)
373      END IF
374      IF (IPRRSP.GT.200) THEN
375         WRITE(LUPRI,'(//A)')'F12 in AO basis in C3FOCK'
376         CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
377      END IF
378C
379      CALL FAOMO(ISYMDM,WRK(KFAO),FDTI,CMO,WRK(KFREE),LFREE)
380      IF (IPRRSP.GT.200) THEN
381         WRITE(LUPRI,'(//A)')'F12 in MO basis in C3FOCK'
382         CALL OUTPUT(FDTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
383      END IF
384C
385C     Second contribution
386C
387      CALL DZERO(WRK(KFAO),2*N2BASX)
388      CALL CDENS1(ISYMB,CMO,ZYMB,WRK(KDAO),WRK(KFREE),LFREE)
389      CALL DSCAL(N2BASX,D4,WRK(KDAO),1)
390      IF (IPRRSP .GT. 200) THEN
391         WRITE(LUPRI,'(//A)')'D1 in C3FOCK'
392         CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
393      END IF
394C
395      ISYMDM = ISYMB
396Chj
397C     Check if density matrix is unsymmetric (IX=0),
398C     symmetric (IX=10), antisymmetric (IX=20), or zero matrix (IX=30)
399C     to threshold THRZER
400C
401      IX = 10 * MATSYM(NBAST,NBAST,WRK(KDAO),THRZER)
402C
403      IF (IX .EQ. 30) THEN
404C        zero density matrix, do nothing !
405         IFCTYP = 0
406      ELSE IF (IX .EQ. 20) THEN
407C        Only exchange if antisymmetric density matrix
408         IFCTYP = IX + 2
409      ELSE
410C        Coulomb+exchange
411         IFCTYP = IX + 3
412      END IF
413      IF (IFCTYP .NE. 0) THEN
414         CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
415     *               DIRFCK,WRK(KFREE),LFREE)
416      END IF
417      IF (IPRRSP .GT. 200) THEN
418         WRITE(LUPRI,'(//A)')'F1 in AO basis in C3FOCK'
419         CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
420      END IF
421C
422      CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
423      IF (IPRRSP .GT. 200) THEN
424         WRITE(LUPRI,'(//A)')'F1 in MO basis in C3FOCK'
425         CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
426      END IF
427C
428      KSYMOP = ISYMB
429      CALL FCKOIN(1,FC,DUMMY,ZYMB,WRK(KDAO),DUMMY)
430      IF (IBEQC .EQ. 1) THEN
431C
432C     If B and C vector identical, all that remains to do is to scale
433C     this term by a factor of two.
434C
435         CALL DSCAL(N2ORBX,D2,WRK(KDAO),1)
436      ELSE
437         CALL OITH1(ISYMC,ZYMC,WRK(KDAO),FDTI,ISYMDM)
438C
439C     Third contribution
440C
441         CALL DZERO(WRK(KFAO),2*N2BASX)
442         CALL CDENS1(ISYMC,CMO,ZYMC,WRK(KDAO),WRK(KFREE),LFREE)
443         CALL DSCAL(N2BASX,D4,WRK(KDAO),1)
444         IF (IPRRSP .GT. 200) THEN
445            WRITE(LUPRI,'(//A)')'D2 in C3FOCK'
446            CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
447         END IF
448C
449         ISYMDM = ISYMC
450Chj
451C     Check if density matrix is unsymmetric (IX=0),
452C     symmetric (IX=10), antisymmetric (IX=20), or zero matrix (IX=30)
453C     to threshold THRZER
454C
455         IX = 10 * MATSYM(NBAST,NBAST,WRK(KDAO),THRZER)
456C
457         IF (IX .EQ. 30) THEN
458C           zero density matrix, do nothing !
459            IFCTYP = 0
460         ELSE IF (IX .EQ. 20) THEN
461C           Only exchange if antisymmetric density matrix
462            IFCTYP = IX + 2
463         ELSE
464C           Coulomb+exchange
465            IFCTYP = IX + 3
466         END IF
467         IF (IFCTYP .NE. 0) THEN
468            CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
469     *                  DIRFCK,WRK(KFREE),LFREE)
470         END IF
471         IF (IPRRSP .GT. 200) THEN
472            WRITE(LUPRI,'(//A)')'F2 in AO basis in C3FOCK'
473            CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
474         END IF
475C
476         CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
477         IF (IPRRSP .GT. 200) THEN
478            WRITE(LUPRI,'(//A)')'F2 in MO basis in C3FOCK'
479            CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
480         END IF
481C
482         KSYMOP = ISYMC
483         CALL FCKOIN(1,FC,DUMMY,ZYMC,WRK(KDAO),DUMMY)
484         CALL OITH1(ISYMB,ZYMB,WRK(KDAO),FDTI,ISYMDM)
485      END IF
486C
487C     All contributions to FDTI now finished, we
488C     scale with 1/2 ( from definition of E3 )
489C
490      CALL DSCAL(NORBT*NORBT,1/D2,FDTI,1)
491C
492      IF (IPRRSP.GT.100) THEN
493         WRITE(LUPRI,'(//A)')'Final result in C3FCKM'
494         CALL OUTPUT(FDTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
495      END IF
496C
497      CALL QEXIT('C3FCKM')
498      RETURN
499      END
500C
501      SUBROUTINE C4FOCK(IBCDEQ,FTTI,VECB,VECC,VECD,
502     *                 ZYMB,ZYMC,ZYMD,KZYVB,KZYVC,KZYVD,
503     *                 ISYMA,ISYMB,ISYMC,ISYMD,
504     *                 CMO,FC,MJWOP,WRK,LWRK)
505C
506C PURPOSE:
507C To calculate Fock matrix (FTTI) with triply one-index transformed
508C integrals to be used in direct cubic RPA.
509C March 1997: tsaue Just one call to SIRFCK, big speedup !
510C July 1997:  kruud. Interface routine; further path depends on available
511C                    memory
512C
513#include "implicit.h"
514C
515      DIMENSION FTTI(N2ORBX)
516      DIMENSION VECB(*),VECC(*),VECD(*)
517      DIMENSION ZYMB(NORBT,NORBT),ZYMC(NORBT,NORBT),ZYMD(NORBT,NORBT)
518      DIMENSION CMO(*),FC(*),WRK(*)
519C
520C  INFDIM : NASHDI
521C  INFINP : DIRFCK
522C  WRKRSP : KSYMOP
523C
524#include "maxorb.h"
525#include "maxash.h"
526#include "priunit.h"
527#include "inforb.h"
528#include "infinp.h"
529#include "infvar.h"
530      DIMENSION MJWOP(2,MAXWOP,8)
531#include "wrkrsp.h"
532#include "infrsp.h"
533C
534      IF (14*N2BASX .GT. LWRK) THEN
535         WRITE (LUPRI,'(A)') ' Due to potential lack of available '//
536     &     'memory, I use a slower C4FOCK routine'
537         CALL C4FCKM(IBCDEQ,FTTI,VECB,VECC,VECD,
538     *               ZYMB,ZYMC,ZYMD,KZYVB,KZYVC,KZYVD,
539     *               ISYMA,ISYMB,ISYMC,ISYMD,
540     *               CMO,FC,MJWOP,WRK,LWRK)
541      ELSE
542         CALL C4FCKO(IBCDEQ,FTTI,VECB,VECC,VECD,
543     *               ZYMB,ZYMC,ZYMD,KZYVB,KZYVC,KZYVD,
544     *               ISYMA,ISYMB,ISYMC,ISYMD,
545     *               CMO,FC,MJWOP,WRK,LWRK)
546      END IF
547C
548      RETURN
549      END
550C
551      SUBROUTINE C4FCKO(IBCDEQ,FTTI,VECB,VECC,VECD,
552     *                  ZYMB,ZYMC,ZYMD,KZYVB,KZYVC,KZYVD,
553     *                  ISYMA,ISYMB,ISYMC,ISYMD,
554     *                  CMO,FC,MJWOP,WRK,LWRK)
555#include "implicit.h"
556#include "dummy.h"
557C
558      DIMENSION FTTI(N2ORBX)
559      DIMENSION VECB(*),VECC(*),VECD(*)
560      DIMENSION ZYMB(NORBT,NORBT),ZYMC(NORBT,NORBT),ZYMD(NORBT,NORBT)
561      DIMENSION CMO(*),FC(*),WRK(*),ISYMDM(7),IFCTYP(7)
562C
563C  INFDIM : NASHDI
564C  INFINP : DIRFCK
565C  WRKRSP : KSYMOP
566C
567#include "maxorb.h"
568#include "maxash.h"
569#include "priunit.h"
570#include "inforb.h"
571#include "infinp.h"
572#include "infvar.h"
573      DIMENSION MJWOP(2,MAXWOP,8)
574#include "wrkrsp.h"
575#include "infrsp.h"
576#include "infpri.h"
577C
578      PARAMETER ( D1 = 1.0D0, D2 = 2.0D0, D6 = 6.0D0 , D8 = 8.0D0 ,
579     &           D12 = 12.0D0)
580C Allocate work space for Fock matrices
581C Corresponding Fock matrices in equations
582C
583C   F-1 <--> F123 + permutations
584C   F-2 <--> F12 + F21
585C   F-3 <--> F13 + F31
586C   F-4 <--> F23 + F32
587C   F-5 <--> F1
588C   F-6 <--> F2
589C   F-7 <--> F3
590C
591      KF1AO = 1
592      KF2AO = KF1AO + N2BASX
593      KF3AO = KF2AO + N2BASX
594      KF4AO = KF3AO + N2BASX
595      KF5AO = KF4AO + N2BASX
596      KF6AO = KF5AO + N2BASX
597      KF7AO = KF6AO + N2BASX
598      KD1AO = KF7AO + N2BASX
599      KD2AO = KD1AO + N2BASX
600      KD3AO = KD2AO + N2BASX
601      KD4AO = KD3AO + N2BASX
602      KD5AO = KD4AO + N2BASX
603      KD6AO = KD5AO + N2BASX
604      KD7AO = KD6AO + N2BASX
605      KFREE = KD7AO + N2BASX
606      LFREE = LWRK  - KFREE
607C
608      IF (IPRRSP.GT.7) WRITE(LUPRI,'(//A,2(/A,I10))')
609     *' Memory allocation in C4FCKO',
610     *' Available workspace: ', LWRK,
611     *' Allocated          : ', KFREE
612C
613C     Unpack the response vectors
614C
615      CALL GTZYMT(1,VECB,KZYVB,ISYMB,ZYMB,MJWOP)
616      CALL GTZYMT(1,VECC,KZYVC,ISYMC,ZYMC,MJWOP)
617      CALL GTZYMT(1,VECD,KZYVD,ISYMD,ZYMD,MJWOP)
618C
619      CALL DZERO(WRK(KF1AO),14*N2BASX)
620C
621C Alternativ construction need to assign KTMP, K..ZYM, K...ZYM
622C
623      KBCZYM  = KF1AO
624      KBCDZYM = KF2AO
625      KTMP    = KF3AO
626      KBDZYM  = KBCZYM
627      KCDZYM  = KBCZYM
628C
629      ISYMDM(1) = ISYMA
630      ISYMDM(2) = MULD2H(ISYMB,ISYMC)
631      ISYMDM(3) = MULD2H(ISYMB,ISYMD)
632      ISYMDM(4) = MULD2H(ISYMC,ISYMD)
633      ISYMDM(5) = ISYMB
634      ISYMDM(6) = ISYMC
635      ISYMDM(7) = ISYMD
636      ISYMBCD   = MULD2H(ISYMDM(2),ISYMD)
637C
638      CALL ZYTRA1(ISYMB,ZYMB,CMO,WRK(KTMP),WRK(KD5AO))
639      CALL ZYTRA1(ISYMC,ZYMC,CMO,WRK(KTMP),WRK(KD6AO))
640      CALL ZYTRA1(ISYMD,ZYMD,CMO,WRK(KTMP),WRK(KD7AO))
641C
642      CALL DZERO(WRK(KBCZYM),3*N2ORBX)
643C
644      CALL ZYMUL2(ISYMB,ISYMD,ZYMB,ZYMD,WRK(KBDZYM))
645      CALL ZYMUL2(ISYMD,ISYMB,ZYMD,ZYMB,WRK(KBDZYM))
646      CALL ZYMUL3(ISYMDM(3),ISYMC,WRK(KBDZYM),ZYMC,WRK(KBCDZYM))
647      CALL ZYTRA2(ISYMDM(3),WRK(KBDZYM),CMO,WRK(KTMP),WRK(KD3AO))
648C
649      CALL DZERO(WRK(KBDZYM),N2ORBX)
650      CALL DZERO(WRK(KTMP),N2BASX)
651      CALL ZYMUL2(ISYMB,ISYMC,ZYMB,ZYMC,WRK(KBCZYM))
652      CALL ZYMUL2(ISYMC,ISYMB,ZYMC,ZYMB,WRK(KBCZYM))
653      CALL ZYMUL3(ISYMDM(2),ISYMD,WRK(KBCZYM),ZYMD,WRK(KBCDZYM))
654      CALL ZYTRA2(ISYMDM(2),WRK(KBCZYM),CMO,WRK(KTMP),WRK(KD2AO))
655C
656      CALL DZERO(WRK(KCDZYM),N2ORBX)
657      CALL DZERO(WRK(KCDZYM),N2ORBX)
658      CALL ZYMUL2(ISYMC,ISYMD,ZYMC,ZYMD,WRK(KCDZYM))
659      CALL ZYMUL2(ISYMD,ISYMC,ZYMD,ZYMC,WRK(KCDZYM))
660      CALL ZYMUL3(ISYMDM(4),ISYMB,WRK(KCDZYM),ZYMB,WRK(KBCDZYM))
661      CALL ZYTRA2(ISYMDM(4),WRK(KCDZYM),CMO,WRK(KTMP),WRK(KD4AO))
662C
663      CALL DZERO(WRK(KTMP),N2BASX)
664      CALL ZYTRA1(ISYMBCD,WRK(KBCDZYM),CMO,WRK(KTMP),WRK(KD1AO))
665C
666C Transpose and scale by 2 to agree with convention used by SIRFCK
667C Further scale by 3 due to permutations. D2AO -> D7AO scaled in one call
668C
669      CALL DSCAL(N2BASX,D8,WRK(KD1AO),1)
670      CALL DSCAL(3*N2BASX,D12,WRK(KD2AO),1)
671      CALL DSCAL(3*N2BASX,D6,WRK(KD5AO),1)
672C
673      IF (IPRRSP.GT.200) THEN
674         WRITE(LUPRI,'(//A)')'D123 in C4FOCK'
675         CALL OUTPUT(WRK(KD1AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
676         WRITE(LUPRI,'(//A)')'D12 in C4FOCK'
677         CALL OUTPUT(WRK(KD2AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
678         WRITE(LUPRI,'(//A)')'D13 in C4FOCK'
679         CALL OUTPUT(WRK(KD3AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
680         WRITE(LUPRI,'(//A)')'D23 in C4FOCK'
681         CALL OUTPUT(WRK(KD4AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
682         WRITE(LUPRI,'(//A)')'D1 in C4FOCK'
683         CALL OUTPUT(WRK(KD5AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
684         WRITE(LUPRI,'(//A)')'D2 in C4FOCK'
685         CALL OUTPUT(WRK(KD6AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
686         WRITE(LUPRI,'(//A)')'D3 in C4FOCK'
687         CALL OUTPUT(WRK(KD7AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
688      END IF
689C
690C Construct Fock matrices in AO basis
691C
692      DO I = 1,7
693        IFCTYP(I) = 3
694      ENDDO
695      NFMATX = 7
696      IF (IBCDEQ.EQ.24) THEN
697         NFMATX    = 5
698         IFCTYP(3) = 0
699         IFCTYP(4) = 0
700         IFCTYP(6) = 0
701         IFCTYP(7) = 0
702      ELSE IF (IBCDEQ.EQ.2) THEN
703         IFCTYP(4) = 0
704         IFCTYP(6) = 0
705      ELSE IF (IBCDEQ.EQ.3) THEN
706         NFMATX    = 6
707         IFCTYP(4) = 0
708         IFCTYP(7) = 0
709      ELSE IF (IBCDEQ.EQ.4) THEN
710         NFMATX    = 6
711         IFCTYP(3) = 0
712         IFCTYP(7) = 0
713      ENDIF
714C
715      CALL DZERO(WRK(KF1AO),7*N2BASX)
716      CALL SIRFCK(WRK(KF1AO),WRK(KD1AO),NFMATX,ISYMDM,IFCTYP,
717     *            DIRFCK,WRK(KFREE),LFREE)
718C
719      IF (IBCDEQ.EQ.24) THEN
720         CALL DCOPY(N2BASX,WRK(KF2AO),1,WRK(KF3AO),1)
721         CALL DCOPY(N2BASX,WRK(KF2AO),1,WRK(KF4AO),1)
722         CALL DCOPY(N2BASX,WRK(KF5AO),1,WRK(KF6AO),1)
723         CALL DCOPY(N2BASX,WRK(KF5AO),1,WRK(KF7AO),1)
724      ELSE IF (IBCDEQ.EQ.2) THEN
725         CALL DCOPY(N2BASX,WRK(KF3AO),1,WRK(KF4AO),1)
726         CALL DCOPY(N2BASX,WRK(KF5AO),1,WRK(KF6AO),1)
727      ELSE IF (IBCDEQ.EQ.3) THEN
728         CALL DCOPY(N2BASX,WRK(KF2AO),1,WRK(KF4AO),1)
729         CALL DCOPY(N2BASX,WRK(KF5AO),1,WRK(KF7AO),1)
730      ELSE IF (IBCDEQ.EQ.4) THEN
731         CALL DCOPY(N2BASX,WRK(KF2AO),1,WRK(KF3AO),1)
732         CALL DCOPY(N2BASX,WRK(KF6AO),1,WRK(KF7AO),1)
733      END IF
734C
735      IF (IPRRSP.GT.200) THEN
736         WRITE(LUPRI,'(//A)')'F123 in AO basis in C4FOCK'
737         CALL OUTPUT(WRK(KF1AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
738         WRITE(LUPRI,'(//A)')'F12 in AO basis in C4FOCK'
739         CALL OUTPUT(WRK(KF2AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
740         WRITE(LUPRI,'(//A)')'F13 in AO basis in C4FOCK'
741         CALL OUTPUT(WRK(KF3AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
742         WRITE(LUPRI,'(//A)')'F23 in AO basis in C4FOCK'
743         CALL OUTPUT(WRK(KF4AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
744         WRITE(LUPRI,'(//A)')'F1 in AO basis in C4FOCK'
745         CALL OUTPUT(WRK(KF5AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
746         WRITE(LUPRI,'(//A)')'F2 in AO basis in C4FOCK'
747         CALL OUTPUT(WRK(KF6AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
748         WRITE(LUPRI,'(//A)')'F3 in AO basis in C4FOCK'
749         CALL OUTPUT(WRK(KF7AO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
750      END IF
751C
752C Transform Fock matrices to MO basis
753C Reuse memory used for density matrices
754C
755      KF2MO = KD1AO
756      KF3MO = KF2MO + N2ORBX
757      KF4MO = KF3MO + N2ORBX
758      KF5MO = KF4MO + N2ORBX
759      KF6MO = KF5MO + N2ORBX
760      KF7MO = KF6MO + N2ORBX
761C
762      CALL FAOMO(ISYMDM(1),WRK(KF1AO),FTTI,CMO,WRK(KFREE),LFREE)
763      CALL FAOMO(ISYMDM(2),WRK(KF2AO),WRK(KF2MO),CMO,WRK(KFREE),LFREE)
764      CALL FAOMO(ISYMDM(3),WRK(KF3AO),WRK(KF3MO),CMO,WRK(KFREE),LFREE)
765      CALL FAOMO(ISYMDM(4),WRK(KF4AO),WRK(KF4MO),CMO,WRK(KFREE),LFREE)
766      CALL FAOMO(ISYMDM(5),WRK(KF5AO),WRK(KF5MO),CMO,WRK(KFREE),LFREE)
767      CALL FAOMO(ISYMDM(6),WRK(KF6AO),WRK(KF6MO),CMO,WRK(KFREE),LFREE)
768      CALL FAOMO(ISYMDM(7),WRK(KF7AO),WRK(KF7MO),CMO,WRK(KFREE),LFREE)
769C
770      IF (IPRRSP.GT.200) THEN
771         WRITE(LUPRI,'(//A)')'F123 in MO basis in C4FOCK'
772         CALL OUTPUT(FTTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
773         WRITE(LUPRI,'(//A)')'F12 in MO basis in C4FOCK'
774         CALL OUTPUT(WRK(KF2MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
775         WRITE(LUPRI,'(//A)')'F13 in MO basis in C4FOCK'
776         CALL OUTPUT(WRK(KF3MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
777         WRITE(LUPRI,'(//A)')'F23 in MO basis in C4FOCK'
778         CALL OUTPUT(WRK(KF4MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
779         WRITE(LUPRI,'(//A)')'F1 in MO basis in C4FOCK'
780         CALL OUTPUT(WRK(KF5MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
781         WRITE(LUPRI,'(//A)')'F2 in MO basis in C4FOCK'
782         CALL OUTPUT(WRK(KF6MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
783         WRITE(LUPRI,'(//A)')'F3 in MO basis in C4FOCK'
784         CALL OUTPUT(WRK(KF7MO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
785      END IF
786C
787C Evaluate FTTI see formula
788C
789C 3*F3 + [k3,Fc] --> F-7
790C 3*F2 + [k2,Fc] --> F-6
791C 3*F1 + [k1,Fc] --> F-5
792C
793      KSYMOP = ISYMD
794      CALL FCKOIN(1,FC,DUMMY,ZYMD,WRK(KF7MO),DUMMY)
795      KSYMOP = ISYMC
796      CALL FCKOIN(1,FC,DUMMY,ZYMC,WRK(KF6MO),DUMMY)
797      KSYMOP = ISYMB
798      CALL FCKOIN(1,FC,DUMMY,ZYMB,WRK(KF5MO),DUMMY)
799C
800C F23 + [k2,F-7] --> F-4
801C F-4 + [k3,F-6] --> F-4
802C
803      CALL OITH1(ISYMC,ZYMC,WRK(KF7MO),WRK(KF4MO),ISYMDM(7))
804      CALL OITH1(ISYMD,ZYMD,WRK(KF6MO),WRK(KF4MO),ISYMDM(6))
805C
806C F13 + [k1,F-7] --> F-3
807C F-3 + [k3,F-5] --> F-3
808C
809      CALL OITH1(ISYMB,ZYMB,WRK(KF7MO),WRK(KF3MO),ISYMDM(7))
810      CALL OITH1(ISYMD,ZYMD,WRK(KF5MO),WRK(KF3MO),ISYMDM(5))
811C
812C F12 + [k2,F-5] --> F-2
813C F-2 + [k1,F-6] --> F-2
814C
815      CALL OITH1(ISYMC,ZYMC,WRK(KF5MO),WRK(KF2MO),ISYMDM(5))
816      CALL OITH1(ISYMB,ZYMB,WRK(KF6MO),WRK(KF2MO),ISYMDM(6))
817C
818C F123 + [k1,F-4] --> FTTI
819C FTTI + [k2,F-3] --> FTTI
820C FTTI + [k3,F-2] --> FTTI
821C
822      CALL OITH1(ISYMB,ZYMB,WRK(KF4MO),FTTI,ISYMDM(4))
823      CALL OITH1(ISYMC,ZYMC,WRK(KF3MO),FTTI,ISYMDM(3))
824      CALL OITH1(ISYMD,ZYMD,WRK(KF2MO),FTTI,ISYMDM(2))
825C
826C Scale with 1/6 ( from definition of E4, without minus sign )
827C
828      CALL DSCAL(NORBT*NORBT,1/D6,FTTI,1)
829C
830      IF (IPRRSP.GT.100) THEN
831         WRITE(LUPRI,'(//A)')'Final result in C4FOCK'
832         CALL OUTPUT(FTTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
833      END IF
834C
835      RETURN
836      END
837C
838      SUBROUTINE C4FCKM(IBCDEQ,FTTI,VECB,VECC,VECD,
839     *                 ZYMB,ZYMC,ZYMD,KZYVB,KZYVC,KZYVD,
840     *                 ISYMA,ISYMB,ISYMC,ISYMD,
841     *                 CMO,FC,MJWOP,WRK,LWRK)
842C
843C PURPOSE:
844C To calculate Fock matrix (FTTI) with triply one-index transformed
845C integrals to be used in direct cubic RPA.
846C July  1997: kruud: Memory efficient version of C4FOCK
847C
848#include "implicit.h"
849#include "dummy.h"
850C
851      DIMENSION FTTI(N2ORBX)
852      DIMENSION VECB(*),VECC(*),VECD(*)
853      DIMENSION ZYMB(NORBT,NORBT),ZYMC(NORBT,NORBT),ZYMD(NORBT,NORBT)
854      DIMENSION CMO(*),FC(*),WRK(*)
855C
856C  INFDIM : NASHDI
857C  INFINP : DIRFCK
858C  WRKRSP : KSYMOP
859C
860#include "maxorb.h"
861#include "maxash.h"
862#include "priunit.h"
863#include "inforb.h"
864#include "infinp.h"
865#include "infvar.h"
866      DIMENSION MJWOP(2,MAXWOP,8)
867#include "wrkrsp.h"
868#include "infrsp.h"
869#include "infpri.h"
870C
871      PARAMETER ( D1 = 1.0D0, D2 = 2.0D0, D6 = 6.0D0, D3 = 3.0D0)
872C
873C Allocate work space for Fock matrices
874C Corresponding Fock matrices in equations
875C
876C   F-1 <--> F123 + permutations
877C   F-2 <--> F12 + F21
878C   F-3 <--> F13 + F31
879C   F-4 <--> F23 + F32
880C   F-5 <--> F1
881C   F-6 <--> F2
882C   F-7 <--> F3
883C
884      IFCTYP = 3
885      NFMATX = 1
886C
887      KFAO  = 1
888      KDAO  = KFAO  + N2BASX
889      KFRE1 = KDAO  + N2BASX
890      KFRE2 = KFRE1 + N2BASX
891      KFREE = KFRE2 + N2BASX
892      LFREE = LWRK  - KFREE
893      IF (LFREE.LT.0) CALL ERRWRK('C4FCKM',LFREE,LWRK)
894C
895      IF (IPRRSP.GT.7) WRITE(LUPRI,'(//A,2(/A,I10))')
896     *' Memory allocation in C4FCKM',
897     *' Available workspace: ', LWRK,
898     *' Allocated          : ', KFREE
899C
900      KFREE = KFRE1
901      LFREE = LWRK - KFREE
902C
903C     Unpack the response vectors
904C
905      CALL GTZYMT(1,VECB,KZYVB,ISYMB,ZYMB,MJWOP)
906      CALL GTZYMT(1,VECC,KZYVC,ISYMC,ZYMC,MJWOP)
907      CALL GTZYMT(1,VECD,KZYVD,ISYMD,ZYMD,MJWOP)
908C
909C     F123 contribution
910C
911      CALL DZERO(WRK(KFAO),4*N2BASX)
912      CALL CDENS3(ISYMB,ISYMC,ISYMD,CMO,ZYMB,ZYMC,ZYMD,
913     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
914      CALL CDENS3(ISYMB,ISYMD,ISYMC,CMO,ZYMB,ZYMD,ZYMC,
915     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
916      CALL CDENS3(ISYMC,ISYMB,ISYMD,CMO,ZYMC,ZYMB,ZYMD,
917     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
918      CALL CDENS3(ISYMC,ISYMD,ISYMB,CMO,ZYMC,ZYMD,ZYMB,
919     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
920      CALL CDENS3(ISYMD,ISYMB,ISYMC,CMO,ZYMD,ZYMB,ZYMC,
921     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
922      CALL CDENS3(ISYMD,ISYMC,ISYMB,CMO,ZYMD,ZYMC,ZYMB,
923     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
924C
925      CALL DSCAL(N2BASX,D2,WRK(KDAO),1)
926      IF (IPRRSP.GT.200) THEN
927         WRITE(LUPRI,'(//A)')'D123 in C4FOCK'
928         CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
929      END IF
930C
931      ISYMDM = ISYMA
932      CALL DZERO(WRK(KFAO),N2BASX)
933      CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
934     *            DIRFCK,WRK(KFREE),LFREE)
935      IF (IPRRSP.GT.200) THEN
936         WRITE(LUPRI,'(//A)')'F123 in AO basis in C4FOCK'
937         CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
938      END IF
939C
940      CALL FAOMO(ISYMDM,WRK(KFAO),FTTI,CMO,WRK(KFREE),LFREE)
941      IF (IPRRSP.GT.200) THEN
942         WRITE(LUPRI,'(//A)')'F123 in MO basis in C4FOCK'
943         CALL OUTPUT(FTTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
944      END IF
945C
946C     F12 matrix
947C
948      CALL DZERO(WRK(KFAO),4*N2BASX)
949      CALL CDENS2(ISYMB,ISYMC,CMO,ZYMB,ZYMC,
950     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
951      CALL CDENS2(ISYMC,ISYMB,CMO,ZYMC,ZYMB,
952     *            WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
953      CALL DSCAL(N2BASX,D6,WRK(KDAO),1)
954      IF (IPRRSP .GT. 200) THEN
955         WRITE(LUPRI,'(//A)')'D12 in C4FOCK'
956         CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
957      END IF
958C
959      ISYMDM = MULD2H(ISYMB,ISYMC)
960      CALL DZERO(WRK(KFAO),N2BASX)
961      CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
962     *            DIRFCK,WRK(KFREE),LFREE)
963      IF (IPRRSP .GT. 200) THEN
964         WRITE(LUPRI,'(//A)')'F12 in AO basis in C4FOCK'
965         CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
966      END IF
967C
968      CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
969      IF (IPRRSP .GT. 200) THEN
970         WRITE(LUPRI,'(//A)')'F12 in MO basis in C4FOCK'
971         CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
972      END IF
973      IF (IBCDEQ .EQ. 24) THEN
974         CALL DSCAL(N2BASX,D3,WRK(KDAO),1)
975      END IF
976      CALL OITH1(ISYMD,ZYMD,WRK(KDAO),FTTI,ISYMDM)
977C
978C     F13 matrix
979C
980      IF (IBCDEQ .NE. 24) THEN
981         CALL DZERO(WRK(KFAO),4*N2BASX)
982         CALL CDENS2(ISYMB,ISYMD,CMO,ZYMB,ZYMD,
983     *               WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
984         CALL CDENS2(ISYMD,ISYMB,CMO,ZYMD,ZYMB,
985     *               WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
986         CALL DSCAL(N2BASX,D6,WRK(KDAO),1)
987         IF (IPRRSP .GT. 200) THEN
988            WRITE(LUPRI,'(//A)')'D13 in C4FOCK'
989            CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
990         END IF
991C
992         ISYMDM = MULD2H(ISYMB,ISYMD)
993         CALL DZERO(WRK(KFAO),N2BASX)
994         CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
995     *        DIRFCK,WRK(KFREE),LFREE)
996         IF (IPRRSP .GT. 200) THEN
997            WRITE(LUPRI,'(//A)')'F13 in AO basis in C4FOCK'
998            CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
999         END IF
1000C
1001         CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
1002         IF (IPRRSP .GT. 200) THEN
1003            WRITE(LUPRI,'(//A)')'F13 in MO basis in C4FOCK'
1004            CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1005         END IF
1006         CALL OITH1(ISYMC,ZYMC,WRK(KDAO),FTTI,ISYMDM)
1007C
1008C     F23 contribution
1009C
1010         CALL DZERO(WRK(KFAO),4*N2BASX)
1011         CALL CDENS2(ISYMC,ISYMD,CMO,ZYMC,ZYMD,
1012     *        WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
1013         CALL CDENS2(ISYMD,ISYMC,CMO,ZYMD,ZYMC,
1014     *        WRK(KDAO),WRK(KFAO),WRK(KFRE1),WRK(KFRE2))
1015         CALL DSCAL(N2BASX,D6,WRK(KDAO),1)
1016         IF (IPRRSP .GT. 20) THEN
1017            WRITE(LUPRI,'(//A)')'D23 in C4FOCK'
1018            CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1019         END IF
1020C
1021         ISYMDM = MULD2H(ISYMC,ISYMD)
1022         CALL DZERO(WRK(KFAO),N2BASX)
1023         CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
1024     *        DIRFCK,WRK(KFREE),LFREE)
1025         IF (IPRRSP .GT. 200) THEN
1026            WRITE(LUPRI,'(//A)')'F23 in AO basis in C4FOCK'
1027            CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1028         END IF
1029C
1030         CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
1031         IF (IPRRSP .GT. 200) THEN
1032            WRITE(LUPRI,'(//A)')'F23 in MO basis in C4FOCK'
1033            CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1034         END IF
1035         CALL OITH1(ISYMB,ZYMB,WRK(KDAO),FTTI,ISYMDM)
1036      END IF
1037C
1038C     F1 contribution
1039C
1040      CALL DZERO(WRK(KFAO),2*N2BASX)
1041      CALL CDENS1(ISYMB,CMO,ZYMB,WRK(KDAO),
1042     *            WRK(KFREE),LFREE)
1043      CALL DSCAL(N2BASX,D6,WRK(KDAO),1)
1044      IF (IPRRSP .GT. 200) THEN
1045         WRITE(LUPRI,'(//A)')'D1 in C4FOCK'
1046         CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1047      END IF
1048C
1049      ISYMDM = ISYMB
1050      CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
1051     *            DIRFCK,WRK(KFREE),LFREE)
1052      IF (IPRRSP .GT. 200) THEN
1053         WRITE(LUPRI,'(//A)')'F1 in AO basis in C4FOCK'
1054         CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1055      END IF
1056C
1057      CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
1058      IF (IPRRSP .GT. 200) THEN
1059         WRITE(LUPRI,'(//A)')'F1 in MO basis in C4FOCK'
1060         CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1061      END IF
1062C
1063      KSYMOP = ISYMB
1064      CALL FCKOIN(1,FC,DUMMY,ZYMB,WRK(KDAO),DUMMY)
1065      CALL DZERO(WRK(KFAO),N2BASX)
1066      CALL OITH1(ISYMD,ZYMD,WRK(KDAO),WRK(KFAO),ISYMDM)
1067      IF (IBCDEQ .EQ. 24) THEN
1068         CALL DSCAL(N2ORBX,D3,WRK(KFAO),1)
1069      END IF
1070      CALL OITH1(ISYMC,ZYMC,WRK(KFAO),FTTI,MULD2H(ISYMB,ISYMD))
1071C
1072      CALL DZERO(WRK(KFAO),N2BASX)
1073      CALL OITH1(ISYMC,ZYMC,WRK(KDAO),WRK(KFAO),ISYMDM)
1074      IF (IBCDEQ .EQ. 24) THEN
1075         CALL DSCAL(N2ORBX,D3,WRK(KFAO),1)
1076      END IF
1077      CALL OITH1(ISYMD,ZYMD,WRK(KFAO),FTTI,MULD2H(ISYMB,ISYMC))
1078C
1079C     F2 contribution
1080C
1081      IF (IBCDEQ .NE. 24) THEN
1082         CALL DZERO(WRK(KFAO),2*N2BASX)
1083         CALL CDENS1(ISYMC,CMO,ZYMC,WRK(KDAO),
1084     *        WRK(KFREE),LFREE)
1085         CALL DSCAL(N2BASX,D6,WRK(KDAO),1)
1086         IF (IPRRSP .GT. 200) THEN
1087            WRITE(LUPRI,'(//A)')'D2 in C4FOCK'
1088            CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1089         END IF
1090C
1091         ISYMDM = ISYMC
1092         CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
1093     *        DIRFCK,WRK(KFREE),LFREE)
1094         IF (IPRRSP .GT. 200) THEN
1095            WRITE(LUPRI,'(//A)')'F2 in AO basis in C4FOCK'
1096            CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1097         END IF
1098C
1099         CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
1100         IF (IPRRSP .GT. 200) THEN
1101            WRITE(LUPRI,'(//A)')'F2 in MO basis in C4FOCK'
1102            CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1103         END IF
1104C
1105         KSYMOP = ISYMC
1106         CALL FCKOIN(1,FC,DUMMY,ZYMC,WRK(KDAO),DUMMY)
1107         CALL DZERO(WRK(KFAO),N2BASX)
1108         CALL OITH1(ISYMD,ZYMD,WRK(KDAO),WRK(KFAO),ISYMDM)
1109         CALL OITH1(ISYMB,ZYMB,WRK(KFAO),FTTI,MULD2H(ISYMC,ISYMD))
1110C
1111         CALL DZERO(WRK(KFAO),N2BASX)
1112         CALL OITH1(ISYMB,ZYMB,WRK(KDAO),WRK(KFAO),ISYMDM)
1113         CALL OITH1(ISYMD,ZYMD,WRK(KFAO),FTTI,MULD2H(ISYMB,ISYMC))
1114C
1115C     F3 contribution
1116C
1117         CALL CDENS1(ISYMD,CMO,ZYMD,WRK(KDAO),WRK(KFREE),LFREE)
1118         CALL DSCAL(N2BASX,D6,WRK(KDAO),1)
1119         IF (IPRRSP .GT. 200) THEN
1120            WRITE(LUPRI,'(//A)')'D3 in C4FOCK'
1121            CALL OUTPUT(WRK(KDAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1122         END IF
1123C
1124         ISYMDM = ISYMD
1125         CALL SIRFCK(WRK(KFAO),WRK(KDAO),NFMATX,ISYMDM,IFCTYP,
1126     *        DIRFCK,WRK(KFREE),LFREE)
1127         IF (IPRRSP .GT. 200) THEN
1128            WRITE(LUPRI,'(//A)')'F3 in AO basis in C4FOCK'
1129            CALL OUTPUT(WRK(KFAO),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI)
1130         END IF
1131C
1132         CALL FAOMO(ISYMDM,WRK(KFAO),WRK(KDAO),CMO,WRK(KFREE),LFREE)
1133         IF (IPRRSP .GT. 200) THEN
1134            WRITE(LUPRI,'(//A)')'F3 in MO basis in C4FOCK'
1135            CALL OUTPUT(WRK(KDAO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1136         END IF
1137C
1138         KSYMOP = ISYMD
1139         CALL FCKOIN(1,FC,DUMMY,ZYMD,WRK(KDAO),DUMMY)
1140         CALL DZERO(WRK(KFAO),N2BASX)
1141         CALL OITH1(ISYMC,ZYMC,WRK(KDAO),WRK(KFAO),ISYMDM)
1142         CALL OITH1(ISYMB,ZYMB,WRK(KFAO),FTTI,MULD2H(ISYMC,ISYMD))
1143C
1144         CALL DZERO(WRK(KFAO),N2BASX)
1145         CALL OITH1(ISYMB,ZYMB,WRK(KDAO),WRK(KFAO),ISYMDM)
1146         CALL OITH1(ISYMC,ZYMC,WRK(KFAO),FTTI,MULD2H(ISYMB,ISYMD))
1147      END IF
1148C
1149C     All contributions now gathered in FTTI, and we finish by
1150C     scaling with 1/6 ( from definition of E4, without minus sign )
1151C
1152      CALL DSCAL(NORBT*NORBT,1/D6,FTTI,1)
1153C
1154      IF (IPRRSP.GT.100) THEN
1155         WRITE(LUPRI,'(//A)')'Final result in C4FOCK'
1156         CALL OUTPUT(FTTI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI)
1157      END IF
1158C
1159      RETURN
1160      END
1161      SUBROUTINE FAOMO(ISYMF,FAO,FMO,CMO,WRK,LWRK)
1162C
1163#include "implicit.h"
1164C
1165      DIMENSION FAO(N2BASX),FMO(N2ORBX),CMO(*),WRK(*)
1166C
1167#include "maxorb.h"
1168#include "maxash.h"
1169#include "inforb.h"
1170#include "infinp.h"
1171#include "wrkrsp.h"
1172C
1173      CALL DZERO(FMO,N2ORBX)
1174      DO 100 ISYM=1,NSYM
1175         JSYM   = MULD2H(ISYM,ISYMF)
1176         NORBI  = NORB(ISYM)
1177         NORBJ  = NORB(JSYM)
1178         IF (NORBI.EQ.0 .OR. NORBJ.EQ.0) GO TO 100
1179         CALL AUTPV(ISYM,JSYM,CMO(ICMO(ISYM)+1),CMO(ICMO(JSYM)+1),
1180     *              FAO,NBAS,NBAST,FMO,NORB,NORBT,WRK,LWRK)
1181 100  CONTINUE
1182      RETURN
1183      END
1184