      SUBROUTINE Set_version_hfodd_lipcorr
      !
      USE hfodd_sizes
      !
      CHARACTER NAMMOD*16
      INTEGER   MODUVE,MODSET
      !
      COMMON /VERMOD/ NAMMOD(NDMODU),MODUVE(NDMODU),MODSET(NDMODU)
      !
      NAMMOD(14) = 'hfodd_lipcorr   '
      MODUVE(14) = 54
      MODSET(14) = 14
      !
      END SUBROUTINE Set_version_hfodd_lipcorr
      !
!=======================================================================
!
MODULE OVRHFB
!
!=======================================================================
!     MODULE COMPUTING THE OVERLAP BETWEEN TWO HFB WAVE FUNCTIONS, BASED
!     ON PFAFFIANS, ACCORDING TO EQS. (7) FOR EVEN SYSTEMS,  OR (19) FOR
!     ODD ONES, OF BERTSCH AND ROBLEDO, PRL 108, 042505 (2012).
!=======================================================================
!     THE RUSSIAN CONVENTION HAS TO BE REMOVED FOR  THE LOWER COMPONENTS
!     OF THE HFB WAVE FUNCTIONS.
!=======================================================================
!     THE ARRAYS AMATRG/BMATRG (AMATLF/BMATLF) ARE THE A AND B  MATRICES
!     FOR THE RIGHT (LEFT) HFB STATE.  FOR AN ODD PARTICLE NUMBER, THESE
!     ARRAYS  CONTAIN ONLY THE  "EVEN" HFB VACUUM ON TOP OF WHICH THE QP
!     EXCITATIONS ARE BUILT. THE BLOCKED QP ARE STORED IN ARRAYS ABLOLF,
!     BBLOLF, ABLORG, BBLORG (q, p, q', p' IN EQ. (19) RESPECTIVELY).
!=======================================================================
!     THE PFAFFIAN IS COMPUTED VIA THE ROUTINE SKPF10 FROM  THE WIMMER'S
!     LIBRARY (ACM TRANS. MATH. SOFTWARE 38, 30 (2012)).
!=======================================================================
!     FROM VERSION 34, THIS MODULE ALSO INCLUDES THE ROUTINE OVRONI THAT
!     COMPUTE THE OVERLAP BETWEEN HFB STATES THROUGH THE ONISHI FORMULA.
!     THE  ONISHI  FORMULA  CONSIDERED INVOLVES A PREFACTOR AS WELL AS A
!     A NORMALIZATION,  IN CONTRAST TO THE EXPRESSION  USUALLY  GIVEN IN
!     THE LITERATURE.  THESE FACTORS COME FROM THE FACT THAT WE COMPUTE
!     IT FROM THE "PFAFFIAN" FORM OF HFB STATES,  AND NOT THEIR THOULESS
!     FORM.
!     THE MODULE ALSO INCLUSES AN EXTRA ROUTINE,  OVRMAT,  COMPUTING THE
!     MATRIX OF OVERLAPS BETWEEN THE QUASIPARTICLE STATES.
!=======================================================================
!     COMMENTS FROM DENSHF:
!     COMPONENTS OF THE SINGLE-PARTICLE (FOR IPAHFB=0) OR  QUASIPARTICLE
!     (FOR IPAHFB>0) WAVE FUNCTIONS  ARE  COLLECTED  IN  THE  SPACE-SPIN
!     REPRESENTATION IN THE "WARIGH" ARRAY. ALL THESE COMPONETS HAVE THE
!     FIRST INDEX SPANNING THE HO BASIS (1:LDBASE) AND THE  THIRD  INDEX
!     SPANNING THE TWO SPIN VALUES (1:NDSPIN). DIFFERENT  COMPONETS  ARE
!     STORED CONSECUTIVELY IN THE SECOND INDEX, WHICH  SPANS  VALUES  OF
!     (1:N*NDSTAT) WITH VALUES OF "N" ADJUSTED TO THE  CURRENT  SYMMETRY
!     CONDITIONS, AND GIVEN BY:
!     N=M   (FOR IPNMIX=0) AND
!     N=2*M (FOR IPNMIX=1), WITH INTEGERS M LISTED BELOW.
!=======================================================================
!     DENOTING THE COMPONENTS BY:
!                   TU=TIME-UP   (SIMPLEX=+i,IREVER=0)
!                   TD=TIME-DOWN (SIMPLEX=-i,IREVER=1)
!                   QU=UPPER HFB COMPONENT  (B*-ARRAYS)
!                   QD=LOWER HFB COMPONENT  (A*-ARRAYS),  WE HAVE:
!     FOR IPAHFB=1, ISIMPY=1, MREVER=0: M=2,(QU&TU,QD&TD)
!     FOR IPAHFB=1, ISIMPY=1, MREVER=1: M=4,(QU&TU,QU&TD,QD&TD,QD&TU)
!     FOR IPAHFB=1, ISIMPY=0, MREVER=1: M=4,(QU&TU,QU&TD,QD&TU,QD&TD)
!=======================================================================
!
USE BLOSAV
USE hfodd_sizes
USE F77_PFAPACK
USE SAVRIG
IMPLICIT NONE
!
CONTAINS
!
!=======================================================================
!=======================================================================
!
SUBROUTINE OVRPFA(LDBASE,LUPPER,IBLOCL,IBLOCR,IDBLOL,IDBLOR,IBLQPL,    &
                  IBLQPR,PREPFA,WALEFT,WARIGH,BTAPLF,OVERLP,ICHARG)
!
INTEGER, INTENT(IN) :: LDBASE,LUPPER,IBLOCR,IBLOCL,IDBLOL,IDBLOR,ICHARG
INTEGER, DIMENSION(0:1,1:NDBLOC), INTENT(IN) :: IBLQPL,IBLQPR
COMPLEX, DIMENSION(2), INTENT(IN) :: PREPFA
COMPLEX, DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN),                      &
                                         INTENT(IN) :: WALEFT, WARIGH
COMPLEX, DIMENSION(LUPPER+NDBLOC,LUPPER+NDBLOC),                       &
                                         INTENT(IN) :: BTAPLF
!
COMPLEX, DIMENSION(2), INTENT(OUT) :: OVERLP
!
!=======================================================================
!     THIS SUBROUTINE COMPUTES THE OVERLAP BETWEEN HFB STATES.
!
!     BTAPLF AND BTAPRG ARE PRODUCTS B^TA FOR THE LEFT AND RIGHT STATES,
!     RESPECTIVELY.  THEY ENTER  AS  THE DIAGONAL BLOCKS OF THE PFAFFIAN
!     MATRIX. AS THEY DO NOT DEPEND ON THE SYMMETRY TRANSFORMATION, THEY
!     ARE PRE-COMPUTED IN PRPFAF AND FEED THE PRESENT ROUTINE.
!=======================================================================
!
INTEGER :: IALLOC,ISTATE,NUBLOC,MUBLOC,MDBLOL,IBLOPB,IBLOPA,MDBLOR,    &
           LDPFAF,INFOPF,JSTATE
COMPLEX :: C_ZERO,C_UNIT,SIGNPF
INTEGER, DIMENSION(:), ALLOCATABLE :: IPFWRK
REAL,    DIMENSION(:), ALLOCATABLE :: RPFWRK
COMPLEX, DIMENSION(2) :: PFAFFI
COMPLEX, DIMENSION(:), ALLOCATABLE :: PFAWRK
COMPLEX, DIMENSION(:,:), ALLOCATABLE :: BMATLF,BMATRG,BAPRRG,ABLOLF,ABLORG,   &
                                        BBPROD,BLQRPR,BRQLPR,QLQRPR,          &
                                        PFAMAT,BTAPRG,AMATRG,EAMATR,          &
                                        BBLORG,BPPRRG,QPPRRG
