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