Energy decomposition analysis (#2906)

src/CMakeLists.txt

@@ -128,6 +128,7 @@ list(
   ec_env_types.F
   ec_methods.F
   ec_orth_solver.F
+  ed_analysis.F
   efield_tb_methods.F
   efield_utils.F
   eip_environment.F

src/common/bibliography.F

@@ -90,7 +90,7 @@ MODULE bibliography
                             Behler2007, Behler2011, Schran2020a, Schran2020b, &
                             Rycroft2009, Thomas2015, Brehm2018, Brehm2020, Shigeta2001, Heinecke2016, &
                             Brehm2021, Bussy2021a, Bussy2021b, Ditler2021, Ditler2022, Mattiat2019, Mattiat2022, &
-                            Belleflamme2023, Knizia2013, Musaelian2023
+                            Belleflamme2023, Knizia2013, Musaelian2023, Eriksen2020
 
 CONTAINS
 
@@ -4841,6 +4841,18 @@ SUBROUTINE add_all_references()
                          "ER"), &
                          DOI="10.1038/s41467-023-36329-y")
 
+      CALL add_reference(key=Eriksen2020, ISI_record=s2a( &
+                         "TY JOUR", &
+                         "PT J", &
+                         "AU Eriksen, JJ", &
+                         "TI Mean-Field density matrix decompositions", &
+                         "SO The Journal of Chemical Physics", &
+                         "PY 2020", &
+                         "VL 153", &
+                         "AR 214109", &
+                         "ER"), &
+                         DOI="10.1063/5.0030764")
+
    END SUBROUTINE add_all_references
 
 END MODULE bibliography

src/ewald_methods_tb.F

@@ -94,7 +94,7 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
       LOGICAL, INTENT(in)                                :: calculate_forces
       TYPE(virial_type), POINTER                         :: virial
       LOGICAL, INTENT(in)                                :: use_virial
-      TYPE(atprop_type), POINTER                         :: atprop
+      TYPE(atprop_type), OPTIONAL, POINTER               :: atprop
 
       CHARACTER(len=*), PARAMETER                        :: routineN = 'tb_spme_evaluate'
 
@@ -196,65 +196,67 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
       CALL rs_grid_set_box(grid_spme, rs=rpot)
 
       ! Atomic Stress
-      IF (ASSOCIATED(atprop)) THEN
-         IF (atprop%stress .AND. use_virial) THEN
-
-            CALL rs_grid_zero(rpot)
-            rhob_g%cc = phi_g%cc*green%p3m_charge%cr
-            CALL pw_transfer(rhob_g, rhob_r)
-            CALL rs_pw_transfer(rpot, rhob_r, pw2rs)
-            ipart = 0
-            DO
-               CALL set_list(particle_set, npart, center, p1, rden, ipart)
-               IF (p1 == 0) EXIT
-               ! calculate function on small boxes
-               CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
-                              is_shell=.FALSE., unit_charge=.TRUE.)
-               CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
-               atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*mcharge(p1)*fint*dvols
-               atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*mcharge(p1)*fint*dvols
-               atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*mcharge(p1)*fint*dvols
-            END DO
-
-            NULLIFY (phib_g)
-            ALLOCATE (phib_g)
-            CALL pw_pool_create_pw(pw_pool, phib_g, &
-                                   use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
-            ffa = (0.5_dp/alpha)**2
-            ffb = 1.0_dp/fourpi
-            DO i = 1, 3
-               DO ig = grid_spme%first_gne0, grid_spme%ngpts_cut_local
-                  phib_g%cc(ig) = ffb*dphi_g(i)%cc(ig)*(1.0_dp + ffa*grid_spme%gsq(ig))
-                  phib_g%cc(ig) = phib_g%cc(ig)*green%influence_fn%cc(ig)
+      IF (PRESENT(atprop)) THEN
+         IF (ASSOCIATED(atprop)) THEN
+            IF (atprop%stress .AND. use_virial) THEN
+
+               CALL rs_grid_zero(rpot)
+               rhob_g%cc = phi_g%cc*green%p3m_charge%cr
+               CALL pw_transfer(rhob_g, rhob_r)
+               CALL rs_pw_transfer(rpot, rhob_r, pw2rs)
+               ipart = 0
+               DO
+                  CALL set_list(particle_set, npart, center, p1, rden, ipart)
+                  IF (p1 == 0) EXIT
+                  ! calculate function on small boxes
+                  CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
+                                 is_shell=.FALSE., unit_charge=.TRUE.)
+                  CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
+                  atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*mcharge(p1)*fint*dvols
+                  atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*mcharge(p1)*fint*dvols
+                  atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*mcharge(p1)*fint*dvols
                END DO
