program export_integrals_ao
  call run
end

subroutine write_2d(name,A)
  implicit none
  character*(*), intent(in)      :: name
  double precision, intent(in)   :: A(ao_num,ao_num)
  integer                        :: i, j
  integer                        :: iunit
  integer                        :: getunitandopen

  iunit = getunitandopen(name,'w')
  do j=1,ao_num
    do i=1,ao_num
      if (A(i,j) /= 0.d0) then
        write (iunit,*) i,j, A(i,j)
      endif
    enddo
  enddo
  close(iunit)
end

subroutine run
  use map_module
  implicit none
  BEGIN_DOC
! Program to import integrals in the AO basis.
!
! one-electron integrals : format is : i j value
! two-electron integrals : format is : i j k l value
!  Dirac's notation is used : <ij|kl> is <r1 r2|r1 r2>
!
! These files are read:
!
! E.qp : Contains the nuclear repulsion energy
!
! T.qp : kinetic energy integrals
!
! S.qp : overlap matrix
!
! P.qp : pseudopotential integrals
!
! V.qp : electron-nucleus potential
!
! W.qp : electron repulsion integrals
!
  END_DOC

  integer                        :: iunit
  integer                        :: getunitandopen

  integer                        :: i,j,k,l
  double precision               :: integral
  double precision, allocatable  :: A(:,:)

  integer                        :: n_integrals
  integer(key_kind)              :: idx

  double precision, external     :: get_ao_two_e_integral

  allocate (A(ao_num,ao_num))
  call ezfio_set_nuclei_nuclear_repulsion(nuclear_repulsion)

  call write_2d('T.qp',ao_kinetic_integrals)
  call write_2d('S.qp',ao_overlap)
  call write_2d('P.qp',ao_pseudo_integrals)
  call write_2d('V.qp',ao_integrals_n_e)

  iunit = getunitandopen('E.qp','w')
  write(iunit,*) nuclear_repulsion
  close(iunit)

  iunit = getunitandopen('W.qp','w')
  do l=1,ao_num
    do k=1,ao_num
      do j=1,ao_num
        do i=1,ao_num
          integral = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
          if (integral /= 0.d0) then
            write (iunit,'(4(I5,2X),E22.15)') i,j,k,l, integral
          endif
        enddo
      enddo
    enddo
  enddo
  close(iunit)

end