!
!===============================================================================
!     Module for Adiabatic Time-Dependent Hartree-Fock-Bogoliubov (ATDHFB)
!     and Adiabatic Time-Dependent Hartree-Fock (ATDHF) calculations
!     ATDHF(B) 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) = 19
      MODSET(17) = 17
      !
      END SUBROUTINE Set_version_hfodd_adiabatic


!-------------------------------------------------------------------------------
! Module containing data from the 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
      ! == ATDHFB
      USE ALLQUZ
      USE SAVQUA
      use MAP_PP
      use MAP_PM
      use PAIDEL
      USE PP_DEN
      USE PD_DEN
      USE PDTDEN

      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                                                            &
     &       /SPIPHA/ PHASPI(1:NDBASE,0:NDREVE,0:NDSPIN,1:NDTWCE)
      COMMON                                                            &
     &       /CFIPRI/ NFIPRI
      ! --- USED FOR TESTING
      COMMON                                                            &
     &       /DENEXP/ EXPAUX(NDXHRM,NDYHRM,NDZHRM,NDTWDD)
      COMMON                                                            &
     &       /SIZSPE/ SPESIZ(1:4*NDBASE)
      COMMON                                                            &
     &       /EIGLIM/ EIGMIN,EIGMAX
      common                                                            &
     &       /pri_flag/ idsi_flag
      ! --- common block added by ATDHF ---
      COMMON                                                            &
     &       /SVOCCU/ SV_OCC(1:4*NDSTAT,NDISOS)
      COMMON                                                            &
     &       /PRSING/ V_SING(1:4*NDSTAT)
      COMMON                                                            &
     &       /SWSING/ I_SING
      COMMON                                                            &
     &       /NPADBC/ NPTOTL(0:NDISOS)

     ! --- ATDHFB ---
      COMMON                                                            &
     &       /ADBMOD/ IADBMD
      COMMON                                                            &
     &       /ADBEFE/ E_FERMI(0:NDISOS)
      COMMON                                                            &
     &       /QPNUMS/ NUMBQP(0:NDREVE,0:NDISOS)
!      COMMON                                                            &
!     &       /QUAPAR/ ENQUAP(1:2*NDSTAT,0:NDREVE),                      &
!     &                V2QUAP(1:2*NDSTAT,0:NDREVE)
!      COMMON                                                            &
!     &       /QUAISO/ EQPISO(1:2*NDSTAT,0:NDREVE,0:NDISOS),             &
!     &                VQPISO(1:2*NDSTAT,0:NDREVE,0:NDISOS)
      COMMON                                                            &
     &       /QPLEVS/ QPENER(1:4*NDSTAT,0:NDISOS),                      &
     &                QPOCCU(1:4*NDSTAT,0:NDISOS),                      &
     &                QPPARR(1:4*NDSTAT,0:NDISOS),                      &
     &                QPEEQU(1:4*NDSTAT,0:NDISOS),                      &
     &                QPDEQU(1:4*NDSTAT,0:NDISOS),                      &
     &                QPCORR(1:4*NDSTAT,0:NDISOS)
      COMMON                                                            &
     &       /LIMOCC/ NOTOCC(1:2*NDSTAT,0:NDREVE),                      &
     &       /OCCISO/ NOCISO(1:2*NDSTAT,0:NDREVE,0:NDISOS)

      COMMON                                                            &
     &       /DENTOC/ DENSIC(NDXHRM,NDYHRM,NDZHRM),                     &
     &                DENCHC(NDXHRM,NDYHRM,NDZHRM)
      COMPLEX                                                           &
     &          DENSIC,DENCHC,DENSIU
     ! -- ADIABATIC SWITCHES --
      COMMON                                                            &
     &       /ADBRUN/ IABNTE,IABCOL,IABICH,IABMOD,IFAMON
      COMMON                                                            &
     &       /ADBPRC/ DABPRS,ETAFAM,XCHIMX
     ! -- FLAGS --
      COMMON /ISWITS/ IGYUSD
      endmodule hfodd_gs_data

!===============================================================================
! Call HFODD subroutines:
!      DENSHF (ICH = 10), SKFILD (ICH = 1), INTEGH (ICH = 2), TO_COUL (ICH = 3)
!      SKPAIR (ICH = 4),  INTEGD (ICH = 5), GOGENE,GOGPAI,REGENE(ICH = 6)
! Parameters of these subroutines should be pushed into:
! DSIPARAI & DSIPARAR & SIPARAC
!
!      XX20 (ICH = 20) : prepare matrix elements of density operator
!
      subroutine HFODD_SUBS(ICHARG, ICH_IN)
      !use hfodd_gs_data, only : I_SING, V_SING, V_OCCU
      use hfodd_gs_data

      ! Integer INPUT parameters of HFODD SUBS
      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,NEWCOU,           &
     &                    ITWCEN,LTWCEN,LDTWCE,ICOTYP,IDOTHC,INDCOE,    &
     &                    IDEFIN,IDEFIP,IGOGPA,IREGPA,ISEPPA,ICOUPA,    &
     &                    ISREGA,ISREGP,N3LORD,I_SLOW,ISGOGA,ISGOGP,    &
     &                                                IN_FIX,IZ_FIX


     ! Character INPUT parameters of HFODD SUBS
      common /DSIPARAC/ NAMEPN
      CHARACTER NAMEPN*8

      ! Real INPUT parameters of HFODD SUBS
      common /DSIPARAR/                                                 &
     &             AROTLT,ALINLT,PRINIT,XOLDEV,XOLDOD,EFER2N,EFER2P,    &
     &                                  SLOWEV,SLOWPA,DELFIN,DELFIP

      COMMON /QBDIMS/ NCONFX(0:1)

      integer ISTATE, IBASE, ISPIN
      COMPLEX EKECOD, EKECOE, EKESCA, EKEVEC, EKEGOG
      integer IBASEI, IBASEJ
      COMPLEX SU
      integer NCONQB
      integer nich, ich_in, ICH

      ICH = ICH_IN
      ! --- For ATDHFB, density should be calculated in two steps.
      ! --- if IPAHFB = 1, positive ICH means the first step, put zeros to density;
      ! ---                negative ICH means the second step
      if(ich .ge. 0) nich = 1
      if(ich .lt. 0) nich = 2
      ICH = abs(ICH)
      IF(IPAHFB.eq.0 .and. NICH.eq.2) stop 'FATAL ERROR: for no pairing case, USE ICH >= 0'

      ! CHECK whether GOGNY is used
      IGYUSD=ISGOGA

      ! 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

      ! Two center : only one center basis is used
      LTWCEN = 1
      INDCOE = 1

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

      ! In SKFILD, no mixure from previous iteration
      XOLDEV = 0.d0
      XOLDOD = 0.d0
      XOLDPA = 0.d0

! --- call DENSHF
      IF (ICH .EQ. 10) THEN
         IF(NICH .eq. 1) 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

! --- call SKFILD
      IF (ICH .EQ. 1) THEN
         CALL ZEFILD(0)
         CALL ZEFILD(1)
         CALL SKFILD(XOLDEV,XOLDOD,IPNMIX,ITWCEN,LTWCEN)
      ENDIF

! --- call INTEGH
      IF (ICH .EQ. 2) THEN
         !CALL SPZERO
         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,                &
     &                                    LTWCEN,ITWCEN,                &
     &                                           INDCOE)
      ENDIF

! --- calculate Time-odd Coulomb field
      IF (ICH .EQ. 3) THEN
         IF(ICHARG.NE.1) STOP 'error in TO_COUL: Coulomb requires ICHARG = 1'
         IF(ICOUDI.NE.2 .OR. ICOUEX.NE.2) return

         EKECOD=0.d0
         EKECOE=0.d0
         EKESCA=0.d0
         EKEVEC=0.d0

         if(IPAHFB.eq.0 .OR. IPAHFB.eq.1) then

            IF (ICOUDI.EQ.2.OR.ICOUEX.EQ.2) THEN
!
                CALL COUENE(ISIMPY,                 &
     &                      ICOTYP,ICOUDI,ICOUEX,   &
     &                                    IDOTHC,   &
     &                      EKECOD,EKESCA,EKEVEC)
!
                IF (ICOUEX.EQ.2) EKECOE=EKESCA+EKEVEC
!
            END IF

         endif

      ENDIF

! --- call SKPAIR, calculate pairing fields
      IF (ICH .EQ. 4) THEN

          CALL SKPAIR(XOLDPA)

      END IF

! --- call INTEGD
      IF (ICH .EQ. 5) THEN

          CALL INTEGD(ITPNMX,ISIMPY,ICHARG,                 &
     &                DELFIN,DELFIP,IDEFIN,IDEFIP,          &
     &                              JETACM,KETA_R,          &
     &                IGOGPA,IREGPA,ISEPPA,ICOUPA,          &
     &                              NEWGOG,NEWCOU)

      END IF

! --- calculate Gogny mean-field & pairing field
! --- call GOGENE, GOGPAI, REGENE
      IF (ICH .EQ. 6) THEN
         !IF(I_GOGA.EQ.0) return
!
         EKEGOG=0.d0

         if(IPAHFB.eq.0 .OR. IPAHFB.eq.1) then

            IF (I_GOGA.GE.1.AND.NEWGOG.EQ.0) THEN
!
                CALL GOGENE(I_SLOW,SLOWEV,ISGOGA,I_GOGA,ISIMPY,EKEGOG)
!
            END IF

            IF ((IGOGPA.EQ.1.OR.IGOGPA.EQ.2).AND.NEWGOG.EQ.0) THEN
!
                CALL GOGPAI(I_SLOW,SLOWPA,ISGOGP,IGOGPA,ISIMPY,IN_FIX,IZ_FIX)
!
            END IF

            IF (I_REGA.GE.1.AND.NEWGOG.EQ.0.OR.                          &
           &    I_GOGA.GE.1.AND.NEWGOG.EQ.1) THEN
!
                CALL REGENE(I_SLOW,SLOWEV,SLOWPA,                        &
           &                I_REGA,IREGPA,ISREGA,ISREGP,                 &
           &                I_GOGA,IGOGPA,ISGOGA,ISGOGP,                 &
           &                N3LORD,ISIMPY,IN_FIX,IZ_FIX,EKEREG,NEWGOG)
!
            END IF
         endif

      ENDIF


