      SUBROUTINE Set_version_hfodd_twocen
      !
      USE hfodd_sizes
      !
      CHARACTER NAMMOD*16
      INTEGER   MODUVE,MODSET
      !
      COMMON /VERMOD/ NAMMOD(NDMODU),MODUVE(NDMODU),MODSET(NDMODU)
      !
      NAMMOD(18) = 'hfodd_twocen    '
      MODUVE(18) = 14
      MODSET(18) = 18
      !
      END SUBROUTINE Set_version_hfodd_twocen
      !
      module hfodd_twocen
C=======================================================================
C
      USE hfodd_sizes
C
C=======================================================================
C
      CONTAINS
C
C=======================================================================
      SUBROUTINE INDGEN(ITWCEN)
C=======================================================================
C
      USE hfodd_sizes
C=======================================================================
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *      /TWCEDI/  LDTWCE
      COMMON
     *      /HERMIN/  INTWHE(1:NDTWCE,1:NDTWCE,1:NDTWCE,1:NDTWCE)
C
                       IND2HR(       1)=1
      IF (ITWCEN.EQ.2) IND2HR(ITWCEN  )=2
      IF (ITWCEN.EQ.2) IND2HR(ITWCEN+1)=1
C
                       IND2HL(       1)=1
      IF (ITWCEN.EQ.2) IND2HL(ITWCEN  )=2
      IF (ITWCEN.EQ.2) IND2HL(ITWCEN+1)=2

C
      ICOUNT=0
C
      DO I=1,ITWCEN
         DO J=1,ITWCEN
            DO K=1,ITWCEN
               DO L=1,ITWCEN
C
                  ICOUNT=ICOUNT+1
C
                  IF (ICOUNT.GT.NDTWHE)
     *                STOP ' INCREASE NDTWHE IN INDGEN'
C
                  IND4HI(ICOUNT)=I
                  IND4HJ(ICOUNT)=J
                  IND4HK(ICOUNT)=K
                  IND4HL(ICOUNT)=L
C
                  INTWHE(I,J,K,L) = ICOUNT
C
               END DO
            END DO
         END DO
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_DEFINT(NOSACT,KARTEZ,NGAUSS,ITWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          PHERMI(ND2MAX+1),DHERMI(ND2MAX+1)
      DIMENSION
     *          XHERMI(NDGAUS),WHERMI(NDGAUS)
C
      COMMON
     *       /TWCINT/ WGHERM(1:NDGAUS,1:NDKART),
     *                PTHERM(1:NDGAUS,1:NDKART)
C
      COMMON
     *       /FACHER/ HERFAC(0:ND2MAX)
      COMMON
     *       /CFIPRI/ NFIPRI
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *       /SCALNG/ HOMSCA(1:NDKART,1:NDTWCE)
      COMMON
     *       /CENPOS/ CENCOR(1:NDKART,1:NDTWCE)
C
      COMMON
     *       /TWCHER/ HERLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DHELE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DDHLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                HERRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DHERI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DDHRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                HERLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
      COMMON
     *       /TWCEXP/ HEREXP(0:ND2MAX,1:NDGAUS,1:NDKART)

C=======================================================================
C
C          THIS SUBROUTINE CALCULATES HARMONIC OSCILLATOR WAVE FUNCTIONS
C          AND THEIR FIRST AND SECOND DERIVATIVES, AS WELL  AS  PREPARES
C          POINTS AND WEIGHTS FOR THE GAUSS-HERMITE INTEGRATION.
C
C          EXPONENTIAL FACTORS ARE    N  O  T   INCLUDED  NEITHER IN THE
C          WAVE FUNCTIONS NOR IN THE DERIVATIVES.
C
           IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_DEFINT'
C
C=======================================================================
C          EVEN IF THE BASIS CONTAINS THE WAVE FUNCTIONS ONLY UP TO  THE
C          NUMBER OF OSCILLATOR QUANTA EQUAL TO "NOSACT", THE CODE NEEDS
C          THE WAVE FUNCTIONS UP TO THE MAXIMUM NUMBER OF QUANTA  EQUAL:
C=======================================================================
C
      NMAIN2=2*NOSACT+2
C
      IF (NGAUSS.GT.NDGAUS) THEN
C
         WRITE(NFIPRI,'(//,1X,''NGAUSS='',I2,'' .GT. NDGAUS='',I2,
     *                         ''   STOP IN DEFINT'',//)') NGAUSS,NDGAUS
         STOP ' NGAUSS.GT.NDGAUS  IN  DEFINT'
C
      END IF
C
C=======================================================================
C     DEFINING ZEROS AND WEIGHTS FOR THE GAUSS-HERMITE  INTEGRATION
C=======================================================================
C
      CALL HERMIT(NGAUSS,XHERMI,WHERMI,NDGAUS)

C
C=======================================================================
C     DEFINING THE HERMITE POLYNOMIALS FOR THE TRANSFORMATION OF THE
C     PRODUCT OF TWO OF THEM INTO A SUM
C=======================================================================
C
      DO IZEROS = 1,NGAUSS
C
         WGHERM(IZEROS,KARTEZ)=WHERMI(IZEROS)
         PTHERM(IZEROS,KARTEZ)=XHERMI(IZEROS)
C
         XZER = XHERMI(IZEROS)
C
         CALL D_HERM(XZER,NMAIN2,PHERMI,DHERMI,ND2MAX+1)
C
         DO NOSCIL = 0,NMAIN2
C
           HEREXP(NOSCIL,IZEROS,KARTEZ)= PHERMI(NOSCIL+1)/
     *                                     HERFAC(NOSCIL)
         END DO
C
      END DO
C
C=======================================================================
C         DEFINING HERMITE POLYNOMIALS WHEN ONLY TWO ARE
C         INVOLVED (DENSITY NON-DEPENDENT HAMILTONIAN TERMS AND NORM
C         NORM OVERLAP MATRIX FOR THE GENERALIZED EIGENVALUE PROBLEM)
C=======================================================================
C
      DO INDBLO=1,NDTWBL
C
          I=IND2HL(INDBLO)
          J=IND2HR(INDBLO)
C
          HOSCAI=HOMSCA(KARTEZ,I)
          HOSCAJ=HOMSCA(KARTEZ,J)
C
          HOSHII=CENCOR(KARTEZ,I)
          HOSHIJ=CENCOR(KARTEZ,J)
C
          SUMOSL=HOSCAI**2+HOSCAJ**2
          SUMCEN=HOSCAI**2*HOSHII+HOSCAJ**2*HOSHIJ
C
          FACTO1=SQRT(2.D0/SUMOSL)
          FACTO2=SUMCEN/SUMOSL
C
          FACNE2(KARTEZ,INDBLO)=           SQRT(2.D0/SUMOSL)
     *      *EXP( -0.5D0*HOSCAI**2*HOSCAJ**2*(HOSHII-HOSHIJ)
     *                                          **2/SUMOSL )
C
C=======================================================================
C     LEFT WAVE FUNCTION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = HOSCAI*(XHERMI(IZEROS)*FACTO1+FACTO2-HOSHII)
C
              CALL D_HERM(XZER,NMAIN2,PHERMI,DHERMI,ND2MAX+1)
C
              DO NOSCIL = 0,NMAIN2
C
                  HERLE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                     PHERMI(NOSCIL+1)
     *                                   *     SQRT(HOSCAI)
     *                                   / HERFAC( NOSCIL )
C
                  DHELE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAI)
     *                                   / HERFAC( NOSCIL )
C
                  DDHLE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAI)
     *                                   / HERFAC( NOSCIL )
C
              END DO
C
C=======================================================================
C
C        DERIVATIVES OF THE "HERMITE * EXPONENT" PRODUCTS  (SEE COMMENTS
C        IN DEVHER)  ARE  MULTIPLIED  BY  THE  APPROPRIATE POWERS OF THE
C        DIMENSIONED STRETCHING FACTORS "HOMSCA". IN THIS WAY THE  FIRST
C        DERIVATIVE HAS THE CORRECT DIMENSION OF FERMI**(-3/2). ALSO THE
C        SECOND DERIVATIVE HAS THE CORRECT DIMENSION OF FERMI**(-5/2).
C
C=======================================================================
C
              DO NOSCIL = 0,NMAIN2
C
                  DHELE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                    DHELE2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                               *HOSCAI
C
                  DDHLE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                   DDHLE2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                        *HOSCAI*HOSCAI
              END DO
C
          END DO
C
C=======================================================================
C     RIGHT WAVE FUNCTION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS

              XZER = HOSCAJ*(XHERMI(IZEROS)*FACTO1+FACTO2-HOSHIJ)
C
              CALL D_HERM(XZER,NMAIN2,PHERMI,DHERMI,ND2MAX+1)
C
              DO NOSCIL = 0,NMAIN2
C
                  HERRI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                     PHERMI(NOSCIL+1)
     *                                   *     SQRT(HOSCAJ)
     *                                   / HERFAC( NOSCIL )
C
                  DHERI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                        (DHERMI(NOSCIL+1)
     *                                 - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAJ)
     *                                   / HERFAC( NOSCIL )
C
                  DDHRI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAJ)
     *                                   / HERFAC( NOSCIL )
C
              END DO
C
              DO NOSCIL = 0,NMAIN2
C
                 DHERI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                    DHERI2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                               *HOSCAJ
C
                 DDHRI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                   DDHRI2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                        *HOSCAJ*HOSCAJ
              END DO
C
          END DO
C
      END DO
C
C=======================================================================
C         DEFINING HERMITE POLYNOMIALS WHEN FOUR OF THEM ARE
C         INVOLVED (DENSITY DEPENDENT HAMILTONIAN TERMS)
C=======================================================================
C
      DO INDBLO=1,NDTWHE
C
          I=IND4HI(INDBLO)
          J=IND4HJ(INDBLO)
          K=IND4HK(INDBLO)
          L=IND4HL(INDBLO)
