!
!===============================================================================
!     Module for Adiabatic Time-Dependent Hartree-Fock (ATDHF) calculations
!     ATDHF calculation is triggered if IADBAT = 1
!===============================================================================
      SUBROUTINE Set_version_hfodd_adiabatic
      !
      USE hfodd_sizes
      !
      CHARACTER NAMMOD*16
      INTEGER   MODUVE,MODSET
      !
      COMMON /VERMOD/ NAMMOD(NDMODU),MODUVE(NDMODU),MODSET(NDMODU)
      !
      NAMMOD(17) = 'hfodd_adiabatic '
      MODUVE(17) = 6
      MODSET(17) = 17
      !
      END SUBROUTINE Set_version_hfodd_adiabatic

!===============================================================================
! Call DENSHF (ICH = 0), SKFILD (ICH = 1), INTEGH (ICH = 2)
! Parameters of these subroutines should be pushed into DSIPARAI & DSIPARAR & SIPARAC
!
      subroutine DENSHF_SKFILD_INTEGH(ICHARG, ICH)

      ! Integer INPUT parameters of DENSHF & SKFILD & INTEGH
      common /DSIPARAI/                                                 &
     &                    ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,    &
     &                    IPAHFB,MREVER,MIN_QP,IPNMIX,ITPNMX,ITIREP,    &
     &                    IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,INIKAR,    &
     &                    ISAWAV,IKERNE,IMOVAX,ITILAX,ITISAX,NO_ORB,    &
     &                    IROTAT,MPAHFB,IPOTHO,ISCHIF,IOCONT,NMUCON,    &
     &                    NSICON,I_YUKA,I_GOGA,I_REGA,I_SEPA,ITWOLI,    &
     &                    JETACM,KETA_R,ICOUDI,ICOUEX,LIPKIN,LIPKIP,    &
     &                    IF_RPA,ILIPON,ILIPOP,NEWGOG,ITWCEN,LTWCEN

      common /DSIPARAC/ NAMEPN
      CHARACTER NAMEPN*8
      ! Real INPUT parameters of DENSHF & SKFILD & INTEGH
      common /DSIPARAR/                                                 &
     &             AROTLT,ALINLT,PRINIT,XOLDEV,XOLDOD,EFER2N,EFER2P

      ! There is a security check in DENSHF, requiring ITPNMX = ICHARG.
      ! Searching string 'ICHARG<>ITPNMX IN DENSHF' in main source file
      ! for further details
      ITPNMX = ICHARG

      if(ICHARG .eq. 0) NAMEPN = 'NEUTRONS'
      if(ICHARG .eq. 1) NAMEPN = ' PROTONS'

      XOLDEV = 0.d0
      XOLDOD = 0.d0

! ---
      IF (ICH .EQ. 0) THEN
         CALL ZEDENS(ICHARG)
         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,         &
     &                                           ITWCEN,LTWCEN)
      ENDIF
! ---
      IF (ICH .EQ. 1) THEN
         CALL ZEFILD(0)
         CALL ZEFILD(1)
         CALL SKFILD(XOLDEV,XOLDOD,IPNMIX)
      ENDIF

! ---
      IF (ICH .EQ. 2) THEN
         CALL INTEGH(IMOVAX,ITILAX,ITISAX,NO_ORB,                       &
     &               IPNMIX,ITPNMX,IROTAT,ICHARG,ISIMPY,                &
     &                      MPAHFB,IPOTHO,ISCHIF,IOCONT,                &
     &                                    NMUCON,NSICON,                &
     &                      I_YUKA,I_GOGA,I_REGA,I_SEPA,                &
     &                             ITWOLI,JETACM,ALINLT,                &
     &                                    KETA_R,AROTLT,                &
     &                                    ICOUDI,ICOUEX,                &
     &                                    LIPKIN,LIPKIP,                &
     &                                    EFER2N,EFER2P,                &
     &                                           IF_RPA,                &
     &                                    ILIPON,ILIPOP,                &
     &                                           NEWGOG,                &
     &                                           ITWCEN)
      ENDIF

      endsubroutine DENSHF_SKFILD_INTEGH

!===============================================================================
! This subroutine should be called before atdhf_implementation
! It pushes input parameters of DENSHF, SKFILD, INTEGH into common block
! DSIPARAI(for integer), DSIPARAR(for real).
!
      subroutine PUSH_DSI_PARA(                                         &
     &                    ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,    &
     &                    IPAHFB,MREVER,MIN_QP,IPNMIX,ITPNMX,ITIREP,    &
     &                    IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,INIKAR,    &
     &                    ISAWAV,IKERNE,IMOVAX,ITILAX,ITISAX,NO_ORB,    &
     &                    IROTAT,MPAHFB,IPOTHO,ISCHIF,IOCONT,NMUCON,    &
     &                    NSICON,I_YUKA,I_GOGA,I_REGA,I_SEPA,ITWOLI,    &
     &                    JETACM,KETA_R,ICOUDI,ICOUEX,LIPKIN,LIPKIP,    &
     &                    IF_RPA,ILIPON,ILIPOP,NEWGOG,ITWCEN,LTWCEN,    &
     &                                                       NAMEPN,    &
     &             AROTLT,ALINLT,PRINIT,XOLDEV,XOLDOD,EFER2N,EFER2P)

      !use ADBATC, only : IADB_DEBUG
      PARAMETER (NUMB_DSI_PARA_INT = 48, NUMB_DSI_PARA_REAL = 7)
      COMMON /DSIPARAI/ DSI_PARA_INT(NUMB_DSI_PARA_INT)
      COMMON /DSIPARAR/ DSI_PARA_REAL(NUMB_DSI_PARA_REAL)
      common /DSIPARAC/ DSI_PARA_CHAR
      CHARACTER NAMEPN*8, DSI_PARA_CHAR*8
      integer DSI_PARA_INT
      real    DSI_PARA_REAL

      DSI_PARA_INT = (/                                                 &
     &                    ISIMTX,JSIMTY,ISIMTZ,ISIGNY,ISIMPY,ISIQTY,    &
     &                    IPAHFB,MREVER,MIN_QP,IPNMIX,ITPNMX,ITIREP,    &
     &                    IDEVAR,ITERUN,ISYMDE,INIROT,INIINV,INIKAR,    &
     &                    ISAWAV,IKERNE,IMOVAX,ITILAX,ITISAX,NO_ORB,    &
     &                    IROTAT,MPAHFB,IPOTHO,ISCHIF,IOCONT,NMUCON,    &
     &                    NSICON,I_YUKA,I_GOGA,I_REGA,I_SEPA,ITWOLI,    &
     &                    JETACM,KETA_R,ICOUDI,ICOUEX,LIPKIN,LIPKIP,    &
     &                    IF_RPA,ILIPON,ILIPOP,NEWGOG,ITWCEN,LTWCEN/)

      DSI_PARA_REAL = (/AROTLT,ALINLT,PRINIT,XOLDEV,XOLDOD,EFER2N,EFER2P/)

      DSI_PARA_CHAR = NAMEPN

      if(1 .eq. 1) then
         write(*,*) 'Checking the parameters pushed into DSIPARAI & DSIPARAR & DSI_PARA_CHAR'
         write(*,'(6i5)') DSI_PARA_INT
         write(*,'(7f12.5)') DSI_PARA_REAL
         write(*,'(a)') DSI_PARA_CHAR
      endif

      endsubroutine PUSH_DSI_PARA

