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