diff --git a/src/core_ae.F b/src/core_ae.F
index c8d92955ae..61034f3b9d 100644
--- a/src/core_ae.F
+++ b/src/core_ae.F
@@ -78,9 +78,11 @@ CONTAINS
 !> \param sac_ae ...
 !> \param nimages ...
 !> \param cell_to_index ...
+!> \param atcore ...
 ! **************************************************************************************************
    SUBROUTINE build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
-                            qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
+                            qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
+                            nimages, cell_to_index, atcore)
 
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p
       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
@@ -95,6 +97,8 @@ CONTAINS
          POINTER                                         :: sab_orb, sac_ae
       INTEGER, INTENT(IN)                                :: nimages
       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
+      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
+         OPTIONAL                                        :: atcore
 
       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_core_ae'
 
@@ -106,9 +110,10 @@ CONTAINS
       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
                                                             npgfb, nsgfa, nsgfb
       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
-      LOGICAL                                            :: dokp, found
-      REAL(KIND=dp)                                      :: alpha_c, core_charge, core_radius, dab, &
-                                                            dac, dbc, f0, rab2, rac2, rbc2, zeta_c
+      LOGICAL                                            :: doat, dokp, found
+      REAL(KIND=dp)                                      :: alpha_c, atk0, atk1, core_charge, &
+                                                            core_radius, dab, dac, dbc, f0, rab2, &
+                                                            rac2, rbc2, zeta_c
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ff
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: habd, work
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: hab, pab, verf, vnuc
@@ -119,6 +124,7 @@ CONTAINS
       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
       TYPE(all_potential_type), POINTER                  :: all_potential
+      REAL(KIND=dp), DIMENSION(SIZE(particle_set))       :: at_thread
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: h_block, p_block, rpgfa, rpgfb, sphi_a, &
                                                             sphi_b, zeta, zetb
       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
@@ -142,9 +148,10 @@ CONTAINS
       nkind = SIZE(atomic_kind_set)
       natom = SIZE(particle_set)
 
+      doat = PRESENT(atcore)
       dokp = (nimages > 1)
 
-      IF (calculate_forces) THEN
+      IF (calculate_forces .OR. doat) THEN
          IF (SIZE(matrix_p, 1) == 2) THEN
             DO img = 1, nimages
                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
@@ -156,6 +163,7 @@ CONTAINS
       END IF
 
       force_thread = 0.0_dp
+      at_thread = 0.0_dp
       pv_thread = 0.0_dp
 
       ALLOCATE (basis_set_list(nkind))
@@ -185,7 +193,7 @@ CONTAINS
 !$OMP SHARED  (ap_iterator, basis_set_list, calculate_forces, use_virial, &
 !$OMP          matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
 !$OMP          sab_orb, sac_ae, nthread, ncoset, nkind, cell_to_index, &
-!$OMP          slot, ldsab,  maxnset, ldai, nder, maxl, maxco, dokp, locks, natom) &
+!$OMP          slot, ldsab,  maxnset, ldai, nder, maxl, maxco, dokp, doat, locks, natom) &
 !$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
 !$OMP          first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
 !$OMP          zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
@@ -193,9 +201,9 @@ CONTAINS
 !$OMP          sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
 !$OMP          rac, dac, rbc, rab2, rac2, rbc2, dbc, na_plus, nb_plus, verf, vnuc, &
 !$OMP          set_radius_a,  core_radius, rpgfa, force_a, force_b, mepos, &
-!$OMP          habd, f0, katom, cellind, img, nij, ff, &
+!$OMP          atk0, atk1, habd, f0, katom, cellind, img, nij, ff, &
 !$OMP          sgp_potential, all_potential, hash, hash1, hash2, iatom8) &
-!$OMP REDUCTION (+ : pv_thread, force_thread )
+!$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )
 
 !$OMP SINGLE
 !$    ALLOCATE (locks(nlock))
@@ -212,7 +220,7 @@ CONTAINS
 
       ALLOCATE (hab(ldsab, ldsab, maxnset*maxnset), work(ldsab, ldsab))
       ALLOCATE (verf(ldai, ldai, 2*maxl + nder + 1), vnuc(ldai, ldai, 2*maxl + nder + 1), ff(0:2*maxl + nder))
-      IF (calculate_forces) THEN
+      IF (calculate_forces .OR. doat) THEN
          ALLOCATE (pab(maxco, maxco, maxnset*maxnset))
       END IF
 
