C -*- Fortran -*-
C
C  DELCD Debug
C  DELCR compare with Dalton integrals
C  DELCX all integrals via the 4-index transformation
C
C integral generator - no symmetry, kinetic, Hamilton and
C bi-electronic integrals
C
C..FILE 'common_syst.h'
C
C..FILE 'common_flow.h'
C..FILE 'constants.h'
C
C..FILE 'common_spheri.h'
C
C..FILE 'common_locat.h'
C
C..FILE 'common_polyta.h'
C
C..FILE 'common_tabaux.h'
C
C..FILE 'common_integ.h'
 
 
C..FILE 'dimens_intu.h'
C
C $Log: intcal.r,v $
C Revision 1.31  2011/03/18 15:32:41  edrisse
C Include option for penetration integrals PEN_INT,
C as bielectronic integrals (aa|aa) and (aa|bb)
C
C Revision 1.30  2011/03/13 21:19:11  reinh
C DIPY in the dipole integrals corrected (was already corrected but not
C included in the sources)
C
C Revision 1.29  2011/03/13 00:15:35  reinh
C Cartesian bielectronic integrals in the same order as
C in Dalton
C
C Revision 1.28  2011/01/31 01:10:17  reinh
C Some corrections in the cartesian part, which should not affect the calculations in solid harmonics.
C
C Revision 1.27  2009/03/17 22:40:49  reinh
C Option LKINET added for having only kinetic-energy integrals
C
C Revision 1.26  2008/10/06 19:43:21  reinh
C spherical harmonics for g-functions added
C
C Revision 1.25  2008/04/30 11:40:45  reinh
C reprogrammed completely the Gamma function
C
C Revision 1.24  2008/04/29 23:52:50  reinh
C now we are at the same position as with the CRYSTAL Gamma function
C
C Revision 1.23  2008/04/29 23:29:03  reinh
C nearly finished, the replacement of BUXTAB/GAMF88
C
C Revision 1.22  2008/04/29 18:06:55  reinh
C own Gamma table routine
C
C Revision 1.21  2008/04/28 18:02:30  reinh
C Do not touch BUXTAB !
C
C Revision 1.20  2008/04/17 14:23:35  reinh
C Added dipole integrals
C
C Revision 1.19  2008/04/16 13:28:07  reinh
C more polite error message
C
C Revision 1.18  2008/04/16 13:25:54  reinh
C no f functions calculatable
C
C Revision 1.17  2008/04/16 13:04:23  reinh
C litte y removed
C
C Revision 1.16  2008/04/16 13:01:24  reinh
C The working version, finally !!
C
C Revision 1.15  2008/04/15 19:32:39  reinh
C Add Hamilton and kinetic to write to HAMILTO
C
C Revision 1.14  2008/04/15 19:19:56  reinh
C This sould be the final version, still puzling about the f
C
C Revision 1.13  2008/04/14 21:30:13  reinh
C seems working, only problems for the (f0 f0 | f3 f3) integrals
C
C Revision 1.12  2008/04/14 08:45:44  reinh
C working for (ss|ss), (ss|cd) and (ab|ss)
C
C Revision 1.11  2008/04/11 16:27:33  reinh
C only one :log
C
C Revision 1.10  2008/04/11 16:26:35  reinh
C Still at (ss|ss), but closer to the goal
C
C Revision 1.9  2008/04/11 00:33:21  reinh
C works still for s only
C
C
C $Id: intcal.r,v 1.31 2011/03/18 15:32:41 edrisse Exp $
C
      PROGRAM INTCAL
C
      INCLUDE 'param.h'                                                 
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
C integrals in hermitian Gaussians                                      
      COMMON /INTEGH/ R000(LHMAX)                                       
C..   INCLUDE 'common_integ.h'                                          
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
C
      CHARACTER*8 PNAME
      LOGICAL LREAD,LN8,LN5
      CHARACTER*3 CLSTRG(0:NLMAX-1,NLMAX+NLMAX+1),CBSTR(NBASM)
      CHARACTER*1 CSTYP(0:5)
 
      DATA CSTYP /'S','P','D','F','G','H'/
 
      INCLUDE 'compiler_stamp'                                          
      X=CPTIME(3)
      PNAME='INTCAL  '
      CALL DATING(PNAME,1)
      WRITE(6,*)
      WRITE(6,*) '           I N T C A L '
      WRITE(6,*) '       integral generator '
      WRITE(6,*) '          molecular case '
      WRITE(6,*) '     P.Reinhardt, Paris 6/2002, 4/2008 '
      WRITE(6,*)
      WRITE(6,*)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ORTHO - programs for ab-initio calculations in localized orbitals
C Copyright (C) 2008 Peter Reinhardt (Univ. Paris VI, France)
C
C This program is free software: you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation, either version 3 of the License, or
C     (at your option) any later version.
C
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with this program.  If not, see <http://www.gnu.org/licenses/>.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      CLSTRG(0,1)='S  '
      CLSTRG(1,1)='PX '
      CLSTRG(1,2)='PY '
      CLSTRG(1,3)='PZ '
      CLSTRG(2,1)='D-2'
      CLSTRG(2,2)='D-1'
      CLSTRG(2,3)='D 0'
      CLSTRG(2,4)='D 1'
      CLSTRG(2,5)='D 2'
      CLSTRG(3,1)='F-3'
      CLSTRG(3,2)='F-2'
      CLSTRG(3,3)='F-1'
      CLSTRG(3,4)='F 0'
      CLSTRG(3,5)='F 1'
      CLSTRG(3,6)='F 2'
      CLSTRG(3,7)='F 3'
      CLSTRG(4,1)='G-4'
      CLSTRG(4,2)='G-3'
      CLSTRG(4,3)='G-2'
      CLSTRG(4,4)='G-1'
      CLSTRG(4,5)='G 0'
      CLSTRG(4,6)='G 1'
      CLSTRG(4,7)='G 2'
      CLSTRG(4,8)='G 3'
      CLSTRG(4,9)='G 4'
C
C numerical value - worse for comparison with Dalton
c$$$      PI=3.141592653589793238D0
c$$$      TP25=34.98683665524972567962
 
      PI=2.D0*ACOS(0.D0)
      TP25=2.D0*(PI**2.5D0)
      ONE=1.D0
C
      CALL RDINP
C
C reading options
C
      WRITE(6,*)
      WRITE(6,*) ' INTEGRAL THRESHOLD IS ',THRESH
      WRITE(6,*)
 
      IF (LFORM) THEN
       WRITE(6,*) ' BIELECTRONIC INTEGRALS ARE WRITTEN FORMATTED '
      ELSE
       WRITE(6,*) ' BIELECTRONIC INTEGRALS ARE WRITTEN UNFORMATTED '
      END IF
      WRITE(6,*)
 
      IF (LKINET) THEN
       WRITE(6,*) ' Calculating only kinetic-energy integrals '
       LMONO=.FALSE.
       LBIEL=.FALSE.
      END IF
 
C
      IF (LMONO) THEN
       WRITE(6,*) ' OVERLAP and KINETIC ENERGY INTEGRALS'
       WRITE(6,*) ' ELECTRON-NUCLEUS INTERACTION INTEGRALS '
      END IF
 
      IF (LBIEL) THEN
       WRITE(6,*) ' ELECTRON-ELECTRON REPULSION INTEGRALS '
      END IF
 
      IF (LCART) THEN
       WRITE(6,*) ' INTEGRALS ARE OVER CARTESIAN FUNCTIONS (6D, 10F) '
      ELSE
       WRITE(6,*) ' INTEGRALS ARE OVER SPHERICAL HARMONICS (5D, 7F) '
      END IF
C
C ----------------------------------------------------------------
C
      IWFILE=22
      IUNITR=26
      ANTOAU=0.529177247D0
      PI=2.D0*ACOS(0.0D0)
C
C read information on system
C
      OPEN(UNIT=IUNITR,FILE='SYSTEM.ORTHO',STATUS='OLD',
     -     FORM='FORMATTED',ERR=901)
      READ(IUNITR,*) NATOM
      WRITE(6,*)
      NSHLT=0
      NBAS=1
      LMAX=0
      NST=1
      ISHL=0
      DO IAT=1,NATOM
