      MODULE SWORD_DATA
C
C Module containing parameters / data required throughout the
C sword program.
C
C Joanne Carr, July 2010
C
      USE precisn, ONLY : wp ! for specifying the kind of reals
      USE params, ONLY : MAXSYM
      IMPLICIT NONE
      SAVE ! for the former common block data

      INTEGER, PARAMETER :: IWH1=(MAXSYM*(MAXSYM+1))/2
      INTEGER, PARAMETER :: IWH2=(IWH1*(IWH1+1))/2 ! 666 when maxsym=8 ; used to dimension the IWHERE array

C---- ILINE is the number of integrals to be printed per line
      INTEGER, PARAMETER :: ILINE=3

      INTEGER, PARAMETER :: IWRITE=6 ! unit number for printing

      INTEGER, PARAMETER :: NOUNIT=444 ! unit for writing total numbers of basis functions for NO runs

      REAL(KIND=wp), PARAMETER :: VSMALL=1.0D-07 ! used as a tolerance on overlap matrix elements

C Data from the ARRAYS and INDS common blocks.  Note that variables set but not used have been deleted.
      INTEGER :: IPFLG, NSYMT
      INTEGER, DIMENSION(MAXSYM,MAXSYM) :: IJPAIR
      INTEGER, DIMENSION(4) :: INTCNT ! dimensioned by the fundamental number of types of 2-electron integrals
      INTEGER, DIMENSION(IWH2,4) :: IWHERE ! second dimension is as described for INTCNT
      INTEGER, DIMENSION(MAXSYM) :: NBFT
C      COMMON /ARRAYS/ NBFT, NZBF, NBFSYM, MT, IWHERE, NSYMT, IPFLG,
C     &                INTCNT
C      COMMON /INDS  / IND, IJPAIR

      END MODULE SWORD_DATA

!*==ord77.spg  processed by SPAG 6.56Rc at 11:57 on 13 May 2010
      PROGRAM ORD77
C***********************************************************************
C
C     Main Program for ORD77 - The Sweden energy integrals ordering
C                              module.
C
C     This section simply allocates dynamic core* and invokes the
C     driver ORDDRV. The main comments on the porting of the code
C     are given in ORDDRV.
C * JMC this section no longer "allocates dynamic core"
C
C     The original code used a 'compiled-in' array for dynamic
C     core so that this is a new section.
C
C     Present version is machine portable
C
C***********************************************************************
      USE sword_data, ONLY : IWRITE
      IMPLICIT NONE
C
C*** Start of declarations rewritten by SPAG
C*** End of declarations rewritten by SPAG
C
      WRITE(IWRITE,10)
 10   FORMAT(//,' PROGRAM SWORD (4 Mar 96) : ',//)
C
C---- Call main routine
C
      CALL ORDDRV
C
      STOP
      END PROGRAM ORD77
!*==orddrv.spg  processed by SPAG 6.56Rc at 11:57 on 13 May 2010
      SUBROUTINE ORDDRV
C***********************************************************************
C
C     Molecule-Sweden Atomic Integrals Ordering Program
C
C     Version: Recoded in Fortran 77 in Rome Italy
C
C     This code takes integrals calculated by MOLECULE and orders
C     them for use by the SCF, MCSCF and TRMO codes. For scattering
C     calculations this module reads the tail integrals computed by
C     GAUSTAIL and subtracts these from the integrals over all space
C     that are computed by MOLECULE.
C
C     Normalization: The Alchemy II version of MOLECULE produces
C                    integrals which are over normalized contracted
C                    Gaussians but un-normalized symmetry adapted
C                    basis functions. In the OT section of Alchemy II
C                    the symmetry adapted overlap matrix is used to
C                    normalize integrals over symmetry adapted basis
C                    basis functions. Sweden did not use this convention
C                    leaving it up to the variational principle in the
C                    SCF program to produce orthonormal MOs; this is
C                    all that is required by the chemistry in fact !
C
C                    In order that we can use Alchemy II produced SCF
C                    (and MCSCF and NO-CI) molecular orbitals, the
C                    integral normalization step has been introduced
C                    into this code.
C
C     Version History:
C
C                   i) Obtained from Jan Almlof at Minnesota Super
C                      Computer Institute, December 1989.
C                  ii) Comment cards added by Charles J Gillan during
C                      visit to Riken in Japan, January 1990.
C                 iii) Implemented on Cray X-MP at Rutherford under
C                      Unicos, April 1990.
C                  iv) Modified to read and subtract tail integrals.
C                      Feb+March 1995.
C                   v) Substantial code re-write with removal of some
C                      redundant options eg. we will never need to
C                      use this code to copy already ordered integrals.
C                      The new code works on IBM 3090 and probably on
C                      other 32 bit machines.
C
C       Logical unit     Contents
C       ------------     --------
C
C             2           Header records defining basis set;
C                         All one electron 'energy' integrals;
C                         Two electron integrals of symmetry blocks IIII
C            81           Two electron integrals blocks of form IJIJ
C            82           Two electron integrals blocks of form IIJJ
C            83           Two electron integrals blocks of form IJKL
C
C     Modus Operandi:
C
C     It is important to realize that the MOLECULE code writes out the
C     indices for ordering the integrals. This is really a glorified
C     GATHER/SCATTER code which writes data structures that can be read
C     by the SWEDEN four index transform code.
C
C***********************************************************************
      USE precisn, ONLY : wp ! for specifying the kind of reals                        
      USE CONSTS, ONLY : IZERO, MINUS7
      USE params, ONLY : LRECL=>LBUF, CISYLIS=>NTYPI, IPD2H, C8STARS, 
     &                   CBLANK, CORD1, CORD2, CTWOTAILS, MAXSYM
      USE global_utils, ONLY : utils_get_bytesizes, COPYD
      USE sword_data, ONLY : NOUNIT, IWRITE, IWHERE, NSYMT, NBFT,
     &                       INTCNT, IJPAIR, IPFLG, IWH2
      IMPLICIT NONE
C
C*** Start of declarations rewritten by SPAG
C
C Local variables
C
      CHARACTER(LEN=192) :: CLABEL
      INTEGER :: I, I1ONLY, ICASE, IDOSPH, IERR, II, IJ, IJI, IJK, 
     &           IJKLTST, IJMAX, IJS, IJTST, IS, IST, ITAPE, ITIME, 
     &           IU1EL, IUNIT, IUORD, J, JEND, K, KK, KL, 
     &           KLS, KLTST, L, LEND, MQF, MQI, 
     &           N1, NBF, NBFA, NBFB, NBTOT, NFTTAIL, NINT, 
     &           NINTS, NREQ, NX, lreal, linteg
      INTEGER, DIMENSION(4) :: IU2EL=(/2, 81, 82, 83/) ! dimensioned for the types of 2-electron integral
      INTEGER, DIMENSION(4) :: IUDELTA=(/180, 181, 182, 183/)
      REAL(KIND=wp) :: POTNUC
      REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: XINT, XNORM
      LOGICAL :: ZTAIL, ZDELTA, LSWORD
C
C*** End of declarations rewritten by SPAG
C
      NAMELIST /INPUTO/ IPFLG, IU2EL, IU1EL, I1ONLY, iuord, NFTTAIL, 
     &   ZTAIL, ZDELTA
C
      CALL WRITE_SVN_INFO(IWRITE)
      WRITE(IWRITE,1000)

      IUNIT=18 ! JMC initialization of a local variable, replacing a data statement
C
C---- Define some constants which are entered into the common blocks (JMC: now variables in the sword_data module)
C     and/or the namelist default definitions
C
      IPFLG=0
      IU1EL=2
      I1ONLY=0
      NFTTAIL=96
      IUORD=0
      ZTAIL=.FALSE.
      ZDELTA = .FALSE.
C
C.... By setting IWHERE so large we will generate errors in the
C     ordering section of the code if there is a mistake.
C
C     If the code is ported to a machine with more than 1 024 Mw
C     then this will have to be changed.
C
      DO I=1, IWH2
C JMC         IWHERE(I,1)=1024000000
C JMC         IWHERE(I,2)=1024000000
C JMC         IWHERE(I,3)=1024000000
C JMC         IWHERE(I,4)=1024000000
C JMC changing this to use the largest integer because if there's an error filling the iwhere array, 
C then when it is used (with things added to it) to provide an index used to access another array then 
C a code compiled with bounds checking should stop with an error.
         IWHERE(I,1)=HUGE(1)
         IWHERE(I,2)=HUGE(1)
         IWHERE(I,3)=HUGE(1)
         IWHERE(I,4)=HUGE(1)
      END DO
