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