C
          HOSCAI=HOMSCA(KARTEZ,I)
          HOSCAJ=HOMSCA(KARTEZ,J)
          HOSCAK=HOMSCA(KARTEZ,K)
          HOSCAL=HOMSCA(KARTEZ,L)
C
          HOSHII=CENCOR(KARTEZ,I)
          HOSHIJ=CENCOR(KARTEZ,J)
          HOSHIK=CENCOR(KARTEZ,K)
          HOSHIL=CENCOR(KARTEZ,L)
C
          SUMOSL=HOSCAI**2+HOSCAJ**2+HOSCAK**2+HOSCAL**2
          SUMCEN=HOSCAI**2*HOSHII+HOSCAJ**2*HOSHIJ
     *          +HOSCAK**2*HOSHIK+HOSCAL**2*HOSHIL
C
          FACLE1=HOSCAJ*SQRT(2.D0/SUMOSL)
          FACLE2=HOSCAJ*(SUMCEN/SUMOSL-HOSHIJ)
C
          FACRI1=HOSCAI*SQRT(2.D0/SUMOSL)
          FACRI2=HOSCAI*(SUMCEN/SUMOSL-HOSHII)
C
C=======================================================================
C     LEFT WAVE FUNCTION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = XHERMI(IZEROS)*FACLE1+FACLE2
C
              CALL D_HERM(XZER,NMAIN2,PHERMI,DHERMI,ND2MAX+1)
C
              DO NOSCIL = 0,NMAIN2
C
                   HERLE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                     PHERMI(NOSCIL+1)
     *                                   *     SQRT(HOSCAJ)
     *                                   / HERFAC( NOSCIL )
C
                  DHELE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAJ)
     *                                   / HERFAC( NOSCIL )
C
                  DDHLE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAJ)
     *                                   / HERFAC( NOSCIL )
C
              END DO
C
C=======================================================================
C
C        DERIVATIVES OF THE "HERMITE * EXPONENT" PRODUCTS  (SEE COMMENTS
C        IN DEVHER)  ARE  MULTIPLIED  BY  THE  APPROPRIATE POWERS OF THE
C        DIMENSIONED STRETCHING FACTORS "HOMSCA". IN THIS WAY THE  FIRST
C        DERIVATIVE HAS THE CORRECT DIMENSION OF FERMI**(-3/2). ALSO THE
C        SECOND DERIVATIVE HAS THE CORRECT DIMENSION OF FERMI**(-5/2).
C
C=======================================================================
C
              DO NOSCIL = 0,NMAIN2
C
                  DHELE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                    DHELE2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                               *HOSCAJ
C
                  DDHLE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                   DDHLE2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                        *HOSCAJ*HOSCAJ
              END DO
C
          END DO
C
C=======================================================================
C     RIGHT WAVE FUNCTION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS

              XZER = XHERMI(IZEROS)*FACRI1+FACRI2
C
              CALL D_HERM(XZER,NMAIN2,PHERMI,DHERMI,ND2MAX+1)
C
              DO NOSCIL = 0,NMAIN2
C
                  HERRI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                     PHERMI(NOSCIL+1)
     *                                   *     SQRT(HOSCAI)
     *                                   / HERFAC( NOSCIL )
C
                  DHERI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAI)
     *                                   / HERFAC( NOSCIL )
C
                  DDHRI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAI)
     *                                   / HERFAC( NOSCIL )
C
              END DO
C
              DO NOSCIL = 0,NMAIN2
C
                 DHERI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                    DHERI2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                               *HOSCAI
C
                 DDHRI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                   DDHRI2(NOSCIL,IZEROS,KARTEZ,INDBLO)
     *                                        *HOSCAI*HOSCAI
              END DO
C
          END DO
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_COPHER(NOSACT,KARTEZ,NGAUSS,ITWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMPLEX
     *          CHELE2,CDHLE2,CDDLE2,CHERI2,CDHRI2,CDDRI2,
     *          CHELE4,CDHLE4,CDDLE4,CHERI4,CDHRI4,CDDRI4
C
      COMMON
     *       /TWCHER/ HERLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DHELE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DDHLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                HERRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DHERI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                DDHRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                HERLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE)
      COMMON
     *       /CTWCHE/ CHELE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                CDHLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                CDDLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                CHERI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                CDHRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                CDDRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
     *                CHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                CDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                CDDLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                CHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                CDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
     *                CDDRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE)
C
C=======================================================================
C
C          THIS SUBROUTINE COPIES THE HARMONIC OSCILLATOR WAVE FUNCTIONS
C          AND THEIR FIRST AND SECOND DERIVATIVES TO COMPLEX ARRAYS.
C
C=======================================================================
C
           IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_COPHER'
C
C=======================================================================
C
      NMAIN2=2*NOSACT+2
C
      DO INDBLO = 1, NDTWBL
         DO IZEROS = 1,NGAUSS
            DO NOSCIL = 0,NMAIN2
