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