LCOV - code coverage report
Current view: top level - src - qs_energy_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:2fce0f8) Lines: 167 180 92.8 %
Date: 2024-12-21 06:28:57 Functions: 3 3 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Utility subroutine for qs energy calculation
      10             : !> \par History
      11             : !>      none
      12             : !> \author MK (29.10.2002)
      13             : ! **************************************************************************************************
      14             : MODULE qs_energy_utils
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      16             :    USE atprop_types,                    ONLY: atprop_array_add,&
      17             :                                               atprop_array_init,&
      18             :                                               atprop_type
      19             :    USE cp_control_types,                ONLY: dft_control_type
      20             :    USE cp_control_utils,                ONLY: read_ddapc_section
      21             :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      22             :                                               dbcsr_create,&
      23             :                                               dbcsr_p_type,&
      24             :                                               dbcsr_release,&
      25             :                                               dbcsr_set
      26             :    USE et_coupling,                     ONLY: calc_et_coupling
      27             :    USE et_coupling_proj,                ONLY: calc_et_coupling_proj
      28             :    USE hartree_local_methods,           ONLY: Vh_1c_gg_integrals
      29             :    USE hartree_local_types,             ONLY: ecoul_1center_type
      30             :    USE input_section_types,             ONLY: section_vals_get,&
      31             :                                               section_vals_get_subs_vals,&
      32             :                                               section_vals_type
      33             :    USE kinds,                           ONLY: dp
      34             :    USE message_passing,                 ONLY: mp_para_env_type
      35             :    USE mulliken,                        ONLY: atom_trace
      36             :    USE post_scf_bandstructure_methods,  ONLY: post_scf_bandstructure
      37             :    USE pw_env_types,                    ONLY: pw_env_get,&
      38             :                                               pw_env_type
      39             :    USE pw_methods,                      ONLY: pw_axpy,&
      40             :                                               pw_scale
      41             :    USE pw_pool_types,                   ONLY: pw_pool_type
      42             :    USE pw_types,                        ONLY: pw_r3d_rs_type
      43             :    USE qs_core_hamiltonian,             ONLY: core_matrices
      44             :    USE qs_dispersion_types,             ONLY: qs_dispersion_type
      45             :    USE qs_energy_types,                 ONLY: qs_energy_type
      46             :    USE qs_environment_types,            ONLY: get_qs_env,&
      47             :                                               qs_environment_type
      48             :    USE qs_integrate_potential,          ONLY: integrate_v_core_rspace,&
      49             :                                               integrate_v_rspace
      50             :    USE qs_kind_types,                   ONLY: qs_kind_type
      51             :    USE qs_ks_atom,                      ONLY: update_ks_atom
      52             :    USE qs_ks_methods,                   ONLY: qs_ks_update_qs_env
      53             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      54             :    USE qs_linres_module,                ONLY: linres_calculation_low
      55             :    USE qs_local_rho_types,              ONLY: local_rho_type
      56             :    USE qs_rho0_ggrid,                   ONLY: integrate_vhg0_rspace
      57             :    USE qs_rho_atom_types,               ONLY: rho_atom_type,&
      58             :                                               zero_rho_atom_integrals
      59             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      60             :                                               qs_rho_type
      61             :    USE qs_scf,                          ONLY: scf
      62             :    USE qs_tddfpt2_methods,              ONLY: tddfpt
      63             :    USE qs_vxc,                          ONLY: qs_xc_density
      64             :    USE qs_vxc_atom,                     ONLY: calculate_vxc_atom
      65             :    USE tip_scan_methods,                ONLY: tip_scanning
      66             :    USE xas_methods,                     ONLY: xas
      67             :    USE xas_tdp_methods,                 ONLY: xas_tdp
      68             :    USE xc_derivatives,                  ONLY: xc_functionals_get_needs
      69             :    USE xc_rho_cflags_types,             ONLY: xc_rho_cflags_type
      70             : #include "./base/base_uses.f90"
      71             : 
      72             :    IMPLICIT NONE
      73             : 
      74             :    PRIVATE
      75             : 
      76             : ! *** Global parameters ***
      77             : 
      78             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_utils'
      79             : 
      80             :    PUBLIC :: qs_energies_properties
      81             : 
      82             : CONTAINS
      83             : 
      84             : ! **************************************************************************************************
      85             : !> \brief Refactoring of qs_energies_scf. Moves computation of properties
      86             : !>        into separate subroutine
      87             : !> \param qs_env ...
      88             : !> \param calc_forces ...
      89             : !> \par History
      90             : !>      05.2013 created [Florian Schiffmann]
      91             : ! **************************************************************************************************
      92             : 
      93       81356 :    SUBROUTINE qs_energies_properties(qs_env, calc_forces)
      94             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      95             :       LOGICAL, INTENT(IN)                                :: calc_forces
      96             : 
      97             :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_properties'
      98             : 
      99             :       INTEGER                                            :: handle, natom
     100             :       LOGICAL                                            :: do_et, do_et_proj, &
     101             :                                                             do_post_scf_bandstructure, do_tip_scan
     102             :       REAL(KIND=dp)                                      :: ekts
     103             :       TYPE(atprop_type), POINTER                         :: atprop
     104             :       TYPE(dft_control_type), POINTER                    :: dft_control
     105             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     106             :       TYPE(pw_env_type), POINTER                         :: pw_env
     107             :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
     108             :       TYPE(qs_energy_type), POINTER                      :: energy
     109             :       TYPE(section_vals_type), POINTER                   :: input, post_scf_bands_section, &
     110             :                                                             proj_section, rest_b_section, &
     111             :                                                             tip_scan_section
     112             : 
     113       20339 :       NULLIFY (atprop, energy, pw_env)
     114       20339 :       CALL timeset(routineN, handle)
     115             : 
     116             :       ! atomic energies using Mulliken partition
     117             :       CALL get_qs_env(qs_env, &
     118             :                       dft_control=dft_control, &
     119             :                       input=input, &
     120             :                       atprop=atprop, &
     121             :                       energy=energy, &
     122             :                       v_hartree_rspace=v_hartree_rspace, &
     123             :                       para_env=para_env, &
     124       20339 :                       pw_env=pw_env)
     125       20339 :       IF (atprop%energy) THEN
     126         140 :          CALL qs_energies_mulliken(qs_env)
     127         140 :          CALL get_qs_env(qs_env, natom=natom)
     128             :          IF (.NOT. dft_control%qs_control%semi_empirical .AND. &
     129         140 :              .NOT. dft_control%qs_control%xtb .AND. &
     130             :              .NOT. dft_control%qs_control%dftb) THEN
     131             :             ! Nuclear charge correction
     132          36 :             CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
     133          36 :             IF (.NOT. ASSOCIATED(atprop%ateb)) THEN
     134           8 :                CALL atprop_array_init(atprop%ateb, natom)
     135             :             END IF
     136             :             ! Kohn-Sham Functional corrections
     137          36 :             CALL ks_xc_correction(qs_env)
     138             :          END IF
     139         140 :          CALL atprop_array_add(atprop%atener, atprop%ateb)
     140         140 :          CALL atprop_array_add(atprop%atener, atprop%ateself)
     141         140 :          CALL atprop_array_add(atprop%atener, atprop%atexc)
     142         140 :          CALL atprop_array_add(atprop%atener, atprop%atecoul)
     143         140 :          CALL atprop_array_add(atprop%atener, atprop%atevdw)
     144         140 :          CALL atprop_array_add(atprop%atener, atprop%ategcp)
     145         140 :          CALL atprop_array_add(atprop%atener, atprop%atecc)
     146         140 :          CALL atprop_array_add(atprop%atener, atprop%ate1c)
     147             :          ! entropic energy
     148         140 :          ekts = energy%kts/REAL(natom, KIND=dp)/REAL(para_env%num_pe, KIND=dp)
     149        7504 :          atprop%atener(:) = atprop%atener(:) + ekts
     150             :       END IF
     151             : 
     152             :       ! ET coupling - projection-operator approach
     153       20339 :       NULLIFY (proj_section)
     154             :       proj_section => &
     155       20339 :          section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING%PROJECTION")
     156       20339 :       CALL section_vals_get(proj_section, explicit=do_et_proj)
     157       20339 :       IF (do_et_proj) THEN
     158          10 :          CALL calc_et_coupling_proj(qs_env)
     159             :       END IF
     160             : 
     161             :       ! **********  Calculate the electron transfer coupling elements********
     162       20339 :       do_et = .FALSE.
     163       20339 :       do_et = dft_control%qs_control%et_coupling_calc
     164       20339 :       IF (do_et) THEN
     165           0 :          qs_env%et_coupling%energy = energy%total
     166           0 :          qs_env%et_coupling%keep_matrix = .TRUE.
     167           0 :          qs_env%et_coupling%first_run = .TRUE.
     168           0 :          CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.)
     169           0 :          qs_env%et_coupling%first_run = .FALSE.
     170           0 :          IF (dft_control%qs_control%ddapc_restraint) THEN
     171           0 :             rest_b_section => section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING%DDAPC_RESTRAINT_B")
     172             :             CALL read_ddapc_section(qs_control=dft_control%qs_control, &
     173           0 :                                     ddapc_restraint_section=rest_b_section)
     174             :          END IF
     175           0 :          CALL scf(qs_env=qs_env)
     176           0 :          qs_env%et_coupling%keep_matrix = .TRUE.
     177             : 
     178           0 :          CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.)
     179           0 :          CALL calc_et_coupling(qs_env)
     180             :       END IF
     181             : 
     182             :       !Properties
     183       20339 :       IF (dft_control%do_xas_calculation) THEN
     184          42 :          CALL xas(qs_env, dft_control)
     185             :       END IF
     186             : 
     187       20339 :       IF (dft_control%do_xas_tdp_calculation) THEN
     188          50 :          CALL xas_tdp(qs_env)
     189             :       END IF
     190             : 
     191             :       ! Compute Linear Response properties as post-scf
     192       20339 :       IF (.NOT. qs_env%linres_run) THEN
     193       19839 :          CALL linres_calculation_low(qs_env)
     194             :       END IF
     195             : 
     196       20339 :       IF (dft_control%tddfpt2_control%enabled) THEN
     197        1054 :          CALL tddfpt(qs_env, calc_forces)
     198             :       END IF
     199             : 
     200             :       ! post-SCF bandstructure calculation from higher level methods
     201       20339 :       NULLIFY (post_scf_bands_section)
     202       20339 :       post_scf_bands_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%BANDSTRUCTURE")
     203       20339 :       CALL section_vals_get(post_scf_bands_section, explicit=do_post_scf_bandstructure)
     204       20339 :       IF (do_post_scf_bandstructure) THEN
     205          34 :          CALL post_scf_bandstructure(qs_env, post_scf_bands_section)
     206             :       END IF
     207             : 
     208             :       ! tip scan
     209       20339 :       NULLIFY (tip_scan_section)
     210       20339 :       tip_scan_section => section_vals_get_subs_vals(input, "PROPERTIES%TIP_SCAN")
     211       20339 :       CALL section_vals_get(tip_scan_section, explicit=do_tip_scan)
     212       20339 :       IF (do_tip_scan) THEN
     213           0 :          CALL tip_scanning(qs_env, tip_scan_section)
     214             :       END IF
     215             : 
     216       20339 :       CALL timestop(handle)
     217             : 
     218       20339 :    END SUBROUTINE qs_energies_properties
     219             : 
     220             : ! **************************************************************************************************
     221             : !> \brief   Use a simple Mulliken-like energy decomposition
     222             : !> \param qs_env ...
     223             : !> \date    07.2011
     224             : !> \author  JHU
     225             : !> \version 1.0
     226             : ! **************************************************************************************************
     227         140 :    SUBROUTINE qs_energies_mulliken(qs_env)
     228             : 
     229             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     230             : 
     231             :       INTEGER                                            :: ispin, natom, nspin
     232         140 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: atcore
     233             :       TYPE(atprop_type), POINTER                         :: atprop
     234             :       TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), &
     235         140 :          TARGET                                          :: core_mat
     236         140 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_h, matrix_ks, rho_ao
     237         140 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: math, matp
     238             :       TYPE(dft_control_type), POINTER                    :: dft_control
     239             :       TYPE(qs_rho_type), POINTER                         :: rho
     240             : 
     241         140 :       CALL get_qs_env(qs_env=qs_env, atprop=atprop)
     242         140 :       IF (atprop%energy) THEN
     243         140 :          CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_h=matrix_h, rho=rho)
     244         140 :          CALL qs_rho_get(rho, rho_ao=rho_ao)
     245             :          ! E = 0.5*Tr(H*P+F*P)
     246        7504 :          atprop%atener = 0._dp
     247         140 :          nspin = SIZE(rho_ao)
     248         284 :          DO ispin = 1, nspin
     249             :             CALL atom_trace(matrix_h(1)%matrix, rho_ao(ispin)%matrix, &
     250         144 :                             0.5_dp, atprop%atener)
     251             :             CALL atom_trace(matrix_ks(ispin)%matrix, rho_ao(ispin)%matrix, &
     252         284 :                             0.5_dp, atprop%atener)
     253             :          END DO
     254             :          !
     255         140 :          CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     256             :          IF (.NOT. dft_control%qs_control%semi_empirical .AND. &
     257         140 :              .NOT. dft_control%qs_control%xtb .AND. &
     258             :              .NOT. dft_control%qs_control%dftb) THEN
     259          36 :             CALL get_qs_env(qs_env=qs_env, natom=natom)
     260         108 :             ALLOCATE (atcore(natom))
     261         160 :             atcore = 0.0_dp
     262         108 :             ALLOCATE (core_mat(1))
     263          36 :             ALLOCATE (core_mat(1)%matrix)
     264          36 :             CALL dbcsr_create(core_mat(1)%matrix, template=matrix_h(1)%matrix)
     265          36 :             CALL dbcsr_copy(core_mat(1)%matrix, matrix_h(1)%matrix)
     266          36 :             CALL dbcsr_set(core_mat(1)%matrix, 0.0_dp)
     267          36 :             math(1:1, 1:1) => core_mat(1:1)
     268          36 :             matp(1:nspin, 1:1) => rho_ao(1:nspin)
     269          36 :             CALL core_matrices(qs_env, math, matp, .FALSE., 0, atcore=atcore)
     270         160 :             atprop%atener = atprop%atener + 0.5_dp*atcore
     271          76 :             DO ispin = 1, nspin
     272             :                CALL atom_trace(core_mat(1)%matrix, rho_ao(ispin)%matrix, &
     273          76 :                                -0.5_dp, atprop%atener)
     274             :             END DO
     275          36 :             DEALLOCATE (atcore)
     276          36 :             CALL dbcsr_release(core_mat(1)%matrix)
     277          36 :             DEALLOCATE (core_mat(1)%matrix)
     278          36 :             DEALLOCATE (core_mat)
     279             :          END IF
     280             :       END IF
     281             : 
     282         280 :    END SUBROUTINE qs_energies_mulliken
     283             : 
     284             : ! **************************************************************************************************
     285             : !> \brief ...
     286             : !> \param qs_env ...
     287             : ! **************************************************************************************************
     288          36 :    SUBROUTINE ks_xc_correction(qs_env)
     289             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     290             : 
     291             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ks_xc_correction'
     292             : 
     293             :       INTEGER                                            :: handle, iatom, ispin, natom, nspins
     294             :       LOGICAL                                            :: gapw, gapw_xc
     295             :       REAL(KIND=dp)                                      :: eh1, exc1
     296          36 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     297             :       TYPE(atprop_type), POINTER                         :: atprop
     298          36 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, rho_ao, xcmat
     299          36 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p
     300             :       TYPE(dft_control_type), POINTER                    :: dft_control
     301          36 :       TYPE(ecoul_1center_type), DIMENSION(:), POINTER    :: ecoul_1c
     302             :       TYPE(local_rho_type), POINTER                      :: local_rho_set
     303             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     304             :       TYPE(pw_env_type), POINTER                         :: pw_env
     305             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     306             :       TYPE(pw_r3d_rs_type)                               :: xc_den
     307          36 :       TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:)    :: vtau, vxc
     308             :       TYPE(pw_r3d_rs_type), POINTER                      :: v_hartree_rspace
     309             :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     310          36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     311             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     312             :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     313          36 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
     314             :       TYPE(section_vals_type), POINTER                   :: xc_fun_section, xc_section
     315             :       TYPE(xc_rho_cflags_type)                           :: needs
     316             : 
     317          36 :       CALL timeset(routineN, handle)
     318             : 
     319          36 :       CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env, atprop=atprop)
     320             : 
     321          36 :       IF (atprop%energy) THEN
     322             : 
     323          36 :          nspins = dft_control%nspins
     324          36 :          xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
     325          36 :          xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     326          36 :          needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
     327          36 :          gapw = dft_control%qs_control%gapw
     328          36 :          gapw_xc = dft_control%qs_control%gapw_xc
     329             : 
     330             :          ! Nuclear charge correction
     331          36 :          CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
     332          36 :          IF (gapw .OR. gapw_xc) THEN
     333             :             CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, &
     334             :                             rho_atom_set=rho_atom_set, ecoul_1c=ecoul_1c, &
     335          12 :                             natom=natom, para_env=para_env)
     336          12 :             CALL zero_rho_atom_integrals(rho_atom_set)
     337          12 :             CALL calculate_vxc_atom(qs_env, .FALSE., exc1)
     338          12 :             IF (gapw) THEN
     339           8 :                CALL Vh_1c_gg_integrals(qs_env, eh1, ecoul_1c, local_rho_set, para_env, tddft=.FALSE.)
     340           8 :                CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
     341             :                CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., &
     342           8 :                                           local_rho_set=local_rho_set, atener=atprop%ateb)
     343             :             END IF
     344             :          END IF
     345             : 
     346          36 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     347          36 :          CALL auxbas_pw_pool%create_pw(xc_den)
     348         148 :          ALLOCATE (vxc(nspins))
     349          76 :          DO ispin = 1, nspins
     350          76 :             CALL auxbas_pw_pool%create_pw(vxc(ispin))
     351             :          END DO
     352          36 :          IF (needs%tau .OR. needs%tau_spin) THEN
     353          36 :             ALLOCATE (vtau(nspins))
     354          24 :             DO ispin = 1, nspins
     355          48 :                CALL auxbas_pw_pool%create_pw(vtau(ispin))
     356             :             END DO
     357             :          END IF
     358             : 
     359          36 :          IF (gapw_xc) THEN
     360           4 :             CALL get_qs_env(qs_env, rho_xc=rho_struct, dispersion_env=dispersion_env)
     361             :          ELSE
     362          32 :             CALL get_qs_env(qs_env, rho=rho_struct, dispersion_env=dispersion_env)
     363             :          END IF
     364          36 :          IF (needs%tau .OR. needs%tau_spin) THEN
     365             :             CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
     366          12 :                                xc_den=xc_den, vxc=vxc, vtau=vtau)
     367             :          ELSE
     368             :             CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
     369          24 :                                xc_den=xc_den, vxc=vxc)
     370             :          END IF
     371          36 :          CALL get_qs_env(qs_env, rho=rho_struct)
     372          36 :          CALL qs_rho_get(rho_struct, rho_ao=rho_ao)
     373          36 :          CALL get_qs_env(qs_env, natom=natom, matrix_s=matrix_s)
     374          36 :          CALL atprop_array_init(atprop%atexc, natom)
     375         148 :          ALLOCATE (xcmat(nspins))
     376          76 :          DO ispin = 1, nspins
     377          40 :             ALLOCATE (xcmat(ispin)%matrix)
     378          40 :             CALL dbcsr_create(xcmat(ispin)%matrix, template=matrix_s(1)%matrix)
     379          40 :             CALL dbcsr_copy(xcmat(ispin)%matrix, matrix_s(1)%matrix)
     380          40 :             CALL dbcsr_set(xcmat(ispin)%matrix, 0.0_dp)
     381          40 :             CALL pw_scale(vxc(ispin), -0.5_dp)
     382          40 :             CALL pw_axpy(xc_den, vxc(ispin))
     383          40 :             CALL pw_scale(vxc(ispin), vxc(ispin)%pw_grid%dvol)
     384             :             CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vxc(ispin), hmat=xcmat(ispin), &
     385          40 :                                     calculate_forces=.FALSE., gapw=(gapw .OR. gapw_xc))
     386          76 :             IF (needs%tau .OR. needs%tau_spin) THEN
     387          12 :                CALL pw_scale(vtau(ispin), -0.5_dp*vtau(ispin)%pw_grid%dvol)
     388             :                CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vtau(ispin), &
     389             :                                        hmat=xcmat(ispin), calculate_forces=.FALSE., &
     390          12 :                                        gapw=(gapw .OR. gapw_xc), compute_tau=.TRUE.)
     391             :             END IF
     392             :          END DO
     393          36 :          IF (gapw .OR. gapw_xc) THEN
     394             :             ! remove one-center potential matrix part
     395          12 :             CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
     396          12 :             CALL update_ks_atom(qs_env, xcmat, matrix_p, forces=.FALSE., kscale=-0.5_dp)
     397          12 :             CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
     398          12 :             CALL atprop_array_init(atprop%ate1c, natom)
     399          48 :             atprop%ate1c = 0.0_dp
     400          48 :             DO iatom = 1, natom
     401             :                atprop%ate1c(iatom) = atprop%ate1c(iatom) + &
     402          48 :                                      rho_atom_set(iatom)%exc_h - rho_atom_set(iatom)%exc_s
     403             :             END DO
     404          12 :             IF (gapw) THEN
     405           8 :                CALL get_qs_env(qs_env=qs_env, ecoul_1c=ecoul_1c)
     406          32 :                DO iatom = 1, natom
     407             :                   atprop%ate1c(iatom) = atprop%ate1c(iatom) + &
     408             :                                         ecoul_1c(iatom)%ecoul_1_h - ecoul_1c(iatom)%ecoul_1_s + &
     409          32 :                                         ecoul_1c(iatom)%ecoul_1_z - ecoul_1c(iatom)%ecoul_1_0
     410             :                END DO
     411             :             END IF
     412             :          END IF
     413          76 :          DO ispin = 1, nspins
     414          40 :             CALL atom_trace(xcmat(ispin)%matrix, rho_ao(ispin)%matrix, 1.0_dp, atprop%atexc)
     415          40 :             CALL dbcsr_release(xcmat(ispin)%matrix)
     416          76 :             DEALLOCATE (xcmat(ispin)%matrix)
     417             :          END DO
     418          36 :          DEALLOCATE (xcmat)
     419             : 
     420          36 :          CALL auxbas_pw_pool%give_back_pw(xc_den)
     421          76 :          DO ispin = 1, nspins
     422          76 :             CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
     423             :          END DO
     424          36 :          IF (needs%tau .OR. needs%tau_spin) THEN
     425          24 :             DO ispin = 1, nspins
     426          48 :                CALL auxbas_pw_pool%give_back_pw(vtau(ispin))
     427             :             END DO
     428             :          END IF
     429             : 
     430             :       END IF
     431             : 
     432          36 :       CALL timestop(handle)
     433             : 
     434          72 :    END SUBROUTINE ks_xc_correction
     435             : 
     436             : END MODULE qs_energy_utils

Generated by: LCOV version 1.15