diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f
index a0f9f897..66b9deb2 100644
--- a/src/hartree_fock/hf_energy.irp.f
+++ b/src/hartree_fock/hf_energy.irp.f
@@ -22,13 +22,39 @@ END_PROVIDER
  HF_energy = nuclear_repulsion
  HF_two_electron_energy = 0.d0
  HF_one_electron_energy = 0.d0
- do j=1,ao_num
-   do i=1,ao_num
-    HF_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) &
-                                       +ao_two_e_integral_beta(i,j)  * SCF_density_matrix_ao_beta(i,j) )
-    HF_one_electron_energy += ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
-   enddo
- enddo
+  if (is_periodic) then
+    complex*16 :: hf_1e_tmp, hf_2e_tmp
+    hf_1e_tmp = (0.d0,0.d0)
+    hf_2e_tmp = (0.d0,0.d0)
+    do j=1,ao_num
+      do i=1,ao_num
+        hf_2e_tmp += 0.5d0 * ( ao_two_e_integral_alpha_complex(i,j) * SCF_density_matrix_ao_alpha_complex(j,i) &
+                              +ao_two_e_integral_beta_complex(i,j)  * SCF_density_matrix_ao_beta_complex(j,i) )
+        hf_1e_tmp += ao_one_e_integrals_complex(i,j) * (SCF_density_matrix_ao_alpha_complex(j,i) &
+                                                      + SCF_density_matrix_ao_beta_complex (j,i) )
+      enddo
+    enddo
+    if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then
+      print*,'HF_2e energy should be real:',irp_here
+      stop -1
+    else
+      HF_two_electron_energy = dble(hf_2e_tmp)
+    endif
+    if (dabs(dimag(hf_1e_tmp)).gt.1.d-10) then
+      print*,'HF_1e energy should be real:',irp_here
+      stop -1
+    else
+      HF_one_electron_energy = dble(hf_1e_tmp)
+    endif
+  else
+    do j=1,ao_num
+      do i=1,ao_num
+        HF_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) &
+                                           +ao_two_e_integral_beta(i,j)  * SCF_density_matrix_ao_beta(i,j) )
+        HF_one_electron_energy += ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
+      enddo
+    enddo
+  endif
  HF_energy += HF_two_electron_energy + HF_one_electron_energy
 END_PROVIDER
 
diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f
index 9c4d54e7..b59f921b 100644
--- a/src/scf_utils/fock_matrix.irp.f
+++ b/src/scf_utils/fock_matrix.irp.f
@@ -164,8 +164,8 @@ BEGIN_PROVIDER [ double precision, SCF_energy ]
    do j=1,ao_num
      do i=1,ao_num
        scf_e_tmp += 0.5d0 * (                                          &
-           (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_alpha_complex(i,j) ) *  SCF_density_matrix_ao_alpha_complex(i,j) +&
-           (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_beta_complex (i,j) ) *  SCF_density_matrix_ao_beta_complex (i,j) )
+           (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_alpha_complex(i,j) ) *  SCF_density_matrix_ao_alpha_complex(j,i) +&
+           (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_beta_complex (i,j) ) *  SCF_density_matrix_ao_beta_complex (j,i) )
      enddo
    enddo
    !TODO: add check for imaginary part? (should be zero)