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 8870 : 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 8870 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
143 8870 : IF (dft_control%qs_control%semi_empirical) THEN
144 0 : CPABORT("Linear response not available with SE methods")
145 8870 : ELSEIF (dft_control%qs_control%dftb) THEN
146 0 : CPABORT("Linear response not available with DFTB")
147 8870 : ELSEIF (dft_control%qs_control%xtb) THEN
148 110 : CALL apply_op_2_xtb(qs_env, p_env)
149 : ELSE
150 8760 : CALL apply_op_2_dft(qs_env, p_env)
151 8760 : CALL apply_hfx(qs_env, p_env)
152 8760 : CALL apply_xc_admm(qs_env, p_env)
153 8760 : IF (dft_control%do_admm) CALL p_env_finish_kpp1(qs_env, p_env)
154 : END IF
155 :
156 18724 : DO ispin = 1, SIZE(c0)
157 9854 : 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 18724 : ncol=ncol, alpha=1.0_dp, beta=1.0_dp)
162 : END DO
163 :
164 8870 : END SUBROUTINE apply_op_2
165 :
166 : ! **************************************************************************************************
167 : !> \brief ...
168 : !> \param qs_env ...
169 : !> \param p_env ...
170 : ! **************************************************************************************************
171 8760 : 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 8760 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
184 : TYPE(cp_logger_type), POINTER :: logger
185 8760 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: k1mat, matrix_s, rho1_ao, rho_ao
186 8760 : 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 8760 : 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 8760 : 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 8760 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho_r, tau1_r, v_rspace_new, &
201 8760 : 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 8760 : 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 8760 : CALL timeset(routineN, handle)
211 :
212 8760 : NULLIFY (auxbas_pw_pool, pw_env, v_rspace_new, para_env, rho1_r, &
213 8760 : v_xc, rho1_ao, rho_ao, poisson_env, input, rho, dft_control, &
214 8760 : logger, rho1_g, v_xc_tau)
215 8760 : logger => cp_get_default_logger()
216 :
217 8760 : energy_hartree = 0.0_dp
218 8760 : energy_hartree_1c = 0.0_dp
219 :
220 8760 : CPASSERT(ASSOCIATED(p_env%kpp1))
221 8760 : CPASSERT(ASSOCIATED(p_env%kpp1_env))
222 8760 : 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 8760 : dft_control=dft_control)
234 :
235 8760 : gapw = dft_control%qs_control%gapw
236 8760 : gapw_xc = dft_control%qs_control%gapw_xc
237 8760 : lr_triplet = linres_control%lr_triplet
238 :
239 8760 : rho1 => p_env%rho1
240 8760 : rho1_xc => p_env%rho1_xc
241 8760 : CPASSERT(ASSOCIATED(rho1))
242 8760 : IF (gapw_xc) THEN
243 496 : CPASSERT(ASSOCIATED(rho1_xc))
244 : END IF
245 :
246 8760 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
247 8760 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
248 :
249 8760 : nspins = SIZE(p_env%kpp1)
250 8760 : lrigpw = dft_control%qs_control%lrigpw
251 8760 : 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 8760 : IF (.NOT. ASSOCIATED(kpp1_env%v_ao)) THEN
259 1154 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
260 1154 : CALL dbcsr_allocate_matrix_set(kpp1_env%v_ao, nspins)
261 2440 : DO ispin = 1, nspins
262 1286 : ALLOCATE (kpp1_env%v_ao(ispin)%matrix)
263 : CALL dbcsr_copy(kpp1_env%v_ao(ispin)%matrix, matrix_s(1)%matrix, &
264 2440 : name="kpp1%v_ao-"//ADJUSTL(cp_to_string(ispin)))
265 : END DO
266 : END IF
267 :
268 8760 : IF (dft_control%do_admm) THEN
269 1920 : xc_section => admm_env%xc_section_primary
270 : ELSE
271 6840 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
272 : END IF
273 :
274 : ! gets the tmp grids
275 8760 : CPASSERT(ASSOCIATED(pw_env))
276 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
277 8760 : poisson_env=poisson_env)
278 8760 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
279 8760 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
280 :
281 8760 : IF (gapw .OR. gapw_xc) &
282 2164 : CALL prepare_gapw_den(qs_env, p_env%local_rho_set, do_rho0=(.NOT. gapw_xc))
283 :
284 : ! *** calculate the hartree potential on the total density ***
285 8760 : CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
286 :
287 8760 : CALL qs_rho_get(rho1, rho_g=rho1_g)
288 8760 : CALL pw_copy(rho1_g(1), rho1_tot_gspace)
289 9744 : DO ispin = 2, nspins
290 9744 : CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
291 : END DO
292 8760 : IF (gapw) THEN
293 1668 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
294 1668 : IF (ASSOCIATED(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
295 0 : CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho1_tot_gspace)
296 : END IF
297 : END IF
298 :
299 8760 : IF (.NOT. (nspins == 1 .AND. lr_triplet)) THEN
300 : CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
301 : energy_hartree, &
302 8760 : v_hartree_gspace)
303 8760 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
304 : END IF
305 :
306 8760 : CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
307 :
308 : ! *** calculate the xc potential ***
309 8760 : NULLIFY (rho1a)
310 8760 : IF (gapw_xc) THEN
311 496 : rho0 => rho_xc
312 496 : rho1a => rho1_xc
313 : ELSE
314 8264 : rho0 => rho
315 8264 : rho1a => rho1
316 : END IF
317 :
318 8760 : deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
319 8760 : NULLIFY (v_xc_tau)
320 8760 : IF (deriv2_analytic) THEN
321 8700 : CALL qs_rho_get(rho1a, rho_r=rho1_r, tau_r=tau1_r)
322 8700 : CALL get_qs_env(qs_env, xcint_weights=weights)
323 : CALL qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, weights, auxbas_pw_pool, &
324 8700 : lr_triplet, v_xc, v_xc_tau)
325 8700 : IF (gapw .OR. gapw_xc) THEN
326 2164 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
327 2164 : rho1_atom_set => p_env%local_rho_set%rho_atom_set
328 : CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
329 2164 : do_triplet=lr_triplet)
330 : END IF
331 : ELSE
332 60 : CALL qs_fxc_fdiff(ks_env, rho0, rho1a, xc_section, 6, lr_triplet, v_xc, v_xc_tau)
333 60 : CPASSERT((.NOT. gapw) .AND. (.NOT. gapw_xc))
334 : END IF
335 :
336 8760 : v_rspace_new => v_xc
337 8760 : NULLIFY (v_xc)
338 :
339 8760 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
340 18504 : DO ispin = 1, nspins
341 9744 : CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
342 18504 : IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
343 : END DO
344 :
345 : ! ADMM Correction
346 8760 : IF (dft_control%do_admm) THEN
347 1920 : IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
348 1112 : IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
349 164 : CPASSERT(.NOT. lr_triplet)
350 164 : CALL get_qs_env(qs_env, xcint_weights=weights)
351 164 : xc_section_aux => admm_env%xc_section_aux
352 164 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
353 164 : CALL qs_rho_get(rho_aux, rho_r=rho_r)
354 3772 : ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
355 : CALL xc_prep_2nd_deriv(kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm, &
356 : rho_r, auxbas_pw_pool, weights, &
357 164 : xc_section=xc_section_aux)
358 : END IF
359 : END IF
360 : END IF
361 :
362 : !-------------------------------!
363 : ! Add both hartree and xc terms !
364 : !-------------------------------!
365 18504 : DO ispin = 1, nspins
366 9744 : CALL dbcsr_set(kpp1_env%v_ao(ispin)%matrix, 0.0_dp)
367 :
368 9744 : IF (gapw_xc) THEN
369 : ! XC and Hartree are integrated separatedly
370 : ! XC uses the soft basis set only
371 :
372 512 : IF (nspins == 1) THEN
373 :
374 480 : IF (.NOT. (lr_triplet)) THEN
375 480 : CALL pw_scale(v_rspace_new(1), 2.0_dp)
376 480 : IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1), 2.0_dp)
377 : END IF
378 480 : CALL qs_rho_get(rho1, rho_ao=rho1_ao)
379 : ! remove kpp1_env%v_ao and work directly on k_p_p1 ?
380 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
381 : pmat=rho1_ao(ispin), &
382 : hmat=kpp1_env%v_ao(ispin), &
383 : qs_env=qs_env, &
384 480 : calculate_forces=.FALSE., gapw=gapw_xc)
385 :
386 480 : IF (ASSOCIATED(v_xc_tau)) THEN
387 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
388 : pmat=rho1_ao(ispin), &
389 : hmat=kpp1_env%v_ao(ispin), &
390 : qs_env=qs_env, &
391 : compute_tau=.TRUE., &
392 0 : calculate_forces=.FALSE., gapw=gapw_xc)
393 : END IF
394 :
395 : ! add hartree only for SINGLETS
396 480 : IF (.NOT. lr_triplet) THEN
397 480 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp, 0.0_dp)
398 :
399 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
400 : pmat=rho_ao(ispin), &
401 : hmat=kpp1_env%v_ao(ispin), &
402 : qs_env=qs_env, &
403 480 : calculate_forces=.FALSE., gapw=gapw)
404 : END IF
405 : ELSE
406 : ! remove kpp1_env%v_ao and work directly on k_p_p1 ?
407 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
408 : pmat=rho_ao(ispin), &
409 : hmat=kpp1_env%v_ao(ispin), &
410 : qs_env=qs_env, &
411 32 : calculate_forces=.FALSE., gapw=gapw_xc)
412 :
413 32 : IF (ASSOCIATED(v_xc_tau)) THEN
414 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
415 : pmat=rho_ao(ispin), &
416 : hmat=kpp1_env%v_ao(ispin), &
417 : qs_env=qs_env, &
418 : compute_tau=.TRUE., &
419 0 : calculate_forces=.FALSE., gapw=gapw_xc)
420 : END IF
421 :
422 32 : CALL pw_copy(v_hartree_rspace, v_rspace_new(ispin))
423 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
424 : pmat=rho_ao(ispin), &
425 : hmat=kpp1_env%v_ao(ispin), &
426 : qs_env=qs_env, &
427 32 : calculate_forces=.FALSE., gapw=gapw)
428 : END IF
429 :
430 : ELSE
431 :
432 9232 : IF (nspins == 1) THEN
433 7296 : IF (.NOT. (lr_triplet)) THEN
434 7296 : CALL pw_scale(v_rspace_new(1), 2.0_dp)
435 7296 : IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1), 2.0_dp)
436 : END IF
437 : ! add hartree only for SINGLETS
438 : !IF (res_etype == tddfpt_singlet) THEN
439 7296 : IF (.NOT. lr_triplet) THEN
440 7296 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(1), 2.0_dp)
441 : END IF
442 : ELSE
443 1936 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin), 1.0_dp)
444 : END IF
445 :
446 9232 : IF (lrigpw) THEN
447 72 : IF (ASSOCIATED(v_xc_tau)) &
448 0 : CPABORT("metaGGA-functionals not supported with LRI!")
449 :
450 72 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
451 72 : CALL get_qs_env(qs_env, nkind=nkind)
452 216 : DO ikind = 1, nkind
453 43008 : lri_v_int(ikind)%v_int = 0.0_dp
454 : END DO
455 : CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
456 72 : lri_v_int, .FALSE., "LRI_AUX")
457 216 : DO ikind = 1, nkind
458 85800 : CALL para_env%sum(lri_v_int(ikind)%v_int)
459 : END DO
460 144 : ALLOCATE (k1mat(1))
461 72 : k1mat(1)%matrix => kpp1_env%v_ao(ispin)%matrix
462 72 : IF (lri_env%exact_1c_terms) THEN
463 : CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), k1mat(1)%matrix, &
464 0 : rho_ao(ispin)%matrix, qs_env, .FALSE., "ORB")
465 : END IF
466 72 : CALL calculate_lri_ks_matrix(lri_env, lri_v_int, k1mat, atomic_kind_set)
467 72 : DEALLOCATE (k1mat)
468 : ELSE
469 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
470 : pmat=rho_ao(ispin), &
471 : hmat=kpp1_env%v_ao(ispin), &
472 : qs_env=qs_env, &
473 9160 : calculate_forces=.FALSE., gapw=gapw)
474 :
475 9160 : IF (ASSOCIATED(v_xc_tau)) THEN
476 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
477 : pmat=rho_ao(ispin), &
478 : hmat=kpp1_env%v_ao(ispin), &
479 : qs_env=qs_env, &
480 : compute_tau=.TRUE., &
481 196 : calculate_forces=.FALSE., gapw=gapw)
482 : END IF
483 : END IF
484 :
485 : END IF
486 :
487 18504 : CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, kpp1_env%v_ao(ispin)%matrix)
488 : END DO
489 :
490 8760 : IF (gapw) THEN
491 1668 : IF (.NOT. ((nspins == 1 .AND. lr_triplet))) THEN
492 : CALL Vh_1c_gg_integrals(qs_env, energy_hartree_1c, &
493 : p_env%hartree_local%ecoul_1c, &
494 : p_env%local_rho_set, &
495 1668 : para_env, tddft=.TRUE., core_2nd=.TRUE.)
496 :
497 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
498 : calculate_forces=.FALSE., &
499 1668 : local_rho_set=p_env%local_rho_set)
500 : END IF
501 : ! *** Add single atom contributions to the KS matrix ***
502 : ! remap pointer
503 1668 : ns = SIZE(p_env%kpp1)
504 1668 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
505 1668 : ns = SIZE(rho_ao)
506 1668 : psmat(1:ns, 1:1) => rho_ao(1:ns)
507 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
508 1668 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
509 7092 : ELSEIF (gapw_xc) THEN
510 496 : ns = SIZE(p_env%kpp1)
511 496 : ksmat(1:ns, 1:1) => p_env%kpp1(1:ns)
512 496 : ns = SIZE(rho_ao)
513 496 : psmat(1:ns, 1:1) => rho_ao(1:ns)
514 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
515 496 : rho_atom_external=p_env%local_rho_set%rho_atom_set)
516 : END IF
517 :
518 : ! KG embedding, contribution of kinetic energy functional to kernel
519 8760 : IF (dft_control%qs_control%do_kg .AND. .NOT. (lr_triplet .OR. gapw .OR. gapw_xc)) THEN
520 16 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
521 :
522 10 : CALL qs_rho_get(rho1, rho_ao=rho1_ao)
523 10 : alpha = 1.0_dp
524 :
525 : ekin_mol = 0.0_dp
526 10 : CALL get_qs_env(qs_env, kg_env=kg_env)
527 : CALL kg_ekin_subset(qs_env=qs_env, &
528 : ks_matrix=p_env%kpp1, &
529 : ekin_mol=ekin_mol, &
530 : calc_force=.FALSE., &
531 : do_kernel=.TRUE., &
532 10 : pmat_ext=rho1_ao)
533 : END IF
534 : END IF
535 :
536 8760 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
537 8760 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
538 18504 : DO ispin = 1, nspins
539 18504 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
540 : END DO
541 8760 : DEALLOCATE (v_rspace_new)
542 8760 : IF (ASSOCIATED(v_xc_tau)) THEN
543 392 : DO ispin = 1, nspins
544 392 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
545 : END DO
546 196 : DEALLOCATE (v_xc_tau)
547 : END IF
548 :
549 8760 : CALL timestop(handle)
550 :
551 8760 : END SUBROUTINE apply_op_2_dft
552 :
553 : ! **************************************************************************************************
554 : !> \brief ...
555 : !> \param qs_env ...
556 : !> \param p_env ...
557 : ! **************************************************************************************************
558 110 : SUBROUTINE apply_op_2_xtb(qs_env, p_env)
559 : TYPE(qs_environment_type), POINTER :: qs_env
560 : TYPE(qs_p_env_type) :: p_env
561 :
562 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_2_xtb'
563 :
564 : INTEGER :: atom_a, handle, iatom, ikind, is, ispin, &
565 : na, natom, natorb, nkind, ns, nsgf, &
566 : nspins
567 : INTEGER, DIMENSION(25) :: lao
568 : INTEGER, DIMENSION(5) :: occ
569 : LOGICAL :: lr_triplet
570 110 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
571 110 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1
572 110 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
573 110 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
574 110 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_p1, matrix_s
575 : TYPE(dft_control_type), POINTER :: dft_control
576 : TYPE(linres_control_type), POINTER :: linres_control
577 : TYPE(mp_para_env_type), POINTER :: para_env
578 110 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
579 : TYPE(pw_env_type), POINTER :: pw_env
580 110 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
581 : TYPE(qs_kpp1_env_type), POINTER :: kpp1_env
582 : TYPE(qs_rho_type), POINTER :: rho, rho1
583 : TYPE(xtb_atom_type), POINTER :: xtb_kind
584 :
585 110 : CALL timeset(routineN, handle)
586 :
587 110 : CPASSERT(ASSOCIATED(p_env%kpp1_env))
588 110 : CPASSERT(ASSOCIATED(p_env%kpp1))
589 110 : kpp1_env => p_env%kpp1_env
590 :
591 110 : rho1 => p_env%rho1
592 110 : CPASSERT(ASSOCIATED(rho1))
593 :
594 : CALL get_qs_env(qs_env=qs_env, &
595 : pw_env=pw_env, &
596 : para_env=para_env, &
597 : rho=rho, &
598 : linres_control=linres_control, &
599 110 : dft_control=dft_control)
600 :
601 110 : CALL qs_rho_get(rho, rho_ao=rho_ao)
602 :
603 110 : lr_triplet = linres_control%lr_triplet
604 110 : CPASSERT(.NOT. lr_triplet)
605 :
606 110 : nspins = SIZE(p_env%kpp1)
607 :
608 220 : DO ispin = 1, nspins
609 220 : CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
610 : END DO
611 :
612 110 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
613 : ! Mulliken charges
614 106 : CALL get_qs_env(qs_env, particle_set=particle_set, matrix_s_kp=matrix_s)
615 106 : natom = SIZE(particle_set)
616 106 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
617 106 : CALL qs_rho_get(rho1, rho_ao_kp=matrix_p1)
618 530 : ALLOCATE (mcharge(natom), charges(natom, 5))
619 318 : ALLOCATE (mcharge1(natom), charges1(natom, 5))
620 8066 : charges = 0.0_dp
621 8066 : charges1 = 0.0_dp
622 106 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
623 106 : nkind = SIZE(atomic_kind_set)
624 106 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
625 424 : ALLOCATE (aocg(nsgf, natom))
626 7536 : aocg = 0.0_dp
627 318 : ALLOCATE (aocg1(nsgf, natom))
628 7536 : aocg1 = 0.0_dp
629 106 : CALL ao_charges(matrix_p, matrix_s, aocg, para_env)
630 106 : CALL ao_charges(matrix_p1, matrix_s, aocg1, para_env)
631 370 : DO ikind = 1, nkind
632 264 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
633 264 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
634 264 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
635 2120 : DO iatom = 1, na
636 1486 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
637 8916 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
638 5782 : DO is = 1, natorb
639 4032 : ns = lao(is) + 1
640 4032 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
641 5518 : charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
642 : END DO
643 : END DO
644 : END DO
645 106 : DEALLOCATE (aocg, aocg1)
646 1592 : DO iatom = 1, natom
647 8916 : mcharge(iatom) = SUM(charges(iatom, :))
648 9022 : mcharge1(iatom) = SUM(charges1(iatom, :))
649 : END DO
650 : ! Coulomb Kernel
651 106 : CALL xtb_coulomb_hessian(qs_env, p_env%kpp1, charges1, mcharge1, mcharge)
652 : !
653 212 : DEALLOCATE (charges, mcharge, charges1, mcharge1)
654 : END IF
655 :
656 110 : CALL timestop(handle)
657 :
658 220 : END SUBROUTINE apply_op_2_xtb
659 :
660 : ! **************************************************************************************************
661 : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
662 : !> \param qs_env ...
663 : !> \param p_env ...
664 : !> \par History
665 : !> * 11.2019 adapted from tddfpt_apply_hfx
666 : ! **************************************************************************************************
667 18672 : SUBROUTINE apply_hfx(qs_env, p_env)
668 : TYPE(qs_environment_type), POINTER :: qs_env
669 : TYPE(qs_p_env_type) :: p_env
670 :
671 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_hfx'
672 :
673 : INTEGER :: handle, ispin, nspins
674 : LOGICAL :: do_hfx
675 : REAL(KIND=dp) :: alpha
676 : TYPE(cp_logger_type), POINTER :: logger
677 9336 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h1_mat, matrix_s, rho1_ao, work
678 : TYPE(dft_control_type), POINTER :: dft_control
679 : TYPE(section_vals_type), POINTER :: hfx_section, input
680 :
681 9336 : CALL timeset(routineN, handle)
682 :
683 9336 : logger => cp_get_default_logger()
684 :
685 : CALL get_qs_env(qs_env=qs_env, &
686 : input=input, &
687 : matrix_s=matrix_s, &
688 9336 : dft_control=dft_control)
689 9336 : nspins = dft_control%nspins
690 :
691 9336 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
692 9336 : CALL section_vals_get(hfx_section, explicit=do_hfx)
693 :
694 9336 : IF (do_hfx) THEN
695 :
696 3824 : IF (dft_control%do_admm) THEN
697 2052 : IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
698 0 : CPABORT("ADMM: Linear Response needs purification_method=none")
699 : END IF
700 2052 : IF (dft_control%admm_control%scaling_model /= do_admm_exch_scaling_none) THEN
701 0 : CPABORT("ADMM: Linear Response needs scaling_model=none")
702 : END IF
703 2052 : IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN
704 0 : CPABORT("ADMM: Linear Response needs admm_method=basis_projection")
705 : END IF
706 : !
707 2052 : rho1_ao => p_env%p1_admm
708 2052 : h1_mat => p_env%kpp1_admm
709 : ELSE
710 1772 : rho1_ao => p_env%p1
711 1772 : h1_mat => p_env%kpp1
712 : END IF
713 :
714 3824 : NULLIFY (work)
715 3824 : CALL dbcsr_allocate_matrix_set(work, nspins)
716 8024 : DO ispin = 1, nspins
717 4200 : ALLOCATE (work(ispin)%matrix)
718 4200 : CALL dbcsr_create(work(ispin)%matrix, template=h1_mat(ispin)%matrix)
719 4200 : CALL dbcsr_copy(work(ispin)%matrix, h1_mat(ispin)%matrix)
720 8024 : CALL dbcsr_set(work(ispin)%matrix, 0.0_dp)
721 : END DO
722 :
723 3824 : CALL hfx_matrix(work, rho1_ao, qs_env, hfx_section)
724 :
725 3824 : alpha = 2.0_dp
726 3824 : IF (nspins == 2) alpha = 1.0_dp
727 :
728 8024 : DO ispin = 1, nspins
729 8024 : CALL dbcsr_add(h1_mat(ispin)%matrix, work(ispin)%matrix, 1.0_dp, alpha)
730 : END DO
731 :
732 3824 : CALL dbcsr_deallocate_matrix_set(work)
733 :
734 : END IF
735 :
736 9336 : CALL timestop(handle)
737 :
738 9336 : END SUBROUTINE apply_hfx
739 :
740 : ! **************************************************************************************************
741 : !> \brief Add the hfx contributions to the Hamiltonian
742 : !>
743 : !> \param matrix_ks ...
744 : !> \param rho_ao ...
745 : !> \param qs_env ...
746 : !> \param hfx_sections ...
747 : !> \param external_x_data ...
748 : !> \param ex ...
749 : !> \note
750 : !> Simplified version of subroutine hfx_ks_matrix()
751 : ! **************************************************************************************************
752 3824 : SUBROUTINE hfx_matrix(matrix_ks, rho_ao, qs_env, hfx_sections, external_x_data, ex)
753 : TYPE(dbcsr_p_type), DIMENSION(:), TARGET :: matrix_ks, rho_ao
754 : TYPE(qs_environment_type), POINTER :: qs_env
755 : TYPE(section_vals_type), POINTER :: hfx_sections
756 : TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, TARGET :: external_x_data
757 : REAL(KIND=dp), OPTIONAL :: ex
758 :
759 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_matrix'
760 :
761 : INTEGER :: handle, irep, ispin, mspin, n_rep_hf, &
762 : nspins
763 : LOGICAL :: distribute_fock_matrix, &
764 : hfx_treat_lsd_in_core, &
765 : s_mstruct_changed
766 : REAL(KIND=dp) :: eh1, ehfx
767 3824 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp
768 : TYPE(dft_control_type), POINTER :: dft_control
769 3824 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
770 : TYPE(mp_para_env_type), POINTER :: para_env
771 :
772 3824 : CALL timeset(routineN, handle)
773 :
774 3824 : NULLIFY (dft_control, para_env, matrix_ks_kp, rho_ao_kp, x_data)
775 :
776 : CALL get_qs_env(qs_env=qs_env, &
777 : dft_control=dft_control, &
778 : para_env=para_env, &
779 : s_mstruct_changed=s_mstruct_changed, &
780 3824 : x_data=x_data)
781 :
782 3824 : IF (PRESENT(external_x_data)) x_data => external_x_data
783 :
784 3824 : CPASSERT(dft_control%nimages == 1)
785 3824 : nspins = dft_control%nspins
786 :
787 3824 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
788 : CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
789 3824 : i_rep_section=1)
790 :
791 3824 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
792 3824 : distribute_fock_matrix = .TRUE.
793 :
794 3824 : mspin = 1
795 3824 : IF (hfx_treat_lsd_in_core) mspin = nspins
796 :
797 3824 : matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
798 3824 : rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
799 :
800 7648 : DO irep = 1, n_rep_hf
801 3824 : ehfx = 0.0_dp
802 :
803 7648 : IF (x_data(irep, 1)%do_hfx_ri) THEN
804 : CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, &
805 : rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
806 170 : nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
807 :
808 : ELSE
809 :
810 7308 : DO ispin = 1, mspin
811 : CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
812 3654 : s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
813 7308 : ehfx = ehfx + eh1
814 : END DO
815 :
816 : END IF
817 : END DO
818 :
819 : ! Export energy
820 3824 : IF (PRESENT(ex)) ex = ehfx
821 :
822 3824 : CALL timestop(handle)
823 :
824 3824 : END SUBROUTINE hfx_matrix
825 :
826 : ! **************************************************************************************************
827 : !> \brief ...
828 : !> \param qs_env ...
829 : !> \param p_env ...
830 : ! **************************************************************************************************
831 18672 : SUBROUTINE apply_xc_admm(qs_env, p_env)
832 : TYPE(qs_environment_type), POINTER :: qs_env
833 : TYPE(qs_p_env_type) :: p_env
834 :
835 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_xc_admm'
836 :
837 : CHARACTER(LEN=default_string_length) :: basis_type
838 : INTEGER :: handle, ispin, ns, nspins
839 : INTEGER, DIMENSION(2, 3) :: bo
840 : LOGICAL :: gapw, lsd
841 : REAL(KIND=dp) :: alpha
842 : TYPE(admm_type), POINTER :: admm_env
843 : TYPE(dbcsr_p_type) :: xcmat
844 9336 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
845 9336 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
846 : TYPE(dft_control_type), POINTER :: dft_control
847 : TYPE(linres_control_type), POINTER :: linres_control
848 : TYPE(mp_para_env_type), POINTER :: para_env
849 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
850 9336 : POINTER :: sab_aux_fit
851 9336 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_aux_g
852 : TYPE(pw_env_type), POINTER :: pw_env
853 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
854 18672 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_aux_r, tau_pw, v_xc, v_xc_tau
855 : TYPE(pw_r3d_rs_type), POINTER :: weights
856 18672 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
857 : TYPE(section_vals_type), POINTER :: xc_fun_section, xc_section
858 : TYPE(task_list_type), POINTER :: task_list
859 : TYPE(xc_rho_cflags_type) :: needs
860 : TYPE(xc_rho_set_type) :: rho1_set
861 :
862 9336 : CALL timeset(routineN, handle)
863 :
864 9336 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
865 :
866 9336 : IF (dft_control%do_admm) THEN
867 2052 : IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
868 : ! nothing to do
869 : ELSE
870 1182 : CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
871 1182 : CPASSERT(.NOT. dft_control%qs_control%lrigpw)
872 1182 : CPASSERT(.NOT. linres_control%lr_triplet)
873 :
874 1182 : nspins = dft_control%nspins
875 :
876 : ! AUX basis contribution
877 1182 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
878 1182 : CPASSERT(ASSOCIATED(pw_env))
879 1182 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
880 1182 : NULLIFY (tau_pw)
881 : ! calculate the xc potential
882 1182 : lsd = (nspins == 2)
883 1182 : CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s)
884 1182 : ALLOCATE (xcmat%matrix)
885 1182 : CALL dbcsr_create(xcmat%matrix, template=matrix_s(1)%matrix)
886 :
887 1182 : CALL get_qs_env(qs_env, admm_env=admm_env)
888 1182 : gapw = admm_env%do_gapw
889 :
890 1182 : NULLIFY (weights)
891 1182 : CALL get_qs_env(qs_env=qs_env, xcint_weights=weights)
892 :
893 1182 : CALL qs_rho_get(p_env%rho1_admm, rho_r=rho1_aux_r, rho_g=rho1_aux_g)
894 1182 : xc_section => admm_env%xc_section_aux
895 11820 : bo = rho1_aux_r(1)%pw_grid%bounds_local
896 : ! create the place where to store the argument for the functionals
897 : CALL xc_rho_set_create(rho1_set, bo, &
898 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
899 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
900 1182 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
901 :
902 1182 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
903 1182 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
904 :
905 : ! calculate the arguments needed by the functionals
906 : CALL xc_rho_set_update(rho1_set, rho1_aux_r, rho1_aux_g, tau_pw, needs, &
907 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
908 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
909 1182 : auxbas_pw_pool)
910 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, p_env%kpp1_env%deriv_set_admm, p_env%kpp1_env%rho_set_admm, &
911 : rho1_aux_r, rho1_aux_g, tau_pw, auxbas_pw_pool, weights, gapw=.FALSE., &
912 1182 : xc_section=xc_section)
913 1182 : IF (ASSOCIATED(v_xc_tau)) THEN
914 0 : CPABORT("Meta-GGA ADMM functionals not yet supported!")
915 : END IF
916 1182 : CALL xc_rho_set_release(rho1_set)
917 :
918 1182 : basis_type = "AUX_FIT"
919 1182 : CALL get_qs_env(qs_env, para_env=para_env)
920 1182 : CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
921 1182 : IF (admm_env%do_gapw) THEN
922 : CALL prepare_gapw_den(qs_env, local_rho_set=p_env%local_rho_set_admm, &
923 274 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
924 274 : rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
925 274 : rho1_atom_set => p_env%local_rho_set_admm%rho_atom_set
926 : CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
927 274 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
928 274 : basis_type = "AUX_FIT_SOFT"
929 274 : task_list => admm_env%admm_gapw_env%task_list
930 : END IF
931 :
932 1182 : alpha = 1.0_dp
933 1182 : IF (nspins == 1) alpha = 2.0_dp
934 :
935 2462 : DO ispin = 1, nspins
936 1280 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
937 1280 : CALL dbcsr_copy(xcmat%matrix, matrix_s(1)%matrix)
938 1280 : CALL dbcsr_set(xcmat%matrix, 0.0_dp)
939 : CALL integrate_v_rspace(v_rspace=v_xc(ispin), hmat=xcmat, qs_env=qs_env, &
940 : calculate_forces=.FALSE., basis_type=basis_type, &
941 1280 : task_list_external=task_list)
942 2462 : CALL dbcsr_add(p_env%kpp1_admm(ispin)%matrix, xcmat%matrix, 1.0_dp, alpha)
943 : END DO
944 :
945 1182 : IF (admm_env%do_gapw) THEN
946 274 : CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
947 274 : ns = SIZE(p_env%kpp1_admm)
948 274 : ksmat(1:ns, 1:1) => p_env%kpp1_admm(1:ns)
949 274 : psmat(1:ns, 1:1) => p_env%p1_admm(1:ns)
950 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE., tddft=.TRUE., &
951 : rho_atom_external=p_env%local_rho_set_admm%rho_atom_set, &
952 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
953 : oce_external=admm_env%admm_gapw_env%oce, &
954 274 : sab_external=sab_aux_fit)
955 : END IF
956 :
957 2462 : DO ispin = 1, nspins
958 2462 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
959 : END DO
960 1182 : DEALLOCATE (v_xc)
961 1182 : CALL dbcsr_deallocate_matrix(xcmat%matrix)
962 :
963 : END IF
964 : END IF
965 :
966 9336 : CALL timestop(handle)
967 :
968 205392 : END SUBROUTINE apply_xc_admm
969 :
970 : END MODULE qs_linres_kernel
|
