      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) = 55
      MODSET(18) = 18
      !
      END SUBROUTINE Set_version_hfodd_twocen
      !
      module hfodd_twocen
C=======================================================================
C
      USE hfodd_sizes
C
C=======================================================================
C
      CONTAINS
C
C=======================================================================
C
      SUBROUTINE TWC_INDGEN(ITWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *      /HERMIN/  INTWHE(1:NDTWCE,1:NDTWCE,1:NDTWCE,1:NDTWCE)
      COMMON
     *      /INFAST/  IDCOPY(1:NDTWDD,0:1),IDCOCO(1:NDTWDD,0:1)
C=======================================================================
C     THIS SUBROUTINES SETS THE INDICES NEEDED FOR THE TWO-CENTER BASIS
C     LOOPS. TWO DIFFERENT TYPES: WHEN TWO WAVE FUNCTIONS ARE INVOLVED
C     OR WHEN THERE ARE FOUR OF THEM.
C=======================================================================
C     TWO WAVE FUNCTIONS INVOLVED. USED TO:
C     - DETERMINE THE QUANTUM NUMBERS OF THE DIFFERENT CENTRE BLOCKS
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
C=======================================================================
C     FOUR WAVE FUNCTIONS INVOLVED. USED TO:
C     - STORE DIFFERENT DENSITIES IN FUNCTION OF THE DIFFERENT
C       GAUSSIAN FACTORS.
C     - OBTAIN MEAN FIELD MATRIX ELEMENTS.
C     - SUM DIFFERENT TERMS IN THE ENERGY DENSITY INTEGRAL.
C
C     ORDER OF THE INDICES WHEN USED IN INTEGH:  < I | O_KL | J >
C=======================================================================
C
      IF (ITWCEN.EQ.2) THEN
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     DEFINING INDICES FOR DENSITY-DEPENDENT TERMS
C=======================================================================
C
          IND4HK(NDTWHE+1*ISWTWC)=2
          IND4HK(NDTWHE+2*ISWTWC)=2
          IND4HK(NDTWHE+3*ISWTWC)=2
          IND4HK(NDTWHE+4*ISWTWC)=1
          IND4HK(NDTWHE+5*ISWTWC)=1
          IND4HK(NDTWHE+6*ISWTWC)=2
! FOR PAIRING DENSITY-DEPENDENT TERMS
          IND4HK(NDTWHE+7*ISWTWC)=1
          IND4HK(NDTWHE+8*ISWTWC)=1
C
          IND4HL(NDTWHE+1*ISWTWC)=2
          IND4HL(NDTWHE+2*ISWTWC)=1
          IND4HL(NDTWHE+3*ISWTWC)=2
          IND4HL(NDTWHE+4*ISWTWC)=1
          IND4HL(NDTWHE+5*ISWTWC)=1
          IND4HL(NDTWHE+6*ISWTWC)=1
! FOR PAIRING DENSITY-DEPENDENT TERMS
          IND4HL(NDTWHE+7*ISWTWC)=2
          IND4HL(NDTWHE+8*ISWTWC)=2
C
      ELSEIF (ITWCEN.EQ.1) THEN
C
           IND4HI(:)=1
           IND4HJ(:)=1
           IND4HK(:)=1
           IND4HL(:)=1
C
           INTWHE(:,:,:,:)=1
C
      END IF
C
C=======================================================================
C     HERE WE DEFINE INDICES FOR COPYING DENSITIES TO AVOID CALLING
C     DENSHF SIXTEEN TIMES. TWO INDICES:
C        IDCOPY: LABEL OF THE ORIGINAL DENSITY. IF 0 DENSHF IS CALLED
C        IDCOCO: DETERMINE WHETHER TO COPY THE ORIGINAL OR ITS COMPLEX
C                CONJUGATE
C=======================================================================
C
      IDCOPY(:,:)=0
      IDCOCO(:,:)=0
C
      IF (ITWCEN.EQ.2) THEN
C
          IDCOPY(ISWTWC*3 +(1-ISWTWC),0)= 2
          IDCOPY(ISWTWC*7 +(1-ISWTWC),0)= 6
          IDCOPY(ISWTWC*9 +(1-ISWTWC),0)= 5
          IDCOPY(ISWTWC*10+(1-ISWTWC),0)= 6
          IDCOPY(ISWTWC*11+(1-ISWTWC),0)= 7
          IDCOPY(ISWTWC*12+(1-ISWTWC),0)= 8
          IDCOPY(ISWTWC*15+(1-ISWTWC),0)=14
C
          IDCOCO(ISWTWC*3 +(1-ISWTWC),0)= 1
          IDCOCO(ISWTWC*7 +(1-ISWTWC),0)= 1
          IDCOCO(ISWTWC*15+(1-ISWTWC),0)= 1
C
          IDCOPY(ISWTWC*9 +(1-ISWTWC),1)= 5
          IDCOPY(ISWTWC*10+(1-ISWTWC),1)= 6
          IDCOPY(ISWTWC*11+(1-ISWTWC),1)= 7
          IDCOPY(ISWTWC*12+(1-ISWTWC),1)= 8
      END IF
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:NDTWDD),IND4HL(1:NDTWDD),
     *                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:NDTWDD),
     *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
      COMMON
     *       /TWCEXP/ HEREXP(0:ND2MAX,1:NDGAUS,1:NDKART)
      COMMON
     *       /TWCPMU/ PNTTWC(1:NDGAUS,1:NDKART,1:NDTWBL)
      COMMON
     *       /TWCFAC/ QUADSUM(1:NDKART,1:NDTWBL),
     *                QUADCEN(1:NDKART,1:NDTWBL)
      COMMON
     *       /TWPOIN/ XZERLE(1:NDGAUS,1:NDKART,1:NDTWDD),
     *                XZERRI(1:NDGAUS,1:NDKART,1:NDTWDD),
     *                PNTWC4(1:NDGAUS,1:NDKART,1:NDTWHE)
C
C=======================================================================
      IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_DEFINT'
C=======================================================================
C
      NMAIN2=2*NOSACT+2
C
      IF (NGAUSS.GT.NDGAUS) THEN
C
         WRITE(NFIPRI,'(//,1X,''NGAUSS='',I2,'' .GT. NDGAUS='',I2,
     *                   ''   STOP IN TWC_DEFINT'',//)') NGAUSS,NDGAUS
         STOP ' NGAUSS.GT.NDGAUS  IN  TWC_DEFINT'
C
      END IF
C
C=======================================================================
C     COMPUTING HERMITE POLYNOMIALS EVALUATED IN THE COMMON LATTICE
C     OF INTEGRATION. TO BE USED WITH THE COEFFICIENTS CF4H FROM DEVHER
C=======================================================================
C
      CALL HERMIT(NGAUSS,XHERMI,WHERMI,NDGAUS)
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     COMPUTING HERMITE POLYNOMIALS EVALUATED IN THE COMMON LATTICE
C     OF INTEGRATION AFTER THE CHANGE OF VARIABLE WHEN ONLY TWO WAVE
C     FUNCTIONS ARE ACTING. TO BE USED FOR COMPUTING MATRIX ELEMENTS OF
C     ONE-BODY OPERATORS (NORM OVERLAP, L, S, ETC)
C=======================================================================
C
      DO INDBLO=1,NDTWBL
C
          HOSCAI=HOMSCA(KARTEZ,IND2HL(INDBLO))
          HOSCAJ=HOMSCA(KARTEZ,IND2HR(INDBLO))
C
          HOSHII=CENCOR(KARTEZ,IND2HL(INDBLO))
          HOSHIJ=CENCOR(KARTEZ,IND2HR(INDBLO))
C
          SUMOSL=HOSCAI**2+HOSCAJ**2
          SUMCEN=HOSCAI**2*HOSHII+HOSCAJ**2*HOSHIJ
C
          FACTO1=SQRT(2.D0/SUMOSL)
          FACTO2=SUMCEN/SUMOSL
C
C=======================================================================
C     LEFT WAVE FUNCTION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = HOSCAI*(XHERMI(IZEROS)*FACTO1+FACTO2-HOSHII)
              PNTTWC(IZEROS,KARTEZ,INDBLO)=XHERMI(IZEROS)*FACTO1+FACTO2
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)*HOSCAI
     *                                   / HERFAC( NOSCIL )
C
                  DDHLE2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAI)
     *                                   *HOSCAI*HOSCAI
     *                                   / HERFAC( NOSCIL )
C
              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)*HOSCAJ
     *                                   / HERFAC( NOSCIL )
C
                  DDHRI2(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAJ)
     *                                   * HOSCAJ*HOSCAJ
     *                                   / HERFAC( NOSCIL )
C
              END DO
C
          END DO
C
C=======================================================================
C         SAVING COEFFICIENTS TO BE USED LATER
C=======================================================================
C
          FACNE2(KARTEZ,INDBLO)= SQRT(2.D0/SUMOSL)
     *     * EXP( -0.5D0*HOSCAI**2*HOSCAJ**2*(HOSHII-HOSHIJ)**2/SUMOSL )
C
          QUADSUM(KARTEZ,INDBLO)=SQRT(2.D0/SUMOSL)
C
          QUADCEN(KARTEZ,INDBLO)=SUMCEN/SUMOSL
C
      END DO
C
C=======================================================================
C     COMPUTING HERMITE POLYNOMIALS EVALUATED IN THE COMMON LATTICE
C     OF INTEGRATION AFTER THE CHANGE OF VARIABLE WHEN FOUR WAVE
C     FUNCTIONS ARE ACTING. TO BE USED FOR OBTAINING THE MEAN FIELDS AND
C     THE ENERGY OF THE FUNCTIONAL
C=======================================================================
C
      DO INDBLO=1,NDTWHE
C
          HOSCAI=HOMSCA(KARTEZ,IND4HI(INDBLO))
          HOSCAJ=HOMSCA(KARTEZ,IND4HJ(INDBLO))
          HOSCAK=HOMSCA(KARTEZ,IND4HK(INDBLO))
          HOSCAL=HOMSCA(KARTEZ,IND4HL(INDBLO))
C
          HOSHII=CENCOR(KARTEZ,IND4HI(INDBLO))
          HOSHIJ=CENCOR(KARTEZ,IND4HJ(INDBLO))
          HOSHIK=CENCOR(KARTEZ,IND4HK(INDBLO))
          HOSHIL=CENCOR(KARTEZ,IND4HL(INDBLO))
C
          SUMOSL=HOSCAI**2+HOSCAJ**2+HOSCAK**2+HOSCAL**2
          SUMCEN=HOSCAI**2*HOSHII+HOSCAJ**2*HOSHIJ
     *          +HOSCAK**2*HOSHIK+HOSCAL**2*HOSHIL
C
          FACTO1=SQRT(2.D0/SUMOSL)
          FACTO2=SUMCEN/SUMOSL
C
C=======================================================================
C     LEFT WAVE FUNCTION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS
C
              XZER = HOSCAK*(XHERMI(IZEROS)*FACTO1+FACTO2-HOSHIK)
              PNTWC4(IZEROS,KARTEZ,INDBLO)=XHERMI(IZEROS)*FACTO1+FACTO2
              XZERLE(IZEROS,KARTEZ,INDBLO) = XZER
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(HOSCAK)
     *                                   / HERFAC( NOSCIL )
C
                  DHELE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAK)*HOSCAK
     *                                   / HERFAC( NOSCIL )
C
                  DDHLE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAK)
     *                                   *HOSCAK*HOSCAK
     *                                   / HERFAC( NOSCIL )
C
              END DO
C
          END DO
C
C=======================================================================
C     RIGHT WAVE FUNCTION
C=======================================================================
C
          DO IZEROS = 1,NGAUSS

              XZER = HOSCAL*(XHERMI(IZEROS)*FACTO1+FACTO2-HOSHIL)
              XZERRI(IZEROS,KARTEZ,INDBLO) = XZER
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(HOSCAL)
     *                                   / HERFAC( NOSCIL )
C
                  DHERI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAL) *HOSCAL
     *                                   / HERFAC( NOSCIL )
C
                  DDHRI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAL)
     *                                   * HOSCAL*HOSCAL
     *                                   / HERFAC( NOSCIL )
C
              END DO
C
          END DO
C
C=======================================================================
C         SAVING COEFFICIENTS TO BE USED LATER
C=======================================================================
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
      END DO
C
C=======================================================================
C        TERMS NEEDED FOR DENSITY-DEPENDENT INTERACTION
C=======================================================================
C
      IF (ITWCEN.EQ.2) THEN !NOT INDENTED
C
      DO INDBLO=1,NDTWHE
C
C=======================================================================
C        CONDITIONALS TO COMPUTE SPECIFIC POLYNOMIALS
C=======================================================================
C
         IF (INDBLO.NE.1 .AND.INDBLO.NE. 2.AND.
     *       INDBLO.NE.8 .AND.INDBLO.NE.16     ) CYCLE
C
                                          IAGAIN=1
         IF (INDBLO.EQ.1.OR.INDBLO.EQ.16) IAGAIN=3
C
         DO IAGAI2=1,IAGAIN
C
            IF (INDBLO.EQ.1) THEN
                IF (IAGAI2.EQ.1 ) THEN
                    INDE_C= 2
                    INDE_D= 2
                    INDNEW=17
                ELSEIF (IAGAI2.EQ.2) THEN
                    INDE_C= 2
                    INDE_D= 1
                    INDNEW=18
                ELSEIF (IAGAI2.EQ.3) THEN
                    INDE_C= 1
                    INDE_D= 2
                    INDNEW=23
                END IF
            ELSEIF (INDBLO.EQ.2 ) THEN
                INDE_C= 2
                INDE_D= 2
                INDNEW=19
            ELSEIF (INDBLO.EQ.8 ) THEN
                INDE_C= 1
                INDE_D= 1
                INDNEW=20
            ELSEIF (INDBLO.EQ.16) THEN
                IF (IAGAI2.EQ.1) THEN
                    INDE_C= 1
                    INDE_D= 1
                    INDNEW=21
                ELSEIF(IAGAI2.EQ.2) THEN
                    INDE_C= 2
                    INDE_D= 1
                    INDNEW=22
                ELSEIF(IAGAI2.EQ.3) THEN
                    INDE_C= 1
                    INDE_D= 2
                    INDNEW=24
                END IF
            END IF
C
            HOSCAC=HOMSCA(KARTEZ,INDE_C)
            HOSHIC=CENCOR(KARTEZ,INDE_C)
            HOSCAD=HOMSCA(KARTEZ,INDE_D)
            HOSHID=CENCOR(KARTEZ,INDE_D)
C
            DO IZEROS = 1,NGAUSS
C
               XZER=HOSCAC*(PNTWC4(IZEROS,KARTEZ,INDBLO)-HOSHIC)
               XZERLE(IZEROS,KARTEZ,INDNEW) = XZER
C
               CALL D_HERM(XZER,NMAIN2,PHERMI,DHERMI,ND2MAX+1)
C
               DO NOSCIL = 0,NMAIN2
C
                  HERLE4(NOSCIL,IZEROS,KARTEZ,INDNEW) =
     *                                     PHERMI(NOSCIL+1)
     *                                   *     SQRT(HOSCAC)
     *                                   / HERFAC( NOSCIL )
C
                  DHELE4(NOSCIL,IZEROS,KARTEZ,INDNEW) =
     *                                         (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAC)*HOSCAC
     *                                   / HERFAC( NOSCIL )
C
                  DDHLE4(NOSCIL,IZEROS,KARTEZ,INDNEW) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAC)
     *                                   *HOSCAC*HOSCAC
     *                                   / HERFAC( NOSCIL )
C
               END DO
C
            END DO
C
            DO IZEROS = 1,NGAUSS
C
               XZER=HOSCAD*(PNTWC4(IZEROS,KARTEZ,INDBLO)-HOSHID)
               XZERRI(IZEROS,KARTEZ,INDNEW) = XZER
C
               CALL D_HERM(XZER,NMAIN2,PHERMI,DHERMI,ND2MAX+1)
C
               DO NOSCIL = 0,NMAIN2
C
                  HERRI4(NOSCIL,IZEROS,KARTEZ,INDNEW) =
     *                                     PHERMI(NOSCIL+1)
     *                                   *     SQRT(HOSCAD)
     *                                   / HERFAC( NOSCIL )
C
                  DHERI4(NOSCIL,IZEROS,KARTEZ,INDNEW) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                                   * SQRT(HOSCAD)*HOSCAD
     *                                   / HERFAC( NOSCIL )
C
                  DDHRI4(NOSCIL,IZEROS,KARTEZ,INDNEW) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                                   * SQRT(HOSCAD)
     *                                   *HOSCAD*HOSCAD
     *                                   / HERFAC( NOSCIL )
C
               END DO
C
            END DO
C
         END DO !IGAIN2
C
      END DO !INDBLO
C
      END IF !ITWCEN NOT INDENTED
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_DEVMUL(NOSACT,NGAUSS,KARTEZ,ITWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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:NDTWDD),
     *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
      COMMON
     *       /TWCINT/ WGHERM(1:NDGAUS,1:NDKART),
     *                PTHERM(1:NDGAUS,1:NDKART)
      COMMON
     *       /POWNMU/ NMUPOW
      COMMON
     *       /TWCPMU/ PNTTWC(1:NDGAUS,1:NDKART,1:NDTWBL)
      COMMON
     *     /COMU2C/ CMUTWC(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART,1:NDTWBL)
      COMMON
     *     /COMUSH/ CMUSH1(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                1:NDTWBL),
     *              CMUSH2(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                1:NDTWBL)
C
C=======================================================================
C     THIS SUBROUTINE IS THE TWO-CENTER EQUIVALENT OF DEVMUL.
C=======================================================================
C
      IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_DEVMUL'
C
      IF (NMUPOW.GT.NDMULT) STOP ' INCREASE NDMULT IN TWC_DEVMUL'
C=======================================================================
C
      DO LTWCEN=1,NDTWBL
C
         DO N=0,NOSACT
            DO M=0,NOSACT
               DO K=0,NMUPOW
C
                  RESULT=0.0D0
C
                  DO IGAUSS=1,NGAUSS
C
                     RESULT=RESULT+  WGHERM(IGAUSS,KARTEZ)
     *                     *HERLE2(N,IGAUSS,KARTEZ,LTWCEN)
     *                     *HERRI2(M,IGAUSS,KARTEZ,LTWCEN)
     *                     *  PNTTWC(IGAUSS,KARTEZ,LTWCEN)**K
C
                  END DO
C
                 CMUTWC(K,N,M,KARTEZ,LTWCEN)=RESULT
     *                         *FACNE2(KARTEZ,LTWCEN)
C
               END DO
            END DO
         END DO
C
      END DO
C
C=======================================================================
C     HERE ARE COMPUTED THE COEFFICIENTS FOR THE SHIFTED MULTIPOLE
C     DEFORMATIONS (EACH CENTER)
C=======================================================================
C
      DO LTWCEN=1,NDTWBL
C
         DO N=0,NOSACT
            DO M=0,NOSACT
               DO K=0,NMUPOW
C
                  RESUL1=0.0D0
                  RESUL2=0.0D0
C
                  DO IGAUSS=1,NGAUSS
C
                     RESUL1 =             RESUL1 + WGHERM(IGAUSS,KARTEZ)
     *       *                            HERLE2(N,IGAUSS,KARTEZ,LTWCEN)
     *       *                            HERRI2(M,IGAUSS,KARTEZ,LTWCEN)
     *       *        (PNTTWC(IGAUSS,KARTEZ,LTWCEN)-CENCOR(KARTEZ,1))**K
C
                     RESUL2 =             RESUL2 + WGHERM(IGAUSS,KARTEZ)
     *       *                            HERLE2(N,IGAUSS,KARTEZ,LTWCEN)
     *       *                            HERRI2(M,IGAUSS,KARTEZ,LTWCEN)
     *       *   (PNTTWC(IGAUSS,KARTEZ,LTWCEN)-CENCOR(KARTEZ,NDTWCE))**K
C
                  END DO
C
                  CMUSH1(K,N,M,KARTEZ,LTWCEN)=RESUL1
     *                        * FACNE2(KARTEZ,LTWCEN)
C
                  CMUSH2(K,N,M,KARTEZ,LTWCEN)=RESUL2
     *                        * FACNE2(KARTEZ,LTWCEN)
C
               END DO
            END DO
         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:NDTWDD),
     *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD)
      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:NDTWDD),
     *                CDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                CDDLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                CHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                CDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                CDDRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD)
C
C=======================================================================
       IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_COPHER'
C=======================================================================
C
      NMAIN2=2*NOSACT+2
C
      DO INDBLO = 1, NDTWDD
         DO IZEROS = 1,NGAUSS
            DO NOSCIL = 0,NMAIN2