!
!=======================================================================
!      INITIALISATION
!=======================================================================
!
C_ZERO=CMPLX(0.0D0,0.0D0)
C_UNIT=CMPLX(1.0D0,0.0D0)
MDBLOL=LUPPER+IDBLOL
MDBLOR=LUPPER+IDBLOR
!
ALLOCATE (BMATLF(1:2*LDBASE,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('BMATLF','OVRPFA')
ALLOCATE (BMATRG(1:2*LDBASE,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('BMATRG','OVRPFA')
ALLOCATE (AMATRG(1:2*LDBASE,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('AMATRG','OVRPFA')
ALLOCATE (EAMATR(1:LUPPER,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('EAMATR','OVRPFA')
ALLOCATE (BAPRRG(1:LUPPER,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('BAPRRG','OVRPFA')
ALLOCATE (BBPROD(1:LUPPER,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('BBPROD','OVRPFA')
ALLOCATE (BTAPRG(1:MDBLOR,1:MDBLOR),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('BTAPRG','OVRPFA')
!
BMATLF(:,:)=C_ZERO
BMATRG(:,:)=C_ZERO
AMATRG(:,:)=C_ZERO
EAMATR(:,:)=C_ZERO
BBPROD(:,:)=C_ZERO
BTAPRG(:,:)=C_ZERO
BAPRRG(:,:)=C_ZERO
!
IF (IBLOCL.EQ.1) THEN
!
    ALLOCATE (ABLOLF(1:IDBLOL,1:2*LDBASE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('ABLOLF','OVRPFA')
    ALLOCATE (BRQLPR(1:LUPPER,1:IDBLOL),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BRQLPR','OVRPFA')
!
    ABLOLF(:,:)=C_ZERO
    BRQLPR(:,:)=C_ZERO
!
END IF
!
IF (IBLOCR.EQ.1) THEN
!
    ALLOCATE (ABLORG(1:IDBLOR,1:2*LDBASE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('ABLORG','OVRPFA')
    ALLOCATE (BLQRPR(1:LUPPER,1:IDBLOR),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BLQRPR','OVRPFA')
    ALLOCATE (BBLORG(1:IDBLOR,1:2*LDBASE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BBLORG','OVRPFA')
    ALLOCATE (BPPRRG(1:LUPPER,1:IDBLOR),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BPPRRG','OVRPFA')
    ALLOCATE (QPPRRG(1:IDBLOR,1:IDBLOR),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('QPPRRG','OVRPFA')
!
    ABLORG(:,:)=C_ZERO
    BLQRPR(:,:)=C_ZERO
    BBLORG(:,:)=C_ZERO
    BPPRRG(:,:)=C_ZERO
    QPPRRG(:,:)=C_ZERO
!
END IF
!
IF ((IBLOCL.EQ.1).AND.(IBLOCR.EQ.1)) THEN
!
    ALLOCATE (QLQRPR(1:IDBLOL,1:IDBLOR),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('QLQRPR','OVRPFA')
!
    QLQRPR(:,:)=C_ZERO
!
END IF
!
!=======================================================================
!     COLLECTING THE B MATRICES FOR THE LEFT AND RIGHT STATES.
!     FOR A BLOCKED HFB STATE,  THE BLOCKED  QP ARE COLLECTED SEPARATELY
!     IN BBLOLF, ABLOLF, AND THE COMPONENTS ARE FLIPPED BACK IN B AND A.
!=======================================================================
!
DO ISTATE=1,LUPPER
!
    JSTATE=ISTATE+LUPPER
!
    BMATLF(1:LDBASE         ,ISTATE)= CONJG(WALEFT(1:LDBASE,ISTATE,0))
    BMATLF(LDBASE+1:2*LDBASE,ISTATE)= CONJG(WALEFT(1:LDBASE,ISTATE,1))
    BMATRG(1:LDBASE         ,ISTATE)= CONJG(WARIGH(1:LDBASE,ISTATE,0))
    BMATRG(LDBASE+1:2*LDBASE,ISTATE)= CONJG(WARIGH(1:LDBASE,ISTATE,1))
    AMATRG(1:LDBASE         ,ISTATE)= CONJG(WARIGH(1:LDBASE,JSTATE,1))
    AMATRG(LDBASE+1:2*LDBASE,ISTATE)=-CONJG(WARIGH(1:LDBASE,JSTATE,0))
!
END DO
!
IF (IBLOCL.EQ.1) THEN
!
    MUBLOC=0
    DO NUBLOC=1,NDBLOC
       IBLOPB=IBLQPL(0,NUBLOC)
       IBLOPA=IBLQPL(1,NUBLOC)
!
       IF (IBLOPB.GT.0) THEN
           BMATLF(1:LDBASE         ,IBLOPB)= WALEFT(1:LDBASE,IBLOPA,1)
           BMATLF(LDBASE+1:2*LDBASE,IBLOPB)=-WALEFT(1:LDBASE,IBLOPA,0)
!
           MUBLOC=MUBLOC+1
           ABLOLF(MUBLOC,1:LDBASE         )=WALEFT(1:LDBASE,IBLOPB,0)
           ABLOLF(MUBLOC,LDBASE+1:2*LDBASE)=WALEFT(1:LDBASE,IBLOPB,1)
       END IF
!
   END DO
!
END IF
!
IF (IBLOCR.EQ.1) THEN
!
    MUBLOC=0
    DO NUBLOC=1,NDBLOC
       IBLOPB=IBLQPR(0,NUBLOC)
       IBLOPA=IBLQPR(1,NUBLOC)
!
       IF (IBLOPB.GT.0) THEN
           BMATRG(1:LDBASE         ,IBLOPB)= WARIGH(1:LDBASE,IBLOPA,1)
           BMATRG(LDBASE+1:2*LDBASE,IBLOPB)=-WARIGH(1:LDBASE,IBLOPA,0)
           AMATRG(1:LDBASE         ,IBLOPB)= WARIGH(1:LDBASE,IBLOPB,0)
           AMATRG(LDBASE+1:2*LDBASE,IBLOPB)= WARIGH(1:LDBASE,IBLOPB,1)

           MUBLOC=MUBLOC+1
           BBLORG(MUBLOC,1:LDBASE         )= WARIGH(1:LDBASE,IBLOPA,1)
           BBLORG(MUBLOC,LDBASE+1:2*LDBASE)=-WARIGH(1:LDBASE,IBLOPA,0)
           ABLORG(MUBLOC,1:LDBASE         )= WARIGH(1:LDBASE,IBLOPB,0)
           ABLORG(MUBLOC,LDBASE+1:2*LDBASE)= WARIGH(1:LDBASE,IBLOPB,1)
       END IF
!
    END DO
!
    IF(MUBLOC.NE.IDBLOR) STOP 'MUBLOC NE IDBLOR IN PRPFAF'
!
END IF
!
!=======================================================================
!     BUILDING THE PFAFFIAN MATRIX FOR HFB VACUUM
!=======================================================================
!
IF ((IBLOCL.EQ.0).AND.(IBLOCR.EQ.0)) THEN
!
    LDPFAF=2*LUPPER
    ALLOCATE (PFAMAT(1:LDPFAF,1:LDPFAF),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('PFAMAT','OVRPFA')
    PFAMAT(:,:)=C_ZERO
!
!   B^T A (precomputed)
    PFAMAT(1:LUPPER,1:LUPPER)=BTAPLF(1:LUPPER,1:LUPPER)
!
!   B'^+ B
    CALL ZGEMM('C','N',LUPPER,LUPPER,2*LDBASE,C_UNIT,                  &
               BMATRG,2*LDBASE,BMATLF,2*LDBASE,C_ZERO,BBPROD,LUPPER)
    PFAMAT(LUPPER+1:LDPFAF,1:LUPPER)=-BBPROD(1:LUPPER,1:LUPPER)
    PFAMAT(1:LUPPER,LUPPER+1:LDPFAF)= TRANSPOSE(                       &
                                      BBPROD(1:LUPPER,1:LUPPER))
!
!     A^+B^* (computed in PRPFAF before version 43).
!
    CALL ZGEMM('T','N',LUPPER,LUPPER,2*LDBASE,C_UNIT,                  &
             BMATRG,2*LDBASE,AMATRG,2*LDBASE,C_ZERO,BAPRRG,LUPPER)
!
    BTAPRG(1:LUPPER,1:LUPPER)=                                         &
                        TRANSPOSE(CONJG(BAPRRG(1:LUPPER,1:LUPPER)))
!
    BTAPRG=0.5_8*(BTAPRG-TRANSPOSE(BTAPRG))
!
    PFAMAT(LUPPER+1:LDPFAF,LUPPER+1:LDPFAF)=BTAPRG(1:LUPPER,1:LUPPER)
!
!=======================================================================
!     BUILDING THE PFAFFIAN MATRIX FOR STATES WITH (MULTI)QP EXCITATION.
!     NOTE: THE DIMENSION OF THE DIAGONAL BLOCKS AUTOMATICALLY ADJUSTS
!           TO THE CASE SINCE LDBLOX=LUPPER WITOUT BLOCKING (X=L,R).
!=======================================================================
!
ELSE ! IBLOCR=1 AND/OR IBLOCL=1
!
    LDPFAF=MDBLOL+MDBLOR
    ALLOCATE (PFAMAT(1:LDPFAF,1:LDPFAF),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('PFAMAT','OVRPFA')
    PFAMAT(:,:)=C_ZERO
!
!   B^T A / B^T p^+ / -p^* B / q^*p^+ (precomputed)
    PFAMAT(1:MDBLOL,1:MDBLOL)=BTAPLF(1:MDBLOL,1:MDBLOL)
!
!   B'^+ B
    CALL ZGEMM('C','N',LUPPER,LUPPER,2*LDBASE,C_UNIT,                  &
               BMATRG,2*LDBASE,BMATLF,2*LDBASE,C_ZERO,BBPROD,LUPPER)
    PFAMAT(MDBLOL+IDBLOR+1:LDPFAF,1:LUPPER)=-BBPROD(1:LUPPER,1:LUPPER)
    PFAMAT(1:LUPPER,MDBLOL+IDBLOR+1:LDPFAF)=                           &
                                   TRANSPOSE(BBPROD(1:LUPPER,1:LUPPER))
!
!  B^T q'^T
    IF (IBLOCR.EQ.1) THEN
        CALL ZGEMM('T','T',LUPPER,IDBLOR,2*LDBASE,C_UNIT,              &
                   BMATLF,2*LDBASE,ABLORG,IDBLOR,C_ZERO,BLQRPR,LUPPER)
        PFAMAT(1:LUPPER,MDBLOL+1:MDBLOL+IDBLOR)=                       &
                                              BLQRPR(1:LUPPER,1:IDBLOR)
        PFAMAT(MDBLOL+1:MDBLOL+IDBLOR,1:LUPPER)=                       &
                                   -TRANSPOSE(BLQRPR(1:LUPPER,1:IDBLOR))
    END IF
!
!   B'^+ q^+
    IF (IBLOCL.EQ.1) THEN
        CALL ZGEMM('C','C',LUPPER,IDBLOL,2*LDBASE,C_UNIT,              &
                   BMATRG,2*LDBASE,ABLOLF,IDBLOL,C_ZERO,BRQLPR,LUPPER)
        PFAMAT(MDBLOL+IDBLOR+1:LDPFAF,LUPPER+1:MDBLOL)=                &
                                             -BRQLPR(1:LUPPER,1:IDBLOL)
        PFAMAT(LUPPER+1:MDBLOL,MDBLOL+IDBLOR+1:LDPFAF)=                &
                                    TRANSPOSE(BRQLPR(1:LUPPER,1:IDBLOL))
    END IF
!
!   q q'^+
    IF ((IBLOCL.EQ.1).AND.(IBLOCR.EQ.1)) THEN
        CALL ZGEMM('N','C',IDBLOL,IDBLOR,2*LDBASE,C_UNIT,              &
                   ABLOLF,IDBLOL,ABLORG,IDBLOR,C_ZERO,QLQRPR,IDBLOL)
        PFAMAT(LUPPER+1:LUPPER+IDBLOL,MDBLOL+1:MDBLOL+IDBLOR)=         &
                                      CONJG(QLQRPR(1:IDBLOL,1:IDBLOR))
        PFAMAT(MDBLOL+1:MDBLOL+IDBLOR,LUPPER+1:LUPPER+IDBLOL)=         &
                           -TRANSPOSE(CONJG(QLQRPR(1:IDBLOL,1:IDBLOR)))
    END IF
!
!   A'^+ B'^* / p'B'^* / -B'^+p'^T / p'q'^T (precomputed)
    CALL ZGEMM('T','N',LUPPER,LUPPER,2*LDBASE,C_UNIT,                  &
             BMATRG,2*LDBASE,AMATRG,2*LDBASE,C_ZERO,BAPRRG,LUPPER)
!
    BTAPRG(IDBLOR+1:MDBLOR,IDBLOR+1:MDBLOR)=TRANSPOSE(CONJG(           &
                                           BAPRRG(1:LUPPER,1:LUPPER)))
    CALL ZGEMM('T','C',LUPPER,IDBLOR,2*LDBASE,C_UNIT,                  &
               BMATRG,2*LDBASE,BBLORG,IDBLOR,C_ZERO,BPPRRG,LUPPER)
!
    BTAPRG(IDBLOR+1:MDBLOR,1:IDBLOR)=                                  &
                                 -CONJG(BPPRRG(1:LUPPER,1:IDBLOR))
    BTAPRG(1:IDBLOR,IDBLOR+1:MDBLOR)=                                  &
                        TRANSPOSE(CONJG(BPPRRG(1:LUPPER,1:IDBLOR)))

    CALL ZGEMM('N','T',IDBLOR,IDBLOR,2*LDBASE,C_UNIT,                  &
               ABLORG,IDBLOR,BBLORG,IDBLOR,C_ZERO,QPPRRG,IDBLOR)
!
    BTAPRG(1:IDBLOR,1:IDBLOR)=TRANSPOSE(QPPRRG(1:IDBLOR,1:IDBLOR))

    BTAPRG=0.5_8*(BTAPRG-TRANSPOSE(BTAPRG))
!
    PFAMAT(MDBLOL+1:LDPFAF,MDBLOL+1:LDPFAF)= BTAPRG(1:MDBLOR,1:MDBLOR)
!
!   Deallocating
    IF (IBLOCL.EQ.1) DEALLOCATE (ABLOLF,BRQLPR)
!
    IF (IBLOCR.EQ.1) DEALLOCATE (ABLORG,BBLORG,BLQRPR,BPPRRG,QPPRRG)
!
    IF ((IBLOCL.EQ.1).AND.(IBLOCR.EQ.1)) DEALLOCATE (QLQRPR)
!
END IF
!
DEALLOCATE (BMATLF,BMATRG,BBPROD,BTAPRG,AMATRG,EAMATR)
!
!=======================================================================
!     COMPUTING THE NORM/OVERLAP
!=======================================================================
!
LDPFAF=SIZE(PFAMAT,1)
ALLOCATE (IPFWRK(LDPFAF),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('IPFWRK','OVRPFA')
ALLOCATE (PFAWRK(LDPFAF),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('PFAWRK','OVRPFA')
ALLOCATE (RPFWRK(LDPFAF),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('RPFWRK','OVRPFA')
!
CALL SKPF10('U','P',LDPFAF,PFAMAT,LDPFAF,PFAFFI,IPFWRK,PFAWRK,LDPFAF,  &
                    RPFWRK,INFOPF)
!
IF(INFOPF.NE.0) THEN
    WRITE(*,*) 'INFO PFAFFIAN =',INFOPF
    STOP
END IF
!
SIGNPF=(-C_UNIT)**(LUPPER*(LUPPER-1)/2) * (-C_UNIT)**(IDBLOL*(IDBLOL-1)/2)
OVERLP(1)=SIGNPF*PREPFA(1)*PFAFFI(1)
OVERLP(2)=PREPFA(2)+PFAFFI(2)
!
IF (AIMAG(OVERLP(2)).GT.1.0D-9) STOP 'OVERLP(2) CANNOT BE IMAGINARY'
!
DEALLOCATE (PFAMAT,IPFWRK,PFAWRK,RPFWRK)
!
!=======================================================================
!
END SUBROUTINE OVRPFA
!
!=======================================================================
!=======================================================================
!
SUBROUTINE PRPFAF(LDBASE,LUPPER,IBLOCR,IDBLOR,IBLQPR,WARIGH,BTAPLF,    &
                                                            BTAPRG)
!
INTEGER, INTENT(IN) :: LDBASE,LUPPER,IBLOCR,IDBLOR
INTEGER, DIMENSION(0:1,1:NDBLOC), INTENT(IN) :: IBLQPR
COMPLEX, DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN), INTENT(IN) :: WARIGH
!
COMPLEX, DIMENSION(LUPPER+NDBLOC,LUPPER+NDBLOC),                       &
                                         INTENT(OUT) :: BTAPLF,BTAPRG
!
!=======================================================================
!     PREPARE THE COMPUTATION OF THE OVERLAP BETWEEN HFB STATES  BY PRE-
!     COMPUTING QUANTITIES THAT DO NOT DEPEND ON THE  SYMMETRY OPERATOR.
!     THE ROUTINE COMPUTES B^TA=BTAPLF AND A^+B^*=BTAPRG FOR THE LEFT OR
!     RIGHT STATES, (IDIAGO=1), I.E. THE DIAGONAL BLOCKS OF THE PFAFFIAN
!     MATRIX.
!     THE CASE IDIAGO=0 IS TREATED VIA TWO SEPARATE CALLS OF PRPFAF, THE
!     SECOND ONE WITH WARIGH=WALEFT AND VARIABLE IXXXXR=IXXXXL.
!=======================================================================
!
INTEGER :: IALLOC,ISTATE,JSTATE,NUBLOC,MUBLOC,IBLOPB,IBLOPA,MDBLOR
COMPLEX :: C_ZERO,C_UNIT
COMPLEX, DIMENSION(:,:), ALLOCATABLE :: AMATRG,BMATRG,BAPRRG,ABLORG,   &
                                        BBLORG,BPPRRG,QPPRRG
!
!=======================================================================
!     INITIALISATION
!=======================================================================
!
C_ZERO=CMPLX(0.0D0,0.0D0)
C_UNIT=CMPLX(1.0D0,0.0D0)
BTAPLF(:,:)=C_ZERO
BTAPRG(:,:)=C_ZERO
!
!=======================================================================
!     ALLOCATION
!=======================================================================
!
MDBLOR=LUPPER+IDBLOR
!
ALLOCATE (AMATRG(1:2*LDBASE,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('AMATRG','PRPFAF')
ALLOCATE (BMATRG(1:2*LDBASE,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('BMATRG','PRPFAF')
ALLOCATE (BAPRRG(1:LUPPER,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('BAPRRG','PRPFAF')
!
AMATRG(:,:)=C_ZERO
BMATRG(:,:)=C_ZERO
BAPRRG(:,:)=C_ZERO
!
IF (IBLOCR.EQ.1) THEN
!
    ALLOCATE (ABLORG(1:IDBLOR,1:2*LDBASE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('ABLORG','PRPFAF')
    ALLOCATE (BBLORG(1:IDBLOR,1:2*LDBASE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BBLORG','PRPFAF')
    ALLOCATE (BPPRRG(1:LUPPER,1:IDBLOR),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BPPRRG','PRPFAF')
    ALLOCATE (QPPRRG(1:IDBLOR,1:IDBLOR),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('QPPRRG','PRPFAF')
!
    ABLORG(:,:)=C_ZERO
    BBLORG(:,:)=C_ZERO
    BPPRRG(:,:)=C_ZERO
    QPPRRG(:,:)=C_ZERO
!
END IF
!
!=======================================================================
!     COLLECTING THE  B AND A MATRICES FOR THE RIGHT STATE, AND REMOVING
!     THE RUSSIAN CONVENTION.  FOR A BLOCKED HFB STATE,  THE BLOCKED  QP
!     ARE COLLECTED SEPARATELY IN  BBLORG AND ABLORG  AND THE COMPONENTS
!     ARE FLIPPED BACK IN B AND A.
!=======================================================================
!
DO ISTATE=1,LUPPER
    JSTATE=LUPPER+ISTATE
!
    BMATRG(1:LDBASE         ,ISTATE)= CONJG(WARIGH(1:LDBASE,ISTATE,0))
    BMATRG(LDBASE+1:2*LDBASE,ISTATE)= CONJG(WARIGH(1:LDBASE,ISTATE,1))
    AMATRG(1:LDBASE         ,ISTATE)= CONJG(WARIGH(1:LDBASE,JSTATE,1))
    AMATRG(LDBASE+1:2*LDBASE,ISTATE)=-CONJG(WARIGH(1:LDBASE,JSTATE,0))
!
END DO
!
IF (IBLOCR.EQ.1) THEN
!
    MUBLOC=0
    DO NUBLOC=1,NDBLOC
       IBLOPB=IBLQPR(0,NUBLOC)
       IBLOPA=IBLQPR(1,NUBLOC)
!
       IF (IBLOPB.GT.0) THEN
           MUBLOC=MUBLOC+1

           BMATRG(1:LDBASE         ,IBLOPB)= WARIGH(1:LDBASE,IBLOPA,1)
           BMATRG(LDBASE+1:2*LDBASE,IBLOPB)=-WARIGH(1:LDBASE,IBLOPA,0)
           AMATRG(1:LDBASE         ,IBLOPB)= WARIGH(1:LDBASE,IBLOPB,0)
           AMATRG(LDBASE+1:2*LDBASE,IBLOPB)= WARIGH(1:LDBASE,IBLOPB,1)

           BBLORG(MUBLOC,1:LDBASE         )= WARIGH(1:LDBASE,IBLOPA,1)
           BBLORG(MUBLOC,LDBASE+1:2*LDBASE)=-WARIGH(1:LDBASE,IBLOPA,0)
           ABLORG(MUBLOC,1:LDBASE         )= WARIGH(1:LDBASE,IBLOPB,0)
           ABLORG(MUBLOC,LDBASE+1:2*LDBASE)= WARIGH(1:LDBASE,IBLOPB,1)
       END IF
!
    END DO
!
    IF(MUBLOC.NE.IDBLOR) STOP 'MUBLOC NE IDBLOR IN PRPFAF'
!
END IF
!
!=======================================================================
!     THE PRODUCT B^T A IS COMPUTED AND STORED IN BTAPLF.
!=======================================================================
!
CALL ZGEMM('T','N',LUPPER,LUPPER,2*LDBASE,C_UNIT,                      &
           BMATRG,2*LDBASE,AMATRG,2*LDBASE,C_ZERO,BAPRRG,LUPPER)
!
BTAPLF(1:LUPPER,1:LUPPER)=BAPRRG(1:LUPPER,1:LUPPER)
!
!=======================================================================
!     THE EXTRA ROWS AND COLUMNS  ENTERING THE PFAFFIAN  MATRIX  IN THE
!     CASE OF AN EXCITED QP STATE ARE COMPUTED.
!=======================================================================
!
IF (IBLOCR.EQ.1) THEN
!
! B^T p^+
    CALL ZGEMM('T','C',LUPPER,IDBLOR,2*LDBASE,C_UNIT,                  &
               BMATRG,2*LDBASE,BBLORG,IDBLOR,C_ZERO,BPPRRG,LUPPER)
!
    BTAPLF(1:LUPPER,LUPPER+1:MDBLOR)=           BPPRRG(1:LUPPER,1:IDBLOR)
    BTAPLF(LUPPER+1:MDBLOR,1:LUPPER)=-TRANSPOSE(BPPRRG(1:LUPPER,1:IDBLOR))
!
! q^*p^+
    CALL ZGEMM('N','T',IDBLOR,IDBLOR,2*LDBASE,C_UNIT,                  &
               ABLORG,IDBLOR,BBLORG,IDBLOR,C_ZERO,QPPRRG,IDBLOR)
!
    BTAPLF(LUPPER+1:MDBLOR,LUPPER+1:MDBLOR)=                           &
                                       CONJG(QPPRRG(1:IDBLOR,1:IDBLOR))
!
END IF
!
!=======================================================================
!     A^+B^*=BTAPRG, BUILD USING THE RESULTS FOR B^TA=BTAPLF (IDIAGO=1).
!=======================================================================
!
IF (IBLOCR.EQ.0) THEN
!
    BTAPRG(1:LUPPER,1:LUPPER)=                                         &
                        TRANSPOSE(CONJG(BAPRRG(1:LUPPER,1:LUPPER)))
ELSE
    BTAPRG(IDBLOR+1:MDBLOR,IDBLOR+1:MDBLOR)=                           &
                        TRANSPOSE(CONJG(BAPRRG(1:LUPPER,1:LUPPER)))
    BTAPRG(IDBLOR+1:MDBLOR,1:IDBLOR)=                                  &
                                 -CONJG(BPPRRG(1:LUPPER,1:IDBLOR))
    BTAPRG(1:IDBLOR,IDBLOR+1:MDBLOR)=                                  &
                        TRANSPOSE(CONJG(BPPRRG(1:LUPPER,1:IDBLOR)))
    BTAPRG(1:IDBLOR,1:IDBLOR)=TRANSPOSE(QPPRRG(1:IDBLOR,1:IDBLOR))
!
    DEALLOCATE (ABLORG,BBLORG,BPPRRG,QPPRRG)
!
END IF
!
DEALLOCATE (AMATRG,BMATRG,BAPRRG)
!
!=======================================================================
!          ENFORCING THE SKEW-SYMMETRY FOR PRECISION
!=======================================================================
!
BTAPLF=0.5_8*(BTAPLF-TRANSPOSE(BTAPLF))
BTAPRG=0.5_8*(BTAPRG-TRANSPOSE(BTAPRG))
!
!=======================================================================
!
END SUBROUTINE PRPFAF
!
!=======================================================================
!=======================================================================
!
SUBROUTINE OVRONI(LDBASE,LUPPER,PREPFA,DETOVR,WALEFT,WARIGH,OVERLP)
!
INTEGER, INTENT(IN) :: LDBASE,LUPPER
COMPLEX, DIMENSION(2), INTENT(IN) :: DETOVR,PREPFA
COMPLEX, DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN),                      &
                                         INTENT(IN) :: WALEFT, WARIGH
!
COMPLEX, DIMENSION(2), INTENT(OUT) :: OVERLP
!
!=======================================================================
!     THIS SUBROUTINE COMPUTES THE  OVERLAP BETWEEN HFB STATES USING THE
!     ONISHI FORMULA:
!
!       <WALEFT|WARIGH> = N*S*SQRT(DET[(B')^T B^*])*SQRT(DET(OVRLAP))
!
!     WHERE S = (-1)^(LUPPER*(LUPPER-1)/2)  IS A PHASE FACTOR, AND N THE
!     NORMALIZATION PRECOMPUTED BY CALLING THE SUBROUTINE WITH N = 1 AND
!     STORED IN PREPFA. DET(OVRLAP) IS STORED IN DETOVR.
!=======================================================================
!
INTEGER :: IALLOC,ISTATE,JSTATE,IBASE,ISPIN
REAL    :: RICOND,DPOWER
COMPLEX :: C_ZERO,C_UNIT,UNIT_I,ONIMOD,ONIPHS,SIGNPF
COMPLEX, DIMENSION(2) :: DETWRK,ONIPRE,ONIOVR
INTEGER, DIMENSION(:), ALLOCATABLE :: IAUXPR
COMPLEX, DIMENSION(:), ALLOCATABLE :: AUXDPR
COMPLEX, DIMENSION(:,:), ALLOCATABLE :: ONIMAT
!
!=======================================================================
!      INITIALISATION
!=======================================================================
!
C_ZERO=CMPLX(0.0D0,0.0D0)
C_UNIT=CMPLX(1.0D0,0.0D0)
UNIT_I=CMPLX(1.0D0,1.0D0)
!
ALLOCATE (ONIMAT(1:LUPPER,1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('ONIMAT','OVRONI')
!
ONIMAT(:,:)=C_ZERO
!
ALLOCATE (IAUXPR(1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('IAUXPR','OVRONI')
ALLOCATE (AUXDPR(1:LUPPER),STAT=IALLOC)
IF (IALLOC.NE.0) CALL NOALLO('AUXDPR','OVRONI')
!
!=======================================================================
!      COMPUTING THE PREFACTOR : SQRT(DET[(B')^T B^*])
!=======================================================================
!
DO ISTATE=1,LUPPER
   DO JSTATE=1,LUPPER
!
      DO IBASE=1,LDBASE
         DO ISPIN=0,1
!
            ONIMAT(ISTATE,JSTATE)=ONIMAT(ISTATE,JSTATE)                &
                                 +CONJG(WALEFT(IBASE,ISTATE,ISPIN))    &
                                       *WARIGH(IBASE,JSTATE,ISPIN)
!
         END DO
      END DO
!
   END DO
END DO
!
CALL ZGECO(ONIMAT,LUPPER,LUPPER,IAUXPR,RICOND,AUXDPR)
CALL ZGEDI(ONIMAT,LUPPER,LUPPER,IAUXPR,DETWRK,AUXDPR,10)
!
IF (ABS(RICOND).LT.1.0D-12)                                            &
   WRITE(*,*)  'BE CAREFUL - THE PRE-ONISHI MATRIX MAY BE SINGULAR'
!
!DPOWER=REAL(DETWRK(2))
!IF (ABS(MOD(DPOWER,2.0D0)).GT.1.0D-9) THEN
!    IF (DPOWER.LT.0.0D0) THEN
!       DETWRK(2)=DETWRK(2)+C_UNIT
!       DETWRK(1)=DETWRK(1)*0.1D0
!    ELSE
!       DETWRK(2)=DETWRK(2)-C_UNIT
!       DETWRK(1)=DETWRK(1)*10.0D0
!    END IF
!END IF
!
!ONIPRE(1)=DETWRK(1)
!ONIMOD=ABS(ONIPRE(1))
!ONIPHS=ATAN2(AIMAG(ONIPRE(1)),REAL(ONIPRE(1)))
!ONIPRE(1)=SQRT(ONIMOD)*EXP(UNIT_I*ONIPHS/2.0D0)
!ONIPRE(2)=DETWRK(2)/2.0D0
!
ONIPRE(1)=SQRT(DETWRK(1))
ONIPRE(2)=DETWRK(2)/2.0D0
!
DEALLOCATE (ONIMAT,IAUXPR,AUXDPR)
!
!=======================================================================
!      COMPUTING THE MAIN FACTOR SQRT(DET(OVRLAP)) THE DETERMINANT BEING
!      AN INPUT OF THE ROUTINE
!=======================================================================
!
DETWRK(:)=DETOVR(:)
!
!DPOWER=REAL(DETOVR(2))
!IF (ABS(MOD(DPOWER,2.0D0)).GT.1.0D-9) THEN
!    IF (DPOWER.LT.0.0D0) THEN
!        DETWRK(2)=DETOVR(2)+C_UNIT
!        DETWRK(1)=DETOVR(1)*0.1D0
!    ELSE
!        DETWRK(2)=DETOVR(2)-C_UNIT
!        DETWRK(1)=DETOVR(1)*10.0D0
!    END IF
!END IF
!
!ONIOVR(1)=DETWRK(1)
!ONIMOD=ABS(ONIOVR(1))
!ONIPHS=ATAN2(AIMAG(ONIOVR(1)),REAL(ONIOVR(1)))
!ONIOVR(1)=SQRT(ONIMOD)*EXP(UNIT_I*ONIPHS/2.0D0)
!ONIOVR(2)=DETWRK(2)/2.0D0
ONIOVR(1)=SQRT(DETWRK(1))
ONIOVR(2)=DETWRK(2)/2.0D0
!
!=======================================================================
!      COMPUTING THE NORM/OVERLAP
!=======================================================================
!
!SIGNPF=(-C_UNIT)**(LUPPER*(LUPPER-1)/2)
!
OVERLP(1)=PREPFA(1)*ONIPRE(1)*ONIOVR(1)
OVERLP(2)=PREPFA(2)+ONIPRE(2)+ONIOVR(2)
!
IF (AIMAG(OVERLP(2)).GT.1.0D-9) STOP 'OVERLP(2) CANNOT BE IMAGINARY'
!
END SUBROUTINE OVRONI
!
!=======================================================================
!=======================================================================
!
SUBROUTINE OVRMAT(LDBASE,LUPPER,WALEFT,WARIGH,OVRLAP)
!
INTEGER, INTENT(IN) :: LDBASE,LUPPER
COMPLEX, DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN),                      &
                                         INTENT(IN) :: WALEFT, WARIGH
!
COMPLEX, DIMENSION(1:4*NDSTAT,1:4*NDSTAT), INTENT(OUT) :: OVRLAP
!
!=======================================================================
!     THIS SUBROUTINE COMPUTES THE MATRIX OF OVERLAPS BETWEEN THE QUASI-
!     PARTICLES OF LEFT AND RIGHT HFB STATES:
!            (A')^T A^* + (B')^T B^*
!=======================================================================
!
INTEGER :: ISTATE,ISTALO,JSTATE,JSTALO,ISPIN,IBASE
!
DO ISTATE=1,LUPPER
   ISTALO=ISTATE+LUPPER
!
   DO JSTATE=1,LUPPER
      JSTALO=JSTATE+LUPPER
!
      DO IBASE=1,LDBASE
         DO ISPIN=0,1
!
            OVRLAP(ISTATE,JSTATE)=OVRLAP(ISTATE,JSTATE)                &
                                 +CONJG(WALEFT(IBASE,ISTATE,ISPIN))    &
                                       *WARIGH(IBASE,JSTATE,ISPIN)     &
                                 +CONJG(WALEFT(IBASE,ISTALO,ISPIN))    &
                                       *WARIGH(IBASE,JSTALO,ISPIN)
!
         END DO
      END DO
!
   END DO
!
END DO
!
END SUBROUTINE OVRMAT
!
!=======================================================================
!=======================================================================
!
END MODULE OVRHFB
!
!=======================================================================
!=======================================================================
!=======================================================================
!=======================================================================
!=======================================================================
!=======================================================================
!
SUBROUTINE INIADD
!
!=======================================================================
!     MANAGE THE ALLOCATION OF THE ARRAYS USED FOR THE BLOCKING.
!=======================================================================
!
USE BLOSAV
USE hfodd_sizes
IMPLICIT NONE
INTEGER::IALLOC
COMPLEX::C_ZERO
!=======================================================================
!ALLOCATE SPACE TO ARRAY
!=======================================================================
!
 IALLOC=0
 C_ZERO=CMPLX(0.0D0,0.0D0)
!
 IF (.NOT.ALLOCATED(BLOLAB)) THEN
    ALLOCATE (BLOLAB(0:NDISOS,1:NDBLOC),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BLOLAB','INIADD')
    BLOLAB(:,:)%IR=0
    BLOLAB(:,:)%N=0
 END IF
!
 IF (.NOT.ALLOCATED(WLABBR)) THEN
    ALLOCATE (WLABBR(0:NDISOS,1:NDBLOC),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('WLABBR','INIADD')
    WLABBR(:,:)=-1
 END IF
!
 IF (.NOT.ALLOCATED(WLABAR)) THEN
    ALLOCATE (WLABAR(0:NDISOS,1:NDBLOC),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('WLABAR','INIADD')
    WLABAR(:,:)=-1
 END IF

 IF (.NOT.ALLOCATED(WLABBL)) THEN
    ALLOCATE (WLABBL(0:NDISOS,1:NDBLOC),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('WLABBL','INIADD')
    WLABBL(:,:)=-1
 END IF
!
 IF (.NOT.ALLOCATED(WLABAL)) THEN
    ALLOCATE (WLABAL(0:NDISOS,1:NDBLOC),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('WLABAL','INIADD')
    WLABAL(:,:)=-1
 END IF

!
 IF (.NOT.ALLOCATED(LPRODR)) THEN
    ALLOCATE (LPRODR(0:NDISOS),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('LPRODR','INIADD')
    LPRODR(:)=0
 END IF
!
 IF (.NOT.ALLOCATED(LPRODL)) THEN
    ALLOCATE (LPRODL(0:NDISOS),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('LPRODL','INIADD')
    LPRODL(:)=0
 END IF
!
 IF (.NOT.ALLOCATED(LDBLOL)) THEN
    ALLOCATE (LDBLOL(0:NDISOS),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('LDBLOL','INIADD')
    LPRODL(:)=0
 END IF
!
 RETURN
END SUBROUTINE INIADD
!
!
!=======================================================================
!=======================================================================
!=======================================================================
!=======================================================================
!=======================================================================
!=======================================================================
!
MODULE O2AVRG
!
!=======================================================================
!      MODULE COMPUTING <o^2> FOR ONE-BODY OPERATORS O VARIOUS FORMS
!      THE ROUTINES AVRXSD, AVRXSO, EXCXSD, EXCXSO, PAIXSD, PAIXSO ARE
!      USED IN LIPCORR ONLY
!=======================================================================
!
CONTAINS
!=======================================================================
SUBROUTINE AVRXSD(HAUXIL,AOBSLI,AKBSLI,RHO_PP)
!=======================================================================
!Calculate Tr[x.rho], x=[x(0),0;0,x(1)]
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 REAL   ::AOBSLI
 COMPLEX::AKBSLI
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::HAUXIL,RHO_PP
!
! INTEGER::INCX,INCY
! COMPLEX::RES1,RES2,ZDOTU
 INTEGER::IBRA,IKET
!
 INTEGER::LDBASE
 COMPLEX::RESULT,C_ZERO
!=======================================================================
 COMMON /DIMENS/ LDBASE
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 RESULT=C_ZERO
!
 DO IBRA=1,LDBASE
    DO IKET=1,LDBASE
       RESULT=RESULT+HAUXIL(IBRA,IKET,0)*RHO_PP(IKET,IBRA,0)&
                    +HAUXIL(IBRA,IKET,1)*RHO_PP(IKET,IBRA,1)

    END DO
 END DO
 AOBSLI=REAL(RESULT)
 AKBSLI=RESULT
 !INCX=1; INCY=1
 !RES1=ZDOTU(NDBASE*NDBASE,HAUXIL(:,:,0),INCX,TRANSPOSE(RHO_PP(:,:,0)),INCY)
 !RES2=ZDOTU(NDBASE*NDBASE,HAUXIL(:,:,1),INCX,TRANSPOSE(RHO_PP(:,:,0)),INCY)
 !AOBSLI=REAL(RES1+RES2)
 !AKBSLI=RES1+RES2
!=======================================================================
 RETURN
END SUBROUTINE AVRXSD
!=======================================================================
SUBROUTINE AVRXSO(HAUXIL,AOBSLI,AKBSLI,RHO_PM)
!=======================================================================
! Calculate Tr[x.rho], x=[0,x(0);x(1),0]
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 REAL   ::AOBSLI
 COMPLEX::AKBSLI
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::HAUXIL,RHO_PM
!
! INTEGER::INCX,INCY
! COMPLEX::RES1,RES2,ZDOTU
 INTEGER::IBRA,IKET
 COMPLEX::RESULT,C_ZERO
!
 INTEGER::LDBASE
 COMMON /DIMENS/ LDBASE
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 RESULT=C_ZERO
!
 DO IBRA=1,LDBASE
    DO IKET=1,LDBASE
       RESULT=RESULT+HAUXIL(IBRA,IKET,0)*RHO_PM(IKET,IBRA,1)&
                    +HAUXIL(IBRA,IKET,1)*RHO_PM(IKET,IBRA,0)
    END DO
 END DO
 AOBSLI=REAL(RESULT)
 AKBSLI=RESULT
 !INCX=1; INCY=1
 !RES1=ZDOTU(NDBASE*NDBASE,HAUXIL(:,:,0),INCX,TRANSPOSE(RHO_PM(:,:,1)),INCY)
 !RES2=ZDOTU(NDBASE*NDBASE,HAUXIL(:,:,1),INCX,TRANSPOSE(RHO_PM(:,:,0)),INCY)
 !AOBSLI=REAL(RES1+RES2)
 !AKBSLI=RES1+RES2
!=======================================================================
 RETURN
END SUBROUTINE AVRXSO
!=======================================================================
SUBROUTINE EXCXSD(IMSIGN,MREVER,ISIMPY,JETKET,EOBSSQ,EKBSSQ,XMAT,RPP,RPM,APP,APM)
!=======================================================================
!Tr[(x.rho)^2], x=[x(0),0;0,x(1)]
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::IMSIGN,MREVER,ISIMPY,JETKET
 REAL   ::EOBSSQ
 COMPLEX::EKBSSQ
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::XMAT,RPP
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE),OPTIONAL::RPM,APP,APM
!
 INTEGER::IREVER
 COMPLEX::C_ZERO,RESULT,C_UNIT,CTWO
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::AUXILI
!
! INTEGER::INCX,INCY
! COMPLEX::RES1,RES2,ZDOTU
 INTEGER::IBRA,IKET
!
 INTEGER::LDBASE
 COMMON /DIMENS/ LDBASE
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 IF(JETKET.EQ.2) THEN
                    APP(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
    IF(ISIMPY.EQ.0) APM(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
 END IF
!=======================================================================
!X.RHO
!=======================================================================
 C_UNIT=CMPLX(1.0D0,0.0D0)
 DO IREVER=0,MREVER
    AUXILI(:,:,IREVER) = C_ZERO
    CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                       C_UNIT,XMAT(:,:,IREVER),NDBASE, &
                       RPP(:,:,IREVER),NDBASE,C_ZERO,AUXILI(:,:,IREVER),NDBASE)
 END DO
!=======================================================================
!TR[(X.RHO)^2] AND 2(X.RHO.X)_{IJ}
!======================================================================
 RESULT=C_ZERO
! INCX=1; INCY=1
 DO IREVER=0,MREVER
!    RESULT=RESULT+ZDOTU(NDBASE*NDBASE,AUXILI(:,:,IREVER),INCX,TRANSPOSE(AUXILI(:,:,IREVER)),INCY)
    DO IBRA=1,LDBASE
       DO IKET=1,LDBASE
          RESULT=RESULT+AUXILI(IBRA,IKET,IREVER)*AUXILI(IKET,IBRA,IREVER)
       END DO
    END DO
    IF (JETKET.EQ.2) THEN
        CTWO=CMPLX(-2.0D0,0.0D0)
        CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                           CTWO,AUXILI(:,:,IREVER),NDBASE, &
                           XMAT(:,:,IREVER),NDBASE,C_ZERO,APP(:,:,IREVER),NDBASE)
    END IF
 END DO
 !
 IF (MREVER.EQ.0) THEN
    RESULT=RESULT+CONJG(RESULT)*(IMSIGN**2)
    IF(JETKET.EQ.2) APP(:,:,1)=CONJG(APP(:,:,0))*(IMSIGN**2)
 END IF
 !
 EOBSSQ = REAL(RESULT)
 EKBSSQ = RESULT
 !
!=======================================================================
!ADDITIONAL TERMS WHEN SIMPLEX IS NOT CONSERVED
!=======================================================================
 IF (ISIMPY.NE.1) THEN
!=======================================================================
!X.RHO
!=======================================================================
    C_UNIT=CMPLX(1.0D0,0.0D0); AUXILI = C_ZERO
    CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                       C_UNIT,XMAT(:,:,0),NDBASE, &
                       RPM(:,:,0),NDBASE,C_ZERO,AUXILI(:,:,0),NDBASE)
    CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                       C_UNIT,XMAT(:,:,1),NDBASE, &
                       RPM(:,:,1),NDBASE,C_ZERO,AUXILI(:,:,1),NDBASE)
!=======================================================================
!TR[(X.RHO)^2] AND (X.RHO.X)_{IJ}
!=======================================================================
!    INCX=1; INCY=1
!    RESULT=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,0),INCX,TRANSPOSE(AUXILI(:,:,1)),INCY)
    RESULT=C_ZERO
    DO IBRA=1,LDBASE
       DO IKET=1,LDBASE
          RESULT=RESULT+AUXILI(IBRA,IKET,0)*AUXILI(IKET,IBRA,1)
       END DO
    END DO
    IF (JETKET.EQ.2) THEN
        IF (MREVER.EQ.0) STOP ' WRONG MREVER IN EXCXSD'
        CTWO=CMPLX(-2.0D0,0.0D0)
        CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                           CTWO,AUXILI(:,:,0),NDBASE, &
                           XMAT(:,:,1),NDBASE,C_ZERO,APM(:,:,0),NDBASE)
        CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                           CTWO,AUXILI(:,:,1),NDBASE, &
                           XMAT(:,:,0),NDBASE,C_ZERO,APM(:,:,1),NDBASE)
    END IF
    !
    RESULT=2*RESULT
    !
    EOBSSQ = EOBSSQ + REAL(RESULT)
    EKBSSQ = EKBSSQ + RESULT
    !
 END IF
!======================================================================
 RETURN
 END SUBROUTINE EXCXSD
!=======================================================================
SUBROUTINE EXCXSO(IMSIGN,MREVER,ISIMPY,JETKET,EOBSSQ,EKBSSQ,XMAT,RPP,RPM,APP,APM)
!=======================================================================
! Tr[(x.rho)^2], x=[0,x(0);x(1),0]
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::IMSIGN,MREVER,ISIMPY,JETKET
 REAL   ::EOBSSQ
 COMPLEX::EKBSSQ
 COMPLEX::C_ZERO,RESULT,C_UNIT,CTWO
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::XMAT,RPP
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE),OPTIONAL::RPM,APP,APM
!
 INTEGER::IREVER
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::AUXILI
!
! INTEGER::INCX,INCY
! COMPLEX::RES1,RES2,ZDOTU
 INTEGER::IBRA,IKET
!
 INTEGER::LDBASE
 COMMON /DIMENS/ LDBASE
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
  IF(JETKET.EQ.2)THEN
                    APP(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
    IF(ISIMPY.EQ.0) APM(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
 END IF
!=======================================================================
!X.RHO
!=======================================================================
 C_UNIT=CMPLX(1.0D0,0.0D0); AUXILI = C_ZERO
 DO IREVER=0,MREVER
    CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                       C_UNIT,XMAT(:,:,1-IREVER),NDBASE, &
                       RPP(:,:,IREVER),NDBASE,C_ZERO,AUXILI(:,:,IREVER),NDBASE)
 END DO
 IF (MREVER.EQ.0) THEN
    AUXILI(:,:,1)=CONJG(AUXILI(:,:,0))*IMSIGN
 END IF
!=======================================================================
!TR[(X.RHO)^2] AND 2(X.RHO.X)
!=======================================================================
! INCX=1; INCY=1
! RESULT=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,1),INCX,TRANSPOSE(AUXILI(:,:,0)),INCY)
    RESULT=C_ZERO
    DO IBRA=1,LDBASE
       DO IKET=1,LDBASE
          RESULT=RESULT+AUXILI(IBRA,IKET,1)*AUXILI(IKET,IBRA,0)
       END DO
    END DO
 IF (JETKET.EQ.2) THEN
     CTWO=CMPLX(-2.0D0,0.0D0)
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        CTWO,AUXILI(:,:,1),NDBASE, &
                        XMAT(:,:,1),NDBASE,C_ZERO,APP(:,:,0),NDBASE)
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        CTWO,AUXILI(:,:,0),NDBASE, &
                        XMAT(:,:,0),NDBASE,C_ZERO,APP(:,:,1),NDBASE)
 END IF
 !
 EOBSSQ = 2*REAL(RESULT)
 EKBSSQ = 2*RESULT
 !
!=======================================================================
!ADDITIONAL TERMS WHEN SIMPLEX IS NOT CONSERVED
!=======================================================================
 IF (ISIMPY.NE.1) THEN
!=======================================================================
!X.RHO
!=======================================================================
     C_UNIT=CMPLX(1.0D0,0.0D0); AUXILI = C_ZERO
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,1),NDBASE, &
                        RPM(:,:,0),NDBASE,C_ZERO,AUXILI(:,:,0),NDBASE)
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,0),NDBASE, &
                        RPM(:,:,1),NDBASE,C_ZERO,AUXILI(:,:,1),NDBASE)
!=======================================================================
!TR[(X.RHO)^2] AND (X.RHO.X)_{IJ}
!=======================================================================
!     INCX=1; INCY=1
!     RES1=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,0),INCX,TRANSPOSE(AUXILI(:,:,0)),INCY)
!     RES2=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,1),INCX,TRANSPOSE(AUXILI(:,:,1)),INCY)
!     RESULT=RES1+RES2
     RESULT=C_ZERO
     DO IBRA=1,LDBASE
        DO IKET=1,LDBASE
           RESULT=RESULT+AUXILI(IBRA,IKET,0)*AUXILI(IKET,IBRA,0)&
                        +AUXILI(IBRA,IKET,1)*AUXILI(IKET,IBRA,1)
        END DO
     END DO
     IF (JETKET.EQ.2) THEN
         CTWO=CMPLX(-2.0D0,0.0D0)
         CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                            CTWO,AUXILI(:,:,1),NDBASE, &
                            XMAT(:,:,0),NDBASE,C_ZERO,APM(:,:,0),NDBASE)
         CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                            CTWO,AUXILI(:,:,0),NDBASE, &
                            XMAT(:,:,1),NDBASE,C_ZERO,APM(:,:,1),NDBASE)
     END IF
     !
     EOBSSQ = EOBSSQ + REAL(RESULT)
     EKBSSQ = EKBSSQ + RESULT
     !
 END IF
!=======================================================================
 RETURN
END SUBROUTINE EXCXSO
!=======================================================================
SUBROUTINE PAIXSD(MREVER,ISIMPY,JETKET,POBSSQ,PKBSSQ,XMAT,KPM1,KPM2,KPP1,KPP2,PPM,PPP)
!=======================================================================
! Tr[(X.kappa)(X.kappa')*] x=[x(0),0;0,x(1)]
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,JETKET
 REAL   ::POBSSQ
 COMPLEX::PKBSSQ
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::XMAT,KPM1,KPM2
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE),OPTIONAL::KPP1,KPP2,PPM,PPP
!
 INTEGER::IREVER
 COMPLEX::C_ZERO,RESULT,C_UNIT,CTWO
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::AUXILI
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::AUXIL2
!
! INTEGER::INCX,INCY
! COMPLEX::RES1,RES2,ZDOTU
 INTEGER::IBRA,IKET
!
 INTEGER::LDBASE
 COMMON /DIMENS/ LDBASE
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 IF(JETKET.EQ.2)THEN
                    PPM(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
    IF(ISIMPY.EQ.0) PPP(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
 END IF
!=======================================================================
!X.KAPPA AND X.KAPPA'
!=======================================================================
 C_UNIT=CMPLX(1.0D0,0.0D0); AUXILI = C_ZERO
 CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,0),NDBASE, &
                    KPM1(:,:,0),NDBASE,C_ZERO,AUXILI(:,:,0),NDBASE)
 C_UNIT=CMPLX(-1.0D0,0.0D0)
 CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,1),NDBASE, &
                    KPM1(:,:,0),NDBASE,C_ZERO,AUXILI(:,:,1),NDBASE)
 C_UNIT=CMPLX(1.0D0,0.0D0); AUXIL2 = C_ZERO
 CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,0),NDBASE, &
                    KPM2(:,:,0),NDBASE,C_ZERO,AUXIL2(:,:,0),NDBASE)
 C_UNIT=CMPLX(-1.0D0,0.0D0)
 CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,1),NDBASE, &
                    KPM2(:,:,0),NDBASE,C_ZERO,AUXIL2(:,:,1),NDBASE)
!=======================================================================
!TR[(X.KAPPA)(X.KAPPA')*] AND 2[X.KAPPA.X']
!=======================================================================
! INCX=1; INCY=1
! RES1=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,0),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,1))),INCY)
! RES2=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,1),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,0))),INCY)
! RESULT=RES1+RES2
 RESULT=C_ZERO
 DO IBRA=1,LDBASE
    DO IKET=1,LDBASE
       RESULT=RESULT+AUXILI(IBRA,IKET,0)*CONJG(AUXIL2(IKET,IBRA,1))&
                    +AUXILI(IBRA,IKET,1)*CONJG(AUXIL2(IKET,IBRA,0))
    END DO
 END DO
 IF (JETKET.EQ.2) THEN
     CTWO=CMPLX(2.0D0,0.0D0)
     CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                        CTWO,AUXILI(:,:,1),NDBASE, &
                        XMAT(:,:,0),NDBASE,C_ZERO,PPM(:,:,1),NDBASE)
     CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                        CTWO,AUXILI(:,:,0),NDBASE, &
                        XMAT(:,:,1),NDBASE,C_ZERO,PPM(:,:,0),NDBASE)
 END IF
!=======================================================================
!    attention: before version 11 and 278h, the two IREVER components
!               above where exchanged, which did not conform  to  the
!               standard of defining arrays PLINPM and PROTPM.
!               This bug was corrected on 06/12/2016
!=======================================================================
 !
 POBSSQ = -REAL(RESULT)
 PKBSSQ = -RESULT
 !
!=======================================================================
!ADDITIONAL TERMS WHEN SIMPLEX IS NOT CONSERVED
!=======================================================================
 IF (ISIMPY.NE.1) THEN
!=======================================================================
!X.KAPPA
!=======================================================================
     C_UNIT=CMPLX(1.0D0,0.0D0)
     AUXILI = C_ZERO
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,0),NDBASE, &
                        KPP1(:,:,0),NDBASE,C_ZERO,AUXILI(:,:,0),NDBASE)
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,1),NDBASE, &
                        KPP1(:,:,1),NDBASE,C_ZERO,AUXILI(:,:,1),NDBASE)
     AUXIL2 = C_ZERO
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,0),NDBASE, &
                        KPP2(:,:,0),NDBASE,C_ZERO,AUXIL2(:,:,0),NDBASE)
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,1),NDBASE, &
                        KPP2(:,:,1),NDBASE,C_ZERO,AUXIL2(:,:,1),NDBASE)
!=======================================================================
!TR[(X.KAPPA)(X.KAPPA')*] AND 2[X.KAPPA.X]
!=======================================================================
!     INCX=1; INCY=1
!     RES1=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,0),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,0))),INCY)
!     RES2=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,1),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,1))),INCY)
!     RESULT=RES1+RES2
     RESULT=C_ZERO
     DO IBRA=1,LDBASE
        DO IKET=1,LDBASE
           RESULT=RESULT+AUXILI(IBRA,IKET,0)*CONJG(AUXIL2(IKET,IBRA,0))&
                        +AUXILI(IBRA,IKET,1)*CONJG(AUXIL2(IKET,IBRA,1))
        END DO
     END DO
     IF (JETKET.EQ.2) THEN
         CTWO=CMPLX(2.0D0,0.0D0)
         CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                            CTWO,AUXILI(:,:,0),NDBASE, &
                            XMAT(:,:,0),NDBASE,C_ZERO,PPP(:,:,0),NDBASE)
         CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                            CTWO,AUXILI(:,:,1),NDBASE, &
                            XMAT(:,:,1),NDBASE,C_ZERO,PPP(:,:,1),NDBASE)
     END IF
     !
     POBSSQ = POBSSQ - REAL(RESULT)
     PKBSSQ = PKBSSQ - RESULT
     !
 END IF
!=======================================================================
 RETURN
END SUBROUTINE PAIXSD
!=======================================================================
SUBROUTINE PAIXSO(MREVER,ISIMPY,JETKET,POBSSQ,PKBSSQ,XMAT,KPM1,KPM2,KPP1,KPP2,PPM,PPP)
!=======================================================================
! Tr[(X.kappa)(X.kappa')*] x=[0,x(0);x(1),0]
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,JETKET
 REAL   ::POBSSQ
 COMPLEX::PKBSSQ
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::XMAT,KPM1,KPM2
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE),OPTIONAL::KPP1,KPP2,PPM,PPP
!
 INTEGER::IREVER
 COMPLEX::C_ZERO,RESULT,C_UNIT,CTWO
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::AUXILI
 COMPLEX,DIMENSION(1:NDBASE,1:NDBASE,0:NDREVE)::AUXIL2
!
! INTEGER::INCX,INCY
! COMPLEX::RES1,RES2,ZDOTU
 INTEGER::IBRA,IKET
!
 INTEGER::LDBASE
 COMMON /DIMENS/ LDBASE
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 IF(JETKET.EQ.2)THEN
                    PPM(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
    IF(ISIMPY.EQ.0) PPP(1:NDBASE,1:NDBASE,0:NDREVE)=C_ZERO
 END IF
!=======================================================================
!X.KAPPA AND X.KAPPA'
!=======================================================================
 C_UNIT=CMPLX(1.0D0,0.0D0); AUXILI = C_ZERO
 CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,0),NDBASE, &
                    KPM1(:,:,1),NDBASE,C_ZERO,AUXILI(:,:,0),NDBASE)
 C_UNIT=CMPLX(-1.0D0,0.0D0)
 CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,1),NDBASE, &
                    KPM1(:,:,1),NDBASE,C_ZERO,AUXILI(:,:,1),NDBASE)
 C_UNIT=CMPLX(1.0D0,0.0D0); AUXIL2 = C_ZERO
 CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,0),NDBASE, &
                    KPM2(:,:,1),NDBASE,C_ZERO,AUXIL2(:,:,0),NDBASE)
 C_UNIT=CMPLX(-1.0D0,0.0D0)
 CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                    C_UNIT,XMAT(:,:,1),NDBASE, &
                    KPM2(:,:,1),NDBASE,C_ZERO,AUXIL2(:,:,1),NDBASE)
