LCOV - code coverage report
Current view: top level - src - ed_analysis.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.3 % 598 576
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculate Energy Decomposition analysis
      10              : !> \par History
      11              : !>      07.2023 created [JGH]
      12              : !> \author JGH
      13              : ! **************************************************************************************************
      14              : MODULE ed_analysis
      15              :    USE admm_types,                      ONLY: admm_type
      16              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      17              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      18              :                                               gto_basis_set_type
      19              :    USE bibliography,                    ONLY: Eriksen2020,&
      20              :                                               cite_reference
      21              :    USE cell_types,                      ONLY: cell_type
      22              :    USE core_ae,                         ONLY: build_erfc
      23              :    USE cp_control_types,                ONLY: dft_control_type
      24              :    USE cp_dbcsr_api,                    ONLY: &
      25              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_info, &
      26              :         dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_symmetric
      27              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot,&
      28              :                                               dbcsr_reserve_diag_blocks
      29              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_plus_fm_fm_t,&
      30              :                                               cp_dbcsr_sm_fm_multiply
      31              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      32              :                                               cp_fm_struct_release,&
      33              :                                               cp_fm_struct_type
      34              :    USE cp_fm_types,                     ONLY: &
      35              :         cp_fm_create, cp_fm_get_diag, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, &
      36              :         cp_fm_release, cp_fm_set_all, cp_fm_set_element, cp_fm_set_submatrix, cp_fm_to_fm, &
      37              :         cp_fm_type
      38              :    USE hartree_local_methods,           ONLY: Vh_1c_gg_integrals
      39              :    USE hartree_local_types,             ONLY: ecoul_1center_type
      40              :    USE hfx_admm_utils,                  ONLY: tddft_hfx_matrix
      41              :    USE iao_analysis,                    ONLY: iao_calculate_dmat,&
      42              :                                               iao_charges,&
      43              :                                               iao_wfn_analysis,&
      44              :                                               print_center_spread
      45              :    USE iao_types,                       ONLY: iao_env_type,&
      46              :                                               iao_set_default
      47              :    USE input_constants,                 ONLY: do_admm_aux_exch_func_none
      48              :    USE input_section_types,             ONLY: section_vals_get,&
      49              :                                               section_vals_get_subs_vals,&
      50              :                                               section_vals_type,&
      51              :                                               section_vals_val_get
      52              :    USE kinds,                           ONLY: dp
      53              :    USE mathconstants,                   ONLY: pi
      54              :    USE message_passing,                 ONLY: mp_comm_type,&
      55              :                                               mp_para_env_type
      56              :    USE mulliken,                        ONLY: mulliken_charges
      57              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      58              :    USE particle_methods,                ONLY: get_particle_set
      59              :    USE particle_types,                  ONLY: particle_type
      60              :    USE physcon,                         ONLY: angstrom
      61              :    USE pw_env_types,                    ONLY: pw_env_get,&
      62              :                                               pw_env_type
      63              :    USE pw_methods,                      ONLY: pw_axpy,&
      64              :                                               pw_scale,&
      65              :                                               pw_transfer
      66              :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      67              :    USE pw_poisson_types,                ONLY: pw_poisson_type
      68              :    USE pw_pool_types,                   ONLY: pw_pool_type
      69              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      70              :                                               pw_r3d_rs_type
      71              :    USE qs_collocate_density,            ONLY: calculate_rho_core
      72              :    USE qs_core_energies,                ONLY: calculate_ecore_alpha,&
      73              :                                               calculate_ecore_overlap,&
      74              :                                               calculate_ecore_self
      75              :    USE qs_core_hamiltonian,             ONLY: core_matrices
      76              :    USE qs_dispersion_pairpot,           ONLY: calculate_dispersion_pairpot
      77              :    USE qs_dispersion_types,             ONLY: qs_dispersion_type
      78              :    USE qs_energy_types,                 ONLY: qs_energy_type
      79              :    USE qs_environment_types,            ONLY: get_qs_env,&
      80              :                                               qs_environment_type
      81              :    USE qs_gcp_method,                   ONLY: calculate_gcp_pairpot
      82              :    USE qs_gcp_types,                    ONLY: qs_gcp_type
      83              :    USE qs_integrate_potential,          ONLY: integrate_v_core_rspace,&
      84              :                                               integrate_v_gaussian_rspace,&
      85              :                                               integrate_v_rspace
      86              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      87              :                                               qs_kind_type
      88              :    USE qs_ks_atom,                      ONLY: update_ks_atom
      89              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      90              :    USE qs_local_rho_types,              ONLY: local_rho_type
      91              :    USE qs_mo_methods,                   ONLY: calculate_subspace_eigenvalues
      92              :    USE qs_mo_types,                     ONLY: deallocate_mo_set,&
      93              :                                               duplicate_mo_set,&
      94              :                                               get_mo_set,&
      95              :                                               mo_set_type
      96              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      97              :    USE qs_rho0_ggrid,                   ONLY: integrate_vhg0_rspace
      98              :    USE qs_rho_atom_types,               ONLY: rho_atom_type,&
      99              :                                               zero_rho_atom_integrals
     100              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
     101              :                                               qs_rho_type
     102              :    USE qs_vxc,                          ONLY: qs_xc_density
     103              :    USE qs_vxc_atom,                     ONLY: calculate_vxc_atom
     104              :    USE xc_derivatives,                  ONLY: xc_functionals_get_needs
     105              :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
     106              : #include "./base/base_uses.f90"
     107              : 
     108              :    IMPLICIT NONE
     109              :    PRIVATE
     110              : 
     111              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ed_analysis'
     112              : 
     113              :    PUBLIC ::  edmf_analysis
     114              : 
     115              : ! **************************************************************************************************
     116              : 
     117              : CONTAINS
     118              : 
     119              : ! **************************************************************************************************
     120              : !> \brief ...
     121              : !> \param qs_env ...
     122              : !> \param input_section ...
     123              : !> \param unit_nr ...
     124              : ! **************************************************************************************************
     125           58 :    SUBROUTINE edmf_analysis(qs_env, input_section, unit_nr)
     126              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     127              :       TYPE(section_vals_type), POINTER                   :: input_section
     128              :       INTEGER, INTENT(IN)                                :: unit_nr
     129              : 
     130              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'edmf_analysis'
     131              : 
     132              :       INTEGER                                            :: handle, iatom, ikind, iorb, ispin, jorb, &
     133              :                                                             nao, natom, nimages, nkind, no, norb, &
     134              :                                                             nref, nspin
     135              :       INTEGER, DIMENSION(2)                              :: nocc
     136           58 :       INTEGER, DIMENSION(:), POINTER                     :: refbas_blk_sizes
     137              :       LOGICAL :: detailed_ener, do_hfx, ewald_correction, explicit, gapw, gapw_xc, &
     138              :          ref_orb_canonical, skip_localize, uniform_occupation, uocc
     139              :       REAL(KIND=dp)                                      :: ateps, checksum, e1, e2, e_pot, ealpha, &
     140              :                                                             ecc, egcp, ehfx, ekts, evdw, focc, &
     141              :                                                             sum_energy
     142           58 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: amval, atcore, ate1h, ate1xc, atecc, &
     143           58 :                                                             ateks, atener, atewald, odiag
     144           58 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: atdet, atmul, mcharge, mweight
     145           58 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: bcenter
     146           58 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: occupation_numbers
     147              :       TYPE(cell_type), POINTER                           :: cell
     148              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     149              :       TYPE(cp_fm_type)                                   :: cvec, cvec2, smo
     150           58 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: c_iao_coef, ciao, fij_mat, orb_weight, &
     151           58 :                                                             rotmat
     152              :       TYPE(cp_fm_type), POINTER                          :: cloc, moref
     153              :       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
     154              :       TYPE(dbcsr_p_type)                                 :: dve_mat
     155           58 :       TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:)      :: exc_mat, ks_mat, vhxc_mat
     156           58 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: core_mat, matrix_h, matrix_hfx, &
     157           58 :                                                             matrix_ks, matrix_p, matrix_s, matrix_t
     158           58 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: math, matp
     159              :       TYPE(dbcsr_type)                                   :: dkmat, dmat
     160              :       TYPE(dbcsr_type), POINTER                          :: smat
     161              :       TYPE(dft_control_type), POINTER                    :: dft_control
     162           58 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: ref_basis_set_list
     163              :       TYPE(gto_basis_set_type), POINTER                  :: refbasis
     164              :       TYPE(iao_env_type)                                 :: iao_env
     165           58 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos, mos_loc
     166              :       TYPE(mp_comm_type)                                 :: group
     167              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     168           58 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     169              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     170              :       TYPE(qs_energy_type), POINTER                      :: energy
     171              :       TYPE(qs_gcp_type), POINTER                         :: gcp_env
     172           58 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     173              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     174              :       TYPE(qs_rho_type), POINTER                         :: rho
     175              :       TYPE(section_vals_type), POINTER                   :: hfx_sections, input
     176              : 
     177           58 :       CALL section_vals_get(input_section, explicit=explicit)
     178           58 :       IF (.NOT. explicit) RETURN
     179              : 
     180           30 :       CALL timeset(routineN, handle)
     181              : 
     182           30 :       IF (unit_nr > 0) THEN
     183              :          WRITE (UNIT=unit_nr, FMT="(/,4(T2,A))") &
     184           15 :             "!-----------------------------------------------------------------------------!", &
     185           15 :             "!                        ENERGY DECOMPOSITION ANALYSIS                        !", &
     186           15 :             "!                    Janus J Eriksen, JCP 153 214109 (2020)                   !", &
     187           30 :             "!-----------------------------------------------------------------------------!"
     188              :       END IF
     189           30 :       CALL cite_reference(Eriksen2020)
     190              : 
     191              :       ! input options
     192           30 :       CALL section_vals_val_get(input_section, "REFERENCE_ORB_CANONICAL", l_val=ref_orb_canonical)
     193           30 :       CALL section_vals_val_get(input_section, "DETAILED_ENERGY", l_val=detailed_ener)
     194           30 :       CALL section_vals_val_get(input_section, "SKIP_LOCALIZATION", l_val=skip_localize)
     195           30 :       CALL section_vals_val_get(input_section, "EWALD_ALPHA_PARAMETER", r_val=ealpha)
     196           30 :       ewald_correction = (ealpha /= 0.0_dp)
     197              :       ! k-points?
     198           30 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     199           30 :       nimages = dft_control%nimages
     200           30 :       IF (nimages > 1) THEN
     201            0 :          IF (unit_nr > 0) THEN
     202              :             WRITE (UNIT=unit_nr, FMT="(T2,A)") &
     203            0 :                "K-Points: MAO's determined and analyzed using Gamma-Point only."
     204              :          END IF
     205              :       END IF
     206           30 :       gapw = dft_control%qs_control%gapw
     207           30 :       gapw_xc = dft_control%qs_control%gapw_xc
     208           30 :       IF (ewald_correction) THEN
     209            2 :          IF (gapw .OR. gapw_xc) THEN
     210            0 :             CALL cp_warn(__LOCATION__, "Ewald Correction for GAPW and GAPW_XC not available.")
     211            0 :             CPABORT("Ewald Correction")
     212              :          END IF
     213              :       END IF
     214           30 :       CALL get_qs_env(qs_env, input=input)
     215           30 :       hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
     216           30 :       CALL section_vals_get(hfx_sections, explicit=do_hfx)
     217              : 
     218           30 :       CALL get_qs_env(qs_env, mos=mos)
     219           30 :       nspin = dft_control%nspins
     220          122 :       ALLOCATE (mos_loc(nspin))
     221              : 
     222              :       ! do we have a uniform occupation
     223           30 :       uocc = .TRUE.
     224           62 :       DO ispin = 1, nspin
     225           32 :          CALL get_mo_set(mos(ispin), uniform_occupation=uniform_occupation)
     226           62 :          IF (.NOT. uniform_occupation) uocc = .FALSE.
     227              :       END DO
     228           30 :       IF (unit_nr > 0) THEN
     229           15 :          IF (uocc) THEN
     230              :             WRITE (UNIT=unit_nr, FMT="(T2,A)") &
     231           14 :                "MO's have a uniform occupation pattern"
     232              :          ELSE
     233              :             WRITE (UNIT=unit_nr, FMT="(T2,A)") &
     234            1 :                "MO's have varying occupations"
     235              :          END IF
     236              :       END IF
     237              : 
     238              :       ! perform IAO analysis
     239           30 :       CALL iao_set_default(iao_env)
     240           30 :       iao_env%do_iao = .TRUE.
     241           30 :       iao_env%do_charges = .TRUE.
     242           30 :       IF (skip_localize) THEN
     243            2 :          iao_env%do_bondorbitals = .FALSE.
     244            2 :          iao_env%do_center = .FALSE.
     245              :       ELSE
     246           28 :          iao_env%do_bondorbitals = .TRUE.
     247           28 :          iao_env%do_center = .TRUE.
     248              :       END IF
     249           30 :       iao_env%eps_occ = 1.0E-4_dp
     250           30 :       CALL get_qs_env(qs_env, cell=cell)
     251           36 :       iao_env%pos_periodic = .NOT. ALL(cell%perd == 0)
     252           30 :       no = 0
     253           30 :       nocc = 0
     254           62 :       DO ispin = 1, nspin
     255           32 :          CALL duplicate_mo_set(mos_loc(ispin), mos(ispin))
     256           32 :          CALL get_mo_set(mos_loc(ispin), nmo=norb)
     257           32 :          no = MAX(no, norb)
     258           32 :          nocc(ispin) = norb
     259           62 :          IF (ref_orb_canonical) THEN
     260           32 :             CALL get_mo_set(mos_loc(ispin), mo_coeff=moref)
     261           32 :             CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
     262              :             CALL calculate_subspace_eigenvalues(moref, matrix_ks(ispin)%matrix, &
     263           32 :                                                 do_rotation=.TRUE.)
     264              :          END IF
     265              :       END DO
     266          120 :       ALLOCATE (bcenter(5, no, nspin))
     267         1154 :       bcenter = 0.0_dp
     268              :       CALL iao_wfn_analysis(qs_env, iao_env, unit_nr, &
     269           30 :                             c_iao_coef=c_iao_coef, mos=mos_loc, bond_centers=bcenter)
     270              : 
     271              :       ! output bond centers
     272           30 :       IF (iao_env%do_bondorbitals .AND. iao_env%do_center) THEN
     273           28 :          CALL print_center_spread(bcenter, nocc, iounit=unit_nr)
     274              :       END IF
     275              : 
     276              :       ! Calculate orbital rotation matrix
     277           30 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     278           30 :       smat => matrix_s(1)%matrix
     279          122 :       ALLOCATE (rotmat(nspin))
     280           62 :       DO ispin = 1, nspin
     281           32 :          CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
     282           32 :          CALL get_mo_set(mos(ispin), mo_coeff=moref)
     283           32 :          CALL cp_fm_get_info(cloc, nrow_global=nao, ncol_global=norb)
     284           32 :          CALL cp_fm_create(smo, cloc%matrix_struct)
     285           32 :          NULLIFY (fm_struct)
     286              :          CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
     287           32 :                                   template_fmstruct=cloc%matrix_struct)
     288           32 :          CALL cp_fm_create(rotmat(ispin), fm_struct)
     289           32 :          CALL cp_fm_struct_release(fm_struct)
     290              :          ! ROTMAT = Cref(T)*S*Cloc
     291           32 :          CALL cp_dbcsr_sm_fm_multiply(smat, cloc, smo, ncol=norb)
     292              :          CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, moref, &
     293           32 :                             smo, 0.0_dp, rotmat(ispin))
     294           94 :          CALL cp_fm_release(smo)
     295              :       END DO
     296              : 
     297              :       ! calculate occupation matrix
     298           30 :       IF (.NOT. uocc) THEN
     299            6 :          ALLOCATE (fij_mat(nspin))
     300            4 :          DO ispin = 1, nspin
     301            2 :             CALL cp_fm_create(fij_mat(ispin), rotmat(ispin)%matrix_struct)
     302            2 :             CALL cp_fm_set_all(fij_mat(ispin), 0.0_dp)
     303            2 :             CALL cp_fm_create(smo, rotmat(ispin)%matrix_struct)
     304              :             ! fii = f
     305            2 :             CALL get_mo_set(mos(ispin), nmo=norb, occupation_numbers=occupation_numbers)
     306           54 :             DO iorb = 1, norb
     307           54 :                CALL cp_fm_set_element(fij_mat(ispin), iorb, iorb, occupation_numbers(iorb))
     308              :             END DO
     309              :             ! fij = U(T)*f*U
     310              :             CALL parallel_gemm('N', 'N', norb, norb, norb, 1.0_dp, fij_mat(ispin), &
     311            2 :                                rotmat(ispin), 0.0_dp, smo)
     312              :             CALL parallel_gemm('T', 'N', norb, norb, norb, 1.0_dp, rotmat(ispin), &
     313            2 :                                smo, 0.0_dp, fij_mat(ispin))
     314            6 :             CALL cp_fm_release(smo)
     315              :          END DO
     316              :       END IF
     317              : 
     318              :       ! localized orbitals in IAO basis => CIAO
     319           92 :       ALLOCATE (ciao(nspin))
     320           62 :       DO ispin = 1, nspin
     321           32 :          CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
     322           32 :          CALL cp_fm_get_info(cloc, nrow_global=nao, ncol_global=norb)
     323           32 :          CALL cp_fm_get_info(c_iao_coef(ispin), ncol_global=nref)
     324           32 :          CALL cp_fm_create(smo, cloc%matrix_struct)
     325           32 :          NULLIFY (fm_struct)
     326              :          CALL cp_fm_struct_create(fm_struct, nrow_global=nref, &
     327           32 :                                   template_fmstruct=cloc%matrix_struct)
     328           32 :          CALL cp_fm_create(ciao(ispin), fm_struct)
     329           32 :          CALL cp_fm_struct_release(fm_struct)
     330              :          ! CIAO = A(T)*S*C
     331           32 :          CALL cp_dbcsr_sm_fm_multiply(smat, cloc, smo, ncol=norb)
     332              :          CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, c_iao_coef(ispin), &
     333           32 :                             smo, 0.0_dp, ciao(ispin))
     334           94 :          CALL cp_fm_release(smo)
     335              :       END DO
     336              : 
     337              :       ! Reference basis set
     338           30 :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
     339           30 :       nkind = SIZE(qs_kind_set)
     340          148 :       ALLOCATE (ref_basis_set_list(nkind))
     341           88 :       DO ikind = 1, nkind
     342           58 :          qs_kind => qs_kind_set(ikind)
     343           58 :          NULLIFY (ref_basis_set_list(ikind)%gto_basis_set)
     344           58 :          NULLIFY (refbasis)
     345           58 :          CALL get_qs_kind(qs_kind=qs_kind, basis_set=refbasis, basis_type="MIN")
     346           88 :          IF (ASSOCIATED(refbasis)) ref_basis_set_list(ikind)%gto_basis_set => refbasis
     347              :       END DO
     348           30 :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, natom=natom)
     349           90 :       ALLOCATE (refbas_blk_sizes(natom))
     350              :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=refbas_blk_sizes, &
     351           30 :                             basis=ref_basis_set_list)
     352              : 
     353              :       ! Atomic orbital weights
     354           92 :       ALLOCATE (orb_weight(nspin))
     355           90 :       ALLOCATE (mcharge(natom, 1))
     356           62 :       DO ispin = 1, nspin
     357           32 :          CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc)
     358           32 :          NULLIFY (fm_struct)
     359              :          CALL cp_fm_struct_create(fm_struct, nrow_global=natom, &
     360           32 :                                   template_fmstruct=cloc%matrix_struct)
     361           32 :          CALL cp_fm_create(orb_weight(ispin), fm_struct)
     362           32 :          CALL cp_fm_struct_release(fm_struct)
     363           62 :          CALL cp_fm_set_all(orb_weight(ispin), 0.0_dp)
     364              :       END DO
     365              :       !
     366           30 :       CALL dbcsr_get_info(smat, distribution=dbcsr_dist)
     367              :       CALL dbcsr_create(matrix=dmat, name="DMAT", &
     368              :                         dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
     369           30 :                         row_blk_size=refbas_blk_sizes, col_blk_size=refbas_blk_sizes)
     370           30 :       CALL dbcsr_reserve_diag_blocks(dmat)
     371              :       !
     372           30 :       NULLIFY (fm_struct)
     373              :       CALL cp_fm_struct_create(fm_struct, ncol_global=1, &
     374           30 :                                template_fmstruct=ciao(1)%matrix_struct)
     375           30 :       CALL cp_fm_create(cvec, fm_struct)
     376           30 :       CALL cp_fm_create(cvec2, fm_struct)
     377           30 :       CALL cp_fm_struct_release(fm_struct)
     378              :       !
     379           62 :       DO ispin = 1, nspin
     380              :          CALL get_mo_set(mos_loc(ispin), &
     381              :                          mo_coeff=cloc, nmo=norb, &
     382           32 :                          occupation_numbers=occupation_numbers)
     383           62 :          IF (uocc) THEN
     384              :             ! uniform occupation
     385          158 :             DO iorb = 1, norb
     386          128 :                CALL cp_fm_to_fm(ciao(ispin), cvec, ncol=1, source_start=iorb, target_start=1)
     387          128 :                focc = occupation_numbers(iorb)
     388          128 :                CALL dbcsr_set(dmat, 0.0_dp)
     389          128 :                CALL cp_dbcsr_plus_fm_fm_t(dmat, cvec, ncol=1, alpha=focc)
     390          128 :                CALL iao_charges(dmat, mcharge(:, 1))
     391              :                CALL cp_fm_set_submatrix(orb_weight(ispin), mcharge, start_row=1, &
     392          128 :                                         start_col=iorb, n_rows=natom, n_cols=1)
     393          560 :                checksum = SUM(mcharge(:, 1))
     394          158 :                IF (ABS(focc - checksum) > 1.E-6_dp) THEN
     395            0 :                   CALL cp_warn(__LOCATION__, "Sum of atomic orbital weights is incorrect")
     396            0 :                   IF (unit_nr > 0) THEN
     397              :                      WRITE (UNIT=unit_nr, FMT="(T2,A,F10.6,T40,A,F10.6)") &
     398            0 :                         "Orbital occupation:", focc, &
     399            0 :                         "Sum of atomic orbital weights:", checksum
     400              :                   END IF
     401              :                END IF
     402              :             END DO
     403              :          ELSE
     404              :             ! non-diagonal occupation matrix
     405            6 :             ALLOCATE (odiag(norb))
     406            2 :             CALL cp_fm_get_diag(fij_mat(ispin), odiag)
     407           54 :             DO iorb = 1, norb
     408           52 :                IF (odiag(iorb) < 1.E-8_dp) CYCLE
     409           44 :                CALL dbcsr_set(dmat, 0.0_dp)
     410           44 :                CALL cp_fm_to_fm(ciao(ispin), cvec, ncol=1, source_start=iorb, target_start=1)
     411         1188 :                DO jorb = 1, norb
     412         1144 :                   CALL cp_fm_get_element(fij_mat(ispin), iorb, jorb, focc)
     413         1144 :                   IF (focc < 1.E-12_dp) CYCLE
     414          486 :                   CALL cp_fm_to_fm(ciao(ispin), cvec2, ncol=1, source_start=jorb, target_start=1)
     415         1674 :                   CALL cp_dbcsr_plus_fm_fm_t(dmat, cvec, cvec2, 1, alpha=focc, symmetry_mode=1)
     416              :                END DO
     417           44 :                CALL iao_charges(dmat, mcharge(:, 1))
     418          396 :                checksum = SUM(mcharge(:, 1))
     419           44 :                focc = odiag(iorb)
     420           44 :                IF (ABS(focc - checksum) > 1.E-6_dp) THEN
     421            0 :                   CALL cp_warn(__LOCATION__, "Sum of atomic orbital weights is incorrect")
     422            0 :                   IF (unit_nr > 0) THEN
     423              :                      WRITE (UNIT=unit_nr, FMT="(T2,A,F10.6,T40,A,F10.6)") &
     424            0 :                         "Orbital occupation:", focc, &
     425            0 :                         "Sum of atomic orbital weights:", checksum
     426              :                   END IF
     427              :                END IF
     428          396 :                mcharge(:, 1) = mcharge(:, 1)/focc
     429              :                CALL cp_fm_set_submatrix(orb_weight(ispin), mcharge, start_row=1, &
     430           54 :                                         start_col=iorb, n_rows=natom, n_cols=1)
     431              :             END DO
     432            2 :             DEALLOCATE (odiag)
     433              :          END IF
     434              :       END DO
     435           30 :       DEALLOCATE (mcharge)
     436           30 :       CALL dbcsr_release(dmat)
     437           30 :       CALL cp_fm_release(cvec)
     438           30 :       CALL cp_fm_release(cvec2)
     439              : 
     440           30 :       CALL get_qs_env(qs_env, para_env=para_env)
     441           30 :       group = para_env
     442              :       ! energy arrays
     443          180 :       ALLOCATE (atener(natom), ateks(natom), atecc(natom), atcore(natom))
     444          120 :       ALLOCATE (ate1xc(natom), ate1h(natom), atewald(natom))
     445          136 :       atener = 0.0_dp
     446          136 :       ateks = 0.0_dp
     447          136 :       atecc = 0.0_dp
     448          136 :       atcore = 0.0_dp
     449          136 :       ate1xc = 0.0_dp
     450          136 :       ate1h = 0.0_dp
     451          136 :       atewald = 0.0_dp
     452           30 :       IF (detailed_ener) THEN
     453           30 :          ALLOCATE (atdet(natom, 10), atmul(natom, 10), amval(natom))
     454          246 :          atdet = 0.0_dp
     455          246 :          atmul = 0.0_dp
     456              :       END IF
     457              :       ! atom dependent density matrix
     458           30 :       CALL dbcsr_create(dkmat, template=smat)
     459           30 :       CALL dbcsr_copy(dkmat, smat)
     460           30 :       CALL dbcsr_set(dkmat, 0.0_dp)
     461              :       ! KS matrix + correction
     462           30 :       CALL get_qs_env(qs_env, matrix_h=matrix_h, matrix_ks=matrix_ks, kinetic=matrix_t)
     463          244 :       ALLOCATE (ks_mat(nspin), core_mat(1), vhxc_mat(nspin), exc_mat(1))
     464           62 :       DO ispin = 1, nspin
     465           32 :          ALLOCATE (ks_mat(ispin)%matrix)
     466           32 :          CALL dbcsr_create(ks_mat(ispin)%matrix, template=matrix_h(1)%matrix)
     467           32 :          CALL dbcsr_copy(ks_mat(ispin)%matrix, matrix_h(1)%matrix)
     468           32 :          CALL dbcsr_add(ks_mat(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, 1.0_dp)
     469           32 :          CALL dbcsr_scale(ks_mat(ispin)%matrix, 0.5_dp)
     470              :          !
     471           32 :          ALLOCATE (vhxc_mat(ispin)%matrix)
     472           32 :          CALL dbcsr_create(vhxc_mat(ispin)%matrix, template=smat)
     473           32 :          CALL dbcsr_copy(vhxc_mat(ispin)%matrix, smat)
     474           62 :          CALL dbcsr_set(vhxc_mat(ispin)%matrix, 0.0_dp)
     475              :       END DO
     476              :       !
     477           30 :       ALLOCATE (core_mat(1)%matrix)
     478           30 :       CALL dbcsr_create(core_mat(1)%matrix, template=matrix_h(1)%matrix)
     479           30 :       CALL dbcsr_copy(core_mat(1)%matrix, matrix_h(1)%matrix)
     480           30 :       CALL dbcsr_set(core_mat(1)%matrix, 0.0_dp)
     481              :       !
     482           30 :       ALLOCATE (exc_mat(1)%matrix)
     483           30 :       CALL dbcsr_create(exc_mat(1)%matrix, template=smat)
     484           30 :       CALL dbcsr_copy(exc_mat(1)%matrix, smat)
     485           30 :       CALL dbcsr_set(exc_mat(1)%matrix, 0.0_dp)
     486              :       !
     487           30 :       CALL vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
     488              :       !
     489           30 :       CALL get_qs_env(qs_env, rho=rho)
     490           30 :       CALL qs_rho_get(rho, rho_ao=matrix_p)
     491           30 :       math(1:1, 1:1) => core_mat(1:1)
     492           30 :       matp(1:nspin, 1:1) => matrix_p(1:nspin)
     493           30 :       CALL core_matrices(qs_env, math, matp, .FALSE., 0, atcore=atcore)
     494           30 :       CALL group%sum(atcore)
     495              :       !
     496           30 :       IF (ewald_correction) THEN
     497            2 :          ALLOCATE (dve_mat%matrix)
     498            2 :          CALL dbcsr_create(dve_mat%matrix, template=matrix_h(1)%matrix)
     499            2 :          CALL dbcsr_copy(dve_mat%matrix, matrix_h(1)%matrix)
     500            2 :          CALL dbcsr_set(dve_mat%matrix, 0.0_dp)
     501            2 :          CALL vh_ewald_correction(qs_env, ealpha, dve_mat, atewald)
     502            2 :          CALL group%sum(atewald)
     503              :       END IF
     504              :       !
     505           62 :       DO ispin = 1, nspin
     506           32 :          CALL dbcsr_add(ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix, 1.0_dp, 1.0_dp)
     507           62 :          CALL dbcsr_add(ks_mat(ispin)%matrix, core_mat(1)%matrix, 1.0_dp, -0.5_dp)
     508              :       END DO
     509              :       !
     510           30 :       IF (detailed_ener .AND. do_hfx) THEN
     511            6 :          ALLOCATE (matrix_hfx(nspin))
     512            4 :          DO ispin = 1, nspin
     513            2 :             ALLOCATE (matrix_hfx(nspin)%matrix)
     514            2 :             CALL dbcsr_create(matrix_hfx(ispin)%matrix, template=smat)
     515            2 :             CALL dbcsr_copy(matrix_hfx(ispin)%matrix, smat)
     516            4 :             CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
     517              :          END DO
     518            2 :          CALL get_hfx_ks_matrix(qs_env, matrix_hfx)
     519              :       END IF
     520              :       ! Loop over spins and atoms
     521           62 :       DO ispin = 1, nspin
     522           32 :          CALL get_mo_set(mos_loc(ispin), mo_coeff=cloc, nmo=norb)
     523           96 :          ALLOCATE (mweight(1, norb))
     524          144 :          DO iatom = 1, natom
     525              :             CALL cp_fm_get_submatrix(orb_weight(ispin), mweight, start_row=iatom, &
     526          112 :                                      start_col=1, n_rows=1, n_cols=norb)
     527          112 :             IF (uocc) THEN
     528           96 :                CALL iao_calculate_dmat(cloc, dkmat, mweight(1, :), .FALSE.)
     529              :             ELSE
     530           16 :                CALL iao_calculate_dmat(cloc, dkmat, mweight(1, :), fij_mat(ispin))
     531              :             END IF
     532          112 :             CALL dbcsr_dot(dkmat, ks_mat(ispin)%matrix, ecc)
     533          112 :             ateks(iatom) = ateks(iatom) + ecc
     534              :             !
     535          112 :             IF (detailed_ener) THEN
     536           18 :                CALL dbcsr_dot(dkmat, matrix_t(1)%matrix, e1)
     537           18 :                atdet(iatom, 1) = atdet(iatom, 1) + e1
     538           18 :                CALL dbcsr_dot(dkmat, matrix_h(1)%matrix, e2)
     539           18 :                atdet(iatom, 2) = atdet(iatom, 2) + (e2 - e1)
     540           18 :                IF (do_hfx) THEN
     541            6 :                   CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 0.5_dp)
     542            6 :                   CALL dbcsr_dot(dkmat, matrix_hfx(ispin)%matrix, ehfx)
     543            6 :                   atdet(iatom, 5) = atdet(iatom, 5) + ehfx
     544            6 :                   CALL dbcsr_scale(matrix_hfx(ispin)%matrix, 2.0_dp)
     545              :                ELSE
     546           12 :                   ehfx = 0.0_dp
     547              :                END IF
     548           18 :                CALL dbcsr_dot(dkmat, exc_mat(1)%matrix, e1)
     549           18 :                atdet(iatom, 3) = atdet(iatom, 3) + ecc - e2 - e1 - ehfx
     550           18 :                atdet(iatom, 4) = atdet(iatom, 4) + e1
     551              :             END IF
     552          256 :             IF (ewald_correction) THEN
     553            6 :                CALL dbcsr_dot(dkmat, dve_mat%matrix, ecc)
     554            6 :                atewald(iatom) = atewald(iatom) + ecc
     555              :             END IF
     556              :             !
     557              :          END DO
     558           94 :          DEALLOCATE (mweight)
     559              :       END DO
     560              :       !
     561           30 :       IF (detailed_ener) THEN
     562              :          ! Mulliken
     563            6 :          CALL get_qs_env(qs_env, rho=rho, para_env=para_env)
     564            6 :          CALL qs_rho_get(rho, rho_ao=matrix_p)
     565              :          ! kinetic energy
     566           12 :          DO ispin = 1, nspin
     567              :             CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_t(1)%matrix, &
     568            6 :                                   para_env, amval)
     569           24 :             atmul(1:natom, 1) = atmul(1:natom, 1) + amval(1:natom)
     570              :             CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_h(1)%matrix, &
     571            6 :                                   para_env, amval)
     572           24 :             atmul(1:natom, 2) = atmul(1:natom, 2) + amval(1:natom)
     573           24 :             atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
     574              :             CALL mulliken_charges(matrix_p(ispin)%matrix, ks_mat(ispin)%matrix, &
     575            6 :                                   para_env, amval)
     576           24 :             atmul(1:natom, 3) = atmul(1:natom, 3) + amval(1:natom)
     577            6 :             IF (do_hfx) THEN
     578              :                CALL mulliken_charges(matrix_p(ispin)%matrix, matrix_hfx(ispin)%matrix, &
     579            2 :                                      para_env, amval)
     580            8 :                atmul(1:natom, 5) = atmul(1:natom, 5) + 0.5_dp*amval(1:natom)
     581            8 :                atmul(1:natom, 3) = atmul(1:natom, 3) - 0.5_dp*amval(1:natom)
     582              :             END IF
     583              :             CALL mulliken_charges(matrix_p(ispin)%matrix, exc_mat(1)%matrix, &
     584            6 :                                   para_env, amval)
     585           24 :             atmul(1:natom, 4) = atmul(1:natom, 4) + amval(1:natom)
     586           30 :             atmul(1:natom, 3) = atmul(1:natom, 3) - amval(1:natom)
     587              :          END DO
     588           24 :          atmul(1:natom, 2) = atmul(1:natom, 2) - atmul(1:natom, 1)
     589           24 :          atmul(1:natom, 3) = atmul(1:natom, 3) + 0.5_dp*atmul(1:natom, 2)
     590           24 :          atmul(1:natom, 2) = 0.5_dp*atmul(1:natom, 2) + 0.5_dp*atcore(1:natom)
     591              :       END IF
     592              :       !
     593           30 :       CALL dbcsr_release(dkmat)
     594           62 :       DO ispin = 1, nspin
     595           32 :          CALL dbcsr_release(ks_mat(ispin)%matrix)
     596           32 :          CALL dbcsr_release(vhxc_mat(ispin)%matrix)
     597           32 :          DEALLOCATE (ks_mat(ispin)%matrix, vhxc_mat(ispin)%matrix)
     598           62 :          CALL deallocate_mo_set(mos_loc(ispin))
     599              :       END DO
     600           30 :       CALL dbcsr_release(core_mat(1)%matrix)
     601           30 :       DEALLOCATE (core_mat(1)%matrix)
     602           30 :       CALL dbcsr_release(exc_mat(1)%matrix)
     603           30 :       DEALLOCATE (exc_mat(1)%matrix)
     604           30 :       DEALLOCATE (ks_mat, core_mat, vhxc_mat, exc_mat)
     605           30 :       DEALLOCATE (mos_loc)
     606           30 :       DEALLOCATE (refbas_blk_sizes)
     607           30 :       DEALLOCATE (ref_basis_set_list)
     608           30 :       IF (detailed_ener .AND. do_hfx) THEN
     609            4 :          DO ispin = 1, nspin
     610            2 :             CALL dbcsr_release(matrix_hfx(ispin)%matrix)
     611            4 :             DEALLOCATE (matrix_hfx(ispin)%matrix)
     612              :          END DO
     613            2 :          DEALLOCATE (matrix_hfx)
     614              :       END IF
     615           30 :       IF (ewald_correction) THEN
     616            2 :          CALL dbcsr_release(dve_mat%matrix)
     617            2 :          DEALLOCATE (dve_mat%matrix)
     618              :       END IF
     619           30 :       CALL cp_fm_release(orb_weight)
     620           30 :       CALL cp_fm_release(ciao)
     621           30 :       CALL cp_fm_release(rotmat)
     622           30 :       CALL cp_fm_release(c_iao_coef)
     623           30 :       IF (.NOT. uocc) THEN
     624            2 :          CALL cp_fm_release(fij_mat)
     625              :       END IF
     626              : 
     627              :       ! KS energy
     628          136 :       atener(1:natom) = ateks(1:natom)
     629              :       ! 1/2 of VPP contribution Tr[VPP(K)*P]
     630          136 :       atener(1:natom) = atener(1:natom) + 0.5_dp*atcore(1:natom)
     631              :       ! core energy corrections
     632           30 :       CALL group%sum(atecc)
     633          136 :       atener(1:natom) = atener(1:natom) + atecc(1:natom)
     634           48 :       IF (detailed_ener) atdet(1:natom, 6) = atdet(1:natom, 6) + atecc(1:natom)
     635              :       ! one center terms (GAPW)
     636           30 :       CALL group%sum(ate1xc)
     637           30 :       CALL group%sum(ate1h)
     638          136 :       atener(1:natom) = atener(1:natom) + ate1xc(1:natom) + ate1h(1:natom)
     639           30 :       IF (detailed_ener) THEN
     640            6 :          IF (gapw .OR. gapw_xc) atdet(1:natom, 8) = atdet(1:natom, 8) + ate1xc(1:natom)
     641            0 :          IF (gapw) atdet(1:natom, 9) = atdet(1:natom, 9) + ate1h(1:natom)
     642              :       END IF
     643              :       ! core correction
     644          136 :       atecc(1:natom) = 0.0_dp
     645           30 :       CALL calculate_ecore_overlap(qs_env, para_env, .FALSE., atecc=atecc)
     646           30 :       CALL group%sum(atecc)
     647          136 :       atener(1:natom) = atener(1:natom) + atecc(1:natom)
     648           48 :       IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
     649           30 :       IF (ewald_correction) THEN
     650            8 :          atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
     651              :       END IF
     652              :       ! e self
     653          136 :       atecc(1:natom) = 0.0_dp
     654           30 :       CALL calculate_ecore_self(qs_env, atecc=atecc)
     655           30 :       CALL group%sum(atecc)
     656          136 :       atener(1:natom) = atener(1:natom) + atecc(1:natom)
     657           48 :       IF (detailed_ener) atdet(1:natom, 7) = atdet(1:natom, 7) + atecc(1:natom)
     658           30 :       IF (ewald_correction) THEN
     659            8 :          atewald(1:natom) = atewald(1:natom) - atecc(1:natom)
     660            2 :          CALL calculate_ecore_alpha(qs_env, ealpha, atewald)
     661              :       END IF
     662              :       ! vdW pair-potentials
     663          136 :       atecc(1:natom) = 0.0_dp
     664           30 :       CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
     665           30 :       CALL calculate_dispersion_pairpot(qs_env, dispersion_env, evdw, .FALSE., atevdw=atecc)
     666              :       ! Pair potential gCP energy
     667           30 :       CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env)
     668           30 :       IF (ASSOCIATED(gcp_env)) THEN
     669           30 :          CALL calculate_gcp_pairpot(qs_env, gcp_env, egcp, .FALSE., ategcp=atecc)
     670              :       END IF
     671           30 :       CALL group%sum(atecc)
     672          136 :       atener(1:natom) = atener(1:natom) + atecc(1:natom)
     673           48 :       IF (detailed_ener) atdet(1:natom, 10) = atdet(1:natom, 10) + atecc(1:natom)
     674              :       ! distribute the entropic energy
     675           30 :       CALL get_qs_env(qs_env, energy=energy)
     676           30 :       ekts = energy%kts/REAL(natom, KIND=dp)
     677          136 :       atener(1:natom) = atener(1:natom) + ekts
     678              :       ! 0.5 Vpp(at)*D + 0.5 * Vpp*D(at)
     679           30 :       IF (detailed_ener) THEN
     680           24 :          atdet(1:natom, 3) = atdet(1:natom, 3) + 0.5_dp*atdet(1:natom, 2)
     681           24 :          atdet(1:natom, 2) = 0.5_dp*atdet(1:natom, 2) + 0.5_dp*atcore(1:natom)
     682              :       END IF
     683              :       !
     684           30 :       IF (detailed_ener) THEN
     685            6 :          IF (unit_nr > 0) THEN
     686            3 :             WRITE (unit_nr, FMT="(/,T2,A)") "Detailed IAO Atomic Energy Components "
     687            3 :             CALL write_atdet(unit_nr, atdet(:, 1), "Kinetic Energy")
     688            3 :             CALL write_atdet(unit_nr, atdet(:, 2), "Short-Range Core and/or PP Energy")
     689            3 :             CALL write_atdet(unit_nr, atdet(:, 3), "Hartree Energy (long-ranged part)")
     690            3 :             CALL write_atdet(unit_nr, atdet(:, 5), "Exact Exchange Energy")
     691            3 :             CALL write_atdet(unit_nr, atdet(:, 4), "Exchange-Correlation Energy")
     692            3 :             CALL write_atdet(unit_nr, atdet(:, 6), "Atomic Core Hartree: Int(nc V(n+nc) dx")
     693            3 :             CALL write_atdet(unit_nr, atdet(:, 7), "Core Self Energy Corrections")
     694            3 :             IF (gapw) THEN
     695            0 :                CALL write_atdet(unit_nr, atdet(:, 9), "GAPW: 1-center Hartree Energy")
     696              :             END IF
     697            3 :             IF (gapw .OR. gapw_xc) THEN
     698            0 :                CALL write_atdet(unit_nr, atdet(:, 8), "GAPW: 1-center XC Energy")
     699              :             END IF
     700            3 :             IF (ABS(evdw) > 1.E-14_dp) THEN
     701            0 :                CALL write_atdet(unit_nr, atdet(:, 10), "vdW Pairpotential Energy")
     702              :             END IF
     703           12 :             DO iatom = 1, natom
     704          102 :                atecc(iatom) = SUM(atdet(iatom, 1:10)) - atener(iatom)
     705              :             END DO
     706            3 :             CALL write_atdet(unit_nr, atecc(:), "Missing Energy")
     707              :             !
     708            3 :             WRITE (unit_nr, FMT="(/,T2,A)") "Detailed Mulliken Atomic Energy Components "
     709            3 :             CALL write_atdet(unit_nr, atmul(:, 1), "Kinetic Energy")
     710            3 :             CALL write_atdet(unit_nr, atmul(:, 2), "Short-Range Core and/or PP Energy")
     711            3 :             CALL write_atdet(unit_nr, atmul(:, 3), "Hartree Energy (long-ranged part)")
     712            3 :             CALL write_atdet(unit_nr, atmul(:, 5), "Exact Exchange Energy")
     713            3 :             CALL write_atdet(unit_nr, atmul(:, 4), "Exchange-Correlation Energy")
     714            3 :             CALL write_atdet(unit_nr, atdet(:, 6), "Atomic Core Hartree: Int(nc V(n+nc) dx")
     715            3 :             CALL write_atdet(unit_nr, atdet(:, 7), "Core Self Energy Corrections")
     716            3 :             IF (gapw) THEN
     717            0 :                CALL write_atdet(unit_nr, atdet(:, 9), "GAPW: 1-center Hartree Energy")
     718              :             END IF
     719            3 :             IF (gapw .OR. gapw_xc) THEN
     720            0 :                CALL write_atdet(unit_nr, atdet(:, 8), "GAPW: 1-center XC Energy")
     721              :             END IF
     722            3 :             IF (ABS(evdw) > 1.E-14_dp) THEN
     723            0 :                CALL write_atdet(unit_nr, atdet(:, 10), "vdW Pairpotential Energy")
     724              :             END IF
     725              :          END IF
     726              :       END IF
     727              : 
     728           30 :       IF (unit_nr > 0) THEN
     729           15 :          e_pot = energy%total
     730           15 :          ateps = 1.E-6_dp
     731           15 :          CALL write_atener(unit_nr, particle_set, atener, "Atomic Energy Decomposition")
     732           68 :          sum_energy = SUM(atener(:))
     733           15 :          checksum = ABS(e_pot - sum_energy)
     734              :          WRITE (UNIT=unit_nr, FMT="((T2,A,T56,F25.13))") &
     735           15 :             "Potential energy (Atomic):", sum_energy, &
     736           15 :             "Potential energy (Total) :", e_pot, &
     737           30 :             "Difference               :", checksum
     738           15 :          CPASSERT((checksum < ateps*ABS(e_pot)))
     739              :          !
     740           15 :          IF (ewald_correction) THEN
     741            1 :             WRITE (UNIT=unit_nr, FMT="(/,(T2,A,F10.3,A))") "Universal Ewald Parameter: ", ealpha, " [Angstrom]"
     742            1 :             CALL write_atener(unit_nr, particle_set, atewald, "Change in Atomic Energy Decomposition")
     743            4 :             sum_energy = SUM(atewald(:))
     744              :             WRITE (UNIT=unit_nr, FMT="((T2,A,T56,F25.13))") &
     745            1 :                "Total Change in Potential energy:", sum_energy
     746              :          END IF
     747              :       END IF
     748              : 
     749           30 :       IF (detailed_ener) THEN
     750            6 :          DEALLOCATE (atdet, atmul, amval)
     751              :       END IF
     752              : 
     753           30 :       IF (unit_nr > 0) THEN
     754              :          WRITE (UNIT=unit_nr, FMT="(/,T2,A)") &
     755           15 :             "!--------------------------- END OF ED ANALYSIS ------------------------------!"
     756              :       END IF
     757           30 :       DEALLOCATE (bcenter)
     758           30 :       DEALLOCATE (atener, ateks, atecc, atcore, ate1xc, ate1h, atewald)
     759              : 
     760           30 :       CALL timestop(handle)
     761              : 
     762          266 :    END SUBROUTINE edmf_analysis
     763              : 
     764              : ! **************************************************************************************************
     765              : !> \brief ...
     766              : !> \param qs_env ...
     767              : !> \param vhxc_mat ...
     768              : !> \param exc_mat ...
     769              : !> \param atecc ...
     770              : !> \param ate1xc ...
     771              : !> \param ate1h ...
     772              : ! **************************************************************************************************
     773           30 :    SUBROUTINE vhxc_correction(qs_env, vhxc_mat, exc_mat, atecc, ate1xc, ate1h)
     774              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     775              :       TYPE(dbcsr_p_type), DIMENSION(:)                   :: vhxc_mat, exc_mat
     776              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: atecc, ate1xc, ate1h
     777              : 
     778              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'vhxc_correction'
     779              : 
     780              :       INTEGER                                            :: handle, iatom, ispin, natom, nspins
     781              :       LOGICAL                                            :: do_admm_corr, gapw, gapw_xc
     782              :       REAL(KIND=dp)                                      :: eh1, exc1
     783              :       TYPE(admm_type), POINTER                           :: admm_env
     784           30 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     785           30 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p
     786              :       TYPE(dft_control_type), POINTER                    :: dft_control
     787           30 :       TYPE(ecoul_1center_type), DIMENSION(:), POINTER    :: ecoul_1c
     788              :       TYPE(local_rho_type), POINTER                      :: local_rho_set
     789              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     790              :       TYPE(pw_env_type), POINTER                         :: pw_env
     791              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     792              :       TYPE(pw_r3d_rs_type)                               :: xc_den
     793           30 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: vtau, vxc
     794              :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
     795              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     796           30 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     797              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     798              :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     799           30 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
     800              :       TYPE(section_vals_type), POINTER                   :: xc_fun_section, xc_section
     801              :       TYPE(xc_rho_cflags_type)                           :: needs
     802              : 
     803           30 :       CALL timeset(routineN, handle)
     804              : 
     805           30 :       CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env)
     806              : 
     807           30 :       nspins = dft_control%nspins
     808           30 :       gapw = dft_control%qs_control%gapw
     809           30 :       gapw_xc = dft_control%qs_control%gapw_xc
     810           30 :       do_admm_corr = .FALSE.
     811           30 :       IF (dft_control%do_admm) THEN
     812            0 :          IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) do_admm_corr = .TRUE.
     813              :       END IF
     814              :       IF (do_admm_corr) THEN
     815            0 :          CALL get_qs_env(qs_env, admm_env=admm_env)
     816            0 :          xc_section => admm_env%xc_section_aux
     817              :       ELSE
     818           30 :          xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
     819              :       END IF
     820           30 :       xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     821           30 :       needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
     822              : 
     823           30 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     824           30 :       CALL auxbas_pw_pool%create_pw(xc_den)
     825          122 :       ALLOCATE (vxc(nspins))
     826           62 :       DO ispin = 1, nspins
     827           62 :          CALL auxbas_pw_pool%create_pw(vxc(ispin))
     828              :       END DO
     829           30 :       IF (needs%tau .OR. needs%tau_spin) THEN
     830           24 :          ALLOCATE (vtau(nspins))
     831           16 :          DO ispin = 1, nspins
     832           38 :             CALL auxbas_pw_pool%create_pw(vtau(ispin))
     833              :          END DO
     834              :       END IF
     835              : 
     836              :       ! Nuclear charge correction
     837           30 :       CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
     838           30 :       CALL integrate_v_core_rspace(v_hartree_rspace, qs_env, atecc=atecc)
     839           30 :       IF (gapw .OR. gapw_xc) THEN
     840              :          CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, &
     841              :                          rho_atom_set=rho_atom_set, ecoul_1c=ecoul_1c, &
     842            8 :                          natom=natom, para_env=para_env)
     843            8 :          CALL zero_rho_atom_integrals(rho_atom_set)
     844            8 :          CALL calculate_vxc_atom(qs_env, .FALSE., exc1)
     845            8 :          IF (gapw) THEN
     846            6 :             CALL Vh_1c_gg_integrals(qs_env, eh1, ecoul_1c, local_rho_set, para_env, tddft=.FALSE.)
     847            6 :             CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
     848              :             CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., &
     849            6 :                                        local_rho_set=local_rho_set, atener=atecc)
     850              :          END IF
     851              :       END IF
     852              : 
     853            8 :       IF (gapw_xc) THEN
     854            2 :          CALL get_qs_env(qs_env, rho_xc=rho_struct, dispersion_env=dispersion_env)
     855              :       ELSE
     856           28 :          CALL get_qs_env(qs_env, rho=rho_struct, dispersion_env=dispersion_env)
     857              :       END IF
     858              :       !
     859           30 :       CPASSERT(.NOT. do_admm_corr)
     860              :       !
     861           30 :       IF (needs%tau .OR. needs%tau_spin) THEN
     862              :          CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
     863            8 :                             xc_den=xc_den, vxc=vxc, vtau=vtau)
     864              :       ELSE
     865              :          CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
     866           22 :                             xc_den=xc_den, vxc=vxc)
     867              :       END IF
     868           62 :       DO ispin = 1, nspins
     869           32 :          CALL pw_scale(vxc(ispin), -0.5_dp)
     870           32 :          CALL pw_axpy(xc_den, vxc(ispin))
     871           32 :          CALL pw_scale(vxc(ispin), vxc(ispin)%pw_grid%dvol)
     872              :          CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vxc(ispin), &
     873              :                                  hmat=vhxc_mat(ispin), calculate_forces=.FALSE., &
     874           32 :                                  gapw=(gapw .OR. gapw_xc))
     875           62 :          IF (needs%tau .OR. needs%tau_spin) THEN
     876            8 :             CALL pw_scale(vtau(ispin), -0.5_dp*vtau(ispin)%pw_grid%dvol)
     877              :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vtau(ispin), &
     878              :                                     hmat=vhxc_mat(ispin), calculate_forces=.FALSE., &
     879            8 :                                     compute_tau=.TRUE., gapw=(gapw .OR. gapw_xc))
     880              :          END IF
     881              :       END DO
     882           30 :       CALL pw_scale(xc_den, xc_den%pw_grid%dvol)
     883              :       CALL integrate_v_rspace(qs_env=qs_env, v_rspace=xc_den, &
     884              :                               hmat=exc_mat(1), calculate_forces=.FALSE., &
     885           30 :                               gapw=(gapw .OR. gapw_xc))
     886              : 
     887           30 :       IF (gapw .OR. gapw_xc) THEN
     888              :          ! remove one-center potential matrix part
     889            8 :          CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
     890            8 :          CALL update_ks_atom(qs_env, vhxc_mat, matrix_p, forces=.FALSE., kscale=-0.5_dp)
     891              :          !
     892           32 :          DO iatom = 1, natom
     893              :             ate1xc(iatom) = ate1xc(iatom) + &
     894           32 :                             rho_atom_set(iatom)%exc_h - rho_atom_set(iatom)%exc_s
     895              :          END DO
     896            8 :          IF (gapw) THEN
     897           24 :             DO iatom = 1, natom
     898              :                ate1h(iatom) = ate1h(iatom) + &
     899              :                               ecoul_1c(iatom)%ecoul_1_h - ecoul_1c(iatom)%ecoul_1_s + &
     900           24 :                               ecoul_1c(iatom)%ecoul_1_z - ecoul_1c(iatom)%ecoul_1_0
     901              :             END DO
     902              :          END IF
     903              :       END IF
     904              : 
     905           30 :       CALL auxbas_pw_pool%give_back_pw(xc_den)
     906           62 :       DO ispin = 1, nspins
     907           62 :          CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
     908              :       END DO
     909           30 :       IF (needs%tau .OR. needs%tau_spin) THEN
     910           16 :          DO ispin = 1, nspins
     911           38 :             CALL auxbas_pw_pool%give_back_pw(vtau(ispin))
     912              :          END DO
     913              :       END IF
     914              : 
     915           30 :       CALL timestop(handle)
     916              : 
     917           60 :    END SUBROUTINE vhxc_correction
     918              : 
     919              : ! **************************************************************************************************
     920              : !> \brief ...
     921              : !> \param qs_env ...
     922              : !> \param ealpha ...
     923              : !> \param vh_mat ...
     924              : !> \param atewd ...
     925              : ! **************************************************************************************************
     926            2 :    SUBROUTINE vh_ewald_correction(qs_env, ealpha, vh_mat, atewd)
     927              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     928              :       REAL(KIND=dp), INTENT(IN)                          :: ealpha
     929              :       TYPE(dbcsr_p_type)                                 :: vh_mat
     930              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: atewd
     931              : 
     932              :       CHARACTER(len=*), PARAMETER :: routineN = 'vh_ewald_correction'
     933              : 
     934              :       INTEGER                                            :: handle, ikind, ispin, natom, nkind
     935              :       REAL(KIND=dp)                                      :: eps_core_charge, rhotot, zeff
     936              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: alpha, atcore, atecc, ccore
     937            2 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     938              :       TYPE(dbcsr_p_type)                                 :: e_mat
     939            2 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p
     940              :       TYPE(dft_control_type), POINTER                    :: dft_control
     941              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     942            2 :          POINTER                                         :: sab_orb, sac_ae
     943            2 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     944              :       TYPE(pw_c1d_gs_type)                               :: rho_tot_ewd_g, v_hewd_gspace
     945              :       TYPE(pw_env_type), POINTER                         :: pw_env
     946              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     947              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     948              :       TYPE(pw_r3d_rs_type)                               :: rho_tot_ewd_r, v_hewd_rspace
     949            2 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
     950              :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
     951            2 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     952              :       TYPE(qs_rho_type), POINTER                         :: rho
     953              : 
     954            2 :       CALL timeset(routineN, handle)
     955              : 
     956            2 :       natom = SIZE(atewd)
     957            8 :       ALLOCATE (atecc(natom), atcore(natom))
     958              :       CALL get_qs_env(qs_env=qs_env, nkind=nkind, qs_kind_set=qs_kind_set, &
     959            2 :                       dft_control=dft_control)
     960            8 :       ALLOCATE (alpha(nkind), ccore(nkind))
     961            6 :       DO ikind = 1, nkind
     962            4 :          CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
     963            4 :          alpha(ikind) = ealpha
     964            6 :          ccore(ikind) = zeff*(ealpha/pi)**1.5_dp
     965              :       END DO
     966              :       !
     967            2 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
     968              :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     969            2 :                       poisson_env=poisson_env)
     970              :       !
     971            2 :       CALL auxbas_pw_pool%create_pw(v_hewd_gspace)
     972            2 :       CALL auxbas_pw_pool%create_pw(v_hewd_rspace)
     973            2 :       CALL auxbas_pw_pool%create_pw(rho_tot_ewd_g)
     974            2 :       CALL auxbas_pw_pool%create_pw(rho_tot_ewd_r)
     975              :       rhotot = 0.0_dp
     976            2 :       CALL calculate_rho_core(rho_tot_ewd_r, rhotot, qs_env, calpha=alpha, ccore=ccore)
     977              :       ! Get the total density in g-space [ions + electrons]
     978            2 :       CALL get_qs_env(qs_env=qs_env, rho=rho)
     979            2 :       CALL qs_rho_get(rho, rho_r=rho_r)
     980            4 :       DO ispin = 1, dft_control%nspins
     981            4 :          CALL pw_axpy(rho_r(ispin), rho_tot_ewd_r)
     982              :       END DO
     983            2 :       CALL pw_transfer(rho_tot_ewd_r, rho_tot_ewd_g)
     984              :       !
     985            2 :       CALL pw_poisson_solve(poisson_env, density=rho_tot_ewd_g, vhartree=v_hewd_gspace)
     986            2 :       CALL pw_transfer(v_hewd_gspace, v_hewd_rspace)
     987            2 :       CALL pw_scale(v_hewd_rspace, v_hewd_rspace%pw_grid%dvol)
     988            8 :       atecc = 0.0_dp
     989            2 :       CALL integrate_v_gaussian_rspace(v_hewd_rspace, qs_env, alpha, ccore, atecc)
     990            8 :       atewd(1:natom) = atewd(1:natom) + atecc(1:natom)
     991              :       !
     992            8 :       atecc = 0.0_dp
     993            2 :       CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
     994            2 :       CALL integrate_v_core_rspace(v_hartree_rspace, qs_env, atecc=atecc)
     995            8 :       atewd(1:natom) = atewd(1:natom) - atecc(1:natom)
     996              :       !
     997            2 :       CALL pw_axpy(v_hartree_rspace, v_hewd_rspace, -1.0_dp)
     998            2 :       CALL dbcsr_set(vh_mat%matrix, 0.0_dp)
     999              :       CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hewd_rspace, hmat=vh_mat, &
    1000            2 :                               calculate_forces=.FALSE.)
    1001              :       !
    1002            2 :       CALL dbcsr_scale(vh_mat%matrix, 0.5_dp)
    1003              :       !
    1004              :       ! delta erfc
    1005            2 :       CALL qs_rho_get(rho, rho_ao=matrix_p)
    1006            2 :       eps_core_charge = dft_control%qs_control%eps_core_charge
    1007            2 :       ALLOCATE (e_mat%matrix)
    1008            2 :       CALL dbcsr_create(e_mat%matrix, template=vh_mat%matrix)
    1009            2 :       CALL dbcsr_copy(e_mat%matrix, vh_mat%matrix)
    1010              :       !
    1011              :       CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
    1012            2 :                       particle_set=particle_set, sab_orb=sab_orb, sac_ppl=sac_ae)
    1013            2 :       CALL dbcsr_set(e_mat%matrix, 0.0_dp)
    1014            8 :       atcore = 0.0_dp
    1015              :       CALL build_erfc(e_mat, matrix_p, qs_kind_set, atomic_kind_set, particle_set, &
    1016            2 :                       alpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore)
    1017            8 :       atewd(1:natom) = atewd(1:natom) + 0.5_dp*atcore(1:natom)
    1018            2 :       CALL dbcsr_add(vh_mat%matrix, e_mat%matrix, 1.0_dp, 0.5_dp)
    1019            6 :       DO ikind = 1, nkind
    1020              :          CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha(ikind), &
    1021            6 :                           ccore_charge=ccore(ikind))
    1022              :       END DO
    1023            2 :       CALL dbcsr_set(e_mat%matrix, 0.0_dp)
    1024            8 :       atcore = 0.0_dp
    1025              :       CALL build_erfc(e_mat, matrix_p, qs_kind_set, atomic_kind_set, particle_set, &
    1026            2 :                       alpha, ccore, eps_core_charge, sab_orb, sac_ae, atcore)
    1027            8 :       atewd(1:natom) = atewd(1:natom) - 0.5_dp*atcore(1:natom)
    1028            2 :       CALL dbcsr_add(vh_mat%matrix, e_mat%matrix, 1.0_dp, -0.5_dp)
    1029              :       !
    1030            2 :       CALL dbcsr_release(e_mat%matrix)
    1031            2 :       DEALLOCATE (e_mat%matrix)
    1032            2 :       CALL auxbas_pw_pool%give_back_pw(v_hewd_gspace)
    1033            2 :       CALL auxbas_pw_pool%give_back_pw(v_hewd_rspace)
    1034            2 :       CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_g)
    1035            2 :       CALL auxbas_pw_pool%give_back_pw(rho_tot_ewd_r)
    1036            2 :       DEALLOCATE (atecc, atcore)
    1037            2 :       DEALLOCATE (alpha, ccore)
    1038              : 
    1039            2 :       CALL timestop(handle)
    1040              : 
    1041            4 :    END SUBROUTINE vh_ewald_correction
    1042              : 
    1043              : ! **************************************************************************************************
    1044              : !> \brief ...
    1045              : !> \param qs_env ...
    1046              : !> \param matrix_hfx ...
    1047              : ! **************************************************************************************************
    1048            2 :    SUBROUTINE get_hfx_ks_matrix(qs_env, matrix_hfx)
    1049              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1050              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_hfx
    1051              : 
    1052              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'get_hfx_ks_matrix'
    1053              : 
    1054              :       INTEGER                                            :: handle
    1055            2 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p
    1056              :       TYPE(qs_rho_type), POINTER                         :: rho
    1057              : 
    1058            2 :       CALL timeset(routineN, handle)
    1059              : 
    1060            2 :       CALL get_qs_env(qs_env, rho=rho)
    1061            2 :       CALL qs_rho_get(rho, rho_ao=matrix_p)
    1062              :       CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, update_energy=.FALSE., &
    1063            2 :                             recalc_integrals=.FALSE.)
    1064              : 
    1065            2 :       CALL timestop(handle)
    1066              : 
    1067            2 :    END SUBROUTINE get_hfx_ks_matrix
    1068              : ! **************************************************************************************************
    1069              : !> \brief Write the atomic coordinates, types, and energies
    1070              : !> \param iounit ...
    1071              : !> \param particle_set ...
    1072              : !> \param atener ...
    1073              : !> \param label ...
    1074              : !> \date    05.06.2023
    1075              : !> \author  JGH
    1076              : !> \version 1.0
    1077              : ! **************************************************************************************************
    1078           16 :    SUBROUTINE write_atener(iounit, particle_set, atener, label)
    1079              : 
    1080              :       INTEGER, INTENT(IN)                                :: iounit
    1081              :       TYPE(particle_type), DIMENSION(:)                  :: particle_set
    1082              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: atener
    1083              :       CHARACTER(LEN=*), INTENT(IN)                       :: label
    1084              : 
    1085              :       INTEGER                                            :: i, natom
    1086              : 
    1087           16 :       IF (iounit > 0) THEN
    1088           16 :          WRITE (UNIT=iounit, FMT="(/,T2,A)") TRIM(label)
    1089              :          WRITE (UNIT=iounit, FMT="(T4,A,T30,A,T42,A,T54,A,T69,A)") &
    1090           16 :             "Atom  Element", "X", "Y", "Z", "Energy[a.u.]"
    1091           16 :          natom = SIZE(atener)
    1092           72 :          DO i = 1, natom
    1093           56 :             WRITE (UNIT=iounit, FMT="(I6,T12,A2,T24,3F12.6,F21.10)") i, &
    1094           56 :                TRIM(ADJUSTL(particle_set(i)%atomic_kind%element_symbol)), &
    1095          296 :                particle_set(i)%r(1:3)*angstrom, atener(i)
    1096              :          END DO
    1097           16 :          WRITE (UNIT=iounit, FMT="(A)") ""
    1098              :       END IF
    1099              : 
    1100           16 :    END SUBROUTINE write_atener
    1101              : 
    1102              : ! **************************************************************************************************
    1103              : !> \brief Write the one component of atomic energies
    1104              : !> \param iounit ...
    1105              : !> \param atener ...
    1106              : !> \param label ...
    1107              : !> \date    18.08.2023
    1108              : !> \author  JGH
    1109              : !> \version 1.0
    1110              : ! **************************************************************************************************
    1111           45 :    SUBROUTINE write_atdet(iounit, atener, label)
    1112              :       INTEGER, INTENT(IN)                                :: iounit
    1113              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: atener
    1114              :       CHARACTER(LEN=*), INTENT(IN)                       :: label
    1115              : 
    1116              :       INTEGER                                            :: natom
    1117              : 
    1118           45 :       IF (iounit > 0) THEN
    1119           45 :          natom = SIZE(atener)
    1120           45 :          WRITE (UNIT=iounit, FMT="(T2,A)") TRIM(label)
    1121           45 :          WRITE (UNIT=iounit, FMT="(5G16.8)") atener(1:natom)
    1122              :       END IF
    1123              : 
    1124           45 :    END SUBROUTINE write_atdet
    1125              : 
    1126              : END MODULE ed_analysis
        

Generated by: LCOV version 2.0-1