C atomic positions are in a.u.
       READ(IUNITR,*) NUCKEL,NSHLAT,(POS(J,IAT),J=1,3)
       NZ(IAT)=NUCKEL
       NSH(IAT)=NSHLAT
       NSHLT=NSHLT+NSHLAT
 9931  FORMAT('  Atom No ',I3,' at ',3F20.12,/,'   Nuclear charge ',I3
     $      ,'  No of shells ',I5)
       WRITE(6,9931) IAT,(POS(J,IAT),J=1,3),NUCKEL,NSHLAT
       DO ISH=1,NSHLAT
        ISHL=ISHL+1
        NZENTR(ISHL)=IAT
        NPST(ISHL)=NST
        READ(IUNITR,*) ITYPE,NPRIMM
 9932   FORMAT('    Shell No ',I3,' of type ',A1,'  with  ',I4,' primiti
     $ves ')
C folded 1 (fixf $Revision: 1.3 $) 
        WRITE(6,9932) ISH,CSTYP(ITYPE),NPRIMM
        IL(ISHL)=ITYPE
        NPRIM(ISHL)=NPRIMM
        LMAX=MAX(LMAX,ITYPE)
        IF (LCART) THEN
         NBAS=NBAS+(ITYPE+1)*(ITYPE+2)/2
        ELSE
         DO I=1,ITYPE+ITYPE+1
          CBSTR(NBAS)=CLSTRG(ITYPE,I)
          NBAS=NBAS+1
         END DO
        END IF
        DO III=1,NPRIMM
         READ(IUNITR,*) EXX(NST),COEFF(NST)
         NST=NST+1
        END DO
       END DO
      END DO
      NST=NST-1
      NBAS=NBAS-1
      CLOSE(IUNITR)
C
C write information
C
      WRITE(6,*)
CD    WRITE(6,*) ' atomic information '
CD    DO I=1,NATOM
CD     WRITE(6,9901) I,NSH(I),NZ(I),(POS(J,I),J=1,3)
 9901  FORMAT(' ATOM ',I4,': NSH ',I4,' NZ ',I4,' at (a.u.) ',3F20.12)
CD    END DO
CD    WRITE(6,*)
CD    WRITE(6,*) ' information per shell'
CD    DO I=1,NSHLT
CD     WRITE(6,9902) I,NPRIM(I),IL(I),NZENTR(I),NPST(I)
 9902  FORMAT(' SHELL ',I4,'  NPRIM',I4,'    L',I2,'  NZENTR',I4
     $      ,'  NPST ',I4)
CD    END DO
CD    WRITE(6,*)
CD    WRITE(6,*) ' information per primitive'
CD    DO I=1,NST
CD     WRITE(6,9903) I,EXX(I),COEFF(I)
 9903  FORMAT(' primitive ',I4,' EXPONENT ',G20.12,' COEFF ',G20.12)
CD    END DO
C
      WRITE(6,*)
      WRITE(6,*) ' READ INFORMATION ON SYSTEM '
      WRITE(6,9946) NATOM,NATMX,NSHLT,NSHLMX,NBAS,NBASM,LMAX,NLMAX
 9946 FORMAT('             ACTUAL      MAXIMUM ',/,
     -     '   NATOMS   ',2I8,/,
     -     '   NSHLT    ',2I8,/,
     -     '   NBAS     ',2I8,/,
     -     '   LMAX     ',2I8)
      IF (NBAS.GT.NBASM) THEN
       WRITE(6,*) ' NBAS (',NBAS,' ) > NBASM (',NBASM,' )'
       STOP 'INCREASE NBASM'
      END IF
      IF (NSHLT.GT.NSHLMX) THEN
       WRITE(6,*) ' NSHLT (',NSHLT,' ) > NSHLMX (',NSHLMX,' )'
       STOP 'INCREASE NSHLMX'
      END IF
 
c$$$      IF (LMAX.GT.2) THEN
c$$$       WRITE(6,*)
c$$$       WRITE(6,*) ' we cannot calculate with f functions or beynd '
c$$$       WRITE(6,*) '  ... there is still an error for '
c$$$     $      ,'the integrals (f0 f0 | f3 f3) '
c$$$       WRITE(6,*)
c$$$       WRITE(6,*) ' and g-functions are not yet implemented '
c$$$       WRITE(6,*)
c$$$       STOP ' sorry '
c$$$      END IF
C
C information on basis sets
C
c$$$      IF (.NOT.LCART) THEN
c$$$       WRITE(6,*)
c$$$       IBAS=0
c$$$       ISHL=0
c$$$       DO IAT=1,NATOM
c$$$        DO IISHL=1,NSH(IAT)
c$$$         ISHL=ISHL+1
c$$$         L=IL(ISHL)
c$$$         DO I=1,L+L+1
c$$$          IBAS=IBAS+1
c$$$          WRITE(6,9226) IBAS,CBSTR(IBAS),NZENTR(ISHL)
c$$$         END DO
c$$$        END DO
c$$$        END DO
c$$$       END IF
 
 9226 FORMAT(' AO No ',I5,' (',A3,') is on atom ',I4)
      WRITE(6,*)
 
C
C prepare tables for the recursions
      WRITE(6,*)
      WRITE(6,*) ' PREPARING TABLES  '
      CALL PREPA
C
C normalizing all primitives
C we suppose that contractions ARE already normalized to 1
C
      WRITE(6,*) ' NORMALIZING PRIMITIVES  '
      WRITE(6,*)
      DO ISHLA=1,NSHLT
       LA=IL(ISHLA)
C the number and start address for the primitives
       NPRIMA=NPRIM(ISHLA)
       NPSTA =NPST (ISHLA)
       IF (NPRIMA.EQ.1) THEN
        EXA=EXX(NPSTA)
        CALL NORM_PRIM(EXA,LA,XNORM)
        COEFF(NPSTA)=XNORM
       ELSE
        DO IIA=1,NPRIMA
         IA=NPSTA+IIA-1
         EXA=EXX(IA)
         CALL NORM_PRIM(EXA,LA,XNORM)
         COEFF(IA)=COEFF(IA)*XNORM
        END DO
       END IF
      END DO
      CALL TIMING('PREP')
 
      IF (LMONO.OR.LKINET) THEN
       CALL MONINT
       CALL TIMING('ONE ')
      END IF
 
      IF (LBIEL) THEN
       CALL TWOINT
       CALL TIMING('TWO ')
      END IF
C that were the mono-electronic integrals, now the bi-electronic ones
 
      X=CPTIME(4)
      CALL DATING(PNAME,2)
      STOP
 901  CONTINUE
      WRITE(6,*) ' NO FILE <SYSTEM.ORTHO> FOUND, PLEASE PROVIDE '
      STOP
 902  CONTINUE
      WRITE(6,*) ' NO FILE <VECTOR> FOUND, PLEASE PROVIDE '
      STOP
 903  CONTINUE
      WRITE(6,*) ' NO FILE <POINTS> FOUND, PLEASE PROVIDE '
      STOP
C
      END
C
      SUBROUTINE RDINP
      INCLUDE 'param.h'                                                 
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
      CHARACTER*4 KEYW(2),STR3
      CHARACTER*6 KEYOPT(10),STR6
      CHARACTER*80 LINE
      DATA KEYW /'*INT','*END'/
      DATA KEYOPT /'NOMONO','THRESH','FORMAT','NOBIEL','CARTES',
     -             'MULINT','KINONL','XXXXXX','XXXXXX','XXXXXX'/
C
C DEFAULTS
      THRESH=1.D-12
      LFORM=.FALSE.
      LMONO=.TRUE.
      LBIEL=.TRUE.
      LCART=.FALSE.
      LKINET=.FALSE.
C
      IOINP=83
      OPEN(UNIT=IOINP,FILE='INPUT.INT',ERR=2217,FORM='FORMATTED',
     -     STATUS='OLD')
C first, look for keyword '*INT'
 1    CONTINUE
      READ(IOINP,'(A4)',END=22,ERR=921) STR3
      IF (STR3.EQ.KEYW(1)) THEN
 11    CONTINUE
       READ(IOINP,'(A6)',END=23,ERR=921) STR6
C       WRITE(6,*) ' READ STRING: |',STR6,'|'
C keyword *END
       IF (STR6(1:4).EQ.KEYW(2)) GO TO 2
       IF (STR6.EQ.KEYOPT(1)) THEN
        LMONO=.FALSE.
       ELSE IF (STR6.EQ.KEYOPT(2)) THEN
        READ(IOINP,*) THRESH
       ELSE IF (STR6.EQ.KEYOPT(3)) THEN
        LFORM=.TRUE.
       ELSE IF (STR6.EQ.KEYOPT(4)) THEN
        LBIEL=.FALSE.
       ELSE IF (STR6.EQ.KEYOPT(5)) THEN
        LCART=.TRUE.
       ELSE IF (STR6.EQ.KEYOPT(6)) THEN
        READ(IOINP,*) MULINT
        WRITE(6,*) '  option MULINT without action, beyond '
     -        ,'dipole not yet implemented '
       ELSE IF (STR6.EQ.KEYOPT(7)) THEN
        LKINET=.TRUE.
       ELSE
        WRITE(6,*) ' POSSIBLE OPTIONS IN THE BLOCK *INT ... *END ARE:'
        WRITE(6,*)
        WRITE(6,*) '              ',KEYOPT(1)
        WRITE(6,*) '              ',KEYOPT(2)
        WRITE(6,*) '              ',KEYOPT(3)
        WRITE(6,*) '              ',KEYOPT(4)
        WRITE(6,*) '              ',KEYOPT(5)
        WRITE(6,*) '              ',KEYOPT(6)
        WRITE(6,*) '              ',KEYOPT(7)
        STOP ' CHOOSE CORRECT OPTION '
       END IF
       GO TO 11
      END IF
      GO TO 1
 2    CONTINUE
      CLOSE(IOINP)
      RETURN
 921  CONTINUE
      CLOSE(IOINP)
      STOP ' ERROR IN INPUT'
 23   CONTINUE
      CLOSE(IOINP)
      WRITE(6,*) ' YOU SHOULD TERMINATE THE BLOCK <*INT> BY <*END> '
      RETURN
 22   CONTINUE
      WRITE(6,*) ' NO KEYWORD <*INT>, USING THE DEFAULT VALUES'
      RETURN
 2217 CONTINUE
      WRITE(6,*) ' NO FILE <INPUT.INT>, USING THE DEFAULT VALUES'
      RETURN
      END
C
      SUBROUTINE NORM_PRIM(EXPON,LVAL,FACNOR)
      INCLUDE 'param.h'                                                 
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
      COMMON /POLYTA/ FAC(NLMAX+NLMAX+2),FACC(NLMAX+1)                  
C..   INCLUDE 'common_polyta.h'                                         
C
C here we normalize a primitive Gaussian function
C
C sqrt[
C l-dependent:
C         2^l*(2a)^(l+3/2)/(pi^3/2 (2l-1)!!)
C m-dependent: (not used here)
C          [(l-|m|)!/(l+|m|)!] * (2-delta_m0)
C       ]
C
      XVAL=DBLE(LVAL)
      PIF=FACC(LVAL+1)*(PI**1.5D0)
CD    WRITE(6,*) ' FACC, PI ',FACC(LVAL+1),PI
      PIF=(2.D0**XVAL)/PIF
      CCC=(2.D0*EXPON)**(XVAL+1.5D0)
      FACNOR=SQRT(PIF*CCC)
CD    WRITE(6,*) ' NORM_PRIM: L,EXPON,FACNOR ',LVAL,EXPON,FACNOR
      RETURN
      END
C
      SUBROUTINE RMCMUR(R000,F,RPC,NMAX)
      INCLUDE 'param.h'                                                 
 
      DIMENSION RPC(3),F(NMAX+1),R000(*)
 
      DIMENSION RU(2,0:12),RV(2,0:12),RT(2,0:12)
C
C nmax = l_a + l_b ; also von 0 .. 8 ; l_max=4
C
C in this routine we perform the recursion
C   from R(000;klm) to R(tuv;000)
C
C the result is communicated via the common block /INTEG/
C NUMHER is the number of hermitian integrals produced
C
C no need for recursion for <s|H|s>
      IF (NMAX.EQ.0) THEN
       R000(1)=F(1)
       RETURN
      END IF
 
      INDX=0
      IIT=1
      IIU=1
      IIV=1
      IIT2=2
      IIU2=2
      IIV2=2
      XPC=RPC(1)
      YPC=RPC(2)
      ZPC=RPC(3)
      DO 2 I=0,NMAX
       RR=F(I+1)
       RT(1,I)=RR
       RT(2,I)=0.D0
       RU(1,I)=RR
       RU(2,I)=0.D0
       RV(1,I)=RR
       RV(2,I)=0.D0
CD     WRITE (6,*) 'I, RV(1,I) = F(I+1) ', I,RV(1,I)
 2    CONTINUE
      DO 11 IT=0,NMAX
       DO 111 I=0,NMAX
         RU(IIU,I)=RT(IIT,I)
 111   CONTINUE
       DO 12 IU=0,NMAX-IT
        DO 121 I=0,NMAX
         RV(IIV,I)=RU(IIU,I)
 121    CONTINUE
        DO 13 IV=0,NMAX-IT-IU
         INDX=INDX+1
C we fill R000 as result of the operation
         R000(INDX)=RV(IIV,0)
         DO 131 K=0,NMAX-IT-IU-IV-1
          RV(IIV2,K)=RV(IIV2,K+1)*DBLE(IV)+ZPC*RV(IIV,K+1)
 131     CONTINUE
         IDUM=IIV
         IIV=IIV2
         IIV2=IDUM
 13     CONTINUE
        DO 122 K=0,NMAX-IT-IU-1
         RU(IIU2,K)=RU(IIU2,K+1)*DBLE(IU)+YPC*RU(IIU,K+1)
 122    CONTINUE
        IDUM=IIU
        IIU=IIU2
        IIU2=IDUM
 12    CONTINUE
       DO 112 K=0,NMAX-IT-1
        RT(IIT2,K)=RT(IIT2,K+1)*DBLE(IT)+XPC*RT(IIT,K+1)
 112   CONTINUE
       IDUM=IIT
       IIT=IIT2
       IIT2=IDUM
 11   CONTINUE
C
CD    WRITE(6,*)
CD    WRITE(6,*) ' we produced ',INDX,' hermitian intermediates'
CD    WRITE(6,*) ' NMAX = ',NMAX
C
      RETURN
      END
C
C we prepare arrays, the Gamma table, the transformation Cart -> spherical
C
      SUBROUTINE PREPA
      INCLUDE 'param.h'                                                 
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
 
      CALL FILSPH
      CALL FFCAL
C Table for the Gamma function
 
      MXABCD=LMAX*4+1
      CALL GAMTAB(MXABCD)
 
      RETURN
      END
C
      SUBROUTINE FFCAL
      INCLUDE 'param.h'                                                 
      COMMON /POLYTA/ FAC(NLMAX+NLMAX+2),FACC(NLMAX+1)                  
C..   INCLUDE 'common_polyta.h'                                         
C
C the (2l-1)!! and the (l+|m|)!
C first the (2l-1)!!
C
      FAC(1)=1.D0
      FAC(2)=1.D0
      DO I=3,NLMAX+NLMAX+2
       XMLT=DBLE(I-1)
       FAC(I)=XMLT*FAC(I-1)
      END DO
C
      L=1
      FACC(1)=1.D0
CD    WRITE(6,*) ' L, FACC(L) ',L-1,FACC(L)
      DO L=2,NLMAX+1
       FACC(L)=FACC(L-1)*DBLE(L+L-3)
CD     WRITE(6,*) ' L, FACC(L) ',L-1,FACC(L)
      END DO
C
CD    DO L=1,NLMAX
CD     WRITE(6,*) ' L, FAC(L), FACC(L) ',L-1,FAC(L),FACC(L)
CD    END DO
CD    DO L=1+NLMAX,NLMAX+NLMAX
CD     WRITE(6,*) ' L, FAC(L)          ',L-1,FAC(L)
CD    END DO
C
      RETURN
      END
C
      SUBROUTINE FILSPH
      INCLUDE 'param.h'                                                 
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      COMMON /POLYTA/ FAC(NLMAX+NLMAX+2),FACC(NLMAX+1)                  
C..   INCLUDE 'common_polyta.h'                                         
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
C
      INCR=1
      NCART(0)=1
      DO L=1,LMAX
       INCR=INCR+1
       NCART(L)=NCART(L-1)+INCR
      END DO
 
      do l=0,lmax
       do m=1,(l+1)*(l+2)/2
        carmul(l,m)=1.d0
       end do
      end do
C proper normalization factors for cartesian d functions
      carmul(2,2)=sqrt(3.d0)
      carmul(2,4)=sqrt(3.d0)
      carmul(2,5)=sqrt(3.d0)
C proper normalization factors for cartesian f functions
      carmul(3,2)=sqrt(5.d0)
      carmul(3,3)=sqrt(5.d0)
      carmul(3,5)=sqrt(5.d0)
      carmul(3,6)=sqrt(15.d0)
      carmul(3,7)=sqrt(5.d0)
      carmul(3,8)=sqrt(5.d0)
      carmul(3,9)=sqrt(5.d0)
C
      DO I=0,LMAX
       DO J=1,NCART(I)
        DO K=1,I+I+1
         SPHCAR(I,J,K)=0.D0
        END DO
       END DO
      END DO
 
      IF (LMAX.GE.0) THEN
C s
C cartesian:
C    1
C    1
C   000
       SPHCAR(0,1,1)=1.D0
       IF (LMAX.GE.1) THEN
C p  spherical harmonics:
C m   -1   0   1
C      x   z   y
C cartesian:
C    1   2   3
C    z   y   x
C   001 010 100
C
C we put them in the order x, y, z
C for both cases
C
        SPHCAR(1,1,3)=1.D0
        SPHCAR(1,2,2)=1.D0
        SPHCAR(1,3,1)=1.D0
        IF (LMAX.GE.2) THEN
C d spherical harmonics:
C m   -2  -1   0      1     2
C     xy, yz, 3z2-r2, xz, x2-y2
C
C cartesian:
C     1   2   3   4   5   6
C    zz  yz  yy  xz  xy  xx
C   002 011 020 101 110 200
C
         SPHCAR(2,5,1)= SQRT(3.D0)
         SPHCAR(2,2,2)= SQRT(3.D0)
c$$$         SPHCAR(2,1,3)= SQRT(3.D0)*2.D0/SQRT(12.D0)
c$$$         SPHCAR(2,3,3)=-SQRT(3.D0)/SQRT(12.D0)
c$$$         SPHCAR(2,6,3)=-SQRT(3.D0)/SQRT(12.D0)
         SPHCAR(2,1,3)= 1.D0
         SPHCAR(2,3,3)=-0.5D0
         SPHCAR(2,6,3)=-0.5D0
 
         SPHCAR(2,4,4)= SQRT(3.D0)
         SPHCAR(2,6,5)= SQRT(3.D0)/2.D0
         SPHCAR(2,3,5)=-SQRT(3.D0)/2.D0
         IF (LMAX.GE.3) THEN
C
C f
C m  -3    -2       -1         0       1              2        3
C    33s   32s      31s       30       31c           32c      33c
C 3xxy-yyy xyz 4yzz-xxy-yyy 5zzz-3zrr 4xzz-xxx-xyy z(xx-yy) xxx-3xyy
C
C     1   2   3   4   5   6  7   8   9  10
C   zzz yzz yyz yyy xzz xyz xyy xxz xxy xxx
C   003 012 021 030 102 111 120 201 210 300
C
          SPHCAR(3, 9,1) = 3.D0*SQRT(10.D0)/4.D0
          SPHCAR(3, 4,1) =     -SQRT(10.D0)/4.D0
          SPHCAR(3, 6,2) =  SQRT(15.D0)
          SPHCAR(3, 2,3) =  SQRT(6.D0)
          SPHCAR(3, 9,3) = -SQRT(6.D0)/4.D0
          SPHCAR(3, 4,3) = -SQRT(6.D0)/4.D0
C this is the convention of Helgaker et al, table 6.3
C but leading to negative monoelectronic integrals
C with respect to Dalton
c$$$          SPHCAR(3, 1,4) =  1.D0
c$$$          SPHCAR(3, 8,4) = -1.5D0
c$$$          SPHCAR(3, 3,4) = -1.5D0
C this is the convention used as well in Dalton
          SPHCAR(3, 1,4) =  -1.D0
          SPHCAR(3, 8,4) =   1.5D0
          SPHCAR(3, 3,4) =   1.5D0
          SPHCAR(3, 5,5) = SQRT(6.D0)
          SPHCAR(3,10,5) =-SQRT(6.D0)/4.D0
          SPHCAR(3, 7,5) =-SQRT(6.D0)/4.D0
          SPHCAR(3, 8,6) = SQRT(15.D0)/2.D0
          SPHCAR(3, 3,6) =-SQRT(15.D0)/2.D0
C this is the convention of Helgaker et al, table 6.3
c$$$          SPHCAR(3, 7,7) = -3.D0*SQRT(10.D0)/4.D0
c$$$          SPHCAR(3,10,7) =  SQRT(10.D0)/4.D0
C this is the convention used as well in Dalton
          SPHCAR(3, 7,7) =  3.D0*SQRT(10.D0)/4.D0
          SPHCAR(3,10,7) = -SQRT(10.D0)/4.D0
          IF (LMAX.GE.4) THEN
C
C g functions
C m  -4 -3 -2 -1 0 1 2 3 4
C
C cartesian
C   1    2    3    4    5    6   7    8    9   10   11   12   13   14   15
C zzzz yzzz yyzz yyyz yyyy xzzz xyzz xyyz xyyy xxzz xxyz xxyy xxxz xxxy xxxx
C 004  013  022  031  040  103  112  121  130  202  211  220  301  310  400
C
C we took the Dalton table which is printed with .PRINT with 9
           SPHCAR(4,14,1)= SQRT(35.D0)/2.D0
           SPHCAR(4,9, 1)=-SQRT(35.D0)/2.D0
 
           SPHCAR(4, 4,2)=-1.D0*SQRT(35.D0/8.D0)
           SPHCAR(4,11,2)= 3.D0*SQRT(35.D0/8.D0)
 
           SPHCAR(4,14,3)=     -1.D0*sqrt(5.D0)/2.D0
           SPHCAR(4, 9,3)=     -1.D0*sqrt(5.D0)/2.D0
           SPHCAR(4, 7,3)=      6.D0*sqrt(5.D0)/2.D0
 
           SPHCAR(4,11,4)=      -1.D0/sqrt(8.d0/45.d0)
           SPHCAR(4, 4,4)=      -1.D0/sqrt(8.d0/45.d0)
           SPHCAR(4, 2,4)=      4.D0/3.D0/sqrt(8.d0/45.d0)
 
C the same as above, the sign is opposite to that of Helgaker, table 6.3
           SPHCAR(4,15,5)=      -1.D0/8.D0*3.D0
           SPHCAR(4,12,5)=      -1.D0/4.D0*3.D0
           SPHCAR(4,10,5)=       1.D0*3.D0
           SPHCAR(4, 5,5)=      -1.D0/8.D0*3.D0
           SPHCAR(4, 3,5)=       1.D0*3.D0
           SPHCAR(4, 1,5)=      -1.D0/3.D0*3.D0
 
           SPHCAR(4,13,6)=      -1.D0/sqrt(8.d0/45.d0)
           SPHCAR(4, 8,6)=      -1.D0/sqrt(8.d0/45.d0)
           SPHCAR(4, 6,6)=       4.D0/3.D0/sqrt(8.d0/45.d0)
 
           SPHCAR(4,15,7)=      -1.D0*sqrt(5.D0)/4
           SPHCAR(4,10,7)=       6.D0*sqrt(5.D0)/4
           SPHCAR(4, 5,7)=       1.D0*sqrt(5.D0)/4
           SPHCAR(4, 3,7)=      -6.D0*sqrt(5.D0)/4
C the same as above, the sign is opposite to that of Helgaker, table 6.3
           SPHCAR(4,13,8)=      -1.D0/3.D0*sqrt(35.D0)*3.D0/4.D0*sqrt(2.
     $D0)
C folded 1 (fixf $Revision: 1.3 $) 
           SPHCAR(4, 8,8)=       1.D0*sqrt(35.D0)*3.D0/4.D0*sqrt(2.D0)
 
           SPHCAR(4,15,9)=       -1.D0/6.d0*sqrt(35.D0)*3.D0/4.D0
           SPHCAR(4,12,9)=        1.D0*sqrt(35.D0)*3.D0/4.D0
           SPHCAR(4, 5,9)=       -1.D0/6.d0*sqrt(35.D0)*3.D0/4.D0
 
          END IF
         END IF
        END IF
       END IF
      END IF
 
      DO I=0,LMAX
       DO K=1,I+I+1
        DO J=1,NCART(I)
         IF (ABS(SPHCAR(I,J,K)).GT.1.D-9) THEN
          WRITE(6,9901) I,J,K,SPHCAR(I,J,K)
 9901     FORMAT('     L= ',I3,' CART = ',I3,' SPH= ',I3,' : ',F20.12)
         end if
        END DO
       END DO
       WRITE(6,*)
      END DO
 
      write(6,*)
 
      RETURN
      END
C
      SUBROUTINE SPHHAM(XINTC,XINTS)
      INCLUDE 'param.h'                                                 
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      DIMENSION XINTC(NCMAX,NCMAX)
      DIMENSION XINTS(NSPMAX,NSPMAX)
C
C transformation from cartesian coordinates on spherical harmonics
C
      DO I=1,LA+LA+1
       DO J=1,LB+LB+1
        XINTS(I,J)=0.D0
       END DO
      END DO
C
      DO MA=1,LA+LA+1
       DO I=1,NCART(LA)
        XMA=SPHCAR(LA,I,MA)
        IF (ABS(XMA).GT.1.D-12) THEN
         DO MB=1,LB+LB+1
          DO J=1,NCART(LB)
           XMB=SPHCAR(LB,J,MB)
           XINTS(MA,MB) =XINTS(MA,MB)+XMA*XMB*XINTC(I,J)
          END DO
         END DO
        END IF
       END DO
      END DO
C
      RETURN
      END
 
      SUBROUTINE MONINT
      INCLUDE 'param.h'                                                 
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
C integrals in hermitian Gaussians                                      
      COMMON /INTEGH/ R000(LHMAX)                                       
C..   INCLUDE 'common_integ.h'                                          
C static quantities for generating the tables and the interpolation form
C of the incomplete gamma function                                      
      PARAMETER (NTABG=20000)                                           
      COMMON/TABAUX/TABGAM(NTABG),NVALS,NSEGM,NSEGM2,NSEGM3,NSEGM4      
     -     ,NSEGM5,NSEGM6                                               
C..   INCLUDE 'common_tabaux.h'                                         
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      COMMON /POLYTA/ FAC(NLMAX+NLMAX+2),FACC(NLMAX+1)                  
C..   INCLUDE 'common_polyta.h'                                         
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
C
C integrals in cartesians                                               
      COMMON /INTEGC/ POTC(NCMAX,NCMAX),OVERL(NCMAX,NCMAX)              
     $     ,EKIN(NCMAX,NCMAX),DIPXC(NCMAX,NCMAX),DIPYC(NCMAX,NCMAX)     
     $     ,DIPZC(NCMAX,NCMAX),NUMCAR                                   
C integrals in spherical harmonics                                      
      COMMON /INTEGS/ POTS(NSPMAX,NSPMAX),OVERS(NSPMAX,NSPMAX)          
     $     ,EKINS(NSPMAX,NSPMAX),DIPXS(NSPMAX,NSPMAX),DIPYS(NSPMAX      
     $     ,NSPMAX),DIPZS(NSPMAX,NSPMAX)                                
C final one-electron integrals                                          
      COMMON /INTU/ VPOT(NBASM,NBASM),TKIN(NBASM,NBASM),S(NBASM,NBASM)  
     $     ,DIPX(NBASM,NBASM),DIPY(NBASM,NBASM),DIPZ(NBASM,NBASM)       
C..   INCLUDE 'dimens_intu.h'                                           
 
      DIMENSION XPA(3),XPB(3),XPC(3)
      DIMENSION SS(3,-1:NLMAX+2,-1:NLMAX)
      DIMENSION DD(3,0:NLMAX,0:NLMAX)
 
      do i=1,3
       do j=-1,lmax+2
        do k=-1,lmax
         ss(i,j,k)=0.d0
        end do
       end do
      end do
 
CR    LOGICAL LERR
CR    DIMENSION SREAD(NBASM,NBASM),TREAD(NBASM,NBASM),VREAD(NBASM,NBASM)
 
C we read the DALTON integrals for comparison
C
CR    DO I=1,NBAS
CR     DO J=1,NBAS
CR      SREAD(I,J)=0.D0
CR      TREAD(I,J)=0.D0
CR      VREAD(I,J)=0.D0
CR     END DO
CR    END DO
 
CR    OPEN(UNIT=15,File='OVERLAP',STATUS='OLD',FORM='FORMATTED')
 107  CONTINUE
CR    READ(15,*,IOSTAT=KK) I,J,XDUM
CR    IF (KK.NE.0) GO TO 108
CR    SREAD(I,J)=XDUM
CR    SREAD(J,I)=XDUM
CR    GO TO 107
 108  CONTINUE
CR    CLOSE(15)
 
CR    OPEN(UNIT=15,File='KINETIC',STATUS='OLD',FORM='FORMATTED')
 207  CONTINUE
CR    READ(15,*,IOSTAT=KK) I,J,XDUM
CR    IF (KK.NE.0) GO TO 208
CR    TREAD(I,J)=XDUM
CR    TREAD(J,I)=XDUM
CR    GO TO 207
 208  CONTINUE
CR    CLOSE(15)
 
CR    OPEN(UNIT=15,File='HAMILTO',STATUS='OLD',FORM='FORMATTED')
 307  CONTINUE
CR    READ(15,*,IOSTAT=KK) I,J,XDUM
CR    IF (KK.NE.0) GO TO 308
CR    VREAD(I,J)=XDUM
CR    VREAD(J,I)=XDUM
CR    GO TO 307
 308  CONTINUE
CR    CLOSE(15)
C
CR    DO I=1,NBAS
CR     DO J=1,NBAS
CR      VREAD(I,J)=VREAD(I,J)-TREAD(I,J)
CR     END DO
CR    END DO
C
CR    DO I=1,NBAS
CR     DO J=I,NBAS
CR      WRITE(6,9777) I,J,SREAD(I,J),TREAD(I,J),VREAD(I,J)
CR     END DO
CR    END DO
C
      DO I=1,NBAS
       DO J=1,NBAS
        VPOT(I,J)=0.D0
        TKIN(I,J)=0.D0
        S   (I,J)=0.D0
        DIPX(I,J)=0.D0
        DIPY(I,J)=0.D0
        DIPZ(I,J)=0.D0
       END DO
      END DO
 
      TWOPI=2.D0*PI
C
C the loop over atomic shells
C
C the real work
 
      ISTA=1
      DO ISHLA=1,NSHLT
CD     WRITE(6,*) ' TREATING SHELL A = ',ISHLA
C the atom
       IATA=NZENTR(ISHLA)
C get the position
       DO I=1,3
        XA(I)=POS(I,IATA)
       END DO
C the angular momentum
       LA=IL(ISHLA)
C the number and start address for the primitives
       NPRIMA=NPRIM(ISHLA)
       NPSTA =NPST (ISHLA)
 
       ISTB=ISTA
       DO ISHLB=ISHLA,NSHLT
CD      WRITE(6,*) '  TREATING SHELL B = ',ISHLB
        IATB=NZENTR(ISHLB)
C get the angular momentum, the position, the primitives and the
C coefficients
        DO I=1,3
         XB(I)=POS(I,IATB)
         XAB(I)=XA(I)-XB(I)
        END DO
        LB=IL(ISHLB)
C the number and start address for the primitives
        NPRIMB=NPRIM(ISHLB)
        NPSTB =NPST (ISHLB)
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Overlap and Kinetic energy integrals
C
CD      WRITE(6,*)
CD      WRITE(6,*) ' ------------------------------------------------'
CD      WRITE(6,*)
C
C
C set cartesian integrals to zero
C
        DO I=1,NCART(LA)
         DO J=1,NCART(LB)
          OVERL(I,J)=0.D0
          EKIN (I,J)=0.D0
          POTC (I,J)=0.D0
          DIPXC(I,J)=0.D0
          DIPYC(I,J)=0.D0
          DIPZC(I,J)=0.D0
         END DO
        END DO
 
        LAB=LA+LB
C
        IA=NPSTA
        DO IIA=1,NPRIMA
         EXA=EXX(IA)
         CCA=COEFF(IA)
         IA=IA+1
 
         IB=NPSTB
         DO IIB=1,NPRIMB
          EXB=EXX(IB)
          COEFFP=CCA*COEFF(IB)
          IB=IB+1
C
          P=EXA+EXB
          XMU=EXA*EXB/P
          RAB=0.D0
          DO I=1,3
           XP(I)=(EXA*XA(I)+EXB*XB(I))/P
           XPA(I)=-EXB/P*XAB(I)
           XPB(I)= EXA/P*XAB(I)
           RAB=RAB+XAB(I)*XAB(I)
           SS(I,0,0)=SQRT(PI/P)
           SS(I,1,0)=XPA(I)*SS(I,0,0)
           SS(I,0,1)=XPB(I)*SS(I,0,0)
          END DO
          EMU=EXP(-XMU*RAB)
          COEFFP=COEFFP*EMU
          DO I=1,3
           DO L=2,LA+LB+2
            JJ=0
            IF (L.LE.LA+2) THEN
             SS(I,L,0)=XPA(I)*SS(I,L-1,0)+0.5D0/P*DBLE(L-1)*SS(I,L-2,0)
            END IF
            IF (L.LE.LB) THEN
             SS(I,0,L)=XPB(I)*SS(I,0,L-1)+0.5D0/P*DBLE(L-1)*SS(I,0,L-2)
            END IF
C horizontal recurrence - we create all intermediates i+j=l
            DO JJ=1,L-1
             II=L-JJ
             IF (II.LE.LA+2.AND.JJ.LE.LB) THEN
              YYY=DBLE(II)*SS(I,II-1,JJ-1)
              IF (JJ.GE.2) THEN
               YYY=YYY+DBLE(JJ-1)*SS(I,II,JJ-2)
              END IF
              SS(I,II,JJ)=XPB(I)*SS(I,II,JJ-1)+0.5D0/P*YYY
             END IF
            END DO
           END DO
          END DO
C
          DO I=1,3
           DO L2=0,LB
            DO L1=0,LA
             DD(I,L1,L2)=2.D0*EXA*(2.D0*EXA*SS(I,L1+2,L2)-DBLE(2*L1+1)
     $            *SS(I,L1,L2))
             IF (L1.GE.2) THEN
              DD(I,L1,L2)=DD(I,L1,L2) + DBLE(L1*(L1-1))*SS(I,L1-2,L2)
             END IF
            END DO
           END DO
          END DO
C
          INDX1=0
          DO I=0,LA
           DO J=0,LA-I
            INDX1=INDX1+1
            INDX2=0
            DO II=0,LB
             DIPXH=XP(1)*SS(1,I,II)+0.5D0/P*(DBLE(I)*SS(1,I-1,II)
     $            +DBLE(II-1)*SS(1,I,II-1))
             DO JJ=0,LB-II
              INDX2=INDX2+1
              K =LA-I-J
              KK=LB-II-JJ
 
              OVERL(INDX1,INDX2)=OVERL(INDX1,INDX2)+ COEFFP*SS(1,I,II)
     $             *SS(2,J,JJ)*SS(3,K,KK)
 
              DIPXH=XP(1)*SS(1,I,II)+0.5D0/P*(DBLE(I)*SS(1,I-1,II)
     $            +DBLE(II)*SS(1,I,II-1))
              DIPYH=XP(2)*SS(2,J,JJ)+0.5D0/P*(DBLE(J)*SS(2,J-1,JJ)
     $             +DBLE(JJ)*SS(2,J,JJ-1))
              DIPZH=XP(3)*SS(3,K,KK)+0.5D0/P*(DBLE(K)*SS(3,K-1,KK)
     $             +DBLE(KK)*SS(3,K,KK-1))
 
              DIPXC(INDX1,INDX2)=DIPXC(INDX1,INDX2) +
     -             COEFFP*DIPXH*SS(2,J,JJ)*SS(3,K,KK)
              DIPYC(INDX1,INDX2)=DIPYC(INDX1,INDX2) +
     -             COEFFP*SS(1,I,II)*DIPYH*SS(3,K,KK)
              DIPZC(INDX1,INDX2)=DIPZC(INDX1,INDX2) +
     -             COEFFP*SS(1,I,II)*SS(2,J,JJ)*DIPZH
 
              EKIN(INDX1,INDX2)=EKIN(INDX1,INDX2)  -COEFFP*0.5D0*(DD(1,I
     $             ,II)*SS(2,J,JJ)*SS(3,K,KK)+SS(1,I,II)*DD(2,J,JJ)*SS(3
     $             ,K,KK)+SS(1,I,II)*SS(2,J,JJ)*DD(3,K,KK))
CD            WRITE(6,*) ' INDX1, INDX2, S (',I,J,K,') (',II,JJ,KK,') '
CD   $             ,OVERL(INDX1,INDX2),EKIN(INDX1,INDX2)
             END DO
            END DO
           END DO
          END DO
C
C loop over atomic positions
C we do not need el-n integrals in the case LKINET
          IF (.NOT.LKINET) THEN
          PIKALL=COEFFP*TWOPI/P
          DO IAT=1,NATOM
 
           IF (NZ(IAT).NE.0) THEN
            DO I=1,3
             XC(I)=POS(I,IAT)
            END DO
CD          WRITE(6,*) ' POSITION OF NUCLEUS: ',(XC(J),J=1,3)
C
C adapt quantities to common block PQCON
C
            RRPC=0.D0
            DO I=1,3
             XPC(I)=XP(I)-XC(I)
             RRPC=RRPC+XPC(I)*XPC(I)
            END DO
            ARG=P*RRPC
C
CD    WRITE(6,*)
CD    WRITE(6,*) ' entering BUXG for the incomplete Gamma function'
CD    WRITE(6,*) ' P * RRPC = ',PRRPC
CD    WRITE(6,*) ' Prefactors = ',EMU,COEFFP,XMU
CD    WRITE(6,*)
 
 
            NMAX=LAB+1
            DO L=1,NMAX
             F(L)=0.D0
            END DO
            CALL GAMF(F,ARG,NMAX)
C
            AMUL=1.D0
            I=1
CD    WRITE(6,*) ' I, F(I) = ',I,F(I)
            DO I=2,LAB+1
             AMUL=-(P+P)*AMUL
             F(I)=F(I)*AMUL
CD     WRITE(6,*) ' I, F(I) = ',I,F(I)
CD     WRITE(6,*)
C
            END DO
C
C that were the preparations
C
            CALL RMCMUR(R000,F,XPC,LAB)
            PIK=DBLE(NZ(IAT))*PIKALL
C we accumulate the e-n interaction for each nucleus in POTC
            CALL PHIREC(POTC,R000,XAB,EXA,EXB,PIK,LA,LB)
C
           END IF
          END DO
C close IF of LKINET
          END IF
C close loops over primitives EXA and EXB
         END DO
        END DO
C
C if cartesian integrals are demanded we have all here
C we come from the recursion with xx, xy, xz, yy, yz, zz and so on
C
        IF (LCART) THEN
         DO IMA=1,NCART(LA)
          IALPH=ISTA+ICARCAR(LA,IMA)-1
          DO IMB=1,NCART(LB)
           IBETA=ISTB+ICARCAR(LB,IMB)-1
           XMULT=CARMUL(LA,IMA)*CARMUL(LB,IMB)
           TKIN(IALPH,IBETA)=EKIN (IMA,IMB)*XMULT
           S   (IALPH,IBETA)=OVERL(IMA,IMB)*XMULT
           VPOT(IALPH,IBETA)=POTC (IMA,IMB)*XMULT
           DIPX(IALPH,IBETA)=DIPXC(IMA,IMB)*XMULT
           DIPY(IALPH,IBETA)=DIPYC(IMA,IMB)*XMULT
           DIPZ(IALPH,IBETA)=DIPZC(IMA,IMB)*XMULT
CD         WRITE(6,*)
CD         WRITE(6,9777) IALPH,IBETA,S(IALPH,IBETA),TKIN(IALPH,IBETA)
CD   $          ,VPOT(IALPH,IBETA)
CD         WRITE(6,*)
          END DO
         END DO
         ISTB=ISTB+NCART(LB)
        ELSE
C
C transformation onto spherical harmonics
C
         CALL SPHHAM(EKIN ,EKINS)
         IF (.NOT.LKINET) THEN
          CALL SPHHAM(OVERL,OVERS)
          CALL SPHHAM(POTC ,POTS )
          CALL SPHHAM(DIPXC,DIPXS)
          CALL SPHHAM(DIPYC,DIPYS)
          CALL SPHHAM(DIPZC,DIPZS)
         END IF
 
         DO IMA=1,1+LA+LA
          IALPH=ISTA+IMA-1
          DO IMB=1,1+LB+LB
           IBETA=ISTB+IMB-1
         IF (.NOT.LKINET) THEN
           TKIN(IALPH,IBETA)=EKINS(IMA,IMB)
           S   (IALPH,IBETA)=OVERS(IMA,IMB)
           VPOT(IALPH,IBETA)=POTS (IMA,IMB)
           DIPX(IALPH,IBETA)=DIPXS(IMA,IMB)
           DIPY(IALPH,IBETA)=DIPYS(IMA,IMB)
           DIPZ(IALPH,IBETA)=DIPZS(IMA,IMB)
         END IF
CD         WRITE(6,*)
CD         WRITE(6,9777) IALPH,IBETA,S(IALPH,IBETA),TKIN(IALPH,IBETA)
CD   $          ,VPOT(IALPH,IBETA)
 9777      FORMAT(' I=',I3,'  J=',I3,'  S=',F18.11,'  T=',F18.11,'  V='
     $          ,F18.11)
CD         WRITE(6,*)
          END DO
         END DO
         ISTB=ISTB+LB+LB+1
        END IF
C
C close loops over shells
C
       END DO
       IF (LCART) THEN
        ISTA=ISTA+NCART(LA)
       ELSE
        ISTA=ISTA+LA+LA+1
       END IF
c$$$       CALL TIMING('SH A')
      END DO
C
CR    WRITE(6,*)
CR    WRITE(6,*) ' DIFFERENCE TO READ VALUES: ( > 1.D-8 ) '
CR    WRITE(6,*)
CR    LERR=.FALSE.
CR    DO I=1,NBAS
CR     DO J=I,NBAS
CR      IF (ABS(VREAD(I,J)-VPOT(I,J)).GT.1.D-8) THEN
CR       WRITE(6,9877) I,J,VPOT(I,J)-VREAD(I,J),VPOT(I,J),VREAD(I,J)
CR       LERR=.TRUE.
CR      END IF
CR      IF (I.EQ.50.AND.J.EQ.51) WRITE(6,*) ' VVV ',I,J,VPOT(I,J)
CR      IF (I.EQ.50.AND.J.EQ.51) WRITE(6,*) ' TTT ',I,J,TKIN(I,J)
CR      IF (ABS(SREAD(I,J)-S(I,J)).GT.1.D-8) THEN
CR       WRITE(6,9878) I,J,S(I,J)-SREAD(I,J),S(I,J),SREAD(I,J)
CR       LERR=.TRUE.
CR      END IF
CR      IF (ABS(TREAD(I,J)-TKIN(I,J)).GT.1.D-8) THEN
CR       WRITE(6,9879) I,J,TKIN(I,J)-TREAD(I,J),TKIN(I,J),TREAD(I,J)
CR       LERR=.TRUE.
CR      END IF
CR     END DO
CR    END DO
CR    IF (.NOT.LERR) THEN
CR     WRITE(6,*) ' NO DIFFERENCE TO REPORT !!! '
CR    END IF
 9877 FORMAT('V I',I3,' J',I3,' DIFF=',F18.11,' C=',F18.11,' R='
     $     ,F18.11)
 9878 FORMAT('S I',I3,' J',I3,' DIFF=',F18.11,' C=',F18.11,' R='
     $     ,F18.11)
 9879 FORMAT('T I',I3,' J',I3,' DIFF=',F18.11,' C=',F18.11,' R='
     $     ,F18.11)
 
C write the integrals to file
C--SKIPCR
 
      WRITE(6,*) ' WRITING FILE <KINETIC> '
      OPEN(UNIT=16,FILE='KINETIC',FORM='FORMATTED',STATUS='UNKNOWN')
      DO I=1,NBAS
       DO J=I,NBAS
          IF (ABS(TKIN(I,J)).GT.THRESH)
     -       WRITE(16,'(2I5,F22.12)') J,I,TKIN(I,J)
       END DO
      END DO
      CLOSE(16)
 
      IF (.NOT.LKINET) THEN
       WRITE(6,*) ' WRITING FILE <HAMILTO> '
       WRITE(6,*) ' WRITING FILE <OVERLAP> '
       WRITE(6,*) ' WRITING FILE <DIPOL_X> '
       WRITE(6,*) ' WRITING FILE <DIPOL_Y> '
       WRITE(6,*) ' WRITING FILE <DIPOL_Z> '
       OPEN(UNIT=14,FILE='OVERLAP',FORM='FORMATTED',STATUS='UNKNOWN')
       OPEN(UNIT=15,FILE='HAMILTO',FORM='FORMATTED',STATUS='UNKNOWN')
       OPEN(UNIT=17,FILE='DIPOL_X',FORM='FORMATTED',STATUS='UNKNOWN')
       OPEN(UNIT=18,FILE='DIPOL_Y',FORM='FORMATTED',STATUS='UNKNOWN')
       OPEN(UNIT=19,FILE='DIPOL_Z',FORM='FORMATTED',STATUS='UNKNOWN')
       DO I=1,NBAS
        DO J=I,NBAS
         IF (ABS(S(I,J)).GT.THRESH)
     -        WRITE(14,'(2I5,F22.12)') J,I,S(I,J)
         VDUM=VPOT(I,J)+TKIN(I,J)
         IF (ABS(VDUM).GT.THRESH)
     -        WRITE(15,'(2I5,F22.12)') J,I,VDUM
         IF (ABS(TKIN(I,J)).GT.THRESH)
     -        WRITE(16,'(2I5,F22.12)') J,I,TKIN(I,J)
         IF (ABS(DIPX(I,J)).GT.THRESH)
     -        WRITE(17,'(2I5,F22.12)') J,I,DIPX(I,J)
         IF (ABS(DIPY(I,J)).GT.THRESH)
     -        WRITE(18,'(2I5,F22.12)') J,I,DIPY(I,J)
         IF (ABS(DIPZ(I,J)).GT.THRESH)
     -        WRITE(19,'(2I5,F22.12)') J,I,DIPZ(I,J)
        END DO
       END DO
       CLOSE(14)
       CLOSE(15)
       CLOSE(17)
       CLOSE(18)
       CLOSE(19)
      END IF
 
C//ENDCR
 
      RETURN
      END
 
      SUBROUTINE TWOINT
      INCLUDE 'param.h'                                                 
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
C integrals in hermitian Gaussians                                      
      COMMON /INTEGH/ R000(LHMAX)                                       
C..   INCLUDE 'common_integ.h'                                          
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      COMMON /POLYTA/ FAC(NLMAX+NLMAX+2),FACC(NLMAX+1)                  
C..   INCLUDE 'common_polyta.h'                                         
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
C
      DIMENSION XPA(3),XPB(3)
      DIMENSION SS(3,0:NLMAX+2,0:NLMAX)
      DIMENSION DD(3,0:NLMAX,0:NLMAX)
      LOGICAL LTAKE
      LOGICAL LOCANC
C
C with f functions we have at most 2401 possibilities in spherical
C harmonics, but 10000 for cartesian functions for the (ff|ff)
C integrals
C
      DIMENSION HHHBUF(NCMAX*NCMAX*NCMAX*NCMAX)
      IF (LFORM) THEN
       OPEN(UNIT=14,FILE='AOTWOINT_FORMATTED',STATUS='UNKNOWN',FORM
     $      ='FORMATTED')
      ELSE
       CALL WOPEN(4,1)
      END IF
 
      XDUMP=0.D0
      XCALC=0.D0
      XSSSS=0.D0
      XSSCD=0.D0
      XABSS=0.D0
      XABCD=0.D0
 
      INDXA=0
      DO ISHLA=1,NSHLT
CD     WRITE(6,*) ' TREATING SHELL A = ',ISHLA
C the atom
       IATA=NZENTR(ISHLA)
       DO I=1,3
        XA(I)=POS(I,IATA)
       END DO
       LA=IL(ISHLA)
       NPRIMA=NPRIM(ISHLA)
       NPSTA =NPST (ISHLA)
       IF (LCART) THEN
        NINDA=NCART(LA)
       ELSE
        NINDA=LA+LA+1
       END IF
 
       INDXB=INDXA
       DO ISHLB=ISHLA,NSHLT
CD     WRITE(6,*) '  TREATING SHELL B = ',ISHLB
C the atom
        IATB=NZENTR(ISHLB)
        DO I=1,3
         XB(I)=POS(I,IATB)
        END DO
        LB=IL(ISHLB)
        NPRIMB=NPRIM(ISHLB)
        NPSTB =NPST (ISHLB)
        IF (LCART) THEN
         NINDB=NCART(LB)
        ELSE
         NINDB=LB+LB+1
        END IF
 
        LAB=LA+LB
 
        INDXC=INDXA
        DO ISHLC=ISHLA,NSHLT
CD     WRITE(6,*) '   TREATING SHELL C = ',ISHLC
C the atom
         IATC=NZENTR(ISHLC)
         DO I=1,3
          XC(I)=POS(I,IATC)
         END DO
         LC=IL(ISHLC)
         NPRIMC=NPRIM(ISHLC)
         NPSTC =NPST (ISHLC)
         IF (LCART) THEN
          NINDC=NCART(LC)
         ELSE
          NINDC=LC+LC+1
         END IF
 
         LABC=LAB+LC
 
         INDXD=INDXC
 
         DO ISHLD=ISHLC,NSHLT
CD        WRITE(6,*) '    TREATING SHELL D = ',ISHLA,ISHLB,ISHLC,ISHLD
CD     WRITE(6,*) ' LA LB LC LD = ',LA,LB,LC,LD
C the atom
          IATD=NZENTR(ISHLD)
          DO I=1,3
           XD(I)=POS(I,IATD)
          END DO
          LD=IL(ISHLD)
          NPRIMD=NPRIM(ISHLD)
          NPSTD =NPST (ISHLD)
          IF (LCART) THEN
           NINDD=NCART(LD)
          ELSE
           NINDD=LD+LD+1
          END IF
CD     WRITE(6,*) ' IATA IATB IATC IATD = ',IATA,IATB,IATC,IATD
 
          LABCD=LABC+LD
          LCD=LC+LD
 
          IF (LABCD.EQ.0) THEN
C this is the simplest of the integrals
           CALL SSSS(HHH)
           I=INDXA+1
           J=INDXB+1
           K=INDXC+1
           L=INDXD+1
           IF (LOCANC(I,J,K,L)) THEN
            XMULT=FACTOR(I,J,K,L)
            HHH=XMULT*HHH
CD                WRITE(6,9901) I,J,K,L,HHH
            XSSSS=XSSSS+1.D0
            IF (ABS(HHH).GE.THRESH) THEN
             XDUMP=XDUMP+1.D0
             IF (LFORM) THEN
              WRITE(14,'(4I5,E24.16)') I,J,K,L,HHH
 9901         FORMAT(' wrote integral ',4I5,E24.16)
             ELSE
              CALL WADD(I,J,K,L,HHH)
             END IF
            END IF
           END IF
          ELSE
           IF (LAB.EQ.0) THEN
            IF (LABCD.EQ.0) STOP '  why 1 '
C this is for integrals (ss|cd)
            CALL SSCD(HHHBUF)
C distribute the integrals
            INDX=0
            DO IC=1,NINDC
             DO ID=1,NINDD
              INDX=INDX+1
              I=INDXA+1
              J=INDXB+1
              K=INDXC+IC
              L=INDXD+ID
CD        WRITE(6,*) ' I J K L: ',I,J,K,L
              IF (LOCANC(I,J,K,L)) THEN
               HHH=HHHBUF(INDX)
               XMULT=FACTOR(I,J,K,L)
               HHH=XMULT*HHHBUF(INDX)
CD                WRITE(6,9901) I,J,K,L,HHH
               XSSCD=XSSCD+1.D0
               IF (ABS(HHH).GE.THRESH) THEN
                XDUMP=XDUMP+1.D0
                IF (LFORM) THEN
                 WRITE(14,'(4I5,E24.16)') I,J,K,L,HHH
                ELSE
                 CALL WADD(I,J,K,L,HHH)
                END IF
               END IF
              END IF
             END DO
            END DO
           ELSE
            IF (LCD.EQ.0) THEN
             IF (LABCD.EQ.0) STOP '  why 2 '
C this is for integrals (ab|ss)
             CALL ABSS(HHHBUF)
C distribute the integrals
             INDX=0
             DO IA=1,NINDA
              DO IB=1,NINDB
               INDX=INDX+1
               I=INDXA+IA
               J=INDXB+IB
               IF (I.LE.J) THEN
C as I and K come from different shells, the outer loop ensures that
C   I < K and K <= L
                K=INDXC+1
                L=INDXD+1
CD        WRITE(6,*) ' I J K L: ',I,J,K,L
                IF (LOCANC(I,J,K,L)) THEN
                 HHH=HHHBUF(INDX)
                 XMULT=FACTOR(I,J,K,L)
                 HHH=XMULT*HHHBUF(INDX)
CD                WRITE(6,9901) I,J,K,L,HHH
                 XABSS=XABSS+1.D0
                 IF (ABS(HHH).GE.THRESH) THEN
                  XDUMP=XDUMP+1.D0
                  IF (LFORM) THEN
                   WRITE(14,'(4I5,E24.16)') I,J,K,L,HHH
                  ELSE
                   CALL WADD(I,J,K,L,HHH)
                  END IF
                 END IF
                END IF
               END IF
              END DO
             END DO
            ELSE
             IF (LABCD.EQ.0) STOP '  why 3 '
C the general case with the double recursion
 
C calculate the integrals over this quadruplet of shells
             CALL CALBI(HHHBUF)
 
CD        WRITE(6,*) ' INDXA INDXB INDXC INDXD: ',INDXA,INDXB,INDXC
CD   $         ,INDXD
 
             INDX=0
             DO IA=1,NINDA
              DO IB=1,NINDB
               DO IC=1,NINDC
                DO ID=1,NINDD
                 INDX=INDX+1
                 I=INDXA+IA
                 J=INDXB+IB
                 K=INDXC+IC
                 L=INDXD+ID
CD        WRITE(6,*) ' I J K L: ',I,J,K,L
                 IF (LOCANC(I,J,K,L)) THEN
                  HHH=HHHBUF(INDX)
                  XMULT=FACTOR(I,J,K,L)
                  HHH=XMULT*HHHBUF(INDX)
CD                WRITE(6,9901) I,J,K,L,HHH
                  XABCD=XABCD+1.D0
                  IF (ABS(HHH).GE.THRESH) THEN
                   XDUMP=XDUMP+1.D0
                   IF (LFORM) THEN
                    WRITE(14,'(4I5,E24.16)') I,J,K,L,HHH
                   ELSE
                    CALL WADD(I,J,K,L,HHH)
                   END IF
                  END IF
                 END IF
                END DO
               END DO
              END DO
             END DO
C close IF (abss) ...
            END IF
C close IF (sscd) ...
           END IF
C close IF (ssss) ...
          END IF
 
          INDXD=INDXD+NINDD
         END DO
         INDXC=INDXC+NINDC
        END DO
        INDXB=INDXB+NINDB
       END DO
       INDXA=INDXA+NINDA
      END DO
 
      WRITE(6,*) ' Bielectronic integrals : '
      XCALC=XSSSS+XSSCD+XABSS+XABCD
      WRITE(6,9941) XSSSS,XSSCD,XABSS,XABCD,XCALC,XDUMP
 
 9941 FORMAT('  we calculated  ',F23.0,' (ss|ss) integrals ',/
     $     ,'  we calculated  ',F23.0,' (ss|cd) integrals ',/,
     $     '  we calculated  ',F23.0,' (ab|ss) integrals ',/
     $     ,'  we calculated  ',F23.0,' (ab|cd) integrals ',/
     $     ,'  we calculated in total ',F15.0,' integrals ',/
     $     ,'  of these we stored ',F19.0,' integrals ')
 
      IF (LFORM) THEN
       CLOSE(UNIT=14)
      ELSE
       CALL WCLOS(1)
      END IF
 
      RETURN
      END
 
C
      FUNCTION FACTOR(I,J,K,L)
      IMPLICIT NONE
      DOUBLE PRECISION FACTOR
      INTEGER I,J,K,L,IDUM
      LOGICAL LIJ,LKL,LIK,LJL
C
C canonical ordering:
C
C we do not need to test or exchange, as we have the function LOCANC
c$$$C i<j,k<l, i<k; if i=k then j<l
c$$$      IF (I.GT.J) THEN
c$$$       IDUM=J
c$$$       J=I
c$$$       I=IDUM
c$$$      END IF
c$$$      IF (K.GT.L) THEN
c$$$       IDUM=L
c$$$       L=K
c$$$       K=IDUM
c$$$      END IF
c$$$      IF (I.GT.K) THEN
c$$$       IDUM=I
c$$$       I=K
c$$$       K=IDUM
c$$$       IDUM=J
c$$$       J=L
c$$$       L=IDUM
c$$$      ELSE IF (I.EQ.K) THEN
c$$$       IF (J.GT.L) THEN
c$$$        IDUM=J
c$$$        J=L
c$$$        L=IDUM
c$$$       END IF
c$$$      END IF
C
      LIJ=.FALSE.
      LIK=.FALSE.
      LJL=.FALSE.
      LKL=.FALSE.
      IF (I.EQ.J) LIJ=.TRUE.
      IF (K.EQ.L) LKL=.TRUE.
      IF (I.EQ.K) LIK=.TRUE.
      IF (J.EQ.L) LJL=.TRUE.
C
C (aa|aa) 1/8
C (aa|ab) 1/2
C (aa|bb) 1/4
C (ab|ab) 1/2
C (aa|bc) 1/2
C (ab|bc) 1
C (ab|cd) 1
C
 
      FACTOR=1.00D0
C (aa|aa), (aa|bb), (aa|ab), (aa|bc)
      IF (LIJ) FACTOR=FACTOR*0.5D0
C (aa|aa), (aa|bb), (ab|bb), (ab|cc)
      IF (LKL) FACTOR=FACTOR*0.5D0
C (aa|aa)
      IF (LIK.AND.LIJ.AND.LKL) FACTOR=FACTOR*0.5D0
C (ab|ab)
      IF (LIK.AND.LJL.AND..NOT.LIJ) FACTOR=0.50D0
C
      RETURN
      END
C
      SUBROUTINE SSSS(HVAL)
      INCLUDE 'param.h'                                                 
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
 
      DIMENSION XPQ(3)
 
      IONE=1
 
      HVAL=0.D0
 
C this is the direct formula for the (ss|ss) integrals
C
C common quantities
 
      RRAB=0.D0
      RRCD=0.D0
      DO I=1,3
       XAB(I)=XA(I)-XB(I)
       XCD(I)=XC(I)-XD(I)
       RRAB=RRAB+XAB(I)*XAB(I)
       RRCD=RRCD+XCD(I)*XCD(I)
      END DO
 
C loop over contractions
      INDXA=NPSTA
      DO IPRIMA=1,NPRIMA
       EXA=EXX(INDXA)
       COEFFA=COEFF(INDXA)
       INDXA=INDXA+1
 
       INDXB=NPSTB
       DO IPRIMB=1,NPRIMB
        EXB=EXX(INDXB)
        COEFFB=COEFF(INDXB)
        INDXB=INDXB+1
 
        P=EXA+EXB
        DO I=1,3
         XP(I)=(EXA*XA(I)+EXB*XB(I))/P
        END DO
 
        XMU=EXA*EXB/P
 
        CAB=COEFFA*COEFFB
        EMU=EXP(-XMU*RRAB)
        CAB=CAB*EMU
 
        INDXC=NPSTC
        DO IPRIMC=1,NPRIMC
         EXC=EXX(INDXC)
         COEFFC=COEFF(INDXC)
         INDXC=INDXC+1
 
         CABC=CAB*COEFFC
 
         INDXD=NPSTD
         DO IPRIMD=1,NPRIMD
          EXD=EXX(INDXD)
          COEFFD=COEFF(INDXD)
          INDXD=INDXD+1
 
          CABCD=CABC*COEFFD
 
          Q=EXC+EXD
          XNU=EXC*EXD/Q
          ENU=EXP(-XNU*RRCD)
          CABCD=CABCD*ENU
 
          ALPHA=P*Q/(P+Q)
 
          RRPQ=0.D0
          DO I=1,3
           XQ(I)=(EXC*XC(I)+EXD*XD(I))/Q
           XPQ(I)=XP(I)-XQ(I)
           RRPQ=RRPQ+XPQ(I)*XPQ(I)
          END DO
C the arguments for the gamma function
          PRRPQ=ALPHA*RRPQ
          CALL GAMF(F,PRRPQ,IONE)
 
CD        IONE=1
CD        WRITE(6,*) ' 1, F(1) = ',IONE,F(1)
C
          BCON=TP25/(P*Q*SQRT(P+Q))
          PIK=CABCD*BCON
          HVAL=HVAL+PIK*F(1)
C close the four loops over the four contractions (or shells)
         END DO
        END DO
       END DO
      END DO
C
      RETURN
      END
C
      SUBROUTINE PHIREC(CARTS,R000,RQ,EXA,EXB,PIK,LA,LB)
      INCLUDE 'param.h'                                                 
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
 
      DIMENSION RQ(3)
 
C Hermite integrals, input
      DIMENSION R000(*)
C Cartesian integrals, output
      DIMENSION CARTS(NCMAX,NCMAX)
C
C a bunch of Hermite integrals R(tuv;000) is transformed to
C a set of cartesian integrals
C
      DIMENSION PT(2,-1:9,-1:9,-1:9),PTT(2,-1:9,-1:9,-1:9)
      DIMENSION PU(2,-1:9,-1:9),PUU(2,-1:9,-1:9)
      DIMENSION PV(2,-1:9),PVV(2,-1:9)
C
      LAB=LA+LB
      P=EXA+EXB
 
C on input we have EXA, EXB, LA and LB, and the R000
C and the F array, and RQ
 
      Q=(EXA*EXB)/P
c$$$      PIK=COEFFP*2.D0*PI/P
 
CD    WRITE(6,*) ' PIK = ',PIK
C
C for <s|H|s> we have the exact result without recursion
C
      IF (LAB.EQ.0) THEN
       CARTS(1,1)=CARTS(1,1)-PIK*R000(1)
CD     WRITE (6,9901)
CD     WRITE(6,'(2I6,2I3,F20.10)') IONE,IZERO,IZERO,IZERO,-PIK*R000(1)
CD     IZERO=0
CD     IONE=1
CD     WRITE (6,9902)
CD     WRITE(6,'(4I6,5I3,3F15.10)')
CD   -             IONE,IONE,IONE,IZERO,IZERO,IZERO,IZERO,IZERO,
CD   -             IZERO,CARTS(1,1)
       RETURN
      END IF
 
      ITT=1
      IUU=1
      IVV=1
      IITT=2
      IIUU=2
      IIVV=2
      IT=1
      IU=1
      IV=1
      IIT=2
      IIU=2
      IIV=2
      TP=0.5D0/P
      QB=-Q/EXB
      QA=-Q/EXA
C
C well, the minus sign is not really elegant but seemed necessary here
C
      QX=-RQ(1)
      QY=-RQ(2)
      QZ=-RQ(3)
C
      PIS=SQRT(PI/P)
      IND=0
CD    WRITE(6,*)
CD    WRITE(6,*) 'TP,PIS',TP,PIS
CD    WRITE(6,*)
CD    WRITE (6,9901)
 9901 FORMAT(//' R(0;TUV) -- Hermite Gaussians',/,80('='),/,
     -  '      I     T  U  V      R(0;TUV)',/)
      DO I=0,LAB
       DO J=0,LAB-I
        DO K=0,LAB-I-J
         IND=IND+1
         PT(IT,I,J,K)=R000(IND)
CD       WRITE(6,'(2I6,2I3,F20.10)') IND,I,J,K,R000(IND)
        END DO
       END DO
      END DO
CD    WRITE (6,9902)
 9902 FORMAT(//' P(000;IJKLMN) -- cartesian Gaussians',/,
     -     80('='),/,' P = potential ',
     $     /,80('='),/,'     I    LEFT  RIGHT  I  J  K  L  M  N',4X
     $     ,'P(000;IJKLMN)',/)
C
C Rekursion in T
C
      INDI=1
      DO I=0,LA
       IIMAX=LAB-I
       DO III=0,IIMAX
        DO JJJ=0,IIMAX-III
         DO KKK=0,IIMAX-III-JJJ
          PTT(ITT,III,JJJ,KKK)=PT(IT,III,JJJ,KKK)
         END DO
        END DO
       END DO
 
       INDJ=1
       DO II=0,LB
        INDII=INDI-1
C
C Rekursion in U
C
        JMAX=LAB-I-II
        DO JJJ=0,JMAX
         DO KKK=0,JMAX-JJJ
          PU(IU,JJJ,KKK)=PTT(ITT,0,JJJ,KKK)
         END DO
        END DO
        DO J=0,LA-I
         INDII=INDII+1
         INDJJ=INDJ-1
         JJMAX=LAB-I-II-J
         DO JJJ=0,JJMAX
          DO KKK=0,JJMAX-JJJ
           PUU(IUU,JJJ,KKK)=PU(IU,JJJ,KKK)
          END DO
         END DO
         DO JJ=0,LB-II
          INDJJ=INDJJ+1
          KMAX=LAB-I-II-J-JJ
          DO KKK=0,KMAX
           PV(IV,KKK)=PUU(IUU,0,KKK)
          END DO
C
C Rekursion in V
C
          DO K=0,LA-I-J
           KKMAX=LAB-I-J-K-II-JJ
           DO KKK=0,KKMAX
            PVV(IVV,KKK)=PV(IV,KKK)
           END DO
           DO KK=0,KKMAX
            IF (I+J+K.EQ.LA) THEN
             IF (II+JJ+KK.EQ.LB) THEN
              CARTS(INDII,INDJJ)=CARTS(INDII,INDJJ)-PVV(IVV,0)*PIK
CD            WRITE(6,'(3I6,5I3,3F15.10)')
CD   -             INDII,INDJJ,I,J,K,II,JJ,KK,CARTS(INDII,INDJJ)
             END IF
            END IF
C
C (00v;ijkiijjkk) -> (00v-1;ijkiijjkk)
C
            DO KKK=0,LAB-I-J-K-II-JJ-KK-1
             PVV(IIVV,KKK)=TP*PVV(IVV,KKK+1)+QZ*QB*PVV(IVV,KKK)+
     -            DBLE(KKK)*PVV(IVV,KKK-1)
            END DO
            IDUM=IIVV
            IIVV=IVV
            IVV=IDUM
           END DO
C
C (00v;ijkiijj0) -> (00v-1;ijkiijj0)
C
           DO KKK=0,LAB-I-J-K-II-JJ-1
            PV(IIV,KKK)=TP*PV(IV,KKK+1)-
     -           QZ*QA*PV(IV,KKK)+DBLE(KKK)*PV(IV,KKK-1)
           END DO
           IDUM=IIV
           IIV=IV
           IV=IDUM
          END DO
C
C (0uv;ij0iijj0) -> (0u-1v;ij0iijj0)
C
          DO JJJ=0,LAB-I-II-J-JJ-1
           DO KKK=0,LAB-JJJ-1
            PUU(IIUU,JJJ,KKK)=TP*PUU(IUU,JJJ+1,KKK)+
     -           QY*QB*PUU(IUU,JJJ,KKK)+DBLE(JJJ)*PUU(IUU,JJJ-1,KKK)
           END DO
          END DO
          IDUM=IIUU
          IIUU=IUU
          IUU=IDUM
         END DO
C
C (0uv;ij0ii00) -> (0u-1v;ij0ii00)
C
         DO JJJ=0,LAB-I-II-J-1
          DO KKK=0,LAB-JJJ-1
           PU(IIU,JJJ,KKK)=TP*PU(IU,JJJ+1,KKK)-
     -          QY*QA*PU(IU,JJJ,KKK)+DBLE(JJJ)*PU(IU,JJJ-1,KKK)
          END DO
         END DO
         IDUM=IIU
         IIU=IU
         IU=IDUM
        END DO
C
C (tuv;i00ii00) -> (t-1uv;i00ii00)
C
        DO III=0,LAB-I-II-1
         DO JJJ=0,LAB-III-1
          DO KKK=0,LAB-III-JJJ-1
           PTT(IITT,III,JJJ,KKK)=TP*PTT(ITT,III+1,JJJ,KKK)+
     -          QX*QB*PTT(ITT,III,JJJ,KKK)+DBLE(III)*PTT(ITT,III-1,JJJ
     $          ,KKK)
          END DO
         END DO
        END DO
        INDJ=INDJ+LB-II+1
        IDUM=IITT
        IITT=ITT
        ITT=IDUM
       END DO
C
C (tuv;i00000) -> (t-1uv;i00000)
C
       DO III=0,LAB-I-1
        DO JJJ=0,LAB-III-1
         DO KKK=0,LAB-III-JJJ-1
          PT(IIT,III,JJJ,KKK)=TP*PT(IT,III+1,JJJ,KKK)-
     -         QX*QA*PT(IT,III,JJJ,KKK)+
     -         DBLE(III)*PT(IT,III-1,JJJ,KKK)
         END DO
        END DO
       END DO
       INDI=INDI+LA-I+1
       IDUM=IIT
       IIT=IT
       IT=IDUM
      END DO
C
C that was all
C
      RETURN
      END
 
C
      SUBROUTINE SSCD(HHHBUF)
      INCLUDE 'param.h'                                                 
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
C integrals in hermitian Gaussians                                      
      COMMON /INTEGH/ R000(LHMAX)                                       
C..   INCLUDE 'common_integ.h'                                          
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
 
      DIMENSION HHHBUF(*)
      DIMENSION XPQ(3)
      DIMENSION RCART(NCMAX,NCMAX)
      DIMENSION RSPHER(NSPMAX,NSPMAX)
C
C the case (ss|cd)
C
C for the first couple we have no recursion, only for the second
C one we do the McMurchie-Davidson and the Hermite->cartesian->spherical
C harmonics
C
      LCD=LC+LD
 
      LCFIN=NCART(LC)
      LDFIN=NCART(LD)
 
      IF (LCART) THEN
       NUMINT=LCFIN*LDFIN
      ELSE
       NUMINT=(2*LC+1)*(2*LD+1)
      END IF
 
      DO I=1,NUMINT
       HHHBUF(I)=0.D0
      END DO
C
      DO I=1,LCFIN
       DO J=1,LDFIN
        RCART(I,J)=0.D0
       END DO
      END DO
 
C common quantities
 
      RRAB=0.D0
      RRCD=0.D0
      DO I=1,3
       XAB(I)=XA(I)-XB(I)
       XCD(I)=XC(I)-XD(I)
       RRAB=RRAB+XAB(I)*XAB(I)
       RRCD=RRCD+XCD(I)*XCD(I)
      END DO
 
C loop over contractions
      INDXA=NPSTA
      DO IPRIMA=1,NPRIMA
       EXA=EXX(INDXA)
       COEFFA=COEFF(INDXA)
       INDXA=INDXA+1
 
       INDXB=NPSTB
       DO IPRIMB=1,NPRIMB
        EXB=EXX(INDXB)
        COEFFB=COEFF(INDXB)
        INDXB=INDXB+1
 
        P=EXA+EXB
        DO I=1,3
         XP(I)=(EXA*XA(I)+EXB*XB(I))/P
        END DO
 
        XMU=EXA*EXB/P
 
        CAB=COEFFA*COEFFB
        EMU=EXP(-XMU*RRAB)
        CAB=CAB*EMU
 
        INDXC=NPSTC
        DO IPRIMC=1,NPRIMC
         EXC=EXX(INDXC)
         COEFFC=COEFF(INDXC)
         INDXC=INDXC+1
 
         CABC=CAB*COEFFC
 
         INDXD=NPSTD
         DO IPRIMD=1,NPRIMD
          EXD=EXX(INDXD)
          COEFFD=COEFF(INDXD)
          INDXD=INDXD+1
 
          CABCD=CABC*COEFFD
 
          Q=EXC+EXD
          XNU=EXC*EXD/Q
          ENU=EXP(-XNU*RRCD)
          CABCD=CABCD*ENU
 
          ALPHA=P*Q/(P+Q)
 
          RRPQ=0.D0
          DO I=1,3
           XQ(I)=(EXC*XC(I)+EXD*XD(I))/Q
           XPQ(I)=XP(I)-XQ(I)
           RRPQ=RRPQ+XPQ(I)*XPQ(I)
          END DO
C the arguments for the gamma function
          PRRPQ=ALPHA*RRPQ
 
          NMAX=LCD+1
          CALL GAMF(F,PRRPQ,NMAX)
 
CD        IONE=1
CD        WRITE(6,*) ' 1, F(1) = ',IONE,F(1)
C
          BCON=TP25/(P*Q*SQRT(P+Q))
          PIK=-CABCD*BCON
CD         I=1
CD    WRITE(6,*) ' I, F(I) = ',I,F(I)
          AMUL=1.D0
          ACON=ALPHA+ALPHA
          DO I=2,LCD+1
           AMUL=-ACON*AMUL
           F(I)=F(I)*AMUL
CD     WRITE(6,*) ' I, F(I) = ',I,F(I)
CD     WRITE(6,*)
C
          END DO
C
          CALL RMCMUR(R000,F,XPQ,LCD)
C
C add the correct sign to the R000s
C
          INDX=0
          ASIGN=1.D0
          DO I=0,LCD
           BSIGN=ASIGN
           DO J=0,LCD-I
            CSIGN=BSIGN
            DO K=0,LCD-I-J
             INDX=INDX+1
             R000(INDX)=CSIGN*R000(INDX)
             CSIGN=-CSIGN
            END DO
            BSIGN=-BSIGN
           END DO
           ASIGN=-ASIGN
          END DO
C
C the common factor the cartesian recursion is PIK
C
          CALL PHIREC(RCART,R000,XCD,EXC,EXD,PIK,LC,LD)
C
C close the loop over the contractions
         END DO
        END DO
       END DO
      END DO
C
C
C for cartesians we just have to copy the integrals to HHHBUF
C
      IF (LCART) THEN
       INDX=0
       DO I=1,LCFIN
        II=ICARCAR(LC,I)
        DO J=1,LDFIN
         JJ=ICARCAR(LD,J)
         xmult=carmul(lc,ii)*carmul(ld,jj)
         INDX=INDX+1
         HHHBUF(INDX)=RCART(II,JJ)*xmult
        END DO
       END DO
      ELSE
C otherwise we transform the block on spherical harmonics
C
C the routine transforms via the angular momenta LA and LB
C
       LA=LC
       LB=LD
       CALL SPHHAM(RCART,RSPHER)
C restore the momenta
       LA=0
       LB=0
       INDX=0
       DO II=1,2*LC+1
        DO JJ=1,2*LD+1
         INDX=INDX+1
         HHHBUF(INDX)=RSPHER(II,JJ)
        END DO
       END DO
      END IF
C
      RETURN
      END
C
      SUBROUTINE ABSS(HHHBUF)
      INCLUDE 'param.h'                                                 
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
C integrals in hermitian Gaussians                                      
      COMMON /INTEGH/ R000(LHMAX)                                       
C..   INCLUDE 'common_integ.h'                                          
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
 
      DIMENSION HHHBUF(*)
      DIMENSION XPQ(3)
      DIMENSION RCART(NCMAX,NCMAX)
      DIMENSION RSPHER(NSPMAX,NSPMAX)
C
C the case (ab|ss)
C
      LAB=LA+LB
 
      LAFIN=NCART(LA)
      LBFIN=NCART(LB)
 
      IF (LCART) THEN
       NUMINT=LAFIN*LBFIN
      ELSE
       NUMINT=(2*LA+1)*(2*LB+1)
      END IF
 
      DO I=1,NUMINT
       HHHBUF(I)=0.D0
      END DO
 
C
C we should set RCART2 as well to zero as we accumulate for LAB=0
c or LCD=0 through PHIREC
C
      DO I=1,LAFIN
       DO J=1,LBFIN
        RCART(I,J)=0.D0
       END DO
      END DO
 
C common quantities
 
      RRAB=0.D0
      RRCD=0.D0
      DO I=1,3
       XAB(I)=XA(I)-XB(I)
       XCD(I)=XC(I)-XD(I)
       RRAB=RRAB+XAB(I)*XAB(I)
       RRCD=RRCD+XCD(I)*XCD(I)
      END DO
 
C loop over contractions
      INDXA=NPSTA
      DO IPRIMA=1,NPRIMA
       EXA=EXX(INDXA)
       COEFFA=COEFF(INDXA)
       INDXA=INDXA+1
 
       INDXB=NPSTB
       DO IPRIMB=1,NPRIMB
        EXB=EXX(INDXB)
        COEFFB=COEFF(INDXB)
        INDXB=INDXB+1
 
        P=EXA+EXB
        DO I=1,3
         XP(I)=(EXA*XA(I)+EXB*XB(I))/P
        END DO
 
        XMU=EXA*EXB/P
 
        CAB=COEFFA*COEFFB
        EMU=EXP(-XMU*RRAB)
        CAB=CAB*EMU
 
        INDXC=NPSTC
        DO IPRIMC=1,NPRIMC
         EXC=EXX(INDXC)
         COEFFC=COEFF(INDXC)
         INDXC=INDXC+1
 
         CABC=CAB*COEFFC
 
         INDXD=NPSTD
         DO IPRIMD=1,NPRIMD
          EXD=EXX(INDXD)
          COEFFD=COEFF(INDXD)
          INDXD=INDXD+1
 
          CABCD=CABC*COEFFD
 
          Q=EXC+EXD
          XNU=EXC*EXD/Q
          ENU=EXP(-XNU*RRCD)
          CABCD=CABCD*ENU
 
          ALPHA=P*Q/(P+Q)
 
          RRPQ=0.D0
          DO I=1,3
           XQ(I)=(EXC*XC(I)+EXD*XD(I))/Q
           XPQ(I)=XP(I)-XQ(I)
           RRPQ=RRPQ+XPQ(I)*XPQ(I)
          END DO
C the arguments for the gamma function
          PRRPQ=ALPHA*RRPQ
 
          NMAX=LAB+1
          CALL GAMF(F,PRRPQ,NMAX)
 
CD        IONE=1
CD        WRITE(6,*) ' 1, F(1) = ',IONE,F(1)
C
          BCON=TP25/(P*Q*SQRT(P+Q))
          PIK=-CABCD*BCON
C
CD         I=1
CD    WRITE(6,*) ' I, F(I) = ',I,F(I)
          AMUL=1.D0
          ACON=ALPHA+ALPHA
          DO I=2,LAB+1
           AMUL=-ACON*AMUL
           F(I)=F(I)*AMUL
CD     WRITE(6,*) ' I, F(I) = ',I,F(I)
CD     WRITE(6,*)
C
          END DO
C
          CALL RMCMUR(R000,F,XPQ,LAB)
C
C in this case we do not need to take care of the sign t'+u'+v'
C
          CALL PHIREC(RCART,R000,XAB,EXA,EXB,PIK,LA,LB)
C close the four loops over the four contractions (or shells)
         END DO
        END DO
       END DO
      END DO
C
C For cartesians we just have to copy the integrals to HHHBUF
      IF (LCART) THEN
       INDX=0
       DO I=1,LAFIN
        II=ICARCAR(LA,I)
        DO J=1,LBFIN
         JJ=ICARCAR(LB,J)
         xmult=carmul(la,ii)*carmul(lb,jj)
         INDX=INDX+1
         HHHBUF(INDX)=RCART(II,JJ)*xmult
        END DO
       END DO
      ELSE
C
C we have to transform the accumulated Cartesian integrals
C on spherical harmonics if it were not integrals (ss|ss)
C
       CALL SPHHAM(RCART,RSPHER)
       INDX=0
       DO II=1,2*LA+1
        DO JJ=1,2*LB+1
         INDX=INDX+1
         HHHBUF(INDX)=RSPHER(II,JJ)
        END DO
       END DO
      END IF
C
      RETURN
      END
C
      SUBROUTINE CALBI(HHHBUF)
      INCLUDE 'param.h'                                                 
      COMMON /LOCAT/ EXA,EXB,EXC,EXD,XA(3),XB(3),XAB(3),XC(3),XD(3)     
     $     ,XCD(3),P,Q,ALPHA,XP(3),XQ(3),F(4*NLMAX+1),LA,LB,LC,LD       
     $     ,NPRIMA,NPRIMB,NPRIMC,NPRIMD,NPSTA,NPSTB,NPSTC,NPSTD         
C..   INCLUDE 'common_locat.h'                                          
      COMMON /SYST/ NATOM,NBAS,NOCC,NVIRT                               
      COMMON /BAS/ EXX(NPRIMX),COEFF(NPRIMX),POS(3,NATMX),NZ(NATMX)     
     $     ,NSH(NATMX),NPRIM(NSHLMX),IL(NSHLMX),NZENTR(NSHLMX)          
     $     ,NPST(NSHLMX),NSHLT,LMAX                                     
C..   INCLUDE 'common_syst.h'                                           
C integrals in hermitian Gaussians                                      
      COMMON /INTEGH/ R000(LHMAX)                                       
C..   INCLUDE 'common_integ.h'                                          
      COMMON /SPHERI/ SPHCAR(0:NLMAX,NCMAX,NSPMAX),NCART(0:NLMAX)       
     $     ,CARMUL(0:NLMAX,(NLMAX+1)*(NLMAX+2)/2)                       
C..   INCLUDE 'common_spheri.h'                                         
      COMMON /FLOW/ THRESH,MULINT,LFORM,LMONO,LBIEL,LCART,LKINET        
      LOGICAL LFORM,LMONO,LBIEL,LCART,LKINET                            
C..   INCLUDE 'common_flow.h'                                           
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
 
      DIMENSION HHHBUF(*)
      DIMENSION XPQ(3),XTMP(3)
      DIMENSION RCART(NCMAX,NCMAX)
C
      DIMENSION RCARTH(NCMAX,NCMAX,165)
      DIMENSION RCARTF(NCMAX,NCMAX,NCMAX,NCMAX)
      DIMENSION RCACCU(NCMAX,NCMAX,NCMAX,NCMAX)
      DIMENSION RSPHER(NSPMAX,NSPMAX)
      DIMENSION RSAVEH(0:20,0:20,0:20)
 
C we try to program Reine, Tellgren, Helgaker PCCP 9 (2007) 4771
C programming eq 67 and 96, the McMurchie scheme
C
C well, it is more the lecture notes of the
C      Sostrup Slot Summer School
C
C for s functions we have straightforward
C     (ss|ss) = 2/(pq sqrt(p+q)) pi^5/2 * F_0(alpha RRPQ)
C
C for (ss|cd) we do the cartesion and spherical recursion on C and D
C for (ab|ss) we do the cartesion and spherical recursion on A and B
C   from the R000(tuv) in these couples
C
C in general we create all R000(t+t',u+u',v+v')
C then we extract all the R000(t+t',u+u',v+v') for a common t',u',v'
C into an array R000(tuv) and do the cartesian recursion on that
C for all t',u',v' done we get R000(t'u'v') for a common cartesian
C couple (jik,lmn) and do the recursion for (i'j'k',l'm'n')
C
C ----------------------------------------------------------------
C
C how many integrals in this batch ?
C limits in cartesian indices
C save LA and LB
c$$$      WRITE(6,*) ' entering CALBI LA,LB,LC,LD = ',LA,LB,LC,LD
C
C save LA and LB for the transformation cartesian -> spherical
C
      LAS=LA
      LBS=LB
 
      LAB=LA+LB
      LCD=LC+LD
      LABCD=LA+LB+LC+LD
 
      LAFIN=NCART(LA)
      LBFIN=NCART(LB)
      LCFIN=NCART(LC)
      LDFIN=NCART(LD)
 
      IF (LCART) THEN
       NUMINT=LAFIN*LBFIN*LCFIN*LDFIN
      ELSE
       NUMINT=(2*LA+1)*(2*LB+1)*(2*LC+1)*(2*LD+1)
      END IF
 
      DO I=1,NUMINT
       HHHBUF(I)=0.D0
      END DO
C
C we should set RCART2 as well to zero as we accumulate for LAB=0
c or LCD=0 through PHIREC
C
      NCOUPL=NPRIMA*NPRIMB*NPRIMC*NPRIMD
CD    WRITE(6,*) NCOUPL,' integrals (ab|cd) ',NUMINT
 
C set the array for the accumulation to zero
      DO I=1,LAFIN
       DO J=1,LBFIN
        DO K=1,LCFIN
         DO L=1,LDFIN
          RCACCU(I,J,K,L)=0.D0
         END DO
        END DO
       END DO
      END DO
 
C common quantities
 
      RRAB=0.D0
      RRCD=0.D0
      DO I=1,3
       XAB(I)=XA(I)-XB(I)
       XCD(I)=XC(I)-XD(I)
       RRAB=RRAB+XAB(I)*XAB(I)
       RRCD=RRCD+XCD(I)*XCD(I)
      END DO
 
C loop over contractions
      INDXA=NPSTA
      DO IPRIMA=1,NPRIMA
       EXA=EXX(INDXA)
       COEFFA=COEFF(INDXA)
       INDXA=INDXA+1
 
       INDXB=NPSTB
       DO IPRIMB=1,NPRIMB
        EXB=EXX(INDXB)
        COEFFB=COEFF(INDXB)
        INDXB=INDXB+1
 
        P=EXA+EXB
        DO I=1,3
         XP(I)=(EXA*XA(I)+EXB*XB(I))/P
        END DO
 
        XMU=EXA*EXB/P
 
        CAB=COEFFA*COEFFB
        EMU=EXP(-XMU*RRAB)
        CAB=CAB*EMU
 
        INDXC=NPSTC
        DO IPRIMC=1,NPRIMC
         EXC=EXX(INDXC)
         COEFFC=COEFF(INDXC)
         INDXC=INDXC+1
 
         CABC=CAB*COEFFC
 
         INDXD=NPSTD
         DO IPRIMD=1,NPRIMD
          EXD=EXX(INDXD)
          COEFFD=COEFF(INDXD)
          INDXD=INDXD+1
 
          CABCD=CABC*COEFFD
 
          Q=EXC+EXD
          XNU=EXC*EXD/Q
          ENU=EXP(-XNU*RRCD)
          CABCD=CABCD*ENU
 
          ALPHA=P*Q/(P+Q)
 
          RRPQ=0.D0
          DO I=1,3
           XQ(I)=(EXC*XC(I)+EXD*XD(I))/Q
           XPQ(I)=XP(I)-XQ(I)
           RRPQ=RRPQ+XPQ(I)*XPQ(I)
          END DO
C the arguments for the gamma function
          PRRPQ=ALPHA*RRPQ
 
          NMAX=LABCD+1
          CALL GAMF(F,PRRPQ,NMAX)
 
CD        IONE=1
CD        WRITE(6,*) ' 1, F(1) = ',IONE,F(1)
C
          BCON=TP25/(P*Q*SQRT(P+Q))
          PIK=-CABCD*BCON
C
C for integrals (ss|ss) we have already all
C
CD         I=1
CD    WRITE(6,*) ' I, F(I) = ',I,F(I)
          AMUL=1.D0
          ACON=ALPHA+ALPHA
          DO I=2,LABCD+1
           AMUL=-ACON*AMUL
           F(I)=F(I)*AMUL
CD     WRITE(6,*) ' I, F(I) = ',I,F(I)
CD     WRITE(6,*)
C
          END DO
C
C now we have to do the McMurchie-Davidson reccurrence
C we should calculate R_{t+t',u+u',v+v'}(alpha,RPQ) via recursion
C
C t + u + v =l_a+l_b
C t'+ u'+ v'=l_c+l_d -> R_{T,U,V} with T+U+V=l_a+l_b+l_c+l_d
C
C thus we have to do the same recurrence in RMCMUR, but for higher
C indices
C
C the difference to the 1-electronic integrals comes in the recurrence
C of the hermite -> cartesian transformation.
C
          CALL RMCMUR(R000,F,XPQ,LABCD)
C
C we have the array R(t+t',u+u',v+v') which we have to transform to
C cartesian Gaussians Phi(ijk,lmn,IJK,LMN), perhaps in 2 steps:
C R(t+t',u+u',v+v') -> PhiH(t,u,v,IJK,LMN) -> Phi(ijk,lmn,IJK,LMN)
C
C Hermite,Hermite
C these we transform via recursion to cartesian integrals
C
C if l_a+l_b=0 we have to recurse only on C and D
C
C in the same sense, if lc+ld=0, the recursion runs only over A and B
C
C we have to transform A+B and C+D
C
          INDX=0
          DO I=0,LABCD
           DO J=0,LABCD-I
            DO K=0,LABCD-I-J
             INDX=INDX+1
             RSAVEH(I,J,K)=R000(INDX)
            END DO
           END DO
          END DO
C
C we transform C and D first
C
          INDX2=0
          DO I=0,LAB
           DO J=0,LAB-I
            DO K=0,LAB-I-J
             INDX2=INDX2+1
C
C for these indices tuv we have to get all the t' u' v' with t'+u'+v'.le
C .lc+ld; we have to take care of the sign (-1)^(t'+u'+v')
C we fill these in R000, adapt the common locat and transform for C and
C D
C
             INDX=0
             ASIGN=1.D0
             DO II=0,LCD
              BSIGN=ASIGN
              DO JJ=0,LCD-II
               CSIGN=BSIGN
               DO KK=0,LCD-II-JJ
                INDX=INDX+1
                R000(INDX)=RSAVEH(I+II,J+JJ,K+KK)
                CSIGN=-CSIGN
               END DO
               BSIGN=-BSIGN
              END DO
              ASIGN=-ASIGN
             END DO
 
C set the cartesians to zero
             DO II=1,LCFIN
              DO JJ=1,LDFIN
               RCARTH(II,JJ,INDX2)=0.D0
              END DO
             END DO
 
C recursion for transforming C and D on cartesians
 
             DO III=1,3
              XTMP(III)=-XCD(III)
             END DO
 
             CALL PHIREC(RCARTH(1,1,INDX2),R000,XTMP,EXC,EXD,-PIK,LC,LD)
            END DO
           END DO
          END DO
C
C all integrals are half-transformed
C
C the second halftransformation
C
          DO I=1,LCFIN
           DO J=1,LDFIN
C
C the sign is +1
C
            INDX=0
            DO II=0,LAB
             DO JJ=0,LAB-II
              DO KK=0,LAB-II-JJ
               INDX=INDX+1
               R000(INDX)=RCARTH(I,J,INDX)
              END DO
             END DO
            END DO
C
C initialize to zero
            DO II=1,LAFIN
             DO JJ=1,LBFIN
              RCART(II,JJ)=0.D0
             END DO
            END DO
C for the second recursion in PHIREC we use PIK
C
C the sign of the integrals is determined by lc+ld
C
C the sign seems to be (-1)^LCD
            IF (MOD(LCD,2).EQ.0) THEN
             XFACT= 1.D0
            ELSE
             XFACT=-1.D0
            END IF
 
            CALL PHIREC(RCART,R000,XAB,EXA,EXB,XFACT,LA,LB)
 
            DO II=1,LAFIN
             DO JJ=1,LBFIN
              RCACCU(II,JJ,I,J)=RCACCU(II,JJ,I,J)+RCART(II,JJ)
             END DO
            END DO
 
           END DO
          END DO
C
C close the four loops over the four contractions (or shells)
         END DO
        END DO
       END DO
      END DO
C
C For cartesians we just have to copy the integrals to HHHBUF
      IF (LCART) THEN
       INDX=0
       DO I=1,LAFIN
        II=ICARCAR(LA,I)
        DO J=1,LBFIN
         JJ=ICARCAR(LB,J)
         DO K=1,LCFIN
          KK=ICARCAR(LC,K)
          DO L=1,LDFIN
           LL=ICARCAR(LD,L)
C normalize correctly
           xmult=carmul(la,ii)*carmul(lb,jj)*carmul(lc,kk)*carmul(ld,ll)
           INDX=INDX+1
           HHHBUF(INDX)=RCACCU(II,JJ,KK,LL)*xmult
          END DO
         END DO
        END DO
       END DO
      ELSE
C
C we may use the cartesian array RCARTF, as we have less spherical than
C cartesian integrals
C
C save LA, LB
       LAS=LA
       LBS=LB
       DO I=1,LCFIN
        DO J=1,LDFIN
C transform A and B
         CALL SPHHAM(RCACCU(1,1,I,J),RSPHER)
         DO II=1,2*LA+1
          DO JJ=1,2*LB+1
           RCARTF(I,J,II,JJ)=RSPHER(II,JJ)
          END DO
         END DO
        END DO
       END DO
 
C adapt LA and LB in common locat
       LA=LC
       LB=LD
       INDX=0
       DO I=1,2*LAS+1
        DO J=1,2*LBS+1
C transform C and D
         CALL SPHHAM(RCARTF(1,1,I,J),RSPHER)
         DO II=1,2*LC+1
          DO JJ=1,2*LD+1
           INDX=INDX+1
           HHHBUF(INDX)=RSPHER(II,JJ)
          END DO
         END DO
        END DO
       END DO
C restore LA and LB
       LA=LAS
       LB=LBS
 
      END IF
C
      RETURN
      END
C
      FUNCTION LOCANC(I,J,K,L)
      IMPLICIT NONE
      INTEGER I,J,K,L
      LOGICAL LOCANC
C
C test of canonical ordering
C
      LOCANC=.TRUE.
C canonical ordering:
C
C i<j,k<l, i<k; if i=k then j<l
      IF (I.GT.J) THEN
       LOCANC=.FALSE.
      END IF
      IF (K.GT.L) THEN
       LOCANC=.FALSE.
      END IF
      IF (I.GT.K) THEN
       LOCANC=.FALSE.
      ELSE
       IF (I.EQ.K) THEN
        IF (J.GT.L) THEN
         LOCANC=.FALSE.
        END IF
       END IF
      END IF
C
      RETURN
      END
C
      SUBROUTINE GAMTAB(MXABCD)
      INCLUDE 'param.h'                                                 
C static quantities for generating the tables and the interpolation form
C of the incomplete gamma function                                      
      PARAMETER (NTABG=20000)                                           
      COMMON/TABAUX/TABGAM(NTABG),NVALS,NSEGM,NSEGM2,NSEGM3,NSEGM4      
     -     ,NSEGM5,NSEGM6                                               
C..   INCLUDE 'common_tabaux.h'                                         
 
      DIMENSION FF(50)
 
 
      RANGE=12.D0
      GINC=0.1D0
C
C tabulation of the incomplete gamma function F_n(x) for 0 < X < RANGE,
C in intervals GINC
C
C how many entries do we have ?
C per order
c$$$      NSEGM=RANGE/GINC+0.5D0
 
      NSEGM=121
      NSEGM2=2*NSEGM
      NSEGM3=3*NSEGM
      NSEGM4=4*NSEGM
      NSEGM5=5*NSEGM
      NSEGM6=6*NSEGM
C
C total: for the interpolation we need 6 more orders
      NVALS=NSEGM*(MXABCD+7)
      WRITE(6,*) ' GAMTAB : preparing ',NVALS,' entries for MXABCD = '
     $     ,MXABCD
 
      IF (NVALS.GT.NTABG) STOP
     $     ' too high a momentum for the Gamma table '
 
C zero argument: F_n(0) = 1/(2n+1)
C
C we need the table up to J_max+6 for the interpolation
C
      DO I=1,MXABCD+7
       INDX=NSEGM*(I-1)+1
       TABGAM(INDX)=1.D0/DBLE(2*I-1)
c$$$       WRITE(6,*) 'F_n(0) at position ',INDX,' for N= ',I-1
      END DO
 
C fill the highest order for downward recursion
C
C we start at INDX
      INDX=(MXABCD+6)*NSEGM+1
      NTWON=2*(MXABCD+6)
      DO I=2,NSEGM
C
C we had already the first one
C
       INDX=INDX+1
       X=DBLE(I-1)*GINC
C
C formula F_n(x)=sum_{k=0}^{\infty} {{(-x)^k}\over{k! (2n+2k+1)}}
C
       K=0
       NTWOK=0
       TERM=1.D0/DBLE(NTWON+1)
       XK=1.D0
 
       TABGAM(INDX)=TERM
 
 100   CONTINUE
       K=K+1
       NTWOK=NTWOK+2
       XK=-(X/DBLE(K))*XK
       DENOM=DBLE(NTWON+NTWOK+1)
       TERM=XK/DENOM
 
       TABGAM(INDX)=TABGAM(INDX)+TERM
 
       IF (ABS(TERM).GT.1.D-15) THEN
        GO TO 100
       END IF
C
C we continue with the downward recursion
C     F_n(x) = {{2xF_{n+1}(x)+exp(-x)}\over{(2n+1)}}
C
       FN1=TABGAM(INDX)
 
       TWOX=2.D0*X
       EXPFAC=EXP(-X)
       INDX2=INDX
       DO N=MXABCD+5,0,-1
        INDX2=INDX2-NSEGM
        FN=(TWOX*FN1+EXPFAC)/DBLE(2*N+1)
        TABGAM(INDX2)=FN
        FN1=FN
       END DO
 
C close loop over X
      END DO
 
      RETURN
      END
 
      SUBROUTINE GAMF(F,ARG,JMAX)
      INCLUDE 'param.h'                                                 
C static quantities for generating the tables and the interpolation form
C of the incomplete gamma function                                      
      PARAMETER (NTABG=20000)                                           
      COMMON/TABAUX/TABGAM(NTABG),NVALS,NSEGM,NSEGM2,NSEGM3,NSEGM4      
     -     ,NSEGM5,NSEGM6                                               
C..   INCLUDE 'common_tabaux.h'                                         
      COMMON /CONST/ PI,ANTOAU,TP25,ONE                                 
C..   INCLUDE 'constants.h'                                             
 
      DIMENSION F(JMAX)
      DIMENSION FF(50)
C
C for ARG < RANGE we do the interpolation in a polynomial of order 6
C in fact RANGE is 12, in steps of 0.1, as in the paper of McMurchie and
C Davidson, J Comp Phys 26 (1978) 218
C
      IF (ARG.LE.12.D0) THEN
       INDX=10.D0*ARG+1.5D0
       ARGTAB=DBLE(INDX-1)*0.1D0
       DELTAX=ARG-ARGTAB
       INDX=INDX+(JMAX-1)*NSEGM
 
C JMAX is the highest momentum + 1
 
       F(JMAX) = TABGAM(INDX)
     -      - DELTAX*(TABGAM(INDX + NSEGM)
     -      - DELTAX*(TABGAM(INDX + NSEGM2)/2.D0
     -      - DELTAX*(TABGAM(INDX + NSEGM3)/6.D0
     -      - DELTAX*(TABGAM(INDX + NSEGM4)/24.D0
     -      - DELTAX*(TABGAM(INDX + NSEGM5)/120.D0
     -      - DELTAX* TABGAM(INDX + NSEGM6)/720.D0)))))
C downward recursion
       EXPFAC=EXP(-ARG)
       TWOARG=2.D0*ARG
       DO J=JMAX-2,0,-1
        F(J+1)=(TWOARG*F(J+2)+EXPFAC)/DBLE(2*J+1)
       END DO
 
      ELSE
       IF (ARG.LE.30.D0) THEN
        EXPFAC=EXP(-ARG)
        SQPIFAC=0.5D0*SQRT(PI/ARG)
        IF (ARG.LE.15.D0) THEN
C between 12 and 15
         G=((-0.3811559346D0/ARG + 0.321180909D0)/ARG -
     $        0.2473631686D0)/ARG + 0.4999489092D0
        ELSE
         IF (ARG.LE.18.D0) THEN
C between 15 and 18
          G= (0.24642845D0/ARG -0.24249438D0)/ARG + 0.4998436875D0
         ELSE
          IF (ARG.LE.24.D0) THEN
C between 18 and 24
           G=  -0.2152832D0/ARG+0.499093162D0
C between 24 and 30
          ELSE
           G=0.490D0
          END IF
         END IF
        END IF
 
        F(1)=SQPIFAC-EXPFAC/ARG*G
 
C upward recursion
        TWOARG=2.D0*ARG
        DO J=1,JMAX-1
         F(J+1)=(DBLE(2*J-1)*F(J)-EXPFAC)/DBLE(TWOARG)
        END DO
 
       ELSE
C aymptotic formula for ARG > 30
        SQPIFAC=0.5D0*SQRT(PI/ARG)
        F(1)=SQPIFAC
C upward recursion
        EXPFAC=EXP(-ARG)
        TWOARG=2.D0*ARG
        DO J=1,JMAX-1
         F(J+1)=(DBLE(2*J-1)*F(J)-EXPFAC)/DBLE(TWOARG)
        END DO
 
       END IF
      END IF
 
      RETURN
      END
 
      FUNCTION ICARCAR(L,I)
      IMPLICIT NONE
      INTEGER ICARCAR,I,L,LL
C      IF (L.NE.1) THEN
       LL=(L+1)*(L+2)/2
       ICARCAR=LL-I+1
C      ELSE
C       ICARCAR=I
C      END IF
      RETURN
      END