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