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