qp_plugins_scemama/devel/svdwf/psiSVD_pt2_v0.irp.f

program psiSVD_pt2_v0

  implicit none

  BEGIN_DOC
  ! perturbative approach to build psi_postsvd
  END_DOC

  read_wf = .True.
  TOUCH read_wf

  PROVIDE N_int

  call run()
end


subroutine run

  USE OMP_LIB

  implicit none

  integer(bit_kind)             :: det1(N_int,2), det2(N_int,2)
  integer                       :: degree, i_state

  integer                       :: i, j, k, l, m, n
  double precision              :: x, y, h12

  double precision, allocatable :: Uref(:,:), Dref(:), Vtref(:,:), Aref(:,:), Vref(:,:)

  integer                       :: rank_max
  double precision              :: E0, overlop, Ept2
  double precision, allocatable :: H0(:,:)
  double precision, allocatable :: eigvec0(:,:), eigval0(:), coeff_psi(:), coeff_tmp(:)

  integer                       :: ii, jj, ia, ib
  double precision, allocatable :: Hdiag(:), Hkl_save(:,:), Hkl_1d(:), Hkl_tmp(:,:), Hdiag_tmp(:)
  double precision, allocatable :: H0_1d(:), H0_tmp(:,:)

  integer                       :: na_new, nb_new, ind_new, ind_gs
  double precision              :: ctmp, coeff_new
  double precision, allocatable :: epsil(:), epsil_energ(:), check_ov(:)

  double precision, allocatable :: Uezfio(:,:,:), Dezfio(:,:), Vezfio(:,:,:)

  integer                       :: n_selected, n_toselect, n_tmp, na_max, nb_max
  integer,          allocatable :: numalpha_selected(:), numbeta_selected(:)
  integer,          allocatable :: numalpha_toselect(:), numbeta_toselect(:)
  integer,          allocatable :: numalpha_tmp(:), numbeta_tmp(:)

  integer                       :: cantor_pairing_ij, cantor_pairing_new
  integer, allocatable          :: cantor_pairing(:), cantor_pairing_tmp(:)

  double precision              :: t_beg, t_end
  integer(kind=8)               :: W_tbeg, W_tend, W_tbeg_it, W_tend_it, W_ir
  real(kind=8)                  :: W_tot_time, W_tot_time_it
  integer                       :: nb_taches

  !$OMP PARALLEL
  nb_taches = OMP_GET_NUM_THREADS()
  !$OMP END PARALLEL

  call SYSTEM_CLOCK(COUNT=W_tbeg, COUNT_RATE=W_ir)

  i_state = 1

  det1(:,1) = psi_det_alpha_unique(:,1)
  det2(:,1) = psi_det_alpha_unique(:,1)
  det1(:,2) = psi_det_beta_unique(:,1)
  det2(:,2) = psi_det_beta_unique(:,1)
  call get_excitation_degree_spin(det1(1,1),det2(1,1),degree,N_int)
  call get_excitation_degree(det1,det2,degree,N_int)
  call i_H_j(det1, det2, N_int, h12)

  ! ---------------------------------------------------------------------------------------
  !                      construct the initial CISD matrix

  print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'
  print *, ' CISD matrix:', n_det_alpha_unique,'x',n_det_beta_unique
  print *, ' N det      :', N_det
  print *, ' ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~'

  allocate( Aref(n_det_alpha_unique,n_det_beta_unique) )
  Aref(:,:) = 0.d0
  do k = 1, N_det
    i         = psi_bilinear_matrix_rows(k)
    j         = psi_bilinear_matrix_columns(k)
    Aref(i,j) = psi_bilinear_matrix_values(k,i_state)
  enddo

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



  ! ---------------------------------------------------------------------------------------
  !                                  perform a Full SVD

  allocate( Uref(n_det_alpha_unique,n_det_alpha_unique)     )
  allocate( Dref(min(n_det_alpha_unique,n_det_beta_unique)) )
  allocate( Vtref(n_det_beta_unique,n_det_beta_unique)      )
  call cpu_time(t_beg)
  call svd_s(Aref, size(Aref,1), Uref, size(Uref,1), Dref, Vtref &
              , size(Vtref,1), n_det_alpha_unique, n_det_beta_unique)
  call cpu_time(t_end)
  print *, " SVD is performed after (min)", (t_end-t_beg)/60.

  allocate( Vref(n_det_beta_unique,n_det_beta_unique) )
  do l = 1, n_det_beta_unique
    do i = 1, n_det_beta_unique
      Vref(i,l) = Vtref(l,i)
    enddo
  enddo
  deallocate( Vtref )
  deallocate( Aref )

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

  !  ***   PARAMETERS   ***  !
  na_max = n_det_alpha_unique
  nb_max = n_det_beta_unique
  !  ***     *****      ***  !

  print *, ' na_max  = ', na_max
  print *, ' nb_max  = ', nb_max

  ! ---------------------------------------------------------------------------------------
  !                            initial wavefunction: psi_0

  n_selected = 1
  allocate(numalpha_selected(n_selected), numbeta_selected(n_selected), cantor_pairing(n_selected))

  numalpha_selected(1) = 1
  numbeta_selected (1) = 1
  cantor_pairing   (1) = 4 !int( 0.5*(1+1)*(1+1+1) ) + 1

  allocate( coeff_psi(n_selected) )
  coeff_psi(1) = 1.d0

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



  ! ---------------------------------------------------------------------------------------
  !         construnc the initial basis to select phi_1 from the FSVD

  n_toselect = na_max * nb_max - n_selected
  print *, ' toselect = ', n_toselect
  print *, ' to trun  = ', n_det_alpha_unique*n_det_beta_unique - na_max*nb_max

  allocate( numalpha_toselect(n_toselect) , numbeta_toselect(n_toselect) )
  k = 0
  do i = 1, na_max
    do j = 1, nb_max

      cantor_pairing_ij = int( 0.5*(i+j)*(i+j+1) ) + j
      if( ANY(cantor_pairing .eq. cantor_pairing_ij) ) cycle

      k = k + 1
      numalpha_toselect(k) = i
      numbeta_toselect (k) = j

    enddo
  enddo
  if( k.ne.n_toselect ) then
    print *, " error in chosing vectors toselect"
    print *, " n_toselect =", n_toselect
    print *, " k          =", k
    stop
  endif

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



  ! ---------------------------------------------------------------------------------------
  !                        read < u_k v_l | H | u_k v_l > for all vectors 

  allocate( Hdiag(n_toselect) , H0(n_selected,n_selected) )
 
  open( unit=11, FILE="klHkl_v0.dat", ACTION="READ")

    read(11,*) i, i, E0
    H0(1,1) = E0

    do i = 1, n_toselect
      read(11,*) ia, ib, ctmp 
      if( (numalpha_toselect(i).ne.ia) .or. (numbeta_toselect(i).ne.ib) ) then
        print *, ' error in reading klHkl_v0 '
        print *, ia, ib
        print *, numalpha_toselect(i), numbeta_toselect(i)
        stop
      endif
      Hdiag(i) = ctmp
    enddo

  close(11)

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

  E0 = E0 + nuclear_repulsion
  print*, ' space dimen     = ', n_selected
  print*, ' ground state E0 = ', E0

  na_new = 1
  nb_new = 1

  !________________________________________________________________________________________________________
  !
  !                              increase the size of psi0 iteratively
  !________________________________________________________________________________________________________

  !  ***   PARAMETERS   ***  !
  rank_max = na_max * nb_max
  !rank_max = 50 * 50
  !  ***     *****      ***  !

  if( rank_max .gt. (na_max*nb_max) ) then
    print *, " rank_max should be less then na_max x nb_max"
    stop
  endif


  allocate( Hkl_save(n_toselect,n_selected) )

  do while( n_selected .lt. rank_max )

    call SYSTEM_CLOCK(COUNT=W_tbeg_it, COUNT_RATE=W_ir) 

    print*, '                         '
    print*, '       new iteration     '

    if( n_toselect .lt. 1 ) then

      print*, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
      print*, ' no more vectors to construct a new basis '
      print*, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
      stop

    else

      ! ---------------------------------------------------------------------------------------
      !                                select a new vector

      allocate( Hkl_1d(n_toselect)  )
      call const_Hkl_1d(na_new, nb_new, na_max, nb_max, n_toselect, Uref, Vref, numalpha_toselect, numbeta_toselect, Hkl_1d)
      Hkl_save(1:n_toselect,n_selected) = Hkl_1d(1:n_toselect)
      deallocate( Hkl_1d )

      ! choose the best vector
      allocate( epsil(n_toselect) , epsil_energ(n_toselect) ) 
      do ii = 1, n_toselect

        ctmp = 0.d0 
        do l = 1, n_selected
          ctmp = ctmp + coeff_psi(l) * Hkl_save(ii,l)
        enddo
        epsil(ii) = ctmp * ctmp / ( E0 - (Hdiag(ii)+nuclear_repulsion) )

        epsil_energ(ii) = epsil(ii)
        epsil(ii)       = dabs( epsil(ii) ) 
      enddo

      ind_new = MAXLOC( epsil, DIM=1 )

      ept2 = epsil_energ(ind_new) 
      if( ept2 .gt. 0.d0 ) then
        print *, ' ept2 > 0 !!!!!!!!!! '
        print *, na_new, nb_new, ept2
        stop
      endif

      na_new             = numalpha_toselect(ind_new) 
      nb_new             = numbeta_toselect (ind_new) 
      cantor_pairing_new = int( 0.5 * (na_new+nb_new) * (na_new+nb_new+1) ) + nb_new

      print *, ' best vector', na_new, nb_new, ept2
      deallocate(epsil,epsil_energ)

      ! new coefficient
      coeff_new = 0.d0
      do l = 1, n_selected
        coeff_new += coeff_psi(l) * Hkl_save(ind_new,l)
      enddo
      coeff_new = coeff_new / ( E0 - (Hdiag(ind_new)+nuclear_repulsion) )
      print *, ' new coeff = ', coeff_new

      ! < psi_old | H | delta_psi > 
      allocate( H0_1d(n_selected) )
      call const_H0_1d(na_new, nb_new, na_max, nb_max, n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0_1d)

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



      ! ---------------------------------------------------------------------------------------
      !                                      new psi

      allocate( numalpha_tmp(n_selected), numbeta_tmp(n_selected), coeff_tmp(n_selected) )
      allocate( cantor_pairing_tmp(n_selected) )
      allocate( H0_tmp(n_selected,n_selected)  )

      coeff_tmp         (:) = coeff_psi        (:)
      numalpha_tmp      (:) = numalpha_selected(:)
      numbeta_tmp       (:) = numbeta_selected (:)
      cantor_pairing_tmp(:) = cantor_pairing   (:)
      H0_tmp          (:,:) = H0             (:,:)

      deallocate( numalpha_selected, numbeta_selected, coeff_psi, cantor_pairing, H0 )

      n_tmp      = n_selected
      n_selected = n_selected + 1

      allocate( numalpha_selected(n_selected) , numbeta_selected(n_selected) , coeff_psi(n_selected) )
      allocate( cantor_pairing(n_selected) )
      allocate( H0(n_selected,n_selected)  )
      H0(:,:) = 0.d0

      do l = 1, n_tmp
        coeff_psi        (l) = coeff_tmp         (l)
        numalpha_selected(l) = numalpha_tmp      (l)
        numbeta_selected (l) = numbeta_tmp       (l)
        cantor_pairing   (l) = cantor_pairing_tmp(l)
      enddo
      H0(1:n_tmp,1:n_tmp) = H0_tmp(1:n_tmp,1:n_tmp)

      deallocate( numalpha_tmp, numbeta_tmp, coeff_tmp, cantor_pairing_tmp, H0_tmp )

      coeff_psi        (n_selected) = coeff_new
      numalpha_selected(n_selected) = na_new
      numbeta_selected (n_selected) = nb_new
      cantor_pairing   (n_selected) = cantor_pairing_new

      H0(1:n_tmp,n_selected) = H0_1d(1:n_tmp)
      H0(n_selected,1:n_tmp) = H0_1d(1:n_tmp)
      deallocate( H0_1d )
      H0(n_selected,n_selected) = Hdiag(ind_new) 

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



      ! ---------------------------------------------------------------------------------------
      !                                      new energy

      allocate( eigvec0(n_selected,n_selected), eigval0(n_selected) )
      call lapack_diag(eigval0, eigvec0, H0, n_selected, n_selected)

      ! get the postsvd ground state
      allocate( check_ov(n_selected) )
      do l = 1, n_selected
        overlop = 0.d0   
        do i = 1, n_selected
          ia = numalpha_selected(i)
          ib = numbeta_selected (i)
          if( ia .eq. ib ) overlop = overlop + eigvec0(i,l) * Dref(ia)
          !overlop = overlop + eigvec0(i,l) * coeff_psi(i)
        enddo
        check_ov(l) = dabs(overlop)
      enddo
      ind_gs       = MAXLOC( check_ov, DIM=1 )
      overlop      = check_ov(ind_gs)
      E0           = eigval0(ind_gs)+nuclear_repulsion
      coeff_psi(:) = eigvec0(:,ind_gs)

      deallocate( check_ov, eigval0, eigvec0 )

      print*, ' space dimen = ', n_selected
      print*, ' diag energy = ', E0
      print*, ' overlop     = ', overlop
      print*, ' index       = ', ind_gs
 
      ! ---------------------------------------------------------------------------------------


  write(2110, '( 3(I5,3X), 3(F15.8,3X) )') n_selected, na_new, nb_new, ept2, E0, overlop


      ! ---------------------------------------------------------------------------------------
      !                      remove selected pair | na_new   nb_new >

      allocate( numalpha_tmp(n_toselect), numbeta_tmp(n_toselect), Hdiag_tmp(n_toselect) )
      numalpha_tmp(:) = numalpha_toselect(:)
      numbeta_tmp (:) = numbeta_toselect (:)
      Hdiag_tmp   (:) = Hdiag            (:)

      ii = n_selected - 1
      allocate( Hkl_tmp(n_toselect,ii) )
      Hkl_tmp(1:n_toselect,1:ii) = Hkl_save(1:n_toselect,1:ii)

      deallocate( numalpha_toselect , numbeta_toselect, Hkl_save, Hdiag )

      n_tmp      = n_toselect
      n_toselect = n_toselect - 1
      print*, ' rank to select  = ', n_toselect

      allocate(numalpha_toselect(n_toselect), numbeta_toselect(n_toselect), Hkl_save(n_toselect,n_selected))
      allocate(Hdiag(n_toselect))

      Hkl_save = 0.d0
      l = 0
      do k = 1, n_tmp

        ia = numalpha_tmp(k) 
        ib = numbeta_tmp (k)
        cantor_pairing_ij = int( 0.5*(ia+ib)*(ia+ib+1) ) + ib
        if( ANY(cantor_pairing .eq. cantor_pairing_ij) ) cycle

        l = l + 1
        numalpha_toselect(l) = numalpha_tmp(k)
        numbeta_toselect (l) = numbeta_tmp (k)
        Hdiag            (l) = Hdiag_tmp   (k)

        Hkl_save(l,1:ii)     = Hkl_tmp(k,1:ii)

      enddo
      if( l .ne. n_toselect) then
        print *, " error in updating to select vectors"
        print *, " l          = ", l
        print *, " n_toselect = ", n_toselect
        stop
      endif

      deallocate( numalpha_tmp , numbeta_tmp , Hkl_tmp, Hdiag_tmp )

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


    endif 

    call SYSTEM_CLOCK(COUNT=W_tend_it, COUNT_RATE=W_ir)
    W_tot_time_it = real(W_tend_it-W_tbeg_it, kind=8) / real(W_ir, kind=8)
    print*, " "
    print*, "  elapsed time (min) = ", W_tot_time_it/60.d0
 
  end do
  !________________________________________________________________________________________________________
  !________________________________________________________________________________________________________



  ! ***************************************************************************************************
  ! save to ezfion
  !allocate( Uezfio(n_det_alpha_unique,rank0,1), Dezfio(rank0,1), Vezfio(n_det_beta_unique,rank0,1) )
  !do l = 1, rank0
  !  Dezfio(l,1)   = coeff_psi(l)
  !  Uezfio(:,l,1) = U0(:,l)
  !  Vezfio(:,l,1) = V0(:,l)
  !enddo
  !call ezfio_set_spindeterminants_n_det(N_det)
  !call ezfio_set_spindeterminants_n_states(N_states)
  !call ezfio_set_spindeterminants_n_det_alpha(n_det_alpha_unique)
  !call ezfio_set_spindeterminants_n_det_beta(n_det_beta_unique)
  !call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows)
  !call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns)
  !call ezfio_set_spindeterminants_psi_coef_matrix_values(psi_bilinear_matrix_values)

  !call ezfio_set_spindeterminants_n_svd_coefs(rank0)
  !call ezfio_set_spindeterminants_psi_svd_alpha(Uezfio)
  !call ezfio_set_spindeterminants_psi_svd_beta(Vezfio )
  !call ezfio_set_spindeterminants_psi_svd_coefs(Dezfio)
  !deallocate( Uezfio, Dezfio, Vezfio )
  ! ***************************************************************************************************

  call SYSTEM_CLOCK(COUNT=W_tend, COUNT_RATE=W_ir)
  W_tot_time   = real(W_tend - W_tbeg, kind=8) / real(W_ir, kind=8)
  print *, ' ___________________________________________________________________'
  print *, ' '
  print *, " Execution avec ", nb_taches, " threads"
  print *, " total elapsed time (min) = ", W_tot_time/60.d0
  print *, ' ___________________________________________________________________'



  deallocate( Dref )
  deallocate( Uref, Vref )
  
  deallocate( psi_coef )
  deallocate( numalpha_selected, numbeta_selected, numalpha_toselect, numbeta_toselect )
  deallocate( H0, Hdiag, Hkl_save )