!-- XX20
      IF (ICH .EQ. 20) THEN
         ! ---- calculate matrix elements of density operator ----
         if(IPAHFB.eq.0 .and. .true.) then
            DEN_PP(:,:,:,ICHARG) = 0; DEN_PM(:,:,:,ICHARG) = 0;
            DO IBASEI=1,LDBASE
               DO IBASEJ=1,LDBASE
                  DO ISTATE=1,NUMBSP(0,ICHARG)+NUMBSP(1,ICHARG)
                     SU = V_OCCU(ISTATE)**2*V_SING(ISTATE)
                     DEN_PP(IBASEI,IBASEJ,0,ICHARG) = DEN_PP(IBASEI,IBASEJ,0,ICHARG) &
                    &   + WAVOCC(IBASEI,ISTATE,0)*SU*conjg(WAVOCC(IBASEJ,ISTATE,0))
                     DEN_PP(IBASEI,IBASEJ,1,ICHARG) = DEN_PP(IBASEI,IBASEJ,1,ICHARG) &
                    &   + WAVOCC(IBASEI,ISTATE,1)*SU*conjg(WAVOCC(IBASEJ,ISTATE,1))
                     DEN_PM(IBASEI,IBASEJ,0,ICHARG) = DEN_PM(IBASEI,IBASEJ,0,ICHARG) &
                    &   + WAVOCC(IBASEI,ISTATE,0)*SU*conjg(WAVOCC(IBASEJ,ISTATE,1))
                     DEN_PM(IBASEI,IBASEJ,1,ICHARG) = DEN_PM(IBASEI,IBASEJ,1,ICHARG) &
                    &   + WAVOCC(IBASEI,ISTATE,1)*SU*conjg(WAVOCC(IBASEJ,ISTATE,0))
                  ENDDO
               ENDDO
            ENDDO
         endif

         ! ---- calculate matrix elements of density operator and pairing tensor ----
         if(IPAHFB.eq.1 .and. .true.) then
            IF(NICH .eq. 1) THEN
               DEN_PP(:,:,:,ICHARG) = 0; DEN_PM(:,:,:,ICHARG) = 0;
               CHI_PP(:,:,:,ICHARG) = 0; CHI_PM(:,:,:,ICHARG) = 0;
            ENDIF
            NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
            NCONQB = NCONFX(ICHARG)
            DO IBASEI=1,LDBASE
               DO IBASEJ=1,LDBASE
                  DO ISTATE=1,NCONQB
                     SU = V_SING(ISTATE)
                     DEN_PP(IBASEI,IBASEJ,0,ICHARG) = DEN_PP(IBASEI,IBASEJ,0,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,0)*SU*conjg(BWAQUZ(IBASEJ,ISTATE,0))
                     DEN_PP(IBASEI,IBASEJ,1,ICHARG) = DEN_PP(IBASEI,IBASEJ,1,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,1)*SU*conjg(BWAQUZ(IBASEJ,ISTATE,1))
                     DEN_PM(IBASEI,IBASEJ,0,ICHARG) = DEN_PM(IBASEI,IBASEJ,0,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,0)*SU*conjg(BWAQUZ(IBASEJ,ISTATE,1))
                     DEN_PM(IBASEI,IBASEJ,1,ICHARG) = DEN_PM(IBASEI,IBASEJ,1,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,1)*SU*conjg(BWAQUZ(IBASEJ,ISTATE,0))

                     CHI_PP(IBASEI,IBASEJ,0,ICHARG) = CHI_PP(IBASEI,IBASEJ,0,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,0)*SU*conjg(AWAQUZ(IBASEJ,ISTATE,0))
                     CHI_PP(IBASEI,IBASEJ,1,ICHARG) = CHI_PP(IBASEI,IBASEJ,1,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,1)*SU*conjg(AWAQUZ(IBASEJ,ISTATE,1))
                     CHI_PM(IBASEI,IBASEJ,0,ICHARG) = CHI_PM(IBASEI,IBASEJ,0,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,0)*SU*conjg(AWAQUZ(IBASEJ,ISTATE,1))
                     CHI_PM(IBASEI,IBASEJ,1,ICHARG) = CHI_PM(IBASEI,IBASEJ,1,ICHARG) &
                    &   + BWAQUZ(IBASEI,ISTATE,1)*SU*conjg(AWAQUZ(IBASEJ,ISTATE,0))
                  ENDDO
               ENDDO
            ENDDO
         endif

      ENDIF
!-- XX20

      endsubroutine HFODD_SUBS

!===============================================================================
! This subroutine should be called before atdhf_implementation
! It pushes input parameters of HFODD subs used by ATDHF(B) into common block
! DSIPARAI(for integer), DSIPARAR(for real), and DSIPARAC(for character).
!
      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,NEWCOU,           &
     &                    ITWCEN,LTWCEN,LDTWCE,ICOTYP,IDOTHC,INDCOE,    &
     &                    IDEFIN,IDEFIP,IGOGPA,IREGPA,ISEPPA,ICOUPA,    &
     &                    ISREGA,ISREGP,N3LORD,I_SLOW,ISGOGA,ISGOGP,    &
     &                                                IN_FIX,IZ_FIX,    &
     &                                                       NAMEPN,    &
     &             AROTLT,ALINLT,PRINIT,XOLDEV,XOLDOD,EFER2N,EFER2P,    &
     &                                  SLOWEV,SLOWPA,DELFIN,DELFIP     )

      PARAMETER (NUMB_DSI_PARA_INT = 67, NUMB_DSI_PARA_REAL = 11)
      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,NEWCOU,           &
     &                    ITWCEN,LTWCEN,LDTWCE,ICOTYP,IDOTHC,INDCOE,    &
     &                    IDEFIN,IDEFIP,IGOGPA,IREGPA,ISEPPA,ICOUPA,    &
     &                    ISREGA,ISREGP,N3LORD,I_SLOW,ISGOGA,ISGOGP,    &
     &                                                IN_FIX,IZ_FIX    /)

      DSI_PARA_REAL = (/AROTLT,ALINLT,PRINIT,XOLDEV,XOLDOD,EFER2N,EFER2P,&
     &                                       SLOWEV,SLOWPA,DELFIN,DELFIP/)

      DSI_PARA_CHAR = NAMEPN

      endsubroutine PUSH_DSI_PARA

MODULE ADBATC
!===============================================================================
! %%- Data for ATDHF calculations
! ## parameters, provided by the main call?
!     -- max_iteration: maximal iterations
!     -- atd_precision: numerical tolerance for ATDHF(B) calculation
!     -- xchi : mixing parameter for ATDHF(B) iterations
!     -- ecut_adb: energy cutoff
!     -- M, N: dimension of particle/hole configurations
!
! ## data blocks
!     # ATDHF
!     -- elp, elh: particle/hole single particle energy
!     -- 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 \rho1
!     -- U, V: eigenvector of \rho1
!     -- mass: collective mass
!     -- vec : eigenvector of rho1
!     -- rho1_old: rho1 in previous iteration
! ##
      implicit none

      integer  max_iteration
      real     atd_precision
      real     xchi
      real     ecut_adb
      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(:,:,:)
      complex, allocatable :: rho1_old(:,:,:)

! ##
!     # ATDHFB
!     -- F_MATS: matrix element of i*F
!     -- E_MATS: matrix element of EUV
!     -- D_CORR: matrix element of Z
!     -- H_CORR: matrix element of E1
!     -- D_CORR_old: D_CORR in previous iteration
!     -- NCONQX: length of the above arrays
!     -- NCONFX: dimension of q.p. space after truncation
! ##
      complex, allocatable :: F_MATS(:,:,:)
      complex, allocatable :: E_MATS(:,:,:)
      complex, allocatable :: D_CORR(:,:,:)
      complex, allocatable :: H_CORR(:,:,:)
      complex, allocatable :: D_CORR_old(:,:,:)
      integer  NCONQX, NCONFX
      COMMON /QBDIMS/ NCONFX(0:1)

! ##
!     # Array in which HF(B) density are saved
      !complex, allocatable :: DE_RHO_SAV(:,:,:,:,:)
      !complex, allocatable :: PD_RHO_SAV(:,:,:,:,:)
      !complex, allocatable :: PP_RHO_SAV(:,:,:,:,:)
      !complex, allocatable :: DE_SPI_SAV(:,:,:,:,:,:)

      complex, allocatable :: DE_RHO_SAV(:,:,:,:,:)
      complex, allocatable :: DE_TAU_SAV(:,:,:,:,:)
      complex, allocatable :: DE_LPR_SAV(:,:,:,:,:)
      complex, allocatable :: DE_DIV_SAV(:,:,:,:,:)
      !complex, allocatable :: DE_DIJ_SAV(:,:,:,:,:)
      complex, allocatable :: PD_RHO_SAV(:,:,:,:,:)
      complex, allocatable :: PD_TAU_SAV(:,:,:,:,:)
      complex, allocatable :: PD_LPR_SAV(:,:,:,:,:)
      complex, allocatable :: PP_RHO_SAV(:,:,:,:,:)
      complex, allocatable :: PP_TAU_SAV(:,:,:,:,:)
      complex, allocatable :: PP_LPR_SAV(:,:,:,:,:)
      complex, allocatable :: DE_SCU_SAV(:,:,:,:,:,:,:)
      complex, allocatable :: DE_DES_SAV(:,:,:,:,:,:,:)
      complex, allocatable :: PD_SCU_SAV(:,:,:,:,:,:,:)
      complex, allocatable :: PP_SCU_SAV(:,:,:,:,:,:,:)

! --- Auxiliary arrays ...
!     -- NUMBHE: number of hole state
!     -- NUMBPE: number of particle state
!     -- NUMBCS: number of total states under cutoff 'ecut_adb'
      integer  NUMBCS(0:1)
      integer  NUMBHE(0:1), NUMBPE(0:1)
      complex, allocatable :: BAVOCC(:,:,:)

! --- Debug mode flags
      integer  IADB_DEBUG
      real adj_rho, adj_gamma
      integer ifammode
      integer i_svd_on, imat_hh
      integer iset_num

!---- FAM coefficient ----
      real fam_eta

!---- Lables for the collective variables
      character*2 axis_str(0:2)

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

      integer i
      real mass_old, mass_res
      ! == IADB_DEBUG
      integer ICHARG

! --- Initialization
      call read_input_para
      call atdhfb_data_init(0)

! --- Cranking approximation
      call IB_mass

! --- Main iteration
      mass_old = -9999

      H_CORR = 0.d0
      call E1_to_Z
      do i = 1, max_iteration

          if(i .gt. 1) then
             D_CORR = xchi*D_CORR + (1-xchi)*D_CORR_old
          endif
          D_CORR_old = D_CORR
          call Z_to_E1
          call E1_to_Z

          call ATDHFB_mass
          mass_res = mass_old - mass
          if(abs(mass_old - mass) < atd_precision .and. i .gt. 15) exit
          if(IADB_DEBUG .eq. 1 .or. IADB_DEBUG .eq. -1) then
             write(*,'(a,i4,a,f12.5,a,f12.5,23X,1H*)') &
             &'# iter = ', i, ' mass = ', mass, ' residual ',mass_old - mass
          endif
          mass_old = mass
      enddo
      if(i > max_iteration) then
         write(*,'(a,i4,a,f12.5,11X,1H*)') &
         &'* Calculation terminated after ', i-1, ' steps, residual is ', mass_res
      endif

! --- Generating output & clean up
      call atdhfb_output
      !call atdhfb_data_init(1)

      endsubroutine atdhfb_implementation

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

      integer ich
      integer INDXMU, INDXNU
      integer NUMBUX, NUMBVX
      integer ICHARG
      integer NCONQB,NDIMUE,NDIMVE
      integer NCONQN,NCONQP

! ---  Set the mass term to 0
      if(ifammode .eq. 0) HBMRPA = 0.d0

! --- Determine the dimension of arrays according to Quasiparticle configuration
      NCONQN = NUMBQP(0,0) + NUMBQP(1,0)
      NCONQP = NUMBQP(0,1) + NUMBQP(1,1)
      NCONQX = max(NCONQN,NCONQP)
      !NCONQX = 2*LDBASE
      !NCONFX = NCONQX
      NCONFX(0) = NCONQN
      NCONFX(1) = NCONQP

