1#define FCKTRA_DEBUG -1
2!
3!  Dalton, a molecular electronic structure program
4!  Copyright (C) by the authors of Dalton.
5!
6!  This program is free software; you can redistribute it and/or
7!  modify it under the terms of the GNU Lesser General Public
8!  License version 2.1 as published by the Free Software Foundation.
9!
10!  This program is distributed in the hope that it will be useful,
11!  but WITHOUT ANY WARRANTY; without even the implied warranty of
12!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13!  Lesser General Public License for more details.
14!
15!  If a copy of the GNU LGPL v2.1 was not distributed with this
16!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
17!
18!
19
20!  FILE: DALTON/sirius/sirfcktra.F
21!  (c) Copyright Hans Joergen Aa. Jensen, hjj@sdu.dk (2016)
22
23! ===================================================================
24!     sirfcktra.F section 0: Sirius interface
25! ===================================================================
26      subroutine SIR_INTOPEN_FCKTRA()
27!
28!  Written by Hans Joergen Aa. Jensen October 2014
29!
30      implicit none
31#include "priunit.h"
32
33#if FCKTRA_DEBUG > 0
34      write(LUPRI,*) 'SIR_INTOPEN_FCKTRA called'
35#endif
36
37      end
38
39      subroutine SIR_FCKTRA_CTL(TYPE_TEXT, ITRSIR, THR_MOTWOINT, CMO,
40     &   WORK, LWORK)
41!
42!  Written by Hans Joergen Aa. Jensen Apr 2017
43!
44      implicit none
45      integer ITRSIR, LWORK
46      real*8  THR_MOTWOINT, CMO(*), WORK(LWORK)
47#include "priunit.h"
48
49! Used from include files:
50!  gnrinf.h : PARCAL
51!  inforb.h : NNORBX, NORBT, ???
52!  inftra.h : IPRTRA, LBUF
53!  inftap.h : LUINTM, LBINTM
54!  cbihr2.h : IFTHRS
55!  ccom.h   : THRS
56#include "iratdef.h"
57#include "maxaqn.h"
58#include "thrzer.h"
59#include "gnrinf.h"
60#include "inforb.h"
61#include "inftra.h"
62#include "inftap.h"
63#include "cbihr2.h"
64#include "ccom.h"
65
66      character*4 TYPE_TEXT
67      integer JTRSIR, NTRLVL, IOSVAL, MX_NDMAT
68      integer KFREE, LFREE, KFRSAV, KUCMO, KDMAT, KFMAT, I, J, K
69      integer MISH_TEST(8), MASH_TEST(8), IDAMAX, L_H2CD
70      integer, allocatable :: ICDTRA(:), ITRTYP(:)
71      character*8 TABLE(3), LAB123(3), INT_LABEL
72      data        TABLE /'MUL_2INT','DRC_2INT','END_2INT'/
73      data        LAB123/'********','********','********'/
74      logical     FEXIST, OLDDX, FNDLAB
75
76C --- start of executable code
77
78      call QENTER('SIR_FCKTRA_CTL')
79
80      JTRSIR = abs(ITRSIR)
81
82#if FCKTRA_DEBUG > 0
83      IPRTRA = MAX(IPRTRA,FCKTRA_DEBUG)
84      write(lupri,*) 'DEBUG: IPRTRA set to',IPRTRA
85#endif
86
87      IF (IPRTRA .GT. 0) THEN
88         CALL TITLER('.FCKTRA AO->MO integral transformation','*',200)
89         WRITE(LUPRI,'(A,I4,1P,D10.2)')
90     &   ' Sirius integral transformation level & screening threshold:',
91     &   JTRSIR,THR_MOTWOINT
92         FLUSH(LUPRI)
93      END IF
94
95      ! call SIR_INTOPEN_FCKTRA()
96
97C     Test if MO integrals already available
98
99      if (TYPE_TEXT .EQ. 'COUL') then
100         INT_LABEL = TABLE(1)
101         L_H2CD = NNORBT
102      else if (TYPE_TEXT .EQ. 'EXCH') then
103         INT_LABEL = TABLE(2)
104         L_H2CD = N2ORBT
105      else
106         call QUIT('Unknown TYPE_TEXT : '//TYPE_TEXT)
107      end if
108      LBUF = IRAT*L_H2CD
109
110      REWIND(LUINTM)
111      IF ( FNDLAB(INT_LABEL,LUINTM) ) THEN ! true iff existing MOTWOINT is created by .FCKTRA
112         REWIND(LUINTM)
113         READ  (LUINTM)
114         READ  (LUINTM) J,NTRLVL ! originally LBINTM, JTRLVL
115         IF (NTRLVL .EQ. 3 .AND. JTRSIR .EQ. 2) NTRLVL = -1 ! (aa|ii) and (ai|ai) missing
116         IF (NTRLVL .EQ. 2 .AND. JTRSIR .EQ. 3) NTRLVL =  3
117         IF (NTRLVL .EQ. 5 .AND. JTRSIR .EQ. 6) NTRLVL =  6
118         IF (NTRLVL .EQ. 6 .AND. JTRSIR .EQ. 4) NTRLVL = -1 ! (gg|ii) and (gi|gi) missing
119         IF (NTRLVL .EQ. 6 .AND. JTRSIR .EQ. 5) NTRLVL = -1 ! (gg|gi) missing
120         IF (NTRLVL .GE. JTRSIR) THEN
121C           level OK, now read CMO matrix from LUINTM
122C           and subtract from input CMO matrix
123            READ  (LUINTM, IOSTAT=IOSVAL) WORK(1:NCMOT)
124         IF (IOSVAL .EQ. 0) THEN
125            WORK(1:NCMOT) = WORK(1:NCMOT) - CMO(1:NCMOT)
126            I = IDAMAX(NCMOT,WORK,1)
127            READ (LUINTM) MISH_TEST(1:8), MASH_TEST(1:8)
128            MISH_TEST(1:8) = MISH_TEST(1:8) - NISH(1:8)
129            MASH_TEST(1:8) = MASH_TEST(1:8) - NASH(1:8)
130            J = 0
131            DO K = 1,8
132               J = J + ABS(MISH_TEST(K)) + ABS(MASH_TEST(K))
133            END DO
134            IF (ABS(WORK(I)) .LE. THRZER .AND. J .EQ. 0) THEN
135               IF (IPRTRA.GE.0) WRITE (LUPRI,'(/A)')
136     &            ' SIR_FCKTRA_CTL: Integral transformation abandoned,'
137     &          //' the required MO integrals are already available.'
138               GO TO 9999
139            END IF
140         END IF ! IOSVAL test
141         END IF ! NTRLVL test
142      END IF
143
144      IF (THR_MOTWOINT .LT. THRS) THEN
145         WRITE(LUPRI,'(/A,1P,D10.2,A,D10.2)')
146     &   '.FCKTRA INFO: AO screening threshold raised from',
147     &   THR_MOTWOINT,' to the general integral threshold',THRS
148         THR_MOTWOINT = THRS
149         FLUSH(LUPRI)
150      END IF
151
152      allocate(ICDTRA(NNORBX))
153
154      KFREE  = 1
155      LFREE  = LWORK
156      KFRSAV = KFREE
157
158!     translate sirius integral level code to FCKTRA and NEWTRA integral level code
159      call N_TRALVL(JTRSIR, NTRLVL)
160      IF (JTRSIR .GE. 1 .AND. JTRSIR .LE. 4) NTRLVL = NTRLVL + 1
161      ! we need third order, because USEDRC not defined for FCKTRA
162#if FCKTRA_DEBUG > 0
163      write(LUPRI,*) 'DEBUG: SIR_FCKTRA_CTL called, debug level',
164     &   FCKTRA_DEBUG
165      write(LUPRI,*) 'DEBUG: ITRSIR, NTRLVL, LWORK:',ITRSIR,NTRLVL,LWORK
166      flush(LUPRI)
167#endif
168
169!     call N_TRASET(sirntra.F) in order to be able to call N_TRALIM(sirntra.F)
170!     which sets :
171!        ICDTRA: index array for C,D distributions (**|CD) to calculate
172!        ITRTYP(1:NORBT) = number of integral indices in which this orbital
173!                          enters (i.e. 0,1,2,3, or 4)
174      call N_TRASET(-1,LWORK)
175#if FCKTRA_DEBUG > 0
176      IPRTRA = MAX(IPRTRA,FCKTRA_DEBUG)
177      write(lupri,*) 'DEBUG: IPRTRA set to',IPRTRA
178#endif
179      allocate(ITRTYP(NORBT))
180      call N_TRALIM(NTRLVL,ICDTRA,ITRTYP)
181      deallocate ( ITRTYP )
182
183      REWIND(LUINTM)
184      WRITE (LUINTM) (0, I=1,8) ! info not used for FCKTRA
185      WRITE (LUINTM) LBINTM,JTRSIR,NSYM,NORB,NBAS
186      WRITE (LUINTM) CMO(1:NCMOT)
187      WRITE (LUINTM) NISH(1:8), NASH(1:8)
188
189      call GETDAT(LAB123(2),LAB123(3))   ! place date in LAB123(2) and time in LAB123(3)
190      WRITE (LUINTM) LAB123,INT_LABEL
191      WRITE (LUINTM) NNORBX, ICDTRA
192!
193!...  Open direct access file for final MO integrals
194!     (the ICDTRA array which links a C,D to a record has been saved on LUINTM)
195!
196      if (TYPE_TEXT .eq. 'COUL' ) then
197
198         CALL GPINQ('MO2INT_FCKTRA','EXIST',FEXIST)
199         IF (FEXIST) THEN
200            IF (LUMINT .LE. 0) CALL GPOPEN(LUMINT,
201     &         'MO2INT_FCKTRA','OLD','DIRECT',' ',LBUF,OLDDX)
202            CALL GPCLOSE(LUMINT,'DELETE')
203         END IF
204         LUMINT = -1
205         CALL GPOPEN(LUMINT,'MO2INT_FCKTRA','NEW',
206     &      'DIRECT','UNFORMATTED',LBUF,OLDDX)
207      else if (TYPE_TEXT .eq. 'EXCH') then
208
209         CALL GPINQ('DRC2INT_FCKTRA','EXIST',FEXIST)
210         IF (FEXIST) THEN
211            IF (LUMINT .LE. 0) CALL GPOPEN(LUMINT,
212     &         'DRC2INT_FCKTRA','OLD','DIRECT',' ',LBUF,OLDDX)
213            CALL GPCLOSE(LUMINT,'DELETE')
214         END IF
215         LUMINT = -1
216         CALL GPOPEN(LUMINT,'DRC2INT_FCKTRA','NEW',
217     &      'DIRECT','UNFORMATTED',LBUF,OLDDX)
218      else
219         call QUIT('Unknown TYPE_TEXT : '//TYPE_TEXT)
220      end if
221
222
223      CALL FCKTRA_CTL(TYPE_TEXT, ICDTRA, THR_MOTWOINT, CMO, WORK, LWORK)
224
225
226      WRITE (LUINTM) LAB123,TABLE(3)
227      REWIND(LUINTM)
228      CALL GPCLOSE(LUINTM, 'KEEP')
229
230      deallocate ( ICDTRA )
231 9999 call QEXIT('SIR_FCKTRA_CTL')
232      return
233      end
234
235! ===================================================================
236!     sirfcktra.F section 1: calculate MO integrals
237! ===================================================================
238      subroutine FCKTRA_CTL( TYPE_TEXT, ICDTRA, THR_MOTWOINT, CMO,
239     &   WORK, LWORK)
240!
241!  Written by Hans Joergen Aa. Jensen 2016
242!
243      implicit none
244      character*4 TYPE_TEXT
245      integer   ICDTRA(1:NNORBX), LWORK
246      real*8    THR_MOTWOINT, CMO(*), WORK(LWORK)
247#include "priunit.h"
248
249! Used from include files:
250!  gnrinf.h : PARCAL
251!  inforb.h : NNORBX, NORBT, ???
252!  inftra.h : IPRTRA, LBUF
253!  inftap.h : LUSUPM
254!  cbihr2.h : IFTHRS
255!  ccom.h   : THRS
256#include "iratdef.h"
257#include "maxaqn.h"
258#include "gnrinf.h"
259#include "inforb.h"
260#include "inftra.h"
261#include "inftap.h"
262#include "cbihr2.h"
263#include "ccom.h"
264
265      integer MX_NDMAT, LUSUPM_save, IFTHRS_save
266      integer KFREE, LFREE, KFRSAV, KUCMO, KDMAT, KFMAT, I, J
267      integer MISH_TEST(8), MASH_TEST(8), IDAMAX
268      integer, allocatable :: ICD_NODE(:)
269      character*8 TABLE(2),LAB123(3)
270      data        TABLE /'MUL_2INT','END INTM'/
271      data        LAB123/'********','********','********'/
272      logical     FEXIST, OLDDX, FNDLAB
273
274C --- start of executable code
275
276      call QENTER('FCKTRA_CTL')
277
278#if FCKTRA_DEBUG > 0
279      IPRTRA = MAX(IPRTRA,FCKTRA_DEBUG)
280      write(lupri,*) 'FCKTRA_DEBUG: IPRTRA set to',IPRTRA
281      write(lupri,*) 'FCKTRA_DEBUG: THR_MOTWOINT ',THR_MOTWOINT
282      flush(lupri)
283#endif
284
285      IF (IPRTRA .GT. 0) THEN
286!        WRITE(LUPRI,'(A,1P,D10.2)')
287!    &   ' Transformation screening threshold:', THR_MOTWOINT
288         FLUSH(LUPRI)
289      END IF
290
291      IFTHRS_save = IFTHRS
292      IFTHRS = NINT(-LOG10(THR_MOTWOINT))
293
294
295      KFREE  = 1
296      LFREE  = LWORK
297      KFRSAV = KFREE
298
299
300      if (FCKTRA_TYPE .eq. 1 .AND. PARCAL) then
301
302#ifdef VAR_MPI
303         CALL FCKTRA_DISTRIBUTED_MASTER( TYPE_TEXT, ICDTRA, IFTHRS,
304     &      CMO, WORK, LWORK)
305
306#else
307         CALL QUIT('Parallel FCKTRA only possible if compiled with MPI')
308#endif
309      else ! not PARCAL
310
311         call MEMGET2('REAL','UCMO',  KUCMO,NBAST*NORBT,
312     &      WORK,KFREE,LFREE)
313
314         MX_NDMAT = LFREE/(4*N2BASX) ! use half of free memory for DMAT and FMAT
315         MX_NDMAT = MIN(NNORBX,MX_NDMAT)
316         call MEMGET2('REAL','DMAT',  KDMAT,MX_NDMAT*N2BASX,
317     &      WORK,KFREE,LFREE)
318         call MEMGET2('REAL','FMAT',  KFMAT,MX_NDMAT*N2BASX,
319     &      WORK,KFREE,LFREE)
320
321         call UPKCMO(CMO,WORK(KUCMO))
322
323         allocate(ICD_NODE(NNORBX))
324         ICD_NODE(:) = 0
325
326         LUSUPM_save = LUSUPM
327         if (FCKTRA_TYPE .eq. 3) LUSUPM = -9999 ! use RDSUPM
328         if (FCKTRA_TYPE .eq. 4) LUSUPM = -1    ! do not use RDSUPM
329
330         call FCKTRA_FCK_DISTRIBUTED(TYPE_TEXT, ICDTRA, ICD_NODE,
331     &      WORK(KUCMO), MX_NDMAT,
332     &      WORK(KDMAT), WORK(KFMAT), WORK(KFREE), LFREE)
333
334         LUSUPM = LUSUPM_save
335         deallocate ( ICD_NODE )
336
337         call MEMREL('after FCKTRA_FCK_DISTRIBUTED',WORK,
338     &      1,KFRSAV,KFREE,LFREE)
339
340      end if !  if (FCKTRA_TYPE .eq. 1 .AND. PARCAL) then ... else ...
341
342      IFTHRS = IFTHRS_save
343      call QEXIT('FCKTRA_CTL')
344      return
345      end
346
347      subroutine FCKTRA_FCK_DISTRIBUTED( TYPE_TEXT, ICDTRA,
348     &   ICD_NODE, UCMO, MX_NDMAT, DMAT, FMAT, WORK, LWORK)
349!
350!  Written by Hans Joergen Aa. Jensen 2016
351!
352      implicit none
353      character*4 TYPE_TEXT
354      integer ICDTRA(NNORBX), ICD_NODE(NNORBX)
355      integer MX_NDMAT, LWORK
356      real*8  DMAT(N2BASX,MX_NDMAT), FMAT(N2BASX,MX_NDMAT)
357      real*8  UCMO(NBAST,NORBT), WORK(LWORK)
358#include "priunit.h"
359
360! Used from include files:
361!  infpri.h : IPRFCK
362!  inftra.h : IPRTRA
363!  infinp.h : DIRFCK, ADDSRI
364!  inforb.h : N2BASX, NORBT,NNORBX, ...
365!  infind.h : ISMO(:)
366!  inftap.h : LUMINT
367!  dftcom.h : HFXFAC
368!  cbihrs.h : NOSSUP, ONLY_J
369!  infpar.h : MYNUM, MASTER
370!  mtags.h  : MTAG7
371#include "maxorb.h"
372#include "maxash.h"
373#include "infpri.h"
374#include "inftra.h"
375#include "infinp.h"
376#include "inforb.h"
377#include "infind.h"
378#include "inftap.h"
379#include "dftcom.h"
380#include "cbihrs.h"
381#include "infpar.h"
382#include "mtags.h"
383
384      integer ICD, IC, ID, NDMAT, L_H2CD
385      integer, allocatable :: ICD_REC(:), IFCTYP(:), ISYMDM(:)
386      logical FEXIST, OLDDX, FCKTRA_TIMING
387      real*8  TCPU0,TCPU1,TCPU2,TWAL0,TWAL1,TWAL2
388      real*8  TCPU_FCK,TCPU_WR,TWAL_FCK,TWAL_WR
389
390      logical NOSSUP_SAVE, ONLY_J_SAVE, ADDSRI_save
391      logical ABCD_EQ_BACD
392      real*8  HFXFAC_SAVE
393
394      allocate ( ICD_REC(MX_NDMAT) )
395      allocate ( IFCTYP(MX_NDMAT) )
396      allocate ( ISYMDM(MX_NDMAT) )
397
398      FCKTRA_TIMING = IPRTRA .gt. 0  .or. FCKTRA_DEBUG .ge. 0
399      IF (FCKTRA_TIMING) THEN
400         call GETTIM(TCPU0,TWAL0)
401         TCPU_FCK = 0.0D0
402         TCPU_WR  = 0.0D0
403         TWAL_FCK = 0.0D0
404         TWAL_WR  = 0.0D0
405      END IF
406#if FCKTRA_DEBUG > 0
407      write(LUPRI,*) 'FCKTRA_FCK_DISTRIBUTED called'
408      write(LUPRI,*) 'MYNUM, MX_NDMAT, LWORK:',MYNUM,MX_NDMAT,LWORK
409#if FCKTRA_DEBUG > 5
410      write(LUPRI,*) 'Test output of UCMO:'
411      call OUTPUT(UCMO,1,NBAST,1,NORBT,NBAST,NORBT,-1,LUPRI)
412#endif
413      flush(lupri)
414#endif
415
416!     if srdft, only long range integrals
417      ADDSRI_save = ADDSRI
418      ADDSRI = .FALSE.
419
420!     for RDSUPM calls
421      NOSSUP_SAVE = NOSSUP
422      NOSSUP      = .true.
423      ONLY_J_SAVE = ONLY_J
424
425      HFXFAC_SAVE = HFXFAC
426
427      if (TYPE_TEXT(1:4) .eq. 'COUL') then ! Coulomb/Mulliken integrals
428         HFXFAC      = 0.0D0
429         ONLY_J      = .true.
430         ABCD_EQ_BACD= .true.
431         IFCTYP(1:MX_NDMAT) = 11 ! symmetric density matrix, only Coulomb
432      else if (TYPE_TEXT(1:4) .eq. 'EXCH') then ! Exchange/Dirac integrals
433         HFXFAC      = 2.0D0
434         ABCD_EQ_BACD= .false.
435         IFCTYP(1:MX_NDMAT) =  2 ! non-symmetric density matrix, only exchange
436         if (LUSUPM .ne. -1) then
437            call quit(
438     &      'ERROR: FCKTRA "EXCH" integrals not yet for "SUPMAT"')
439         end if
440      else if (TYPE_TEXT(1:4) .eq. 'ANTI') then ! antisymmetrized exchange integrals
441         HFXFAC      = 2.0D0
442         ABCD_EQ_BACD= .false.
443         IFCTYP(1:MX_NDMAT) =  3 ! non-symmetric density matrix, Coulomb and exchange
444         CALL QUIT('FCKTRA ERROR: '//
445     &   'Antisymmetrized Dirac integrals not implemented yet.')
446         if (LUSUPM .ne. -1) then
447            call quit(
448     &      'ERROR: FCKTRA "EXCH" integrals not yet for "SUPMAT"')
449         end if
450      else
451         write(LUPRI,'(///2A)') 'FCKTRA ERROR: '//
452     &   'Unknown MO integral type: ',TYPE_TEXT
453         CALL QUIT('FCKTRA ERROR: '//
454     &   'Requested MO integral type not implemented.')
455      end if
456
457
458      NDMAT = 0
459      ICD   = 0
460      do IC = 1, NORBT
461         do ID = 1, IC
462            ICD = ICD + 1
463#if FCKTRA_DEBUG > 4
464            write (lupri,*) 'mynum, icd, icdtra, icd_node',
465     &      mynum, icd, icdtra(icd), icd_node(icd) ; flush(lupri)
466#endif
467
468            if (ICDTRA(ICD) .gt. 0 .and.
469     &          ICD_NODE(ICD) .eq. MYNUM) then
470
471               NDMAT = NDMAT + 1
472               ISYMDM(NDMAT)  = MULD2H(ISMO(IC), ISMO(ID))
473               ICD_REC(NDMAT) = ICDTRA(ICD)
474               call FCKTRA_CD_DMAT(UCMO(1,IC),UCMO(1,ID),
475     &            DMAT(1,NDMAT),ABCD_EQ_BACD)
476#if FCKTRA_DEBUG > 2
477               write(lupri,*)
478     &           'mynum, NDMAT, IC, ID, ICD, ICDTRA, ISYMDM:',
479     &            mynum,NDMAT,IC,ID,ICD,ICDTRA(ICD), ISYMDM(NDMAT)
480#if FCKTRA_DEBUG > 5
481               write(lupri,*) 'DMAT no. NDMAT'
482               call output(DMAT(1,NDMAT),1,NBAST,1,NBAST,
483     &            NBAST,NBAST,-1,LUPRI)
484#endif
485               flush(lupri)
486#endif
487            end if
488
489            if (NDMAT .eq. MX_NDMAT .or.
490     &         (NDMAT .gt. 0 .and. ICD .eq. NNORBX) ) then
491               call DZERO(FMAT,NDMAT*N2BASX)
492               if (FCKTRA_TIMING) call GETTIM(TCPU1,TWAL1)
493
494               if (LUSUPM .ne. -1) then
495#if FCKTRA_DEBUG > 1
496                  write (lupri,*) 'calling RDSUPM, NDMAT',NDMAT
497                  flush(lupri)
498#endif
499                  call RDSUPM(NDMAT,FMAT,DMAT,ISYMDM,WORK,LWORK)
500               else
501#if FCKTRA_DEBUG > 2
502                  IPRFCK = MAX(IPRFCK,FCKTRA_DEBUG)
503                  write (lupri,*)'calling SIRFCK; NDMAT, IPRFCK, DIRFCK'
504     &            ,NDMAT, IPRFCK, DIRFCK
505                  flush(lupri)
506#endif
507                  call SIRFCK(FMAT,DMAT,NDMAT,ISYMDM,IFCTYP,DIRFCK,
508     &                     WORK,LWORK)
509               end if
510               if (FCKTRA_TIMING) then
511                  call GETTIM(TCPU2,TWAL2)
512#if FCKTRA_DEBUG > 2
513                  write (lupri,*) 'CPU and wall time',
514     &               TCPU2-TCPU1,TWAL2-TWAL1
515                  flush(lupri)
516#endif
517                  TCPU_FCK = TCPU_FCK + TCPU2 - TCPU1
518                  TWAL_FCK = TWAL_FCK + TWAL2 - TWAL1
519               end if
520
521               ! (**|CD) in AO basis now in FMAT
522               if (ABCD_EQ_BACD) then
523                  call FCKTRA_AB_TO_MO(NDMAT,UCMO,FMAT,DMAT)
524                  ! (**|CD) in MO basis now in DMAT
525                  CALL FCKTRA_AB_PACK(NDMAT,DMAT,FMAT,ISYMDM)
526                  ! packed (**|CD) in MO basis now in FMAT
527                  L_H2CD = NNORBT
528               else
529                  call FCKTRA_AB_TO_MO(NDMAT,UCMO,FMAT,FMAT)
530                  ! (**|CD) in MO basis now in FMAT
531                  call quit('packing abcd .ne. bacd not implemented')
532                  L_H2CD = N2ORBT
533               end if
534#if FCKTRA_DEBUG > 2
535               write (lupri,*) 'Back from FCKTRA_AB_TO_MO'; flush(lupri)
536#endif
537               if (MYNUM .ne. 0) then
538                  call MPIXSEND(MYNUM,1,'INTEGER',MASTER,MTAG7)
539                  call MPIXSEND(NDMAT,1,'INTEGER',MASTER,MTAG7)
540                  call MPIXSEND(ICD_REC,NDMAT,'INTEGER',MASTER,MTAG7)
541                  call MPIXSEND(FMAT,NDMAT*L_H2CD,'DOUBLE',MASTER,MTAG7)
542#if FCKTRA_DEBUG > 1
543                  write(lupri,*)'done to master, mtag7',master,mtag7
544                  flush(lupri)
545#endif
546               else
547#if FCKTRA_DEBUG > 1
548                  write(lupri,*) 'writing to MOTWOINT, ndmat',ndmat
549                  flush(lupri)
550#endif
551                  call FCKTRA_WR_LUMINT(LUMINT,NDMAT,FMAT,
552     &               L_H2CD,ICD_REC)
553               end if
554               if (FCKTRA_TIMING) then
555                  call GETTIM(TCPU1,TWAL1)
556                  TCPU_WR  = TCPU_WR  + TCPU1 - TCPU2
557                  TWAL_WR  = TWAL_WR  + TWAL1 - TWAL2
558               end if
559               NDMAT = 0
560            end if
561         end do ! ID = 1, IC
562      end do ! IC = 1, NORBT
563
564      if (FCKTRA_TIMING) then
565         call GETTIM(TCPU1,TWAL1)
566         write(lupri,'(A,2F20.3)') 'CPU and WALL time (s) AO Fock',
567     &      TCPU_FCK,TWAL_FCK
568         write(lupri,'(A,2F20.3)') 'CPU and WALL time (s) AO->MO ',
569     &      TCPU_WR,TWAL_WR
570         TCPU2 = TCPU1 - TCPU0 - TCPU_FCK - TCPU_WR
571         TWAL2 = TWAL1 - TWAL0 - TWAL_FCK - TWAL_WR
572         write(lupri,'(A,2F20.3)') 'CPU and WALL time (s) rest   ',
573     &      TCPU2,TWAL2
574         flush(lupri)
575      end if
576
577      HFXFAC = HFXFAC_SAVE
578      NOSSUP = NOSSUP_SAVE
579      ONLY_J = ONLY_J_SAVE
580      ADDSRI = ADDSRI_save
581
582      deallocate ( ICD_REC )
583      deallocate ( IFCTYP )
584      deallocate ( ISYMDM )
585
586      return
587      end
588! ===================================================================
589      subroutine FCKTRA_CD_DMAT(CMO_C,CMO_D,DMAT_CD,SYM_CD)
590!
591!  Written by Hans Joergen Aa. Jensen 2016
592!
593      implicit none
594      real*8  CMO_C(NBAST), CMO_D(NBAST), DMAT_CD(NBAST,NBAST)
595      logical SYM_CD
596
597! Used from include files:
598!  inforb.h : NBAST
599#include "inforb.h"
600
601      integer JA, JB
602
603      do JB = 1, NBAST
604         do JA = 1,NBAST
605            DMAT_CD(JA,JB) = CMO_C(JA)*CMO_D(JB)
606         end do
607      end do
608
609      if (SYM_CD) then ! symmetrize
610         call DGETSI(NBAST,DMAT_CD,DMAT_CD)
611      end if
612
613      end
614! ===================================================================
615      subroutine FCKTRA_AB_TO_MO(NDMAT,UCMO,H2CD_AO,H2CD_MO)
616!
617!  Written by Hans Joergen Aa. Jensen 2016-2017
618!
619!  Transform (ab|CD) to (AB|CD); ABCD MO indices, ab AO indices
620!
621!  H2CD_AO contains NDMAT "CD" records of (ab|CD) on input,
622!  H2CD_MO contains NDMAT "CD" records of (AB|CD) on output.
623!
624!  Note: H2CD_AO and H2CD_MO may be equivalenced
625!        (i.e. share memory), as N2ORBX .le. N2BASX
626!
627!
628      implicit none
629      integer NDMAT
630      real*8  UCMO(NBAST,NORBT)
631      real*8  H2CD_AO(N2BASX,NDMAT), H2CD_MO(N2ORBX,NDMAT)
632      real*8, allocatable :: H2CD_TMP(:)
633#include "priunit.h"
634
635! Used from include files:
636!  inforb.h : NBAST,NORBT,N2BASX,N2ORBX, ...
637#include "inforb.h"
638
639      integer IDMAT
640
641      allocate (H2CD_TMP(N2BASX))
642
643! TO DO: this should be done on GPU if available /hjaaj Aug. 2016
644      do IDMAT = 1, NDMAT
645         call DGEMM('T','N',NORBT,NBAST,NBAST, 1.0D0,
646     &                 UCMO,NBAST,
647     &                 H2CD_AO(1,IDMAT),NBAST, 0.0D0,
648     &                 H2CD_TMP,NORBT)
649         call DGEMM('N','N',NORBT,NORBT,NBAST, 1.0D0,
650     &                 H2CD_TMP,NORBT,
651     &                 UCMO,NBAST, 0.0D0,
652     &                 H2CD_MO(1,IDMAT),NORBT)
653#if FCKTRA_DEBUG > 5
654            write (LUPRI,*) 'H2CD_MO, record',IDMAT
655            call OUTPUT(H2CD_MO(1,IDMAT),1,NORBT,1,NORBT,
656     &         NORBT,NORBT,-1,LUPRI)
657#endif
658      end do
659
660      deallocate(H2CD_TMP)
661
662      end
663! ===================================================================
664      subroutine FCKTRA_AB_PACK(NDMAT,H2CD_MO,H2CD_MO_PK,ICDSYM)
665!
666!  Written by Hans Joergen Aa. Jensen 2016-2017
667!
668!  Purpose: Pack Mulliken MO integrals before writing them to disk,
669!           using (AB|CD) = (BA|CD) and (AB|CD)=0 if sym(AB) .ne. sym(CD).
670!
671!  H2CD_MO contains NDMAT "CD" records of (AB|CD) on input,
672!  packed into H2CD_MO_PK on output.
673!
674!
675      implicit none
676      integer NDMAT, ICDSYM(NDMAT)
677      real*8  H2CD_MO(N2ORBX,NDMAT), H2CD_MO_PK(NNORBT,NDMAT)
678
679#include "priunit.h"
680
681! Used from include files:
682!  inforb.h : N2ORBX,NNORBT, NSYM, ...
683#include "inforb.h"
684
685      integer IDMAT
686      real*8, allocatable :: H2CD_TMP(:)
687
688      if (NSYM > 1) allocate (H2CD_TMP(NNORBX))
689
690      do IDMAT = 1, NDMAT
691         if (NSYM > 1) then
692            call DGETSP(NORBT, H2CD_MO(1,IDMAT), H2CD_TMP)
693            call TRDPAK( H2CD_TMP, H2CD_MO_PK(1,IDMAT),
694     &         NORB,IORB,NORBT,ICDSYM(IDMAT),1)
695         else
696            call DGETSP(NORBT, H2CD_MO(1,IDMAT), H2CD_MO_PK(1,IDMAT))
697         end if
698#if FCKTRA_DEBUG > 5
699         write (LUPRI,*) 'IDMAT, sym',
700     &      IDMAT,ICDSYM(IDMAT)
701         if (ICDSYM(idmat) == 1) then
702            write (LUPRI,*) 'Packed H2CD_MO'
703            call OUTPKB(H2CD_MO_PK(1,IDMAT),NORB,NSYM,-1,LUPRI)
704         end if
705#endif
706      end do
707
708      if (NSYM > 1) deallocate(H2CD_TMP)
709
710      end
711
712! ===================================================================
713      subroutine FCKTRA_WR_LUMINT(LUMINT,NDMAT,H2CD,L_H2CD,ICD_REC)
714!
715!  Written by Hans Joergen Aa. Jensen 2016
716!
717      implicit none
718      integer LUMINT, NDMAT, L_H2CD, ICD_REC(NDMAT), IOS
719      real*8  H2CD(L_H2CD,NDMAT)
720
721#include "priunit.h"
722
723      integer IDMAT
724
725      do IDMAT = 1, NDMAT
726         ! write NDMAT (**|CD) MO integral distributions to file
727         write(LUMINT, rec=ICD_REC(IDMAT), iostat=IOS)
728     &      H2CD(1:L_H2CD,IDMAT)
729         if (IOS .ne. 0) then
730            call QENTER('FCKTRA_WR_LUMINT')
731            write(lupri,*) 'IDMAT, rec, IOS',IDMAT,ICD_REC(IDMAT),IOS
732            write(lupri,*) 'LUMINT, L_H2CD',LUMINT,L_H2CD
733            call QUIT('IOS .ne. 0')
734         end if
735      end do
736
737      end
738
739! ===================================================================
740!     sirfcktra.F section 2: FCKTRA type 1 routines (if VAR_MPI)
741! ===================================================================
742#ifdef VAR_MPI
743      subroutine FCKTRA_DISTRIBUTED_MASTER( TYPE_TEXT, ICDTRA, IFTHRS,
744     &   CMO, WORK, LWORK)
745!
746!  Written by Hans Joergen Aa. Jensen 2016
747!
748      implicit none
749      character*4 TYPE_TEXT
750      integer ICDTRA(NNORBX), IFTHRS, LWORK
751      real*8  CMO(NCMOT), WORK(LWORK)
752
753#include "priunit.h"
754
755#include "mpif.h"
756
757! Used from include files:
758!  gnrinf.h : CHIVAL
759!  inforb.h : N2BASX, NORBT,NNORBX, ...
760!  infind.h : ISMO(:)
761!  inftra.h : IPRTRA
762!  inftap.h : LUMINT
763!  infpar.h : MYNUM
764!  iprtyp.h : defined parallel calculation types
765!  mtags.h  : MTAG7
766
767#include "maxorb.h"
768#include "maxash.h"
769#include "gnrinf.h"
770#include "inforb.h"
771#include "infind.h"
772#include "inftra.h"
773#include "inftap.h"
774#include "infpar.h"
775#include "iprtyp.h"
776#include "mtags.h"
777
778      integer, allocatable :: ICD_NODE(:)
779      integer IPRTYP, IPRINT_slaves, ICD, NCDTRA, IWHO, NWHO
780      integer KFREE, LFREE
781      integer NDMAT, IDIST, N2GAB, L_H2CD
782
783      call qenter('FCKTRA_DISTRIBUTED_MASTER')
784
785      if (TYPE_TEXT .eq. 'COUL') then
786         L_H2CD = NNORBT
787      else if (TYPE_TEXT .eq. 'EXCH') then
788         L_H2CD = N2ORBT
789      else
790         call QUIT('Unknown TYPE_TEXT '//TYPE_TEXT)
791      end if
792
793!     start nodes on the distributed 2-electron integral transformation task
794
795      IPRTYP = CALL_FCKTRA_DISTRIBUTED ! get code number from iprtyp.h
796      call MPIXBCAST(IPRTYP,1,'INTEGER',MASTER)
797      IPRINT_slaves = IPRTRA ! TODO set print level in input
798      call MPIXBCAST(IPRINT_slaves,1,'INTEGER',MASTER)
799
800!     Now nodes are in subroutine FCKTRA_DISTRIBUTED_NODE.
801!     First make sure basis set information is known to nodes
802
803      KFREE = 1
804      LFREE = LWORK
805      CALL HER_sendreceive_molinfo(IPRINT_slaves,WORK,KFREE,LFREE)
806
807!     Transfer integral transformation information to nodes
808
809      allocate( ICD_NODE(NNORBX) )
810      ICD_NODE(:) = 0
811      NCDTRA = 0
812      DO ICD = 1,NNORBX
813         IF (ICDTRA(ICD) .EQ. 0) CYCLE
814         ICD_NODE(ICD) = MOD(NCDTRA,NODTOT) + 1
815         NCDTRA = NCDTRA + 1
816      END DO
817
818#if FCKTRA_DEBUG > 2
819      write (*,*) 'FCKTRA_DISTRIBUTED_MASTER has started slaves'
820      write (lupri,*) 'FCKTRA_DISTRIBUTED_MASTER has started slaves'
821      write (lupri,*) 'IPRTYP,IPRINT_slaves=',IPRTYP,IPRINT_slaves
822      write (lupri,*) 'ICD_NODE is'
823      write (lupri,'(10I4,5X,10I4)') ICD_NODE(1:NNORBX)
824      flush(lupri)
825#endif
826      call MPIXBCAST(TYPE_TEXT,4,'CHARACTER',MASTER)
827      call MPIXBCAST(IFTHRS,1,'INTEGER',MASTER)
828      call MPIXBCAST(CHIVAL,1,'DOUBLE',MASTER)
829      call MPIXBCAST(ICDTRA,NNORBX,'INTEGER',MASTER)
830      call MPIXBCAST(ISMO,NORBT,'INTEGER',MASTER)
831      call MPIXBCAST(ICD_NODE,NNORBX,'INTEGER',MASTER)
832      call MPIXBCAST(CMO,NCMOT,'DOUBLE',MASTER)
833
834      IDIST = 0
835  100 CONTINUE
836
837      IWHO = -1
838      call MPIXRECV(NWHO,1,'INTEGER',IWHO,MTAG7)
839      call MPIXRECV(NDMAT,1,'INTEGER',NWHO,MTAG7)
840      IF (NDMAT .LE. 0) THEN
841         write(LUPRI,*)
842     &   'ERROR: master received ndmat.le.0, ndmat, nwho',ndmat,nwho
843         call quit('NDMAT .le. 0')
844      END IF
845      ! we reuse ICD_NODE to store ICD_REC
846      ! CALL MPIXRECV(ICD_REC,NDMAT,'INTEGER',NWHO,MTAG7)
847      CALL MPIXRECV(ICD_NODE,NDMAT,'INTEGER',NWHO,MTAG7)
848      IF (NDMAT*L_H2CD .GT. LWORK) CALL QUIT('NDMAT*L_H2CD .gt. lwork')
849      CALL MPIXRECV(WORK,NDMAT*L_H2CD,'DOUBLE',NWHO,MTAG7)
850
851      call FCKTRA_WR_LUMINT(LUMINT,NDMAT,WORK,L_H2CD,ICD_NODE)
852
853#if FCKTRA_DEBUG > 1
854      write (lupri,*) 'Received from nwho =',nwho
855      write (lupri,*) 'Received and saved ndmat =',ndmat
856      write(lupri,'(A,(20I5))') 'The CD records:',ICD_NODE(1:NDMAT)
857#endif
858
859      IDIST = IDIST + NDMAT
860      IF (IDIST .GT. NCDTRA) THEN
861         write(lupri,*) 'ERROR: IDIST .gt. NCDTRA',IDIST,NCDTRA
862         call quit('IDIST .gt. NCDTRA, uha')
863      else if (IDIST .LT. NCDTRA) then
864         go to 100
865      end if
866
867      deallocate( ICD_NODE )
868      call qexit('FCKTRA_DISTRIBUTED_MASTER')
869      return
870      end subroutine FCKTRA_DISTRIBUTED_MASTER
871! ===================================================================
872      subroutine FCKTRA_DISTRIBUTED_NODE(WORK,LWORK,IPRINT)
873!
874!  Written by Hans Joergen Aa. Jensen 2016
875!
876      implicit none
877      integer LWORK, IPRINT
878      real*8  WORK(LWORK)
879
880      integer, allocatable :: ICDTRA(:), ICD_NODE(:)
881      integer ICD, NCDTRA, MX_NDMAT
882      integer KFREE, LFREE, KDMAT, KFMAT, KUCMO
883      logical PARHER_save, SLAVE_save, DIRFCK_save
884      character*4 TYPE_TEXT
885#include "priunit.h"
886
887! Used from include files:
888!  gnrinf.h : PARCAL, CHIVAL
889!  infinp.h : DIRFCK
890!  infind.h : ISMO(:)
891!  inforb.h : NORBT, NNORBX, NCMOT
892!  infpar.h : MYNUM,SLAVE,PARHER,?
893!  inftap.h : LUSUPM
894!  cbihr2.h : IFTHRS
895#include "maxorb.h"
896#include "maxash.h"
897#include "gnrinf.h"
898#include "infinp.h"
899#include "infind.h"
900#include "inforb.h"
901#include "infpar.h"
902#include "inftap.h"
903#include "cbihr2.h"
904
905!  mpif.h   : for MPI
906#include "mpif.h"
907
908!     First make sure basis set information is known to this node
909
910      KFREE = 1
911      LFREE = LWORK
912      CALL HER_sendreceive_molinfo(IPRINT,WORK,KFREE,LFREE)
913
914#if FCKTRA_DEBUG > 1
915      write (lupri,*)
916     &'FCKTRA_DISTRIBUTED_NODE has started, MYNYM =',MYNUM
917      write (lupri,*) 'NORBT, NNORBX, NCMOT',
918     &                 NORBT, NNORBX, NCMOT
919      flush(lupri)
920#endif
921
922!    We check that orbital information has been transferred to slaves,
923!    and then we call SETORB to define inforb.h and infind.h (we need ISMO(:))
924
925      IF (NORBT .LE. 0) THEN
926         CALL QUIT(
927     &   'ERROR: Orbital information not transferred to slaves.')
928      END IF
929      CALL SETORB
930
931!    Some set up for this task (AO to MO transformation)
932
933      KFREE = 1
934      LFREE = LWORK
935
936      SLAVE_save  = SLAVE
937      SLAVE  = .FALSE. ! Be sure the Fock matrices are NOT constructed in parallel,
938
939      PARHER_save = PARHER
940      PARHER = .FALSE. ! be sure the Fock matrices are NOT constructed in parallel
941                       ! we are parallelizing on the distributions
942
943      DIRFCK_save = DIRFCK
944      DIRFCK = .TRUE.  ! NEVER read from AOTWOINT when parallel
945
946      LUSUPM = -1      ! Do not use AO integrals from file AO2_JINT
947
948
949      allocate( ICDTRA(NNORBX) )
950      allocate( ICD_NODE(NNORBX) )
951
952      call MEMGET2('REAL','UCMO',  KUCMO,NBAST*NORBT,
953     &      WORK,KFREE,LFREE)
954
955      call MPIXBCAST(TYPE_TEXT,4,'CHARACTER',MASTER)
956      call MPIXBCAST(IFTHRS,1,'INTEGER',MASTER)            ! accuracy of MO integrals (screening)
957      call MPIXBCAST(CHIVAL,1,'DOUBLE',MASTER)
958      call MPIXBCAST(ICDTRA,NNORBX,'INTEGER',MASTER)
959      call MPIXBCAST(ISMO,NORBT,'INTEGER',MASTER)
960      call MPIXBCAST(ICD_NODE,NNORBX,'INTEGER',MASTER)
961
962      call MPIXBCAST(WORK(KFREE),NCMOT,'DOUBLE',MASTER)
963      call UPKCMO(WORK(KFREE),WORK(KUCMO))
964
965      NCDTRA = 0
966      DO ICD = 1,NNORBX
967         IF (ICD_NODE(ICD) .EQ. MYNUM) NCDTRA = NCDTRA + 1
968      END DO
969
970      IF (NCDTRA .EQ. 0) THEN
971      ! nothing to do for me this time ...
972         GO TO 1000
973      END IF
974
975      MX_NDMAT = LFREE/(4*N2BASX) ! use up to half of free memory for DMAT and FMAT
976      MX_NDMAT = MIN(NCDTRA,MX_NDMAT)
977
978      IF (IPRINT .GT. 0) THEN
979         write(lupri,'(A,2I12)') 'FCKTRA # of distributions and'//
980     &   ' # of distributions per load',NCDTRA, MX_NDMAT
981         flush(lupri)
982      END IF
983
984      IF (MX_NDMAT .EQ. 0) THEN
985         CALL QUIT('MX_NDMAT = 0, too little memory')
986      END IF
987      call MEMGET2('REAL','DMAT', KDMAT,MX_NDMAT*N2BASX,
988     &     WORK,KFREE,LFREE)
989      call MEMGET2('REAL','FMAT', KFMAT,MX_NDMAT*N2BASX,
990     &     WORK,KFREE,LFREE)
991
992#if FCKTRA_DEBUG > 2
993      write(lupri,*) 'NCDTRA,MX_NDMAT,NNORBX =',NCDTRA,MX_NDMAT,NNORBX
994      write (lupri,'(10I4,5X,10I4)') ICD_NODE(1:NNORBX)
995      write(lupri,*) 'chival',chival
996      write (lupri,*) mynum,'node; ICD_NODE is'
997      write (lupri,*) MYNUM,' is calling fcktra_fck_distributed'
998      flush(lupri)
999#endif
1000
1001      CALL  FCKTRA_FCK_DISTRIBUTED( TYPE_TEXT,
1002     &      ICDTRA, ICD_NODE, WORK(KUCMO), MX_NDMAT,
1003     &      WORK(KDMAT), WORK(KFMAT), WORK(KFREE), LFREE)
1004
1005
1006 1000 CONTINUE
1007
1008      deallocate( ICDTRA )
1009      deallocate( ICD_NODE )
1010
1011      PARHER = PARHER_save
1012      SLAVE  = SLAVE_save
1013      DIRFCK = DIRFCK_save
1014
1015      return
1016      end subroutine FCKTRA_DISTRIBUTED_NODE
1017#endif   /* VAR_MPI */
1018
1019! ===================================================================
1020!     sirfcktra.F section 3: read MO integrals
1021! ===================================================================
1022      subroutine FCKTRA_NXTH2M(IC,ID,H2CD,NEEDTP,WORK,IDIST)
1023!
1024!  Written by Hans Joergen Aa. Jensen 2016
1025!  This version is interface routine for .FCKTRA integral transformation.
1026!
1027! Purpose:
1028!    Read next Mulliken two-electron integral distribution (**|cd)
1029!    where |cd) distribution is needed according to NEEDTP(ITYPCD)
1030!
1031! Usage:
1032! 1) Set IDIST = 0 before first call.
1033!    Do **NOT** change IDIST in calling routine
1034!    until last distribution has been read (signalled by IDIST .eq. -1)
1035! 2) if IDIST<0 on input, then IC, ID must be set to correspond to
1036!    the distribution wanted in H2CD
1037!
1038
1039      implicit none
1040#include "priunit.h"
1041
1042      integer  :: IC, ID, NEEDTP(-4:6), IDIST
1043      real*8   :: H2CD(NORBT,NORBT), WORK(*)
1044
1045! Used from include files:
1046!  inforb.h : NORBT, NNORBT, NNORBX, ???
1047!  infind.h : ISMO(:)
1048!  inftra.h : LBUF
1049!  inftap.h : LUINTM, LUMINT
1050#include "maxorb.h"
1051#include "maxash.h"
1052#include "dummy.h"
1053#include "iratdef.h"
1054#include "inforb.h"
1055#include "infind.h"
1056#include "inftra.h"
1057#include "inftap.h"
1058#include "orbtypdef.h"
1059
1060      integer       :: MMORBX, ICD, ICD_1, ICDSYM
1061      integer       :: ITYPC, ITYPD, ITYPCD
1062      integer, save :: ICD_SAVE
1063      integer, save, allocatable :: ICDTRA(:)
1064      logical       :: OLDDX
1065
1066#if FCKTRA_DEBUG > 3
1067      write(LUPRI,*) 'FCKTRA_NXTH2M called, IC, ID, IDIST =',
1068     &   IC, ID, IDIST
1069      flush(LUPRI)
1070#endif
1071
1072      if (IDIST .le. 0) then
1073         if (allocated(ICDTRA)) deallocate(ICDTRA)
1074         allocate(ICDTRA(1:NNORBX))
1075         if (LUINTM .le. 0) THEN
1076            CALL GPOPEN(LUINTM,FNINTM,'OLD',' ',
1077     &                  'UNFORMATTED',IDUMMY,.FALSE.)
1078         END IF
1079         LBUF = IRAT*NNORBT
1080         IF (LUMINT .LE. 0) THEN
1081            CALL GPOPEN(LUMINT,'MO2INT_FCKTRA','OLD','DIRECT',' ',
1082     &         LBUF,OLDDX)
1083         END IF
1084         rewind (LUINTM)
1085         call MOLLAB('MUL_2INT',LUINTM,LUPRI)
1086         read (LUINTM) MMORBX, ICDTRA(1:MMORBX)
1087         if (NNORBX .ne. MMORBX) then
1088            call QENTER('FCKTRA_NXTH2M')
1089            write (LUPRI,'(//A,2I10)')  'FATAL ERROR '//
1090     &      'NNORBX on MO2INT_FCKTRA .ne. NNORBX',MMORBX,NNORBX
1091            call QUIT(
1092     &      'NNORBX on MO2INT_FCKTRA .ne. NNORBX in FCKTRA_NXTH2M')
1093         end if
1094         ICD_SAVE = 0
1095      end if
1096
1097      if (IDIST .lt. 0) then
1098
1099         if (IC .le. 0 .or. IC .gt. NORBT .and.
1100     &       ID .le. 0 .or. ID .gt. NORBT) then
1101            write (LUPRI,*) 'FCKTRA_NXTH2M: invalid input IC, ID',IC,ID
1102            call quit('FCKTRA_NXTH2M: invalid input IC, ID')
1103         end if
1104         if (IC .ge. ID) then
1105            ICD = IC*(IC-1)/2 + ID
1106         else
1107            ICD = ID*(ID-1)/2 + IC
1108         end if
1109         IDIST = ICDTRA(ICD)
1110         if (IDIST .le. 0) then
1111            call QUIT('IC,ID record no. .le. 0 in FCKTRA_NXTH2M')
1112         end if
1113         deallocate(ICDTRA)
1114
1115      else ! IDIST .ge. 0
1116
1117  100    continue
1118         IDIST = -1
1119
1120         ICD_1 = ICD_SAVE + 1
1121         do ICD = ICD_1,NNORBX
1122            if (ICDTRA(ICD) .gt. 0) then
1123               IDIST    = ICDTRA(ICD)
1124               ICD_SAVE = ICD
1125               exit
1126            end if
1127         end do
1128         if (IDIST .eq. -1) then
1129            ! finished
1130            deallocate(ICDTRA)
1131            return
1132         end if
1133
1134         do IC = 1, NORBT
1135            ICD = IC*(IC+1)/2
1136            if (ICD .lt. ICD_SAVE) cycle
1137            ID = ICD_SAVE - ICD + IC
1138            exit
1139         end do
1140         ITYPC  = IOBTYP(IC)
1141         ITYPD  = IOBTYP(ID)
1142         ITYPCD = IDBTYP(ITYPC,ITYPD)
1143         if ( NEEDTP(ITYPCD) .eq. 0 ) go to 100
1144
1145      end if
1146
1147         ICDSYM = MULD2H( ISMO(IC), ISMO(ID) )
1148         read ( LUMINT, rec=IDIST ) WORK(1:NNORBT)
1149         call TRDPAK(WORK(1+NNORBT),WORK,NORB,IORB,NORBT,
1150     &      ICDSYM,-1)
1151         call DSPTGE(NORBT,WORK(1+NNORBT),H2CD)
1152
1153#if FCKTRA_DEBUG > 9
1154         write (LUPRI,*) 'FCKTRA_NXTH2M, rec,IC,ID',IDIST,IC,ID
1155         write (LUPRI,*) 'Unpacked H2CD, ICDSYM',ICDSYM
1156         call output(H2CD,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI)
1157         if (ICDSYM.eq.1) then ! only print some of them, of ICDSYM.eq.1
1158            write (LUPRI,*) 'Packed H2CDPK'
1159            call outpkb(WORK,NORB,NSYM,-1,LUPRI)
1160         end if
1161#endif
1162
1163      return
1164      end
1165
1166      subroutine FCKTRA_NXTH2D(IC,ID,H2CD,NEEDTP,IDIST)
1167!
1168!  Written by Hans Joergen Aa. Jensen 2016
1169!  This version is interface routine for .FCKTRA integral transformation.
1170!
1171! Purpose:
1172!    Read next Dirac two-electron integral distribution <**|cd> = (*c|*d)
1173!    where (cd) distribution is needed according to NEEDTP(ITYPCD)
1174!
1175!
1176
1177      implicit none
1178#include "priunit.h"
1179
1180      integer  IC, ID, NEEDTP(*), IDIST
1181      real*8   H2CD(NORBT,NORBT)
1182
1183! Used from include files:
1184!  inforb.h : NORBT, ???
1185#include "inforb.h"
1186
1187      call quit('ERROR: FCKTRA_NXTH2D is not implemented')
1188
1189
1190      end
1191! -- end of DALTON/sirius/sirfcktra.F --
1192