C
               IF (INDBLO.LE.NDTWBL) THEN
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)
               END IF
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=======================================================================
C
      SUBROUTINE TWC_MORDEN(INDCOE,ITPNMX,IGRAIN,IPAIRI)
C=======================================================================
      USE hfodd_sizes
      USE HE_DEN
      USE HETDEN
      USE PD_DEN
      USE PP_DEN
      USE PDTDEN
C=======================================================================
      COMMON
     *      /INFAST/  IDCOPY(1:NDTWDD,0:1),IDCOCO(1:NDTWDD,0:1)
C
C=======================================================================
C     HERE WE COPY DIFFERENT DENSITIES FROM DENSHF INTO DIFFERENT CENTER
C     INDICES TAKING ADVANTAGE OF THE SYMMETRIES WHEN WAVE FUNCTIONS
C     ARE EVALUATED IN NEW COMMON LATTICE OF INTEGRATION
C=======================================================================
C
      IF (IDCOCO(INDCOE,IPAIRI).EQ.0) THEN
C
       DE_RHO(:,:,:,ITPNMX,INDCOE)=
     * DE_RHO(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_TAU(:,:,:,ITPNMX,INDCOE)=
     * DE_TAU(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_LPR(:,:,:,ITPNMX,INDCOE)=
     * DE_LPR(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_DIV(:,:,:,ITPNMX,INDCOE)=
     * DE_DIV(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       DE_SPI(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_SPI(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_KIS(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_KIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_LPS(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_LPS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_GRR(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_GRR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_CUR(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_CUR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_KIF(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_KIF(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_ROS(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_ROS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_ROC(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_ROC(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       DE_SCU(:,:,:,:,:,ITPNMX,INDCOE)=
     *            DE_SCU(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_DES(:,:,:,:,:,ITPNMX,INDCOE)=
     *            DE_DES(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       IF (IGRAIN.EQ.1) THEN
C
       DE_DIJ(:,:,:,ITPNMX,INDCOE)=
     *            DE_DIJ(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       DE_DIS(:,:,:,:,ITPNMX,INDCOE)=
     *            DE_DIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       END IF
C
C=======================================================================
C      PAIRING DENSITIES
C=======================================================================
C
       PD_RHO(:,:,:,ITPNMX,INDCOE)=
     * PD_RHO(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       PD_TAU(:,:,:,ITPNMX,INDCOE)=
     * PD_TAU(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       PD_LPR(:,:,:,ITPNMX,INDCOE)=
     * PD_LPR(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))

C
       PD_GRR(:,:,:,:,ITPNMX,INDCOE)=
     * PD_GRR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       PD_SCU(:,:,:,:,:,ITPNMX,INDCOE)=
     * PD_SCU(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       PP_RHO(:,:,:,ITPNMX,INDCOE)=
     * PP_RHO(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       PP_TAU(:,:,:,ITPNMX,INDCOE)=
     * PP_TAU(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       PP_LPR(:,:,:,ITPNMX,INDCOE)=
     * PP_LPR(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       PP_GRR(:,:,:,:,ITPNMX,INDCOE)=
     * PP_GRR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
       PP_SCU(:,:,:,:,:,ITPNMX,INDCOE)=
     * PP_SCU(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       IF (IGRAIN.EQ.1) THEN
C
       PD_DIS(:,:,:,:,ITPNMX,INDCOE)=
     * PD_DIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       PP_DIS(:,:,:,:,ITPNMX,INDCOE)=
     *PD_DIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI))
C
       END IF
C
      !COPIES THE COMPLEX CONJUGATE
      ! NOT APPLICABLE TO PAIRING DENSITIES
      ELSEIF (IDCOCO(INDCOE,IPAIRI).EQ.1) THEN
C
       DE_RHO(:,:,:,ITPNMX,INDCOE)=
     *            CONJG(DE_RHO(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_TAU(:,:,:,ITPNMX,INDCOE)=
     *            CONJG(DE_TAU(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_LPR(:,:,:,ITPNMX,INDCOE)=
     *            CONJG(DE_LPR(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_DIV(:,:,:,ITPNMX,INDCOE)=
     *            CONJG(DE_DIV(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       DE_SPI(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_SPI(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_KIS(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_KIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_LPS(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_LPS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_GRR(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_GRR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_CUR(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_CUR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_KIF(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_KIF(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_ROS(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_ROS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_ROC(:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_ROC(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       DE_SCU(:,:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_SCU(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       DE_DES(:,:,:,:,:,ITPNMX,INDCOE)=
     *  CONJG(DE_DES(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       IF (IGRAIN.EQ.1) THEN
           DE_DIJ(:,:,:,ITPNMX,INDCOE)=
     *            CONJG(DE_DIJ(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
           DE_DIS(:,:,:,:,ITPNMX,INDCOE)=
     *            CONJG(DE_DIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       END IF
C
C=======================================================================
C      PAIRING DENSITIES
C=======================================================================
C
       PD_RHO(:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PD_RHO(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       PD_TAU(:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PD_TAU(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       PD_LPR(:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PD_LPR(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))

C
       PD_GRR(:,:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PD_GRR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       PD_SCU(:,:,:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PD_SCU(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       PP_RHO(:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PP_RHO(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       PP_TAU(:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PP_TAU(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       PP_LPR(:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PP_LPR(:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       PP_GRR(:,:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PP_GRR(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
       PP_SCU(:,:,:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PP_SCU(:,:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       IF (IGRAIN.EQ.1) THEN
C
       PD_DIS(:,:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PD_DIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       PP_DIS(:,:,:,:,ITPNMX,INDCOE)=
     *          CONJG(PD_DIS(:,:,:,:,ITPNMX,IDCOPY(INDCOE,IPAIRI)))
C
       END IF
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWCCOP(ICHARG,LTWCEN)
C=======================================================================
      USE MAT_PP
      USE MAT_PM
      USE TWCEPP
      USE TWCEPM
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C=======================================================================
C
      IALLOC=0
C
C=======================================================================
      IF (.NOT.ALLOCATED(TWC_PP).AND.LTWCEN.GE.1) THEN
          ALLOCATE (TWC_PP(1:2*NDBASE,1:2*NDBASE,0:NDREVE,
     *                                           0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWC_PP','TWCCOP')
      END IF
C=======================================================================
      IF (.NOT.ALLOCATED(TWC_PM).AND.LTWCEN.GE.1) THEN
          ALLOCATE (TWC_PM(1:2*NDBASE,1:2*NDBASE,0:NDREVE,
     *                                           0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWC_PM','TWCCOP')
      END IF
C=======================================================================
C          THIS SUBROUTINE COPIES THE CALCULATED MATRIX ELEMENTS OF THE
C          MEAN-FIELD INTO THE TWO-CENTER STRUCTURE.
C=======================================================================
C
      IF (LTWCEN.EQ.1) THEN

          TWC_PP=CMPLX(0.0D0,0.0D0)
          TWC_PM=CMPLX(0.0D0,0.0D0)
C
          MBRA=0
          MKET=0
C
      END IF
C
      IF (LTWCEN.EQ.2) THEN
C
          MBRA=LDBATW(1)
          MKET=LDBATW(1)
C
      END IF
C
      IF (LTWCEN.EQ.3) THEN
C
          MBRA=LDBATW(1)
          MKET=0
C
      END IF
C
C=======================================================================
C     SETTING THE MATRIX ELEMENTS
C=======================================================================
C
      DO IKET=1,LDBATW(IND2HR(LTWCEN))
         DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
            TWC_PP(MBRA+IBRA,MKET+IKET,0,ICHARG)=BIG_PP(IBRA,IKET,0)
            TWC_PP(MBRA+IBRA,MKET+IKET,1,ICHARG)=BIG_PP(IBRA,IKET,1)
C
            TWC_PM(MBRA+IBRA,MKET+IKET,0,ICHARG)=BIG_PM(IBRA,IKET,0)
            TWC_PM(MBRA+IBRA,MKET+IKET,1,ICHARG)=BIG_PM(IBRA,IKET,1)
C
         END DO
      END DO
C
      IF (LTWCEN.EQ.3) THEN
C
          MBRA=0
          MKET=LDBATW(1)
C
          DO IKET=1,LDBATW(IND2HL(LTWCEN))
             DO IBRA=1,LDBATW(IND2HR(LTWCEN))
C
                TWC_PP(MBRA+IBRA,MKET+IKET,0,ICHARG)=
     *                                CONJG(BIG_PP(IKET,IBRA,0))
                TWC_PP(MBRA+IBRA,MKET+IKET,1,ICHARG)=
     *                                CONJG(BIG_PP(IKET,IBRA,1))
C
                TWC_PM(MBRA+IBRA,MKET+IKET,0,ICHARG)=
     *                                CONJG(BIG_PM(IKET,IBRA,1))
                TWC_PM(MBRA+IBRA,MKET+IKET,1,ICHARG)=
     *                                CONJG(BIG_PM(IKET,IBRA,0))
C
             END DO
          END DO
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      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+6)
      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:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *       /SCALNG/ HOMSCA(1:NDKART,1:NDTWCE)
      COMMON
     *       /CFIPRI/ NFIPRI
      COMMON
     *       /CENPOS/ CENCOR(1:NDKART,1:NDTWCE)
      COMMON
     *       /TWINTF/ FACNE2(1:NDKART,1:NDTWBL),
     *                FACNE4(1:NDKART,1:NDTWHE)
C
C=======================================================================
C     YEP, YOU PROBLABLY HAVE ALREADY FIGURED IT OUT... THIS DOES THE
C     SAME AS DEVHER. PRETTY ORIGINAL, ISN'T IT?
C     AS WE WILL NEED HIGHER ORDER DERIVATIVES, WE COMPUTE EVERYTHING
C     AGAIN, BUT THIS WILL BE SIMPLIFIED IN THE NEAR FUTURE.
C=======================================================================
C
      IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_DEVHER'
C
C=======================================================================
C
      NMAIN1=2*NOSACT+6
      NGAUSS=2*NOSACT+7
C
      CF2H00(:,:,:,KARTEZ,:)=0.0D0
      CF2H01(:,:,:,KARTEZ,:)=0.0D0
      CF2H02(:,:,:,KARTEZ,:)=0.0D0
C
      CF4H00(:,:,:,KARTEZ,:)=0.0D0
      CF4H01(:,:,:,KARTEZ,:)=0.0D0
      CF4H11(:,:,:,KARTEZ,:)=0.0D0
      CF4H02(:,:,:,KARTEZ,:)=0.0D0
C
C=======================================================================
C
      CALL HERMIT(NGAUSS,XHERMI,WHERMI,NGAUSS)
      CALL NORFAC(HERFAC,NMAIN1)
C
      DO INDBLO=1,NDTWBL
C
         HOSCAI=HOMSCA(KARTEZ,IND2HL(INDBLO))
         HOSCAJ=HOMSCA(KARTEZ,IND2HR(INDBLO))
C
         HOSHII=CENCOR(KARTEZ,IND2HL(INDBLO))
         HOSHIJ=CENCOR(KARTEZ,IND2HR(INDBLO))

         SUMOSL=HOSCAI**2+HOSCAJ**2
         SUMCEN=HOSCAI**2*HOSHII+HOSCAJ**2*HOSHIJ
C
         FACTO1=SQRT(2.D0/SUMOSL)
         FACTO2=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,NMAIN1+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,NMAIN1+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,NMAIN1+1)
C
            DO NOSCIL = 0,NMAIN1
C
                HERAUX(NOSCIL,IZEROS)= PHERMI(NOSCIL+1)/
     *                                   HERFAC(NOSCIL)
C
            END DO
C
         END DO
C
C=======================================================================
C        MAIN LOOPS IN QUANTM NUMBERS N AND M
C=======================================================================
C
         DO M=0,NOSACT
            DO N=0,NOSACT
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 00
C=======================================================================
C
               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
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 01
C=======================================================================
C
               DO K=0,N+M+1
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
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 02
C=======================================================================
C
               DO K=0,N+M+2
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
C
              END DO
          END DO
C
      END DO
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS WHEN 4 WAVE FUNCTIONS ARE USED.
C        BEHOLD! HERE ARE COMPUTED THE "OUTER WAVE FUNCTIONS" OF THE
C        SANDWICHES IN INTEGH, DON'T CONFUSE THEM WITH THOSE OBTAINED
C        IN DEFINT, WHICH WILL BE USED IN DENSHF.
C=======================================================================
C
      DO INDBLO=1,NDTWHE
C
         HOSCAI=HOMSCA(KARTEZ,IND4HI(INDBLO))
         HOSCAJ=HOMSCA(KARTEZ,IND4HJ(INDBLO))
         HOSCAK=HOMSCA(KARTEZ,IND4HK(INDBLO))
         HOSCAL=HOMSCA(KARTEZ,IND4HL(INDBLO))
C
         HOSHII=CENCOR(KARTEZ,IND4HI(INDBLO))
         HOSHIJ=CENCOR(KARTEZ,IND4HJ(INDBLO))
         HOSHIK=CENCOR(KARTEZ,IND4HK(INDBLO))
         HOSHIL=CENCOR(KARTEZ,IND4HL(INDBLO))
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=HOSCAI*SQRT(2.D0/SUMOSL)
         FACLE2=HOSCAI*(SUMCEN/SUMOSL-HOSHII)
C
         FACRI1=HOSCAJ*SQRT(2.D0/SUMOSL)
         FACRI2=HOSCAJ*(SUMCEN/SUMOSL-HOSHIJ)
C
         DO IZEROS = 1,NGAUSS
C
            XZER = XHERMI(IZEROS)*FACLE1+FACLE2
C
            CALL D_HERM(XZER,NMAIN1,PHERMI,DHERMI,NMAIN1+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,NMAIN1+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
         DO N=0,NOSACT
            DO M=0,NOSACT
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 00
C=======================================================================
C
               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
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 01
C=======================================================================
C
               DO K=0,N+M+1
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
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS 11
C=======================================================================
C
               DO K=0,N+M+2
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
C
            END DO
         END DO
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_NORMAT(IEIGCU,EIGCUT)
C=======================================================================
      USE hfodd_sizes
      USE OVM_PP
      USE OVM_PM
      USE NOREIG
      USE BIGOVE
C=======================================================================
      DIMENSION
     *          AUXTW0(1:2*NDBASE,1:2*NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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
     *       /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,NATDIM
      COMMON
     *       /EGVOVM/ EGVNOR(1:4*NDBASE)
      COMMON
     *       /CFIPRI/ NFIPRI
      COMPLEX
     *          CELNOR,C_ZERO
      COMPLEX
     *          RESULT,AUXTW0
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT

C=======================================================================
      ALLOCATABLE CELNOR(:)
C=======================================================================
C
      C_ZERO=CMPLX(0.0D0,0.0D0)
      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,LDBATW(IND2HL(LTWCEN))
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBATW(IND2HR(LTWCEN))
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+LDBATW(1),IKET+LDBATW(1))
     *                             = RESULT
                   IF(LTWCEN.EQ.3) AUXTW0(IBRA+LDBATW(1),IKET)
     *                             = RESULT
                   IF(LTWCEN.EQ.3) AUXTW0(IKET,IBRA+LDBATW(1))
     *                             = RESULT
            END DO
         END DO
C
      END DO
C=======================================================================
      IF (.NOT.ALLOCATED(OVMAPP)) THEN
          ALLOCATE (OVMAPP(1:2*NDBASE,1:2*NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('OVMAPP','NORMAT')
      END IF
      IF (.NOT.ALLOCATED(OVMAPM)) THEN
          ALLOCATE (OVMAPM(1:2*NDBASE,1:2*NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('OVMAPM','NORMAT')
      END IF
C=======================================================================
      DO IBRA=1,LDBTOT
         DO IKET=1,LDBTOT
C
         IF (IBRA.LE.LDBATW(1)) MY=NYVECT(IBRA,IND2HL(1))
         IF (IBRA.GT.LDBATW(1)) MY=NYVECT(IBRA-LDBATW(1),IND2HL(NDTWCE))
C
         IF (IKET.LE.LDBATW(1)) NY=NYVECT(IKET,IND2HR(1))
         IF (IKET.GT.LDBATW(1)) NY=NYVECT(IKET-LDBATW(1),IND2HR(NDTWCE))
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)=AUXTW0(IBRA,IKET)*IPHAPM(MY,NY,0)
         OVMAPM(IBRA,IKET,1)=AUXTW0(IBRA,IKET)*IPHAMP(MY,NY,0)
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
      IF (.NOT.ALLOCATED(BIGNOR)) THEN
          ALLOCATE (BIGNOR(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('BIGNOR','NORMAT')
      END IF
C=======================================================================
      DO IBRA=1,2*LDBTOT
         DO IKET=1,2*LDBTOT
C
            IF (IBRA.GE.IKET) THEN
C
                NCOUNT=IBRA+((4*LDBTOT-IKET)*(IKET-1))/2
C
                IF (IBRA.LE.LDBTOT.AND.IKET.LE.LDBTOT)
     *
     *              CELNOR(NCOUNT)=OVMAPP(IBRA,IKET,0)
C
                IF (IBRA.GT.LDBTOT.AND.IKET.LE.LDBTOT)
     *
     *              CELNOR(NCOUNT)=OVMAPM(IBRA-LDBTOT,IKET,1)
C
                IF (IBRA.GT.LDBTOT.AND.IKET.GT.LDBTOT)
     *
     *              CELNOR(NCOUNT)=OVMAPP(IBRA-LDBTOT,IKET-LDBTOT,1)
C
            END IF
C
         END DO
      END DO
C
      BIGNOR=0.0D0
      DO IBRA=1,2*LDBTOT
         DO IKET=1,2*LDBTOT
C
                IF (IBRA.LE.LDBTOT.AND.IKET.LE.LDBTOT)
     *
     *              BIGNOR(IBRA,IKET)=OVMAPP(IBRA,IKET,0)

                IF (IBRA.GT.LDBTOT.AND.IKET.LE.LDBTOT)
     *
     *              BIGNOR(IBRA,IKET)=OVMAPM(IBRA-LDBTOT,IKET,1)

                IF (IBRA.LE.LDBTOT.AND.IKET.GT.LDBTOT)
     *
     *              BIGNOR(IBRA,IKET)=OVMAPM(IBRA,IKET-LDBTOT,0)

                IF (IBRA.GT.LDBTOT.AND.IKET.GT.LDBTOT)
     *
     *          BIGNOR(IBRA,IKET)=OVMAPP(IBRA-LDBTOT,IKET-LDBTOT,1)
C
         END DO
      END DO
c
      CALL DIAMAT(CELNOR,EGVNOR,EWTNOR,2*LDBTOT,4*NDBASE,2*LDBTOT)
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=2*LDBTOT,1,-1
C
          IF(EGVNOR(IBAS).GT.EIGCUT) ICUTOF=IBAS
C
      END DO

      IF (ICUTOF.EQ.0) STOP 'ICUTOF=0, NO STATES IN THE BASIS :('
C
      ICOUNT=0
C
      DO INDNAT=ICUTOF,2*LDBTOT
      ICOUNT=ICOUNT+1
         DO IBASE=1,2*LDBTOT
            EWTAUX(IBASE,ICOUNT)=
     *      EWTNOR(IBASE,INDNAT)/SQRT(EGVNOR(INDNAT))
         END DO
      END DO
C
      NATDIM=ICOUNT
C
      WRITE(NFIPRI,'(79(1H*),/,1H*,77X,1H*,/,
     *                   1H*,1X,''CUT-OFF IN THE SMALLEST NORM '',
     *          ''EIGENVALUE            ='',E12.3,12X,1H*)')
     * EIGCUT

      WRITE(NFIPRI,'(1H*,1X,''NUMBER OF STATES IN NATURAL BASIS'',
     *                          ''                  ='',I12,12X,1H*,/,
     *                                                    1H*,77X,1H*)')
     * NATDIM
C
C=======================================================================
C
      RETURN
      END
C
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)
C
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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,LDBATW(IND2HL(LTWCEN))
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMLIN(IBRA,IKET)= FACINP * IPHASE(MY,NY) * 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
C=======================================================================
C
      RETURN
      END
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)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
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,LDBATW(IND2HL(LTWCEN))
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMLIN(IBRA,IKET)= FACINP * IPHASE(MY,NY) * 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=======================================================================
C
      RETURN
      END
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
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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,LDBATW(IND2HL(LTWCEN))
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMLIN(IBRA,IKET)= FACINP * IPHASE(MY,NY) * 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
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_INTKIN(HAMKIN,HBMACT,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
C
      DIMENSION
     *          HAMKIN(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
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,LDBATW(IND2HL(LTWCEN))
C
            MX=NXVECT(IBRA,IND2HL(LTWCEN))
            MY=NYVECT(IBRA,IND2HL(LTWCEN))
            MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
            DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
               NX=NXVECT(IKET,IND2HR(LTWCEN))
               NY=NYVECT(IKET,IND2HR(LTWCEN))
               NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
               HAMKIN(IBRA,IKET) = HBMACT * FACTOR
     *       * (
     *           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(HAMANG,FACINP,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
C
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=-PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
         MX=NXVECT(IBRA,IND2HL(LTWCEN))
         MY=NYVECT(IBRA,IND2HL(LTWCEN))
         MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
         DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
            NX=NXVECT(IKET,IND2HR(LTWCEN))
            NY=NYVECT(IKET,IND2HR(LTWCEN))
            NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
            HAMANG(IBRA,IKET) = FACINP * FACTOR *
     *                                 CF2H00(0,MX,NX,1,LTWCEN) *
     *      (
     *        CF2H01(0,MZ,NZ,3,LTWCEN) *
     *      ( CF2H00(1,MY,NY,2,LTWCEN)*QUADSUM(2,LTWCEN)/SQRT(2.D0) +
     *                 CF2H00(0,MY,NY,2,LTWCEN)*QUADCEN(2,LTWCEN) )
     *      - CF2H01(0,MY,NY,2,LTWCEN) *
     *      ( CF2H00(1,MZ,NZ,3,LTWCEN)*QUADSUM(3,LTWCEN)/SQRT(2.D0) +
     *        CF2H00(0,MZ,NZ,3,LTWCEN)*QUADCEN(3,LTWCEN) )
     *      )
C
         END DO
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_LY(HAMANG,FACINP,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=-PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
         MX=NXVECT(IBRA,IND2HL(LTWCEN))
         MY=NYVECT(IBRA,IND2HL(LTWCEN))
         MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
         DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
            NX=NXVECT(IKET,IND2HR(LTWCEN))
            NY=NYVECT(IKET,IND2HR(LTWCEN))
            NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
            HAMANG(IBRA,IKET) = FACINP * FACTOR *
     *                                 CF2H00(0,MY,NY,2,LTWCEN) *
     *      (
     *        CF2H01(0,MX,NX,1,LTWCEN) *
     *      ( CF2H00(1,MZ,NZ,3,LTWCEN)*QUADSUM(3,LTWCEN)/SQRT(2.D0) +
     *                 CF2H00(0,MZ,NZ,3,LTWCEN)*QUADCEN(3,LTWCEN) )
     *      - CF2H01(0,MZ,NZ,3,LTWCEN) *
     *      ( CF2H00(1,MX,NX,1,LTWCEN)*QUADSUM(1,LTWCEN)/SQRT(2.D0) +
     *        CF2H00(0,MX,NX,1,LTWCEN)*QUADCEN(1,LTWCEN) )
     *      )
C
         END DO
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_LZ(HAMANG,FACINP,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=-PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
         MX=NXVECT(IBRA,IND2HL(LTWCEN))
         MY=NYVECT(IBRA,IND2HL(LTWCEN))
         MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
         DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
            NX=NXVECT(IKET,IND2HR(LTWCEN))
            NY=NYVECT(IKET,IND2HR(LTWCEN))
            NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
            HAMANG(IBRA,IKET) = FACINP * FACTOR *
     *                                 CF2H00(0,MZ,NZ,3,LTWCEN) *
     *      (
     *        CF2H01(0,MY,NY,2,LTWCEN) *
     *      ( CF2H00(1,MX,NX,1,LTWCEN)*QUADSUM(1,LTWCEN)/SQRT(2.D0) +
     *                 CF2H00(0,MX,NX,1,LTWCEN)*QUADCEN(1,LTWCEN) )
     *      - CF2H01(0,MX,NX,1,LTWCEN) *
     *      ( CF2H00(1,MY,NY,2,LTWCEN)*QUADSUM(2,LTWCEN)/SQRT(2.D0) +
     *        CF2H00(0,MY,NY,2,LTWCEN)*QUADCEN(2,LTWCEN) )
     *      )
C
            END DO
         END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_SX(HAMANG,FACINP,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)*
     *       FACINP*0.5D0
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
         MX=NXVECT(IBRA,IND2HL(LTWCEN))
         MY=NYVECT(IBRA,IND2HL(LTWCEN))
         MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
         DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
            NX=NXVECT(IKET,IND2HR(LTWCEN))
            NY=NYVECT(IKET,IND2HR(LTWCEN))
            NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
            HAMANG(IBRA,IKET) =  FACTOR
     *                        * CF2H00(0,MX,NX,1,LTWCEN)
     *                        * CF2H00(0,MY,NY,2,LTWCEN)
     *                        * CF2H00(0,MZ,NZ,3,LTWCEN)
C
            END DO
         END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
      SUBROUTINE TWC_SY(HAMANG,FACINP,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)*
     *       FACINP*0.5D0
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
         MX=NXVECT(IBRA,IND2HL(LTWCEN))
         MY=NYVECT(IBRA,IND2HL(LTWCEN))
         MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
         DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
            NX=NXVECT(IKET,IND2HR(LTWCEN))
            NY=NYVECT(IKET,IND2HR(LTWCEN))
            NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
            HAMANG(IBRA,IKET) = FACTOR
     *                        * CF2H00(0,MX,NX,1,LTWCEN)
     *                        * CF2H00(0,MY,NY,2,LTWCEN)
     *                        * CF2H00(0,MZ,NZ,3,LTWCEN)
C
            END DO
         END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_SZ(HAMANG,FACINP,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          HAMANG(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                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
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      PIFAC=(4.0D0*ATAN(1.00D0))**(3.D0/4.D0)
      FACTOR=PIFAC*FACNE2(1,LTWCEN)*FACNE2(2,LTWCEN)*FACNE2(3,LTWCEN)*
     *       FACINP*0.5D0
C
C=======================================================================
C
      HAMANG(:,:)=0.D0
C
      DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
         MX=NXVECT(IBRA,IND2HL(LTWCEN))
         MY=NYVECT(IBRA,IND2HL(LTWCEN))
         MZ=NZVECT(IBRA,IND2HL(LTWCEN))
C
         DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
            NX=NXVECT(IKET,IND2HR(LTWCEN))
            NY=NYVECT(IKET,IND2HR(LTWCEN))
            NZ=NZVECT(IKET,IND2HR(LTWCEN))
C
            HAMANG(IBRA,IKET) = FACTOR
     *                        * CF2H00(0,MX,NX,1,LTWCEN)
     *                        * CF2H00(0,MY,NY,2,LTWCEN)
     *                        * CF2H00(0,MZ,NZ,3,LTWCEN)
C
            END DO
         END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
      SUBROUTINE TWCSIZ_NATBAS(ICHARG,LDMEFI,IPAHFB,EFERMI,ECUTOF,
     *                         LAMCUT,INSIZB,IDSIZB,IFIBLB,ISABLB,
     *                                       IREMQB,INUMQB,JNUMQB)
C=======================================================================
      USE hfodd_sizes
      USE OVM_PP
      USE OVM_PM
      USE TWCEPP
      USE TWCEPM
      USE NOREIG
      USE TWCWAV
      USE BIGOVE
      USE TWCDEL
      USE TWCQUZ
      USE BIGFRA
      USE BIGHAM
C=======================================================================
      COMPLEX
     *          C_ZERO,C_UNIT
      COMPLEX
     *          CELMTS,CELHFB
      COMPLEX
     *          HAMAUX,CCWEWT
      COMPLEX
     *          COMHA2,HAMCUT,HAMCU2,
     *          COMDEL,DELCUT,DELCU2
      COMPLEX
     *          EWTNAT,WAVTWC
      COMPLEX
     *          EQUNAT
      COMPLEX
     *          AWAAUX,BWAAUX,AWANAT,BWANAT
      COMPLEX
     *          RESAUX
      DIMENSION
     *          SPENAT(1:4*NDBASE),EQUTWC(1:8*NDBASE)
      DIMENSION
     *          AWAAUX(1:4*NDBASE,1:2*NDSTAT),
     *          BWAAUX(1:4*NDBASE,1:2*NDSTAT)
      DIMENSION
     *          AWANAT(1:4*NDBASE,1:2*NDSTAT),
     *          BWANAT(1:4*NDBASE,1:2*NDSTAT)
      DIMENSION
     *          IDSTAT(0:NDREVE)
      DIMENSION
     *          EORDER(1:8*NDBASE),NORDER(1:8*NDBASE),
     *          NSTATE(1:8*NDBASE)
      DIMENSION
     *          KSTATE(1:2*NDSTAT,0:NDREVE)
      DIMENSION
     *          INSIZB(1:NDBLOC),IDSIZB(1:NDBLOC)
      DIMENSION
     *          IREMQB(1:NDBLOC),INUMQB(1:NDBLOC),JNUMQB(1:NDBLOC)
      DIMENSION
     *          CCWEWT(1:4*NDBASE,1:4*NDBASE)
C=======================================================================
      COMMON
     *       /NATCUT/ ICUTOF,NATDIM
      COMMON
     *       /SPNUMS/ NUMBSP(0:NDREVE,0:NDISOS)
      COMMON
     *       /TWCDIA/ TWCSPE(1:2*NDSTAT,0:NDISOS)
      COMMON
     *       /SCALNG/ HOMSCA(1:NDKART,1:NDTWCE)
      COMMON
     *       /CENPOS/ CENCOR(1:NDKART,1:NDTWCE)
      COMMON
     *       /EGVOVM/ EGVNOR(1:4*NDBASE)
      COMMON
     *       /PARBAS/ IBAPAR(1:NDBASE,1:NDTWCE)
      COMMON
     *       /QPNUMS/ NUMBQP(0:NDREVE,0:NDISOS)
      COMMON
     *       /XEQUIV/ EEQUIV(1:2*NDSTAT,0:NDREVE),
     *                DEQUIV(1:2*NDSTAT,0:NDREVE)
      COMMON
     *       /QUAPAR/ ENQUAP(1:2*NDSTAT,0:NDREVE),
     *                V2QUAP(1:2*NDSTAT,0:NDREVE)
      COMMON
     *       /SIMPQU/ QPSIMP(1:2*NDSTAT,0:NDREVE)
      COMMON
     *       /CFIPRI/ NFIPRI
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C=======================================================================
      ALLOCATABLE COMHA2(:,:),
     *            COMDEL(:,:)
C
      ALLOCATABLE HAMCUT(:,:),HAMCU2(:,:),
     *            DELCUT(:,:),DELCU2(:,:)

      ALLOCATABLE CELMTS(:), CELHFB(:)
C
      ALLOCATABLE EWTNAT(:,:),EQUNAT(:,:)
C
      ALLOCATABLE HAMAUX(:,:)
C
      ALLOCATABLE WAVTWC(:,:)
C=======================================================================
      C_ZERO=CMPLX(0.0D0,0.0D0)
      C_UNIT=CMPLX(1.0D0,0.0D0)
      ONE_16=1.0D0
C=======================================================================
C     BUILDING THE PROPER SIMPLEX STRUCTURE OF THE HAMILTONIAN
C=======================================================================
      IALLOC=0
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
      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
      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
      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=======================================================================
C     BUILDING THE MATRIX ELEMENTS OF THE CONSTRAINTS ON THE FRAGMENTS
C=======================================================================
C
      CALL TWC_FRGMAT(ICHARG)
C
C=======================================================================
C
      COMHAM(:,:)=C_ZERO
C
      DO IBRA=1,2*LDBTOT
         DO IKET=1,2*LDBTOT
C
            IF (IBRA.LE.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
                COMHAM(IBRA,IKET)=TWC_PP(IBRA,IKET,0,ICHARG)
C
            ELSEIF (IBRA.GT.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
                COMHAM(IBRA,IKET)=TWC_PM(IBRA-LDBTOT,IKET,1,ICHARG)
C
            ELSEIF (IBRA.LE.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
                COMHAM(IBRA,IKET)=TWC_PM(IBRA,IKET-LDBTOT,0,ICHARG)
C
            ELSEIF (IBRA.GT.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
                COMHAM(IBRA,IKET)=
     *          TWC_PP(IBRA-LDBTOT,IKET-LDBTOT,1,ICHARG)
C
            END IF
C
            COMHAM(IBRA,IKET)=COMHAM(IBRA,IKET)-FRAMAT(IBRA,IKET)
C
         END DO
      END DO
C
C=======================================================================
C     HF CASE
C=======================================================================
C
      IALLOC=0
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
      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
      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
      IF (.NOT.ALLOCATED(TWCOCC)) THEN
          ALLOCATE (TWCOCC(1:2*NDBASE,1:2*NDSTAT,0:NDREVE,0:NDISOS),
     *                                                 STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWCOCC','TWCSIZ_NATBAS')
      END IF
C=======================================================================
C
      HAMAUX(:,:)=C_ZERO
      HAMCUT(:,:)=C_ZERO
C
      CALL ZGEMM('C','N',NATDIM,2*LDBTOT,2*LDBTOT,
     *            C_UNIT,EWTAUX,4*NDBASE,
     *                          COMHAM,4*NDBASE,
     *                   C_ZERO,HAMAUX,4*NDBASE)

C
      CALL ZGEMM('N','N',NATDIM,NATDIM,2*LDBTOT,
     *            C_UNIT,HAMAUX,4*NDBASE,
     *                          EWTAUX,4*NDBASE,
     *                   C_ZERO,HAMCUT,4*NDBASE)
C
C=======================================================================
C     DIAGONALIZING H IN THE "NATURAL BASIS"
C=======================================================================
C
C=======================================================================
C     INSTEAD OF COMPUTING THE CANONICAL BASIS, WE DIAGONALIZE THE
C     P-H CHANNEL TO COMPUTE THE SINGLE-PARTICLE ENERGIES AND THE
C     DIFFERENT SINGLE-PARTICLE PROPERTIES.
C=======================================================================
C
      CELMTS(:)=C_ZERO
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
C=======================================================================
C     INSTEAD OF COMPUTING THE CANONICAL BASIS, WE DIAGONALIZE THE
C     P-H CHANNEL TO COMPUTE THE SINGLE-PARTICLE ENERGIES AND THE
C     DIFFERENT SINGLE-PARTICLE PROPERTIES.
C=======================================================================
C
      LDMEFI=NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
C
      IF (LDMEFI.GT.NATDIM) STOP 'LDMEFI.GT.NATDIM IN TWCSIZ_NATBAS'

      CALL DIAMAT(CELMTS,SPENAT,EWTNAT,NATDIM,4*NDBASE,NATDIM)
C
C=======================================================================
C     GOING BACK TO THE HO BASIS
C=======================================================================
C
      WAVTWC(:,:)=C_ZERO
C
      CALL ZGEMM('N','N',2*LDBTOT,NATDIM,NATDIM,
     *                   C_UNIT,EWTAUX,4*NDBASE,
     *                          EWTNAT,4*NDBASE,
     *                   C_ZERO,WAVTWC,4*NDBASE)
C
C=======================================================================
C     SAVING ENERGIES AND AMPLITUDES OF MANY BODY WAVE FUNCTION
C=======================================================================
C
      DO IALLST=1,LDMEFI
C
         TWCSPE(IALLST,ICHARG)=SPENAT(IALLST)
C
         DO IBASE=1,LDBTOT
            TWCOCC(IBASE,IALLST,0,ICHARG)=WAVTWC(IBASE       ,IALLST)
            TWCOCC(IBASE,IALLST,1,ICHARG)=WAVTWC(IBASE+LDBTOT,IALLST)
         END DO
C
      END DO
C
C=======================================================================
C     HFB CASE:
C
C     1. SUBSTRACTION THE FERMI ENERGY TERM FROM THE MEAN FIELD
C     2. TRANSFORMATION P-H AND P-P INTO THE NATURAL BASIS
C
C          | U     0   | | H       D   |  | U^dagger     0   |
C          | 0     U^* | | -D^*   -H^* |  | 0            U^T |
C
C     3. DIAGONALIZATION OF THE MATRIX
C
C    |  U   H  U^dagger       U   D  U^T  |   |  HAMCUT       DELCUT  |
C    |                                    | = |                       |
C    | -U^* D  U^dagger      -U^* H  U^T  |   |  DELCU2       HAMCU2  |
C
C=======================================================================
C
      IF (IPAHFB.EQ.1) THEN
C
C=======================================================================
C
      IALLOC=0
      IF (.NOT.ALLOCATED(COMHA2)) THEN
          ALLOCATE (COMHA2(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('COMHA2','TWCSIZ_NATBAS')
      END IF
      IF (.NOT.ALLOCATED(HAMCU2)) THEN
          ALLOCATE (HAMCU2(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('HAMCU2','TWCSIZ_NATBAS')
      END IF
      IF (.NOT.ALLOCATED(COMDEL)) THEN
          ALLOCATE (COMDEL(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('COMDEL','TWCSIZ_NATBAS')
      END IF
      IF (.NOT.ALLOCATED(DELCUT)) THEN
          ALLOCATE (DELCUT(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('DELCUT','TWCSIZ_NATBAS')
      END IF
      IF (.NOT.ALLOCATED(DELCU2)) THEN
          ALLOCATE (DELCU2(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('DELCU2','TWCSIZ_NATBAS')
      END IF
      IF (.NOT.ALLOCATED(CELHFB)) THEN
          ALLOCATE (CELHFB(1:(8*NDBASE+1)*4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('CELHFB','TWCSIZ_NATBAS')
      END IF
      IF (.NOT.ALLOCATED(EQUNAT)) THEN
          ALLOCATE (EQUNAT(1:8*NDBASE,1:8*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('EQUNAT','TWCSIZ_NATBAS')
      END IF
      IF (.NOT.ALLOCATED(AWATWC)) THEN
          ALLOCATE (AWATWC(1:2*NDBASE,1:2*NDSTAT,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('AWATWC','HFBSIZ')
      END IF
      IF (.NOT.ALLOCATED(BWATWC)) THEN
          ALLOCATE (BWATWC(1:2*NDBASE,1:2*NDSTAT,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('BWATWC','HFBSIZ')
      END IF
C
C=======================================================================
C
      COMHA2(:,:)=-CONJG(COMHAM(:,:))
C
C=======================================================================
C   TRANSFORM H
C=======================================================================
C
      HAMAUX(:,:)=C_ZERO
      HAMCUT(:,:)=C_ZERO
C
      CALL ZGEMM('C','N',NATDIM,2*LDBTOT,2*LDBTOT,
     *                     C_UNIT,EWTAUX,4*NDBASE,
     *                            COMHAM,4*NDBASE,
     *                     C_ZERO,HAMAUX,4*NDBASE)

C
      CALL ZGEMM('N','N',NATDIM,NATDIM,2*LDBTOT,
     *                   C_UNIT,HAMAUX,4*NDBASE,
     *                          EWTAUX,4*NDBASE,
     *                   C_ZERO,HAMCUT,4*NDBASE)
C
C=======================================================================
C   TRANSFORM -H^*+LAMBDA
C=======================================================================
C
      HAMAUX(:,:)=C_ZERO
      HAMCU2(:,:)=C_ZERO
C
      CCWEWT=CONJG(EWTAUX)
C
      CALL ZGEMM('T','N',NATDIM,2*LDBTOT,2*LDBTOT,
     *                C_UNIT,EWTAUX,4*NDBASE,
     *                              COMHA2,4*NDBASE,
     *                       C_ZERO,HAMAUX,4*NDBASE)
C
      CALL ZGEMM('N','N',NATDIM,NATDIM,2*LDBTOT,
     *                C_UNIT,HAMAUX,4*NDBASE,
     *                              CCWEWT,4*NDBASE,
     *                       C_ZERO,HAMCU2,4*NDBASE)
C
C=======================================================================
C   TRANSFORM D
C=======================================================================
C
      COMDEL(:,:)=C_ZERO
      DELCUT(:,:)=C_ZERO
      HAMAUX(:,:)=C_ZERO
C
      DO IBRA=1,2*LDBTOT
         DO IKET=1,2*LDBTOT
C
                IF (IBRA.LE.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
                    COMDEL(IBRA,IKET)=TWDEPP(IBRA,IKET,0,ICHARG)
C
                ELSEIF (IBRA.GT.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
                    COMDEL(IBRA,IKET)=TWDEPM(IBRA-LDBTOT,IKET,1,ICHARG)
C
                ELSEIF (IBRA.LE.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
                    COMDEL(IBRA,IKET)=TWDEPM(IBRA,IKET-LDBTOT,0,ICHARG)
C
                ELSEIF (IBRA.GT.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
                    COMDEL(IBRA,IKET)=
     *              TWDEPP(IBRA-LDBTOT,IKET-LDBTOT,1,ICHARG)
C
                END IF
C
         END DO
      END DO
C
      CALL ZGEMM('C','N',NATDIM,2*LDBTOT,2*LDBTOT,
     *                C_UNIT,EWTAUX,4*NDBASE,
     *                              COMDEL,4*NDBASE,
     *                       C_ZERO,HAMAUX,4*NDBASE)
C
      CALL ZGEMM('N','N',NATDIM,NATDIM,2*LDBTOT,
     *                C_UNIT,HAMAUX,4*NDBASE,
     *                              CCWEWT,4*NDBASE,
     *                       C_ZERO,DELCUT,4*NDBASE)
C
C=======================================================================
C   OBTAINING LOWER DIAGONAL PAIRING TERM
C=======================================================================
C
      DELCU2(:,:)=C_ZERO
C
      DO IBRA=1,NATDIM
         DO IKET=1,NATDIM
C
            DELCU2(IBRA,IKET)=CONJG(DELCUT(IKET,IBRA))
C
         END DO
      END DO
C
C=======================================================================
C   PUTTING ALL TOGETHER IN CELHFB
C=======================================================================
C
      DO IBRA=1,2*NATDIM
         DO IKET=1,2*NATDIM
C
            IF (IBRA.GE.IKET) THEN
C
                NCOUNT=IBRA+((4*NATDIM-IKET)*(IKET-1))/2
C
                IF (IBRA.LE.NATDIM.AND.IKET.LE.NATDIM)
     *              CELHFB(NCOUNT)=HAMCUT(IBRA,IKET)
C
                IF (IBRA.GT.NATDIM.AND.IKET.LE.NATDIM)
     *              CELHFB(NCOUNT)=DELCU2(IBRA-NATDIM,IKET)
C
                IF (IBRA.GT.NATDIM.AND.IKET.GT.NATDIM)
     *              CELHFB(NCOUNT)=HAMCU2(IBRA-NATDIM,IKET-NATDIM)
C
            END IF
C
         END DO
      END DO
C
      DO IBRA=1,2*NATDIM
C
         IKET=IBRA
C
         NCOUNT=IBRA+((4*NATDIM-IKET)*(IKET-1))/2
C
         IF (IBRA.LE.NATDIM)
     *
     *       CELHFB(NCOUNT)=CELHFB(NCOUNT)-EFERMI
C
         IF (IBRA.GT.NATDIM)
     *
     *       CELHFB(NCOUNT)=CELHFB(NCOUNT)+EFERMI
C
      END DO
C
C=======================================================================
C   HERE WE GO!
C=======================================================================
C
      LDQUAP=NATDIM
C
      CALL DIAMAT(CELHFB,EQUTWC,EQUNAT,2*NATDIM,8*NDBASE,LDQUAP)
C
      DO IBASE=1,NATDIM
         DO ISTATE=1,LDQUAP
C
            BWANAT(IBASE,ISTATE)=EQUNAT(IBASE,ISTATE)
            AWANAT(IBASE,ISTATE)=EQUNAT(IBASE+NATDIM,ISTATE)
C
         END DO
      END DO
C
C=======================================================================
C   GOING BACK TO THE HO BASIS
C=======================================================================
C
      AWAAUX(:,:)=C_ZERO
      BWAAUX(:,:)=C_ZERO
C
      CALL ZGEMM('N','N',2*LDBTOT,LDQUAP,NATDIM,
     *                   C_UNIT,EWTAUX,4*NDBASE,
     *                          BWANAT,4*NDBASE,
     *                   C_ZERO,BWAAUX,4*NDBASE)
C
      CALL ZGEMM('N','N',2*LDBTOT,LDQUAP,NATDIM,
     *                   C_UNIT,CCWEWT,4*NDBASE,
     *                          AWANAT,4*NDBASE,
     *                   C_ZERO,AWAAUX,4*NDBASE)
C
C
C=======================================================================
C         DETERMINING THE ORDERED QUASIPARTICLE SPECTRUM
C=======================================================================
C
      IORDER=0
C
      DO ISTATE=1,LDQUAP
C
         IORDER=IORDER+1
         NORDER(IORDER)=IORDER
         EORDER(IORDER)=-EQUTWC(ISTATE)
         NSTATE(IORDER)=ISTATE
C
      END DO
C
      LDORDE=IORDER
C
      CALL ORDER1(EORDER,NORDER,LDORDE)
C
C=======================================================================
C         LOOP OVER ORDERED QUASIPARTICLE STATES
C=======================================================================
C
      INCORR=0
C
      IDSTAT(0)=0
      IDSTAT(1)=0
C
      DO IALLST=1,LDORDE
C
         ISTATE=NSTATE(NORDER(IALLST))
C
C=======================================================================
C        CALCULATING THE NORMS OF THE "FIRST" COMPONENTS
C=======================================================================
C
         RESAUX=0.0D0
C
         DO IBASE=1,NATDIM
C
            RESAUX=RESAUX+CONJG(BWANAT(IBASE,ISTATE))
     *                   *      BWANAT(IBASE,ISTATE)
C
         END DO

         VORDER=REAL(RESAUX)
C
C=======================================================================
C        CHECKING VALUES OF VORDER - JUST IN CASE. VALUES WHICH ARE OUT
C        OF BOUNDS UP TO A SMALL MACHINE-PRECISION NUMBER "EMACHI"  ARE
C        REPLACED BY THE VALUES AT BOUNDS.
C=======================================================================
C
         EMACHI=1.0D-10
C
         IF (VORDER.GT.1+EMACHI.OR.VORDER.LT. -EMACHI) THEN
C
             WRITE(NFIPRI,'(/,1X,20(1H/),
     *             '' ERROR IN NORMALIZATION OF'',      9X,20(1H/),/,
     *             1X,20(1H/),
     *             '' QUASIPARTICLES IN NATBAS.'',      9X,20(1H/),/,
     *             1X,20(1H/),
     *             '' IF THIS IS A ROUND-OFF PROBLEM'', 4X,20(1H/),/,
     *             1X,20(1H/),
     *             '' YOU MAY CHANGE VALUE OF EMACHI.'',3X,20(1H/),/,
     *             1X,20(1H/),'' ISTATE='',I5,''        '',2X,
     *                                           12X,20(1H/),/,
     *             1X,20(1H/),'' VORDER='',D25.18,
     *                                            2X,20(1H/),/)')
     *
     *             ISTATE,VORDER
C
             STOP 'ERROR IN NORMALIZATION OF QUASIPARTICLES IN HFBSIZ'
C
         END IF
C
         IF (VORDER.GT.1) VORDER=1
         IF (VORDER.LT.0) VORDER=0
C
C=======================================================================
C
         HORDER=EORDER(IALLST)*(1-2*VORDER)
C
                          CORDER=HORDER+EFERMI
         IF (LAMCUT.EQ.1) CORDER=HORDER
C
         IF (CORDER.GT.ECUTOF) GO TO 1
C
C=======================================================================
C         WE PUT ALL THE QUASI-PARTICLES IN THE SAME BLOCK
C=======================================================================
C
          IREVER=0
C
C=======================================================================
C
         IDSTAT(IREVER)=IDSTAT(IREVER)+1
C
         IF (IDSTAT(IREVER).GT.2*NDSTAT) THEN
C
             INCORR=1
             GO TO 1
C
         END IF
C
         KSTATE(IDSTAT(IREVER),IREVER)=ISTATE
C
C=======================================================================
C
         DEQUIV(IDSTAT(IREVER),IREVER)=
     *                            SQRT(ABS(EORDER(IALLST)**2-HORDER**2))
C
         EEQUIV(IDSTAT(IREVER),IREVER)=HORDER+EFERMI
C
C=======================================================================
C
         ENQUAP(IDSTAT(IREVER),IREVER)=EORDER(IALLST)
C
         V2QUAP(IDSTAT(IREVER),IREVER)=VORDER
C

         QPSIMP(IDSTAT(IREVER),IREVER)=0.0D0
C
    1    CONTINUE
C
      END DO
C
C=======================================================================
C         DEFINING NUMBERS OF QUASIPARTICLES KEPT AFTER THE CUT-OFF
C=======================================================================
C
C      DO IREVER=0,NDREVE
C
         NUMBQP(0,ICHARG)=IDSTAT(0)
C
C      END DO
C
C=======================================================================
C         STOP IN CASE OF AN INCORRECT PHASE SPACE
C=======================================================================
C
      IF (INCORR.EQ.1) THEN
C
          WRITE(NFIPRI,'(/,                              1X,20(1H/),
     *              '' TOO LARGE PHASE-SPACE SIZE IN NATBAS'',
     *                                                   1X,20(1H/),/,
     *                 1X,20(1H/),
     *              '' ICHARG='',I2,3X,''2NDSTAT='',I4, 13X,20(1H/),/,
     *                 1X,20(1H/),
     *              '' IREVER='',I2,3X,'' REQUIRES NUMBQP='',I4,
     *                                                   4X,20(1H/),/,
     *                 1X,20(1H/),
     *              '' IREVER='',I2,3X,'' REQUIRES NUMBQP='',I4,
     *                                                   4X,20(1H/),/,
     *                 1X,20(1H/),
     *              '' INCREASE NDSTAT UP TO NUMBQP OR MORE'',
     *                                                   1X,20(1H/),/)')
     *
     *    ICHARG,2*NDSTAT,(IREVER,NUMBQP(IREVER,ICHARG),IREVER=0,NDREVE)
C
          STOP ' TOO LARGE PHASE-SPACE SIZE IN HFBSIZ'
C
      END IF
C
      BWATWC=C_ZERO
      AWATWC=C_ZERO
C
      DO IREVER=0,NDREVE
         DO ISTATE=1,NUMBQP(IREVER,ICHARG)
C
            JSTATE=KSTATE(ISTATE,IREVER)
C
            LSTATE=ISTATE+IREVER*NUMBQP(0,ICHARG)
C
            DO IBASE=1,LDBTOT
C
               BWATWC(IBASE,LSTATE,0)=BWAAUX(IBASE       ,JSTATE)
               BWATWC(IBASE,LSTATE,1)=BWAAUX(IBASE+LDBTOT,JSTATE)
C
               AWATWC(IBASE,LSTATE,0)=AWAAUX(IBASE       ,JSTATE)
               AWATWC(IBASE,LSTATE,1)=AWAAUX(IBASE+LDBTOT,JSTATE)
C
            END DO
C
         END DO
      END DO
C
      IREMQB(:)=0
      INUMQB(:)=0
      JNUMQB(:)=0
C
      DO NUBLOC=1,NDBLOC
C
C               WRITE (*,'(A8,3I3,3X,4I4,/,20X,4I4)') 'HFBSIZ-A',
C    *                          NUBLOC,0,ICHARG,
C    *                          INSIZB(NUBLOC),IDSIZB(NUBLOC)
C
          IF (IABS(IDSIZB(NUBLOC)).EQ.1) THEN
C
             CALL BLOSIZ(ICHARG,INSIZB(NUBLOC),IDSIZB(NUBLOC),
     *                          INUMQB(NUBLOC),IREMQB(NUBLOC),
     *                                  IFIBLB,ISABLB,NUBLOC,2)
C
             DO MUBLOC=1,NUBLOC-1
C
                IF (INUMQB(MUBLOC).EQ.INUMQB(NUBLOC).AND.
     *              IREMQB(MUBLOC).EQ.IREMQB(NUBLOC))
     *
     *          STOP ' QUASIPARTICLES BLOCKED TWICE IN TWCSIZ_NATBAS'
C
             END DO
C
          END IF
C
      END DO
C
      DO NUBLOC=1,NDBLOC
C
          IF (IABS(IDSIZB(NUBLOC)).EQ.1) THEN
C
             JSTATE=KSTATE(INUMQB(NUBLOC),IREMQB(NUBLOC))
C
             LSTATE=INUMQB(NUBLOC)+IREMQB(NUBLOC)*NUMBQP(0,ICHARG)
C
C=======================================================================
C         HERE THE INVERTED QUASIPARTICLE REPLACES THE STANDARD ONE
C=======================================================================
C
             DO IBASE=1,LDBTOT
C
                BWATWC(IBASE,LSTATE,0)=
     *                              CONJG(AWAAUX(IBASE       ,JSTATE))
                BWATWC(IBASE,LSTATE,1)=
     *                              CONJG(AWAAUX(IBASE+LDBTOT,JSTATE))
C
                AWATWC(IBASE,LSTATE,0)=
     *                              CONJG(BWAAUX(IBASE       ,JSTATE))
                AWATWC(IBASE,LSTATE,1)=
     *                              CONJG(BWAAUX(IBASE+LDBTOT,JSTATE))
C
             END DO
C
             V2QUAP(INUMQB(NUBLOC),IREMQB(NUBLOC))=1.0D0-
     *       V2QUAP(INUMQB(NUBLOC),IREMQB(NUBLOC))
C
             JNUMQB(NUBLOC)=INUMQB(NUBLOC)
C
          END IF
C
      END DO
C
C=======================================================================
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_BIGDEN(ICHARG)
C=======================================================================
      USE hfodd_sizes
      USE TWCEPP
      USE TWCEPM
      USE BIGDEN
C=======================================================================
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C=======================================================================
C     HERE WE BUILD THE DENSITY MATRIX IN TWO-CENTRE HARMONIC
C     OSCILLATOR BASIS. THAT WILL BE USED TO:
C     1. OBTAIN EXPECTATION VALUES OF ONE-BODY OPERATORS
C     2. OTHER QUANTUM STUFF
C=======================================================================
C
      IALLOC=0
      IF (.NOT.ALLOCATED(BIGTWC)) THEN
          ALLOCATE (BIGTWC(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('BIGTWC','TWC_BIGDEN')
      END IF
C=======================================================================
C
      DO IBRA=1,2*LDBTOT
         DO IKET=1,2*LDBTOT
C
            IF (IBRA.LE.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
                BIGTWC(IBRA,IKET)= TWC_PP(IBRA,IKET,0,ICHARG)
C
            ELSEIF (IBRA.GT.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
                BIGTWC(IBRA,IKET)=TWC_PM(IBRA-LDBTOT,IKET,1,ICHARG)
C
            ELSEIF (IBRA.LE.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
                BIGTWC(IBRA,IKET)=TWC_PM(IBRA,IKET-LDBTOT,0,ICHARG)
C
            ELSEIF (IBRA.GT.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
                BIGTWC(IBRA,IKET)=
     *             TWC_PP(IBRA-LDBTOT,IKET-LDBTOT,1,ICHARG)
C
            END IF
C
         END DO
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_INTMUL(ELEMUL,LAMBDA,MIU,IMAREA,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          ELEMUL(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /DATYLM/ COEYLM(0:NDMULT,0:NDMULT),
     *
     *                NREYLM(0:NDMULT,0:NDMULT),
     *                KREYLM(0:NDMULT,0:NDMULT,1:NDTERM,0:NDKART),
     *
     *                NIMYLM(0:NDMULT,1:NDMULT),
     *                KIMYLM(0:NDMULT,1:NDMULT,1:NDTERM,0:NDKART)
      COMMON
     *       /OURUNI/ QUNITS(0:NDMULT,0:NDMULT)
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *     /COMU2C/ CMUTWC(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART,1:NDTWBL)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      NIU=IABS(MIU)
C
      IF (MIU.GE.0) THEN
          JPHASE=IMAREA
      ELSE
          JPHASE=(-1)**MIU
      END IF
C
      ELEMUL(:,:)=0.0D0
C
      IF (IMAREA.EQ.1) THEN
C
C=======================================================================
C         CALCULATING THE REAL PART
C=======================================================================
C
          DO NUTERM=1,NREYLM(LAMBDA,NIU)
C
             DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
                NX_KET=NXVECT(IKET,IND2HR(LTWCEN))
                NY_KET=NYVECT(IKET,IND2HR(LTWCEN))
                NZ_KET=NZVECT(IKET,IND2HR(LTWCEN))
C
                DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                   NX_BRA=NXVECT(IBRA,IND2HL(LTWCEN))
                   NY_BRA=NYVECT(IBRA,IND2HL(LTWCEN))
                   NZ_BRA=NZVECT(IBRA,IND2HL(LTWCEN))
C
                   ELEMUL(IBRA,IKET) =
     *             ELEMUL(IBRA,IKET) + KREYLM(LAMBDA,NIU,NUTERM,0)
     *
     *      * CMUTWC(KREYLM(LAMBDA,NIU,NUTERM,1),NX_BRA,NX_KET,1,LTWCEN)
     *      * CMUTWC(KREYLM(LAMBDA,NIU,NUTERM,2),NY_BRA,NY_KET,2,LTWCEN)
     *      * CMUTWC(KREYLM(LAMBDA,NIU,NUTERM,3),NZ_BRA,NZ_KET,3,LTWCEN)
     *      * JPHASE
C
                END DO
             END DO
C
          END DO
C
      ELSE
C
          IF (MIU.EQ.0) STOP 'NO IMAGINARY PART FOR MIU=0 IN TWC_INTMUL'
C
C=======================================================================
C         CALCULATING THE IMAGINARY PART
C=======================================================================
C
          DO NUTERM=1,NIMYLM(LAMBDA,NIU)
C
             DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
                NX_KET=NXVECT(IKET,IND2HR(LTWCEN))
                NY_KET=NYVECT(IKET,IND2HR(LTWCEN))
                NZ_KET=NZVECT(IKET,IND2HR(LTWCEN))
C
                DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                   NX_BRA=NXVECT(IBRA,IND2HL(LTWCEN))
                   NY_BRA=NYVECT(IBRA,IND2HL(LTWCEN))
                   NZ_BRA=NZVECT(IBRA,IND2HL(LTWCEN))
C
                   ELEMUL(IBRA,IKET) =
     *             ELEMUL(IBRA,IKET) + KIMYLM(LAMBDA,NIU,NUTERM,0)
     *
     *      * CMUTWC(KIMYLM(LAMBDA,NIU,NUTERM,1),NX_BRA,NX_KET,1,LTWCEN)
     *      * CMUTWC(KIMYLM(LAMBDA,NIU,NUTERM,2),NY_BRA,NY_KET,2,LTWCEN)
     *      * CMUTWC(KIMYLM(LAMBDA,NIU,NUTERM,3),NZ_BRA,NZ_KET,3,LTWCEN)
     *      * JPHASE
C
                END DO
             END DO
C
          END DO
C
      END IF
C
C=======================================================================
C        NORMALIZATION FACTOR AND PHASE
C=======================================================================
C
      DO IKET=1,LDBATW(IND2HR(LTWCEN))
         NY_KET=NYVECT(IKET,IND2HR(LTWCEN))
         DO IBRA=1,LDBATW(IND2HL(LTWCEN))
            NY_BRA=NYVECT(IBRA,IND2HL(LTWCEN))
C
            ELEMUL(IBRA,IKET) = ELEMUL(IBRA,IKET)
     *                        * QUNITS(LAMBDA,NIU)
     *                        * COEYLM(LAMBDA,NIU)
C
         END DO
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_MULFRA(ELEMUL,LAMBDA,MIU,IMAREA,LTWCEN,ICENTE)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          ELEMUL(1:NDBASE,1:NDBASE)
C
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /DATYLM/ COEYLM(0:NDMULT,0:NDMULT),
     *
     *                NREYLM(0:NDMULT,0:NDMULT),
     *                KREYLM(0:NDMULT,0:NDMULT,1:NDTERM,0:NDKART),
     *
     *                NIMYLM(0:NDMULT,1:NDMULT),
     *                KIMYLM(0:NDMULT,1:NDMULT,1:NDTERM,0:NDKART)
      COMMON
     *       /OURUNI/ QUNITS(0:NDMULT,0:NDMULT)
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *     /COMUSH/ CMUSH1(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                1:NDTWBL),
     *              CMUSH2(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART,
     *                                                1:NDTWBL)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      NIU=IABS(MIU)
C
      IF (MIU.GE.0) THEN
          JPHASE=IMAREA
      ELSE
          JPHASE=(-1)**MIU
      END IF
C
      ELEMUL(:,:)=0.0D0
C
      IF (IMAREA.EQ.1) THEN
C
C=======================================================================
C         CALCULATING THE REAL PART
C=======================================================================
C
          DO NUTERM=1,NREYLM(LAMBDA,NIU)
C
             DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
                NX_KET=NXVECT(IKET,IND2HR(LTWCEN))
                NY_KET=NYVECT(IKET,IND2HR(LTWCEN))
                NZ_KET=NZVECT(IKET,IND2HR(LTWCEN))
C
                DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                   NX_BRA=NXVECT(IBRA,IND2HL(LTWCEN))
                   NY_BRA=NYVECT(IBRA,IND2HL(LTWCEN))
                   NZ_BRA=NZVECT(IBRA,IND2HL(LTWCEN))
C
                   IF (ICENTE.EQ.1) THEN
C
                   ELEMUL(IBRA,IKET) =
     *             ELEMUL(IBRA,IKET) + KREYLM(LAMBDA,NIU,NUTERM,0)
     *
     *      * CMUSH1(KREYLM(LAMBDA,NIU,NUTERM,1),NX_BRA,NX_KET,1,LTWCEN)
     *      * CMUSH1(KREYLM(LAMBDA,NIU,NUTERM,2),NY_BRA,NY_KET,2,LTWCEN)
     *      * CMUSH1(KREYLM(LAMBDA,NIU,NUTERM,3),NZ_BRA,NZ_KET,3,LTWCEN)
     *      * JPHASE
C
                   ELSEIF (ICENTE.EQ.2) THEN
C
                   ELEMUL(IBRA,IKET) =
     *             ELEMUL(IBRA,IKET) + KREYLM(LAMBDA,NIU,NUTERM,0)
     *
     *      * CMUSH2(KREYLM(LAMBDA,NIU,NUTERM,1),NX_BRA,NX_KET,1,LTWCEN)
     *      * CMUSH2(KREYLM(LAMBDA,NIU,NUTERM,2),NY_BRA,NY_KET,2,LTWCEN)
     *      * CMUSH2(KREYLM(LAMBDA,NIU,NUTERM,3),NZ_BRA,NZ_KET,3,LTWCEN)
     *      * JPHASE
C
                   END IF
C
                END DO
             END DO
C
          END DO
C
      ELSE
C
          IF (MIU.EQ.0) STOP 'NO IMAGINARY PART FOR MIU=0 IN TWC_INTMUL'
C
C=======================================================================
C         CALCULATING THE IMAGINARY PART
C=======================================================================
C
          DO NUTERM=1,NIMYLM(LAMBDA,NIU)
C
             DO IKET=1,LDBATW(IND2HR(LTWCEN))
C
                NX_KET=NXVECT(IKET,IND2HR(LTWCEN))
                NY_KET=NYVECT(IKET,IND2HR(LTWCEN))
                NZ_KET=NZVECT(IKET,IND2HR(LTWCEN))
C
                DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                   NX_BRA=NXVECT(IBRA,IND2HL(LTWCEN))
                   NY_BRA=NYVECT(IBRA,IND2HL(LTWCEN))
                   NZ_BRA=NZVECT(IBRA,IND2HL(LTWCEN))
C
                   IF (ICENTE.EQ.1) THEN
C
                   ELEMUL(IBRA,IKET) =
     *             ELEMUL(IBRA,IKET) + KIMYLM(LAMBDA,NIU,NUTERM,0)
     *      * CMUSH1(KIMYLM(LAMBDA,NIU,NUTERM,1),NX_BRA,NX_KET,1,LTWCEN)
     *      * CMUSH1(KIMYLM(LAMBDA,NIU,NUTERM,2),NY_BRA,NY_KET,2,LTWCEN)
     *      * CMUSH1(KIMYLM(LAMBDA,NIU,NUTERM,3),NZ_BRA,NZ_KET,3,LTWCEN)
     *      * JPHASE
                   ELSEIF (ICENTE.EQ.2) THEN
C
                   ELEMUL(IBRA,IKET) =
     *             ELEMUL(IBRA,IKET) + KIMYLM(LAMBDA,NIU,NUTERM,0)
     *      * CMUSH2(KIMYLM(LAMBDA,NIU,NUTERM,1),NX_BRA,NX_KET,1,LTWCEN)
     *      * CMUSH2(KIMYLM(LAMBDA,NIU,NUTERM,2),NY_BRA,NY_KET,2,LTWCEN)
     *      * CMUSH2(KIMYLM(LAMBDA,NIU,NUTERM,3),NZ_BRA,NZ_KET,3,LTWCEN)
     *      * JPHASE
C
                   END IF
C
                END DO
             END DO
C
          END DO
C
      END IF
C
C=======================================================================
C        NORMALIZATION FACTOR AND PHASE
C=======================================================================
C
      DO IKET=1,LDBATW(IND2HR(LTWCEN))
         NY_KET=NYVECT(IKET,IND2HR(LTWCEN))
         DO IBRA=1,LDBATW(IND2HL(LTWCEN))
            NY_BRA=NYVECT(IBRA,IND2HL(LTWCEN))
C
            ELEMUL(IBRA,IKET) = ELEMUL(IBRA,IKET)
     *                        * QUNITS(LAMBDA,NIU)
     *                        * COEYLM(LAMBDA,NIU)
C
         END DO
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_BIGSIM(AUXILI,LTWCEN,IMAREA,KARTEZ)
C=======================================================================
      USE hfodd_sizes
      USE BIGOPP
      USE BIGOPM
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
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      DIMENSION
     *          AUXILI(1:NDBASE,1:NDBASE)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C=======================================================================
      IALLOC=0
      IF (.NOT.ALLOCATED(TWCOPP)) THEN
          ALLOCATE (TWCOPP(1:2*NDBASE,1:2*NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWCOPP','TWC_BIGSIM')
      END IF
      IF (.NOT.ALLOCATED(TWCOPM)) THEN
          ALLOCATE (TWCOPM(1:2*NDBASE,1:2*NDBASE,0:NDREVE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWCOPM','TWC_BIGSIM')
      END IF
C=======================================================================
C    MATRIX ELEMENTS OF ONE-BODY OPERATORS FOR DIFFERENT CENTERS ARE PUT
C    IN DIFFERENT SIMPLEX ARRAYS, IN THE STYLE OF TWCCOP.
C=======================================================================

      IF (LTWCEN.EQ.1) THEN
          MBRA=0
          MKET=0
C
      ELSEIF (LTWCEN.EQ.2) THEN
C
          MBRA=LDBATW(1)
          MKET=LDBATW(1)
C
      ELSEIF (LTWCEN.EQ.3) THEN
C
          MBRA=LDBATW(1)
          MKET=0
C
      END IF
C
      DO IKET=1,LDBATW(IND2HR(LTWCEN))
         DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
	        NY=NYVECT(IBRA,IND2HL(LTWCEN))
	        MY=NYVECT(IKET,IND2HR(LTWCEN))
C
              TWCOPP(MBRA+IBRA,MKET+IKET,0)=
     *         AUXILI(IBRA,IKET)*IPHAPP(NY,MY,KARTEZ)
              TWCOPP(MBRA+IBRA,MKET+IKET,1)=
     *         AUXILI(IBRA,IKET)*IPHAMM(NY,MY,KARTEZ)
              TWCOPM(MBRA+IBRA,MKET+IKET,0)=
     *         AUXILI(IBRA,IKET)*IPHAPM(NY,MY,KARTEZ)
              TWCOPM(MBRA+IBRA,MKET+IKET,1)=
     *         AUXILI(IBRA,IKET)*IPHAMP(NY,MY,KARTEZ)
C
         END DO
      END DO
C
      IF (LTWCEN.EQ.3) THEN
C
          MBRA=0
          MKET=LDBATW(1)
C
          IF (IMAREA.EQ.1) THEN
C
          DO IKET=1,LDBATW(IND2HL(LTWCEN))
             DO IBRA=1,LDBATW(IND2HR(LTWCEN))
C
             TWCOPP(MBRA+IBRA,MKET+IKET,0)=TWCOPP(MKET+IKET,MBRA+IBRA,0)
             TWCOPP(MBRA+IBRA,MKET+IKET,1)=TWCOPP(MKET+IKET,MBRA+IBRA,1)
             TWCOPM(MBRA+IBRA,MKET+IKET,0)=TWCOPM(MKET+IKET,MBRA+IBRA,1)
             TWCOPM(MBRA+IBRA,MKET+IKET,1)=TWCOPM(MKET+IKET,MBRA+IBRA,0)
C
             END DO
          END DO
C
          ELSEIF (IMAREA.EQ.-1) THEN
C
          DO IKET=1,LDBATW(IND2HL(LTWCEN))
             DO IBRA=1,LDBATW(IND2HR(LTWCEN))
C
            TWCOPP(MBRA+IBRA,MKET+IKET,0)=-TWCOPP(MKET+IKET,MBRA+IBRA,0)
            TWCOPP(MBRA+IBRA,MKET+IKET,1)=-TWCOPP(MKET+IKET,MBRA+IBRA,1)
            TWCOPM(MBRA+IBRA,MKET+IKET,0)=-TWCOPM(MKET+IKET,MBRA+IBRA,1)
            TWCOPM(MBRA+IBRA,MKET+IKET,1)=-TWCOPM(MKET+IKET,MBRA+IBRA,0)
C
             END DO
          END DO
C
          ELSE
C
              STOP ' IMAREA WRONG IN TWC_BIGSIM'
C
          END IF
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
      SUBROUTINE TWC_BIGOBO
C=======================================================================
      USE hfodd_sizes
      USE TWCOBO
      USE BIGOPP
      USE BIGOPM
C=======================================================================
      COMPLEX
     *          C_ZERO
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C=======================================================================
      C_ZERO=CMPLX(0.0D0,0.0D0)
C=======================================================================
      IALLOC=0
      IF (.NOT.ALLOCATED(BIGOBO)) THEN
          ALLOCATE (BIGOBO(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('BIGOBO','TWC_BIGOBO')
      END IF
C=======================================================================
C     THE ARRAYS OBTAINED IN BIGSIM ARE COMBINED IN BIGOBO, THIS IS TO
C     COMPUTE THE TRACE EXPLICITLY WITH THE DENSITY MATRIX.
C=======================================================================
C
      BIGOBO=C_ZERO
C
      DO IBRA=1,2*LDBTOT
         DO IKET=1,2*LDBTOT
C
            IF (IBRA.LE.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
               BIGOBO(IBRA,IKET)=TWCOPP(IBRA,IKET,0)
C
            ELSEIF (IBRA.GT.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
               BIGOBO(IBRA,IKET)=TWCOPP(IBRA-LDBTOT,IKET-LDBTOT,1)
C
            ELSEIF (IBRA.LE.LDBTOT.AND.IKET.GT.LDBTOT) THEN
C
               BIGOBO(IBRA,IKET)=TWCOPM(IBRA,IKET-LDBTOT,0)
C
            ELSEIF (IBRA.GT.LDBTOT.AND.IKET.LE.LDBTOT) THEN
C
               BIGOBO(IBRA,IKET)=TWCOPM(IBRA-LDBTOT,IKET,1)
C
            END IF
C
         END DO
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_INTCON(ICHARG,NMUCON,HAMCON,LTWCEN,IF_RPA)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      ALLOCATABLE ELEMUL(:,:)
C=======================================================================
      DIMENSION
     *          HAMCON(1:NDBASE,1:NDBASE)
      DIMENSION
     *          CMUL_T(0:NDMULT,-NDMULT:NDMULT)
      DIMENSION
     *          CMUL_V(0:NDMULT,-NDMULT:NDMULT),
     *          QMUL_V(0:NDMULT,-NDMULT:NDMULT)
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
     *       /QMULTI/ QMUL_N(0:NDMULT,-NDMULT:NDMULT),
     *                QMUL_P(0:NDMULT,-NDMULT:NDMULT),
     *                QMUL_T(0:NDMULT,-NDMULT:NDMULT)
      COMMON
     *       /QCNSTR/ STIFFQ(0:NDMULT,-NDMULT:NDMULT),
     *                QASKED(0:NDMULT,-NDMULT:NDMULT),
     *                IFLAGQ(0:NDMULT,-NDMULT:NDMULT)
      COMMON
     *       /VCNSTR/ STIFFV(0:NDMULT,-NDMULT:NDMULT),
     *                VASKED(0:NDMULT,-NDMULT:NDMULT),
     *                IFLAGV(0:NDMULT,-NDMULT:NDMULT)
      COMMON
     *       /GCNSTR/ STIFFG(0:NDMULT),
     *                GASKED(0:NDMULT),
     *                IFLAGG(0:NDMULT)
      COMMON
     *       /QLASTR/ GALMUQ(0:NDMULT,-NDMULT:NDMULT),
     *                QLINEA(0:NDMULT,-NDMULT:NDMULT),
     *                IFLALQ(0:NDMULT,-NDMULT:NDMULT)
      COMMON
     *       /VLASTR/ GALMUV(0:NDMULT,-NDMULT:NDMULT),
     *                VLINEA(0:NDMULT,-NDMULT:NDMULT),
     *                IFLALV(0:NDMULT,-NDMULT:NDMULT)
      COMMON
     *       /OURUNI/ QUNITS(0:NDMULT,0:NDMULT)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)

C
C=======================================================================
C
      IALLOC=0
C
C=======================================================================
      ALLOCATE (ELEMUL(1:NDBASE,1:NDBASE),STAT=IALLOC)
      IF (IALLOC.NE.0) CALL NOALLO('ELEMUL','TWC_INTCON')
C
C=======================================================================
C
      QMUL_V(:,:)=QMUL_N(:,:)-QMUL_P(:,:)
C
      HAMCON(:,:)=0.0D0
      ELEMUL(:,:)=0.0D0
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS THAT IN THE CONSTRAINT TERMS
C        MULTIPLY THE MULTIPOLE MOMENTS
C=======================================================================
C
      DO LAMBDA=0,NMUCON
C
         RESULT=0.0D0
C
         DO MIU=-LAMBDA,LAMBDA
            NIU=IABS(MIU)
            RESULT=RESULT+(QMUL_T(LAMBDA,MIU)/QUNITS(LAMBDA,NIU))**2
         END DO
C
         GMUL_T=SQRT(RESULT)*QUNITS(LAMBDA,0)
C
         DO MIU=-LAMBDA,LAMBDA
C
            NIU=IABS(MIU)
C
            CMUL_T(LAMBDA,MIU)=0.0D0
C
            IF (IFLAGQ(LAMBDA,MIU).EQ.1) THEN
C
                IF (IF_RPA.NE.1)
     *              CMUL_T(LAMBDA,MIU)=CMUL_T(LAMBDA,MIU)
     *                +2.0D0*(QMUL_T(LAMBDA,MIU)-QASKED(LAMBDA,MIU))
     *                      * STIFFQ(LAMBDA,MIU)
C
                IF (ABS(IFLALQ(LAMBDA,MIU)).EQ.1)
     *
     *              CMUL_T(LAMBDA,MIU)=CMUL_T(LAMBDA,MIU)
     *                                -GALMUQ(LAMBDA,MIU)
C
            END IF
C
            IF (IFLAGG(LAMBDA).EQ.1.AND.GMUL_T.GT.0)
     *
     *          CMUL_T(LAMBDA,MIU)=CMUL_T(LAMBDA,MIU)
     *                +2.0D0*(GMUL_T            -GASKED(LAMBDA))
     *                      * STIFFG(LAMBDA)
     *                      * QMUL_T(LAMBDA,MIU)/GMUL_T
     *                      * QUNITS(LAMBDA, 0 )**2
     *                      / QUNITS(LAMBDA,NIU)**2
C
            CMUL_V(LAMBDA,MIU)=0.0D0
C
            IF (IFLAGV(LAMBDA,MIU).EQ.1) THEN
C
                IF (IF_RPA.NE.1)
     *              CMUL_V(LAMBDA,MIU)=CMUL_V(LAMBDA,MIU)
     *                +2.0D0*(QMUL_V(LAMBDA,MIU)-VASKED(LAMBDA,MIU))
     *                      * STIFFV(LAMBDA,MIU)
C
                IF (ABS(IFLALV(LAMBDA,MIU)).EQ.1)
     *
     *              CMUL_V(LAMBDA,MIU)=CMUL_V(LAMBDA,MIU)
     *                                -GALMUV(LAMBDA,MIU)
C
            END IF
C
         END DO
C
      END DO
C
C=======================================================================
C        HERE WE CALCULATE CONTRIBUTIONS FROM REAL AND IMAGINARY
C        PARTS OF THE MULTIPOLE MOMENTS
C=======================================================================
C
      DO LAMBDA=0,NMUCON
         DO MIU=-LAMBDA,LAMBDA
C
            IF (IFLAGQ(LAMBDA,MIU).EQ.1.OR.
     *          IFLAGV(LAMBDA,MIU).EQ.1.OR.IFLAGG(LAMBDA).EQ.1) THEN
C
                IF (MIU.LT.0) IMAREA=-1
                IF (MIU.GE.0) IMAREA= 1
C
                CALL TWC_INTMUL(ELEMUL,LAMBDA,MIU,IMAREA,LTWCEN)
C
                DO IKET=1,LDBATW(IND2HR(LTWCEN))
                   DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                      IF (IFLAGQ(LAMBDA,MIU).EQ.1)
     *                    HAMCON(IBRA,IKET)=HAMCON(IBRA,IKET)
     *                                     +ELEMUL(IBRA,IKET)
     *                                     *CMUL_T(LAMBDA,MIU)
C
                      IF (IFLAGV(LAMBDA,MIU).EQ.1) THEN
C
                          IF (ICHARG.EQ.0)
     *                    HAMCON(IBRA,IKET)=HAMCON(IBRA,IKET)
     *                                     +ELEMUL(IBRA,IKET)
     *                                     *CMUL_V(LAMBDA,MIU)
                          IF (ICHARG.EQ.1)
     *                    HAMCON(IBRA,IKET)=HAMCON(IBRA,IKET)
     *                                     -ELEMUL(IBRA,IKET)
     *                                     *CMUL_V(LAMBDA,MIU)
                      END IF
C
                   END DO
                END DO
C
             END IF
C
         END DO
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_SPAVER(AVRVAL,ICHARG,IMAREA,IHERMI)
C=======================================================================
      USE ALLWAV
      USE TWCWAV
      USE TWCOBO
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          AVRVAL(1:2*NDSTAT)
C
      COMPLEX
     *          RESULT,RESUL0
      COMPLEX
     *          WAVTWC
C=======================================================================
      COMMON
     *       /SPNUMS/ NUMBSP(0:NDREVE,0:NDISOS)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C=======================================================================
      ALLOCATABLE WAVTWC(:,:)
C=======================================================================
C     THIS SUBROUTINE CALCULATES SINGLE-PARTICLE AVERAGE VALUES OF
C     AN OPERATOR IN THE TWO-CENTRE CASE WHICH IS:
C        1. REAL OR IMAGINARY IN THE SIMPLEX BASIS (IMAREA=+1/-1)
C        2. TIME-EVEN OR TIME-ODD (ITIPAR=+1/-1)
C        3. HERMITIAN OR ANTI-HERMITIAN (IHERMI=+1/-1)
C     INPUT ARRAY BIGOBO SHOULD CONTAIN
C        1. REAL OR IMAGINARY PARTS OF MATRIX ELEMENTS (IMAREA=+1/-1)
C     RESULTING AVRVAL VECTOR THEN CONTAINS
C        1. REAL OR IMAGINARY PARTS OF AVERAGE VALUES (IHERMI=+1/-1)
C=======================================================================
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

      LSTATE=NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
      DO ISTATE=1,LSTATE
C
         DO IBASE=1,LDBTOT
            WAVTWC(IBASE       ,ISTATE)=TWCOCC(IBASE,ISTATE,0,ICHARG)
            WAVTWC(IBASE+LDBTOT,ISTATE)=TWCOCC(IBASE,ISTATE,1,ICHARG)
         END DO
      END DO

      DO ISTATE=1,LSTATE
C
         RESULT=CMPLX(0.0D0,0.0D0)
C
            DO IBRA=1,2*LDBTOT
C
               RESUL0=CMPLX(0.0D0,0.0D0)
C
               DO IKET=1,2*LDBTOT
                  RESUL0=RESUL0+BIGOBO(IBRA,IKET)*WAVTWC(IKET,ISTATE)
               END DO
C
               RESULT=RESULT+RESUL0*CONJG(WAVTWC(IBRA,ISTATE))
C
            END DO
C
         IF (IMAREA.EQ.-1) RESULT=RESULT*CMPLX(0.0D0,1.0D0)
C
         IF (IHERMI.EQ.1) THEN
             AVRVAL(ISTATE)=REAL(RESULT)
         ELSE
             AVRVAL(ISTATE)=AIMAG(RESULT)
         END IF
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_TRUHER(NOSACT,KARTEZ,NGAUSS)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          PHERMI(ND2MAX+1),DHERMI(ND2MAX+1)
      DIMENSION
     *          XHERMI(NDGAUS),WHERMI(NDGAUS)
C
      COMMON
     *       /FACHER/ HERFAC(0:ND2MAX)
      COMMON
     *       /CFIPRI/ NFIPRI
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
      COMMON
     *       /SCALNG/ HOMSCA(1:NDKART,1:NDTWCE)
      COMMON
     *       /CENPOS/ CENCOR(1:NDKART,1:NDTWCE)
C
      COMMON
     *       /TRUHER/ XZER_I(1:NDGAUS,1:NDKART,1:NDTWBL),
     *                XZER_J(1:NDGAUS,1:NDKART,1:NDTWBL)
      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:NDTWDD),
     *                DHELE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHLE4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                HERRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DHERI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD),
     *                DDHRI4(0:ND2MAX,1:NDGAUS,1:NDKART,1:NDTWDD)
C
C=======================================================================
C     THIS SUBROUTINE COMPUTES THE HERMITE POLYNOMIALS IN THE TRUE
C     SPATIAL COORDINATES, INSTEAD OF USING THE COMMON INTEGRATION
C     LATTICE. THIS ALLOWS TO SAVE IN THE REVIEW FILE THE PROPER
C     DENSITIES IN SPACE TO BE PLOTTED.
C
C     MUST BE CALLED AT THE END OF THE CALCULATION AS IT WILL DESTROY
C     THE INFORMATION SAVED IN THE DENSITY ARRAYS.
C=======================================================================
C
      HERLE4(:,:,KARTEZ,:)=0.0D0
      DHELE4(:,:,KARTEZ,:)=0.0D0
      DDHLE4(:,:,KARTEZ,:)=0.0D0
      HERRI4(:,:,KARTEZ,:)=0.0D0
      DHERI4(:,:,KARTEZ,:)=0.0D0
      DDHRI4(:,:,KARTEZ,:)=0.0D0
C
      NMAIN2=2*NOSACT+2
C
      IF (NGAUSS.GT.NDGAUS) THEN
C
         WRITE(NFIPRI,'(//,1X,''NGAUSS='',I2,'' .GT. NDGAUS='',I2,
     *                   ''   STOP IN TWC_DEFINT'',//)') NGAUSS,NDGAUS
         STOP ' NGAUSS.GT.NDGAUS  IN  TWC_DEFINT'
C
      END IF
C
      CALL HERMIT(NGAUSS,XHERMI,WHERMI,NDGAUS)
C
      DO INDBLO=1,NDTWBL
C
         HOSCAI=HOMSCA(KARTEZ,IND2HL(INDBLO))
         HOSCAJ=HOMSCA(KARTEZ,IND2HR(INDBLO))
C
         HOSHII=CENCOR(KARTEZ,IND2HL(INDBLO))
         HOSHIJ=CENCOR(KARTEZ,IND2HR(INDBLO))
C
C=======================================================================
C     LEFT WAVE FUNCTION
C=======================================================================
C
         DO IZEROS = 1,NGAUSS
C
            XZER = HOSCAI*(XHERMI(IZEROS)-HOSHII)
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(HOSCAI)
     *                                   / HERFAC( NOSCIL)
                DHELE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                             * SQRT(HOSCAI) *HOSCAI
     *                                   / HERFAC( NOSCIL )
C
                DDHLE4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                      * SQRT(HOSCAI) *HOSCAI*HOSCAI
     *                                   / HERFAC( NOSCIL )
C
            END DO
C
C=======================================================================
C      SAVING THE POINTS TO COMPUTE THE EXPONENTIAL PART OF THE LEFT WF
C=======================================================================
C
            XZER_I(IZEROS,KARTEZ,INDBLO)=HOSCAI*(XHERMI(IZEROS)-HOSHII)
C
         END DO
C
C=======================================================================
C     RIGHT WAVE FUNCTION
C=======================================================================
C
         DO IZEROS = 1,NGAUSS


            XZER = HOSCAJ*(XHERMI(IZEROS)-HOSHIJ)
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(HOSCAJ)
     *                                   / HERFAC( NOSCIL )
C
                DHERI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                                          (DHERMI(NOSCIL+1)
     *                                   - XZER*PHERMI(NOSCIL+1))
     *                             * SQRT(HOSCAJ) *HOSCAJ
     *                                   / HERFAC( NOSCIL )
C
                DDHRI4(NOSCIL,IZEROS,KARTEZ,INDBLO) =
     *                               (XZER*XZER-2*NOSCIL-1)
     *                                   * PHERMI(NOSCIL+1)
     *                      * SQRT(HOSCAJ) *HOSCAJ*HOSCAJ
     *                                   / HERFAC( NOSCIL )
            END DO
C
C
C=======================================================================
C      SAVING THE POINTS TO COMPUTE THE EXPONENTIAL PART OF THE RIGHT WF
C=======================================================================
C
            XZER_J(IZEROS,KARTEZ,INDBLO)=HOSCAJ*(XHERMI(IZEROS)-HOSHIJ)
C
         END DO
C
      END DO
C
C
C=======================================================================
C     THE SECOND ITERATION REFERS TO THE SECOND CENTER FOR BOTH
C     WAVE FUNCTIONS, WICH CORRESPONDS TO THE LAST INDEX IN DENSHF.
C     THIS MUST BE DONE TO USE THE PROPER ARRAYS INSIDE DENSHF
C=======================================================================
C
      HERLE4(:,:,:,NDTWHE)=HERLE4(:,:,:,NDTWCE)
      DHELE4(:,:,:,NDTWHE)=DHELE4(:,:,:,NDTWCE)
      DDHLE4(:,:,:,NDTWHE)=DDHLE4(:,:,:,NDTWCE)
      HERRI4(:,:,:,NDTWHE)=HERRI4(:,:,:,NDTWCE)
      DHERI4(:,:,:,NDTWHE)=DHERI4(:,:,:,NDTWCE)
      DDHRI4(:,:,:,NDTWHE)=DDHRI4(:,:,:,NDTWCE)
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_TRUEXP
C=======================================================================
      USE hfodd_sizes
C=======================================================================
C
      COMMON
     *       /TRUEXP/ TRUEXP(NDXHRM,NDYHRM,NDZHRM,NDTWBL)
      COMMON
     *       /GAUMAX/ NXHERM,NYHERM,NZHERM
      COMMON
     *       /TRUHER/ XZER_I(1:NDGAUS,1:NDKART,1:NDTWBL),
     *                XZER_J(1:NDGAUS,1:NDKART,1:NDTWBL)
C
C=======================================================================
C     IN THE SPIRIT OF DEFEXP HERE WE COMPUTE THE EXPONENTIAL PART OF
C     THE PRODUCT OF TWO HO WAVE FUNCTIONS, BUT IN THE REAL POSITIONS
C     IN SPACE. JUST FOR SAVING IN REVIEW FILE.
C=======================================================================
C
      DO INDBLO=1,NDTWBL
C
      DO IZ=1,NZHERM
         DO IY=1,NYHERM
            DO IX=1,NXHERM
C
         TRUEXP(IX,IY,IZ,INDBLO)=EXP(-0.5D0*(XZER_I(IX,1,INDBLO)**2
     *                                      +XZER_I(IY,2,INDBLO)**2
     *                                      +XZER_I(IZ,3,INDBLO)**2))
     *                          *EXP(-0.5D0*(XZER_J(IX,1,INDBLO)**2
     *                                      +XZER_J(IY,2,INDBLO)**2
     *                                      +XZER_J(IZ,3,INDBLO)**2))
C
            END DO
         END DO
      END DO
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C=======================================================================
C
      SUBROUTINE TWC_REVDEN(ISIMTX,JSIMTY,ISIMTZ,
     *            ISIGNY,ISIMPY,ISIQTY,MPAHFB,MREVER,
     *                                 MIN_QP,IPNMIX,
     *                   ITIREP,PRINIT,IDEVAR,ITERUN,
     *            ISYMDE,INIROT,INIINV,INIKAR,ITWCEN)
C=======================================================================
      USE hfodd_sizes
      USE WAVR_L
      USE SAVRIG
      USE HE_DEN
      USE TRUDEN
C=======================================================================
      COMMON
     *       /GAUMAX/ NXHERM,NYHERM,NZHERM
      COMMON
     *       /AXILIZ/ NOSCIL,IF_THO,IAXIAP
      COMMON
     *       /TRUEXP/ TRUEXP(NDXHRM,NDYHRM,NDZHRM,NDTWBL)
      COMMON
     *       /CFIPRI/ NFIPRI
C
      CHARACTER
     *          NAMEPN*8
      LOGICAL
     *          PRINIT
C=======================================================================
C     THIS SUBROUTINE, CALLED BEFORE REVIEW, COMPUTES THE DENSITIES IN
C     THE REAL SPACE COORDINATES. STEPS INVOLVED:
C     1. COMPUTATION OF THE HERMITE POLYNOMIALS IN THE GOOD LATTICE.
C     2. COMPUTATION OF THE EXPONENTIAL FACTORS IN THE GOOD LATTICE.
C     3. COMPUTATION OF DENSITIES USING AGAIN DENSHF. CALLED 3 TIMES.
C     4. SUM OF 4 CONTRIBUTIONS FROM DIFFERENT CENTERS.
C     5. SAVING THE RESULT IN THE FIRST INDEX OF THE DENSITIES.
C=======================================================================
      IF (ITWCEN.NE.2) STOP ' WRONG ITWCEN IN TWC_REVDEN'
C=======================================================================
      IALLER=0
C
      IF(ALLOCATED(REVRHO)) THEN
         DEALLOCATE(REVRHO,REVTAU,REVLPR,REVDIV)
      END IF
C
      IALLOC=0
C
      ALLOCATE(REVRHO(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:NDISOS),
     *         REVTAU(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:NDISOS),
     *         REVLPR(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:NDISOS),
     *         REVDIV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:NDISOS),
     *         STAT=IALLOC)
C
      IF (IALLOC.NE.0) THEN
          WRITE(NFIPRI,'(''ERROR ALLOCATING ARRAYS 1 IN TWC_REVDEN'')')
          IALLER=1
      END IF
C
C=======================================================================
C
      IF(ALLOCATED(REVKIS)) THEN
         DEALLOCATE(REVKIS,REVGRR,REVLPS,REVROS,REVROC,
     *              REVCUR,REVKIF,REVSPI)
      END IF
C
      IALLOC=0
C
      ALLOCATE(REVKIS(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         REVGRR(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         REVLPS(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         REVROS(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         REVROC(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         REVCUR(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         REVKIF(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         REVSPI(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,0:NDISOS),
     *         STAT=IALLOC)
C
      IF (IALLOC.NE.0) THEN
          WRITE(NFIPRI,'(''ERROR ALLOCATING ARRAYS 2 IN TWC_REVDEN'')')
          IALLER=1
      END IF
C
C=======================================================================
C
      IF(ALLOCATED(REVSCU)) THEN
         DEALLOCATE(REVSCU,REVDES)
      END IF
C
      IALLOC=0
C
      ALLOCATE(REVSCU(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,1:NDKART,
     *                                                    0:NDISOS),
     *         REVDES(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,1:NDKART,
     *                                                    0:NDISOS),
     *         STAT=IALLOC)
C
      IF (IALLOC.NE.0) THEN
          WRITE(NFIPRI,'(''ERROR ALLOCATING ARRAYS 3 IN TWC_REVDEN'')')
          IALLER=1
      END IF

C=======================================================================
      NOSACT=NOSCIL
C
      NGAUSS=NXHERM
      KARTEZ=1
C
      CALL TWC_TRUHER(NOSACT,KARTEZ,NGAUSS)
      CALL TWC_COPHER(NOSACT,KARTEZ,NGAUSS,ITWCEN)
C
      NGAUSS=NYHERM
      KARTEZ=2
C
      CALL TWC_TRUHER(NOSACT,KARTEZ,NGAUSS)
      CALL TWC_COPHER(NOSACT,KARTEZ,NGAUSS,ITWCEN)
C
      NGAUSS=NZHERM
      KARTEZ=3
C
      CALL TWC_TRUHER(NOSACT,KARTEZ,NGAUSS)
      CALL TWC_COPHER(NOSACT,KARTEZ,NGAUSS,ITWCEN)
C
      CALL TWC_TRUEXP

      CALL ZEDENS(0)
      CALL ZEDENS(1)
C
      DO ICHARG=0,NDISOS

      IF (ICHARG.EQ.0)  THEN
          NAMEPN='NEUTRONS'
          ITPNMX=0
      ELSEIF (ICHARG.EQ.1) THEN
          NAMEPN='PROTONS'
          ITPNMX=1
      END IF

C
C=======================================================================
C     CALLING DENSHF TO COMPUTE IN THE PROPER LATTICE
C=======================================================================
C
      ISAWAV=1
      IKERNE=0
C
C=======================================================================
C     CALLING DENSHF FOR I=J=K=L=1 (SAME AS IN ONE-CENTRE OPTION)
C=======================================================================
C
      CALL DENSHF(ISIMTX,JSIMTY,ISIMTZ,
     *            ISIGNY,ISIMPY,ISIQTY,MPAHFB,MREVER,ICHARG,
     *                                 MIN_QP,IPNMIX,ITPNMX,
     *                   ITIREP,NAMEPN,PRINIT,IDEVAR,ITERUN,
     *            ISYMDE,INIROT,INIINV,INIKAR,ISAWAV,IKERNE,
     *                                             ITWCEN,1)
C
C=======================================================================
C     CALLING DENSHF FOR I=J=K=L=2
C=======================================================================
C
      CALL DENSHF(ISIMTX,JSIMTY,ISIMTZ,
     *            ISIGNY,ISIMPY,ISIQTY,MPAHFB,MREVER,ICHARG,
     *                                 MIN_QP,IPNMIX,ITPNMX,
     *                   ITIREP,NAMEPN,PRINIT,IDEVAR,ITERUN,
     *            ISYMDE,INIROT,INIINV,INIKAR,ISAWAV,IKERNE,
     *                                        ITWCEN,16)
C
       IKERNE=1
       ISAWAV=0
C
       IALLOC=0
C
       IF (.NOT.ALLOCATED(WALEFT)) THEN
           ALLOCATE (WALEFT(1:NDBASE,1:2*NDSTAT,0:NDSPIN),
     *                                                   STAT=IALLOC)
           IF (IALLOC.NE.0) CALL NOALLO('WALEFT','HFODD ')
       END IF
C
C=======================================================================
C     SELECTING THE COEFFICIENTS FOR THE REST OF DENSITIES
C=======================================================================
C
       WALEFT(:,:,:)=SARIGH(:,:,:,ICHARG,2)
       WARIGH(:,:,:)=SARIGH(:,:,:,ICHARG,1)

       CALL DENSHF(ISIMTX,JSIMTY,ISIMTZ,
     *            ISIGNY,ISIMPY,ISIQTY,MPAHFB,MREVER,ICHARG,
     *                                 MIN_QP,IPNMIX,ITPNMX,
     *                   ITIREP,NAMEPN,PRINIT,IDEVAR,ITERUN,
     *            ISYMDE,INIROT,INIINV,INIKAR,ISAWAV,IKERNE,
     *                                        ITWCEN,3)
C
      END DO
C
C=======================================================================
C     SAVING DENSITIES
C=======================================================================
C
      IND1=1
      IND2=1
      IN22=1
      IND3=1
C
      IF (ITWCEN.EQ.2) THEN
C
          IND1=1
          IND2=2
          IN22=16
          IND3=3
C
      END IF
      DO ICHARG=0,NDISOS
C
      DO IX=1,NXHERM
         DO IY=1,NYHERM
            DO IZ=1,NZHERM
C
               REVRHO(IX,IY,IZ,ICHARG)=
     *               DE_RHO(IX,IY,IZ,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +             DE_RHO(IX,IY,IZ,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +     2.0d0 * DE_RHO(IX,IY,IZ,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
               REVTAU(IX,IY,IZ,ICHARG)=
     *               DE_TAU(IX,IY,IZ,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +             DE_TAU(IX,IY,IZ,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +     2.0d0 * DE_TAU(IX,IY,IZ,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
               REVLPR(IX,IY,IZ,ICHARG)=
     *               DE_LPR(IX,IY,IZ,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +             DE_LPR(IX,IY,IZ,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +     2.0d0 * DE_LPR(IX,IY,IZ,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C

C
               REVDIV(IX,IY,IZ,ICHARG)=
     *               DE_DIV(IX,IY,IZ,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +             DE_DIV(IX,IY,IZ,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +     2.0d0 * DE_DIV(IX,IY,IZ,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                  DO K=1,NDKART
C
                    REVKIS(IX,IY,IZ,K,ICHARG)=
     *              DE_KIS(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_KIS(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_KIS(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                    REVGRR(IX,IY,IZ,K,ICHARG)=
     *              DE_GRR(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_GRR(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_GRR(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                    REVLPS(IX,IY,IZ,K,ICHARG)=
     *              DE_LPS(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_LPS(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_LPS(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                    REVROS(IX,IY,IZ,K,ICHARG)=
     *              DE_ROS(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_ROS(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_ROS(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                    REVROC(IX,IY,IZ,K,ICHARG)=
     *              DE_ROC(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_ROC(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_ROC(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                    REVCUR(IX,IY,IZ,K,ICHARG)=
     *              DE_CUR(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_CUR(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_CUR(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                    REVKIF(IX,IY,IZ,K,ICHARG)=
     *              DE_KIF(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_KIF(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_KIF(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                    REVSPI(IX,IY,IZ,K,ICHARG)=
     *              DE_SPI(IX,IY,IZ,K,ICHARG,IND1)*TRUEXP(IX,IY,IZ,IND1)
     * +            DE_SPI(IX,IY,IZ,K,ICHARG,IN22)*TRUEXP(IX,IY,IZ,IND2)
     * +    2.0d0 * DE_SPI(IX,IY,IZ,K,ICHARG,IND3)*TRUEXP(IX,IY,IZ,IND3)
C
                     DO L=1,NDKART
C
                        REVSCU(IX,IY,IZ,K,L,ICHARG)=
     *                  DE_SCU(IX,IY,IZ,K,L,ICHARG,IND1)
     *                *            TRUEXP(IX,IY,IZ,IND1)
     * +                DE_SCU(IX,IY,IZ,K,L,ICHARG,IN22)
     *                *            TRUEXP(IX,IY,IZ,IND2)
     * +        2.0d0 * DE_SCU(IX,IY,IZ,K,L,ICHARG,IND3)
     *                *            TRUEXP(IX,IY,IZ,IND3)
C
                        REVDES(IX,IY,IZ,K,L,ICHARG)=
     *                  DE_DES(IX,IY,IZ,K,L,ICHARG,IND1)
     *               *             TRUEXP(IX,IY,IZ,IND1)
     * +                DE_DES(IX,IY,IZ,K,L,ICHARG,IN22)
     *               *             TRUEXP(IX,IY,IZ,IND2)
     * +       2.0d0 *  DE_DES(IX,IY,IZ,K,L,ICHARG,IND3)
     *               *             TRUEXP(IX,IY,IZ,IND3)
C
                     END DO !L
C
                  END DO !K
C
               END DO !IX
            END DO !IY
         END DO !IZ
C
      END DO !ICHARG
C
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_INTRAD(HAMRAD,KARTEZ,NUMORD,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      PARAMETER (NDMORD=3,NDTORD=10)
C=======================================================================
      DIMENSION
     *          HAMRAD(1:NDBASE,1:NDBASE)
      DIMENSION
     *          N_KART(1:NDKART)
C
      COMMON
     *       /MATORD/ NUTORD(0:NDMORD),
     *                MUTORD(1:NDTORD,0:NDMORD),
     *                IUTORD(1:NDTORD,0:NDMORD,1:NDKART)
      COMMON
     *     /COMU2C/ CMUTWC(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART,1:NDTWBL)
      COMMON
     *       /BASISO/ NXVECT(1:NDBASE,1:NDTWCE),
     *                NYVECT(1:NDBASE,1:NDTWCE),
     *                NZVECT(1:NDBASE,1:NDTWCE)
      COMMON
     *       /RADUNI/ RUNITS(0:NDMORD)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C=======================================================================
C      THIS SUBROUTINE CALCULATES:
C      FOR KARTEZ>0, MATRIX ELEMENTS OF X**2, Y**2, OR Z**2
C      FOR KARTEZ=0, MATRIX ELEMENTS OF R**(2*NUMORD)
C=======================================================================
C
      HAMRAD(:,:)=0.0D0
C
      ITWBRA=IND2HL(LTWCEN)
      ITWKET=IND2HL(LTWCEN)
C
      IF (KARTEZ.GT.0) THEN
C
          N_KART(1)=0
          N_KART(2)=0
          N_KART(3)=0
          N_KART(KARTEZ)=2
C
C=======================================================================
C        CALCULATING THE MATRIX ELEMENT MULTIPLIED BY THE UNITS DEFINED
C        ONCE FOR ALL IN THE SUBROUTINE "DEFUNI" AND BY THE PHASE.
C=======================================================================
C
          DO IKET=1,LDBATW(IND2HR(LTWCEN))
             DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                HAMRAD(IBRA,IKET) =
     *
     *             CMUTWC(N_KART(1),NXVECT(IBRA,ITWBRA),
     *                              NXVECT(IKET,ITWKET),1,LTWCEN)
     *           * CMUTWC(N_KART(2),NYVECT(IBRA,ITWBRA),
     *                              NYVECT(IKET,ITWKET),2,LTWCEN)
     *           * CMUTWC(N_KART(3),NZVECT(IBRA,ITWBRA),
     *                              NZVECT(IKET,ITWKET),3,LTWCEN)
     *           * RUNITS(1)
C
             END DO
          END DO
C
      ELSE
C
          DO IKET=1,LDBATW(IND2HR(LTWCEN))
             DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                DO LUTORD=1,NUTORD(NUMORD)
C
                   HAMRAD(IBRA,IKET) = HAMRAD(IBRA,IKET)
     *           + MUTORD(       LUTORD,NUMORD  )
     *           * CMUTWC(IUTORD(LUTORD,NUMORD,1),NXVECT(IBRA,ITWBRA),
     *                           NXVECT(IKET,ITWKET),1,LTWCEN)
     *           * CMUTWC(IUTORD(LUTORD,NUMORD,2),NYVECT(IBRA,ITWBRA),
     *                           NYVECT(IKET,ITWKET),2,LTWCEN)
     *           * CMUTWC(IUTORD(LUTORD,NUMORD,3),NZVECT(IBRA,ITWBRA),
     *                           NZVECT(IKET,ITWKET),3,LTWCEN)
     *           * RUNITS(NUMORD)
C
                END DO
C
             END DO
          END DO
C
      END IF
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWDECO(ICHARG,LTWCEN)
C=======================================================================
      USE PAIDEL
      USE TWCDEL
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C=======================================================================
C
      IALLOC=0
C
C=======================================================================
      IF (.NOT.ALLOCATED(TWDEPP).AND.LTWCEN.GE.1) THEN
          ALLOCATE (TWDEPP(1:2*NDBASE,1:2*NDBASE,0:NDREVE,
     *                                           0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWDEPP','TWDECO')
      END IF
C=======================================================================
      IF (.NOT.ALLOCATED(TWDEPM).AND.LTWCEN.GE.1) THEN
          ALLOCATE (TWDEPM(1:2*NDBASE,1:2*NDBASE,0:NDREVE,
     *                                           0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('TWDEPM','TWDECO')
      END IF
C=======================================================================
C          THIS SUBROUTINE COPIES THE CALCULATED MATRIX ELEMENTS OF THE
C          PAIRING FIELD INTO THE TWO-CENTER STRUCTURE IN THE STYLE OF
C          TWCCOP.
C=======================================================================
C
      IF (LTWCEN.EQ.1) THEN

          TWDEPP=CMPLX(0.0D0,0.0D0)
          TWDEPM=CMPLX(0.0D0,0.0D0)
C
          MBRA=0
          MKET=0
C
      END IF
C
      IF (LTWCEN.EQ.2) THEN
C
          MBRA=LDBATW(1)
          MKET=LDBATW(1)
C
      END IF
C
      IF (LTWCEN.EQ.3) THEN
C
          MBRA=LDBATW(1)
          MKET=0
C
      END IF
C
C
C=======================================================================
C     SETTING THE MATRIX ELEMENTS
C=======================================================================
C
      IF (LTWCEN.LT.3) THEN
C
      DO IKET=1,LDBATW(IND2HR(LTWCEN))
         DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
            TWDEPP(MBRA+IBRA,MKET+IKET,0,ICHARG)= HAMIDE(IBRA,IKET,1)
            TWDEPP(MBRA+IBRA,MKET+IKET,1,ICHARG)= HAMIDE(IBRA,IKET,1)
C
C
            TWDEPM(MBRA+IBRA,MKET+IKET,1,ICHARG)= HAMIDE(IBRA,IKET,0)
            TWDEPM(MBRA+IBRA,MKET+IKET,0,ICHARG)=-HAMIDE(IKET,IBRA,0)
C
         END DO
      END DO
C
      ELSEIF (LTWCEN.EQ.3) THEN
C
          DO IKET=1,LDBATW(IND2HR(LTWCEN))
              DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
               TWDEPP(MBRA+IBRA,MKET+IKET,0,ICHARG)=HAMIDE(IBRA,IKET,1)
               TWDEPP(MBRA+IBRA,MKET+IKET,1,ICHARG)=HAMIDE(IBRA,IKET,1)
C
               TWDEPM(MBRA+IBRA,MKET+IKET,1,ICHARG)= HAMIDE(IBRA,IKET,0)
               TWDEPM(MBRA+IBRA,MKET+IKET,0,ICHARG)=-HAMIDE(IBRA,IKET,0)
C
             END DO
          END DO
C
          DO IKET=1,LDBATW(IND2HL(LTWCEN))
              DO IBRA=1,LDBATW(IND2HR(LTWCEN))

C
               TWDEPP(MKET+IBRA,MBRA+IKET,0,ICHARG)=-HAMIDE(IKET,IBRA,1)
               TWDEPP(MKET+IBRA,MBRA+IKET,1,ICHARG)=-HAMIDE(IKET,IBRA,1)
C
               TWDEPM(MKET+IBRA,MBRA+IKET,1,ICHARG)= HAMIDE(IKET,IBRA,0)
               TWDEPM(MKET+IBRA,MBRA+IKET,0,ICHARG)=-HAMIDE(IKET,IBRA,0)
C
             END DO
          END DO
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_FRGMAT(ICHARG)
C=======================================================================
      USE BIGOVE
      USE BIGDEN
      USE BIGFRA
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMPLEX
     *          C_ZERO
C=======================================================================
      COMMON
     *       /FRCONS/ FRGSTF(0:NDISOS,1:NDTWCE),
     *                FRGASK(0:NDISOS,1:NDTWCE),
     *                IFLFRG(0:NDISOS,1:NDTWCE),
     *                CLMFRA(0:NDISOS,1:NDTWCE),
     *                AVEFRA(0:NDISOS,1:NDTWCE)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C=======================================================================
C
      C_ZERO=CMPLX(0.0D0,0.0D0)
C
C=======================================================================
C
      IALLOC=0
      IF (.NOT.ALLOCATED(FRAMAT)) THEN
          ALLOCATE (FRAMAT(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('FRAMAT','TWC_FRGMAT')
      END IF
C
C=======================================================================
C     OBTAINING MATRIX ELEMENTS TO INSERT IN THE ROUTHIAN
C=======================================================================
C
      FRAMAT(:,:)=C_ZERO
C
      CONSFA=CLMFRA(ICHARG,1)
      CONSFB=CLMFRA(ICHARG,NDTWCE)
C
      IF (IFLFRG(ICHARG,1).EQ.0.OR.IFLFRG(ICHARG,1).EQ.3)
     *    CONSFA = 0.0D0
      IF (IFLFRG(ICHARG,NDTWCE).EQ.0.OR.IFLFRG(ICHARG,NDTWCE).EQ.3)
     *    CONSFB = 0.0D0
C
      IF (CONSFA.LT.1E-15.AND.CONSFB.LT.1E-15) GO TO 5000
C
      DO IBRA=1,2*LDBTOT
C
         DO IKET=1,2*LDBTOT
C
            IF (IBRA.LE.LDBATW(1))                            IBRACEN=1
            IF (IBRA.GT.LDBATW(1).AND.IBRA.LE.LDBTOT)         IBRACEN=2
            IF (IBRA.GT.LDBTOT.AND.IBRA.LE.LDBTOT+LDBATW(1))  IBRACEN=1
            IF (IBRA.GT.LDBTOT+LDBATW(1))                     IBRACEN=2
C
            IF (IKET.LE.LDBATW(1))                            IKETCEN=1
            IF (IKET.GT.LDBATW(1).AND.IKET.LE.LDBTOT)         IKETCEN=2
            IF (IKET.GT.LDBTOT.AND.IKET.LE.LDBTOT+LDBATW(1))  IKETCEN=1
            IF (IKET.GT.LDBTOT+LDBATW(1))                     IKETCEN=2
C
C
            IF (IBRACEN.EQ.1.AND.IKETCEN.EQ.1) THEN

                FRAMAT(IBRA,IKET)=CONSFA*BIGNOR(IBRA,IKET)
C
            ELSEIF (IBRACEN.NE.IKETCEN) THEN
C
                FRAMAT(IBRA,IKET)=(CONSFA+CONSFB)/2
     *                            *BIGNOR(IBRA,IKET)
C
            ELSEIF (IBRACEN.EQ.2.AND.IKETCEN.EQ.2) THEN
C
                FRAMAT(IBRA,IKET)=CONSFB*BIGNOR(IBRA,IKET)
C
            END IF
C
         END DO
C
      END DO
C
C=======================================================================
C
 5000 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_FRGENE(ICHARG,ITWCEN,ENFRCO)
C=======================================================================
      USE BIGOVE
      USE BIGDEN
      USE BIGFRA
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMPLEX
     *          C_ZERO
C
      COMPLEX
     *          RESULA,RESULB,RESULT
C=======================================================================
      COMMON
     *       /FRCONS/ FRGSTF(0:NDISOS,1:NDTWCE),
     *                FRGASK(0:NDISOS,1:NDTWCE),
     *                IFLFRG(0:NDISOS,1:NDTWCE),
     *                CLMFRA(0:NDISOS,1:NDTWCE),
     *                AVEFRA(0:NDISOS,1:NDTWCE)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C=======================================================================
C
      C_ZERO=CMPLX(0.0D0,0.0D0)
C
C=======================================================================
C     OBTAINING AVERAGE PARTICLE NUMBER IN BOTH FRAGMENTS
C=======================================================================
C
      RESULA=C_ZERO
      RESULB=C_ZERO
C
      DO IBRA=1,2*LDBTOT
C
         DO IKET=1,2*LDBTOT
C
            IF (IKET.LE.LDBATW(1).OR.IKET.GT.LDBTOT.
     *                        AND.IKET.LE.LDBTOT+LDBATW(1)) THEN
C
                RESULA=RESULA+BIGNOR(IBRA,IKET)*BIGTWC(IKET,IBRA)
C
            ELSE
C
                RESULB=RESULB+BIGNOR(IBRA,IKET)*BIGTWC(IKET,IBRA)
C
            END IF
C
         END DO
C
      END DO
C
      AVEFRA(ICHARG,1)     = REAL(RESULA)
      AVEFRA(ICHARG,NDTWCE)= REAL(RESULB)
C
      ICHECK=0
C
      DO LTWCEN=1,ITWCEN
C
         IF (IFLFRG(ICHARG,LTWCEN).EQ.1.OR.IFLFRG(ICHARG,LTWCEN).EQ.2)
     *       ICHECK=ICHECK+1
C
      END DO
C
      IF (ICHECK.GT.0) THEN
C
          RESULT=C_ZERO
C
          DO LTWCEN=1,ITWCEN
C
             RESULT=RESULT-CLMFRA(ICHARG,LTWCEN)*AVEFRA(ICHARG,LTWCEN)
C
          END DO
C
          ENFRCO=REAL(RESULT/2.0D0)
C
C=======================================================================
C               UPDATE THE LAGRANGE MULTIPLIERS
C=======================================================================
C
          DO LTWCEN=1,ITWCEN
C
             IF (IFLFRG(ICHARG,LTWCEN).EQ.1.OR.
     *           IFLFRG(ICHARG,LTWCEN).EQ.2) THEN
C
                 CLMFRA(ICHARG,LTWCEN)=CLMFRA(ICHARG,LTWCEN)
     *             - (AVEFRA(ICHARG,LTWCEN) - FRGASK(ICHARG,LTWCEN))
     *             *  FRGSTF(ICHARG,LTWCEN)*2.0D0
C
             ELSE
C
                 CLMFRA(ICHARG,LTWCEN)=0.0D0
C
             END IF
C
          END DO
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_PRIFRA(NMOMUL)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMMON
     *       /FRCONS/ FRGSTF(0:NDISOS,1:NDTWCE),
     *                FRGASK(0:NDISOS,1:NDTWCE),
     *                IFLFRG(0:NDISOS,1:NDTWCE),
     *                CLMFRA(0:NDISOS,1:NDTWCE),
     *                AVEFRA(0:NDISOS,1:NDTWCE)
      COMMON
     *       /CFIPRI/ NFIPRI
      COMMON
     *       /QMUTWC/ QFRA_N(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE),
     *                QFRA_P(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE),
     *                QFRA_T(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE)
      CHARACTER
     *          PRIMOM(-NDMULT:NDMULT)*15
      CHARACTER
     *          PRIMIU*2
C
C=======================================================================
C         THIS SUBROUTINE PRINTS THE PROPERTIES OF THE FRAGMENTS
C         IN THE TWO-CENTER BASIS FORMALISM.
C         AT THE MOMENT IS ONLY THE NUMBER OF PARTICLES.
C=======================================================================
C
      WRITE(NFIPRI,'(79(1H*),/,                        1H*,77X,1H*,/,
     *             1H*,26X,'' PROPERTIES OF FRAGMENTS'',   27X,1H*,/,
     *                                                 1H*,77X,1H*,/,
     *                                                     79(1H*),/,
     *                                                 1H*,77X,1H*,/,
     *             1H*, 2X,''NEUTRONS:(FR1)='',F13.6,4X,
     *                              ''(FR2)='',F13.6,4X,
     *                              ''(TOT)='',F13.6,       1X,1H*)')
     *
     *       AVEFRA(0,1),AVEFRA(0,NDTWCE),AVEFRA(0,1)+AVEFRA(0,NDTWCE)
C
      WRITE(NFIPRI,'(1H*,77X,1H*,/,
     *             1H*, 2X,''PROTONS :(FR1)='',F13.6,4X,
     *                              ''(FR2)='',F13.6,4X,
     *                              ''(TOT)='',F13.6,       1X,1H*)')
     *
     *       AVEFRA(1,1),AVEFRA(1,NDTWCE),AVEFRA(1,1)+AVEFRA(1,NDTWCE)
C
C=======================================================================
C     MULTIPOLDE DEFORMATIONS
C=======================================================================
C
      DO LTWCEN=1,NDTWCE
c
      IF (LTWCEN.EQ.1) THEN
      WRITE(NFIPRI,'(1H*,77X,1H*)')
      WRITE(NFIPRI,'(79(1H*),/,                        1H*,77X,1H*,/,
     *             1H*,26X,''MULTIP. MOMENTS 1ST FRAG'',   27X,1H*,/,
     *                                                 1H*,77X,1H*,/,
     *                                                     79(1H*))')
      ELSEIF(LTWCEN.EQ.2) THEN
      WRITE(NFIPRI,'(79(1H*),/,                        1H*,77X,1H*,/,
     *             1H*,26X,''MULTIP. MOMENTS 2ND FRAG'',   27X,1H*,/,
     *                                                 1H*,77X,1H*,/,
     *                                                     79(1H*))')
      END IF
C
      DO MIU=-NDMULT,NDMULT
         PRIMOM(MIU)='  .............'
      END DO
C
      DO LAMBDA=0,MIN(NMOMUL,NDMULT)
C
         DO MIU=-LAMBDA,LAMBDA
C
            QMOMEN = QFRA_T(LAMBDA,MIU,LTWCEN)
C
            WRITE(PRIMIU,'(I2)') MIU
            IF (PRIMIU(1:2).EQ.' 0') PRIMIU(1:2)='0 '
            IF (PRIMIU(1:1).EQ.' ' ) PRIMIU(1:1)='+'
C
            IF (ABS(QMOMEN).LT.1.0E-10) THEN
C
                WRITE(PRIMOM(MIU),'(2X,A1,I1,A2,''='',4X,''ZERO'')')
     *                                          'Q',LAMBDA,PRIMIU
            ELSE
C
               IF (ABS(QMOMEN).LT.1.0E-3.OR.QMOMEN.GT.1000.OR.
     *                                      QMOMEN.LT.-100) THEN
C
                   WRITE(PRIMOM(MIU),'(2X,A1,I1,A2,''='',1PE8.1)')
     *                                 'Q',LAMBDA,PRIMIU,QMOMEN
               ELSE
C
                   WRITE(PRIMOM(MIU),'(2X,A1,I1,A2,''='',F8.4)')
     *                                'Q',LAMBDA,PRIMIU,QMOMEN
               END IF
C
            END IF
C
         END DO
C
C=======================================================================
C
         WRITE(NFIPRI,'(1H*,77X,1H*,/,1H*,5A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=0,4)
         WRITE(NFIPRI,'(1H*,     15X,4A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=-1,-4,-1)
C
C=======================================================================
C
         IF (LAMBDA.EQ.5) THEN
         WRITE(NFIPRI,'(1H*,77X,1H*,/,1H*,60X, A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=5,5)
         WRITE(NFIPRI,'(1H*,         60X, A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=-5,-5)
         END IF
C
C=======================================================================
C
         IF (LAMBDA.EQ.6) THEN
         WRITE(NFIPRI,'(1H*,77X,1H*,/,1H*,45X,2A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=5,6)
         WRITE(NFIPRI,'(1H*,         45X,2A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=-5,-6,-1)
         END IF
C
C=======================================================================
C
         IF (LAMBDA.EQ.7) THEN
         WRITE(NFIPRI,'(1H*,77X,1H*,/,1H*,30X,3A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=5,7)
         WRITE(NFIPRI,'(1H*,         30X,3A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=-5,-7,-1)
         END IF
C
C=======================================================================
C
         IF (LAMBDA.EQ.8) THEN
         WRITE(NFIPRI,'(1H*,77X,1H*,/,1H*,15X,4A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=5,8)
         WRITE(NFIPRI,'(1H*,         15X,4A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=-5,-8,-1)
         END IF
C
C=======================================================================
C
         IF (LAMBDA.EQ.9) THEN
         WRITE(NFIPRI,'(1H*,77X,1H*,/,1H*,    5A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=5,9)
         WRITE(NFIPRI,'(1H*,             5A15,2X,1H*)')
     *                                        (PRIMOM(MIU),MIU=-5,-9,-1)
         END IF
C
      END DO
C
      WRITE(NFIPRI,'(1H*,77X,1H*,/,79(1H*),/)')
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_CANQUZ(ENECAN)
C=======================================================================
      USE hfodd_sizes
      USE BIGHAM
      USE BIGDEN
C=======================================================================
      COMPLEX
     *          EWAVEF,CELMTS
      COMPLEX
     *          C_ZERO,C_UNIT
      COMPLEX
     *          AUXILI
      COMPLEX
     *          HAMCAN
      COMPLEX
     *          ENESUM
C=======================================================================
      ALLOCATABLE EWAVEF(:,:),CELMTS(:)
      ALLOCATABLE AUXILI(:,:)
C
      DIMENSION
     *          V2OCCU(1:4*NDBASE)
      DIMENSION
     *          HAMCAN(1:4*NDBASE,1:4*NDBASE)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C=======================================================================
      IALLOC=0
      IF (.NOT.ALLOCATED(EWAVEF)) THEN
          ALLOCATE (EWAVEF(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('EWAVEF','TWC_CANQUZ')
      END IF
      IF (.NOT.ALLOCATED(CELMTS)) THEN
          ALLOCATE (CELMTS(1:((4*NDBASE+1)*2*NDBASE)),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('CELMTS','TWC_CANQUZ')
      END IF
C=======================================================================
      C_ZERO=CMPLX(0.0D0,0.0D0)
      C_UNIT=CMPLX(1.0D0,0.0D0)
C=======================================================================
C
      ENESUM=0.0D0
C
      DO IBRA=1,2*LDBTOT
         DO IKET=1,2*LDBTOT
C
            IF (IBRA.GE.IKET) THEN

                NCOUNT=IBRA+((4*LDBTOT-IKET)*(IKET-1))/2
                CELMTS(NCOUNT)=-BIGTWC(IBRA,IKET)
C
            END IF
C
         END DO
      END DO
C
      EWAVEF=C_ZERO

      CALL DIAMAT(CELMTS,V2OCCU,EWAVEF,2*LDBTOT,4*NDBASE,2*LDBTOT)
C
      DO I_BASE=1,2*LDBTOT
C
         V2OCCU(I_BASE)=-V2OCCU(I_BASE)
C
      END DO
C
      ALLOCATE (AUXILI(1:4*NDBASE,1:4*NDBASE),STAT=IALLOC)
      IF (IALLOC.NE.0) CALL NOALLO('AUXILI','TWC_CANQUZ')
C
      CALL ZGEMM('C','N',2*LDBTOT,2*LDBTOT,2*LDBTOT,C_UNIT,EWAVEF,
     *                   4*NDBASE,COMHAM,4*NDBASE,C_ZERO,AUXILI,
     *                   4*NDBASE)
C
      CALL ZGEMM('N','N',2*LDBTOT,2*LDBTOT,2*LDBTOT,C_UNIT,AUXILI,
     *                   4*NDBASE,EWAVEF,4*NDBASE,C_ZERO,HAMCAN,
     *                   4*NDBASE)
C
      DO I_BASE=1,2*LDBTOT
      ENESUM=ENESUM+V2OCCU(I_BASE)*HAMCAN(I_BASE,I_BASE)
      END DO

      ENECAN=REAL(ENESUM)
C
      DEALLOCATE(AUXILI,CELMTS)
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_FRAGPAR(INNUMB,IZNUMB,LTWCEN,IZFRAG,INFRAG)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMMON
     *       /FRCONS/ FRGSTF(0:NDISOS,1:NDTWCE),
     *                FRGASK(0:NDISOS,1:NDTWCE),
     *                IFLFRG(0:NDISOS,1:NDTWCE),
     *                CLMFRA(0:NDISOS,1:NDTWCE),
     *                AVEFRA(0:NDISOS,1:NDTWCE)
      COMMON
     *       /CFIPRI/ NFIPRI
C
C=======================================================================
C         THIS SUBROUTINE COMPUTES THE PARTICLE NUMBERS FOR THE
C         NILSSON HAMILTONIAN
C=======================================================================
C
      IF (LTWCEN.EQ.1) THEN
C
          IF (IFLFRG(0,LTWCEN).EQ.1.OR.IFLFRG(0,LTWCEN).EQ.3) THEN
C
              INFRAG=INT(FRGASK(0,LTWCEN))
C
          ELSE
C
              INFRAG=INT(INNUMB/2)
C
          END IF
C
          IF (IFLFRG(1,LTWCEN).EQ.1.OR.IFLFRG(1,LTWCEN).EQ.3) THEN
C
              IZFRAG=INT(FRGASK(1,LTWCEN))
C
          ELSE
C
              IZFRAG=INT(IZNUMB/2)
C
          END IF
C
      ELSEIF (LTWCEN.EQ.2) THEN
C
          IF (IFLFRG(0,LTWCEN).EQ.1.OR.IFLFRG(0,LTWCEN).EQ.3) THEN
C
              INFRAG=INT(FRGASK(0,LTWCEN))
C
          ELSE
C
              IF (IFLFRG(0,LTWCEN-1).EQ.1.
     *         OR.IFLFRG(0,LTWCEN-1).EQ.3) THEN
C
                   INFRAG=INNUMB-INT(FRGASK(0,LTWCEN-1))
C
              ELSE
C
                  INFRAG=INNUMB-INT(INNUMB/2)
C
              END IF
C
          END IF
C
          IF (IFLFRG(1,LTWCEN).EQ.1.OR.IFLFRG(1,LTWCEN).EQ.3) THEN
C
              IZFRAG=INT(FRGASK(1,LTWCEN))
C
          ELSE
C
              IF (IFLFRG(1,LTWCEN-1).EQ.1.
     *         OR.IFLFRG(1,LTWCEN-1).EQ.3) THEN
C
                   IZFRAG=IZNUMB-INT(FRGASK(1,LTWCEN-1))
C
              ELSE
C
                  IZFRAG=IZNUMB-INT(IZNUMB/2)
C
              END IF
C
          END IF
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_SAVFIL(IPACOP,ICHARG,LTWCEN)
C=======================================================================
      USE MAT_PP
      USE MAT_PM
      USE PAIDEL
      USE MAFTWC
C=======================================================================
      USE hfodd_sizes
C=======================================================================
C
      IALLOC=0
c
      IF (IPACOP.EQ.0) THEN
C
          IF (.NOT.ALLOCATED(FITWPP)) THEN
              ALLOCATE (FITWPP(1:NDBASE,1:NDBASE,0:NDREVE,
     *                                  1:NDTWBL,0:NDISOS),STAT=IALLOC)
              IF (IALLOC.NE.0) CALL NOALLO('FITWPP','TWC_SAVFIL')
              FITWPP=CMPLX(0.0D0,0.0D0)
          END IF
C
          IF (.NOT.ALLOCATED(FITWPM)) THEN
              ALLOCATE (FITWPM(1:NDBASE,1:NDBASE,0:NDREVE,
     *                                  1:NDTWBL,0:NDISOS),STAT=IALLOC)
              IF (IALLOC.NE.0) CALL NOALLO('FITWPM','TWC_SAVFIL')
              FITWPM=CMPLX(0.0D0,0.0D0)
          END IF
C
          FITWPP(:,:,:,LTWCEN,ICHARG)=BIG_PP(:,:,:)
          FITWPM(:,:,:,LTWCEN,ICHARG)=BIG_PM(:,:,:)
C
      ELSEIF (IPACOP.EQ.1) THEN
C
          IF (.NOT.ALLOCATED(FITWDE)) THEN
              ALLOCATE (FITWDE(1:NDBASE,1:NDBASE,0:NDREVE,
     *                                  1:NDTWBL,0:NDISOS),STAT=IALLOC)
             IF (IALLOC.NE.0) CALL NOALLO('FITWDE','TWC_SAVFIL')
             FITWDE=CMPLX(0.0D0,0.0D0)
          END IF
C
          FITWDE(:,:,:,LTWCEN,ICHARG)=HAMIDE(:,:,:)
C
      END IF
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_FRGMOM(NMOMUL,QMOFRA)
C=======================================================================
      USE TWCOBO
      USE BIGDEN
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      COMPLEX
     *          C_ZERO
      COMPLEX
     *          RESULT,RESULA,RESULB
C=======================================================================
      DIMENSION
     *          QMOFRA(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE)
C=======================================================================
      ALLOCATABLE ELEMUL(:,:)
C=======================================================================
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
C
C=======================================================================
C
      IALLOC=0
C
C=======================================================================
      ALLOCATE (ELEMUL(1:NDBASE,1:NDBASE),STAT=IALLOC)
      IF (IALLOC.NE.0) CALL NOALLO('ELEMUL','TWC_FRGMOM')
C=======================================================================
      C_ZERO=CMPLX(0.0D0,0.0D0)
C
      QMOFRA=0.0D0
C
      DO ICENTE=1,NDTWCE
C
         DO LAMBDA=0,NMOMUL
            DO MIU=-LAMBDA,LAMBDA
C
               IF (MIU.LT.0) IMAREA=-1
               IF (MIU.GE.0) IMAREA= 1
C
               DO LTWCEN=1,NDTWBL
C
                  KARTEZ=0
                  CALL TWC_MULFRA(ELEMUL,LAMBDA,MIU,
     *                         IMAREA,LTWCEN,ICENTE)
                  CALL TWC_BIGSIM(ELEMUL,LTWCEN,IMAREA,KARTEZ)
C
               END DO
C
               CALL TWC_BIGOBO
C
               RESULA=C_ZERO
               RESULB=C_ZERO
C
               DO IBRA=1,2*LDBTOT
                  DO IKET=1,2*LDBTOT
C
                      IF (IKET.LE.LDBATW(1).OR.IKET.GT.LDBTOT.
     *                            AND.IKET.LE.LDBTOT+LDBATW(1)) THEN
C
                          RESULA=RESULA+BIGOBO(IBRA,IKET)
     *                                 *BIGTWC(IKET,IBRA)
C
                      ELSE
C
                          RESULB=RESULB+BIGOBO(IBRA,IKET)
     *                                 *BIGTWC(IKET,IBRA)
C
                      END IF
C
                  END DO
               END DO
C
               IF (ICENTE.EQ.1) RESULT=RESULA
               IF (ICENTE.EQ.2) RESULT=RESULB
C
               QMOFRA(LAMBDA,MIU,ICENTE)=REAL(RESULT)
C
            END DO
         END DO
C
      END DO
C
C=======================================================================
C
      DEALLOCATE (ELEMUL)
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_FRACON(NMUCON,HAMCON,LTWCEN)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      ALLOCATABLE ELEMUL(:,:)
C=======================================================================
      DIMENSION
     *          HAMCON(1:NDBASE,1:NDBASE)
      DIMENSION
     *          CMUL_T(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE)
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
     *       /QMUTWC/ QFRA_N(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE),
     *                QFRA_P(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE),
     *                QFRA_T(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE)
      COMMON
     *       /QTWCNS/ STQTWC(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE),
     *                QTWCAS(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE),
     *                IFLQTW(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE),
     *                GALMTW(0:NDMULT,-NDMULT:NDMULT,1:NDTWCE)
      COMMON
     *       /OURUNI/ QUNITS(0:NDMULT,0:NDMULT)
      COMMON
     *       /DIMTWC/ LDBATW(1:NDTWCE),LDBTOT
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C
C=======================================================================
C
      IALLOC=0
C
C=======================================================================
      ALLOCATE (ELEMUL(1:NDBASE,1:NDBASE),STAT=IALLOC)
      IF (IALLOC.NE.0) CALL NOALLO('ELEMUL','TWC_INTCON')
C
C=======================================================================
C
      HAMCON(:,:)=0.0D0
      ELEMUL(:,:)=0.0D0
C
C=======================================================================
C        CALCULATING THE COEFFICIENTS THAT IN THE CONSTRAINT TERMS
C        MULTIPLY THE MULTIPOLE MOMENTS
C=======================================================================
C
      DO ICENTE=1,NDTWCE
C
         DO LAMBDA=0,NMUCON
            DO MIU=-LAMBDA,LAMBDA
C
               NIU=IABS(MIU)
               CMUL_T(LAMBDA,MIU,ICENTE)=0.0D0
C
               IF (IFLQTW(LAMBDA,MIU,ICENTE).EQ.1) THEN
C
                   CMUL_T(LAMBDA,MIU,ICENTE)=
     *                                  2.0D0*STQTWC(LAMBDA,MIU,ICENTE)
     *      *      (QFRA_T(LAMBDA,MIU,ICENTE)-QTWCAS(LAMBDA,MIU,ICENTE))
     *                                       -GALMTW(LAMBDA,MIU,ICENTE)
C
               END IF
C
            END DO
         END DO
C
C=======================================================================
C        HERE WE CALCULATE CONTRIBUTIONS FROM REAL AND IMAGINARY
C        PARTS OF THE MULTIPOLE MOMENTS
C=======================================================================
C
         DO LAMBDA=0,NMUCON
            DO MIU=-LAMBDA,LAMBDA
C
               IF (IFLQTW(LAMBDA,MIU,ICENTE).EQ.1) THEN
C
                   IF (LTWCEN.EQ.1.AND.ICENTE.EQ.1) FACTOR=1.0D0
                   IF (LTWCEN.EQ.2.AND.ICENTE.EQ.1) FACTOR=0.0D0
                   IF (LTWCEN.EQ.3.AND.ICENTE.EQ.1) FACTOR=0.5D0
C
                   IF (LTWCEN.EQ.1.AND.ICENTE.EQ.2) FACTOR=0.0D0
                   IF (LTWCEN.EQ.2.AND.ICENTE.EQ.2) FACTOR=1.0D0
                   IF (LTWCEN.EQ.3.AND.ICENTE.EQ.2) FACTOR=0.5D0
C
                   IF (MIU.LT.0) IMAREA=-1
                   IF (MIU.GE.0) IMAREA= 1
C
                   CALL TWC_MULFRA(ELEMUL,LAMBDA,MIU,IMAREA,
     *                                        LTWCEN,ICENTE)
C
                   DO IKET=1,LDBATW(IND2HR(LTWCEN))
                      DO IBRA=1,LDBATW(IND2HL(LTWCEN))
C
                          HAMCON(IBRA,IKET)= HAMCON(IBRA,IKET)
     *                              +        ELEMUL(IBRA,IKET)
     *                              *CMUL_T(LAMBDA,MIU,ICENTE)
     *                              *                  FACTOR
C
                      END DO
                   END DO
C
                END IF
C
            END DO
         END DO
C
      END DO
C
C=======================================================================
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE TWC_ROTBAS(ITWCEN,LDTWCE)
C=======================================================================
      USE hfodd_sizes
C=======================================================================
      DIMENSION
     *          ROBAMA(1:NDKART,1:NDKART,1:NDTWCE),
     *          BDBAMA(1:NDKART,1:NDKART,1:NDTWCE),
     *          DVBAMA(1:NDKART,1:NDTWCE),
     *          DSBAMA(1:NDTWCE)
      DIMENSION
     *          ACOMAT(1:NDKART,1:NDKART),
     *          BCOMAT(1:NDKART)
      DIMENSION
     *          EIGVAL(1:NDKART)
      DIMENSION
     *          WORK(37),IWORK(18)
      DIMENSION
     *          BPCOMT(1:NDKART)
      COMMON
     *       /ROTFRG/ ALFFRA(1:NDTWCE),BETFRA(1:NDTWCE),
     *                GAMFRA(1:NDTWCE),IROFRG(1:NDTWCE)
      COMMON
     *       /SCALNG/ HOMSCA(NDKART,NDTWCE)
      COMMON
     *       /CENPOS/ CENCOR(1:NDKART,1:NDTWCE)
      COMMON
     *       /TWCEIN/ IND4HI(1:NDTWHE),IND4HJ(1:NDTWHE),
     *                IND4HK(1:NDTWDD),IND4HL(1:NDTWDD),
     *                IND2HR(1:NDTWBL),IND2HL(1:NDTWBL)
C=======================================================================
      IF (NDKART.NE.3) STOP 'NDKART MUST BE 3 IN TWC_ROTBAS'
      IF (ITWCEN.NE.2) STOP 'ITWCEN MUST BE 2 IN TWC_ROTBAS'
C=======================================================================
      ICHECK=0
C
      DO LTWCEN=1,ITWCEN
C
         ICHECK=ICHECK+IROFRG(LTWCEN)
C
         IF (IROFRG(LTWCEN).EQ.0) THEN
C
             ALFFRA(LTWCEN)=0.0D0
             BETFRA(LTWCEN)=0.0D0
             GAMFRA(LTWCEN)=0.0D0
C
         END IF
C
      END DO
C
      IF (ICHECK.EQ.0) GO TO 6000
C
C=======================================================================
C                         PART I: DEFINITIONS
C
C     -ROBAMA CONTAINS THE ROTATION MATRIX TIMES THE OSCILLATOR LENGTH
C
C         ROBAMA_MUNU^I=R_MUNU(ALPHA,BETA,GAMMA)^I*HOMSCA_NU^I
C
C     -BDBAMA, NEEDED TO BUILD THE QUADRATIC FORM OF THE GAUSSIAN, IS
C
C                           BDBAMA_MUNU^I=
C         SUM_LAMBDA ROBAMA_LAMBDA MU^I ROBAMA_LAMBDA NU^I
C
C     -DVBAMA, NEEDED TO BUILD THE QUADRATIC FORM OF THE GAUSSIAN, IS
C
C             DVBAMA_MU^I=SUM_NU BDBAMA_MUNU^I R_0 NU^I
C
C     -DSBAMA, NEEDED TO BUILD THE QUADRATIC FORM OF THE GAUSSIAN, IS
C
C           DSBAMA=SUM_MUNU BDBAMA_MUNU^I R_0 MU^I R_0 NU^I
C
C=======================================================================
C
      DO LTWCEN=1,ITWCEN
C
         ROBAMA(1,1,LTWCEN)=
     *       COS(BETFRA(LTWCEN))*COS(GAMFRA(LTWCEN))
         ROBAMA(1,2,LTWCEN)=
     *       SIN(ALFFRA(LTWCEN))*SIN(BETFRA(LTWCEN))*COS(GAMFRA(LTWCEN))
     *      -COS(ALFFRA(LTWCEN))*SIN(GAMFRA(LTWCEN))
         ROBAMA(1,3,LTWCEN)=
     *       COS(ALFFRA(LTWCEN))*SIN(BETFRA(LTWCEN))*COS(GAMFRA(LTWCEN))
     *      +SIN(ALFFRA(LTWCEN))*SIN(GAMFRA(LTWCEN))
C
         ROBAMA(2,1,LTWCEN)=
     *       COS(BETFRA(LTWCEN))*SIN(GAMFRA(LTWCEN))
         ROBAMA(2,2,LTWCEN)=
     *       SIN(ALFFRA(LTWCEN))*SIN(BETFRA(LTWCEN))*SIN(GAMFRA(LTWCEN))
     *      +COS(ALFFRA(LTWCEN))*COS(GAMFRA(LTWCEN))
         ROBAMA(2,3,LTWCEN)=
     *       COS(ALFFRA(LTWCEN))*SIN(BETFRA(LTWCEN))*SIN(GAMFRA(LTWCEN))
     *      -SIN(ALFFRA(LTWCEN))*COS(GAMFRA(LTWCEN))
C
C
         ROBAMA(3,1,LTWCEN)=
     *      -SIN(BETFRA(LTWCEN))
         ROBAMA(3,2,LTWCEN)=
     *       SIN(ALFFRA(LTWCEN))*COS(BETFRA(LTWCEN))
         ROBAMA(3,3,LTWCEN)=
     *       COS(ALFFRA(LTWCEN))*COS(BETFRA(LTWCEN))
C
         DO MU=1,NDKART
            DO NU=1,NDKART
C
               ROBAMA(MU,NU,LTWCEN)=
     *         ROBAMA(MU,NU,LTWCEN)*HOMSCA(NU,LTWCEN)
C
            END DO
         END DO
C
C   !TENSOR BDBAMA
C
         DO MU=1,NDKART
            DO NU=1,NDKART
C
               BDBAMA(MU,NU,LTWCEN)=0.0D0
               DO LAMBDA=1,NDKART
C
               BDBAMA(MU,NU,LTWCEN)= BDBAMA(MU,NU,LTWCEN)
     *                         + ROBAMA(MU,LAMBDA,LTWCEN)
     *                         * ROBAMA(NU,LAMBDA,LTWCEN)
C
               END DO
C
            END DO
         END DO
C
C   !VECTOR DVBAMA
C
         DO MU=1,NDKART
C
            DVBAMA(MU,LTWCEN)=0.0D0
            DO NU=1,NDKART
C
               DVBAMA(MU,LTWCEN)= DVBAMA(MU,   LTWCEN)
     *                         +  BDBAMA(MU,NU,LTWCEN)
     *                         *  CENCOR(   NU,LTWCEN)
C
            END DO
C
         END DO
C
C   !SCALAR DSBAMA
C
         DSBAMA(LTWCEN)=0.0D0
         DO MU=1,NDKART
            DO NU=1,NDKART
C
               DSBAMA(LTWCEN)= DSBAMA(LTWCEN)
     *                 + BDBAMA(MU,NU,LTWCEN)
     *                 * CENCOR(MU,   LTWCEN)
     *                 * CENCOR(   NU,LTWCEN)
C
            END DO
         END DO
C
      END DO
C
C   !COEFFICIENTS OF THE QUADRATIC FORM
C
      DO LTWCEN=1,LDTWCE
C
         I=IND4HI(LTWCEN)
         J=IND4HJ(LTWCEN)
         K=IND4HK(LTWCEN)
         L=IND4HL(LTWCEN)
C
C
C   !MATRIX A (QUADRATIC TERM)
C
         DO MU=1,NDKART
            DO NU=1,NDKART
C
               ACOMAT(MU,NU)=0.5D0*(BDBAMA(MU,NU,I)+BDBAMA(MU,NU,J)
     *                 +            BDBAMA(MU,NU,K)+BDBAMA(MU,NU,L))
C
            END DO
         END DO
C
C   !MATRIX B (LINEAR TERM)
C
         DO MU=1,NDKART
C
            BCOMAT(MU)=-DVBAMA(MU,I)-DVBAMA(MU,J)
     *                 -DVBAMA(MU,K)-DVBAMA(MU,L)
C
         END DO
C   !COEFFICIENT C
C
         CCOMAT=0.5D0*(DSBAMA(I)+DSBAMA(J)+DSBAMA(K)+DSBAMA(L))
C
C=======================================================================
C                         PART II: REMOVING CROSS TERMS
C      DIAGONALIZING MATRIX A TO DEFINE A NEW SET OF COORDINATES
C=======================================================================
C
         CALL DSYEVD('V','L',3,ACOMAT,3,EIGVAL,WORK,37,IWORK,18,INFO)
         IF (INFO.NE.0) THEN
             STOP 'ERROR IN DSYEVD IN TWC_ROTBAS'
         END IF
C
C=======================================================================
C                         PART III: DEFINE NEW LINEAR TERMS
C                              AND THE NEW SCALAR TERM
C=======================================================================
C
         DO MU=1,NDKART
C
            BPCOMT(MU)=0.0D0
            DO NU=1,NDKART
               BPCOMT(MU)=BPCOMT(MU)+BCOMAT(NU)*ACOMAT(NU,MU)
            END DO
C
         END DO
C
         CPCOMT=0.0D0
         DO MU=1,NDKART
            CPCOMT=CPCOMT-BPCOMT(MU)**2/(4*EIGVAL(MU))
         END DO
C
      END DO
C

 6000 CONTINUE
C
      RETURN
      END
C
C=======================================================================
C
      END MODULE hfodd_twocen
