diff --git a/src/trexio/EZFIO.cfg b/src/trexio/EZFIO.cfg
new file mode 100644
index 00000000..8606e908
--- /dev/null
+++ b/src/trexio/EZFIO.cfg
@@ -0,0 +1,54 @@
+[backend]
+type: integer
+doc: Back-end used in TREXIO. 0: HDF5, 1:Text
+interface: ezfio, ocaml, provider
+default: 0
+
+[trexio_file]
+type: character*(256)
+doc: Name of the exported TREXIO file
+interface: ezfio, ocaml, provider
+default: None
+
+[export_rdm]
+type: logical
+doc: If True, export two-body reduced density matrix
+interface: ezfio, ocaml, provider
+default: False
+
+[export_ao_one_e_ints]
+type: logical
+doc: If True, export one-electron integrals in AO basis
+interface: ezfio, ocaml, provider
+default: False
+
+[export_mo_one_e_ints]
+type: logical
+doc: If True, export one-electron integrals in MO basis
+interface: ezfio, ocaml, provider
+default: False
+
+[export_ao_two_e_ints]
+type: logical
+doc: If True, export two-electron integrals in AO basis
+interface: ezfio, ocaml, provider
+default: False
+
+[export_ao_two_e_ints_cholesky]
+type: logical
+doc: If True, export Cholesky-decomposed two-electron integrals in AO basis
+interface: ezfio, ocaml, provider
+default: False
+
+[export_mo_two_e_ints]
+type: logical
+doc: If True, export two-electron integrals in MO basis
+interface: ezfio, ocaml, provider
+default: False
+
+[export_mo_two_e_ints_cholesky]
+type: logical
+doc: If True, export Cholesky-decomposed two-electron integrals in MO basis
+interface: ezfio, ocaml, provider
+default: False
+
diff --git a/src/trexio/README.rst b/src/trexio/README.rst
new file mode 100644
index 00000000..7d7304c6
--- /dev/null
+++ b/src/trexio/README.rst
@@ -0,0 +1,6 @@
+======
+trexio
+======
+
+Module for handling TREXIO files.
+See https://github.com/trex-coe/trexio
diff --git a/src/trexio/export_trexio.irp.f b/src/trexio/export_trexio.irp.f
new file mode 100644
index 00000000..3ae0dcb4
--- /dev/null
+++ b/src/trexio/export_trexio.irp.f
@@ -0,0 +1,7 @@
+program export_trexio_prog
+  implicit none
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call export_trexio
+end
+
diff --git a/src/trexio/export_trexio_routines.irp.f b/src/trexio/export_trexio_routines.irp.f
new file mode 100644
index 00000000..d69e7a70
--- /dev/null
+++ b/src/trexio/export_trexio_routines.irp.f
@@ -0,0 +1,604 @@
+subroutine export_trexio
+  use trexio
+  implicit none
+  BEGIN_DOC
+  !     Exports the wave function in TREXIO format
+  END_DOC
+
+  integer(trexio_t)              :: f ! TREXIO file handle
+  integer(trexio_exit_code)      :: rc
+  double precision, allocatable  :: factor(:)
+
+  print *, 'TREXIO file : '//trim(trexio_filename)
+  print *, ''
+
+  call system('cp '//trim(trexio_filename)//' '//trim(trexio_filename)//'.bak')
+  if (backend == 0) then
+    f = trexio_open(trexio_filename, 'u', TREXIO_HDF5, rc)
+  else if (backend == 1) then
+    f = trexio_open(trexio_filename, 'u', TREXIO_TEXT, rc)
+  endif
+  if (f == 0_8) then
+    print *, 'Unable to open TREXIO file for writing'
+    print *, 'rc = ', rc
+    stop -1
+  endif
+  call ezfio_set_trexio_trexio_file(trexio_filename)
+
+! ------------------------------------------------------------------------------
+
+! Electrons
+! ---------
+
+  print *, 'Electrons'
+
+  rc = trexio_write_electron_up_num(f, elec_alpha_num)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_electron_dn_num(f, elec_beta_num)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+
+! Nuclei
+! ------
+
+  print *, 'Nuclei'
+
+  rc = trexio_write_nucleus_num(f, nucl_num)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_nucleus_charge(f, nucl_charge)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_nucleus_coord(f, nucl_coord_transp)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_nucleus_label(f, nucl_label, 32)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_nucleus_repulsion(f, nuclear_repulsion)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+
+! Pseudo-potentials
+! -----------------
+
+  if (do_pseudo) then
+
+    print *, 'ECP'
+    integer                       :: num
+
+    num = 0
+    do k=1,pseudo_klocmax
+      do i=1,nucl_num
+        if (pseudo_dz_k(i,k) /= 0.d0) then
+          num = num+1
+        end if
+      end do
+    end do
+
+    do l=0,pseudo_lmax
+      do k=1,pseudo_kmax
+        do i=1,nucl_num
+          if (pseudo_dz_kl(i,k,l) /= 0.d0) then
+            num = num+1
+          end if
+        end do
+      end do
+    end do
+
+    integer, allocatable          :: ang_mom(:), nucleus_index(:), power(:), lmax(:)
+    double precision, allocatable :: exponent(:), coefficient(:)
+
+    allocate(ang_mom(num), nucleus_index(num), exponent(num), coefficient(num), power(num), &
+            lmax(nucl_num) )
+
+    do i=1,nucl_num
+      lmax(i) = -1
+      do l=0,pseudo_lmax
+        do k=1,pseudo_kmax
+          if (pseudo_dz_kl_transp(k,l,i) /= 0.d0) then
+            lmax(i) = max(lmax(i), l)
+          end if
+        end do
+      end do
+    end do
+
+    j = 0
+    do i=1,nucl_num
+      do k=1,pseudo_klocmax
+        if (pseudo_dz_k_transp(k,i) /= 0.d0) then
+          j = j+1
+          ang_mom(j) = lmax(i)+1
+          nucleus_index(j) = i
+          exponent(j) = pseudo_dz_k_transp(k,i)
+          coefficient(j) = pseudo_v_k_transp(k,i)
+          power(j) = pseudo_n_k_transp(k,i)
+        end if
+      end do
+
+      do l=0,lmax(i)
+        do k=1,pseudo_kmax
+          if (pseudo_dz_kl_transp(k,l,i) /= 0.d0) then
+            j = j+1
+            ang_mom(j) = l
+            nucleus_index(j) = i
+            exponent(j) = pseudo_dz_kl_transp(k,l,i)
+            coefficient(j) = pseudo_v_kl_transp(k,l,i)
+            power(j) = pseudo_n_kl_transp(k,l,i)
+          end if
+        end do
+      end do
+    end do
+
+
+    lmax(:) = lmax(:)+1
+    rc = trexio_write_ecp_max_ang_mom_plus_1(f, lmax)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ecp_z_core(f, int(nucl_charge_remove))
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ecp_num(f, num)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ecp_ang_mom(f, ang_mom)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ecp_nucleus_index(f, nucleus_index)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ecp_exponent(f, exponent)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ecp_coefficient(f, coefficient)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ecp_power(f, power)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+  endif
+
+
+! Basis
+! -----
+
+  print *, 'Basis'
+
+
+  rc = trexio_write_basis_type(f, 'Gaussian', len('Gaussian'))
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_basis_prim_num(f, prim_num)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+   rc = trexio_write_basis_shell_num(f, shell_num)
+   call trexio_assert(rc, TREXIO_SUCCESS)
+
+   rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
+   call trexio_assert(rc, TREXIO_SUCCESS)
+
+   rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
+   call trexio_assert(rc, TREXIO_SUCCESS)
+
+   allocate(factor(shell_num))
+   if (ao_normalized) then
+     factor(1:shell_num) = shell_normalization_factor(1:shell_num)
+   else
+     factor(1:shell_num) = 1.d0
+   endif
+   rc = trexio_write_basis_shell_factor(f, factor)
+   call trexio_assert(rc, TREXIO_SUCCESS)
+
+   deallocate(factor)
+
+  rc = trexio_write_basis_shell_index(f, shell_index)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_basis_exponent(f, prim_expo)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_basis_coefficient(f, prim_coef)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  allocate(factor(prim_num))
+  if (primitives_normalized) then
+    factor(1:prim_num) = prim_normalization_factor(1:prim_num)
+  else
+    factor(1:prim_num) = 1.d0
+  endif
+  rc = trexio_write_basis_prim_factor(f, factor)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+  deallocate(factor)
+
+
+! Atomic orbitals
+! ---------------
+
+  print *, 'AOs'
+
+  rc = trexio_write_ao_num(f, ao_num)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_ao_cartesian(f, 1)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_ao_shell(f, ao_shell)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  integer :: i, pow0(3), powA(3), j, k, l, nz
+  double precision :: normA, norm0, C_A(3), overlap_x, overlap_z, overlap_y, c
+  nz=100
+
+  C_A(1) = 0.d0
+  C_A(2) = 0.d0
+  C_A(3) = 0.d0
+
+  allocate(factor(ao_num))
+  if (ao_normalized) then
+    do i=1,ao_num
+      l = ao_first_of_shell(ao_shell(i))
+      factor(i) = (ao_coef_normalized(i,1)+tiny(1.d0))/(ao_coef_normalized(l,1)+tiny(1.d0))
+    enddo
+  else
+    factor(:) = 1.d0
+  endif
+  rc = trexio_write_ao_normalization(f, factor)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+  deallocate(factor)
+
+! One-e AO integrals
+! ------------------
+
+  if (export_ao_one_e_ints) then
+    print *, 'AO one-e integrals'
+
+    rc = trexio_write_ao_1e_int_overlap(f,ao_overlap)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ao_1e_int_kinetic(f,ao_kinetic_integrals)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_ao_1e_int_potential_n_e(f,ao_integrals_n_e)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    if (do_pseudo) then
+      rc = trexio_write_ao_1e_int_ecp(f, ao_pseudo_integrals_local + ao_pseudo_integrals_non_local)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    endif
+
+    rc = trexio_write_ao_1e_int_core_hamiltonian(f,ao_one_e_integrals)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+  end if
+
+! Two-e AO integrals
+! ------------------
+
+  if (export_ao_two_e_ints) then
+    print *, 'AO two-e integrals'
+    PROVIDE ao_two_e_integrals_in_map
+
+    integer(8), parameter :: BUFSIZE=100000_8
+    double precision :: eri_buffer(BUFSIZE), integral
+    integer(4) :: eri_index(4,BUFSIZE)
+    integer(8) :: icount, offset
+
+    double precision, external :: get_ao_two_e_integral
+
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,ao_num
+      do k=1,ao_num
+        do j=l,ao_num
+          do i=k,ao_num
+            if (i==j .and. k<l) cycle
+            if (i<j) cycle
+            integral = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+              call trexio_assert(rc, TREXIO_SUCCESS)
+              offset += icount
+              icount = 0_8
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    end if
+  end if
+
+! Two-e AO integrals - Cholesky
+! -----------------------------
+
+  integer(4) :: chol_index(3,BUFSIZE)
+  double precision :: chol_buffer(BUFSIZE)
+
+  if (export_ao_two_e_ints_cholesky) then
+    print *, 'AO two-e integrals Cholesky'
+
+    rc = trexio_write_ao_2e_int_eri_cholesky_num(f, cholesky_ao_num)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    icount = 0_8
+    offset = 0_8
+    do k=1,cholesky_ao_num
+     do j=1,ao_num
+      do i=1,ao_num
+         integral = cholesky_ao(i,j,k)
+         if (integral == 0.d0) cycle
+         icount += 1_8
+         chol_buffer(icount) = integral
+         chol_index(1,icount) = i
+         chol_index(2,icount) = j
+         chol_index(3,icount) = k
+         if (icount == BUFSIZE) then
+           rc = trexio_write_ao_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
+           call trexio_assert(rc, TREXIO_SUCCESS)
+           offset += icount
+           icount = 0_8
+         end if
+      end do
+     end do
+    end do
+
+    if (icount > 0_8) then
+      rc = trexio_write_ao_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    end if
+  end if
+
+
+
+! Molecular orbitals
+! ------------------
+
+  print *, 'MOs'
+
+  rc = trexio_write_mo_type(f, mo_label, len(trim(mo_label)))
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_mo_num(f, mo_num)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_mo_coefficient(f, mo_coef)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  if ( (trim(mo_label) == 'Canonical').and. &
+       (export_mo_two_e_ints_cholesky.or.export_mo_two_e_ints) ) then
+    rc = trexio_write_mo_energy(f, fock_matrix_diag_mo)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+  endif
+
+  rc = trexio_write_mo_class(f, mo_class, len(mo_class(1)))
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+! One-e MO integrals
+! ------------------
+
+  if (export_mo_one_e_ints) then
+    print *, 'MO one-e integrals'
+
+    rc = trexio_write_mo_1e_int_kinetic(f,mo_kinetic_integrals)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_write_mo_1e_int_potential_n_e(f,mo_integrals_n_e)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    if (do_pseudo) then
+      rc = trexio_write_mo_1e_int_ecp(f,mo_pseudo_integrals_local)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    endif
+
+    rc = trexio_write_mo_1e_int_core_hamiltonian(f,mo_one_e_integrals)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+  end if
+
+! Two-e MO integrals
+! ------------------
+
+  if (export_mo_two_e_ints) then
+    print *, 'MO two-e integrals'
+    PROVIDE mo_two_e_integrals_in_map
+
+    double precision, external :: mo_two_e_integral
+
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,mo_num
+      do k=1,mo_num
+        do j=l,mo_num
+          do i=k,mo_num
+            if (i==j .and. k<l) cycle
+            if (i<j) cycle
+            integral = mo_two_e_integral(i,j,k,l)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+              call trexio_assert(rc, TREXIO_SUCCESS)
+              offset += icount
+              icount = 0_8
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount > 0_8) then
+      rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    end if
+  end if
+
+! Two-e MO integrals - Cholesky
+! -----------------------------
+
+  if (export_mo_two_e_ints_cholesky) then
+    print *, 'MO two-e integrals Cholesky'
+
+    rc = trexio_write_mo_2e_int_eri_cholesky_num(f, cholesky_ao_num)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    icount = 0_8
+    offset = 0_8
+    do k=1,cholesky_ao_num
+     do j=1,mo_num
+      do i=1,mo_num
+         integral = cholesky_mo(i,j,k)
+         if (integral == 0.d0) cycle
+         icount += 1_8
+         chol_buffer(icount) = integral
+         chol_index(1,icount) = i
+         chol_index(2,icount) = j
+         chol_index(3,icount) = k
+         if (icount == BUFSIZE) then
+           rc = trexio_write_mo_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
+           call trexio_assert(rc, TREXIO_SUCCESS)
+           offset += icount
+           icount = 0_8
+         end if
+      end do
+     end do
+    end do
+
+    if (icount > 0_8) then
+      rc = trexio_write_mo_2e_int_eri_cholesky(f, offset, icount, chol_index, chol_buffer)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    end if
+  end if
+
+
+! One-e RDM
+! ---------
+
+  rc = trexio_write_rdm_1e(f,one_e_dm_mo)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_rdm_1e_up(f,one_e_dm_mo_alpha_average)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+  rc = trexio_write_rdm_1e_dn(f,one_e_dm_mo_beta_average)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+
+! Two-e RDM
+! ---------
+
+  if (export_rdm) then
+    PROVIDE two_e_dm_mo
+    print *, 'Two-e RDM'
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,mo_num
+      do k=1,mo_num
+        do j=1,mo_num
+          do i=1,mo_num
+            integral = two_e_dm_mo(i,j,k,l)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
+              call trexio_assert(rc, TREXIO_SUCCESS)
+              offset += icount
+              icount = 0_8
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    end if
+  end if
+
+
+! ------------------------------------------------------------------------------
+
+  ! Determinants
+  ! ------------
+
+    integer*8, allocatable :: det_buffer(:,:,:)
+    double precision, allocatable :: coef_buffer(:,:)
+    integer :: nint
+
+!    rc = trexio_read_determinant_int64_num(f, nint)
+!    call trexio_assert(rc, TREXIO_SUCCESS)
+    nint = N_int
+    if (nint /= N_int) then
+       stop 'Problem with N_int'
+    endif
+    allocate ( det_buffer(nint, 2, BUFSIZE), coef_buffer(BUFSIZE, n_states) )
+
+    icount = 0_8
+    offset = 0_8
+    rc = trexio_write_state_num(f, n_states)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+
+    rc = trexio_set_state (f, 0)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+    do k=1,n_det
+       icount += 1_8
+       det_buffer(1:nint, 1:2, icount) = psi_det(1:N_int, 1:2, k)
+       coef_buffer(icount,1:N_states) = psi_coef(k,1:N_states)
+       if (icount == BUFSIZE) then
+         call trexio_assert(rc, TREXIO_SUCCESS)
+         rc = trexio_write_determinant_list(f, offset, icount, det_buffer)
+         call trexio_assert(rc, TREXIO_SUCCESS)
+         do i=1,N_states
+           rc = trexio_set_state (f, i-1)
+           call trexio_assert(rc, TREXIO_SUCCESS)
+           rc = trexio_write_determinant_coefficient(f, offset, icount, coef_buffer(1,i))
+         end do
+         rc = trexio_set_state (f, 0)
+         offset += icount
+         icount = 0_8
+       end if
+    end do
+
+  if (icount >= 0_8) then
+       call trexio_assert(rc, TREXIO_SUCCESS)
+       rc = trexio_write_determinant_list(f, offset, icount, det_buffer)
+       call trexio_assert(rc, TREXIO_SUCCESS)
+       do i=1,N_states
+         rc = trexio_set_state (f, i-1)
+         call trexio_assert(rc, TREXIO_SUCCESS)
+         rc = trexio_write_determinant_coefficient(f, offset, icount, coef_buffer(1,i))
+       end do
+       rc = trexio_set_state (f, 0)
+  end if
+
+  deallocate ( det_buffer, coef_buffer )
+
+  rc = trexio_close(f)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+end
+
+
+! -*- mode: f90 -*-
diff --git a/src/trexio/import_trexio_determinants.irp.f b/src/trexio/import_trexio_determinants.irp.f
new file mode 100644
index 00000000..1759bb94
--- /dev/null
+++ b/src/trexio/import_trexio_determinants.irp.f
@@ -0,0 +1,79 @@
+program import_determinants_ao
+  call run
+end
+
+subroutine run
+  use trexio
+  use map_module
+  implicit none
+  BEGIN_DOC
+! Program to import determinants from TREXIO
+  END_DOC
+
+  integer(trexio_t)              :: f ! TREXIO file handle
+  integer(trexio_exit_code)      :: rc
+
+  integer :: m
+
+  double precision, allocatable :: coef_buffer(:,:)
+  integer*8       , allocatable :: det_buffer(:,:,:)
+
+  f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
+  if (f == 0_8) then
+    print *, 'Unable to open TREXIO file for reading'
+    print *, 'rc = ', rc
+    stop -1
+  endif
+
+
+
+  ! Determinants
+  ! ------------
+
+  integer :: nint, nstates
+  integer :: bufsize
+
+  rc = trexio_read_state_num(f, nstates)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+
+!  rc = trexio_read_determinant_int64_num(f, nint)
+!  call trexio_assert(rc, TREXIO_SUCCESS)
+  nint = N_int
+  if (nint /= N_int) then
+     stop 'Problem with N_int'
+  endif
+
+  integer*8 :: offset, icount
+
+  rc = trexio_read_determinant_num(f, bufsize)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+  print *, 'N_det = ', bufsize
+
+  allocate ( det_buffer(nint, 2, bufsize), coef_buffer(bufsize, n_states) )
+
+
+  offset = 0_8
+  icount = bufsize
+
+  rc = trexio_read_determinant_list(f, offset, icount, det_buffer)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+  if (icount /= bufsize) then
+      print *, 'error: bufsize /= N_det: ', bufsize, icount
+      stop -1
+  endif
+
+  do m=1,nstates
+    rc = trexio_set_state(f, m-1)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+    rc = trexio_read_determinant_coefficient(f, offset, icount, coef_buffer(1,m))
+    call trexio_assert(rc, TREXIO_SUCCESS)
+    if (icount /= bufsize) then
+        print *, 'error: bufsize /= N_det for state', m, ':', icount, bufsize
+        stop -1
+    endif
+  enddo
+
+  call save_wavefunction_general(bufsize,nstates,det_buffer,size(coef_buffer,1),coef_buffer)
+
+
+end
diff --git a/src/trexio/import_trexio_integrals.irp.f b/src/trexio/import_trexio_integrals.irp.f
new file mode 100644
index 00000000..9f9ad9d6
--- /dev/null
+++ b/src/trexio/import_trexio_integrals.irp.f
@@ -0,0 +1,146 @@
+program import_integrals_ao
+  use trexio
+  implicit none
+  integer(trexio_t)              :: f ! TREXIO file handle
+  integer(trexio_exit_code)      :: rc
+
+  f = trexio_open(trexio_filename, 'r', TREXIO_AUTO, rc)
+  if (f == 0_8) then
+    print *, 'Unable to open TREXIO file for reading'
+    print *, 'rc = ', rc
+    stop -1
+  endif
+
+  call run(f)
+  rc = trexio_close(f)
+  call trexio_assert(rc, TREXIO_SUCCESS)
+end
+
+subroutine run(f)
+  use trexio
+  use map_module
+  implicit none
+  BEGIN_DOC
+! Program to import integrals from TREXIO
+  END_DOC
+
+  integer(trexio_t), intent(in)  :: f ! TREXIO file handle
+  integer(trexio_exit_code)      :: rc
+
+  integer ::i,j,k,l
+  integer(8) :: m, n_integrals
+  double precision :: integral
+
+  integer(key_kind), allocatable   :: buffer_i(:)
+  real(integral_kind), allocatable :: buffer_values(:)
+
+
+  double precision, allocatable :: A(:,:)
+  double precision, allocatable :: V(:)
+  integer         , allocatable :: Vi(:,:)
+  double precision              :: s
+
+  if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
+    rc = trexio_read_nucleus_repulsion(f, s)
+    call trexio_assert(rc, TREXIO_SUCCESS)
+    if (rc /= TREXIO_SUCCESS) then
+      print *, irp_here, rc
+      print *, 'Error reading nuclear repulsion'
+      stop -1
+    endif
+    call ezfio_set_nuclei_nuclear_repulsion(s)
+    call ezfio_set_nuclei_io_nuclear_repulsion('Read')
+  endif
+
+  ! AO integrals
+  ! ------------
+
+  allocate(A(ao_num, ao_num))
+
+
+  if (trexio_has_ao_1e_int_overlap(f) == TREXIO_SUCCESS) then
+    rc = trexio_read_ao_1e_int_overlap(f, A)
+    if (rc /= TREXIO_SUCCESS) then
+      print *, irp_here
+      print *, 'Error reading AO overlap'
+      stop -1
+    endif
+    call ezfio_set_ao_one_e_ints_ao_integrals_overlap(A)
+    call ezfio_set_ao_one_e_ints_io_ao_integrals_overlap('Read')
+  endif
+
+  if (trexio_has_ao_1e_int_kinetic(f) == TREXIO_SUCCESS) then
+    rc = trexio_read_ao_1e_int_kinetic(f, A)
+    if (rc /= TREXIO_SUCCESS) then
+      print *, irp_here
+      print *, 'Error reading AO kinetic integrals'
+      stop -1
+    endif
+    call ezfio_set_ao_one_e_ints_ao_integrals_kinetic(A)
+    call ezfio_set_ao_one_e_ints_io_ao_integrals_kinetic('Read')
+  endif
+
+!  if (trexio_has_ao_1e_int_ecp(f) == TREXIO_SUCCESS) then
+!    rc = trexio_read_ao_1e_int_ecp(f, A)
+!    if (rc /= TREXIO_SUCCESS) then
+!      print *, irp_here
+!      print *, 'Error reading AO ECP local integrals'
+!      stop -1
+!    endif
+!    call ezfio_set_ao_one_e_ints_ao_integrals_pseudo(A)
+!    call ezfio_set_ao_one_e_ints_io_ao_integrals_pseudo('Read')
+!  endif
+
+  if (trexio_has_ao_1e_int_potential_n_e(f) == TREXIO_SUCCESS) then
+    rc = trexio_read_ao_1e_int_potential_n_e(f, A)
+    if (rc /= TREXIO_SUCCESS) then
+      print *, irp_here
+      print *, 'Error reading AO potential N-e integrals'
+      stop -1
+    endif
+    call ezfio_set_ao_one_e_ints_ao_integrals_n_e(A)
+    call ezfio_set_ao_one_e_ints_io_ao_integrals_n_e('Read')
+  endif
+
+  deallocate(A)
+
+  ! AO 2e integrals
+  ! ---------------
+  PROVIDE ao_integrals_map
+
+  integer*4 :: BUFSIZE
+  BUFSIZE=ao_num**2
+  allocate(buffer_i(BUFSIZE), buffer_values(BUFSIZE))
+  allocate(Vi(4,BUFSIZE), V(BUFSIZE))
+
+  integer*8 :: offset, icount
+
+  offset = 0_8
+  icount = BUFSIZE
+  rc = TREXIO_SUCCESS
+  do while (icount == size(V))
+    rc = trexio_read_ao_2e_int_eri(f, offset, icount, Vi, V)
+    do m=1,icount
+      i = Vi(1,m)
+      j = Vi(2,m)
+      k = Vi(3,m)
+      l = Vi(4,m)
+      integral = V(m)
+      call two_e_integrals_index(i, j, k, l, buffer_i(m) )
+      buffer_values(m) = integral
+    enddo
+    call insert_into_ao_integrals_map(int(icount,4),buffer_i,buffer_values)
+    offset = offset + icount
+    if (rc /= TREXIO_SUCCESS) then
+        exit
+    endif
+  end do
+  n_integrals = offset
+
+  call map_sort(ao_integrals_map)
+  call map_unique(ao_integrals_map)
+
+  call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map)
+  call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read')
+
+end
diff --git a/src/trexio/qp_import_trexio.py b/src/trexio/qp_import_trexio.py
new file mode 100755
index 00000000..de8d1269
--- /dev/null
+++ b/src/trexio/qp_import_trexio.py
@@ -0,0 +1,415 @@
+#!/usr/bin/env python3
+"""
+convert TREXIO file to EZFIO
+
+Usage:
+    qp_import_trexio [-o EZFIO_DIR] FILE
+
+Options:
+    -o --output=EZFIO_DIR    Produced directory
+                             by default is FILE.ezfio
+
+"""
+
+import sys
+import os
+import trexio
+import numpy as np
+from functools import reduce
+from ezfio import ezfio
+from docopt import docopt
+
+
+try:
+    QP_ROOT = os.environ["QP_ROOT"]
+    QP_EZFIO = os.environ["QP_EZFIO"]
+except KeyError:
+    print("Error: QP_ROOT environment variable not found.")
+    sys.exit(1)
+else:
+    sys.path = [QP_EZFIO + "/Python",
+                QP_ROOT + "/install/resultsFile",
+                QP_ROOT + "/install",
+                QP_ROOT + "/scripts"] + sys.path
+
+
+def generate_xyz(l):
+
+    def create_z(x,y,z):
+       return (x, y, l-(x+y))
+
+    def create_y(accu,x,y,z):
+       if y == 0:
+          result = [create_z(x,y,z)] + accu
+       else:
+          result = create_y([create_z(x,y,z)] + accu , x, y-1, z)
+       return result
+
+    def create_x(accu,x,y,z):
+       if x == 0:
+          result = create_y([], x,y,z) + accu
+       else:
+          xnew = x-1
+          ynew = l-xnew
+          result = create_x(create_y([],x,y,z) + accu , xnew, ynew, z)
+       return result
+
+    result = create_x([], l, 0, 0)
+    result.reverse()
+    return result
+
+
+
+def write_ezfio(trexio_filename, filename):
+
+    try:
+        trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_TEXT)
+    except:
+        trexio_file = trexio.File(trexio_filename,mode='r',back_end=trexio.TREXIO_HDF5)
+
+    ezfio.set_file(filename)
+    ezfio.set_trexio_trexio_file(trexio_filename)
+
+    print("Nuclei\t\t...\t", end=' ')
+
+    charge = [0.]
+    if trexio.has_nucleus(trexio_file):
+        charge = trexio.read_nucleus_charge(trexio_file)
+        ezfio.set_nuclei_nucl_num(len(charge))
+        ezfio.set_nuclei_nucl_charge(charge)
+
+        coord = trexio.read_nucleus_coord(trexio_file)
+        coord = np.transpose(coord)
+        ezfio.set_nuclei_nucl_coord(coord)
+
+        label = trexio.read_nucleus_label(trexio_file)
+        nucl_num = trexio.read_nucleus_num(trexio_file)
+
+        # Transformt H1 into H
+        import re
+        p = re.compile(r'(\d*)$')
+        label = [p.sub("", x).capitalize() for x in label]
+        ezfio.set_nuclei_nucl_label(label)
+
+    else:
+        ezfio.set_nuclei_nucl_num(1)
+        ezfio.set_nuclei_nucl_charge([0.])
+        ezfio.set_nuclei_nucl_coord([0.,0.,0.])
+        ezfio.set_nuclei_nucl_label(["X"])
+
+    print("OK")
+
+
+    print("Electrons\t...\t", end=' ')
+
+    try:
+        num_beta = trexio.read_electron_dn_num(trexio_file)
+    except:
+        num_beta = sum(charge)//2
+
+    try:
+        num_alpha = trexio.read_electron_up_num(trexio_file)
+    except:
+        num_alpha = sum(charge) - num_beta
+
+    if num_alpha == 0:
+        print("\n\nError: There are zero electrons in the TREXIO file.\n\n")
+        sys.exit(1)
+    ezfio.set_electrons_elec_alpha_num(num_alpha)
+    ezfio.set_electrons_elec_beta_num(num_beta)
+
+    print("OK")
+
+    print("Basis\t\t...\t", end=' ')
+
+    shell_num = 0
+    try:
+        basis_type = trexio.read_basis_type(trexio_file)
+
+        if basis_type.lower() not in ["gaussian", "slater"]:
+            raise TypeError
+
+        shell_num   = trexio.read_basis_shell_num(trexio_file)
+        prim_num    = trexio.read_basis_prim_num(trexio_file)
+        ang_mom     = trexio.read_basis_shell_ang_mom(trexio_file)
+        nucl_index  = trexio.read_basis_nucleus_index(trexio_file)
+        exponent    = trexio.read_basis_exponent(trexio_file)
+        coefficient = trexio.read_basis_coefficient(trexio_file)
+        shell_index = trexio.read_basis_shell_index(trexio_file)
+        ao_shell    = trexio.read_ao_shell(trexio_file)
+
+        ezfio.set_basis_basis("Read from TREXIO")
+        ezfio.set_basis_shell_num(shell_num)
+        ezfio.set_basis_prim_num(prim_num)
+        ezfio.set_basis_shell_ang_mom(ang_mom)
+        ezfio.set_basis_basis_nucleus_index([ x+1 for x in nucl_index ])
+        ezfio.set_basis_prim_expo(exponent)
+        ezfio.set_basis_prim_coef(coefficient)
+
+        nucl_shell_num = []
+        prev = None
+        m = 0
+        for i in ao_shell:
+            if i != prev:
+               m += 1
+               if prev is None or nucl_index[i] != nucl_index[prev]:
+                    nucl_shell_num.append(m)
+                    m = 0
+            prev = i
+        assert (len(nucl_shell_num) == nucl_num)
+
+        shell_prim_num = []
+        prev = shell_index[0]
+        count = 0
+        for i in shell_index:
+            if i != prev:
+               shell_prim_num.append(count)
+               count = 0
+            count += 1
+            prev = i
+        shell_prim_num.append(count)
+
+        assert (len(shell_prim_num) == shell_num)
+
+        ezfio.set_basis_shell_prim_num(shell_prim_num)
+        ezfio.set_basis_shell_index([x+1 for x in shell_index])
+        ezfio.set_basis_nucleus_shell_num(nucl_shell_num)
+
+
+        shell_factor = trexio.read_basis_shell_factor(trexio_file)
+        prim_factor  = trexio.read_basis_prim_factor(trexio_file)
+
+        print("OK")
+    except:
+        print("None")
+        ezfio.set_ao_basis_ao_cartesian(True)
+
+    print("AOS\t\t...\t", end=' ')
+
+    try:
+        cartesian = trexio.read_ao_cartesian(trexio_file)
+    except:
+        cartesian = True
+
+    if not cartesian:
+        raise TypeError('Only cartesian TREXIO files can be converted')
+
+    ao_num = trexio.read_ao_num(trexio_file)
+    ezfio.set_ao_basis_ao_num(ao_num)
+
+    if shell_num > 0:
+        ao_shell    = trexio.read_ao_shell(trexio_file)
+        at = [ nucl_index[i]+1 for i in ao_shell ]
+        ezfio.set_ao_basis_ao_nucl(at)
+
+        num_prim0 = [ 0 for i in range(shell_num) ]
+        for i in shell_index:
+            num_prim0[i] += 1
+
+        coef = {}
+        expo = {}
+        for i,c in enumerate(coefficient):
+            idx = shell_index[i]
+            if idx in coef:
+              coef[idx].append(c)
+              expo[idx].append(exponent[i])
+            else:
+              coef[idx] = [c]
+              expo[idx] = [exponent[i]]
+
+        coefficient = []
+        exponent    = []
+        power_x     = []
+        power_y     = []
+        power_z     = []
+        num_prim    = []
+
+        for i in range(shell_num):
+            for x,y,z in generate_xyz(ang_mom[i]):
+                power_x.append(x)
+                power_y.append(y)
+                power_z.append(z)
+                coefficient.append(coef[i])
+                exponent.append(expo[i])
+                num_prim.append(num_prim0[i])
+
+        assert (len(coefficient) == ao_num)
+        ezfio.set_ao_basis_ao_power(power_x + power_y + power_z)
+        ezfio.set_ao_basis_ao_prim_num(num_prim)
+
+        prim_num_max = max( [ len(x) for x in coefficient ] )
+
+        for i in range(ao_num):
+            coefficient[i] += [0. for j in range(len(coefficient[i]), prim_num_max)]
+            exponent   [i] += [0. for j in range(len(exponent[i]), prim_num_max)]
+
+        coefficient = reduce(lambda x, y: x + y, coefficient, [])
+        exponent    = reduce(lambda x, y: x + y, exponent   , [])
+
+        coef = []
+        expo = []
+        for i in range(prim_num_max):
+            for j in range(i, len(coefficient), prim_num_max):
+                coef.append(coefficient[j])
+                expo.append(exponent[j])
+
+#        ezfio.set_ao_basis_ao_prim_num_max(prim_num_max)
+        ezfio.set_ao_basis_ao_coef(coef)
+        ezfio.set_ao_basis_ao_expo(expo)
+        ezfio.set_ao_basis_ao_basis("Read from TREXIO")
+
+    print("OK")
+
+
+    #                _
+    # |\/|  _   _   |_)  _.  _ o  _
+    # |  | (_) _>   |_) (_| _> | _>
+    #
+
+    print("MOS\t\t...\t", end=' ')
+
+    labels = { "Canonical" : "Canonical",
+               "RHF" : "Canonical",
+               "BOYS" : "Localized",
+               "ROHF" : "Canonical",
+               "UHF" : "Canonical",
+               "Natural": "Natural" }
+    try:
+      label = labels[trexio.read_mo_type(trexio_file)]
+    except:
+      label = "None"
+    ezfio.set_mo_basis_mo_label(label)
+
+    try:
+      clss = trexio.read_mo_class(trexio_file)
+      core     = [ i for i in clss if i.lower() == "core" ]
+      inactive = [ i for i in clss if i.lower() == "inactive" ]
+      active   = [ i for i in clss if i.lower() == "active" ]
+      virtual  = [ i for i in clss if i.lower() == "virtual" ]
+      deleted  = [ i for i in clss if i.lower() == "deleted" ]
+    except trexio.Error:
+      pass
+
+    try:
+      mo_num = trexio.read_mo_num(trexio_file)
+      ezfio.set_mo_basis_mo_num(mo_num)
+
+      MoMatrix = trexio.read_mo_coefficient(trexio_file)
+      ezfio.set_mo_basis_mo_coef(MoMatrix)
+
+      mo_occ = [ 0. for i in range(mo_num) ]
+      for i in range(num_alpha):
+         mo_occ[i] += 1.
+      for i in range(num_beta):
+         mo_occ[i] += 1.
+      ezfio.set_mo_basis_mo_occ(mo_occ)
+    except:
+      pass
+
+    print("OK")
+
+
+    print("Pseudos\t\t...\t", end=' ')
+
+    ezfio.set_pseudo_do_pseudo(False)
+
+    if trexio.has_ecp_ang_mom(trexio_file):
+        ezfio.set_pseudo_do_pseudo(True)
+        max_ang_mom_plus_1 = trexio.read_ecp_max_ang_mom_plus_1(trexio_file)
+        z_core = trexio.read_ecp_z_core(trexio_file)
+        ang_mom = trexio.read_ecp_ang_mom(trexio_file)
+        nucleus_index = trexio.read_ecp_nucleus_index(trexio_file)
+        exponent = trexio.read_ecp_exponent(trexio_file)
+        coefficient = trexio.read_ecp_coefficient(trexio_file)
+        power = trexio.read_ecp_power(trexio_file)
+
+        lmax = max( max_ang_mom_plus_1 ) - 1
+        ezfio.set_pseudo_pseudo_lmax(lmax)
+        ezfio.set_pseudo_nucl_charge_remove(z_core)
+
+        prev_center = None
+        ecp = {}
+        for i in range(len(ang_mom)):
+            center = nucleus_index[i]
+            if center != prev_center:
+               ecp[center] = { "lmax": max_ang_mom_plus_1[center],
+                               "zcore": z_core[center],
+                               "contr": {} }
+               for j in range(max_ang_mom_plus_1[center]+1):
+                    ecp[center]["contr"][j] = []
+
+            ecp[center]["contr"][ang_mom[i]].append( (coefficient[i], power[i], exponent[i]) )
+            prev_center = center
+
+        ecp_loc = {}
+        ecp_nl  = {}
+        kmax    = 0
+        klocmax    = 0
+        for center in ecp:
+            ecp_nl [center] = {}
+            for k in ecp[center]["contr"]:
+                if k == ecp[center]["lmax"]:
+                    ecp_loc[center] = ecp[center]["contr"][k]
+                    klocmax = max(len(ecp_loc[center]), klocmax)
+                else:
+                    ecp_nl [center][k] = ecp[center]["contr"][k]
+                    kmax = max(len(ecp_nl [center][k]), kmax)
+
+        ezfio.set_pseudo_pseudo_klocmax(klocmax)
+        ezfio.set_pseudo_pseudo_kmax(kmax)
+
+        pseudo_n_k = [[0  for _ in range(nucl_num)] for _ in range(klocmax)]
+        pseudo_v_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)]
+        pseudo_dz_k = [[0. for _ in range(nucl_num)] for _ in range(klocmax)]
+        pseudo_n_kl = [[[0  for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
+        pseudo_v_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
+        pseudo_dz_kl = [[[0. for _ in range(nucl_num)] for _ in range(kmax)] for _ in range(lmax+1)]
+        for center in ecp_loc:
+            for k in range( len(ecp_loc[center]) ):
+                v, n, dz = ecp_loc[center][k]
+                pseudo_n_k[k][center] = n
+                pseudo_v_k[k][center] = v
+                pseudo_dz_k[k][center] = dz
+
+        ezfio.set_pseudo_pseudo_n_k(pseudo_n_k)
+        ezfio.set_pseudo_pseudo_v_k(pseudo_v_k)
+        ezfio.set_pseudo_pseudo_dz_k(pseudo_dz_k)
+
+        for center in ecp_nl:
+            for l in range( len(ecp_nl[center]) ):
+                for k in range( len(ecp_nl[center][l]) ):
+                    v, n, dz = ecp_nl[center][l][k]
+                    pseudo_n_kl[l][k][center] = n
+                    pseudo_v_kl[l][k][center] = v
+                    pseudo_dz_kl[l][k][center] = dz
+
+        ezfio.set_pseudo_pseudo_n_kl(pseudo_n_kl)
+        ezfio.set_pseudo_pseudo_v_kl(pseudo_v_kl)
+        ezfio.set_pseudo_pseudo_dz_kl(pseudo_dz_kl)
+
+
+    print("OK")
+
+
+
+
+def get_full_path(file_path):
+    file_path = os.path.expanduser(file_path)
+    file_path = os.path.expandvars(file_path)
+    return file_path
+
+
+if __name__ == '__main__':
+    ARGUMENTS = docopt(__doc__)
+
+    FILE = get_full_path(ARGUMENTS['FILE'])
+    trexio_filename = FILE
+
+    if ARGUMENTS["--output"]:
+        EZFIO_FILE = get_full_path(ARGUMENTS["--output"])
+    else:
+        EZFIO_FILE = "{0}.ezfio".format(FILE)
+
+    write_ezfio(trexio_filename, EZFIO_FILE)
+    sys.stdout.flush()
+
diff --git a/src/trexio/trexio_file.irp.f b/src/trexio/trexio_file.irp.f
new file mode 100644
index 00000000..c9897748
--- /dev/null
+++ b/src/trexio/trexio_file.irp.f
@@ -0,0 +1,20 @@
+BEGIN_PROVIDER [ character*(1024), trexio_filename ]
+ implicit none
+ BEGIN_DOC
+ ! Name of the TREXIO file
+ END_DOC
+ character*(1024) :: prefix
+
+ trexio_filename = trexio_file
+
+ if (trexio_file == 'None') then
+    prefix = trim(ezfio_work_dir)//trim(ezfio_filename)
+    if (backend == 0) then
+      trexio_filename = trim(prefix)//'.h5'
+    else if (backend == 1) then
+      trexio_filename = trim(prefix)
+    endif
+ endif
+END_PROVIDER
+
+
diff --git a/src/trexio/trexio_module.F90 b/src/trexio/trexio_module.F90
new file mode 100644
index 00000000..acd08492
--- /dev/null
+++ b/src/trexio/trexio_module.F90
@@ -0,0 +1 @@
+#include "trexio_f.f90"