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 Routines to calculate EXX in RPA and energy correction methods
10 : !> \par History
11 : !> 07.2020 separated from mp2.F [F. Stein, code by Jan Wilhelm]
12 : !> 06.2022 EXX contribution to the forces [A. Bussy]
13 : !> 03.2023 Generalized for energy correction methods
14 : !> \author Jan Wilhelm, Frederick Stein, Augustin Bussy, Fabian Belleflamme
15 : ! **************************************************************************************************
16 : MODULE hfx_exx
17 : USE admm_methods, ONLY: admm_projection_derivative
18 : USE admm_types, ONLY: admm_env_create,&
19 : admm_env_release,&
20 : admm_type,&
21 : get_admm_env
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
24 : dbcsr_copy,&
25 : dbcsr_create,&
26 : dbcsr_p_type,&
27 : dbcsr_release,&
28 : dbcsr_scale,&
29 : dbcsr_set,&
30 : dbcsr_type
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 : copy_fm_to_dbcsr,&
33 : dbcsr_allocate_matrix_set,&
34 : dbcsr_deallocate_matrix_set
35 : USE cp_log_handling, ONLY: cp_get_default_logger,&
36 : cp_logger_get_default_unit_nr,&
37 : cp_logger_type
38 : USE hfx_admm_utils, ONLY: create_admm_xc_section,&
39 : tddft_hfx_matrix
40 : USE hfx_derivatives, ONLY: derivatives_four_center
41 : USE hfx_energy_potential, ONLY: integrate_four_center
42 : USE hfx_ri, ONLY: hfx_ri_update_forces,&
43 : hfx_ri_update_ks
44 : USE hfx_types, ONLY: hfx_type
45 : USE input_constants, ONLY: do_admm_aux_exch_func_none
46 : USE input_section_types, ONLY: section_vals_create,&
47 : section_vals_duplicate,&
48 : section_vals_get,&
49 : section_vals_get_subs_vals,&
50 : section_vals_release,&
51 : section_vals_set_subs_vals,&
52 : section_vals_type,&
53 : section_vals_val_get
54 : USE kinds, ONLY: dp
55 : USE machine, ONLY: m_walltime
56 : USE message_passing, ONLY: mp_para_env_type
57 : USE parallel_gemm_api, ONLY: parallel_gemm
58 : USE pw_env_types, ONLY: pw_env_get,&
59 : pw_env_type
60 : USE pw_methods, ONLY: pw_scale
61 : USE pw_pool_types, ONLY: pw_pool_type
62 : USE pw_types, ONLY: pw_r3d_rs_type
63 : USE qs_energy_types, ONLY: qs_energy_type
64 : USE qs_environment_types, ONLY: get_qs_env,&
65 : qs_environment_type
66 : USE qs_integrate_potential, ONLY: integrate_v_rspace
67 : USE qs_ks_reference, ONLY: ks_ref_potential
68 : USE qs_rho_types, ONLY: qs_rho_get,&
69 : qs_rho_type
70 : USE qs_vxc, ONLY: qs_vxc_create
71 : USE task_list_types, ONLY: task_list_type
72 : USE virial_types, ONLY: virial_type
73 :
74 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
75 :
76 : #include "./base/base_uses.f90"
77 :
78 : IMPLICIT NONE
79 :
80 : PRIVATE
81 :
82 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_exx'
83 :
84 : PUBLIC :: calculate_exx, add_exx_to_rhs, calc_exx_admm_xc_contributions, exx_pre_hfx, exx_post_hfx
85 :
86 : CONTAINS
87 :
88 : ! **************************************************************************************************
89 : !> \brief ...
90 : !> \param qs_env ...
91 : !> \param unit_nr ...
92 : !> \param hfx_sections ...
93 : !> \param x_data ...
94 : !> \param do_gw ...
95 : !> \param do_admm ...
96 : !> \param calc_forces ...
97 : !> \param reuse_hfx ...
98 : !> \param do_im_time ...
99 : !> \param E_ex_from_GW ...
100 : !> \param E_admm_from_GW ...
101 : !> \param t3 ...
102 : ! **************************************************************************************************
103 436 : SUBROUTINE calculate_exx(qs_env, unit_nr, hfx_sections, x_data, &
104 : do_gw, do_admm, calc_forces, reuse_hfx, do_im_time, &
105 : E_ex_from_GW, E_admm_from_GW, t3)
106 : TYPE(qs_environment_type), POINTER :: qs_env
107 : INTEGER, INTENT(IN) :: unit_nr
108 : TYPE(section_vals_type), POINTER :: hfx_sections
109 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
110 : LOGICAL, INTENT(IN) :: do_gw, do_admm, calc_forces, reuse_hfx, &
111 : do_im_time
112 : REAL(KIND=dp), INTENT(IN) :: E_ex_from_GW, E_admm_from_GW(2), t3
113 :
114 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_exx'
115 :
116 : INTEGER :: handle, i, irep, ispin, mspin, n_rep_hf, &
117 : nspins
118 : LOGICAL :: calc_ints, hfx_treat_lsd_in_core, &
119 : use_virial
120 : REAL(KIND=dp) :: eh1, ehfx, t1, t2, tf1, tf2
121 218 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_ks_aux_fit, rho_ao, &
122 218 : rho_ao_aux_fit, rho_ao_resp
123 218 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_2d, rho_ao_2d
124 : TYPE(dft_control_type), POINTER :: dft_control
125 : TYPE(mp_para_env_type), POINTER :: para_env
126 : TYPE(qs_energy_type), POINTER :: energy
127 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
128 : TYPE(section_vals_type), POINTER :: input
129 : TYPE(virial_type), POINTER :: virial
130 :
131 218 : CALL timeset(routineN, handle)
132 :
133 218 : t1 = m_walltime()
134 :
135 218 : NULLIFY (input, para_env, matrix_ks, matrix_ks_aux_fit, rho, rho_ao, virial, &
136 218 : dft_control, rho_aux_fit, rho_ao_aux_fit)
137 :
138 218 : CALL exx_pre_hfx(hfx_sections, x_data, reuse_hfx)
139 :
140 : CALL get_qs_env(qs_env=qs_env, &
141 : input=input, &
142 : para_env=para_env, &
143 : energy=energy, &
144 : rho=rho, &
145 : matrix_ks=matrix_ks, &
146 : virial=virial, &
147 218 : dft_control=dft_control)
148 218 : CALL qs_rho_get(rho, rho_ao=rho_ao)
149 :
150 218 : IF (do_admm) THEN
151 : CALL get_admm_env(qs_env%admm_env, &
152 : matrix_ks_aux_fit=matrix_ks_aux_fit, &
153 42 : rho_aux_fit=rho_aux_fit)
154 42 : CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
155 :
156 42 : IF (qs_env%admm_env%do_gapw) THEN
157 0 : CPABORT("ADMM EXX only implmented with GPW")
158 : END IF
159 : END IF
160 :
161 218 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
162 : CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
163 218 : i_rep_section=1)
164 :
165 : ! put matrix_ks to zero
166 464 : DO i = 1, SIZE(matrix_ks)
167 246 : CALL dbcsr_set(matrix_ks(i)%matrix, 0.0_dp)
168 464 : IF (do_admm) THEN
169 48 : CALL dbcsr_set(matrix_ks_aux_fit(i)%matrix, 0.0_dp)
170 : END IF
171 : END DO
172 :
173 : ! take the exact exchange energy from GW or calculate it
174 218 : IF (do_gw) THEN
175 :
176 54 : IF (calc_forces) CPABORT("Not implemented")
177 :
178 54 : IF (qs_env%mp2_env%ri_g0w0%update_xc_energy) THEN
179 24 : CALL remove_exc_energy(energy)
180 24 : energy%total = energy%total + E_ex_from_GW
181 24 : energy%ex = E_ex_from_GW
182 24 : IF (do_admm) THEN
183 4 : energy%total = energy%total + E_admm_from_GW(1) + E_admm_from_GW(2)
184 4 : energy%exc = E_admm_from_GW(1)
185 4 : energy%exc_aux_fit = E_admm_from_GW(2)
186 : END IF
187 24 : t2 = m_walltime()
188 :
189 24 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Time=', t2 - t1 + t3
190 24 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX energy = ', energy%ex
191 24 : IF (do_admm .AND. unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') &
192 2 : 'EXX ADMM XC correction = ', E_admm_from_GW(1) + E_admm_from_GW(2)
193 : END IF
194 :
195 : ELSE
196 :
197 164 : CALL remove_exc_energy(energy)
198 :
199 164 : nspins = dft_control%nspins
200 164 : mspin = 1
201 164 : IF (hfx_treat_lsd_in_core) mspin = nspins
202 :
203 164 : calc_ints = .TRUE.
204 164 : IF (reuse_hfx) calc_ints = .FALSE.
205 164 : IF (calc_forces .AND. do_im_time) calc_ints = .FALSE.
206 :
207 164 : ehfx = 0.0_dp
208 164 : IF (do_admm) THEN
209 34 : matrix_ks_2d(1:nspins, 1:1) => matrix_ks_aux_fit(1:nspins)
210 34 : rho_ao_2d(1:nspins, 1:1) => rho_ao_aux_fit(1:nspins)
211 : ELSE
212 130 : matrix_ks_2d(1:nspins, 1:1) => matrix_ks(1:nspins)
213 130 : rho_ao_2d(1:nspins, 1:1) => rho_ao(1:nspins)
214 : END IF
215 :
216 328 : DO irep = 1, n_rep_hf
217 :
218 328 : IF (x_data(irep, 1)%do_hfx_ri) THEN
219 : CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, &
220 : rho_ao=rho_ao_2d, geometry_did_change=calc_ints, nspins=nspins, &
221 8 : hf_fraction=x_data(irep, 1)%general_parameter%fraction)
222 : ELSE
223 :
224 312 : DO ispin = 1, mspin
225 :
226 : CALL integrate_four_center(qs_env, x_data, matrix_ks_2d, eh1, &
227 : rho_ao_2d, hfx_sections, para_env, &
228 156 : calc_ints, irep, .TRUE., ispin=ispin)
229 312 : ehfx = ehfx + eh1
230 : END DO
231 : END IF
232 : END DO
233 :
234 : ! include the EXX contribution to the total energy
235 164 : energy%ex = ehfx
236 164 : energy%total = energy%total + energy%ex
237 :
238 164 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
239 164 : IF (use_virial) THEN
240 0 : virial%pv_calculate = .TRUE.
241 0 : virial%pv_fock_4c = 0.0_dp
242 : END IF
243 :
244 164 : IF (calc_forces) THEN
245 58 : tf1 = m_walltime()
246 :
247 : !Note: no need to remove xc forces: they are not even calculated in the first place
248 58 : NULLIFY (rho_ao_resp)
249 :
250 116 : DO irep = 1, n_rep_hf
251 :
252 116 : IF (x_data(irep, 1)%do_hfx_ri) THEN
253 :
254 : CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
255 : x_data(irep, 1)%general_parameter%fraction, &
256 : rho_ao=rho_ao_2d, rho_ao_resp=rho_ao_resp, &
257 4 : use_virial=use_virial)
258 : ELSE
259 :
260 : CALL derivatives_four_center(qs_env, rho_ao_2d, rho_ao_resp, &
261 : hfx_sections, para_env, irep, &
262 54 : use_virial, external_x_data=x_data)
263 :
264 : END IF
265 :
266 : END DO !irep
267 :
268 58 : tf2 = m_walltime()
269 58 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Force Time=', tf2 - tf1
270 : END IF
271 :
272 164 : IF (use_virial) THEN
273 0 : virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
274 0 : virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
275 : END IF
276 :
277 : ! ADMM XC correction
278 164 : IF (do_admm) THEN
279 :
280 : CALL calc_exx_admm_xc_contributions(qs_env, matrix_ks, matrix_ks_aux_fit, x_data, &
281 : energy%exc, energy%exc_aux_fit, calc_forces, &
282 34 : use_virial)
283 :
284 : ! ADMM overlap forces
285 34 : IF (calc_forces) CALL admm_projection_derivative(qs_env, matrix_ks_aux_fit, rho_ao)
286 :
287 34 : energy%total = energy%total + energy%exc_aux_fit
288 34 : energy%total = energy%total + energy%exc
289 :
290 : END IF
291 :
292 164 : IF (use_virial) THEN
293 0 : virial%pv_calculate = .FALSE.
294 : END IF
295 :
296 164 : t2 = m_walltime()
297 :
298 164 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Time=', t2 - t1 + t3
299 164 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX energy = ', energy%ex
300 164 : IF (do_admm .AND. unit_nr > 0) THEN
301 17 : WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX ADMM XC correction = ', energy%exc + energy%exc_aux_fit
302 : END IF
303 : END IF
304 :
305 218 : CALL exx_post_hfx(qs_env, x_data, reuse_hfx)
306 :
307 218 : CALL timestop(handle)
308 :
309 218 : END SUBROUTINE calculate_exx
310 :
311 : ! **************************************************************************************************
312 : !> \brief Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian
313 : !> \param rhs ...
314 : !> \param qs_env ...
315 : !> \param ext_hfx_section ...
316 : !> \param x_data ...
317 : !> \param recalc_integrals ...
318 : !> \param do_admm ...
319 : !> \param do_ec ...
320 : !> \param do_exx ...
321 : !> \param reuse_hfx ...
322 : ! **************************************************************************************************
323 230 : SUBROUTINE add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, &
324 : recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
325 :
326 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: rhs
327 : TYPE(qs_environment_type), POINTER :: qs_env
328 : TYPE(section_vals_type), POINTER :: ext_hfx_section
329 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
330 : LOGICAL, INTENT(IN), OPTIONAL :: recalc_integrals, do_admm, do_ec, &
331 : do_exx, reuse_hfx
332 :
333 : CHARACTER(LEN=*), PARAMETER :: routineN = 'add_exx_to_rhs'
334 :
335 : INTEGER :: handle, ispin, nao, nao_aux, nspins, &
336 : unit_nr
337 : LOGICAL :: calc_ints, do_hfx, my_do_ec, my_do_exx, &
338 : my_recalc_integrals
339 : REAL(dp) :: dummy_real1, dummy_real2
340 : TYPE(admm_type), POINTER :: admm_env
341 : TYPE(cp_logger_type), POINTER :: logger
342 230 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_work, matrix_ks, matrix_s_aux, &
343 230 : rho_ao, rho_ao_aux, work_admm
344 : TYPE(dbcsr_type) :: dbcsr_tmp
345 : TYPE(dft_control_type), POINTER :: dft_control
346 : TYPE(pw_env_type), POINTER :: pw_env
347 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
348 : TYPE(pw_r3d_rs_type) :: vh_rspace
349 230 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vadmm_rspace, vtau_rspace, vxc_rspace
350 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
351 : TYPE(section_vals_type), POINTER :: hfx_section
352 : TYPE(task_list_type), POINTER :: task_list_aux_fit
353 :
354 230 : NULLIFY (dbcsr_work, matrix_ks, rho, hfx_section, pw_env, vadmm_rspace, vtau_rspace, vxc_rspace, &
355 230 : auxbas_pw_pool, dft_control, rho_ao, rho_aux_fit, rho_ao_aux, work_admm, &
356 230 : matrix_s_aux, admm_env, task_list_aux_fit)
357 :
358 460 : logger => cp_get_default_logger()
359 230 : IF (logger%para_env%is_source()) THEN
360 124 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
361 : ELSE
362 : unit_nr = -1
363 : END IF
364 :
365 230 : CALL timeset(routineN, handle)
366 :
367 230 : my_recalc_integrals = .FALSE.
368 230 : IF (PRESENT(recalc_integrals)) my_recalc_integrals = recalc_integrals
369 230 : my_do_ec = .FALSE.
370 230 : IF (PRESENT(do_ec)) my_do_ec = do_ec
371 230 : my_do_exx = .TRUE.
372 230 : IF (PRESENT(do_exx)) my_do_exx = do_exx
373 :
374 : ! Strategy: we take the ks_matrix, remove the current xc contribution, and then add the RPA HF one
375 230 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks, rho=rho, pw_env=pw_env, dft_control=dft_control)
376 230 : nspins = dft_control%nspins
377 :
378 : ! do_exx ; subtract XC and EX
379 : ! do_dcdft: subtract EX
380 :
381 230 : CALL dbcsr_allocate_matrix_set(dbcsr_work, nspins)
382 468 : DO ispin = 1, nspins
383 238 : ALLOCATE (dbcsr_work(ispin)%matrix)
384 238 : CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix)
385 468 : CALL dbcsr_set(dbcsr_work(ispin)%matrix, 0.0_dp)
386 : END DO
387 :
388 230 : IF (dft_control%do_admm) THEN
389 80 : CALL get_qs_env(qs_env, admm_env=admm_env)
390 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux, task_list_aux_fit=task_list_aux_fit, &
391 80 : rho_aux_fit=rho_aux_fit)
392 80 : CALL dbcsr_allocate_matrix_set(work_admm, nspins)
393 164 : DO ispin = 1, nspins
394 84 : ALLOCATE (work_admm(ispin)%matrix)
395 84 : CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
396 84 : CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
397 164 : CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
398 : END DO
399 80 : CALL dbcsr_create(dbcsr_tmp, template=matrix_ks(1)%matrix)
400 80 : CALL dbcsr_copy(dbcsr_tmp, matrix_ks(1)%matrix)
401 :
402 80 : nao = admm_env%nao_orb
403 80 : nao_aux = admm_env%nao_aux_fit
404 80 : CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
405 : END IF
406 :
407 230 : CALL qs_rho_get(rho, rho_ao=rho_ao)
408 :
409 : ! Remove the standard XC + HFX contribution
410 230 : CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, dummy_real1, dummy_real2)
411 :
412 : ! Only remove standard XC and/or HFX
413 : ! (due to different requirements for RHS)
414 230 : IF (my_do_exx) THEN
415 :
416 60 : DO ispin = 1, nspins
417 60 : CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix)
418 : END DO
419 :
420 60 : DO ispin = 1, nspins
421 34 : CALL pw_scale(vxc_rspace(ispin), -1.0_dp)
422 : CALL integrate_v_rspace(v_rspace=vxc_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, &
423 34 : calculate_forces=.FALSE.)
424 60 : IF (ASSOCIATED(vtau_rspace)) THEN
425 0 : CALL pw_scale(vtau_rspace(ispin), -1.0_dp)
426 : CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, &
427 0 : calculate_forces=.FALSE., compute_tau=.TRUE.)
428 : END IF
429 : END DO
430 :
431 : END IF ! do_exx
432 :
433 : ! Remove standard HFX
434 230 : IF (my_do_exx .OR. my_do_ec) THEN
435 :
436 198 : hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
437 198 : CALL section_vals_get(hfx_section, explicit=do_hfx)
438 198 : IF (do_hfx) THEN
439 100 : IF (dft_control%do_admm) THEN
440 :
441 : !note: factor -1.0 taken care of later, after HFX ADMM contribution is taken
442 56 : IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
443 76 : DO ispin = 1, nspins
444 : CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=work_admm(ispin), &
445 : qs_env=qs_env, calculate_forces=.FALSE., basis_type="AUX_FIT", &
446 76 : task_list_external=task_list_aux_fit)
447 : END DO
448 : END IF
449 :
450 56 : CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., my_recalc_integrals)
451 :
452 116 : DO ispin = 1, nspins
453 60 : CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
454 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
455 60 : 0.0_dp, admm_env%work_aux_orb)
456 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
457 60 : 0.0_dp, admm_env%work_orb_orb)
458 60 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_tmp, keep_sparsity=.TRUE.)
459 116 : CALL dbcsr_add(dbcsr_work(ispin)%matrix, dbcsr_tmp, 1.0_dp, -1.0_dp)
460 : END DO
461 : ELSE
462 88 : DO ispin = 1, nspins
463 88 : CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp)
464 : END DO
465 44 : CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., my_recalc_integrals)
466 88 : DO ispin = 1, nspins
467 88 : CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp)
468 : END DO
469 : END IF
470 : END IF !do_hfx
471 :
472 : END IF ! do_exx
473 :
474 : ! Clean
475 230 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
476 230 : CALL auxbas_pw_pool%give_back_pw(vh_rspace)
477 468 : DO ispin = 1, nspins
478 238 : CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
479 468 : IF (ASSOCIATED(vtau_rspace)) THEN
480 16 : CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
481 : END IF
482 : END DO
483 230 : DEALLOCATE (vxc_rspace)
484 230 : IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
485 :
486 230 : IF (dft_control%do_admm) THEN
487 164 : DO ispin = 1, nspins
488 164 : IF (ASSOCIATED(vadmm_rspace)) THEN
489 56 : CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
490 : END IF
491 : END DO
492 80 : IF (ASSOCIATED(vadmm_rspace)) DEALLOCATE (vadmm_rspace)
493 : END IF
494 :
495 : !Add the HF contribution from RI_RPA/EC_ENV to the ks matrix
496 230 : CALL exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx)
497 :
498 230 : calc_ints = .TRUE.
499 230 : IF (reuse_hfx) calc_ints = .FALSE.
500 230 : IF (do_admm) THEN
501 68 : DO ispin = 1, nspins
502 68 : CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
503 : END DO
504 :
505 : CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., calc_ints, ext_hfx_section, &
506 32 : x_data)
507 :
508 : !ADMM XC correction
509 : CALL calc_exx_admm_xc_contributions(qs_env, dbcsr_work, work_admm, x_data, dummy_real1, &
510 32 : dummy_real2, .FALSE., .FALSE.)
511 :
512 68 : DO ispin = 1, nspins
513 36 : CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
514 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
515 36 : 0.0_dp, admm_env%work_aux_orb)
516 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
517 36 : 0.0_dp, admm_env%work_orb_orb)
518 36 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_tmp, keep_sparsity=.TRUE.)
519 68 : CALL dbcsr_add(dbcsr_work(ispin)%matrix, dbcsr_tmp, 1.0_dp, 1.0_dp)
520 : END DO
521 :
522 : ELSE
523 198 : CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., calc_ints, ext_hfx_section, x_data)
524 : END IF
525 230 : CALL exx_post_hfx(qs_env, x_data, reuse_hfx)
526 :
527 : !Update the RHS
528 468 : DO ispin = 1, nspins
529 468 : CALL dbcsr_add(rhs(ispin)%matrix, dbcsr_work(ispin)%matrix, 1.0_dp, 1.0_dp)
530 : END DO
531 :
532 230 : CALL dbcsr_deallocate_matrix_set(dbcsr_work)
533 230 : IF (dft_control%do_admm) THEN
534 80 : CALL dbcsr_deallocate_matrix_set(work_admm)
535 80 : CALL dbcsr_release(dbcsr_tmp)
536 : END IF
537 :
538 230 : CALL timestop(handle)
539 :
540 230 : END SUBROUTINE add_exx_to_rhs
541 :
542 : ! **************************************************************************************************
543 : !> \brief get the ADMM XC section from the ri_rpa/ec_env type if available, create and store them otherwise
544 : !> \param qs_env ...
545 : !> \param x_data ...
546 : !> \param xc_section_aux ...
547 : !> \param xc_section_primary ...
548 : ! **************************************************************************************************
549 74 : SUBROUTINE get_exx_admm_xc_sections(qs_env, x_data, xc_section_aux, xc_section_primary)
550 : TYPE(qs_environment_type), POINTER :: qs_env
551 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
552 : TYPE(section_vals_type), POINTER :: xc_section_aux, xc_section_primary
553 :
554 : INTEGER :: natom
555 : TYPE(admm_type), POINTER :: qs_admm_env, tmp_admm_env
556 : TYPE(dft_control_type), POINTER :: dft_control
557 : TYPE(mp_para_env_type), POINTER :: para_env
558 : TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_empty, xc_section, &
559 : xc_section_empty
560 :
561 74 : NULLIFY (qs_admm_env, tmp_admm_env, para_env, xc_section, xc_section_empty, xc_fun_empty, &
562 74 : xc_fun, dft_control)
563 :
564 74 : IF (ASSOCIATED(qs_env%mp2_env)) THEN
565 26 : IF (ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_aux) .AND. &
566 : ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_primary)) THEN
567 12 : xc_section_aux => qs_env%mp2_env%ri_rpa%xc_section_aux
568 12 : xc_section_primary => qs_env%mp2_env%ri_rpa%xc_section_primary
569 : END IF
570 48 : ELSEIF (qs_env%energy_correction) THEN
571 48 : IF (ASSOCIATED(qs_env%ec_env%xc_section_aux) .AND. &
572 : ASSOCIATED(qs_env%ec_env%xc_section_primary)) THEN
573 42 : xc_section_aux => qs_env%ec_env%xc_section_aux
574 42 : xc_section_primary => qs_env%ec_env%xc_section_primary
575 : END IF
576 : END IF
577 :
578 74 : IF (.NOT. ASSOCIATED(xc_section_aux) .OR. .NOT. ASSOCIATED(xc_section_primary)) THEN
579 :
580 20 : CALL get_qs_env(qs_env, admm_env=qs_admm_env, natom=natom, para_env=para_env, dft_control=dft_control)
581 20 : CPASSERT(ASSOCIATED(qs_admm_env))
582 :
583 : !create XC section with XC_FUNCITONAL NONE (aka empty XC_FUNCTIONAL section)
584 20 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
585 20 : xc_fun => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
586 20 : CALL section_vals_duplicate(xc_section, xc_section_empty)
587 20 : CALL section_vals_create(xc_fun_empty, xc_fun%section)
588 20 : CALL section_vals_set_subs_vals(xc_section_empty, "XC_FUNCTIONAL", xc_fun_empty)
589 :
590 : CALL admm_env_create(tmp_admm_env, dft_control%admm_control, qs_admm_env%mos_aux_fit, &
591 20 : para_env, natom, qs_admm_env%nao_aux_fit)
592 :
593 : CALL create_admm_xc_section(x_data=x_data, xc_section=xc_section_empty, &
594 20 : admm_env=tmp_admm_env)
595 :
596 20 : CALL section_vals_duplicate(tmp_admm_env%xc_section_aux, xc_section_aux)
597 20 : CALL section_vals_duplicate(tmp_admm_env%xc_section_primary, xc_section_primary)
598 :
599 20 : IF (ASSOCIATED(qs_env%mp2_env)) THEN
600 14 : qs_env%mp2_env%ri_rpa%xc_section_aux => xc_section_aux
601 14 : qs_env%mp2_env%ri_rpa%xc_section_primary => xc_section_primary
602 : END IF
603 20 : IF (qs_env%energy_correction) THEN
604 6 : qs_env%ec_env%xc_section_aux => xc_section_aux
605 6 : qs_env%ec_env%xc_section_primary => xc_section_primary
606 : END IF
607 :
608 20 : CALL section_vals_release(xc_section_empty)
609 20 : CALL section_vals_release(xc_fun_empty)
610 20 : CALL admm_env_release(tmp_admm_env)
611 :
612 : END IF
613 :
614 74 : END SUBROUTINE get_exx_admm_xc_sections
615 :
616 : ! **************************************************************************************************
617 : !> \brief Calculate the RI_RPA%HF / EC_ENV%HF ADMM XC contributions to the KS matrices and the respective energies
618 : !> \param qs_env ...
619 : !> \param matrix_prim ...
620 : !> \param matrix_aux ...
621 : !> \param x_data ...
622 : !> \param exc ...
623 : !> \param exc_aux_fit ...
624 : !> \param calc_forces ...
625 : !> \param use_virial ...
626 : ! **************************************************************************************************
627 222 : SUBROUTINE calc_exx_admm_xc_contributions(qs_env, matrix_prim, matrix_aux, x_data, exc, exc_aux_fit, &
628 : calc_forces, use_virial)
629 : TYPE(qs_environment_type), POINTER :: qs_env
630 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_prim, matrix_aux
631 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
632 : REAL(dp), INTENT(INOUT) :: exc, exc_aux_fit
633 : LOGICAL, INTENT(IN) :: calc_forces, use_virial
634 :
635 : INTEGER :: ispin, nspins
636 : REAL(dp), DIMENSION(3, 3) :: pv_loc
637 : TYPE(admm_type), POINTER :: admm_env
638 74 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao, rho_ao_aux_fit
639 : TYPE(dft_control_type), POINTER :: dft_control
640 : TYPE(pw_env_type), POINTER :: pw_env
641 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
642 74 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_dummy, v_rspace, v_rspace_aux_fit
643 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
644 : TYPE(section_vals_type), POINTER :: xc_section_aux, xc_section_primary
645 : TYPE(virial_type), POINTER :: virial
646 :
647 74 : NULLIFY (xc_section_aux, xc_section_primary, rho, rho_aux_fit, v_dummy, v_rspace, v_rspace_aux_fit, &
648 74 : auxbas_pw_pool, pw_env, rho_ao, rho_ao_aux_fit, dft_control, admm_env)
649 :
650 74 : CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, rho=rho, admm_env=admm_env, virial=virial)
651 74 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
652 :
653 74 : nspins = dft_control%nspins
654 74 : CALL qs_rho_get(rho, rho_ao=rho_ao)
655 74 : CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
656 74 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
657 :
658 74 : CALL get_exx_admm_xc_sections(qs_env, x_data, xc_section_aux, xc_section_primary)
659 :
660 74 : IF (use_virial) virial%pv_xc = 0.0_dp
661 : CALL qs_vxc_create(qs_env%ks_env, rho_struct=rho_aux_fit, xc_section=xc_section_aux, &
662 74 : vxc_rho=v_rspace_aux_fit, vxc_tau=v_dummy, exc=exc_aux_fit)
663 74 : IF (use_virial) THEN
664 0 : virial%pv_exc = virial%pv_exc - virial%pv_xc
665 0 : virial%pv_virial = virial%pv_virial - virial%pv_xc
666 : END IF
667 :
668 74 : IF (.NOT. dft_control%admm_control%aux_exch_func == do_admm_aux_exch_func_none) THEN
669 74 : IF (use_virial) pv_loc = virial%pv_virial
670 158 : DO ispin = 1, nspins
671 84 : CALL pw_scale(v_rspace_aux_fit(ispin), v_rspace_aux_fit(ispin)%pw_grid%dvol)
672 : CALL integrate_v_rspace(v_rspace=v_rspace_aux_fit(ispin), hmat=matrix_aux(ispin), &
673 : pmat=rho_ao_aux_fit(ispin), qs_env=qs_env, &
674 : basis_type="AUX_FIT", calculate_forces=calc_forces, &
675 158 : task_list_external=qs_env%admm_env%task_list_aux_fit)
676 : END DO
677 74 : IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
678 : END IF
679 :
680 74 : IF (ASSOCIATED(v_rspace_aux_fit)) THEN
681 158 : DO ispin = 1, nspins
682 158 : CALL auxbas_pw_pool%give_back_pw(v_rspace_aux_fit(ispin))
683 : END DO
684 74 : DEALLOCATE (v_rspace_aux_fit)
685 : END IF
686 74 : IF (ASSOCIATED(v_dummy)) THEN
687 0 : DO ispin = 1, nspins
688 0 : CALL auxbas_pw_pool%give_back_pw(v_dummy(ispin))
689 : END DO
690 0 : DEALLOCATE (v_dummy)
691 : END IF
692 :
693 74 : IF (use_virial) virial%pv_xc = 0.0_dp
694 : CALL qs_vxc_create(qs_env%ks_env, rho_struct=rho, xc_section=xc_section_primary, &
695 74 : vxc_rho=v_rspace, vxc_tau=v_dummy, exc=exc)
696 74 : IF (use_virial) THEN
697 0 : virial%pv_exc = virial%pv_exc - virial%pv_xc
698 0 : virial%pv_virial = virial%pv_virial - virial%pv_xc
699 : END IF
700 :
701 74 : IF (.NOT. dft_control%admm_control%aux_exch_func == do_admm_aux_exch_func_none) THEN
702 74 : IF (use_virial) pv_loc = virial%pv_virial
703 158 : DO ispin = 1, nspins
704 84 : CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
705 : CALL integrate_v_rspace(v_rspace=v_rspace(ispin), hmat=matrix_prim(ispin), &
706 : pmat=rho_ao(ispin), qs_env=qs_env, &
707 158 : calculate_forces=calc_forces)
708 : END DO
709 74 : IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
710 : END IF
711 :
712 74 : IF (ASSOCIATED(v_rspace)) THEN
713 158 : DO ispin = 1, nspins
714 158 : CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
715 : END DO
716 74 : DEALLOCATE (v_rspace)
717 : END IF
718 74 : IF (ASSOCIATED(v_dummy)) THEN
719 0 : DO ispin = 1, nspins
720 0 : CALL auxbas_pw_pool%give_back_pw(v_dummy(ispin))
721 : END DO
722 0 : DEALLOCATE (v_dummy)
723 : END IF
724 :
725 74 : END SUBROUTINE calc_exx_admm_xc_contributions
726 :
727 : ! **************************************************************************************************
728 : !> \brief Prepare the external x_data for integration. Simply change the HFX fraction in case the
729 : !> qs_env%x_data is reused
730 : !> \param ext_hfx_section ...
731 : !> \param x_data ...
732 : !> \param reuse_hfx ...
733 : ! **************************************************************************************************
734 612 : SUBROUTINE exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx)
735 : TYPE(section_vals_type), POINTER :: ext_hfx_section
736 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
737 : LOGICAL :: reuse_hfx
738 :
739 : INTEGER :: irep, n_rep_hf
740 : REAL(dp) :: frac
741 :
742 498 : IF (.NOT. reuse_hfx) RETURN
743 :
744 114 : CALL section_vals_get(ext_hfx_section, n_repetition=n_rep_hf)
745 :
746 228 : DO irep = 1, n_rep_hf
747 114 : CALL section_vals_val_get(ext_hfx_section, "FRACTION", r_val=frac, i_rep_section=irep)
748 342 : x_data(irep, :)%general_parameter%fraction = frac
749 : END DO
750 :
751 : END SUBROUTINE exx_pre_hfx
752 :
753 : ! **************************************************************************************************
754 : !> \brief Revert back to the proper HFX fraction in case qs_env%x_data is reused
755 : !> \param qs_env ...
756 : !> \param x_data ...
757 : !> \param reuse_hfx ...
758 : ! **************************************************************************************************
759 612 : SUBROUTINE exx_post_hfx(qs_env, x_data, reuse_hfx)
760 : TYPE(qs_environment_type), POINTER :: qs_env
761 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
762 : LOGICAL :: reuse_hfx
763 :
764 : INTEGER :: irep, n_rep_hf
765 : REAL(dp) :: frac
766 : TYPE(section_vals_type), POINTER :: input, qs_hfx_section
767 :
768 498 : IF (.NOT. reuse_hfx) RETURN
769 :
770 114 : CALL get_qs_env(qs_env, input=input)
771 114 : qs_hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
772 114 : CALL section_vals_get(qs_hfx_section, n_repetition=n_rep_hf)
773 :
774 228 : DO irep = 1, n_rep_hf
775 114 : CALL section_vals_val_get(qs_hfx_section, "FRACTION", r_val=frac, i_rep_section=irep)
776 342 : x_data(irep, :)%general_parameter%fraction = frac
777 : END DO
778 :
779 : END SUBROUTINE exx_post_hfx
780 :
781 : ! **************************************************************************************************
782 : !> \brief ...
783 : !> \param energy ...
784 : ! **************************************************************************************************
785 188 : ELEMENTAL SUBROUTINE remove_exc_energy(energy)
786 : TYPE(qs_energy_type), INTENT(INOUT) :: energy
787 :
788 : ! Remove the Exchange-correlation energy contributions from the total energy
789 : energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
790 188 : energy%exc_aux_fit + energy%exc1_aux_fit)
791 :
792 188 : energy%exc = 0.0_dp
793 188 : energy%exc1 = 0.0_dp
794 188 : energy%exc_aux_fit = 0.0_dp
795 188 : energy%exc1_aux_fit = 0.0_dp
796 188 : energy%ex = 0.0_dp
797 :
798 188 : END SUBROUTINE
799 :
800 : END MODULE hfx_exx
801 :
|
