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