-               IF (grid_spme%have_g0) phib_g%cc(1) = 0.0_dp
-               DO j = 1, i
-                  nd = 0
-                  nd(j) = 1
-                  CALL pw_copy(phib_g, rhob_g)
-                  CALL pw_derive(rhob_g, nd)
-                  rhob_g%cc = rhob_g%cc*green%p3m_charge%cr
-                  CALL pw_transfer(rhob_g, rhob_r)
-                  CALL rs_pw_transfer(rpot, rhob_r, pw2rs)
-
-                  ipart = 0
-                  DO
-                     CALL set_list(particle_set, npart, center, p1, rden, ipart)
-                     IF (p1 == 0) EXIT
-                     ! calculate function on small boxes
-                     CALL get_patch(particle_set, delta, green, p1, rhos, &
-                                    is_core=.FALSE., is_shell=.FALSE., unit_charge=.TRUE.)
-                     ! integrate box and potential
-                     CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
-                     atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fint*dvols*mcharge(p1)
-                     IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fint*dvols*mcharge(p1)
+
+               NULLIFY (phib_g)
+               ALLOCATE (phib_g)
+               CALL pw_pool_create_pw(pw_pool, phib_g, &
+                                      use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE)
+               ffa = (0.5_dp/alpha)**2
+               ffb = 1.0_dp/fourpi
+               DO i = 1, 3
+                  DO ig = grid_spme%first_gne0, grid_spme%ngpts_cut_local
+                     phib_g%cc(ig) = ffb*dphi_g(i)%cc(ig)*(1.0_dp + ffa*grid_spme%gsq(ig))
+                     phib_g%cc(ig) = phib_g%cc(ig)*green%influence_fn%cc(ig)
                   END DO
+                  IF (grid_spme%have_g0) phib_g%cc(1) = 0.0_dp
+                  DO j = 1, i
+                     nd = 0
+                     nd(j) = 1
+                     CALL pw_copy(phib_g, rhob_g)
+                     CALL pw_derive(rhob_g, nd)
+                     rhob_g%cc = rhob_g%cc*green%p3m_charge%cr
+                     CALL pw_transfer(rhob_g, rhob_r)
+                     CALL rs_pw_transfer(rpot, rhob_r, pw2rs)
+
+                     ipart = 0
+                     DO
+                        CALL set_list(particle_set, npart, center, p1, rden, ipart)
+                        IF (p1 == 0) EXIT
+                        ! calculate function on small boxes
+                        CALL get_patch(particle_set, delta, green, p1, rhos, &
+                                       is_core=.FALSE., is_shell=.FALSE., unit_charge=.TRUE.)
+                        ! integrate box and potential
+                        CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
+                        atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fint*dvols*mcharge(p1)
+                        IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fint*dvols*mcharge(p1)
+                     END DO
 
+                  END DO
                END DO
-            END DO
-            CALL pw_pool_give_back_pw(pw_pool, phib_g)
-            DEALLOCATE (phib_g)
+               CALL pw_pool_give_back_pw(pw_pool, phib_g)
+               DEALLOCATE (phib_g)
 
+            END IF
          END IF
       END IF
 
@@ -367,7 +369,7 @@ SUBROUTINE tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial
          POINTER                                         :: n_list
       TYPE(virial_type), POINTER                         :: virial
       LOGICAL, INTENT(IN)                                :: use_virial
-      TYPE(atprop_type), POINTER                         :: atprop
+      TYPE(atprop_type), OPTIONAL, POINTER               :: atprop
 
       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'tb_ewald_overlap'
 
@@ -404,10 +406,12 @@ SUBROUTINE tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial
                   pfr = -dfr*mcharge(iatom)*mcharge(jatom)
                END IF
                CALL virial_pair_force(virial%pv_virial, -pfr, rij, rij)
-               IF (ASSOCIATED(atprop)) THEN
-                  IF (atprop%stress) THEN
-                     CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*pfr, rij, rij)
-                     CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp*pfr, rij, rij)
+               IF (PRESENT(atprop)) THEN
+                  IF (ASSOCIATED(atprop)) THEN
+                     IF (atprop%stress) THEN
+                        CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp*pfr, rij, rij)
+                        CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp*pfr, rij, rij)
+                     END IF
                   END IF
                END IF
             END IF

src/force_env_methods.F

@@ -69,6 +69,7 @@ MODULE force_env_methods
         use_eip_force, use_embed, use_fist_force, use_mixed_force, use_nnp_force, use_prog_name, &
         use_pwdft_force, use_qmmm, use_qmmmx, use_qs_force
    USE force_env_utils,                 ONLY: rescale_forces,&