!=======================================================================
!TR[(X.KAPPA)(X.KAPPA')*] AND 2[X.KAPPA.X']
!=======================================================================
! INCX=1; INCY=1
! RES1=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,0),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,0))),INCY)
! RES2=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,1),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,1))),INCY)
! RESULT=RES1+RES2
 RESULT=C_ZERO
 DO IBRA=1,LDBASE
    DO IKET=1,LDBASE
       RESULT=RESULT+AUXILI(IBRA,IKET,0)*CONJG(AUXIL2(IKET,IBRA,0))&
                    +AUXILI(IBRA,IKET,1)*CONJG(AUXIL2(IKET,IBRA,1))
    END DO
 END DO
 IF (JETKET.EQ.2) THEN
     CTWO=CMPLX(2.0D0,0.0D0)
     CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                        CTWO,AUXILI(:,:,1),NDBASE, &
                        XMAT(:,:,0),NDBASE,C_ZERO,PPM(:,:,1),NDBASE)
     CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                        CTWO,AUXILI(:,:,0),NDBASE, &
                        XMAT(:,:,1),NDBASE,C_ZERO,PPM(:,:,0),NDBASE)
 END IF
!=======================================================================
!    attention: before version 11 and 278h, the two IREVER components
!               above where exchanged, which did not conform  to  the
!               standard of defining arrays PLINPM and PROTPM.
!               This bug was corrected on 06/12/2016
!=======================================================================
 !
 POBSSQ = -REAL(RESULT)
 PKBSSQ = -RESULT
 !
!=======================================================================
!ADDITIONAL TERMS WHEN SIMPLEX IS NOT CONSERVED
!=======================================================================
 IF (ISIMPY.NE.1) THEN
!=======================================================================
!X.KAPPA AND X.KAPPA'
!=======================================================================
     C_UNIT=CMPLX(1.0D0,0.0D0)
     AUXILI = C_ZERO
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,1),NDBASE, &
                        KPP1(:,:,0),NDBASE,C_ZERO,AUXILI(:,:,0),NDBASE)
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,0),NDBASE, &
                        KPP1(:,:,1),NDBASE,C_ZERO,AUXILI(:,:,1),NDBASE)
     AUXIL2 = C_ZERO
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,1),NDBASE, &
                        KPP2(:,:,0),NDBASE,C_ZERO,AUXIL2(:,:,0),NDBASE)
     CALL ZGEMM('N','N',LDBASE,LDBASE,LDBASE, &
                        C_UNIT,XMAT(:,:,0),NDBASE, &
                        KPP2(:,:,1),NDBASE,C_ZERO,AUXIL2(:,:,1),NDBASE)
!=======================================================================
!TR[(X.KAPPA)(X.KAPPA')*] AND 2[X.KAPPA.X]
!=======================================================================
!     INCX=1; INCY=1
!     RES1=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,1),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,0))),INCY)
!     RES2=ZDOTU(NDBASE*NDBASE,AUXILI(:,:,0),INCX,CONJG(TRANSPOSE(AUXIL2(:,:,1))),INCY)
!     RESULT=RES1+RES2
     RESULT=C_ZERO
     DO IBRA=1,LDBASE
        DO IKET=1,LDBASE
           RESULT=RESULT+AUXILI(IBRA,IKET,1)*CONJG(AUXIL2(IKET,IBRA,0))&
                        +AUXILI(IBRA,IKET,0)*CONJG(AUXIL2(IKET,IBRA,1))
        END DO
     END DO
     IF (JETKET.EQ.2) THEN
         CTWO=CMPLX(2.0D0,0.0D0)
         CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                            CTWO,AUXILI(:,:,1),NDBASE, &
                            XMAT(:,:,0),NDBASE,C_ZERO,PPP(:,:,0),NDBASE)
         CALL ZGEMM('N','T',LDBASE,LDBASE,LDBASE, &
                            CTWO,AUXILI(:,:,0),NDBASE, &
                            XMAT(:,:,1),NDBASE,C_ZERO,PPP(:,:,1),NDBASE)
     END IF
     !
     POBSSQ = POBSSQ - REAL(RESULT)
     PKBSSQ = PKBSSQ - RESULT
     !
 END IF
!=======================================================================
 RETURN
END SUBROUTINE PAIXSO
!=======================================================================
END MODULE O2AVRG
!
!=======================================================================
!=======================================================================
!
!
MODULE OFFDIA
!
!=======================================================================
!      MODULE CONTAINING MANY ROUTINES USED IN THE NEXT ROUTINES BELOW
!      AND THAT COMPUTES MASS AND INERTIA RENORMALIZATION
!=======================================================================
!
!routines to calculate transformation matrix, transformed wave function or energy kernal
USE hfodd_sizes
COMPLEX,ALLOCATABLE::VRVAUX(:,:),URUAUX(:,:)
!
CONTAINS
SUBROUTINE WLRMAT(IHFB,LDBASE,LUPPER,WALEFT,WARIGH,OVRV,OVRU)
!CALCULATE OVRLAP MATRIX FROM WL, WR
IMPLICIT NONE
 INTEGER::IHFB,LDBASE,LUPPER
 COMPLEX,DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN)::WALEFT,WARIGH
 COMPLEX,DIMENSION(1:LUPPER,1:LUPPER)::OVRV
 COMPLEX,DIMENSION(1:LUPPER,1:LUPPER),OPTIONAL::OVRU
!
 COMPLEX::C_ZERO,C_UNIT
!
 C_ZERO=CMPLX(0.0D0,0.0D0)
 C_UNIT=CMPLX(1.0D0,0.0D0)
!
 CALL ZGEMM('C','N',LUPPER,LUPPER,LDBASE,C_UNIT,&
             WALEFT(1:LDBASE,1:LUPPER,0),LDBASE,&
             WARIGH(1:LDBASE,1:LUPPER,0),LDBASE,C_ZERO,OVRV,LUPPER)
!
 CALL ZGEMM('C','N',LUPPER,LUPPER,LDBASE,C_UNIT,&
             WALEFT(1:LDBASE,1:LUPPER,1),LDBASE,&
             WARIGH(1:LDBASE,1:LUPPER,1),LDBASE,C_UNIT,OVRV,LUPPER)
!
 IF(IHFB.GE.1)THEN
    CALL ZGEMM('C','N',LUPPER,LUPPER,LDBASE,C_UNIT,&
                WALEFT(1:LDBASE,1+LUPPER:2*LUPPER,0),LDBASE,&
                WARIGH(1:LDBASE,1+LUPPER:2*LUPPER,0),LDBASE,C_ZERO,OVRU,LUPPER)
!
    CALL ZGEMM('C','N',LUPPER,LUPPER,LDBASE,C_UNIT,&
                WALEFT(1:LDBASE,1+LUPPER:2*LUPPER,1),LDBASE,&
                WARIGH(1:LDBASE,1+LUPPER:2*LUPPER,1),LDBASE,C_UNIT,OVRU,LUPPER)
 ENDIF
!
 RETURN
END SUBROUTINE WLRMAT
!=======================================================================
SUBROUTINE MAT_WR(IHFB,LDBASE,LUPPER,OVRLAP,WARIGH)
IMPLICIT NONE
 INTEGER::IHFB,LDBASE,LUPPER
 COMPLEX,DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN)::WARIGH
 COMPLEX,DIMENSION(1:LUPPER,1:LUPPER)::OVRLAP
!
 INTEGER::ISPIN
 COMPLEX::C_UNIT,C_ZERO
 COMPLEX,ALLOCATABLE::WARNEW(:,:)
!
 C_ZERO=CMPLX(0.0D0,0.0D0)
 C_UNIT=CMPLX(1.0D0,0.0D0)
 ALLOCATE(WARNEW(1:LDBASE,1:LUPPER))

 DO ISPIN=0,NDSPIN
    CALL ZGEMM('N','N',LDBASE,LUPPER,LUPPER,C_UNIT,&
             WARIGH(1:LDBASE,1:LUPPER,ISPIN),LDBASE,&
             OVRLAP(1:LUPPER,1:LUPPER),      LUPPER,&
             C_ZERO,WARNEW,LDBASE)
    WARIGH(1:LDBASE,1:LUPPER,ISPIN)=WARNEW(1:LDBASE,1:LUPPER)
!
    IF(IHFB.GE.1)THEN
       CALL ZGEMM('N','N',LDBASE,LUPPER,LUPPER,C_UNIT,&
                WARIGH(1:LDBASE,1+LUPPER:2*LUPPER,ISPIN),LDBASE,&
                OVRLAP(1:LUPPER,1       :  LUPPER),      LUPPER,&
                C_ZERO,WARNEW,LDBASE)
       WARIGH(1:LDBASE,1+LUPPER:2*LUPPER,ISPIN)=WARNEW(1:LDBASE,1:LUPPER)
    ENDIF
 ENDDO
 DEALLOCATE(WARNEW)
 RETURN
END SUBROUTINE MAT_WR
!=======================================================================
SUBROUTINE RMATCN(MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,LTOTAL,INROTA,ICHARG,&
                  RMATRX,RMAINV,WALEFT,WARIGH)
!=======================================================================
! RETURN ROTATION MATRIX FOR UPPER AND LOWER BLOCK IN CANONICAL
! BASIS, THE LATTER WILL BE INVERSED LATER
!=======================================================================
USE CANBAS
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,LTOTAL,INROTA,ICHARG
 COMPLEX,DIMENSION(1:LUPPER,1:LUPPER)::RMATRX, RMAINV
 COMPLEX,DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN)::WALEFT,WARIGH
!
 INTEGER:: JSTATE,ISTATE,IREVER,ISPIN,IBASE
 COMPLEX::C_ZERO!,AVRCMB,AVRCMA
!
 INTEGER::NUMBQP
 COMPLEX::PHASPI
 COMMON /SPIPHA/ PHASPI(1:NDBASE,0:NDREVE,0:NDSPIN,1:NDTWCE)
 COMMON /QPNUMS/ NUMBQP(0:NDREVE,0:NDISOS)
!=======================================================================
 CALL CPUTIM('RMATCN',1)
!=======================================================================
 IF(ISIMPY.EQ.0) STOP 'NEED ISIMPY=1 IN RMATCN'
 C_ZERO=CMPLX(0.0D0,0.0D0)
 WALEFT(:,:,:)=C_ZERO
!=======================================================================
 JSTATE=0
 DO IREVER=0,MREVER
    DO ISTATE=1,LTIMUP
       JSTATE=JSTATE+1
       DO ISPIN=0,NDSPIN
!D FOR UPPER BLOCK
          WALEFT(1:LDBASE,JSTATE,ISPIN)=&
          WAVCAN(1:LDBASE,ISTATE,IREVER,ICHARG)*PHASPI(1:LDBASE,IREVER,ISPIN,1)
!D FOR LOWER BLOCK
          WALEFT(1:LDBASE,JSTATE+LUPPER,ISPIN)=(2*ISPIN-1)*&
          CONJG(WAVCAN(1:LDBASE,ISTATE,1-IREVER,ICHARG)*PHASPI(1:LDBASE,1-IREVER,1-ISPIN,1))
       END DO
    END DO
 END DO
! RD
 WARIGH(:,:,:)=WALEFT(:,:,:)
 IF(INROTA.NE.0) CALL ROTWAV(LTOTAL,WARIGH)
! DRD
 RMATRX(:,:)=C_ZERO
 RMAINV(:,:)=C_ZERO
!
 CALL WLRMAT(1,LDBASE,LUPPER,WALEFT,WARIGH,RMATRX,RMAINV)

!(rotated) wavefunction
 JSTATE=0
 DO IREVER=0,MREVER
    DO ISTATE=1,LTIMUP
       JSTATE=JSTATE+1
       DO ISPIN=0,NDSPIN
!
          WALEFT(1:LDBASE,JSTATE,ISPIN)&
         =V_CAN(ISTATE,IREVER,ICHARG)*WALEFT(1:LDBASE,JSTATE,ISPIN)
          WALEFT(1:LDBASE,JSTATE+LUPPER,ISPIN)&
         =U_CAN(ISTATE,1-IREVER,ICHARG)*WALEFT(1:LDBASE,JSTATE+LUPPER,ISPIN)
!
          WARIGH(1:LDBASE,JSTATE,ISPIN)&
         =V_CAN(ISTATE,IREVER,ICHARG)*WARIGH(1:LDBASE,JSTATE,ISPIN)
          WARIGH(1:LDBASE,JSTATE+LUPPER,ISPIN)&
         =U_CAN(ISTATE,1-IREVER,ICHARG)*WARIGH(1:LDBASE,JSTATE+LUPPER,ISPIN)
          END DO
    END DO
 END DO

!=======================================================================
 CALL CPUTIM('RMATCN',0)
!=======================================================================
 RETURN
END SUBROUTINE RMATCN
!=======================================================================
SUBROUTINE VRVMAT(MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,ICHARG, RMAT_V,RMAT_U,OVRLAP)
!=======================================================================
! calculate V*R_V*V+U*CONJG(TRANSPOSE(R_U))*U matrix in canonical basis.
!=======================================================================
USE CANBAS
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,LDBASE,LUPPER,LTIMUP,INROTA,ICHARG
 COMPLEX,DIMENSION(1:LUPPER,1:LUPPER)::OVRLAP,RMAT_V,RMAT_U
!
 INTEGER::JSTATE,ISTATE,JREVER,IREVER,IBRA,JKET,ISPIN,IBASE
 COMPLEX::C_ZERO,AVRCMP
!
 INTEGER::NUMBQP
 COMPLEX::PHASPI
 COMMON /SPIPHA/ PHASPI(1:NDBASE,0:NDREVE,0:NDSPIN,1:NDTWCE)
 COMMON /QPNUMS/ NUMBQP(0:NDREVE,0:NDISOS)
!=======================================================================
 CALL CPUTIM('VRVMAT',1)
!=======================================================================
 IF(ISIMPY.EQ.0) STOP 'NEED ISIMPY=1 IN VRVMAT'
 C_ZERO=CMPLX(0.0D0,0.0D0)
!=======================================================================
 OVRLAP(:,:)=C_ZERO
! CALCULATE URU+VRV
 DO IREVER=0,MREVER!SIMPLEX BLOCK OF LEFT V,U
    DO JREVER=0,MREVER!SIMPLEX BLOCK OF RIGHT V,U
       DO ISTATE=1,LTIMUP
          DO JSTATE=1,LTIMUP
             IBRA=ISTATE+IREVER*LTIMUP!THE LABEL OF R_U,R_V, IT CAN HAVE OFF DIAGONAL PART SINCE ROTATION CAN BREAK SIMPLEX
             JKET=JSTATE+JREVER*LTIMUP
             AVRCMP=V_CAN(ISTATE,IREVER,ICHARG)*RMAT_V(IBRA,JKET)*V_CAN(JSTATE,JREVER,ICHARG)&
                             +U_CAN(ISTATE,1-IREVER,ICHARG)*CONJG(RMAT_U(JKET,IBRA))*U_CAN(JSTATE,1-JREVER,ICHARG)
             OVRLAP(IBRA,JKET)=AVRCMP
          ENDDO
       ENDDO
    ENDDO
 ENDDO
!=======================================================================
      CALL CPUTIM('VRVMAT',0)
!=======================================================================
 RETURN
END SUBROUTINE VRVMAT
!=======================================================================
SUBROUTINE ROTWLR(MREVER,ISIMPY,ICANON,IDOGOA,IPAHFB,INROTA,IEXTEN,LDBASE,LTIMUP,LUPPER,LTOTAL,ICHARG,OVKERN)
!=======================================================================
!CALCULATE B.OVERLAP^-1,A.OVERLAP^-1 IN ROTATIONAL CASE
!ICANON>1 ONLY WHEN ISIMPY=1
!=======================================================================
USE SAVLEF
USE SAVRIG
USE WAVR_L
!USE PFAFFI, ONLY:NORM_N,VUAUX,LV2AUX,PV2AUX
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,ICANON,IDOGOA,IPAHFB,INROTA,IEXTEN,LDBASE,LTIMUP,LUPPER,LTOTAL,ICHARG
 COMPLEX::OVKERN