!-------------------------------------------------------------------------------
! Module containing data from ground state solver
      module hfodd_gs_data
      USE hfodd_sizes
      USE ALLWAV
      USE WAVR_L
      USE MAT_PP
      USE MAT_PM
      ! == IADB_DEBUG
      USE HE_DEN
      USE VE_FLD
      USE MAD_PP
      USE MAD_PM
      USE WD_FLD
      ! == IADB_DEBUG

      COMMON                                                            &
     &       /SPNUMS/ NUMBSP(0:NDREVE,0:NDISOS)
      COMMON                                                            &
     &       /SIZOCC/ M_OCCU(1:2*NDSTAT,0:NDISOS),                      &
     &                M_OCVA(1:2*NDSTAT,0:NDISOS)
      COMMON                                                            &
     &       /PROCCU/ V_OCCU(1:4*NDSTAT)
      COMMON                                                            &
     &       /DIMENS/ LDBASE
      COMMON                                                            &
     &       /SPLEVS/ SPENER(1:2*NDSTAT,0:NDISOS),                      &
     &                SPOCCU(1:2*NDSTAT,0:NDISOS),                      &
     &                SPALIG(1:2*NDSTAT,0:NDISOS),                      &
     &                SPSPAL(1:2*NDSTAT,0:NDISOS),                      &
     &                SPPARI(1:2*NDSTAT,0:NDISOS),                      &
     &                SPSIMP(1:2*NDSTAT,0:NDISOS)
      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                                                            &
     &       /PLANCK/ HBMASS,HBMRPA,HBMINP
      COMMON                                                            &
     &       /CCPSKY/ CRHO_T,CRHO_S,CRHODT,CRHODS,                      &
     &                CLPR_T,CLPR_S,                                    &
     &                CTAU_T,CTAU_S,                                    &
     &                CSCU_T,CSCU_S,                                    &
     &                CDIV_T,CDIV_S,                                    &
     &                CSPI_T,CSPI_S,CSPIDT,CSPIDS,                      &
     &                CLPS_T,CLPS_S,                                    &
     &                CCUR_T,CCUR_S,                                    &
     &                CKIS_T,CKIS_S,                                    &
     &                CROT_T,CROT_S,                                    &
     &                CSCT_T,CSCT_S,                                    &
     &                CKIT_T,CKIT_S,                                    &
     &                CSPT_T,CSPT_S

      ! --- common block added by ATDHF ---
      COMMON                                                            &
     &       /SVOCCU/ SV_OCC(1:4*NDSTAT,NDISOS)
      endmodule hfodd_gs_data

MODULE ADBATC
!===============================================================================
! %%- Data for ATDHF calculations
! ## parameters, provided by the main call
!     -- max_iteration: maximal iterations
!     -- atd_precision: numerical tolerance for ATDHF calculation
!     -- M, N: dimension of particle/hole configurations
!     -- elp, elh: particle/hole single particle energy
!
! ## data blocks
!     -- rho1: time-odd correction of density
!     -- gamma1: time-odd correction of mean field
!     -- rho_derivative: derivative of density with respect to the collective
!                        variables
!     -- eph: difference of particle/hole energy
!     -- S: sigular value of \tilde\rho1
!     -- U, V: eigenvector of \tilde\rho1
!     -- mass: collective mass
!     -- vec : eigenvector of rho1
! ##
      implicit none

      integer max_iteration
      real atd_precision
      integer M, N
      real,    allocatable :: elp(:,:), elh(:,:)
      complex, allocatable :: rho1(:,:,:)
      complex, allocatable :: gamma1(:,:,:)
      complex, allocatable :: rho_derivative(:,:,:)
      real,    allocatable :: eph(:,:,:)
      real,    allocatable :: S(:,:)
      complex, allocatable :: U(:,:,:), V(:,:,:)
      real  mass
      complex, allocatable :: vec(:,:,:)

! --- Auxiliary arrays ...
!     -- NUMBHE: number of hole state
!     -- NUMBPE: number of particle state
      integer  NUMBHE(0:1), NUMBPE(0:1)
      complex, allocatable :: BAVOCC(:,:,:)
! --- Debug mode flag
      integer  IADB_DEBUG

CONTAINS

!===============================================================================
! --- Implementation of adiabatic time-dependent Hartree-Fock method
      subroutine atdhf_implementation
      ! == IADB_DEBUG
      use hfodd_gs_data
      implicit none

      integer i
      real mass_old
      ! == IADB_DEBUG
      integer ICHARG

! --- Initialization
      call atdhf_data_init(0)