+                                              write_atener,&
                                               write_forces
    USE force_fields_util,               ONLY: get_generic_info
    USE fp_methods,                      ONLY: fp_eval
@@ -413,10 +414,9 @@ RECURSIVE SUBROUTINE force_env_calc_energy_force(force_env, calc_force, &
          CALL force_env%para_env%sum(atprop_env%atener)
          CALL force_env_get(force_env, potential_energy=e_pot)
          IF (output_unit > 0) THEN
-            IF (logger%iter_info%print_level > low_print_level) THEN
-               WRITE (UNIT=output_unit, FMT="(/,T6,A,T15,A)") "Atom", "Potential energy"
-               WRITE (UNIT=output_unit, FMT="(T2,I8,1X,F20.10)") &
-                  (i, atprop_env%atener(i), i=1, SIZE(atprop_env%atener))
+            IF (logger%iter_info%print_level >= low_print_level) THEN
+               CALL cp_subsys_get(subsys=subsys, particles=particles)
+               CALL write_atener(output_unit, particles, atprop_env%atener, "Mulliken Atomic Energies")
             END IF
             sum_energy = accurate_sum(atprop_env%atener(:))
             checksum = ABS(e_pot - sum_energy)

src/force_env_utils.F

@@ -31,6 +31,7 @@ MODULE force_env_utils
    USE molecule_types,                  ONLY: global_constraint_type
    USE particle_list_types,             ONLY: particle_list_type
    USE particle_types,                  ONLY: update_particle_set
+   USE physcon,                         ONLY: angstrom
 #include "./base/base_uses.f90"
 
    IMPLICIT NONE
@@ -42,6 +43,7 @@ MODULE force_env_utils
    PUBLIC :: force_env_shake, &
              force_env_rattle, &
              rescale_forces, &
+             write_atener, &
              write_forces
 
 CONTAINS
@@ -475,4 +477,38 @@ SUBROUTINE write_forces(particles, iw, label, ndigits, total_force, &
 
    END SUBROUTINE write_forces
 
+! **************************************************************************************************
+!> \brief Write the atomic coordinates, types, and energies
+!> \param iounit ...
+!> \param particles ...
+!> \param atener ...
+!> \param label ...
+!> \date    05.06.2023
+!> \author  JGH
+!> \version 1.0
+! **************************************************************************************************
+   SUBROUTINE write_atener(iounit, particles, atener, label)
+
+      INTEGER, INTENT(IN)                                :: iounit
+      TYPE(particle_list_type), POINTER                  :: particles
+      REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: atener
+      CHARACTER(LEN=*), INTENT(IN)                       :: label
+
+      INTEGER                                            :: i, natom
+
+      IF (iounit > 0) THEN
+         WRITE (UNIT=iounit, FMT="(/,T2,A)") TRIM(label)
+         WRITE (UNIT=iounit, FMT="(T4,A,T30,A,T42,A,T54,A,T69,A)") &
+            "Atom  Element", "X", "Y", "Z", "Energy[a.u.]"
+         natom = particles%n_els
+         DO i = 1, natom
+            WRITE (iounit, "(I6,T12,A2,T24,3F12.6,F21.10)") i, &
+               TRIM(ADJUSTL(particles%els(i)%atomic_kind%element_symbol)), &
+               particles%els(i)%r(1:3)*angstrom, atener(i)
+         END DO
+         WRITE (iounit, "(A)") ""
+      END IF
+
+   END SUBROUTINE write_atener
+
 END MODULE force_env_utils

src/hartree_local_methods.F

@@ -8,8 +8,6 @@
 MODULE hartree_local_methods
    USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                               get_atomic_kind
-   USE atprop_types,                    ONLY: atprop_array_init,&
-                                              atprop_type
    USE basis_set_types,                 ONLY: get_gto_basis_set,&
                                               gto_basis_set_type
    USE cell_types,                      ONLY: cell_type
@@ -225,11 +223,11 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
 
       INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, lmax0, &
          lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_s_harm, maxl, maxso, mepos, n1, nat, &
-         natom, nchan_0, nkind, nr, nset, nsotot, nspins, num_pe
+         nchan_0, nkind, nr, nset, nsotot, nspins, num_pe
       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: cg_n_list
       INTEGER, ALLOCATABLE, DIMENSION(:, :, :)           :: cg_list
       INTEGER, DIMENSION(:), POINTER                     :: atom_list, lmax, lmin, npgf
-      LOGICAL                                            :: atenergy, core_charge, l_2nd_local_rho, &
+      LOGICAL                                            :: core_charge, l_2nd_local_rho, &
                                                             my_core_2nd, my_periodic, paw_atom
       REAL(dp)                                           :: back_ch, factor
       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: gexp, sqrtwr
@@ -239,7 +237,6 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
                                                             Vh1_s, vrrad_0, zet
       REAL(dp), DIMENSION(:, :, :), POINTER              :: my_CG, Qlm_gg, Qlm_gg_2nd
       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
-      TYPE(atprop_type), POINTER                         :: atprop
       TYPE(cell_type), POINTER                           :: cell
       TYPE(dft_control_type), POINTER                    :: dft_control
       TYPE(grid_atom_type), POINTER                      :: grid_atom
@@ -262,12 +259,12 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
       NULLIFY (rho0_mpole, rhoz_set)
       NULLIFY (atom_list, grid_atom, harmonics)
       NULLIFY (basis_1c, lmin, lmax, npgf, zet)
-      NULLIFY (gsph, atprop)
+      NULLIFY (gsph)
 
       CALL get_qs_env(qs_env=qs_env, &
                       cell=cell, dft_control=dft_control, &
                       atomic_kind_set=atomic_kind_set, &
-                      qs_kind_set=qs_kind_set, atprop=atprop, &
+                      qs_kind_set=qs_kind_set, &
                       pw_env=pw_env, qs_charges=qs_charges)
 
       CALL pw_env_get(pw_env, poisson_env=poisson_env)
@@ -305,15 +302,6 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
       CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso)
       CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax_0)
 