!
 INTEGER::IALLOC,ISTATE,JSTATE,KSTATE,ISPIN,IBASE,SIGNS,IFLAG,SIGV2,IFGVI2,LBOLD,LBNEW
 INTEGER::LOCCUD(0:NDREVE)
 INTEGER::IAUXDI(1:2*NDSTAT)
 REAL   ::NOTOCC,RICOND,LOGNORM,PHASE,LOGPF,LOGV2,phav2
 COMPLEX::C_ZERO,C_UNIT,C_I,AVRCMP,CPHASE
 COMPLEX::DETWRK(2),DETRMA(2),AUXDIA(1:2*NDSTAT)
 COMPLEX,DIMENSION(1:2*NDSTAT,1:2*NDSTAT)::VU
 COMPLEX,ALLOCATABLE::OVRLAP(:,:),RMATRX(:,:),RMAINV(:,:),EWAVEF(:,:),WARAUX(:,:)
!
 INTEGER::NFIPRI
 REAL   ::ROTATX,ROTATY,ROTATZ
 COMPLEX::PHASPI
 COMMON /CFIPRI/ NFIPRI
 COMMON /ROTATT/ ROTATX,ROTATY,ROTATZ
 COMMON /SPIPHA/ PHASPI(1:NDBASE,0:NDREVE,0:NDSPIN,1:NDTWCE)
!=======================================================================
 IALLOC=1
