program export_trexio
  use trexio
  implicit none
  BEGIN_DOC
  !     Exports the wave function in TREXIO format
  END_DOC

  integer(8)                     :: f ! TREXIO file handle
  integer                        :: rc

  print *, 'TREXIO file : '//trim(trexio_filename)
  print *, ''

  if (trexio_backend == 0) then
    f = trexio_open(trexio_filename, 'w', TREXIO_HDF5)
  else if (trexio_backend == 1) then
    f = trexio_open(trexio_filename, 'w', TREXIO_TEXT)
  endif
  if (f == 0) then
    print *, 'Unable to open TREXIO file for writing'
    stop -1
  endif

! ------------------------------------------------------------------------------

! Nuclei
! ------

  rc = trexio_write_nucleus_num(f, nucl_num)
  call check_success(rc)

  rc = trexio_write_nucleus_charge(f, nucl_charge)
  call check_success(rc)

  rc = trexio_write_nucleus_coord(f, nucl_coord_transp)
  call check_success(rc)

  ! Electrons
  rc = trexio_write_electron_up_num(f, elec_alpha_num)
  call check_success(rc)
!
  rc = trexio_write_electron_dn_num(f, elec_beta_num)
  call check_success(rc)


! Basis
! -----

!  rc = trexio_write_basis_type(f, 'Gaussian')
!  call check_success(rc)

  rc = trexio_write_basis_shell_num(f, sum(Nucl_num_shell_Aos))
  call check_success(rc)


! Atomic orbitals
! ---------------

  rc = trexio_write_ao_num(f, ao_num)
  call check_success(rc)

  rc = trexio_write_ao_cartesian(f, 1)
  call check_success(rc)


! One-e AO integrals
! ------------------

  rc = trexio_write_ao_1e_int_overlap(f,ao_overlap)
  call check_success(rc)

  rc = trexio_write_ao_1e_int_kinetic(f,ao_kinetic_integrals)
  call check_success(rc)

  rc = trexio_write_ao_1e_int_potential_n_e(f,ao_integrals_n_e)
  call check_success(rc)

  if (do_pseudo) then
    rc = trexio_write_ao_1e_int_ecp_local(f,ao_pseudo_integrals_local)
    call check_success(rc)

    rc = trexio_write_ao_1e_int_ecp_non_local(f,ao_pseudo_integrals_non_local)
    call check_success(rc)
  endif

  rc = trexio_write_ao_1e_int_core_hamiltonian(f,ao_one_e_integrals)
  call check_success(rc)


! Molecular orbitals
! ------------------

!  rc = trexio_write_mo_type(f, mo_label)
!  call check_success(rc)

  rc = trexio_write_mo_num(f, mo_num)
  call check_success(rc)

  rc = trexio_write_mo_coef(f, mo_coef)
  call check_success(rc)


! One-e MO integrals
! ------------------

  rc = trexio_write_mo_1e_int_kinetic(f,mo_kinetic_integrals)
  call check_success(rc)

  rc = trexio_write_mo_1e_int_potential_n_e(f,mo_integrals_n_e)
  call check_success(rc)

  if (do_pseudo) then
    rc = trexio_write_mo_1e_int_ecp_local(f,mo_pseudo_integrals_local)
    call check_success(rc)

    rc = trexio_write_mo_1e_int_ecp_non_local(f,mo_pseudo_integrals_non_local)
    call check_success(rc)
  endif

  rc = trexio_write_mo_1e_int_core_hamiltonian(f,one_e_dm_mo)
  call check_success(rc)


! RDM
! ----

!  rc = trexio_write_rdm_one_e(f,one_e_dm_mo)
!  call check_success(rc)
!
!  rc = trexio_write_rdm_one_e_up(f,one_e_dm_mo_alpha_average)
!  call check_success(rc)
!
!  rc = trexio_write_rdm_one_e_dn(f,one_e_dm_mo_beta_average)
!  call check_success(rc)


! ------------------------------------------------------------------------------

  rc = trexio_close(f)
  call check_success(rc)

end


subroutine check_success(rc)
  use trexio
  implicit none
  integer, intent(in)            :: rc
  character*(128)                :: str
  if (rc /= TREXIO_SUCCESS) then
    call trexio_string_of_error(rc,str)
    print *, 'TREXIO Error: ' //trim(str)
    stop -1
  endif
end

! -*- mode: f90 -*-