C
C---- Create a direct access file on unit ISCR
C
C     OPEN(ISCR,ACCESS='DIRECT',FORM='UNFORMATTED',RECL=4095)
C
C---- Read the namelist &INPUTO which controls the run
C
      READ(5,INPUTO)
      IF(iuord.GT.0)iunit=iuord
C
C---- Notify the user about this data now
C
c      WRITE(IWRITE,1100) IPFLG,I1ONLY,IU2EL,IU1EL,iuord,NFTTAIL,ZTAIL
C
      ITAPE=nfttail
      itime=0
      IF(.NOT.ztail)GO TO 201
 202  OPEN(UNIT=ITAPE,FORM='unformatted')
C
C---- Start Reading from ONE electron integrals tape
C
      READ(ITAPE,ERR=201)CLABEL, NSYMT, (NBFT(I),I=1,NSYMT), POTNUC, 
     &                   IDOSPH
      GO TO 203
 201  itape=iu1el
      itime=itime+1
      IF(itime.GT.1)THEN
         WRITE(iwrite,9900)itape
         STOP
      END IF
      GO TO 202
C
C---- The following information is passed to TRMO77 to tell it
C     that it deals with integrals stored in FULL form that is to
C     not packed. The packing option has been disabled here !
C     It was very Cray specific !
C
 203  NX=-NSYMT
C
C---- Write Header on the output unit and summary info to printer
C
      OPEN(UNIT=iunit,FORM='unformatted')
      WRITE(IUNIT)CLABEL, POTNUC, NX, IDOSPH
      WRITE(IUNIT)(NBFT(I),I=1,NSYMT), izero, izero, izero, izero, 
     &            izero, izero, izero
C
      WRITE(IWRITE,16)CLABEL(1:48)
      WRITE(IWRITE,17)POTNUC
      WRITE(IWRITE,15)(I,NBFT(I),I=1,NSYMT)
 
 
! write numbers of basis functions per symmetry to fort.444 for NO runs
      OPEN(UNIT=nounit,STATUS='unknown')
      WRITE(nounit,*)(NBFT(I),I=1,NSYMT)
      CLOSE(nounit)
 
 
      WRITE(IWRITE,4447)
      DO I=1, NSYMT
         WRITE(IWRITE,14)(IPD2H(I,J),J=1,NSYMT)
      END DO
C
C---- Copy over the rest of the header data - this defines the basis set
C     and so on. The routine COPYD has been generalized and is used in
C     the REFORM code too.
C     Header copy takes place from ITAPE ---> IUNIT
C
      LSWORD=.TRUE.
      CALL COPYD(ITAPE,IWRITE,IUNIT,LSWORD)
C
C=======================================================================
C
C     One electron integral ordering section
C
C=======================================================================
C
C---- Place header record on output dataset of ordered integrals
C
      WRITE(IUNIT)C8STARS, CBLANK, CBLANK, CORD1

CCC RUI
cc      WRITE(6,*) 'index transform code   ', (A(I),I=1,NLOWR)
      
C
C---- Compute : the total number of basis function, NBTOT;
C               the total storage for all lower half triangles, N1;
C               pointers into each symmetries lower triangle, NZBF. (JMC set but not used so removed)
C
      NBTOT=0
      N1=0
C
      DO I=1, NSYMT
         N1=N1+NBFT(I)*(NBFT(I)+1)/2
         NBTOT=NBTOT+NBFT(I)
      END DO
C
C---- Make core allocations for standard case (no extra potentials)
C
C     It is useful to remember that:
C
C     INORM points to storage for normalization coefficients for
C                                 each sym. adpt. basis function.
C     IXINT  " "   to  " " "  for buffer used to read un-ordered
C                                 integrals written by MOLECULE
C      IIND  " "   to  " " "  for buffer to read the indices that
C                             accompany the values in IXINT
C      IORD  " "   to  " " "  for the ordered integrals
C      IMAP  " "   to  " " "  for --- not used -----
C       IA2  " "   to  " " "  for --- not used -----
C
C     There are LRECL/2 integrals and LRECL/2 indices per record ! JMC ??? don't believe that...
C
      ALLOCATE(xnorm(maxsym*nbtot))
C
C---- Invoke the subroutine which computes symmetry adapted basis
C     function normalization factors. These are obtained from the
C     overlap integrals.
C
      CALL NORMAL(NSYMT,NBFT,XNORM,IWRITE,N1,lrecl,iu1el,NFTTAIL,ZTAIL)
C
C---- Invoke routine which reads, orders and writes the 1-el integrals
C
      CALL ORD1EL(iu1el,IUNIT,NFTTAIL,ZTAIL,lrecl,N1,NSYMT,NBFT,IWRITE,
     &            XNORM,IPFLG)
C
C---- Normalize the raw basis function amplitudes now - this is
C     only needed in a scattering run which will be indicated
C     by having the tails flag switched on.
C
      IF(ZTAIL)CALL AMPERE(NFTTAIL,IUNIT,IWRITE,NBFT,NSYMT,XNORM)
C
C---- Terminate if only one electron integral ordering is required
C
      WRITE(IWRITE,4400)
      IF(I1ONLY.NE.0)THEN
         WRITE(IWRITE,
     &             FMT='(//,5X,''**** Only one electron integral '',/)')
         WRITE(IWRITE,
     &             FMT='(   5X,''     ordering has been requested'',/)')
         WRITE(IWRITE,
     &             FMT='(   5X,''     Successful completion      '',/)')
         STOP
      END IF
C
C=======================================================================
C
C     Two electron integral ordering section
C
C=======================================================================
C
      ITAPE=IU2EL(1)
C
C---- Write header on the output dataset of ordered integrals
C
      WRITE(IUNIT)C8STARS, CBLANK, CBLANK, CORD2
C
      DO I=1, 4
         INTCNT(I)=0
      END DO
C
      MQI=1
C
C---- All IIII blocks of integrals
C
      DO I=1, NSYMT
         NBF=NBFT(I)
         MQF=NBF+MQI-1
         MQI=MQF+1
         IJ=NBF*(NBF+1)/2
         NINTS=IJ*IJ
         IS=I*(I+1)/2
         IST=IS*(IS+1)/2
         IWHERE(IST,1)=INTCNT(1)
         INTCNT(1)=INTCNT(1)+NINTS
      END DO
C
C---- If there is only one symmetry the jump go straight to
C     the ordering section
C
      IF(NSYMT.LE.1)GO TO 110
C
C---- All IJIJ blocks of integrals followed by IIJJ blocks
C
      DO I=2, NSYMT
         JEND=I-1
         DO J=1, JEND
            IJ=NBFT(I)*NBFT(J)
            NINTS=IJ*IJ
            IJS=I*(I-1)/2+J
            IST=IJS*(IJS-1)/2+IJS
            IWHERE(IST,2)=INTCNT(2)
            INTCNT(2)=INTCNT(2)+NINTS
         END DO
      END DO
C
C.... The IIJJ blocks of integrals
C
      DO I=2, NSYMT
         JEND=I-1
         NBFA=NBFT(I)
         II=NBFA*(NBFA+1)/2
         DO J=1, JEND
            NBFB=NBFT(J)
            KK=NBFB*(NBFB+1)/2
            NINTS=II*KK
            IJI=I*(I+1)/2
            IJK=J*(J+1)/2
            IST=IJI*(IJI-1)/2+IJK
            IWHERE(IST,3)=INTCNT(3)
            INTCNT(3)=INTCNT(3)+NINTS
         END DO
      END DO
C
C---- If the number of symmetries is less than four there are
C     no IJKL blocks. Go straight to the ordering section
C
      IF(NSYMT.GE.4)THEN
C
         DO I=2, NSYMT
            JEND=I-1
            DO J=1, JEND
               DO K=2, I
                  LEND=K-1
                  IF(I.EQ.K)LEND=J
                  DO L=1, LEND
                     IF(I.EQ.K .AND. J.EQ.L)CYCLE
                     IJTST=IPD2H(I,J)
                     KLTST=IPD2H(K,L)
                     IJKLTST=IPD2H(IJTST,KLTST)
                     IF(IJKLTST.NE.1)CYCLE
                     IJS=((I-1)*I)/2+J
                     KLS=((K-1)*K)/2+L
                     IST=IJS*(IJS-1)/2+KLS
                     IJ=NBFT(I)*NBFT(J)
                     KL=NBFT(K)*NBFT(L)
                     NINTS=IJ*KL
                     IWHERE(IST,4)=INTCNT(4)
                     INTCNT(4)=INTCNT(4)+NINTS
                  END DO
               END DO
            END DO
         END DO
      END IF
