Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_scale,&
24 : dbcsr_set
25 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
26 : dbcsr_deallocate_matrix_set
27 : USE cp_log_handling, ONLY: cp_get_default_logger,&
28 : cp_logger_type,&
29 : cp_to_string
30 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
31 : cp_print_key_should_output,&
32 : cp_print_key_unit_nr
33 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
34 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
35 : do_method_gapw,&
36 : do_method_gapw_xc,&
37 : tddfpt_excitations,&
38 : tddfpt_triplet
39 : USE input_section_types, ONLY: section_get_ival,&
40 : section_vals_get,&
41 : section_vals_get_subs_vals,&
42 : section_vals_type,&
43 : section_vals_val_get
44 : USE kahan_sum, ONLY: accurate_sum
45 : USE kinds, ONLY: dp
46 : USE lri_environment_types, ONLY: lri_density_type,&
47 : lri_environment_type,&
48 : lri_kind_type
49 : USE lri_ks_methods, ONLY: calculate_lri_ks_matrix
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE pw_env_types, ONLY: pw_env_get,&
52 : pw_env_type
53 : USE pw_methods, ONLY: pw_axpy,&
54 : pw_copy,&
55 : pw_integrate_function,&
56 : pw_scale,&
57 : pw_transfer
58 : USE pw_poisson_methods, ONLY: pw_poisson_solve
59 : USE pw_poisson_types, ONLY: pw_poisson_type
60 : USE pw_pool_types, ONLY: pw_pool_type
61 : USE pw_types, ONLY: pw_c1d_gs_type,&
62 : pw_r3d_rs_type
63 : USE qs_energy_types, ONLY: allocate_qs_energy,&
64 : deallocate_qs_energy,&
65 : qs_energy_type
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_environment_type
68 : USE qs_gapw_densities, ONLY: prepare_gapw_den
69 : USE qs_integrate_potential, ONLY: integrate_v_rspace,&
70 : integrate_v_rspace_diagonal,&
71 : integrate_v_rspace_one_center
72 : USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type
73 : USE qs_ks_atom, ONLY: update_ks_atom
74 : USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
75 : USE qs_p_env_types, ONLY: qs_p_env_type
76 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
77 : USE qs_rho_atom_types, ONLY: rho_atom_type
78 : USE qs_rho_types, ONLY: qs_rho_get,&
79 : qs_rho_type
80 : USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom
81 : USE xc, ONLY: xc_calc_2nd_deriv,&
82 : xc_prep_2nd_deriv
83 : USE xc_derivative_set_types, ONLY: xc_dset_release
84 : USE xc_rho_set_types, ONLY: xc_rho_set_release
85 : #include "./base/base_uses.f90"
86 :
87 : IMPLICIT NONE
88 :
89 : PRIVATE
90 :
91 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
92 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kpp1_env_methods'
93 :
94 : PUBLIC :: kpp1_create, &
95 : kpp1_calc_k_p_p1, &
96 : kpp1_calc_k_p_p1_fdiff, &
97 : kpp1_did_change, &
98 : kpp1_check_i_alloc, &
99 : calc_kpp1
100 :
101 : CONTAINS
102 :
103 : ! **************************************************************************************************
104 : !> \brief allocates and initializes a kpp1_env
105 : !> \param kpp1_env the environment to initialize
106 : !> \par History
107 : !> 07.2002 created [fawzi]
108 : !> \author Fawzi Mohamed
109 : ! **************************************************************************************************
110 1626 : SUBROUTINE kpp1_create(kpp1_env)
111 : TYPE(qs_kpp1_env_type) :: kpp1_env
112 :
113 1626 : NULLIFY (kpp1_env%v_ao, kpp1_env%rho_set, kpp1_env%deriv_set, &
114 1626 : kpp1_env%rho_set_admm, kpp1_env%deriv_set_admm)
115 1626 : END SUBROUTINE kpp1_create
116 :
117 : ! **************************************************************************************************
118 : !> \brief calculates the k_p_p1 kernel of the perturbation theory
119 : !> \param p_env perturbation environment containing kpp1 kernel and kpp1_env
120 : !> \param qs_env kpp1's qs_env
121 : !> \param rho1 the density that represent the first direction along which
122 : !> you should evaluate the derivatives
123 : !> \param rho1_xc ...
124 : ! **************************************************************************************************
125 1920 : SUBROUTINE kpp1_calc_k_p_p1(p_env, qs_env, rho1, rho1_xc)
126 :
127 : TYPE(qs_p_env_type) :: p_env
128 : TYPE(qs_environment_type), POINTER :: qs_env
129 : TYPE(qs_rho_type), POINTER :: rho1
130 : TYPE(qs_rho_type), OPTIONAL, POINTER :: rho1_xc
131 :
132 : CHARACTER(len=*), PARAMETER :: routineN = 'kpp1_calc_k_p_p1'
133 :
134 : INTEGER :: excitations, handle, res_etype
135 : LOGICAL :: do_excitations, do_triplet, explicit, &
136 : lsd_singlets
137 : TYPE(section_vals_type), POINTER :: input, xc_section
138 :
139 480 : CALL timeset(routineN, handle)
140 :
141 480 : NULLIFY (input)
142 :
143 480 : CPASSERT(ASSOCIATED(rho1))
144 :
145 : CALL get_qs_env(qs_env=qs_env, &
146 480 : input=input)
147 :
148 : CALL section_vals_val_get(input, "DFT%EXCITATIONS", &
149 480 : i_val=excitations)
150 : CALL section_vals_val_get(input, "DFT%TDDFPT%LSD_SINGLETS", &
151 480 : l_val=lsd_singlets)
152 : CALL section_vals_val_get(input, "DFT%TDDFPT%RES_ETYPE", &
153 480 : i_val=res_etype)
154 :
155 480 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
156 480 : IF (excitations == tddfpt_excitations) THEN
157 480 : xc_section => section_vals_get_subs_vals(input, "DFT%TDDFPT%XC")
158 480 : CALL section_vals_get(xc_section, explicit=explicit)
159 480 : IF (.NOT. explicit) THEN
160 0 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
161 : END IF
162 : END IF
163 :
164 480 : do_excitations = (excitations == tddfpt_excitations)
165 480 : do_triplet = (res_etype == tddfpt_triplet)
166 :
167 : CALL calc_kpp1(rho1_xc, rho1, xc_section, .TRUE., &
168 : lsd_singlets, .FALSE., do_excitations, &
169 480 : do_triplet, qs_env, p_env)
170 :
171 480 : CALL timestop(handle)
172 480 : END SUBROUTINE kpp1_calc_k_p_p1
173 :
174 : ! **************************************************************************************************
175 : !> \brief ...
176 : !> \param rho1_xc ...
177 : !> \param rho1 ...
178 : !> \param xc_section ...
179 : !> \param do_tddft ...
180 : !> \param lsd_singlets ...
181 : !> \param lrigpw ...
182 : !> \param do_excitations ...
183 : !> \param do_triplet ...
184 : !> \param qs_env ...
185 : !> \param p_env ...
186 : !> \param calc_forces ...
187 : !> \param calc_virial ...
188 : !> \param virial ...
189 : ! **************************************************************************************************
190 2192 : SUBROUTINE calc_kpp1(rho1_xc, rho1, xc_section, do_tddft, lsd_singlets, lrigpw, &
191 : do_excitations, do_triplet, qs_env, p_env, calc_forces, &
192 : calc_virial, virial)
193 :
194 : TYPE(qs_rho_type), POINTER :: rho1_xc, rho1
195 : TYPE(section_vals_type), POINTER :: xc_section
196 : LOGICAL, INTENT(IN) :: do_tddft, lsd_singlets, lrigpw, &
197 : do_excitations, do_triplet
198 : TYPE(qs_environment_type), POINTER :: qs_env
199 : TYPE(qs_p_env_type) :: p_env
200 : LOGICAL, INTENT(IN), OPTIONAL :: calc_forces, calc_virial
201 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
202 : OPTIONAL :: virial
203 :
204 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_kpp1'
205 :
206 : INTEGER :: handle, ikind, ispin, nkind, ns, nspins, &
207 : output_unit
208 : LOGICAL :: gapw, gapw_xc, lsd, my_calc_forces
209 : REAL(KIND=dp) :: alpha, energy_hartree, energy_hartree_1c
210 2192 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
211 : TYPE(cp_logger_type), POINTER :: logger
212 2192 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: k1mat, rho_ao
213 2192 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
214 : TYPE(lri_density_type), POINTER :: lri_density
215 : TYPE(lri_environment_type), POINTER :: lri_env
216 2192 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
217 : TYPE(mp_para_env_type), POINTER :: para_env
218 : TYPE(pw_c1d_gs_type) :: rho1_tot_gspace
219 2192 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_pw
220 : TYPE(pw_env_type), POINTER :: pw_env
221 : TYPE(pw_poisson_type), POINTER :: poisson_env
222 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
223 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
224 2192 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho1_r_pw, tau1_r, tau1_r_pw, &
225 2192 : v_rspace_new, v_xc, v_xc_tau
226 : TYPE(qs_rho_type), POINTER :: rho
227 2192 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
228 : TYPE(section_vals_type), POINTER :: input, scf_section
229 :
230 2192 : CALL timeset(routineN, handle)
231 :
232 2192 : NULLIFY (v_xc, rho1_g, pw_env, rho1_g_pw, tau1_r_pw)
233 2192 : logger => cp_get_default_logger()
234 :
235 2192 : CPASSERT(ASSOCIATED(p_env%kpp1))
236 2192 : CPASSERT(ASSOCIATED(p_env%kpp1_env))
237 2192 : CPASSERT(ASSOCIATED(rho1))
238 :
239 2192 : nspins = SIZE(p_env%kpp1)
240 2192 : lsd = (nspins == 2)
241 :
242 2192 : my_calc_forces = .FALSE.
243 2192 : IF (PRESENT(calc_forces)) my_calc_forces = calc_forces
244 :
245 : CALL get_qs_env(qs_env, &
246 : pw_env=pw_env, &
247 : input=input, &
248 : para_env=para_env, &
249 2192 : rho=rho)
250 :
251 2192 : CPASSERT(ASSOCIATED(rho1))
252 :
253 2192 : IF (lrigpw) THEN
254 : CALL get_qs_env(qs_env, &
255 : lri_env=lri_env, &
256 : lri_density=lri_density, &
257 0 : atomic_kind_set=atomic_kind_set)
258 : END IF
259 :
260 2192 : gapw = (section_get_ival(input, "DFT%QS%METHOD") == do_method_gapw)
261 2192 : gapw_xc = (section_get_ival(input, "DFT%QS%METHOD") == do_method_gapw_xc)
262 2192 : IF (gapw_xc) THEN
263 0 : CPASSERT(ASSOCIATED(rho1_xc))
264 : END IF
265 :
266 2192 : CALL kpp1_check_i_alloc(p_env%kpp1_env, qs_env, do_excitations, lsd_singlets, do_triplet)
267 :
268 2192 : CALL qs_rho_get(rho, rho_ao=rho_ao)
269 2192 : CALL qs_rho_get(rho1, rho_g=rho1_g)
270 :
271 : ! gets the tmp grids
272 2192 : CPASSERT(ASSOCIATED(pw_env))
273 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
274 2192 : poisson_env=poisson_env)
275 2192 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
276 :
277 2192 : IF (gapw .OR. gapw_xc) &
278 0 : CALL prepare_gapw_den(qs_env, p_env%local_rho_set, do_rho0=(.NOT. gapw_xc))
279 :
280 : ! *** calculate the hartree potential on the total density ***
281 2192 : CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
282 :
283 2192 : CALL pw_copy(rho1_g(1), rho1_tot_gspace)
284 2802 : DO ispin = 2, nspins
285 2802 : CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
286 : END DO
287 2192 : IF (gapw) &
288 0 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
289 :
290 2192 : scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
291 2192 : IF (cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES") &
292 : /= 0) THEN
293 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
294 0 : extension=".scfLog")
295 0 : CALL print_densities(rho1, rho1_tot_gspace, output_unit)
296 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
297 0 : "PRINT%TOTAL_DENSITIES")
298 : END IF
299 :
300 2192 : IF (.NOT. (nspins == 1 .AND. do_excitations .AND. do_triplet)) THEN
301 : BLOCK
302 : TYPE(pw_c1d_gs_type) :: v_hartree_gspace
303 2096 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
304 : CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
305 : energy_hartree, &
306 2096 : v_hartree_gspace)
307 2096 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
308 2096 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
309 : END BLOCK
310 4192 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
311 : END IF
312 :
313 2192 : CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
314 :
315 : ! *** calculate the xc potential ***
316 2192 : IF (gapw_xc) THEN
317 0 : CALL qs_rho_get(rho1_xc, rho_r=rho1_r, tau_r=tau1_r)
318 : ELSE
319 2192 : CALL qs_rho_get(rho1, rho_r=rho1_r, tau_r=tau1_r)
320 : END IF
321 :
322 2192 : IF (nspins == 1 .AND. do_excitations .AND. &
323 : (lsd_singlets .OR. do_triplet)) THEN
324 :
325 96 : lsd = .TRUE.
326 288 : ALLOCATE (rho1_r_pw(2))
327 288 : DO ispin = 1, 2
328 192 : CALL rho1_r_pw(ispin)%create(rho1_r(1)%pw_grid)
329 288 : CALL pw_transfer(rho1_r(1), rho1_r_pw(ispin))
330 : END DO
331 :
332 96 : IF (ASSOCIATED(tau1_r)) THEN
333 0 : ALLOCATE (tau1_r_pw(2))
334 0 : DO ispin = 1, 2
335 0 : CALL tau1_r_pw(ispin)%create(tau1_r(1)%pw_grid)
336 0 : CALL pw_transfer(tau1_r(1), tau1_r_pw(ispin))
337 : END DO
338 : END IF
339 :
340 : ELSE
341 :
342 2096 : rho1_r_pw => rho1_r
343 :
344 2096 : tau1_r_pw => tau1_r
345 :
346 : END IF
347 :
348 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, p_env%kpp1_env%deriv_set, p_env%kpp1_env%rho_set, &
349 : rho1_r_pw, rho1_g_pw, tau1_r_pw, auxbas_pw_pool, xc_section, .FALSE., &
350 : lsd_singlets=lsd_singlets, do_excitations=do_excitations, &
351 : do_triplet=do_triplet, do_tddft=do_tddft, &
352 2192 : compute_virial=calc_virial, virial_xc=virial)
353 :
354 4994 : DO ispin = 1, nspins
355 4994 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
356 : END DO
357 2192 : v_rspace_new => v_xc
358 2192 : IF (SIZE(v_xc) /= nspins) THEN
359 96 : CALL auxbas_pw_pool%give_back_pw(v_xc(2))
360 : END IF
361 2192 : NULLIFY (v_xc)
362 2192 : IF (ASSOCIATED(v_xc_tau)) THEN
363 560 : DO ispin = 1, nspins
364 560 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
365 : END DO
366 222 : IF (SIZE(v_xc_tau) /= nspins) THEN
367 0 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(2))
368 : END IF
369 : END IF
370 :
371 2192 : IF (gapw .OR. gapw_xc) THEN
372 0 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
373 0 : rho1_atom_set => p_env%local_rho_set%rho_atom_set
374 : CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
375 0 : do_tddft=do_tddft, do_triplet=do_triplet)
376 : END IF
377 :
378 2192 : IF (nspins == 1 .AND. do_excitations .AND. &
379 : (lsd_singlets .OR. do_triplet)) THEN
380 288 : DO ispin = 1, SIZE(rho1_r_pw)
381 288 : CALL rho1_r_pw(ispin)%release()
382 : END DO
383 96 : DEALLOCATE (rho1_r_pw)
384 96 : IF (ASSOCIATED(tau1_r_pw)) THEN
385 0 : DO ispin = 1, SIZE(tau1_r_pw)
386 0 : CALL tau1_r_pw(ispin)%release()
387 : END DO
388 0 : DEALLOCATE (tau1_r_pw)
389 : END IF
390 : END IF
391 :
392 2192 : alpha = 1.0_dp
393 2192 : IF (do_excitations .AND. nspins == 1) alpha = 2.0_dp
394 :
395 : !-------------------------------!
396 : ! Add both hartree and xc terms !
397 : !-------------------------------!
398 4994 : DO ispin = 1, nspins
399 2802 : CALL dbcsr_set(p_env%kpp1_env%v_ao(ispin)%matrix, 0.0_dp)
400 :
401 : ! XC and Hartree are integrated separatedly
402 : ! XC uses the soft basis set only
403 2802 : IF (gapw_xc) THEN
404 :
405 0 : IF (do_excitations .AND. nspins == 1) THEN
406 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
407 : pmat=rho_ao(ispin), &
408 : hmat=p_env%kpp1_env%v_ao(ispin), &
409 : qs_env=qs_env, &
410 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
411 :
412 0 : IF (ASSOCIATED(v_xc_tau)) THEN
413 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
414 : pmat=rho_ao(ispin), &
415 : hmat=p_env%kpp1_env%v_ao(ispin), &
416 : qs_env=qs_env, &
417 : compute_tau=.TRUE., &
418 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
419 : END IF
420 :
421 : ! add hartree only for SINGLETS
422 0 : IF (.NOT. do_triplet) THEN
423 0 : CALL pw_copy(v_hartree_rspace, v_rspace_new(1))
424 :
425 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
426 : pmat=rho_ao(ispin), &
427 : hmat=p_env%kpp1_env%v_ao(ispin), &
428 : qs_env=qs_env, &
429 0 : calculate_forces=my_calc_forces, gapw=gapw)
430 : END IF
431 : ELSE
432 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
433 : pmat=rho_ao(ispin), &
434 : hmat=p_env%kpp1_env%v_ao(ispin), &
435 : qs_env=qs_env, &
436 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
437 :
438 0 : IF (ASSOCIATED(v_xc_tau)) THEN
439 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
440 : pmat=rho_ao(ispin), &
441 : hmat=p_env%kpp1_env%v_ao(ispin), &
442 : qs_env=qs_env, &
443 : compute_tau=.TRUE., &
444 0 : calculate_forces=my_calc_forces, gapw=gapw_xc)
445 : END IF
446 :
447 0 : CALL pw_copy(v_hartree_rspace, v_rspace_new(ispin))
448 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
449 : pmat=rho_ao(ispin), &
450 : hmat=p_env%kpp1_env%v_ao(ispin), &
451 : qs_env=qs_env, &
452 0 : calculate_forces=my_calc_forces, gapw=gapw)
453 : END IF
454 :
455 : ELSE
456 :
457 2802 : IF (do_excitations .AND. nspins == 1) THEN
458 :
459 : ! add hartree only for SINGLETS
460 1582 : IF (.NOT. do_triplet) THEN
461 1486 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(1))
462 : END IF
463 : ELSE
464 1220 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin))
465 : END IF
466 :
467 2802 : IF (lrigpw) THEN
468 0 : IF (ASSOCIATED(v_xc_tau)) CPABORT("Meta-GGA functionals not supported with LRI!")
469 :
470 0 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
471 0 : CALL get_qs_env(qs_env, nkind=nkind)
472 0 : DO ikind = 1, nkind
473 0 : lri_v_int(ikind)%v_int = 0.0_dp
474 : END DO
475 : CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
476 0 : lri_v_int, .FALSE., "LRI_AUX")
477 0 : DO ikind = 1, nkind
478 0 : CALL para_env%sum(lri_v_int(ikind)%v_int)
479 : END DO
480 0 : ALLOCATE (k1mat(1))
481 0 : k1mat(1)%matrix => p_env%kpp1_env%v_ao(ispin)%matrix
482 0 : IF (lri_env%exact_1c_terms) THEN
483 : CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, &
484 0 : rho_ao(ispin)%matrix, qs_env, my_calc_forces, "ORB")
485 : END IF
486 0 : CALL calculate_lri_ks_matrix(lri_env, lri_v_int, k1mat, atomic_kind_set)
487 0 : DEALLOCATE (k1mat)
488 : ELSE
489 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
490 : pmat=rho_ao(ispin), &
491 : hmat=p_env%kpp1_env%v_ao(ispin), &
492 : qs_env=qs_env, &
493 2802 : calculate_forces=my_calc_forces, gapw=gapw)
494 :
495 2802 : IF (ASSOCIATED(v_xc_tau)) THEN
496 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
497 : pmat=rho_ao(ispin), &
498 : hmat=p_env%kpp1_env%v_ao(ispin), &
499 : qs_env=qs_env, &
500 : compute_tau=.TRUE., &
501 338 : calculate_forces=my_calc_forces, gapw=gapw)
502 : END IF
503 : END IF
504 : END IF
505 :
506 4994 : CALL dbcsr_add(p_env%kpp1(ispin)%matrix, p_env%kpp1_env%v_ao(ispin)%matrix, 1.0_dp, alpha)
507 : END DO
508 :
509 2192 : IF (gapw) THEN
510 0 : IF (.NOT. (do_excitations .AND. nspins == 1 .AND. do_triplet)) THEN
511 : CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
512 : p_env%hartree_local%ecoul_1c, &
513 : p_env%local_rho_set, &
514 0 : para_env, tddft=.TRUE., core_2nd=.TRUE.)
515 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
516 : calculate_forces=my_calc_forces, &
517 0 : local_rho_set=p_env%local_rho_set)
518 : END IF
519 : ! *** Add single atom contributions to the KS matrix ***
520 : ! remap pointer
521 0 : ns = SIZE(p_env%kpp1)
522 0 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
523 0 : ns = SIZE(rho_ao)
524 0 : psmat(1:ns, 1:1) => rho_ao(1:ns)
525 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=my_calc_forces, tddft=.TRUE., &
526 0 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
527 2192 : ELSEIF (gapw_xc) THEN
528 0 : ns = SIZE(p_env%kpp1)
529 0 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
530 0 : ns = SIZE(rho_ao)
531 0 : psmat(1:ns, 1:1) => rho_ao(1:ns)
532 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=my_calc_forces, tddft=.TRUE., &
533 0 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
534 : END IF
535 :
536 2192 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
537 5090 : DO ispin = 1, SIZE(v_rspace_new)
538 5090 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
539 : END DO
540 2192 : DEALLOCATE (v_rspace_new)
541 2192 : IF (ASSOCIATED(v_xc_tau)) THEN
542 560 : DO ispin = 1, SIZE(v_xc_tau)
543 560 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
544 : END DO
545 222 : DEALLOCATE (v_xc_tau)
546 : END IF
547 :
548 2192 : CALL timestop(handle)
549 2192 : END SUBROUTINE calc_kpp1
550 :
551 : ! **************************************************************************************************
552 : !> \brief calcualtes the k_p_p1 kernel of the perturbation theory with finite
553 : !> differences
554 : !> \param qs_env kpp1's qs_env
555 : !> \param k_p_p1 the sparse matrix that will contain the kernel k_p_p1
556 : !> \param rho the density where to evaluate the derivatives (i.e. p along
557 : !> with with its grid representations, that must be valid)
558 : !> \param rho1 the density that represent the first direction along which
559 : !> you should evaluate the derivatives
560 : !> \param diff the amount of the finite difference step
561 : !> \par History
562 : !> 01.2003 created [fawzi]
563 : !> \author Fawzi Mohamed
564 : !> \note
565 : !> useful for testing purposes.
566 : !> rescale my_diff depending on the norm of rho1?
567 : ! **************************************************************************************************
568 0 : SUBROUTINE kpp1_calc_k_p_p1_fdiff(qs_env, k_p_p1, rho, rho1, &
569 : diff)
570 : TYPE(qs_environment_type), POINTER :: qs_env
571 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: k_p_p1
572 : TYPE(qs_rho_type), POINTER :: rho, rho1
573 : REAL(KIND=dp), INTENT(in), OPTIONAL :: diff
574 :
575 : INTEGER :: ispin, nspins
576 : REAL(KIND=dp) :: my_diff
577 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_2, matrix_s, rho1_ao, rho_ao
578 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g, rho_g
579 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho_r
580 : TYPE(qs_energy_type), POINTER :: qs_energy
581 :
582 0 : NULLIFY (ks_2, matrix_s, qs_energy, rho_ao, rho1_ao, rho_r, rho1_r, rho_g, rho1_g)
583 0 : nspins = SIZE(k_p_p1)
584 0 : my_diff = 1.0e-6_dp
585 0 : IF (PRESENT(diff)) my_diff = diff
586 0 : CALL allocate_qs_energy(qs_energy)
587 :
588 0 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r, rho_g=rho_g)
589 0 : CALL qs_rho_get(rho1, rho_ao=rho1_ao, rho_r=rho1_r, rho_g=rho1_g)
590 0 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
591 :
592 : ! rho = rho0+h/2*rho1
593 0 : my_diff = my_diff/2.0_dp
594 0 : DO ispin = 1, SIZE(k_p_p1)
595 : CALL dbcsr_add(rho_ao(ispin)%matrix, rho1_ao(ispin)%matrix, &
596 0 : alpha_scalar=1.0_dp, beta_scalar=my_diff)
597 0 : CALL pw_axpy(rho1_r(ispin), rho_r(ispin), my_diff)
598 0 : CALL pw_axpy(rho1_g(ispin), rho_g(ispin), my_diff)
599 : END DO
600 :
601 : CALL qs_ks_build_kohn_sham_matrix(qs_env, &
602 : ext_ks_matrix=k_p_p1, &
603 : calculate_forces=.FALSE., &
604 0 : just_energy=.FALSE.)
605 :
606 0 : CALL dbcsr_allocate_matrix_set(ks_2, nspins)
607 0 : DO ispin = 1, nspins
608 0 : ALLOCATE (ks_2(ispin)%matrix)
609 : CALL dbcsr_copy(ks_2(ispin)%matrix, matrix_s(1)%matrix, &
610 0 : name="tmp_ks2-"//ADJUSTL(cp_to_string(ispin)))
611 : END DO
612 :
613 : ! rho = rho0-h/2*rho1
614 0 : my_diff = -2.0_dp*my_diff
615 0 : DO ispin = 1, nspins
616 : CALL dbcsr_add(rho_ao(ispin)%matrix, rho1_ao(ispin)%matrix, &
617 0 : alpha_scalar=1.0_dp, beta_scalar=my_diff)
618 0 : CALL pw_axpy(rho1_r(ispin), rho_r(ispin), my_diff)
619 0 : CALL pw_axpy(rho1_g(ispin), rho_g(ispin), my_diff)
620 : END DO
621 :
622 : CALL qs_ks_build_kohn_sham_matrix(qs_env, &
623 : ext_ks_matrix=ks_2, &
624 : calculate_forces=.FALSE., &
625 0 : just_energy=.FALSE.)
626 :
627 : ! rho = rho0
628 0 : my_diff = -0.5_dp*my_diff
629 0 : DO ispin = 1, nspins
630 : CALL dbcsr_add(rho_ao(ispin)%matrix, rho1_ao(ispin)%matrix, &
631 0 : alpha_scalar=1.0_dp, beta_scalar=my_diff)
632 0 : CALL pw_axpy(rho1_r(ispin), rho_r(ispin), my_diff)
633 0 : CALL pw_axpy(rho1_g(ispin), rho_g(ispin), my_diff)
634 : END DO
635 :
636 : ! k_p_p1=(H(rho0+h/2 rho1)-H(rho0-h/2 rho1))/h
637 0 : DO ispin = 1, nspins
638 : CALL dbcsr_add(k_p_p1(ispin)%matrix, ks_2(ispin)%matrix, &
639 0 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
640 0 : CALL dbcsr_scale(k_p_p1(ispin)%matrix, alpha_scalar=0.5_dp/my_diff)
641 : END DO
642 :
643 0 : CALL dbcsr_deallocate_matrix_set(ks_2)
644 0 : CALL deallocate_qs_energy(qs_energy)
645 0 : END SUBROUTINE kpp1_calc_k_p_p1_fdiff
646 :
647 : ! **************************************************************************************************
648 : !> \brief checks that the intenal storage is allocated, and allocs it if needed
649 : !> \param kpp1_env the environment to check
650 : !> \param qs_env the qs environment this kpp1_env lives in
651 : !> \param do_excitations ...
652 : !> \param lsd_singlets ...
653 : !> \param do_triplet ...
654 : !> \author Fawzi Mohamed
655 : !> \note
656 : !> private routine
657 : ! **************************************************************************************************
658 2192 : SUBROUTINE kpp1_check_i_alloc(kpp1_env, qs_env, do_excitations, lsd_singlets, do_triplet)
659 :
660 : TYPE(qs_kpp1_env_type) :: kpp1_env
661 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
662 : LOGICAL, INTENT(IN) :: do_excitations, lsd_singlets, do_triplet
663 :
664 : INTEGER :: ispin, nspins
665 : TYPE(admm_type), POINTER :: admm_env
666 2192 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
667 : TYPE(dft_control_type), POINTER :: dft_control
668 : TYPE(pw_env_type), POINTER :: pw_env
669 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
670 2192 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_rho_r, my_tau_r, rho_r, tau_r
671 : TYPE(qs_rho_type), POINTER :: rho
672 : TYPE(section_vals_type), POINTER :: admm_xc_section, input, xc_section
673 :
674 : ! ------------------------------------------------------------------
675 :
676 2192 : NULLIFY (pw_env, auxbas_pw_pool, matrix_s, rho, rho_r, admm_env, dft_control, my_rho_r, my_tau_r)
677 :
678 : CALL get_qs_env(qs_env, pw_env=pw_env, &
679 : matrix_s=matrix_s, rho=rho, input=input, &
680 2192 : admm_env=admm_env, dft_control=dft_control)
681 :
682 2192 : CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
683 2192 : nspins = SIZE(rho_r)
684 :
685 2192 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
686 :
687 2192 : IF (.NOT. ASSOCIATED(kpp1_env%v_ao)) THEN
688 272 : CALL dbcsr_allocate_matrix_set(kpp1_env%v_ao, nspins)
689 632 : DO ispin = 1, nspins
690 360 : ALLOCATE (kpp1_env%v_ao(ispin)%matrix)
691 : CALL dbcsr_copy(kpp1_env%v_ao(ispin)%matrix, matrix_s(1)%matrix, &
692 632 : name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin)))
693 : END DO
694 : END IF
695 :
696 2192 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set)) THEN
697 :
698 272 : IF (nspins == 1 .AND. (do_excitations .AND. &
699 : (lsd_singlets .OR. do_triplet))) THEN
700 6 : ALLOCATE (my_rho_r(2))
701 6 : DO ispin = 1, 2
702 4 : CALL auxbas_pw_pool%create_pw(my_rho_r(ispin))
703 6 : CALL pw_axpy(rho_r(1), my_rho_r(ispin), 0.5_dp, 0.0_dp)
704 : END DO
705 2 : IF (dft_control%use_kinetic_energy_density) THEN
706 0 : ALLOCATE (my_tau_r(2))
707 0 : DO ispin = 1, 2
708 0 : CALL auxbas_pw_pool%create_pw(my_tau_r(ispin))
709 0 : CALL pw_axpy(tau_r(1), my_tau_r(ispin), 0.5_dp, 0.0_dp)
710 : END DO
711 : END IF
712 : ELSE
713 270 : my_rho_r => rho_r
714 270 : IF (dft_control%use_kinetic_energy_density) THEN
715 36 : my_tau_r => tau_r
716 : END IF
717 : END IF
718 :
719 272 : IF (dft_control%do_admm) THEN
720 40 : xc_section => admm_env%xc_section_primary
721 : ELSE
722 232 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
723 : END IF
724 :
725 6256 : ALLOCATE (kpp1_env%deriv_set, kpp1_env%rho_set)
726 : CALL xc_prep_2nd_deriv(kpp1_env%deriv_set, kpp1_env%rho_set, &
727 : my_rho_r, auxbas_pw_pool, &
728 272 : xc_section=xc_section, tau_r=my_tau_r)
729 :
730 272 : IF (nspins == 1 .AND. (do_excitations .AND. &
731 : (lsd_singlets .OR. do_triplet))) THEN
732 6 : DO ispin = 1, SIZE(my_rho_r)
733 6 : CALL my_rho_r(ispin)%release()
734 : END DO
735 2 : DEALLOCATE (my_rho_r)
736 2 : IF (ASSOCIATED(my_tau_r)) THEN
737 0 : DO ispin = 1, SIZE(my_tau_r)
738 0 : CALL my_tau_r(ispin)%release()
739 : END DO
740 0 : DEALLOCATE (my_tau_r)
741 : END IF
742 : END IF
743 : END IF
744 :
745 : ! ADMM Correction
746 2192 : IF (dft_control%do_admm) THEN
747 212 : IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
748 92 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
749 24 : CPASSERT(.NOT. do_triplet)
750 24 : admm_xc_section => admm_env%xc_section_aux
751 24 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho)
752 24 : CALL qs_rho_get(rho, rho_r=rho_r)
753 552 : ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
754 : CALL xc_prep_2nd_deriv(kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, &
755 : rho_r, auxbas_pw_pool, &
756 24 : xc_section=admm_xc_section)
757 : END IF
758 : END IF
759 : END IF
760 :
761 2192 : END SUBROUTINE kpp1_check_i_alloc
762 :
763 : ! **************************************************************************************************
764 : !> \brief function to advise of changes either in the grids
765 : !> \param kpp1_env the kpp1_env
766 : !> \par History
767 : !> 11.2002 created [fawzi]
768 : !> \author Fawzi Mohamed
769 : ! **************************************************************************************************
770 1626 : SUBROUTINE kpp1_did_change(kpp1_env)
771 : TYPE(qs_kpp1_env_type) :: kpp1_env
772 :
773 1626 : IF (ASSOCIATED(kpp1_env%deriv_set)) THEN
774 0 : CALL xc_dset_release(kpp1_env%deriv_set)
775 0 : DEALLOCATE (kpp1_env%deriv_set)
776 : NULLIFY (kpp1_env%deriv_set)
777 : END IF
778 1626 : IF (ASSOCIATED(kpp1_env%rho_set)) THEN
779 0 : CALL xc_rho_set_release(kpp1_env%rho_set)
780 0 : DEALLOCATE (kpp1_env%rho_set)
781 : END IF
782 :
783 1626 : END SUBROUTINE kpp1_did_change
784 :
785 : ! **************************************************************************************************
786 : !> \brief ...
787 : !> \param rho1 ...
788 : !> \param rho1_tot_gspace ...
789 : !> \param out_unit ...
790 : ! **************************************************************************************************
791 0 : SUBROUTINE print_densities(rho1, rho1_tot_gspace, out_unit)
792 :
793 : TYPE(qs_rho_type), POINTER :: rho1
794 : TYPE(pw_c1d_gs_type), INTENT(IN) :: rho1_tot_gspace
795 : INTEGER :: out_unit
796 :
797 : REAL(KIND=dp) :: total_rho_gspace
798 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho1_r
799 :
800 0 : NULLIFY (tot_rho1_r)
801 :
802 0 : total_rho_gspace = pw_integrate_function(rho1_tot_gspace, isign=-1)
803 0 : IF (out_unit > 0) THEN
804 0 : CALL qs_rho_get(rho1, tot_rho_r=tot_rho1_r)
805 : WRITE (UNIT=out_unit, FMT="(T3,A,T60,F20.10)") &
806 0 : "KPP1 total charge density (r-space):", &
807 0 : accurate_sum(tot_rho1_r), &
808 0 : "KPP1 total charge density (g-space):", &
809 0 : total_rho_gspace
810 : END IF
811 :
812 0 : END SUBROUTINE print_densities
813 :
814 : END MODULE qs_kpp1_env_methods
|
