BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral_dressing, (ao_num_align,ao_num) ]
  implicit none
  BEGIN_DOC
  ! Dressing of the core hamiltonian in the orthogonal AO basis set
  END_DOC

  integer                        :: i,j,k
  integer                        :: mu, nu, lambda, A
  double precision               :: tmp

  ao_ortho_mono_elec_integral_dressing = 0.d0
  i = idx_dressing
  do mu=1,ao_num
    if (dabs(mo_coef_in_ao_ortho_basis(mu,i)) > 1.d-5) then
      do A=1,nucl_num
        tmp = 0.d0
        do nu=1,ao_num
          tmp += AO_orthoSlaOverlap_matrix(nu,A) * ao_ortho_mono_elec_integral(mu,nu)
        enddo
        ao_ortho_mono_elec_integral_dressing(mu,mu) += cusp_C(A,i) * (AO_orthoSlaH_matrix(mu,A) - tmp)
      enddo
      ao_ortho_mono_elec_integral_dressing(mu,mu) *= 1.d0/mo_coef_in_ao_ortho_basis(mu,i)
    endif
  enddo
END_PROVIDER


BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral, (ao_num_align, ao_num) ]
 implicit none
 BEGIN_DOC
 ! h core in orthogonal AO basis
 END_DOC
  double precision, allocatable ::  T(:,:)
  allocate(T(ao_num,ao_num))
  call dgemm('T','N',ao_num,ao_num,ao_num,1.d0,                    &
      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),      &
      ao_mono_elec_integral, size(ao_mono_elec_integral,1),          &
      0.d0, T, size(T,1))
  call dgemm('N','N',ao_num,ao_num,ao_num,1.d0,                    &
      T, size(T,1),                                                  &
      ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1),      &
      0.d0, ao_ortho_mono_elec_integral, size(ao_ortho_mono_elec_integral,1))
  deallocate(T)
END_PROVIDER


BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_dressing, (ao_num,ao_num) ]
  implicit none
  BEGIN_DOC
  ! Dressing of the core hamiltonian in the AO basis set
  END_DOC
  call ao_ortho_cano_to_ao(ao_ortho_mono_elec_integral_dressing,size(ao_ortho_mono_elec_integral_dressing,1),&
    ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1))
END_PROVIDER

BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_dressing, (mo_tot_num_align,mo_tot_num) ]
  implicit none
  BEGIN_DOC
  ! Dressing of the core hamiltonian in the MO basis set
  END_DOC

  call ao_to_mo(ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1),&
    mo_mono_elec_integral_dressing,size(mo_mono_elec_integral_dressing,1))
END_PROVIDER


BEGIN_PROVIDER [ integer, idx_dressing ]
 implicit none
 BEGIN_DOC
 ! Index of the MO which is being dressed
 END_DOC
 idx_dressing = 1
END_PROVIDER


BEGIN_PROVIDER [ double precision, cusp_corrected_mos, (ao_num_align,mo_tot_num) ]
  implicit none
  BEGIN_DOC
  ! Dressing core hamiltonian in the AO basis set
  END_DOC
  integer :: i,j
  double precision, allocatable :: F(:,:), M(:,:)
  allocate(F(mo_tot_num_align,mo_tot_num),M(ao_num,mo_tot_num))

  logical :: oneshot

!  oneshot = .True. 
  oneshot = .False.

  if (oneshot) then
    cusp_corrected_mos(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num)
    slater_coef(1:nucl_num,1:mo_tot_num) = cusp_C(1:nucl_num,1:mo_tot_num)
    return
  endif

  do idx_dressing=1,mo_tot_num

    if (idx_dressing>1) then
      TOUCH idx_dressing
    endif

    do j=1,mo_tot_num
      do i=1,mo_tot_num
        F(i,j) = Fock_matrix_mo(i,j)
      enddo
    enddo

    do j=1,mo_tot_num
      do i=1,ao_num
        M(i,j) = mo_coef(i,j)
      enddo
    enddo

    integer :: it
    do it=1,128

!      print *,  'C', mo_coef(1:ao_num,1:mo_tot_num)
!      print *,  'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num)
!      print *,  'cAi', cusp_C(1:nucl_num,1:mo_tot_num)
!      print *,  'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
!      print *,  'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num)
!      print *,  'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num)
!      print *,  'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
!      print *,  'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
!      print *,  'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num)
!      print *,  'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
!      print *,  'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num)
      do j=1,mo_tot_num
        do i=1,mo_tot_num
          Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j)
        enddo
      enddo
      do i=1,mo_tot_num
        Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i)
      enddo
      double precision :: conv
      conv = 0.d0
      do j=1,mo_tot_num
        do i=1,mo_tot_num
          if (i==j) cycle
          conv = max(conv,Fock_matrix_mo(i,j))
        enddo
      enddo
      TOUCH Fock_matrix_mo Fock_matrix_diag_mo

      mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num) 
      TOUCH mo_coef
!      print *,  'C', mo_coef(1:ao_num,1:mo_tot_num)
!      print *,  '-----'
      print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv)
      if (conv < 1.d-5) exit

    enddo
    cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing)
    slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing)
  enddo

  idx_dressing = 1
  mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num)
  soft_TOUCH mo_coef idx_dressing slater_coef

END_PROVIDER