@@ -281,7 +289,7 @@ CONTAINS
          NULLIFY (h_block)
          CALL dbcsr_get_block_p(matrix=matrix_h(1, img)%matrix, &
                                 row=irow, col=icol, BLOCK=h_block, found=found)
-         IF (calculate_forces) THEN
+         IF (calculate_forces .OR. doat) THEN
             NULLIFY (p_block)
             CALL dbcsr_get_block_p(matrix=matrix_p(1, img)%matrix, &
                                    row=irow, col=icol, BLOCK=p_block, found=found)
@@ -352,6 +360,9 @@ CONTAINS
                      rbc2 = dbc*dbc
                      nij = jset + (iset - 1)*maxnset
                      ! *** Calculate the GTH pseudo potential forces ***
+                     IF (doat) THEN
+                        atk0 = f0*SUM(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
+                     END IF
                      IF (calculate_forces) THEN
                         na_plus = npgfa(iset)*ncoset(la_max(iset) + nder)
                         nb_plus = npgfb(jset)*ncoset(lb_max(jset))
@@ -396,6 +407,11 @@ CONTAINS
                            alpha_c, core_radius, zeta_c, core_charge, &
                            rab, rab2, rac, rac2, rbc2, hab(:, :, nij), verf, vnuc, ff(0:))
                      END IF
+                     ! calculate atomic contributions
+                     IF (doat) THEN
+                        atk1 = f0*SUM(hab(1:ncoa, 1:ncob, nij)*pab(1:ncoa, 1:ncob, nij))
+                        at_thread(katom) = at_thread(katom) + (atk1 - atk0)
+                     END IF
                   END DO
                END DO
             END DO
@@ -472,6 +488,9 @@ CONTAINS
          END DO
 !$OMP END DO
       END IF
+      IF (doat) THEN
+         atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
+      END IF
 
       IF (calculate_forces .AND. use_virial) THEN
          virial%pv_ppl = virial%pv_ppl + pv_thread
diff --git a/src/core_ppl.F b/src/core_ppl.F
index eac058c654..061df7b0e3 100644
--- a/src/core_ppl.F
+++ b/src/core_ppl.F
@@ -88,10 +88,11 @@ CONTAINS
 !> \param cell_to_index ...
 !> \param basis_type ...
 !> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
+!> \param atcore ...
 ! **************************************************************************************************
    SUBROUTINE build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
                              qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
-                             nimages, cell_to_index, basis_type, deltaR)
+                             nimages, cell_to_index, basis_type, deltaR, atcore)
 
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p
       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
@@ -109,6 +110,8 @@ CONTAINS
       CHARACTER(LEN=*), INTENT(IN)                       :: basis_type
       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
          OPTIONAL                                        :: deltaR
+      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
+         OPTIONAL                                        :: atcore
 
       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_core_ppl'
       INTEGER, PARAMETER                                 :: nexp_max = 30
@@ -126,9 +129,11 @@ CONTAINS
                                                             nct_lpot, npgfa, npgfb, nsgfa, nsgfb
       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
       INTEGER, DIMENSION(nexp_max)                       :: nct_ppl
-      LOGICAL                                            :: do_dR, dokp, ecp_local, ecp_semi_local, &
-                                                            found, lpotextended, only_gaussians
-      REAL(KIND=dp)                                      :: alpha, dab, dac, dbc, f0, ppl_radius
+      LOGICAL                                            :: do_dR, doat, dokp, ecp_local, &
+                                                            ecp_semi_local, found, lpotextended, &
+                                                            only_gaussians
+      REAL(KIND=dp)                                      :: alpha, atk0, atk1, dab, dac, dbc, f0, &
+                                                            ppl_radius
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: qab, work
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: hab2_w, ppl_fwork, ppl_work
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: hab, pab
@@ -143,6 +148,7 @@ CONTAINS
       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
       TYPE(gth_potential_type), POINTER                  :: gth_potential
+      REAL(KIND=dp), DIMENSION(SIZE(particle_set))       :: at_thread
       REAL(KIND=dp), DIMENSION(nexp_max)                 :: alpha_ppl
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval_lpot, h1_1block, h1_2block, &
                                                             h1_3block, h_block, p_block, rpgfa, &