C
C---- Now generate IJPAIR array for ordering use.
C
 110  ijmax=0
      DO I=1, NSYMT
         DO J=1, i-1
            IJPAIR(I,J)=NBFT(I)*NBFT(J)
            IF(ijpair(i,j).GT.ijmax)ijmax=ijpair(i,j)
            IJPAIR(J,I)=IJPAIR(I,J)
         END DO
         IJPAIR(I,I)=NBFT(I)*(NBFT(I)+1)/2
         IF(ijpair(i,i).GT.ijmax)ijmax=ijpair(i,i)
      END DO
C
      J=0
      DO I=1, NBTOT
         J=J+I
      END DO
C
C---- Perform the sorting procedure for each type of integral. The
C     intial statements in the loop to 88888 set up the variables
C     and buffers. The actual ordering is done in routine SRT2EL
C
      DO ICASE=1, 4
C
C....    Case numbers here are 1) IIII  2) IJIJ  3) IIJJ  4) IJKL
C
         IF(ICASE.EQ.4 .AND. NSYMT.LT.4)CYCLE
         IF(ICASE.GT.1 .AND. NSYMT.LT.2)CYCLE
C
         NINT=INTCNT(ICASE)
         IF (NINT .LT. 0) THEN
            WRITE(iwrite,9903) NINT
            STOP
         END IF
         WRITE(IWRITE,87889)CISYLIS(ICASE), NINT
C
C..... Invoke the routine SRT2EL which reads and orders the integrals
C
C      Remember that in Sweden the four types of two electron integrals
C      are stored on different units. A code has been added to split
C      the Alchemy II MOLECULE output to conform to Sweden conventions.
C
         ITAPE=IU2EL(ICASE)
         IF (ZDELTA) THEN
           ITAPE = IUDELTA(ICASE)
         ENDIF
         write(*,*) '***** reading file: case, unit' ,itape, icase
C
         CALL utils_get_bytesizes(linteg, lreal)
C JMC 3 real arrays of dimension nint, ijmax and lrecl, and 11 integer arrays of dim. lrecl are declared in srt2el
         nreq=(lreal*(nint+ijmax+lrecl)+(11*linteg*lrecl))/1000000
         WRITE(iwrite,9902)nreq
         ALLOCATE(xint(nint),stat=ierr)
         IF(ierr.NE.0)THEN
            WRITE(iwrite,9901)
            STOP
         END IF
c
         CALL SRT2EL(xint,NINT,ITAPE,LRECL,IUNIT,ijmax,ICASE,XNORM,
     &               ZTAIL,NFTTAIL,CTWOTAILS(icase),nbtot,zdelta)
c
         DEALLOCATE(xint)
C
         WRITE(IWRITE,4500)CISYLIS(ICASE)
C
      END DO
C
C---- Write terminator string on the output file of integrals
C     and terminate the run itself.
C
      WRITE(IUNIT)MINUS7, MINUS7, MINUS7, MINUS7, MINUS7, MINUS7, 
     &            MINUS7, MINUS7
c
      DEALLOCATE(xnorm)
c
      STOP
