LCOV - code coverage report
Current view: top level - src - kg_correction.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 87.0 % 338 294
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 5 5

            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 Routines for a Kim-Gordon-like partitioning into molecular subunits
      10              : !> \par History
      11              : !>       2012.06 created [Martin Haeufel]
      12              : !> \author Martin Haeufel and Florian Schiffmann
      13              : ! **************************************************************************************************
      14              : MODULE kg_correction
      15              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      16              :    USE cp_control_types,                ONLY: dft_control_type
      17              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      18              :                                               dbcsr_p_type
      19              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot
      20              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      21              :                                               cp_logger_get_default_unit_nr,&
      22              :                                               cp_logger_type
      23              :    USE ec_methods,                      ONLY: create_kernel
      24              :    USE input_constants,                 ONLY: kg_tnadd_atomic,&
      25              :                                               kg_tnadd_embed,&
      26              :                                               kg_tnadd_embed_ri,&
      27              :                                               kg_tnadd_none
      28              :    USE input_section_types,             ONLY: section_vals_type
      29              :    USE kg_environment_types,            ONLY: kg_environment_type
      30              :    USE kinds,                           ONLY: dp
      31              :    USE lri_environment_methods,         ONLY: calculate_lri_densities,&
      32              :                                               lri_kg_rho_update
      33              :    USE lri_environment_types,           ONLY: lri_density_type,&
      34              :                                               lri_environment_type,&
      35              :                                               lri_kind_type
      36              :    USE lri_forces,                      ONLY: calculate_lri_forces
      37              :    USE lri_ks_methods,                  ONLY: calculate_lri_ks_matrix
      38              :    USE message_passing,                 ONLY: mp_para_env_type
      39              :    USE pw_env_types,                    ONLY: pw_env_get,&
      40              :                                               pw_env_type
      41              :    USE pw_methods,                      ONLY: pw_integral_ab,&
      42              :                                               pw_scale
      43              :    USE pw_pool_types,                   ONLY: pw_pool_type
      44              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      45              :                                               pw_r3d_rs_type
      46              :    USE qs_environment_types,            ONLY: get_qs_env,&
      47              :                                               qs_environment_type
      48              :    USE qs_integrate_potential,          ONLY: integrate_v_rspace,&
      49              :                                               integrate_v_rspace_one_center
      50              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      51              :    USE qs_rho_methods,                  ONLY: qs_rho_rebuild,&
      52              :                                               qs_rho_update_rho
      53              :    USE qs_rho_types,                    ONLY: qs_rho_create,&
      54              :                                               qs_rho_get,&
      55              :                                               qs_rho_release,&
      56              :                                               qs_rho_set,&
      57              :                                               qs_rho_type,&
      58              :                                               qs_rho_unset_rho_ao
      59              :    USE qs_vxc,                          ONLY: qs_vxc_create
      60              :    USE virial_types,                    ONLY: virial_type
      61              : #include "./base/base_uses.f90"
      62              : 
      63              :    IMPLICIT NONE
      64              : 
      65              :    PRIVATE
      66              : 
      67              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kg_correction'
      68              : 
      69              :    PUBLIC :: kg_ekin_subset
      70              : 
      71              : CONTAINS
      72              : 
      73              : ! **************************************************************************************************
      74              : !> \brief Calculates the subsystem Hohenberg-Kohn kinetic energy and the forces
      75              : !> \param qs_env ...
      76              : !> \param ks_matrix ...
      77              : !> \param ekin_mol ...
      78              : !> \param calc_force ...
      79              : !> \param do_kernel Contribution of kinetic energy functional to kernel in response calculation
      80              : !> \param pmat_ext Response density used to fold 2nd deriv or to integrate kinetic energy functional
      81              : !> \par History
      82              : !>       2012.06 created [Martin Haeufel]
      83              : !>       2014.01 added atomic potential option [JGH]
      84              : !>       2020.01 Added KG contribution to linear response [fbelle]
      85              : !> \author Martin Haeufel and Florian Schiffmann
      86              : ! **************************************************************************************************
      87          846 :    SUBROUTINE kg_ekin_subset(qs_env, ks_matrix, ekin_mol, calc_force, do_kernel, pmat_ext)
      88              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      89              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix
      90              :       REAL(KIND=dp), INTENT(out)                         :: ekin_mol
      91              :       LOGICAL, INTENT(IN)                                :: calc_force, do_kernel
      92              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
      93              :          POINTER                                         :: pmat_ext
      94              : 
      95              :       LOGICAL                                            :: lrigpw
      96              :       TYPE(dft_control_type), POINTER                    :: dft_control
      97              :       TYPE(kg_environment_type), POINTER                 :: kg_env
      98              : 
      99          846 :       CALL get_qs_env(qs_env, kg_env=kg_env, dft_control=dft_control)
     100          846 :       lrigpw = dft_control%qs_control%lrigpw
     101              : 
     102          846 :       IF (kg_env%tnadd_method == kg_tnadd_embed) THEN
     103          644 :          IF (lrigpw) THEN
     104           20 :             CALL kg_ekin_embed_lri(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
     105              :          ELSE
     106              :             CALL kg_ekin_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force, &
     107          624 :                                do_kernel, pmat_ext)
     108              :          END IF
     109          202 :       ELSE IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
     110              :          CALL kg_ekin_ri_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force, &
     111           20 :                                do_kernel, pmat_ext)
     112          182 :       ELSE IF (kg_env%tnadd_method == kg_tnadd_atomic) THEN
     113          144 :          CALL kg_ekin_atomic(qs_env, ks_matrix, ekin_mol)
     114           38 :       ELSE IF (kg_env%tnadd_method == kg_tnadd_none) THEN
     115           38 :          ekin_mol = 0.0_dp
     116              :       ELSE
     117            0 :          CPABORT("Unknown KG embedding method")
     118              :       END IF
     119              : 
     120          846 :    END SUBROUTINE kg_ekin_subset
     121              : 
     122              : ! **************************************************************************************************
     123              : !> \brief ...
     124              : !> \param qs_env ...
     125              : !> \param kg_env ...
     126              : !> \param ks_matrix ...
     127              : !> \param ekin_mol ...
     128              : !> \param calc_force ...
     129              : !> \param do_kernel Contribution of kinetic energy functional to kernel in response calculation
     130              : !> \param pmat_ext Response density used to fold 2nd deriv or to integrate kinetic energy functional
     131              : ! **************************************************************************************************
     132         1248 :    SUBROUTINE kg_ekin_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force, do_kernel, pmat_ext)
     133              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     134              :       TYPE(kg_environment_type), POINTER                 :: kg_env
     135              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix
     136              :       REAL(KIND=dp), INTENT(out)                         :: ekin_mol
     137              :       LOGICAL, INTENT(IN)                                :: calc_force, do_kernel
     138              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
     139              :          POINTER                                         :: pmat_ext
     140              : 
     141              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'kg_ekin_embed'
     142              : 
     143              :       INTEGER                                            :: handle, iounit, ispin, isub, nspins
     144              :       LOGICAL                                            :: use_virial
     145              :       REAL(KIND=dp)                                      :: alpha, ekin_imol
     146              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: xcvirial
     147              :       TYPE(cp_logger_type), POINTER                      :: logger
     148          624 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: density_matrix
     149              :       TYPE(dft_control_type), POINTER                    :: dft_control
     150          624 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho1_g
     151              :       TYPE(pw_env_type), POINTER                         :: pw_env
     152              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     153          624 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r, rho_r, tau1_r, vxc_rho, vxc_tau
     154              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     155              :       TYPE(qs_rho_type), POINTER                         :: old_rho, rho1, rho_struct
     156              :       TYPE(section_vals_type), POINTER                   :: xc_section
     157              :       TYPE(virial_type), POINTER                         :: virial
     158              : 
     159          624 :       CALL timeset(routineN, handle)
     160              : 
     161          624 :       logger => cp_get_default_logger()
     162          624 :       iounit = cp_logger_get_default_unit_nr(logger)
     163              : 
     164          624 :       NULLIFY (ks_env, dft_control, old_rho, pw_env, rho_struct, virial, vxc_rho, vxc_tau)
     165              : 
     166              :       CALL get_qs_env(qs_env, &
     167              :                       ks_env=ks_env, &
     168              :                       rho=old_rho, &
     169              :                       dft_control=dft_control, &
     170              :                       virial=virial, &
     171          624 :                       pw_env=pw_env)
     172          624 :       nspins = dft_control%nspins
     173          624 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     174          624 :       use_virial = use_virial .AND. calc_force
     175              : 
     176              :       ! Kernel potential in response calculation (no forces calculated at this point)
     177              :       ! requires spin-factor
     178              :       ! alpha = 2 closed-shell
     179              :       ! alpha = 1 open-shell
     180          624 :       alpha = 1.0_dp
     181          624 :       IF (do_kernel .AND. .NOT. calc_force .AND. nspins == 1) alpha = 2.0_dp
     182              : 
     183          624 :       NULLIFY (auxbas_pw_pool)
     184          624 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     185              : 
     186              :       ! get the density matrix
     187          624 :       CALL qs_rho_get(old_rho, rho_ao=density_matrix)
     188              :       ! allocate and initialize the density
     189          624 :       ALLOCATE (rho_struct)
     190          624 :       CALL qs_rho_create(rho_struct)
     191              :       ! set the density matrix to the blocked matrix
     192          624 :       CALL qs_rho_set(rho_struct, rho_ao=density_matrix) ! blocked_matrix
     193          624 :       CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
     194              :       ! full density kinetic energy term
     195          624 :       CALL qs_rho_update_rho(rho_struct, qs_env)
     196              :       ! get blocked density that has been put on grid
     197          624 :       CALL qs_rho_get(rho_struct, rho_r=rho_r)
     198              : 
     199              :       ! If external density associated then it is needed either for
     200              :       ! 1) folding of second derivative while partially integrating, or
     201              :       ! 2) integration of response forces
     202          624 :       NULLIFY (rho1)
     203          624 :       IF (PRESENT(pmat_ext)) THEN
     204           58 :          ALLOCATE (rho1)
     205           58 :          CALL qs_rho_create(rho1)
     206           58 :          CALL qs_rho_set(rho1, rho_ao=pmat_ext)
     207           58 :          CALL qs_rho_rebuild(rho1, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
     208           58 :          CALL qs_rho_update_rho(rho1, qs_env)
     209              :       END IF
     210              : 
     211              :       ! XC-section pointing to kinetic energy functional in KG environment
     212              :       NULLIFY (xc_section)
     213          624 :       xc_section => kg_env%xc_section_kg
     214              : 
     215          624 :       ekin_imol = 0.0_dp
     216              : 
     217              :       ! calculate xc potential or kernel
     218          624 :       IF (do_kernel) THEN
     219              :          ! derivation wrt to rho_struct and evaluation at rho_struct
     220          142 :          IF (use_virial) virial%pv_xc = 0.0_dp
     221           46 :          CALL qs_rho_get(rho1, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
     222              :          CALL create_kernel(qs_env, &
     223              :                             vxc=vxc_rho, &
     224              :                             vxc_tau=vxc_tau, &
     225              :                             rho=rho_struct, &
     226              :                             rho1_r=rho1_r, &
     227              :                             rho1_g=rho1_g, &
     228              :                             tau1_r=tau1_r, &
     229              :                             xc_section=xc_section, &
     230              :                             compute_virial=use_virial, &
     231           46 :                             virial_xc=virial%pv_xc)
     232              :       ELSE
     233              :          CALL qs_vxc_create(ks_env=ks_env, &
     234              :                             rho_struct=rho_struct, &
     235              :                             xc_section=xc_section, &
     236              :                             vxc_rho=vxc_rho, &
     237              :                             vxc_tau=vxc_tau, &
     238          578 :                             exc=ekin_imol)
     239              :       END IF
     240              : 
     241          624 :       IF (ASSOCIATED(vxc_tau)) THEN
     242            0 :          CPABORT(" KG with meta-kinetic energy functionals not implemented")
     243              :       END IF
     244              : 
     245              :       ! Integrate xc-potential with external density for outer response forces
     246          624 :       IF (PRESENT(pmat_ext) .AND. .NOT. do_kernel) THEN
     247           12 :          CALL qs_rho_get(rho1, rho_ao=density_matrix, rho_r=rho1_r)
     248              :          ! Direct volume term of virial
     249              :          ! xc-potential is unscaled
     250           12 :          IF (use_virial) THEN
     251            8 :             ekin_imol = 0.0_dp
     252           16 :             DO ispin = 1, nspins
     253           16 :                ekin_imol = ekin_imol + pw_integral_ab(rho1_r(ispin), vxc_rho(ispin))
     254              :             END DO
     255              :          END IF
     256              :       END IF
     257              : 
     258         1248 :       DO ispin = 1, nspins
     259         1248 :          CALL pw_scale(vxc_rho(ispin), alpha*vxc_rho(ispin)%pw_grid%dvol)
     260              :       END DO
     261              : 
     262         1248 :       DO ispin = 1, nspins
     263              :          CALL integrate_v_rspace(v_rspace=vxc_rho(ispin), &
     264              :                                  pmat=density_matrix(ispin), hmat=ks_matrix(ispin), &
     265          624 :                                  qs_env=qs_env, calculate_forces=calc_force)
     266         1248 :          CALL auxbas_pw_pool%give_back_pw(vxc_rho(ispin))
     267              :       END DO
     268          624 :       DEALLOCATE (vxc_rho)
     269          624 :       ekin_mol = -ekin_imol
     270          624 :       xcvirial(1:3, 1:3) = 0.0_dp
     271          624 :       IF (use_virial) THEN
     272          208 :          xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) - virial%pv_xc(1:3, 1:3)
     273              :       END IF
     274              : 
     275              :       ! loop over all subsets
     276         1928 :       DO isub = 1, kg_env%nsubsets
     277              :          ! calculate the densities for the given blocked density matrix
     278              :          ! pass the subset task_list
     279              :          CALL qs_rho_update_rho(rho_struct, qs_env, &
     280         1304 :                                 task_list_external=kg_env%subset(isub)%task_list)
     281              :          ! Same for external (response) density if present
     282         1304 :          IF (PRESENT(pmat_ext)) THEN
     283              :             CALL qs_rho_update_rho(rho1, qs_env, &
     284          116 :                                    task_list_external=kg_env%subset(isub)%task_list)
     285              :          END IF
     286              : 
     287         1304 :          ekin_imol = 0.0_dp
     288         1304 :          NULLIFY (vxc_rho)
     289              : 
     290              :          ! calculate  Hohenberg-Kohn kinetic energy of the density
     291              :          ! corresponding to the remaining molecular block(s)
     292              :          ! info per block in rho_struct now
     293              : 
     294              :          ! calculate xc-potential or kernel
     295         1304 :          IF (do_kernel) THEN
     296          284 :             IF (use_virial) virial%pv_xc = 0.0_dp
     297           92 :             CALL qs_rho_get(rho1, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
     298              :             CALL create_kernel(qs_env, &
     299              :                                vxc=vxc_rho, &
     300              :                                vxc_tau=vxc_tau, &
     301              :                                rho=rho_struct, &
     302              :                                rho1_r=rho1_r, &
     303              :                                rho1_g=rho1_g, &
     304              :                                tau1_r=tau1_r, &
     305              :                                xc_section=xc_section, &
     306              :                                compute_virial=use_virial, &
     307           92 :                                virial_xc=virial%pv_xc)
     308              :          ELSE
     309              :             CALL qs_vxc_create(ks_env=ks_env, &
     310              :                                rho_struct=rho_struct, &
     311              :                                xc_section=xc_section, &
     312              :                                vxc_rho=vxc_rho, &
     313              :                                vxc_tau=vxc_tau, &
     314         1212 :                                exc=ekin_imol)
     315              :          END IF
     316              : 
     317              :          ! Integrate with response density for outer response forces
     318         1304 :          IF (PRESENT(pmat_ext) .AND. .NOT. do_kernel) THEN
     319           24 :             CALL qs_rho_get(rho1, rho_ao=density_matrix)
     320              :             ! Direct volume term of virial
     321              :             ! xc-potential is unscaled
     322           24 :             IF (use_virial) THEN
     323           16 :                ekin_imol = 0.0_dp
     324           32 :                DO ispin = 1, nspins
     325           32 :                   ekin_imol = ekin_imol + pw_integral_ab(rho1_r(ispin), vxc_rho(ispin))
     326              :                END DO
     327              :             END IF
     328              :          END IF
     329              : 
     330         2608 :          DO ispin = 1, nspins
     331         1304 :             CALL pw_scale(vxc_rho(ispin), -alpha*vxc_rho(ispin)%pw_grid%dvol)
     332              : 
     333              :             CALL integrate_v_rspace(v_rspace=vxc_rho(ispin), &
     334              :                                     pmat=density_matrix(ispin), &
     335              :                                     hmat=ks_matrix(ispin), &
     336              :                                     qs_env=qs_env, &
     337              :                                     calculate_forces=calc_force, &
     338         1304 :                                     task_list_external=kg_env%subset(isub)%task_list)
     339              :             ! clean up vxc_rho
     340         1304 :             CALL auxbas_pw_pool%give_back_pw(vxc_rho(ispin))
     341         2608 :             IF (ASSOCIATED(vxc_tau)) THEN
     342            0 :                CALL pw_scale(vxc_tau(ispin), -alpha*vxc_tau(ispin)%pw_grid%dvol)
     343              :                CALL integrate_v_rspace(v_rspace=vxc_tau(ispin), &
     344              :                                        pmat=density_matrix(ispin), &
     345              :                                        hmat=ks_matrix(ispin), &
     346              :                                        qs_env=qs_env, &
     347              :                                        compute_tau=.TRUE., &
     348              :                                        calculate_forces=calc_force, &
     349            0 :                                        task_list_external=kg_env%subset(isub)%task_list)
     350              :                ! clean up vxc_rho
     351            0 :                CALL auxbas_pw_pool%give_back_pw(vxc_tau(ispin))
     352              :             END IF
     353              :          END DO
     354         1304 :          DEALLOCATE (vxc_rho)
     355              : 
     356         1304 :          ekin_mol = ekin_mol + ekin_imol
     357              : 
     358         1928 :          IF (use_virial) THEN
     359          416 :             xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) + virial%pv_xc(1:3, 1:3)
     360              :          END IF
     361              : 
     362              :       END DO
     363              : 
     364          624 :       IF (use_virial) THEN
     365          208 :          virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3)
     366              :       END IF
     367              : 
     368              :       ! clean up rho_struct
     369          624 :       CALL qs_rho_unset_rho_ao(rho_struct)
     370          624 :       CALL qs_rho_release(rho_struct)
     371          624 :       DEALLOCATE (rho_struct)
     372          624 :       IF (PRESENT(pmat_ext)) THEN
     373           58 :          CALL qs_rho_unset_rho_ao(rho1)
     374           58 :          CALL qs_rho_release(rho1)
     375           58 :          DEALLOCATE (rho1)
     376              :       END IF
     377              : 
     378          624 :       CALL timestop(handle)
     379              : 
     380          624 :    END SUBROUTINE kg_ekin_embed
     381              : 
     382              : ! **************************************************************************************************
     383              : !> \brief ...
     384              : !> \param qs_env ...
     385              : !> \param kg_env ...
     386              : !> \param ks_matrix ...
     387              : !> \param ekin_mol ...
     388              : !> \param calc_force ...
     389              : ! **************************************************************************************************
     390           20 :    SUBROUTINE kg_ekin_embed_lri(qs_env, kg_env, ks_matrix, ekin_mol, calc_force)
     391              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     392              :       TYPE(kg_environment_type), POINTER                 :: kg_env
     393              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix
     394              :       REAL(KIND=dp), INTENT(out)                         :: ekin_mol
     395              :       LOGICAL                                            :: calc_force
     396              : 
     397              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'kg_ekin_embed_lri'
     398              : 
     399              :       INTEGER                                            :: color, handle, iatom, ikind, imol, &
     400              :                                                             ispin, isub, natom, nkind, nspins
     401              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atomlist
     402              :       LOGICAL                                            :: use_virial
     403              :       REAL(KIND=dp)                                      :: ekin_imol
     404              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: xcvirial
     405           20 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     406           20 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: density_matrix, ksmat
     407           20 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: pmat
     408              :       TYPE(dft_control_type), POINTER                    :: dft_control
     409              :       TYPE(lri_density_type), POINTER                    :: lri_density
     410              :       TYPE(lri_environment_type), POINTER                :: lri_env
     411           20 :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: lri_v_int
     412              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     413              :       TYPE(pw_env_type), POINTER                         :: pw_env
     414              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     415           20 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: vxc_rho, vxc_tau
     416              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     417              :       TYPE(qs_rho_type), POINTER                         :: old_rho, rho_struct
     418              :       TYPE(virial_type), POINTER                         :: virial
     419              : 
     420           20 :       CALL timeset(routineN, handle)
     421              : 
     422           20 :       NULLIFY (vxc_rho, vxc_tau, old_rho, rho_struct, ks_env)
     423              : 
     424           20 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     425              : 
     426              :       ! get set of molecules, natom, dft_control, pw_env
     427              :       CALL get_qs_env(qs_env, &
     428              :                       ks_env=ks_env, &
     429              :                       rho=old_rho, &
     430              :                       natom=natom, &
     431              :                       dft_control=dft_control, &
     432              :                       virial=virial, &
     433              :                       para_env=para_env, &
     434           20 :                       pw_env=pw_env)
     435              : 
     436           20 :       nspins = dft_control%nspins
     437           20 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     438            0 :       use_virial = use_virial .AND. calc_force
     439              : 
     440           20 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     441              : 
     442              :       ! get the density matrix
     443           20 :       CALL qs_rho_get(old_rho, rho_ao=density_matrix)
     444              :       ! allocate and initialize the density
     445           20 :       ALLOCATE (rho_struct)
     446           20 :       CALL qs_rho_create(rho_struct)
     447              :       ! set the density matrix to the blocked matrix
     448           20 :       CALL qs_rho_set(rho_struct, rho_ao=density_matrix) ! blocked_matrix
     449           20 :       CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
     450              : 
     451           20 :       CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density, nkind=nkind)
     452           20 :       IF (lri_env%exact_1c_terms) THEN
     453            0 :          CPABORT(" KG with LRI and exact one-center terms not implemented")
     454              :       END IF
     455           60 :       ALLOCATE (atomlist(natom))
     456           40 :       DO ispin = 1, nspins
     457           20 :          lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     458           80 :          DO ikind = 1, nkind
     459        14660 :             lri_v_int(ikind)%v_int = 0.0_dp
     460           60 :             IF (calc_force) THEN
     461           46 :                lri_v_int(ikind)%v_dadr = 0.0_dp
     462           46 :                lri_v_int(ikind)%v_dfdr = 0.0_dp
     463              :             END IF
     464              :          END DO
     465              :       END DO
     466              : 
     467              :       ! full density kinetic energy term
     468          120 :       atomlist = 1
     469           20 :       CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
     470              :       ekin_imol = 0.0_dp
     471              :       CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
     472           20 :                          vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
     473           20 :       IF (ASSOCIATED(vxc_tau)) THEN
     474            0 :          CPABORT(" KG with meta-kinetic energy functionals not implemented")
     475              :       END IF
     476           40 :       DO ispin = 1, nspins
     477           20 :          CALL pw_scale(vxc_rho(ispin), vxc_rho(ispin)%pw_grid%dvol)
     478           20 :          lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     479           20 :          CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, lri_v_int, calc_force, "LRI_AUX")
     480           40 :          CALL auxbas_pw_pool%give_back_pw(vxc_rho(ispin))
     481              :       END DO
     482           20 :       DEALLOCATE (vxc_rho)
     483           20 :       ekin_mol = -ekin_imol
     484           20 :       xcvirial(1:3, 1:3) = 0.0_dp
     485           20 :       IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) - virial%pv_xc(1:3, 1:3)
     486              : 
     487              :       ! loop over all subsets
     488           60 :       DO isub = 1, kg_env%nsubsets
     489          240 :          atomlist = 0
     490          240 :          DO iatom = 1, natom
     491          200 :             imol = kg_env%atom_to_molecule(iatom)
     492          200 :             color = kg_env%subset_of_mol(imol)
     493          240 :             IF (color == isub) atomlist(iatom) = 1
     494              :          END DO
     495           40 :          CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
     496              : 
     497              :          ekin_imol = 0.0_dp
     498              :          ! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s)
     499              :          CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, &
     500           40 :                             vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol)
     501           40 :          ekin_mol = ekin_mol + ekin_imol
     502              : 
     503           80 :          DO ispin = 1, nspins
     504           40 :             CALL pw_scale(vxc_rho(ispin), -vxc_rho(ispin)%pw_grid%dvol)
     505           40 :             lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     506              :             CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, &
     507              :                                                lri_v_int, calc_force, &
     508           40 :                                                "LRI_AUX", atomlist=atomlist)
     509              :             ! clean up vxc_rho
     510           80 :             CALL auxbas_pw_pool%give_back_pw(vxc_rho(ispin))
     511              :          END DO
     512           40 :          DEALLOCATE (vxc_rho)
     513              : 
     514          100 :          IF (use_virial) THEN
     515            0 :             xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) + virial%pv_xc(1:3, 1:3)
     516              :          END IF
     517              : 
     518              :       END DO
     519              : 
     520           20 :       IF (use_virial) THEN
     521            0 :          virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3)
     522              :       END IF
     523              : 
     524           20 :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
     525           40 :       ALLOCATE (ksmat(1))
     526           40 :       DO ispin = 1, nspins
     527           20 :          lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     528           60 :          DO ikind = 1, nkind
     529        29300 :             CALL para_env%sum(lri_v_int(ikind)%v_int)
     530              :          END DO
     531           20 :          ksmat(1)%matrix => ks_matrix(ispin)%matrix
     532           40 :          CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set)
     533              :       END DO
     534           20 :       IF (calc_force) THEN
     535            2 :          pmat(1:nspins, 1:1) => density_matrix(1:nspins)
     536            2 :          CALL calculate_lri_forces(lri_env, lri_density, qs_env, pmat, atomic_kind_set)
     537              :       END IF
     538           20 :       DEALLOCATE (atomlist, ksmat)
     539              : 
     540              :       ! clean up rho_struct
     541           20 :       CALL qs_rho_unset_rho_ao(rho_struct)
     542           20 :       CALL qs_rho_release(rho_struct)
     543           20 :       DEALLOCATE (rho_struct)
     544              : 
     545           20 :       CALL timestop(handle)
     546              : 
     547           60 :    END SUBROUTINE kg_ekin_embed_lri
     548              : 
     549              : ! **************************************************************************************************
     550              : !> \brief ...
     551              : !> \param qs_env ...
     552              : !> \param kg_env ...
     553              : !> \param ks_matrix ...
     554              : !> \param ekin_mol ...
     555              : !> \param calc_force ...
     556              : !> \param do_kernel Contribution of kinetic energy functional to kernel in response calculation
     557              : !> \param pmat_ext Response density used to fold 2nd deriv or to integrate kinetic energy functional
     558              : ! **************************************************************************************************
     559           20 :    SUBROUTINE kg_ekin_ri_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force, &
     560              :                                do_kernel, pmat_ext)
     561              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     562              :       TYPE(kg_environment_type), POINTER                 :: kg_env
     563              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix
     564              :       REAL(KIND=dp), INTENT(out)                         :: ekin_mol
     565              :       LOGICAL                                            :: calc_force, do_kernel
     566              :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
     567              :          POINTER                                         :: pmat_ext
     568              : 
     569              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'kg_ekin_ri_embed'
     570              : 
     571              :       INTEGER                                            :: color, handle, iatom, ikind, imol, &
     572              :                                                             iounit, ispin, isub, natom, nkind, &
     573              :                                                             nspins
     574           20 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atomlist
     575           20 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     576              :       LOGICAL                                            :: use_virial
     577              :       REAL(KIND=dp)                                      :: alpha, ekin_imol
     578              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: xcvirial
     579           20 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     580              :       TYPE(cp_logger_type), POINTER                      :: logger
     581           20 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: density_matrix, ksmat
     582           20 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: pmat
     583              :       TYPE(dft_control_type), POINTER                    :: dft_control
     584              :       TYPE(lri_density_type), POINTER                    :: lri_density, lri_rho1
     585              :       TYPE(lri_environment_type), POINTER                :: lri_env, lri_env1
     586           20 :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: lri_v_int
     587              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     588           20 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho1_g
     589              :       TYPE(pw_env_type), POINTER                         :: pw_env
     590              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     591           20 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho1_r, tau1_r, vxc_rho, vxc_tau
     592              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     593              :       TYPE(qs_rho_type), POINTER                         :: rho, rho1, rho_struct
     594              :       TYPE(section_vals_type), POINTER                   :: xc_section
     595              :       TYPE(virial_type), POINTER                         :: virial
     596              : 
     597           20 :       CALL timeset(routineN, handle)
     598              : 
     599           20 :       logger => cp_get_default_logger()
     600           20 :       iounit = cp_logger_get_default_unit_nr(logger)
     601              : 
     602              :       CALL get_qs_env(qs_env, &
     603              :                       ks_env=ks_env, &
     604              :                       rho=rho, &
     605              :                       natom=natom, &
     606              :                       nkind=nkind, &
     607              :                       dft_control=dft_control, &
     608              :                       virial=virial, &
     609              :                       para_env=para_env, &
     610           20 :                       pw_env=pw_env)
     611              : 
     612           20 :       nspins = dft_control%nspins
     613           20 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     614            0 :       use_virial = use_virial .AND. calc_force
     615              : 
     616              :       ! Kernel potential in response calculation (no forces calculated at this point)
     617              :       ! requires spin-factor
     618              :       ! alpha = 2 closed-shell
     619              :       ! alpha = 1 open-shell
     620           20 :       alpha = 1.0_dp
     621           20 :       IF (do_kernel .AND. .NOT. calc_force .AND. nspins == 1) alpha = 2.0_dp
     622              : 
     623           20 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     624              : 
     625              :       ! get the density matrix
     626           20 :       CALL qs_rho_get(rho, rho_ao=density_matrix)
     627              :       ! allocate and initialize the density
     628              :       NULLIFY (rho_struct)
     629           20 :       ALLOCATE (rho_struct)
     630           20 :       CALL qs_rho_create(rho_struct)
     631              :       ! set the density matrix to the blocked matrix
     632           20 :       CALL qs_rho_set(rho_struct, rho_ao=density_matrix)
     633           20 :       CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
     634              : 
     635           20 :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
     636           20 :       ALLOCATE (cell_to_index(1, 1, 1))
     637           20 :       cell_to_index(1, 1, 1) = 1
     638           20 :       lri_env => kg_env%lri_env
     639           20 :       lri_density => kg_env%lri_density
     640              : 
     641           20 :       NULLIFY (pmat)
     642          100 :       ALLOCATE (pmat(nspins, 1))
     643           40 :       DO ispin = 1, nspins
     644           40 :          pmat(ispin, 1)%matrix => density_matrix(ispin)%matrix
     645              :       END DO
     646              :       CALL calculate_lri_densities(lri_env, lri_density, qs_env, pmat, cell_to_index, &
     647           20 :                                    rho_struct, atomic_kind_set, para_env, response_density=.FALSE.)
     648           20 :       kg_env%lri_density => lri_density
     649              : 
     650           20 :       DEALLOCATE (pmat)
     651              : 
     652           20 :       IF (PRESENT(pmat_ext)) THEN
     653              :          ! If external density associated then it is needed either for
     654              :          ! 1) folding of second derivative while partially integrating, or
     655              :          ! 2) integration of response forces
     656              :          NULLIFY (rho1)
     657            0 :          ALLOCATE (rho1)
     658            0 :          CALL qs_rho_create(rho1)
     659            0 :          CALL qs_rho_set(rho1, rho_ao=pmat_ext)
     660            0 :          CALL qs_rho_rebuild(rho1, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.)
     661              : 
     662            0 :          lri_env1 => kg_env%lri_env1
     663            0 :          lri_rho1 => kg_env%lri_rho1
     664              :          ! calculate external density as LRI-densities
     665            0 :          NULLIFY (pmat)
     666            0 :          ALLOCATE (pmat(nspins, 1))
     667            0 :          DO ispin = 1, nspins
     668            0 :             pmat(ispin, 1)%matrix => pmat_ext(ispin)%matrix
     669              :          END DO
     670              :          CALL calculate_lri_densities(lri_env1, lri_rho1, qs_env, pmat, cell_to_index, &
     671            0 :                                       rho1, atomic_kind_set, para_env, response_density=.FALSE.)
     672            0 :          kg_env%lri_rho1 => lri_rho1
     673            0 :          DEALLOCATE (pmat)
     674              : 
     675              :       END IF
     676              : 
     677              :       ! XC-section pointing to kinetic energy functional in KG environment
     678              :       NULLIFY (xc_section)
     679           20 :       xc_section => kg_env%xc_section_kg
     680              : 
     681              :       ! full density kinetic energy term
     682           20 :       ekin_imol = 0.0_dp
     683           20 :       NULLIFY (vxc_rho, vxc_tau)
     684              : 
     685              :       ! calculate xc potential or kernel
     686           20 :       IF (do_kernel) THEN
     687              :          ! kernel total
     688              :          ! derivation wrt to rho_struct and evaluation at rho_struct
     689            0 :          CALL qs_rho_get(rho1, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
     690              :          CALL create_kernel(qs_env, &
     691              :                             vxc=vxc_rho, &
     692              :                             vxc_tau=vxc_tau, &
     693              :                             rho=rho_struct, &
     694              :                             rho1_r=rho1_r, &
     695              :                             rho1_g=rho1_g, &
     696              :                             tau1_r=tau1_r, &
     697            0 :                             xc_section=xc_section)
     698              :       ELSE
     699              :          ! vxc total
     700              :          CALL qs_vxc_create(ks_env=ks_env, &
     701              :                             rho_struct=rho_struct, &
     702              :                             xc_section=xc_section, &
     703              :                             vxc_rho=vxc_rho, &
     704              :                             vxc_tau=vxc_tau, &
     705           20 :                             exc=ekin_imol)
     706              : 
     707              :       END IF
     708              : 
     709           20 :       IF (ASSOCIATED(vxc_tau)) THEN
     710            0 :          CPABORT(" KG with meta-kinetic energy functionals not implemented")
     711              :       END IF
     712              : 
     713           40 :       DO ispin = 1, nspins
     714           20 :          CALL pw_scale(vxc_rho(ispin), alpha*vxc_rho(ispin)%pw_grid%dvol)
     715              : 
     716           20 :          IF (PRESENT(pmat_ext) .AND. .NOT. do_kernel) THEN
     717              :             ! int w/ pmat_ext
     718            0 :             lri_v_int => lri_rho1%lri_coefs(ispin)%lri_kinds
     719              :          ELSE
     720              :             ! int w/ rho_ao
     721           20 :             lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     722              :          END IF
     723           20 :          CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, lri_v_int, calc_force, "LRI_AUX")
     724           40 :          CALL auxbas_pw_pool%give_back_pw(vxc_rho(ispin))
     725              :       END DO
     726              : 
     727           20 :       DEALLOCATE (vxc_rho)
     728           20 :       ekin_mol = -ekin_imol
     729           20 :       xcvirial(1:3, 1:3) = 0.0_dp
     730           20 :       IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) - virial%pv_xc(1:3, 1:3)
     731              : 
     732              :       ! loop over all subsets
     733           60 :       ALLOCATE (atomlist(natom))
     734           60 :       DO isub = 1, kg_env%nsubsets
     735          240 :          atomlist = 0
     736          240 :          DO iatom = 1, natom
     737          200 :             imol = kg_env%atom_to_molecule(iatom)
     738          200 :             color = kg_env%subset_of_mol(imol)
     739          240 :             IF (color == isub) atomlist(iatom) = 1
     740              :          END DO
     741              :          ! update ground-state density
     742           40 :          CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist)
     743              : 
     744              :          ! Same for external (response) density if present
     745           40 :          IF (PRESENT(pmat_ext)) THEN
     746              :             ! update response density
     747            0 :             CALL lri_kg_rho_update(rho1, qs_env, lri_env1, lri_rho1, atomlist)
     748              :          END IF
     749              : 
     750           40 :          ekin_imol = 0.0_dp
     751              :          ! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s)
     752           40 :          NULLIFY (vxc_rho, vxc_tau)
     753              : 
     754              :          ! calculate xc potential or kernel
     755           40 :          IF (do_kernel) THEN
     756              :             ! subsys kernel
     757            0 :             CALL qs_rho_get(rho1, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
     758              :             CALL create_kernel(qs_env, &
     759              :                                vxc=vxc_rho, &
     760              :                                vxc_tau=vxc_tau, &
     761              :                                rho=rho_struct, &
     762              :                                rho1_r=rho1_r, &
     763              :                                rho1_g=rho1_g, &
     764              :                                tau1_r=tau1_r, &
     765            0 :                                xc_section=xc_section)
     766              :          ELSE
     767              : 
     768              :             ! subsys xc-potential
     769              :             CALL qs_vxc_create(ks_env=ks_env, &
     770              :                                rho_struct=rho_struct, &
     771              :                                xc_section=xc_section, &
     772              :                                vxc_rho=vxc_rho, &
     773              :                                vxc_tau=vxc_tau, &
     774           40 :                                exc=ekin_imol)
     775              :          END IF
     776           40 :          ekin_mol = ekin_mol + ekin_imol
     777              : 
     778           80 :          DO ispin = 1, nspins
     779           40 :             CALL pw_scale(vxc_rho(ispin), -alpha*vxc_rho(ispin)%pw_grid%dvol)
     780              : 
     781           40 :             IF (PRESENT(pmat_ext) .AND. .NOT. do_kernel) THEN
     782              :                ! int w/ pmat_ext
     783            0 :                lri_v_int => lri_rho1%lri_coefs(ispin)%lri_kinds
     784              :             ELSE
     785              :                ! int w/ rho_ao
     786           40 :                lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     787              :             END IF
     788              : 
     789              :             CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, &
     790              :                                                lri_v_int, calc_force, &
     791           40 :                                                "LRI_AUX", atomlist=atomlist)
     792              :             ! clean up vxc_rho
     793           80 :             CALL auxbas_pw_pool%give_back_pw(vxc_rho(ispin))
     794              :          END DO
     795           40 :          DEALLOCATE (vxc_rho)
     796              : 
     797           60 :          IF (use_virial) THEN
     798            0 :             xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) + virial%pv_xc(1:3, 1:3)
     799              :          END IF
     800              : 
     801              :       END DO
     802              : 
     803           20 :       IF (use_virial) THEN
     804            0 :          virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3)
     805              :       END IF
     806              : 
     807           40 :       ALLOCATE (ksmat(1))
     808           40 :       DO ispin = 1, nspins
     809           20 :          ksmat(1)%matrix => ks_matrix(ispin)%matrix
     810           40 :          IF (PRESENT(pmat_ext) .AND. .NOT. do_kernel) THEN
     811              :             ! KS int with rho_ext"
     812            0 :             lri_v_int => lri_rho1%lri_coefs(ispin)%lri_kinds
     813            0 :             DO ikind = 1, nkind
     814            0 :                CALL para_env%sum(lri_v_int(ikind)%v_int)
     815              :             END DO
     816            0 :             CALL calculate_lri_ks_matrix(lri_env1, lri_v_int, ksmat, atomic_kind_set)
     817              :          ELSE
     818              :             ! KS int with rho_ao"
     819           20 :             lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
     820           60 :             DO ikind = 1, nkind
     821        29300 :                CALL para_env%sum(lri_v_int(ikind)%v_int)
     822              :             END DO
     823           20 :             CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set)
     824              :          END IF
     825              : 
     826              :       END DO
     827           20 :       IF (calc_force) THEN
     828              : 
     829            2 :          NULLIFY (pmat)
     830            8 :          ALLOCATE (pmat(nspins, 1))
     831              : 
     832            2 :          IF (PRESENT(pmat_ext) .AND. .NOT. do_kernel) THEN
     833              :             ! Forces with rho_ext
     834            0 :             DO ispin = 1, nspins
     835            0 :                pmat(ispin, 1)%matrix => pmat_ext(ispin)%matrix
     836              :             END DO
     837            0 :             CALL calculate_lri_forces(lri_env1, lri_rho1, qs_env, pmat, atomic_kind_set)
     838              :          ELSE
     839              :             ! Forces with rho_ao
     840            4 :             DO ispin = 1, nspins
     841            4 :                pmat(ispin, 1)%matrix => density_matrix(ispin)%matrix
     842              :             END DO
     843            2 :             CALL calculate_lri_forces(lri_env, lri_density, qs_env, pmat, atomic_kind_set)
     844              :          END IF
     845              : 
     846            2 :          DEALLOCATE (pmat)
     847              : 
     848              :       END IF
     849           20 :       DEALLOCATE (atomlist, ksmat)
     850              : 
     851              :       ! clean up rho_struct
     852           20 :       CALL qs_rho_unset_rho_ao(rho_struct)
     853           20 :       CALL qs_rho_release(rho_struct)
     854           20 :       DEALLOCATE (rho_struct)
     855           20 :       IF (PRESENT(pmat_ext)) THEN
     856            0 :          CALL qs_rho_unset_rho_ao(rho1)
     857            0 :          CALL qs_rho_release(rho1)
     858            0 :          DEALLOCATE (rho1)
     859              :       END IF
     860           20 :       DEALLOCATE (cell_to_index)
     861              : 
     862           20 :       CALL timestop(handle)
     863              : 
     864           40 :    END SUBROUTINE kg_ekin_ri_embed
     865              : 
     866              : ! **************************************************************************************************
     867              : !> \brief ...
     868              : !> \param qs_env ...
     869              : !> \param ks_matrix ...
     870              : !> \param ekin_mol ...
     871              : ! **************************************************************************************************
     872          144 :    SUBROUTINE kg_ekin_atomic(qs_env, ks_matrix, ekin_mol)
     873              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     874              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix
     875              :       REAL(KIND=dp), INTENT(out)                         :: ekin_mol
     876              : 
     877              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'kg_ekin_atomic'
     878              : 
     879              :       INTEGER                                            :: handle, ispin, nspins
     880          144 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: density_matrix, tnadd_matrix
     881              :       TYPE(kg_environment_type), POINTER                 :: kg_env
     882              :       TYPE(qs_rho_type), POINTER                         :: rho
     883              : 
     884          144 :       NULLIFY (rho, kg_env, density_matrix, tnadd_matrix)
     885              : 
     886          144 :       CALL timeset(routineN, handle)
     887          144 :       CALL get_qs_env(qs_env, kg_env=kg_env, rho=rho)
     888              : 
     889          144 :       nspins = SIZE(ks_matrix)
     890              :       ! get the density matrix
     891          144 :       CALL qs_rho_get(rho, rho_ao=density_matrix)
     892              :       ! get the tnadd matrix
     893          144 :       tnadd_matrix => kg_env%tnadd_mat
     894              : 
     895          144 :       ekin_mol = 0.0_dp
     896          288 :       DO ispin = 1, nspins
     897          144 :          CALL dbcsr_dot(tnadd_matrix(1)%matrix, density_matrix(ispin)%matrix, ekin_mol)
     898              :          CALL dbcsr_add(ks_matrix(ispin)%matrix, tnadd_matrix(1)%matrix, &
     899          288 :                         alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     900              :       END DO
     901              :       ! definition is inverted (see qs_ks_methods)
     902          144 :       ekin_mol = -ekin_mol
     903              : 
     904          144 :       CALL timestop(handle)
     905              : 
     906          144 :    END SUBROUTINE kg_ekin_atomic
     907              : 
     908              : END MODULE kg_correction
        

Generated by: LCOV version 2.0-1