C
               CHELE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                HERLE2(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDHLE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DHELE2(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDDLE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DDHLE2(NOSCIL,IZEROS,KARTEZ,INDBLO)
C
               CHERI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                HERRI2(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDHRI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DHERI2(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDDRI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DDHRI2(NOSCIL,IZEROS,KARTEZ,INDBLO)
C
            END DO
         END DO
      END DO
C
      DO INDBLO = 1, NDTWHE
         DO IZEROS = 1,NGAUSS
            DO NOSCIL = 0,NMAIN2
C
               CHELE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                HERLE4(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDHLE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DHELE4(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDDLE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DDHLE4(NOSCIL,IZEROS,KARTEZ,INDBLO)
C
               CHERI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                HERRI4(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDHRI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DHERI4(NOSCIL,IZEROS,KARTEZ,INDBLO)
               CDDRI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                DDHRI4(NOSCIL,IZEROS,KARTEZ,INDBLO)
C
            END DO
         END DO
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
      SUBROUTINE TWCCOP(ISIMPY,LTWCEN)
C=======================================================================
      USE MAT_PP
      USE MAT_PM
      USE TWCEPP
      USE TWCEPM
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMMON
     *       /DIMENS/ LDBASE
C
C=======================================================================
C
      IALLOC=0
C
C=======================================================================
      IF (.NOT.ALLOCATED(TWC_PP)) THEN
          ALLOCATE (TWC_PP(1:2*NDBASE,1:2*NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWC_PP','TWCCOP')
      END IF
C=======================================================================
      IF (.NOT.ALLOCATED(TWC_PM).AND.LTWCEN.GE.1.AND.ISIMPY.NE.1) THEN
          ALLOCATE (TWC_PM(1:2*NDBASE,1:2*NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWC_PM','TWCCOP')
      END IF
C=======================================================================
C
C=======================================================================
C          THIS SUBROUTINE COPIES THE CALCULATED MATRIX ELEMENTS OF THE
C          MEAN-FIELD INTO THE TWO-CENTER HAMILTONIAN
C=======================================================================
C
      IF (LTWCEN.EQ.1) THEN
C
C=======================================================================
C          ATTENTION: THE TWO LINES BELOW SHOULD ONLY BE INCLUDED IF THE
C                     TWO-CENTER CROSS-TERM IS TO BE SET TO ZERO.
C=======================================================================
C
          TWC_PP=CMPLX(0.0D0,0.0D0)
          TWC_PM=CMPLX(0.0D0,0.0D0)
C
C=======================================================================
C
          MBRA=0
          MKET=0
C
      END IF
C
      IF (LTWCEN.EQ.2) THEN
C
          MBRA=LDBASE
          MKET=LDBASE
C
      END IF
C
      IF (LTWCEN.EQ.3) THEN
C
          MBRA=LDBASE
          MKET=0
C
      END IF
C
C=======================================================================
C     SETTING THE MATRIX ELEMENTS
C=======================================================================
C
      DO IKET=1,LDBASE
         DO IBRA=1,LDBASE
C
            TWC_PP(MBRA+IBRA,MKET+IKET,0)=BIG_PP(IBRA,IKET,0)
            TWC_PP(MBRA+IBRA,MKET+IKET,1)=BIG_PP(IBRA,IKET,1)
C
         END DO
      END DO
C

      IF (ISIMPY.NE.1) THEN
C
          DO IKET=1,LDBASE
             DO IBRA=1,LDBASE
C
                TWC_PM(MBRA+IBRA,MKET+IKET,0)=BIG_PM(IBRA,IKET,0)
                TWC_PM(MBRA+IBRA,MKET+IKET,1)=BIG_PM(IBRA,IKET,1)
C
             END DO
          END DO
C
      END IF
C
C=======================================================================
C     FOR THE INITIAL STEP, THE NILSSON HAMILTONIAN ONLY HAS DIAGONAL
C     BLOCKS, SO WE PUT THE OFF-DIAGONAL TERMS TO ZERO.
C=======================================================================
C
      IF (LTWCEN.EQ.3) THEN
C
          MBRA=0
          MKET=LDBASE
C
          DO IKET=1,LDBASE
             DO IBRA=1,LDBASE
C
                TWC_PP(MBRA+IBRA,MKET+IKET,0)=CONJG(BIG_PP(IKET,IBRA,0))
                TWC_PP(MBRA+IBRA,MKET+IKET,1)=CONJG(BIG_PP(IKET,IBRA,1))
C
             END DO
          END DO
C
          IF (ISIMPY.NE.1) THEN
C
              DO IKET=1,LDBASE
                 DO IBRA=1,LDBASE
C
                TWC_PM(MBRA+IBRA,MKET+IKET,0)=CONJG(BIG_PM(IKET,IBRA,0))
                TWC_PM(MBRA+IBRA,MKET+IKET,1)=CONJG(BIG_PM(IKET,IBRA,1))
C
                 END DO
              END DO
C
          END IF
C
      END IF
C
C=======================================================================
C
      RETURN
      END

      SUBROUTINE TWC_DEVHER(NOSACT,KARTEZ,ITWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
C
      DIMENSION
     *          PHERMI(2*NOSACT+7),DHERMI(2*NOSACT+7)
      DIMENSION
     *          XHERMI(2*NOSACT+7),WHERMI(2*NOSACT+7)
      DIMENSION
     *          HERLEF(0:NOSACT,1:2*NOSACT+7),
     *          DH1LEF(0:NOSACT,1:2*NOSACT+7),
     *          DH2LEF(0:NOSACT,1:2*NOSACT+7)
      DIMENSION
     *          HERRIG(0:NOSACT,1:2*NOSACT+7),
     *          DH2RIG(0:NOSACT,1:2*NOSACT+7),
     *          DH1RIG(0:NOSACT,1:2*NOSACT+7)
      DIMENSION
     *          HERFAC(0:2*NOSACT+7)
      DIMENSION
     *           HERAUX(0:2*NOSACT+6,1:2*NOSACT+7)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
      COMMON
     *       /COEFH4/ CF4H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWHE),
     *                CF4H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWHE),
     *                CF4H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWHE),
     *                CF4H11(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWHE)

      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)

      COMMON
     *       /SCALNG/ HOMSCA(1:NDKART,1:NDTWCE)
      COMMON
     *       /CFIPRI/ NFIPRI
      COMMON
     *       /CENPOS/ CENCOR(1:NDKART,1:NDTWCE)
      COMMON
     *       /TWCFAC/ QUADSUM(1:NDKART,1:NDTWBL),
     *                QUADCEN(1:NDKART,1:NDTWBL)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
C=======================================================================
C
      IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_DEVHER'
C
C=======================================================================
C
      NMAIN1=2*NOSACT+6
      NGAUSS=2*NOSACT+7
C
C=======================================================================
C
      DO INDBLO=1,NDTWBL
C
          I=IND2HL(INDBLO)
          J=IND2HR(INDBLO)
C
          HOSCAI=HOMSCA(KARTEZ,I)
          HOSCAJ=HOMSCA(KARTEZ,J)
C
          HOSHII=CENCOR(KARTEZ,I)
          HOSHIJ=CENCOR(KARTEZ,J)
C
C=======================================================================
C     DEFINING ZEROS AND WEIGHTS FOR THE GAUSS-HERMITE  INTEGRATION
C=======================================================================
C
          CALL HERMIT(NGAUSS,XHERMI,WHERMI,NGAUSS)
C
C=======================================================================
C     DEFINING HERMITE POLYNOMIAL NORMALIZATION CONSTANTS AND OTHERS
C=======================================================================
C
          CALL NORFAC(HERFAC,NMAIN1)
C
          SUMOSL=HOSCAI**2+HOSCAJ**2
          SUMCEN=HOSCAI**2*HOSHII+HOSCAJ**2*HOSHIJ
C
          FACTO1=SQRT(2.D0/SUMOSL)
          FACTO2=SUMCEN/SUMOSL
C
          QUADSUM(KARTEZ,INDBLO)=SQRT(2.D0/SUMOSL)
          QUADCEN(KARTEZ,INDBLO)=SUMCEN/SUMOSL
C
C=======================================================================
C     DEFINING VALUES OF THE LEFT HERMITE POLYNOMIALS AT GAUSS ZEROS
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = HOSCAI*(XHERMI(IZEROS)*FACTO1+FACTO2-HOSHII)
C
              CALL D_HERM(XZER,NMAIN1,PHERMI,DHERMI,ND1MAX+1)
C
              DO NOSCIL = 0,NOSACT
C
                  HERLEF(NOSCIL,IZEROS)= PHERMI(NOSCIL+1)/HERFAC(NOSCIL)
C
                  DH1LEF(NOSCIL,IZEROS) = (DHERMI(NOSCIL+1)
     *                             - XZER*PHERMI(NOSCIL+1))
     *                                   / HERFAC( NOSCIL )
C
                  DH2LEF(NOSCIL,IZEROS) = (XZER*XZER-2*NOSCIL-1)
     *                                        * PHERMI(NOSCIL+1)
     *                                        / HERFAC( NOSCIL )
C
              END DO
C
          END DO
C
C=======================================================================
C     DEFINING VALUES OF THE RIGHT HERMITE POLYNOMIALS AT GAUSS ZEROS
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = HOSCAJ*(XHERMI(IZEROS)*FACTO1+FACTO2-HOSHIJ)
C
              CALL D_HERM(XZER,NMAIN1,PHERMI,DHERMI,ND1MAX+1)
C
              DO NOSCIL = 0,NOSACT
C
                  HERRIG(NOSCIL,IZEROS)= PHERMI(NOSCIL+1)/HERFAC(NOSCIL)
C
                  DH1RIG(NOSCIL,IZEROS) = (DHERMI(NOSCIL+1)
     *                             - XZER*PHERMI(NOSCIL+1))
     *                                   / HERFAC( NOSCIL )
C
                  DH2RIG(NOSCIL,IZEROS) = (XZER*XZER-2*NOSCIL-1)
     *                                        * PHERMI(NOSCIL+1)
     *                                        / HERFAC( NOSCIL )
C
              END DO
C
          END DO
C
C=======================================================================
C        DEFINING THE NEW HERMITE POLYNOMIALS OF THE EXPANSION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = XHERMI(IZEROS)
C
              CALL D_HERM(XZER,NMAIN1,PHERMI,DHERMI,ND1MAX+1)
C
              DO NOSCIL = 0,NMAIN1
C
                  HERAUX(NOSCIL,IZEROS)= PHERMI(NOSCIL+1)/
     *                                     HERFAC(NOSCIL)
C
              END DO
C
          END DO
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 00
C=======================================================================
C
          DO M=0,NOSACT
              DO N=0,NOSACT
                  DO K=0,N+M
C
                      RESULT=0.0D0
C
                      DO IGAUSS=1,NGAUSS
C
                          RESULT=RESULT+WHERMI(IGAUSS)
     *                               *HERLEF(M,IGAUSS)
     *                               *HERRIG(N,IGAUSS)
     *                               *HERAUX(K,IGAUSS)
C
                      END DO
C
                      CF2H00(K,M,N,KARTEZ,INDBLO) = RESULT
     *                               * SQRT(HOSCAI*HOSCAJ)
C
                  END DO
              END DO
          END DO
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 01
C=======================================================================
C
          DO M=0,NOSACT
              DO N=0,NOSACT
                  DO K=0,N+M
C
                      RESULT=0.0D0
C
                      DO IGAUSS=1,NGAUSS
C
                          RESULT=RESULT+WHERMI(IGAUSS)
     *                               *HERLEF(M,IGAUSS)
     *                               *DH1RIG(N,IGAUSS)
     *                               *HERAUX(K,IGAUSS)
C
                      END DO
C
                      CF2H01(K,M,N,KARTEZ,INDBLO) = RESULT
     *                               * SQRT(HOSCAI*HOSCAJ)
     *                               *             HOSCAJ
C
                  END DO
              END DO
          END DO
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 02
C=======================================================================
C
          DO M=0,NOSACT
              DO N=0,NOSACT
                  DO K=0,N+M
C
                      RESULT=0.0D0
C
                      DO IGAUSS=1,NGAUSS
C
                          RESULT=RESULT+WHERMI(IGAUSS)
     *                               *HERLEF(M,IGAUSS)
     *                               *DH2RIG(N,IGAUSS)
     *                               *HERAUX(K,IGAUSS)
C
                      END DO
C
                      CF2H02(K,M,N,KARTEZ,INDBLO) = RESULT
     *                               * SQRT(HOSCAI*HOSCAJ)
     *                               *      HOSCAJ*HOSCAJ
C
                  END DO
              END DO
          END DO
C
      END DO
C=======================================================================
C        CALCULATING THE COEFFICIENTS WHEN 4 WAVE FUNCTIONS ARE USED
C=======================================================================
      DO INDBLO=1,NDTWHE
C
          I=IND4HI(INDBLO)
          J=IND4HJ(INDBLO)
          K=IND4HK(INDBLO)
          L=IND4HL(INDBLO)
C
          HOSCAI=HOMSCA(KARTEZ,I)
          HOSCAJ=HOMSCA(KARTEZ,J)
          HOSCAK=HOMSCA(KARTEZ,K)
          HOSCAL=HOMSCA(KARTEZ,L)
C
          HOSHII=CENCOR(KARTEZ,I)
          HOSHIJ=CENCOR(KARTEZ,J)
          HOSHIK=CENCOR(KARTEZ,K)
          HOSHIL=CENCOR(KARTEZ,L)
C
C=======================================================================
C     DEFINING VALUES OF THE LEFT HERMITE POLYNOMIALS AT GAUSS ZEROS
C=======================================================================
C
          SUMOSL=HOSCAI**2+HOSCAJ**2+HOSCAK**2+HOSCAL**2
          SUMCEN=HOSCAI**2*HOSHII+HOSCAJ**2*HOSHIJ
     *          +HOSCAK**2*HOSHIK+HOSCAL**2*HOSHIL
C
          FACLE1=HOSCAJ*SQRT(2.D0/SUMOSL)
          FACLE2=HOSCAJ*(SUMCEN/SUMOSL-HOSHIJ)
C
          FACRI1=HOSCAI*SQRT(2.D0/SUMOSL)
          FACRI2=HOSCAI*(SUMCEN/SUMOSL-HOSHII)
C
          FACNE4(KARTEZ,INDBLO)=         SQRT(2.D0/SUMOSL)
     *   * EXP( -0.5D0*((HOSCAI*HOSHII)**2+(HOSCAJ*HOSHIJ)**2
     *                 +(HOSCAK*HOSHIK)**2+(HOSCAL*HOSHIL)**2
     *                                   - SUMCEN**2/SUMOSL))
C
          DO IZEROS = 1,NGAUSS
C
              XZER = XHERMI(IZEROS)*FACLE1+FACLE2
C
              CALL D_HERM(XZER,NMAIN1,PHERMI,DHERMI,
     *                                     ND1MAX+1)
C
              DO NOSCIL = 0,NOSACT
C
                  HERLEF(NOSCIL,IZEROS)= PHERMI(NOSCIL+1)/
     *                                     HERFAC(NOSCIL)
C
                  DH1LEF(NOSCIL,IZEROS) = (DHERMI(NOSCIL+1)
     *                             - XZER*PHERMI(NOSCIL+1))
     *                                   / HERFAC( NOSCIL )
C
                  DH2LEF(NOSCIL,IZEROS) = (XZER*XZER-2*NOSCIL-1)
     *                                        * PHERMI(NOSCIL+1)
     *                                        / HERFAC( NOSCIL )
              END DO
C
          END DO
C
C=======================================================================
C        DEFINING VALUES OF THE RIGHT HERMITE POLYNOMIALS AT GAUSS ZEROS
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = XHERMI(IZEROS)*FACRI1+FACRI2
C
              CALL D_HERM(XZER,NMAIN1,PHERMI,DHERMI,
     *                                           ND1MAX+1)
C
              DO NOSCIL = 0,NOSACT
C
                  HERRIG(NOSCIL,IZEROS)= PHERMI(NOSCIL+1)/
     *                                           HERFAC(NOSCIL)
C
                  DH1RIG(NOSCIL,IZEROS) = (DHERMI(NOSCIL+1)
     *                             - XZER*PHERMI(NOSCIL+1))
     *                                   / HERFAC( NOSCIL )
C
                  DH2RIG(NOSCIL,IZEROS) = (XZER*XZER-2*NOSCIL-1)
     *                                        * PHERMI(NOSCIL+1)
     *                                        / HERFAC( NOSCIL )
C
              END DO
C
          END DO
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 00
C=======================================================================
C
          DO N=0,NOSACT
              DO M=0,NOSACT
                  DO K=0,N+M
C
                      RESULT=0.0D0
C
                      DO IGAUSS=1,NGAUSS
C
                          RESULT=RESULT+WHERMI(IGAUSS)
     *                               *HERLEF(M,IGAUSS)
     *                               *HERRIG(N,IGAUSS)
     *                               *HERAUX(K,IGAUSS)
C
                      END DO
C
                      CF4H00(K,M,N,KARTEZ,INDBLO)=RESULT
     *                              * SQRT(HOSCAI*HOSCAJ)
     *                            * FACNE4(KARTEZ,INDBLO)
C
                  END DO
              END DO
          END DO
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 01
C=======================================================================
C
          DO N=0,NOSACT
              DO M=0,NOSACT
                  DO K=0,N+M
C
                      RESULT=0.0D0
C
                      DO IGAUSS=1,NGAUSS
C
                          RESULT=RESULT+WHERMI(IGAUSS)
     *                               *HERLEF(M,IGAUSS)
     *                               *DH1RIG(N,IGAUSS)
     *                               *HERAUX(K,IGAUSS)
C
                      END DO
C
                      CF4H01(K,M,N,KARTEZ,INDBLO)=RESULT
     *                            *   SQRT(HOSCAI*HOSCAJ)
     *                            *               HOSCAJ
     *                            * FACNE4(KARTEZ,INDBLO)
C
                  END DO
              END DO
          END DO
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 11
C=======================================================================
C
          DO N=0,NOSACT
              DO M=0,NOSACT
                  DO K=0,N+M
C
                      RESULT=0.0D0
C
                      DO IGAUSS=1,NGAUSS
C
                          RESULT=RESULT+WHERMI(IGAUSS)
     *                               *DH1LEF(M,IGAUSS)
     *                               *DH1RIG(N,IGAUSS)
     *                               *HERAUX(K,IGAUSS)
C
                      END DO
C
                      CF4H11(K,M,N,KARTEZ,INDBLO)=RESULT
     *                              * SQRT(HOSCAI*HOSCAJ)
     *                              *      HOSCAI*HOSCAJ
     *                            * FACNE4(KARTEZ,INDBLO)
C
                  END DO
              END DO
          END DO
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C=======================================================================
C
      SUBROUTINE NORMAT(IEIGCU,EIGCUT)
C=======================================================================
      USE hfodd_sizes
      USE OVM_PP
      USE OVM_PM
      USE NOREIG
C=======================================================================
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /AUXTWC/ AUXTW0(1:NDBASE,1:NDBASE),
     *                AUXTW1(1:NDBASE,1:NDBASE)
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
      COMMON
     *       /NATCUT/ ICUTOF
      COMMON
     *       /EGVOVM/ EGVNOR(1:4*NDBASE)
      COMPLEX
     *          CELNOR
C=======================================================================
      ALLOCATABLE CELNOR(:)
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
C
C=======================================================================
C
      AUXTW0(:,:)=0.D0
C
      DO LTWCEN=1,NDTWBL
C
         FACTOR=PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
         DO IBRA=1,LDBASE
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBASE
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               RESULT =                    FACTOR
     *                  * CF2H00(0,MX,NX,1,LTWCEN)
     *                  * CF2H00(0,MY,NY,2,LTWCEN)
     *                  * CF2H00(0,MZ,NZ,3,LTWCEN)
C
                   IF(LTWCEN.EQ.1) AUXTW0(IBRA,IKET)=RESULT
                   IF(LTWCEN.EQ.2) AUXTW0(IBRA+LDBASE,IKET+LDBASE)
     *                             = RESULT
                   IF(LTWCEN.EQ.3) AUXTW0(IBRA+LDBASE,IKET       )
     *                             = RESULT
                   IF(LTWCEN.EQ.3) AUXTW0(IKET,IBRA+LDBASE       )
     *                             = RESULT
            END DO
         END DO
C
      END DO
C=======================================================================
      IF (.NOT.ALLOCATED(OVMAPP)) THEN
          ALLOCATE (OVMAPP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('OVMAPP','NORMAT')
      END IF
      IF (.NOT.ALLOCATED(OVMAPM)) THEN
          ALLOCATE (OVMAPM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('OVMAPM','NORMAT')
      END IF
C=======================================================================
      DO IBRA=1,2*LDBASE
         DO IKET=1,2*LDBASE
C
          IF (IBRA.LE.LDBASE) ICENBR=1
          IF (IBRA.GT.LDBASE) ICENBR=NDTWCE
          IF (IKET.LE.LDBASE) ICENKE=1
          IF (IKET.GT.LDBASE) ICENKE=NDTWCE
C
          MY=NYVECT(IBRA,IND2HL(ICENBR))
          NY=NYVECT(IKET,IND2HR(ICENKE))
C
          OVMAPP(IBRA,IKET,0)=AUXTW0(IBRA,IKET)*IPHAPP(MY,NY,0)
          OVMAPP(IBRA,IKET,1)=AUXTW0(IBRA,IKET)*IPHAMM(MY,NY,0)
          OVMAPM(IBRA,IKET,0)=0.D0
          OVMAPM(IBRA,IKET,1)=0.D0
C
         END DO
      END DO
C=======================================================================
C     DIAGONALIZING THE OVERLAP MATRIX
C=======================================================================
      IF (.NOT.ALLOCATED(CELNOR)) THEN
          ALLOCATE (CELNOR(1:(4*NDBASE+1)*2*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('CELNOR','NORMAT')
      END IF
C
      IF (.NOT.ALLOCATED(EWTNOR)) THEN
          ALLOCATE (EWTNOR(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('EWTNOR','NORMAT')
      END IF
      IF (.NOT.ALLOCATED(EWTAUX)) THEN
          ALLOCATE (EWTAUX(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('EWTAUX','NORMAT')
      END IF
C=======================================================================
      DO IBRA=1,4*LDBASE
         DO IKET=1,4*LDBASE
C
            IF (IBRA.GE.IKET) THEN
C
                NCOUNT=IBRA+((8*LDBASE-IKET)*(IKET-1))/2
C
                IF (IBRA.LE.2*LDBASE.AND.IKET.LE.2*LDBASE)
     *
     *              CELNOR(NCOUNT)=OVMAPP(IBRA,IKET,0)
C
                IF (IBRA.GT.2*LDBASE.AND.IKET.LE.2*LDBASE)
     *
     *              CELNOR(NCOUNT)=OVMAPM(IBRA-2*LDBASE,IKET,1)
C
                IF (IBRA.GT.2*LDBASE.AND.IKET.GT.2*LDBASE)
     *
     *              CELNOR(NCOUNT)=OVMAPP(IBRA-2*LDBASE,IKET-2*LDBASE,1)
C
            END IF
C
         END DO
      END DO
C
      CALL DIAMAT(CELNOR,EGVNOR,EWTNOR,4*LDBASE,4*NDBASE,4*LDBASE)
C
C=======================================================================
C     INDEX FOR THE FIRST STATE TO BE TAKEN INTO ACCOUNT
C=======================================================================
C
      ICUTOF=0
C
      IF (IEIGCU.EQ.0) EIGCUT=-EGVNOR(1)
C
      DO IBAS=4*LDBASE,1,-1
C
          IF(EGVNOR(IBAS).GT.EIGCUT) ICUTOF=IBAS
C
      END DO
C
      ICOUNT=0
C
      DO INDNAT=ICUTOF,4*LDBASE
      ICOUNT=ICOUNT+1
         DO IBASE=1,4*LDBASE
         EWTAUX(IBASE,ICOUNT)=EWTNOR(IBASE,INDNAT)/SQRT(EGVNOR(INDNAT))
         END DO
      END DO
C
C=======================================================================
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_PX(IPHASE,HAMLIN,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          IPHASE(0:NDYMAX,0:NDYMAX)
      DIMENSION
     *          HAMLIN(1:NDBASE,1:NDBASE)
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
C
C=======================================================================
C
      HAMLIN(:,:)=0.D0
C
      DO LTWCEN=1,NDTWBL
      FACTOR=-PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
         DO IBRA=1,LDBASE
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBASE
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMLIN(IBRA,IKET)= FACINP * IPHASE(NY,MY) * FACTOR
     *                                  * CF2H01(0,MX,NX,1,LTWCEN)
     *                                  * CF2H00(0,MY,NY,2,LTWCEN)
     *                                  * CF2H00(0,MZ,NZ,3,LTWCEN)
C
            END DO
C
         END DO
C
      END DO
C=======================================================================
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_PY(IPHASE,HAMLIN,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          IPHASE(0:NDYMAX,0:NDYMAX)
      DIMENSION
     *          HAMLIN(1:NDBASE,1:NDBASE)
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
C
C=======================================================================
C
      HAMLIN(:,:)=0.D0
C
      DO LTWCEN=1,NDTWBL
      FACTOR=-PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
         DO IBRA=1,LDBASE
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBASE
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMLIN(IBRA,IKET)= FACINP * IPHASE(NY,MY) * FACTOR
     *                                  * CF2H00(0,MX,NX,1,LTWCEN)
     *                                  * CF2H01(0,MY,NY,2,LTWCEN)
     *                                  * CF2H00(0,MZ,NZ,3,LTWCEN)
C
            END DO
C
         END DO
C
      END DO
C=======================================================================
C
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_PZ(IPHASE,HAMLIN,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          IPHASE(0:NDYMAX,0:NDYMAX)
      DIMENSION
     *          HAMLIN(1:NDBASE,1:NDBASE)
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
C
C=======================================================================
C
      HAMLIN(:,:)=0.D0
C
      DO LTWCEN=1,NDTWBL
      FACTOR=-PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
         DO IBRA=1,LDBASE
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBASE
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMLIN(IBRA,IKET)= FACINP * IPHASE(NY,MY) * FACTOR
     *                                  * CF2H00(0,MX,NX,1,LTWCEN)
     *                                  * CF2H00(0,MY,NY,2,LTWCEN)
     *                                  * CF2H01(0,MZ,NZ,3,LTWCEN)
C
            END DO
C
         END DO
C
      END DO
C=======================================================================
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_INTKIN(HAMKIN,HBMACT,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMKIN(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
C
C=======================================================================
C
      HAMKIN(:,:)=0.D0
C
      FACTOR=-PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
         DO IBRA=1,LDBASE
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBASE
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMKIN(IBRA,IKET) = HBMACT * FACTOR
     *                           * IPHAPP(MY,NY,0)
     *       * (
     *           CF2H02(0,MX,NX,1,LTWCEN)
     *       *   CF2H00(0,MY,NY,2,LTWCEN)
     *       *   CF2H00(0,MZ,NZ,3,LTWCEN)
     *
     *       +   CF2H00(0,MX,NX,1,LTWCEN)
     *       *   CF2H02(0,MY,NY,2,LTWCEN)
     *       *   CF2H00(0,MZ,NZ,3,LTWCEN)
     *
     *       +   CF2H00(0,MX,NX,1,LTWCEN)
     *       *   CF2H00(0,MY,NY,2,LTWCEN)
     *       *   CF2H02(0,MZ,NZ,3,LTWCEN)
     *          )
C
            END DO
C
         END DO
C
C=======================================================================
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_LX(IPHASE,HAMANG,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
C
      DIMENSION
     *          IPHASE(0:NDYMAX,0:NDYMAX)
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
      COMMON
     *       /TWCFAC/ QUADSUM(1:NDKART,1:NDTWBL),
     *                QUADCEN(1:NDKART,1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=-PIFAC*FACNE2(1,NDTWBL)*FACNE2(2,NDTWBL)*FACNE2(3,NDTWBL)
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBASE
C
         MX=NXVECT(IBRA,IND2HL(NDTWBL))
         MY=NYVECT(IBRA,IND2HL(NDTWBL))
         MZ=NZVECT(IBRA,IND2HL(NDTWBL))
C
         DO IKET=1,LDBASE
C
            NX=NXVECT(IKET,IND2HR(NDTWBL))
            NY=NYVECT(IKET,IND2HR(NDTWBL))
            NZ=NZVECT(IKET,IND2HR(NDTWBL))
C
            HAMANG(IBRA,IKET) = FACINP * FACTOR * IPHASE(MY,NY) *
     *                                 CF2H00(0,MX,NX,1,NDTWBL) *
     *      (
     *        CF2H01(0,MZ,NZ,3,NDTWBL) *
     *      ( CF2H00(1,MY,NY,2,NDTWBL)*QUADSUM(2,NDTWBL)/SQRT(2.D0) +
     *                 CF2H00(0,MY,NY,2,NDTWBL)*QUADCEN(2,NDTWBL) )
     *      - CF2H01(0,MY,NY,2,NDTWBL) *
     *      ( CF2H00(1,MZ,NZ,3,NDTWBL)*QUADSUM(3,NDTWBL)/SQRT(2.D0) +
     *        CF2H00(0,MZ,NZ,3,NDTWBL)*QUADCEN(3,NDTWBL) )
     *      )
C
         END DO
C
      END DO
C
C=======================================================================
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_LY(IPHASE,HAMANG,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          IPHASE(0:NDYMAX,0:NDYMAX)
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
      COMMON
     *       /TWCFAC/ QUADSUM(1:NDKART,1:NDTWBL),
     *                QUADCEN(1:NDKART,1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=-PIFAC*FACNE2(1,NDTWBL)*FACNE2(2,NDTWBL)*FACNE2(3,NDTWBL)
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBASE
C
         MX=NXVECT(IBRA,IND2HL(NDTWBL))
         MY=NYVECT(IBRA,IND2HL(NDTWBL))
         MZ=NZVECT(IBRA,IND2HL(NDTWBL))
C
         DO IKET=1,LDBASE
C
            NX=NXVECT(IKET,IND2HR(NDTWBL))
            NY=NYVECT(IKET,IND2HR(NDTWBL))
            NZ=NZVECT(IKET,IND2HR(NDTWBL))
C
            HAMANG(IBRA,IKET) = FACINP * FACTOR * IPHASE(MY,NY) *
     *                                 CF2H00(0,MY,NY,2,NDTWBL) *
     *      (
     *        CF2H01(0,MX,NX,1,NDTWBL) *
     *      ( CF2H00(1,MZ,NZ,3,NDTWBL)*QUADSUM(3,NDTWBL)/SQRT(2.D0) +
     *                 CF2H00(0,MZ,NZ,3,NDTWBL)*QUADCEN(3,NDTWBL) )
     *      - CF2H01(0,MZ,NZ,3,NDTWBL) *
     *      ( CF2H00(1,MX,NX,1,NDTWBL)*QUADSUM(1,NDTWBL)/SQRT(2.D0) +
     *        CF2H00(0,MX,NX,1,NDTWBL)*QUADCEN(1,NDTWBL) )
     *      )
C
         END DO
C
      END DO
C
C=======================================================================
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_LZ(IPHASE,HAMANG,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          IPHASE(0:NDYMAX,0:NDYMAX)
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
      COMMON
     *       /TWCFAC/ QUADSUM(1:NDKART,1:NDTWBL),
     *                QUADCEN(1:NDKART,1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=-PIFAC*FACNE2(1,NDTWBL)*FACNE2(2,NDTWBL)*FACNE2(3,NDTWBL)
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBASE
C
         MX=NXVECT(IBRA,IND2HL(NDTWBL))
         MY=NYVECT(IBRA,IND2HL(NDTWBL))
         MZ=NZVECT(IBRA,IND2HL(NDTWBL))
C
         DO IKET=1,LDBASE
C
            NX=NXVECT(IKET,IND2HR(NDTWBL))
            NY=NYVECT(IKET,IND2HR(NDTWBL))
            NZ=NZVECT(IKET,IND2HR(NDTWBL))
C
            HAMANG(IBRA,IKET) = FACINP * FACTOR * IPHASE(MY,NY) *
     *                                 CF2H00(0,MZ,NZ,3,NDTWBL) *
     *      (
     *        CF2H01(0,MY,NY,2,NDTWBL) *
     *      ( CF2H00(1,MX,NX,1,NDTWBL)*QUADSUM(1,NDTWBL)/SQRT(2.D0) +
     *                 CF2H00(0,MX,NX,1,NDTWBL)*QUADCEN(1,NDTWBL) )
     *      - CF2H01(0,MX,NX,1,NDTWBL) *
     *      ( CF2H00(1,MY,NY,2,NDTWBL)*QUADSUM(2,NDTWBL)/SQRT(2.D0) +
     *        CF2H00(0,MY,NY,2,NDTWBL)*QUADCEN(2,NDTWBL) )
     *      )
C
            END DO
         END DO
C
C=======================================================================
      END SUBROUTINE
C=======================================================================
C

      SUBROUTINE TWC_SX(HAMANG,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=PIFAC*FACNE2(1,NDTWBL)*FACNE2(2,NDTWBL)*FACNE2(3,NDTWBL)*
     *       FACINP*0.5D0
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBASE
C
         MX=NXVECT(IBRA,IND2HL(NDTWBL))
         MY=NYVECT(IBRA,IND2HL(NDTWBL))
         MZ=NZVECT(IBRA,IND2HL(NDTWBL))
C
         DO IKET=1,LDBASE
C
            NX=NXVECT(IKET,IND2HR(NDTWBL))
            NY=NYVECT(IKET,IND2HR(NDTWBL))
            NZ=NZVECT(IKET,IND2HR(NDTWBL))
C
            HAMANG(IBRA,IKET) = IPHAPM(MY,NY,1) * FACTOR
     *                        * CF2H00(0,MX,NX,1,NDTWBL)
     *                        * CF2H00(0,MY,NY,2,NDTWBL)
     *                        * CF2H00(0,MZ,NZ,3,NDTWBL)
C
            END DO
         END DO
C
C=======================================================================
      END SUBROUTINE
C
C=======================================================================
      SUBROUTINE TWC_SY(HAMANG,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=PIFAC*FACNE2(1,NDTWBL)*FACNE2(2,NDTWBL)*FACNE2(3,NDTWBL)*
     *       FACINP*0.5D0
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBASE
C
         MX=NXVECT(IBRA,IND2HL(NDTWBL))
         MY=NYVECT(IBRA,IND2HL(NDTWBL))
         MZ=NZVECT(IBRA,IND2HL(NDTWBL))
C
         DO IKET=1,LDBASE
C
            NX=NXVECT(IKET,IND2HR(NDTWBL))
            NY=NYVECT(IKET,IND2HR(NDTWBL))
            NZ=NZVECT(IKET,IND2HR(NDTWBL))
C
            HAMANG(IBRA,IKET) = IPHAPP(MY,NY,2) * FACTOR
     *                        * CF2H00(0,MX,NX,1,NDTWBL)
     *                        * CF2H00(0,MY,NY,2,NDTWBL)
     *                        * CF2H00(0,MZ,NZ,3,NDTWBL)
C
            END DO
         END DO
C
C=======================================================================
      END SUBROUTINE
C=======================================================================
C
      SUBROUTINE TWC_SZ(HAMANG,FACINP)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
      COMMON
     *       /DIMENS/ LDBASE
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
     *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
      COMMON
     *       /COEFH2/ CF2H00(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H01(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL),
     *                CF2H02(0:ND2MAX,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                  1:NDTWBL)
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=PIFAC*FACNE2(1,NDTWBL)*FACNE2(2,NDTWBL)*FACNE2(3,NDTWBL)*
     *       FACINP*0.5D0
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBASE
C
         MX=NXVECT(IBRA,IND2HL(NDTWBL))
         MY=NYVECT(IBRA,IND2HL(NDTWBL))
         MZ=NZVECT(IBRA,IND2HL(NDTWBL))
C
         DO IKET=1,LDBASE
C
            NX=NXVECT(IKET,IND2HR(NDTWBL))
            NY=NYVECT(IKET,IND2HR(NDTWBL))
            NZ=NZVECT(IKET,IND2HR(NDTWBL))
C
            HAMANG(IBRA,IKET) = IPHAPM(MY,NY,3) * FACTOR
     *                        * CF2H00(0,MX,NX,1,NDTWBL)
     *                        * CF2H00(0,MY,NY,2,NDTWBL)
     *                        * CF2H00(0,MZ,NZ,3,NDTWBL)
C
            END DO
         END DO
C
C=======================================================================
      END SUBROUTINE
C
C=======================================================================
      SUBROUTINE TWCSIZ_NATBAS(ICHARG)
C=======================================================================
      USE hfodd_sizes
      USE OVM_PP
      USE OVM_PM
      USE TWCEPP
      USE TWCEPM
      USE NOREIG
      USE TWCWAV
C=======================================================================
      COMPLEX
     *          C_ZERO,C_UNIT
      COMPLEX
     *          CELMTS
      COMPLEX
     *          COMHAM,HAMAUX,HAMCUT
      COMPLEX
     *          EWTNAT,WAVTWC
      DIMENSION
     *          SPENAT(1:4*NDBASE)
      COMMON
     *       /DIMENS/ LDBASE
      COMMON
     *       /NATCUT/ ICUTOF
      COMMON
     *       /SPNUMS/ NUMBSP(0:NDREVE,0:NDISOS)
      COMMON
     *       /TWCDIA/ TWCSPE(1:4*NDSTAT,0:NDISOS)
C=======================================================================
      ALLOCATABLE COMHAM(:,:)
      ALLOCATABLE CELMTS(:)
C
      ALLOCATABLE HAMAUX(:,:),HAMCUT(:,:)
C
      ALLOCATABLE EWTNAT(:,:),WAVTWC(:,:)
C=======================================================================
      C_ZERO=CMPLX(0.0D0,0.0D0)
      C_UNIT=CMPLX(1.0D0,0.0D0)
C=======================================================================
C     BUILDING THE PROPER SIMPLEX STRUCTURE OF THE HAMILTONIAN
C=======================================================================
      IF (.NOT.ALLOCATED(COMHAM)) THEN
          ALLOCATE (COMHAM(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('COMHAM','TWCSIZ_NATBAS')
      END IF
C=======================================================================
C
      DO IBRA=1,4*LDBASE
         DO IKET=1,4*LDBASE
C
                IF (IBRA.LE.2*LDBASE.AND.IKET.LE.2*LDBASE) THEN
C
                    COMHAM(IBRA,IKET)=TWC_PP(IBRA,IKET,0)
C
                ELSEIF (IBRA.GT.2*LDBASE.AND.IKET.LE.2*LDBASE) THEN
C
                    COMHAM(IBRA,IKET)=TWC_PM(IBRA-2*LDBASE,IKET,1)
C
                ELSEIF (IBRA.LE.2*LDBASE.AND.IKET.GT.2*LDBASE) THEN
C
                    COMHAM(IBRA,IKET)=TWC_PM(IBRA,IKET-2*LDBASE,0)
C
                ELSEIF (IBRA.GT.2*LDBASE.AND.IKET.GT.2*LDBASE) THEN
C
                COMHAM(IBRA,IKET)=TWC_PP(IBRA-2*LDBASE,IKET-2*LDBASE,1)
C
                END IF
C
         END DO
      END DO
C
C=======================================================================
C     DIAGONALIZING H IN THE "NATURAL BASIS"
C=======================================================================
      IF (.NOT.ALLOCATED(CELMTS)) THEN
          ALLOCATE (CELMTS(1:(4*NDBASE+1)*2*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('CELMTS','TWCSIZ_NATBAS')
      END IF
C
      IF (.NOT.ALLOCATED(EWTNAT)) THEN
          ALLOCATE (EWTNAT(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('EWTNAT','TWCSIZ_NATBAS')
      END IF
C=======================================================================
      NATDIM=4*LDBASE-ICUTOF+1
C=======================================================================
      IF (.NOT.ALLOCATED(HAMAUX)) THEN
          ALLOCATE (HAMAUX(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('HAMAUX','TWCSIZ_NATBAS')
      END IF
C
      IF (.NOT.ALLOCATED(HAMCUT)) THEN
          ALLOCATE (HAMCUT(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('HAMCUT','TWCSIZ_NATBAS')
      END IF
C
      IF (.NOT.ALLOCATED(EWTAUX)) THEN
          ALLOCATE (EWTAUX(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('EWTAUX','TWCSIZ_NATBAS')
      END IF
C=======================================================================
      CALL ZGEMM('C','N',NATDIM,4*LDBASE,4*LDBASE,
     *                C_UNIT,EWTAUX,4*NDBASE,
     *                              COMHAM,4*NDBASE,
     *                       C_ZERO,HAMAUX,4*NDBASE)
C
      CALL ZGEMM('N','N',NATDIM,NATDIM,4*LDBASE,
     *                C_UNIT,HAMAUX,4*NDBASE,
     *                              EWTAUX,4*NDBASE,
     *                       C_ZERO,HAMCUT,4*NDBASE)
C
      DO IBRA=1,NATDIM
         DO IKET=1,NATDIM
C
            IF (IBRA.GE.IKET) THEN
C
                NCOUNT=IBRA+((2*NATDIM-IKET)*(IKET-1))/2
                CELMTS(NCOUNT)=HAMCUT(IBRA,IKET)
C
            END IF
C
         END DO
      END DO
C
      CALL DIAMAT(CELMTS,SPENAT,EWTNAT,NATDIM,4*NDBASE,NATDIM)
C=======================================================================
C     GOING BACK TO THE HO BASIS
C=======================================================================
      IF (.NOT.ALLOCATED(WAVTWC)) THEN
          ALLOCATE (WAVTWC(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('WAVTWC','TWCSIZ_NATBAS')
      END IF
C=======================================================================
      CALL ZGEMM('N','N',4*LDBASE,NATDIM,NATDIM,
     *                C_UNIT,EWTAUX,4*NDBASE,
     *                              EWTNAT,4*NDBASE,
     *                       C_ZERO,WAVTWC,4*NDBASE)
C=======================================================================
C     SAVING ENERGIES AND AMPLITUDES OF MANY BODY WAVE FUNCTION
C=======================================================================
      IF (.NOT.ALLOCATED(TWCOCC)) THEN
          ALLOCATE (TWCOCC(1:4*NDBASE,1:2*NDSTAT,0:NDREVE,0:NDISOS),
     *                                                 STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWCOCC','TWCSIZ')
      END IF
C=======================================================================
C
      DO IALLST=1,NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
C
         TWCSPE(IALLST,ICHARG)=SPENAT(IALLST)
C
         DO IBASE=1,2*LDBASE
            TWCOCC(IBASE,IALLST,0,ICHARG)=WAVTWC(IBASE         ,IALLST)
            TWCOCC(IBASE,IALLST,1,ICHARG)=WAVTWC(IBASE+2*LDBASE,IALLST)
         END DO
      END DO
C
C=======================================================================
      END SUBROUTINE

C=======================================================================
C
C       SUBROUTINE TWC_LX_TEST
C=======================================================================
C       USE hfodd_sizes
C=======================================================================
C       COMMON
C      *       /DIMENS/ LDBASE
C       COMMON
C      *       /WEIGH2/ WGHER2(1:NDGAUS,1:NDKART),
C      *                WGHER4(1:NDGAUS,1:NDKART,1:NDTWHE)
C
C       COMMON
C      *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
C      *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
C      *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C       COMMON
C      *       /TWCHER/ HERLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DHELE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DDHLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                HERRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DHERI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DDHRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                HERLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE)
C
C       COMMON
C      *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
C      *                NYVECT(1:NDBASE,1:NDTWCE),
C      *                NZVECT(1:NDBASE,1:NDTWCE)
C       COMMON
C      *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL)
C       COMMON
C      *       /GAUMAX/ NXHERM,NYHERM,NZHERM
C
C       COMMON
C      *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
C       COMMON
C      *       /XYZINT/ VARINT(1:NDGAUS,1:NDKART,1:NDTWBL)
C       DIMENSION
C      *          HAMANG(1:NDBASE,1:NDBASE)
C=======================================================================
C
C
C       HAMANG(:,:)=0.0D0
C       do ltwcen=1,ndtwbl
C
C       DO IBRA=1,LDBASE
C
C          MX=NXVECT(IBRA,IND2HL(LTWCEN))
C          MY=NYVECT(IBRA,IND2HL(LTWCEN))
C          MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
C          DO IKET=1,LDBASE
C
C             NX=NXVECT(IKET,IND2HR(LTWCEN))
C             NY=NYVECT(IKET,IND2HR(LTWCEN))
C             NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
C             RESULT=0.D0
C
C             DO IX=1,NXHERM
C                DO IY=1,NYHERM
C                   DO IZ=1,NZHERM
C
C                      RESULT =   RESULT +
C      *                          WGHER2(IX,1)*WGHER2(IY,2)*WGHER2(IZ,3)
C      *              *FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C      *              *     HERLE2(MX,IX,1,LTWCEN)*HERRI2(NX,IX,1,LTWCEN)
C      *              *    (HERLE2(MY,IY,2,LTWCEN)*HERRI2(NY,IY,2,LTWCEN)
C      *              *     HERLE2(MZ,IZ,3,LTWCEN)*DHERI2(NZ,IZ,3,LTWCEN)
C      *              *                               VARINT(IY,2,LTWCEN)
C      *              -     HERLE2(MZ,IZ,3,LTWCEN)*HERRI2(NZ,IZ,3,LTWCEN)
C      *              *     HERLE2(MY,IY,2,LTWCEN)*DHERI2(NY,IY,2,LTWCEN)
C      *              *                              VARINT(IZ,3,LTWCEN))
C
C                   END DO
C                END DO
C             END DO
C
C             PRINT*, LTWCEN,IBRA,IKET, RESULT
C          END DO
C       END DO
C       end do
C
C       END SUBROUTINE
C=======================================================================
C
C       SUBROUTINE TWC_LY_TEST
C=======================================================================
C       USE hfodd_sizes
C=======================================================================
C       COMMON
C      *       /DIMENS/ LDBASE
C       COMMON
C      *       /WEIGH2/ WGHER2(1:NDGAUS,1:NDKART),
C      *                WGHER4(1:NDGAUS,1:NDKART,1:NDTWHE)
C
C       COMMON
C      *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
C      *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
C      *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C       COMMON
C      *       /TWCHER/ HERLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DHELE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DDHLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                HERRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DHERI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DDHRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                HERLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE)
C
C       COMMON
C      *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
C      *                NYVECT(1:NDBASE,1:NDTWCE),
C      *                NZVECT(1:NDBASE,1:NDTWCE)
C       COMMON
C      *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL)
C       COMMON
C      *       /GAUMAX/ NXHERM,NYHERM,NZHERM
C
C       COMMON
C      *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
C       COMMON
C      *       /XYZINT/ VARINT(1:NDGAUS,1:NDKART,1:NDTWBL)
C       DIMENSION
C      *          HAMANG(1:NDBASE,1:NDBASE),HAUXPM(1:NDBASE,1:NDBASE)
C=======================================================================
C
C
C       HAMANG(:,:)=0.0D0
C       DO LTWCEN=1,NDTWBL
C
C       DO IBRA=1,LDBASE
C
C          MX=NXVECT(IBRA,IND2HL(LTWCEN))
C          MY=NYVECT(IBRA,IND2HL(LTWCEN))
C          MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
C          DO IKET=1,LDBASE
C
C             NX=NXVECT(IKET,IND2HR(LTWCEN))
C             NY=NYVECT(IKET,IND2HR(LTWCEN))
C             NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
C             RESULT=0.D0
C
C             DO IX=1,NXHERM
C                DO IY=1,NYHERM
C                   DO IZ=1,NZHERM
C
C                      RESULT =                                   RESULT
C      *              +            WGHER2(IX,1)*WGHER2(IY,2)*WGHER2(IZ,3)
C      *              *FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C      *              *     HERLE2(MY,IY,2,LTWCEN)*HERRI2(NY,IY,2,LTWCEN)
C      *              *    (HERLE2(MZ,IZ,3,LTWCEN)*HERRI2(NZ,IZ,3,LTWCEN)
C      *              *     HERLE2(MX,IX,1,LTWCEN)*DHERI2(NX,IX,1,LTWCEN)
C      *              *                               VARINT(IZ,3,LTWCEN)
C      *              -     HERLE2(MX,IX,1,LTWCEN)*HERRI2(NX,IX,1,LTWCEN)
C      *              *     HERLE2(MZ,IZ,3,LTWCEN)*DHERI2(NZ,IZ,3,LTWCEN)
C      *              *                              VARINT(IX,1,LTWCEN))
C
C                   END DO
C                END DO
C             END DO
C
C           PRINT*, LTWCEN,IBRA,IKET,RESULT
C          END DO
C       END DO
C       END DO
C
C       END SUBROUTINE
C=======================================================================
C
C       SUBROUTINE TWC_LZ_TEST
C=======================================================================
C       USE hfodd_sizes
C=======================================================================
C       COMMON
C      *       /DIMENS/ LDBASE
C       COMMON
C      *       /WEIGH2/ WGHER2(1:NDGAUS,1:NDKART),
C      *                WGHER4(1:NDGAUS,1:NDKART,1:NDTWHE)
C
C       COMMON
C      *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
C      *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
C      *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C       COMMON
C      *       /TWCHER/ HERLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DHELE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DDHLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                HERRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DHERI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                DDHRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C      *                HERLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C      *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE)
C
C       COMMON
C      *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
C      *                NYVECT(1:NDBASE,1:NDTWCE),
C      *                NZVECT(1:NDBASE,1:NDTWCE)
C       COMMON
C      *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL)
C       COMMON
C      *       /GAUMAX/ NXHERM,NYHERM,NZHERM
C
C       COMMON
C      *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
C      *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
C       COMMON
C      *       /XYZINT/ VARINT(1:NDGAUS,1:NDKART,1:NDTWBL)
C       DIMENSION
C      *          HAMANG(1:NDBASE,1:NDBASE),HAUXPM(1:NDBASE,1:NDBASE)
C=======================================================================
C
C
C       HAMANG(:,:)=0.0D0
C
C
C       DO LTWCEN=1,NDTWBL
C       DO IBRA=1,LDBASE
C
C          MX=NXVECT(IBRA,IND2HL(LTWCEN))
C          MY=NYVECT(IBRA,IND2HL(LTWCEN))
C          MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
C          DO IKET=1,LDBASE
C
C             NX=NXVECT(IKET,IND2HR(LTWCEN))
C             NY=NYVECT(IKET,IND2HR(LTWCEN))
C             NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
C             RESULT=0.D0
C
C             DO IX=1,NXHERM
C                DO IY=1,NYHERM
C                   DO IZ=1,NZHERM
C
C                      RESULT =                                   RESULT
C      *              +            WGHER2(IX,1)*WGHER2(IY,2)*WGHER2(IZ,3)
C      *              *FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C      *              *     HERLE2(MZ,IZ,3,LTWCEN)*HERRI2(NZ,IZ,3,LTWCEN)
C      *              *    (HERLE2(MX,IX,1,LTWCEN)*HERRI2(NX,IX,1,LTWCEN)
C      *              *     HERLE2(MY,IY,2,LTWCEN)*DHERI2(NY,IY,2,LTWCEN)
C      *              *                               VARINT(IX,1,LTWCEN)
C      *              -     HERLE2(MY,IY,2,LTWCEN)*HERRI2(NY,IY,2,LTWCEN)
C      *              *     HERLE2(MX,IX,1,LTWCEN)*DHERI2(NX,IX,1,LTWCEN)
C      *              *                              VARINT(IY,2,LTWCEN))
C
C                  END DO
C               END DO
C            END DO
C
C            PRINT*, LTWCEN,IBRA,IKET, RESULT
C         END DO
C      END DO
C      END DO
C
C      END SUBROUTINE
C=======================================================================
C
C=======================================================================
C
C      SUBROUTINE TWC_INTKIN_TEST(HAMKIN,HBMASS,LTWCEN)
C=======================================================================
C      USE hfodd_sizes
C=======================================================================
C      COMMON
C     *       /DIMENS/ LDBASE
C      COMMON
C     *       /WEIGH2/ WGHER2(1:NDGAUS,1:NDKART),
C     *                WGHER4(1:NDGAUS,1:NDKART,1:NDTWHE)
C
C      COMMON
C     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
C     *                IND4HK(1:NDTWHE),IND4HL(1:NDTWHE),
C     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C      COMMON
C     *       /TWCHER/ HERLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C     *                DHELE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C     *                DDHLE2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C     *                HERRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C     *                DHERI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C     *                DDHRI2(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWBL),
C     *                HERLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C     *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C     *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C     *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C     *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE),
C     *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWHE)
C
C      COMMON
C     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
C     *                NYVECT(1:NDBASE,1:NDTWCE),
C     *                NZVECT(1:NDBASE,1:NDTWCE)
C      COMMON
C     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL)
C      COMMON
C     *       /GAUMAX/ NXHERM,NYHERM,NZHERM
C
C      COMMON
C     *       /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),
C     *                IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),
C     *                IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),
C     *                IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
C      DIMENSION
C     *          HAMKIN(1:NDBASE,1:NDBASE)
C=======================================================================
C
C
C      DO IBRA=1,LDBASE
CC
C         MX=NXVECT(IBRA,IND2HL(LTWCEN))
C         MY=NYVECT(IBRA,IND2HL(LTWCEN))
C         MZ=NZVECT(IBRA,IND2HL(LTWCEN))
CC
C         DO IKET=1,LDBASE
CC
C            NX=NXVECT(IKET,IND2HR(LTWCEN))
C            NY=NYVECT(IKET,IND2HR(LTWCEN))
C            NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
C            RESULX=0.D0
C            RESULY=0.D0
C            RESULZ=0.D0
CC
C            RESKIX=0.D0
C            RESKIY=0.D0
C            RESKIZ=0.D0
C
CC
CC=======================================================================
CC     X DIRECTION
CC=======================================================================
CC
C            DO IX=1,NXHERM
CC
C               RESULX=RESULX+WGHER2(IX,1)*FACNE2(1,LTWCEN)
C     *                             *HERLE2(MX,IX,1,LTWCEN)
C     *                             *HERRI2(NX,IX,1,LTWCEN)
C               RESKIX=RESKIX+WGHER2(IX,1)*FACNE2(1,LTWCEN)
C     *                             *HERLE2(MX,IX,1,LTWCEN)
C     *                             *DDHRI2(NX,IX,1,LTWCEN)
CC
C            END DO
CC
C            DO IY=1,NYHERM
CC
C               RESULY=RESULY+WGHER2(IY,2)*FACNE2(2,LTWCEN)
C     *                             *HERLE2(MY,IY,2,LTWCEN)
C     *                             *HERRI2(NY,IY,2,LTWCEN)
C               RESKIY=RESKIY+WGHER2(IY,2)*FACNE2(2,LTWCEN)
C     *                             *HERLE2(MY,IY,2,LTWCEN)
C     *                             *DDHRI2(NY,IY,2,LTWCEN)
CC
C            END DO
CC
C            DO IZ=1,NZHERM
CC
C               RESULZ=RESULZ+WGHER2(IZ,3)*FACNE2(3,LTWCEN)
C     *                             *HERLE2(MZ,IZ,3,LTWCEN)
C     *                             *HERRI2(NZ,IZ,3,LTWCEN)
C               RESKIZ=RESKIZ+WGHER2(IZ,3)*FACNE2(3,LTWCEN)
C     *                             *HERLE2(MZ,IZ,3,LTWCEN)
C     *                             *DDHRI2(NZ,IZ,3,LTWCEN)
CC
C            END DO
CC
C            HAMKIN(IBRA,IKET)= -HBMASS*IPHAPP(MY,NY,0)
C     *                         *(RESKIX*RESULY*RESULZ
C     *                          +RESULX*RESKIY*RESULZ
C     *                          +RESULX*RESULY*RESKIZ)
CC
C         END DO
C      END DO
CC
CC     SIMPLEX BLOCKS???
CC
C      END SUBROUTINE
C=======================================================================
C
C
C=======================================================================
C
C      SUBROUTINE TWCSIZ(ICHARG)
C=======================================================================
C      USE TWCWAV
C      USE TWCEPP
C      USE TWCEPM
C      USE TWCEIZ
C      USE OVM_PP
C      USE OVM_PM
C      USE OVMSIM
C=======================================================================
C      USE hfodd_sizes
C=======================================================================
C      COMPLEX
C     *          CELMTS
C      COMPLEX
C     *          AUXDIC
C
C      ALLOCATABLE AUXDIC(:)
C      ALLOCATABLE AUXDID(:)
C=======================================================================
C      ALLOCATABLE CELMTS(:,:)
C=======================================================================
C      COMMON
C     *       /TWCDIA/ TWCSPE(1:4*NDSTAT,0:NDISOS),
C     *                SPTWCZ(1:4*NDBASE)
C      COMMON
C     *       /EIGLIM/ EIGMIN,EIGMAX
C      COMMON
C     *       /SPNUMS/ NUMBSP(0:NDREVE,0:NDISOS)
C      COMMON
C     *       /DIMENS/ LDBASE
C      COMMON
C     *       /CFIPRI/ NFIPRI
C
C=======================================================================
C
C      CALL CPUTIM('TWCSIZ',1)
C
C=======================================================================
C
C      IALLOC=0
C
C=======================================================================
C      IF (.NOT.ALLOCATED(TWCOCC)) THEN
C          ALLOCATE (TWCOCC(1:4*NDBASE,1:2*NDSTAT,0:NDREVE,0:NDISOS),
C     *                                                 STAT=IALLOC)
C          IF (IALLOC.NE.0) CALL NOALLO('TWCOCC','TWCSIZ')
C      END IF
C=======================================================================
C      IF (.NOT.ALLOCATED(CELMTS)) THEN
C          ALLOCATE (CELMTS(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
C          IF (IALLOC.NE.0) CALL NOALLO('CELMTS','TWCSIZ')
C      END IF
C=======================================================================
C      IF (.NOT.ALLOCATED(OVERME)) THEN
C         ALLOCATE (OVERME(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
C          IF (IALLOC.NE.0) CALL NOALLO('OVERME','TWCSIZ')
C      END IF
C=======================================================================
C
C         CELMTS(:,:)=0.D0
C         OVERME(:,:)=0.D0
C
C=======================================================================
C        THIS SUBROUTINE PERFORMS THE HF DIAGONALIZATION IN THE CASE  OF
C        NO CONSERVED SYMMETRIES AND TWO-CENTER BASIS
C=======================================================================
C        CONSTRUCTING THE HF HAMILTONIAN MATRIX, WHICH HAS THE FOLLOWING
C        CENTER 1 (1),  CENTER 2 (2),  SIMPLEX=+i (U),AND SIMPLEX=-i (D)
C        STRUCTURE OF THE 4x4 LDBASE BLOCKS:
C=======================================================================
C
C         -------------------------------------------------------------
C         |         SIMPLEX-UP          |         SIMPLEX-DOWN        |
C         |-----------------------------|-----------------------------|
C         |   CENTER 1   |   CENTER 2   |   CENTER 1   |   CENTER 2   |
C         -------------------------------------------------------------
C
C         -------------------------------------------------------------
C         |    1U-1U     |    1U-2U     |    1U-1D     |    1U-2D     |
C         |--------------|--------------|--------------|--------------|
C         |    2U-1U     |    2U-2U     |    2U-1D     |    2U-2D     |
C         |--------------|--------------|--------------|--------------|
C         |    1D-1U     |    1D-2U     |    1D-1D     |    1D-2D     |
C         |--------------|--------------|--------------|--------------|
C         |    2D-1U     |    2D-2U     |    2D-1D     |    2D-2D     |
C         -------------------------------------------------------------
C
C=======================================================================
C
C      DO IBRA=1,4*LDBASE
C         DO IKET=1,4*LDBASE
C
C                IF (IBRA.LE.2*LDBASE.AND.IKET.LE.2*LDBASE) THEN
C
C                    CELMTS(IBRA,IKET)=TWC_PP(IBRA,IKET,0)
C                    OVERME(IBRA,IKET)=OVMAPP(IBRA,IKET,0)
C
C                ELSEIF (IBRA.GT.2*LDBASE.AND.IKET.LE.2*LDBASE) THEN
C
C                    CELMTS(IBRA,IKET)=TWC_PM(IBRA-2*LDBASE,IKET,1)
C                    OVERME(IBRA,IKET)=OVMAPM(IBRA-2*LDBASE,IKET,1)
C
C                ELSEIF (IBRA.LE.2*LDBASE.AND.IKET.GT.2*LDBASE) THEN
C
C                    CELMTS(IBRA,IKET)=TWC_PM(IBRA,IKET-2*LDBASE,0)
C                    OVERME(IBRA,IKET)=OVMAPM(IBRA,IKET-2*LDBASE,0)
C
C                ELSEIF (IBRA.GT.2*LDBASE.AND.IKET.GT.2*LDBASE) THEN
C
C                CELMTS(IBRA,IKET)=TWC_PP(IBRA-2*LDBASE,IKET-2*LDBASE,1)
C               OVERME(IBRA,IKET)=OVMAPP(IBRA-2*LDBASE,IKET-2*LDBASE,1)
C
C                END IF
C
C         END DO
C      END DO
C
C=======================================================================
C
C      IALLOC=0
C      INFORM=0
C
C      LDACTU=4*LDBASE
C
C      LDWORK=2*LDACTU-1
C      IF (.NOT.ALLOCATED(AUXDIC)) THEN
C      ALLOCATE(AUXDIC(1:LDWORK),STAT=IALLOC)
C      IF (IALLOC.NE.0) CALL NOALLO('AUXDIC','TWCSIZ')
C      END IF
C
C      LRWORK=3*LDACTU-2
C      IF (.NOT.ALLOCATED(AUXDID)) THEN
C      ALLOCATE(AUXDID(1:LRWORK),STAT=IALLOC)
C      IF (IALLOC.NE.0) CALL NOALLO('AUXDID','TWCSIZ')
C      END IF
C
C      CALL ZHEGV(1,'V','L',4*LDBASE,CELMTS,4*NDBASE,OVERME,4*NDBASE,
C     *           SPTWCZ,AUXDIC,LDWORK,AUXDID,INFORM)
C
C      IF (INFORM.NE.0) THEN
C
C          WRITE(NFIPRI,'(//,1X,20(1H=),
C     *                   ''  ERROR IN ZHEGV, INFORM='',
C     *                   I4,1X,20(1H=),/,
C     *                      1X,79(1H=))')
C     *
C     *          INFORM
C          STOP ' ERROR IN ZHEGV'
C
C       END IF
C
C      DEALLOCATE(AUXDIC,AUXDID)
C
C=======================================================================
C
C      DO IALLST=1,NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
C         TWCSPE(IALLST,ICHARG)=SPTWCZ(IALLST)
C      END DO
C
C=======================================================================
C
C      DO IALLST=1,NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
C         DO IBASE=1,2*LDBASE
C            TWCOCC(IBASE,IALLST,0,ICHARG)=CELMTS(IBASE         ,IALLST)
C            TWCOCC(IBASE,IALLST,1,ICHARG)=CELMTS(IBASE+2*LDBASE,IALLST)
C         END DO
C      END DO
C
C
C=======================================================================
C      DEALLOCATE (CELMTS)
C=======================================================================
C
C      CALL CPUTIM('TWCSIZ',0)
C
C=======================================================================
C
C      RETURN
C      END
C=======================================================================
C
      END MODULE hfodd_twocen
