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