Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 1688 : SUBROUTINE kpp1_create(kpp1_env)
98 : TYPE(qs_kpp1_env_type) :: kpp1_env
99 :
100 1688 : NULLIFY (kpp1_env%v_ao, kpp1_env%rho_set, kpp1_env%deriv_set, &
101 1688 : kpp1_env%rho_set_admm, kpp1_env%deriv_set_admm)
102 1688 : 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(qs_rho_type), POINTER :: rho
152 1784 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
153 : TYPE(section_vals_type), POINTER :: input, scf_section
154 :
155 1784 : CALL timeset(routineN, handle)
156 :
157 1784 : NULLIFY (v_xc, rho1_g, pw_env, rho1_g_pw, tau1_r_pw)
158 1784 : logger => cp_get_default_logger()
159 :
160 1784 : CPASSERT(ASSOCIATED(p_env%kpp1))
161 1784 : CPASSERT(ASSOCIATED(p_env%kpp1_env))
162 1784 : CPASSERT(ASSOCIATED(rho1))
163 :
164 1784 : nspins = SIZE(p_env%kpp1)
165 1784 : lsd = (nspins == 2)
166 :
167 1784 : my_calc_forces = .FALSE.
168 1784 : IF (PRESENT(calc_forces)) my_calc_forces = calc_forces
169 :
170 : CALL get_qs_env(qs_env, &
171 : pw_env=pw_env, &
172 : input=input, &
173 : para_env=para_env, &
174 1784 : rho=rho)
175 :
176 1784 : CPASSERT(ASSOCIATED(rho1))
177 :
178 1784 : IF (lrigpw) THEN
179 : CALL get_qs_env(qs_env, &
180 : lri_env=lri_env, &
181 : lri_density=lri_density, &
182 0 : atomic_kind_set=atomic_kind_set)
183 : END IF
184 :
185 1784 : gapw = (section_get_ival(input, "DFT%QS%METHOD") == do_method_gapw)
186 1784 : gapw_xc = (section_get_ival(input, "DFT%QS%METHOD") == do_method_gapw_xc)
187 1784 : IF (gapw_xc) THEN
188 0 : CPASSERT(ASSOCIATED(rho1_xc))
189 : END IF
190 :
191 1784 : CALL kpp1_check_i_alloc(p_env%kpp1_env, qs_env, do_triplet)
192 :
193 1784 : CALL qs_rho_get(rho, rho_ao=rho_ao)
194 1784 : CALL qs_rho_get(rho1, rho_g=rho1_g)
195 :
196 : ! gets the tmp grids
197 1784 : CPASSERT(ASSOCIATED(pw_env))
198 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
199 1784 : poisson_env=poisson_env)
200 1784 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
201 :
202 1784 : IF (gapw .OR. gapw_xc) &
203 0 : CALL prepare_gapw_den(qs_env, p_env%local_rho_set, do_rho0=(.NOT. gapw_xc))
204 :
205 : ! *** calculate the hartree potential on the total density ***
206 1784 : CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
207 :
208 1784 : CALL pw_copy(rho1_g(1), rho1_tot_gspace)
209 2334 : DO ispin = 2, nspins
210 2334 : CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
211 : END DO
212 1784 : IF (gapw) THEN
213 0 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
214 0 : IF (ASSOCIATED(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
215 0 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho1_tot_gspace)
216 : END IF
217 : END IF
218 :
219 1784 : scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
220 1784 : IF (cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES") &
221 : /= 0) THEN
222 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
223 0 : extension=".scfLog")
224 0 : CALL print_densities(rho1, rho1_tot_gspace, output_unit)
225 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
226 0 : "PRINT%TOTAL_DENSITIES")
227 : END IF
228 :
229 1784 : IF (.NOT. (nspins == 1 .AND. do_triplet)) THEN
230 : BLOCK
231 : TYPE(pw_c1d_gs_type) :: v_hartree_gspace
232 1784 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
233 : CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
234 : energy_hartree, &
235 1784 : v_hartree_gspace)
236 1784 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
237 1784 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
238 : END BLOCK
239 3568 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
240 : END IF
241 :
242 1784 : CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
243 :
244 : ! *** calculate the xc potential ***
245 1784 : IF (gapw_xc) THEN
246 0 : CALL qs_rho_get(rho1_xc, rho_r=rho1_r, tau_r=tau1_r)
247 : ELSE
248 1784 : CALL qs_rho_get(rho1, rho_r=rho1_r, tau_r=tau1_r)
249 : END IF
250 :
251 1784 : IF (nspins == 1 .AND. do_triplet) THEN
252 :
253 0 : lsd = .TRUE.
254 0 : ALLOCATE (rho1_r_pw(2))
255 0 : DO ispin = 1, 2
256 0 : CALL rho1_r_pw(ispin)%create(rho1_r(1)%pw_grid)
257 0 : CALL pw_transfer(rho1_r(1), rho1_r_pw(ispin))
258 : END DO
259 :
260 0 : IF (ASSOCIATED(tau1_r)) THEN
261 0 : ALLOCATE (tau1_r_pw(2))
262 0 : DO ispin = 1, 2
263 0 : CALL tau1_r_pw(ispin)%create(tau1_r(1)%pw_grid)
264 0 : CALL pw_transfer(tau1_r(1), tau1_r_pw(ispin))
265 : END DO
266 : END IF
267 :
268 : ELSE
269 :
270 1784 : rho1_r_pw => rho1_r
271 :
272 1784 : tau1_r_pw => tau1_r
273 :
274 : END IF
275 :
276 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, p_env%kpp1_env%deriv_set, p_env%kpp1_env%rho_set, &
277 : rho1_r_pw, rho1_g_pw, tau1_r_pw, auxbas_pw_pool, xc_section, .FALSE., &
278 : do_excitations=.TRUE., do_triplet=do_triplet, &
279 1784 : compute_virial=calc_virial, virial_xc=virial)
280 :
281 4118 : DO ispin = 1, nspins
282 4118 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
283 : END DO
284 1784 : v_rspace_new => v_xc
285 1784 : IF (SIZE(v_xc) /= nspins) THEN
286 0 : CALL auxbas_pw_pool%give_back_pw(v_xc(2))
287 : END IF
288 1784 : NULLIFY (v_xc)
289 1784 : IF (ASSOCIATED(v_xc_tau)) THEN
290 616 : DO ispin = 1, nspins
291 616 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
292 : END DO
293 244 : IF (SIZE(v_xc_tau) /= nspins) THEN
294 0 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(2))
295 : END IF
296 : END IF
297 :
298 1784 : IF (gapw .OR. gapw_xc) THEN
299 0 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
300 0 : rho1_atom_set => p_env%local_rho_set%rho_atom_set
301 : CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
302 0 : do_triplet=do_triplet)
303 : END IF
304 :
305 1784 : IF (nspins == 1 .AND. do_triplet) THEN
306 0 : DO ispin = 1, SIZE(rho1_r_pw)
307 0 : CALL rho1_r_pw(ispin)%release()
308 : END DO
309 0 : DEALLOCATE (rho1_r_pw)
310 0 : IF (ASSOCIATED(tau1_r_pw)) THEN
311 0 : DO ispin = 1, SIZE(tau1_r_pw)
312 0 : CALL tau1_r_pw(ispin)%release()
313 : END DO
314 0 : DEALLOCATE (tau1_r_pw)
315 : END IF
316 : END IF
317 :
318 550 : alpha = 1.0_dp
319 1234 : IF (nspins == 1) alpha = 2.0_dp
320 :
321 : !-------------------------------!
322 : ! Add both hartree and xc terms !
323 : !-------------------------------!
324 4118 : DO ispin = 1, nspins
325 2334 : CALL dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix, 0.0_dp)
326 :
327 : ! XC and Hartree are integrated separatedly
328 : ! XC uses the soft basis set only
329 2334 : IF (gapw_xc) THEN
330 :
331 0 : IF (nspins == 1) THEN
332 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
333 : pmat=rho_ao(ispin), &
334 : hmat=p_env%kpp1_env%v_ao(ispin), &
335 : qs_env=qs_env, &
336 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
337 :
338 0 : IF (ASSOCIATED(v_xc_tau)) THEN
339 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
340 : pmat=rho_ao(ispin), &
341 : hmat=p_env%kpp1_env%v_ao(ispin), &
342 : qs_env=qs_env, &
343 : compute_tau=.TRUE., &
344 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
345 : END IF
346 :
347 : ! add hartree only for SINGLETS
348 0 : IF (.NOT. do_triplet) THEN
349 0 : CALL pw_copy(v_hartree_rspace, v_rspace_new(1))
350 :
351 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
352 : pmat=rho_ao(ispin), &
353 : hmat=p_env%kpp1_env%v_ao(ispin), &
354 : qs_env=qs_env, &
355 0 : calculate_forces=my_calc_forces, gapw=gapw)
356 : END IF
357 : ELSE
358 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
359 : pmat=rho_ao(ispin), &
360 : hmat=p_env%kpp1_env%v_ao(ispin), &
361 : qs_env=qs_env, &
362 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
363 :
364 0 : IF (ASSOCIATED(v_xc_tau)) THEN
365 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
366 : pmat=rho_ao(ispin), &
367 : hmat=p_env%kpp1_env%v_ao(ispin), &
368 : qs_env=qs_env, &
369 : compute_tau=.TRUE., &
370 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
371 : END IF
372 :
373 0 : CALL pw_copy(v_hartree_rspace, v_rspace_new(ispin))
374 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
375 : pmat=rho_ao(ispin), &
376 : hmat=p_env%kpp1_env%v_ao(ispin), &
377 : qs_env=qs_env, &
378 0 : calculate_forces=my_calc_forces, gapw=gapw)
379 : END IF
380 :
381 : ELSE
382 :
383 2334 : IF (nspins == 1) THEN
384 :
385 : ! add hartree only for SINGLETS
386 1234 : IF (.NOT. do_triplet) THEN
387 1234 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(1))
388 : END IF
389 : ELSE
390 1100 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin))
391 : END IF
392 :
393 2334 : IF (lrigpw) THEN
394 0 : IF (ASSOCIATED(v_xc_tau)) CPABORT("Meta-GGA functionals not supported with LRI!")
395 :
396 0 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
397 0 : CALL get_qs_env(qs_env, nkind=nkind)
398 0 : DO ikind = 1, nkind
399 0 : lri_v_int(ikind)%v_int = 0.0_dp
400 : END DO
401 : CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
402 0 : lri_v_int, .FALSE., "LRI_AUX")
403 0 : DO ikind = 1, nkind
404 0 : CALL para_env%sum(lri_v_int(ikind)%v_int)
405 : END DO
406 0 : ALLOCATE (k1mat(1))
407 0 : k1mat(1)%matrix => p_env%kpp1_env%v_ao(ispin)%matrix
408 0 : IF (lri_env%exact_1c_terms) THEN
409 : CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, &
410 0 : rho_ao(ispin)%matrix, qs_env, my_calc_forces, "ORB")
411 : END IF
412 0 : CALL calculate_lri_ks_matrix(lri_env, lri_v_int, k1mat, atomic_kind_set)
413 0 : DEALLOCATE (k1mat)
414 : ELSE
415 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
416 : pmat=rho_ao(ispin), &
417 : hmat=p_env%kpp1_env%v_ao(ispin), &
418 : qs_env=qs_env, &
419 2334 : calculate_forces=my_calc_forces, gapw=gapw)
420 :
421 2334 : IF (ASSOCIATED(v_xc_tau)) THEN
422 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
423 : pmat=rho_ao(ispin), &
424 : hmat=p_env%kpp1_env%v_ao(ispin), &
425 : qs_env=qs_env, &
426 : compute_tau=.TRUE., &
427 372 : calculate_forces=my_calc_forces, gapw=gapw)
428 : END IF
429 : END IF
430 : END IF
431 :
432 4118 : CALL dbcsr_add(p_env%kpp1(ispin)%matrix, p_env%kpp1_env%v_ao(ispin)%matrix, 1.0_dp, alpha)
433 : END DO
434 :
435 1784 : IF (gapw) THEN
436 0 : IF (.NOT. (nspins == 1 .AND. do_triplet)) THEN
437 : CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
438 : p_env%hartree_local%ecoul_1c, &
439 : p_env%local_rho_set, &
440 0 : para_env, tddft=.TRUE., core_2nd=.TRUE.)
441 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
442 : calculate_forces=my_calc_forces, &
443 0 : local_rho_set=p_env%local_rho_set)
444 : END IF
445 : ! *** Add single atom contributions to the KS matrix ***
446 : ! remap pointer
447 0 : ns = SIZE(p_env%kpp1)
448 0 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
449 0 : ns = SIZE(rho_ao)
450 0 : psmat(1:ns, 1:1) => rho_ao(1:ns)
451 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=my_calc_forces, tddft=.TRUE., &
452 0 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
453 1784 : ELSEIF (gapw_xc) THEN
454 0 : ns = SIZE(p_env%kpp1)
455 0 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
456 0 : ns = SIZE(rho_ao)
457 0 : psmat(1:ns, 1:1) => rho_ao(1:ns)
458 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=my_calc_forces, tddft=.TRUE., &
459 0 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
460 : END IF
461 :
462 1784 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
463 4118 : DO ispin = 1, SIZE(v_rspace_new)
464 4118 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
465 : END DO
466 1784 : DEALLOCATE (v_rspace_new)
467 1784 : IF (ASSOCIATED(v_xc_tau)) THEN
468 616 : DO ispin = 1, SIZE(v_xc_tau)
469 616 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
470 : END DO
471 244 : DEALLOCATE (v_xc_tau)
472 : END IF
473 :
474 1784 : CALL timestop(handle)
475 1784 : END SUBROUTINE calc_kpp1
476 :
477 : ! **************************************************************************************************
478 : !> \brief checks that the intenal storage is allocated, and allocs it if needed
479 : !> \param kpp1_env the environment to check
480 : !> \param qs_env the qs environment this kpp1_env lives in
481 : !> \param do_triplet ...
482 : !> \author Fawzi Mohamed
483 : !> \note
484 : !> private routine
485 : ! **************************************************************************************************
486 1784 : SUBROUTINE kpp1_check_i_alloc(kpp1_env, qs_env, do_triplet)
487 :
488 : TYPE(qs_kpp1_env_type) :: kpp1_env
489 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
490 : LOGICAL, INTENT(IN) :: do_triplet
491 :
492 : INTEGER :: ispin, nspins
493 : TYPE(admm_type), POINTER :: admm_env
494 1784 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
495 : TYPE(dft_control_type), POINTER :: dft_control
496 : TYPE(pw_env_type), POINTER :: pw_env
497 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
498 1784 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_rho_r, my_tau_r, rho_r, tau_r
499 : TYPE(qs_rho_type), POINTER :: rho
500 : TYPE(section_vals_type), POINTER :: admm_xc_section, input, xc_section
501 :
502 : ! ------------------------------------------------------------------
503 :
504 1784 : NULLIFY (pw_env, auxbas_pw_pool, matrix_s, rho, rho_r, admm_env, dft_control, my_rho_r, my_tau_r)
505 :
506 : CALL get_qs_env(qs_env, pw_env=pw_env, &
507 : matrix_s=matrix_s, rho=rho, input=input, &
508 1784 : admm_env=admm_env, dft_control=dft_control)
509 :
510 1784 : CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
511 1784 : nspins = SIZE(rho_r)
512 :
513 1784 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
514 :
515 1784 : IF (.NOT. ASSOCIATED(kpp1_env%v_ao)) THEN
516 272 : CALL dbcsr_allocate_matrix_set(kpp1_env%v_ao, nspins)
517 630 : DO ispin = 1, nspins
518 358 : ALLOCATE (kpp1_env%v_ao(ispin)%matrix)
519 : CALL dbcsr_copy(kpp1_env%v_ao(ispin)%matrix, matrix_s(1)%matrix, &
520 630 : name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin)))
521 : END DO
522 : END IF
523 :
524 1784 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set)) THEN
525 :
526 272 : IF (nspins == 1 .AND. do_triplet) THEN
527 0 : ALLOCATE (my_rho_r(2))
528 0 : DO ispin = 1, 2
529 0 : CALL auxbas_pw_pool%create_pw(my_rho_r(ispin))
530 0 : CALL pw_axpy(rho_r(1), my_rho_r(ispin), 0.5_dp, 0.0_dp)
531 : END DO
532 0 : IF (dft_control%use_kinetic_energy_density) THEN
533 0 : ALLOCATE (my_tau_r(2))
534 0 : DO ispin = 1, 2
535 0 : CALL auxbas_pw_pool%create_pw(my_tau_r(ispin))
536 0 : CALL pw_axpy(tau_r(1), my_tau_r(ispin), 0.5_dp, 0.0_dp)
537 : END DO
538 : END IF
539 : ELSE
540 272 : my_rho_r => rho_r
541 272 : IF (dft_control%use_kinetic_energy_density) THEN
542 40 : my_tau_r => tau_r
543 : END IF
544 : END IF
545 :
546 272 : IF (dft_control%do_admm) THEN
547 40 : xc_section => admm_env%xc_section_primary
548 : ELSE
549 232 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
550 : END IF
551 :
552 6256 : ALLOCATE (kpp1_env%deriv_set, kpp1_env%rho_set)
553 : CALL xc_prep_2nd_deriv(kpp1_env%deriv_set, kpp1_env%rho_set, &
554 : my_rho_r, auxbas_pw_pool, &
555 272 : xc_section=xc_section, tau_r=my_tau_r)
556 :
557 272 : IF (nspins == 1 .AND. do_triplet) THEN
558 0 : DO ispin = 1, SIZE(my_rho_r)
559 0 : CALL my_rho_r(ispin)%release()
560 : END DO
561 0 : DEALLOCATE (my_rho_r)
562 0 : IF (ASSOCIATED(my_tau_r)) THEN
563 0 : DO ispin = 1, SIZE(my_tau_r)
564 0 : CALL my_tau_r(ispin)%release()
565 : END DO
566 0 : DEALLOCATE (my_tau_r)
567 : END IF
568 : END IF
569 : END IF
570 :
571 : ! ADMM Correction
572 1784 : IF (dft_control%do_admm) THEN
573 212 : IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
574 92 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
575 24 : CPASSERT(.NOT. do_triplet)
576 24 : admm_xc_section => admm_env%xc_section_aux
577 24 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho)
578 24 : CALL qs_rho_get(rho, rho_r=rho_r)
579 552 : ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
580 : CALL xc_prep_2nd_deriv(kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, &
581 : rho_r, auxbas_pw_pool, &
582 24 : xc_section=admm_xc_section)
583 : END IF
584 : END IF
585 : END IF
586 :
587 1784 : END SUBROUTINE kpp1_check_i_alloc
588 :
589 : ! **************************************************************************************************
590 : !> \brief function to advise of changes either in the grids
591 : !> \param kpp1_env the kpp1_env
592 : !> \par History
593 : !> 11.2002 created [fawzi]
594 : !> \author Fawzi Mohamed
595 : ! **************************************************************************************************
596 1688 : SUBROUTINE kpp1_did_change(kpp1_env)
597 : TYPE(qs_kpp1_env_type) :: kpp1_env
598 :
599 1688 : IF (ASSOCIATED(kpp1_env%deriv_set)) THEN
600 0 : CALL xc_dset_release(kpp1_env%deriv_set)
601 0 : DEALLOCATE (kpp1_env%deriv_set)
602 : NULLIFY (kpp1_env%deriv_set)
603 : END IF
604 1688 : IF (ASSOCIATED(kpp1_env%rho_set)) THEN
605 0 : CALL xc_rho_set_release(kpp1_env%rho_set)
606 0 : DEALLOCATE (kpp1_env%rho_set)
607 : END IF
608 :
609 1688 : END SUBROUTINE kpp1_did_change
610 :
611 : ! **************************************************************************************************
612 : !> \brief ...
613 : !> \param rho1 ...
614 : !> \param rho1_tot_gspace ...
615 : !> \param out_unit ...
616 : ! **************************************************************************************************
617 0 : SUBROUTINE print_densities(rho1, rho1_tot_gspace, out_unit)
618 :
619 : TYPE(qs_rho_type), POINTER :: rho1
620 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho1_tot_gspace
621 : INTEGER :: out_unit
622 :
623 : REAL(KIND=dp) :: total_rho_gspace
624 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho1_r
625 :
626 0 : NULLIFY (tot_rho1_r)
627 :
628 0 : total_rho_gspace = pw_integrate_function(rho1_tot_gspace, isign=-1)
629 0 : IF (out_unit > 0) THEN
630 0 : CALL qs_rho_get(rho1, tot_rho_r=tot_rho1_r)
631 : WRITE (UNIT=out_unit, FMT="(T3,A,T60,F20.10)") &
632 0 : "KPP1 total charge density (r-space):", &
633 0 : accurate_sum(tot_rho1_r), &
634 0 : "KPP1 total charge density (g-space):", &
635 0 : total_rho_gspace
636 : END IF
637 :
638 0 : END SUBROUTINE print_densities
639 :
640 : END MODULE qs_kpp1_env_methods
|
