C  5 Jul 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  9 dec 03 - mws - synch common block runopt
C 13 Jun 01 - DGF - some code relocated to new source module solib
C 25 Mar 00 - dgf - readcsf: fix spin projection quantities
C 16 feb 00 - mws - mimeit: argument of integer type treated correctly
C 21 dec 99 - dgf - changes from 7/jan/98 to 4/sep/98 brought online
C  4 sep 98 - dgf - add dynamic load balancing for FF calcs (CIMAT)
C  7 aug 98 - dgf - new parallel: split recsize instead of the # of rec
C  5 jul 98 - dgf - implement a 2nd strategy for FF updating of large AS
C  5 jul 98 - dgf - deport some critical bit routines to unport
C  4 jul 98 - dgf - purge the core-active form-factor ballast
C 17 jan 98 - dgf - generalise code for arbitrary spin multiplicities
C                   and arbitrary number of CI runs
C  7 jan 98 - dgf - rid spntrn of the bothersome arrays
C  6 jan 98 - dgf - change CSF data order to allow for 64 bit machines
C                   parlst eluminated, sort merged with keptrk
C                   several common blocks (SPINFO,flags) redefined
C                   parallelise keptrk to distribute sorted FF's storage
C                   remove spin-function generation code
C 14 oct 97 - dgf,trf - many changes to allow for 10 active orbitals
C                   interface Gamess CSFs directly with SOC code
C 16 JUL 97 - MWS - CFSORT,CIMAT,CSFGEN,GETCSF,SORT,SPCSF: MOVE DECLARES
C  2 APR 97 - DGF - ADAPT TOM FURLANI'S FORM FACTOR CODE FOR GAMESS
C
C*MODULE SOFFAC  *DECK A2B2
      SUBROUTINE A2B2(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE GENERATES AND STORES THE TWO PARTICLE FORMFACTOR
C *** WHEN THERE ARE TWO DISCOS - ONE OF EACH SPIN;
C             <A(S),B(S)//A(T),B(T)>
C *** A - ALPHA, B - TRIPLET; S - SINGLET, T - TRIPLET
C
C *** PICK UP THE DISCOINCIDENT ORBITALS IN THE SINGLET AND TRIPLET
C     DETERMINANTS
C
      NORB1U = IAND(MSKNAO,NDIS1)
      NORB2U = IAND(MSKNAO,NDIS2)
C
C *** DETERMINE PHASE
C
      CALL PHASE2(NORB1U,NORB2U,NCOIN,KOUNT1,NAO)
C
C *** FIND ORBITAL NUMBER OF DISCOINCIDENT BETA ORBITALS
C
      LOC1U=IFINDEM(NORB1U)
      LOC2U=IFINDEM(NORB2U)
C
C *** NOW GET THE ALPHA CONTRIBUTION
C
      NSDET1 = ISHFT(NDIS1,-NAO)
      NSDET2 = ISHFT(NDIS2,-NAO)
      NORB1L = IAND(MSKNAO,NSDET1)
      NORB2L = IAND(MSKNAO,NSDET2)
      LOC1L=IFINDEM(NORB1L)
      LOC2L=IFINDEM(NORB2L)
C
C *** PLACE ALPHA ORBITALS IN THEIR PROPER POSITION IN ORDER TO PICKUP
C *** COINCIDENT ORBITALS
C
      NORB1L = ISHFT(NORB1L,NAO)
      NORB2L = ISHFT(NORB2L,NAO)
      CALL PHASE2(NORB1L,NORB2L,NCOIN,KOUNT2,NAOD)
C
C *** DETERMINE PHASE CHANGE, IF KOUNT IS ODD THEN CHANGE SIGN
C
      KOUNT = KOUNT1+KOUNT2
      COEF=CPROD
      IF (MOD(KOUNT,2) .EQ. 1) COEF = -COEF
      COEF = -HALF*COEF
C
C *** FORM THE FOUR INDICES FOR THE COULOMB INTEGRALS
C *** <LOC1L,LOC1U/L12/LOC2L,LOC2U> AND <LOC1U,LOC1L/L12/LOC2U,LOC2L>
C
C *** STORE THE FORMFACTOR FOR IJKL
C *** STORE THE FORMFACTOR FOR JILK
C
      CALL STOR4DX(LOC1L,LOC1U,LOC2L,LOC2U,COEF,-COEF,D2Z,JFLAGZ)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK AAAB
      SUBROUTINE AAAB(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE CALCULATES FORM FACTORS FOR DETERMINANTS WHICH
C *** HAVE TWO DISCOS AND THE DISCOS HAVE THE FOLLOWING SPIN
C ***             <A(S),A(S)//A(T),B(T)>
C
C *** PICK UP THE DISCOINCIDENT ORBITALS, TRIPLET FIRST
C
      NORB2U = IAND(MSKNAO,NDIS2)
      NORB2L = IAND(MSKNAO,ISHFT(NDIS2,-NAO))
      LOC2L=IFINDEM(NORB2L)
      LOC2U=IFINDEM(NORB2U)
C
C *** NOW SINGLET
C
      NTORB = IAND(MSKNAO,ISHFT(NDIS1,-NAO))
      CALL FIND2(NTORB,LOC1U,LOC1L)
C
C *** PUT 1 BIT IN PROPER LOCATION
C
      NORB2L = ISHFT(NORB2L,NAO)
      NORB1L = ISHFT(1,(NAOD-LOC1L))
      NORB1U = ISHFT(1,(NAOD-LOC1U))
C
C *** CALCULATE THE PHASE CHANGE
C
      CALL PHASE2(NORB1L,NORB2L,NCOIN,KOUNT1,NAOD)
      CALL PHASE2(NORB1U,NORB2U,NCOIN,KOUNT2,NAOD)
      KOUNT = KOUNT1 + KOUNT2
      COEF = CPROD
      IF (MOD(KOUNT,2) .EQ. 1) COEF = -COEF
      COEF1 = HALF*COEF
C
C *** FORM THE FOUR INDICES FOR THE COULOMB INTEGRALS
C *** <LOC1L,LOC1U/L12/LOC2L,LOC2U> AND <LOC1U,LOC1L/L12/LOC2U,LOC2L>
C
C *** STORE THE FORM FACTORS FOR THE COULOMB INTEGRALS IJKL,JILK
C
      CALL STOR4DX(LOC1L,LOC1U,LOC2L,LOC2U,COEF,COEF1,D2Z,JFLAGZ)
C
C *** NOW DO THE EXCHANGE TERMS <1L,1U/L12/2U,2L> AND <1U,1L/L12/2L,2U>
C     IJLK,JIKL
C
      COEF = -COEF
      COEF1 = HALF*COEF
      CALL STOR4DX(LOC1L,LOC1U,LOC2U,LOC2L,COEF1,COEF,D2Z,JFLAGZ)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ABAA
      SUBROUTINE ABAA(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE CALCULATES FORM FACTORS FOR DETERMINANTS WHICH
C *** HAVE TWO DISCOS AND THE DISCOS HAVE THE FOLLOWING SPIN
C ***             <A(S),B(S)//A(T),A(T)>
C
C *** PICK UP THE DISCOINCIDENT ORBITALS, SINGLET FIRST
C
      NORB1U = IAND(MSKNAO,NDIS1)
      NORB1L = IAND(MSKNAO,ISHFT(NDIS1,-NAO))
      LOC1L=IFINDEM(NORB1L)
      LOC1U=IFINDEM(NORB1U)
C
C *** NOW TRIPLET
C
      NTORB = IAND(MSKNAO,ISHFT(NDIS2,-NAO))
      CALL FIND2(NTORB,LOC2U,LOC2L)
C
C *** PUT 1 BIT IN PROPER LOCATION
C
      NORB1L = ISHFT(NORB1L,NAO)
      NORB2L = ISHFT(1,(NAOD-LOC2L))
      NORB2U = ISHFT(1,(NAOD-LOC2U))
C
C *** CALCULATE THE PHASE CHANGE
C
      CALL PHASE2(NORB1L,NORB2L,NCOIN,KOUNT1,NAOD)
      CALL PHASE2(NORB1U,NORB2U,NCOIN,KOUNT2,NAOD)
      KOUNT = KOUNT1 + KOUNT2
C
C *** FORM THE FOUR INDICES FOR THE COULOMB INTEGRALS
C *** <LOC1L,LOC1U/L12/LOC2L,LOC2U> AND <LOC1U,LOC1L/L12/LOC2U,LOC2L>
C
C *** STORE THE FORM FACTORS FOR THE COULOMB INTEGRALS IJKL,JILK
C
      COEF = CPROD
      IF (MOD(KOUNT,2) .EQ. 1) COEF = -COEF
      COEF1 = HALF*COEF
      CALL STOR4DX(LOC1L,LOC1U,LOC2L,LOC2U,COEF,COEF1,D2Z,JFLAGZ)
C
C *** NOW DO THE EXCHANGE TERMS <1L,1U/L12/2U,2L> AND <1U,1L/L12/2L,2U>
C     IJLK,JIKL
C
      COEF = -COEF
      COEF1 = HALF*COEF
      CALL STOR4DX(LOC1L,LOC1U,LOC2U,LOC2L,COEF,COEF1,D2Z,JFLAGZ)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ABBB
      SUBROUTINE ABBB(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE CALCULATES FORM FACTORS FOR DETERMINANTS WHICH
C *** HAVE TWO DISCOS AND THE DISCOS HAVE THE FOLLOWING SPIN
C ***             <A(S),B(S)//B(T),B(T)>
C
C *** PICK UP THE DISCOINCIDENT ORBITALS, SINGLET FIRST
C
      NORB1U = IAND(MSKNAO,NDIS1)
      NORB1L = IAND(MSKNAO,ISHFT(NDIS1,-NAO))
      LOC1L=IFINDEM(NORB1L)
      LOC1U=IFINDEM(NORB1U)
C
C *** NOW TRIPLET
C
      NTORB = IAND(MSKNAO,NDIS2)
      CALL FIND2(NTORB,LOC2U,LOC2L)
      NORB2L = ISHFT(1,(NAO-LOC2L))
      NORB2U = ISHFT(1,(NAO-LOC2U))
      NORB1L = ISHFT(NORB1L,NAO)
C
C *** CALCULATE THE PHASE CHANGE
C
      CALL PHASE2(NORB1U,NORB2U,NCOIN,KOUNT1,NAO)
      CALL PHASE2(NORB1L,NORB2L,NCOIN,KOUNT2,NAOD)
      KOUNT = KOUNT1 + KOUNT2
C
C *** FORM THE FOUR INDICES FOR THE COULOMB INTEGRALS
C *** <LOC1L,LOC1U/L12/LOC2L,LOC2U> AND <LOC1U,LOC1L/L12/LOC2U,LOC2L>
C
C *** STORE THE FORM FACTORS FOR THE COULOMB INTEGRALS IJKL,JILK
C
      COEF = CPROD
      IF (MOD(KOUNT,2) .EQ. 1) COEF = -COEF
      COEF1 = HALF*COEF
      CALL STOR4DX(LOC1L,LOC1U,LOC2L,LOC2U,COEF1,COEF,D2Z,JFLAGZ)
C
C *** NOW DO THE EXCHANGE TERMS <1L,1U/L12/2U,2L> AND <1U,1L/L12/2L,2U>
C     IJLK,JIKL
C
      COEF = -COEF
      COEF1 = HALF*COEF
      CALL STOR4DX(LOC1L,LOC1U,LOC2U,LOC2L,COEF1,COEF,D2Z,JFLAGZ)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ALLSME
      SUBROUTINE ALLSME(ALLB,NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(THREE2=1.5D+00)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*)
      LOGICAL ALLB
C
C *** THIS ROUTINE CALCULATES FORM FACTORS FOR DETERMINANTS WHICH
C *** HAVE TWO DISCOS AND ALL DISCOS HAVE THE SAME SPIN
C *** <A(S),A(S)//A(T),A(T)> OR <B(S),B(S)//B(T),B(T)>
C
C *** PICK UP DISCOINCIDENT ORBITALS, IF ALLB = TRUE,
C     THEN ALL BETA SPINS
C
      IF (ALLB) THEN
         NORB1 = IAND(MSKNAO,NDIS1)
         NORB2 = IAND(MSKNAO,NDIS2)
      ELSE
         NORB1 = IAND(MSKNAO,ISHFT(NDIS1,-NAO))
         NORB2 = IAND(MSKNAO,ISHFT(NDIS2,-NAO))
      ENDIF
      CALL FIND2(NORB1,LOC1U,LOC1L)
      CALL FIND2(NORB2,LOC2U,LOC2L)
C
C *** CALCULATE THE PHASE CHANGE
C
      IF (ALLB) THEN
         NUM = NAO
         NORB1L = ISHFT(1,(NUM-LOC1L))
         NORB2L = ISHFT(1,(NUM-LOC2L))
         NORB1U = ISHFT(1,(NUM-LOC1U))
         NORB2U = ISHFT(1,(NUM-LOC2U))
      ELSE
         NUM = NAO + NAO
         NORB1L = ISHFT(1,(NUM-LOC1L))
         NORB2L = ISHFT(1,(NUM-LOC2L))
         NORB1U = ISHFT(1,(NUM-LOC1U))
         NORB2U = ISHFT(1,(NUM-LOC2U))
      ENDIF
      CALL PHASE2(NORB1L,NORB2L,NCOIN,KOUNT1,NUM)
      CALL PHASE2(NORB1U,NORB2U,NCOIN,KOUNT2,NUM)
      KOUNT = KOUNT1 + KOUNT2
      COEF = CPROD
      IF (MOD(KOUNT,2) .EQ. 1) COEF = -COEF
C
C *** FORM FACTOR DEPENDS ON ORBITAL SPINS
C
      IF (ALLB) THEN
         COEF1 = -THREE2*COEF
         COEF2 = -COEF1
      ELSE
         COEF1 = THREE2*COEF
         COEF2 = -COEF1
      ENDIF
C
C *** CALCULATE THE COULOMB FOUR INDICES IJKL AND JILK FOR
C *** <LOC1L,LOC1U/L12/LOC2L,LOC2U> AND <LOC1U,LOC1L/L12/LOC2U,LOC2L>
C
C *** STORE FORM FACTORS, FIRST COULOMB IJKL,JILK
C
      CALL STOR4DX(LOC1L,LOC1U,LOC2L,LOC2U,COEF1,COEF1,D2Z,JFLAGZ)
C
C *** CALCULATE FOUR INDICES FOR THE EXCHANGE TYPE INTEGRALS
C *** <LOC1L,LOC1U/L12/LOC2U,LOC2L> AND <LOC1U,LOC1L/L12/LOC2L,LOC2U>
C
C *** NOW EXCHANGE IJLK,JIKL
C
      CALL STOR4DX(LOC1L,LOC1U,LOC2U,LOC2L,COEF2,COEF2,D2Z,JFLAGZ)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ALPHFF
      SUBROUTINE ALPHFF(NALPHA,KOUNTA,GAMMA,SPNDT1,COEF,LOC1,LOC2,D2Z,
     *           JFLAGZ,LOCC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER GAMMA,SPNDT1
      DIMENSION JFLAGZ(*),LOCC(*),D2Z(*)
C
C *** THIS ROUTINE GENERATES AND STORES 2 PARTICLE FORM FACTORS
C *** FOR THE CASE OF ONE DISCOINCIDENCE AND ALPHA COINCIDENT
C *** ACTIVE ORBITALS. IT IS CALLED FROM SUBROUTINE ONED2P
C
C
C *** FIND THE ORBITAL NUMBERS FOR THE COINCIDENT ALPHA ORBITALS
C
      CALL FNDALL(NALPHA,LOCC,KOUNTA)
C
C *** GET THE FORM FACTORS FOR THE COULOMB AND EXCHANGE INTEGRALS
C
      CALL GETFFA(COEF,C1,C2,C3,C4,GAMMA,SPNDT1)
C
C *** LOOP OVER THE OCCUPIED ALPHA ORBITALS,
C *** L12 IS DEFINED AS THE -GAMMA COMPONENT OF (R12 X P1)/(R12**3)
C *** WHERE R12 IS THE DISTANCE AND DIRECTION BETWEEN ELECTRONS 1 & 2
C *** AND P1 IS THE MOMENTUM OPERATOR FOR ELECTRON 1.
C *** NOTE: (R12**3) REFERS TO MAGNITUDE ONLY
C
      DO 100 K=1,KOUNTA
        I = LOCC(K)
C
C ***   FORM THE FOUR INDICES FOR <I,LOC1//I,LOC2> , <LOC1,I//LOC2,I>
C ***   <I,LOC1//LOC2,I> AND <LOC1,I//I,LOC2>
C ***   LOCC(K) - ORBITAL INDEX OF THE KTH COINCIDENT ALPHA ORBITAL
C
C ***   STORE THE COULOMB FORM FACTORS IJKL,JILK
C
        CALL STOR4DX(I,LOC1,I,LOC2,C1,C2,D2Z,JFLAGZ)
C
C ***   STORE THE EXCHANGE FORM FACTORS, THERE ARE NONE FOR BETA DISCOS
C ***   AND ALPHA COINCIDENT
C       IJLK,JIKL
C
        IF ( (GAMMA .EQ. 0) .AND. (SPNDT1 .EQ. 0) ) THEN
        ELSE
           CALL STOR4DX(I,LOC1,LOC2,I,C3,C4,D2Z,JFLAGZ)
        ENDIF
100   CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK BBAB
      SUBROUTINE BBAB(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE CALCULATES FORM FACTORS FOR DETERMINANTS WHICH
C *** HAVE TWO DISCOS AND THE DISCOS HAVE THE FOLLOWING SPIN
C ***             <B(S),B(S)//A(T),B(T)>
C
C *** PICK UP THE DISCOINCIDENT ORBITALS, TRIPLET FIRST
C
      NORB2U = IAND(MSKNAO,NDIS2)
      NORB2L = IAND(MSKNAO,ISHFT(NDIS2,-NAO))
      LOC2L=IFINDEM(NORB2L)
      LOC2U=IFINDEM(NORB2U)
C
C *** NOW SINGLET
C
      NTORB = IAND(MSKNAO,NDIS1)
      CALL FIND2(NTORB,LOC1U,LOC1L)
      NORB1L = ISHFT(1,(NAO-LOC1L))
      NORB1U = ISHFT(1,(NAO-LOC1U))
      NORB2L = ISHFT(NORB2L,NAO)
C
C *** CALCULATE THE PHASE CHANGE
C
      CALL PHASE2(NORB1U,NORB2U,NCOIN,KOUNT1,NAO)
      CALL PHASE2(NORB1L,NORB2L,NCOIN,KOUNT2,NAOD)
      KOUNT = KOUNT1 + KOUNT2
      COEF = CPROD
      IF (MOD(KOUNT,2) .EQ. 1) COEF = -COEF
      COEF1 = HALF*COEF
C
C *** FORM THE FOUR INDICES FOR THE COULOMB INTEGRALS
C *** <LOC1L,LOC1U/L12/LOC2L,LOC2U> AND <LOC1U,LOC1L/L12/LOC2U,LOC2L>
C
C
C *** STORE THE FORM FACTORS FOR THE COULOMB INTEGRALS IJKL,JILK
C
      CALL STOR4DX(LOC1L,LOC1U,LOC2L,LOC2U,COEF1,COEF,D2Z,JFLAGZ)
C
C *** NOW DO THE EXCHANGE TERMS <1L,1U/L12/2U,2L> AND <1U,1L/L12/2L,2U>
C     IJLK,JIKL
C
      COEF = -COEF
      COEF1 = HALF*COEF
      CALL STOR4DX(LOC1L,LOC1U,LOC2U,LOC2L,COEF,COEF1,D2Z,JFLAGZ)
C
      RETURN
      END
C
C*MODULE SOFFAC  *DECK BETAFF
      SUBROUTINE BETAFF(NBETA,KOUNTB,GAMMA,SPNDT1,COEF,LOC1,LOC2,D2Z,
     *   JFLAGZ,LOCC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER GAMMA,SPNDT1
      DIMENSION JFLAGZ(*),LOCC(*),D2Z(*)
C
C *** THIS ROUTINE GENERATES AND STORES 2 PARTICLE FORM FACTORS
C *** FOR THE CASE OF ONE DISCOINCIDENCE AND BETA COINCIDENT
C *** ACTIVE ORBITALS. IT IS CALLED FROM SUBROUTINE ONED2P
C
C
C *** FIND THE ORBITAL NUMBERS FOR THE COINCIDENT BETA ORBITALS
C
      CALL FNDALL(NBETA,LOCC,KOUNTB)
C
C *** GET THE FORM FACTORS FOR THE COULOMB AND EXCHANGE INTEGRALS
C
      CALL GETFFB(COEF,C1,C2,C3,C4,GAMMA,SPNDT1)
C
C *** LOOP OVER THE OCCUPIED BETA ORBITALS,
C
      DO 100 K=1,KOUNTB
        I = LOCC(K)
C
C ***   FORM THE FOUR INDICES FOR <I,LOC1//I,LOC2> , <LOC1,I//LOC2,I>
C ***   <I,LOC1//LOC2,I> AND <LOC1,I//I,LOC2>
C ***   LOCC(K) - ORBITAL INDEX OF THE KTH COINCIDENT BETA ORBITAL
C
C ***   STORE THE COULOMB FORM FACTORS IJKL,JILK
C
        CALL STOR4DX(I,LOC1,I,LOC2,C1,C2,D2Z,JFLAGZ)
C
C ***   STORE THE EXCHANGE FORM FACTORS, THERE ARE NONE FOR ALPHA DISCOS
C ***   AND BETA COINCIDENT
C       IJLK,JIKL
C
        IF ( (GAMMA .EQ. 0) .AND. (SPNDT1 .EQ. 1) ) THEN
        ELSE
           CALL STOR4DX(I,LOC1,LOC2,I,C3,C4,D2Z,JFLAGZ)
        ENDIF
100   CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK CIMAT
      SUBROUTINE CIMAT(ISODA,IBUFII,IBUFJJ,GAMMA,D1Z,D2Z,IDX1Z,IDX2Z,
     *                 IFLAGZ,JFLAGZ,LOCA,LOCB,FBIN1Z,FBIN2Z,IBINP,
     *                 IBINP1,SAMEMUL,SAVDSK,SLOWFF,ACTION,ICI1,ICI2,
     *                 LDARISO,IBTYP1)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MAXCP=4096)
C
      LOGICAL OUT,OUT1,SAMEMUL,SAVDSK,SLOWFF,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ISODA(*),IBUFII(*),IBUFJJ(*),GAMMA(*),IDX1Z(*),IDX2Z(*),
     *          D1Z(*),D2Z(*),IFLAGZ(*),JFLAGZ(*),LOCA(*),LOCB(*),
     *          FBIN1Z(IBNLEN,*),FBIN2Z(IBNLEN,*),IBINP(*),IBINP1(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      COMMON /SOGUG/  CP(MAXCP),NUNIQ,LIOBP,ISODAF,NSODA,LSTREC,JSODAF,
     *                JSODA,NRECJ
      COMMON /SOBUF/  IBNLEN,IBN2,MAXGAM,MAXFL,MAXFL2,ND1FZ,ND2FZ
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /FMCOM / X(1)
      DATA DEBUG,DBGME/8HDEBUG   ,8HCIMAT   /,USUAL/8HNORMAL  /
C
C *** THIS ROUTINE FORMS THE ONE AND TWO PARTICLE FORM FACTORS FROM THE
C *** SINGLET AND TRIPLET CSFS GENERATED IN SUBROUTINE CSFGEN
C
C     DGF: restrictions: # Slater determinants per CSF 2**10
C     iobp buffer length (roughly # CSFs x NAO) 2**22
C
      MASK12 = 2**12 - 1
      IF(NWDVAR.EQ.1) THEN
         LEXP=52
         LEFCSF=32
      ELSE
         LEXP=20
         LEFCSF=16
      ENDIF
      MASKL=2**LEXP-1
C
      OUT=MASWRK.AND.(NPRINT.EQ.5.OR.EXETYP.EQ.DEBUG.OR.EXETYP.EQ.DBGME)
      IF(OUT) WRITE(IW,*) 'THE MAX # OF UNIQUE COEFS(CSFDET) IS',NUNIQ
C
C *** INITIALIZE VARIABLES
C
      CALL VICLR(IDX2Z,1,NAO4)
      CALL VICLR(IDX1Z,1,NAO2)
      CALL VCLR(D2Z,1,NAO4)
      CALL VCLR(D1Z,1,NAO2)
      CALL VICLR(JFLAGZ,1,MAXFL2)
      CALL VICLR(IFLAGZ,1,MAXFL)
      CALL VICLR(IBINP,1,NAO4+NAO2)
C
C *** GET THE CSFS
C
      CALL READCSF(ICI1,ICI2,X(LIOBP),IBUFII,IBUFJJ,X(JSODA),
     *             SAMEMUL,OUT)
      NOPNSI=IBUFII(1)
      NMATE=0
      NOPNSJ=IBUFJJ(1)
      NXTHWI=1
      KFIG2=1
C
      IF(IBTYP1.EQ.1) CALL DDI_DLBNEXT(NEXT)
      DO 950 KFIG1=1,NOPNSI
        IOLDHW=NXTHWI+1
        IF (KFIG1.LT.NOPNSI) NXTHWI=IBUFII(IOLDHW+1)
        IC1=IOLDHW
        IF(GOPARR) THEN
          IF(IBTYP1.EQ.1) THEN
C           dynamic load balancing
            IF(KFIG1-1.NE.NEXT) GOTO 875
            CALL DDI_DLBNEXT(NEXT)
          ELSE
C           static load balancing
            IF(MOD(KFIG1,NPROC).NE.ME) GOTO 875
          ENDIF
        ENDIF
        NXTHWJ=1
C
        KF2TMP=KFIG2
        DO 850 LL=1,NOPNSJ
          JOLDHW=NXTHWJ+1
          IF (LL.LT.NOPNSJ) NXTHWJ=IBUFJJ(JOLDHW+1)
          IC2=JOLDHW
C
          IC1T=IC1
          IC2T=IC2
          CALL COMPAR(IC1,IC2,IBUFII,IBUFJJ,NDT1,NDT2,MASKL,IZRTST)
          IF(IZRTST.EQ.0) THEN
C
C *** IF compar returns izrtst.eq.0 THEN THERE WERE
C *** LESS THAN 3 DISCOINCIDENCES BETWEEN OPEN SHELL DISTRIBUTIONS
C *** II AND JJ, OTHERWISE izrtst=1 AND THE CI MATRIX ELEMENTS
C *** BETWEEN THE TWO OPEN SHELL DISTRIBUTIONS ARE ZERO
C
            ND1FZ=0
            ND2FZ=0
            NDT2=IAND(IBUFJJ(IC2),MASK12)
C
            CALL SLATER(NDT1,NDT2,IC1,IC2,IBUFII,IBUFJJ,CP,CP,D1Z,D2Z,
     *                  IFLAGZ,JFLAGZ,LOCA,LOCB,MASKL,OUT)
C
            IFZERO=MAX(ND1FZ,ND2FZ)
            IF (IFZERO.NE.0) THEN
              IF(ND1FZ.GT.MAXFL.OR.ND2FZ.GT.MAXFL2) THEN
                WRITE(IW,*) 'SEVERE THUNDERSTORM WARNING',ND1FZ,
     *                      ND2FZ,MAXFL,MAXFL2
                CALL ABRT
              ENDIF
C
              KF12=IOR(ISHFT(KFIG1,LEFCSF),KFIG2)
              NMATE=NMATE+1
              NOFF = 0
              CALL KEPTRK(ISODA,IDX1Z,IFLAGZ,D1Z,ND1FZ,KF12,NOFF,
     *                    GAMMA,FBIN1Z,FBIN1Z,IBINP,NAO2,IBINP1,
     *                    SAVDSK,SLOWFF,LDARISO)
              NOFF = NOFF + NAO2
              CALL KEPTRK(ISODA,IDX2Z,JFLAGZ,D2Z,ND2FZ,KF12,NOFF,
     *                    GAMMA,FBIN2Z,FBIN2Z,IBINP,NAO4,IBINP1,
     *                    SAVDSK,SLOWFF,LDARISO)
C
            ENDIF
            IC2=IC2+NDT2+2
            IF(OUT) WRITE(IW,7776)NDT1,NDT2,IC1,IC2
            IF(OUT) WRITE(IW,4972)KFIG1,KFIG2
C           IF ((KFIG1.EQ.KFIG2).AND.((KFIG1+10).GE.NCONMX)) THEN
C              IF(OUT) WRITE(IW,4975)KFIG1,KFIG2
C           ENDIF
            IC1=IC1+NDT1+2
            NDT1=IAND(IBUFII(IC1),MASK12)
            IF(OUT) WRITE(IW,7776)NDT1,NDT2
            IC2=IC2T
          ENDIF
          IC1=IC1T
          KFIG2=KFIG2+1
  850   CONTINUE
        KFIG2=KF2TMP
  875   CONTINUE
  950 CONTINUE
      IF(IBTYP1.EQ.1) CALL DDI_DLBRESET
C
C *** NOW BE SURE ALL THE BINS ARE EMPTY
C
      IFF1Z=0
      IFF2Z=0
C     call ddi_sync(353)
      NOFF = 0
      CALL EMPTY(ISODA,IBINP,IDX1Z,FBIN1Z,NAO2,NOFF,GAMMA,IBINP1,SAVDSK,
     *           SLOWFF,IFF1Z,LDARISO)
      NOFF = NOFF + NAO2
      CALL EMPTY(ISODA,IBINP,IDX2Z,FBIN2Z,NAO4,NOFF,GAMMA,IBINP1,SAVDSK,
     *           SLOWFF,IFF2Z,LDARISO)
C
C     for savdsk=true empty does an internal check if the number of the
C     sorted FFs is equal to the number of pregenerated. If the code
C     reach this point, the result is positive.
C
      OUT1=NPRINT.EQ.5.OR.EXETYP.EQ.DEBUG.OR.EXETYP.EQ.DBGME
      IF(GOPARR) THEN
         IF(OUT1) THEN
            WRITE(IW,*) 'NODE ',ME,'ACCUMULATED:'
            WRITE(IW,9000) IFF1Z,IFF2Z,IFF1Z+IFF2Z,MS
         ENDIF
         CALL DDI_GSUMI(2302,IFF1Z,1)
         CALL DDI_GSUMI(2303,IFF2Z,1)
      ENDIF
      IF(MASWRK) THEN
         WRITE(IW,8990)
         WRITE(IW,9000) IFF1Z,IFF2Z,IFF1Z+IFF2Z,MS
      ENDIF
      NSODAM=0
      DO 960 I=1,NAO4+NAO2
         IF(IBINP(I).NE.0) NSODAM=NSODAM+1
  960 CONTINUE
C     currently nsodam need not be in the common block
      IF(OUT) WRITE(IW,*) 'THE NUMBER OF ACTUAL INDICES WRITTEN TO DA',
     *                    NSODAM
      IF(OUT1) WRITE(IW,*) 'IBINP ON NODE',ME
      IF(OUT1) WRITE(IW,2349)(IBINP(KKK),KKK=1,NAO4+NAO2)
C
C *** WRITE OUT THE RESULTS
C
      IF(OUT1.AND.((.NOT.SAVDSK).OR.IBINP1(1).GE.0))
     *   CALL OUTFB(ISODA,IBINP,GAMMA)
      IF(OUT) WRITE(IW,7666)NMATE
      MAXGM=-1
      DO 1000 I=1,NAO4+NAO2
         IF(MAXGM.LT.IBINP(I)) MAXGM=IBINP(I)
 1000 CONTINUE
      IF(SAVDSK) THEN
         IF(IBINP1(1).GE.0) THEN
            IF(MAXGAM.NE.MAXGM) THEN
               WRITE(IW,*) 'CONFUSION IN CIMAT, MAXGAM',MAXGAM,MAXGM
               CALL ABRT
            ENDIF
         ELSE
            IF(MAXGAM.LT.MAXGM) THEN
               WRITE(IW,*) 'INSUFFICENT GUESS OF MAXGAM',MAXGAM,MAXGM
               CALL ABRT
            ELSE
               MAXGAM=MAXGM
            ENDIF
         ENDIF
      ELSE
         IF(MAXGM.LT.MAXGAM) THEN
            IF(MASWRK.OR.OUT1) WRITE(IW,9800) MAXGAM,MAXGM,ME
C           it can be continued because the buffer is big enough
            IF(ACTION.NE.USUAL) THEN
               WRITE(IW,9830)
            ELSE
Cnb            maxgam=maxgm
            ENDIF
C
C       this will be helpful for the second ms (if any)
C       since then it would be a perfect fit (hopefully?)
C       for the runs with action=save one would then have a maxgam which
C       will not be known to the next run with action=read if set
C       maxgam=maxgm. A good solution is use savdsk=.t. in this case.
C
         ELSE IF(MAXGM.NE.MAXGAM) THEN
               WRITE(IW,*) 'MAXGAM IS TOO SMALL ON NODE',ME
C           it must be stopped because the buffer is too small.
                 CALL ABRT
              ELSE
                 WRITE(IW,9810) MAXGAM,ME
         ENDIF
      ENDIF
C
      RETURN
 2349 FORMAT(9I8)
 4972 FORMAT(//,T20,'KFIG1=',I10,5X,'KFIG2=',I10)
 7666 FORMAT(/1X,'THE NUMBER OF PRODUCTIVE CSF PAIRS NMATE= ',I10)
 7776 FORMAT(3X,'NDT1:',I5,3X,'NDT2:',I5,3X,'IC1:',I5,3X,'IC2:',I5)
 8990 FORMAT(/1X,'ALL NODES TOGETHER GENERATED AND SORTED:',/)
 9000 FORMAT(I12,' 1-PARTICLE ACTIVE ORBITAL FFS,',/,
     *       I12,' 2-PARTICLE ACTIVE ORBITAL FFS,',/,
     *       1X,61(1H-)/
     *       I12,' TOTAL FOR MS=',I2,/)
 9800 FORMAT(/1X,'PREDICTED SORTING BUFFER LENGTH',I7,
     *         ' IS GREATER THAN NEEDED',I7,' ON NODE',I3,/)
 9810 FORMAT(/1X,'PREDICTED SORTING BUFFER LENGTH',I7,
     *       ' IS OPTIMUM ON NODE',I3,/)
 9830 FORMAT(/
     *  1X,'FYI: SECOND MS (IF ANY) CANNOT BENEFIT FROM THE KNOWN'/
     *  1X,'MAXGAM BECAUSE YOU ARE SAVING FORM FACTORS AND NOT ',/,
     *  1X,'USING DSKSAV=.T., CONTINUING...',/)
      END
C
C*MODULE SOFFAC  *DECK COMPAR
      SUBROUTINE COMPAR(IC1,IC2,IBUFII,IBUFJJ,NDT1,NDT2,MASKL,IZRTST)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION IBUFII(*),IBUFJJ(*)
C
C *** THIS ROUTINE COMPARES ORBITAL OCCUPATION SCHEMES BETWEEN
C *** SINGLET AND TRIPLET CSFS
C
      MASK12 = 2**12 - 1
C
C *** 1ST CONFIGURATION
C
      NDT1 = IAND(MASK12,IBUFII(IC1))
      ICFIG = IAND(MASKL,ISHFT(IBUFII(IC1),-12))
C
C *** 2ND CONFIGURATION
C
      NDT2 = IAND(MASK12,IBUFJJ(IC2))
      JCFIG = IAND(MASKL,ISHFT(IBUFJJ(IC2),-12))
C
C *** COMPARE ORBITAL OCCUPATION SCHEMES,
C *** FIND THE NUMBER OF DISCOINCIDENCES
C
      NCMPR = IEOR(ICFIG,JCFIG)
      NDISCO=KOUNTBIT(NCMPR,NAOD)
C
C *** IF THERE ARE MORE THEN 2 DISCOINCIDENCES (NDISCO > 5) return 1
C
      IZRTST = 1
      IF (NDISCO .LE. 5) IZRTST = 0
C
C *** DIVIDE BY TWO TO GET THE NUMBER OF SPATIAL DISCOINCIDENCES
C
C     NORBS = NDISCO/2
      RETURN
      END
C
C*MODULE SOFFAC  *DECK EMPTY
      SUBROUTINE EMPTY(ISODA,IBINP,INDXL,STOREL,NUMBR,NOFF,SORTED,IBINP1
     *                ,SAVDSK,SLOWFF,IFF,LDARISO)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MAXCP=4096)
      LOGICAL SAVDSK,SLOWFF
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SOGUG/  CP(MAXCP),NUNIQ,LIOBP,ISODAF,NSODA,LSTREC,JSODAF,
     *                JSODA,NRECJ
      COMMON /SOBUF/  IBNLEN,IBN2,MAXGAM,MAXFL,MAXFL2,ND1FZ,ND2FZ
C
      DIMENSION ISODA(*),IBINP(*),INDXL(*),STOREL(IBNLEN,*),SORTED(*),
     *          IBINP1(*)
C
C *** THIS ROUTINE MAKES SURE THAT ALL THE SORTING BINS ARE EMPTY.
C *** CALLED FROM PARLST
C
C *** VARIABLES
C
C       IBINP  - POINTER ARRAY WHICH KEEPS TRACK OF THE NUMBER OF FORM
C                FOR INDEX I WRITTEN TO DISK
C       INDXL  - ARRAY WHICH KEEPS TRACK OF THE NUMBER OF FORM FACTORS
C                CURRENTLY HELD IN ARRAY STOREL FOR A GIVEN PAIR (FOUR)
C                INDEX
C       STOREL - HOLDS THE FORM FACTORS ACCORDING TO PAIR OR FOUR INDEX
C
      DO 100 I=1,NUMBR
         NSOFAR = IBINP(NOFF+I)
         NLEFT = INDXL(I)
         IF (NLEFT.NE.0) THEN
           IF((.NOT.SAVDSK).OR.IBINP1(1).GE.0) THEN
             IF (NSOFAR+NLEFT .GT. MAXGAM) THEN
               WRITE(IW,*) '*** ERROR IN EMPTY: INCREASE RECLEN',
     *                     NSOFAR+NLEFT,MAXGAM
               CALL ABRT
             ENDIF
             IF(SLOWFF) THEN
               IF (NSOFAR.NE.0) CALL RAREAD(ISODAF,ISODA,SORTED,NSOFAR,
     *                                      NOFF+I+LSTREC,0)
               CALL DCOPY(NLEFT,STOREL(1,I),1,SORTED(NSOFAR+1),1)
               NPHYSOFF=0
               CALL RAWRITE(ISODAF,ISODA,SORTED,NSOFAR+NLEFT,
     *                      NSOFAR+NLEFT,NOFF+I+LSTREC,NPHYSOFF)
             ELSE
               NPHYSOFF=NSOFAR/LDARISO
               CALL RAWRITE(ISODAF,ISODA,STOREL(1,I),NSOFAR+NLEFT,
     *                      NLEFT,NOFF+I+LSTREC,NPHYSOFF)
             ENDIF
           ENDIF
           IBINP(NOFF+I) = NSOFAR + NLEFT
           INDXL(I) = 0
C         double check if the precalculated sizes are equal to saved
           IF(SAVDSK.AND.IBINP1(1).GE.0.AND.IBINP(NOFF+I).NE.
     *        IBINP1(NOFF+I)) THEN
               WRITE(IW,*) 'IBINP SIZE MISMATCH',IBINP(NOFF+I),
     *                     IBINP1(NOFF+I),I,NOFF
               CALL ABRT
            ENDIF
         ENDIF
         IFF=IFF+IBINP(NOFF+I)
C
 100  CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK FNDALL
      SUBROUTINE FNDALL(NORB,LOC,NN)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION LOC(*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
C
C *** THIS ROUTINE RETURNS THE LOCATION OF ALL ORBITALS IN NORB,
C *** I.E., THE ORBITAL NUMBERS
C
      MASK1 = 1
      NTEMP = NORB
      J = 0
      DO 100 I=NAO,1,-1
        IF (IAND(MASK1,NTEMP) .EQ. 1) THEN
          J = J + 1
          LOC(J) = I
        ENDIF
        NTEMP = ISHFT(NTEMP,-1)
100   CONTINUE
      IF (J .NE. NN) THEN
        WRITE(IW,'(A)')' ERROR - J NE NN IN FNDALL'
        CALL ABRT
      ENDIF
C
C *** SWITCH THE ORDER OF THE ORBITALS
C
      K=J
      DO 200 I=1,J/2
         NTEMP = LOC(I)
         LOC(I) = LOC(K)
         LOC(K) = NTEMP
         K = K - 1
200   CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK GETFFA
      SUBROUTINE GETFFA(COEF,C1,C2,C3,C4,GAMMA,SPNDT1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(ZERO=0.0D+00,HALF=0.5D+00,THREE2=1.5D+00)
      INTEGER GAMMA,SPNDT1
C
C *** THIS ROUTINE ASSIGNS VALUES TO THE COULOMB AND EXCHANGE FORM
C *** FACTORS WHEN THE TWO COINCIDENT ORBITALS HAVE ALPHA SPIN
C
        IF (GAMMA .EQ. 0) THEN
           IF (SPNDT1 .EQ. 0) THEN
C
C ***        NOTE THERE IS NO EXCHANGE TERM FOR GAMMA = 0 WITH BOTH
C ***        ORBITALS BETA (FOR ALPHA COINCIDENT ORBITALS ONLY)
C
              C1 = -HALF*COEF
              C2 = HALF*COEF
              C3 = ZERO
              C4 = ZERO
           ELSE
              C1 = THREE2*COEF
              C2 = C1
              C3 = -THREE2*COEF
              C4 = C3
           ENDIF
        ELSE IF (GAMMA .EQ. 1) THEN
           C1 = COEF
           C2 = HALF*COEF
           C3 = -C2
           C4 = -COEF
        ELSE
           C1 = COEF
           C2 = HALF*COEF
           C3 = -C1
           C4 = -C2
        ENDIF
      RETURN
      END
C
C*MODULE SOFFAC  *DECK GETFFB
      SUBROUTINE GETFFB(COEF,C1,C2,C3,C4,GAMMA,SPNDT1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(ZERO=0.0D+00,HALF=0.5D+00,THREE2=1.5D+00)
      INTEGER GAMMA,SPNDT1
C
C *** THIS ROUTINE ASSIGNS VALUES TO THE COULOMB AND EXCHANGE FORM
C *** FACTORS WHEN THE TWO COINCIDENT ORBITALS HAVE BETA SPIN
C
C
        IF (GAMMA .EQ. 0) THEN
C
           IF (SPNDT1 .EQ. 0) THEN
              C1 = -THREE2*COEF
              C2 = C1
              C3 = -C1
              C4 = C3
           ELSE
C
C ***         NOTE THERE IS NO EXCHANGE TERM FOR GAMMA = 0 WITH BOTH
C ***         ORBITALS ALPHA (FOR BETA COINCIDENT ORBITALS ONLY)
C
              C1 = HALF*COEF
              C2 = -C1
              C3 = ZERO
              C4 = ZERO
           ENDIF
        ELSE IF (GAMMA .EQ. 1) THEN
           C1 = COEF
           C2 = HALF*COEF
           C3 = -C1
           C4 = -C2
        ELSE
           C1 = COEF
           C2 = HALF*COEF
           C3 = -C2
           C4 = -C1
        ENDIF
      RETURN
      END
C
C*MODULE SOFFAC  *DECK KEPTRK
      SUBROUTINE KEPTRK(ISODA,INDXL,NFLAG,STOR,NUMBR,KF12,NOFF,SORTED,
     *             STOREL,ISTORL,IBINP,MMM,IBINP1,SAVDSK,SLOWFF,LDARISO)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MAXCP=4096,ZERO=0.0D+00)
      LOGICAL SAVDSK,SLOWFF
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /SOGUG/  CP(MAXCP),NUNIQ,LIOBP,ISODAF,NSODA,LSTREC,JSODAF,
     *                JSODA,NRECJ
      COMMON /SOBUF/  IBNLEN,IBN2,MAXGAM,MAXFL,MAXFL2,ND1FZ,ND2FZ
      COMMON /TRNMOM/ OPERR,SYMTOL,NOCC,NUMCI,NFZC,NUMVEC,ICI,MXRT,NSTAT
C
      DIMENSION ISODA(*),INDXL(*),STOR(*),NFLAG(*),STOREL(IBNLEN,*),
     *          ISTORL(IBN2,*),IBINP(*),SORTED(*),IBINP1(*)
C
C     DATA TOL/1.0D-08/,zero/0.0d+00/
C
C *** THIS ROUTINE STORES THE ONE AND TWO PARTICLE FORM FACTORS IN A
C *** PACKED MANNER, I.E., WITH THEIR CSF LABEL, ETC.
C
C    NFLAG  - OUTPUT FROM SLATER, ARRAY WHICH CONTAINS ORBITAL INDICES
C    STOR   - OUTPUT FORM SLATER, ARRAY WHICH CONTAINS FORM FACTOR FOR
C             ORBITAL INDEX IN NFLAG
C    INDXL  - KEEPS TRACK OF THE NUMBER OF FORM FACTORS STORED THUS
C             FAR FOR A GIVEN PAIR OR FOUR INDEX
C
C     small is determined by the Clebsch-Gordon coefs (eventually)
C     (through coefficients in front of determinants)
C
      SMALL=SYMTOL*1.0D-04
      DO 100 K=1,NUMBR
         INDEX = NFLAG(K)
         COEF = STOR(INDEX)
         STOR(INDEX) = ZERO
         IF (ABS(COEF).GE.SMALL) THEN
            IF(INDEX.GT.MMM) THEN
               WRITE(IW,*) 'TOO LARGE INDEX',INDEX,MMM
               CALL ABRT
            ENDIF
            NFF = INDXL(INDEX)
            NFF = NFF + 1
C           ISTORL((2*NFF-1),INDEX) = CSFIDX
            ISTORL((NFF-1)*NWDVAR+1,INDEX) = KF12
            NFF = NFF + 1
            STOREL(NFF,INDEX) = COEF
C
C *** CHECK TO SEE IF THE BIN IS FULL, IF IT IS THEN READ IN THE RECORD
C *** CORRESPONDING TO THAT INDEX, ADD IN THE ADDITIONAL FORM FACTORS
C *** AND WRITE IT BACK OUT
C
            IF (NFF .GE. IBNLEN-1) THEN
              NSOFAR = IBINP(NOFF+INDEX)
              IF((.NOT.SAVDSK).OR.IBINP1(1).GE.0) THEN
C
C             ie if doing real sorting, not just FF calculation
C
                IF ((NSOFAR+NFF) .GT. MAXGAM) THEN
                  WRITE(IW,*) '*** ERROR KEPTRK: INCREASE RECLEN',
     *                        NSOFAR+NFF,MAXGAM
                  CALL ABRT
                ENDIF
C           lenrec is the maximum
C           number of words for this record, although only first
C           NSOFAR+NFF are meaningful numbers for the first rawrit on
C           this record. The rest is a balast to reserve
C           space on the file for future rawrit's to that record
                IF(SAVDSK) THEN
                  LENREC=IBINP1(NOFF+INDEX)
                ELSE
                  LENREC=MAXGAM
                ENDIF
C
C          implement two different strategies
C          one (ibnlen.lt.ldar) reads the existing record,
C          appends the new chunk in RAM and writes off all to disk.
C          The other uses chunk size multiple to the number
C          of physical records. The new chunk is written directly to
C          disk with a shift in physical records.
C
                IF(SLOWFF) THEN
                  IF(NSOFAR.NE.0) CALL RAREAD(ISODAF,ISODA,SORTED,
     *                                    NSOFAR,NOFF+INDEX+LSTREC,0)
                  CALL DCOPY(NFF,STOREL(1,INDEX),1,SORTED(NSOFAR+1),1)
                  NPHYSOFF=0
                  CALL RAWRITE(ISODAF,ISODA,SORTED,LENREC,NSOFAR+NFF,
     *                         NOFF+INDEX+LSTREC,NPHYSOFF)
                ELSE
                  NPHYSOFF=NSOFAR/LDARISO
                  IF(NFF.NE.IBNLEN) THEN
                    WRITE(IW,*) 'QUEER NFF',NFF,IBNLEN
                    CALL ABRT
                  ENDIF
                  CALL RAWRITE(ISODAF,ISODA,STOREL(1,INDEX),LENREC,NFF,
     *                         NOFF+INDEX+LSTREC,NPHYSOFF)
                ENDIF
              ENDIF
C              IBINP(NOFF+INDEX) = IBINP(NOFF+INDEX) + ibnlen
              IBINP(NOFF+INDEX) = IBINP(NOFF+INDEX) + NFF
              INDXL(INDEX) = 0
            ELSE
              INDXL(INDEX) = NFF
            ENDIF
C 200    continue
         ENDIF
  100 CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ONED2P
      SUBROUTINE ONED2P(LOC1,LOC2,GAMMA,NCOIN,COEF,SPNDT1,D2Z,JFLAGZ,
     *                  LOCA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*),LOCA(*)
      INTEGER GAMMA,SPNDT1
C
C *** THIS ROUTINE FORMS THE TWO PARTICLE FORM FACTORS WHEN
C *** THERE IS ONE DISCOINCIDENCE. CALLED FROM ONEDIS
C
C *** BEGIN BY PICKING UP THE COINCIDENT ACTIVE ALPHA AND BETA ORBITALS
C
      NALPHA = ISHFT(NCOIN,-NAO)
      NBETA  = IAND(NCOIN,MSKNAO)
C
C *** COUNT THE NUMBER OF COINCIDENT ALPHA AND BETA SPINS
C
      KOUNTA=KOUNTBIT(NALPHA,NAO)
      KOUNTB=KOUNTBIT(NBETA,NAO)
C
C *** GET FORM FACTORS FOR ALPHA COINCIDENT ORBITALS
C
      IF (KOUNTA .GT. 0) CALL ALPHFF(NALPHA,KOUNTA,GAMMA,SPNDT1,COEF,
     *   LOC1,LOC2,D2Z,JFLAGZ,LOCA)
C
C *** GET FORM FACTORS FOR BETA COINCIDENT ORBITALS
C
      IF (KOUNTB .GT. 0) CALL BETAFF(NBETA,KOUNTB,GAMMA,SPNDT1,COEF,
     *   LOC1,LOC2,D2Z,JFLAGZ,LOCA)
C
C *** GET THE FORM FACTORS INVOLVING CORE AND ACTIVE ORBITALS
C
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ONEDIS
      SUBROUTINE ONEDIS(CPROD,NDET1,NDET2,ORDET,D1Z,D2Z,IFLAGZ,JFLAGZ,
     *                  LOCA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D1Z(*),D2Z(*),IFLAGZ(*),JFLAGZ(*),LOCA(*)
C
      INTEGER ORDET,GAMMA,SPNDT1,SPNDT2
C
C *** THIS ROUTINE GENERATES THE ONE AND TWO PARTICLE FORM FACTORS
C *** WHEN THERE IS ONE DISCOINCIDENT ORBITAL
C
C *** BEGIN BY PICKING UP THE DISCOINCIDENT ORBITAL IN DETERMINANT 1
C *** (NDET1 - SINGLET) AND IN NDET2 (TRIPLET)
C
      NDIS1 = IAND(NDET1,ORDET)
      NDIS2 = IAND(NDET2,ORDET)
C
C *** FIND THE COINCIDENT ORBITALS
C
      NCOIN = IAND(NDET1,NDET2)
C
C *** FIND THE SPIN OF THE DISCOINCIDENT ORBITALS
C *** SPNDT1 = 0 BETA SPIN, = 1 ALPHA SPIN
C
      SPNDT1 = 0
      IF (NDIS1 .GT. MSKNAO) SPNDT1 = 1
      SPNDT2 = 0
      IF (NDIS2 .GT. MSKNAO) SPNDT2 = 1
C
C *** DETERMINE WHICH SPIN OPERATOR WILL CONNECT THE DISCO ORBITALS.
C *** <A//A> OR <B//B> GAMMA = 0 (S(Z)); <A//B> GAMMA = 1 (S(+));
C *** <B//A> GAMMA = -1 (S(-))
C
      IF (SPNDT1.EQ.SPNDT2) THEN
         GAMMA = 0
      ELSE IF (SPNDT1 .EQ. 0) THEN
         GAMMA = -1
      ELSE
         GAMMA = 1
      ENDIF
C
C *** FIND PHASE CHANGE THAT OCCURS WHEN PUTTING THE DETERMINANTS
C *** INTO MAXIMUM COINCIDENCE
C
      CALL PHASE2(NDIS1,NDIS2,NCOIN,KOUNT,NAOD)
      COEF = CPROD
      IF (MOD(KOUNT,2) .EQ. 1) COEF = -COEF
C
C *** FIND ORBITAL INDICES
C
      NORB1 = NDIS1
      IF (SPNDT1 .EQ. 1) NORB1 = ISHFT(NORB1,-NAO)
      NORB2 = NDIS2
      IF (SPNDT2 .EQ. 1) NORB2 = ISHFT(NORB2,-NAO)
      LOC1=IFINDEM(NORB1)
      LOC2=IFINDEM(NORB2)
C
C *** STORE THE ONE PARTICLE FORM FACTOR
C
      CALL SONEP(LOC1,LOC2,COEF,SPNDT1,SPNDT2,D1Z,IFLAGZ)
C
C *** NOW DO THE TWO ELECTRON INTEGRALS
C *** ONED2P - ONE DISCO 2 PARTICLE FORM FACTORS
C
      CALL ONED2P(LOC1,LOC2,GAMMA,NCOIN,COEF,SPNDT1,D2Z,JFLAGZ,LOCA)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK OUTFB
      SUBROUTINE OUTFB(ISODA,IBINP,GAMMA)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION ISODA(*),IBINP(*),GAMMA(*)
C
C *** THIS ROUTINE PRINTS OUT THE FORM FACTORS SORTED BY INDEX
C
C *** VARIABLES
C
C     IBINP - POINTER ARRAY WHICH KEEPS TRACK OF THE NUMBER OF FORM
C             FOR INDEX I WRITTEN TO DISK
C
C *** DO ONE PARTICLE FORM FACTORS
C
      NOFF = 0
      IF(GOPARR) THEN
         WRITE(IW,9010) ME,'1P'
      ELSE
         WRITE(IW,9000) '1P'
      ENDIF
      CALL PSORT(ISODA,IBINP,NAO2,NOFF,GAMMA,GAMMA)
C
C *** DO TWO PARTICLE FORM FACTORS
C
      NOFF = NOFF + NAO2
      IF(GOPARR) THEN
         WRITE(IW,9010) ME,'2P'
      ELSE
         WRITE(IW,9000) '2P'
      ENDIF
      CALL PSORT(ISODA,IBINP,NAO4,NOFF,GAMMA,GAMMA)
      RETURN
 9000 FORMAT(//T10,'*** ',A3,
     *       ' FORM FACTORS SORTED BY ORBITAL INDEX ***',//)
 9010 FORMAT(//I3,':',T10,'*** ',A3,
     *       ' FORM FACTORS SORTED BY ORBITAL INDEX ***',//)
      END
C
C*MODULE SOFFAC  *DECK PHASE2
      SUBROUTINE PHASE2(NORB1,NORB2,NCOIN,KOUNTT,NUM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C *** THIS ROUTINE DETERMINES THE NUMBER OF COINCIDENT ORBITALS (KOUNT)
C *** BETWEEN THE TWO DISCOINCIDENT ORBITALS NORB1 AND NORB2.
C *** NCOIN - CONTAINS THE COINCIDENT ORBITALS. IF KOUNT IS EVEN NO
C     PHASE CHANGE OCCURS WHEN PUTTING THE ORBITALS INTO MAXIMUM
C     COINCIDENCE.
C
      IF (NORB1 .GT. NORB2) THEN
        N = NORB1 - 2*NORB2
      ELSE
        N = NORB2 - 2*NORB1
      ENDIF
      NBITS = IAND(NCOIN,N)
C
C *** COUNT THE 1 BITS (COINCIDENT ORBITALS) IN NBITS
C
      KOUNTT=KOUNTBIT(NBITS,NUM)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK PSORT
      SUBROUTINE PSORT(ISODA,IBINP,NUMBR,NOFF,SORTED,ISORTD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MAXCP=4096)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /SOGUG/  CP(MAXCP),NUNIQ,LIOBP,ISODAF,NSODA,LSTREC,JSODAF,
     *                JSODA,NRECJ
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DIMENSION ISODA(*),IBINP(*),SORTED(*),ISORTD(*)
C      EQUIVALENCE (SORTED,ISORTD)
C
C *** THIS ROUTINE PRINTS SORTED FORM FACTORS, IT IS CALLED FROM OUTFB
C
      DO 100 I=1,NUMBR
         IF(GOPARR) THEN
            WRITE(IW,9000) ME,I
         ELSE
            WRITE(IW,9005) I
         ENDIF
         NSOFAR = IBINP(NOFF+I)
         IF (NSOFAR .NE. 0) THEN
            CALL RAREAD(ISODAF,ISODA,SORTED,NSOFAR,NOFF+I+LSTREC,0)
            DO 50 J=1,NSOFAR,2
C        2 means a pair of (real,one 64 bit int or two 32 bit ints)
               IF(GOPARR) THEN
                  WRITE(IW,9010) ME,ISORTD((J-1)*NWDVAR+1),SORTED(J+1)
               ELSE
                  WRITE(IW,9020)    ISORTD((J-1)*NWDVAR+1),SORTED(J+1)
               ENDIF
  50        CONTINUE
         ENDIF
 100  CONTINUE
      RETURN
 9000 FORMAT(/,I3,':',T20,'*** ORBITAL INDEX = ',I5,' ***')
 9005 FORMAT(/,T20,'*** ORBITAL INDEX = ',I5,' ***')
 9010 FORMAT(I3,':','LABEL = ',O22,5X,'FF = ',E22.15)
 9020 FORMAT(2X,'LABEL = ',O22,5X,'FF = ',E22.15)
      END
C
C*MODULE SOFFAC  *DECK SLATER
      SUBROUTINE SLATER(NDT1,NDT2,IC1,IC2,IBUFII,IBUFJJ,CP1,CP3,D1Z,D2Z,
     *                  IFLAGZ,JFLAGZ,LOCA,LOCB,MASKL,OUT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION IBUFII(*),IBUFJJ(*),CP1(*),CP3(*),D1Z(*),D2Z(*),
     *          IFLAGZ(*),JFLAGZ(*),LOCA(*),LOCB(*)
      LOGICAL OUT
C
C *** THIS ROUTINE GENERATES THE ONE AND TWO PARTCILE FORM FACTORS
C *** FOR A PAIR OF SLATER DETERMINANTS
C
      MASK12 = 2**12 - 1
      NOFF1 = IC1 + 1
      NOFF2 = IC2 + 1
C
C *** LOOP OVER TRIPLET DETERMINANTS
C
      DO 8000 JJ=1,NDT2
        IDET2 = IBUFJJ(NOFF2+JJ)
        NDET2 = IAND(MASKL,ISHFT(IDET2,-12))
        NORM2 = IAND(MASK12,IDET2)
C
C ***   LOOP OVER SINGLET DETERMINENTS
C
        DO 7000 II=1,NDT1
          IDET1 = IBUFII(NOFF1+II)
C
C ***     FIND THE DISCOINCIDENT ORBITALS
C
C ***     STRIP NORMALIZATION OFF DETERMINANTS 1 AND 2
C
          NDET1 = IAND(MASKL,ISHFT(IDET1,-12))
          LORDET = IEOR(NDET1,NDET2)
C
C ***     COUNT THE DISCOINCIDENT BITS = 2*(NUMBER OF DISCO'S IN ORDET)
C
          NDISCO=KOUNTBIT(LORDET,NAOD)
C
C ***  GOTO TO WORK IF THERE ARE LESS THAN 3 DISCOINCIDENCES (NDISCO<5)
C
          IF (NDISCO .GE. 5) GOTO 7000
C
C ***     BEGIN BY FORMING THE COEFFICIENT PRODUCT
C
          NORM1 = IAND(MASK12,IDET1)
          CPROD = CP1(NORM1) * CP3(NORM2)
C
C ***     DETERMINE WHICH SUBROUTINE TO BRANCH TO
C
          IF(OUT) WRITE(IW,9000)NDISCO,CPROD,NDET1,NDET2
          IF (NDISCO .EQ. 0) THEN
             CALL ZDISCO(CPROD,NDET1,D1Z,D2Z,IFLAGZ,JFLAGZ,LOCA,LOCB)
          ELSE IF (NDISCO .LE. 2) THEN
             CALL ONEDIS(CPROD,NDET1,NDET2,LORDET,D1Z,D2Z,IFLAGZ,JFLAGZ,
     *                   LOCA)
          ELSE
             CALL TWODIS(CPROD,NDET1,NDET2,LORDET,D2Z,JFLAGZ)
          ENDIF
C
7000    CONTINUE
8000  CONTINUE
      RETURN
9000  FORMAT(//,' NDISCO,CPROD,NDET1,NDET2 FROM SLATER =',
     *       /,I5,2X,E22.14,2X,O11.11,2X,O11.11)
      END
C
C*MODULE SOFFAC  *DECK SOFFAC
      SUBROUTINE SOFFAC(ISODA,IBUFII,IBUFJJ,GAMMA,FBIN1Z,FBIN2Z,IBINP,
     *                  IBP1Z,IBP2Z,IFLAGZ,JFLAGZ,LOCA,LOCB,D1Z,D2Z,
     *                  IBINP1,SAMEMUL,SAVDSK,SLOWFF,ACTION,ICI1,ICI2,
     *                  LDARISO,IBTYP1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MAXCP=4096)
      LOGICAL SAMEMUL,SAVDSK,SLOWFF,GOPARR,DSKWRK,MASWRK
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SOGUG/  CP(MAXCP),NUNIQ,LIOBP,ISODAF,NSODA,LSTREC,JSODAF,
     *                JSODA,NRECJ
      COMMON /SOBUF/  IBNLEN,IBN2,MAXGAM,MAXFL,MAXFL2,ND1FZ,ND2FZ
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
C
      DIMENSION ISODA(*),IBUFII(*),IBUFJJ(*),GAMMA(*),FBIN1Z(*),
     *          FBIN2Z(*),IBINP(*),IBP1Z(*),IBP2Z(*),IFLAGZ(*),JFLAGZ(*)
     *         ,LOCA(*),LOCB(*),D1Z(*),D2Z(*),IBINP1(*)
C
      DATA USUAL/8HNORMAL  /
C     -----------------------------------------------------------------
C     THIS PROGRAM GENERATES THE ONE AND TWO PARTICLE FORM
C     FACTORS FOR THE CALCULATION OF SPIN ORBIT COUPLING
C     CONSTANTS BETWEEN SINGLET AND TRIPLET ELECTRONIC STATES
C                  TRF 9/21/84
C     -----------------------------------------------------------------
C
C     if(maswrk) WRITE(IW,1001)
C
      IF(SAVDSK) IBINP1(1)=-1
C     this is a hidden argument to cimat, if ibinp(1)<0 then the call
C     of cimat is preliminary (no data written to disk etc)
C     in the second call to cimat ibinp1 inside of cimat will be ibinp
C     from the first call to cimat, i.e. ibinp1(1)>=0
C
      IF(.NOT.SAVDSK.AND.ACTION.EQ.USUAL)
     *   CALL RAOPEN2(ISODAF,ISODA,1,NSODA,MAXGAM,LDARISO,NPRINT)
C
      CALL CIMAT(ISODA,IBUFII,IBUFJJ,GAMMA,D1Z,D2Z,IBP1Z,IBP2Z,IFLAGZ,
     *           JFLAGZ,LOCA,LOCB,FBIN1Z,FBIN2Z,IBINP,IBINP1,
     *           SAMEMUL,SAVDSK,SLOWFF,ACTION,ICI1,ICI2,LDARISO,IBTYP1)
      CALL TIMIT(1)
C
      IF(SAVDSK) THEN
         IF(MASWRK) WRITE(IW,9010)
C     in the first call to CIMAT, save the record lenths in ibinp1
C     then provide this to the second call to cimat
C     note that the order of ibinp and ibinp1 is switched now
C     and ibinp(1) is >= 0
         IF(ACTION.EQ.USUAL)
     *      CALL RAOPEN2(ISODAF,ISODA,1,NSODA,MAXGAM,LDARISO,NPRINT)
         CALL CIMAT(ISODA,IBUFII,IBUFJJ,GAMMA,D1Z,D2Z,IBP1Z,IBP2Z,
     *              IFLAGZ,JFLAGZ,LOCA,LOCB,FBIN1Z,FBIN2Z,IBINP1,IBINP,
     *              SAMEMUL,SAVDSK,SLOWFF,ACTION,ICI1,ICI2,LDARISO,
     *              IBTYP1)
         CALL TIMIT(1)
      ENDIF
C
      RETURN
C1001 FORMAT(/1X,'Calculating the form factors... ',/)
 9010 FORMAT(/1X,'DOING THE SECOND FF CALCULATION...',/)
      END
C
C*MODULE SOFFAC  *DECK SONEP
      SUBROUTINE SONEP(LOC1,LOC2,COEF,SPNDT1,SPNDT2,D1Z,IFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      COMMON /SOBUF/  IBNLEN,IBN2,MAXGAM,MAXFL,MAXFL2,ND1FZ,ND2FZ
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D1Z(*),IFLAGZ(*)
      INTEGER SPNDT1,SPNDT2
C
C *** THIS ROUTINE STORE THE ONE PARTICLE FORM FACTOR FOR ONE DISCO
C
C *** FORM THE PAIR INDEX
C
      IJ = NAO*(LOC1-1) + LOC2
      C = COEF*HALF
      IF ( (SPNDT1.EQ.0) .AND. (SPNDT2.EQ.0) ) C = -C
      CALL STORE2(IJ,D1Z,ND1FZ,IFLAGZ,C)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK STORE2
      SUBROUTINE STORE2(INDX,STOR,NUMBR,KFLAG,COEF)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION STOR(*),KFLAG(*)
C
C *** THIS ROUTINE STORES THE TWO PARTICLE FORMFACTOR IN ARRAY STOR
C
C
C *** VARIABLES
C
C     NUMBR - KEEPS TRACK OF THE NUMBER OF UNIQUE FORM FACTORS
C             STORED SO FAR
C     KfLAG  - LIST OF UNIQUE FOUR INDEXES STORED SO FAR
C
      IF (STOR(INDX) .NE. 0.0D+00) THEN
        STOR(INDX) = STOR(INDX) + COEF
      ELSE
        NUMBR = NUMBR + 1
        KFLAG(NUMBR) = INDX
C       STOR(INDX) = STOR(INDX) + COEF
        STOR(INDX) = COEF
      ENDIF
      RETURN
      END
C
C*MODULE SOFFAC  *DECK TWODIS
      SUBROUTINE TWODIS(CPROD,NDET1,NDET2,ORDET,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D2Z(*),JFLAGZ(*)
      INTEGER ORDET,SPINS1,SPINS2
      LOGICAL ALLB
C
C *** THIS ROUTINE GENERATES TWO PARTICLE FORMFACTORS
C
C *** BEGIN BY DETERMINING THE TYPE OF DISCOINCIDENCE, SPIN, SPATIAL OR
C *** BOTH. FIRST PICK UP THE TWO DISCO ORBITALS IN DET1 (SINGLET) AND
C *** DET2 (TRIPLET)
C
      NDIS1 = IAND(NDET1,ORDET)
      NDIS2 = IAND(NDET2,ORDET)
C
C *** FIND THE COINCIDENT ORBITALS
C
      NCOIN = IAND(NDET1,NDET2)
C
C *** RIGHT SHIFT NDIS1 AND NDIS2 BY NAO POSITIONS, IF THE NUMBER OF
C *** BITS LEFT IS 0 THEN DISCO SPINS ARE BETA; 1= ONE OF EACH;
C     2= 2 ALPHA
C
      NTEMP1 = ISHFT(NDIS1,-NAO)
      NTEMP2 = ISHFT(NDIS2,-NAO)
      SPINS1=KOUNTBIT(NTEMP1,NAO)
      SPINS2=KOUNTBIT(NTEMP2,NAO)
C
C *** DETERMINE WHICH SUBROUTINE TO BRANCH TO, IF THERE ARE TWO SPIN
C *** DISCOS THEN GO GET A NEW PAIR OF DETERMINANTS
C
      IF (ABS(SPINS2-SPINS1) .EQ. 2 ) RETURN
C
      IF (SPINS1 .EQ. 0) THEN
          IF (SPINS2 .EQ. 0) THEN
             ALLB = .TRUE.
             CALL ALLSME(ALLB,NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
          ELSE
             CALL BBAB(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
          ENDIF
      ELSE IF (SPINS1 .EQ. 1) THEN
          IF (SPINS2 .EQ.0) THEN
             CALL ABBB(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
          ELSE IF (SPINS2 .EQ. 1) THEN
             CALL A2B2(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
          ELSE
             CALL ABAA(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
          ENDIF
      ELSE
          IF (SPINS2 .EQ. 1) THEN
             CALL AAAB(NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
          ELSE
             ALLB = .FALSE.
             CALL ALLSME(ALLB,NDIS1,NDIS2,NCOIN,CPROD,D2Z,JFLAGZ)
          ENDIF
      ENDIF
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ZDAAAA
      SUBROUTINE ZDAAAA(COEF,KOUNTA,LOCA,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(THREE2=1.5D+00)
      DIMENSION LOCA(*),D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE STORES TWO PARTICLE FORM FACTORS FOR INTEGRALS OF THE
C *** FOLLOWING TYPE <AA//AA> (ZERO DISCOS).  CALLED FROM ZROD2P
C
      C1 = THREE2*COEF
      C2 = -C1
C
C *** LOOP OVER THE ALPHA SPINS
C
      DO 200 M=1,KOUNTA
         I = LOCA(M)
         DO 100 N=M+1,KOUNTA
            J = LOCA(N)
C
C ***   FORM THE FOUR INDEX FOR BOTH COULOMB INTEGRALS AND BOTH EXCHANGE
C
C ***   STORE THE FORM FACTORS FOR THE COULOMB INTEGRALS IJIJ,JIJI
C
            CALL STOR4DX(I,J,I,J,C1,C1,D2Z,JFLAGZ)
C
C ***       STORE THE FORM FACTORS FOR THE EXCHANGE INTEGRALS IJJI,JIIJ
C
            CALL STOR4DX(I,J,J,I,C2,C2,D2Z,JFLAGZ)
  100   CONTINUE
  200 CONTINUE
C
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ZDABAB
      SUBROUTINE ZDABAB(COEF,KOUNTA,KOUNTB,LOCA,LOCB,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      DIMENSION LOCA(*),LOCB(*),D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE STORES TWO PARTICLE FORM FACTORS FOR INTEGRALS OF THE
C *** FOLLOWING TYPE <AB//AB> (ZERO DISCOS).  CALLED FROM ZROD2P
C
      C1 = -HALF*COEF
      C2 = -C1
C
C *** LOOP OVER THE ALPHA AND BETA SPINS
C
      DO 200 IA=1,KOUNTA
         I = LOCA(IA)
         DO 100 JB=1,KOUNTB
            J = LOCB(JB)
C
C ***       FORM THE FOUR INDEX FOR BOTH COULOMB INTEGRALS (THERE ARE NO
C ***       EXCHANGE INTEGRALS - TWO SPIN DISCO <AB//BA>
C
C ***       STORE THE FORM FACTORS FOR THE COULOMB INTEGRALS IJIJ,JIJI
C
            CALL STOR4DX(I,J,I,J,C1,C2,D2Z,JFLAGZ)
  100   CONTINUE
  200 CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ZDBBBB
      SUBROUTINE ZDBBBB(COEF,KOUNTB,LOCB,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(THREE2=1.5D+00)
      DIMENSION LOCB(*),D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE STORES TWO PARTICLE FORM FACTORS FOR INTEGRALS OF THE
C *** FOLLOWING TYPE <BB//BB> (ZERO DISCOS).  CALLED FROM ZROD2P
C
      C1 = -THREE2*COEF
      C2 = -C1
C
C *** LOOP OVER THE BETA SPINS
C
      DO 200 M=1,KOUNTB
         I = LOCB(M)
         DO 100 N=M+1,KOUNTB
            J = LOCB(N)
C
C ***   FORM THE FOUR INDEX FOR BOTH COULOMB INTEGRALS AND BOTH EXCHANGE
C
C ***   STORE THE FORM FACTORS FOR THE COULOMB INTEGRALS IJIJ,JIJI
C
            CALL STOR4DX(I,J,I,J,C1,C1,D2Z,JFLAGZ)
C
C ***       STORE THE FORM FACTORS FOR THE EXCHANGE INTEGRALS IJJI,JIIJ
C
            CALL STOR4DX(I,J,J,I,C2,C2,D2Z,JFLAGZ)
  100   CONTINUE
  200 CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ZDISCO
      SUBROUTINE ZDISCO(CPROD,NDET1,D1Z,D2Z,IFLAGZ,JFLAGZ,LOCA,LOCB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION D1Z(*),D2Z(*),IFLAGZ(*),JFLAGZ(*),LOCA(*),LOCB(*)
C
C *** THIS ROUTINE GENERATES THE ONE AND TWO PARTICLE FORM FACTORS
C *** WHEN THERE ARE NO DISCOINCIDENT ORBITALS
C
C *** GET THE ALPHA AND BETA ORBITALS
C
      NBETA = IAND(NDET1,MSKNAO)
      NALPHA = IAND(ISHFT(NDET1,-NAO),MSKNAO)
C
C *** COUNT THE NUMBER OF ALPHA AND BETA COINCIDENT ORBITALS
C
      KOUNTA=KOUNTBIT(NALPHA,NAO)
      KOUNTB=KOUNTBIT(NBETA,NAO)
C
C *** FIND THE ORBITAL INDICES OF THE COINCIDENT ORBITALS
C
      CALL FNDALL(NALPHA,LOCA,KOUNTA)
      CALL FNDALL(NBETA,LOCB,KOUNTB)
C
C *** STORE THE ONE PARTICLE FORM FACTORS
C
      CALL ZROD1P(CPROD,KOUNTA,KOUNTB,LOCA,LOCB,D1Z,IFLAGZ)
C
C *** STORE THE TWO PARTICLE FORM FACTORS
C
      CALL ZROD2P(CPROD,KOUNTA,KOUNTB,LOCA,LOCB,D2Z,JFLAGZ)
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ZROD1P
      SUBROUTINE ZROD1P(CPROD,KOUNTA,KOUNTB,LOCA,LOCB,D1Z,IFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(HALF=0.5D+00)
      COMMON /SOBUF/  IBNLEN,IBN2,MAXGAM,MAXFL,MAXFL2,ND1FZ,ND2FZ
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION LOCA(*),LOCB(*),D1Z(*),IFLAGZ(*)
C
C *** THIS ROUTINE STORES ONE PARTICLE FORM FACTORS FOR DETERMINANTS
C *** WHICH HAVE ZERO DISCOINCIDENT ORBITALS. CALLED FROM ZDISCO
C
      COEF = CPROD
      C1 = HALF*COEF
C
C *** LOOP OVER ALPHA SPIN ORBITALS FIRST
C
      DO 100 I=1,KOUNTA
C
C ***   PAIR INDEX
C
        IJ = NAO*(LOCA(I)-1) + LOCA(I)
C
C ***   KEEP A LIST OF INDICES (IN HCOREL) STORED THUS FAR
C
        CALL STORE2(IJ,D1Z,ND1FZ,IFLAGZ,C1)
  100 CONTINUE
C
C *** NOW DO THE BETA ORBITALS
C
      C2 = -C1
      DO 200 I=1,KOUNTB
        IJ = NAO*(LOCB(I)-1) + LOCB(I)
        CALL STORE2(IJ,D1Z,ND1FZ,IFLAGZ,C2)
  200 CONTINUE
      RETURN
      END
C
C*MODULE SOFFAC  *DECK ZROD2P
      SUBROUTINE ZROD2P(CPROD,KOUNTA,KOUNTB,LOCA,LOCB,D2Z,JFLAGZ)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION LOCA(*),LOCB(*),D2Z(*),JFLAGZ(*)
C
C *** THIS ROUTINE CALCULATES TWO PARTICLE FORM FACTORS FOR DETERMINANTS
C *** WHICH HAVE ZERO DISCOINCIDENT ORBITALS. CALLED FROM ZDISCO
C
C *** THERE ARE 4 CLASSES OF TWO E- INTEGRALS.  THE FIRST THREE INVOLVE
C *** ONLY ACTIVE ORBITALS THEY ARE <A(S),A(S)//A(T),A(T)>, <BB//BB>,
C *** AND <AB//AB> (A = ALPHA, B = BETA). THE LAST CLASS ARE ONES IN
C *** WHICH THERE ARE TWO CORE ORBITALS AND TWO ACTIVE ORBITALS
C *** NOTE: // = L12 WHERE
C ***
C *** L12 IS DEFINED AS THE -GAMMA COMPONENT OF (R12 X P1)/(R12**3)
C *** WHERE R12 IS THE DISTANCE AND DIRECTION BETWEEN ELECTRONS 1 & 2
C *** AND P1 IS THE MOMENTUM OPERATOR FOR ELECTRON 1.
C *** NOTE: (R12**3) REFERS TO MAGNITUDE ONLY
C
C *** FIRST CALCULATE FORM FACTORS FOR INTEGRALS OF THE ALL ACTIVE FORM
C ***                     <I(A),J(A)//I(A),J(A)>
C
      COEF = CPROD
      CALL ZDAAAA(COEF,KOUNTA,LOCA,D2Z,JFLAGZ)
C
C *** NOW DO <BB//BB>
C
      CALL ZDBBBB(COEF,KOUNTB,LOCB,D2Z,JFLAGZ)
C
C *** DO <AB//AB>
C
      CALL ZDABAB(COEF,KOUNTA,KOUNTB,LOCA,LOCB,D2Z,JFLAGZ)
C
C *** FINISH UP WITH THE CORE ORBITALS <CORE,ALPHA//CORE,ALPHA>
C
      RETURN
      END
C*MODULE SOFFAC  *DECK READCSF
      SUBROUTINE READCSF(ICI1,ICI2,IOBP,IBUFII,IBUFJJ,JSOD,SAMEMUL,OUT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MAXCP=4096)
      LOGICAL SAMEMUL,OUT
      DIMENSION IBUFII(*),IBUFJJ(*),IOBP(2,*),JSOD(*)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /SOGUG/  CP(MAXCP),NUNIQ,LIOBP,ISODAF,NSODA,LSTREC,JSODAF,
     *                JSODA,NRECJ
C
      IF(SAMEMUL) THEN
         IF(MS.NE.1) THEN
            MS1=0
         ELSE
            MS1=1
         ENDIF
         IF(MS.NE.-1) THEN
            MS2=0
         ELSE
            MS2=1
         ENDIF
      ELSE
         MS1=0
         MS2=1-MS
      ENDIF
C
      CALL RAREAD(JSODAF,JSOD,IBUFII,(IOBP(MS1+1,ICI1)-1)/NWDVAR+1,
     *            (ICI1-1)*NRECJ+MS1+1,1)
      CALL RAREAD(JSODAF,JSOD,IBUFJJ,(IOBP(MS2+1,ICI2)-1)/NWDVAR+1,
     *            (ICI2-1)*NRECJ+MS2+1,1)
C
      IF(OUT) THEN
         WRITE(IW,9783) IOBP(MS1+1,ICI1)
         WRITE(IW,9776)(IBUFII(KMP),KMP=1,IOBP(MS1+1,ICI1))
         WRITE(IW,9784) IOBP(MS2+1,ICI2)
         WRITE(IW,9776)(IBUFJJ(KMP),KMP=1,IOBP(MS2+1,ICI2))
      ENDIF
C
      RETURN
 9776 FORMAT(//,1X,5(O21))
 9783 FORMAT(////,T20,'LIST OF IBUFII (SINGLET CSFS) ',I10)
 9784 FORMAT(//,T20,'LIST OF IBUFJJ (TRIPLET CSFS) ',I10)
      END
C
C*MODULE SOFFAC  *DECK STOR4DX
      SUBROUTINE STOR4DX(I,J,K,L,CIJKL,CJILK,STOR,KFLAG)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION STOR(*),KFLAG(*)
      COMMON /SOBUF/  IBNLEN,IBN2,MAXGAM,MAXFL,MAXFL2,ND1FZ,ND2FZ
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
C
C     stor4dx is fourdx+store2 (each twice) in the old code
C     it calculates the four index and stores the form factor (2P)
C     either coulomb or exchange
C
C *** THIS ROUTINE CALCULATES A FOUR INDEX
C
C     call FOURDX(I,J,K,L,IJKL)
C     call fourdx(j,i,l,k,jilk)
      IJKL = NAO2*((NAO*(I-1) + J) -1) + NAO*(K-1) + L
      JILK = NAO2*((NAO*(J-1) + I) -1) + NAO*(L-1) + K
C
C     call store2(ijkl,stor,ND2FZ,kflag,cijkl)
      IF (STOR(IJKL) .NE. 0.0D+00) THEN
        STOR(IJKL) = STOR(IJKL) + CIJKL
      ELSE
        ND2FZ = ND2FZ + 1
        KFLAG(ND2FZ) = IJKL
        STOR(IJKL) = CIJKL
      ENDIF
C     call store2(jilk,stor,ND2FZ,kflag,cjilk)
      IF (STOR(JILK) .NE. 0.0D+00) THEN
        STOR(JILK) = STOR(JILK) + CJILK
      ELSE
        ND2FZ = ND2FZ + 1
        KFLAG(ND2FZ) = JILK
        STOR(JILK) = CJILK
      ENDIF
      RETURN
      END
C*MODULE SOFFAC  *DECK KOUNTBIT
      FUNCTION KOUNTBIT(NWORD,NUM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C *** THIS ROUTINE COUNTS THE NUMBER 1 BITS IN NWORD
C
      DIMENSION KKK(0:63)
      DATA KKK/0,1,1,2, 1,2,2,3, 1,2,2,3, 2,3,3,4,
     *         1,2,2,3, 2,3,3,4, 2,3,3,4, 3,4,4,5,
     *         1,2,2,3, 2,3,3,4, 2,3,3,4, 3,4,4,5,
     *         2,3,3,4, 3,4,4,5, 3,4,4,5, 4,5,5,6/
      KOUNTBIT = 0
      NTMP = NWORD
      NUM1 = (NUM-1)/12+1
      DO 100 IK=1,NUM1
       KOUNTBIT=KOUNTBIT+KKK(IAND(NTMP,63))+KKK(IAND(ISHFT(NTMP,-6),63))
       NTMP = ISHFT(NTMP,-12)
100   CONTINUE
      RETURN
      END
C*MODULE SOFFAC  *DECK IFINDEM
      FUNCTION IFINDEM(NORB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION KKK(0:63)
      DATA KKK/-1,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,
     *          4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,
     *          5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,
     *          4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0/
C
C     find the first "set" bit from the low end counting from the high
C     end, i.e. for 00100010 the answer is 7
C
C *** THIS ROUTINE RETURNS THE LOCATION OF THE DISCOINCIDENT ORBITAL
C *** IN NORB, I.E., ITS ORBITAL NUMBER
C
      NTEMP = NORB
      IFINDEM = NAO
  100 CONTINUE
         LOCC=KKK(IAND(NTEMP,63))
         IF(LOCC.LT.0) THEN
            IF(NTEMP.EQ.0) THEN
C              do not accept "not found"
               WRITE(IW,*) ' ERROR - COULD NOT FIND ORBITAL IN FINDEM'
               CALL ABRT
            ELSE
               NTEMP = ISHFT(NTEMP,-6)
               IFINDEM=IFINDEM-6
            ENDIF
         ELSE
            IFINDEM=IFINDEM-LOCC
            RETURN
         ENDIF
      GOTO 100
C     return
      END
C*MODULE SOFFAC  *DECK FIND2
      SUBROUTINE FIND2(NORB,LOC1,LOC2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /SPINFO/ NCORE,MS,NAO,NAOD,NAO2,NAO4,MSKNAO
      DIMENSION KKK1(0:15),KKK2(0:15)
      DATA KKK1/-1,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0/
      DATA KKK2/-1,-1,-1,1,-1,2,2,1,-1,3,3,1,3,2,2,1/
C
C     find the first and the second "set" bit from the low end counting
C     from the high end, i.e. for 00100010 the answer is 7 and 3
C
C *** THIS ROUTINE RETURNS THE LOCATION OF THE DISCOINCIDENT ORBITAL
C *** IN NORB, I.E., ITS ORBITAL NUMBER
C
      NTEMP = NORB
      LOC1 = NAO
      LOC2 = NAO
      IDONE=0
  100 CONTINUE
         I=IAND(NTEMP,15)
         LOCC1=KKK1(I)
         LOCC2=KKK2(I)
         IF(NTEMP.EQ.0) THEN
C           do not accept "not found"
            WRITE(IW,*) ' ERROR - COULD NOT FIND ORBITAL IN FIND2'
            CALL ABRT
         ELSE
            NTEMP = ISHFT(NTEMP,-4)
         ENDIF
         IF(LOCC1.GE.0) THEN
            IF(IDONE.NE.0) THEN
               LOC2=LOC2-LOCC1
               RETURN
            ELSE
               LOC1=LOC1-LOCC1
               IF(LOCC2.GE.0) THEN
                  LOC2=LOC2-LOCC2
                  RETURN
               ENDIF
               IDONE=1
            ENDIF
         ENDIF
         IF(IDONE.EQ.0) LOC1=LOC1-4
         LOC2=LOC2-4
      GOTO 100
C     return
      END
