Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_r3d_rs_type
44 : USE qs_environment_types, ONLY: get_qs_env,&
45 : qs_environment_type
46 : USE qs_fxc, ONLY: qs_fxc_create
47 : USE qs_integrate_potential, ONLY: integrate_v_rspace
48 : USE qs_kpp1_env_methods, ONLY: calc_kpp1
49 : USE qs_linres_types, ONLY: linres_control_type
50 : USE qs_p_env_methods, ONLY: p_env_finish_kpp1
51 : USE qs_p_env_types, ONLY: qs_p_env_type
52 : USE qs_rho_atom_types, ONLY: rho_atom_type
53 : USE qs_rho_types, ONLY: qs_rho_get,&
54 : qs_rho_type
55 : USE task_list_types, ONLY: task_list_type
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : ! *** Public subroutines ***
63 : PUBLIC :: build_dm_response
64 : PUBLIC :: 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 12380 : 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 12380 : nspins = SIZE(dm, 1)
88 :
89 26198 : DO ispin = 1, nspins
90 13818 : CALL dbcsr_set(dm(ispin)%matrix, 0.0_dp)
91 13818 : 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 26198 : symmetry_mode=1)
98 : END DO
99 :
100 12380 : 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, qs_env, p_env, &
161 1784 : 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 :: 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 : TYPE(pw_env_type), POINTER :: pw_env
285 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
286 188 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
287 : TYPE(qs_rho_type), POINTER :: rho_aux
288 188 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
289 : TYPE(section_vals_type), POINTER :: xc_section
290 : TYPE(task_list_type), POINTER :: task_list_aux_fit
291 :
292 188 : CALL timeset(routineN, handle)
293 :
294 188 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
295 :
296 188 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
297 92 : CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
298 92 : CPASSERT(.NOT. dft_control%qs_control%gapw)
299 92 : CPASSERT(.NOT. dft_control%qs_control%gapw_xc)
300 92 : CPASSERT(.NOT. dft_control%qs_control%lrigpw)
301 92 : CPASSERT(.NOT. linres_control%lr_triplet)
302 92 : IF (.NOT. ASSOCIATED(p_env%kpp1_admm)) THEN
303 0 : CPABORT("kpp1_admm has to be associated if ADMM kernel calculations are requested")
304 : END IF
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 : ! calculate the xc potential
313 92 : NULLIFY (v_xc, v_xc_tau)
314 92 : CALL get_qs_env(qs_env=qs_env, admm_env=admm_env)
315 92 : xc_section => admm_env%xc_section_aux
316 92 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux)
317 : CALL qs_fxc_create(qs_env, rho_aux, p_env%rho1_admm, rho0_atom_set, xc_section, &
318 : .FALSE., v_xc, v_xc_tau, rho1_atom_set, &
319 92 : compute_virial=calc_virial, virial_xc=virial)
320 :
321 : NULLIFY (work_hmat%matrix)
322 92 : ALLOCATE (work_hmat%matrix)
323 92 : CALL dbcsr_copy(work_hmat%matrix, p_env%kpp1_admm(1)%matrix)
324 :
325 92 : alpha = 1.0_dp
326 92 : IF (nspins == 1) alpha = 2.0_dp
327 :
328 92 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux)
329 92 : CALL get_admm_env(admm_env, task_list_aux_fit=task_list_aux_fit)
330 :
331 92 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux, ncol_global=nao)
332 220 : DO ispin = 1, nspins
333 128 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
334 128 : CALL dbcsr_set(work_hmat%matrix, 0.0_dp)
335 : CALL integrate_v_rspace(v_rspace=v_xc(ispin), hmat=work_hmat, qs_env=qs_env, &
336 : calculate_forces=my_calc_forces, basis_type="AUX_FIT", &
337 128 : task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
338 128 : IF (ASSOCIATED(v_xc_tau)) THEN
339 0 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
340 : CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), hmat=work_hmat, qs_env=qs_env, &
341 : compute_tau=.TRUE., &
342 : calculate_forces=my_calc_forces, basis_type="AUX_FIT", &
343 0 : task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
344 : END IF
345 220 : CALL dbcsr_add(p_env%kpp1_admm(ispin)%matrix, work_hmat%matrix, 1.0_dp, alpha)
346 :
347 : END DO
348 :
349 92 : CALL dbcsr_release(work_hmat%matrix)
350 92 : DEALLOCATE (work_hmat%matrix)
351 :
352 92 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
353 92 : CPASSERT(ASSOCIATED(pw_env))
354 92 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
355 :
356 220 : DO ispin = 1, nspins
357 220 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
358 : END DO
359 92 : DEALLOCATE (v_xc)
360 92 : IF (ASSOCIATED(v_xc_tau)) THEN
361 0 : DO ispin = 1, nspins
362 0 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
363 : END DO
364 0 : DEALLOCATE (v_xc_tau)
365 : END IF
366 :
367 : END IF
368 :
369 188 : CALL timestop(handle)
370 :
371 188 : END SUBROUTINE apply_xc_admm_ao
372 :
373 : END MODULE qs_2nd_kernel_ao
|
