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