SUBROUTINE INITIAL(TFIX,L,EOLD,ENEW)
USE PARS
USE HRHOMAT
USE GAUSS
USE gaussquad
USE EIGENS
IMPLICIT NONE
 INTEGER:: TFIX,L
 REAL(DP)::EOLD,ENEW
 INTEGER:: I
 REAL(DP)::ALPHA
!
 EOLD=0.0D0;ENEW=0.0D0;RHOS(:,:)=0.0D0
 EIGVAL(:)=0.0D0;EIGVAL_OLD(:)=0.0D0
! 
 CALL ANLS;  ALPHA=L!CALCULTAE LOG(ANL)VALUE
 CALL cpquad(GUSEX,ALPHA,TRIM(QUAD),WX(1:GUSEX),XMESH(1:GUSEX))!FIND MESH POINTS
 CALL TPLUSU(L) !CALCULATE ONE-BODY MATRIX ELEMENTS
 CALL V2TOTAL(L)!CALCULATE TWO-BODY MATRIX ELEMENTS
!
 DO I= 0,TFIX-1 !INITIALIZE RHO MATRIX
   RHOS(I,I)=1.0D0
 ENDDO

END SUBROUTINE INITIAL

SUBROUTINE TMEAN(ISPINC,TUES)!TR((T+U)*RHO)
USE PARS
USE TVUMAT
USE HRHOMAT
IMPLICIT NONE
 INTEGER::I,J,ISPINC
 REAL(DP)::TUES,SUMS
 REAL(DP)::TEMP
 SUMS=0.0D0
 DO I=0,NUSE
    SUMS=SUMS+TUMAT(I,I)*RHOS(I,I)
 ENDDO
 TUES=SUMS!WHEN ONLY TREATING UPPER PART, THE 2 FACTOR ARE ADDED IN ENERGY_NEW SUBROUTINE
END SUBROUTINE TMEAN

SUBROUTINE MEANFIELD!BUILD MEANFIELD FROM VMAT AND RHOS, THEN BUILD H
USE PARS
USE HRHOMAT
USE TVUMAT
IMPLICIT NONE
 INTEGER::BRA1,BRA2,KET1,KET2
 REAL(DP)::SUMS
 V1MAT(:,:)=0.0D0;HMAT(:,:)=0.0D0
 DO BRA1=0,NUSE
!    V1MAT(BRA1,BRA1)=1.0D0
    DO KET1=0,NUSE
       SUMS=0.0D0
       DO BRA2=0,NUSE
          DO KET2=0,NUSE
             SUMS=SUMS+VMAT(BRA2,BRA1,KET2,KET1)*RHOS(KET2,BRA2)
             !SUMS=SUMS+8.0D0*RHOS(KET2,BRA2)
          ENDDO
       ENDDO
       V1MAT(BRA1,KET1)=SUMS
    ENDDO
 ENDDO
 HMAT(:,:)=TUMAT(:,:)+V1MAT(:,:)
END SUBROUTINE MEANFIELD

SUBROUTINE RHOMAT(TFIX,ISPINC,RMIX)!CONSTRUCT RHO MATRIX
USE PARS
USE HRHOMAT
USE EIGENS
IMPLICIT NONE
 INTEGER::TFIX,ISPINC
 REAL(DP)::RMIX
 REAL(DP)::RHOS_NEW(0:NMAX,0:NMAX)
 INTEGER::I,IORDER,BRA,KET
 REAL(DP)::SUMS
!CONSTRUCT RHOS
! RHOS_OLD(:,:)=RHOS(:,:)
! write(*,*)rmix!EIGVEC(0,0)*HMAT(0,0)
 CALL REORDR!determine the order of eigenvectors
 RHOS_NEW(:,:)=0.0D0
 DO BRA=0,NUSE
    DO KET=0,NUSE
       SUMS=0.0D0      
       DO I=0,TFIX-1
          !IORDER=EIGORD(I)
          IORDER=I
          SUMS=SUMS+EIGVEC(BRA,IORDER)*EIGVEC(KET,IORDER)
       ENDDO
!       RHOS(BRA,KET)=SUMS
        RHOS_NEW(BRA,KET)=SUMS!*(ISPINC+1)!|PHI_i> HAS SPIN LABEL WHEN ONLY TREATING UPPER BLOCK
    ENDDO
 ENDDO
! RHOS(:,:)=RMIX*RHOS(:,:)+(1-RMIX)*RHOS_OLD(:,:)
 RHOS(:,:)=RMIX*RHOS_NEW(:,:)+(1-RMIX)*RHOS(:,:)
END SUBROUTINE RHOMAT