!=======================================================================
 ALLOCATE (OVRLAP(1:2*NDSTAT,1:2*NDSTAT),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('OVRLAP','ROTWLR')
!
 IF(ICANON.EQ.2)THEN
    ALLOCATE (RMATRX(1:2*NDSTAT,1:2*NDSTAT),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('RMATRX','ROTWLR')
    ALLOCATE (RMAINV(1:2*NDSTAT,1:2*NDSTAT),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('RMAINV','ROTWLR')
 ENDIF
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 C_UNIT=CMPLX(1.0D0,0.0D0)
 C_I=CMPLX(0.0D0,1.0D0)
 LOCCUD(:)=LTIMUP
!=======================================================================
!(BRB+ARA) IN HFB CASE
!=======================================================================
 IF(IPAHFB.EQ.1)THEN

    IF(ICANON.LT.2)THEN
!=======================================================================
!PUT WAVE FUNCTION INTO WALEFT AND WARIGH
!=======================================================================
       IF(ICANON.EQ.1)THEN
          CALL CANWAV(MREVER,ISIMPY,WALEFT,ICHARG,LDBASE,LOCCUD)
          WARIGH(:,:,:)=WALEFT(:,:,:)
       ELSE
          WALEFT(:,:,:)=SARIGH(:,:,:,ICHARG,1)
          WARIGH(:,:,:)=SARIGH(:,:,:,ICHARG,1)
       ENDIF
!=======================================================================
!ROTATE WARIGH
!=======================================================================
       IF(INROTA.NE.0) CALL ROTWAV(LTOTAL,WARIGH)
!=======================================================================
!/CALCULATE OVERLAP MATRIX
!=======================================================================
       OVRLAP(:,:)=C_ZERO
       IAUXDI(:)=0
       AUXDIA(:)=0.0D0
!
       ALLOCATE(VRVAUX(1:LUPPER,1:LUPPER))
       ALLOCATE(URUAUX(1:LUPPER,1:LUPPER))
!
       CALL WLRMAT(IPAHFB,LDBASE,LUPPER,WALEFT,WARIGH,VRVAUX,URUAUX)
       OVRLAP(1:LUPPER,1:LUPPER)=VRVAUX(1:LUPPER,1:LUPPER)+URUAUX(1:LUPPER,1:LUPPER)
!
       DEALLOCATE(VRVAUX,URUAUX)
!=======================================================================
! INVERTING THE OVERLAP MATRIX
!=======================================================================
       CALL ZGECO(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,RICOND,AUXDIA)
       CALL ZGEDI(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,DETWRK,AUXDIA,01)
       IF(IDOGOA.NE.0)THEN!GOA NEED THE KERNAL
          IF(MREVER.EQ.0) STOP 'CALCULATE KERNAL WITH PFAFFIAN IN HFB CASE NEED MREVER=1 IN ROTWLR'
          !USE PFAFFIAN METHOD TO CALCULATE THE KERNAL
           IF(INROTA.NE.0)THEN
              STOP 'ROTWLR NEEDS TO BE ADAPTED SINCE LIPCOR VERSION 16'
!              ALLOCATE(VUAUX(1:2*NDSTAT,1:2*NDSTAT,0:NDISOS))
!              CALL NORM_N(0,IPAHFB,LUPPER,LDBASE,WALEFT,WALEFT,VUAUX(:,:,ICHARG),LV2AUX(ICHARG),PV2AUX(ICHARG),OVKERN)
!              CALL NORM_N(1,IPAHFB,LUPPER,LDBASE,WALEFT,WARIGH,VUAUX(:,:,ICHARG),LV2AUX(ICHARG),PV2AUX(ICHARG),OVKERN)
!in hfb case, pfaffian is the overlap; in hf case, determinant of d^+ d is the overlap
!              DEALLOCATE(VUAUX)
           ELSE!the kernal is 1 if there is no rotation
              OVKERN=1.0d0
           ENDIF
       ENDIF
!
       IF (ABS(RICOND).LT.1.0D-12) THEN
          WRITE(NFIPRI,*)'BE CAREFUL - R MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
          WRITE(NFIPRI,*)'INROTA=',INROTA,' RICOND=',RICOND
       END IF
!
    ELSEIF(ICANON.EQ.2)THEN
!=======================================================================
! CONSTRUT R AND R^-1 MATRIX UNDER CANONICAL BASIS
! ALSO CALCULATE WAVE FUNCTION AND ROTATED WAVE FUNCTION
!=======================================================================
       CALL RMATCN(MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,LTOTAL,INROTA,ICHARG,&
                   RMATRX(1:LUPPER,1:LUPPER),RMAINV(1:LUPPER,1:LUPPER),&
                   WALEFT,WARIGH)
!
       IAUXDI(:)=0
       AUXDIA(:)=0.0D0
       CALL ZGECO(RMAINV,2*NDSTAT,LUPPER,IAUXDI,RICOND,AUXDIA)
       CALL ZGEDI(RMAINV,2*NDSTAT,LUPPER,IAUXDI,DETRMA,AUXDIA,11)
!
       IF (ABS(RICOND).LT.1.0D-12) THEN
          WRITE(NFIPRI,*)'BE CAREFUL - R MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
          WRITE(NFIPRI,*)'INROTA=',INROTA,' RICOND=',RICOND
       END IF
!=======================================================================
!(BRB+AR^-1A) IN HFB CASE
!=======================================================================
       CALL VRVMAT(MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,ICHARG,&
                   RMATRX(1:LUPPER,1:LUPPER),RMAINV(1:LUPPER,1:LUPPER),&
                   OVRLAP(1:LUPPER,1:LUPPER))
!=======================================================================
! INVERTING THE OVERLAP MATRIX
!=======================================================================
       IAUXDI(:)=0
       AUXDIA(:)=0.0D0
       CALL ZGECO(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,RICOND,AUXDIA)
       CALL ZGEDI(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,DETWRK,AUXDIA,11)
!
       IF (ABS(RICOND).LT.1.0D-12) THEN
          WRITE(NFIPRI,*)'BE CAREFUL - R MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
          WRITE(NFIPRI,*)'INROTA=',INROTA,' RICOND=',RICOND
       END IF
!
       OVKERN=DETWRK(1)*10.0D0**DETWRK(2)*(DETRMA(1)*10.0D0**DETRMA(2))

    ELSE
          STOP 'WRONG ICANO IN RENINE'

    ENDIF
!=======================================================================
!MULTIPLYING THE RIGHT WAVE FUNCTIONS BY THE INVERTED OVERLAP MATRIX
!=======================================================================
    CALL MAT_WR(IPAHFB,LDBASE,LUPPER,OVRLAP(1:LUPPER,1:LUPPER),WARIGH)
 ELSE
!=======================================================================
!CALCULATE B(BRB+ARA)^-1, A(BRB+ARA)^-1 IN HF CASE
!=======================================================================
    WALEFT(:,:,:)=SARIGH(:,:,:,ICHARG,1)
    WARIGH(:,:,:)=SARIGH(:,:,:,ICHARG,1)
    IF(INROTA.NE.0) CALL ROTWAV(LTOTAL,WARIGH)
!
    CALL WLRMAT(IPAHFB,LDBASE,LTOTAL,WALEFT,WARIGH,OVRLAP(1:LTOTAL,1:LTOTAL))
!
!=======================================================================
! INVERTING THE OVERLAP MATRIX
!=======================================================================
    CALL ZGECO(OVRLAP,2*NDSTAT,LTOTAL,IAUXDI,RICOND,AUXDIA)
!
    IF(IDOGOA.GE.1)THEN
       CALL ZGEDI(OVRLAP,2*NDSTAT,LTOTAL,IAUXDI,DETWRK,AUXDIA,11)
       OVKERN=DETWRK(1)*10.0D0**DETWRK(2)
    ELSE
       CALL ZGEDI(OVRLAP,2*NDSTAT,LTOTAL,IAUXDI,DETWRK,AUXDIA,01)
    ENDIF
!
    IF (ABS(RICOND).LT.1.0D-12) THEN
        WRITE(NFIPRI,*) 'BE CAREFULL - OVERLAP MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
        WRITE(NFIPRI,*) 'INROTA=',INROTA,' RICOND=',RICOND
    END IF
!=======================================================================
!MULTIPLYING THE RIGHT WAVE FUNCTIONS BY THE INVERTED OVERLAP MATRIX
!=======================================================================
    CALL MAT_WR(IPAHFB,LDBASE,LTOTAL,OVRLAP(1:LTOTAL,1:LTOTAL),WARIGH)
 END IF ! End IF (IPAHFB.EQ.1)
!=======================================================================
 DEALLOCATE(OVRLAP)
 IF(ICANON.EQ.2)DEALLOCATE(RMATRX,RMAINV)
!=======================================================================
END SUBROUTINE ROTWLR
!=======================================================================
SUBROUTINE LMATCN(MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,LTOTAL,INSHIF,ICHARG,&
                  LMATRX,LMAINV,WALEFT,WARIGH)
!=======================================================================
! RETURN ROTATION MATRIX FOR UPPER AND LOWER BLOCK IN CANONICAL
! BASIS, THE LATTER WILL BE INVERSED LATER
!=======================================================================
USE CANBAS
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,LTOTAL,INSHIF,ICHARG
 COMPLEX,DIMENSION(1:LUPPER,1:LUPPER)::LMATRX, LMAINV
 COMPLEX,DIMENSION(1:NDBASE,1:4*NDSTAT,0:NDSPIN)::WALEFT,WARIGH
!
 INTEGER:: JSTATE,ISTATE,IREVER,ISPIN,IBASE
 COMPLEX::C_ZERO,AVRCMB,AVRCMA
!
 INTEGER::NUMBQP
 COMPLEX::PHASPI
 COMMON /SPIPHA/ PHASPI(1:NDBASE,0:NDREVE,0:NDSPIN,1:NDTWCE)
 COMMON /QPNUMS/ NUMBQP(0:NDREVE,0:NDISOS)
!=======================================================================
 CALL CPUTIM('LMATCN',1)
!=======================================================================
 IF(ISIMPY.EQ.0) STOP 'NEED ISIMPY=1 IN LMATCN'
 C_ZERO=CMPLX(0.0D0,0.0D0)
 WALEFT(:,:,:)=C_ZERO
!=======================================================================
 JSTATE=0
 DO IREVER=0,MREVER
    DO ISTATE=1,LTIMUP
       JSTATE=JSTATE+1
       DO ISPIN=0,NDSPIN
!D FOR UPPER BLOCK
          WALEFT(1:LDBASE,JSTATE,ISPIN)=&
          WAVCAN(1:LDBASE,ISTATE,IREVER,ICHARG)*PHASPI(1:LDBASE,IREVER,ISPIN,1)
!D FOR LOWER BLOCK
          WALEFT(1:LDBASE,JSTATE+LUPPER,ISPIN)=(2*ISPIN-1)*&
          CONJG(WAVCAN(1:LDBASE,ISTATE,1-IREVER,ICHARG)*PHASPI(1:LDBASE,1-IREVER,1-ISPIN,1))
       END DO
    END DO
 END DO
! RD
 WARIGH(:,:,:)=WALEFT(:,:,:)
 IF(INSHIF.NE.0) CALL SHIWAV(LTOTAL,WARIGH)
! DRD
 LMATRX(:,:)=C_ZERO
 LMAINV(:,:)=C_ZERO
!
 CALL WLRMAT(1,LDBASE,LUPPER,WALEFT,WARIGH,LMATRX,LMAINV)
!(shifted) wavefunction
 JSTATE=0
 DO IREVER=0,MREVER
    DO ISTATE=1,LTIMUP
       JSTATE=JSTATE+1
       DO ISPIN=0,NDSPIN
!
          WALEFT(1:LDBASE,JSTATE,ISPIN)&
         =V_CAN(ISTATE,IREVER,ICHARG)*WALEFT(1:LDBASE,JSTATE,ISPIN)
          WALEFT(1:LDBASE,JSTATE+LUPPER,ISPIN)&
         =U_CAN(ISTATE,1-IREVER,ICHARG)*WALEFT(1:LDBASE,JSTATE+LUPPER,ISPIN)
!
          WARIGH(1:LDBASE,JSTATE,ISPIN)&
         =V_CAN(ISTATE,IREVER,ICHARG)*WARIGH(1:LDBASE,JSTATE,ISPIN)
          WARIGH(1:LDBASE,JSTATE+LUPPER,ISPIN)&
         =U_CAN(ISTATE,1-IREVER,ICHARG)*WARIGH(1:LDBASE,JSTATE+LUPPER,ISPIN)
          END DO
    END DO
 END DO

!=======================================================================
 CALL CPUTIM('LMATCN',0)
!=======================================================================
 RETURN
END SUBROUTINE LMATCN
!=======================================================================
SUBROUTINE LINWLR(MREVER,ISIMPY,ICANON,IDOGOA,IPAHFB,INSHIF,IEXTEN,LDBASE,LTIMUP,LUPPER,LTOTAL,ICHARG,OVKERN)
!=======================================================================
!CALCULATE B.OVERLAP^-1,A.OVERLAP^-1 IN ROTATIONAL CASE
!ICANON>1 ONLY WHEN ISIMPY=1
!=======================================================================
USE SAVLEF
USE SAVRIG
USE WAVR_L
!USE PFAFFI, ONLY:NORM_N,VUAUX,LV2AUX,PV2AUX
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,ICANON,IDOGOA,IPAHFB,INSHIF,IEXTEN,LDBASE,LTIMUP,LUPPER,LTOTAL,ICHARG
 COMPLEX::OVKERN
!
 INTEGER::IALLOC,ISTATE,JSTATE,KSTATE,ISPIN,IBASE
 INTEGER::LOCCUD(0:NDREVE)
 INTEGER::IAUXDI(1:2*NDSTAT)
 REAL   ::RICOND
 COMPLEX::C_ZERO,C_UNIT,AVRCMP
 COMPLEX::DETWRK(2),DETRMA(2),AUXDIA(1:2*NDSTAT)
 COMPLEX,ALLOCATABLE::OVRLAP(:,:),LMATRX(:,:),LMAINV(:,:),EWAVEF(:,:),WARAUX(:,:)
!
 INTEGER::NFIPRI
 COMMON /CFIPRI/ NFIPRI
!=======================================================================
 IALLOC=1
!=======================================================================
 ALLOCATE (OVRLAP(1:2*NDSTAT,1:2*NDSTAT),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('OVRLAP','LINWLR')
!
 IF(ICANON.EQ.2)THEN
    ALLOCATE (LMATRX(1:2*NDSTAT,1:2*NDSTAT),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('LMATRX','LINWLR')
    ALLOCATE (LMAINV(1:2*NDSTAT,1:2*NDSTAT),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('LMAINV','LINWLR')
 ENDIF
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 C_UNIT=CMPLX(1.0D0,0.0D0)
 LOCCUD(:)=LTIMUP
!=======================================================================
!(BRB+ARA) IN HFB CASE
!=======================================================================
 IF(IPAHFB.EQ.1)THEN

    IF(ICANON.LT.2)THEN
!=======================================================================
!PUT WAVE FUNCTION INTO WALEFT AND WARIGH
!=======================================================================
       IF(ICANON.EQ.1)THEN
          CALL CANWAV(MREVER,ISIMPY,WALEFT,ICHARG,LDBASE,LOCCUD)
          WARIGH(:,:,:)=WALEFT(:,:,:)
       ELSE
          WALEFT(:,:,:)=SARIGH(:,:,:,ICHARG,1)
          WARIGH(:,:,:)=SARIGH(:,:,:,ICHARG,1)
       ENDIF
!=======================================================================
!ROTATE WARIGH
!=======================================================================
       IF(INSHIF.NE.0) CALL SHIWAV(LTOTAL,WARIGH)
!=======================================================================
!CALCULATE OVERLAP MATRIX
!=======================================================================
       OVRLAP(:,:)=C_ZERO
       IAUXDI(:)=0
       AUXDIA(:)=0.0D0
!
       DO JSTATE=1,LUPPER
          DO ISTATE=1,LUPPER
             AVRCMP=C_ZERO
             DO ISPIN=0,NDSPIN
                DO IBASE=1,LDBASE
                   AVRCMP=AVRCMP&
            +CONJG(WALEFT(IBASE,ISTATE,ISPIN))*WARIGH(IBASE,JSTATE,ISPIN)&
            +CONJG(WALEFT(IBASE,ISTATE+LUPPER,ISPIN))*WARIGH(IBASE,JSTATE+LUPPER,ISPIN)
                END DO
             END DO
             OVRLAP(ISTATE,JSTATE)=AVRCMP
          END DO
       END DO
!=======================================================================
! INVERTING THE OVERLAP MATRIX
!=======================================================================
       CALL ZGECO(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,RICOND,AUXDIA)
       CALL ZGEDI(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,DETWRK,AUXDIA,01)
       IF(IDOGOA.NE.0)THEN!GOA NEED THE KERNAL
          IF(MREVER.EQ.0) STOP 'CALCULATE KERNAL WITH PFAFFIAN IN HFB CASE NEED MREVER=1 IN LINWLR'
          !USE PFAFFIAN METHOD TO CALCULATE THE KERNAL
           IF(INSHIF.NE.0)THEN
              STOP 'ROTWLR NEEDS TO BE ADAPTED SINCE LIPCOR VERSION 16'
!              ALLOCATE(VUAUX(1:2*NDSTAT,1:2*NDSTAT,0:NDISOS))
!              CALL NORM_N(0,IPAHFB,LUPPER,LDBASE,WALEFT,WALEFT,VUAUX(:,:,ICHARG),LV2AUX(ICHARG),PV2AUX(ICHARG),OVKERN)
!              CALL NORM_N(1,IPAHFB,LUPPER,LDBASE,WALEFT,WARIGH,VUAUX(:,:,ICHARG),LV2AUX(ICHARG),PV2AUX(ICHARG),OVKERN)
!in hfb case, pfaffian is the overlap; in hf case, determinant of d^+ d is the overlap
!              DEALLOCATE(VUAUX)
           ELSE
              OVKERN=1.0d0
           ENDIF
       ENDIF
!
       IF (ABS(RICOND).LT.1.0D-12) THEN
          WRITE(NFIPRI,*)'BE CAREFUL - R MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
          WRITE(NFIPRI,*)'INSHIF=',INSHIF,' RICOND=',RICOND
       END IF
!
    ELSEIF(ICANON.EQ.2)THEN
!=======================================================================
! CONSTRUT R AND R^-1 MATRIX UNDER CANONICAL BASIS
! ALSO CALCULATE WAVE FUNCTION AND ROTATED WAVE FUNCTION
!=======================================================================
       CALL LMATCN(MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,LTOTAL,INSHIF,ICHARG,&
                   LMATRX(1:LUPPER,1:LUPPER),LMAINV(1:LUPPER,1:LUPPER),&
                   WALEFT,WARIGH)
!
       IAUXDI(:)=0
       AUXDIA(:)=0.0D0
       CALL ZGECO(LMAINV,2*NDSTAT,LUPPER,IAUXDI,RICOND,AUXDIA)
       CALL ZGEDI(LMAINV,2*NDSTAT,LUPPER,IAUXDI,DETRMA,AUXDIA,11)
!
       IF (ABS(RICOND).LT.1.0D-12) THEN
          WRITE(NFIPRI,*)'BE CAREFUL - R MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
          WRITE(NFIPRI,*)'INSHIF=',INSHIF,' RICOND=',RICOND
       END IF
!=======================================================================
!(BRB+AR^-1A) IN HFB CASE
!=======================================================================
       CALL VRVMAT(MREVER,ISIMPY,LDBASE,LTIMUP,LUPPER,ICHARG,&
                   LMATRX(1:LUPPER,1:LUPPER),LMAINV(1:LUPPER,1:LUPPER),&
                   OVRLAP(1:LUPPER,1:LUPPER))
! write(*,*)ovrlap(1,1),ovrlap(1,1+ltimup),'ovrlap'
!=======================================================================
! INVERTING THE OVERLAP MATRIX
!=======================================================================
       IAUXDI(:)=0
       AUXDIA(:)=0.0D0
       CALL ZGECO(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,RICOND,AUXDIA)
       CALL ZGEDI(OVRLAP,2*NDSTAT,LUPPER,IAUXDI,DETWRK,AUXDIA,11)
!
       IF (ABS(RICOND).LT.1.0D-12) THEN
          WRITE(NFIPRI,*)'BE CAREFUL - R MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
          WRITE(NFIPRI,*)'INSHIF=',INSHIF,' RICOND=',RICOND
       END IF
!
       OVKERN=DETWRK(1)*10.0D0**DETWRK(2)*(DETRMA(1)*10.0D0**DETRMA(2))
    ELSE
          STOP 'WRONG ICANO IN RENINE'

    ENDIF
!=======================================================================
!MULTIPLYING THE RIGHT WAVE FUNCTIONS BY THE INVERTED OVERLAP MATRIX
!=======================================================================
    ALLOCATE(EWAVEF(1:NDBASE,1:2*NDSTAT))
    ALLOCATE(WARAUX(1:NDBASE,1:2*NDSTAT))
!
    DO ISPIN=0,NDSPIN
       EWAVEF(:,:)=C_ZERO
       WARAUX(:,:)=C_ZERO
!UPPER BLOCK
       DO ISTATE=1,LUPPER
          CALL ZCOPY(NDBASE,WARIGH(1,ISTATE,ISPIN),1,EWAVEF(1,ISTATE),1)
       END DO
       CALL ZGEMM('N','N',LDBASE,LUPPER,LUPPER,&
                   C_UNIT,EWAVEF,NDBASE,OVRLAP,2*NDSTAT,&
                   C_ZERO,WARAUX,NDBASE)
       DO ISTATE=1,LUPPER
          CALL ZCOPY(NDBASE,WARAUX(1,ISTATE),1,WARIGH(1,ISTATE,ISPIN),1)
       END DO
!LOWER BLOCK
       DO ISTATE=1,LUPPER
          CALL ZCOPY(NDBASE,WARIGH(1,ISTATE+LUPPER,ISPIN),1,EWAVEF(1,ISTATE),1)
       END DO
       CALL ZGEMM('N','N',LDBASE,LUPPER,LUPPER,&
                   C_UNIT,EWAVEF,NDBASE,OVRLAP,2*NDSTAT,&
                   C_ZERO,WARAUX,NDBASE)
       DO ISTATE=1,LUPPER
          CALL ZCOPY(NDBASE, WARAUX(1,ISTATE),1,WARIGH(1,ISTATE+LUPPER,ISPIN),1)
       END DO
    END DO
!
    DEALLOCATE(EWAVEF)
    DEALLOCATE(WARAUX)
 ELSE
!=======================================================================
!CALCULATE B(BRB+ARA)^-1, A(BRB+ARA)^-1 IN HF CASE
!=======================================================================
    WALEFT(:,:,:)=SARIGH(:,:,:,ICHARG,1)
    WARIGH(:,:,:)=SARIGH(:,:,:,ICHARG,1)
    IF(INSHIF.NE.0) CALL SHIWAV(LTOTAL,WARIGH)
!
    DO JSTATE=1,LTOTAL
       DO ISTATE=1,LTOTAL
          AVRCMP=C_ZERO
          DO ISPIN=0,NDSPIN
             DO IBASE=1,LDBASE
                AVRCMP=AVRCMP+CONJG(WALEFT(IBASE,ISTATE,ISPIN))*WARIGH(IBASE,JSTATE,ISPIN)
             END DO
          END DO
          OVRLAP(ISTATE,JSTATE)=AVRCMP
       END DO
    END DO
!=======================================================================
! INVERTING THE OVERLAP MATRIX
!=======================================================================
    CALL ZGECO(OVRLAP,2*NDSTAT,LTOTAL,IAUXDI,RICOND,AUXDIA)
    IF(IDOGOA.GE.1)THEN
       CALL ZGEDI(OVRLAP,2*NDSTAT,LTOTAL,IAUXDI,DETWRK,AUXDIA,11)
       OVKERN=DETWRK(1)*10.0D0**DETWRK(2)
    ELSE
       CALL ZGEDI(OVRLAP,2*NDSTAT,LTOTAL,IAUXDI,DETWRK,AUXDIA,01)
    ENDIF
    IF (ABS(RICOND).LT.1.0D-12) THEN
        WRITE(NFIPRI,*) 'BE CAREFULL - OVERLAP MATRIX ,MAY BE SINGULAR IN RENINE FOR:'
        WRITE(NFIPRI,*) 'INSHIF=',INSHIF,' RICOND=',RICOND
    END IF
!=======================================================================
!MULTIPLYING THE RIGHT WAVE FUNCTIONS BY THE INVERTED OVERLAP MATRIX
!=======================================================================
    ALLOCATE (WARAUX(1:2*NDSTAT,0:NDSPIN),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('WARAUX','LINWRL')

    DO JSTATE=1,LDBASE
       DO ISTATE=1,LTOTAL
          DO ISPIN=0,NDSPIN
             WARAUX(ISTATE,ISPIN)=C_ZERO
             DO KSTATE=1,LTOTAL
                WARAUX(ISTATE,ISPIN)=WARAUX(ISTATE,ISPIN)&
          +WARIGH(JSTATE,KSTATE,ISPIN)*OVRLAP(KSTATE,ISTATE)
             END DO
          END DO
       END DO
!
       DO ISTATE=1,LTOTAL
          DO ISPIN=0,NDSPIN
              WARIGH(JSTATE,ISTATE,ISPIN)=WARAUX(ISTATE,ISPIN)
          END DO
       END DO
    END DO
    DEALLOCATE(WARAUX)
 END IF ! End IF (IPAHFB.EQ.1)
!=======================================================================
 DEALLOCATE(OVRLAP)
 IF(ICANON.EQ.2)DEALLOCATE(LMATRX,LMAINV)
!=======================================================================
END SUBROUTINE LINWLR
!=======================================================================
SUBROUTINE H_KERN(ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,IPAHFB,MREVER,&
                  ICHARG,MIN_QP,IPNMIX,ITIREP,&
                  NAMEPN,PRINIT,IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,INIKAR,ISAWAV,IKERNE,&
                  ISHIFY,JMULMO,NUMCOU,BOUCOU,IN_FIX,IZ_FIX,&
                  ISKYRM,EKEKIN,EPAIKI,DROTSQ,DKOTSQ,EKECOD,EKECOX,EKESKY)
!=======================================================================
!CALCULATE RHO, KINETIC, PAIRING AND COULOMB ENERGY
!=======================================================================
IMPLICIT NONE
!=======================================================================
 CHARACTER::NAMEPN*8
 LOGICAL::PRINIT
 INTEGER::ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,IPAHFB,MREVER,&
          ICHARG,MIN_QP,IPNMIX,ITPNMX,ITIREP,&
          IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,INIKAR,ISAWAV,IKERNE,&
          ISHIFY,JMULMO,NUMCOU,IN_FIX,IZ_FIX,ISKYRM
 REAL   ::BOUCOU,ENEKIN,DROTSQ
 COMPLEX::EKEKIN,DKOTSQ,EPAIKI,EKECOD,EKECOX,EKESKY
!
 INTEGER::LDPNMX
 REAL   ::COMULT,EKIN_N,EKIN_P,EKIN_T,EPOT_N,EPOT_P,EPOT_T,ESUM_N,ESUM_P,ESUM_T,EPAI_N,EPAI_P,EPAI_T,&
          EREA_N,EREA_P,EREA_T,ELIP_N,ELIP_P,ELIP_T,ECOULD,ECOULE,ECOULT,ECOULS,ECOULV, EMULCO,&
          EMUSLO,EMUREA,ESIFCO,ESISLO,ESIREA,ESPICO,ESPSLO,ESPREA,ENREAR,ECORCM,ECOR_R,EEVEW0,&
          EODDW0,ENE_W0,ENEVEN,ENEODD,ENESKY,ESTABN,ETOTSP,ETOTFU
 COMPLEX::QMUT_N,QMUT_P,QMUT_T,QMUL_N,QMUL_P,QMUL_T
 COMMON /QMUTTI/ QMUT_N(0:NDMULT,-NDMULT:NDMULT),QMUT_P(0:NDMULT,-NDMULT:NDMULT),QMUT_T(0:NDMULT,-NDMULT:NDMULT)
 COMMON /QMULTI/ QMUL_N(0:NDMULT,-NDMULT:NDMULT),QMUL_P(0:NDMULT,-NDMULT:NDMULT),QMUL_T(0:NDMULT,-NDMULT:NDMULT)
 COMMON /COEMUL/ COMULT(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART)
 COMMON /ALLENE/ EKIN_N,EKIN_P,EKIN_T,EPOT_N,EPOT_P,EPOT_T,ESUM_N,ESUM_P,ESUM_T,&
                 EPAI_N,EPAI_P,EPAI_T,EREA_N,EREA_P,EREA_T,ELIP_N,ELIP_P,ELIP_T,&
                 ECOULD,ECOULE,ECOULT,ECOULS,ECOULV,&
                 EMULCO,EMUSLO,EMUREA,ESIFCO,ESISLO,ESIREA,ESPICO,ESPSLO,ESPREA,&
                 ENREAR,ECORCM,ECOR_R,&
                 EEVEW0,EODDW0,ENE_W0,ENEVEN,ENEODD,ENESKY,ESTABN,ETOTSP,ETOTFU
!
 COMPLEX::C_ZERO
!=======================================================================
 ITPNMX=ICHARG
 C_ZERO=CMPLX(0.0D0,0.0D0)
!=======================================================================
!              DENSITY
!=======================================================================
 CALL ZEDENS(ITPNMX)
 CALL DENSHF(ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,IPAHFB,MREVER,&
             ICHARG,MIN_QP,IPNMIX,ITPNMX,ITIREP,&
             NAMEPN,PRINIT,IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,INIKAR,&
             ISAWAV,IKERNE,1,1)
 IF(MREVER.EQ.0) CALL DBLING(ITPNMX,1)
!=======================================================================
!              KINETIC ENERGY
!=======================================================================
! CALL DENMAC(MREVER,ICHARG,ISIMPY,IPAHFB,WALEFT,WARIGH,1)!CALL BEFORE THIS SUBROUTINE
 CALL EKINET(ENEKIN,DROTSQ,EKEKIN,DKOTSQ,1)
!=======================================================================
!              PAIRING ENERGY
!=======================================================================
 IF (IPAHFB.EQ.1) THEN
    CALL EPAIRI(IN_FIX,IZ_FIX,ITPNMX,EPAIKI)
    IF (MREVER.EQ.0) EPAIKI=EPAIKI*4.0d0
 ELSE
    EPAIKI=C_ZERO
 END IF
!
 IF(ICHARG.EQ.1)THEN
! if(.false.)then
!=======================================================================
!              CALCULATING THE COULOMB ENERGY
!=======================================================================
    CALL MOMETS(ISIMPY,ISIGNY,ISIQTY,QMUL_P,QMUT_P,COMULT,JMULMO,ISHIFY,1)
    CALL COUMAT(NUMCOU,BOUCOU,ISIMPY,IKERNE)
    CALL COULOD(ISIMPY,EKECOD)
    CALL TRUCHD
    CALL COULOE(EKECOX)
 ENDIF
!
 IF(ISKYRM.EQ.1)THEN
!=======================================================================
!              CALCULATING THE SKYRME ENERGY
!=======================================================================
       LDPNMX=1
       CALL TRUTOD(1,1)
       CALL ESKYRM(ENESKY,ENEVEN,ENEODD,&
                   ENREAR,ENE_W0,EEVEW0,EODDW0,EKESKY,LDPNMX,0,1)
 ENDIF
 RETURN
END SUBROUTINE H_KERN
!=======================================================================
END MODULE OFFDIA
!
!=======================================================================
!=======================================================================
!
!
!
!=======================================================================
!      ROUTINES CALLED IN MAIN PROGRAM:
!       - CANWAV: BUILDS ROTATED CANONICAL WAVE FUNCTION
!       - RENINE: CALCULATES THE INERTIA RENORMALIZATION FACTORS
!       - LINAVR: CALCULATES LINEAR MOMENTA
!       - RENMAS: CALCULATES THE MASS RENORMALIZATION FACTORS
!=======================================================================
!
!
!
!=======================================================================
!ROUTINES CALLED IN PROGRAM MAIN
SUBROUTINE CANWAV(MREVER,ISIMPY,WALEFT,ICHARG,LDBASE,LOCCUD)
!=======================================================================
! use u_can, v_can ,wavcan to build A=DU, B=DV
!=======================================================================
USE CANBAS
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 INTEGER::MREVER,ISIMPY,ICHARG,LDBASE,LTIMUP,LUPPER,LSTATE,LTOTAL
 INTEGER::LOCCUD(0:NDREVE)
!
 INTEGER::JSTATE,IREVER,ISTATE,ISPIN
 COMPLEX::C_ZERO
 COMPLEX::WALEFT(1:NDBASE,1:4*NDSTAT,0:NDSPIN)
!
 INTEGER::NUMBQP
 COMPLEX::PHASPI
 COMMON /SPIPHA/ PHASPI(1:NDBASE,0:NDREVE,0:NDSPIN,1:NDTWCE)
 COMMON /QPNUMS/ NUMBQP(0:NDREVE,0:NDISOS)
!=======================================================================
 CALL CPUTIM('CANWAV',1)
!=======================================================================
 IF(ISIMPY.EQ.0) STOP 'NOT IMPLEMENT YET'
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 WALEFT(:,:,:)=C_ZERO
!=======================================================================
 JSTATE=0
! B
! do istate=1,loccud(0)
! write(*,*)v_can(istate,0,icharg)**2+u_can(istate,1,icharg)**2,&
!           v_can(istate,0,icharg)**2+u_can(istate,0,icharg)**2,istate
! enddo
! stop
! V_CAN(:,1,:)=-V_CAN(:,0,:)
! U_CAN(:,1,:)=U_CAN(:,0,:)
 !WAVCAN(:,:,1,:)=CONJG(WAVCAN(:,:,0,:))
 DO IREVER=0,MREVER
    DO ISTATE=1,LOCCUD(IREVER)
          JSTATE=JSTATE+1
          DO ISPIN=0,NDSPIN
             WALEFT(1:LDBASE,JSTATE,ISPIN)=V_CAN(ISTATE,IREVER,ICHARG)&
           * WAVCAN(1:LDBASE,ISTATE,IREVER,ICHARG)* PHASPI(1:LDBASE,IREVER,ISPIN,1)
          END DO
    END DO
 END DO
! A
 DO IREVER=1,(1-MREVER),-1
     DO ISTATE=1,LOCCUD(IREVER)
          JSTATE=JSTATE+1
          DO ISPIN=0,NDSPIN
             WALEFT(1:LDBASE,JSTATE,ISPIN)=U_CAN(ISTATE,1-IREVER,ICHARG) &
          * (2*ISPIN-1)*CONJG(WAVCAN(1:LDBASE,ISTATE,IREVER,ICHARG)* PHASPI(1:LDBASE,IREVER,1-ISPIN,1))
          END DO
     END DO
 END DO
!=======================================================================
      CALL CPUTIM('CANWAV',0)
!=======================================================================
 RETURN
END SUBROUTINE CANWAV
!=======================================================================
SUBROUTINE ROTAVR(ISIMPY,ISIGNY,ISIQTY,&
                  IPAHFB,ICHARG,KETA_R,IROTAT,MREVER,CORROT,IKERNE,&
                  DROTSQ,EROTSQ,PROTSQ,TROTSQ,AROTLI,PROTLI,PROTKI,&
                  DKOTSQ,EKOTSQ,PKOTSQ,TKOTSQ,AKOTLI,PKOTLI,PKOTKI)
!=======================================================================
USE MAT_PP
USE MAT_PM
USE MAD_PP
USE MAD_PM
USE KAP_PM
USE KAP2PM
USE KAP_PP
USE KAP2PP
USE MAP_PP
USE MAP_PM
USE AROSTO
USE PROSTO
USE O2AVRG
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 LOGICAL::CORROT
 INTEGER::ISIMPY,ISIGNY,ISIQTY,IPAHFB,ICHARG,KETA_R,IROTAT,MREVER,IKERNE
 REAL   ::DROTSQ(0:NDKART),EROTSQ(0:NDKART),PROTSQ(0:NDKART),TROTSQ(0:NDKART),&
          AROTLI(0:NDKART),PROTLI(0:NDKART),PROTKI(0:NDKART)
 COMPLEX::DKOTSQ(0:NDKART),EKOTSQ(0:NDKART),PKOTSQ(0:NDKART),TKOTSQ(0:NDKART),&
          AKOTLI(0:NDKART),PKOTLI(0:NDKART),PKOTKI(0:NDKART)
!
 INTEGER::IALLOC,IBRA,KARTEZ,IMSIGN
 REAL   ::FACTOR
 COMPLEX:: RESULT,C_ZERO,UNIT_I
 REAL   ,ALLOCATABLE:: HAUXIL(:,:)
 COMPLEX,ALLOCATABLE:: HAUXIJ(:,:,:)
 COMPLEX,ALLOCATABLE:: RHO_PP(:,:,:),RHO_PM(:,:,:),KAI_PM(:,:,:),KAI2PM(:,:,:),KAI_PP(:,:,:),KAI2PP(:,:,:)
!
 INTEGER::IPHAPP,IPHAPM,IPHAMP,IPHAMM,LDBASE,ITESTA,ITESTO,NOITER,NUMITE
 COMMON /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),&
                 IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
 COMMON /DIMENS/ LDBASE
 COMMON /ITERTS/ ITESTA,ITESTO,NOITER,NUMITE
!=======================================================================
 CALL CPUTIM('ROTAVR',1)
 IALLOC=0
 !
 C_ZERO=CMPLX(0.0D0,0.0D0)
 UNIT_I=CMPLX(0.0D0,1.0D0)
 FACTOR=1.0D0
!=======================================================================
!allocate space for 2[J.rho.J] and 2[J.kAPPa.J] matrix
!=======================================================================
 IF(KETA_R.EQ.2)THEN
    IF(.NOT.ALLOCATED(AROTPP))THEN
       ALLOCATE(AROTPP(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('AROTPP','ROTAVR')
    ENDIF
    IF(.NOT.ALLOCATED(PROTPM) .AND. IPAHFB.EQ.1)THEN
       ALLOCATE(PROTPM(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('PROTPM','ROTAVR')
    ENDIF
!
    IF(ISIMPY.NE.1)THEN
       IF(.NOT.ALLOCATED(AROTPM))THEN
          ALLOCATE(AROTPM(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('AROTPM','ROTAVR')
       ENDIF
       IF(.NOT.ALLOCATED(PROTPP) .AND. IPAHFB.EQ.1)THEN
          ALLOCATE(PROTPP(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('PROTPP','ROTAVR')
       ENDIF
    ENDIF
 ENDIF
!=======================================================================
!allocate space for rho and kappa in case they are not allocate yet
!=======================================================================
 IF(.NOT.ALLOCATED(BIG_PP)) THEN
    ALLOCATE(BIG_PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BIG_PP','ROTAVR')
 ENDIF
 IF(.NOT.ALLOCATED(DEN_PP)) THEN
    ALLOCATE(DEN_PP(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('DEN_PP','ROTAVR')
 ENDIF
!
 IF(.NOT.ALLOCATED(CHI_PM) .AND. IPAHFB.EQ.1)THEN
    ALLOCATE(CHI_PM(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS), STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('CHI_PM','ROTAVR')
 ENDIF
!
 IF(ISIMPY.NE.1)THEN
    IF(.NOT.ALLOCATED(BIG_PM))THEN
        ALLOCATE (BIG_PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
        IF (IALLOC.NE.0) CALL NOALLO('BIG_PM','ROTAVR')
    END IF
    IF(.NOT.ALLOCATED(DEN_PM))THEN
        ALLOCATE (DEN_PM(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS),STAT=IALLOC)
        IF (IALLOC.NE.0) CALL NOALLO('DEN_PM','ROTAVR')
    ENDIF
!
    IF (.NOT.ALLOCATED(CHI_PP) .AND. IPAHFB.EQ.1)THEN
       ALLOCATE (CHI_PP(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('CHI_PP','ROTAVR')
    END IF
 ENDIF
!=======================================================================
!allocate work space
!=======================================================================
 ALLOCATE (HAUXIL(1:NDBASE,1:NDBASE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('HAUXIL','ROTAVR')
 ALLOCATE (HAUXIJ(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('HAUXIJ','ROTAVR')
!
 ALLOCATE (RHO_PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('RHO_PP','ROTAVR')
!
 IF(ISIMPY.NE.1)THEN
    ALLOCATE (RHO_PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('RHO_PM','ROTAVR')
 ENDIF
!
 IF (IPAHFB.EQ.1.AND.KETA_R.GT.0) THEN
    ALLOCATE (KAI_PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('KAI_PM','ROTAVR')
    ALLOCATE (KAI2PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('KAI2PM','ROTAVR')
!
    IF(ISIMPY.NE.1)THEN
       ALLOCATE (KAI_PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('KAI_PP','ROTAVR')
       ALLOCATE (KAI2PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('KAI2PP','ROTAVR')
    ENDIF
 ENDIF
 IF (KETA_R.NE.2)THEN
     IF(.NOT.ALLOCATED(BIG_PP)) STOP 'BIG_PP NOT ALLOCATED IN ROTAVR'
     IF(.NOT.ALLOCATED(BIG_PM).AND.ISIMPY.NE.1) STOP 'BIG_PM NOT ALLOCATED IN ROTAVR'
     IF(.NOT.ALLOCATED(PAI_PM).AND.IPAHFB.EQ.1.AND.KETA_R.GT.0) STOP 'PAI_PM NOT ALLOCATED IN ROTAVR'
     IF(.NOT.ALLOCATED(PAI2PM).AND.IPAHFB.EQ.1.AND.KETA_R.GT.0.AND.IKERNE.EQ.1) STOP 'PAI2PM NOT ALLOCATED IN ROTAVR'
     IF(.NOT.ALLOCATED(PAI_PP).AND.IPAHFB.EQ.1.AND.KETA_R.GT.0.AND.ISIMPY.NE.1) STOP 'PAI_PP NOT ALLOCATED IN ROTAVR'
     IF(.NOT.ALLOCATED(PAI2PP).AND.IPAHFB.EQ.1.AND.KETA_R.GT.0.AND.ISIMPY.NE.1.AND.IKERNE.EQ.1) STOP 'PAI2PP NOT ALLOCATED IN ROTAVR'
 ENDIF
!=======================================================================
!    ZEROING ARRAYS OF AVERAGE VALUES
!=======================================================================
 DROTSQ(:)=0.0D0
 EROTSQ(:)=0.0D0
 PROTSQ(:)=0.0D0
 TROTSQ(:)=0.0D0
 AROTLI(:)=0.0D0
 PROTLI(:)=0.0D0
 PROTKI(:)=0.0D0
!
 DKOTSQ(:)=C_ZERO
 EKOTSQ(:)=C_ZERO
 PKOTSQ(:)=C_ZERO
 TKOTSQ(:)=C_ZERO
 AKOTLI(:)=C_ZERO
 PKOTLI(:)=C_ZERO
 PKOTKI(:)=C_ZERO
!
 IF(KETA_R.EQ.2)THEN
                    AROTPP(:,:,:,:,ICHARG)=C_ZERO
    IF(ISIMPY.NE.1) AROTPM(:,:,:,:,ICHARG)=C_ZERO
!
    IF(IPAHFB.EQ.1)THEN
                        PROTPM(:,:,:,:,ICHARG)=C_ZERO
       IF (ISIMPY.NE.1) PROTPP(:,:,:,:,ICHARG)=C_ZERO
    ENDIF
 ENDIF
!
                 RHO_PP(:,:,:)=C_ZERO
 IF(ISIMPY.NE.1) RHO_PM(:,:,:)=C_ZERO
 IF (IPAHFB.EQ.1.AND.KETA_R.GT.0) THEN
    KAI_PM(:,:,:)=C_ZERO
    KAI2PM(:,:,:)=C_ZERO
     IF(ISIMPY.NE.1) THEN
        KAI_PP(:,:,:)=C_ZERO
        KAI2PP(:,:,:)=C_ZERO
     ENDIF
 ENDIF
!=======================================================================
!    store rho and kAPPa into matrix
!=======================================================================
 IF (KETA_R.EQ.2) THEN !when calculating matrix element,use mixed rho and kAPPa
    IF(ITESTA.EQ.NUMITE) THEN!initial den, chi at the first iteration
                        DEN_PP(:,:,:,ICHARG)= BIG_PP(:,:,:)
       IF (ISIMPY.NE.1) DEN_PM(:,:,:,ICHARG)= BIG_PM(:,:,:)
       IF (IPAHFB.EQ.1) THEN
                           CHI_PM(:,:,:,ICHARG)=PAI_PM(:,:,:)
          IF (ISIMPY.NE.1) CHI_PP(:,:,:,ICHARG)=PAI_PP(:,:,:)
       ENDIF
    ENDIF
!
                     RHO_PP(:,:,:)=DEN_PP(:,:,:,ICHARG)
    IF (ISIMPY.NE.1) RHO_PM(:,:,:)=DEN_PM(:,:,:,ICHARG)
    IF (IPAHFB.EQ.1) THEN
       KAI_PM(:,:,:)=CHI_PM(:,:,:,ICHARG)
       KAI2PM(:,:,:)=CHI_PM(:,:,:,ICHARG) ! ANALYZE THAT
       IF (ISIMPY.NE.1)THEN
          KAI_PP(:,:,:)=CHI_PP(:,:,:,ICHARG)
          KAI2PP(:,:,:)=CHI_PP(:,:,:,ICHARG)
       ENDIF
    ENDIF
!
 ELSE!when there is no correction or inside RENINE, use rho and kAPPa from DENMAC and PAIMAC
                     RHO_PP(:,:,:)=BIG_PP(:,:,:)
    IF (ISIMPY.NE.1) RHO_PM(:,:,:)=BIG_PM(:,:,:)
    IF (IPAHFB.EQ.1.AND.KETA_R.GT.0) THEN
        KAI_PM(:,:,:)=PAI_PM(:,:,:)
        IF (IKERNE.EQ.1) THEN
            KAI2PM(:,:,:)=PAI2PM(:,:,:)
        ELSE
            KAI2PM(:,:,:)=PAI_PM(:,:,:)
        ENDIF
        IF (ISIMPY.NE.1) THEN
            KAI_PP(:,:,:)=PAI_PP(:,:,:)
            IF (IKERNE.EQ.1) THEN
                KAI2PP(:,:,:)=PAI2PP(:,:,:)
            ELSE
                KAI2PP(:,:,:)=PAI_PP(:,:,:)
            ENDIF
        ENDIF
    ENDIF
 ENDIF
!=======================================================================
!                         X-DIRECTION
!=======================================================================
!
!=======================================================================
![Jx]
!=======================================================================
 CALL INT_LX(IPHAPM(0,0,0),HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)= UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)=-UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!
 CALL INT_SX(HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIJ(1:NDBASE,1:NDBASE,0)+HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIJ(1:NDBASE,1:NDBASE,1)+HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[Jx.rho]
!=======================================================================
 IF (ISIMPY.NE.1) THEN
! IF (IROTAT.EQ.1 .AND. ISIMPY.NE.1 .AND. ISIQTY.NE.1) THEN
     CALL AVRXSO(HAUXIJ,AROTLI(1),AKOTLI(1),RHO_PM)
 END IF
!
 IF (CORROT) THEN
    IMSIGN=1
!=======================================================================
!Tr[(Jx.rho)^2] and 2[Jx.rho.Jx]
!=======================================================================

    IF(ISIMPY.EQ.1 .AND. KETA_R.NE.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(1),EKOTSQ(1),HAUXIJ,RHO_PP)
    ELSE IF(ISIMPY.EQ.0 .AND. KETA_R.NE.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(1),EKOTSQ(1),HAUXIJ,RHO_PP,RPM=RHO_PM)
    ELSE IF(ISIMPY.EQ.1 .AND.  KETA_R.EQ.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(1),EKOTSQ(1),HAUXIJ,RHO_PP,APP=AROTPP(:,:,1,:,ICHARG))
    ELSE IF(ISIMPY.EQ.0 .AND.  KETA_R.EQ.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(1),EKOTSQ(1),HAUXIJ,RHO_PP,RPM=RHO_PM,&
                   APP=AROTPP(:,:,1,:,ICHARG),APM=AROTPM(:,:,1,:,ICHARG))
    ELSE
       stop 'wrong ISIMPY or KETA_R in rotavr'
    END IF
!=======================================================================
!Tr[(Jx.kAPPa)(Jx.kAPPa')*] and 2[Jx.kAPPa.Jx]
!=======================================================================
    IF (IPAHFB.EQ.1.AND.KETA_R.GT.0) THEN
       IF(ISIMPY.EQ.1 .AND. KETA_R.NE.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(1),PKOTSQ(1),HAUXIJ,KAI_PM,KAI2PM)
       ELSE IF(ISIMPY.EQ.0 .AND. KETA_R.NE.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(1),PKOTSQ(1),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP)
       ELSE IF(ISIMPY.EQ.1 .AND.  KETA_R.EQ.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(1),PKOTSQ(1),HAUXIJ,KAI_PM,KAI2PM,PPM=PROTPM(:,:,1,:,ICHARG))
       ELSE IF(ISIMPY.EQ.0 .AND.  KETA_R.EQ.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(1),PKOTSQ(1),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP,&
                      PPM=PROTPM(:,:,1,:,ICHARG),PPP=PROTPP(:,:,1,:,ICHARG))
       ELSE
         stop 'wrong ISIMPY or KETA_R in rotavr'
       END IF
    END IF

!=======================================================================
![Jx^2]=[Lx^2]+2[Lx.Sx]+[Sx^2]
!=======================================================================
    CALL INTLX2(HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIL(1:NDBASE,1:NDBASE)
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIL(1:NDBASE,1:NDBASE)
!
    CALL INT_LX(IPHAPP(0,0,1),HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIJ(1:NDBASE,1:NDBASE,0)+UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIJ(1:NDBASE,1:NDBASE,1)-UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!
    DO IBRA=1,LDBASE
       HAUXIJ(IBRA,IBRA,0)=HAUXIJ(IBRA,IBRA,0)+CMPLX(0.25D0,0.0D0)
       HAUXIJ(IBRA,IBRA,1)=HAUXIJ(IBRA,IBRA,1)+CMPLX(0.25D0,0.0D0)
    END DO
!=======================================================================
!Tr[Jx^2.rho]
!=======================================================================
    CALL AVRXSD(HAUXIJ,DROTSQ(1),DKOTSQ(1),RHO_PP)
 END IF
!=======================================================================
!                         Z-DIRECTION
!=======================================================================
!
!=======================================================================
![Jz]
!=======================================================================
 CALL INT_LZ(IPHAPM(0,0,0),HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)=UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)=-UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!
 CALL INT_SZ(HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)= HAUXIJ(1:NDBASE,1:NDBASE,0)+UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)= HAUXIJ(1:NDBASE,1:NDBASE,1)-UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[Jz.rho]
!=======================================================================
 IF (ISIMPY.NE.1) THEN
! IF (IROTAT.EQ.1 .AND. ISIMPY.NE.1 .AND. ISIQTY.NE.1) THEN
     CALL AVRXSO(HAUXIJ,AROTLI(3),AKOTLI(3),RHO_PM)
 END IF
!
 IF (CORROT) THEN
    IMSIGN=1
!=======================================================================
!Tr[(Jz.rho)^2] and 2[Jz.rho.Jz]
!=======================================================================

    IF(ISIMPY.EQ.1 .AND. KETA_R.NE.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(3),EKOTSQ(3),HAUXIJ,RHO_PP)
    ELSE IF(ISIMPY.EQ.0 .AND. KETA_R.NE.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(3),EKOTSQ(3),HAUXIJ,RHO_PP,RPM=RHO_PM)
    ELSE IF(ISIMPY.EQ.1 .AND.  KETA_R.EQ.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(3),EKOTSQ(3),HAUXIJ,RHO_PP,APP=AROTPP(:,:,3,:,ICHARG))
    ELSE IF(ISIMPY.EQ.0 .AND.  KETA_R.EQ.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(3),EKOTSQ(3),HAUXIJ,RHO_PP,RPM=RHO_PM,&
                   APP=AROTPP(:,:,3,:,ICHARG), APM=AROTPM(:,:,3,:,ICHARG))
    ELSE
       stop 'wrong ISIMPY or KETA_R in rotavr'
    END IF
!=======================================================================
!Tr[(Jz.kAPPa)(Jz.kAPPa')*] and 2[Jz.kAPPa.Jz]
!=======================================================================
    IF (IPAHFB.EQ.1.AND.KETA_R.GT.0) THEN
       IF(ISIMPY.EQ.1 .AND. KETA_R.NE.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(3),PKOTSQ(3),HAUXIJ,KAI_PM,KAI2PM)
       ELSEIF(ISIMPY.EQ.0 .AND. KETA_R.NE.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(3),PKOTSQ(3),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP)
       ELSEIF(ISIMPY.EQ.1 .AND.  KETA_R.EQ.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(3),PKOTSQ(3),HAUXIJ,KAI_PM,KAI2PM,PPM=PROTPM(:,:,3,:,ICHARG))
       ELSEIF(ISIMPY.EQ.0 .AND.  KETA_R.EQ.2)then
          CALL PAIXSO(MREVER,ISIMPY,KETA_R,PROTSQ(3),PKOTSQ(3),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP,&
                      PPM=PROTPM(:,:,3,:,ICHARG),PPP=PROTPP(:,:,3,:,ICHARG))
       ELSE
         stop 'wrong ISIMPY or KETA_R in rotavr'
       ENDIF
    END IF
!=======================================================================
![Jz^2]
!=======================================================================
    CALL INTLZ2(HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIL(1:NDBASE,1:NDBASE)
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIL(1:NDBASE,1:NDBASE)
!
    CALL INT_LZ(IPHAPP(0,0,3),HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIJ(1:NDBASE,1:NDBASE,0)-HAUXIL(1:NDBASE,1:NDBASE)
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIJ(1:NDBASE,1:NDBASE,1)-HAUXIL(1:NDBASE,1:NDBASE)
!
    DO IBRA=1,LDBASE
       HAUXIJ(IBRA,IBRA,0)=HAUXIJ(IBRA,IBRA,0)+CMPLX(0.25D0,0.0D0)
       HAUXIJ(IBRA,IBRA,1)=HAUXIJ(IBRA,IBRA,1)+CMPLX(0.25D0,0.0D0)
    END DO
!=======================================================================
!Tr[(Jz^2.rho]
!=======================================================================
    CALL AVRXSD(HAUXIJ,DROTSQ(3),DKOTSQ(3),RHO_PP)
 END IF
!=======================================================================
!                         Y-DIRECTION
!=======================================================================
!
!=======================================================================
![Jy]
!=======================================================================
 CALL INT_LY(IPHAPP(0,0,0),HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)=UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)=UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!
 CALL INT_SY(HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIJ(1:NDBASE,1:NDBASE,0)+HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIJ(1:NDBASE,1:NDBASE,1)-HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[Jy.rho]
!=======================================================================
!  IF (IROTAT.EQ.1 .AND. ((ISIMPY.EQ.1 .AND. ISIGNY.NE.1) .OR. (ISIMPY.NE.1 .AND. ISIQTY.NE.1))) THEN
 CALL AVRXSD(HAUXIJ,AROTLI(2),AKOTLI(2),RHO_PP)
!  END IF
!
 IF (CORROT) THEN
    IMSIGN=-1
!=======================================================================
!Tr[(Jy.rho)^2] and 2[Jy.rho.Jy]
!=======================================================================
    IF(ISIMPY.EQ.1 .AND. KETA_R.NE.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(2),EKOTSQ(2),HAUXIJ,RHO_PP)
    ELSEIF(ISIMPY.EQ.0 .AND. KETA_R.NE.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(2),EKOTSQ(2),HAUXIJ,RHO_PP,RPM=RHO_PM)
    ELSEIF(ISIMPY.EQ.1 .AND.  KETA_R.EQ.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(2),EKOTSQ(2),HAUXIJ,RHO_PP,APP=AROTPP(:,:,2,:,ICHARG))
    ELSEIF(ISIMPY.EQ.0 .AND.  KETA_R.EQ.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,KETA_R,EROTSQ(2),EKOTSQ(2),HAUXIJ,RHO_PP,RPM=RHO_PM,&
                   APP=AROTPP(:,:,2,:,ICHARG), APM=AROTPM(:,:,2,:,ICHARG))
    ELSE
       stop 'wrong ISIMPY or KETA_R in rotavr'
    ENDIF
!=======================================================================
!Tr[(Jy.kAPPa)(Jy.kAPPa')*] and 2[Jy.kAPPa.Jy]
!=======================================================================
    IF (IPAHFB.EQ.1.AND.KETA_R.GT.0) THEN
       IF(ISIMPY.EQ.1 .AND. KETA_R.NE.2)then
          CALL PAIXSD(MREVER,ISIMPY,KETA_R,PROTSQ(2),PKOTSQ(2),HAUXIJ,KAI_PM,KAI2PM)
       ELSEIF(ISIMPY.EQ.0 .AND. KETA_R.NE.2)then
          CALL PAIXSD(MREVER,ISIMPY,KETA_R,PROTSQ(2),PKOTSQ(2),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP)
       ELSEIF(ISIMPY.EQ.1 .AND.  KETA_R.EQ.2)then
          CALL PAIXSD(MREVER,ISIMPY,KETA_R,PROTSQ(2),PKOTSQ(2),HAUXIJ,KAI_PM,KAI2PM,PPM=PROTPM(:,:,2,:,ICHARG))
       ELSEIF(ISIMPY.EQ.0 .AND.  KETA_R.EQ.2)then
          CALL PAIXSD(MREVER,ISIMPY,KETA_R,PROTSQ(2),PKOTSQ(2),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP,&
                      PPM=PROTPM(:,:,2,:,ICHARG),PPP=PROTPP(:,:,2,:,ICHARG))
       ELSE
         stop 'wrong ISIMPY or KETA_R in rotavr'
       ENDIF
    END IF
!=======================================================================
![Jy^2]
!=======================================================================
    CALL INTLY2(HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIL(1:NDBASE,1:NDBASE)
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIL(1:NDBASE,1:NDBASE)
!
    CALL INT_LY(IPHAPP(0,0,2),HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIJ(1:NDBASE,1:NDBASE,0)+UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)

    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIJ(1:NDBASE,1:NDBASE,1)-UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!
    DO IBRA=1,LDBASE
       HAUXIJ(IBRA,IBRA,0)=HAUXIJ(IBRA,IBRA,0)+CMPLX(0.25D0,0.0D0)
       HAUXIJ(IBRA,IBRA,1)=HAUXIJ(IBRA,IBRA,1)+CMPLX(0.25D0,0.0D0)
    END DO
!=======================================================================
!Tr[(Jy^2.rho]
!=======================================================================
    CALL AVRXSD(HAUXIJ,DROTSQ(2),DKOTSQ(2),RHO_PP)
 END IF
!
!=======================================================================
!    HERE ALL TERMS ARE COMBINED INTO THE AVERAGE OF THE SQUARE OF
!    THE TOTAL ANGULAR MOMENTUM.
!=======================================================================
 DO KARTEZ=1,NDKART
    TROTSQ(KARTEZ)=DROTSQ(KARTEZ)-EROTSQ(KARTEZ)&
                  +PROTSQ(KARTEZ)&
                  +AROTLI(KARTEZ)**2
    TKOTSQ(KARTEZ)=DKOTSQ(KARTEZ)-EKOTSQ(KARTEZ)&
                  +PKOTSQ(KARTEZ)&
                  +AKOTLI(KARTEZ)**2
 END DO
!
 DO KARTEZ=1,NDKART
    DROTSQ(0)=DROTSQ(0)+DROTSQ(KARTEZ)
    EROTSQ(0)=EROTSQ(0)+EROTSQ(KARTEZ)
    PROTSQ(0)=PROTSQ(0)+PROTSQ(KARTEZ)
    TROTSQ(0)=TROTSQ(0)+TROTSQ(KARTEZ)
    AROTLI(0)=AROTLI(0)+AROTLI(KARTEZ)**2
    PROTLI(0)=PROTLI(0)+PROTLI(KARTEZ)**2
    PROTKI(0)=PROTKI(0)+PROTKI(KARTEZ)**2
!
    DKOTSQ(0)=DKOTSQ(0)+DKOTSQ(KARTEZ)
    EKOTSQ(0)=EKOTSQ(0)+EKOTSQ(KARTEZ)
    PKOTSQ(0)=PKOTSQ(0)+PKOTSQ(KARTEZ)
    TKOTSQ(0)=TKOTSQ(0)+TKOTSQ(KARTEZ)
    AKOTLI(0)=AKOTLI(0)+AKOTLI(KARTEZ)**2
    PKOTLI(0)=PKOTLI(0)+PKOTLI(KARTEZ)**2
    PKOTKI(0)=PKOTKI(0)+PKOTKI(KARTEZ)**2
 END DO
!=======================================================================
!Deallocate work space
!=======================================================================
 DEALLOCATE (HAUXIL)
 DEALLOCATE (HAUXIJ)
 DEALLOCATE (RHO_PP)
 IF(ISIMPY.NE.1) DEALLOCATE (RHO_PM)
 IF(IPAHFB.EQ.1.AND.KETA_R.GT.0) DEALLOCATE (KAI_PM,KAI2PM)
 IF(IPAHFB.EQ.1.AND.KETA_R.GT.0.AND.ISIMPY.NE.1) DEALLOCATE (KAI_PP,KAI2PP)
!=======================================================================
 CALL CPUTIM('ROTAVR',0)
 RETURN
END SUBROUTINE ROTAVR
!=======================================================================
SUBROUTINE RENINE(MIN_QP,IRENIN,IDOGOA,&
                  KSIMTX,JSIMTY,KSIMTZ,KSIGNY,KSIMPY,KSIQTY,&
                         IPAHFB,IROTAT,ITIREP,KREVER,IPNMIX,&
                  IDEVAR,ITERUN,NMUCOU,ISHIFY,NUMCOU,BOUCOU,&
                  ISYMDE,INIROT,INIINV,INIKAR,ISAWAV,IKERNE,&
                  NOSCIL,NUMETA,NFICOU,FILCOU,ICOULI,ICOULO,&
                         IN_FIX,IZ_FIX,ICANRP,I_SLOW,SLOWRP)
!=======================================================================
 USE SAVLEF
 USE SAVRIG
 USE PD_DEN
 USE WAVR_L
 USE HCOULO
 USE OFFDIA, ONLY:RMATCN,ROTWLR,VRVMAT,H_KERN
!=======================================================================
 USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 CHARACTER::FILCOU*68
 INTEGER::MIN_QP,IRENIN,IDOGOA,&
          KSIMTX,JSIMTY,KSIMTZ,KSIGNY,KSIMPY,KSIQTY,&
                 IPAHFB,IROTAT,ITIREP,KREVER,IPNMIX,&
          IDEVAR,ITERUN,NMUCOU,ISHIFY,NUMCOU,ISYMDE,&
          INIROT,INIINV,INIKAR,ISAWAV,IKERNE,&
          NOSCIL,NUMETA,NFICOU,ICOULI,ICOULO,&
          IN_FIX,IZ_FIX,ICANRP,I_SLOW
 REAL   ::BOUCOU,SLOWRP
!=======================================================================
 CHARACTER::NAMEPN*8,FCOUOF*68,CLAB*10
 LOGICAL::PRINIT,CORROT
 INTEGER::IALLOC,JMULMO,JDKART,ICANON,IEXTEN,IFULLD,INROTA,ICHARG,IDKART,KETA_R,&
          LTIMUP,LUPPER,LTOTAL,LDPNMX,MREVER,ISIMPY,ISIGNY,ISIQTY,INPCOU,I_SUCC,&
          ITRUNC,IDIGTA,IMASKN,IREINE,ISKYRM,ISIMTX,ISIMTZ,NLAB
 REAL   ::TSHIFT,H2PARA,A_PARA,REINER,GOAMAS,BAREKE
 COMPLEX::C_ZERO,C_UNIT,EKESKY,BAREKK
 INTEGER,DIMENSION(0:NDISOS)::LSTEMP,LUTEMP,LTTEMP,LDTEMP
 REAL   ,DIMENSION(0:NDKART)::DROTSQ,EROTSQ,PROTSQ,TROTSQ,AROTLI,PROTLI,PROTKI
 COMPLEX,DIMENSION(0:NDSHIF,1:NDKART,0:NDISOS)::AKKERN
 COMPLEX,DIMENSION(0:NDSHIF,0:NDISOS)::OVKERN,EKKERN,EPKERN,DKKERN
 COMPLEX,DIMENSION(0:NDSHIF)::SKKERN,CDKERN,CXKERN,ETKERN,JTKERN
 COMPLEX,DIMENSION(0:NDKART)::DKOTSQ,EKOTSQ,PKOTSQ,TKOTSQ,AKOTLI,PKOTLI,PKOTKI
 COMPLEX,ALLOCATABLE::COUHPP(:,:),COUHPM(:,:)
!=======================================================================
 INTEGER::LDBASE,LDTOTS,LDSTAT,LDUPPE,LDTIMU,NFIPRI,ITESTA,ITESTO,NOITER,NUMITE
 REAL   ::ALPROT,BETROT,GAMROT,ROTATX,ROTATY,ROTATZ,HBMASS,HBMRPA,HBMINP,ROTREN,&
          GOAREN,GOAREI,ROTRIN
!
 REAL   ::COMULT
 COMPLEX::QMUT_N,QMUT_P,QMUT_T,QMUL_N,QMUL_P,QMUL_T
 REAL   ,DIMENSION(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART)::COM_TMP
 COMPLEX,DIMENSION(0:NDMULT,-NDMULT:NDMULT)::QTN_TMP,QTP_TMP,QTT_TMP,QLN_TMP,QLP_TMP,QLT_TMP
!
 COMMON /QMUTTI/ QMUT_N(0:NDMULT,-NDMULT:NDMULT),QMUT_P(0:NDMULT,-NDMULT:NDMULT),QMUT_T(0:NDMULT,-NDMULT:NDMULT)
 COMMON /QMULTI/ QMUL_N(0:NDMULT,-NDMULT:NDMULT),QMUL_P(0:NDMULT,-NDMULT:NDMULT),QMUL_T(0:NDMULT,-NDMULT:NDMULT)
 COMMON /COEMUL/ COMULT(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART)
!
! REAL   ::COEYLM,QUNITS
! INTEGER::NREYLM,KREYLM,NIMYLM,KIMYLM
! REAL   ,DIMENSION(0:NDMULT,0:NDMULT)::COETMP,QUNTMP
! INTEGER,DIMENSION(0:NDMULT,0:NDMULT)::NRETMP
! INTEGER,DIMENSION(0:NDMULT,0:NDMULT,1:NDTERM,0:NDKART)::KRETMP
! INTEGER,DIMENSION(0:NDMULT,1:NDMULT,1:NDTERM,0:NDKART)::KIMTMP
! INTEGER,DIMENSION(0:NDMULT,1:NDMULT)::NIMTMP
! 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 /RANGLE/ ALPROT,BETROT,GAMROT
 COMMON /DIMENS/ LDBASE
 COMMON /DIMSTA/ LDTOTS(0:NDISOS),LDSTAT(0:NDISOS),LDUPPE(0:NDISOS),LDTIMU(0:NDISOS)
 COMMON /ROTATT/ ROTATX,ROTATY,ROTATZ
 COMMON /PLANCK/ HBMASS,HBMRPA,HBMINP
 COMMON /RENROT/ ROTREN(NDKART)
 COMMON /RENGOA/ GOAREN(NDKART),GOAREI(NDKART)
 COMMON /CFIPRI/ NFIPRI
 COMMON /ITERTS/ ITESTA,ITESTO,NOITER,NUMITE
 COMMON /RINROT/ ROTRIN(NDKART)
!=======================================================================
 CALL CPUTIM('RENINE',1)
!=======================================================================
 IALLOC=0
 C_ZERO=CMPLX(0.0D0,0.0D0)
 C_UNIT=CMPLX(1.0D0,0.0D0)
!=======================================================================
 IF (.NOT.ALLOCATED(WARIGH)) THEN
     ALLOCATE (WARIGH(1:NDBASE,1:4*NDSTAT,0:NDSPIN),STAT=IALLOC)
     IF (IALLOC.NE.0) CALL NOALLO('WARIGH','RENINE')
     WARIGH(:,:,:)=C_ZERO
 END IF
 IF (.NOT.ALLOCATED(WALEFT)) THEN
     ALLOCATE (WALEFT(1:NDBASE,1:4*NDSTAT,0:NDSPIN),STAT=IALLOC)
     IF (IALLOC.NE.0) CALL NOALLO('WALEFT','RENINE')
     WALEFT(:,:,:)=C_ZERO
 END IF
!=======================================================================
 ALLOCATE (COUHPP(1:NDBASE,1:NDBASE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('COUHPP','RENINE')
 ALLOCATE (COUHPM(1:NDBASE,1:NDBASE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('COUHPM','RENINE')
!=======================================================================
!
!=======================================================================
!         THIS SUBROUTINE CALCULATES THE INERTIA RENORMALIZATION FACTORS
!=======================================================================
 NAMEPN='DUMMY   '
 PRINIT=.FALSE.
 JMULMO=NMUCOU
 CORROT=.TRUE.
 KETA_R=0
!=======================================================================
!         STORING THE COULOMB MATRIX ELEMENTS AND SO ON.
!         THESE WILL BE RETURNED AFTER THE INERTIA RENORMALIZATION.
!=======================================================================
                  COUHPP(:,:)=HPPCOU(:,:)
 IF (KSIMPY.NE.1) COUHPM(:,:)=HPMCOU(:,:)
!
 COM_TMP=COMULT
 QTN_TMP=QMUT_N
 QTP_TMP=QMUT_P
 QTT_TMP=QMUT_T
 QLN_TMP=QMUL_N
 QLP_TMP=QMUL_P
 QLT_TMP=QMUL_T
! COETMP=COEYLM
! NRETMP=NREYLM
! KRETMP=KREYLM
! NIMTMP=NIMYLM
! KIMTMP=KIMYLM
! QUNTMP=QUNITS

 LSTEMP(:)=LDTOTS(:)
 LUTEMP(:)=LDUPPE(:)
 LTTEMP(:)=LDTIMU(:)
 LDTEMP(:)=LDSTAT(:)
!go to 1
!=======================================================================
!STORE SYMMETRY PARAMETERS,
!THESE VALUES MAY BE CHANGED DEPENDING ON ROTATION AXIS
!=======================================================================
 ISIMPY=KSIMPY
 ISIGNY=KSIGNY
 ISIQTY=KSIQTY
 ISIMTX=KSIMTX
 ISIMTZ=KSIMTZ
 MREVER=KREVER
!
 JDKART = 0
 ICANON = ABS(ICANRP)!calculate under canonical basis
 IFULLD = 0
 IEXTEN = 0!extend mrever=0 wave function to mrever=1
 IF(ICANRP .EQ. -1) IFULLD=1!use full canonical basis when ICANRP=-1
 IF(ABS(ROTATX).GT. 1.0D-8)THEN!WHEN ROTATING ALONG X OR Z AXIX, SYMEMTRY WOULD BE BROKEN
    ISIMPY = 0
    ISIGNY = 0
    ISIMTZ = 0
!
    IF(KREVER.EQ.0)THEN!After rotate along x or z axis ,wave function would have off diagonal part, both simplex block is need.
       IEXTEN=1
    ENDIF
 ENDIF
 IF(ABS(ROTATZ).GT. 1.0D-8)THEN!WHEN ROTATING ALONG X OR Z AXIX, SYMEMTRY WOULD BE BROKEN
    ISIMPY = 0
    ISIGNY = 0
    ISIMTX = 0
!
    IF(KREVER.EQ.0)THEN!After rotate along x or z axis ,wave function would have off diagonal part, both simplex block is need.
       IEXTEN=1
    ENDIF
 ENDIF
 IF(IEXTEN.EQ.1) MREVER=1!use both simpliex block
 IF(IPAHFB.NE.1) ICANON=0!in HF and HFB+canonical case one dosen't need to transform into canonical basis
!=======================================================================
!CHECK INPUT PARAMTER
!=======================================================================
 IF (KREVER.EQ.0 .AND. ICANON.EQ.0 .AND.(ABS(ROTATX) .GT. 1.0d-8 .OR. ABS(ROTATZ) .GT. 1.0d-8))&
    STOP ' ROTATION ALONG X OR Z WITH MREVER=0 IN RENINE'
 IF ((ABS(ROTATX).GT. 1.0d-8).AND.(ABS(ROTATY).GT. 1.0d-8).OR.&
     (ABS(ROTATY).GT. 1.0d-8).AND.(ABS(ROTATZ).GT. 1.0d-8).OR.&
     (ABS(ROTATZ).GT. 1.0d-8).AND.(ABS(ROTATX).GT. 1.0d-8))&
    STOP ' ROTATION ALONG MULTIPLE DIRECTION NOT IMPLEMENTED YET in RENINE'
 IF ((ABS(ROTATX).LE. 1.0d-8).AND.(ABS(ROTATY).LE. 1.0d-8).AND.(ABS(ROTATZ).LE. 1.0d-8))&
    STOP ' NO ROTATION ANGLE DETECTED'
! IF (ICANON.GT.0 .AND. KSIMPY.NE.1) STOP 'ICANON>0 WITH ISIMPY =0 NOT IMPLINENTED YET in RENINE'
! IF (ICANON.EQ.0 .AND. IFULLD.EQ.1) STOP 'ONLY FULL CANONICAL BASIS ARE STORED'
 IF (IRENIN.GT.NDSHIF) STOP ' IRENIN.GT.NDSHIF IN RENINE'
! IF (IEXTEN.EQ.1 .AND. (KREVER.EQ.1 .OR. ICANON.EQ.0)) STOP 'IEXTEN=1 needs KREVER=0 and ICANON>0 in RENINE'
!=======================================================================
!    TREAT COULOMB PART WHEN SYMMETRY IS BROKEN INSIDE THIS ROUTINE
!=======================================================================
 IF(ISIMPY.NE.KSIMPY .OR. ISIGNY.NE.KSIGNY .OR. ISIQTY.NE.KSIQTY )THEN
    CALL MULCOU(ISIMPY,ISIGNY,ISIQTY,NMUCOU)
    FCOUOF=TRIM(FILCOU)//'OFD'
    IF(ITESTA.EQ.NUMITE)THEN!DO CALCULATION ONLY AT THE FIRST ITERATION
      ICOULI=0
      ICOULO=1
      CALL PRECOU(NOSCIL,NUMCOU,NUMETA,BOUCOU,NFICOU,FCOUOF,ICOULI,ICOULO)
    ELSE!READ FROM FILE
      INPCOU=1
      I_SUCC=1
      OPEN(NFICOU,FILE=FCOUOF,STATUS='OLD',FORM='UNFORMATTED')
      CALL RECOUL(NFICOU,INPCOU,I_SUCC,NUMCOU,NUMETA,BOUCOU)
      CLOSE(NFICOU)
    ENDIF
 ENDIF
!=======================================================================
!         HERE BEGIN THE LOOP THE ROTATION ANGLES.
!=======================================================================
 DO INROTA=IRENIN,0,-1
    IF (ABS(ROTATY) .GT. 1.0d-8) THEN
     ALPROT=0.0d0
     BETROT=INROTA*ROTATY
     GAMROT=0.0d0
    ELSEIF (ABS(ROTATX) .GT. 1.0d-8) THEN
     ALPROT=ATAN(1.0d0)*2.0d0
     BETROT=INROTA*ROTATX
     GAMROT=-ALPROT
    ELSE
     ALPROT=0.0d0
     BETROT=0.0d0
     GAMROT=INROTA*ROTATZ
    ENDIF
!
    DO ICHARG=0,NDISOS
       IF(IFULLD.EQ.1)THEN!USE FULL CANONICAL BASIS
          LDTIMU(ICHARG)=LDBASE
          LDUPPE(ICHARG)=LDBASE*(1+MREVER)
          LDSTAT(ICHARG)=LDBASE*(2+MREVER)
          LDTOTS(ICHARG)=LDBASE*2*(1+MREVER)
       ENDIF
!
       IF(IEXTEN.EQ.1)THEN!modify orbits number countings
          LDUPPE(ICHARG)=LDTIMU(ICHARG)*2
          LDSTAT(ICHARG)=LDTIMU(ICHARG)*3
          LDTOTS(ICHARG)=LDTIMU(ICHARG)*4
       ENDIF
!
       LTOTAL=LDTOTS(ICHARG)
       LUPPER=LDUPPE(ICHARG)
       LTIMUP=LDTIMU(ICHARG)
       IF (LTOTAL.GT.4*NDSTAT)STOP ' LTOTAL.GT.4*NDSTAT IN RENINE'
!=======================================================================
!FETCHING THE LEFT WAVE FUNCTIONS FOR THE GIVEN CHARGE
!AND CALCULATE ROTATED RIGHT WAVE FUNCTIONS MULTIPLYING THE INVERTED OVERLAP MATRIX
!USE 'KSIMPY' ,since symmetry in wave function is broken after this routine is called
!USE 'mrever' since both blocks are already stored, by change lduppe and so on, the extension is done
!=======================================================================
       CALL ROTWLR(MREVER,KSIMPY,ICANON,IDOGOA,IPAHFB,INROTA,IEXTEN,&
                   LDBASE,LTIMUP,LUPPER,LTOTAL,ICHARG,OVKERN(INROTA,ICHARG))
!=======================================================================
!        HERE BEGIN CALCULATIONS OF DENSITIES AND OBSERVABLES
!=======================================================================
!=======================================================================
!        RHO MATRIX
!=======================================================================
       CALL DENMAC(MREVER,ICHARG,ISIMPY,IPAHFB,WALEFT,WARIGH,1)
!=======================================================================
!        KAPPA MATRIX
!=======================================================================
       IF (IDOGOA.EQ.0) THEN
          IF (IPAHFB.EQ.1 .AND. CORROT) CALL PAIMAC(ICHARG,ISIMPY,WALEFT,WARIGH,IKERNE)
!=======================================================================
!        TERMS FOR <|J^2|BETA>
!=======================================================================
          CALL ROTAVR(ISIMPY,ISIGNY,ISIQTY,&
                      IPAHFB,ICHARG,KETA_R,IROTAT,MREVER,CORROT,IKERNE,&
                      DROTSQ,EROTSQ,PROTSQ,TROTSQ,AROTLI,PROTLI,PROTKI,&
                      DKOTSQ,EKOTSQ,PKOTSQ,TKOTSQ,AKOTLI,PKOTLI,PKOTKI)
!
          DKKERN(INROTA,ICHARG)=TKOTSQ(JDKART)
          AKKERN(INROTA,1,ICHARG)=AKOTLI(1)
          AKKERN(INROTA,2,ICHARG)=AKOTLI(2)
          AKKERN(INROTA,3,ICHARG)=AKOTLI(3)
      ELSE
          DKKERN(INROTA,ICHARG)=C_ZERO
          AKKERN(INROTA,1,ICHARG)=C_ZERO
          AKKERN(INROTA,2,ICHARG)=C_ZERO
          AKKERN(INROTA,3,ICHARG)=C_ZERO
      END IF
!=======================================================================
!        <|H|BETA>
!=======================================================================
!goto 1
      ISKYRM=ICHARG
      CALL H_KERN(ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,&
                  IPAHFB,MREVER,ICHARG,MIN_QP,IPNMIX,ITIREP,NAMEPN,&
                  PRINIT,IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,&
                  INIKAR,ISAWAV,IKERNE,ISHIFY,JMULMO,NUMCOU,&
                  BOUCOU,IN_FIX,IZ_FIX,ISKYRM,&
                  EKKERN(INROTA,ICHARG),EPKERN(INROTA,ICHARG),&
                  BAREKE,BAREKK,CDKERN(INROTA),CXKERN(INROTA),SKKERN(INROTA))
!goto 1
    END DO!END ICHARG=0,NDISOS
!goto 1
!=======================================================================
!        SKYREM PART OF <|H|Ri>
!=======================================================================
!       CALL TRUTOD
!       CALL ESKYRM(ENESKY,ENEVEN,ENEODD,&
!                   ENREAR,ENE_W0,EEVEW0,EODDW0,EKESKY,LDPNMX,0,1)
!       SKKERN(INSHIF)=EKESKY
!=======================================================================
!        <|H|BETA> and <|J^2|BETA>
!=======================================================================
    ETKERN(INROTA)=EKKERN(INROTA,0)+EKKERN(INROTA,1)&
                  +EPKERN(INROTA,0)+EPKERN(INROTA,1)&
                  +CDKERN(INROTA)+CXKERN(INROTA)+SKKERN(INROTA)
    JTKERN(INROTA)=DKKERN(INROTA,0) +DKKERN(INROTA,1)&
                +2*AKKERN(INROTA,1,0)*AKKERN(INROTA,1,1)&
                +2*AKKERN(INROTA,2,0)*AKKERN(INROTA,2,1)
!
! write(*,*)abs(etkern(inrota)),abs(jtkern(inrota)),abs(OVKERN(INROTA,0)*OVKERN(INROTA,1)),inrota
! NLAB=INT((ROTATX+ROTATY+ROTATZ)/3.14*50)
! WRITE(CLAB,'(I3)')NLAB
! OPEN(1,FILE='pf'//trim(adjustl(clab))//'.dat')
! write(1,*)etkern(inrota)-(EPKERN(INROTA,0)+EPKERN(INROTA,1)),&
!           EPKERN(INROTA,0)+EPKERN(INROTA,1),OVKERN(INROTA,0)*OVKERN(INROTA,1)
! close(1)
! stop
 END DO!END INROTA=IRENIN,0
!=======================================================================
!        BELOW WE DEFINE THE LIPKIN PARAMETERS, WHICH ARE POSITIVE
!=======================================================================
 IF (IDOGOA.EQ.0) THEN
    REINER=(ETKERN(IRENIN)-ETKERN(0))/(JTKERN(IRENIN)-JTKERN(0))
 ELSE
    TSHIFT=IRENIN**2*(ROTATX**2+ROTATY**2+ROTATZ**2)
!
    H2PARA=-2*   (ETKERN(IRENIN)  -ETKERN(0)       )/TSHIFT
    A_PARA=-2*LOG(OVKERN(IRENIN,0)*OVKERN(IRENIN,1))/TSHIFT
!
    GOAMAS=H2PARA/(2*A_PARA**2)
    IF (MREVER.EQ.0) GOAMAS=GOAMAS/4
! MREVER.EQ.0 with shifting or rotating is not allowed yet
! But this option is still kept inside RENMAS and RENINE.
    REINER=GOAMAS
 END IF
!=======================================================================
!   Only keep a given digital of lipkin parameter to accelate converge
!=======================================================================
 ITRUNC=1
 IDIGTA=8
 IF(ITRUNC.EQ.1)THEN
    IMASKN=10**IDIGTA
    IREINE=REINER*IMASKN
    REINER=IREINE/(IMASKN*1.0d0)
 ENDIF
! IF(REINER.GE.0.1d0) REINER=5.0d-2
! IF(REINER.LE.0.0d0) REINER=0.0d0
!=======================================================================
!        BELOW WE DEFINE THE COEFFICIENTS OF <J^2>, WHICH ARE NEGATIVE
!=======================================================================
 DO IDKART = 1,NDKART
    IF (ROTRIN(IDKART).LE. 0.0d0)THEN
        ROTREN(IDKART)=-ROTRIN(IDKART)
    ELSE
        IF (I_SLOW.NE.1) THEN
            ROTREN(IDKART)=ROTREN(IDKART)*SLOWRP-REINER*(1-SLOWRP)
!           ROTREN(IDKART)=-ROTRIN(IDKART)
        ELSE
            ROTREN(IDKART)=-REINER
        ENDIF
    ENDIF
 ENDDO
!
 GOAREI(1)=-GOAMAS
 GOAREI(2)=-GOAMAS
 GOAREI(3)=-GOAMAS
!=======================================================================
! Rstore quantities temporarily changed inside the routine
!=======================================================================
 LDTOTS(:)=LSTEMP(:)
 LDUPPE(:)=LUTEMP(:)
 LDTIMU(:)=LTTEMP(:)
 LDSTAT(:)=LDTEMP(:)
!
!goto 1
 IF(ISIMPY.NE.KSIMPY .OR. ISIGNY.NE.KSIGNY .OR. ISIQTY.NE.KSIQTY )THEN
    CALL MULCOU(KSIMPY,KSIGNY,KSIQTY,NMUCOU)
    OPEN(NFICOU,FILE=FILCOU,STATUS='OLD',FORM='UNFORMATTED')
    INPCOU=1
    I_SUCC=1
    CALL RECOUL(NFICOU,INPCOU,I_SUCC,NUMCOU,NUMETA,BOUCOU)
    CLOSE(NFICOU)
 ENDIF
! CALL MOMETS(ISIMPY,ISIGNY,ISIQTY,QMUL_P,QMUT_P,COMULT,JMULMO,ISHIFY)
                 HPPCOU(:,:)=COUHPP(:,:)
 IF (ISIMPY.NE.1)HPMCOU(:,:)=COUHPM(:,:)
! COEYLM=COETMP
! NREYLM=NRETMP
! KREYLM=KRETMP
! NIMYLM=NIMTMP
! KIMYLM=KIMTMP
! QUNITS=QUNTMP
 COMULT=COM_TMP
 QMUT_N=QTN_TMP
 QMUT_P=QTP_TMP
 QMUT_T=QTT_TMP
 QMUL_N=QLN_TMP
 QMUL_P=QLP_TMP
 QMUL_T=QLT_TMP
1 continue
!=======================================================================
 CALL CPUTIM('RENINE',0)
 RETURN
 END SUBROUTINE RENINE
!=======================================================================
SUBROUTINE LINAVR(ISIMPY,ISIGNY,ISIQTY,&
                  IPAHFB,ICHARG,JETACM,IROTAT,MREVER,COR_CM,IKERNE,&
                  DLINSQ,ELINSQ,PLINSQ,TLINSQ,ALINLI,PLINLI,PLINKI,&
                  DKINSQ,EKINSQ,PKINSQ,TKINSQ,AKINLI,PKINLI,PKINKI)
!=======================================================================
USE MAT_PP
USE MAT_PM
USE MAD_PP
USE MAD_PM
USE KAP_PM
USE KAP2PM
USE KAP_PP
USE KAP2PP
USE MAP_PP
USE MAP_PM
USE ALISTO
USE PLISTO
USE O2AVRG
! USE OSQUARE
!=======================================================================
USE hfodd_sizes
!=======================================================================
IMPLICIT NONE
!=======================================================================
 LOGICAL::COR_CM
 INTEGER::ISIMPY,ISIGNY,ISIQTY,IPAHFB,ICHARG,JETACM,IROTAT,MREVER,IKERNE
 REAL   ::DLINSQ(0:NDKART),ELINSQ(0:NDKART),PLINSQ(0:NDKART),TLINSQ(0:NDKART),&
          ALINLI(0:NDKART),PLINLI(0:NDKART),PLINKI(0:NDKART)
 COMPLEX::DKINSQ(0:NDKART),EKINSQ(0:NDKART),PKINSQ(0:NDKART),TKINSQ(0:NDKART),&
          AKINLI(0:NDKART),PKINLI(0:NDKART),PKINKI(0:NDKART)
!
 INTEGER::IALLOC,IBRA,KARTEZ,IMSIGN
 REAL   ::FACTOR
 COMPLEX:: RESULT,C_ZERO,UNIT_I
 REAL   ,ALLOCATABLE:: HAUXIL(:,:)
 COMPLEX,ALLOCATABLE:: HAUXIJ(:,:,:)
 COMPLEX,ALLOCATABLE:: RHO_PP(:,:,:),RHO_PM(:,:,:),KAI_PM(:,:,:),KAI2PM(:,:,:),KAI_PP(:,:,:),KAI2PP(:,:,:)
!
 INTEGER::IPHAPP,IPHAPM,IPHAMP,IPHAMM,LDBASE,ITESTA,ITESTO,NOITER,NUMITE
 COMMON /T_PHAS/ IPHAPP(0:NDYMAX,0:NDYMAX,0:NDKART),IPHAPM(0:NDYMAX,0:NDYMAX,0:NDKART),&
                 IPHAMP(0:NDYMAX,0:NDYMAX,0:NDKART),IPHAMM(0:NDYMAX,0:NDYMAX,0:NDKART)
 COMMON /DIMENS/ LDBASE
 COMMON /ITERTS/ ITESTA,ITESTO,NOITER,NUMITE
!=======================================================================
 CALL CPUTIM('LINAVR',1)
 IALLOC=0
 !
 C_ZERO=CMPLX(0.0D0,0.0D0)
 UNIT_I=CMPLX(0.0D0,1.0D0)
 FACTOR=1.0D0
!=======================================================================
!allocate space for 2[J.rho.J] and 2[J.kappa.J] matrix
!=======================================================================
 IF(JETACM.EQ.2)THEN
    IF(.NOT.ALLOCATED(ALINPP))THEN
       ALLOCATE(ALINPP(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('ALINPP','LINAVR')
    ENDIF
    IF(.NOT.ALLOCATED(PLINPM) .AND. IPAHFB.EQ.1)THEN
       ALLOCATE(PLINPM(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('PLINPM','LINAVR')
    ENDIF
!
    IF(ISIMPY.NE.1)THEN
       IF(.NOT.ALLOCATED(ALINPM))THEN
          ALLOCATE(ALINPM(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('ALINPM','LINAVR')
       ENDIF
       IF(.NOT.ALLOCATED(PLINPP) .AND. IPAHFB.EQ.1)THEN
          ALLOCATE(PLINPP(NDBASE,NDBASE,1:NDKART,0:NDREVE,0:NDISOS),STAT=IALLOC)
          IF (IALLOC.NE.0) CALL NOALLO('PLINPP','LINAVR')
       ENDIF
    ENDIF
 ENDIF
!=======================================================================
!allocate space for rho and kappa in case they are not allocated yet
!=======================================================================
 IF(.NOT.ALLOCATED(BIG_PP)) THEN
    ALLOCATE(BIG_PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('BIG_PP','LINAVR')
 ENDIF
 IF(.NOT.ALLOCATED(DEN_PP)) THEN
    ALLOCATE(DEN_PP(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('DEN_PP','LINAVR')
 ENDIF
!
 IF(.NOT.ALLOCATED(CHI_PM) .AND. IPAHFB.EQ.1)THEN
    ALLOCATE(CHI_PM(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS), STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('CHI_PM','LINAVR')
 ENDIF
!
 IF(ISIMPY.NE.1)THEN
    IF(.NOT.ALLOCATED(BIG_PM))THEN
        ALLOCATE (BIG_PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
        IF (IALLOC.NE.0) CALL NOALLO('BIG_PM','LINAVR')
    END IF
    IF(.NOT.ALLOCATED(DEN_PM))THEN
        ALLOCATE (DEN_PM(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS),STAT=IALLOC)
        IF (IALLOC.NE.0) CALL NOALLO('DEN_PM','LINAVR')
    ENDIF
!
    IF (.NOT.ALLOCATED(CHI_PP) .AND. IPAHFB.EQ.1)THEN
       ALLOCATE (CHI_PP(1:NDBASE,1:NDBASE,0:NDREVE,0:NDISOS),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('CHI_PP','LINAVR')
    END IF
 ENDIF
!=======================================================================
!allocate work space
!=======================================================================
 ALLOCATE (HAUXIL(1:NDBASE,1:NDBASE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('HAUXIL','LINAVR')
 ALLOCATE (HAUXIJ(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('HAUXIJ','LINAVR')
!
 ALLOCATE (RHO_PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('RHO_PP','LINAVR')
!
 IF(ISIMPY.NE.1)THEN
    ALLOCATE (RHO_PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('RHO_PM','LINAVR')
 ENDIF
!
 IF (IPAHFB.EQ.1.AND.JETACM.GT.0) THEN
    ALLOCATE (KAI_PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('KAI_PM','LINAVR')
    ALLOCATE (KAI2PM(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
    IF (IALLOC.NE.0) CALL NOALLO('KAI2PM','LINAVR')
!
    IF(ISIMPY.NE.1)THEN
       ALLOCATE (KAI_PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('KAI_PP','LINAVR')
       ALLOCATE (KAI2PP(1:NDBASE,1:NDBASE,0:NDREVE),STAT=IALLOC)
       IF (IALLOC.NE.0) CALL NOALLO('KAI2PP','LINAVR')
    ENDIF
 ENDIF
 IF (JETACM.NE.2)THEN
     IF(.NOT.ALLOCATED(BIG_PP)) STOP 'BIG_PP NOT ALLOCATED IN LINAVR'
     IF(.NOT.ALLOCATED(BIG_PM).AND.ISIMPY.NE.1) STOP 'BIG_PM NOT ALLOCATED IN LINAVR'
     IF(.NOT.ALLOCATED(PAI_PM).AND.IPAHFB.EQ.1.AND.JETACM.GT.0) STOP 'PAI_PM NOT ALLOCATED IN LINAVR'
     IF(.NOT.ALLOCATED(PAI2PM).AND.IPAHFB.EQ.1.AND.JETACM.GT.0.AND.IKERNE.EQ.1) STOP 'PAI2PM NOT ALLOCATED IN LINAVR'
     IF(.NOT.ALLOCATED(PAI_PP).AND.IPAHFB.EQ.1.AND.JETACM.GT.0.AND.ISIMPY.NE.1) STOP 'PAI_PP NOT ALLOCATED IN LINAVR'
     IF(.NOT.ALLOCATED(PAI2PP).AND.IPAHFB.EQ.1.AND.JETACM.GT.0.AND.ISIMPY.NE.1.AND.IKERNE.EQ.1) STOP 'PAI2PP NOT ALLOCATED IN LINAVR'
 ENDIF
!=======================================================================
!    ZEROING ARRAYS OF AVERAGE VALUES
!=======================================================================
 DLINSQ(:)=0.0D0
 ELINSQ(:)=0.0D0
 PLINSQ(:)=0.0D0
 TLINSQ(:)=0.0D0
 ALINLI(:)=0.0D0
 PLINLI(:)=0.0D0
 PLINKI(:)=0.0D0
!
 DKINSQ(:)=C_ZERO
 EKINSQ(:)=C_ZERO
 PKINSQ(:)=C_ZERO
 TKINSQ(:)=C_ZERO
 AKINLI(:)=C_ZERO
 PKINLI(:)=C_ZERO
 PKINKI(:)=C_ZERO
!
 IF(JETACM.EQ.2)THEN
                    ALINPP(:,:,:,:,ICHARG)=C_ZERO
    IF(ISIMPY.NE.1) ALINPM(:,:,:,:,ICHARG)=C_ZERO
!
    IF(IPAHFB.EQ.1)THEN
                        PLINPM(:,:,:,:,ICHARG)=C_ZERO
       IF (ISIMPY.NE.1) PLINPP(:,:,:,:,ICHARG)=C_ZERO
    ENDIF
 ENDIF
!
                 RHO_PP(:,:,:)=C_ZERO
 IF(ISIMPY.NE.1) RHO_PM(:,:,:)=C_ZERO
 IF (IPAHFB.EQ.1.AND.JETACM.GT.0) THEN
    KAI_PM(:,:,:)=C_ZERO
    KAI2PM(:,:,:)=C_ZERO
     IF(ISIMPY.NE.1) THEN
        KAI_PP(:,:,:)=C_ZERO
        KAI2PP(:,:,:)=C_ZERO
     ENDIF
 ENDIF
!
 IF (JETACM.EQ.2) THEN !when calculating matrix element,use mixed rho and kappa
    IF(ITESTA.EQ.NUMITE) THEN!initial den, chi at the first iteration
                        DEN_PP(:,:,:,ICHARG)= BIG_PP(:,:,:)
       IF (ISIMPY.NE.1) DEN_PM(:,:,:,ICHARG)= BIG_PM(:,:,:)
       IF (IPAHFB.EQ.1) THEN
                           CHI_PM(:,:,:,ICHARG)=PAI_PM(:,:,:)
          IF (ISIMPY.NE.1) CHI_PP(:,:,:,ICHARG)=PAI_PP(:,:,:)
       ENDIF
    ENDIF
!
                     RHO_PP(:,:,:)=DEN_PP(:,:,:,ICHARG)
    IF (ISIMPY.NE.1) RHO_PM(:,:,:)=DEN_PM(:,:,:,ICHARG)
    IF (IPAHFB.EQ.1) THEN
       KAI_PM(:,:,:)=CHI_PM(:,:,:,ICHARG)
       KAI2PM(:,:,:)=CHI_PM(:,:,:,ICHARG) ! ANALYZE THAT
       IF (ISIMPY.NE.1)THEN
          KAI_PP(:,:,:)=CHI_PP(:,:,:,ICHARG)
          KAI2PP(:,:,:)=CHI_PP(:,:,:,ICHARG)
       ENDIF
    ENDIF
!
 ELSE!when there is no correction or inside RENMAS, use rho and kappa from DENMAC and PAIMAC
                     RHO_PP(:,:,:)=BIG_PP(:,:,:)
    IF (ISIMPY.NE.1) RHO_PM(:,:,:)=BIG_PM(:,:,:)
    IF (IPAHFB.EQ.1.AND.JETACM.GT.0) THEN
        KAI_PM(:,:,:)=PAI_PM(:,:,:)
        IF (IKERNE.EQ.1) THEN
            KAI2PM(:,:,:)=PAI2PM(:,:,:)
        ELSE
            KAI2PM(:,:,:)=PAI_PM(:,:,:)
        ENDIF
        IF (ISIMPY.NE.1) THEN
            KAI_PP(:,:,:)=PAI_PP(:,:,:)
            IF (IKERNE.EQ.1) THEN
                KAI2PP(:,:,:)=PAI2PP(:,:,:)
            ELSE
                KAI2PP(:,:,:)=PAI_PP(:,:,:)
            ENDIF
        ENDIF
    ENDIF
 ENDIF
!=======================================================================
!                         X-DIRECTION
!=======================================================================
!
!=======================================================================
![Px]
!=======================================================================
 CALL INT_PX(IPHAPP(0,0,0),HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)= UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)= UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[Px.rho]
!=======================================================================
 CALL AVRXSD(HAUXIJ,ALINLI(1),AKINLI(1),RHO_PP)
!
 IF (COR_CM) THEN
    IMSIGN=-1
!=======================================================================
!Tr[(Px.rho)^2] and 2[Px.rho.Px]
!=======================================================================

    IF(ISIMPY.EQ.1 .AND. JETACM.NE.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(1),EKINSQ(1),HAUXIJ,RHO_PP)
    ELSEIF(ISIMPY.EQ.0 .AND. JETACM.NE.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(1),EKINSQ(1),HAUXIJ,RHO_PP,RPM=RHO_PM)
    ELSEIF(ISIMPY.EQ.1 .AND.  JETACM.EQ.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(1),EKINSQ(1),HAUXIJ,RHO_PP,APP=ALINPP(:,:,1,:,ICHARG))
    ELSEIF(ISIMPY.EQ.0 .AND.  JETACM.EQ.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(1),EKINSQ(1),HAUXIJ,RHO_PP,RPM=RHO_PM,&
                   APP=ALINPP(:,:,1,:,ICHARG),APM=ALINPM(:,:,1,:,ICHARG))
    ELSE
       stop 'wrong ISIMPY or JETACM in LINAVR'
    ENDIF
!=======================================================================
!Tr[(Px.kappa)(Px.kappa')*] and 2[Px.kappa.Px]
!=======================================================================
    IF (IPAHFB.EQ.1.AND.JETACM.GT.0) THEN
       IF(ISIMPY.EQ.1 .AND. JETACM.NE.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(1),PKINSQ(1),HAUXIJ,KAI_PM,KAI2PM)
       ELSEIF(ISIMPY.EQ.0 .AND. JETACM.NE.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(1),PKINSQ(1),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP)
       ELSEIF(ISIMPY.EQ.1 .AND.  JETACM.EQ.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(1),PKINSQ(1),HAUXIJ,KAI_PM,KAI2PM,PPM=PLINPM(:,:,1,:,ICHARG))
       ELSEIF(ISIMPY.EQ.0 .AND.  JETACM.EQ.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(1),PKINSQ(1),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP,&
                      PPM=PLINPM(:,:,1,:,ICHARG),PPP=PLINPP(:,:,1,:,ICHARG))
       ELSE
         stop 'wrong ISIMPY or JETACM in LINAVR'
       ENDIF
    END IF

!=======================================================================
![Px^2]
!=======================================================================
    CALL INTPX2(IPHAPP(0,0,0),HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIL(1:NDBASE,1:NDBASE)
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[Px^2.rho]
!=======================================================================
    CALL AVRXSD(HAUXIJ,DLINSQ(1),DKINSQ(1),RHO_PP)
 END IF
!=======================================================================
!                         Z-DIRECTION
!=======================================================================
!
!=======================================================================
![Pz]
!=======================================================================
 CALL INT_PZ(IPHAPP(0,0,0),HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)=UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)=UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[Pz.rho]
!=======================================================================
 CALL AVRXSD(HAUXIJ,ALINLI(3),AKINLI(3),RHO_PP)
!
 IF (COR_CM) THEN
    IMSIGN=-1
!=======================================================================
!Tr[(Pz.rho)^2] and 2[Pz.rho.Pz]
!=======================================================================

    IF(ISIMPY.EQ.1 .AND. JETACM.NE.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(3),EKINSQ(3),HAUXIJ,RHO_PP)
    ELSEIF(ISIMPY.EQ.0 .AND. JETACM.NE.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(3),EKINSQ(3),HAUXIJ,RHO_PP,RPM=RHO_PM)
    ELSEIF(ISIMPY.EQ.1 .AND.  JETACM.EQ.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(3),EKINSQ(3),HAUXIJ,RHO_PP,APP=ALINPP(:,:,3,:,ICHARG))
    ELSEIF(ISIMPY.EQ.0 .AND.  JETACM.EQ.2)then
       CALL EXCXSD(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(3),EKINSQ(3),HAUXIJ,RHO_PP,RPM=RHO_PM,&
                   APP=ALINPP(:,:,3,:,ICHARG), APM=ALINPM(:,:,3,:,ICHARG))
    ELSE
       stop 'wrong ISIMPY or JETACM in LINAVR'
    ENDIF
!=======================================================================
!Tr[(Pz.kappa)(Pz.kappa')*] and 2[Pz.kappa.Pz]
!=======================================================================
    IF (IPAHFB.EQ.1.AND.JETACM.GT.0) THEN
       IF(ISIMPY.EQ.1 .AND. JETACM.NE.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(3),PKINSQ(3),HAUXIJ,KAI_PM,KAI2PM)
       ELSEIF(ISIMPY.EQ.0 .AND. JETACM.NE.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(3),PKINSQ(3),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP)
       ELSEIF(ISIMPY.EQ.1 .AND.  JETACM.EQ.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(3),PKINSQ(3),HAUXIJ,KAI_PM,KAI2PM,PPM=PLINPM(:,:,3,:,ICHARG))
       ELSEIF(ISIMPY.EQ.0 .AND.  JETACM.EQ.2)then
          CALL PAIXSD(MREVER,ISIMPY,JETACM,PLINSQ(3),PKINSQ(3),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP,&
                      PPM=PLINPM(:,:,3,:,ICHARG),PPP=PLINPP(:,:,3,:,ICHARG))
       ELSE
         stop 'wrong ISIMPY or JETACM in LINAVR'
       ENDIF
    END IF
!=======================================================================
![Pz^2]
!=======================================================================
    CALL INTPZ2(IPHAPP(0,0,0),HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIL(1:NDBASE,1:NDBASE)
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[(Pz^2.rho]
!=======================================================================
    CALL AVRXSD(HAUXIJ,DLINSQ(3),DKINSQ(3),RHO_PP)
 END IF
!=======================================================================
!                         Y-DIRECTION
!=======================================================================
!
!=======================================================================
![Py]
!=======================================================================
 CALL INT_PY(IPHAPM(0,0,0),HAUXIL,FACTOR)
 HAUXIJ(1:NDBASE,1:NDBASE,0)= UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
 HAUXIJ(1:NDBASE,1:NDBASE,1)=-UNIT_I*HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[Py.rho]
!=======================================================================
 IF (ISIMPY.NE.1 ) THEN
 CALL AVRXSO(HAUXIJ,ALINLI(2),AKINLI(2),RHO_PM)
  END IF
!
 IF (COR_CM) THEN
    IMSIGN=1
!=======================================================================
!Tr[(Py.rho)^2] and 2[Py.rho.Py]
!=======================================================================
    IF(ISIMPY.EQ.1 .AND. JETACM.NE.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(2),EKINSQ(2),HAUXIJ,RHO_PP)
    ELSEIF(ISIMPY.EQ.0 .AND. JETACM.NE.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(2),EKINSQ(2),HAUXIJ,RHO_PP,RPM=RHO_PM)
    ELSEIF(ISIMPY.EQ.1 .AND.  JETACM.EQ.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(2),EKINSQ(2),HAUXIJ,RHO_PP,APP=ALINPP(:,:,2,:,ICHARG))
    ELSEIF(ISIMPY.EQ.0 .AND.  JETACM.EQ.2)then
       CALL EXCXSO(IMSIGN,MREVER,ISIMPY,JETACM,ELINSQ(2),EKINSQ(2),HAUXIJ,RHO_PP,RPM=RHO_PM,&
                   APP=ALINPP(:,:,2,:,ICHARG), APM=ALINPM(:,:,2,:,ICHARG))
    ELSE
       stop 'wrong ISIMPY or JETACM in LINAVR'
    ENDIF
!=======================================================================
!Tr[(Py.kappa)(Py.kappa')*] and 2[Py.kappa.Py]
!=======================================================================
    IF (IPAHFB.EQ.1.AND.JETACM.GT.0) THEN
       IF(ISIMPY.EQ.1 .AND. JETACM.NE.2)then
          CALL PAIXSO(MREVER,ISIMPY,JETACM,PLINSQ(2),PKINSQ(2),HAUXIJ,KAI_PM,KAI2PM)
       ELSEIF(ISIMPY.EQ.0 .AND. JETACM.NE.2)then
          CALL PAIXSO(MREVER,ISIMPY,JETACM,PLINSQ(2),PKINSQ(2),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP)
       ELSEIF(ISIMPY.EQ.1 .AND.  JETACM.EQ.2)then
          CALL PAIXSO(MREVER,ISIMPY,JETACM,PLINSQ(2),PKINSQ(2),HAUXIJ,KAI_PM,KAI2PM,PPM=PLINPM(:,:,2,:,ICHARG))
       ELSEIF(ISIMPY.EQ.0 .AND.  JETACM.EQ.2)then
          CALL PAIXSO(MREVER,ISIMPY,JETACM,PLINSQ(2),PKINSQ(2),HAUXIJ,KAI_PM,KAI2PM,KPP1=KAI_PP,KPP2=KAI2PP,&
                      PPM=PLINPM(:,:,2,:,ICHARG),PPP=PLINPP(:,:,2,:,ICHARG))
       ELSE
         stop 'wrong ISIMPY or JETACM in LINAVR'
       ENDIF
    END IF
!=======================================================================
![Py^2]
!=======================================================================
    CALL INTPY2(IPHAPP(0,0,0),HAUXIL,FACTOR)
    HAUXIJ(1:NDBASE,1:NDBASE,0)=HAUXIL(1:NDBASE,1:NDBASE)
    HAUXIJ(1:NDBASE,1:NDBASE,1)=HAUXIL(1:NDBASE,1:NDBASE)
!=======================================================================
!Tr[(Py^2.rho]
!=======================================================================
    CALL AVRXSD(HAUXIJ,DLINSQ(2),DKINSQ(2),RHO_PP)
 END IF
!
!=======================================================================
!    HERE ALL TERMS ARE COMBINED INTO THE AVERAGE OF THE SQUARE OF
!    THE TOTAL LINEAR MOMENTUM.
!=======================================================================
 DO KARTEZ=1,NDKART
    TLINSQ(KARTEZ)=DLINSQ(KARTEZ)-ELINSQ(KARTEZ)+PLINSQ(KARTEZ) +ALINLI(KARTEZ)**2
    TKINSQ(KARTEZ)=DKINSQ(KARTEZ)-EKINSQ(KARTEZ)+PKINSQ(KARTEZ) +AKINLI(KARTEZ)**2
 END DO

 DO KARTEZ=1,NDKART
    DLINSQ(0)=DLINSQ(0)+DLINSQ(KARTEZ)
    ELINSQ(0)=ELINSQ(0)+ELINSQ(KARTEZ)
    PLINSQ(0)=PLINSQ(0)+PLINSQ(KARTEZ)
    TLINSQ(0)=TLINSQ(0)+TLINSQ(KARTEZ)
    ALINLI(0)=ALINLI(0)+ALINLI(KARTEZ)**2
    PLINLI(0)=PLINLI(0)+PLINLI(KARTEZ)**2
    PLINKI(0)=PLINKI(0)+PLINKI(KARTEZ)**2
!
    DKINSQ(0)=DKINSQ(0)+DKINSQ(KARTEZ)
    EKINSQ(0)=EKINSQ(0)+EKINSQ(KARTEZ)
    PKINSQ(0)=PKINSQ(0)+PKINSQ(KARTEZ)
    TKINSQ(0)=TKINSQ(0)+TKINSQ(KARTEZ)
    AKINLI(0)=AKINLI(0)+AKINLI(KARTEZ)**2
    PKINLI(0)=PKINLI(0)+PKINLI(KARTEZ)**2
    PKINKI(0)=PKINKI(0)+PKINKI(KARTEZ)**2
 END DO
!=======================================================================
!Deallocate work space
!=======================================================================
 DEALLOCATE (HAUXIL)
 DEALLOCATE (HAUXIJ)
 DEALLOCATE (RHO_PP)
 IF(ISIMPY.NE.1) DEALLOCATE (RHO_PM)
 IF(IPAHFB.EQ.1.AND.JETACM.GT.0) DEALLOCATE (KAI_PM,KAI2PM)
 IF(IPAHFB.EQ.1.AND.JETACM.GT.0.AND.ISIMPY.NE.1) DEALLOCATE (KAI_PP,KAI2PP)
!=======================================================================
 CALL CPUTIM('LINAVR',0)
 RETURN
END SUBROUTINE LINAVR
!=======================================================================
SUBROUTINE RENMAS(MIN_QP,IRENMA,IDOGOA, &
                  KSIMTX,JSIMTY,KSIMTZ,KSIGNY,KSIMPY,KSIQTY,&
                         IPAHFB,IROTAT,ITIREP,KREVER,IPNMIX,&
                  IDEVAR,ITERUN,NMUCOU,ISHIFY,NUMCOU,BOUCOU,&
                  ISYMDE,INIROT,INIINV,INIKAR,ISAWAV,IKERNE,&
                  NOSCIL,NUMETA,NFICOU,FILCOU,ICOULI,ICOULO,&
                         IN_FIX,IZ_FIX,ICANTP,I_SLOW,SLOWTP)
!=======================================================================
USE SAVLEF
USE SAVRIG
USE PD_DEN
USE WAVR_L
USE HCOULO
USE OFFDIA, ONLY:LMATCN,LINWLR,VRVMAT,H_KERN
!=======================================================================
USE hfodd_sizes
!=======================================================================
 CHARACTER::FILCOU*68
 CHARACTER::NAMEPN*8,FCOUOF*68
 LOGICAL::PRINIT,COR_CM
 COMPLEX::EKESKY,EKEKIN,EKECOD,EKECOX
!=======================================================================
 COMPLEX::C_ZERO,C_UNIT
 COMPLEX::AVRCMP,BUFCMP
 COMPLEX::EPAIKN,EPAIKP,BAREKK
!=======================================================================
 COMPLEX,ALLOCATABLE::COUHPP(:,:),COUHPM(:,:)
!=======================================================================
 INTEGER,DIMENSION(NDKART)::IREMAS
 COMPLEX,DIMENSION(2)::DETWRK
 COMPLEX,DIMENSION(0:NDSHIF,1:NDKART,0:NDISOS)::AKKERN
 COMPLEX,DIMENSION(0:NDSHIF,0:NDISOS)::OVKERN,EKKERN,EPKERN,DKKERN
 COMPLEX,DIMENSION(0:NDSHIF)::SKKERN,CDKERN,CXKERN,ETKERN
 REAL   ,DIMENSION(0:NDKART)::DLINSQ,ELINSQ,PLINSQ,TLINSQ,ALINLI,PLINLI,PLINKI
 COMPLEX,DIMENSION(0:NDKART)::DKINSQ,EKINSQ,PKINSQ,TKINSQ,AKINLI,PKINLI,PKINKI
 COMPLEX,DIMENSION(0:NDSHIF,1:NDKART)::ETKERS,DKKERS,OVKERS
 REAL   ,DIMENSION(NDKART)::REMASS,GOAMAS
!
 COMPLEX::QMUT_N,QMUT_P,QMUT_T,QMUL_N,QMUL_P,QMUL_T
 REAL   ,DIMENSION(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART)::COM_TMP
 COMPLEX,DIMENSION(0:NDMULT,-NDMULT:NDMULT)::QTN_TMP,QTP_TMP,QTT_TMP,QLN_TMP,QLP_TMP,QLT_TMP
!
 COMMON /QMUTTI/ QMUT_N(0:NDMULT,-NDMULT:NDMULT),QMUT_P(0:NDMULT,-NDMULT:NDMULT),QMUT_T(0:NDMULT,-NDMULT:NDMULT)
 COMMON /QMULTI/ QMUL_N(0:NDMULT,-NDMULT:NDMULT),QMUL_P(0:NDMULT,-NDMULT:NDMULT),QMUL_T(0:NDMULT,-NDMULT:NDMULT)
 COMMON /COEMUL/ COMULT(0:NDMULT,0:NDOSCI,0:NDOSCI,1:NDKART)
!
 COMMON /SVECTO/ BSHIFT(NDKART)
 COMMON /DIMENS/ LDBASE
 COMMON /DIMSTA/ LDTOTS(0:NDISOS),LDSTAT(0:NDISOS),LDUPPE(0:NDISOS),LDTIMU(0:NDISOS)
 COMMON /DISTAN/ DISTAX,DISTAY,DISTAZ
 COMMON /PLANCK/ HBMASS,HBMRPA,HBMINP
 COMMON /RENORM/ HBMREN(NDKART)
 COMMON /RENGOA/ GOAREN(NDKART),GOAREI(NDKART)
 COMMON /CFIPRI/ NFIPRI
 COMMON /ITERTS/ ITESTA,ITESTO,NOITER,NUMITE
!=======================================================================
 CALL CPUTIM('RENMAS',1)
!=======================================================================
 IALLOC=0
 !IF(MREVER.EQ.0)STOP 'MREVER=0 in RENMAS'
!=======================================================================
 IF (.NOT.ALLOCATED(WARIGH)) THEN
     ALLOCATE (WARIGH(1:NDBASE,1:4*NDSTAT,0:NDSPIN),STAT=IALLOC)
     IF (IALLOC.NE.0) CALL NOALLO('WARIGH','RENMAS')
     WARIGH(:,:,:)=CMPLX(0.0D0,0.0D0)
 END IF
 IF (.NOT.ALLOCATED(WALEFT)) THEN
     ALLOCATE (WALEFT(1:NDBASE,1:4*NDSTAT,0:NDSPIN),STAT=IALLOC)
     IF (IALLOC.NE.0) CALL NOALLO('WALEFT','RENMAS')
     WALEFT(:,:,:)=CMPLX(0.0D0,0.0D0)
 END IF
!=======================================================================
 ALLOCATE (COUHPP(1:NDBASE,1:NDBASE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('COUHPP','RENMAS')
 ALLOCATE (COUHPM(1:NDBASE,1:NDBASE),STAT=IALLOC)
 IF (IALLOC.NE.0) CALL NOALLO('COUHPM','RENMAS')
!=======================================================================
!         THIS SUBROUTINE CALCULATES THE MASS RENORMALIZATION FACTORS
!=======================================================================
 C_ZERO=CMPLX(0.0D0,0.0D0)
 C_UNIT=CMPLX(1.0D0,0.0D0)
 SIF=1.0D-8
!
 NAMEPN='DUMMY   '
 PRINIT=.FALSE.
 JMULMO=NMUCOU
 COR_CM=.TRUE.
 JETACM=0
!=======================================================================
!         STORING THE COULOMB MATRIX ELEMENTS. THESE WILL BE RETURNED
!         AFTER THE MASS RENORMALIZATION.
!=======================================================================
                  COUHPP(:,:)=HPPCOU(:,:)
 IF (KSIMPY.NE.1) COUHPM(:,:)=HPMCOU(:,:)
 COM_TMP=COMULT
 QTN_TMP=QMUT_N
 QTP_TMP=QMUT_P
 QTT_TMP=QMUT_T
 QLN_TMP=QMUL_N
 QLP_TMP=QMUL_P
 QLT_TMP=QMUL_T
!
 ICANON=ABS(ICANTP)
 IEXTEN=0
!=======================================================================
!STORE SYMMETRY PARAMETERS,
!THESE VALUES MAY BE CHANGED DEPENDING ON SHIFTING DIRECTION
!=======================================================================
 ISIMPY=KSIMPY
 ISIGNY=KSIGNY
 ISIQTY=KSIQTY
 ISIMTX=KSIMTX
 ISIMTZ=KSIMTZ
 MREVER=KREVER
!
 ISIQTY = 0
 IF(ABS(DISTAX).GT. SIF) THEN
    ISIGNY = 0
    ISIMTZ = 0
 ENDIF
!
 IF(ABS(DISTAY).GT. SIF) THEN
    ISIMPY = 0
    ISIMTX = 0
    ISIMTZ = 0
    IF(KREVER.EQ.0) IEXTEN=1!when breaking time reversal, the stored wave function should be modified to include all simplex blocks
 ENDIF
!
 IF(ABS(DISTAZ).GT. SIF) THEN
    ISIGNY = 0
    ISIMTX = 0
 ENDIF
!
 IF(IEXTEN.EQ.1) MREVER=1
 IF(IPAHFB.NE.1) ICANON=0!
!
!=======================================================================
!CHECK INPUT PARAMTER
!=======================================================================
 IF (KREVER.EQ.0 .AND. ICANON.EQ.0 .AND. ABS(DISTAY).GT. 1.0d-8)&
    STOP 'SHIFTING ALONG Y WITH MREVER=0 IN RENMAS'
! IF (ICANON.GT.0 .AND. KSIMPY.NE.1) STOP 'ICANON>0 WITH ISIMPY =0 NOT IMPLINENTED YET in RENMAS'
 IF (IRENMA.GT.NDSHIF) STOP ' IRENMA.GT.NDSHIF IN RENMAS'
! IF (IEXTEN.EQ.1 .AND. (KREVER.EQ.1 .OR. ICANON.EQ.0)) STOP 'IEXTEN=1 needs KREVER=0 and ICANON>0 in RENMAS'
! IF(ICANON.GT.1)STOP 'ICANON>1 in RENMAS'
!=======================================================================
!    TREAT COULOMB PART WHEN SYMMETRY IS BROKEN INSIDE THIS ROUTINE
!=======================================================================
 IF(ISIMPY.NE.KSIMPY .OR. ISIGNY.NE.KSIGNY .OR. ISIQTY.NE.KSIQTY )THEN
    CALL MULCOU(ISIMPY,ISIGNY,ISIQTY,NMUCOU)
    FCOUOF=TRIM(FILCOU)//'OFD'
    IF(ITESTA.EQ.NUMITE)THEN!DO CALCULATION ONLY AT THE FIRST ITERATION
      ICOULI=0
      ICOULO=1
      CALL PRECOU(NOSCIL,NUMCOU,NUMETA,BOUCOU,NFICOU,FCOUOF,ICOULI,ICOULO)
    ELSE!READ FROM FILE
      INPCOU=1
      I_SUCC=1
      OPEN(NFICOU,FILE=FCOUOF,STATUS='OLD',FORM='UNFORMATTED')
      CALL RECOUL(NFICOU,INPCOU,I_SUCC,NUMCOU,NUMETA,BOUCOU)
      CLOSE(NFICOU)
    ENDIF
 ENDIF
!
 ETKERS(:,:)=C_ZERO
 DKKERS(:,:)=C_ZERO
 OVKERS(:,:)=C_ZERO
!=======================================================================
!         HERE BEGIN THE LOOP THE SHIFT VECTORS.
!=======================================================================!
 IF (IRENMA.GT.NDSHIF) STOP ' IRENMA.GT.NDSHIF IN RENMAS'
 MDKART=NDKART!UPPER LIMIT OF DIRECTION LOOP
 DO INSHIF=IRENMA,0,-1!SHIFT DISTANCE
    IF (INSHIF.EQ.0) MDKART=1!WHEN THERE IS NO SHIFT,ONLY DO 1 TIME
    DO KARTEZ=1,MDKART!X,Y,Z DIRECTION
     BSHIFT(:)=0.0d0
     IF (KARTEZ.EQ.1) BSHIFT(1)=INSHIF*DISTAX
     IF (KARTEZ.EQ.2) BSHIFT(2)=INSHIF*DISTAY
     IF (KARTEZ.EQ.3) BSHIFT(3)=INSHIF*DISTAZ
     IF ((INSHIF.EQ.0).OR.(ABS(BSHIFT(KARTEZ)).GT.SIF)) THEN
       DO ICHARG=0,NDISOS
          IF(IEXTEN.EQ.1)THEN
             LDUPPE(ICHARG)=LDTIMU(ICHARG)*2
             LDSTAT(ICHARG)=LDTIMU(ICHARG)*3
             LDTOTS(ICHARG)=LDTIMU(ICHARG)*4
          ENDIF
          LTIMUP=LDTIMU(ICHARG)
          LUPPER=LDUPPE(ICHARG)
          LTOTAL=LDTOTS(ICHARG)
!=======================================================================
!FETCHING THE LEFT WAVE FUNCTIONS FOR THE GIVEN CHARGE
!AND CALCULATE ROTATED RIGHT WAVE FUNCTIONS MULTIPLYING THE INVERTED OVERLAP MATRIX
!USE 'KSIMPY' HERE,SINCE CHANGE IN WAVEFUNCTION HAPPENS AFTER THIS ROUTINE IS CALLED
!USE 'mrever' since both blocks are already stored, by change lduppe and so on, the extension is done
!=======================================================================
          CALL LINWLR(MREVER,KSIMPY,ICANON,IDOGOA,IPAHFB,INSHIF,IEXTEN,&
                      LDBASE,LTIMUP,LUPPER,LTOTAL,ICHARG,OVKERN(INSHIF,ICHARG))
!=======================================================================
!        RHO MATRIX
!=======================================================================
          CALL DENMAC(MREVER,ICHARG,ISIMPY,IPAHFB,WALEFT,WARIGH,1)
!=======================================================================
!        KAPPA MATRIX
!=======================================================================
          IF (IDOGOA.EQ.0) THEN
             IF (IPAHFB.EQ.1 .AND. COR_CM) CALL PAIMAC(ICHARG,ISIMPY,WALEFT,WARIGH,IKERNE)
!=======================================================================
!        TERMS FOR <|J^2|Ri>
!=======================================================================
             CALL LINAVR(ISIMPY,ISIGNY,ISIQTY,&
                         IPAHFB,ICHARG,JETACM,IROTAT,MREVER,COR_CM,IKERNE,&
                         DLINSQ,ELINSQ,PLINSQ,TLINSQ,ALINLI,PLINLI,PLINKI,&
                         DKINSQ,EKINSQ,PKINSQ,TKINSQ,AKINLI,PKINLI,PKINKI)
!
             DKKERN(INSHIF,ICHARG)=TKINSQ(0)
             AKKERN(INSHIF,1,ICHARG)=AKINLI(1)
             AKKERN(INSHIF,2,ICHARG)=AKINLI(2)
             AKKERN(INSHIF,3,ICHARG)=AKINLI(3)
          ELSE
             DKKERN(INSHIF,ICHARG)=C_ZERO
             AKKERN(INSHIF,1,ICHARG)=C_ZERO
             AKKERN(INSHIF,2,ICHARG)=C_ZERO
             AKKERN(INSHIF,3,ICHARG)=C_ZERO
          END IF

!=======================================================================
!        <|H|Ri>
!=======================================================================
          ISKYRM=ICHARG
          CALL H_KERN(ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,&
                      IPAHFB,MREVER,ICHARG,MIN_QP,IPNMIX,ITIREP,NAMEPN,&
                      PRINIT,IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,&
                      INIKAR,ISAWAV,IKERNE,ISHIFY,JMULMO,NUMCOU,&
                      BOUCOU,IN_FIX,IZ_FIX,ISKYRM,&
                      EKKERN(INSHIF,ICHARG),EPKERN(INSHIF,ICHARG),&
                      BAREKE,BAREKK,CDKERN(INSHIF),CXKERN(INSHIF),SKKERN(INSHIF))
       ENDDO!END ICHARG
!=======================================================================
!        SKYREM PART OF <|H|Ri>
!=======================================================================
!       CALL TRUTOD
!       CALL ESKYRM(ENESKY,ENEVEN,ENEODD,&
!                   ENREAR,ENE_W0,EEVEW0,EODDW0,EKESKY,LDPNMX,0,1)
!       SKKERN(INSHIF)=EKESKY
!=======================================================================
!        <|H|Ri> and <|J^2|Ri>
!=======================================================================
       ETKERS(INSHIF,KARTEZ)=EKKERN(INSHIF,0)+EKKERN(INSHIF,1)&
                            +EPKERN(INSHIF,0)+EPKERN(INSHIF,1)&
                            +CDKERN(INSHIF)+CXKERN(INSHIF)+SKKERN(INSHIF)
       DKKERS(INSHIF,KARTEZ)= DKKERN(INSHIF,0)+DKKERN(INSHIF,1)&
               +2*AKKERN(INSHIF,1,0)*AKKERN(INSHIF,1,1)&
               +2*AKKERN(INSHIF,2,0)*AKKERN(INSHIF,2,1)&
               +2*AKKERN(INSHIF,3,0)*AKKERN(INSHIF,3,1)
       OVKERS(INSHIF,KARTEZ)= OVKERN(INSHIF,0)*OVKERN(INSHIF,1)
!
!     write(*,*)'e',etkers(inshif,kartez),inshif,kartez
!     write(*,*)'p',dkkers(inshif,kartez),inshif,kartez
     ENDIF! IF ((BSHIFT(KARTEZ).NE.0.0D0).OR.(INSHIF.EQ.0))
    ENDDO! End KARTEZ
 END DO ! End INSHIF
!=======================================================================
!        BELOW WE DEFINE THE LIPKIN PARAMETERS, WHICH ARE POSITIVE
!=======================================================================
 REMASS(:)=0.0d0
 GOAMAS(:)=0.0d0
 ABX=ABS(DISTAX)
 ABY=ABS(DISTAY)
 ABZ=ABS(DISTAZ)

 IF(ABX.GT.SIF)THEN!SHIFT ALONG X AXIS EXISTS
!
    IF (IDOGOA.EQ.0) THEN
       REMASS(1)=REAL((ETKERS(IRENMA,1)-ETKERS(0,1))/(DKKERS(IRENMA,1)-DKKERS(0,1)))
    ELSE
       TSHIFT=REAL(IRENMA**2)*DISTAX**2
       H2PARA=-2*REAL(ETKERS(IRENMA,1)-ETKERS(0,1)) /TSHIFT
       A_PARA=-2*LOG(REAL(OVKERS(IRENMA,1)))/TSHIFT
       GOAMAS(1)=H2PARA/(2*A_PARA**2)
!       REMASS(1)=GOAMAS(1)
    END IF
 ENDIF
!
 IF(ABY.GT.SIF)THEN!SHIFT ALONG Y AXIS EXISTS
!
    IF (IDOGOA.EQ.0) THEN
       REMASS(2)=REAL((ETKERS(IRENMA,2)-ETKERS(0,1))/(DKKERS(IRENMA,2)-DKKERS(0,1)))
    ELSE
       TSHIFT=REAL(IRENMA**2)*DISTAY**2
       H2PARA=-2*REAL(ETKERS(IRENMA,2)-ETKERS(0,1)) /TSHIFT
       A_PARA=-2*LOG(REAL(OVKERS(IRENMA,2)))/TSHIFT
       GOAMAS(2)=H2PARA/(2*A_PARA**2)
!       REMASS(2)=GOAMAS(2)
    END IF
 ENDIF
!
 IF(ABZ.GT.SIF)THEN!SHIFT ALONG Y AXIS EXISTS

!
    IF (IDOGOA.EQ.0) THEN
       REMASS(3)=REAL((ETKERS(IRENMA,3)-ETKERS(0,1))/(DKKERS(IRENMA,3)-DKKERS(0,1)))
    ELSE
       TSHIFT=REAL(IRENMA**2)*DISTAZ**2
       H2PARA=-2*REAL(ETKERS(IRENMA,3)-ETKERS(0,1)) /TSHIFT
       A_PARA=-2*LOG(REAL(OVKERS(IRENMA,3)))/TSHIFT
       GOAMAS(3)=H2PARA/(2*A_PARA**2)
!       REMASS(3)=GOAMAS(3)
    END IF
 ENDIF

 IF (MREVER.EQ.0) GOAMAS(:)=GOAMAS(:)/4
! MREVER.EQ.0 with shifting or rotating is not allowed yet
! But this option is still kept inside RENMAS and RENINE.
 IF (IDOGOA.EQ.1) REMASS(:)=GOAMAS(:)
!
 IF(ABX.LT.SIF) THEN!NO SHIFT ALONG X AXIS
    IF(DISTAY.GT.SIF) THEN
       REMASS(1)=REMASS(2)
       GOAMAS(1)=GOAMAS(2)
    ELSE
       REMASS(1)=REMASS(3)
       GOAMAS(1)=GOAMAS(3)
    ENDIF
 ENDIF
!
 IF(ABY.LT.SIF) THEN!NO SHIFT ALONG Y AXIS
    IF(DISTAX.GT.SIF) THEN
       REMASS(2)=REMASS(1)
       GOAMAS(2)=GOAMAS(1)
     ELSE
       REMASS(2)=REMASS(3)
       GOAMAS(2)=GOAMAS(3)
     ENDIF
 ENDIF
!
 IF(ABZ.LT.SIF) THEN!NO SHIFT ALONG Z AXIS
    IF(DISTAX.GT.SIF) THEN
       REMASS(3)=REMASS(1)
       GOAMAS(3)=GOAMAS(1)
     ELSE
       REMASS(3)=REMASS(2)
       GOAMAS(3)=GOAMAS(2)
     ENDIF
 ENDIF
!=======================================================================
!   Only keep a given digital of lipkin parameter to accelate converge
!=======================================================================
 ITRUNC=1
 IDIGTA=8
 IF(ITRUNC.EQ.1)THEN
    IMASK=10**IDIGTA
    IREMAS(:)=NINT(REMASS(:)*IMASK)
    REMASS(:)=IREMAS(:)/(IMASK*1.0d0)
 ENDIF
!
!=======================================================================
!       BELOW WE DEFINE THE COEFFICIENTS OF <P^2>, WHICH ARE NEGATIVE
!=======================================================================
 IF (I_SLOW.NE.1) THEN
     HBMREN(1)=HBMREN(1)*SLOWTP-REMASS(1)*(1-SLOWTP)
     HBMREN(2)=HBMREN(2)*SLOWTP-REMASS(2)*(1-SLOWTP)
     HBMREN(3)=HBMREN(3)*SLOWTP-REMASS(3)*(1-SLOWTP)
 ELSE
     HBMREN(1)=-REMASS(1)
     HBMREN(2)=-REMASS(2)
     HBMREN(3)=-REMASS(3)
 ENDIF
 !
 GOAREN(1)=-GOAMAS(1)
 GOAREN(2)=-GOAMAS(2)
 GOAREN(3)=-GOAMAS(3)
 !write(*,*)'renmas',remass(1),goamas(1)
!=======================================================================
!    RESTORE COULOMB PART WHEN SYMMETRY IS BROKEN INSIDE THIS ROUTINE
!=======================================================================
 IF(ISIMPY.NE.KSIMPY .OR. ISIGNY.NE.KSIGNY .OR. ISIQTY.NE.KSIQTY )THEN
    CALL MULCOU(KSIMPY,KSIGNY,KSIQTY,NMUCOU)
    OPEN(NFICOU,FILE=FILCOU,STATUS='OLD',FORM='UNFORMATTED')
    INPCOU=1
    I_SUCC=1
    CALL RECOUL(NFICOU,INPCOU,I_SUCC,NUMCOU,NUMETA,BOUCOU)
    CLOSE(NFICOU)
 ENDIF
                  HPPCOU(:,:)=COUHPP(:,:)
 IF (ISIMPY.NE.1) HPMCOU(:,:)=COUHPM(:,:)
 COMULT=COM_TMP
 QMUT_N=QTN_TMP
 QMUT_P=QTP_TMP
 QMUT_T=QTT_TMP
 QMUL_N=QLN_TMP
 QMUL_P=QLP_TMP
 QMUL_T=QLT_TMP
!=======================================================================
 CALL CPUTIM('RENMAS',0)
 RETURN
END SUBROUTINE RENMAS
!=======================================================================