!     DO ICHARG=0,NDISOS
!        if(IADB_DEBUG .eq. 1) then
!           print*, ' --- AT the beginning of ATDHF ---'
!           print*, 'ICHARG = ', ICHARG
!           print*, 'Maxval of DEN_PP:', maxval(abs(DEN_PP(:,:,:,ICHARG)))
!           print*, 'Maxval of DEN_PM:', maxval(abs(DEN_PM(:,:,:,ICHARG)))
!           print*, ' --- --- --- ---'
!           print*, 'NDXHRM, NDYHRM, NDZHRM', NDXHRM, NDYHRM, NDZHRM
!           print*, 'Maxval of DE_RHO:', maxval(abs(DE_RHO(:,:,:,ICHARG)))
!           print*, 'Maxloc of DE_RHO:', maxloc(abs(DE_RHO(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_TAU:', maxval(abs(DE_TAU(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_LPR:', maxval(abs(DE_LPR(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_DIV:', maxval(abs(DE_DIV(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_SPI:', maxval(abs(DE_SPI(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_KIS:', maxval(abs(DE_KIS(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_GRR:', maxval(abs(DE_GRR(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_LPS:', maxval(abs(DE_LPS(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_ROS:', maxval(abs(DE_ROS(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_ROC:', maxval(abs(DE_ROC(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_CUR:', maxval(abs(DE_CUR(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_KIF:', maxval(abs(DE_KIF(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_SCU:', maxval(abs(DE_SCU(:,:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_DES:', maxval(abs(DE_DES(:,:,:,:,:,ICHARG)))
!        endif
!     ENDDO ! For ICHARG

      if(IADB_DEBUG .eq. 1) then
         print*, ' --- AT the beginning of ATDHF ---'
         print*, 'Maxval of VE_MAS:', maxval(abs(VE_MAS(:,:,:,0:1,:)))
         print*, 'Maxval of VE_CEN:', maxval(abs(VE_CEN(:,:,:,0:1,:)))
         print*, 'Maxval of VE_SPT:', maxval(abs(VE_SPT(:,:,:,0:1)))
         print*, 'Maxval of VE_KIS:', maxval(abs(VE_KIS(:,:,:,:,0:1)))
         print*, 'Maxval of VE_SPI:', maxval(abs(VE_SPI(:,:,:,:,0:1)))
         print*, 'Maxval of VE_CUR:', maxval(abs(VE_CUR(:,:,:,:,0:1)))
         print*, 'Maxval of VE_KIT:', maxval(abs(VE_KIT(:,:,:,:,0:1)))
         print*, 'Maxval of VE_SOR:', maxval(abs(VE_SOR(:,:,:,:,:,0:1)))
         print*, 'Maxval of WD_CEN:', maxval(abs(WD_CEN(:,:,:,0:1)))
         print*, 'Maxval of WD_TAU:', maxval(abs(WD_TAU(:,:,:,0:1)))
         print*, 'Maxval of WD_SOR:', maxval(abs(WD_SOR(:,:,:,:,:,0:1)))
      endif

! --- Cranking approximation
      call cranking_mass

! --- Main iteration
      mass_old = -9999
      !rho1 = rho_derivative/eph
      !rho1 = 1.d0/eph * rho_derivative
      gamma1 = 0.d0
      call gamma1_to_rho1
      do i = 1, max_iteration
          call rho1_to_gamma1
          call gamma1_to_rho1
          !call broyden_mix_rho1(i)
          call calculate_mass
          if(abs(mass_old - mass) < atd_precision .and. i .gt. 5) exit
          if(IADB_DEBUG .eq. 1) then
             write(*,'(a,i4,a,f12.5,a,f12.5)') &
             &'# iter = ', i, ' mass = ', mass, ' residual ',mass_old - mass
          endif
          mass_old = mass
      enddo
      if(i > max_iteration) stop 'Iteration failed in ATDHF_IMPLEMENTATION!'

      if(IADB_DEBUG .eq. 1) then
         write(*,'(a,i5,f12.5)') 'Calculation converged after :',i,abs(mass_old - mass)
      endif
! --- Generating output & clean up
      call atdhf_output
      call atdhf_data_init(1)

      endsubroutine atdhf_implementation

!===============================================================================
! initialization of ATDHF calculations, allocating and preparing necessary arrays
      subroutine atdhf_data_init(ich)
      use hfodd_gs_data
      implicit none

      integer ich
      integer ISTATP, ISTATH
      integer NUMBPX, NUMBHX
      integer ICHARG, ISTATE, iccount, jccount

      IADB_DEBUG = 1

! --- Setup iteration precision & maximum
      atd_precision = 1.d-6
      max_iteration = 199

! --- Determine the dimension of arrays according to PH configuration
      call phsize(NUMBPX, NUMBHX)
      M = NUMBPX
      N = NUMBHX

      if(IADB_DEBUG .eq. 1 .and. ICH .eq. 0) then
         write(*,*) 'checking dimension of PH space'
         write(*,'(2i5)') M, N
         write(*,*) 'for neutron and proton:'
         write(*,'(4i5)') NUMBHE(0),NUMBPE(0),NUMBHE(1),NUMBPE(1)
      endif

      if(ich == 0) then
          ! --- allocations
          allocate(elp(M,0:NDISOS), elh(N,0:NDISOS))
          allocate(rho1(M,N,0:NDISOS), gamma1(M,N,0:NDISOS))
          allocate(S(N,0:NDISOS), U(M,N,0:NDISOS), V(N,N,0:NDISOS))
          allocate(rho_derivative(M,N,0:NDISOS), eph(M,N,0:NDISOS))
          allocate(BAVOCC(1:NDBASE,1:2*NDSTAT,0:NDREVE))
          allocate(vec(M+N,N,0:NDISOS))

          ! --- reset skfild, remove the effect of linear mixing
          Call DENSHF_SKFILD_INTEGH(-1,1)

          ! --- extract particle/hole energy
          DO ICHARG=0,NDISOS
             iccount = 0; jccount = 0
             DO ISTATE=1,NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
                if(M_OCVA(ISTATE,ICHARG) .eq. 1) then
                   iccount = iccount + 1
                   elh(iccount,ICHARG) = SPENER(ISTATE,ICHARG)
                else
                   jccount = jccount + 1
                   elp(jccount,ICHARG) = SPENER(ISTATE,ICHARG)
                endif
             ENDDO
             if(iccount .ne. NUMBHE(ICHARG)) stop 'ERROR in extracting PH energy...'
          ENDDO

          ! --- Prepare the single particle energy difference
          eph = 0.d0
          DO ICHARG=0,NDISOS
             do ISTATP = 1, NUMBPE(ICHARG)
                do ISTATH = 1, NUMBHE(ICHARG)
                   eph(ISTATP,ISTATH,ICHARG) = elp(ISTATP,ICHARG) - elh(ISTATH,ICHARG)
                   if( abs(eph(ISTATP,ISTATH,ICHARG)) .lt. 1.0d-6 ) then
                      stop 'Unsolved Crossing detected in EPH!'
                   endif
                enddo
             enddo
          ENDDO

          if(IADB_DEBUG .eq. 1) then
             print*, 'checking value of EPH, the abs maxval and minval are:'
             write(*,'(2f12.5)') maxval(abs(eph)), minval(abs(eph))
             write(*,'(a,f12.5)') 'for neutron', maxval(abs(eph(1:NUMBPE(0),1:NUMBHE(0),0)))
             write(*,'(a,f12.5)') 'for proton ', maxval(abs(eph(1:NUMBPE(1),1:NUMBHE(1),1)))
          endif

          ! --- calculate the derivative of density
          call calculate_rho_derivative

          ! --- Switch off time-even parameters
          HBMRPA = 0.d0
          if(IADB_DEBUG .eq. 1 .and. .true.) then
             CRHO_T = 0.d0; CRHO_S = 0.d0; CRHODT = 0.d0; CRHODS = 0.d0
             CLPR_T = 0.d0; CLPR_S = 0.d0;
             CTAU_T = 0.d0; CTAU_S = 0.d0;
             CSCU_T = 0.d0; CSCU_S = 0.d0;
             CDIV_T = 0.d0; CDIV_S = 0.d0;
             CSCT_T = 0.d0; CSCT_S = 0.d0;
             if(IADB_DEBUG .eq. 1) print*, 'C_ parameters of the time-even terms are set to 0'
          endif

          if(IADB_DEBUG .eq. 1 .and. .false.) then
             CSPI_T = 0.d0; CSPI_S = 0.d0; CSPIDT = 0.d0; CSPIDS = 0.d0;
             CLPS_T = 0.d0; CLPS_S = 0.d0;
             CCUR_T = 0.d0; CCUR_S = 0.d0;
             CKIS_T = 0.d0; CKIS_S = 0.d0;
             CROT_T = 0.d0; CROT_S = 0.d0;
             CKIT_T = 0.d0; CKIT_S = 0.d0;
             CSPT_T = 0.d0; CSPT_S = 0.d0;
             if(IADB_DEBUG .eq. 1) print*, 'C_ parameters of the time-odd terms are set to 0'
          endif
      endif

      if(ich == 1) then
         ! --- recover WAVOCC & V_OCCU, and everything ...
         DO ICHARG=0,NDISOS
            WAVOCC = SAVOCC(:,:,:,ICHARG)
            V_OCCU = SV_OCC(:,ICHARG)
            Call DENSHF_SKFILD_INTEGH(ICHARG,0)
         ENDDO
         Call DENSHF_SKFILD_INTEGH(-1,1)
         DO ICHARG=0,NDISOS
            Call DENSHF_SKFILD_INTEGH(ICHARG,2)
         ENDDO

         ! --- free memory space...
         deallocate(rho1, gamma1, rho_derivative, eph, S, U, V)
         deallocate(elp,elh,vec,BAVOCC)
      endif

      endsubroutine atdhf_data_init

!===============================================================================
! after applying SVD to rho1, fed DENSHF with eigen vector of rho1, extracting
! Gamma1 by calling SKFILD and INTEGH
      subroutine rho1_to_gamma1
      use hfodd_gs_data, only : NDISOS
      implicit none

      integer ICHARG, NDIMPE, NDIMHE

      vec = (0.d0, 0.d0)

      DO ICHARG=0,NDISOS
         NDIMPE = NUMBPE(ICHARG)
         NDIMHE = NUMBHE(ICHARG)
         if(IADB_DEBUG .eq. 1) then
            print*, ICHARG,'PH dim:', NDIMPE, NDIMHE
            write(*,'(a,2f12.5)') 'rho1(1,1,ICHARG) = ', rho1(1,1,ICHARG)
            write(*,'(a,2f12.5)') 'rho1(3,4,ICHARG) = ', rho1(3,4,ICHARG)
            call PRINT_MATRIX_Z( 'rho1', NUMBPE(ICHARG), NUMBHE(ICHARG), rho1(:,:,ICHARG), M )
         endif

         ! Performing singular value decomposition to rho1
         !call zcallsvd(M,N,rho1(1,1,ICHARG),U(1,1,ICHARG),V(1,1,ICHARG),S(1,ICHARG))
         call zcallsvd(NDIMPE,NDIMHE,rho1(1,1,ICHARG),  M,              &
     &                 U(1,1,ICHARG), M, V(1,1,ICHARG), N,              &
     &                 S(1,ICHARG), N)
         if(IADB_DEBUG .eq. 1) then
            write(*,'(a,i5,a,i5)') 'Singular value:',ICHARG,' Dimension',NDIMHE
            write(*,'(4f12.5)') S(:,ICHARG)
         endif

         ! Preparing eigenvectors of rho1 with positive energy
         vec(1:NDIMPE,1:NDIMHE,ICHARG) = U(1:NDIMPE,1:NDIMHE,ICHARG)
         vec(NDIMPE+1:NDIMPE+NDIMHE,1:NDIMHE,ICHARG) = V(1:NDIMHE,1:NDIMHE,ICHARG)
      ENDDO

      if(IADB_DEBUG .eq. 1) then
         print*, 'checking value of U, the abs maxval and minval are:'
         write(*,'(2f12.5)') maxval(abs(U)), minval(abs(U))
         write(*,'(a,f12.5)') 'for neutron', maxval(abs(U(1:NUMBPE(0),1:NUMBHE(0),0)))
         write(*,'(a,f12.5)') 'for proton ', maxval(abs(U(1:NUMBPE(1),1:NUMBHE(1),1)))
         print*, 'checking value of V, the abs maxval and minval are:'
         write(*,'(2f12.5)') maxval(abs(V)), minval(abs(V))
         write(*,'(a,f12.5)') 'for neutron', maxval(abs(V(1:NUMBPE(0),1:NUMBHE(0),0)))
         write(*,'(a,f12.5)') 'for proton ', maxval(abs(V(1:NUMBPE(1),1:NUMBHE(1),1)))
         print*, 'checking value of vec, the abs maxval and minval are:'
         write(*,'(2f12.5)') maxval(abs(vec)), minval(abs(vec))
      endif
      ! Feeding DENSHF with vec to get the Gamma1
      call fed_DENSHF

      endsubroutine rho1_to_gamma1

!===============================================================================
      subroutine gamma1_to_rho1
      use hfodd_gs_data, only : NDISOS
      implicit none
      real*8 adj_para
      integer ICHARG, NDIMHE, NDIMPE

      if(IADB_DEBUG .eq. 1) then
         write(*,'(a)') ' checking gamma1, abs maxval and minval'
         write(*,'(2f12.5)') maxval(abs(gamma1)), minval(abs(gamma1))
         write(*,'(a,f12.5)') 'for neutron', maxval(abs(gamma1(1:NUMBPE(0),1:NUMBHE(0),0)))
         write(*,'(a,f12.5)') 'for proton ', maxval(abs(gamma1(1:NUMBPE(1),1:NUMBHE(1),1)))
      endif

      !rho1 = 1.0/eph * (rho_derivative - gamma1)
      rho1 = 0.d0
      DO ICHARG=0,NDISOS
         NDIMHE = NUMBHE(ICHARG)
         NDIMPE = NUMBPE(ICHARG)
         rho1(1:NDIMPE,1:NDIMHE,ICHARG) = 1.d0/eph(1:NDIMPE,1:NDIMHE,ICHARG) * &
        & (rho_derivative(1:NDIMPE,1:NDIMHE,ICHARG) - gamma1(1:NDIMPE,1:NDIMHE,ICHARG))
      ENDDO

      if(IADB_DEBUG .eq. 1 .and. .false.) then
         adj_para = 1.d0
         rho1 = 1.0/eph * (rho_derivative - adj_para*gamma1)
         print*, 'A coefficient is multiplied to gamma1 ...', adj_para
      endif

      if(IADB_DEBUG .eq. 1) then
         print*, '## in gamma1torho1 ##'
         print*, 'gamma1(1,1,0)', gamma1(1,1,0)
         print*, 'gamma1(1,1,1)',gamma1(1,1,1)
         print*, 'rho1(1,1,0)', rho1(1,1,0)
         print*, 'rho1(1,1,1)',rho1(1,1,1)
         print*, 'gamma1(3,4,0)', gamma1(3,4,0)
         print*, 'gamma1(3,4,1)',gamma1(3,4,1)
         print*, 'rho1(3,4,0)', rho1(3,4,0)
         print*, 'rho1(3,4,1)',rho1(3,4,1)
      endif

      endsubroutine gamma1_to_rho1

!===============================================================================
      subroutine calculate_mass
      use hfodd_gs_data, only : NDISOS
      implicit none

      integer ISTATP, ISTATH
      integer ICHARG

      complex resu
      complex resu_half
      real mass_isos(0:NDISOS)

      !if(IADB_DEBUG .eq. 1) print*, 'Warning: Is mass real or imaginary?'
      mass = 0.d0
      ! -- mass = Re[Tr(rho1*rho_derivative)]
      DO ICHARG = 0, NDISOS
         resu = 0.d0
         resu_half = 0.d0
         do ISTATP = 1, NUMBPE(ICHARG)
            do ISTATH = 1, NUMBHE(ICHARG)
               resu = resu + rho1(ISTATP,ISTATH,ICHARG) * conjg(rho_derivative(ISTATP,ISTATH,ICHARG)) &
                    &  + conjg(rho1(ISTATP,ISTATH,ICHARG)) * rho_derivative(ISTATP,ISTATH,ICHARG)
               resu_half = resu_half &
                    &  + rho1(ISTATP,ISTATH,ICHARG) * conjg(rho_derivative(ISTATP,ISTATH,ICHARG))
            enddo
         enddo
         mass_isos(ICHARG) = real(resu)
         mass = mass + mass_isos(ICHARG)
         if(IADB_DEBUG .eq. 1) then
            write(*,'(a)') ' checking rho1 in calculate_mass, abs maxval and minval'
            write(*,'(2f12.5)') maxval(abs(rho1)), minval(abs(rho1))
            write(*,'(a,i5,f12.5)') 'N/P', ICHARG,maxval(abs(rho1(1:NUMBPE(ICHARG),&
            &                                     1:NUMBHE(ICHARG),ICHARG)))
            write(*,'(a)') ' real/imag part of collective mass:'
            write(*,'(i5,2f12.5)') ICHARG, real(resu), imag(resu)
            write(*,'(a)') ' real/imag part of resu_half:'
            write(*,'(i5,2e20.10)') ICHARG, real(resu_half), imag(resu_half)
         endif
      ENDDO

      endsubroutine calculate_mass

!===============================================================================
! calculating matrix elements of rho_derivative in PH space
      subroutine calculate_rho_derivative
      use hfodd_gs_data
      implicit none

      real    FACTOR
      real,   allocatable :: HAUXPP(:,:), HAUXPM(:,:), HAUXMP(:,:), HAUXMM(:,:)
      real,   allocatable :: HAUXTT(:,:)
      integer IBASEI, IBASEJ, ICHARG, ISTATP, ISTATH
      complex RESUPP, RESUPM, RESUMP, RESUMM
      integer NDIMHE

      allocate(HAUXPP(1:NDBASE,1:NDBASE),HAUXPM(1:NDBASE,1:NDBASE))
      allocate(HAUXMP(1:NDBASE,1:NDBASE),HAUXMM(1:NDBASE,1:NDBASE))
      allocate(HAUXTT(1:NDBASE,1:NDBASE))

      ! rho_derivative(ph) = <p|jy|h>

      ! e.g., for orbital parts
      FACTOR = 1.d0
      CALL INT_SY(HAUXTT,FACTOR)
      CALL INT_LY(IPHAPP(0,0,0),HAUXPP,FACTOR)
      CALL INT_LY(IPHAPM(0,0,0),HAUXPM,FACTOR)
      CALL INT_LY(IPHAMP(0,0,0),HAUXMP,FACTOR)
      CALL INT_LY(IPHAMM(0,0,0),HAUXMM,FACTOR)
      HAUXPP = HAUXPP + HAUXTT
      HAUXPM = HAUXPM + HAUXTT
      HAUXMP = HAUXMP + HAUXTT
      HAUXMM = HAUXMM + HAUXTT

      rho_derivative = 0.d0
      DO ICHARG=0,NDISOS
         WAVOCC = SAVOCC(:,:,:,ICHARG)
         NDIMHE = NUMBHE(ICHARG)
         DO ISTATP = 1, NUMBPE(ICHARG)
            DO ISTATH = 1, NUMBHE(ICHARG)
               RESUPP = 0; RESUPM = 0; RESUMP = 0; RESUMM = 0;
               DO IBASEI=1,LDBASE
                  DO IBASEJ=1,LDBASE
                     RESUPP = RESUPP + conjg( WAVOCC(IBASEI,ISTATP+NDIMHE,0) )   &
     &                               * HAUXPP(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,0)
                     RESUPM = RESUPM + conjg( WAVOCC(IBASEI,ISTATP+NDIMHE,0) )   &
     &                               * HAUXPM(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,1)
                     RESUMP = RESUMP + conjg( WAVOCC(IBASEI,ISTATP+NDIMHE,1) )   &
     &                               * HAUXMP(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,0)
                     RESUMM = RESUMM + conjg( WAVOCC(IBASEI,ISTATP+NDIMHE,1) )   &
     &                               * HAUXMM(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,1)
                  ENDDO
               ENDDO
               rho_derivative(ISTATP,ISTATH,ICHARG) = RESUPP + RESUPM + RESUMP + RESUMM
            ENDDO
         ENDDO
      ENDDO ! For ICHARG

      rho_derivative = rho_derivative * (0.d0,1.d0)

      if(IADB_DEBUG .eq. 1) then
         print*, 'checking value of rho_derivative, the abs maxval and minval are:'
         write(*,'(2f12.5)') maxval(abs(rho_derivative)), minval(abs(rho_derivative))
         print*, 'rho_derivative(1,1,0)', rho_derivative(1,1,0)
         print*, 'rho_derivative(1,1,1)', rho_derivative(1,1,1)
         print*, 'rho_derivative(3,4,0)', rho_derivative(3,4,0)
         print*, 'rho_derivative(3,4,1)', rho_derivative(3,4,1)
         print*, 'eph(1,1,0)', eph(1,1,0)
         print*, 'eph(1,1,1)', eph(1,1,1)
         print*, 'eph(3,4,0)', eph(3,4,0)
         print*, 'eph(3,4,1)', eph(3,4,1)
         DO ICHARG=0,NDISOS
            print*, 'ICHARG = ', ICHARG
            call PRINT_MATRIX_Z( 'rho_derivative', &
                & NUMBPE(ICHARG), NUMBHE(ICHARG), rho_derivative(:,:,ICHARG), M )
            call PRINT_MATRIX( 'EPH', &
                & NUMBPE(ICHARG), NUMBHE(ICHARG), eph(:,:,ICHARG), M )
         ENDDO
         print*, 'BIG_PP & BIG_PM before ATDHF calculation:'
         write(*,*) 'MAXVAL of BIG_PP', maxval(abs(BIG_PP))
         write(*,*) 'MAXVAL of BIG_PM', maxval(abs(BIG_PM))
         call PRINT_MATRIX_Z( 'BIG_PP(0)', LDBASE, LDBASE, BIG_PP(1,1,0), NDBASE )
         call PRINT_MATRIX_Z( 'BIG_PP(1)', LDBASE, LDBASE, BIG_PP(1,1,1), NDBASE )
         call PRINT_MATRIX_Z( 'BIG_PM(0)', LDBASE, LDBASE, BIG_PM(1,1,0), NDBASE )
         call PRINT_MATRIX_Z( 'BIG_PM(1)', LDBASE, LDBASE, BIG_PM(1,1,1), NDBASE )
      endif

      deallocate(HAUXPP,HAUXPM,HAUXMP,HAUXMM,HAUXTT)
      endsubroutine calculate_rho_derivative

!===============================================================================
! Calculating the collective inertial under cranking approximation
! Wroking in PH spapce
      subroutine cranking_mass
      implicit none

      if(IADB_DEBUG .eq. 1) print*, '-- Cranking approximation -- '
      !rho1 = rho_derivative/eph
      gamma1 = 0.d0
      call gamma1_to_rho1
      call calculate_mass
      print*, 'Collective mass for rotation along y-axis is:'
      print*, mass
      stop 'ATDHF calculation ends here -- 20221205'

      endsubroutine cranking_mass

!===============================================================================
      subroutine atdhf_output
      implicit none

      ! generating output of ATDHF calculations
      print*, 'Collective mass for rotation along y-axis is:'
      print*, mass

      endsubroutine atdhf_output

!===============================================================================
! Determining the dimension of particle/hole space ...
!
      subroutine phsize(NUMBPX, NUMBHX)
      use hfodd_gs_data
      implicit none

      integer ICHARG, ISTATE
      integer num_total, iccount
      integer NUMBPX, NUMBHX

      DO ICHARG=0,NDISOS
         iccount = 0
         DO ISTATE=1,NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
            if(M_OCVA(ISTATE,ICHARG) .eq. 1) iccount = iccount + 1
            if(IADB_DEBUG .eq. 1) print*,'M_OCVA', ISTATE, ICHARG, M_OCVA(ISTATE,ICHARG)
         ENDDO
         NUMBHE(ICHARG) = iccount
         NUMBPE(ICHARG) = NUMBSP(0,ICHARG) + NUMBSP(1,ICHARG) - NUMBHE(ICHARG)
      ENDDO

      num_total = NUMBSP(0,0) + NUMBSP(1,0) + NUMBSP(0,1) + NUMBSP(1,1)
      NUMBHX = max(NUMBHE(0), NUMBHE(1))
      NUMBPX = max(NUMBPE(0), NUMBPE(1))
      endsubroutine phsize

!===============================================================================
! ##  Feed DENSHF with the eigen vector of RHO1,                              ##
! ##  Then, extract matrix elements of mean field from INTEGH                 ##
!-------------------------------------------------------------------------------
      subroutine fed_DENSHF
      use hfodd_gs_data
      implicit none

      integer ISTATE, ISTATH, ISTATP, ISTATV
      complex RESU0,  RESU1,  RESU2,  RESU3
      real    FAC
      integer IBASE,  IBASEI, IBASEJ
      integer ICHARG
      integer NDIMPE, NDIMHE
      integer ITEMPX
      !###### IADB_DEBUG ######
      complex RESUPP, RESUPM, RESUMP, RESUMM

      complex, allocatable :: RESUXX(:,:)
      complex c_one, czero

      c_one = (1.d0,0.d0)
      czero = (0.d0,0.d0)

! ---
!     The positive and negative eigenvalues of RHO1 are put into the array V_OCCU consecutively
!     =>  |-S_1,+S_1|-S_2,+S_2|...|-S_N,+S_N|
!     In S, only positive parts are stored

      DO ICHARG=0,NDISOS
!         print*, 'replacing is not preformed ..'
!         WAVOCC = SAVOCC(:,:,:,ICHARG)
!         V_OCCU = 0!SV_OCC(:,ICHARG)
!         goto 333
! --- replace V_OCCU by r_i
         V_OCCU = 0.d0
         DO ISTATH = 1, NUMBHE(ICHARG)
            V_OCCU(2*ISTATH - 1) = -S(ISTATH,ICHARG)
            V_OCCU(2*ISTATH    ) = +S(ISTATH,ICHARG)
         ENDDO
         !V_OCCU = 0.d0
! --- replacing WAVOCC
         BAVOCC = 0.d0
         WAVOCC = SAVOCC(:,:,:,ICHARG)

         if(IADB_DEBUG .eq. 1) then
            print*, 'checking value of WAVOCC before replacing:', ICHARG
            write(*,'(2f12.5)') maxval(abs(WAVOCC)), minval(abs(WAVOCC))
         endif
         ! -- $ B_{h' i} = \sum_{\mu} \varphi_{\mu h'} A_{\mu i} $
         DO ISTATH = 1, NUMBHE(ICHARG)
            DO IBASE=1,LDBASE
               RESU0 = 0.d0; RESU1 = 0.d0; RESU2 = 0.d0; RESU3 = 0.d0;
               DO ISTATE=1,NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
                  IF(ISTATE .LE. NUMBHE(ICHARG)) THEN
                     FAC = -1.d0
                     ISTATV = ISTATE + NUMBPE(ICHARG)
                  ELSE
                     FAC = +1.d0
                     ISTATV = ISTATE - NUMBHE(ICHARG)
                  ENDIF
                  RESU0 = RESU0 + WAVOCC(IBASE,ISTATE,0) * VEC(ISTATV,ISTATH,ICHARG) * FAC
                  RESU1 = RESU1 + WAVOCC(IBASE,ISTATE,1) * VEC(ISTATV,ISTATH,ICHARG) * FAC
                  RESU2 = RESU2 + WAVOCC(IBASE,ISTATE,0) * VEC(ISTATV,ISTATH,ICHARG)
                  RESU3 = RESU3 + WAVOCC(IBASE,ISTATE,1) * VEC(ISTATV,ISTATH,ICHARG)
               ENDDO
               BAVOCC(IBASE,2*ISTATH - 1,0) = RESU0
               BAVOCC(IBASE,2*ISTATH - 1,1) = RESU1
               BAVOCC(IBASE,2*ISTATH    ,0) = RESU2
               BAVOCC(IBASE,2*ISTATH    ,1) = RESU3
            ENDDO
         ENDDO
         WAVOCC = BAVOCC

         if(IADB_DEBUG .eq. 1) then
            print*, 'ICHARG=', ICHARG, NUMBSP(0,ICHARG), NUMBSP(1,ICHARG)
            print*, 'checking value of V_OCCU, the abs maxval and minval are:'
            ITEMPX = NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
            write(*,'(2f12.5)') maxval(abs(V_OCCU(1:ITEMPX))), minval(abs(V_OCCU(1:ITEMPX)))
            print*, 'checking value of WAVOCC, the abs maxval and minval are:'
            write(*,'(2f12.5)') maxval(abs(WAVOCC)), minval(abs(WAVOCC))
            write(*,'(a,i5,f12.5)') 'N/P', ICHARG, maxval(abs(WAVOCC(1:LDBASE,&
            &                                     1:NUMBHE(ICHARG)*2,0:1)))
            print*, 'checking value of VEC, the abs maxval and minval are:'
            write(*,'(a,i5,f12.5)') 'N/P', ICHARG, maxval(abs(VEC(:,:,ICHARG)))
         endif

 333     if(maxval(abs(WAVOCC)) .lt. 1.d-6) print*, 'call DENSHF directly ...'

         Call DENSHF_SKFILD_INTEGH(ICHARG,0)

         if(IADB_DEBUG .eq. 1 .and. .true.) then
            DE_RHO = 0.d0;
            DE_LPR = 0.d0;
            DE_TAU = 0.d0;
            DE_SCU = 0.d0;
            DE_DIV = 0.d0;
            DE_DES = 0.d0;
            DE_GRR = 0.d0;
            print*, 'DE_RHO etc. (EVEN) are set to 0'
         endif

         if(IADB_DEBUG .eq. 1 .and. .false.) then
            DE_SPI = 0.d0;
            DE_KIS = 0.d0;
            DE_LPS = 0.d0;
            DE_ROS = 0.d0;
            DE_ROC = 0.d0;
            DE_CUR = 0.d0;
            DE_KIF = 0.d0;
            print*, 'DE_SPI etc. (ODD) are set to 0'
         endif

!        if(IADB_DEBUG .eq. 1) then
!           print*, ' --- after calling DENSHF ---'
!           print*, 'ICHARG = ', ICHARG
!           print*, 'Maxval of DEN_PP:', maxval(abs(DEN_PP(:,:,:,ICHARG)))
!           print*, 'Maxval of DEN_PM:', maxval(abs(DEN_PM(:,:,:,ICHARG)))
!           print*, ' --- --- --- ---'
!           print*, 'NDXHRM, NDYHRM, NDZHRM', NDXHRM, NDYHRM, NDZHRM
!           print*, 'Maxval of DE_RHO:', maxval(abs(DE_RHO(:,:,:,ICHARG)))
!           print*, 'Maxloc of DE_RHO:', maxloc(abs(DE_RHO(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_TAU:', maxval(abs(DE_TAU(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_LPR:', maxval(abs(DE_LPR(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_DIV:', maxval(abs(DE_DIV(:,:,:,ICHARG)))
!           print*, 'Maxval of DE_SPI:', maxval(abs(DE_SPI(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_KIS:', maxval(abs(DE_KIS(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_GRR:', maxval(abs(DE_GRR(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_LPS:', maxval(abs(DE_LPS(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_ROS:', maxval(abs(DE_ROS(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_ROC:', maxval(abs(DE_ROC(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_CUR:', maxval(abs(DE_CUR(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_KIF:', maxval(abs(DE_KIF(:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_SCU:', maxval(abs(DE_SCU(:,:,:,:,:,ICHARG)))
!           print*, 'Maxval of DE_DES:', maxval(abs(DE_DES(:,:,:,:,:,ICHARG)))
!        endif
      ENDDO ! For ICHARG

      Call DENSHF_SKFILD_INTEGH(-1,1)
      !print*, 'VE_CEN is set to 0 manually, to see How it affects BIG_PP'
      !VE_CEN = 0.d0
      if(IADB_DEBUG .eq. 1) then
         print*, ' --- after calling SKFILD ---'
         print*, 'Maxval of VE_MAS:', maxval(abs(VE_MAS(:,:,:,0:1,:)))
         print*, 'Maxval of VE_CEN:', maxval(abs(VE_CEN(:,:,:,0:1,:)))
         print*, 'Maxval of VE_SPT:', maxval(abs(VE_SPT(:,:,:,0:1)))
         print*, 'Maxval of VE_KIS:', maxval(abs(VE_KIS(:,:,:,:,0:1)))
         print*, 'Maxval of VE_SPI:', maxval(abs(VE_SPI(:,:,:,:,0:1)))
         print*, 'Maxval of VE_CUR:', maxval(abs(VE_CUR(:,:,:,:,0:1)))
         print*, 'Maxval of VE_KIT:', maxval(abs(VE_KIT(:,:,:,:,0:1)))
         print*, 'Maxval of VE_SOR:', maxval(abs(VE_SOR(:,:,:,:,:,0:1)))
         print*, 'Maxval of WD_CEN:', maxval(abs(WD_CEN(:,:,:,0:1)))
         print*, 'Maxval of WD_TAU:', maxval(abs(WD_TAU(:,:,:,0:1)))
         print*, 'Maxval of WD_SOR:', maxval(abs(WD_SOR(:,:,:,:,:,0:1)))
      endif

!      stop

! --- extract matrix elements of Gamma1
! -   $ \Gamma_{ph} = sum_{ij}\varphi_{ip}^* \Gamma_{ij} \varphi_{jh} $
      gamma1 = (0.d0,0.d0)
      DO ICHARG = 0, NDISOS
         if(IADB_DEBUG .eq. 1) then
            call PRINT_MAXVAL_Z( '## Active PRINT_MAXVAL_Z', M, N, gamma1, M, -1 )
         endif

         Call DENSHF_SKFILD_INTEGH(ICHARG,2)

         if(IADB_DEBUG .eq. 1) then
            call PRINT_MAXVAL_Z( '## Deactive PRINT_MAXVAL_Z', M, N, gamma1, M, -2 )
         endif

         if(IADB_DEBUG .eq. 1) then
            print*, 'ICHARG=', ICHARG
            print*, 'checking value of BIG_PP, the abs maxval and minval are:'
            write(*,'(2f12.5)') maxval(abs(BIG_PP)), minval(abs(BIG_PP))
            print*, 'checking value of BIG_PM, the abs maxval and minval are:'
            write(*,'(2f12.5)') maxval(abs(BIG_PM)), minval(abs(BIG_PM))
            print*, 'checking value of SAVOCC, the abs maxval and minval are:'
            write(*,'(2f12.5)') maxval(abs(SAVOCC(:,:,:,ICHARG))), minval(abs(SAVOCC(:,:,:,ICHARG)))
         endif
         WAVOCC = SAVOCC(:,:,:,ICHARG)
         NDIMPE = NUMBPE(ICHARG)
         NDIMHE = NUMBHE(ICHARG)
         allocate(RESUXX(NDIMPE,LDBASE))
         !RESUXX = (0.d0,0.d0)
         !gamma1(:,:,ICHARG) = czero
         if(IADB_DEBUG .eq. 1) then
            print*, 'at beginning, RESUXX(1,1):', RESUXX(1,1)
         endif
         ! -- WF(+) * GAMMA(++) * WF(+)
         call ZGEMM('C', 'N', NDIMPE, LDBASE, LDBASE,                          &
     &              c_one, WAVOCC(1,NDIMHE+1,0), NDBASE, BIG_PP(1,1,0), NDBASE,&
     &              czero, RESUXX, NDIMPE)
         call ZGEMM('N', 'N', NDIMPE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,1,0), NDBASE,              &
     &              czero, gamma1(1,1,ICHARG), M)
         if(IADB_DEBUG .eq. 1) then
            print*, '#### special check ####'
            print*, 'WAVOCC(1,NDIMHE+1,0):', WAVOCC(1,NDIMHE+1,0)
            print*, 'maxval WAVOCC(:,NDIMHE+1,:):', maxval(abs(WAVOCC(:,NDIMHE+1,:)))
            print*, 'BIG_PP(1,1,0):', BIG_PP(1,1,0)
            print*, 'checking value of WAVOCC, the abs maxval and minval are:'
            write(*,'(2f12.5)') maxval(abs(WAVOCC)), minval(abs(WAVOCC))
            print*, 'NDIMPE, LDBASE, LDBASE, NDBASE', NDIMPE, LDBASE, LDBASE, NDBASE
            print*, 'RESUXX(1,1):', RESUXX(1,1)
         endif
         if(IADB_DEBUG .eq. 1) then
            print*, 'ICHARG=', ICHARG
            print*, 'checking value of GAMMA(++), the abs maxval and minval are:'
            write(*,'(2f12.5)') maxval(abs(gamma1)), minval(abs(gamma1))
            print*, 'checking value of RESUXX, the abs maxval and minval are:'
            write(*,'(2e12.5)') maxval(abs(RESUXX)), minval(abs(RESUXX))
            write(*,*) NDIMPE, NDIMHE, LDBASE, NDBASE
         endif
         ! -- WF(+) * GAMMA(+-) * WF(-)
         call ZGEMM('C', 'N', NDIMPE, LDBASE, LDBASE,                          &
     &              c_one, WAVOCC(1,NDIMHE+1,0), NDBASE, BIG_PM(1,1,0), NDBASE,&
     &              czero, RESUXX, NDIMPE)
         call ZGEMM('N', 'N', NDIMPE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,1,1), NDBASE,              &
     &              c_one, gamma1(1,1,ICHARG), M)
         ! -- WF(-) * GAMMA(-+) * WF(+)
         call ZGEMM('C', 'N', NDIMPE, LDBASE, LDBASE,                          &
     &              c_one,WAVOCC(1,NDIMHE+1,1), NDBASE, BIG_PM(1,1,1), NDBASE, &
     &              czero, RESUXX, NDIMPE)
         call ZGEMM('N', 'N', NDIMPE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,1,0), NDBASE,              &
     &              c_one, gamma1(1,1,ICHARG), M)
         ! -- WF(-) * GAMMA(--) * WF(-)
         call ZGEMM('C', 'N', NDIMPE, LDBASE, LDBASE,                          &
     &              c_one, WAVOCC(1,NDIMHE+1,1), NDBASE, BIG_PP(1,1,1), NDBASE,&
     &              czero, RESUXX, NDIMPE)
         call ZGEMM('N', 'N', NDIMPE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,1,1), NDBASE,              &
     &              c_one, gamma1(1,1,ICHARG), M)
         deallocate(RESUXX)

         if(IADB_DEBUG .eq. 1) then
            print*, ' ## at the end of fed_DENSHF ##'
            print*, 'ICHARG = ', ICHARG
            print*, 'gamma1(1,1,ICHARG):', gamma1(1,1,ICHARG)
            print*, 'gamma1(3,4,ICHARG):', gamma1(3,4,ICHARG)
            call PRINT_MATRIX_Z( 'gamma1', NUMBPE(ICHARG), NUMBHE(ICHARG), gamma1(:,:,ICHARG), M )
            write(*,*) 'MAXVAL of BIG_PP', maxval(abs(BIG_PP))
            write(*,*) 'MAXVAL of BIG_PM', maxval(abs(BIG_PM))
            call PRINT_MATRIX_Z( 'BIG_PP(0)', LDBASE, LDBASE, BIG_PP(1,1,0), NDBASE )
            call PRINT_MATRIX_Z( 'BIG_PP(1)', LDBASE, LDBASE, BIG_PP(1,1,1), NDBASE )
            call PRINT_MATRIX_Z( 'BIG_PM(0)', LDBASE, LDBASE, BIG_PM(1,1,0), NDBASE )
            call PRINT_MATRIX_Z( 'BIG_PM(1)', LDBASE, LDBASE, BIG_PM(1,1,1), NDBASE )
         endif
      ENDDO

!      stop

      endsubroutine fed_DENSHF

!---------------------------------------------------------------------
      subroutine broyden_mix_rho1(iter)
      use hfodd_gs_data, only : NDISOS
      implicit none

      integer iter
      integer nddim, ndim
      integer nstart, info, n_iter
      real W0_BRO, WN_BRO
      real alphamix, dr2
      integer i, ich
      real, save, allocatable :: vvin(:), vvout(:)

      integer ICHARG, ISTATP, ISTATH

      nstart = 1
      alphamix = 0.5d0
      W0_BRO=0.01D0
      WN_BRO=1.00D0
      n_iter = 7

      ich = 0
      DO ICHARG = 0, NDISOS
         DO ISTATP = 1, NUMBPE(ICHARG)
            DO ISTATH = 1, NUMBHE(ICHARG)
               ich = ich + 2
            ENDDO
         ENDDO
      ENDDO

      nddim = ich
      ndim = nddim
      print*, 'ndim', ndim

      IF( .NOT.ALLOCATED(vvin) .and. iter .eq. 1) then
         allocate(vvin(ndim), vvout(ndim))
         print*, 'vvin and vvout have been allocated'
      endif

      ich = 0
      DO ICHARG = 0, NDISOS
         DO ISTATP = 1, NUMBPE(ICHARG)
            DO ISTATH = 1, NUMBHE(ICHARG)
               ich = ich + 1
               vvin(ich) = real( rho1(ISTATP, ISTATH, ICHARG) )
               ich = ich + 1
               vvin(ich) = aimag( rho1(ISTATP, ISTATH, ICHARG) )
            ENDDO
         ENDDO
      ENDDO

      ! % new vector
      call rho1_to_gamma1
      call gamma1_to_rho1

      ich = 0
      DO ICHARG = 0, NDISOS
         DO ISTATP = 1, NUMBPE(ICHARG)
            DO ISTATH = 1, NUMBHE(ICHARG)
               ich = ich + 1
               vvout(ich) = real( rho1(ISTATP, ISTATH, ICHARG) )
               ich = ich + 1
               vvout(ich) = aimag( rho1(ISTATP, ISTATH, ICHARG) )
            ENDDO
         ENDDO
      ENDDO

      !vvout = vvout - vvin


      print*, 'calling broyden'
      call  broyden(nddim,ndim,vvout,vvin,alphamix,dr2,iter,         &
     &                            n_iter,W0_BRO,WN_BRO,nstart,info)
      print*, 'INFO of Broyden', info

      ich = 0
      DO ICHARG = 0, NDISOS
         DO ISTATP = 1, NUMBPE(ICHARG)
            DO ISTATH = 1, NUMBHE(ICHARG)
               ich = ich + 1
               rho1(ISTATP, ISTATH, ICHARG) = vvin(ich)
               ich = ich + 1
               rho1(ISTATP, ISTATH, ICHARG) = rho1(ISTATP, ISTATH, ICHARG) + vvin(ich)*(0.d0,1.d0)
            ENDDO
         ENDDO
      ENDDO

      print*, ' --  Broyden mixing -- ', iter, dr2
      !if(iter .ge. 3 .and. dr2 .lt. 0.001) stop
      endsubroutine broyden_mix_rho1


END MODULE ADBATC

!===============================================================================
! --------------------------------------------------------------------%%--------
! --    Applying SVD to a complex matrix RHO1 (M by N)
! --    RHO1 (M by N) = U(M by N) * diag[S(N)] * V(N by N); note: M >= N
! --------------------------------------------------------------------%%--------
      subroutine zcallsvd( M, N, in_A, LDA, U, LDU, V, LDV, S, LDS)
      implicit none

      INTEGER    LWMAX
      PARAMETER  ( LWMAX = 16 )
      integer IWORK0(LWMAX)
      complex CWORK0(LWMAX)
      real RWORK0(LWMAX)

      integer M
      integer N
      complex in_A(LDA,LDV)
      integer LDA
      complex U(LDU,LDV)
      integer LDU
      complex V(LDV,LDV)
      integer LDV
      real S(LDS)
      integer LDS
      complex A(LDA,LDV)

      integer NUMRANK
      integer, dimension(:), allocatable:: IWORK
      integer LIWORK
      complex, dimension(:), allocatable:: CWORK
      integer LWORK
      real, dimension(:), allocatable:: RWORK
      integer LRWORK
      integer INFO

      if(LDA .ne. LDU) stop 'Dimension error: LDA /= LDU'
      if(LDV .ne. LDS) stop 'Dimension error: LDV /= LDS'
      if(LDU .lt. LDV) stop 'LDU should larger than LDV'

! Prevent input matrix from being altered
      A = in_A

! Query the optimal workspace.
      LIWORK = -1
      LWORK = -1
      LRWORK = -1
! --------------------------------------------------------------------72--------

      call zgesvdq( 'A',     'P',     'N',     'U',    'V',             &
     &               M,       N,       A,       LDA,                    &
     &               S,       U,       LDU,     V,      LDV,            &
     &               NUMRANK,                                           &
     &               IWORK0,  LIWORK,                                   &
     &               CWORK0,  LWORK,                                    &
     &               RWORK0,  LRWORK,                                   &
     &               INFO )

! --------------------------------------------------------------------72--------

      LIWORK =   IWORK0( 1 )
      LWORK  =   INT( CWORK0( 1 ) )
      LRWORK =   INT( RWORK0( 1 ) )
      allocate( IWORK(LIWORK), CWORK(LWORK), RWORK(LRWORK) )

! Compute SVD.
! --------------------------------------------------------------------72--------

      call zgesvdq( 'A',     'P',     'N',     'U',    'V',             &
     &               M,       N,       A,       LDA,                    &
     &               S,       U,       LDU,     V,      LDV,            &
     &               NUMRANK,                                           &
     &               IWORK,   LIWORK,                                   &
     &               CWORK,   LWORK,                                    &
     &               RWORK,   LRWORK,                                   &
     &               INFO )
! --------------------------------------------------------------------72--------

! -- Clean up --
      deallocate( IWORK, CWORK, RWORK)

      if(INFO /= 0) stop 'Fatal Error in ZCALLSVD'
      endsubroutine

!===============================================================================
      SUBROUTINE PRINT_MATRIX( DESC, M, N, A, LDA )
      CHARACTER*(*)    DESC
      INTEGER          M, N, LDA
      DOUBLE PRECISION A( LDA, * )
!
      INTEGER          I, J
      INTEGER ich_i, ich_j
      ich_i = M / 4
      ich_j = N / 4
!
      WRITE(*,*)
      WRITE(*,*) DESC
      DO I = 1, M, ich_i
         WRITE(*,9998) ( A( I, J ), J = 1, N, ich_j )
      END DO
!
 9998 FORMAT( 11(:,1X,E10.2) )
      RETURN
      END
!===============================================================================
      SUBROUTINE PRINT_MATRIX_Z( DESC, M, N, A, LDA )
      CHARACTER*(*)    DESC
      INTEGER          M, N, LDA
      complex A( LDA, * )
!
      INTEGER          I, J
      INTEGER ich_i, ich_j
      ich_i = M / 4
      ich_j = N / 4
!
      WRITE(*,*)
      WRITE(*,*) DESC
      DO I = 1, M, ich_i
         WRITE(*,9998) ( A( I, J ), J = 1, N, ich_j )
      END DO
!
 9998 FORMAT( 11(:,1X,2E10.2) )
      RETURN
      END

!===============================================================================
! print the maxval of the matrix A
      SUBROUTINE PRINT_MAXVAL_Z( DESC, M, N, A, LDA, ich )

      integer, save :: imode

      CHARACTER*(*)    DESC
      INTEGER          M, N, LDA
      complex A( LDA, N )

      if(ich .eq. -1) imode = 9999
      if(ich .eq. -2) imode = 0

      if(imode .eq. 9999) then
         WRITE(*,*) 'printing the maxval of the matrix:'
         WRITE(*,*) DESC
         write(*,*) maxval(abs(A))
      endif

      end
