LCOV - code coverage report
Current view: top level - src - kg_correction.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 294 338 87.0 %
Date: 2024-03-28 07:31:50 Functions: 5 5 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 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_log_handling,                 ONLY: cp_get_default_logger,&
      18             :                                               cp_logger_get_default_unit_nr,&
      19             :                                               cp_logger_type
      20             :    USE dbcsr_api,                       ONLY: dbcsr_add,&
      21             :                                               dbcsr_dot,&
      22             :                                               dbcsr_p_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        1190 :                                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          40 :                                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         624 :    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          40 :    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          10 :          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          20 :    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 1.15