program export_integrals_mo
  call run
end

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

  iunit = getunitandopen(name,'w')
  do j=1,mo_num
    do i=1,mo_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 MO 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:
!
! T_mo.qp : kinetic energy integrals
!
! P_mo.qp : pseudopotential integrals
!
! V_mo.qp : electron-nucleus potential
!
! W_mo.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_mo_two_e_integral

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

  call write_2d('T_mo.qp',mo_kinetic_integrals)
  call write_2d('P_mo.qp',mo_pseudo_integrals)
  call write_2d('V_mo.qp',mo_integrals_n_e)

  iunit = getunitandopen('W_mo.qp','w')
  do l=1,mo_num
    do k=1,mo_num
      do j=1,mo_num
        do i=1,mo_num
          integral = get_mo_two_e_integral(i,j,k,l,mo_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