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 module that builds the second order perturbation kernel
10 : !> kpp1 = delta_rho|_P delta_rho|_P E drho(P1) drho
11 : !> \par History
12 : !> 07.2002 created [fawzi]
13 : !> \author Fawzi Mohamed
14 : ! **************************************************************************************************
15 : MODULE qs_kpp1_env_methods
16 : USE admm_types, ONLY: admm_type,&
17 : get_admm_env
18 : USE atomic_kind_types, ONLY: atomic_kind_type
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
21 : dbcsr_copy,&
22 : dbcsr_p_type,&
23 : dbcsr_set
24 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_type,&
27 : cp_to_string
28 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 : cp_print_key_should_output,&
30 : cp_print_key_unit_nr
31 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
32 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
33 : do_method_gapw,&
34 : do_method_gapw_xc
35 : USE input_section_types, ONLY: section_get_ival,&
36 : section_vals_get_subs_vals,&
37 : section_vals_type
38 : USE kahan_sum, ONLY: accurate_sum
39 : USE kinds, ONLY: dp
40 : USE lri_environment_types, ONLY: lri_density_type,&
41 : lri_environment_type,&
42 : lri_kind_type
43 : USE lri_ks_methods, ONLY: calculate_lri_ks_matrix
44 : USE message_passing, ONLY: mp_para_env_type
45 : USE pw_env_types, ONLY: pw_env_get,&
46 : pw_env_type
47 : USE pw_methods, ONLY: pw_axpy,&
48 : pw_copy,&
49 : pw_integrate_function,&
50 : pw_scale,&
51 : pw_transfer
52 : USE pw_poisson_methods, ONLY: pw_poisson_solve
53 : USE pw_poisson_types, ONLY: pw_poisson_type
54 : USE pw_pool_types, ONLY: pw_pool_type
55 : USE pw_types, ONLY: pw_c1d_gs_type,&
56 : pw_r3d_rs_type
57 : USE qs_environment_types, ONLY: get_qs_env,&
58 : qs_environment_type
59 : USE qs_gapw_densities, ONLY: prepare_gapw_den
60 : USE qs_integrate_potential, ONLY: integrate_v_rspace,&
61 : integrate_v_rspace_diagonal,&
62 : integrate_v_rspace_one_center
63 : USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type
64 : USE qs_ks_atom, ONLY: update_ks_atom
65 : USE qs_p_env_types, ONLY: qs_p_env_type
66 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
67 : USE qs_rho_atom_types, ONLY: rho_atom_type
68 : USE qs_rho_types, ONLY: qs_rho_get,&
69 : qs_rho_type
70 : USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom
71 : USE xc, ONLY: xc_calc_2nd_deriv,&
72 : xc_prep_2nd_deriv
73 : USE xc_derivative_set_types, ONLY: xc_dset_release
74 : USE xc_rho_set_types, ONLY: xc_rho_set_release
75 : #include "./base/base_uses.f90"
76 :
77 : IMPLICIT NONE
78 :
79 : PRIVATE
80 :
81 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
82 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kpp1_env_methods'
83 :
84 : PUBLIC :: kpp1_create, &
85 : kpp1_did_change, &
86 : calc_kpp1
87 :
88 : CONTAINS
89 :
90 : ! **************************************************************************************************
91 : !> \brief allocates and initializes a kpp1_env
92 : !> \param kpp1_env the environment to initialize
93 : !> \par History
94 : !> 07.2002 created [fawzi]
95 : !> \author Fawzi Mohamed
96 : ! **************************************************************************************************
97 1830 : SUBROUTINE kpp1_create(kpp1_env)
98 : TYPE(qs_kpp1_env_type) :: kpp1_env
99 :
100 1830 : NULLIFY (kpp1_env%v_ao, kpp1_env%rho_set, kpp1_env%deriv_set, &
101 1830 : kpp1_env%rho_set_admm, kpp1_env%deriv_set_admm)
102 1830 : END SUBROUTINE kpp1_create
103 :
104 : ! **************************************************************************************************
105 : !> \brief ...
106 : !> \param rho1_xc ...
107 : !> \param rho1 ...
108 : !> \param xc_section ...
109 : !> \param lrigpw ...
110 : !> \param do_triplet ...
111 : !> \param qs_env ...
112 : !> \param p_env ...
113 : !> \param calc_forces ...
114 : !> \param calc_virial ...
115 : !> \param virial ...
116 : ! **************************************************************************************************
117 1784 : SUBROUTINE calc_kpp1(rho1_xc, rho1, xc_section, lrigpw, do_triplet, qs_env, p_env, &
118 : calc_forces, calc_virial, virial)
119 :
120 : TYPE(qs_rho_type), POINTER :: rho1_xc, rho1
121 : TYPE(section_vals_type), POINTER :: xc_section
122 : LOGICAL, INTENT(IN) :: lrigpw, do_triplet
123 : TYPE(qs_environment_type), POINTER :: qs_env
124 : TYPE(qs_p_env_type) :: p_env
125 : LOGICAL, INTENT(IN), OPTIONAL :: calc_forces, calc_virial
126 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
127 : OPTIONAL :: virial
128 :
129 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_kpp1'
130 :
131 : INTEGER :: handle, ikind, ispin, nkind, ns, nspins, &
132 : output_unit
133 : LOGICAL :: gapw, gapw_xc, lsd, my_calc_forces
134 : REAL(KIND=dp) :: alpha, energy_hartree, energy_hartree_1c
135 1784 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
136 : TYPE(cp_logger_type), POINTER :: logger
137 1784 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: k1mat, rho_ao
138 1784 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
139 : TYPE(lri_density_type), POINTER :: lri_density
140 : TYPE(lri_environment_type), POINTER :: lri_env
141 1784 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
142 : TYPE(mp_para_env_type), POINTER :: para_env
143 : TYPE(pw_c1d_gs_type) :: rho1_tot_gspace
144 1784 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_pw
145 : TYPE(pw_env_type), POINTER :: pw_env
146 : TYPE(pw_poisson_type), POINTER :: poisson_env
147 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
148 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
149 1784 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho1_r_pw, tau1_r, tau1_r_pw, &
150 1784 : v_rspace_new, v_xc, v_xc_tau
151 : TYPE(pw_r3d_rs_type), POINTER :: weights
152 : TYPE(qs_rho_type), POINTER :: rho
153 1784 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
154 : TYPE(section_vals_type), POINTER :: input, scf_section
155 :
156 1784 : CALL timeset(routineN, handle)
157 :
158 1784 : NULLIFY (v_xc, rho1_g, pw_env, rho1_g_pw, tau1_r_pw)
159 1784 : logger => cp_get_default_logger()
160 :
161 1784 : CPASSERT(ASSOCIATED(p_env%kpp1))
162 1784 : CPASSERT(ASSOCIATED(p_env%kpp1_env))
163 1784 : CPASSERT(ASSOCIATED(rho1))
164 :
165 1784 : nspins = SIZE(p_env%kpp1)
166 1784 : lsd = (nspins == 2)
167 :
168 1784 : my_calc_forces = .FALSE.
169 1784 : IF (PRESENT(calc_forces)) my_calc_forces = calc_forces
170 :
171 : CALL get_qs_env(qs_env, &
172 : pw_env=pw_env, &
173 : input=input, &
174 : para_env=para_env, &
175 1784 : rho=rho)
176 :
177 1784 : CPASSERT(ASSOCIATED(rho1))
178 :
179 1784 : IF (lrigpw) THEN
180 : CALL get_qs_env(qs_env, &
181 : lri_env=lri_env, &
182 : lri_density=lri_density, &
183 0 : atomic_kind_set=atomic_kind_set)
184 : END IF
185 :
186 1784 : gapw = (section_get_ival(input, "DFT%QS%METHOD") == do_method_gapw)
187 1784 : gapw_xc = (section_get_ival(input, "DFT%QS%METHOD") == do_method_gapw_xc)
188 1784 : IF (gapw_xc) THEN
189 0 : CPASSERT(ASSOCIATED(rho1_xc))
190 : END IF
191 :
192 1784 : CALL kpp1_check_i_alloc(p_env%kpp1_env, qs_env, do_triplet)
193 :
194 1784 : CALL qs_rho_get(rho, rho_ao=rho_ao)
195 1784 : CALL qs_rho_get(rho1, rho_g=rho1_g)
196 :
197 : ! gets the tmp grids
198 1784 : CPASSERT(ASSOCIATED(pw_env))
199 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
200 1784 : poisson_env=poisson_env)
201 1784 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
202 :
203 1784 : IF (gapw .OR. gapw_xc) THEN
204 0 : CALL prepare_gapw_den(qs_env, p_env%local_rho_set, do_rho0=(.NOT. gapw_xc))
205 : END IF
206 :
207 : ! *** calculate the hartree potential on the total density ***
208 1784 : CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
209 :
210 1784 : CALL pw_copy(rho1_g(1), rho1_tot_gspace)
211 2334 : DO ispin = 2, nspins
212 2334 : CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
213 : END DO
214 1784 : IF (gapw) THEN
215 0 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
216 0 : IF (ASSOCIATED(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
217 0 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho1_tot_gspace)
218 : END IF
219 : END IF
220 :
221 1784 : scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
222 1784 : IF (cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES") &
223 : /= 0) THEN
224 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
225 0 : extension=".scfLog")
226 0 : CALL print_densities(rho1, rho1_tot_gspace, output_unit)
227 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
228 0 : "PRINT%TOTAL_DENSITIES")
229 : END IF
230 :
231 1784 : IF (.NOT. (nspins == 1 .AND. do_triplet)) THEN
232 : BLOCK
233 : TYPE(pw_c1d_gs_type) :: v_hartree_gspace
234 1784 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
235 : CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
236 : energy_hartree, &
237 1784 : v_hartree_gspace)
238 1784 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
239 1784 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
240 : END BLOCK
241 3568 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
242 : END IF
243 :
244 1784 : CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
245 :
246 : ! *** calculate the xc potential ***
247 1784 : IF (gapw_xc) THEN
248 0 : CALL qs_rho_get(rho1_xc, rho_r=rho1_r, tau_r=tau1_r)
249 : ELSE
250 1784 : CALL qs_rho_get(rho1, rho_r=rho1_r, tau_r=tau1_r)
251 : END IF
252 :
253 1784 : IF (nspins == 1 .AND. do_triplet) THEN
254 :
255 0 : lsd = .TRUE.
256 0 : ALLOCATE (rho1_r_pw(2))
257 0 : DO ispin = 1, 2
258 0 : CALL rho1_r_pw(ispin)%create(rho1_r(1)%pw_grid)
259 0 : CALL pw_transfer(rho1_r(1), rho1_r_pw(ispin))
260 : END DO
261 :
262 0 : IF (ASSOCIATED(tau1_r)) THEN
263 0 : ALLOCATE (tau1_r_pw(2))
264 0 : DO ispin = 1, 2
265 0 : CALL tau1_r_pw(ispin)%create(tau1_r(1)%pw_grid)
266 0 : CALL pw_transfer(tau1_r(1), tau1_r_pw(ispin))
267 : END DO
268 : END IF
269 :
270 : ELSE
271 :
272 1784 : rho1_r_pw => rho1_r
273 :
274 1784 : tau1_r_pw => tau1_r
275 :
276 : END IF
277 :
278 1784 : NULLIFY (weights)
279 1784 : CALL get_qs_env(qs_env, xcint_weights=weights)
280 :
281 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, p_env%kpp1_env%deriv_set, p_env%kpp1_env%rho_set, &
282 : rho1_r_pw, rho1_g_pw, tau1_r_pw, auxbas_pw_pool, weights, &
283 : xc_section, .FALSE., do_excitations=.TRUE., do_triplet=do_triplet, &
284 1784 : compute_virial=calc_virial, virial_xc=virial)
285 :
286 4118 : DO ispin = 1, nspins
287 4118 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
288 : END DO
289 1784 : v_rspace_new => v_xc
290 1784 : IF (SIZE(v_xc) /= nspins) THEN
291 0 : CALL auxbas_pw_pool%give_back_pw(v_xc(2))
292 : END IF
293 1784 : NULLIFY (v_xc)
294 1784 : IF (ASSOCIATED(v_xc_tau)) THEN
295 616 : DO ispin = 1, nspins
296 616 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
297 : END DO
298 244 : IF (SIZE(v_xc_tau) /= nspins) THEN
299 0 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(2))
300 : END IF
301 : END IF
302 :
303 1784 : IF (gapw .OR. gapw_xc) THEN
304 0 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
305 0 : rho1_atom_set => p_env%local_rho_set%rho_atom_set
306 : CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
307 0 : do_triplet=do_triplet)
308 : END IF
309 :
310 1784 : IF (nspins == 1 .AND. do_triplet) THEN
311 0 : DO ispin = 1, SIZE(rho1_r_pw)
312 0 : CALL rho1_r_pw(ispin)%release()
313 : END DO
314 0 : DEALLOCATE (rho1_r_pw)
315 0 : IF (ASSOCIATED(tau1_r_pw)) THEN
316 0 : DO ispin = 1, SIZE(tau1_r_pw)
317 0 : CALL tau1_r_pw(ispin)%release()
318 : END DO
319 0 : DEALLOCATE (tau1_r_pw)
320 : END IF
321 : END IF
322 :
323 550 : alpha = 1.0_dp
324 1234 : IF (nspins == 1) alpha = 2.0_dp
325 :
326 : !-------------------------------!
327 : ! Add both hartree and xc terms !
328 : !-------------------------------!
329 4118 : DO ispin = 1, nspins
330 2334 : CALL dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix, 0.0_dp)
331 :
332 : ! XC and Hartree are integrated separatedly
333 : ! XC uses the soft basis set only
334 2334 : IF (gapw_xc) THEN
335 :
336 0 : IF (nspins == 1) THEN
337 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
338 : pmat=rho_ao(ispin), &
339 : hmat=p_env%kpp1_env%v_ao(ispin), &
340 : qs_env=qs_env, &
341 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
342 :
343 0 : IF (ASSOCIATED(v_xc_tau)) THEN
344 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
345 : pmat=rho_ao(ispin), &
346 : hmat=p_env%kpp1_env%v_ao(ispin), &
347 : qs_env=qs_env, &
348 : compute_tau=.TRUE., &
349 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
350 : END IF
351 :
352 : ! add hartree only for SINGLETS
353 0 : IF (.NOT. do_triplet) THEN
354 0 : CALL pw_copy(v_hartree_rspace, v_rspace_new(1))
355 :
356 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
357 : pmat=rho_ao(ispin), &
358 : hmat=p_env%kpp1_env%v_ao(ispin), &
359 : qs_env=qs_env, &
360 0 : calculate_forces=my_calc_forces, gapw=gapw)
361 : END IF
362 : ELSE
363 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
364 : pmat=rho_ao(ispin), &
365 : hmat=p_env%kpp1_env%v_ao(ispin), &
366 : qs_env=qs_env, &
367 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
368 :
369 0 : IF (ASSOCIATED(v_xc_tau)) THEN
370 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
371 : pmat=rho_ao(ispin), &
372 : hmat=p_env%kpp1_env%v_ao(ispin), &
373 : qs_env=qs_env, &
374 : compute_tau=.TRUE., &
375 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
376 : END IF
377 :
378 0 : CALL pw_copy(v_hartree_rspace, v_rspace_new(ispin))
379 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
380 : pmat=rho_ao(ispin), &
381 : hmat=p_env%kpp1_env%v_ao(ispin), &
382 : qs_env=qs_env, &
383 0 : calculate_forces=my_calc_forces, gapw=gapw)
384 : END IF
385 :
386 : ELSE
387 :
388 2334 : IF (nspins == 1) THEN
389 :
390 : ! add hartree only for SINGLETS
391 1234 : IF (.NOT. do_triplet) THEN
392 1234 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(1))
393 : END IF
394 : ELSE
395 1100 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin))
396 : END IF
397 :
398 2334 : IF (lrigpw) THEN
399 0 : IF (ASSOCIATED(v_xc_tau)) CPABORT("Meta-GGA functionals not supported with LRI!")
400 :
401 0 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
402 0 : CALL get_qs_env(qs_env, nkind=nkind)
403 0 : DO ikind = 1, nkind
404 0 : lri_v_int(ikind)%v_int = 0.0_dp
405 : END DO
406 : CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
407 0 : lri_v_int, .FALSE., "LRI_AUX")
408 0 : DO ikind = 1, nkind
409 0 : CALL para_env%sum(lri_v_int(ikind)%v_int)
410 : END DO
411 0 : ALLOCATE (k1mat(1))
412 0 : k1mat(1)%matrix => p_env%kpp1_env%v_ao(ispin)%matrix
413 0 : IF (lri_env%exact_1c_terms) THEN
414 : CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, &
415 0 : rho_ao(ispin)%matrix, qs_env, my_calc_forces, "ORB")
416 : END IF
417 0 : CALL calculate_lri_ks_matrix(lri_env, lri_v_int, k1mat, atomic_kind_set)
418 0 : DEALLOCATE (k1mat)
419 : ELSE
420 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
421 : pmat=rho_ao(ispin), &
422 : hmat=p_env%kpp1_env%v_ao(ispin), &
423 : qs_env=qs_env, &
424 2334 : calculate_forces=my_calc_forces, gapw=gapw)
425 :
426 2334 : IF (ASSOCIATED(v_xc_tau)) THEN
427 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
428 : pmat=rho_ao(ispin), &
429 : hmat=p_env%kpp1_env%v_ao(ispin), &
430 : qs_env=qs_env, &
431 : compute_tau=.TRUE., &
432 372 : calculate_forces=my_calc_forces, gapw=gapw)
433 : END IF
434 : END IF
435 : END IF
436 :
437 4118 : CALL dbcsr_add(p_env%kpp1(ispin)%matrix, p_env%kpp1_env%v_ao(ispin)%matrix, 1.0_dp, alpha)
438 : END DO
439 :
440 1784 : IF (gapw) THEN
441 0 : IF (.NOT. (nspins == 1 .AND. do_triplet)) THEN
442 : CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
443 : p_env%hartree_local%ecoul_1c, &
444 : p_env%local_rho_set, &
445 0 : para_env, tddft=.TRUE., core_2nd=.TRUE.)
446 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
447 : calculate_forces=my_calc_forces, &
448 0 : local_rho_set=p_env%local_rho_set)
449 : END IF
450 : ! *** Add single atom contributions to the KS matrix ***
451 : ! remap pointer
452 0 : ns = SIZE(p_env%kpp1)
453 0 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
454 0 : ns = SIZE(rho_ao)
455 0 : psmat(1:ns, 1:1) => rho_ao(1:ns)
456 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=my_calc_forces, tddft=.TRUE., &
457 0 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
458 1784 : ELSE IF (gapw_xc) THEN
459 0 : ns = SIZE(p_env%kpp1)
460 0 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
461 0 : ns = SIZE(rho_ao)
462 0 : psmat(1:ns, 1:1) => rho_ao(1:ns)
463 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=my_calc_forces, tddft=.TRUE., &
464 0 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
465 : END IF
466 :
467 1784 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
468 4118 : DO ispin = 1, SIZE(v_rspace_new)
469 4118 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
470 : END DO
471 1784 : DEALLOCATE (v_rspace_new)
472 1784 : IF (ASSOCIATED(v_xc_tau)) THEN
473 616 : DO ispin = 1, SIZE(v_xc_tau)
474 616 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
475 : END DO
476 244 : DEALLOCATE (v_xc_tau)
477 : END IF
478 :
479 1784 : CALL timestop(handle)
480 1784 : END SUBROUTINE calc_kpp1
481 :
482 : ! **************************************************************************************************
483 : !> \brief checks that the intenal storage is allocated, and allocs it if needed
484 : !> \param kpp1_env the environment to check
485 : !> \param qs_env the qs environment this kpp1_env lives in
486 : !> \param do_triplet ...
487 : !> \author Fawzi Mohamed
488 : !> \note
489 : !> private routine
490 : ! **************************************************************************************************
491 1784 : SUBROUTINE kpp1_check_i_alloc(kpp1_env, qs_env, do_triplet)
492 :
493 : TYPE(qs_kpp1_env_type) :: kpp1_env
494 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
495 : LOGICAL, INTENT(IN) :: do_triplet
496 :
497 : INTEGER :: ispin, nspins
498 : TYPE(admm_type), POINTER :: admm_env
499 1784 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
500 : TYPE(dft_control_type), POINTER :: dft_control
501 : TYPE(pw_env_type), POINTER :: pw_env
502 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
503 1784 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_rho_r, my_tau_r, rho_r, tau_r
504 : TYPE(pw_r3d_rs_type), POINTER :: weights
505 : TYPE(qs_rho_type), POINTER :: rho
506 : TYPE(section_vals_type), POINTER :: admm_xc_section, input, xc_section
507 :
508 : ! ------------------------------------------------------------------
509 :
510 1784 : NULLIFY (pw_env, auxbas_pw_pool, matrix_s, rho, rho_r, admm_env, dft_control, my_rho_r, my_tau_r)
511 :
512 : CALL get_qs_env(qs_env, pw_env=pw_env, &
513 : matrix_s=matrix_s, rho=rho, input=input, &
514 1784 : admm_env=admm_env, dft_control=dft_control)
515 :
516 1784 : NULLIFY (weights)
517 1784 : CALL get_qs_env(qs_env, xcint_weights=weights)
518 :
519 1784 : CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
520 1784 : nspins = SIZE(rho_r)
521 :
522 1784 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
523 :
524 1784 : IF (.NOT. ASSOCIATED(kpp1_env%v_ao)) THEN
525 272 : CALL dbcsr_allocate_matrix_set(kpp1_env%v_ao, nspins)
526 630 : DO ispin = 1, nspins
527 358 : ALLOCATE (kpp1_env%v_ao(ispin)%matrix)
528 : CALL dbcsr_copy(kpp1_env%v_ao(ispin)%matrix, matrix_s(1)%matrix, &
529 630 : name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin)))
530 : END DO
531 : END IF
532 :
533 1784 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set)) THEN
534 :
535 272 : IF (nspins == 1 .AND. do_triplet) THEN
536 0 : ALLOCATE (my_rho_r(2))
537 0 : DO ispin = 1, 2
538 0 : CALL auxbas_pw_pool%create_pw(my_rho_r(ispin))
539 0 : CALL pw_axpy(rho_r(1), my_rho_r(ispin), 0.5_dp, 0.0_dp)
540 : END DO
541 0 : IF (dft_control%use_kinetic_energy_density) THEN
542 0 : ALLOCATE (my_tau_r(2))
543 0 : DO ispin = 1, 2
544 0 : CALL auxbas_pw_pool%create_pw(my_tau_r(ispin))
545 0 : CALL pw_axpy(tau_r(1), my_tau_r(ispin), 0.5_dp, 0.0_dp)
546 : END DO
547 : END IF
548 : ELSE
549 272 : my_rho_r => rho_r
550 272 : IF (dft_control%use_kinetic_energy_density) THEN
551 40 : my_tau_r => tau_r
552 : END IF
553 : END IF
554 :
555 272 : IF (dft_control%do_admm) THEN
556 40 : xc_section => admm_env%xc_section_primary
557 : ELSE
558 232 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
559 : END IF
560 :
561 6256 : ALLOCATE (kpp1_env%deriv_set, kpp1_env%rho_set)
562 : CALL xc_prep_2nd_deriv(kpp1_env%deriv_set, kpp1_env%rho_set, &
563 : my_rho_r, auxbas_pw_pool, weights, &
564 272 : xc_section=xc_section, tau_r=my_tau_r)
565 :
566 272 : IF (nspins == 1 .AND. do_triplet) THEN
567 0 : DO ispin = 1, SIZE(my_rho_r)
568 0 : CALL my_rho_r(ispin)%release()
569 : END DO
570 0 : DEALLOCATE (my_rho_r)
571 0 : IF (ASSOCIATED(my_tau_r)) THEN
572 0 : DO ispin = 1, SIZE(my_tau_r)
573 0 : CALL my_tau_r(ispin)%release()
574 : END DO
575 0 : DEALLOCATE (my_tau_r)
576 : END IF
577 : END IF
578 : END IF
579 :
580 : ! ADMM Correction
581 1784 : IF (dft_control%do_admm) THEN
582 212 : IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
583 92 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
584 24 : CPASSERT(.NOT. do_triplet)
585 24 : admm_xc_section => admm_env%xc_section_aux
586 24 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho)
587 24 : CALL qs_rho_get(rho, rho_r=rho_r)
588 552 : ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
589 : CALL xc_prep_2nd_deriv(kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, &
590 : rho_r, auxbas_pw_pool, weights, &
591 24 : xc_section=admm_xc_section)
592 : END IF
593 : END IF
594 : END IF
595 :
596 1784 : END SUBROUTINE kpp1_check_i_alloc
597 :
598 : ! **************************************************************************************************
599 : !> \brief function to advise of changes either in the grids
600 : !> \param kpp1_env the kpp1_env
601 : !> \par History
602 : !> 11.2002 created [fawzi]
603 : !> \author Fawzi Mohamed
604 : ! **************************************************************************************************
605 1830 : SUBROUTINE kpp1_did_change(kpp1_env)
606 : TYPE(qs_kpp1_env_type) :: kpp1_env
607 :
608 1830 : IF (ASSOCIATED(kpp1_env%deriv_set)) THEN
609 0 : CALL xc_dset_release(kpp1_env%deriv_set)
610 0 : DEALLOCATE (kpp1_env%deriv_set)
611 : NULLIFY (kpp1_env%deriv_set)
612 : END IF
613 1830 : IF (ASSOCIATED(kpp1_env%rho_set)) THEN
614 0 : CALL xc_rho_set_release(kpp1_env%rho_set)
615 0 : DEALLOCATE (kpp1_env%rho_set)
616 : END IF
617 :
618 1830 : END SUBROUTINE kpp1_did_change
619 :
620 : ! **************************************************************************************************
621 : !> \brief ...
622 : !> \param rho1 ...
623 : !> \param rho1_tot_gspace ...
624 : !> \param out_unit ...
625 : ! **************************************************************************************************
626 0 : SUBROUTINE print_densities(rho1, rho1_tot_gspace, out_unit)
627 :
628 : TYPE(qs_rho_type), POINTER :: rho1
629 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho1_tot_gspace
630 : INTEGER :: out_unit
631 :
632 : REAL(KIND=dp) :: total_rho_gspace
633 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho1_r
634 :
635 0 : NULLIFY (tot_rho1_r)
636 :
637 0 : total_rho_gspace = pw_integrate_function(rho1_tot_gspace, isign=-1)
638 0 : IF (out_unit > 0) THEN
639 0 : CALL qs_rho_get(rho1, tot_rho_r=tot_rho1_r)
640 : WRITE (UNIT=out_unit, FMT="(T3,A,T60,F20.10)") &
641 0 : "KPP1 total charge density (r-space):", &
642 0 : accurate_sum(tot_rho1_r), &
643 0 : "KPP1 total charge density (g-space):", &
644 0 : total_rho_gspace
645 : END IF
646 :
647 0 : END SUBROUTINE print_densities
648 :
649 : END MODULE qs_kpp1_env_methods
|