@@ -160,6 +166,7 @@ CONTAINS
 !$    INTEGER, PARAMETER                                 :: nlock = 501
 
       do_dR = PRESENT(deltaR)
+      doat = PRESENT(atcore)
       MARK_USED(int_8)
 
       IF (calculate_forces) THEN
@@ -177,7 +184,7 @@ CONTAINS
          CPASSERT(PRESENT(cell_to_index) .AND. ASSOCIATED(cell_to_index))
       END IF
 
-      IF (calculate_forces) THEN
+      IF (calculate_forces .OR. doat) THEN
          IF (SIZE(matrix_p, 1) == 2) THEN
             DO img = 1, nimages
                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
@@ -188,6 +195,7 @@ CONTAINS
          END IF
       END IF
       force_thread = 0.0_dp
+      at_thread = 0.0_dp
 
       maxder = ncoset(nder)
 
@@ -224,14 +232,14 @@ CONTAINS
 !$OMP SHARED  (ap_iterator, basis_set_list, calculate_forces, use_virial, &
 !$OMP          matrix_h, matrix_p, atomic_kind_set, qs_kind_set, particle_set, &
 !$OMP          sab_orb, sac_ppl, nthread, ncoset, nkind, cell_to_index, &
-!$OMP          ldsab,  maxnset, maxder, do_dR, deltaR, &
+!$OMP          ldsab,  maxnset, maxder, do_dR, deltaR, doat, &
 !$OMP          maxlgto, nder, maxco, dokp, locks, natom) &
 !$OMP PRIVATE (ikind, jkind, iatom, jatom, rab, basis_set_a, basis_set_b, &
 !$OMP          first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
 !$OMP          zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
 !$OMP          zetb, dab, irow, icol, h_block, found, iset, ncoa, &
 !$OMP          sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, hab2, hab2_w, qab, &
-!$OMP          h1_1block, h1_2block, h1_3block, kkind, nseta, &
+!$OMP          atk0, atk1, h1_1block, h1_2block, h1_3block, kkind, nseta, &
 !$OMP          gth_potential, sgp_potential, alpha, cexp_ppl, lpotextended, &
 !$OMP          ppl_radius, nexp_lpot, nexp_ppl, alpha_ppl, alpha_lpot, nct_ppl, &
 !$OMP          nct_lpot, cval_ppl, cval_lpot, rac, dac, rbc, dbc, &
@@ -240,7 +248,7 @@ CONTAINS
 !$OMP          nloc, nrloc, aloc, bloc, n_local, a_local, c_local, &
 !$OMP          slmax, npot, nrpot, apot, bpot, only_gaussians, &
 !$OMP          ldai, hash, hash1, hash2, iatom8) &
-!$OMP REDUCTION (+ : pv_thread, force_thread )
+!$OMP REDUCTION (+ : pv_thread, force_thread, at_thread )
 
 !$OMP SINGLE
 !$    ALLOCATE (locks(nlock))
@@ -258,7 +266,7 @@ CONTAINS
       ALLOCATE (hab(ldsab, ldsab, maxnset, maxnset), work(ldsab, ldsab*maxder))
       ldai = ncoset(2*maxlgto + 2*nder)
       ALLOCATE (ppl_work(ldai, ldai, MAX(maxder, 2*maxlgto + 2*nder + 1)))
-      IF (calculate_forces) THEN
+      IF (calculate_forces .OR. doat) THEN
          ALLOCATE (pab(maxco, maxco, maxnset, maxnset))
          ldai = ncoset(maxlgto)
          ALLOCATE (ppl_fwork(ldai, ldai, maxder))
@@ -350,7 +358,7 @@ CONTAINS
 
          CALL dbcsr_get_block_p(matrix_h(1, img)%matrix, irow, icol, h_block, found)
          CPASSERT(found)
-         IF (calculate_forces) THEN
+         IF (calculate_forces .OR. doat) THEN
             NULLIFY (p_block)
             CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
             IF (ASSOCIATED(p_block)) THEN
@@ -461,6 +469,10 @@ CONTAINS
                      sgfb = first_sgfb(1, jset)
                      IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
                      ! *** Calculate the GTH pseudo potential forces ***
+                     IF (doat) THEN
+                        atk0 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
+                                      pab(1:ncoa, 1:ncob, iset, jset))
+                     END IF
                      IF (calculate_forces) THEN
 
                         force_a(:) = 0.0_dp