! --- Applying Quasiparticle cutoff ---
! ### Because of the structure of QPENER, an error will happen when
! ### ECUT_ABD < CUTOFF
   if(.false.) then
      DO ICHARG=0,NDISOS
         INDXNU = 0
         do INDXMU = 1, 2*LDBASE
            if(QPENER(INDXMU,ICHARG) .ge. ecut_adb-E_FERMI(ICHARG)) EXIT
            INDXNU = INDXNU + 1
         enddo
         NCONFX(ICHARG) = INDXNU
      ENDDO
      NCONQX = max(NCONFX(0),NCONFX(1))
   endif

      if(ich == 0) then
         ! --- allocations
         allocate(E_MATS(NCONQX,NCONQX,0:NDISOS))
         allocate(D_CORR(NCONQX,NCONQX,0:NDISOS), H_CORR(NCONQX,NCONQX,0:NDISOS))
         allocate(F_MATS(NCONQX,NCONQX,0:NDISOS))
         allocate(D_CORR_old(NCONQX,NCONQX,0:NDISOS))
         call RS_even_density(0)

         ! --- reset skfild, SKPAIR, remove the effect of linear mixing
         Call HFODD_SUBS(-1,1)
         Call HFODD_SUBS(-1,4)
         ! --- save time-even/odd densities ---
         call RS_even_density(1)

         ! --- Prepare the quasiparticle energy matrix
         E_MATS = 0.d0
         DO ICHARG=0,NDISOS
            NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
            NCONQB = NCONFX(ICHARG)
            do INDXMU = 1, NCONQB
               do INDXNU = 1, NCONQB
                  E_MATS(INDXMU,INDXNU,ICHARG) = QPENER(INDXMU,ICHARG) + QPENER(INDXNU,ICHARG)
               enddo
            enddo
         ENDDO

         ! --- calculate the derivative of density
         call calculate_F_MATS

      endif

      if(ich == 1) then
         ! --- free memory space...
         deallocate(D_CORR, H_CORR, F_MATS, E_MATS, D_CORR_old)
      endif

      endsubroutine atdhfb_data_init

!===============================================================================
! Calculating the collective inertial under cranking approximation
! Wroking in QP spapce
      subroutine IB_mass
      use hfodd_gs_data
      implicit none

      !Z = iF/EUV
      H_CORR = 0.d0
      call E1_to_Z
      call ATDHFB_mass

          WRITE(NFIPRI,'(79(1H*),/,                1H*,77X,1H*,/,  &
     &  1H*,2X,''Inglis-Belyaev MOMENT OF INERTIA ALONG '',        &
     &                                               A2,''AXIS '', &
     &                         ''IS '',F12.5,'' MeV-1'', 8X,1H*,/, &
     &                                             1H*,77X,1H*)')  &
     &    axis_str(iset_num), mass

      endsubroutine IB_mass

!===============================================================================
! update Z with E1
      subroutine E1_to_Z
      use hfodd_gs_data, only : NDISOS, NUMBQP
      implicit none
      integer ICHARG,NCONQB

      ! Z = 1.0/EUV * (iF - E1)
      D_CORR = 0.d0
      DO ICHARG=0,NDISOS
         NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
         NCONQB = NCONFX(ICHARG)
         D_CORR(1:NCONQB,1:NCONQB,ICHARG) = 1.d0/E_MATS(1:NCONQB,1:NCONQB,ICHARG) * &
        & (F_MATS(1:NCONQB,1:NCONQB,ICHARG) - H_CORR(1:NCONQB,1:NCONQB,ICHARG))
      ENDDO

      endsubroutine E1_to_Z

!===============================================================================
      subroutine ATDHFB_mass
      use hfodd_gs_data, only : NDISOS, NUMBQP
      implicit none

      integer INDXMU,INDXNU
      integer ICHARG,NCONQB

      complex resu
      complex resu_half
      real mass_isos(0:NDISOS)

      mass = 0.d0
      ! -- mass = 1/2 i * Re[Tr(F^* Z - F Z^*)]
      !         = - 1/2 Re[Tr(F_MATS^* Z + F_MATS Z^*)], note F_MATS = i*F
      DO ICHARG = 0, NDISOS
         NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
         NCONQB = NCONFX(ICHARG)
         resu = 0.d0
         resu_half = 0.d0
         do INDXMU = 1, NCONQB
            do INDXNU = 1, NCONQB
               resu = resu + D_CORR(INDXNU,INDXMU,ICHARG) * conjg(F_MATS(INDXMU,INDXNU,ICHARG)) &
                    & + conjg(D_CORR(INDXNU,INDXMU,ICHARG)) * F_MATS(INDXMU,INDXNU,ICHARG)
               resu_half = resu_half &
                    & + D_CORR(INDXNU,INDXMU,ICHARG) * conjg(F_MATS(INDXMU,INDXNU,ICHARG))
            enddo
         enddo
         mass_isos(ICHARG) = -real(resu)*0.5d0
         mass = mass + mass_isos(ICHARG)

      ENDDO

      endsubroutine ATDHFB_mass

!===============================================================================
subroutine calculate_F_MATS
! calculating the matrix element of iF
! TDRRHO : [j,\rho_0],   note TDRRHO = i \dot{\rho}
! TDRKAP : (j \kappa_0 + \kappa j^*), note TDRKAP = i \dot{\kappa}
use hfodd_gs_data
implicit none

real    FACTOR
real,   allocatable :: HAUXPP(:,:), HAUXPM(:,:), HAUXMP(:,:), HAUXMM(:,:)
complex,allocatable :: FAUXPP(:,:), FAUXPM(:,:), FAUXMP(:,:), FAUXMM(:,:)
real,   allocatable :: HAUXTT(:,:)
complex,allocatable :: TDRRHO(:,:,:,:,:),TDRKAP(:,:,:,:,:)
integer IBASEI, IBASEJ, IBASEM, ICHARG, INDXMU, INDXNU, NCONQB
complex RESUPP, RESUPM, RESUMP, RESUMM
integer ISTATE,IREVEI,IREVEJ,JREVEI,JREVEJ
complex c_one, czero, cmone, i_unit
complex,allocatable :: RESUXX(:,:)
complex anti_sym

c_one = (1.d0,0.d0)
czero = (0.d0,0.d0)
cmone = -c_one
i_unit = (0.d0,1.d0)

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))
allocate(FAUXPP(1:NDBASE,1:NDBASE),FAUXPM(1:NDBASE,1:NDBASE))
allocate(FAUXMP(1:NDBASE,1:NDBASE),FAUXMM(1:NDBASE,1:NDBASE))

allocate(TDRRHO(1:NDBASE,1:NDBASE,0:1,0:1,0:NDISOS),TDRKAP(1:NDBASE,1:NDBASE,0:1,0:1,0:NDISOS))

! iset_num: 0 -- JY; 1 -- JZ; 2 -- JX;
axis_str(0) = 'Y-'; axis_str(1) = 'Z-'; axis_str(2) = 'X-'


!rho = B^* B^T; kappa = B^* A^T
DEN_PP = 0; DEN_PM = 0; CHI_PP = 0; CHI_PM = 0;
DO ICHARG=0,NDISOS
   NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
   !NCONQB = NCONFX(ICHARG) !?
   AWAQUZ(:,:,:)=ASVQUA(:,:,:,ICHARG)
   BWAQUZ(:,:,:)=BSVQUA(:,:,:,ICHARG)
   DO IBASEI=1,LDBASE
      DO IBASEJ=1,LDBASE
         DO ISTATE=1,NCONQB
            DEN_PP(IBASEI,IBASEJ,0,ICHARG) = DEN_PP(IBASEI,IBASEJ,0,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,0)*conjg(BWAQUZ(IBASEJ,ISTATE,0))
            DEN_PP(IBASEI,IBASEJ,1,ICHARG) = DEN_PP(IBASEI,IBASEJ,1,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,1)*conjg(BWAQUZ(IBASEJ,ISTATE,1))
            DEN_PM(IBASEI,IBASEJ,0,ICHARG) = DEN_PM(IBASEI,IBASEJ,0,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,0)*conjg(BWAQUZ(IBASEJ,ISTATE,1))
            DEN_PM(IBASEI,IBASEJ,1,ICHARG) = DEN_PM(IBASEI,IBASEJ,1,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,1)*conjg(BWAQUZ(IBASEJ,ISTATE,0))

            CHI_PP(IBASEI,IBASEJ,0,ICHARG) = CHI_PP(IBASEI,IBASEJ,0,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,0)*conjg(AWAQUZ(IBASEJ,ISTATE,0))
            CHI_PP(IBASEI,IBASEJ,1,ICHARG) = CHI_PP(IBASEI,IBASEJ,1,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,1)*conjg(AWAQUZ(IBASEJ,ISTATE,1))
            CHI_PM(IBASEI,IBASEJ,0,ICHARG) = CHI_PM(IBASEI,IBASEJ,0,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,0)*conjg(AWAQUZ(IBASEJ,ISTATE,1))
            CHI_PM(IBASEI,IBASEJ,1,ICHARG) = CHI_PM(IBASEI,IBASEJ,1,ICHARG) &
           &   + BWAQUZ(IBASEI,ISTATE,1)*conjg(AWAQUZ(IBASEJ,ISTATE,0))
         ENDDO
      ENDDO
   ENDDO

   ! -- anti-symmetrize kappa
   DO IBASEI=1,LDBASE
      DO IBASEJ=1,LDBASE
         anti_sym = (CHI_PM(IBASEI,IBASEJ,0,ICHARG) - CHI_PM(IBASEJ,IBASEI,1,ICHARG))/2
         CHI_PM(IBASEI,IBASEJ,0,ICHARG) = + anti_sym
         CHI_PM(IBASEJ,IBASEI,1,ICHARG) = - anti_sym
      ENDDO
   ENDDO

ENDDO ! ICHARG


FACTOR = 1.d0
if(iset_num == 0) then
   ! -- set 0 --
   if(IADB_DEBUG .eq. 1) print*, 'Matrix element of JY is used'
   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)

   FAUXPP = HAUXPP*i_unit + HAUXTT
   FAUXPM = HAUXPM*i_unit
   FAUXMP = HAUXMP*i_unit
   FAUXMM = HAUXMM*i_unit - HAUXTT

else if(iset_num == 1) then
   ! -- set 1 --
   if(IADB_DEBUG .eq. 1) print*, 'Matrix element of JZ is used'
   CALL INT_SZ(HAUXTT,FACTOR)
   CALL INT_LZ(IPHAPP(0,0,0),HAUXPP,FACTOR)
   CALL INT_LZ(IPHAPM(0,0,0),HAUXPM,FACTOR)
   CALL INT_LZ(IPHAMP(0,0,0),HAUXMP,FACTOR)
   CALL INT_LZ(IPHAMM(0,0,0),HAUXMM,FACTOR)

   FAUXPP = HAUXPP*i_unit
   FAUXPM = HAUXPM*i_unit + HAUXTT*i_unit
   FAUXMP = HAUXMP*i_unit - HAUXTT*i_unit
   FAUXMM = HAUXMM*i_unit

else if(iset_num == 2) then
   ! -- set 2 --
   if(IADB_DEBUG .eq. 1) print*, 'Matrix element of JX is used'
   CALL INT_SX(HAUXTT,FACTOR)
   CALL INT_LX(IPHAPP(0,0,0),HAUXPP,FACTOR)
   CALL INT_LX(IPHAPM(0,0,0),HAUXPM,FACTOR)
   CALL INT_LX(IPHAMP(0,0,0),HAUXMP,FACTOR)
   CALL INT_LX(IPHAMM(0,0,0),HAUXMM,FACTOR)

   FAUXPP = HAUXPP*i_unit
   FAUXPM = HAUXPM*i_unit + HAUXTT
   FAUXMP = HAUXMP*i_unit + HAUXTT
   FAUXMM = HAUXMM*i_unit

endif ! iset_num


