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 Routines to calculate 2nd order kernels from a given response density in ao basis
10 : !> linear response scf
11 : !> \par History
12 : !> created 08-2020 [Frederick Stein], Code by M. Iannuzzi
13 : !> \author Frederick Stein
14 : ! **************************************************************************************************
15 : MODULE qs_2nd_kernel_ao
16 : USE admm_types, ONLY: admm_type,&
17 : get_admm_env
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
20 : dbcsr_copy,&
21 : dbcsr_create,&
22 : dbcsr_p_type,&
23 : dbcsr_release,&
24 : dbcsr_set
25 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
26 : dbcsr_allocate_matrix_set,&
27 : dbcsr_deallocate_matrix_set
28 : USE cp_fm_types, ONLY: cp_fm_get_info,&
29 : cp_fm_type
30 : USE hfx_admm_utils, ONLY: tddft_hfx_matrix
31 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
32 : do_admm_basis_projection,&
33 : do_admm_exch_scaling_none,&
34 : do_admm_purify_none
35 : USE input_section_types, ONLY: section_vals_get,&
36 : section_vals_get_subs_vals,&
37 : section_vals_type
38 : USE kinds, ONLY: dp
39 : USE pw_env_types, ONLY: pw_env_get,&
40 : pw_env_type
41 : USE pw_methods, ONLY: pw_scale
42 : USE pw_pool_types, ONLY: pw_pool_type
43 : USE pw_types, ONLY: pw_c1d_gs_type,&
44 : pw_r3d_rs_type
45 : USE qs_environment_types, ONLY: get_qs_env,&
46 : qs_environment_type
47 : USE qs_integrate_potential, ONLY: integrate_v_rspace
48 : USE qs_kpp1_env_methods, ONLY: calc_kpp1
49 : USE qs_ks_types, ONLY: qs_ks_env_type
50 : USE qs_linres_types, ONLY: linres_control_type
51 : USE qs_p_env_methods, ONLY: p_env_finish_kpp1
52 : USE qs_p_env_types, ONLY: qs_p_env_type
53 : USE qs_rho_types, ONLY: qs_rho_get,&
54 : qs_rho_type
55 : USE task_list_types, ONLY: task_list_type
56 : USE xc, ONLY: xc_calc_2nd_deriv
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : ! *** Public subroutines ***
64 : PUBLIC :: build_dm_response, apply_2nd_order_kernel
65 : PUBLIC :: apply_hfx_ao
66 : PUBLIC :: apply_xc_admm_ao
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_2nd_kernel_ao'
69 :
70 : ! **************************************************************************************************
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief This routine builds response density in dbcsr format
76 : !> \param c0 coefficients of unperturbed system (not changed)
77 : !> \param c1 coefficients of response (not changed)
78 : !> \param dm response density matrix
79 : ! **************************************************************************************************
80 10490 : SUBROUTINE build_dm_response(c0, c1, dm)
81 : !
82 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c0, c1
83 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: dm
84 :
85 : INTEGER :: ispin, ncol, nspins
86 :
87 10490 : nspins = SIZE(dm, 1)
88 :
89 22232 : DO ispin = 1, nspins
90 11742 : CALL dbcsr_set(dm(ispin)%matrix, 0.0_dp)
91 11742 : CALL cp_fm_get_info(c0(ispin), ncol_global=ncol)
92 : CALL cp_dbcsr_plus_fm_fm_t(dm(ispin)%matrix, &
93 : matrix_v=c0(ispin), &
94 : matrix_g=c1(ispin), &
95 : ncol=ncol, alpha=2.0_dp, &
96 : keep_sparsity=.TRUE., &
97 22232 : symmetry_mode=1)
98 : END DO
99 :
100 10490 : END SUBROUTINE build_dm_response
101 :
102 : ! **************************************************************************************************
103 : !> \brief Calculate a second order kernel (DFT, HF, ADMM correction) for a given density
104 : !> \param qs_env ...
105 : !> \param p_env perturbation environment containing the correct density matrices p_env%p1, p_env%p1_admm,
106 : !> the kernel will be saved in p_env%kpp1, p_env%kpp1_admm
107 : !> \param recalc_hfx_integrals whether to recalculate the HFX integrals
108 : !> \param calc_forces whether to calculate forces
109 : !> \param calc_virial whether to calculate virials
110 : !> \param virial collect the virial terms from the XC + ADMM parts (terms from integration will be added to pv_virial)
111 : ! **************************************************************************************************
112 3568 : SUBROUTINE apply_2nd_order_kernel(qs_env, p_env, recalc_hfx_integrals, calc_forces, calc_virial, virial)
113 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
114 : TYPE(qs_p_env_type) :: p_env
115 : LOGICAL, INTENT(IN), OPTIONAL :: recalc_hfx_integrals, calc_forces, &
116 : calc_virial
117 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
118 : OPTIONAL :: virial
119 :
120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_2nd_order_kernel'
121 :
122 : INTEGER :: handle, ispin
123 : LOGICAL :: do_hfx, my_calc_forces, my_calc_virial, &
124 : my_recalc_hfx_integrals
125 : TYPE(admm_type), POINTER :: admm_env
126 : TYPE(dft_control_type), POINTER :: dft_control
127 : TYPE(linres_control_type), POINTER :: linres_control
128 : TYPE(section_vals_type), POINTER :: hfx_sections, input, xc_section
129 :
130 1784 : CALL timeset(routineN, handle)
131 :
132 1784 : my_recalc_hfx_integrals = .FALSE.
133 1784 : IF (PRESENT(recalc_hfx_integrals)) my_recalc_hfx_integrals = recalc_hfx_integrals
134 :
135 1784 : my_calc_forces = .FALSE.
136 1784 : IF (PRESENT(calc_forces)) my_calc_forces = calc_forces
137 :
138 1784 : my_calc_virial = .FALSE.
139 1784 : IF (PRESENT(calc_virial)) my_calc_virial = calc_virial
140 :
141 1784 : CALL get_qs_env(qs_env, dft_control=dft_control)
142 :
143 4118 : DO ispin = 1, SIZE(p_env%kpp1)
144 2334 : CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
145 4118 : IF (dft_control%do_admm) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
146 : END DO
147 :
148 : CALL get_qs_env(qs_env=qs_env, &
149 : input=input, &
150 1784 : linres_control=linres_control)
151 :
152 1784 : IF (dft_control%do_admm) THEN
153 212 : CALL get_qs_env(qs_env, admm_env=admm_env)
154 212 : xc_section => admm_env%xc_section_primary
155 : ELSE
156 1572 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
157 : END IF
158 :
159 : CALL calc_kpp1(p_env%rho1_xc, p_env%rho1, xc_section, &
160 : dft_control%qs_control%lrigpw, linres_control%lr_triplet, &
161 1784 : qs_env, p_env, calc_forces=my_calc_forces, calc_virial=my_calc_virial, virial=virial)
162 :
163 : ! hfx section
164 1784 : NULLIFY (hfx_sections)
165 1784 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
166 1784 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
167 1784 : IF (do_hfx) THEN
168 1286 : CALL apply_hfx_ao(qs_env, p_env, my_recalc_hfx_integrals)
169 :
170 1286 : IF (dft_control%do_admm) THEN
171 188 : CALL apply_xc_admm_ao(qs_env, p_env, my_calc_forces, my_calc_virial, virial)
172 188 : CALL p_env_finish_kpp1(qs_env, p_env)
173 : END IF
174 : END IF
175 :
176 1784 : CALL timestop(handle)
177 :
178 1784 : END SUBROUTINE apply_2nd_order_kernel
179 :
180 : ! **************************************************************************************************
181 : !> \brief This routine applies the Hartree-Fock Exchange kernel to a perturbation density matrix considering ADMM
182 : !> \param qs_env the Quickstep environment
183 : !> \param p_env perturbation environment from which p1/p1_admm and kpp1/kpp1_admm are taken
184 : !> \param recalc_integrals whether the integrals are to be recalculated (default: no)
185 : ! **************************************************************************************************
186 1286 : SUBROUTINE apply_hfx_ao(qs_env, p_env, recalc_integrals)
187 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
188 : TYPE(qs_p_env_type), INTENT(IN) :: p_env
189 : LOGICAL, INTENT(IN), OPTIONAL :: recalc_integrals
190 :
191 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_hfx_ao'
192 :
193 : INTEGER :: handle, ispin, nspins
194 : LOGICAL :: my_recalc_integrals
195 : REAL(KIND=dp) :: alpha
196 1286 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h1_mat, rho1, work_hmat
197 : TYPE(dft_control_type), POINTER :: dft_control
198 :
199 1286 : CALL timeset(routineN, handle)
200 :
201 1286 : my_recalc_integrals = .FALSE.
202 1286 : IF (PRESENT(recalc_integrals)) my_recalc_integrals = recalc_integrals
203 :
204 1286 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
205 :
206 1286 : IF (dft_control%do_admm) THEN
207 188 : IF (dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
208 0 : CPABORT("ADMM: Linear Response needs purification_method=none")
209 : END IF
210 188 : IF (dft_control%admm_control%scaling_model /= do_admm_exch_scaling_none) THEN
211 0 : CPABORT("ADMM: Linear Response needs scaling_model=none")
212 : END IF
213 188 : IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN
214 0 : CPABORT("ADMM: Linear Response needs admm_method=basis_projection")
215 : END IF
216 : !
217 : END IF
218 :
219 1286 : nspins = dft_control%nspins
220 :
221 1286 : IF (dft_control%do_admm) THEN
222 188 : rho1 => p_env%p1_admm
223 188 : h1_mat => p_env%kpp1_admm
224 : ELSE
225 1098 : rho1 => p_env%p1
226 1098 : h1_mat => p_env%kpp1
227 : END IF
228 :
229 2888 : DO ispin = 1, nspins
230 1602 : CPASSERT(ASSOCIATED(rho1(ispin)%matrix))
231 2888 : CPASSERT(ASSOCIATED(h1_mat(ispin)%matrix))
232 : END DO
233 :
234 1286 : NULLIFY (work_hmat)
235 1286 : CALL dbcsr_allocate_matrix_set(work_hmat, nspins)
236 2888 : DO ispin = 1, nspins
237 1602 : ALLOCATE (work_hmat(ispin)%matrix)
238 1602 : CALL dbcsr_create(work_hmat(ispin)%matrix, template=rho1(ispin)%matrix)
239 1602 : CALL dbcsr_copy(work_hmat(ispin)%matrix, rho1(ispin)%matrix)
240 2888 : CALL dbcsr_set(work_hmat(ispin)%matrix, 0.0_dp)
241 : END DO
242 :
243 : ! Calculate kernel
244 1286 : CALL tddft_hfx_matrix(work_hmat, rho1, qs_env, .FALSE., my_recalc_integrals)
245 :
246 1286 : alpha = 2.0_dp
247 1286 : IF (nspins == 2) alpha = 1.0_dp
248 :
249 2888 : DO ispin = 1, nspins
250 2888 : CALL dbcsr_add(h1_mat(ispin)%matrix, work_hmat(ispin)%matrix, 1.0_dp, alpha)
251 : END DO
252 :
253 1286 : CALL dbcsr_deallocate_matrix_set(work_hmat)
254 :
255 1286 : CALL timestop(handle)
256 :
257 1286 : END SUBROUTINE apply_hfx_ao
258 :
259 : ! **************************************************************************************************
260 : !> \brief apply the kernel from the ADMM exchange correction
261 : !> \param qs_env ...
262 : !> \param p_env perturbation environment
263 : !> \param calc_forces whether to calculate forces
264 : !> \param calc_virial whether to calculate gradients
265 : !> \param virial collects the virial terms from the XC functional (virial terms from integration are collected in pv_virial)
266 : ! **************************************************************************************************
267 188 : SUBROUTINE apply_xc_admm_ao(qs_env, p_env, calc_forces, calc_virial, virial)
268 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
269 : TYPE(qs_p_env_type) :: p_env
270 : LOGICAL, INTENT(IN), OPTIONAL :: calc_forces, calc_virial
271 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
272 : OPTIONAL :: virial
273 :
274 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_xc_admm_ao'
275 :
276 : INTEGER :: handle, ispin, nao, nao_aux, nspins
277 : LOGICAL :: lsd, my_calc_forces
278 : REAL(KIND=dp) :: alpha
279 : TYPE(admm_type), POINTER :: admm_env
280 : TYPE(dbcsr_p_type) :: work_hmat
281 188 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux
282 : TYPE(dft_control_type), POINTER :: dft_control
283 : TYPE(linres_control_type), POINTER :: linres_control
284 188 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_aux_g
285 : TYPE(pw_env_type), POINTER :: pw_env
286 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
287 188 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_aux_r, tau1_aux_r, v_xc, v_xc_tau
288 : TYPE(qs_ks_env_type), POINTER :: ks_env
289 : TYPE(qs_rho_type), POINTER :: rho_aux
290 : TYPE(section_vals_type), POINTER :: xc_section
291 : TYPE(task_list_type), POINTER :: task_list_aux_fit
292 :
293 188 : CALL timeset(routineN, handle)
294 :
295 188 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
296 :
297 188 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
298 92 : CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
299 92 : CPASSERT(.NOT. dft_control%qs_control%gapw)
300 92 : CPASSERT(.NOT. dft_control%qs_control%gapw_xc)
301 92 : CPASSERT(.NOT. dft_control%qs_control%lrigpw)
302 92 : CPASSERT(.NOT. linres_control%lr_triplet)
303 92 : IF (.NOT. ASSOCIATED(p_env%kpp1_admm)) &
304 0 : CPABORT("kpp1_admm has to be associated if ADMM kernel calculations are requested")
305 :
306 92 : nspins = dft_control%nspins
307 :
308 92 : my_calc_forces = .FALSE.
309 92 : IF (PRESENT(calc_forces)) my_calc_forces = calc_forces
310 :
311 : ! AUX basis contribution
312 92 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
313 92 : CPASSERT(ASSOCIATED(pw_env))
314 92 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
315 92 : NULLIFY (v_xc)
316 : ! calculate the xc potential
317 92 : lsd = (nspins == 2)
318 92 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, admm_env=admm_env)
319 92 : CALL get_admm_env(admm_env, task_list_aux_fit=task_list_aux_fit)
320 :
321 92 : CALL qs_rho_get(p_env%rho1_admm, rho_r=rho1_aux_r, rho_g=rho1_aux_g, tau_r=tau1_aux_r)
322 92 : xc_section => admm_env%xc_section_aux
323 :
324 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, p_env%kpp1_env%deriv_set_admm, &
325 : p_env%kpp1_env%rho_set_admm, &
326 : rho1_aux_r, rho1_aux_g, tau1_aux_r, auxbas_pw_pool, xc_section=xc_section, gapw=.FALSE., &
327 92 : compute_virial=calc_virial, virial_xc=virial)
328 :
329 : NULLIFY (work_hmat%matrix)
330 92 : ALLOCATE (work_hmat%matrix)
331 92 : CALL dbcsr_copy(work_hmat%matrix, p_env%kpp1_admm(1)%matrix)
332 :
333 92 : alpha = 1.0_dp
334 92 : IF (nspins == 1) alpha = 2.0_dp
335 :
336 92 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
337 92 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
338 :
339 92 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
340 220 : DO ispin = 1, nspins
341 128 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
342 128 : CALL dbcsr_set(work_hmat%matrix, 0.0_dp)
343 : CALL integrate_v_rspace(v_rspace=v_xc(ispin), hmat=work_hmat, qs_env=qs_env, &
344 : calculate_forces=my_calc_forces, basis_type="AUX_FIT", &
345 128 : task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
346 128 : IF (ASSOCIATED(v_xc_tau)) THEN
347 0 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
348 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), hmat=work_hmat, qs_env=qs_env, &
349 : compute_tau=.TRUE., &
350 : calculate_forces=my_calc_forces, basis_type="AUX_FIT", &
351 0 : task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
352 : END IF
353 220 : CALL dbcsr_add(p_env%kpp1_admm(ispin)%matrix, work_hmat%matrix, 1.0_dp, alpha)
354 :
355 : END DO
356 :
357 92 : CALL dbcsr_release(work_hmat%matrix)
358 92 : DEALLOCATE (work_hmat%matrix)
359 :
360 220 : DO ispin = 1, nspins
361 220 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
362 : END DO
363 92 : DEALLOCATE (v_xc)
364 :
365 : END IF
366 :
367 188 : CALL timestop(handle)
368 :
369 188 : END SUBROUTINE apply_xc_admm_ao
370 : END MODULE qs_2nd_kernel_ao
|