@@ -579,6 +591,12 @@ CONTAINS
 !$OMP END CRITICAL(type2)
                         END IF
                      END IF
+                     ! calculate atomic contributions
+                     IF (doat) THEN
+                        atk1 = f0*SUM(hab(1:ncoa, 1:ncob, iset, jset)* &
+                                      pab(1:ncoa, 1:ncob, iset, jset))
+                        at_thread(katom) = at_thread(katom) + (atk1 - atk0)
+                     END IF
                   END DO
                END DO
             END DO
@@ -690,7 +708,7 @@ CONTAINS
 
       DEALLOCATE (basis_set_list)
 
-      IF (calculate_forces) THEN
+      IF (calculate_forces .OR. doat) THEN
          ! *** If LSD, then recover alpha density and beta density     ***
          ! *** from the total density (1) and the spin density (2)     ***
          IF (SIZE(matrix_p, 1) == 2) THEN
@@ -714,6 +732,9 @@ CONTAINS
 !$OMP END DO
          DEALLOCATE (atom_of_kind, kind_of)
       END IF
+      IF (doat) THEN
+         atcore(1:natom) = atcore(1:natom) + at_thread(1:natom)
+      END IF
 
       IF (calculate_forces .AND. use_virial) THEN
          virial%pv_ppl = virial%pv_ppl + pv_thread
diff --git a/src/core_ppnl.F b/src/core_ppnl.F
index 222e38e199..c4dd13ae04 100644
--- a/src/core_ppnl.F
+++ b/src/core_ppnl.F
@@ -82,10 +82,11 @@ CONTAINS
 !> \param basis_type ...
 !> \param deltaR Weighting factors of the derivatives wrt. nuclear positions
 !> \param matrix_l ...
+!> \param atcore ...
 ! **************************************************************************************************
    SUBROUTINE build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
                               qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
-                              nimages, cell_to_index, basis_type, deltaR, matrix_l)
+                              nimages, cell_to_index, basis_type, deltaR, matrix_l, atcore)
 
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p
       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
@@ -106,6 +107,8 @@ CONTAINS
          OPTIONAL                                        :: deltaR
       TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
          POINTER                                         :: matrix_l
+      REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
+         OPTIONAL                                        :: atcore
 
       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'build_core_ppnl'
 
@@ -119,9 +122,9 @@ CONTAINS
       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, npgfa, nprj_ppnl, &
                                                             nsgf_seta
       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa
-      LOGICAL                                            :: do_dR, do_gth, do_kp, do_soc, found, &
-                                                            ppnl_present
-      REAL(KIND=dp)                                      :: dac, f0, ppnl_radius
+      LOGICAL                                            :: do_dR, do_gth, do_kp, do_soc, doat, &
+                                                            found, ppnl_present
+      REAL(KIND=dp)                                      :: atk, dac, f0, ppnl_radius
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: radp
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: sab, work
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: ai_work, lab, work_l
@@ -134,6 +137,7 @@ CONTAINS
       TYPE(gth_potential_p_type), DIMENSION(:), POINTER  :: gpotential
       TYPE(clist_type), POINTER                          :: clist
       TYPE(alist_type), POINTER                          :: alist_ac, alist_bc
+      REAL(KIND=dp), DIMENSION(SIZE(particle_set))       :: at_thread
       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: achint, acint, alkint, bchint, bcint, &
                                                             blkint
       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cprj, h_block, l_block_x, l_block_y, &
@@ -156,6 +160,8 @@ CONTAINS
 
       do_dR = .FALSE.
       IF (PRESENT(deltaR)) do_dR = .TRUE.
