program delta_svd_v2

  BEGIN_DOC
  ! !!!
  END_DOC

  implicit none

  read_wf = .true.
  TOUCH read_wf

  call run()

end 



subroutine run()

  implicit none

  integer                       :: ii, ii_i, ii_j, k, l
  integer                       :: na, nb, nsvd
  double precision, allocatable :: U(:,:), V(:,:), S(:)
  double precision, allocatable :: deltaDet_ex(:), deltaDet_ap(:), deltaSVD_kl(:,:)

  na = n_det_alpha_unique
  nb = n_det_beta_unique

  ! read SVD vectors
  call ezfio_get_spindeterminants_n_svd_coefs(nsvd)
  allocate( U(na,nsvd), S(nsvd), V(nb,nsvd) )
  call ezfio_get_spindeterminants_psi_svd_alpha(U)
  call ezfio_get_spindeterminants_psi_svd_beta (V)
  call ezfio_get_spindeterminants_psi_svd_coefs(S)

  !allocate( deltaSVD_kk(nsvd) )
  !call obtenir_deltaSVD_kk(nsvd, U, V, deltaSVD_kk)
  allocate( deltaSVD_kl(nsvd,nsvd) )
  call obtenir_deltaSVD_kl(nsvd, U, V, deltaSVD_kl)

  allocate( deltaDet_ex(N_det) , deltaDet_ap(N_det) )
  do ii = 1, N_det
    ii_i = psi_bilinear_matrix_rows   (ii)
    ii_j = psi_bilinear_matrix_columns(ii)

    deltaDet_ex(ii) = 0.d0
    deltaDet_ap(ii) = 0.d0
    do k = 1, nsvd
      deltaDet_ap(ii) = deltaDet_ap(ii) + U(ii_i,k) * V(ii_j,k) * deltaSVD_kl(k,k)
      do l = 1, nsvd
        deltaDet_ex(ii) = deltaDet_ex(ii) + U(ii_i,k) * V(ii_j,l) * deltaSVD_kl(k,l)
      enddo
    enddo

    print *, deltaDet_ex(ii) , deltaDet_ap(ii)
  enddo

  !call ezfio_set_dmc_dress_dmc_delta_h(delta_Det_ex)

  deallocate( U , S , V )
  deallocate( deltaSVD_kl )
  deallocate( deltaDet_ex , deltaDet_ap )

end





! _______________________________________________________________________________________________________
!
!     [ delta_SVD ]_kl = < dup_SVD ddn_SVD | H_tilde - H | psi >
!                      = \sum_{i,j}  U_{i,k} V_{j,l} < dup_i ddn_j | H_tilde - H | psi >
!
subroutine obtenir_deltaSVD_kl(nsvd, U, V, deltaSVD_kl)

  implicit none

  integer,          intent(in)  :: nsvd
  double precision, intent(in)  :: U(n_det_alpha_unique,nsvd), V(n_det_beta_unique,nsvd)
  double precision, intent(out) :: deltaSVD_kl(nsvd,nsvd)

  integer                       :: i_state
  integer                       :: k, l, ii, jj, ii_i, ii_j, jj_i, jj_j

  integer(bit_kind)             :: det_ii(N_int,2), det_jj(N_int,2)
  double precision              :: hmono, heff, hderiv, hthree, htilde_ij, hij, delta_mat

  double precision, allocatable :: delta_det(:)

 PROVIDE scalar_mu_r_pot_physicist_mo deriv_mu_r_pot_physicist_mo
 PROVIDE three_body_3_index three_body_3_index_exch_12 three_body_3_index_exch_13 three_body_3_index_exch_23
 PROVIDE three_body_5_index three_body_5_index_exch_13 three_body_5_index_exch_32
 PROVIDE three_body_4_index three_body_4_index_exch_12 three_body_4_index_exch_12_part


  i_state = 1

  ! ! !
  det_ii(:,1) = psi_det_alpha_unique(:,1)
  det_ii(:,2) = psi_det_beta_unique (:,1)
  det_jj(:,1) = psi_det_alpha_unique(:,1)
  det_jj(:,2) = psi_det_beta_unique (:,1)
  call htilde_mat(det_ii, det_jj, hmono, heff, hderiv, hthree, htilde_ij)
  call i_H_j(det_ii, det_jj, N_int, hij)
  ! ! !


  ! -------------------------------------------------------------------------
  ! < dup_i ddn_j | H_tilde - H | psi >
  allocate( delta_det(N_det) )
  delta_det(:) = 0.d0

 !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8)                   &
 !$OMP SHARED(N_det, i_state, N_int, psi_bilinear_matrix_rows          &
 !$OMP      , psi_bilinear_matrix_columns, psi_bilinear_matrix_values  &
 !$OMP      , psi_det_alpha_unique, psi_det_beta_unique, delta_det)    &
 !$OMP PRIVATE(ii, ii_i, ii_j, jj, jj_i, jj_j, det_ii, det_jj          &
 !$OMP       , hmono, heff, hderiv, hthree, htilde_ij, hij)
  do ii = 1, N_det
    ii_i = psi_bilinear_matrix_rows   (ii)
    ii_j = psi_bilinear_matrix_columns(ii)

    det_ii(:,1) = psi_det_alpha_unique(:,ii_i)
    det_ii(:,2) = psi_det_beta_unique (:,ii_j)

    do jj = 1, N_det
      jj_i = psi_bilinear_matrix_rows   (jj)
      jj_j = psi_bilinear_matrix_columns(jj)

      det_jj(:,1) = psi_det_alpha_unique(:,jj_i)
      det_jj(:,2) = psi_det_beta_unique (:,jj_j)

      call htilde_mat(det_ii, det_jj, hmono, heff, hderiv, hthree, htilde_ij)
      call i_H_j(det_ii, det_jj, N_int, hij)

      delta_det(ii) = delta_det(ii) + (htilde_ij-hij) * psi_bilinear_matrix_values(jj,i_state) 
    enddo
  enddo
 !$OMP END PARALLEL DO
  ! -------------------------------------------------------------------------


  deltaSVD_kl(:,:) = 0.d0
 !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8)                   &
 !$OMP SHARED(nsvd, N_det, U, V, delta_det, deltaSVD_kl                &
 !$OMP      , psi_bilinear_matrix_rows, psi_bilinear_matrix_columns)   &
 !$OMP PRIVATE(k, l, ii, ii_i, ii_j)
  do l = 1, nsvd
    do k = 1, nsvd
      do ii = 1, N_det
        ii_i = psi_bilinear_matrix_rows   (ii)
        ii_j = psi_bilinear_matrix_columns(ii)

        deltaSVD_kl(k,l) = deltaSVD_kl(k,l) + U(ii_i,k) * V(ii_j,l) * delta_det(ii)
      enddo
    enddo
  enddo
 !$OMP END PARALLEL DO

  deallocate( delta_det )

  return
end subroutine obtenir_deltaSVD_kl
! _______________________________________________________________________________________________________
! _______________________________________________________________________________________________________