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