+      doat = .FALSE.
+      IF (PRESENT(atcore)) doat = .TRUE.
 
       IF (calculate_forces) THEN
          CALL timeset(routineN//"_forces", handle)
@@ -178,7 +184,7 @@ CONTAINS
             CPASSERT(PRESENT(cell_to_index) .AND. ASSOCIATED(cell_to_index))
          END IF
 
-         IF (calculate_forces) THEN
+         IF (calculate_forces .OR. doat) THEN
             IF (SIZE(matrix_p, 1) == 2) THEN
                DO img = 1, nimages
                   CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
@@ -478,23 +484,24 @@ CONTAINS
          ! We now calculate the Hamiltonian matrix elements
 
          force_thread = 0.0_dp
+         at_thread = 0.0_dp
          pv_thread = 0.0_dp
 
 !$OMP PARALLEL &
 !$OMP DEFAULT (NONE) &
 !$OMP SHARED  (do_kp, basis_set, matrix_h, matrix_l, cell_to_index,&
 !$OMP          sab_orb, matrix_p, sap_int, nkind, eps_ppnl, force, &
-!$OMP          do_dR, deltaR, maxder, nder, &
+!$OMP          doat, do_dR, deltaR, maxder, nder, &
 !$OMP          locks, virial, use_virial, calculate_forces, do_soc, natom) &
 !$OMP PRIVATE (ikind, jkind, iatom, jatom, cell_b, rab, &
 !$OMP          slot, iab, atom_a, f0, irow, icol, h_block, &
 !$OMP          l_block_x, l_block_y, l_block_z, &
-!$OMP          r_2block, r_3block, &
+!$OMP          r_2block, r_3block, atk, &
 !$OMP          found,p_block, iac, ibc, alist_ac, alist_bc, acint, bcint, &
 !$OMP          achint, bchint, alkint, blkint, &
 !$OMP          na, np, nb, katom, j, fa, fb, rbc, rac, &
 !$OMP          kkind, kac, kbc, i, img, hash, iatom8) &
-!$OMP REDUCTION (+ : pv_thread, force_thread )
+!$OMP REDUCTION (+ : at_thread, pv_thread, force_thread )
 
 !$OMP SINGLE
 !$       ALLOCATE (locks(nlock))
@@ -557,7 +564,7 @@ CONTAINS
                CALL dbcsr_get_block_p(matrix_h(3, img)%matrix, irow, icol, r_3block, found)
             END IF
 
-            IF (calculate_forces) THEN
+            IF (calculate_forces .OR. doat) THEN
                NULLIFY (p_block)
                CALL dbcsr_get_block_p(matrix_p(1, img)%matrix, irow, icol, p_block, found)
             END IF
@@ -711,6 +718,19 @@ CONTAINS
                               END IF
 
                            END IF
+                           IF (doat) THEN
+                              IF (ASSOCIATED(p_block)) THEN
+                                 katom = alist_ac%clist(kac)%catom
+                                 IF (iatom <= jatom) THEN
+                                    atk = SUM(p_block(1:na, 1:nb)* &
+                                              MATMUL(achint(1:na, 1:np, 1), TRANSPOSE(bcint(1:nb, 1:np, 1))))
+                                 ELSE
+                                    atk = SUM(p_block(1:nb, 1:na)* &
+                                              MATMUL(bchint(1:nb, 1:np, 1), TRANSPOSE(acint(1:na, 1:np, 1))))
+                                 END IF
+                                 at_thread(katom) = at_thread(katom) + f0*atk
+                              END IF
+                           END IF
                            EXIT ! We have found a match and there can be only one single match
                         END IF
                      END DO
@@ -751,7 +771,11 @@ CONTAINS
             virial%pv_virial = virial%pv_virial + pv_thread
          END IF
 
-         IF (calculate_forces) THEN
+         IF (doat) THEN
+            atcore(1:natom) = atcore(1:natom) + at_thread
+         END IF
+
+         IF (calculate_forces .OR. doat) THEN
             ! If LSD, then recover alpha density and beta density
             ! from the total density (1) and the spin density (2)
             IF (SIZE(matrix_p, 1) == 2) THEN
diff --git a/src/ed_analysis.F b/src/ed_analysis.F
index 3c4fdcfb63..694fc807cf 100644
--- a/src/ed_analysis.F
+++ b/src/ed_analysis.F
@@ -71,6 +71,7 @@ MODULE ed_analysis
    USE qs_core_energies,                ONLY: calculate_ecore_alpha,&
                                               calculate_ecore_overlap,&
                                               calculate_ecore_self
+   USE qs_core_hamiltonian,             ONLY: core_matrices
    USE qs_dispersion_pairpot,           ONLY: calculate_dispersion_pairpot
    USE qs_dispersion_types,             ONLY: qs_dispersion_type
    USE qs_energy_types,                 ONLY: qs_energy_type
@@ -137,8 +138,8 @@ CONTAINS
       REAL(KIND=dp)                                      :: ateps, checksum, e1, e2, e_pot, ealpha, &
                                                             ecc, egcp, ehfx, ekts, evdw, focc, &
                                                             sum_energy
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: amval, ate1h, ate1xc, atecc, ateks, &
-                                                            atener, atewald, odiag
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: amval, atcore, ate1h, ate1xc, atecc, &
+                                                            ateks, atener, atewald, odiag
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: atdet, atmul, mcharge, mweight
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: bcenter
       REAL(KIND=dp), DIMENSION(:), POINTER               :: occupation_numbers
@@ -151,8 +152,9 @@ CONTAINS
       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
       TYPE(dbcsr_p_type)                                 :: dve_mat
       TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:)      :: exc_mat, ks_mat, vhxc_mat