end




subroutine const_H0_1d(na_new, nb_new, na_max, nb_max, n_selected, Uref, Vref, numalpha_selected, numbeta_selected, H0_1d)

  implicit none

  integer, intent(in)           :: na_new, nb_new, na_max, nb_max, n_selected
  integer, intent(in)           :: numalpha_selected(n_selected), numbeta_selected(n_selected)
  double precision, intent(in)  :: Uref(n_det_alpha_unique,n_det_alpha_unique)
  double precision, intent(in)  :: Vref(n_det_beta_unique ,n_det_beta_unique)
  double precision, intent(out) :: H0_1d(n_selected)

  integer(bit_kind)             :: det1(N_int,2)
  integer(bit_kind)             :: det2(N_int,2)
  integer                       :: degree, na, nb

  integer                       :: i, j, k, l, ii, jj, m
  double precision              :: h12

  double precision, allocatable :: Hmat_kl(:,:), tmp1(:,:), tmp2(:,:)
  double precision, allocatable :: U1d(:), V1d(:) 
  double precision, allocatable :: Utmp(:,:), Vtmp(:,:)

  double precision              :: ti, tf

  print *, ""
  print *, " start const_H0_1d"
  call wall_time(ti)

  na = n_det_alpha_unique
  nb = n_det_beta_unique

  allocate( U1d(na) , V1d(nb) )
  U1d(1:na) = Uref(1:na,na_new)
  V1d(1:nb) = Vref(1:nb,nb_new)

  allocate( tmp1(na,nb) )
  tmp1 = 0.d0

  do l = 1, nb
    det2(:,2) = psi_det_beta_unique(:,l)
    do j = 1, nb
      det1(:,2) = psi_det_beta_unique(:,j)

      call get_excitation_degree_spin(det1(1,2),det2(1,2),degree,N_int)
      if(degree .gt. 2) cycle

      do k = 1, na
        det2(:,1) = psi_det_alpha_unique(:,k)
        do i = 1, na
          det1(:,1) = psi_det_alpha_unique(:,i)

          call get_excitation_degree(det1,det2,degree,N_int)
          if(degree .gt. 2) cycle

          call i_H_j(det1, det2, N_int, h12)
          if( h12 .eq. 0.d0) cycle

          tmp1(i,j) += h12 * U1d(k) * V1d(l)

        enddo
      enddo
    enddo
  enddo

  deallocate( U1d , V1d )

  ! tmp2(j,m) = sum_i tmp1(i,j) x Uref(i,m)
  allocate( Utmp(na,na_max) )
  Utmp(1:na,1:na_max) = Uref(1:na,1:na_max)

  allocate( tmp2(nb,na_max) )
  call DGEMM('T', 'N', nb, na_max, na, 1.d0,   &
       tmp1, size(tmp1,1), Utmp, size(Utmp,1), &
       0.d0, tmp2, size(tmp2,1)                )
  deallocate( tmp1 )
  deallocate( Utmp )

  ! Hmat_kl(m,n) = sum_j tmp2(j,m) x Vref(j,n)
  allocate( Vtmp(nb,nb_max) )
  Vtmp(1:nb,1:nb_max) = Vref(1:nb,1:nb_max)

  allocate( Hmat_kl(na_max,nb_max) )
  call DGEMM('T', 'N', na_max, nb_max, nb, 1.d0, &
       tmp2, size(tmp2,1), Vtmp, size(Vtmp,1),   &
       0.d0, Hmat_kl, size(Hmat_kl,1)            )
  deallocate( tmp2 )
  deallocate( Vtmp )

  do m = 1, n_selected
    ii       = numalpha_selected(m)
    jj       = numbeta_selected (m)
    H0_1d(m) = Hmat_kl(ii,jj)
  enddo
  deallocate( Hmat_kl )

  call wall_time(tf)
  print *, " end const_H0_1d after (min) ", (tf-ti)/60.
  print *, ""

  return
