diff --git a/src/tools/print_energy.irp.f b/src/tools/print_energy.irp.f
index 4703e7d4..056531a0 100644
--- a/src/tools/print_energy.irp.f
+++ b/src/tools/print_energy.irp.f
@@ -8,7 +8,11 @@ program print_energy
  ! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
  read_wf = .True.
  touch read_wf
- call run
+ if (is_complex) then
+  call run_complex
+ else
+  call run
+ endif
 end
 
 subroutine run
@@ -32,3 +36,25 @@ subroutine run
    print *, E(i)/norm(i)
  enddo
 end
+
+subroutine run_complex
+  implicit none
+  integer :: i
+  complex*16 :: i_h_psi_array(n_states)
+  double precision :: e(n_states)
+  double precision :: norm(n_states)
+ 
+  e(:) = nuclear_repulsion
+  norm(:) = 0.d0
+  do i=1,n_det
+   call i_H_psi_complex(psi_det(1,1,i), psi_det, psi_coef_complex, N_int, N_det, &
+                size(psi_coef_complex,1), N_states, i_H_psi_array)
+   norm(:) += cdabs(psi_coef_complex(i,:))**2
+   E(:) += dble(i_h_psi_array(:) * dconjg(psi_coef_complex(i,:)))
+  enddo
+ 
+  print *, 'Energy:'
+  do i=1,N_states
+    print *, E(i)/norm(i)
+  enddo
+end