-      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_h, matrix_hfx, matrix_ks, &
-                                                            matrix_p, matrix_s, matrix_t
+      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: core_mat, matrix_h, matrix_hfx, &
+                                                            matrix_ks, matrix_p, matrix_s, matrix_t
+      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: math, matp
       TYPE(dbcsr_type)                                   :: dkmat, dmat
       TYPE(dbcsr_type), POINTER                          :: smat
       TYPE(dft_control_type), POINTER                    :: dft_control
@@ -436,11 +438,12 @@ CONTAINS
       CALL cp_fm_release(cvec2)
 
       ! energy arrays
-      ALLOCATE (atener(natom), ateks(natom), atecc(natom))
+      ALLOCATE (atener(natom), ateks(natom), atecc(natom), atcore(natom))
       ALLOCATE (ate1xc(natom), ate1h(natom), atewald(natom))
       atener = 0.0_dp
       ateks = 0.0_dp
       atecc = 0.0_dp
+      atcore = 0.0_dp
       ate1xc = 0.0_dp
       ate1h = 0.0_dp
       atewald = 0.0_dp
@@ -455,7 +458,7 @@ CONTAINS
       CALL dbcsr_set(dkmat, 0.0_dp)
       ! KS matrix + correction
       CALL get_qs_env(qs_env, matrix_h=matrix_h, matrix_ks=matrix_ks, kinetic=matrix_t)
-      ALLOCATE (ks_mat(nspin), vhxc_mat(nspin), exc_mat(1))
+      ALLOCATE (ks_mat(nspin), core_mat(1), vhxc_mat(nspin), exc_mat(1))
       DO ispin = 1, nspin
          ALLOCATE (ks_mat(ispin)%matrix)
          CALL dbcsr_create(ks_mat(ispin)%matrix, template=matrix_h(1)%matrix)
@@ -468,6 +471,12 @@ CONTAINS
          CALL dbcsr_copy(vhxc_mat(ispin)%matrix, smat)
          CALL dbcsr_set(vhxc_mat(ispin)%matrix, 0.0_dp)
       END DO
+      !
+      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)
+      !
       ALLOCATE (exc_mat(1)%matrix)
       CALL dbcsr_create(exc_mat(1)%matrix, template=smat)
       CALL dbcsr_copy(exc_mat(1)%matrix, smat)
@@ -475,6 +484,12 @@ CONTAINS
       !
       CALL vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
       !
+      CALL get_qs_env(qs_env, rho=rho)
+      CALL qs_rho_get(rho, rho_ao=matrix_p)
+      math(1:1, 1:1) => core_mat(1:1)
+      matp(1:nspin, 1:1) => matrix_p(1:nspin)
+      CALL core_matrices(qs_env, math, matp, .FALSE., 0, atcore=atcore)
+      !
       IF (ewald_correction) THEN
          ALLOCATE (dve_mat%matrix)
          CALL dbcsr_create(dve_mat%matrix, template=matrix_h(1)%matrix)
@@ -485,6 +500,7 @@ CONTAINS
       !
       DO ispin = 1, nspin
          CALL dbcsr_add(ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix, 1.0_dp, 1.0_dp)
+         CALL dbcsr_add(ks_mat(ispin)%matrix, core_mat(1)%matrix, 1.0_dp, -0.5_dp)
       END DO
       !
       IF (detailed_ener .AND. do_hfx) THEN
@@ -575,9 +591,11 @@ CONTAINS
          DEALLOCATE (ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix)
          CALL deallocate_mo_set(mos_loc(ispin))
       END DO
+      CALL dbcsr_release(core_mat(1)%matrix)
+      DEALLOCATE (core_mat(1)%matrix)
       CALL dbcsr_release(exc_mat(1)%matrix)
       DEALLOCATE (exc_mat(1)%matrix)