end subroutine const_H0_1d





subroutine const_Hkl_1d(na_new, nb_new, na_max, nb_max, n_toselect, Uref, Vref, numalpha_toselect, numbeta_toselect, Hkl_1d)

  implicit none

  integer, intent(in)           :: na_new, nb_new, na_max, nb_max, n_toselect
  integer, intent(in)           :: numalpha_toselect(n_toselect), numbeta_toselect(n_toselect)
  double precision, intent(in)  :: Uref(n_det_alpha_unique,n_det_alpha_unique)
  double precision, intent(in)  :: Vref(n_det_beta_unique ,n_det_beta_unique)
  double precision, intent(out) :: Hkl_1d(n_toselect)

  integer(bit_kind)             :: det1(N_int,2)
  integer(bit_kind)             :: det2(N_int,2)
  integer                       :: degree, na, nb

  integer                       :: i, j, k, l, ii, jj, m
  double precision              :: h12

  double precision, allocatable :: Hmat_kl(:,:), tmp1(:,:), tmp2(:,:)
  double precision, allocatable :: U1d(:), V1d(:) 
  double precision, allocatable :: Utmp(:,:), Vtmp(:,:)

  double precision              :: ti, tf

  print *, ""
  print *, " start const_Hkl_1d"
  call wall_time(ti)

  na = n_det_alpha_unique
  nb = n_det_beta_unique

  allocate( U1d(na) , V1d(nb) )
  U1d(1:na) = Uref(1:na,na_new)
  V1d(1:nb) = Vref(1:nb,nb_new)

  allocate( tmp1(na,nb) )
  tmp1 = 0.d0

  do l = 1, nb
    det2(:,2) = psi_det_beta_unique(:,l)
    do j = 1, nb
      det1(:,2) = psi_det_beta_unique(:,j)

      call get_excitation_degree_spin(det1(1,2),det2(1,2),degree,N_int)
      if(degree .gt. 2) cycle

      do k = 1, na
        det2(:,1) = psi_det_alpha_unique(:,k)
        do i = 1, na
          det1(:,1) = psi_det_alpha_unique(:,i)

          call get_excitation_degree(det1,det2,degree,N_int)
          if(degree .gt. 2) cycle

          call i_H_j(det1, det2, N_int, h12)
          if( h12 .eq. 0.d0) cycle

          tmp1(i,j) += h12 * U1d(k) * V1d(l)

        enddo
      enddo
    enddo
  enddo

  deallocate( U1d , V1d )

  ! tmp2(j,m) = sum_i tmp1(i,j) x Uref(i,m)
  allocate( Utmp(na,na_max) )
  Utmp(1:na,1:na_max) = Uref(1:na,1:na_max)

  allocate( tmp2(nb,na_max) )
  call DGEMM('T', 'N', nb, na_max, na, 1.d0,   &
       tmp1, size(tmp1,1), Utmp, size(Utmp,1), &
       0.d0, tmp2, size(tmp2,1)                )
  deallocate( tmp1 , Utmp )

  ! Hmat_kl(m,n) = sum_j tmp2(j,m) x Vref(j,n)
  allocate( Vtmp(nb,nb_max) )
  Vtmp(1:nb,1:nb_max) = Vref(1:nb,1:nb_max)

  allocate( Hmat_kl(na_max,nb_max) )
  call DGEMM('T', 'N', na_max, nb_max, nb, 1.d0, &
       tmp2, size(tmp2,1), Vtmp, size(Vtmp,1),   &
       0.d0, Hmat_kl, size(Hmat_kl,1)            )
  deallocate( tmp2 )
  deallocate( Vtmp )

  do m = 1, n_toselect
    ii        = numalpha_toselect(m)
    jj        = numbeta_toselect (m)
    Hkl_1d(m) = Hmat_kl(ii,jj)
  enddo
  deallocate( Hmat_kl )

  call wall_time(tf)
  print *, " end const_Hkl_1d after (min) ", (tf-ti)/60.
  print *, ""

  return
end subroutine const_Hkl_1d