program qmcpack
  implicit none
  BEGIN_DOC
! Generates a file for CHAMP
  END_DOC

  integer :: i,j,k,l, istate

  read_wf = .True.
  TOUCH read_wf
  do j=1,ao_prim_num_max
    do i=1,ao_num
      ao_coef(i,j) = ao_coef(i,j) * ao_coef_normalization_factor(i)
    enddo
  enddo
  call ezfio_set_ao_basis_ao_coef(ao_coef)
  do j=1,mo_num
    do i=1,ao_num
      mo_coef(i,j) *= 1.d0/ao_coef_normalization_factor(i)
    enddo
  enddo
  call save_mos
  call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5')
  call system('$QP_ROOT/src/champ/qp_convert.py '//trim(ezfio_filename))

!  integer :: iunit
!  integer, external :: getUnitAndOpen
!  iunit = getUnitAndOpen(trim(ezfio_filename)//'.H','w')
!
!  double precision, external :: diag_h_mat_elem
!  write(iunit,*) N_states
!  do istate=1,N_states
!    write(iunit,*) istate, psi_energy_with_nucl_rep(istate)
!  enddo
!  write(iunit,*) N_det
!  do k=1,N_det
!    write(iunit,'(I10,X,F22.15)') k, diag_h_mat_elem(psi_det(1,1,k),N_int) + nuclear_repulsion
!  enddo
!
!
!  double precision :: F(N_states)
!  integer(bit_kind), allocatable :: det(:,:,:)
!  double precision , allocatable :: coef(:,:)
!  integer :: ispin
!  double precision :: norm(N_states), hij
!  allocate(det(N_int,2,N_det), coef(N_det,N_states))
!  do j=1,mo_num
!    do i=1,j-1
!      do ispin=1,2
!        call build_singly_excited_wavefunction(j,i,1,det,coef)
!        F = 0.d0
!        do istate=1,N_states
!          norm(istate) = 0.d0
!          do k=1,N_det
!            norm(istate) = norm(istate) + coef(k,istate) * coef(k,istate)
!          enddo
!          if (norm(istate) > 0.d0) then
!            norm(istate) = (1.d0/dsqrt(norm(istate)))
!          endif
!        enddo
!        if (sum(norm(:)) > 0.d0) then
!          
!          do istate = 1,N_states
!            coef(:,istate) = coef(:,istate) * norm(istate)
!          enddo
!          !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,istate,hij) REDUCTION(+:F)
!          do k=1,N_det
!            if (sum(coef(k,:)*coef(k,:))==0.d0) cycle
!            call i_H_j(det(1,1,k), det(1,1,k), N_int, hij)
!            do istate=1,N_states
!              F(istate) = F(istate) + hij*coef(k,istate)*coef(k,istate)
!            enddo
!            do l=1,k-1
!              if (sum(coef(l,:)*coef(l,:))==0.d0) cycle
!              call i_H_j(det(1,1,k), det(1,1,l), N_int, hij)
!              do istate=1,N_states
!                F(istate) = F(istate) + 2.d0*hij*coef(k,istate)*coef(l,istate)
!              enddo
!            enddo
!          enddo
!          !$OMP END PARALLEL DO
!          F(:) = F(:) - psi_energy(:)
!        endif
!        do istate=1,N_states
!          write(iunit,'(I4,X,I4,X,I1,X,I3,X,F22.15)') i, j, ispin, istate, F(istate) 
!        enddo
!      enddo
!    enddo
!  enddo
!
!  deallocate(det,coef)
!
!  close(iunit)

end