diff --git a/devel/import_integrals/export_integrals_ao.irp.f b/devel/import_integrals/export_integrals_ao.irp.f
new file mode 100644
index 0000000..b192df8
--- /dev/null
+++ b/devel/import_integrals/export_integrals_ao.irp.f
@@ -0,0 +1,89 @@
+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
diff --git a/devel/import_integrals/export_integrals_mo.irp.f b/devel/import_integrals/export_integrals_mo.irp.f
new file mode 100644
index 0000000..e4cea9b
--- /dev/null
+++ b/devel/import_integrals/export_integrals_mo.irp.f
@@ -0,0 +1,80 @@
+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
diff --git a/stable/champ/qp_convert.py b/stable/champ/qp_convert.py
index 5c61fbb..0bae35c 100755
--- a/stable/champ/qp_convert.py
+++ b/stable/champ/qp_convert.py
@@ -1,4 +1,4 @@
-#!/usr/bin/env python2
+#!/usr/bin/env python3
 #
 # Modified from the QMCPACK interface developed by @tapplencourt and @abenali
 
diff --git a/stable/mp2/H_apply.irp.f b/stable/mp2/H_apply.irp.f
index 4d6e5a4..471dde5 100644
--- a/stable/mp2/H_apply.irp.f
+++ b/stable/mp2/H_apply.irp.f
@@ -1,15 +1,15 @@
 use bitmasks
-BEGIN_SHELL [ /usr/bin/env python2 ]
+BEGIN_SHELL [ /usr/bin/env python3 ]
 from generate_h_apply import *
 from perturbation import perturbations
 
 s = H_apply("mp2")
 s.set_perturbation("Moller_plesset")
 #s.set_perturbation("epstein_nesbet")
-print s
+print(s)
 
 s = H_apply("mp2_selection")
 s.set_selection_pt2("Moller_Plesset")
-print s
+print(s)
 END_SHELL