diff --git a/src/ed_analysis.F b/src/ed_analysis.F
index 694fc807cf..a5eb85b736 100644
--- a/src/ed_analysis.F
+++ b/src/ed_analysis.F
@@ -582,6 +582,8 @@ CONTAINS
             atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
          END DO
          atmul(1:natom, 2) = atmul(1:natom, 2) - atmul(1:natom, 1)
+         atmul(1:natom, 3) = atmul(1:natom, 3) + 0.5_dp*atmul(1:natom, 2)
+         atmul(1:natom, 2) = 0.5_dp*atmul(1:natom, 2) + 0.5_dp*atcore(1:natom)
       END IF
       !
       CALL dbcsr_release(dkmat)
@@ -623,7 +625,7 @@ CONTAINS
       ! KS energy
       atener(1:natom) = ateks(1:natom)
       ! 1/2 of VPP contribution Tr[VPP(K)*P]
-      CALL group%sum(atcore)  ! maybe not needed
+      CALL group%sum(atcore)
       atener(1:natom) = atener(1:natom) + 0.5_dp*atcore(1:natom)
       ! core energy corrections
       CALL group%sum(atecc)
@@ -673,7 +675,10 @@ CONTAINS
       ekts = energy%kts/REAL(natom, KIND=dp)
       atener(1:natom) = atener(1:natom) + ekts
       ! 0.5 Vpp(at)*D + 0.5 * Vpp*D(at)
-      IF (detailed_ener) atdet(1:natom, 2) = 0.5_dp*atdet(1:natom, 2) + 0.5_dp*atcore(1:natom)
+      IF (detailed_ener) THEN
+         atdet(1:natom, 3) = atdet(1:natom, 3) + 0.5_dp*atdet(1:natom, 2)
+         atdet(1:natom, 2) = 0.5_dp*atdet(1:natom, 2) + 0.5_dp*atcore(1:natom)
+      END IF
       !
       IF (detailed_ener) THEN
          IF (unit_nr > 0) THEN
diff --git a/src/qs_energy_utils.F b/src/qs_energy_utils.F
index 12672382e1..eb866c73a2 100644
--- a/src/qs_energy_utils.F
+++ b/src/qs_energy_utils.F
@@ -40,6 +40,7 @@ MODULE qs_energy_utils
                                               pw_scale
    USE pw_pool_types,                   ONLY: pw_pool_type
    USE pw_types,                        ONLY: pw_r3d_rs_type
+   USE qs_core_hamiltonian,             ONLY: core_matrices
    USE qs_dispersion_types,             ONLY: qs_dispersion_type
    USE qs_energy_types,                 ONLY: qs_energy_type
    USE qs_environment_types,            ONLY: get_qs_env,&
@@ -227,25 +228,54 @@ CONTAINS
 
       TYPE(qs_environment_type), POINTER                 :: qs_env
 
-      INTEGER                                            :: ispin
+      INTEGER                                            :: ispin, natom, nspin
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: atcore
       TYPE(atprop_type), POINTER                         :: atprop
+      TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), &
+         TARGET                                          :: core_mat
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_h, matrix_ks, rho_ao
+      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: math, matp
+      TYPE(dft_control_type), POINTER                    :: dft_control
       TYPE(qs_rho_type), POINTER                         :: rho
 
-      NULLIFY (atprop, matrix_h, matrix_ks, rho, rho_ao)
-      CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_h=matrix_h, &
-                      rho=rho, atprop=atprop)
-      CALL qs_rho_get(rho, rho_ao=rho_ao)
-
+      CALL get_qs_env(qs_env=qs_env, atprop=atprop)
       IF (atprop%energy) THEN
+         CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_h=matrix_h, rho=rho)
+         CALL qs_rho_get(rho, rho_ao=rho_ao)
          ! E = 0.5*Tr(H*P+F*P)
          atprop%atener = 0._dp
-         DO ispin = 1, SIZE(rho_ao)
+         nspin = SIZE(rho_ao)
+         DO ispin = 1, nspin
             CALL atom_trace(matrix_h(1)%matrix, rho_ao(ispin)%matrix, &
                             0.5_dp, atprop%atener)
             CALL atom_trace(matrix_ks(ispin)%matrix, rho_ao(ispin)%matrix, &
                             0.5_dp, atprop%atener)
          END DO
+         !
+         CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
+         IF (.NOT. dft_control%qs_control%semi_empirical .AND. &
+             .NOT. dft_control%qs_control%xtb .AND. &
+             .NOT. dft_control%qs_control%dftb) THEN
+            CALL get_qs_env(qs_env=qs_env, natom=natom)
+            ALLOCATE (atcore(natom))
+            ALLOCATE (core_mat(1))
+            ALLOCATE (core_mat(1)%matrix)
+            CALL dbcsr_create(core_mat(1)%matrix, template=matrix_h(1)%matrix)
+            CALL dbcsr_copy(core_mat(1)%matrix, matrix_h(1)%matrix)
+            CALL dbcsr_set(core_mat(1)%matrix, 0.0_dp)
+            math(1:1, 1:1) => core_mat(1:1)
+            matp(1:nspin, 1:1) => rho_ao(1:nspin)
+            CALL core_matrices(qs_env, math, matp, .FALSE., 0, atcore=atcore)
+            atprop%atener = atprop%atener + 0.5_dp*atcore
+            DO ispin = 1, nspin
+               CALL atom_trace(core_mat(1)%matrix, rho_ao(ispin)%matrix, &
+                               -0.5_dp, atprop%atener)
+            END DO
+            DEALLOCATE (atcore)
+            CALL dbcsr_release(core_mat(1)%matrix)
+            DEALLOCATE (core_mat(1)%matrix)
+            DEALLOCATE (core_mat)
+         END IF
       END IF
 
    END SUBROUTINE qs_energies_mulliken