-      DEALLOCATE (ks_mat, vhxc_mat, exc_mat)
+      DEALLOCATE (ks_mat, core_mat, vhxc_mat, exc_mat)
       DEALLOCATE (mos_loc)
       DEALLOCATE (refbas_blk_sizes)
       DEALLOCATE (ref_basis_set_list)
@@ -604,6 +622,9 @@ CONTAINS
       group = para_env
       ! KS energy
       atener(1:natom) = ateks(1:natom)
+      ! 1/2 of VPP contribution Tr[VPP(K)*P]
+      CALL group%sum(atcore)  ! maybe not needed
+      atener(1:natom) = atener(1:natom) + 0.5_dp*atcore(1:natom)
       ! core energy corrections
       CALL group%sum(atecc)
       atener(1:natom) = atener(1:natom) + atecc(1:natom)
@@ -651,7 +672,9 @@ CONTAINS
       CALL get_qs_env(qs_env, energy=energy)
       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
          IF (unit_nr > 0) THEN
             WRITE (unit_nr, FMT="(/,T2,A)") "Detailed IAO Atomic Energy Components "
@@ -726,7 +749,7 @@ CONTAINS
             "!--------------------------- END OF ED ANALYSIS ------------------------------!"
       END IF
       DEALLOCATE (bcenter)
-      DEALLOCATE (atener, ateks, atecc, ate1xc, ate1h, atewald)
+      DEALLOCATE (atener, ateks, atecc, atcore, ate1xc, ate1h, atewald)
 
       CALL timestop(handle)
 
diff --git a/src/qs_core_hamiltonian.F b/src/qs_core_hamiltonian.F
index 0d6834ddc9..5a6654db70 100644
--- a/src/qs_core_hamiltonian.F
+++ b/src/qs_core_hamiltonian.F
@@ -140,7 +140,7 @@ MODULE qs_core_hamiltonian
 
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian'
 
-   PUBLIC :: build_core_hamiltonian_matrix
+   PUBLIC :: build_core_hamiltonian_matrix, core_matrices
    PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real
 
 CONTAINS
@@ -412,15 +412,17 @@ CONTAINS
 !> \param matrix_p ...
 !> \param calculate_forces ...
 !> \param nder ...
+!> \param atcore ...
 ! **************************************************************************************************
-   SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
+   SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, atcore)
 
       TYPE(qs_environment_type), POINTER                 :: qs_env
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p
       LOGICAL, INTENT(IN)                                :: calculate_forces
       INTEGER, INTENT(IN)                                :: nder
+      REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atcore
 
-      INTEGER                                            :: nimages
+      INTEGER                                            :: natom, nimages
       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
       LOGICAL                                            :: all_present, my_gt_nl, ppl_present, &
                                                             ppnl_present, use_virial
@@ -438,8 +440,11 @@ CONTAINS
       TYPE(virial_type), POINTER                         :: virial
 
       NULLIFY (dft_control)
-      CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control)
+      CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, natom=natom)
       nimages = dft_control%nimages
+      IF (PRESENT(atcore)) THEN
+         CPASSERT(SIZE(atcore) >= natom)
+      END IF
 
       ! check whether a gauge transformed version of the non-local potential part has to be used
       my_gt_nl = .FALSE.
@@ -484,7 +489,8 @@ CONTAINS
       all_present = ASSOCIATED(sac_ae)
       IF (all_present) THEN
          CALL build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
-                            qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
+                            qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
+                            nimages, cell_to_index, atcore=atcore)
       END IF
 
       ! *** compute the ppl contribution to the core hamiltonian ***
@@ -500,12 +506,12 @@ CONTAINS
                ELSE
                   CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
                                       qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
-                                      nimages, cell_to_index, "ORB")
+                                      nimages, cell_to_index, "ORB", atcore=atcore)
                END IF
             ELSE
                CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
                                    qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
-                                   nimages, cell_to_index, "ORB")
+                                   nimages, cell_to_index, "ORB", atcore=atcore)
             END IF
          END IF
       END IF
@@ -517,7 +523,7 @@ CONTAINS
          IF (.NOT. my_gt_nl) THEN
             CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
                                  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
-                                 nimages, cell_to_index, "ORB")
+                                 nimages, cell_to_index, "ORB", atcore=atcore)
          END IF
       END IF