!----------------------------------------------------------------------
! TDRKAP : (j \kappa_0 + \kappa j^*) == <J,K>
!--- <J,K>++ = J++K++ + J+-K-+ + K++J++ + K+-J+-
!--- <J,K>+- = J++K+- + J+-K-- + K++J-+ + K+-J--
!--- <J,K>-+ = J-+K++ + J--K-+ + K-+J++ + K--J+-
!--- <J,K>-- = J-+K+- + J--K-- + K-+J-+ + K--J--
DO ICHARG=0,NDISOS
   ! j^y_{at}\kappa^0_{tb}+\kappa^0_{at}j^{y}_{bt}
   DO IBASEI=1,LDBASE
      DO IBASEJ=1,LDBASE
         RESUPP = 0; RESUPM = 0; RESUMP = 0; RESUMM = 0;
         DO IBASEM=1,LDBASE
            RESUPP = RESUPP                                               &
  &                + FAUXPP(IBASEI,IBASEM)*CHI_PP(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXPM(IBASEI,IBASEM)*CHI_PM(IBASEM,IBASEJ,1,ICHARG) &
  &                + CHI_PP(IBASEI,IBASEM,0,ICHARG)*FAUXPP(IBASEJ,IBASEM) &
  &                + CHI_PM(IBASEI,IBASEM,0,ICHARG)*FAUXPM(IBASEJ,IBASEM)
            RESUPM = RESUPM                                               &
  &                + FAUXPP(IBASEI,IBASEM)*CHI_PM(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXPM(IBASEI,IBASEM)*CHI_PP(IBASEM,IBASEJ,1,ICHARG) &
  &                + CHI_PP(IBASEI,IBASEM,0,ICHARG)*FAUXMP(IBASEJ,IBASEM) &
  &                + CHI_PM(IBASEI,IBASEM,0,ICHARG)*FAUXMM(IBASEJ,IBASEM)
            RESUMP = RESUMP                                               &
  &                + FAUXMP(IBASEI,IBASEM)*CHI_PP(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXMM(IBASEI,IBASEM)*CHI_PM(IBASEM,IBASEJ,1,ICHARG) &
  &                + CHI_PM(IBASEI,IBASEM,1,ICHARG)*FAUXPP(IBASEJ,IBASEM) &
  &                + CHI_PP(IBASEI,IBASEM,1,ICHARG)*FAUXPM(IBASEJ,IBASEM)
            RESUMM = RESUMM                                               &
  &                + FAUXMP(IBASEI,IBASEM)*CHI_PM(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXMM(IBASEI,IBASEM)*CHI_PP(IBASEM,IBASEJ,1,ICHARG) &
  &                + CHI_PM(IBASEI,IBASEM,1,ICHARG)*FAUXMP(IBASEJ,IBASEM) &
  &                + CHI_PP(IBASEI,IBASEM,1,ICHARG)*FAUXMM(IBASEJ,IBASEM)
         ENDDO
         TDRKAP(IBASEI,IBASEJ,0,0,ICHARG) = RESUPP
         TDRKAP(IBASEI,IBASEJ,0,1,ICHARG) = RESUPM
         TDRKAP(IBASEI,IBASEJ,1,0,ICHARG) = RESUMP
         TDRKAP(IBASEI,IBASEJ,1,1,ICHARG) = RESUMM
      ENDDO
   ENDDO

!----------------------------------------------------------------------
! TDRRHO : [j,\rho_0]
!--- [J,R]++ = J++R++ + J+-R-+ - R++J++ - R+-J-+
!--- [J,R]+- = J++R+- + J+-R-- - R++J+- - R+-J--
!--- [J,R]-+ = J-+R++ + J--R-+ - R-+J++ - R--J-+
!--- [J,R]-- = J-+R+- + J--R-- - R-+J+- - R--J--
   ! j^y_{at}\rho^0_{tb}-\rho^0_{at}j^y_{tb}: anti-Hermitian
   DO IBASEI=1,LDBASE
      DO IBASEJ=1,LDBASE
         RESUPP = 0; RESUPM = 0; RESUMP = 0; RESUMM = 0;
         DO IBASEM=1,LDBASE
            RESUPP = RESUPP                                               &
  &                + FAUXPP(IBASEI,IBASEM)*DEN_PP(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXPM(IBASEI,IBASEM)*DEN_PM(IBASEM,IBASEJ,1,ICHARG) &
  &                - DEN_PP(IBASEI,IBASEM,0,ICHARG)*FAUXPP(IBASEM,IBASEJ) &
  &                - DEN_PM(IBASEI,IBASEM,0,ICHARG)*FAUXMP(IBASEM,IBASEJ)
            RESUPM = RESUPM                                               &
  &                + FAUXPP(IBASEI,IBASEM)*DEN_PM(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXPM(IBASEI,IBASEM)*DEN_PP(IBASEM,IBASEJ,1,ICHARG) &
  &                - DEN_PP(IBASEI,IBASEM,0,ICHARG)*FAUXPM(IBASEM,IBASEJ) &
  &                - DEN_PM(IBASEI,IBASEM,0,ICHARG)*FAUXMM(IBASEM,IBASEJ)
            RESUMP = RESUMP                                               &
  &                + FAUXMP(IBASEI,IBASEM)*DEN_PP(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXMM(IBASEI,IBASEM)*DEN_PM(IBASEM,IBASEJ,1,ICHARG) &
  &                - DEN_PM(IBASEI,IBASEM,1,ICHARG)*FAUXPP(IBASEM,IBASEJ) &
  &                - DEN_PP(IBASEI,IBASEM,1,ICHARG)*FAUXMP(IBASEM,IBASEJ)
            RESUMM = RESUMM                                               &
  &                + FAUXMP(IBASEI,IBASEM)*DEN_PM(IBASEM,IBASEJ,0,ICHARG) &
  &                + FAUXMM(IBASEI,IBASEM)*DEN_PP(IBASEM,IBASEJ,1,ICHARG) &
  &                - DEN_PM(IBASEI,IBASEM,1,ICHARG)*FAUXPM(IBASEM,IBASEJ) &
  &                - DEN_PP(IBASEI,IBASEM,1,ICHARG)*FAUXMM(IBASEM,IBASEJ)
         ENDDO
         TDRRHO(IBASEI,IBASEJ,0,0,ICHARG) = RESUPP
         TDRRHO(IBASEI,IBASEJ,0,1,ICHARG) = RESUPM
         TDRRHO(IBASEI,IBASEJ,1,0,ICHARG) = RESUMP
         TDRRHO(IBASEI,IBASEJ,1,1,ICHARG) = RESUMM
      ENDDO
   ENDDO


ENDDO ! ICHARG


!print*, 'Calculating F_MATS'

!------------------------------------------------------------------------
! F = A^{\dagger}\dot{\rho_0}B^* + A^{\dagger}\dot{\kappa_0}A^*
!   - B^{\dagger}\dot{\kappa}^*_0B^* - B^{\dagger}\dot{\rho}^*_0A^*.
! note: AWAQUZ is A^*; BWAQUZ is B^*
! TDRRHO : Time-derivative of density operator
! TDRKAP : Time-derivative of pairing tensor
!------------------------------------------------------------------------

F_MATS = 0.d0
DO ICHARG=0,NDISOS
   NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
   NCONQB = NCONFX(ICHARG)
   AWAQUZ(:,:,:)=ASVQUA(:,:,:,ICHARG)
   BWAQUZ(:,:,:)=BSVQUA(:,:,:,ICHARG)

   allocate(RESUXX(NCONQB,LDBASE))
   do IREVEI = 0,1
      do IREVEJ = 0,1
         ! \kappa^\dagger = -\kappa^*; \rho^T = \rho^*
         ! --- A^{\dagger}\dot{\rho_0}B^* ---
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                                            &
     &              c_one, AWAQUZ(1,1,IREVEI), NDBASE, TDRRHO(1,1,IREVEI,IREVEJ,ICHARG), NDBASE, &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                                            &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,IREVEJ), NDBASE,                           &
     &              c_one, F_MATS(1,1,ICHARG), NCONQX)

        ! --- A^{\dagger}\dot{\kappa_0}A^*
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                                            &
     &              c_one, AWAQUZ(1,1,IREVEI), NDBASE, TDRKAP(1,1,IREVEI,IREVEJ,ICHARG), NDBASE, &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                                            &
     &              c_one, RESUXX, NCONQB, AWAQUZ(1,1,IREVEJ), NDBASE,                           &
     &              c_one, F_MATS(1,1,ICHARG), NCONQX)

        ! --- -B^{\dagger}\dot{\kappa}^*_0B^* = B^{\dagger}\dot{\kappa}^{\dagger}_0B^*
        ! --- Note TDRKAP = i\dot{\kappa},
        ! --- therefore, TDRKAP^{\dagger} = -i\dot{\kappa}^{\dagger}.
        ! --- Yes, an extra minus appears (CMONE in the following block)...
         if(mod(IREVEI+IREVEJ,2).eq.0) then
            JREVEI = IREVEI; JREVEJ = IREVEJ;
         else
            JREVEI = IREVEJ; JREVEJ = IREVEI;
         endif
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                                            &
     &              c_one, BWAQUZ(1,1,IREVEI), NDBASE, TDRKAP(1,1,JREVEI,JREVEJ,ICHARG), NDBASE, &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                                            &
     &              cmone, RESUXX, NCONQB, BWAQUZ(1,1,IREVEJ), NDBASE,                           &
     &              c_one, F_MATS(1,1,ICHARG), NCONQX)

        ! --- -B^{\dagger}\dot{\rho}^*_0A^* = -B^{\dagger}\dot{\rho}^T_0A^*
         if(mod(IREVEI+IREVEJ,2).eq.0) then
            JREVEI = IREVEI; JREVEJ = IREVEJ;
         else
            JREVEI = IREVEJ; JREVEJ = IREVEI;
         endif
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                                            &
     &              c_one, BWAQUZ(1,1,IREVEI), NDBASE, TDRRHO(1,1,JREVEI,JREVEJ,ICHARG), NDBASE, &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                                            &
     &              cmone, RESUXX, NCONQB, AWAQUZ(1,1,IREVEJ), NDBASE,                           &
     &              c_one, F_MATS(1,1,ICHARG), NCONQX)
      enddo
   enddo
   deallocate(RESUXX)

ENDDO ! ICHARG


deallocate(HAUXPP,HAUXPM,HAUXMP,HAUXMM,HAUXTT)
deallocate(FAUXPP,FAUXPM,FAUXMP,FAUXMM)
deallocate(TDRRHO,TDRKAP)

endsubroutine calculate_F_MATS


!===============================================================================
subroutine Z_to_E1
use hfodd_gs_data
implicit none

integer ICHARG, ISTEP
integer IBASE, ISTATE, JSTATE, LSTATE

COMPLEX, ALLOCATABLE :: CELMTS(:)
REAL, ALLOCATABLE :: EIGVAL(:)
COMPLEX, ALLOCATABLE :: EWAVEF(:,:)
COMPLEX, ALLOCATABLE :: UWAQUA(:,:,:), VWAQUA(:,:,:)

integer LDACTU, NDACTU, NUMVEC
integer NCONQB
integer INDXMU, INDXNU, IKET, IBRA, NCOUNT
integer IBASEI, IREVEI, IREVER
COMPLEX RESU
!
complex c_one, czero, cmone
complex,allocatable :: RESUXX(:,:)

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

DO ICHARG = 0, NDISOS

   !print*, 'ICHARG = ', ICHARG

   NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
   NCONQB = NCONFX(ICHARG)
   NDACTU = 2*NCONQB
   LDACTU = NDACTU
   NUMVEC = NDACTU
   ALLOCATE( CELMTS(1:((NDACTU+1)*NDACTU)/2) )
   ALLOCATE( EIGVAL(1:NDACTU) )
   ALLOCATE( EWAVEF(1:NDACTU,1:NDACTU) )

   ALLOCATE(UWAQUA(NDBASE,NDACTU,0:NDREVE),VWAQUA(NDBASE,NDACTU,0:NDREVE))

   CELMTS = 0.d0

   select case (ifammode)
   case(0)
   ! --- linearalized
      DO INDXMU = 1,NCONQB
         DO INDXNU = 1,NCONQB
            IBRA = INDXMU + NCONQB
            IKET = INDXNU
            NCOUNT=IBRA+((2*NDACTU-IKET)*(IKET-1))/2
            CELMTS(NCOUNT)=-conjg( D_CORR(INDXMU,INDXNU,ICHARG) )
         END DO
      END DO

   case(-1)
   ! --- HFB
      DO INDXMU = 1,NCONQB
            IBRA = INDXMU + NCONQB
            IKET = IBRA
            NCOUNT=IBRA+((2*NDACTU-IKET)*(IKET-1))/2
            CELMTS(NCOUNT)= 1.0d0
      END DO

   case(1)
   ! --- FAM
      DO INDXMU = 1,NCONQB
            IBRA = INDXMU + NCONQB
            IKET = IBRA
            NCOUNT=IBRA+((2*NDACTU-IKET)*(IKET-1))/2
            CELMTS(NCOUNT)= 1.0d0
      END DO
      DO INDXMU = 1,NCONQB
         DO INDXNU = 1,NCONQB
            IBRA = INDXMU + NCONQB
            IKET = INDXNU
            NCOUNT=IBRA+((2*NDACTU-IKET)*(IKET-1))/2
            CELMTS(NCOUNT)=-fam_eta*conjg( D_CORR(INDXMU,INDXNU,ICHARG) )
         END DO
      END DO

   case default
      stop 'unknown runmode ..'
   end select


   CALL DIAMAT(CELMTS,EIGVAL,EWAVEF,LDACTU,NDACTU,NUMVEC)


! -- $ B_{h' i} = \sum_{\mu} \varphi_{\mu h'} A_{\mu i} $
! |A B^*|*|U V^*|
! |B A^*|*|V U^*|

   AWAQUZ(:,:,:)=ASVQUA(:,:,:,ICHARG)
   BWAQUZ(:,:,:)=BSVQUA(:,:,:,ICHARG)
! note AWAQUZ -> A^*, BWAQUZ -> B^*
! U = EWAVEF(1:NCONQB,ISTATE), V = EWAVEF(NCONQB+1:2*NCONQB,ISTATE)
! |\sum_{a} A_{ia}*U_{ab} + B^*_{ia}*V_{ab},    \sum_{a} A_{ia}*V^*_{ab} + B^*_{ia}*U^*_{ab}|
! |\sum_{a} B_{ia}*U_{ab} + A^*_{ia}*V_{ab},    \sum_{a} B_{ia}*V^*_{ab} + A^*_{ia}*U^*_{ab}|
   DO ISTATE = 1, NCONQB
      DO IBASEI = 1, LDBASE
         DO IREVEI = 0, 1
            ! \sum_{a} A_{ia}*U_{ab} + B^*_{ia}*V_{ab}
            RESU = 0.d0
            DO INDXMU = 1, NCONQB
               RESU = RESU + conjg(AWAQUZ(IBASEI,INDXMU,IREVEI))*EWAVEF(INDXMU,ISTATE) &
                    &      + BWAQUZ(IBASEI,INDXMU,IREVEI)*EWAVEF(INDXMU+NCONQB,ISTATE)
            ENDDO
            UWAQUA(IBASEI,ISTATE,IREVEI) = RESU

            ! \sum_{a} B_{ia}*U_{ab} + A^*_{ia}*V_{ab}
            RESU = 0.d0
            DO INDXMU = 1, NCONQB
               RESU = RESU + conjg(BWAQUZ(IBASEI,INDXMU,IREVEI))*EWAVEF(INDXMU,ISTATE) &
                    &      + AWAQUZ(IBASEI,INDXMU,IREVEI)*EWAVEF(INDXMU+NCONQB,ISTATE)
            ENDDO
            VWAQUA(IBASEI,ISTATE,IREVEI) = RESU

            ! \sum_{a} A_{ia}*V^*_{ab} + B^*_{ia}*U^*_{ab}
            RESU = 0.d0
            DO INDXMU = 1, NCONQB
               RESU = RESU + conjg(AWAQUZ(IBASEI,INDXMU,IREVEI))*EWAVEF(INDXMU,ISTATE+NCONQB) &
                    &      + BWAQUZ(IBASEI,INDXMU,IREVEI)*EWAVEF(INDXMU+NCONQB,ISTATE+NCONQB)
            ENDDO
            UWAQUA(IBASEI,ISTATE+NCONQB,IREVEI) = RESU

            ! \sum_{a} B_{ia}*V^*_{ab} + A^*_{ia}*U^*_{ab}
            RESU = 0.d0
            DO INDXMU = 1, NCONQB
               RESU = RESU + conjg(BWAQUZ(IBASEI,INDXMU,IREVEI))*EWAVEF(INDXMU,ISTATE+NCONQB) &
                    &      + AWAQUZ(IBASEI,INDXMU,IREVEI)*EWAVEF(INDXMU+NCONQB,ISTATE+NCONQB)
            ENDDO
            VWAQUA(IBASEI,ISTATE+NCONQB,IREVEI) = RESU
         ENDDO
      ENDDO
   ENDDO

   call active_I_SING(1)
   !print*, 'I_SING is activated ...'

   DO ISTEP = 1, -1, -2
      if(IADB_DEBUG .eq. 1) print*, 'ISTEP = ', ISTEP
      AWAQUZ = 0.d0; BWAQUZ = 0.d0
      IF (ISTEP .eq. 1) then ! B -> A'
         AWAQUZ(1:NDBASE,1:NCONQB,0:NDREVE) = VWAQUA(1:NDBASE,1:NCONQB,0:NDREVE)
         BWAQUZ(1:NDBASE,1:NCONQB,0:NDREVE) = UWAQUA(1:NDBASE,1:NCONQB,0:NDREVE)
         V_SING(1:NCONQB) = EIGVAL(1:NCONQB)
      else ! B -> B'
         AWAQUZ(1:NDBASE,1:NCONQB,0:NDREVE) = VWAQUA(1:NDBASE,NCONQB+1:2*NCONQB,0:NDREVE)
         BWAQUZ(1:NDBASE,1:NCONQB,0:NDREVE) = UWAQUA(1:NDBASE,NCONQB+1:2*NCONQB,0:NDREVE)
         V_SING(1:NCONQB) = EIGVAL(NCONQB+1:2*NCONQB)
      ENDIF
      DO IREVER=0,1
         DO ISTATE=1,NUMBQP(IREVER,ICHARG)
            NOTOCC(ISTATE,IREVER) = 0
         ENDDO
      ENDDO
      NOTOCC = 0

      Call HFODD_SUBS(ICHARG,10*ISTEP)
      Call HFODD_SUBS(ICHARG,20*ISTEP)

   ENDDO

   if(ICHARG.eq.1) Call HFODD_SUBS(ICHARG,3)
   Call HFODD_SUBS(ICHARG,6)


   DEALLOCATE(CELMTS,EIGVAL,EWAVEF)
   DEALLOCATE(UWAQUA,VWAQUA)

ENDDO ! ICHARG

! -- ''Special'' treament to the time-even density  --
! -- will be necessary for density-dependent forces --
!-----------------------------------------------------
call RS_even_density(2)
! -- ''Special'' treament to the time-even density  --
!-----------------------------------------------------

Call HFODD_SUBS(-1,6)
!print*, 'Calculating FIELDS'
Call HFODD_SUBS(-1,1)
Call HFODD_SUBS(-1,4)


!print*, 'Calculating E1 ...'
!E_1 = (A^{\dagger}h-B^{\dagger}\Delta^*)B^* + (A^{\dagger}\Delta-B^{\dagger} h^*)A^*.
DO ICHARG = 0, NDISOS

   Call HFODD_SUBS(ICHARG,2)
   Call HFODD_SUBS(ICHARG,5)


   AWAQUZ(:,:,:)=ASVQUA(:,:,:,ICHARG)
   BWAQUZ(:,:,:)=BSVQUA(:,:,:,ICHARG)

   NCONQB = NUMBQP(0,ICHARG) + NUMBQP(1,ICHARG)
   NCONQB = NCONFX(ICHARG)
   allocate(RESUXX(NCONQB,LDBASE))

   if(IADB_DEBUG.eq.1) then
      if(ifammode .eq. -1) print*, 'subtracting Fermi energy ...',E_FERMI(ICHARG)
   endif
   IF(ifammode .eq. -1) then
      DO IBASEI = 1, LDBASE
         BIG_PP(IBASEI,IBASEI,0) = BIG_PP(IBASEI,IBASEI,0) - E_FERMI(ICHARG)
         BIG_PP(IBASEI,IBASEI,1) = BIG_PP(IBASEI,IBASEI,1) - E_FERMI(ICHARG)
      ENDDO
   endif

   if(imat_hh .eq. 0) then
     ! A^{\dagger} h B^*
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, BIG_PP(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,0), NDBASE,              &
     &              czero, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, BIG_PM(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, BIG_PM(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, BIG_PP(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

     ! - B^{\dagger} \Delta^* B^* = B^{\dagger} \Delta^{\dagger} B^*
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, HAMIDE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, HAMIDE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, HAM2DE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, BWAQUZ(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, HAM2DE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, BWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

     ! A^{\dagger} \Delta A^*
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, HAMIDE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, AWAQUZ(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, HAM2DE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, AWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, HAMIDE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, AWAQUZ(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, HAM2DE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, AWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

     ! -B^{\dagger} h^T A^*
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, BIG_PP(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, AWAQUZ(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, BIG_PM(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, AWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, BIG_PM(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, AWAQUZ(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, BIG_PP(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, AWAQUZ(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

!---------------------------------------------------------------------------------
!  --- reproducing HFB results ---
   else if(imat_hh .eq. 3) then
     ALLOCATE(UWAQUA(1:NDBASE,1:2*NDSTAT,0:NDREVE))
     ALLOCATE(VWAQUA(1:NDBASE,1:2*NDSTAT,0:NDREVE))
     UWAQUA = conjg(AWAQUZ) ! A
     VWAQUA = conjg(BWAQUZ) ! B
     ! A^{\dagger} h A
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, BIG_PP(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, UWAQUA(1,1,0), NDBASE,              &
     &              czero, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, BIG_PM(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, UWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, BIG_PM(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, UWAQUA(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, BIG_PP(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, UWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

     ! - B^{\dagger} \Delta^* A = B^{\dagger} \Delta^{\dagger} A
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, HAMIDE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, UWAQUA(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, HAMIDE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, UWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, HAM2DE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, UWAQUA(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'C', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, HAM2DE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, UWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

     ! A^{\dagger} \Delta B
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, HAMIDE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, VWAQUA(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,0), NDBASE, HAM2DE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, VWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, HAMIDE(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, VWAQUA(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'N', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, AWAQUZ(1,1,1), NDBASE, HAM2DE(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              c_one, RESUXX, NCONQB, VWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

     ! -B^{\dagger} h^T B
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, BIG_PP(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, VWAQUA(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,0), NDBASE, BIG_PM(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, VWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, BIG_PM(1,1,0), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, VWAQUA(1,1,0), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)
         call ZGEMM('T', 'T', NCONQB, LDBASE, LDBASE,                          &
     &              c_one, BWAQUZ(1,1,1), NDBASE, BIG_PP(1,1,1), NDBASE,       &
     &              czero, RESUXX, NCONQB)
         call ZGEMM('N', 'N', NCONQB, NCONQB, LDBASE,                          &
     &              cmone, RESUXX, NCONQB, VWAQUA(1,1,1), NDBASE,              &
     &              c_one, H_CORR(1,1,ICHARG), NCONQX)

     DEALLOCATE(UWAQUA,VWAQUA)

   endif
!---------------------------------------------------------------------------------
   DEALLOCATE(RESUXX)
   !print*, 'E1 is done ...'

ENDDO ! ICHARG

if(ifammode .eq. 1) then
   H_CORR = H_CORR/fam_eta
endif

if(ifammode .eq. -1) STOP 'reproducing HFB results ...'

endsubroutine Z_to_E1

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

      integer i
      real mass_old, mass_res
      ! == IADB_DEBUG
      integer ICHARG

! --- Initialization
      call read_input_para
      call atdhf_data_init(0)

! --- Cranking approximation
      call Inglis_mass

! --- Main iteration
      mass_old = -9999

      gamma1 = 0.d0
      call gamma1_to_rho1
      do i = 1, max_iteration

          if(i .gt. 1) then
             rho1 = xchi*rho1 + (1-xchi)*rho1_old
          endif
          rho1_old = rho1
          call rho1_to_gamma1
          call gamma1_to_rho1
          !call broyden_mix_rho1(i)
          call ATDHF_mass
          mass_res = mass_old - mass
          if(abs(mass_old - mass) < atd_precision .and. i .gt. 15) exit
          if(IADB_DEBUG .eq. 1 .or. IADB_DEBUG .eq. -1) then
             write(*,'(a,i4,a,f12.5,a,f12.5,23X,1H*)') &
             &'# iter = ', i, ' mass = ', mass, ' residual ',mass_old - mass
          endif
          mass_old = mass
      enddo
      if(i > max_iteration) then
         write(*,'(a,i4,a,f12.5,11X,1H*)') &
         &'* Calculation terminated after ', i-1, ' steps, residual is ', mass_res
      endif


! --- Generating output & clean up
      call atdhf_output
      call atdhf_data_init(1)

      endsubroutine atdhf_implementation

!===============================================================================
! active the switch I_SING to feed DENSHF with singular values & vectors
      subroutine active_I_SING(ich)
      use hfodd_gs_data
      implicit none

      integer ich

      if(ich .eq. 1) then
         V_OCCU = 0.d0
         I_SING = 1
         V_SING = 1.d0
      endif

      if(ich .eq. 0) I_SING = 0

      endsubroutine active_I_SING

!===============================================================================
      subroutine read_input_para
      use hfodd_gs_data, only : NFIPRI,IABNTE,IABCOL,IABICH,IABMOD,DABPRS,XCHIMX,IFAMON,ETAFAM
      implicit none

      ! ---- ATDHF input -----------------------------------------------------------------------
      IADB_DEBUG = -1    ! 1: Debug mode is actived; -1: only iteration summary is shown
      adj_rho = 1.d0     ! adjustable parameter multiplying to rho1   (DEBUG)
      adj_gamma= 1.d0    ! adjustable parameter multiplying to gamma1 (DEBUG)
      fam_eta = 1.d-6    ! Finite amplitude difference parameter
      ifammode = 0       ! 1: FAM-ATDHF; 0: linearilized ATDHF; -1: HF iteration
      i_svd_on = 0       ! 1: Singular value decomposition is used
      imat_hh = 0        ! 0: PH-, 1: HH-, 2: PP- matrix elements of Hamiltonian (1,2 for DEBUG)
      ecut_adb = 60.0    ! particle-energy cutoff in ATDHF calculation
      iset_num = 0       ! 0 -- JY; 1 -- JZ; 2 -- JX
      xchi = 0.5d0       ! mixing parameter for ATDHF(B) iteration
      ! --- Setup iteration precision & maximum
      atd_precision = 1.d-6
      max_iteration = 099
      ! ---- ATDHF input -----------------------------------------------------------------------

   ! --- set ATDHF(B) Parameters according to datafile,
   !   DABPRS: precision
   !   XCHIMX: mixing parameter for ATDHF(B) iteration
   !   IABNTE: number of iteration
   !   IABICH: rotating axis
   !   IABCOL: type of collective motion
   !   IABMOD: runmode, 0: admin actived; 1: admin disabled
   !   IFAMON: FAM is used (1) or not (0)
   !   ETAFAM: fam_eta

   if(IABMOD.eq.1) then
      max_iteration = IABNTE
      atd_precision = DABPRS
      xchi = XCHIMX
      !IABCOL=1
      iset_num = IABICH
      !
      ifammode = IFAMON
      fam_eta = ETAFAM
   endif

   if(IABMOD.eq.0) then
! --- read ATDHF input parameters from data file
      read(*,*) ifammode, i_svd_on, imat_hh, IADB_DEBUG
      read(*,*) adj_rho,adj_gamma
      read(*,*) fam_eta
      read(*,*) ecut_adb
      read(*,*) iset_num
      read(*,*) xchi
   endif

! --- compatibility check
      if(ifammode.eq.1 .and. i_svd_on.eq.1) stop 'ifammode = 1 requires i_svd_on = 0'

! --- printing the information of ATDHFB/QFAM calculations
      if(ifammode.eq.1) then
           WRITE(NFIPRI,'(79(1H*),/,                1H*,77X,1H*,/,  &
      &  1H*,2X,''FINITE AMPLITUDE METHOD WILL BE USED:  '',        &
      &                                                 ''ETA = '', &
      &                                            E12.5,18X,1H*,/, &
      &                                             1H*,77X,1H*)')  &
      &  fam_eta
      endif

      endsubroutine read_input_para

!===============================================================================
! 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
      integer NDIMPE,NDIMHE

      idsi_flag = IADB_DEBUG

! --- number of single-particle states
      NUMBCS = NPTOTL

! ---  Set the mass term to 0
      if(ifammode .eq. 0) HBMRPA = 0.d0

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

      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))
         allocate(rho1_old(M,N,0:NDISOS))
         call RS_even_density(0)

         ! --- reset skfild, remove the effect of linear mixing
         Call HFODD_SUBS(-1,1)
         !DE_RHO_SAV = DE_RHO
         call RS_even_density(1)

         ! --- extract particle/hole energy
         DO ICHARG=0,NDISOS

            NDIMPE = NUMBPE(ICHARG)
            NDIMHE = NUMBHE(ICHARG)
            elh(1:NDIMHE,ICHARG) = SPENER(1:NDIMHE,ICHARG)
            elp(1:NDIMPE,ICHARG) = SPENER(NDIMHE+1:NDIMHE+NDIMPE,ICHARG)

         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
                     print*, 'ICHARG, ISTATP, ISTATH, elp(ISTATP,ICHARG), elh(ISTATH,ICHARG)'
                     print*, ICHARG, ISTATP, ISTATH, elp(ISTATP,ICHARG), elh(ISTATH,ICHARG)
                     stop 'Unsolved Crossing detected in EPH!'
                  endif
               enddo
            enddo
         ENDDO

         ! --- calculate the derivative of density
         call calculate_rho_derivative

      endif

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

         ! --- free memory space...
         deallocate(rho1, gamma1, rho_derivative, eph, S, U, V)
         deallocate(elp,elh,vec,BAVOCC)
         deallocate(rho1_old)
      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

   if(i_svd_on .eq. 1) then

      vec = (0.d0, 0.d0)

      DO ICHARG=0,NDISOS
         NDIMPE = NUMBPE(ICHARG)
         NDIMHE = NUMBHE(ICHARG)

         ! Performing singular value decomposition to rho1
         call zcallsvd(NDIMPE,NDIMHE,rho1(1,1,ICHARG),  M,              &
     &                 U(1,1,ICHARG), M, V(1,1,ICHARG), N,              &
     &                 S(1,ICHARG), N)


         ! 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) =  &
     &                     CONJG(TRANSPOSE( V(1:NDIMHE,1:NDIMHE,ICHARG) ))


      ENDDO

      if(IADB_DEBUG .eq. 1) print*, 'VEC is devided by the normalization factor 2'
      vec = vec/sqrt(2.d0)

   endif
      ! Feeding DENSHF with vec to get the Gamma1
      call fed_DENSHF

      if(ifammode .eq. 1) then
         gamma1 = gamma1/fam_eta
      endif

      endsubroutine rho1_to_gamma1

!===============================================================================
! update rho1 with gamma1
      subroutine gamma1_to_rho1
      use hfodd_gs_data, only : NDISOS
      implicit none
      integer ICHARG, NDIMHE, NDIMPE

      !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

      endsubroutine gamma1_to_rho1

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

      integer ISTATP, ISTATH
      integer ICHARG

      complex resu
      complex resu_half
      real mass_isos(0:NDISOS)

      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)

      ENDDO

      endsubroutine ATDHF_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(:,:)
      complex,allocatable :: FAUXPP(:,:), FAUXPM(:,:), FAUXMP(:,:), FAUXMM(:,:)
      real,   allocatable :: HAUXTT(:,:)
      integer IBASEI, IBASEJ, ICHARG, ISTATP, ISTATH
      complex RESUPP, RESUPM, RESUMP, RESUMM
      integer NDIMHE
      real    ppres, pmres, mpres, mmres, ttres
      real    res_max_PP,res_max_PM,res_max_MP,res_max_MM,res_max_TT
      complex J_ISO(0:1), J_TOTAL
      complex i_unit

      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))
      allocate(FAUXPP(1:NDBASE,1:NDBASE),FAUXPM(1:NDBASE,1:NDBASE))
      allocate(FAUXMP(1:NDBASE,1:NDBASE),FAUXMM(1:NDBASE,1:NDBASE))

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

      i_unit = (0.d0,1.d0)

      ! iset_num: 0 -- JY; 1 -- JZ; 2 -- JX;
      axis_str(0) = 'Y-'; axis_str(1) = 'Z-'; axis_str(2) = 'X-'

   FACTOR = 1.d0
   if(iset_num == 0) then
      ! -- set 0 --
      if(IADB_DEBUG .eq. 1) print*, 'Matrix element of JY is used'
      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)

      FAUXPP = HAUXPP*i_unit + HAUXTT
      FAUXPM = HAUXPM*i_unit
      FAUXMP = HAUXMP*i_unit
      FAUXMM = HAUXMM*i_unit - HAUXTT

   else if(iset_num == 1) then
      ! -- set 1 --
      if(IADB_DEBUG .eq. 1) print*, 'Matrix element of JZ is used'
      CALL INT_SZ(HAUXTT,FACTOR)
      CALL INT_LZ(IPHAPP(0,0,0),HAUXPP,FACTOR)
      CALL INT_LZ(IPHAPM(0,0,0),HAUXPM,FACTOR)
      CALL INT_LZ(IPHAMP(0,0,0),HAUXMP,FACTOR)
      CALL INT_LZ(IPHAMM(0,0,0),HAUXMM,FACTOR)

      FAUXPP = HAUXPP*i_unit
      FAUXPM = HAUXPM*i_unit + HAUXTT*i_unit
      FAUXMP = HAUXMP*i_unit - HAUXTT*i_unit
      FAUXMM = HAUXMM*i_unit

   else if(iset_num == 2) then
      ! -- set 2 --
      if(IADB_DEBUG .eq. 1) print*, 'Matrix element of JX is used'
      CALL INT_SX(HAUXTT,FACTOR)
      CALL INT_LX(IPHAPP(0,0,0),HAUXPP,FACTOR)
      CALL INT_LX(IPHAPM(0,0,0),HAUXPM,FACTOR)
      CALL INT_LX(IPHAMP(0,0,0),HAUXMP,FACTOR)
      CALL INT_LX(IPHAMM(0,0,0),HAUXMM,FACTOR)

      FAUXPP = HAUXPP*i_unit
      FAUXPM = HAUXPM*i_unit + HAUXTT
      FAUXMP = HAUXMP*i_unit + HAUXTT
      FAUXMM = HAUXMM*i_unit

   endif ! iset_num

      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) )   &
     &                               * FAUXPP(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,0)
                     RESUPM = RESUPM + conjg( WAVOCC(IBASEI,ISTATP+NDIMHE,0) )   &
     &                               * FAUXPM(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,1)
                     RESUMP = RESUMP + conjg( WAVOCC(IBASEI,ISTATP+NDIMHE,1) )   &
     &                               * FAUXMP(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,0)
                     RESUMM = RESUMM + conjg( WAVOCC(IBASEI,ISTATP+NDIMHE,1) )   &
     &                               * FAUXMM(IBASEI,IBASEJ)                     &
     &                               * WAVOCC(IBASEJ,ISTATH,1)
                  ENDDO
               ENDDO
               rho_derivative(ISTATP,ISTATH,ICHARG) = RESUPP + RESUPM + RESUMP + RESUMM
            ENDDO
         ENDDO
      ENDDO ! For ICHARG

      deallocate(HAUXPP,HAUXPM,HAUXMP,HAUXMM,HAUXTT)
      deallocate(FAUXPP,FAUXPM,FAUXMP,FAUXMM)
      endsubroutine calculate_rho_derivative

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

      !rho1 = rho_derivative/eph
      gamma1 = 0.d0
      call gamma1_to_rho1
      call ATDHF_mass


           WRITE(NFIPRI,'(79(1H*),/,                1H*,77X,1H*,/,  &
      &  1H*,2X,''Inglis Formula MOMENT OF INERTIA ALONG '',        &
      &                                               A2,''AXIS '', &
      &                         ''IS '',F12.5,'' MeV-1'', 8X,1H*,/, &
      &                                             1H*,77X,1H*)')  &
      &    axis_str(iset_num), mass

      endsubroutine Inglis_mass

!===============================================================================
      subroutine atdhf_output
      use hfodd_gs_data
      implicit none

      ! generating output of ATDHF calculations

          WRITE(NFIPRI,'(79(1H*),/,                1H*,77X,1H*,/,  &
     &           1H*,2X,''ATDHF MOMENT OF INERTIA ALONG '',        &
     &                                               A2,''AXIS '', &
     &                         ''IS '',F12.5,'' MeV-1'',17X,1H*,/, &
     &                                             1H*,77X,1H*)')  &
     &    axis_str(iset_num), mass
      endsubroutine atdhf_output

!===============================================================================
      subroutine atdhfb_output
      use hfodd_gs_data
      implicit none

      ! generating output of ATDHF calculations

          WRITE(NFIPRI,'(79(1H*),/,                1H*,77X,1H*,/,  &
     &           1H*,2X,''ATDHFB MOMENT OF INERTIA ALONG '',       &
     &                                               A2,''AXIS '', &
     &                         ''IS '',F12.5,'' MeV-1'',16X,1H*,/, &
     &                                             1H*,77X,1H*)')  &
     &    axis_str(iset_num), mass
      endsubroutine atdhfb_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, jccount
      integer NUMBPX, NUMBHX


      DO ICHARG=0,NDISOS
         iccount = 0; jccount = 0
         DO ISTATE=1,NUMBCS(ICHARG)
            if(M_OCVA(ISTATE,ICHARG) .eq. 1) iccount = iccount + 1
            IF(SPENER(ISTATE,ICHARG) .gt. ecut_adb) EXIT
            jccount = jccount + 1
         ENDDO
         NUMBCS(ICHARG) = jccount
         NUMBHE(ICHARG) = iccount
         NUMBPE(ICHARG) = NUMBCS(ICHARG) - NUMBHE(ICHARG)
      ENDDO

      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
      !###### 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

!     occupation number of the time-odd density
      S = sqrt(S)

      DO ICHARG=0,NDISOS

      if(i_svd_on .eq. 1) then

         call active_I_SING(1)
! --- replace V_SING by r_i
         DO ISTATH = 1, NUMBHE(ICHARG)
            V_SING(2*ISTATH - 1) = -1.d0
            V_SING(2*ISTATH    ) = +1.d0
            V_OCCU(2*ISTATH - 1) = S(ISTATH,ICHARG)
            V_OCCU(2*ISTATH    ) = S(ISTATH,ICHARG)
         ENDDO

! --- replacing WAVOCC
         BAVOCC = 0.d0
         WAVOCC = SAVOCC(:,:,:,ICHARG)

         ! -- $ 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,NUMBCS(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

      else ! i_svd_on
         call REPLACE_WAVOCC(ICHARG)
      endif ! i_svd_on

         Call HFODD_SUBS(ICHARG,10)

         ! calculate time-odd Coulomb field
         if(ICHARG.eq.1) then
            Call HFODD_SUBS(ICHARG,3)
         endif
      ENDDO ! For ICHARG

      call RS_even_density(2)

      Call HFODD_SUBS(-1,1)

! --- extract matrix elements of Gamma1
! -   $ \Gamma_{ph} = sum_{ij}\varphi_{ip}^* \Gamma_{ij} \varphi_{jh} $
      gamma1 = (0.d0,0.d0)
      DO ICHARG = 0, NDISOS

         Call HFODD_SUBS(ICHARG,2)

         WAVOCC = SAVOCC(:,:,:,ICHARG)
         NDIMPE = NUMBPE(ICHARG)
         NDIMHE = NUMBHE(ICHARG)


      if(imat_hh .eq. 0) then
         allocate(RESUXX(NDIMPE,LDBASE))
         ! -- 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)
         ! -- 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)
     !---------- HH ------------
     else if (imat_hh .eq. 1) then
         allocate(RESUXX(NDIMHE,LDBASE))
         ! -- WF(+) * GAMMA(++) * WF(+)
         call ZGEMM('C', 'N', NDIMHE, LDBASE, LDBASE,                          &
     &              c_one, WAVOCC(1,1,0), NDBASE, BIG_PP(1,1,0), NDBASE,&
     &              czero, RESUXX, NDIMHE)
         call ZGEMM('N', 'N', NDIMHE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMHE, WAVOCC(1,1,0), NDBASE,              &
     &              czero, gamma1(1,1,ICHARG), M)
         ! -- WF(+) * GAMMA(+-) * WF(-)
         call ZGEMM('C', 'N', NDIMHE, LDBASE, LDBASE,                          &
     &              c_one, WAVOCC(1,1,0), NDBASE, BIG_PM(1,1,0), NDBASE,&
     &              czero, RESUXX, NDIMHE)
         call ZGEMM('N', 'N', NDIMHE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMHE, WAVOCC(1,1,1), NDBASE,              &
     &              c_one, gamma1(1,1,ICHARG), M)
         ! -- WF(-) * GAMMA(-+) * WF(+)
         call ZGEMM('C', 'N', NDIMHE, LDBASE, LDBASE,                          &
     &              c_one, WAVOCC(1,1,1), NDBASE, BIG_PM(1,1,1), NDBASE,&
     &              czero, RESUXX, NDIMHE)
         call ZGEMM('N', 'N', NDIMHE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMHE, WAVOCC(1,1,0), NDBASE,              &
     &              c_one, gamma1(1,1,ICHARG), M)
         ! -- WF(-) * GAMMA(--) * WF(-)
         call ZGEMM('C', 'N', NDIMHE, LDBASE, LDBASE,                          &
     &              c_one, WAVOCC(1,1,1), NDBASE, BIG_PP(1,1,1), NDBASE,&
     &              czero, RESUXX, NDIMHE)
         call ZGEMM('N', 'N', NDIMHE, NDIMHE, LDBASE,                          &
     &              c_one, RESUXX, NDIMHE, WAVOCC(1,1,1), NDBASE,              &
     &              c_one, gamma1(1,1,ICHARG), M)
     !---------- PP ------------
     else if (imat_hh .eq. 2) then
         allocate(RESUXX(NDIMPE,LDBASE))
         ! -- 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, NDIMPE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,NDIMHE+1,0), NDBASE,       &
     &              czero, gamma1(1,1,ICHARG), M)
         ! -- 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, NDIMPE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,NDIMHE+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, NDIMPE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,NDIMHE+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, NDIMPE, LDBASE,                          &
     &              c_one, RESUXX, NDIMPE, WAVOCC(1,NDIMHE+1,1), NDBASE,       &
     &              c_one, gamma1(1,1,ICHARG), M)
     endif ! imat_hh
     !---------- PP & HH END ------------
         deallocate(RESUXX)

      ENDDO

      endsubroutine fed_DENSHF

!--------------------------------------------------------------------

subroutine REPLACE_WAVOCC(ICHARG)
use hfodd_gs_data
implicit none

integer ICHARG
integer ISTATP, ISTATH, ISTATV
integer i,j
integer IBASE, ISTATE

COMPLEX, ALLOCATABLE :: CELMTS(:)
REAL, ALLOCATABLE :: EIGVAL(:)
COMPLEX, ALLOCATABLE :: EWAVEF(:,:)
COMPLEX, ALLOCATABLE :: AA(:,:)

integer LDACTU, NDACTU, NUMVEC
integer NDIMPE, NDIMHE

COMPLEX RESU0, RESU1

   NDIMPE = NUMBPE(ICHARG)
   NDIMHE = NUMBHE(ICHARG)

   NDACTU = NDIMPE + NDIMHE
   LDACTU = NDACTU
   NUMVEC = NDACTU

   ALLOCATE( AA(1:NDACTU,1:NDACTU))
   ALLOCATE( CELMTS(1:((NDACTU+1)*NDACTU)/2) )
   ALLOCATE( EIGVAL(1:NDACTU) )
   ALLOCATE( EWAVEF(1:NDACTU,1:NDACTU) )

   select case (ifammode)
   case (0)
   ! --- for linearalized ATDHF ---
      if(IADB_DEBUG .eq. 1) print*, 'linearalized ATDHF'
      AA = (0.d0, 0.d0)
      do ISTATP = 1, NDIMPE
         do ISTATH = 1, NDIMHE
            AA(ISTATP,ISTATH + NDIMPE) = rho1(ISTATP,ISTATH,ICHARG)
            AA(ISTATH + NDIMPE,ISTATP) = conjg(rho1(ISTATP,ISTATH,ICHARG))
         enddo
      enddo

   case (-1)
   ! --- for HF iteration ---
      if(IADB_DEBUG .eq. 1) print*, 'reproduction of HF iteration in ATDHF'
      AA = (0.d0, 0.d0)
      do ISTATH = 1, NDIMHE
         AA(ISTATH + NDIMPE,ISTATH + NDIMPE) = 1.d0
      enddo

   case (1)
   ! --- for FAM-ATDHF ---
      if(IADB_DEBUG .eq. 1) print*, 'FAM-ATDHF, fam_eta = ', fam_eta
      AA = (0.d0, 0.d0)
      do ISTATH = 1, NDIMHE
         AA(ISTATH + NDIMPE,ISTATH + NDIMPE) = 1.d0
      enddo
      do ISTATP = 1, NDIMPE
         do ISTATH = 1, NDIMHE
            AA(ISTATP,ISTATH + NDIMPE) = AA(ISTATP,ISTATH + NDIMPE) &
                                    &  + fam_eta*rho1(ISTATP,ISTATH,ICHARG)
            AA(ISTATH + NDIMPE,ISTATP) = AA(ISTATH + NDIMPE,ISTATP) &
                                    &  + fam_eta*conjg(rho1(ISTATP,ISTATH,ICHARG))
         enddo
      enddo

   case default
      stop 'unknown runmode ..'
   end select
   ! --------------------------------------------------------------

   DO j = 1, NDACTU
      DO i = j, NDACTU
         CELMTS(i + (j-1)*(2*NDACTU-j)/2) = AA(i,j)
      ENDDO
   ENDDO

   CALL DIAMAT(CELMTS,EIGVAL,EWAVEF,LDACTU,NDACTU,NUMVEC)


   !---------------------------------------------------------------------!
   ! The eigenstates from DIAMAT are sorted ascendingly.                 !
   ! structure: |-s, -v, ..., 0, ..., 0, ..., V, S|                      !
   !            |   NDIMHE   |---------|  NDIMHE  |                      !
   ! In linearalized ATDHF mode, s = S, v = V, S = V;                    !
   ! In FAM-ATDHF mode, s = S - 1, v = V - 1; (since rho_hh = 1)         !
   ! In ATDHF, only those with nonvanishing eigenvalues are needed.      !
   ! Since in HF iteration, only 2*NDIMHE eigenstates are used, the zero !
   ! eigenstates should be removed like this:                            !
   !---------------------------------------------------------------------!
   EIGVAL(  NDIMHE+1:2*NDIMHE) = EIGVAL(  NDACTU-NDIMHE+1:NDACTU)
   EWAVEF(:,NDIMHE+1:2*NDIMHE) = EWAVEF(:,NDACTU-NDIMHE+1:NDACTU)
   ! For security purpose,
   EIGVAL(2*NDIMHE+1:NDACTU) = 0.d0

   call active_I_SING(1)
   DO ISTATH = 1, 2*NDIMHE
      if(EIGVAL(ISTATH).GE.0) then
         V_OCCU(ISTATH) = sqrt( +EIGVAL(ISTATH) )
         V_SING(ISTATH) = +1.d0
      else
         V_OCCU(ISTATH) = sqrt( -EIGVAL(ISTATH) )
         V_SING(ISTATH) = -1.d0
      endif
   ENDDO

   BAVOCC = 0.d0
   WAVOCC = SAVOCC(:,:,:,ICHARG)

   DO ISTATH = 1, 2*NDIMHE
      DO IBASE=1,LDBASE
         RESU0 = 0.d0; RESU1 = 0.d0;
         DO ISTATE=1,NUMBCS(ICHARG)
            IF(ISTATE .LE. NDIMHE) THEN
               ISTATV = ISTATE + NUMBPE(ICHARG)
            ELSE
               ISTATV = ISTATE - NUMBHE(ICHARG)
            ENDIF

            RESU0 = RESU0 + WAVOCC(IBASE,ISTATE,0) * EWAVEF(ISTATV,ISTATH)
            RESU1 = RESU1 + WAVOCC(IBASE,ISTATE,1) * EWAVEF(ISTATV,ISTATH)
         ENDDO
         BAVOCC(IBASE,ISTATH,0) = RESU0
         BAVOCC(IBASE,ISTATH,1) = RESU1
      ENDDO
   ENDDO

   WAVOCC = BAVOCC
   DEALLOCATE(AA,CELMTS,EIGVAL,EWAVEF)

endsubroutine REPLACE_WAVOCC


!--------------------------------------------------------------------
! -- allocate [0], save [1], retrieve [2] time-even densities --
subroutine RS_even_density(ich)
use hfodd_gs_data
implicit none
integer ich
integer i_rs_level
integer LDPNMX, LDTWDD

i_rs_level = 1
IF(IGYUSD.eq.1) i_rs_level = 0

! pending to be modified for p-n mixing & two-center basis
LDPNMX=1
LDTWDD=1

if(i_rs_level .eq. 0) return

if(ich == 0) then
   if(i_rs_level .gt. 0) then
      IF(ALLOCATED(DE_RHO)) ALLOCATE(DE_RHO_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
   endif
   if(i_rs_level .gt. 1) then
      IF(ALLOCATED(DE_TAU)) ALLOCATE(DE_TAU_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(DE_LPR)) ALLOCATE(DE_LPR_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(DE_DIV)) ALLOCATE(DE_DIV_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(DE_SCU)) ALLOCATE(DE_SCU_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,1:NDKART,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(DE_DES)) ALLOCATE(DE_DES_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,1:NDKART,0:LDPNMX,1:LDTWDD))
   endif
   if(i_rs_level .gt. 2) then
      !IF(ALLOCATED(DE_DIJ)) ALLOCATE(DE_DIJ_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PD_RHO)) ALLOCATE(PD_RHO_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PD_TAU)) ALLOCATE(PD_TAU_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PD_LPR)) ALLOCATE(PD_LPR_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PP_RHO)) ALLOCATE(PP_RHO_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PP_TAU)) ALLOCATE(PP_TAU_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PP_LPR)) ALLOCATE(PP_LPR_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PD_SCU)) ALLOCATE(PD_SCU_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,1:NDKART,0:LDPNMX,1:LDTWDD))
      IF(ALLOCATED(PP_SCU)) ALLOCATE(PP_SCU_SAV(1:NDXHRM,1:NDYHRM,1:NDZHRM,1:NDKART,1:NDKART,0:LDPNMX,1:LDTWDD))
   endif
endif

if(ich == 1) then
   ! -- save time-even density --
   if(i_rs_level .gt. 0) then
      IF(ALLOCATED(DE_RHO)) DE_RHO_SAV = DE_RHO
   endif
   if(i_rs_level .gt. 1) then
      IF(ALLOCATED(DE_TAU)) DE_TAU_SAV = DE_TAU
      IF(ALLOCATED(DE_LPR)) DE_LPR_SAV = DE_LPR
      IF(ALLOCATED(DE_DIV)) DE_DIV_SAV = DE_DIV
      IF(ALLOCATED(DE_SCU)) DE_SCU_SAV = DE_SCU
      IF(ALLOCATED(DE_DES)) DE_DES_SAV = DE_DES
   endif
   if(i_rs_level .gt. 2) then
      !IF(ALLOCATED(DE_DIJ)) DE_DIJ_SAV = DE_DIJ
      IF(ALLOCATED(PD_RHO)) PD_RHO_SAV = PD_RHO
      IF(ALLOCATED(PD_TAU)) PD_TAU_SAV = PD_TAU
      IF(ALLOCATED(PD_LPR)) PD_LPR_SAV = PD_LPR
      IF(ALLOCATED(PP_RHO)) PP_RHO_SAV = PP_RHO
      IF(ALLOCATED(PP_TAU)) PP_TAU_SAV = PP_TAU
      IF(ALLOCATED(PP_LPR)) PP_LPR_SAV = PP_LPR
      IF(ALLOCATED(PD_SCU)) PD_SCU_SAV = PD_SCU
      IF(ALLOCATED(PP_SCU)) PP_SCU_SAV = PP_SCU
   endif
endif

if(ich == 2) then
   ! -- retrieve time-even density --
   if(i_rs_level .gt. 0) then
      IF(ALLOCATED(DE_RHO)) DE_RHO = DE_RHO_SAV
   endif
   if(i_rs_level .gt. 1) then
      IF(ALLOCATED(DE_TAU)) DE_TAU = DE_TAU_SAV
      IF(ALLOCATED(DE_LPR)) DE_LPR = DE_LPR_SAV
      IF(ALLOCATED(DE_DIV)) DE_DIV = DE_DIV_SAV
      IF(ALLOCATED(DE_SCU)) DE_SCU = DE_SCU_SAV
      IF(ALLOCATED(DE_DES)) DE_DES = DE_DES_SAV
   endif
   if(i_rs_level .gt. 2) then
      !IF(ALLOCATED(DE_DIJ)) DE_DIJ = DE_DIJ_SAV
      IF(ALLOCATED(PD_RHO)) PD_RHO = imag(PD_RHO) + real(PD_RHO_SAV)
      IF(ALLOCATED(PD_TAU)) PD_TAU = imag(PD_TAU) + real(PD_TAU_SAV)
      IF(ALLOCATED(PD_LPR)) PD_LPR = imag(PD_LPR) + real(PD_LPR_SAV)
      IF(ALLOCATED(PP_RHO)) PP_RHO = imag(PP_RHO) + real(PP_RHO_SAV)
      IF(ALLOCATED(PP_TAU)) PP_TAU = imag(PP_TAU) + real(PP_TAU_SAV)
      IF(ALLOCATED(PP_LPR)) PP_LPR = imag(PP_LPR) + real(PP_LPR_SAV)
      IF(ALLOCATED(PD_SCU)) PD_SCU = PD_SCU*(0.d0,1.d0) + PD_SCU_SAV*(1.d0,0.d0)
      IF(ALLOCATED(PP_SCU)) PP_SCU = PP_SCU*(0.d0,1.d0) + PP_SCU_SAV*(1.d0,0.d0)
   endif
endif

endsubroutine RS_even_density


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
      integer LDA
      integer LDU
      integer LDV
      integer LDS
      complex in_A(LDA,LDV)
      complex U(LDU,LDV)
      complex V(LDV,LDV)
      real S(LDS)
      !complex A(LDA,LDV)
      complex, allocatable :: A(:,:)

      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
      allocate( A(LDA,LDV) )
      A = in_A

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

      call zgesvdq( 'H',     'P',     'N',     'U',    'A',             &
     &               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( 'H',     'P',     'N',     'U',    'A',             &
     &               M,       N,       A,       LDA,                    &
     &               S,       U,       LDU,     V,      LDV,            &
     &               NUMRANK,                                           &
     &               IWORK,   LIWORK,                                   &
     &               CWORK,   LWORK,                                    &
     &               RWORK,   LRWORK,                                   &
     &               INFO )
! --------------------------------------------------------------------72--------

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

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

!===============================================================================
! --------------------------------------------------------------------%%--------
! --    READ/WRITE WAVEFUNCTION & MULTIPOLE MOMENTS FROM/TO FILE
! --------------------------------------------------------------------%%--------

      SUBROUTINE RWADBP(FSTR,ICH,SMUL_V,NDMULT,SAVQUA,SBVQUA,NDBASE,NDSTAT,NDREVE,NDISOS)
      IMPLICIT NONE
      CHARACTER*(*) FSTR
      INTEGER ICH
      INTEGER NDMULT,NDBASE,NDSTAT,NDREVE,NDISOS
      REAL SMUL_V(0:NDMULT,-NDMULT:NDMULT)
      COMPLEX SAVQUA(1:NDBASE,1:2*NDSTAT,0:NDREVE,0:NDISOS)
      COMPLEX SBVQUA(1:NDBASE,1:2*NDSTAT,0:NDREVE,0:NDISOS)

      INTEGER NFIADB
      CHARACTER SFIADB*72

      NFIADB = 9932

      SFIADB=TRIM(FSTR)

      OPEN(NFIADB,FILE=SFIADB,STATUS='OLD',ERR=9,FORM='UNFORMATTED')

      IF(ICH .EQ. 0) THEN
         OPEN(NFIADB,FILE=SFIADB,FORM='UNFORMATTED')
         WRITE (NFIADB) SMUL_V
         WRITE (NFIADB) SAVQUA
         WRITE (NFIADB) SBVQUA
      ELSE
         OPEN(NFIADB,FILE=SFIADB,STATUS='OLD',ERR=9,FORM='UNFORMATTED')
         READ (NFIADB) SMUL_V
         READ (NFIADB) SAVQUA
         READ (NFIADB) SBVQUA
      ENDIF

      CLOSE(NFIADB)

      RETURN

    9 WRITE(*,'(/,1X,26(1H/),2X,''FILE NOT FOUND ON DISC'',              &
     &                                                 2X,26(1H/),/,     &
     &                 1X, 3(1H/),2X,68X,                  2X, 3(1H/),/, &
     &                 1X, 1(1H/),2X,A72,                  2X, 1(1H/),/, &
     &                 1X, 3(1H/),2X,68X,                  2X, 3(1H/),/, &
     &                 1X,26(1H/),2X,''FILE NOT FOUND ON DISC'',         &
     &                                                 2X,26(1H/))')     &
     &             SFIADB

      STOP ' FILE NOT FOUND IN RWADBP'

      ENDSUBROUTINE

! ------------------------------------------------------------------------------
