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