-      atenergy = .FALSE.
-      IF (ASSOCIATED(atprop)) THEN
-         atenergy = atprop%energy
-         IF (atenergy) THEN
-            CALL get_qs_env(qs_env=qs_env, natom=natom)
-            CALL atprop_array_init(atprop%ate1c, natom)
-         END IF
-      END IF
-
       !   Put to 0 the local hartree energy contribution from 1 center integrals
       energy_hartree_1c = 0.0_dp
 
@@ -434,7 +422,7 @@ SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set
 
                CALL Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
                                       grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, & ! core charge for density
-                                      nchan_0, nspins, max_iso_not0, atenergy, atprop%ate1c)
+                                      nchan_0, nspins, max_iso_not0)
 
                IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)
 
@@ -549,12 +537,10 @@ END SUBROUTINE Vh_1c_atom_potential
 !> \param nchan_0 ...
 !> \param nspins ...
 !> \param max_iso_not0 ...
-!> \param atenergy ...
-!> \param ate1c ...
 ! **************************************************************************************************
    SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
                                 grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, &
-                                nchan_0, nspins, max_iso_not0, atenergy, ate1c)
+                                nchan_0, nspins, max_iso_not0)
 
       REAL(dp), INTENT(INOUT)                            :: energy_hartree_1c
       TYPE(ecoul_1center_type), DIMENSION(:), POINTER    :: ecoul_1c
@@ -566,8 +552,6 @@ SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
       REAL(dp), DIMENSION(:), POINTER                    :: rrad_z
       REAL(dp), DIMENSION(:, :), POINTER                 :: Vh1_h, Vh1_s
       INTEGER, INTENT(IN)                                :: nchan_0, nspins, max_iso_not0
-      LOGICAL, INTENT(IN)                                :: atenergy
-      REAL(dp), DIMENSION(:), POINTER                    :: ate1c
 
       INTEGER                                            :: iso, ispin, nr
       REAL(dp)                                           :: ecoul_1_0, ecoul_1_h, ecoul_1_s, &
@@ -607,11 +591,6 @@ SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
       energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0
       energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s
 
-      ! atomic energy contribution
-      IF (atenergy) THEN
-         ate1c(iatom) = ate1c(iatom) + ecoul_1_z - ecoul_1_0
-      END IF
-
    END SUBROUTINE Vh_1c_atom_energy
 
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/hartree_local_types.F

@@ -24,8 +24,10 @@ MODULE hartree_local_types
 ! **************************************************************************************************
    TYPE ecoul_1center_type
       TYPE(rho_atom_coeff), POINTER   :: Vh1_h, Vh1_s
-      REAL(dp)                        :: ecoul_1_h, ecoul_1_s, &
-                                         ecoul_1_z, ecoul_1_0
+      REAL(dp)                        :: ecoul_1_h = 0.0_dp, &
+                                         ecoul_1_s = 0.0_dp, &
+                                         ecoul_1_z = 0.0_dp, &
+                                         ecoul_1_0 = 0.0_dp
    END TYPE ecoul_1center_type
 
 ! **************************************************************************************************