diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py
index dc7d340e..abfea976 100644
--- a/scripts/generate_h_apply.py
+++ b/scripts/generate_h_apply.py
@@ -130,9 +130,15 @@ class H_apply(object):
   !$OMP END PARALLEL
 
   call dsort(H_jj,iorder,N_det)
-  do k=1,N_states
-    psi_coef(iorder(k),k) = 1.d0
-  enddo
+  if (is_complex) then
+    do k=1,N_states
+      psi_coef_complex(iorder(k),k) = (1.d0,0.d0)
+    enddo
+  else
+    do k=1,N_states
+      psi_coef(iorder(k),k) = 1.d0
+    enddo
+  endif
   deallocate(H_jj,iorder)
     """
 
@@ -141,7 +147,11 @@ class H_apply(object):
   if (s2_eig) then
     call make_s2_eigenfunction
   endif
-  SOFT_TOUCH psi_det psi_coef N_det
+  if (is_complex) then
+    SOFT_TOUCH psi_det psi_coef_complex N_det
+  else
+    SOFT_TOUCH psi_det psi_coef N_det
+  endif
 """
     s["printout_now"]   = """write(6,*)  &
        100.*float(i_generator)/float(N_det_generators), '% in ', wall_1-wall_0, 's'"""
diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f
index 816253c5..5fd76493 100644
--- a/src/cis/cis.irp.f
+++ b/src/cis/cis.irp.f
@@ -77,8 +77,13 @@ subroutine run
   endif
 
   call ezfio_set_cis_energy(CI_energy)
-  psi_coef = ci_eigenvectors
-  SOFT_TOUCH psi_coef
+  if (is_complex) then
+    psi_coef_complex = ci_eigenvectors_complex
+    SOFT_TOUCH psi_coef_complex
+  else
+    psi_coef = ci_eigenvectors
+    SOFT_TOUCH psi_coef
+  endif
   call save_wavefunction
 
 end
diff --git a/src/cisd/cisd.irp.f b/src/cisd/cisd.irp.f
index 6c55e2ff..c3c9f821 100644
--- a/src/cisd/cisd.irp.f
+++ b/src/cisd/cisd.irp.f
@@ -61,16 +61,29 @@ subroutine run
   else
    call H_apply_cisd
   endif
-  psi_coef = ci_eigenvectors
-  SOFT_TOUCH psi_coef
+  if (is_complex) then
+    psi_coef_complex = ci_eigenvectors_complex
+    SOFT_TOUCH psi_coef_complex
+  else
+    psi_coef = ci_eigenvectors
+    SOFT_TOUCH psi_coef
+  endif
   call save_wavefunction
   call ezfio_set_cisd_energy(CI_energy)
 
-  do i = 1,N_states
-    k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1)
-    delta_E  = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
-    cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
-  enddo
+  if (is_complex) then
+    do i = 1,N_states
+      k = maxloc(cdabs(psi_coef_sorted_complex(1:N_det,i)),dim=1)
+      delta_E  = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
+      cisdq(i) = CI_energy(i) + delta_E * (1.d0 - cdabs(psi_coef_sorted_complex(k,i))**2)
+    enddo
+  else
+    do i = 1,N_states
+      k = maxloc(dabs(psi_coef_sorted(1:N_det,i)),dim=1)
+      delta_E  = CI_electronic_energy(i) - diag_h_mat_elem(psi_det_sorted(1,1,k),N_int)
+      cisdq(i) = CI_energy(i) + delta_E * (1.d0 - psi_coef_sorted(k,i)**2)
+    enddo
+  endif
   print *,  'N_det = ', N_det
   print*,''
   print*,'******************************'
diff --git a/src/cisd/cisd_routine.irp.f b/src/cisd/cisd_routine.irp.f
index 93b31e7d..e243e113 100644
--- a/src/cisd/cisd_routine.irp.f
+++ b/src/cisd/cisd_routine.irp.f
@@ -20,8 +20,13 @@ subroutine run_cisd
       print*, i ,CI_energy(i) - CI_energy(1)
     enddo
   endif
-  psi_coef = ci_eigenvectors
-  SOFT_TOUCH psi_coef
+  if (is_complex) then
+    psi_coef_complex = ci_eigenvectors_complex
+    SOFT_TOUCH psi_coef_complex
+  else
+    psi_coef = ci_eigenvectors
+    SOFT_TOUCH psi_coef
+  endif
   call save_wavefunction
   call ezfio_set_cisd_energy(CI_energy)
 
diff --git a/src/determinants/h_apply.irp.f b/src/determinants/h_apply.irp.f
index f73d328a..fcce1645 100644
--- a/src/determinants/h_apply.irp.f
+++ b/src/determinants/h_apply.irp.f
@@ -256,7 +256,7 @@ subroutine copy_H_apply_buffer_to_wf
 
     call remove_duplicates_in_psi_det(found_duplicates)
     do k=1,N_states
-      call normalize(psi_coef_complex(1,k),N_det)
+      call normalize_complex(psi_coef_complex(1,k),N_det)
     enddo
     SOFT_TOUCH N_det psi_det psi_coef_complex
   else