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 linres kernel functions
10 : !> \par History
11 : !> created from qs_linres_methods
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_linres_kernel
15 : USE admm_types, ONLY: admm_type,&
16 : get_admm_env
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
21 : dbcsr_copy,&
22 : dbcsr_create,&
23 : dbcsr_deallocate_matrix,&
24 : dbcsr_p_type,&
25 : dbcsr_set
26 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
27 : dbcsr_allocate_matrix_set,&
28 : dbcsr_deallocate_matrix_set
29 : USE cp_fm_types, ONLY: cp_fm_get_info,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_type
33 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
34 : USE hfx_energy_potential, ONLY: integrate_four_center
35 : USE hfx_ri, ONLY: hfx_ri_update_ks
36 : USE hfx_types, ONLY: hfx_type
37 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
38 : do_admm_basis_projection,&
39 : do_admm_exch_scaling_none,&
40 : do_admm_purify_none,&
41 : kg_tnadd_embed
42 : USE input_section_types, ONLY: section_vals_get,&
43 : section_vals_get_subs_vals,&
44 : section_vals_type,&
45 : section_vals_val_get
46 : USE kg_correction, ONLY: kg_ekin_subset
47 : USE kg_environment_types, ONLY: kg_environment_type
48 : USE kinds, ONLY: default_string_length,&
49 : dp
50 : USE lri_environment_types, ONLY: lri_density_type,&
51 : lri_environment_type,&
52 : lri_kind_type
53 : USE lri_ks_methods, ONLY: calculate_lri_ks_matrix
54 : USE message_passing, ONLY: mp_para_env_type
55 : USE mulliken, ONLY: ao_charges
56 : USE particle_types, ONLY: particle_type
57 : USE pw_env_types, ONLY: pw_env_get,&
58 : pw_env_type
59 : USE pw_methods, ONLY: pw_axpy,&
60 : pw_copy,&
61 : pw_scale,&
62 : pw_transfer
63 : USE pw_poisson_methods, ONLY: pw_poisson_solve
64 : USE pw_poisson_types, ONLY: pw_poisson_type
65 : USE pw_pool_types, ONLY: pw_pool_type
66 : USE pw_types, ONLY: pw_c1d_gs_type,&
67 : pw_r3d_rs_type
68 : USE qs_environment_types, ONLY: get_qs_env,&
69 : qs_environment_type
70 : USE qs_fxc, ONLY: qs_fxc_apply,&
71 : qs_fxc_prep
72 : USE qs_gapw_densities, ONLY: prepare_gapw_den
73 : USE qs_integrate_potential, ONLY: integrate_v_rspace,&
74 : integrate_v_rspace_diagonal,&
75 : integrate_v_rspace_one_center
76 : USE qs_kind_types, ONLY: get_qs_kind,&
77 : get_qs_kind_set,&
78 : qs_kind_type
79 : USE qs_kpp1_env_methods, ONLY: kpp1_check_i_alloc
80 : USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type
81 : USE qs_ks_atom, ONLY: update_ks_atom
82 : USE qs_linres_types, ONLY: linres_control_type
83 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
84 : USE qs_p_env_methods, ONLY: p_env_finish_kpp1
85 : USE qs_p_env_types, ONLY: qs_p_env_type
86 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
87 : USE qs_rho_atom_types, ONLY: rho_atom_type
88 : USE qs_rho_methods, ONLY: qs_rho_rebuild,&
89 : qs_rho_update_rho
90 : USE qs_rho_types, ONLY: qs_rho_create,&
91 : qs_rho_get,&
92 : qs_rho_release,&
93 : qs_rho_type
94 : USE task_list_types, ONLY: task_list_type
95 : USE xtb_ehess, ONLY: xtb_coulomb_hessian
96 : USE xtb_types, ONLY: get_xtb_atom_param,&
97 : xtb_atom_type
98 : #include "./base/base_uses.f90"
99 :
100 : IMPLICIT NONE
101 :
102 : PRIVATE
103 :
104 : ! *** Public subroutines ***
105 : PUBLIC :: apply_xc_admm
106 : PUBLIC :: apply_hfx
107 : PUBLIC :: apply_hxc_kernel_kp
108 : PUBLIC :: apply_op_2
109 : PUBLIC :: hfx_matrix
110 :
111 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_kernel'
112 :
113 : ! **************************************************************************************************
114 :
115 : CONTAINS
116 :
117 : ! **************************************************************************************************
118 : !> \brief Apply the periodic GPW Hartree-XC kernel to a K-point AO density response.
119 : !> \param qs_env Ground-state QS environment defining the density and XC kernel.
120 : !> \param kpp1_env Persistent XC derivative data at the ground-state density.
121 : !> \param rho1_ao_kp Real-space image matrices of the spin-resolved density response.
122 : !> \param v1_ao_kp Real-space image matrices of the resulting Hartree-XC potential response.
123 : !>
124 : !> The input is the physical density response, including occupations and irreducible K-point
125 : !> weights. Consequently, unlike the occupied-virtual linear-response path, no closed-shell
126 : !> factor of two is introduced here.
127 : ! **************************************************************************************************
128 198 : SUBROUTINE apply_hxc_kernel_kp(qs_env, kpp1_env, rho1_ao_kp, v1_ao_kp)
129 : TYPE(qs_environment_type), POINTER :: qs_env
130 : TYPE(qs_kpp1_env_type), INTENT(INOUT) :: kpp1_env
131 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: rho1_ao_kp
132 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: v1_ao_kp
133 :
134 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_hxc_kernel_kp'
135 :
136 : INTEGER :: handle, img, ispin, nimages, nspins
137 : LOGICAL :: do_hfx
138 : REAL(KIND=dp) :: energy_hartree
139 : TYPE(admm_type), POINTER :: admm_env
140 198 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: v1_ao_spin
141 198 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho1_work
142 : TYPE(dft_control_type), POINTER :: dft_control
143 : TYPE(pw_c1d_gs_type) :: rho1_tot_gspace, v_hartree_gspace
144 198 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
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 198 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_hxc, v_xc, v_xc_tau
150 : TYPE(qs_rho_type), POINTER :: rho, rho1
151 : TYPE(qs_rho_type), TARGET :: rho1_store
152 198 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
153 : TYPE(section_vals_type), POINTER :: hfx_section, input, xc_section
154 :
155 198 : CALL timeset(routineN, handle)
156 :
157 198 : NULLIFY (admm_env, auxbas_pw_pool, dft_control, hfx_section, input, poisson_env, pw_env, rho, rho1, &
158 198 : rho1_g, rho1_work, rho_atom_set, rho1_atom_set, v1_ao_spin, v_hxc, v_xc, &
159 198 : v_xc_tau, xc_section)
160 : CALL get_qs_env(qs_env=qs_env, admm_env=admm_env, dft_control=dft_control, input=input, &
161 198 : pw_env=pw_env, rho=rho)
162 :
163 198 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
164 : dft_control%qs_control%xtb) THEN
165 0 : CPABORT("The periodic AO Hartree-XC kernel is only available for DFT")
166 : END IF
167 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc .OR. &
168 198 : dft_control%qs_control%lrigpw .OR. dft_control%qs_control%rigpw) THEN
169 0 : CPABORT("The periodic AO Hartree-XC kernel currently requires GPW")
170 : END IF
171 198 : IF (dft_control%do_admm) THEN
172 0 : CPABORT("The periodic AO Hartree-XC kernel currently excludes ADMM")
173 : END IF
174 198 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
175 198 : CALL section_vals_get(hfx_section, explicit=do_hfx)
176 198 : IF (do_hfx) THEN
177 0 : CPABORT("The periodic AO Hartree-XC kernel currently excludes exact exchange")
178 : END IF
179 :
180 198 : nspins = dft_control%nspins
181 198 : nimages = dft_control%nimages
182 594 : CPASSERT(ALL(SHAPE(rho1_ao_kp) == [nspins, nimages]))
183 :
184 198 : IF (.NOT. ASSOCIATED(v1_ao_kp)) THEN
185 198 : CALL dbcsr_allocate_matrix_set(v1_ao_kp, nspins, nimages)
186 446 : DO ispin = 1, nspins
187 54270 : DO img = 1, nimages
188 53824 : ALLOCATE (v1_ao_kp(ispin, img)%matrix)
189 : CALL dbcsr_copy(v1_ao_kp(ispin, img)%matrix, rho1_ao_kp(ispin, img)%matrix, &
190 54072 : name="K-point Hartree-XC response")
191 : END DO
192 : END DO
193 : ELSE
194 0 : CPASSERT(ALL(SHAPE(v1_ao_kp) == [nspins, nimages]))
195 : END IF
196 446 : DO ispin = 1, nspins
197 54270 : DO img = 1, nimages
198 54072 : CALL dbcsr_set(v1_ao_kp(ispin, img)%matrix, 0.0_dp)
199 : END DO
200 : END DO
201 :
202 198 : CALL qs_rho_create(rho1_store)
203 198 : rho1 => rho1_store
204 198 : CALL qs_rho_rebuild(rho1_store, qs_env, rebuild_ao=.TRUE., rebuild_grids=.TRUE.)
205 198 : CALL qs_rho_get(rho1, rho_ao_kp=rho1_work)
206 446 : DO ispin = 1, nspins
207 54270 : DO img = 1, nimages
208 54072 : CALL dbcsr_copy(rho1_work(ispin, img)%matrix, rho1_ao_kp(ispin, img)%matrix)
209 : END DO
210 : END DO
211 198 : CALL qs_rho_update_rho(rho1, qs_env)
212 :
213 198 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
214 198 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set)) THEN
215 2668 : ALLOCATE (kpp1_env%deriv_set, kpp1_env%rho_set)
216 : CALL qs_fxc_prep(qs_env, rho, kpp1_env%rho_set, kpp1_env%deriv_set, &
217 116 : xc_section, pw_env, is_triplet=.FALSE.)
218 : END IF
219 :
220 198 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
221 198 : CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
222 198 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
223 198 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
224 198 : CALL qs_rho_get(rho1, rho_g=rho1_g)
225 198 : CALL pw_copy(rho1_g(1), rho1_tot_gspace)
226 248 : DO ispin = 2, nspins
227 248 : CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
228 : END DO
229 : energy_hartree = 0.0_dp
230 198 : CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, energy_hartree, v_hartree_gspace)
231 198 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
232 198 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
233 :
234 : CALL qs_fxc_apply(qs_env, kpp1_env%deriv_set, kpp1_env%rho_set, rho1, &
235 198 : rho_atom_set, xc_section, .FALSE., v_xc, v_xc_tau, rho1_atom_set)
236 198 : v_hxc => v_xc
237 198 : NULLIFY (v_xc)
238 446 : DO ispin = 1, nspins
239 248 : CALL pw_scale(v_hxc(ispin), v_hxc(ispin)%pw_grid%dvol)
240 248 : CALL pw_axpy(v_hartree_rspace, v_hxc(ispin))
241 248 : v1_ao_spin => v1_ao_kp(ispin, :)
242 : CALL integrate_v_rspace(v_rspace=v_hxc(ispin), hmat_kp=v1_ao_spin, &
243 446 : qs_env=qs_env, calculate_forces=.FALSE.)
244 : END DO
245 198 : IF (ASSOCIATED(v_xc_tau)) THEN
246 0 : DO ispin = 1, nspins
247 0 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
248 0 : v1_ao_spin => v1_ao_kp(ispin, :)
249 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), hmat_kp=v1_ao_spin, &
250 0 : qs_env=qs_env, compute_tau=.TRUE., calculate_forces=.FALSE.)
251 : END DO
252 : END IF
253 198 : NULLIFY (v1_ao_spin)
254 :
255 198 : CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
256 198 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
257 198 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
258 446 : DO ispin = 1, SIZE(v_hxc)
259 446 : CALL auxbas_pw_pool%give_back_pw(v_hxc(ispin))
260 : END DO
261 198 : DEALLOCATE (v_hxc)
262 198 : IF (ASSOCIATED(v_xc_tau)) THEN
263 0 : DO ispin = 1, SIZE(v_xc_tau)
264 0 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
265 : END DO
266 0 : DEALLOCATE (v_xc_tau)
267 : END IF
268 198 : CALL qs_rho_release(rho1_store)
269 :
270 198 : CALL timestop(handle)
271 198 : END SUBROUTINE apply_hxc_kernel_kp
272 :
273 : ! **************************************************************************************************
274 : !> \brief ...
275 : !> \param qs_env ...
276 : !> \param p_env ...
277 : !> \param c0 ...
278 : !> \param Av ...
279 : ! **************************************************************************************************
280 9582 : SUBROUTINE apply_op_2(qs_env, p_env, c0, Av)
281 : !
282 : TYPE(qs_environment_type), POINTER :: qs_env
283 : TYPE(qs_p_env_type) :: p_env
284 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c0
285 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: Av
286 :
287 : INTEGER :: ispin, ncol
288 : TYPE(dft_control_type), POINTER :: dft_control
289 :
290 9582 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
291 9582 : IF (dft_control%qs_control%semi_empirical) THEN
292 0 : CPABORT("Linear response not available with SE methods")
293 9582 : ELSE IF (dft_control%qs_control%dftb) THEN
294 0 : CPABORT("Linear response not available with DFTB")
295 9582 : ELSE IF (dft_control%qs_control%xtb) THEN
296 206 : CALL apply_op_2_xtb(qs_env, p_env)
297 : ELSE
298 9376 : CALL apply_op_2_dft(qs_env, p_env)
299 9376 : CALL apply_hfx(qs_env, p_env)
300 9376 : CALL apply_xc_admm(qs_env, p_env)
301 9376 : IF (dft_control%do_admm) CALL p_env_finish_kpp1(qs_env, p_env)
302 : END IF
303 :
304 20278 : DO ispin = 1, SIZE(c0)
305 10696 : CALL cp_fm_get_info(c0(ispin), ncol_global=ncol)
306 : CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1(ispin)%matrix, &
307 : c0(ispin), &
308 : Av(ispin), &
309 20278 : ncol=ncol, alpha=1.0_dp, beta=1.0_dp)
310 : END DO
311 :
312 9582 : END SUBROUTINE apply_op_2
313 :
314 : ! **************************************************************************************************
315 : !> \brief ...
316 : !> \param qs_env ...
317 : !> \param p_env ...
318 : ! **************************************************************************************************
319 9376 : SUBROUTINE apply_op_2_dft(qs_env, p_env)
320 : TYPE(qs_environment_type), POINTER :: qs_env
321 : TYPE(qs_p_env_type) :: p_env
322 :
323 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_dft'
324 :
325 : INTEGER :: handle, ikind, ispin, nkind, ns, nspins
326 : LOGICAL :: do_onecenter, gapw, gapw_xc, lr_triplet, &
327 : lrigpw
328 : REAL(KIND=dp) :: alpha, ekin_mol, energy_hartree, &
329 : energy_hartree_1c
330 : TYPE(admm_type), POINTER :: admm_env
331 9376 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
332 : TYPE(cp_logger_type), POINTER :: logger
333 9376 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: k1mat, rho1_ao, rho_ao
334 9376 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
335 : TYPE(dft_control_type), POINTER :: dft_control
336 : TYPE(kg_environment_type), POINTER :: kg_env
337 : TYPE(linres_control_type), POINTER :: linres_control
338 : TYPE(lri_density_type), POINTER :: lri_density
339 : TYPE(lri_environment_type), POINTER :: lri_env
340 9376 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
341 : TYPE(mp_para_env_type), POINTER :: para_env
342 : TYPE(pw_c1d_gs_type) :: rho1_tot_gspace, v_hartree_gspace
343 9376 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
344 : TYPE(pw_env_type), POINTER :: pw_env
345 : TYPE(pw_poisson_type), POINTER :: poisson_env
346 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
347 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
348 9376 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace_new, v_xc, v_xc_tau
349 : TYPE(qs_kpp1_env_type), POINTER :: kpp1_env
350 : TYPE(qs_rho_type), POINTER :: rho, rho1, rho1_xc, rho1a, rho_aux, &
351 : rho_xc
352 9376 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
353 : TYPE(section_vals_type), POINTER :: input, xc_section, xc_section_aux
354 :
355 9376 : CALL timeset(routineN, handle)
356 :
357 9376 : NULLIFY (auxbas_pw_pool, pw_env, v_rspace_new, para_env, v_xc, &
358 9376 : rho1_ao, rho_ao, poisson_env, input, rho, dft_control, logger, &
359 9376 : rho1_g, v_xc_tau)
360 9376 : logger => cp_get_default_logger()
361 :
362 9376 : energy_hartree = 0.0_dp
363 9376 : energy_hartree_1c = 0.0_dp
364 :
365 9376 : CPASSERT(ASSOCIATED(p_env%kpp1))
366 9376 : CPASSERT(ASSOCIATED(p_env%kpp1_env))
367 9376 : kpp1_env => p_env%kpp1_env
368 :
369 : CALL get_qs_env(qs_env=qs_env, &
370 : pw_env=pw_env, &
371 : input=input, &
372 : admm_env=admm_env, &
373 : para_env=para_env, &
374 : rho=rho, &
375 : rho_xc=rho_xc, &
376 : linres_control=linres_control, &
377 9376 : dft_control=dft_control)
378 :
379 9376 : gapw = dft_control%qs_control%gapw
380 9376 : gapw_xc = dft_control%qs_control%gapw_xc
381 9376 : do_onecenter = gapw .OR. gapw_xc
382 9376 : lr_triplet = linres_control%lr_triplet
383 :
384 9376 : rho1 => p_env%rho1
385 9376 : rho1_xc => p_env%rho1_xc
386 9376 : CPASSERT(ASSOCIATED(rho1))
387 9376 : IF (gapw_xc) THEN
388 560 : CPASSERT(ASSOCIATED(rho1_xc))
389 : END IF
390 :
391 9376 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
392 9376 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
393 :
394 9376 : nspins = SIZE(p_env%kpp1)
395 9376 : lrigpw = dft_control%qs_control%lrigpw
396 9376 : IF (lrigpw) THEN
397 : CALL get_qs_env(qs_env, &
398 : lri_env=lri_env, &
399 : lri_density=lri_density, &
400 72 : atomic_kind_set=atomic_kind_set)
401 : END IF
402 :
403 9376 : IF (dft_control%do_admm) THEN
404 2120 : xc_section => admm_env%xc_section_primary
405 : ELSE
406 7256 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
407 : END IF
408 :
409 9376 : CALL kpp1_check_i_alloc(kpp1_env, qs_env, xc_section)
410 :
411 : ! gets the tmp grids
412 9376 : CPASSERT(ASSOCIATED(pw_env))
413 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
414 9376 : poisson_env=poisson_env)
415 9376 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
416 9376 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
417 :
418 9376 : IF (gapw .OR. gapw_xc) THEN
419 2486 : CALL prepare_gapw_den(qs_env, p_env%local_rho_set, do_rho0=(.NOT. gapw_xc))
420 : END IF
421 :
422 : ! *** calculate the hartree potential on the total density ***
423 9376 : CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
424 :
425 9376 : CALL qs_rho_get(rho1, rho_g=rho1_g)
426 9376 : CALL pw_copy(rho1_g(1), rho1_tot_gspace)
427 10394 : DO ispin = 2, nspins
428 10394 : CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
429 : END DO
430 9376 : IF (gapw) THEN
431 1926 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
432 1926 : IF (ASSOCIATED(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
433 0 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho1_tot_gspace)
434 : END IF
435 : END IF
436 :
437 9376 : IF (.NOT. (nspins == 1 .AND. lr_triplet)) THEN
438 : CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
439 : energy_hartree, &
440 9376 : v_hartree_gspace)
441 9376 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
442 : END IF
443 :
444 9376 : CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
445 :
446 : ! *** calculate the xc potential ***
447 : NULLIFY (rho1a)
448 9376 : IF (gapw_xc) THEN
449 560 : rho1a => rho1_xc
450 : ELSE
451 8816 : rho1a => rho1
452 : END IF
453 :
454 9376 : NULLIFY (v_xc_tau)
455 9376 : NULLIFY (rho_atom_set, rho1_atom_set)
456 9376 : IF (do_onecenter) THEN
457 2486 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
458 2486 : rho1_atom_set => p_env%local_rho_set%rho_atom_set
459 : END IF
460 : CALL qs_fxc_apply(qs_env, kpp1_env%deriv_set, kpp1_env%rho_set, rho1a, rho_atom_set, &
461 9376 : xc_section, do_onecenter, v_xc, v_xc_tau, rho1_atom_set)
462 :
463 9376 : v_rspace_new => v_xc
464 9376 : NULLIFY (v_xc)
465 :
466 9376 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
467 19770 : DO ispin = 1, nspins
468 10394 : CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
469 19770 : IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
470 : END DO
471 :
472 : ! ADMM Correction
473 9376 : IF (dft_control%do_admm) THEN
474 2120 : IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
475 1274 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
476 0 : CPASSERT(.NOT. lr_triplet)
477 0 : xc_section_aux => admm_env%xc_section_aux
478 0 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
479 0 : ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
480 : CALL qs_fxc_prep(qs_env, rho_aux, kpp1_env%rho_set_admm, kpp1_env%deriv_set_admm, &
481 0 : xc_section_aux, pw_env, is_triplet=.FALSE.)
482 : END IF
483 : END IF
484 : END IF
485 :
486 : !-------------------------------!
487 : ! Add both hartree and xc terms !
488 : !-------------------------------!
489 19770 : DO ispin = 1, nspins
490 10394 : CALL dbcsr_set(kpp1_env%v_ao(ispin)%matrix, 0.0_dp)
491 :
492 10394 : IF (gapw_xc) THEN
493 : ! XC and Hartree are integrated separatedly
494 : ! XC uses the soft basis set only
495 :
496 588 : IF (nspins == 1) THEN
497 :
498 532 : IF (.NOT. (lr_triplet)) THEN
499 532 : CALL pw_scale(v_rspace_new(1), 2.0_dp)
500 532 : IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1), 2.0_dp)
501 : END IF
502 532 : CALL qs_rho_get(rho1, rho_ao=rho1_ao)
503 : ! remove kpp1_env%v_ao and work directly on k_p_p1 ?
504 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
505 : pmat=rho1_ao(ispin), &
506 : hmat=kpp1_env%v_ao(ispin), &
507 : qs_env=qs_env, &
508 532 : calculate_forces=.FALSE., gapw=gapw_xc)
509 :
510 532 : IF (ASSOCIATED(v_xc_tau)) THEN
511 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
512 : pmat=rho1_ao(ispin), &
513 : hmat=kpp1_env%v_ao(ispin), &
514 : qs_env=qs_env, &
515 : compute_tau=.TRUE., &
516 0 : calculate_forces=.FALSE., gapw=gapw_xc)
517 : END IF
518 :
519 : ! add hartree only for SINGLETS
520 532 : IF (.NOT. lr_triplet) THEN
521 532 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp, 0.0_dp)
522 :
523 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
524 : pmat=rho_ao(ispin), &
525 : hmat=kpp1_env%v_ao(ispin), &
526 : qs_env=qs_env, &
527 532 : calculate_forces=.FALSE., gapw=gapw)
528 : END IF
529 : ELSE
530 : ! remove kpp1_env%v_ao and work directly on k_p_p1 ?
531 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
532 : pmat=rho_ao(ispin), &
533 : hmat=kpp1_env%v_ao(ispin), &
534 : qs_env=qs_env, &
535 56 : calculate_forces=.FALSE., gapw=gapw_xc)
536 :
537 56 : IF (ASSOCIATED(v_xc_tau)) THEN
538 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
539 : pmat=rho_ao(ispin), &
540 : hmat=kpp1_env%v_ao(ispin), &
541 : qs_env=qs_env, &
542 : compute_tau=.TRUE., &
543 0 : calculate_forces=.FALSE., gapw=gapw_xc)
544 : END IF
545 :
546 56 : CALL pw_copy(v_hartree_rspace, v_rspace_new(ispin))
547 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
548 : pmat=rho_ao(ispin), &
549 : hmat=kpp1_env%v_ao(ispin), &
550 : qs_env=qs_env, &
551 56 : calculate_forces=.FALSE., gapw=gapw)
552 : END IF
553 :
554 : ELSE
555 :
556 9806 : IF (nspins == 1) THEN
557 7826 : IF (.NOT. (lr_triplet)) THEN
558 7826 : CALL pw_scale(v_rspace_new(1), 2.0_dp)
559 7826 : IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1), 2.0_dp)
560 : END IF
561 : ! add hartree only for SINGLETS
562 : !IF (res_etype == tddfpt_singlet) THEN
563 : IF (.NOT. lr_triplet) THEN
564 7826 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp)
565 : END IF
566 : ELSE
567 1980 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin), 1.0_dp)
568 : END IF
569 :
570 9806 : IF (lrigpw) THEN
571 72 : IF (ASSOCIATED(v_xc_tau)) THEN
572 0 : CPABORT("metaGGA-functionals not supported with LRI!")
573 : END IF
574 :
575 72 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
576 72 : CALL get_qs_env(qs_env, nkind=nkind)
577 216 : DO ikind = 1, nkind
578 43008 : lri_v_int(ikind)%v_int = 0.0_dp
579 : END DO
580 : CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
581 72 : lri_v_int, .FALSE., "LRI_AUX")
582 216 : DO ikind = 1, nkind
583 85800 : CALL para_env%sum(lri_v_int(ikind)%v_int)
584 : END DO
585 144 : ALLOCATE (k1mat(1))
586 72 : k1mat(1)%matrix => kpp1_env%v_ao(ispin)%matrix
587 72 : IF (lri_env%exact_1c_terms) THEN
588 : CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, &
589 0 : rho_ao(ispin)%matrix, qs_env, .FALSE., "ORB")
590 : END IF
591 72 : CALL calculate_lri_ks_matrix(lri_env, lri_v_int, k1mat, atomic_kind_set)
592 72 : DEALLOCATE (k1mat)
593 : ELSE
594 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
595 : pmat=rho_ao(ispin), &
596 : hmat=kpp1_env%v_ao(ispin), &
597 : qs_env=qs_env, &
598 9734 : calculate_forces=.FALSE., gapw=gapw)
599 :
600 9734 : IF (ASSOCIATED(v_xc_tau)) THEN
601 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
602 : pmat=rho_ao(ispin), &
603 : hmat=kpp1_env%v_ao(ispin), &
604 : qs_env=qs_env, &
605 : compute_tau=.TRUE., &
606 254 : calculate_forces=.FALSE., gapw=gapw)
607 : END IF
608 : END IF
609 :
610 : END IF
611 :
612 19770 : CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, kpp1_env%v_ao(ispin)%matrix)
613 : END DO
614 :
615 9376 : IF (gapw) THEN
616 1926 : IF (.NOT. ((nspins == 1 .AND. lr_triplet))) THEN
617 : CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
618 : p_env%hartree_local%ecoul_1c, &
619 : p_env%local_rho_set, &
620 1926 : para_env, tddft=.TRUE., core_2nd=.TRUE.)
621 :
622 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
623 : calculate_forces=.FALSE., &
624 1926 : local_rho_set=p_env%local_rho_set)
625 : END IF
626 : ! *** Add single atom contributions to the KS matrix ***
627 : ! remap pointer
628 1926 : ns = SIZE(p_env%kpp1)
629 1926 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
630 1926 : ns = SIZE(rho_ao)
631 1926 : psmat(1:ns, 1:1) => rho_ao(1:ns)
632 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
633 1926 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
634 7450 : ELSE IF (gapw_xc) THEN
635 560 : ns = SIZE(p_env%kpp1)
636 560 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
637 560 : ns = SIZE(rho_ao)
638 560 : psmat(1:ns, 1:1) => rho_ao(1:ns)
639 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
640 560 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
641 : END IF
642 :
643 : ! KG embedding, contribution of kinetic energy functional to kernel
644 9376 : IF (dft_control%qs_control%do_kg .AND. .NOT. (lr_triplet .OR. gapw .OR. gapw_xc)) THEN
645 16 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
646 :
647 10 : CALL qs_rho_get(rho1, rho_ao=rho1_ao)
648 10 : alpha = 1.0_dp
649 :
650 : ekin_mol = 0.0_dp
651 10 : CALL get_qs_env(qs_env, kg_env=kg_env)
652 : CALL kg_ekin_subset(qs_env=qs_env, &
653 : ks_matrix=p_env%kpp1, &
654 : ekin_mol=ekin_mol, &
655 : calc_force=.FALSE., &
656 : do_kernel=.TRUE., &
657 10 : pmat_ext=rho1_ao)
658 : END IF
659 : END IF
660 :
661 9376 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
662 9376 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
663 19770 : DO ispin = 1, nspins
664 19770 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
665 : END DO
666 9376 : DEALLOCATE (v_rspace_new)
667 9376 : IF (ASSOCIATED(v_xc_tau)) THEN
668 508 : DO ispin = 1, nspins
669 508 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
670 : END DO
671 254 : DEALLOCATE (v_xc_tau)
672 : END IF
673 :
674 9376 : CALL timestop(handle)
675 :
676 9376 : END SUBROUTINE apply_op_2_dft
677 :
678 : ! **************************************************************************************************
679 : !> \brief ...
680 : !> \param qs_env ...
681 : !> \param p_env ...
682 : ! **************************************************************************************************
683 206 : SUBROUTINE apply_op_2_xtb(qs_env, p_env)
684 : TYPE(qs_environment_type), POINTER :: qs_env
685 : TYPE(qs_p_env_type) :: p_env
686 :
687 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_xtb'
688 :
689 : INTEGER :: atom_a, handle, iatom, ikind, is, ispin, &
690 : na, natom, natorb, nkind, ns, nsgf, &
691 : nspins
692 : INTEGER, DIMENSION(25) :: lao
693 : INTEGER, DIMENSION(5) :: occ
694 : LOGICAL :: lr_triplet
695 206 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
696 206 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1
697 206 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
698 206 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmat, rho_ao
699 206 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_p1, matrix_s
700 : TYPE(dft_control_type), POINTER :: dft_control
701 : TYPE(linres_control_type), POINTER :: linres_control
702 : TYPE(mp_para_env_type), POINTER :: para_env
703 206 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
704 : TYPE(pw_env_type), POINTER :: pw_env
705 206 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
706 : TYPE(qs_kpp1_env_type), POINTER :: kpp1_env
707 : TYPE(qs_rho_type), POINTER :: rho, rho1
708 : TYPE(xtb_atom_type), POINTER :: xtb_kind
709 :
710 206 : CALL timeset(routineN, handle)
711 :
712 206 : CPASSERT(ASSOCIATED(p_env%kpp1_env))
713 206 : CPASSERT(ASSOCIATED(p_env%kpp1))
714 206 : kpp1_env => p_env%kpp1_env
715 :
716 206 : rho1 => p_env%rho1
717 206 : CPASSERT(ASSOCIATED(rho1))
718 :
719 : CALL get_qs_env(qs_env=qs_env, &
720 : pw_env=pw_env, &
721 : para_env=para_env, &
722 : rho=rho, &
723 : linres_control=linres_control, &
724 206 : dft_control=dft_control)
725 :
726 206 : CALL qs_rho_get(rho, rho_ao=rho_ao)
727 :
728 206 : lr_triplet = linres_control%lr_triplet
729 206 : CPASSERT(.NOT. lr_triplet)
730 :
731 206 : nspins = SIZE(p_env%kpp1)
732 :
733 508 : DO ispin = 1, nspins
734 508 : CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
735 : END DO
736 :
737 206 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
738 : ! Mulliken charges
739 202 : CALL get_qs_env(qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
740 202 : natom = SIZE(particle_set)
741 202 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
742 202 : CALL qs_rho_get(rho1, rho_ao_kp=matrix_p1)
743 1010 : ALLOCATE (mcharge(natom), charges(natom, 5))
744 606 : ALLOCATE (mcharge1(natom), charges1(natom, 5))
745 202 : charges = 0.0_dp
746 202 : charges1 = 0.0_dp
747 202 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
748 202 : nkind = SIZE(atomic_kind_set)
749 202 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
750 808 : ALLOCATE (aocg(nsgf, natom))
751 202 : aocg = 0.0_dp
752 606 : ALLOCATE (aocg1(nsgf, natom))
753 202 : aocg1 = 0.0_dp
754 202 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
755 202 : CALL ao_charges(matrix_p1, matrix_s, aocg1, para_env)
756 1872 : IF (nspins == 2) aocg1 = 0.5_dp*aocg1
757 704 : DO ikind = 1, nkind
758 502 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
759 502 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
760 502 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
761 3026 : DO iatom = 1, na
762 1820 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
763 10920 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
764 7306 : DO is = 1, natorb
765 4984 : ns = lao(is) + 1
766 4984 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
767 6804 : charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
768 : END DO
769 : END DO
770 : END DO
771 202 : DEALLOCATE (aocg, aocg1)
772 2022 : DO iatom = 1, natom
773 10920 : mcharge(iatom) = SUM(charges(iatom, :))
774 11122 : mcharge1(iatom) = SUM(charges1(iatom, :))
775 : END DO
776 : ! Coulomb Kernel
777 202 : pmat => matrix_p1(:, 1)
778 202 : CALL xtb_coulomb_hessian(qs_env, p_env%kpp1, charges1, mcharge1, mcharge, pmat)
779 : !
780 404 : DEALLOCATE (charges, mcharge, charges1, mcharge1)
781 : END IF
782 :
783 206 : CALL timestop(handle)
784 :
785 412 : END SUBROUTINE apply_op_2_xtb
786 :
787 : ! **************************************************************************************************
788 : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
789 : !> \param qs_env ...
790 : !> \param p_env ...
791 : !> \par History
792 : !> * 11.2019 adapted from tddfpt_apply_hfx
793 : ! **************************************************************************************************
794 19904 : SUBROUTINE apply_hfx(qs_env, p_env)
795 : TYPE(qs_environment_type), POINTER :: qs_env
796 : TYPE(qs_p_env_type) :: p_env
797 :
798 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_hfx'
799 :
800 : INTEGER :: handle, ispin, nspins
801 : LOGICAL :: do_hfx
802 : REAL(KIND=dp) :: alpha
803 : TYPE(cp_logger_type), POINTER :: logger
804 9952 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h1_mat, matrix_s, rho1_ao, work
805 : TYPE(dft_control_type), POINTER :: dft_control
806 : TYPE(section_vals_type), POINTER :: hfx_section, input
807 :
808 9952 : CALL timeset(routineN, handle)
809 :
810 9952 : logger => cp_get_default_logger()
811 :
812 : CALL get_qs_env(qs_env=qs_env, &
813 : input=input, &
814 : matrix_s=matrix_s, &
815 9952 : dft_control=dft_control)
816 9952 : nspins = dft_control%nspins
817 :
818 9952 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
819 9952 : CALL section_vals_get(hfx_section, explicit=do_hfx)
820 :
821 9952 : IF (do_hfx) THEN
822 :
823 4092 : IF (dft_control%do_admm) THEN
824 2252 : IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
825 0 : CPABORT("ADMM: Linear Response needs purification_method=none")
826 : END IF
827 2252 : IF (dft_control%admm_control%scaling_model /= do_admm_exch_scaling_none) THEN
828 0 : CPABORT("ADMM: Linear Response needs scaling_model=none")
829 : END IF
830 2252 : IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN
831 0 : CPABORT("ADMM: Linear Response needs admm_method=basis_projection")
832 : END IF
833 : !
834 2252 : rho1_ao => p_env%p1_admm
835 2252 : h1_mat => p_env%kpp1_admm
836 : ELSE
837 1840 : rho1_ao => p_env%p1
838 1840 : h1_mat => p_env%kpp1
839 : END IF
840 :
841 4092 : NULLIFY (work)
842 4092 : CALL dbcsr_allocate_matrix_set(work, nspins)
843 8588 : DO ispin = 1, nspins
844 4496 : ALLOCATE (work(ispin)%matrix)
845 4496 : CALL dbcsr_create(work(ispin)%matrix, template=h1_mat(ispin)%matrix)
846 4496 : CALL dbcsr_copy(work(ispin)%matrix, h1_mat(ispin)%matrix)
847 8588 : CALL dbcsr_set(work(ispin)%matrix, 0.0_dp)
848 : END DO
849 :
850 4092 : CALL hfx_matrix(work, rho1_ao, qs_env, hfx_section)
851 :
852 4092 : alpha = 2.0_dp
853 4092 : IF (nspins == 2) alpha = 1.0_dp
854 :
855 8588 : DO ispin = 1, nspins
856 8588 : CALL dbcsr_add(h1_mat(ispin)%matrix, work(ispin)%matrix, 1.0_dp, alpha)
857 : END DO
858 :
859 4092 : CALL dbcsr_deallocate_matrix_set(work)
860 :
861 : END IF
862 :
863 9952 : CALL timestop(handle)
864 :
865 9952 : END SUBROUTINE apply_hfx
866 :
867 : ! **************************************************************************************************
868 : !> \brief Add the hfx contributions to the Hamiltonian
869 : !>
870 : !> \param matrix_ks ...
871 : !> \param rho_ao ...
872 : !> \param qs_env ...
873 : !> \param hfx_sections ...
874 : !> \param external_x_data ...
875 : !> \param ex ...
876 : !> \note
877 : !> Simplified version of subroutine hfx_ks_matrix()
878 : ! **************************************************************************************************
879 4092 : SUBROUTINE hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections, external_x_data, ex)
880 : TYPE(dbcsr_p_type), DIMENSION(:), TARGET :: matrix_ks, rho_ao
881 : TYPE(qs_environment_type), POINTER :: qs_env
882 : TYPE(section_vals_type), POINTER :: hfx_sections
883 : TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, TARGET :: external_x_data
884 : REAL(KIND=dp), OPTIONAL :: ex
885 :
886 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_matrix'
887 :
888 : INTEGER :: handle, irep, ispin, mspin, n_rep_hf, &
889 : nspins
890 : LOGICAL :: distribute_fock_matrix, &
891 : hfx_treat_lsd_in_core, &
892 : s_mstruct_changed
893 : REAL(KIND=dp) :: eh1, ehfx
894 4092 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp
895 : TYPE(dft_control_type), POINTER :: dft_control
896 4092 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
897 : TYPE(mp_para_env_type), POINTER :: para_env
898 :
899 4092 : CALL timeset(routineN, handle)
900 :
901 4092 : NULLIFY (dft_control, para_env, matrix_ks_kp, rho_ao_kp, x_data)
902 :
903 : CALL get_qs_env(qs_env=qs_env, &
904 : dft_control=dft_control, &
905 : para_env=para_env, &
906 : s_mstruct_changed=s_mstruct_changed, &
907 4092 : x_data=x_data)
908 :
909 4092 : IF (PRESENT(external_x_data)) x_data => external_x_data
910 :
911 4092 : CPASSERT(dft_control%nimages == 1)
912 4092 : nspins = dft_control%nspins
913 :
914 4092 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
915 : CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
916 4092 : i_rep_section=1)
917 :
918 4092 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
919 4092 : distribute_fock_matrix = .TRUE.
920 :
921 4092 : mspin = 1
922 4092 : IF (hfx_treat_lsd_in_core) mspin = nspins
923 :
924 4092 : matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
925 4092 : rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
926 :
927 8184 : DO irep = 1, n_rep_hf
928 4092 : ehfx = 0.0_dp
929 :
930 8184 : IF (x_data(irep, 1)%do_hfx_ri) THEN
931 : CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, &
932 : rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
933 170 : nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
934 :
935 : ELSE
936 :
937 7844 : DO ispin = 1, mspin
938 : CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
939 3922 : s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
940 7844 : ehfx = ehfx + eh1
941 : END DO
942 :
943 : END IF
944 : END DO
945 :
946 : ! Export energy
947 4092 : IF (PRESENT(ex)) ex = ehfx
948 :
949 4092 : CALL timestop(handle)
950 :
951 4092 : END SUBROUTINE hfx_matrix
952 :
953 : ! **************************************************************************************************
954 : !> \brief ...
955 : !> \param qs_env ...
956 : !> \param p_env ...
957 : ! **************************************************************************************************
958 9952 : SUBROUTINE apply_xc_admm(qs_env, p_env)
959 : TYPE(qs_environment_type), POINTER :: qs_env
960 : TYPE(qs_p_env_type) :: p_env
961 :
962 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_xc_admm'
963 :
964 : CHARACTER(LEN=default_string_length) :: basis_type
965 : INTEGER :: handle, ispin, ns, nspins
966 : REAL(KIND=dp) :: alpha
967 : TYPE(admm_type), POINTER :: admm_env
968 : TYPE(dbcsr_p_type) :: xcmat
969 9952 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
970 9952 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
971 : TYPE(dft_control_type), POINTER :: dft_control
972 : TYPE(linres_control_type), POINTER :: linres_control
973 : TYPE(mp_para_env_type), POINTER :: para_env
974 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
975 9952 : POINTER :: sab_aux_fit
976 : TYPE(pw_env_type), POINTER :: pw_env
977 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
978 9952 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
979 9952 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
980 : TYPE(qs_kpp1_env_type), POINTER :: kpp1_env
981 : TYPE(qs_rho_type), POINTER :: rho_aux
982 9952 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
983 : TYPE(section_vals_type), POINTER :: xc_section
984 : TYPE(task_list_type), POINTER :: task_list
985 :
986 9952 : CALL timeset(routineN, handle)
987 :
988 9952 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
989 :
990 9952 : IF (dft_control%do_admm) THEN
991 2252 : CALL get_qs_env(qs_env, admm_env=admm_env)
992 2252 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
993 : ! nothing to do
994 : ELSE
995 1344 : CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
996 1344 : CPASSERT(.NOT. dft_control%qs_control%lrigpw)
997 1344 : CPASSERT(.NOT. linres_control%lr_triplet)
998 :
999 1344 : nspins = dft_control%nspins
1000 :
1001 : ! AUX basis contribution
1002 1344 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1003 1344 : CPASSERT(ASSOCIATED(pw_env))
1004 1344 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1005 : ! calculate the xc potential
1006 1344 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s)
1007 1344 : ALLOCATE (xcmat%matrix)
1008 1344 : CALL dbcsr_create(xcmat%matrix, template=matrix_s(1)%matrix)
1009 :
1010 1344 : NULLIFY (v_xc, v_xc_tau)
1011 1344 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux)
1012 1344 : xc_section => admm_env%xc_section_aux
1013 1344 : kpp1_env => p_env%kpp1_env
1014 :
1015 1344 : NULLIFY (rho_atom_set, rho1_atom_set)
1016 1344 : basis_type = "AUX_FIT"
1017 1344 : CALL get_qs_env(qs_env, para_env=para_env, qs_kind_set=kind_set)
1018 1344 : CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
1019 1344 : IF (admm_env%do_gapw) THEN
1020 396 : kind_set => admm_env%admm_gapw_env%admm_kind_set
1021 : CALL prepare_gapw_den(qs_env, local_rho_set=p_env%local_rho_set_admm, &
1022 396 : do_rho0=.FALSE., kind_set_external=kind_set)
1023 396 : rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
1024 396 : rho1_atom_set => p_env%local_rho_set_admm%rho_atom_set
1025 396 : basis_type = "AUX_FIT_SOFT"
1026 396 : task_list => admm_env%admm_gapw_env%task_list
1027 : END IF
1028 :
1029 : CALL qs_fxc_apply(qs_env, kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, p_env%rho1_admm, &
1030 : rho_atom_set, xc_section, admm_env%do_gapw, v_xc, v_xc_tau, rho1_atom_set, &
1031 1344 : kind_set_external=kind_set)
1032 1344 : IF (ASSOCIATED(v_xc_tau)) THEN
1033 0 : CPABORT("Meta-GGA ADMM functionals not yet supported!")
1034 : END IF
1035 :
1036 1344 : alpha = 1.0_dp
1037 1344 : IF (nspins == 1) alpha = 2.0_dp
1038 :
1039 2814 : DO ispin = 1, nspins
1040 1470 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
1041 1470 : CALL dbcsr_copy(xcmat%matrix, matrix_s(1)%matrix)
1042 1470 : CALL dbcsr_set(xcmat%matrix, 0.0_dp)
1043 : CALL integrate_v_rspace(v_rspace=v_xc(ispin), hmat=xcmat, qs_env=qs_env, &
1044 : calculate_forces=.FALSE., basis_type=basis_type, &
1045 1470 : task_list_external=task_list)
1046 2814 : CALL dbcsr_add(p_env%kpp1_admm(ispin)%matrix, xcmat%matrix, 1.0_dp, alpha)
1047 : END DO
1048 :
1049 1344 : IF (admm_env%do_gapw) THEN
1050 396 : CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
1051 396 : ns = SIZE(p_env%kpp1_admm)
1052 396 : ksmat(1:ns, 1:1) => p_env%kpp1_admm(1:ns)
1053 396 : psmat(1:ns, 1:1) => p_env%p1_admm(1:ns)
1054 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
1055 : rho_atom_external=p_env%local_rho_set_admm%rho_atom_set, &
1056 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
1057 : oce_external=admm_env%admm_gapw_env%oce, &
1058 396 : sab_external=sab_aux_fit)
1059 : END IF
1060 :
1061 2814 : DO ispin = 1, nspins
1062 2814 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
1063 : END DO
1064 1344 : DEALLOCATE (v_xc)
1065 1344 : CALL dbcsr_deallocate_matrix(xcmat%matrix)
1066 :
1067 : END IF
1068 : END IF
1069 :
1070 9952 : CALL timestop(handle)
1071 :
1072 9952 : END SUBROUTINE apply_xc_admm
1073 :
1074 : END MODULE qs_linres_kernel
|