C
C---- Format Statements
C
 14   FORMAT(8I5)
 15   FORMAT(//' BASIS FUNCTIONS PER SYMMETRY'/(2I5))
 16   FORMAT(/' TITLE HEADER ON INTEGRALS UNIT: ',/,1X,20A)
 17   FORMAT(/' NUCLEAR REPULSION ENERGY =',F15.8)
C
 1000 FORMAT(///,1X,
     &       'A T O M I C  I N T E G R A L  O R D E R I N G  C O D E ',
     &       //,1X,
     &       'Version modified for R-matrix calculations',//)
 1100 FORMAT(/' Namelist &INPUTO Contents : ',//,1X,'I1ONLY = ',I3,
     &       /' IPFLG  = ',I3,/,1X,'IU2EL  = ',4I3,' IU1EL  = ',I3,
     &       ' IUORD ',I3,/1X,'NFTTAIL= ',I3,'   ZTAIL  = ',L7,/)
87889 FORMAT(' In core sort of case ',A,4X,I10)
 4400 FORMAT(/,' One-electron integral ordering completed ')
 4447 FORMAT(//,' GROUP MULTIPLICATION TABLE ',/)
 4500 FORMAT(' Two-electron ordering for ',A,' completed '/)
 9900 FORMAT(/' ***** Error reading unit ',i3,'  *****'/)
 9901 FORMAT(/' Unable to allocate memory for integral storage')
 9902 FORMAT(/' Memory required is ',I10,' Mbytes')
 9903 FORMAT(/' ERROR: integer out of range. NINT= ',I20,/,
     * ' Recompile for longer integers and rerun everything necessary',
     *       //' Stopping')
C
      END SUBROUTINE ORDDRV
!*==ord1el.spg  processed by SPAG 6.56Rc at 11:57 on 13 May 2010
      SUBROUTINE ORD1EL(ITAPE,IOUT,NFTTAIL,ZTAIL,LTRI,NLOWR,NSYMT,NBFT,
     &                  IWRITE,XNORM,IPFLG)
C***********************************************************************
C
C     Orders the one electron integrals removing R-matrix tails if
C     requested by the user.
C
C     Input data:
C             VV Buffer for reading integrals
C         INDEXV Buffer for reading indices for each integral
C              A buffer for holding ordered one electron integrals.
C          ITAPE Logical unit holding the energy integrals written by
C                MOLECULE.
C           IOUT Logical unit for the output of ordered integrals.
C        NFTTAIL Logical unit for the R-matrix tail integrals
C          ZTAIL Flag indicating if there are R-matrix tails to be
C                subtracted.
C           LTRI Number of integrals on each record
C          NLOWR Total number of integrals summed over all lower
C                half triangles of integrals.
C          NSYMT Number of symmetries in the basis function set
C           NBFT Number of basis function per symmetry.
C         IWRITE Logical unit for the printer
C          XNORM Normalization factor for each symmetry adapted basis
C                function
C          IPFLG Print flag
C
C     Output data:
C
C         None - all data is written to disk
C
C     Copyright 1995 (c) Charles J Gillan
C     All rights reserved
C
C***********************************************************************
      USE precisn, ONLY : wp ! for specifying the kind of reals                        
      USE CONSTS, ONLY : XZERO
      USE params, ONLY : CKINET, CONELH, COVERLAP, MAXSYM
      USE global_utils, ONLY : search
      IMPLICIT NONE
C
C*** Start of declarations rewritten by SPAG
C
C Dummy arguments
C
      INTEGER :: IOUT, IPFLG, ITAPE, IWRITE, LTRI, NFTTAIL, NLOWR, NSYMT
      LOGICAL :: ZTAIL
      INTEGER, DIMENSION(*) :: NBFT
      REAL(KIND=wp), DIMENSION(MAXSYM,*) :: XNORM
      INTENT (IN) IOUT, IPFLG, IWRITE, LTRI, NBFT, NLOWR, NSYMT, XNORM, 
     &            ZTAIL
C
C Local variables
C
      REAL(KIND=wp), DIMENSION(nlowr) :: A
C     The ordering here (in CHEAD) actually corresponds to the order in which
C     these integrals are stored on the file of ordered integrals.
      CHARACTER(LEN=8), DIMENSION(3) :: CHEAD=(/COVERLAP,CONELH,CKINET/)
      INTEGER :: I, IEND, IFAIL, IT, J, LSYM, N, KK
      INTEGER, DIMENSION(ltri) :: INDEXV
      REAL(KIND=wp), DIMENSION(ltri) :: VV
C
C*** End of declarations rewritten by SPAG
C
C---- Debug header
C
c        WRITE(IWRITE,1000)
c        WRITE(IWRITE,1010) ZTAIL,NFTTAIL,ITAPE,IOUT,LTRI
c        WRITE(IWRITE,1020) NSYMT,NLOWR,(NBFT(I),I=1,NSYMT)
C
C=======================================================================
C
C     Ordering of energy integrals and, when requested, removal of
C     R-matrix tails.
C
C=======================================================================
C
C---- Loop over the three types of one electron energy integral that
C     we seek to order.
C
      DO IT=1, 3
C
c        WRITE(IWRITE,FMT='(/,5X,''INTS - Looking for '',A,/)') CHEAD(IT)
C
C.... Locate the first record of these one electron integrals written
C     by MOLECULE.
C
         CALL SEARCH(ITAPE,CHEAD(IT),ifail,iwrite)
         IF(ifail.NE.0)STOP 999
C
C..... Initialize the buffer of ordered integrals now
C
         DO I=1, NLOWR
            A(I)=XZERO
         END DO
C
C..... Read all records of integrals for this type and order these.
C
C      The following is a repeat until loop in Fortran where we read
C      records of integrals until the IEND value is negative.
C      JMC changed the GO TO structure to the following do loop...
C
         DO KK = 1, HUGE(1) ! make sure this loop is technically not infinite
            READ(ITAPE)(VV(I),I=1,LTRI), (INDEXV(I),I=1,LTRI), IEND
C
C.... Debug output for the user
C
c        WRITE(IWRITE,FMT='(/,10X,''Record of Integrals: Number = '',I5,/)')
c     *        IEND
c        WRITE(IWRITE,FMT='(10X,I5,1X,F20.10,1X,Z8,1X,I6)')
c     *        (I,VV(I),INDEXV(I),INDEXV(I),I=1,IEND)
C
C.... If we have integrals then order them, otherwise process accordingl
C
            IF(IEND.GT.0)THEN
               DO I=1, IEND
                  A(INDEXV(I))=VV(I)
               END DO
            ELSE IF(IEND.LT.0)THEN
               EXIT
            END IF
         END DO
C
C---- Removal of R-matrix tail integrals
C
         IF(ZTAIL)THEN
C
            CALL SEARCH(NFTTAIL,CHEAD(IT),ifail,iwrite)
            IF(ifail.NE.0)STOP 999
C
C..... Read all records of integrals for this type and order these.
C
C      The following is a repeat until loop in Fortran where we read
C      records of integrals until the IEND value is negative.
C      JMC changed the GO TO structure to the following do loop...
C
            DO KK = 1, HUGE(1) ! make sure this loop is technically not infinite
               READ(NFTTAIL)(VV(I),I=1,LTRI), (INDEXV(I),I=1,LTRI), IEND
               IF(IEND.GE.0)THEN
                  DO I=1, IEND
                     A(INDEXV(I))=A(INDEXV(I))-VV(I)
                  END DO
               ELSE
                  EXIT
               END IF
            END DO
         END IF
C
C---- Normalize the basis functions in each integral now.
C
C     Here we exploit the fact the matrix is ordered and we can build
C     an appropriate loop structure.
C
C     For an R-matrix run the normalization factors will be for the
C     inner sphere - although this only affects continuum functions
C
         N=0
C
         DO LSYM=1, NSYMT
            DO I=1, NBFT(LSYM)
               DO J=1, I
                  N=N+1
                  A(N)=A(N)*XNORM(LSYM,I)*XNORM(LSYM,J)
               END DO
            END DO
         END DO
C
C---- Write the ordered integrals on disk now for use by the four index
C     transform code. Also, print these for the user if required.
C
         WRITE(IOUT)(A(I),I=1,NLOWR)
C
         IF(IPFLG.NE.0)THEN
            WRITE(IWRITE,7000)CHEAD(IT)
            WRITE(IWRITE,7010)(I,A(I),I=1,NLOWR)
         END IF
C
      END DO
C
C---- Subroutine return point for R-matrix case and for bound state.
C
c        WRITE(IWRITE,8000)
C
      RETURN
C
C---- Format Statements
C
 1000 FORMAT(//,5X,'====> ORD1EL - ORDER 1 El. <====',/)
 1010 FORMAT(5X,'Value of tail switch      (ZTAIL) = ',L7,/,5X,
     &       'Logical unit for tails  (NFTTAIL) = ',I5,/,5X,
     &       'Logical unit for raw ints (ITAPE) = ',I5,/,5X,
     &       'Logical unit for ord ints  (IOUT) = ',I5,/,5X,
     &       'Integrals per record       (LTRI) = ',I5,/)
 1020 FORMAT(5X,'Number of symmetries      (NSYMT) = ',I5,/,5X,
     &       'No. of ordered integrals  (NBTOT) = ',I5,/,5X,
     &       'Basis functions per symmetry : ',/,5X,8(I3,1X))
 7000 FORMAT(/,5X,'Ordered (Normalized) integrals : ',A,/)
 7010 FORMAT((5X,4(I4,1X,D13.7,1X)))
 8000 FORMAT(//,5X,'====== END OF ORD1El ===========',/)
C
      END SUBROUTINE ORD1EL
!*==srt2el.spg  processed by SPAG 6.56Rc at 11:57 on 13 May 2010
      SUBROUTINE SRT2EL(a,NINTS,ITAPE2,LRECL,IUNIT,ijmax,ICASE,XNORM,
     &                  ZTAIL,NFTTAIL,ctail,nbtot,zdelta)
C***********************************************************************
C
C     SRT2EL - Sorts the two electron integrals for a given case i.e.
C
C       IIII, IJIJ, IIJJ or IJKL referred to as ICASE = 1, 2, 3 or 4
C
C     Input data:
C           ITAPE Logical unit holding the two electron integrals
C                 for the present case. Note that these begi
C                 with the header card:
C
C                 ********                TWOELSUP
C
C                 and so each case is assumed to be held on a different
C                 logical unit. This routine is called one for
C                 each case.
C              VV Buffer used to read records of un-ordered two
C                 electron integrals from ITAPE.
C          INDEXV Buffer used to read the symmetry and function
C                 labels associated with each integral in VV.
C     IS,JS,KS,LS Buffers which hold the symmetry label for each
C                 atomic basis function in the two electron integrals.
C     IV,JV,KV,LV Buffers which hold the atomic basis function numbers
C                 for each function in the two electron integrals.
C                 Note that (IS,IV) form a pair etc....
C           XNORM Array of basis function normalization factors.
C           ZTAIL Switch which denotes whether, or not, two electron
C                 tail integrals should be handled
C         NFTTAIL Logical unit number for dataset holding the two
C                 electron tail integrals for this value of ICASE.
C
C     Notes:
C
C     The integrals are read from a dataset, there is a different one
C     for each case, ordered in-core and the written to IUNIT. All
C     cases end up on IUNIT and that file is subsequently used by
C     the FJK and TRMO codes.
C
C     This routine has been substantially recoded. In particular
C     the out of core option has been removed. Also basis function
C     normalization has been introduced. In Sweden the symmetry
C     adapted atomic basis functions are not, in general, normalized
C     by MOLECULE. Note that this leads to different MO coefficients
C     at the SCF stage but, of course, the same values for physical
C     observables.
C
C***********************************************************************
      USE precisn, ONLY : wp ! for specifying the kind of reals                        
      USE CONSTS, ONLY : XZERO, MINUS1, MINUS2, MINUS3, MINUS6, IZERO
      USE params, ONLY : CTWOEL, MAXSYM, IPD2H
      USE integer_packing, ONLY : sword_unpack
      USE global_utils, ONLY : INDFUNC, search
      USE sword_data, ONLY : IWRITE, IWHERE, IJPAIR, NBFT, NSYMT, IPFLG
      IMPLICIT NONE
C
C*** Start of declarations rewritten by SPAG
C
C Dummy arguments
C
      CHARACTER(LEN=8) :: CTAIL
      INTEGER :: ICASE, IJMAX, ITAPE2, IUNIT, LRECL, NBTOT, NFTTAIL, 
     &           NINTS
      LOGICAL :: ZTAIL
      LOGICAL :: zdelta
      REAL(KIND=wp), DIMENSION(nints) :: A
      REAL(KIND=wp), DIMENSION(MAXSYM,nbtot) :: XNORM
      INTENT (IN) CTAIL, ICASE, IJMAX, ITAPE2, IUNIT, LRECL, NBTOT, 
     &            NFTTAIL, NINTS, XNORM, ZTAIL
      INTENT (INOUT) A
C
C Local variables
C
      REAL(KIND=wp), DIMENSION(ijmax) :: BUF
      CHARACTER(LEN=8) :: CHEADER
      INTEGER :: I, IEND, IFAIL, II, IJ, IJKLT, IJMX, IJTST, INDX1, 
     &           INDX2, IOFFI, IPASS, ITAPE, ITIME, IWH, IXX, J, K, KL, 
     &           KLMX, KLTST, L, LEND
      INTEGER, DIMENSION(lrecl) :: INDEX1, INDEX2, IOFF, IS, IV, JS, JV, 
     &                             KS, KV, LS, LV
      LOGICAL :: LCONN
      REAL(KIND=wp), DIMENSION(LRECL) :: VV
C
C*** End of declarations rewritten by SPAG
C
C---- The integrals are writen in full storage format.
C
      DO I=1, NINTS
         A(I)=XZERO
      END DO
C
C---- Find the header which marks the start of the two electron
C     integrals. Remember that ITAPE is a different unit for each
C     block of integrals.
C
      ipass=1
      cheader=ctwoel
      itape=itape2
c
 200  INQUIRE(UNIT=itape,OPENED=lconn)
      IF(lconn)THEN
         REWIND itape
      ELSE
         OPEN(UNIT=itape,FORM='unformatted')
      END IF
      CALL SEARCH(ITAPE,CHEADER,ifail,iwrite)
      IF(ifail.NE.0)THEN
         IF(ipass.EQ.1)THEN
            STOP 999
         ELSE
            WRITE(IWRITE,*)' No tail subtraction for this case'
            GO TO 230
         END IF
      END IF
      itime=0
C
C=======================================================================
C
C     Read and order all integrals
C
C=======================================================================
C
C---- Following is a repeat until loop in Fortran. We read records of
C     integrals until IEND=-1. Successive iterations around the loop
C     begin at line 102.
C
 102  READ(ITAPE)(VV(I),I=1,LRECL), (INDEX1(i),index2(i),I=1,LRECL),IEND
      itime=itime+1
C
C.... Debug output for the user
C
c        WRITE(IWRITE,FMT='(/,10X,''Record of Integrals: Number = '',I5,/)')
c     *        IEND
c        WRITE(IWRITE,FMT='(10X,I5,1X,F20.10,1X,Z8,1X,Z8)')
c     *          (I,VV(I),INDEX1(I),INDEX2(I),I=1,IEND)
C
      IF(IEND.LT.0)GO TO 230
C
C....... Unpack indices and form array pointers
C
C        Packed data consists of function numbers, within symmetries,
C        as well as symmetry numbers too.
C
C        This remains machine dependent - JMC hence it has been moved to an interfacing module
C
      DO II=1, IEND
C JMC         IS(II)=IAND(ISHFT(INDEX1(II),-25),15)
C JMC         IV(II)=IAND(ISHFT(INDEX1(II),-15),1023)
C
C JMC         JS(II)=IAND(ISHFT(INDEX1(II),-10),15)
C JMC         JV(II)=IAND(INDEX1(II),1023)
         CALL SWORD_UNPACK(INDEX1(II),IS(II),IV(II),JS(II),JV(II))
C
C JMC         KS(II)=IAND(ISHFT(INDEX2(II),-25),15)
C JMC         KV(II)=IAND(ISHFT(INDEX2(II),-15),1023)
C
C JMC         LS(II)=IAND(ISHFT(INDEX2(II),-10),15)
C JMC         LV(II)=IAND(INDEX2(II),1023)
         CALL SWORD_UNPACK(INDEX2(II),KS(II),KV(II),LS(II),LV(II))
C
C....... Compute the storage pointer now into the triangulation block
C        and from it the offset for storage of the ordered integrals
C
         IWH=INDFUNC(INDFUNC(IS(II),JS(II)),INDFUNC(KS(II),LS(II)))
         IOFF(II)=IWHERE(IWH,ICASE)
c
C..... As for the one electron integrals, we now normalize the
C      symmetry adapted basis functions in each integral. This
C      was not done in the original Sweden work.
C
C      Read the comments in ORD1EL and at the header of the program
C      for a discussion of this process.
C
         VV(II)=VV(II)*XNORM(IS(II),IV(II))*XNORM(JS(II),JV(II))
     &          *XNORM(KS(II),KV(II))*XNORM(LS(II),LV(II))
C
      END DO
c
      IF(ipass.EQ.2)THEN
         DO ii=1, iend
            vv(ii)=-vv(ii)
         END DO
      END IF
C
C.... Now compute the integral index in the ordered list. In D2h there
C     are only four cases to consider. These are defined by the symmetry
C     representations of the integrals.
C
C     The four cases are:
C
C     1) IIII - all the same symmetry
C     2) IJIJ - two symmetries
C     3) IIJJ - two symmetries
C     4) IJKL - all different symmetries
C
      IF(ICASE.EQ.1)THEN
C
         DO II=1, IEND
            INDX1=(INDFUNC(IV(II),JV(II))-1)*IJPAIR(IS(II),IS(II))
     &            +INDFUNC(KV(II),LV(II))+IOFF(II)
            A(INDX1)=A(INDX1)+VV(II)
            INDX2=(INDFUNC(KV(II),LV(II))-1)*IJPAIR(IS(II),IS(II))
     &            +INDFUNC(IV(II),JV(II))+IOFF(II)
            IF(indx1.EQ.indx2)CYCLE
            A(INDX2)=A(INDX2)+VV(II)
         END DO
      ELSE IF(ICASE.EQ.2)THEN
C
         DO II=1, IEND
            INDX1=((IV(II)-1)*NBFT(JS(II))+JV(II)-1)
     &            *IJPAIR(IS(II),JS(II))+(KV(II)-1)*NBFT(JS(II))+LV(II)
     &            +IOFF(II)
            A(INDX1)=A(INDX1)+VV(II)
            INDX2=((KV(II)-1)*NBFT(LS(II))+LV(II)-1)
     &            *IJPAIR(KS(II),LS(II))+(IV(II)-1)*NBFT(LS(II))+JV(II)
     &            +IOFF(II)
            IF(indx1.EQ.indx2)CYCLE
            A(INDX2)=A(INDX2)+VV(II)
         END DO
      ELSE IF(ICASE.EQ.3)THEN
C
         DO II=1, IEND
            INDX1=(INDFUNC(KV(II),LV(II))-1)*IJPAIR(IS(II),JS(II))
     &            +INDFUNC(IV(II),JV(II))+IOFF(II)
            A(INDX1)=A(INDX1)+VV(II)
         END DO
      ELSE IF(ICASE.EQ.4)THEN
C
         DO II=1, IEND
            INDX1=((KV(II)-1)*NBFT(LS(II))+LV(II)-1)
     &            *IJPAIR(IS(II),JS(II))+(IV(II)-1)*NBFT(JS(II))+JV(II)
     &            +IOFF(II)
            A(INDX1)=A(INDX1)+VV(II)
         END DO
      END IF
C
C.... Loop back to read a new record of integrals
C
      GO TO 102
C
C---- The whole of the two el. integrals have now been read.
C     Set up read of tails, if required.
C
 230  IF(ztail .AND. ipass.EQ.1)THEN
         IF (zdelta.eq.0) THEN
           ipass=2
           ITAPE=nfttail
           CHEADER=ctail
           GO TO 200
         END IF
      END IF
C
C=======================================================================
C
C     Take the ordered integrals and place them on the output dataset.
C
C=======================================================================
C
C---- Branch to the section of code that is appropriate to the
C     block of integrals in question.
C
      SELECT CASE(icase)
      CASE(1)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     I I I I  Type 1 integrals ICASE=1.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
c        WRITE(IWRITE,FMT='(/,10X,''Writing out ordered IIII integ.'')')
C
C---- Begin loop over all symmetries.
C
C     If there are no basis functions in a symmetry then we
C     may skip the details of that symmetry.
C
C     As we pass through the list of ordered integrals
C     we keep a counter updated.
C
C     IOFFI is the offset from the start of the array
C     of fully ordered integrals into the present record
C     of the present symmetry. All records written are the
C     same length for any given symmetry.
C
         IOFFI=0
C
         DO I=1, NSYMT
            IF(NBFT(I).EQ.0)CYCLE
            WRITE(IUNIT)MINUS1, I, I, I, I, izero, izero, izero
C
C....... Each symmetry has nbft(i)*(nbft(i)+1)/2 records
C        each of length nbft(i)*(nbft(i)+1)/2 integrals
C        which have to be written. These are ordered one
C        after another in core.
C
            IJMX=IJPAIR(I,I)
            DO IJ=1, IJMX
               DO KL=1, IJMX
                  BUF(KL)=A(IOFFI+KL)
               END DO
               IOFFI=IOFFI+IJMX
C
               IXX=IJMX
C
               IF(IXX.NE.0)THEN
                  IF(IPFLG.EQ.2)CALL PRINTR(IJ,IXX,BUF)
                  WRITE(IUNIT)IJ, IXX, (IZERO,II=1,6)
                  WRITE(IUNIT)(BUF(II),II=1,IXX), 
     &                        (xzero,ii,ii=ixx+1,8-ixx)
               END IF
            END DO
         END DO
C
      CASE(2)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     I J I J  Type 2 integrals ICASE=2.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
c        WRITE(IWRITE,FMT='(/,10X,''Writing out ordered IJIJ integ.'')')
C
C---- Begin double loop over lower triangle of symmetry pairs.
C
C     If there are no basis functions in a symmetry then we
C     may skip the details of that symmetry.
C
C     As we pass through the list of ordered integrals
C     we keep a counter updated.
C
C     IOFFI is the offset from the start of the array
C     of fully ordered integrals into the present record
C     of the present symmetry. All records written are the
C     same length for any given symmetry.
C
         IOFFI=0
C
         DO I=2, NSYMT
            IF(NBFT(I).EQ.0)CYCLE
            DO J=1, I-1
               IF(NBFT(J).EQ.0)CYCLE
               WRITE(IUNIT)MINUS2, I, J, I, J, izero, izero, izero
C
C.......... Each symmetry has nbft(i)*nbft(j) records
C           each of length nbft(i)*nbft(j) integrals
C           which have to be written. These are ordered
C           one after another in core.
C
               IJMX=IJPAIR(I,J)
               DO IJ=1, IJMX
                  DO KL=1, IJMX
                     BUF(KL)=A(IOFFI+KL)
                  END DO
                  IOFFI=IOFFI+IJMX
                  IXX=IJMX
C
                  IF(IXX.NE.0)THEN
                     IF(IPFLG.EQ.2)CALL PRINTR(IJ,IXX,BUF)
                     WRITE(IUNIT)IJ, IXX, (IZERO,II=1,6)
                     WRITE(IUNIT)(BUF(II),II=1,IXX), 
     &                           (xzero,ii,ii=ixx+1,8-ixx)
                  END IF
               END DO
            END DO
         END DO
c
      CASE(3)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     I I J J  Type 3 integrals ICASE=3.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
c        WRITE(IWRITE,FMT='(/,10X,''Writing out ordered IIJJ integ.'')')
C
C---- Begin double loop over lower triangle of symmetry pairs.
C
C     If there are no basis functions in a symmetry then we
C     may skip the details of that symmetry.
C
C     As we pass through the list of ordered integrals
C     we keep a counter updated.
C
C     IOFFI is the offset from the start of the array
C     of fully ordered integrals into the present record
C     of the present symmetry. All records written are the
C     same length for any given symmetry.
C
         IOFFI=0
C
         DO I=2, NSYMT
            IF(NBFT(I).EQ.0)CYCLE
            DO J=1, I-1
               IF(NBFT(J).EQ.0)CYCLE
               WRITE(IUNIT)MINUS3, I, I, J, J, izero, izero, izero
c      WRITE(IWRITE,3050) MINUS3,I,I,J,J,izero,izero,izero
C
C.......... Note index ordering here !!!!
C
C           We have the integrals in the order JJII so
C           that we have nbft(j)*(nbft(j)+1)/2 records
C           each holding nbft(i)*(nbft(i)+1)/2 integrals.
C
               IJMX=IJPAIR(J,J)
               KLMX=IJPAIR(I,I)
               DO IJ=1, IJMX
                  DO KL=1, KLMX
                     BUF(KL)=A(IOFFI+KL)
                  END DO
                  IOFFI=IOFFI+KLMX
                  IXX=KLMX
                  IF(IXX.NE.0)THEN
                     IF(IPFLG.EQ.2)CALL PRINTR(IJ,IXX,BUF)
                     WRITE(IUNIT)IJ, IXX, (IZERO,II=1,8)
                     WRITE(IUNIT)(BUF(II),II=1,IXX), 
     &                           (xzero,ii,ii=ixx+1,8-ixx)
                  END IF
               END DO
            END DO
         END DO
C
      CASE(4)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C     I J K L  Type 6 integrals ICASE=4.
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
c        WRITE(IWRITE,FMT='(/,10X,''Writing out ordered IJKL integ.'')')
C
C---- We descend into a loop over all four symmetries. Of course,
C     we apply the usual canonical scheme
C
C         J < I ; L < K  and IJ < KL
C
C
C     As before IOFFI is a bias counter.
C
         IOFFI=0
C
         DO I=2, NSYMT
            IF(NBFT(I).EQ.0)CYCLE
            DO J=1, I-1
               IF(NBFT(J).EQ.0)CYCLE
               DO K=2, I
                  IF(NBFT(K).EQ.0)CYCLE
                  LEND=K-1
                  IF(K.EQ.I)LEND=J-1
                  DO L=1, LEND
                     IF(NBFT(L).EQ.0)CYCLE
C
C....... We now apply group theory to pick out only the integral
C        blocks that are valid.
C
                     IJTST=IPD2H(I,J)
                     KLTST=IPD2H(K,L)
                     IJKLT=IPD2H(IJTST,KLTST)
                     IF(IJKLT.NE.1)CYCLE
C
                     WRITE(IUNIT)MINUS6, I, J, K, L, izero, izero, izero
                     KLMX=IJPAIR(I,J) ! JMC are these 2 lines correct? 
                     IJMX=IJPAIR(K,L)
C
                     DO IJ=1, IJMX
                        DO KL=1, KLMX
                           BUF(KL)=A(IOFFI+KL)
                        END DO
                        IOFFI=IOFFI+KLMX
C
                        IXX=KLMX
C
                        IF(IXX.NE.0)THEN
                           IF(IPFLG.EQ.2)CALL PRINTR(IJ,IXX,BUF)
                           WRITE(IUNIT)IJ, IXX, (IZERO,II=1,6)
                           WRITE(IUNIT)(BUF(II),II=1,IXX), 
     &                                 (xzero,ii,ii=ixx+1,8-ixx)
                        END IF
                     END DO
                  END DO
               END DO
            END DO
         END DO
C
      END SELECT
      RETURN
C
C---- Format Statements
C
 3010 FORMAT(/,10X,'Integrals sequence number = ',I5,//,10X,
     &       'Symmetry numbers = ',4(I3,1X),/,10X,'Function numbers = ',
     &       4(I3,1X),/,10X,'Integral value   = ',F20.10,/,10X,
     &       'Triangulation    = ',I3,/,10X,'Offset storage   = ',I10,/)
C
 300  FORMAT(' TYPE=',I4,' SYMS=',4I4)
C
 3050 FORMAT(' Integral block type = ',I4,' Symmetry labels =',4I4)
 3055 FORMAT(' Number of records in ordered set = ',I10,/,
     &       ' Integrals per record             = ',I10,/)
C
      END SUBROUTINE SRT2EL
!*==printr.spg  processed by SPAG 6.56Rc at 11:57 on 13 May 2010
      SUBROUTINE PRINTR(IJ,IXX,BUF)
C***********************************************************************
C
C     Utility to unpack and print a vector BUF of length IXX
C
C     In this code BUF will contain ordered integrals. There will
C     be IXX integrals in the buffer corresponding to record number
C     IJ within a symmetry block.
C
C     Author: Charles J Gillan
C
C     Copyright 1995 (c) Charles J Gillan
C     All rights reserved
C
C***********************************************************************
      USE precisn, ONLY : wp ! for specifying the kind of reals                        
      USE sword_data, ONLY : ILINE, IWRITE
      IMPLICIT NONE
C
C*** Start of declarations rewritten by SPAG
C
C Dummy arguments
C
      INTEGER :: IJ, IXX
      REAL(KIND=wp), DIMENSION(*) :: BUF
      INTENT (IN) BUF, IJ, IXX
C
C Local variables
C
      INTEGER :: I, J, KP
      INTEGER, DIMENSION(ILINE) :: IX
      REAL(KIND=wp), DIMENSION(ILINE) :: X
C
C*** End of declarations rewritten by SPAG
C
C---- Banner header
C
      WRITE(IWRITE,1000)IJ, IXX
C
C---- Loop over the vector in steps of six at a time
C
      KP=0
C
      DO I=1, IXX
         KP=KP+1
         X(KP)=BUF(I)
         IX(KP)=I
         IF(KP.EQ.ILINE)THEN
            WRITE(IWRITE,1020)(IX(J),X(J),J=1,KP)
            KP=0
         END IF
      END DO
C
C---- If the array length is a multiple of six then finished else
C     print the remaining few elements - DO control takes care
C     of this.
C
      WRITE(IWRITE,1020)(IX(J),X(J),J=1,KP)
C
      RETURN
C
C---- Format Statements
C
 1000 FORMAT(/,5X,'Printed record of integrals: ',//,5X,
     &       '  Sequence number (IJ) = ',I7,
     &       '  No. of non-zero integrals (IXX) = ',I7,/)
 1020 FORMAT(6(1X,I5,F15.8))
C
 9900 FORMAT(/,10X,'**** Error in PRINTR: ',/)
C
      END SUBROUTINE PRINTR
!*==normal.spg  processed by SPAG 6.56Rc at 11:57 on 13 May 2010
      SUBROUTINE NORMAL(NSYM,NBFT,XNORM,IWRITE,NORDINTS,LTRI,ITAPE,
     &                  NFTTAIL,ZTAIL)
C***********************************************************************
C
C     NORMAL - Computes normalization factors for the symmetry adapted
C              basis functions in each integral.
C
C***********************************************************************
      USE precisn, ONLY : wp ! for specifying the kind of reals                        
      USE CONSTS, ONLY : XZERO, XONE
      USE params, ONLY : COVERLAP, MAXSYM
      USE global_utils, ONLY : search
      USE sword_data, ONLY : VSMALL
      IMPLICIT NONE
C
C*** Start of declarations rewritten by SPAG
C
C Dummy arguments
C
      INTEGER :: ITAPE, IWRITE, LTRI, NFTTAIL, NORDINTS, NSYM
      LOGICAL :: ZTAIL
      INTEGER, DIMENSION(*) :: NBFT
      REAL(KIND=wp), DIMENSION(MAXSYM,*) :: XNORM
      INTENT (IN) IWRITE, LTRI, NBFT, NORDINTS, NSYM, ZTAIL
      INTENT (INOUT) XNORM
C
C Local variables
C
      INTEGER :: I, IFAIL, IFUNCT, INT, J, L, NUMB, KK
      INTEGER, DIMENSION(ltri) :: INDEXV
      REAL(KIND=wp), DIMENSION(nordints) :: OLAP
      REAL(KIND=wp), DIMENSION(ltri) :: XBUF
C
C*** End of declarations rewritten by SPAG
C
C---- Debug banner header
C
c        WRITE(IWRITE,1000)
c        WRITE(IWRITE,1010) NSYM,(L,NBFT(L),L=1,NSYM)
c        WRITE(IWRITE,1020) NORDINTS,LTRI,ITAPE
C
C---- Initialization of ordered integral buffer
C
      DO I=1, NORDINTS
         OLAP(I)=XZERO
      END DO
C
C---- Find the overlap integrals on the dataset written by MOLECULE.
C     Then read and order these.
C
      CALL SEARCH(ITAPE,COVERLAP,ifail,iwrite)
      IF(ifail.NE.0)STOP 999
C JMC replacing GO TO structure with the following do loop
      DO KK = 1, HUGE(1) ! make sure this loop is technically not infinite
         READ(ITAPE)(XBUF(I),I=1,LTRI), (INDEXV(I),I=1,LTRI), NUMB
C
         IF(NUMB.LT.0)EXIT ! i.e. go to the IF(ZTAIL)THEN line
         DO I=1, NUMB
            OLAP(INDEXV(I))=XBUF(I)
         END DO
      END DO
C
C---- Find the overlap tail integrals and do the subtraction sum if
C     necessary.
C
      IF(ZTAIL)THEN
         REWIND nfttail
         CALL SEARCH(NFTTAIL,COVERLAP,ifail,iwrite)
         IF(ifail.NE.0)STOP 999
C JMC replacing GO TO structure with the following do loop
         DO KK = 1, HUGE(1) ! make sure this loop is technically not infinite
            READ(NFTTAIL)(XBUF(I),I=1,LTRI), (INDEXV(I),I=1,LTRI), NUMB
            IF(NUMB.LT.0)EXIT ! i.e. go to the INT=0 line below
            DO I=1, NUMB
               OLAP(INDEXV(I))=OLAP(INDEXV(I))-XBUF(I)
            END DO
         END DO
      END IF
C
C---- Begin loop over the ordered lower half triangles of integrals
C
C     N.B. this is an assumption about the input !
C
C     We copy out the overlap matrix elements at this stage.
C
      INT=0
      DO L=1, NSYM
         IFUNCT=0
         DO I=1, NBFT(L)
            DO J=1, I
               INT=INT+1
               IF(I.EQ.J)THEN
                  IFUNCT=IFUNCT+1
                  XNORM(L,IFUNCT)=OLAP(INT)
                  IF(OLAP(INT).LT.VSMALL)THEN
                     WRITE(IWRITE,9900)
                     WRITE(IWRITE,9910)L, I, OLAP(INT)
                     STOP 999
                  END IF
               END IF
            END DO
         END DO
      END DO
C
C---- Compute the normalization factors
C
      DO L=1, NSYM
         DO I=1, NBFT(L)
            XNORM(L,I)=XONE/SQRT(XNORM(L,I))
         END DO
      END DO
C
c        WRITE(IWRITE,2000)
c        DO L=1,NSYM
c           WRITE(IWRITE,2010) L
c           WRITE(IWRITE,2020) (I,XNORM(L,I),I=1,NBFT(L))
c        END DO
C
      RETURN
C
C---- Format Statements
C
 1000 FORMAT(/,10X,'====> NORMAL - COMPUTE FACTORS <====',//)
 1010 FORMAT(10X,'No. of symmetries is orbital set (NSYM) = ',I3,/,10X,
     &       'Basis functions per symmetry : ',//,(10X,I3,1X,I3))
 1020 FORMAT(/,10X,'Length of ordered integral vector = ',I6,/,10X,
     &       'Number of ints per record (LTRI)  = ',I6,/,10X,
     &       'Unit for overlap ints     (ITAPE) = ',I6,/)
 2000 FORMAT(//,10X,'Normalization factors follow : ',/)
 2010 FORMAT(/,10X,'Symmetry number = ',I3,/)
 2020 FORMAT(10X,I3,1X,F20.10)
 9900 FORMAT(/,10X,'**** Error in NORMAL : ',//)
 9910 FORMAT(/,10X,'One of the sym. adapted basis function ',/,10X,
     &       'overlap integrals is too small. This   ',/,10X,
     &       'is probably an error ',//,10X,'Symmetry = ',I3,1X,
     &       'Function no. = ',I3,/,10X,'Overlap integral = ',F20.10,/)
 9935 FORMAT(/,10X,'Error in an integral index : ',//,10X,
     &       'Within record sequence = ',I5,/,10X,
     &       'INDEXV value is        = ',I10,/,10X,
     &       'Integral itself        = ',F25.13,/)
C
      END SUBROUTINE NORMAL
!*==ampere.spg  processed by SPAG 6.56Rc at 11:57 on 13 May 2010
      SUBROUTINE AMPERE(NFTTAIL,ITAPE,IWRITE,NBFT,NSYM,XNORM)
C***********************************************************************
C
C     Normalizes the basis function amplitudes - raw amplitudes -
C     which are later used to build the R-matrix. The normalized
C     amplitudes are placed on the file of ordered integrals
C
C***********************************************************************
      USE precisn, ONLY : wp ! for specifying the kind of reals                        
      USE params, ONLY : CRAWB, C8STARS, CBLANK, MAXSYM,
     &                   CRAWB_PT1, CRAWB_PT2, CAMP_PTS ! LRM May 2011
      USE global_utils, ONLY : search
      IMPLICIT NONE
C
C*** Start of declarations rewritten by SPAG
C
C Dummy arguments
C
      INTEGER :: ITAPE, IWRITE, NFTTAIL, NSYM
      INTEGER, DIMENSION(*) :: NBFT
      REAL(KIND=wp), DIMENSION(MAXSYM,*) :: XNORM
      INTENT (IN) ITAPE, IWRITE, NBFT, NSYM, XNORM
C
C Local variables
C
      REAL(KIND=wp), ALLOCATABLE, DIMENSION(:) :: BUF
      INTEGER :: ICOUNT, IFAIL, ISYM, JFUNCT, NSYMABF
C
C Change by Laura Moore, May 2011: new variables needed to enable
C calculation of surface amplitudes at a number of radial distances
C BNDRY_AMP_HEADER_NAME contains the headers defined in global_data.f90 
C that are used to find the location of the amplitudes within the input
C and output files.  It is of dimension 3 because only 3 headers have  
C been defined in global_data.f90.
C The code that follows in this subroutine has additional loops over the
C radial distances at which the amplitudes are to be calculated (these 
C loops are indexed by I_RADIAL_PT)     
      INTEGER :: NUM_PTS_NEEDING_BNDRY_AMPS, I_RADIAL_PT
      REAL(KIND=wp) :: DIST_BTN_PTS_NEEDING_BNDRY_AMPS
      CHARACTER(LEN=8), DIMENSION(3) :: BNDRY_AMP_HEADER_NAME
C
C*** End of declarations rewritten by SPAG
C
C---- Compute the number of amplitudes to be read i.e. number of
C     total functions
C
      JFUNCT=0
C
      DO ISYM=1, NSYM
         JFUNCT=JFUNCT+NBFT(ISYM)
      END DO
C
C---- Pick up the amplitudes now from the tail integrals file
C
C Change by Laura Moore, May 2011: to enable read of amplitudes
C at several radial distances.  First read the number of radial 
C distances at which amplitudes are to be calculated and their
C radial spacing.  If header cannot be found then assume default
C values
      CALL SEARCH(NFTTAIL,CAMP_PTS,ifail,iwrite)
      IF(ifail.NE.0)THEN
C Assume default values      
         NUM_PTS_NEEDING_BNDRY_AMPS=1
         DIST_BTN_PTS_NEEDING_BNDRY_AMPS=0.0_wp
         WRITE(IWRITE,8030) NFTTAIL,NUM_PTS_NEEDING_BNDRY_AMPS
      ELSE
         READ(NFTTAIL)NUM_PTS_NEEDING_BNDRY_AMPS
         READ(NFTTAIL)DIST_BTN_PTS_NEEDING_BNDRY_AMPS     
      END IF 
C  An error check on value of NUM_PTS_NEEDING_BNDRY_AMPS
      IF ((NUM_PTS_NEEDING_BNDRY_AMPS.GT.3).OR. 
     &    (NUM_PTS_NEEDING_BNDRY_AMPS.LT.1)) THEN 
         WRITE(IWRITE,9900)
         WRITE(IWRITE,9920) NUM_PTS_NEEDING_BNDRY_AMPS,NFTTAIL
         STOP 9920
      END IF
C
C LRM: Write the number of points at which the surface
C     amplitudes are to be calculated, and the spacing
C     between these points so that these variables can  
C     be read by subsequent programs
      WRITE(ITAPE)C8STARS, CBLANK, CBLANK, CAMP_PTS
      WRITE(ITAPE)NUM_PTS_NEEDING_BNDRY_AMPS, CBLANK, 
     &            CBLANK, CBLANK, CBLANK
      WRITE(ITAPE)DIST_BTN_PTS_NEEDING_BNDRY_AMPS, CBLANK, 
     &            CBLANK, CBLANK, CBLANK
C 
C LRM: Read nsymabf and then allocate buf once    
      CALL SEARCH(NFTTAIL,CRAWB,ifail,iwrite)
      IF(ifail.NE.0)STOP 999
      READ(NFTTAIL)NSYMABF
      ALLOCATE(buf(nsymabf))
C---- Got to make sure that what we have matches the data on file
C
      IF(NSYMABF.NE.JFUNCT)THEN
         WRITE(IWRITE,9900)
         WRITE(IWRITE,9910)JFUNCT, NSYMABF
         STOP 999
      END IF
C
C LRM: Now normalize the amplitudes at each radial distance
C
C     First set up vector of Bndry_Amp_Header_Names used to find the  
C     location of amplitudes at each radial distance within the input
C     and output files
      BNDRY_AMP_HEADER_NAME(1) = CRAWB
      BNDRY_AMP_HEADER_NAME(2) = CRAWB_PT1
      BNDRY_AMP_HEADER_NAME(3) = CRAWB_PT2
C
      DO I_RADIAL_PT = 1, NUM_PTS_NEEDING_BNDRY_AMPS
C
         CALL SEARCH(NFTTAIL,BNDRY_AMP_HEADER_NAME(I_RADIAL_PT),
     &               ifail,iwrite)
         IF(ifail.NE.0)STOP 999
C
         READ(NFTTAIL)NSYMABF
         READ(NFTTAIL)(BUF(ISYM),ISYM=1,NSYMABF)
C
C---- Loop over all symmetries and normalize each amplitude
C
C     We exploit the fact that the amplitudes are ordered already
C
         ICOUNT=0
C
         DO ISYM=1, NSYM
            DO JFUNCT=1, NBFT(ISYM)
               ICOUNT=ICOUNT+1
c               WRITE(IWRITE,2010) ISYM,JFUNCT,BUF(ICOUNT)
               BUF(ICOUNT)=BUF(ICOUNT)*XNORM(ISYM,JFUNCT)
c               WRITE(IWRITE,2020) XNORM(ISYM,JFUNCT),BUF(ICOUNT)
            END DO
         END DO
C
C---- Write the new amplitudes on the file of ordered integrals
C
         WRITE(ITAPE)C8STARS, CBLANK, CBLANK,
     &               BNDRY_AMP_HEADER_NAME(I_RADIAL_PT)
         WRITE(ITAPE)NSYMABF, CBLANK, CBLANK, CBLANK, CBLANK
         WRITE(ITAPE)(BUF(JFUNCT),JFUNCT=1,ICOUNT), CBLANK, CBLANK, 
     &                CBLANK, CBLANK
C
      END DO
C
      DEALLOCATE(buf)
      RETURN
C
C---- Format Statements
C
 2010 FORMAT(/,10X,'Symmetry number = ',I2,' FunctionNo. = ',I3,/,10X,
     &       'Un-normalized amplitude = ',F25.13)
 2020 FORMAT(10X,'Normalization factor    = ',F25.13,/,10X,
     &       'Normalized amplitude    = ',F25.13)
 8030 FORMAT(10X,'**********  WARNING IN AMPERE ********** : ',/,
     &       10X,' Header CAMP_PTS (defined in global_data.f90) not',
     &       /,10X,' found in unit NFTTAIL.',/,10X,' Unit NFTTAIL ',
     &       'is unit ',I3,/,10X,' Continue using ',
     &       ' the following default value: ',/,10X,' Number of ',
     &       ' radial points at which surface amplitudes ',/,10X,
     &       ' are needed = ',I3,/,10X, ' i.e. only calculate surface',
     &       ' amplitudes on the R-matrix ',/,10X,' boundary.  This',
     &       ' is OK for all runs unless you want to ',/,10X, 
     &       ' generate input for subsequent R-matrix incorporating',
     &       /,10X,' time (RMT) calculations',/,10X,
     &       ' ************************************************ ',/)
 9900 FORMAT(/,10X,'Error in AMPERE : ',/)
 9910 FORMAT(10X,'Computed number of amplitudes (JFUNCT) = ',I10,/,10X,
     &       'Number of amplitudes on file (NSYMABF) = ',I10,/)
 9920 FORMAT(10X,' Number of radial points at which surface ',
     &       'amplitudes',/,10X,' are needed = ',I3,
     &       ' but NUM_PTS_NEEDING_BNDRY_AMPS',
     &       /,10X,' should be in the range 1 to 3.',/,10X,' This is',
     &       ' because header names for output data on unit ',
     &       /,10X,' NFTTAIL (unit ',I3,') are only defined in ',/,10X, 
     &       ' global_data.f90 for up to 3 radial points.',/,10X,
     &       ' The default setting is 1 to calculate amplitudes only ',
     &       /,10X,' on the R-matrix boundary - OK for all runs ',/,
     &       10X,' unless you want to generate input for subsequent',
     &       /,10X,' R-matrix incorporating time (RMT) calculations',
     &       /,10X,' in which case NUM_PTS_NEEDING_BNDRY_AMPS should',
     &       /,10X,' be set to 3 in the GAUSTAIL namelist.',/)
C
      END SUBROUTINE AMPERE
C     
C     Write svn info for file to unit
C
      SUBROUTINE WRITE_SVN_INFO(unit)
      IMPLICIT NONE
      integer, INTENT(IN) :: unit
C     
C     Please do not edit text between pair of $'s
C     subversion automatically edits these entries.
C     
C     The double colon (::) specifies fixed width which 
C     is important to work with fortran fixed format files.
C     
C     The format statement has been carefully constructed to 
C     give maximum width between $$'s
C     
C     If edits are not occuring then enable svn keywords with:
C     svn propset svn:keywords "URL Author Date Rev Id" thisfile.f
C     
 9876 FORMAT(/,1X,                                                     '
     &$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
     &'/,1X,'$',64X,'$',/,
     &     1X,'$',16X,'Subversion Revision Information',17X,'$',/,
     &     1X,'$',64X,'$',/,1X,                                        '
     &$Id:: sword.f 730 2011-10-06 15:25:34Z laura                     $
     &',/,1X,'$',64X,'$',/,1X,                                         '
     &$URL:: http://ccpforge.cse.rl.ac.uk/svn/ukrmol-in/tags/release-1#$
     &',/,1X,                                                          '
     &$LastChangedBy:: laura                                           $
     &',/,1X,                                                          '
     &$LastChangedDate:: 2011-10-06 16:25:34 +0100 (Thu, 06 Oct 2011)  $
     &',/,1X,                                                          '
     &$LastChangedRevision:: 730                                       $
     &',/,1X,'$',64X,'$',/,1X,                                         '
     &$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
     &',/)
C      
      WRITE(unit,9876)
C     
      END SUBROUTINE WRITE_SVN_INFO
