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 Utilities for hfx and admm methods
10 : !>
11 : !>
12 : !> \par History
13 : !> refactoring 03-2011 [MI]
14 : !> Made GAPW compatible 12.2019 (A. Bussy)
15 : !> \author MI
16 : ! **************************************************************************************************
17 : MODULE hfx_admm_utils
18 : USE admm_dm_types, ONLY: admm_dm_create
19 : USE admm_methods, ONLY: kpoint_calc_admm_matrices,&
20 : scale_dm
21 : USE admm_types, ONLY: admm_env_create,&
22 : admm_gapw_r3d_rs_type,&
23 : admm_type,&
24 : get_admm_env,&
25 : set_admm_env
26 : USE atomic_kind_types, ONLY: atomic_kind_type
27 : USE basis_set_container_types, ONLY: add_basis_set_to_container
28 : USE basis_set_types, ONLY: copy_gto_basis_set,&
29 : get_gto_basis_set,&
30 : gto_basis_set_type
31 : USE cell_types, ONLY: cell_type,&
32 : plane_distance
33 : USE cp_blacs_env, ONLY: cp_blacs_env_type
34 : USE cp_control_types, ONLY: admm_control_type,&
35 : dft_control_type
36 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
37 : dbcsr_copy,&
38 : dbcsr_create,&
39 : dbcsr_init_p,&
40 : dbcsr_p_type,&
41 : dbcsr_set,&
42 : dbcsr_type,&
43 : dbcsr_type_no_symmetry
44 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
45 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_m_by_n_from_row_template,&
46 : dbcsr_allocate_matrix_set
47 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type
48 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
49 : cp_fm_struct_release,&
50 : cp_fm_struct_type
51 : USE cp_fm_types, ONLY: cp_fm_create,&
52 : cp_fm_get_info,&
53 : cp_fm_type
54 : USE cp_log_handling, ONLY: cp_get_default_logger,&
55 : cp_logger_get_default_io_unit,&
56 : cp_logger_type,&
57 : cp_to_string
58 : USE distribution_1d_types, ONLY: distribution_1d_type
59 : USE distribution_2d_types, ONLY: distribution_2d_type
60 : USE external_potential_types, ONLY: copy_potential
61 : USE hfx_derivatives, ONLY: derivatives_four_center
62 : USE hfx_energy_potential, ONLY: integrate_four_center
63 : USE hfx_pw_methods, ONLY: pw_hfx
64 : USE hfx_ri, ONLY: hfx_ri_update_forces,&
65 : hfx_ri_update_ks
66 : USE hfx_ri_kp, ONLY: hfx_ri_update_forces_kp,&
67 : hfx_ri_update_ks_kp
68 : USE hfx_types, ONLY: hfx_type
69 : USE input_constants, ONLY: &
70 : do_admm_aux_exch_func_bee, do_admm_aux_exch_func_bee_libxc, do_admm_aux_exch_func_default, &
71 : do_admm_aux_exch_func_default_libxc, do_admm_aux_exch_func_none, &
72 : do_admm_aux_exch_func_opt, do_admm_aux_exch_func_opt_libxc, do_admm_aux_exch_func_pbex, &
73 : do_admm_aux_exch_func_pbex_libxc, do_admm_aux_exch_func_sx_libxc, &
74 : do_admm_basis_projection, do_admm_charge_constrained_projection, do_admm_purify_none, &
75 : do_potential_coulomb, do_potential_id, do_potential_long, do_potential_mix_cl, &
76 : do_potential_mix_cl_trunc, do_potential_short, do_potential_truncated, &
77 : xc_funct_no_shortcut, xc_none
78 : USE input_section_types, ONLY: section_vals_duplicate,&
79 : section_vals_get,&
80 : section_vals_get_subs_vals,&
81 : section_vals_get_subs_vals2,&
82 : section_vals_remove_values,&
83 : section_vals_type,&
84 : section_vals_val_get,&
85 : section_vals_val_set
86 : USE kinds, ONLY: dp
87 : USE kpoint_methods, ONLY: kpoint_initialize_mos
88 : USE kpoint_transitional, ONLY: kpoint_transitional_release,&
89 : set_2d_pointer
90 : USE kpoint_types, ONLY: get_kpoint_info,&
91 : kpoint_type
92 : USE libint_2c_3c, ONLY: cutoff_screen_factor
93 : USE mathlib, ONLY: erfc_cutoff
94 : USE message_passing, ONLY: mp_para_env_type
95 : USE molecule_types, ONLY: molecule_type
96 : USE particle_types, ONLY: particle_type
97 : USE paw_proj_set_types, ONLY: get_paw_proj_set,&
98 : paw_proj_set_type
99 : USE pw_env_types, ONLY: pw_env_get,&
100 : pw_env_type
101 : USE pw_poisson_types, ONLY: pw_poisson_type
102 : USE pw_pool_types, ONLY: pw_pool_type
103 : USE pw_types, ONLY: pw_r3d_rs_type
104 : USE qs_energy_types, ONLY: qs_energy_type
105 : USE qs_environment_types, ONLY: get_qs_env,&
106 : qs_environment_type,&
107 : set_qs_env
108 : USE qs_interactions, ONLY: init_interaction_radii
109 : USE qs_kind_types, ONLY: get_qs_kind,&
110 : get_qs_kind_set,&
111 : init_gapw_basis_set,&
112 : init_gapw_nlcc,&
113 : qs_kind_type
114 : USE qs_ks_types, ONLY: qs_ks_env_type
115 : USE qs_local_rho_types, ONLY: local_rho_set_create
116 : USE qs_matrix_pools, ONLY: mpools_get
117 : USE qs_mo_types, ONLY: allocate_mo_set,&
118 : get_mo_set,&
119 : init_mo_set,&
120 : mo_set_type
121 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
122 : release_neighbor_list_sets
123 : USE qs_neighbor_lists, ONLY: atom2d_build,&
124 : atom2d_cleanup,&
125 : build_neighbor_lists,&
126 : local_atoms_type,&
127 : pair_radius_setup,&
128 : write_neighbor_lists
129 : USE qs_oce_methods, ONLY: build_oce_matrices
130 : USE qs_oce_types, ONLY: allocate_oce_set,&
131 : create_oce_set
132 : USE qs_overlap, ONLY: build_overlap_matrix
133 : USE qs_rho_atom_methods, ONLY: init_rho_atom
134 : USE qs_rho_methods, ONLY: qs_rho_rebuild
135 : USE qs_rho_types, ONLY: qs_rho_create,&
136 : qs_rho_get,&
137 : qs_rho_type
138 : USE rt_propagation_types, ONLY: rt_prop_type
139 : USE task_list_methods, ONLY: generate_qs_task_list
140 : USE task_list_types, ONLY: allocate_task_list,&
141 : deallocate_task_list
142 : USE virial_types, ONLY: virial_type
143 : USE xc_adiabatic_utils, ONLY: rescale_xc_potential
144 : #include "./base/base_uses.f90"
145 :
146 : IMPLICIT NONE
147 :
148 : PRIVATE
149 :
150 : ! *** Public subroutines ***
151 : PUBLIC :: hfx_ks_matrix, hfx_admm_init, aux_admm_init, create_admm_xc_section, &
152 : tddft_hfx_matrix, hfx_ks_matrix_kp
153 :
154 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_admm_utils'
155 :
156 : CONTAINS
157 :
158 : ! **************************************************************************************************
159 : !> \brief ...
160 : !> \param qs_env ...
161 : !> \param calculate_forces ...
162 : ! **************************************************************************************************
163 22436 : SUBROUTINE hfx_admm_init(qs_env, calculate_forces)
164 :
165 : TYPE(qs_environment_type), POINTER :: qs_env
166 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
167 :
168 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_admm_init'
169 :
170 : INTEGER :: handle, ispin, n_rep_hf, nao_aux_fit, &
171 : natoms, nelectron, nmo
172 : LOGICAL :: calc_forces, do_kpoints, &
173 : s_mstruct_changed, use_virial
174 : REAL(dp) :: maxocc
175 : TYPE(admm_type), POINTER :: admm_env
176 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
177 : TYPE(cp_fm_struct_type), POINTER :: aux_fit_fm_struct
178 : TYPE(cp_fm_type), POINTER :: mo_coeff_aux_fit
179 11218 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_aux_fit_kp
180 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
181 : TYPE(dft_control_type), POINTER :: dft_control
182 11218 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
183 : TYPE(mp_para_env_type), POINTER :: para_env
184 11218 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
185 : TYPE(qs_ks_env_type), POINTER :: ks_env
186 : TYPE(section_vals_type), POINTER :: hfx_sections, input, xc_section
187 : TYPE(virial_type), POINTER :: virial
188 :
189 11218 : CALL timeset(routineN, handle)
190 :
191 11218 : NULLIFY (admm_env, hfx_sections, mos, mos_aux_fit, para_env, virial, &
192 11218 : mo_coeff_aux_fit, xc_section, ks_env, dft_control, input, &
193 11218 : qs_kind_set, mo_coeff_b, aux_fit_fm_struct, blacs_env)
194 :
195 : CALL get_qs_env(qs_env, &
196 : mos=mos, &
197 : admm_env=admm_env, &
198 : para_env=para_env, &
199 : blacs_env=blacs_env, &
200 : s_mstruct_changed=s_mstruct_changed, &
201 : ks_env=ks_env, &
202 : dft_control=dft_control, &
203 : input=input, &
204 : virial=virial, &
205 11218 : do_kpoints=do_kpoints)
206 :
207 11218 : calc_forces = .FALSE.
208 11218 : IF (PRESENT(calculate_forces)) calc_forces = .TRUE.
209 :
210 11218 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
211 :
212 11218 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
213 11218 : IF (n_rep_hf > 1) &
214 0 : CPABORT("ADMM can handle only one HF section.")
215 :
216 11218 : IF (.NOT. ASSOCIATED(admm_env)) THEN
217 : ! setup admm environment
218 450 : CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
219 450 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
220 450 : CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
221 450 : CALL set_qs_env(qs_env, admm_env=admm_env)
222 450 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
223 : CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
224 450 : admm_env=admm_env)
225 :
226 : ! Initialize the GAPW data types
227 450 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) &
228 92 : CALL init_admm_gapw(qs_env)
229 :
230 : ! ADMM neighbor lists and overlap matrices
231 450 : CALL admm_init_hamiltonians(admm_env, qs_env, "AUX_FIT")
232 :
233 : !The aux_fit task list and densities
234 450 : ALLOCATE (admm_env%rho_aux_fit)
235 450 : CALL qs_rho_create(admm_env%rho_aux_fit)
236 450 : ALLOCATE (admm_env%rho_aux_fit_buffer)
237 450 : CALL qs_rho_create(admm_env%rho_aux_fit_buffer)
238 450 : CALL admm_update_s_mstruct(admm_env, qs_env, "AUX_FIT")
239 450 : IF (admm_env%do_gapw) CALL update_admm_gapw(qs_env)
240 :
241 : !The ADMM KS matrices
242 450 : CALL admm_alloc_ks_matrices(admm_env, qs_env)
243 :
244 : !The aux_fit MOs and derivatives
245 1916 : ALLOCATE (mos_aux_fit(dft_control%nspins))
246 1016 : DO ispin = 1, dft_control%nspins
247 566 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, nelectron=nelectron, maxocc=maxocc)
248 : CALL allocate_mo_set(mo_set=mos_aux_fit(ispin), &
249 : nao=nao_aux_fit, &
250 : nmo=nmo, &
251 : nelectron=nelectron, &
252 : n_el_f=REAL(nelectron, dp), &
253 : maxocc=maxocc, &
254 1016 : flexible_electron_count=dft_control%relax_multiplicity)
255 : END DO
256 450 : admm_env%mos_aux_fit => mos_aux_fit
257 :
258 1016 : DO ispin = 1, dft_control%nspins
259 566 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
260 : CALL cp_fm_struct_create(aux_fit_fm_struct, context=blacs_env, para_env=para_env, &
261 566 : nrow_global=nao_aux_fit, ncol_global=nmo)
262 566 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, mo_coeff_b=mo_coeff_b)
263 566 : IF (.NOT. ASSOCIATED(mo_coeff_aux_fit)) THEN
264 : CALL init_mo_set(mos_aux_fit(ispin), fm_struct=aux_fit_fm_struct, &
265 566 : name="qs_env%mo_aux_fit"//TRIM(ADJUSTL(cp_to_string(ispin))))
266 : END IF
267 566 : CALL cp_fm_struct_release(aux_fit_fm_struct)
268 :
269 1582 : IF (.NOT. ASSOCIATED(mo_coeff_b)) THEN
270 566 : CALL cp_fm_get_info(mos_aux_fit(ispin)%mo_coeff, ncol_global=nmo)
271 566 : CALL dbcsr_init_p(mos_aux_fit(ispin)%mo_coeff_b)
272 566 : CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit_kp)
273 : CALL cp_dbcsr_m_by_n_from_row_template(mos_aux_fit(ispin)%mo_coeff_b, &
274 : template=matrix_s_aux_fit_kp(1, 1)%matrix, &
275 566 : n=nmo, sym=dbcsr_type_no_symmetry)
276 : END IF
277 : END DO
278 :
279 450 : IF (qs_env%requires_mo_derivs) THEN
280 1046 : ALLOCATE (admm_env%mo_derivs_aux_fit(dft_control%nspins))
281 554 : DO ispin = 1, dft_control%nspins
282 308 : CALL get_mo_set(admm_env%mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
283 554 : CALL cp_fm_create(admm_env%mo_derivs_aux_fit(ispin), mo_coeff_aux_fit%matrix_struct)
284 : END DO
285 : END IF
286 :
287 900 : IF (do_kpoints) THEN
288 : BLOCK
289 : TYPE(kpoint_type), POINTER :: kpoints
290 30 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_aux_fit_kp
291 30 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools_aux_fit
292 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct
293 : INTEGER :: ic, ik, ikk, is
294 : INTEGER, PARAMETER :: nwork1 = 4
295 : LOGICAL :: use_real_wfn
296 :
297 30 : NULLIFY (ao_mo_fm_pools_aux_fit, mos_aux_fit_kp)
298 :
299 30 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
300 30 : CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
301 :
302 : !Test combinations of input values. So far, only ADMM2 is availavle
303 30 : IF (.NOT. admm_env%purification_method == do_admm_purify_none) &
304 0 : CPABORT("Only ADMM_PURIFICATION_METHOD NONE implemeted for ADMM K-points")
305 30 : IF (.NOT. (dft_control%admm_control%method == do_admm_basis_projection &
306 : .OR. dft_control%admm_control%method == do_admm_charge_constrained_projection)) &
307 0 : CPABORT("Only BASIS_PROJECTION and CHARGE_CONSTRAINED_PROJECTION implemented for KP")
308 30 : IF (admm_env%do_admms .OR. admm_env%do_admmp .OR. admm_env%do_admmq) THEN
309 14 : IF (use_real_wfn) CPABORT("Only KP-HFX ADMM2 is implemented with REAL wavefunctions")
310 : END IF
311 :
312 30 : CALL kpoint_initialize_mos(kpoints, admm_env%mos_aux_fit, for_aux_fit=.TRUE.)
313 :
314 30 : CALL mpools_get(kpoints%mpools_aux_fit, ao_mo_fm_pools=ao_mo_fm_pools_aux_fit)
315 244 : DO ik = 1, SIZE(kpoints%kp_aux_env)
316 214 : mos_aux_fit_kp => kpoints%kp_aux_env(ik)%kpoint_env%mos
317 214 : ikk = kpoints%kp_range(1) + ik - 1
318 506 : DO ispin = 1, SIZE(mos_aux_fit_kp, 2)
319 1000 : DO ic = 1, SIZE(mos_aux_fit_kp, 1)
320 524 : CALL get_mo_set(mos_aux_fit_kp(ic, ispin), mo_coeff=mo_coeff_aux_fit, mo_coeff_b=mo_coeff_b)
321 :
322 : ! no sparse matrix representation of kpoint MO vectors
323 524 : CPASSERT(.NOT. ASSOCIATED(mo_coeff_b))
324 :
325 786 : IF (.NOT. ASSOCIATED(mo_coeff_aux_fit)) THEN
326 : CALL init_mo_set(mos_aux_fit_kp(ic, ispin), &
327 : fm_pool=ao_mo_fm_pools_aux_fit(ispin)%pool, &
328 : name="kpoints_"//TRIM(ADJUSTL(cp_to_string(ikk)))// &
329 524 : "%mo_aux_fit"//TRIM(ADJUSTL(cp_to_string(ispin))))
330 : END IF
331 : END DO
332 : END DO
333 : END DO
334 :
335 150 : ALLOCATE (admm_env%scf_work_aux_fit(nwork1))
336 :
337 : ! create an ao_ao parallel matrix structure
338 : CALL cp_fm_struct_create(ao_ao_fm_struct, context=blacs_env, para_env=para_env, &
339 : nrow_global=nao_aux_fit, &
340 30 : ncol_global=nao_aux_fit)
341 :
342 150 : DO is = 1, nwork1
343 : CALL cp_fm_create(admm_env%scf_work_aux_fit(is), &
344 : matrix_struct=ao_ao_fm_struct, &
345 150 : name="SCF-WORK_MATRIX-AUX-"//TRIM(ADJUSTL(cp_to_string(is))))
346 : END DO
347 30 : CALL cp_fm_struct_release(ao_ao_fm_struct)
348 :
349 : ! Create and populate the internal ADMM overlap matrices at each KP
350 60 : CALL kpoint_calc_admm_matrices(qs_env, calc_forces)
351 :
352 : END BLOCK
353 : END IF
354 :
355 10768 : ELSE IF (s_mstruct_changed) THEN
356 360 : CALL admm_init_hamiltonians(admm_env, qs_env, "AUX_FIT")
357 360 : CALL admm_update_s_mstruct(admm_env, qs_env, "AUX_FIT")
358 360 : CALL admm_alloc_ks_matrices(admm_env, qs_env)
359 360 : IF (admm_env%do_gapw) CALL update_admm_gapw(qs_env)
360 360 : IF (do_kpoints) CALL kpoint_calc_admm_matrices(qs_env, calc_forces)
361 : END IF
362 :
363 11218 : IF (admm_env%do_gapw .AND. dft_control%do_admm_dm) THEN
364 0 : CPABORT("GAPW ADMM not implemented for MCWEENY or NONE_DM purification.")
365 : END IF
366 :
367 : !ADMMS and ADMMP stress tensors only available for close-shell systesms, because virial cannot
368 : !be scaled by gsi spin component wise
369 11218 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
370 1192 : IF (use_virial .AND. admm_env%do_admms .AND. dft_control%nspins == 2) THEN
371 0 : CPABORT("ADMMS stress tensor is only available for closed-shell systems")
372 : END IF
373 1192 : IF (use_virial .AND. admm_env%do_admmp .AND. dft_control%nspins == 2) THEN
374 0 : CPABORT("ADMMP stress tensor is only available for closed-shell systems")
375 : END IF
376 :
377 11218 : IF (dft_control%do_admm_dm .AND. .NOT. ASSOCIATED(admm_env%admm_dm)) THEN
378 14 : CALL admm_dm_create(admm_env%admm_dm, dft_control%admm_control, nspins=dft_control%nspins, natoms=natoms)
379 : END IF
380 :
381 11218 : CALL timestop(handle)
382 :
383 11218 : END SUBROUTINE hfx_admm_init
384 :
385 : ! **************************************************************************************************
386 : !> \brief Minimal setup routine for admm_env
387 : !> No forces
388 : !> No k-points
389 : !> No DFT correction terms
390 : !> \param qs_env ...
391 : !> \param mos ...
392 : !> \param admm_env ...
393 : !> \param admm_control ...
394 : !> \param basis_type ...
395 : ! **************************************************************************************************
396 4 : SUBROUTINE aux_admm_init(qs_env, mos, admm_env, admm_control, basis_type)
397 :
398 : TYPE(qs_environment_type), POINTER :: qs_env
399 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
400 : TYPE(admm_type), POINTER :: admm_env
401 : TYPE(admm_control_type), POINTER :: admm_control
402 : CHARACTER(LEN=*) :: basis_type
403 :
404 : CHARACTER(LEN=*), PARAMETER :: routineN = 'aux_admm_init'
405 :
406 : INTEGER :: handle, ispin, nao_aux_fit, natoms, &
407 : nelectron, nmo
408 : LOGICAL :: do_kpoints
409 : REAL(dp) :: maxocc
410 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
411 : TYPE(cp_fm_struct_type), POINTER :: aux_fit_fm_struct
412 : TYPE(cp_fm_type), POINTER :: mo_coeff_aux_fit
413 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_aux_fit_kp
414 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
415 : TYPE(dft_control_type), POINTER :: dft_control
416 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_aux_fit
417 : TYPE(mp_para_env_type), POINTER :: para_env
418 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
419 : TYPE(qs_ks_env_type), POINTER :: ks_env
420 :
421 4 : CALL timeset(routineN, handle)
422 :
423 4 : CPASSERT(.NOT. ASSOCIATED(admm_env))
424 :
425 : CALL get_qs_env(qs_env, &
426 : para_env=para_env, &
427 : blacs_env=blacs_env, &
428 : ks_env=ks_env, &
429 : dft_control=dft_control, &
430 4 : do_kpoints=do_kpoints)
431 :
432 4 : CPASSERT(.NOT. do_kpoints)
433 4 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
434 0 : CPABORT("AUX ADMM not possible with GAPW")
435 : END IF
436 :
437 : ! setup admm environment
438 4 : CALL get_qs_env(qs_env, natom=natoms, qs_kind_set=qs_kind_set)
439 4 : CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type=basis_type)
440 : !
441 4 : CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
442 : ! no XC correction used
443 4 : NULLIFY (admm_env%xc_section_aux, admm_env%xc_section_primary)
444 : ! ADMM neighbor lists and overlap matrices
445 4 : CALL admm_init_hamiltonians(admm_env, qs_env, basis_type)
446 4 : NULLIFY (admm_env%rho_aux_fit, admm_env%rho_aux_fit_buffer)
447 : !The ADMM KS matrices
448 4 : CALL admm_alloc_ks_matrices(admm_env, qs_env)
449 : !The aux_fit MOs and derivatives
450 16 : ALLOCATE (mos_aux_fit(dft_control%nspins))
451 8 : DO ispin = 1, dft_control%nspins
452 4 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, nelectron=nelectron, maxocc=maxocc)
453 : CALL allocate_mo_set(mo_set=mos_aux_fit(ispin), nao=nao_aux_fit, nmo=nmo, &
454 : nelectron=nelectron, n_el_f=REAL(nelectron, dp), &
455 8 : maxocc=maxocc, flexible_electron_count=0.0_dp)
456 : END DO
457 4 : admm_env%mos_aux_fit => mos_aux_fit
458 :
459 8 : DO ispin = 1, dft_control%nspins
460 4 : CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
461 : CALL cp_fm_struct_create(aux_fit_fm_struct, context=blacs_env, para_env=para_env, &
462 4 : nrow_global=nao_aux_fit, ncol_global=nmo)
463 4 : CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, mo_coeff_b=mo_coeff_b)
464 4 : IF (.NOT. ASSOCIATED(mo_coeff_aux_fit)) THEN
465 : CALL init_mo_set(mos_aux_fit(ispin), fm_struct=aux_fit_fm_struct, &
466 4 : name="mo_aux_fit"//TRIM(ADJUSTL(cp_to_string(ispin))))
467 : END IF
468 4 : CALL cp_fm_struct_release(aux_fit_fm_struct)
469 :
470 12 : IF (.NOT. ASSOCIATED(mo_coeff_b)) THEN
471 4 : CALL cp_fm_get_info(mos_aux_fit(ispin)%mo_coeff, ncol_global=nmo)
472 4 : CALL dbcsr_init_p(mos_aux_fit(ispin)%mo_coeff_b)
473 4 : CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit_kp)
474 : CALL cp_dbcsr_m_by_n_from_row_template(mos_aux_fit(ispin)%mo_coeff_b, &
475 : template=matrix_s_aux_fit_kp(1, 1)%matrix, &
476 4 : n=nmo, sym=dbcsr_type_no_symmetry)
477 : END IF
478 : END DO
479 :
480 4 : CALL timestop(handle)
481 :
482 8 : END SUBROUTINE aux_admm_init
483 :
484 : ! **************************************************************************************************
485 : !> \brief Sets up the admm_gapw env
486 : !> \param qs_env ...
487 : ! **************************************************************************************************
488 92 : SUBROUTINE init_admm_gapw(qs_env)
489 :
490 : TYPE(qs_environment_type), POINTER :: qs_env
491 :
492 : INTEGER :: ikind, nkind
493 : TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
494 : TYPE(admm_type), POINTER :: admm_env
495 92 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
496 : TYPE(dft_control_type), POINTER :: dft_control
497 : TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, aux_fit_soft_basis, &
498 : orb_basis, soft_basis
499 : TYPE(mp_para_env_type), POINTER :: para_env
500 92 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: admm_kind_set, qs_kind_set
501 : TYPE(section_vals_type), POINTER :: input
502 :
503 92 : NULLIFY (admm_kind_set, aux_fit_basis, atomic_kind_set, aux_fit_soft_basis, &
504 92 : dft_control, input, orb_basis, para_env, qs_kind_set, soft_basis)
505 :
506 : CALL get_qs_env(qs_env, admm_env=admm_env, &
507 : atomic_kind_set=atomic_kind_set, &
508 : dft_control=dft_control, &
509 : input=input, &
510 : para_env=para_env, &
511 92 : qs_kind_set=qs_kind_set)
512 :
513 92 : admm_env%do_gapw = .TRUE.
514 92 : ALLOCATE (admm_env%admm_gapw_env)
515 92 : admm_gapw_env => admm_env%admm_gapw_env
516 92 : NULLIFY (admm_gapw_env%local_rho_set)
517 92 : NULLIFY (admm_gapw_env%admm_kind_set)
518 92 : NULLIFY (admm_gapw_env%task_list)
519 :
520 : !Create a new kind set for the ADMM stuff (paw_proj soft AUX_FIT basis, etc)
521 92 : nkind = SIZE(qs_kind_set)
522 2292 : ALLOCATE (admm_gapw_env%admm_kind_set(nkind))
523 92 : admm_kind_set => admm_gapw_env%admm_kind_set
524 :
525 : !In this new kind set, we want the AUX_FIT basis to be known as ORB, such that GAPW routines work
526 268 : DO ikind = 1, nkind
527 : !copying over simple data of interest from qs_kind_set
528 176 : admm_kind_set(ikind)%name = qs_kind_set(ikind)%name
529 176 : admm_kind_set(ikind)%element_symbol = qs_kind_set(ikind)%element_symbol
530 176 : admm_kind_set(ikind)%natom = qs_kind_set(ikind)%natom
531 176 : admm_kind_set(ikind)%hard_radius = qs_kind_set(ikind)%hard_radius
532 176 : admm_kind_set(ikind)%max_rad_local = qs_kind_set(ikind)%max_rad_local
533 176 : admm_kind_set(ikind)%gpw_type_forced = qs_kind_set(ikind)%gpw_type_forced
534 176 : admm_kind_set(ikind)%ngrid_rad = qs_kind_set(ikind)%ngrid_rad
535 176 : admm_kind_set(ikind)%ngrid_ang = qs_kind_set(ikind)%ngrid_ang
536 :
537 : !copying potentials of interest from qs_kind_set
538 176 : IF (ASSOCIATED(qs_kind_set(ikind)%all_potential)) THEN
539 40 : CALL copy_potential(qs_kind_set(ikind)%all_potential, admm_kind_set(ikind)%all_potential)
540 : END IF
541 176 : IF (ASSOCIATED(qs_kind_set(ikind)%gth_potential)) THEN
542 136 : CALL copy_potential(qs_kind_set(ikind)%gth_potential, admm_kind_set(ikind)%gth_potential)
543 : END IF
544 176 : IF (ASSOCIATED(qs_kind_set(ikind)%sgp_potential)) THEN
545 0 : CALL copy_potential(qs_kind_set(ikind)%sgp_potential, admm_kind_set(ikind)%sgp_potential)
546 : END IF
547 :
548 176 : NULLIFY (orb_basis)
549 176 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis, basis_type="AUX_FIT")
550 176 : CALL copy_gto_basis_set(aux_fit_basis, orb_basis)
551 268 : CALL add_basis_set_to_container(admm_kind_set(ikind)%basis_sets, orb_basis, "ORB")
552 : END DO
553 :
554 : !Create the corresponding soft basis set (and projectors)
555 : CALL init_gapw_basis_set(admm_kind_set, dft_control%qs_control, input, &
556 92 : modify_qs_control=.FALSE.)
557 :
558 : !Make sure the basis and the projectors are well initialized
559 92 : CALL init_interaction_radii(dft_control%qs_control, admm_kind_set)
560 :
561 : !We also init the atomic grids and harmonics
562 92 : CALL local_rho_set_create(admm_gapw_env%local_rho_set)
563 : CALL init_rho_atom(admm_gapw_env%local_rho_set%rho_atom_set, &
564 92 : atomic_kind_set, admm_kind_set, dft_control, para_env)
565 :
566 : !Make sure that any NLCC potential is well initialized
567 92 : CALL init_gapw_nlcc(admm_kind_set)
568 :
569 : !Need to have access to the soft AUX_FIT basis from the qs_env => add it to the qs_kinds
570 268 : DO ikind = 1, nkind
571 176 : NULLIFY (aux_fit_soft_basis)
572 176 : CALL get_qs_kind(admm_kind_set(ikind), basis_set=soft_basis, basis_type="ORB_SOFT")
573 176 : CALL copy_gto_basis_set(soft_basis, aux_fit_soft_basis)
574 268 : CALL add_basis_set_to_container(qs_kind_set(ikind)%basis_sets, aux_fit_soft_basis, "AUX_FIT_SOFT")
575 : END DO
576 :
577 92 : END SUBROUTINE init_admm_gapw
578 :
579 : ! **************************************************************************************************
580 : !> \brief Builds the ADMM nmeighbor lists and overlap matrix on the model of qs_energies_init_hamiltonians()
581 : !> \param admm_env ...
582 : !> \param qs_env ...
583 : !> \param aux_basis_type ...
584 : ! **************************************************************************************************
585 814 : SUBROUTINE admm_init_hamiltonians(admm_env, qs_env, aux_basis_type)
586 :
587 : TYPE(admm_type), POINTER :: admm_env
588 : TYPE(qs_environment_type), POINTER :: qs_env
589 : CHARACTER(len=*) :: aux_basis_type
590 :
591 : CHARACTER(len=*), PARAMETER :: routineN = 'admm_init_hamiltonians'
592 :
593 : INTEGER :: handle, hfx_pot, ikind, nkind
594 : LOGICAL :: do_kpoints, mic, molecule_only
595 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: aux_fit_present, orb_present
596 : REAL(dp) :: eps_schwarz, omega, pdist, roperator, &
597 : subcells
598 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: aux_fit_radius, orb_radius
599 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
600 814 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
601 : TYPE(cell_type), POINTER :: cell
602 814 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_aux_fit_kp, &
603 814 : matrix_s_aux_fit_vs_orb_kp
604 : TYPE(dft_control_type), POINTER :: dft_control
605 : TYPE(distribution_1d_type), POINTER :: distribution_1d
606 : TYPE(distribution_2d_type), POINTER :: distribution_2d
607 : TYPE(gto_basis_set_type), POINTER :: aux_fit_basis_set, orb_basis_set
608 : TYPE(kpoint_type), POINTER :: kpoints
609 814 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
610 814 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
611 : TYPE(mp_para_env_type), POINTER :: para_env
612 814 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
613 814 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
614 : TYPE(qs_ks_env_type), POINTER :: ks_env
615 : TYPE(section_vals_type), POINTER :: hfx_sections, neighbor_list_section
616 :
617 814 : NULLIFY (particle_set, cell, kpoints, distribution_1d, distribution_2d, molecule_set, &
618 814 : atomic_kind_set, dft_control, neighbor_list_section, aux_fit_basis_set, orb_basis_set, &
619 814 : ks_env, para_env, qs_kind_set, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb_kp)
620 :
621 814 : CALL timeset(routineN, handle)
622 :
623 : CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set, cell=cell, kpoints=kpoints, &
624 : local_particles=distribution_1d, distribution_2d=distribution_2d, &
625 : molecule_set=molecule_set, atomic_kind_set=atomic_kind_set, do_kpoints=do_kpoints, &
626 814 : dft_control=dft_control, para_env=para_env, qs_kind_set=qs_kind_set)
627 3256 : ALLOCATE (orb_present(nkind), aux_fit_present(nkind))
628 5698 : ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind), pair_radius(nkind, nkind))
629 2316 : aux_fit_radius(:) = 0.0_dp
630 :
631 814 : molecule_only = .FALSE.
632 814 : IF (dft_control%qs_control%do_kg) molecule_only = .TRUE.
633 814 : mic = molecule_only
634 814 : IF (kpoints%nkp > 0) THEN
635 44 : mic = .FALSE.
636 770 : ELSEIF (dft_control%qs_control%semi_empirical) THEN
637 0 : mic = .TRUE.
638 : END IF
639 :
640 814 : pdist = dft_control%qs_control%pairlist_radius
641 :
642 814 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
643 814 : neighbor_list_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%NEIGHBOR_LISTS")
644 :
645 3944 : ALLOCATE (atom2d(nkind))
646 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
647 814 : molecule_set, molecule_only, particle_set=particle_set)
648 :
649 2316 : DO ikind = 1, nkind
650 1502 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
651 1502 : IF (ASSOCIATED(orb_basis_set)) THEN
652 1502 : orb_present(ikind) = .TRUE.
653 1502 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
654 : ELSE
655 0 : orb_present(ikind) = .FALSE.
656 : END IF
657 :
658 1502 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type=aux_basis_type)
659 2316 : IF (ASSOCIATED(aux_fit_basis_set)) THEN
660 1502 : aux_fit_present(ikind) = .TRUE.
661 1502 : CALL get_gto_basis_set(gto_basis_set=aux_fit_basis_set, kind_radius=aux_fit_radius(ikind))
662 : ELSE
663 0 : aux_fit_present(ikind) = .FALSE.
664 : END IF
665 : END DO
666 :
667 814 : IF (pdist < 0.0_dp) THEN
668 : pdist = MAX(plane_distance(1, 0, 0, cell), &
669 : plane_distance(0, 1, 0, cell), &
670 2 : plane_distance(0, 0, 1, cell))
671 : END IF
672 :
673 : !In case of K-points, we need to add the HFX potential range to sab_aux_fit, because it is used
674 : !to populate AUX density and KS matrices
675 814 : roperator = 0.0_dp
676 814 : IF (do_kpoints) THEN
677 44 : hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
678 44 : CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", i_val=hfx_pot)
679 :
680 : SELECT CASE (hfx_pot)
681 : CASE (do_potential_id)
682 26 : roperator = 0.0_dp
683 : CASE (do_potential_truncated)
684 32 : CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
685 : CASE (do_potential_mix_cl_trunc)
686 6 : CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
687 : CASE (do_potential_short)
688 0 : CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%OMEGA", r_val=omega)
689 0 : CALL section_vals_val_get(hfx_sections, "SCREENING%EPS_SCHWARZ", r_val=eps_schwarz)
690 0 : CALL erfc_cutoff(eps_schwarz, omega, roperator)
691 : CASE DEFAULT
692 44 : CPABORT("HFX potential not available for K-points (NYI)")
693 : END SELECT
694 : END IF
695 :
696 814 : CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, pair_radius, pdist)
697 5334 : pair_radius = pair_radius + cutoff_screen_factor*roperator
698 : CALL build_neighbor_lists(admm_env%sab_aux_fit, particle_set, atom2d, cell, pair_radius, &
699 814 : mic=mic, molecular=molecule_only, subcells=subcells, nlname="sab_aux_fit")
700 : CALL build_neighbor_lists(admm_env%sab_aux_fit_asymm, particle_set, atom2d, cell, pair_radius, &
701 : mic=mic, symmetric=.FALSE., molecular=molecule_only, subcells=subcells, &
702 814 : nlname="sab_aux_fit_asymm")
703 814 : CALL pair_radius_setup(aux_fit_present, orb_present, aux_fit_radius, orb_radius, pair_radius)
704 : CALL build_neighbor_lists(admm_env%sab_aux_fit_vs_orb, particle_set, atom2d, cell, pair_radius, &
705 : mic=mic, symmetric=.FALSE., molecular=molecule_only, subcells=subcells, &
706 814 : nlname="sab_aux_fit_vs_orb")
707 :
708 : CALL write_neighbor_lists(admm_env%sab_aux_fit, particle_set, cell, para_env, neighbor_list_section, &
709 814 : "/SAB_AUX_FIT", "sab_aux_fit", "AUX_FIT_ORBITAL AUX_FIT_ORBITAL")
710 : CALL write_neighbor_lists(admm_env%sab_aux_fit_vs_orb, particle_set, cell, para_env, neighbor_list_section, &
711 814 : "/SAB_AUX_FIT_VS_ORB", "sab_aux_fit_vs_orb", "ORBITAL AUX_FIT_ORBITAL")
712 :
713 814 : CALL atom2d_cleanup(atom2d)
714 :
715 : !The ADMM overlap matrices (initially in qs_core_hamiltonian.F)
716 814 : CALL get_qs_env(qs_env, ks_env=ks_env)
717 :
718 814 : CALL kpoint_transitional_release(admm_env%matrix_s_aux_fit)
719 : CALL build_overlap_matrix(ks_env, matrixkp_s=matrix_s_aux_fit_kp, &
720 : matrix_name="AUX_FIT_OVERLAP", &
721 : basis_type_a=aux_basis_type, &
722 : basis_type_b=aux_basis_type, &
723 814 : sab_nl=admm_env%sab_aux_fit)
724 814 : CALL set_2d_pointer(admm_env%matrix_s_aux_fit, matrix_s_aux_fit_kp)
725 814 : CALL kpoint_transitional_release(admm_env%matrix_s_aux_fit_vs_orb)
726 : CALL build_overlap_matrix(ks_env, matrixkp_s=matrix_s_aux_fit_vs_orb_kp, &
727 : matrix_name="MIXED_OVERLAP", &
728 : basis_type_a=aux_basis_type, &
729 : basis_type_b="ORB", &
730 814 : sab_nl=admm_env%sab_aux_fit_vs_orb)
731 814 : CALL set_2d_pointer(admm_env%matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp)
732 :
733 814 : CALL timestop(handle)
734 :
735 2442 : END SUBROUTINE admm_init_hamiltonians
736 :
737 : ! **************************************************************************************************
738 : !> \brief Updates the ADMM task_list and density based on the model of qs_env_update_s_mstruct()
739 : !> \param admm_env ...
740 : !> \param qs_env ...
741 : !> \param aux_basis_type ...
742 : ! **************************************************************************************************
743 810 : SUBROUTINE admm_update_s_mstruct(admm_env, qs_env, aux_basis_type)
744 :
745 : TYPE(admm_type), POINTER :: admm_env
746 : TYPE(qs_environment_type), POINTER :: qs_env
747 : CHARACTER(len=*) :: aux_basis_type
748 :
749 : CHARACTER(len=*), PARAMETER :: routineN = 'admm_update_s_mstruct'
750 :
751 : INTEGER :: handle
752 : LOGICAL :: skip_load_balance_distributed
753 : TYPE(dft_control_type), POINTER :: dft_control
754 : TYPE(qs_ks_env_type), POINTER :: ks_env
755 :
756 810 : NULLIFY (ks_env, dft_control)
757 :
758 810 : CALL timeset(routineN, handle)
759 :
760 810 : CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
761 :
762 : !The aux_fit task_list
763 810 : skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
764 810 : IF (ASSOCIATED(admm_env%task_list_aux_fit)) CALL deallocate_task_list(admm_env%task_list_aux_fit)
765 810 : CALL allocate_task_list(admm_env%task_list_aux_fit)
766 : CALL generate_qs_task_list(ks_env, admm_env%task_list_aux_fit, &
767 : reorder_rs_grid_ranks=.FALSE., soft_valid=.FALSE., &
768 : basis_type=aux_basis_type, &
769 : skip_load_balance_distributed=skip_load_balance_distributed, &
770 810 : sab_orb_external=admm_env%sab_aux_fit)
771 :
772 : !The aux_fit densities
773 810 : CALL qs_rho_rebuild(admm_env%rho_aux_fit, qs_env=qs_env, admm=.TRUE.)
774 810 : CALL qs_rho_rebuild(admm_env%rho_aux_fit_buffer, qs_env=qs_env, admm=.TRUE.)
775 :
776 810 : CALL timestop(handle)
777 :
778 810 : END SUBROUTINE admm_update_s_mstruct
779 :
780 : ! **************************************************************************************************
781 : !> \brief Update the admm_gapw_env internals to the current qs_env (i.e. atomic positions)
782 : !> \param qs_env ...
783 : ! **************************************************************************************************
784 252 : SUBROUTINE update_admm_gapw(qs_env)
785 :
786 : TYPE(qs_environment_type), POINTER :: qs_env
787 :
788 : CHARACTER(len=*), PARAMETER :: routineN = 'update_admm_gapw'
789 :
790 : INTEGER :: handle, ikind, nkind
791 : LOGICAL :: paw_atom
792 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: aux_present, oce_present
793 : REAL(dp) :: subcells
794 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: aux_radius, oce_radius
795 252 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
796 : TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
797 : TYPE(admm_type), POINTER :: admm_env
798 252 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
799 : TYPE(cell_type), POINTER :: cell
800 : TYPE(dft_control_type), POINTER :: dft_control
801 : TYPE(distribution_1d_type), POINTER :: distribution_1d
802 : TYPE(distribution_2d_type), POINTER :: distribution_2d
803 : TYPE(gto_basis_set_type), POINTER :: aux_fit_basis
804 252 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
805 252 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
806 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
807 252 : POINTER :: sap_oce
808 252 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
809 : TYPE(paw_proj_set_type), POINTER :: paw_proj
810 252 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: admm_kind_set, qs_kind_set
811 : TYPE(qs_ks_env_type), POINTER :: ks_env
812 :
813 252 : NULLIFY (ks_env, qs_kind_set, admm_kind_set, aux_fit_basis, cell, distribution_1d)
814 252 : NULLIFY (distribution_2d, paw_proj, particle_set, molecule_set, admm_env, admm_gapw_env)
815 252 : NULLIFY (dft_control, atomic_kind_set, sap_oce)
816 :
817 252 : CALL timeset(routineN, handle)
818 :
819 : CALL get_qs_env(qs_env, ks_env=ks_env, qs_kind_set=qs_kind_set, admm_env=admm_env, &
820 252 : dft_control=dft_control)
821 252 : admm_gapw_env => admm_env%admm_gapw_env
822 252 : admm_kind_set => admm_gapw_env%admm_kind_set
823 252 : nkind = SIZE(qs_kind_set)
824 :
825 : !Update the task lisft for the AUX_FIT_SOFT basis
826 252 : IF (ASSOCIATED(admm_gapw_env%task_list)) CALL deallocate_task_list(admm_gapw_env%task_list)
827 252 : CALL allocate_task_list(admm_gapw_env%task_list)
828 :
829 : !note: we set soft_valid to .FALSE. want to use AUX_FIT_SOFT and not the normal ORB SOFT basis
830 : CALL generate_qs_task_list(ks_env, admm_gapw_env%task_list, reorder_rs_grid_ranks=.FALSE., &
831 : soft_valid=.FALSE., basis_type="AUX_FIT_SOFT", &
832 : skip_load_balance_distributed=dft_control%qs_control%skip_load_balance_distributed, &
833 252 : sab_orb_external=admm_env%sab_aux_fit)
834 :
835 : !Update the precomputed oce integrals
836 : !a sap_oce neighbor list is required => build it here
837 1008 : ALLOCATE (aux_present(nkind), oce_present(nkind))
838 1540 : aux_present = .FALSE.; oce_present = .FALSE.
839 1008 : ALLOCATE (aux_radius(nkind), oce_radius(nkind))
840 1540 : aux_radius = 0.0_dp; oce_radius = 0.0_dp
841 :
842 770 : DO ikind = 1, nkind
843 518 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis, basis_type="AUX_FIT")
844 518 : IF (ASSOCIATED(aux_fit_basis)) THEN
845 518 : aux_present(ikind) = .TRUE.
846 518 : CALL get_gto_basis_set(aux_fit_basis, kind_radius=aux_radius(ikind))
847 : END IF
848 :
849 : !note: get oce info from admm_kind_set
850 518 : CALL get_qs_kind(admm_kind_set(ikind), paw_atom=paw_atom, paw_proj_set=paw_proj)
851 770 : IF (paw_atom) THEN
852 324 : oce_present(ikind) = .TRUE.
853 324 : CALL get_paw_proj_set(paw_proj, rcprj=oce_radius(ikind))
854 : END IF
855 : END DO
856 :
857 1008 : ALLOCATE (pair_radius(nkind, nkind))
858 1912 : pair_radius = 0.0_dp
859 252 : CALL pair_radius_setup(aux_present, oce_present, aux_radius, oce_radius, pair_radius)
860 :
861 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
862 : distribution_2d=distribution_2d, local_particles=distribution_1d, &
863 252 : particle_set=particle_set, molecule_set=molecule_set)
864 252 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
865 :
866 1274 : ALLOCATE (atom2d(nkind))
867 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
868 252 : molecule_set, .FALSE., particle_set)
869 : CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, &
870 252 : subcells=subcells, operator_type="ABBA", nlname="AUX_PAW-PRJ")
871 252 : CALL atom2d_cleanup(atom2d)
872 :
873 : !actually compute the oce matrices
874 252 : CALL create_oce_set(admm_gapw_env%oce)
875 252 : CALL allocate_oce_set(admm_gapw_env%oce, nkind)
876 :
877 : !always compute the derivative, cheap anyways
878 : CALL build_oce_matrices(admm_gapw_env%oce%intac, calculate_forces=.TRUE., nder=1, &
879 : qs_kind_set=admm_kind_set, particle_set=particle_set, &
880 252 : sap_oce=sap_oce, eps_fit=dft_control%qs_control%gapw_control%eps_fit)
881 :
882 252 : CALL release_neighbor_list_sets(sap_oce)
883 :
884 252 : CALL timestop(handle)
885 :
886 756 : END SUBROUTINE update_admm_gapw
887 :
888 : ! **************************************************************************************************
889 : !> \brief Allocates the various ADMM KS matrices
890 : !> \param admm_env ...
891 : !> \param qs_env ...
892 : ! **************************************************************************************************
893 814 : SUBROUTINE admm_alloc_ks_matrices(admm_env, qs_env)
894 :
895 : TYPE(admm_type), POINTER :: admm_env
896 : TYPE(qs_environment_type), POINTER :: qs_env
897 :
898 : CHARACTER(len=*), PARAMETER :: routineN = 'admm_alloc_ks_matrices'
899 :
900 : INTEGER :: handle, ic, ispin
901 814 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit_dft_kp, &
902 814 : matrix_ks_aux_fit_hfx_kp, &
903 814 : matrix_ks_aux_fit_kp, &
904 814 : matrix_s_aux_fit_kp
905 : TYPE(dft_control_type), POINTER :: dft_control
906 :
907 814 : NULLIFY (dft_control, matrix_s_aux_fit_kp, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp)
908 :
909 814 : CALL timeset(routineN, handle)
910 :
911 814 : CALL get_qs_env(qs_env, dft_control=dft_control)
912 814 : CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit_kp)
913 :
914 814 : CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit)
915 814 : CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit_dft)
916 814 : CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit_hfx)
917 :
918 814 : CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_kp, dft_control%nspins, dft_control%nimages)
919 814 : CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_dft_kp, dft_control%nspins, dft_control%nimages)
920 814 : CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_hfx_kp, dft_control%nspins, dft_control%nimages)
921 :
922 1796 : DO ispin = 1, dft_control%nspins
923 5460 : DO ic = 1, dft_control%nimages
924 3664 : ALLOCATE (matrix_ks_aux_fit_kp(ispin, ic)%matrix)
925 : CALL dbcsr_create(matrix_ks_aux_fit_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, ic)%matrix, &
926 3664 : name="KOHN-SHAM_MATRIX for ADMM")
927 3664 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
928 3664 : CALL dbcsr_set(matrix_ks_aux_fit_kp(ispin, ic)%matrix, 0.0_dp)
929 :
930 3664 : ALLOCATE (matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix)
931 : CALL dbcsr_create(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, 1)%matrix, &
932 3664 : name="KOHN-SHAM_MATRIX for ADMM")
933 3664 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
934 3664 : CALL dbcsr_set(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, 0.0_dp)
935 :
936 3664 : ALLOCATE (matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix)
937 : CALL dbcsr_create(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, 1)%matrix, &
938 3664 : name="KOHN-SHAM_MATRIX for ADMM")
939 3664 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
940 4646 : CALL dbcsr_set(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, 0.0_dp)
941 : END DO
942 : END DO
943 :
944 : CALL set_admm_env(admm_env, &
945 : matrix_ks_aux_fit_kp=matrix_ks_aux_fit_kp, &
946 : matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft_kp, &
947 814 : matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx_kp)
948 :
949 814 : CALL timestop(handle)
950 :
951 814 : END SUBROUTINE admm_alloc_ks_matrices
952 :
953 : ! **************************************************************************************************
954 : !> \brief Add the HFX K-point contribution to the real-space Hamiltonians
955 : !> \param qs_env ...
956 : !> \param matrix_ks ...
957 : !> \param energy ...
958 : !> \param calculate_forces ...
959 : ! **************************************************************************************************
960 248 : SUBROUTINE hfx_ks_matrix_kp(qs_env, matrix_ks, energy, calculate_forces)
961 : TYPE(qs_environment_type), POINTER :: qs_env
962 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
963 : TYPE(qs_energy_type), POINTER :: energy
964 : LOGICAL, INTENT(in) :: calculate_forces
965 :
966 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ks_matrix_kp'
967 :
968 : INTEGER :: handle, img, irep, ispin, n_rep_hf, &
969 : nimages, nspins
970 : LOGICAL :: do_adiabatic_rescaling, &
971 : s_mstruct_changed, use_virial
972 : REAL(dp) :: eh1, ehfx, eold
973 248 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy
974 248 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit_im, matrix_ks_im
975 248 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks_aux_fit_hfx_kp, &
976 248 : matrix_ks_aux_fit_kp, matrix_ks_orb, &
977 248 : rho_ao_orb
978 : TYPE(dft_control_type), POINTER :: dft_control
979 248 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
980 : TYPE(mp_para_env_type), POINTER :: para_env
981 : TYPE(pw_env_type), POINTER :: pw_env
982 : TYPE(pw_poisson_type), POINTER :: poisson_env
983 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
984 : TYPE(qs_rho_type), POINTER :: rho_orb
985 : TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
986 : hfx_sections, input
987 : TYPE(virial_type), POINTER :: virial
988 :
989 248 : CALL timeset(routineN, handle)
990 :
991 248 : NULLIFY (auxbas_pw_pool, dft_control, hfx_sections, input, &
992 248 : para_env, poisson_env, pw_env, virial, matrix_ks_im, &
993 248 : matrix_ks_orb, rho_ao_orb, matrix_h, matrix_ks_aux_fit_kp, &
994 248 : matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx_kp)
995 :
996 : CALL get_qs_env(qs_env=qs_env, &
997 : dft_control=dft_control, &
998 : input=input, &
999 : matrix_h_kp=matrix_h, &
1000 : para_env=para_env, &
1001 : pw_env=pw_env, &
1002 : virial=virial, &
1003 : matrix_ks_im=matrix_ks_im, &
1004 : s_mstruct_changed=s_mstruct_changed, &
1005 248 : x_data=x_data)
1006 :
1007 : ! No RTP
1008 248 : IF (qs_env%run_rtp) CPABORT("No RTP implementation with K-points HFX")
1009 :
1010 : ! No adiabatic rescaling
1011 248 : adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
1012 248 : CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1013 248 : IF (do_adiabatic_rescaling) CPABORT("No adiabatic rescaling implementation with K-points HFX")
1014 :
1015 248 : IF (dft_control%do_admm) THEN
1016 : CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit_kp, &
1017 : matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, &
1018 138 : matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx_kp)
1019 : END IF
1020 :
1021 248 : nspins = dft_control%nspins
1022 248 : nimages = dft_control%nimages
1023 :
1024 248 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1025 378 : IF (use_virial .AND. calculate_forces) virial%pv_fock_4c = 0.0_dp
1026 :
1027 248 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1028 248 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1029 :
1030 : ! *** Initialize the auxiliary ks matrix to zero if required
1031 248 : IF (dft_control%do_admm) THEN
1032 300 : DO ispin = 1, nspins
1033 8236 : DO img = 1, nimages
1034 8098 : CALL dbcsr_set(matrix_ks_aux_fit_kp(ispin, img)%matrix, 0.0_dp)
1035 : END DO
1036 : END DO
1037 : END IF
1038 580 : DO ispin = 1, nspins
1039 12714 : DO img = 1, nimages
1040 12466 : CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
1041 : END DO
1042 : END DO
1043 :
1044 744 : ALLOCATE (hf_energy(n_rep_hf))
1045 :
1046 248 : eold = 0.0_dp
1047 :
1048 496 : DO irep = 1, n_rep_hf
1049 :
1050 : ! fetch the correct matrices for normal HFX or ADMM
1051 248 : IF (dft_control%do_admm) THEN
1052 138 : CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks_orb, rho_aux_fit=rho_orb)
1053 : ELSE
1054 110 : CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_orb, rho=rho_orb)
1055 : END IF
1056 248 : CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb)
1057 :
1058 : ! Finally the real hfx calulation
1059 : ehfx = 0.0_dp
1060 :
1061 248 : IF (.NOT. x_data(irep, 1)%do_hfx_ri) THEN
1062 0 : CPABORT("Only RI-HFX is implemented for K-points")
1063 : END IF
1064 :
1065 : CALL hfx_ri_update_ks_kp(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1066 : rho_ao_orb, s_mstruct_changed, nspins, &
1067 248 : x_data(irep, 1)%general_parameter%fraction)
1068 :
1069 248 : IF (calculate_forces) THEN
1070 : !Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force
1071 46 : IF (dft_control%do_admm) THEN
1072 26 : CALL scale_dm(qs_env, rho_ao_orb, scale_back=.FALSE.)
1073 : END IF
1074 :
1075 : CALL hfx_ri_update_forces_kp(qs_env, x_data(irep, 1)%ri_data, nspins, &
1076 : x_data(irep, 1)%general_parameter%fraction, &
1077 46 : rho_ao_orb, use_virial=use_virial)
1078 :
1079 46 : IF (dft_control%do_admm) THEN
1080 26 : CALL scale_dm(qs_env, rho_ao_orb, scale_back=.TRUE.)
1081 : END IF
1082 : END IF
1083 :
1084 248 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
1085 248 : eh1 = ehfx - eold
1086 248 : CALL pw_hfx(qs_env, eh1, hfx_sections, poisson_env, auxbas_pw_pool, irep)
1087 744 : eold = ehfx
1088 :
1089 : END DO
1090 :
1091 : ! *** Set the total HFX energy
1092 248 : energy%ex = ehfx
1093 :
1094 : ! *** Add Core-Hamiltonian-Matrix ***
1095 580 : DO ispin = 1, nspins
1096 12714 : DO img = 1, nimages
1097 : CALL dbcsr_add(matrix_ks(ispin, img)%matrix, matrix_h(1, img)%matrix, &
1098 12466 : 1.0_dp, 1.0_dp)
1099 : END DO
1100 : END DO
1101 248 : IF (use_virial .AND. calculate_forces) THEN
1102 130 : virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1103 130 : virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1104 10 : virial%pv_calculate = .FALSE.
1105 : END IF
1106 :
1107 : !update the hfx aux_fit matrix
1108 248 : IF (dft_control%do_admm) THEN
1109 300 : DO ispin = 1, nspins
1110 8236 : DO img = 1, nimages
1111 : CALL dbcsr_add(matrix_ks_aux_fit_hfx_kp(ispin, img)%matrix, matrix_ks_aux_fit_kp(ispin, img)%matrix, &
1112 8098 : 0.0_dp, 1.0_dp)
1113 : END DO
1114 : END DO
1115 : END IF
1116 :
1117 248 : CALL timestop(handle)
1118 :
1119 992 : END SUBROUTINE hfx_ks_matrix_kp
1120 :
1121 : ! **************************************************************************************************
1122 : !> \brief Add the hfx contributions to the Hamiltonian
1123 : !>
1124 : !> \param qs_env ...
1125 : !> \param matrix_ks ...
1126 : !> \param rho ...
1127 : !> \param energy ...
1128 : !> \param calculate_forces ...
1129 : !> \param just_energy ...
1130 : !> \param v_rspace_new ...
1131 : !> \param v_tau_rspace ...
1132 : !> \par History
1133 : !> refactoring 03-2011 [MI]
1134 : ! **************************************************************************************************
1135 :
1136 24468 : SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
1137 : just_energy, v_rspace_new, v_tau_rspace)
1138 :
1139 : TYPE(qs_environment_type), POINTER :: qs_env
1140 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
1141 : TYPE(qs_rho_type), POINTER :: rho
1142 : TYPE(qs_energy_type), POINTER :: energy
1143 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
1144 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
1145 :
1146 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ks_matrix'
1147 :
1148 : INTEGER :: handle, img, irep, ispin, mspin, &
1149 : n_rep_hf, nimages, ns, nspins
1150 : LOGICAL :: distribute_fock_matrix, &
1151 : do_adiabatic_rescaling, &
1152 : hfx_treat_lsd_in_core, &
1153 : s_mstruct_changed, use_virial
1154 : REAL(dp) :: eh1, ehfx, ehfxrt, eold
1155 24468 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy
1156 24468 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_1d, matrix_ks_aux_fit, &
1157 24468 : matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_im, matrix_ks_im, rho_ao_1d, rho_ao_resp
1158 24468 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im, matrix_ks_orb, &
1159 24468 : rho_ao_orb
1160 : TYPE(dft_control_type), POINTER :: dft_control
1161 24468 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1162 24468 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1163 : TYPE(mp_para_env_type), POINTER :: para_env
1164 : TYPE(pw_env_type), POINTER :: pw_env
1165 : TYPE(pw_poisson_type), POINTER :: poisson_env
1166 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1167 : TYPE(qs_rho_type), POINTER :: rho_orb
1168 : TYPE(rt_prop_type), POINTER :: rtp
1169 : TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
1170 : hfx_sections, input
1171 : TYPE(virial_type), POINTER :: virial
1172 :
1173 24468 : CALL timeset(routineN, handle)
1174 :
1175 24468 : NULLIFY (auxbas_pw_pool, dft_control, hfx_sections, input, &
1176 24468 : para_env, poisson_env, pw_env, virial, matrix_ks_im, &
1177 24468 : matrix_ks_orb, rho_ao_orb, matrix_h, matrix_h_im, matrix_ks_aux_fit, &
1178 24468 : matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx)
1179 :
1180 : CALL get_qs_env(qs_env=qs_env, &
1181 : dft_control=dft_control, &
1182 : input=input, &
1183 : matrix_h_kp=matrix_h, &
1184 : matrix_h_im_kp=matrix_h_im, &
1185 : para_env=para_env, &
1186 : pw_env=pw_env, &
1187 : virial=virial, &
1188 : matrix_ks_im=matrix_ks_im, &
1189 : s_mstruct_changed=s_mstruct_changed, &
1190 24468 : x_data=x_data)
1191 :
1192 24468 : IF (dft_control%do_admm) THEN
1193 : CALL get_admm_env(qs_env%admm_env, mos_aux_fit=mo_array, matrix_ks_aux_fit=matrix_ks_aux_fit, &
1194 10856 : matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
1195 : ELSE
1196 13612 : CALL get_qs_env(qs_env=qs_env, mos=mo_array)
1197 : END IF
1198 :
1199 24468 : nspins = dft_control%nspins
1200 24468 : nimages = dft_control%nimages
1201 :
1202 24468 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1203 :
1204 24676 : IF (use_virial .AND. calculate_forces) virial%pv_fock_4c = 0.0_dp
1205 :
1206 24468 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1207 24468 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1208 : CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
1209 24468 : i_rep_section=1)
1210 24468 : adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
1211 24468 : CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1212 :
1213 : ! *** Initialize the auxiliary ks matrix to zero if required
1214 24468 : IF (dft_control%do_admm) THEN
1215 23980 : DO ispin = 1, nspins
1216 23980 : CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
1217 : END DO
1218 : END IF
1219 54174 : DO ispin = 1, nspins
1220 83880 : DO img = 1, nimages
1221 59412 : CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
1222 : END DO
1223 : END DO
1224 :
1225 24468 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1226 :
1227 73404 : ALLOCATE (hf_energy(n_rep_hf))
1228 :
1229 24468 : eold = 0.0_dp
1230 :
1231 48992 : DO irep = 1, n_rep_hf
1232 : ! Remember: Vhfx is added, energy is calclulated from total Vhfx,
1233 : ! so energy of last iteration is correct
1234 :
1235 24524 : IF (do_adiabatic_rescaling .AND. hfx_treat_lsd_in_core) &
1236 0 : CPABORT("HFX_TREAT_LSD_IN_CORE not implemented for adiabatically rescaled hybrids")
1237 : ! everything is calculated with adiabatic rescaling but the potential is not added in a first step
1238 24524 : distribute_fock_matrix = .NOT. do_adiabatic_rescaling
1239 :
1240 24524 : mspin = 1
1241 24524 : IF (hfx_treat_lsd_in_core) mspin = nspins
1242 :
1243 : ! fetch the correct matrices for normal HFX or ADMM
1244 24524 : IF (dft_control%do_admm) THEN
1245 10856 : CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_1d, rho_aux_fit=rho_orb)
1246 10856 : ns = SIZE(matrix_ks_1d)
1247 10856 : matrix_ks_orb(1:ns, 1:1) => matrix_ks_1d(1:ns)
1248 : ELSE
1249 13668 : CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_orb, rho=rho_orb)
1250 : END IF
1251 24524 : CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb)
1252 : ! Finally the real hfx calulation
1253 24524 : ehfx = 0.0_dp
1254 :
1255 24524 : IF (x_data(irep, 1)%do_hfx_ri) THEN
1256 :
1257 : CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1258 : mo_array, rho_ao_orb, &
1259 : s_mstruct_changed, nspins, &
1260 1356 : x_data(irep, 1)%general_parameter%fraction)
1261 1356 : IF (dft_control%do_admm) THEN
1262 : !for ADMMS, we need the exchange matrix k(d) for both spins
1263 382 : DO ispin = 1, nspins
1264 : CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
1265 382 : name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1266 : END DO
1267 : END IF
1268 :
1269 : ELSE
1270 :
1271 46348 : DO ispin = 1, mspin
1272 : CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
1273 : para_env, s_mstruct_changed, irep, distribute_fock_matrix, &
1274 23180 : ispin=ispin)
1275 46348 : ehfx = ehfx + eh1
1276 : END DO
1277 : END IF
1278 :
1279 24524 : IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
1280 : !Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force
1281 724 : IF (dft_control%do_admm) THEN
1282 238 : CALL scale_dm(qs_env, rho_ao_orb, scale_back=.FALSE.)
1283 : END IF
1284 724 : NULLIFY (rho_ao_resp)
1285 :
1286 724 : IF (x_data(irep, 1)%do_hfx_ri) THEN
1287 :
1288 : CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
1289 : x_data(irep, 1)%general_parameter%fraction, &
1290 : rho_ao=rho_ao_orb, mos=mo_array, &
1291 : rho_ao_resp=rho_ao_resp, &
1292 50 : use_virial=use_virial)
1293 :
1294 : ELSE
1295 :
1296 : CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
1297 674 : para_env, irep, use_virial)
1298 :
1299 : END IF
1300 :
1301 : !Scale auxiliary density matrix for ADMMP back with 1/gsi(ispin)
1302 724 : IF (dft_control%do_admm) THEN
1303 238 : CALL scale_dm(qs_env, rho_ao_orb, scale_back=.TRUE.)
1304 : END IF
1305 : END IF
1306 :
1307 : !! If required, the calculation of the forces will be done later with adiabatic rescaling
1308 24524 : IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx
1309 :
1310 : ! special case RTP/EMD we have a full complex density and HFX has a contribution from the imaginary part
1311 24524 : ehfxrt = 0.0_dp
1312 24524 : IF (qs_env%run_rtp) THEN
1313 :
1314 414 : CALL get_qs_env(qs_env=qs_env, rtp=rtp)
1315 860 : DO ispin = 1, nspins
1316 860 : CALL dbcsr_set(matrix_ks_im(ispin)%matrix, 0.0_dp)
1317 : END DO
1318 414 : IF (dft_control%do_admm) THEN
1319 : ! matrix_ks_orb => matrix_ks_aux_fit_im
1320 76 : ns = SIZE(matrix_ks_aux_fit_im)
1321 76 : matrix_ks_orb(1:ns, 1:1) => matrix_ks_aux_fit_im(1:ns)
1322 152 : DO ispin = 1, nspins
1323 152 : CALL dbcsr_set(matrix_ks_aux_fit_im(ispin)%matrix, 0.0_dp)
1324 : END DO
1325 : ELSE
1326 : ! matrix_ks_orb => matrix_ks_im
1327 338 : ns = SIZE(matrix_ks_im)
1328 338 : matrix_ks_orb(1:ns, 1:1) => matrix_ks_im(1:ns)
1329 : END IF
1330 :
1331 414 : CALL qs_rho_get(rho_orb, rho_ao_im=rho_ao_1d)
1332 414 : ns = SIZE(rho_ao_1d)
1333 414 : rho_ao_orb(1:ns, 1:1) => rho_ao_1d(1:ns)
1334 :
1335 414 : ehfxrt = 0.0_dp
1336 :
1337 414 : IF (x_data(irep, 1)%do_hfx_ri) THEN
1338 : CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1339 : mo_array, rho_ao_orb, &
1340 : .FALSE., nspins, &
1341 0 : x_data(irep, 1)%general_parameter%fraction)
1342 0 : IF (dft_control%do_admm) THEN
1343 : !for ADMMS, we need the exchange matrix k(d) for both spins
1344 0 : DO ispin = 1, nspins
1345 : CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
1346 0 : name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1347 : END DO
1348 : END IF
1349 :
1350 : ELSE
1351 828 : DO ispin = 1, mspin
1352 : CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
1353 : para_env, .FALSE., irep, distribute_fock_matrix, &
1354 414 : ispin=ispin)
1355 828 : ehfxrt = ehfxrt + eh1
1356 : END DO
1357 : END IF
1358 :
1359 414 : IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
1360 232 : NULLIFY (rho_ao_resp)
1361 :
1362 232 : IF (x_data(irep, 1)%do_hfx_ri) THEN
1363 :
1364 : CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
1365 : x_data(irep, 1)%general_parameter%fraction, &
1366 : rho_ao=rho_ao_orb, mos=mo_array, &
1367 0 : use_virial=use_virial)
1368 :
1369 : ELSE
1370 : CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
1371 232 : para_env, irep, use_virial)
1372 : END IF
1373 : END IF
1374 :
1375 : !! If required, the calculation of the forces will be done later with adiabatic rescaling
1376 414 : IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx + ehfxrt
1377 :
1378 414 : IF (dft_control%rtp_control%velocity_gauge) THEN
1379 0 : CPASSERT(ASSOCIATED(matrix_h_im))
1380 0 : DO ispin = 1, nspins
1381 : CALL dbcsr_add(matrix_ks_im(ispin)%matrix, matrix_h_im(1, 1)%matrix, &
1382 0 : 1.0_dp, 1.0_dp)
1383 : END DO
1384 : END IF
1385 :
1386 : END IF
1387 :
1388 48992 : IF (.NOT. qs_env%run_rtp) THEN
1389 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1390 24110 : poisson_env=poisson_env)
1391 24110 : eh1 = ehfx - eold
1392 24110 : CALL pw_hfx(qs_env, eh1, hfx_sections, poisson_env, auxbas_pw_pool, irep)
1393 24110 : eold = ehfx
1394 : END IF
1395 :
1396 : END DO
1397 :
1398 : ! *** Set the total HFX energy
1399 24468 : energy%ex = ehfx + ehfxrt
1400 :
1401 : ! *** Add Core-Hamiltonian-Matrix ***
1402 54174 : DO ispin = 1, nspins
1403 83880 : DO img = 1, nimages
1404 : CALL dbcsr_add(matrix_ks(ispin, img)%matrix, matrix_h(1, img)%matrix, &
1405 59412 : 1.0_dp, 1.0_dp)
1406 : END DO
1407 : END DO
1408 24468 : IF (use_virial .AND. calculate_forces) THEN
1409 208 : virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1410 208 : virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1411 16 : virial%pv_calculate = .FALSE.
1412 : END IF
1413 :
1414 : !! If we perform adiabatic rescaling we are now able to rescale the xc-potential
1415 24468 : IF (do_adiabatic_rescaling) THEN
1416 : CALL rescale_xc_potential(qs_env, matrix_ks, rho, energy, v_rspace_new, v_tau_rspace, &
1417 44 : hf_energy, just_energy, calculate_forces, use_virial)
1418 : END IF ! do_adiabatic_rescaling
1419 :
1420 : !update the hfx aux_fit matrixIF (dft_control%do_admm) THEN
1421 24468 : IF (dft_control%do_admm) THEN
1422 23980 : DO ispin = 1, nspins
1423 : CALL dbcsr_add(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, &
1424 23980 : 0.0_dp, 1.0_dp)
1425 : END DO
1426 : END IF
1427 :
1428 24468 : CALL timestop(handle)
1429 :
1430 122340 : END SUBROUTINE hfx_ks_matrix
1431 :
1432 : ! **************************************************************************************************
1433 : !> \brief This routine modifies the xc section depending on the potential type
1434 : !> used for the HF exchange and the resulting correction term. Currently
1435 : !> three types of corrections are implemented:
1436 : !>
1437 : !> coulomb: Ex,hf = Ex,hf' + (PBEx-PBEx')
1438 : !> shortrange: Ex,hf = Ex,hf' + (XWPBEX-XWPBEX')
1439 : !> truncated: Ex,hf = Ex,hf' + ( (XWPBEX0-PBE_HOLE_TC_LR) -(XWPBEX0-PBE_HOLE_TC_LR)' )
1440 : !>
1441 : !> with ' denoting the auxiliary basis set and
1442 : !>
1443 : !> PBEx: PBE exchange functional
1444 : !> XWPBEX: PBE exchange hole for short-range potential (erfc(omega*r)/r)
1445 : !> XWPBEX0: PBE exchange hole for standard coulomb potential
1446 : !> PBE_HOLE_TC_LR: PBE exchange hole for longrange truncated coulomb potential
1447 : !>
1448 : !> Above explanation is correct for the deafult case. If a specific functional is requested
1449 : !> for the correction term (cfun), we get
1450 : !> Ex,hf = Ex,hf' + (cfun-cfun')
1451 : !> for all cases of operators.
1452 : !>
1453 : !> \param x_data ...
1454 : !> \param xc_section the original xc_section
1455 : !> \param admm_env the ADMM environment
1456 : !> \par History
1457 : !> 12.2009 created [Manuel Guidon]
1458 : !> 05.2021 simplify for case of no correction [JGH]
1459 : !> \author Manuel Guidon
1460 : ! **************************************************************************************************
1461 474 : SUBROUTINE create_admm_xc_section(x_data, xc_section, admm_env)
1462 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1463 : TYPE(section_vals_type), POINTER :: xc_section
1464 : TYPE(admm_type), POINTER :: admm_env
1465 :
1466 : LOGICAL, PARAMETER :: debug_functional = .FALSE.
1467 : #if defined (__LIBXC)
1468 : REAL(KIND=dp), PARAMETER :: x_factor_c = 0.930525736349100025_dp
1469 : #endif
1470 :
1471 : CHARACTER(LEN=20) :: name_x_func
1472 : INTEGER :: hfx_potential_type, ifun, iounit, nfun
1473 : LOGICAL :: funct_found
1474 : REAL(dp) :: cutoff_radius, hfx_fraction, omega, &
1475 : scale_coulomb, scale_longrange, scale_x
1476 : TYPE(cp_logger_type), POINTER :: logger
1477 : TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_section
1478 :
1479 474 : logger => cp_get_default_logger()
1480 474 : NULLIFY (admm_env%xc_section_aux, admm_env%xc_section_primary)
1481 :
1482 : !! ** Duplicate existing xc-section
1483 474 : CALL section_vals_duplicate(xc_section, admm_env%xc_section_aux)
1484 474 : CALL section_vals_duplicate(xc_section, admm_env%xc_section_primary)
1485 : !** Now modify the auxiliary basis
1486 : !** First remove all functionals
1487 474 : xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL")
1488 :
1489 : !* Overwrite possible shortcut
1490 : CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1491 474 : i_val=xc_funct_no_shortcut)
1492 :
1493 : !** Get number of Functionals in the list
1494 474 : ifun = 0
1495 474 : nfun = 0
1496 384 : DO
1497 858 : ifun = ifun + 1
1498 858 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1499 858 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1500 384 : nfun = nfun + 1
1501 : END DO
1502 :
1503 : ifun = 0
1504 858 : DO ifun = 1, nfun
1505 384 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=1)
1506 384 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1507 858 : CALL section_vals_remove_values(xc_fun)
1508 : END DO
1509 :
1510 474 : IF (ASSOCIATED(x_data)) THEN
1511 466 : hfx_potential_type = x_data(1, 1)%potential_parameter%potential_type
1512 466 : hfx_fraction = x_data(1, 1)%general_parameter%fraction
1513 : ELSE
1514 8 : CPWARN("ADMM requested without a DFT%XC%HF section. It will be ignored for the SCF.")
1515 8 : admm_env%aux_exch_func = do_admm_aux_exch_func_none
1516 : END IF
1517 :
1518 : !in case of no admm exchange corr., no auxiliary exchange functional needed
1519 474 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
1520 : CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1521 94 : i_val=xc_none)
1522 : hfx_fraction = 0.0_dp
1523 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default) THEN
1524 : ! default PBE Functional
1525 : !! ** Add functionals evaluated with auxiliary basis
1526 184 : SELECT CASE (hfx_potential_type)
1527 : CASE (do_potential_coulomb)
1528 : CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", &
1529 184 : l_val=.TRUE.)
1530 : CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1531 184 : r_val=-hfx_fraction)
1532 : CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", &
1533 184 : r_val=0.0_dp)
1534 : CASE (do_potential_short)
1535 6 : omega = x_data(1, 1)%potential_parameter%omega
1536 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1537 6 : l_val=.TRUE.)
1538 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1539 6 : r_val=-hfx_fraction)
1540 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1541 6 : r_val=0.0_dp)
1542 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1543 6 : r_val=omega)
1544 : CASE (do_potential_truncated)
1545 42 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1546 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1547 42 : l_val=.TRUE.)
1548 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1549 42 : r_val=hfx_fraction)
1550 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1551 42 : r_val=cutoff_radius)
1552 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1553 42 : l_val=.TRUE.)
1554 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1555 42 : r_val=0.0_dp)
1556 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1557 42 : r_val=-hfx_fraction)
1558 : CASE (do_potential_long)
1559 2 : omega = x_data(1, 1)%potential_parameter%omega
1560 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1561 2 : l_val=.TRUE.)
1562 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1563 2 : r_val=hfx_fraction)
1564 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1565 2 : r_val=-hfx_fraction)
1566 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1567 2 : r_val=omega)
1568 : CASE (do_potential_mix_cl)
1569 2 : omega = x_data(1, 1)%potential_parameter%omega
1570 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1571 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1572 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1573 2 : l_val=.TRUE.)
1574 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1575 2 : r_val=hfx_fraction*scale_longrange)
1576 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1577 2 : r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1578 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1579 2 : r_val=omega)
1580 : CASE (do_potential_mix_cl_trunc)
1581 2 : omega = x_data(1, 1)%potential_parameter%omega
1582 2 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1583 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1584 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1585 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1586 2 : l_val=.TRUE.)
1587 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1588 2 : r_val=hfx_fraction*(scale_longrange + scale_coulomb))
1589 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1590 2 : r_val=cutoff_radius)
1591 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1592 2 : l_val=.TRUE.)
1593 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1594 2 : r_val=hfx_fraction*scale_longrange)
1595 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1596 2 : r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1597 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1598 2 : r_val=omega)
1599 : CASE DEFAULT
1600 238 : CPABORT("Unknown potential operator!")
1601 : END SELECT
1602 :
1603 : !** Now modify the functionals for the primary basis
1604 238 : xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
1605 : !* Overwrite possible shortcut
1606 : CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1607 238 : i_val=xc_funct_no_shortcut)
1608 :
1609 184 : SELECT CASE (hfx_potential_type)
1610 : CASE (do_potential_coulomb)
1611 184 : ifun = 0
1612 184 : funct_found = .FALSE.
1613 : DO
1614 338 : ifun = ifun + 1
1615 338 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1616 338 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1617 338 : IF (xc_fun%section%name == "PBE") THEN
1618 148 : funct_found = .TRUE.
1619 : END IF
1620 : END DO
1621 184 : IF (.NOT. funct_found) THEN
1622 : CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", &
1623 36 : l_val=.TRUE.)
1624 : CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1625 36 : r_val=hfx_fraction)
1626 : CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", &
1627 36 : r_val=0.0_dp)
1628 : ELSE
1629 : CALL section_vals_val_get(xc_fun_section, "PBE%SCALE_X", &
1630 148 : r_val=scale_x)
1631 148 : scale_x = scale_x + hfx_fraction
1632 : CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1633 148 : r_val=scale_x)
1634 : END IF
1635 : CASE (do_potential_short)
1636 6 : omega = x_data(1, 1)%potential_parameter%omega
1637 6 : ifun = 0
1638 6 : funct_found = .FALSE.
1639 : DO
1640 18 : ifun = ifun + 1
1641 18 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1642 18 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1643 18 : IF (xc_fun%section%name == "XWPBE") THEN
1644 6 : funct_found = .TRUE.
1645 : END IF
1646 : END DO
1647 6 : IF (.NOT. funct_found) THEN
1648 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1649 0 : l_val=.TRUE.)
1650 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1651 0 : r_val=hfx_fraction)
1652 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1653 0 : r_val=0.0_dp)
1654 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1655 0 : r_val=omega)
1656 : ELSE
1657 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1658 6 : r_val=scale_x)
1659 6 : scale_x = scale_x + hfx_fraction
1660 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1661 6 : r_val=scale_x)
1662 : END IF
1663 : CASE (do_potential_long)
1664 2 : omega = x_data(1, 1)%potential_parameter%omega
1665 2 : ifun = 0
1666 2 : funct_found = .FALSE.
1667 : DO
1668 10 : ifun = ifun + 1
1669 10 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1670 10 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1671 10 : IF (xc_fun%section%name == "XWPBE") THEN
1672 0 : funct_found = .TRUE.
1673 : END IF
1674 : END DO
1675 2 : IF (.NOT. funct_found) THEN
1676 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1677 2 : l_val=.TRUE.)
1678 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1679 2 : r_val=-hfx_fraction)
1680 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1681 2 : r_val=hfx_fraction)
1682 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1683 2 : r_val=omega)
1684 : ELSE
1685 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1686 0 : r_val=scale_x)
1687 0 : scale_x = scale_x - hfx_fraction
1688 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1689 0 : r_val=scale_x)
1690 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1691 0 : r_val=scale_x)
1692 0 : scale_x = scale_x + hfx_fraction
1693 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1694 0 : r_val=scale_x)
1695 :
1696 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1697 0 : r_val=omega)
1698 : END IF
1699 : CASE (do_potential_truncated)
1700 42 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1701 42 : ifun = 0
1702 42 : funct_found = .FALSE.
1703 : DO
1704 62 : ifun = ifun + 1
1705 62 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1706 62 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1707 62 : IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
1708 0 : funct_found = .TRUE.
1709 : END IF
1710 : END DO
1711 42 : IF (.NOT. funct_found) THEN
1712 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1713 42 : l_val=.TRUE.)
1714 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1715 42 : r_val=-hfx_fraction)
1716 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1717 42 : r_val=cutoff_radius)
1718 : ELSE
1719 : CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1720 0 : r_val=scale_x)
1721 0 : scale_x = scale_x - hfx_fraction
1722 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1723 0 : r_val=scale_x)
1724 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1725 0 : r_val=cutoff_radius)
1726 : END IF
1727 42 : ifun = 0
1728 42 : funct_found = .FALSE.
1729 : DO
1730 104 : ifun = ifun + 1
1731 104 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1732 104 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1733 104 : IF (xc_fun%section%name == "XWPBE") THEN
1734 0 : funct_found = .TRUE.
1735 : END IF
1736 : END DO
1737 42 : IF (.NOT. funct_found) THEN
1738 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1739 42 : l_val=.TRUE.)
1740 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1741 42 : r_val=hfx_fraction)
1742 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1743 42 : r_val=0.0_dp)
1744 :
1745 : ELSE
1746 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1747 0 : r_val=scale_x)
1748 0 : scale_x = scale_x + hfx_fraction
1749 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1750 0 : r_val=scale_x)
1751 : END IF
1752 : CASE (do_potential_mix_cl_trunc)
1753 2 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1754 2 : omega = x_data(1, 1)%potential_parameter%omega
1755 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1756 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1757 2 : ifun = 0
1758 2 : funct_found = .FALSE.
1759 : DO
1760 6 : ifun = ifun + 1
1761 6 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1762 6 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1763 6 : IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
1764 0 : funct_found = .TRUE.
1765 : END IF
1766 : END DO
1767 2 : IF (.NOT. funct_found) THEN
1768 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1769 2 : l_val=.TRUE.)
1770 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1771 2 : r_val=-hfx_fraction*(scale_coulomb + scale_longrange))
1772 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1773 2 : r_val=cutoff_radius)
1774 :
1775 : ELSE
1776 : CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1777 0 : r_val=scale_x)
1778 0 : scale_x = scale_x - hfx_fraction*(scale_coulomb + scale_longrange)
1779 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1780 0 : r_val=scale_x)
1781 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1782 0 : r_val=cutoff_radius)
1783 : END IF
1784 2 : ifun = 0
1785 2 : funct_found = .FALSE.
1786 : DO
1787 8 : ifun = ifun + 1
1788 8 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1789 8 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1790 8 : IF (xc_fun%section%name == "XWPBE") THEN
1791 2 : funct_found = .TRUE.
1792 : END IF
1793 : END DO
1794 2 : IF (.NOT. funct_found) THEN
1795 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1796 0 : l_val=.TRUE.)
1797 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1798 0 : r_val=hfx_fraction*(scale_coulomb + scale_longrange))
1799 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1800 0 : r_val=-hfx_fraction*scale_longrange)
1801 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1802 0 : r_val=omega)
1803 :
1804 : ELSE
1805 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1806 2 : r_val=scale_x)
1807 2 : scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
1808 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1809 2 : r_val=scale_x)
1810 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1811 2 : r_val=scale_x)
1812 2 : scale_x = scale_x - hfx_fraction*scale_longrange
1813 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1814 2 : r_val=scale_x)
1815 :
1816 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1817 2 : r_val=omega)
1818 : END IF
1819 : CASE (do_potential_mix_cl)
1820 2 : omega = x_data(1, 1)%potential_parameter%omega
1821 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1822 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1823 2 : ifun = 0
1824 2 : funct_found = .FALSE.
1825 : DO
1826 6 : ifun = ifun + 1
1827 6 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1828 6 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1829 6 : IF (xc_fun%section%name == "XWPBE") THEN
1830 2 : funct_found = .TRUE.
1831 : END IF
1832 : END DO
1833 240 : IF (.NOT. funct_found) THEN
1834 : CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1835 0 : l_val=.TRUE.)
1836 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1837 0 : r_val=hfx_fraction*(scale_coulomb + scale_longrange))
1838 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1839 0 : r_val=-hfx_fraction*scale_longrange)
1840 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1841 0 : r_val=omega)
1842 :
1843 : ELSE
1844 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1845 2 : r_val=scale_x)
1846 2 : scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
1847 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1848 2 : r_val=scale_x)
1849 :
1850 : CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1851 2 : r_val=scale_x)
1852 2 : scale_x = scale_x - hfx_fraction*scale_longrange
1853 : CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1854 2 : r_val=scale_x)
1855 :
1856 : CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1857 2 : r_val=omega)
1858 : END IF
1859 : END SELECT
1860 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default_libxc) THEN
1861 : ! default PBE Functional
1862 : !! ** Add functionals evaluated with auxiliary basis
1863 : #if defined (__LIBXC)
1864 2 : SELECT CASE (hfx_potential_type)
1865 : CASE (do_potential_coulomb)
1866 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1867 2 : l_val=.TRUE.)
1868 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1869 2 : r_val=-hfx_fraction)
1870 : CASE (do_potential_short)
1871 2 : omega = x_data(1, 1)%potential_parameter%omega
1872 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1873 2 : l_val=.TRUE.)
1874 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1875 2 : r_val=-hfx_fraction)
1876 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1877 2 : r_val=omega)
1878 : CASE (do_potential_truncated)
1879 0 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1880 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1881 0 : l_val=.TRUE.)
1882 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1883 0 : r_val=hfx_fraction)
1884 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1885 0 : r_val=cutoff_radius)
1886 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1887 0 : l_val=.TRUE.)
1888 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1889 0 : r_val=-hfx_fraction)
1890 : CASE (do_potential_long)
1891 2 : omega = x_data(1, 1)%potential_parameter%omega
1892 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1893 2 : l_val=.TRUE.)
1894 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1895 2 : r_val=hfx_fraction)
1896 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1897 2 : r_val=omega)
1898 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1899 2 : l_val=.TRUE.)
1900 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1901 2 : r_val=-hfx_fraction)
1902 : CASE (do_potential_mix_cl)
1903 2 : omega = x_data(1, 1)%potential_parameter%omega
1904 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1905 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1906 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1907 2 : l_val=.TRUE.)
1908 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1909 2 : r_val=hfx_fraction*scale_longrange)
1910 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1911 2 : r_val=omega)
1912 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1913 2 : l_val=.TRUE.)
1914 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1915 2 : r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1916 : CASE (do_potential_mix_cl_trunc)
1917 2 : omega = x_data(1, 1)%potential_parameter%omega
1918 2 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1919 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1920 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1921 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1922 2 : l_val=.TRUE.)
1923 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1924 2 : r_val=hfx_fraction*(scale_longrange + scale_coulomb))
1925 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1926 2 : r_val=cutoff_radius)
1927 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1928 2 : l_val=.TRUE.)
1929 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1930 2 : r_val=hfx_fraction*scale_longrange)
1931 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1932 2 : r_val=omega)
1933 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1934 2 : l_val=.TRUE.)
1935 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1936 2 : r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1937 : CASE DEFAULT
1938 10 : CPABORT("Unknown potential operator!")
1939 : END SELECT
1940 :
1941 : !** Now modify the functionals for the primary basis
1942 10 : xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
1943 : !* Overwrite possible shortcut
1944 : CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1945 10 : i_val=xc_funct_no_shortcut)
1946 :
1947 2 : SELECT CASE (hfx_potential_type)
1948 : CASE (do_potential_coulomb)
1949 2 : ifun = 0
1950 2 : funct_found = .FALSE.
1951 : DO
1952 4 : ifun = ifun + 1
1953 4 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1954 4 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1955 4 : IF (xc_fun%section%name == "GGA_X_PBE") THEN
1956 0 : funct_found = .TRUE.
1957 : END IF
1958 : END DO
1959 2 : IF (.NOT. funct_found) THEN
1960 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1961 2 : l_val=.TRUE.)
1962 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1963 2 : r_val=hfx_fraction)
1964 : ELSE
1965 : CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
1966 0 : r_val=scale_x)
1967 0 : scale_x = scale_x + hfx_fraction
1968 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1969 0 : r_val=scale_x)
1970 : END IF
1971 : CASE (do_potential_short)
1972 2 : omega = x_data(1, 1)%potential_parameter%omega
1973 2 : ifun = 0
1974 2 : funct_found = .FALSE.
1975 : DO
1976 4 : ifun = ifun + 1
1977 4 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1978 4 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1979 4 : IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
1980 0 : funct_found = .TRUE.
1981 : END IF
1982 : END DO
1983 2 : IF (.NOT. funct_found) THEN
1984 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1985 2 : l_val=.TRUE.)
1986 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1987 2 : r_val=hfx_fraction)
1988 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1989 2 : r_val=omega)
1990 : ELSE
1991 : CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1992 0 : r_val=scale_x)
1993 0 : scale_x = scale_x + hfx_fraction
1994 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1995 0 : r_val=scale_x)
1996 : END IF
1997 : CASE (do_potential_long)
1998 2 : omega = x_data(1, 1)%potential_parameter%omega
1999 2 : ifun = 0
2000 2 : funct_found = .FALSE.
2001 : DO
2002 4 : ifun = ifun + 1
2003 4 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2004 4 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2005 4 : IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2006 0 : funct_found = .TRUE.
2007 : END IF
2008 : END DO
2009 2 : IF (.NOT. funct_found) THEN
2010 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2011 2 : l_val=.TRUE.)
2012 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2013 2 : r_val=-hfx_fraction)
2014 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2015 2 : r_val=omega)
2016 : ELSE
2017 : CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2018 0 : r_val=scale_x)
2019 0 : scale_x = scale_x - hfx_fraction
2020 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2021 0 : r_val=scale_x)
2022 :
2023 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2024 0 : r_val=omega)
2025 : END IF
2026 2 : ifun = 0
2027 2 : funct_found = .FALSE.
2028 : DO
2029 6 : ifun = ifun + 1
2030 6 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2031 6 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2032 6 : IF (xc_fun%section%name == "GGA_X_PBE") THEN
2033 0 : funct_found = .TRUE.
2034 : END IF
2035 : END DO
2036 2 : IF (.NOT. funct_found) THEN
2037 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2038 2 : l_val=.TRUE.)
2039 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2040 2 : r_val=hfx_fraction)
2041 : ELSE
2042 : CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2043 0 : r_val=scale_x)
2044 0 : scale_x = scale_x + hfx_fraction
2045 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2046 0 : r_val=scale_x)
2047 : END IF
2048 : CASE (do_potential_truncated)
2049 0 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
2050 0 : ifun = 0
2051 0 : funct_found = .FALSE.
2052 : DO
2053 0 : ifun = ifun + 1
2054 0 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2055 0 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2056 0 : IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
2057 0 : funct_found = .TRUE.
2058 : END IF
2059 : END DO
2060 0 : IF (.NOT. funct_found) THEN
2061 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
2062 0 : l_val=.TRUE.)
2063 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2064 0 : r_val=-hfx_fraction)
2065 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2066 0 : r_val=cutoff_radius)
2067 :
2068 : ELSE
2069 : CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2070 0 : r_val=scale_x)
2071 0 : scale_x = scale_x - hfx_fraction
2072 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2073 0 : r_val=scale_x)
2074 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2075 0 : r_val=cutoff_radius)
2076 : END IF
2077 0 : ifun = 0
2078 0 : funct_found = .FALSE.
2079 : DO
2080 0 : ifun = ifun + 1
2081 0 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2082 0 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2083 0 : IF (xc_fun%section%name == "GGA_X_PBE") THEN
2084 0 : funct_found = .TRUE.
2085 : END IF
2086 : END DO
2087 0 : IF (.NOT. funct_found) THEN
2088 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2089 0 : l_val=.TRUE.)
2090 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2091 0 : r_val=hfx_fraction)
2092 :
2093 : ELSE
2094 : CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2095 0 : r_val=scale_x)
2096 0 : scale_x = scale_x + hfx_fraction
2097 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2098 0 : r_val=scale_x)
2099 : END IF
2100 : CASE (do_potential_mix_cl_trunc)
2101 2 : cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
2102 2 : omega = x_data(1, 1)%potential_parameter%omega
2103 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
2104 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
2105 2 : ifun = 0
2106 2 : funct_found = .FALSE.
2107 : DO
2108 4 : ifun = ifun + 1
2109 4 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2110 4 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2111 4 : IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
2112 0 : funct_found = .TRUE.
2113 : END IF
2114 : END DO
2115 2 : IF (.NOT. funct_found) THEN
2116 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
2117 2 : l_val=.TRUE.)
2118 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2119 2 : r_val=-hfx_fraction*(scale_coulomb + scale_longrange))
2120 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2121 2 : r_val=cutoff_radius)
2122 :
2123 : ELSE
2124 : CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2125 0 : r_val=scale_x)
2126 0 : scale_x = scale_x - hfx_fraction*(scale_coulomb + scale_longrange)
2127 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2128 0 : r_val=scale_x)
2129 : CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2130 0 : r_val=cutoff_radius)
2131 : END IF
2132 2 : ifun = 0
2133 2 : funct_found = .FALSE.
2134 : DO
2135 6 : ifun = ifun + 1
2136 6 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2137 6 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2138 6 : IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2139 0 : funct_found = .TRUE.
2140 : END IF
2141 : END DO
2142 2 : IF (.NOT. funct_found) THEN
2143 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2144 2 : l_val=.TRUE.)
2145 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2146 2 : r_val=-hfx_fraction*scale_longrange)
2147 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2148 2 : r_val=omega)
2149 :
2150 : ELSE
2151 : CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2152 0 : r_val=scale_x)
2153 0 : scale_x = scale_x - hfx_fraction*scale_longrange
2154 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2155 0 : r_val=scale_x)
2156 :
2157 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2158 0 : r_val=omega)
2159 : END IF
2160 2 : ifun = 0
2161 2 : funct_found = .FALSE.
2162 : DO
2163 8 : ifun = ifun + 1
2164 8 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2165 8 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2166 8 : IF (xc_fun%section%name == "GGA_X_PBE") THEN
2167 0 : funct_found = .TRUE.
2168 : END IF
2169 : END DO
2170 2 : IF (.NOT. funct_found) THEN
2171 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2172 2 : l_val=.TRUE.)
2173 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2174 2 : r_val=hfx_fraction*(scale_coulomb + scale_longrange))
2175 : ELSE
2176 : CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2177 0 : r_val=scale_x)
2178 0 : scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
2179 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2180 0 : r_val=scale_x)
2181 : END IF
2182 : CASE (do_potential_mix_cl)
2183 2 : omega = x_data(1, 1)%potential_parameter%omega
2184 2 : scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
2185 2 : scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
2186 2 : ifun = 0
2187 2 : funct_found = .FALSE.
2188 : DO
2189 4 : ifun = ifun + 1
2190 4 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2191 4 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2192 4 : IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2193 0 : funct_found = .TRUE.
2194 : END IF
2195 : END DO
2196 2 : IF (.NOT. funct_found) THEN
2197 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2198 2 : l_val=.TRUE.)
2199 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2200 2 : r_val=-hfx_fraction*scale_longrange)
2201 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2202 2 : r_val=omega)
2203 :
2204 : ELSE
2205 : CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2206 0 : r_val=scale_x)
2207 0 : scale_x = scale_x - hfx_fraction*scale_longrange
2208 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2209 0 : r_val=scale_x)
2210 :
2211 : CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2212 0 : r_val=omega)
2213 : END IF
2214 2 : ifun = 0
2215 2 : funct_found = .FALSE.
2216 : DO
2217 6 : ifun = ifun + 1
2218 6 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2219 6 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2220 6 : IF (xc_fun%section%name == "GGA_X_PBE") THEN
2221 0 : funct_found = .TRUE.
2222 : END IF
2223 : END DO
2224 12 : IF (.NOT. funct_found) THEN
2225 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2226 2 : l_val=.TRUE.)
2227 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2228 2 : r_val=hfx_fraction*(scale_coulomb + scale_longrange))
2229 : ELSE
2230 : CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2231 0 : r_val=scale_x)
2232 0 : scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
2233 : CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2234 0 : r_val=scale_x)
2235 : END IF
2236 : END SELECT
2237 : #else
2238 : CALL cp_abort(__LOCATION__, "In order use a LibXC-based ADMM "// &
2239 : "exchange correction functionals, you have to compile and link against LibXC!")
2240 : #endif
2241 :
2242 : ! PBEX (always bare form), OPTX and Becke88 functional
2243 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex .OR. &
2244 : admm_env%aux_exch_func == do_admm_aux_exch_func_opt .OR. &
2245 : admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN
2246 120 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2247 90 : name_x_func = 'PBE'
2248 30 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2249 14 : name_x_func = 'OPTX'
2250 16 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN
2251 16 : name_x_func = 'BECKE88'
2252 : END IF
2253 : !primary basis
2254 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_SECTION_PARAMETERS_", &
2255 120 : l_val=.TRUE.)
2256 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", &
2257 120 : r_val=-hfx_fraction)
2258 :
2259 120 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2260 90 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_C", r_val=0.0_dp)
2261 : END IF
2262 :
2263 120 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2264 14 : IF (admm_env%aux_exch_func_param) THEN
2265 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A1", &
2266 0 : r_val=admm_env%aux_x_param(1))
2267 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A2", &
2268 0 : r_val=admm_env%aux_x_param(2))
2269 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%GAMMA", &
2270 0 : r_val=admm_env%aux_x_param(3))
2271 : END IF
2272 : END IF
2273 :
2274 : !** Now modify the functionals for the primary basis
2275 120 : xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2276 : !* Overwrite possible L")
2277 : !* Overwrite possible shortcut
2278 : CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
2279 120 : i_val=xc_funct_no_shortcut)
2280 :
2281 120 : ifun = 0
2282 120 : funct_found = .FALSE.
2283 : DO
2284 210 : ifun = ifun + 1
2285 210 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2286 210 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2287 210 : IF (xc_fun%section%name == TRIM(name_x_func)) THEN
2288 44 : funct_found = .TRUE.
2289 : END IF
2290 : END DO
2291 120 : IF (.NOT. funct_found) THEN
2292 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_SECTION_PARAMETERS_", &
2293 76 : l_val=.TRUE.)
2294 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", &
2295 76 : r_val=hfx_fraction)
2296 76 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2297 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_C", &
2298 48 : r_val=0.0_dp)
2299 28 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2300 14 : IF (admm_env%aux_exch_func_param) THEN
2301 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A1", &
2302 0 : r_val=admm_env%aux_x_param(1))
2303 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A2", &
2304 0 : r_val=admm_env%aux_x_param(2))
2305 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%GAMMA", &
2306 0 : r_val=admm_env%aux_x_param(3))
2307 : END IF
2308 : END IF
2309 :
2310 : ELSE
2311 : CALL section_vals_val_get(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", &
2312 44 : r_val=scale_x)
2313 44 : scale_x = scale_x + hfx_fraction
2314 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", &
2315 44 : r_val=scale_x)
2316 44 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2317 0 : CPASSERT(.NOT. admm_env%aux_exch_func_param)
2318 : END IF
2319 : END IF
2320 :
2321 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc .OR. &
2322 : admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc .OR. &
2323 : admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc .OR. &
2324 : admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
2325 : #if defined(__LIBXC)
2326 12 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc) THEN
2327 2 : name_x_func = 'GGA_X_PBE'
2328 10 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2329 2 : name_x_func = 'GGA_X_OPTX'
2330 8 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
2331 2 : name_x_func = 'GGA_X_B88'
2332 6 : ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc) THEN
2333 6 : name_x_func = 'LDA_X'
2334 : END IF
2335 : !primary basis
2336 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_SECTION_PARAMETERS_", &
2337 12 : l_val=.TRUE.)
2338 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE", &
2339 12 : r_val=-hfx_fraction)
2340 :
2341 12 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2342 2 : IF (admm_env%aux_exch_func_param) THEN
2343 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_A", &
2344 0 : r_val=admm_env%aux_x_param(1))
2345 : ! LibXC rescales the second parameter of the OPTX functional (see documentation there)
2346 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_B", &
2347 0 : r_val=admm_env%aux_x_param(2)/x_factor_c)
2348 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_GAMMA", &
2349 0 : r_val=admm_env%aux_x_param(3))
2350 : END IF
2351 : END IF
2352 :
2353 : !** Now modify the functionals for the primary basis
2354 12 : xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2355 : !* Overwrite possible L")
2356 : !* Overwrite possible shortcut
2357 : CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
2358 12 : i_val=xc_funct_no_shortcut)
2359 :
2360 12 : ifun = 0
2361 12 : funct_found = .FALSE.
2362 : DO
2363 24 : ifun = ifun + 1
2364 24 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2365 24 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2366 24 : IF (xc_fun%section%name == TRIM(name_x_func)) THEN
2367 0 : funct_found = .TRUE.
2368 : END IF
2369 : END DO
2370 12 : IF (.NOT. funct_found) THEN
2371 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_SECTION_PARAMETERS_", &
2372 12 : l_val=.TRUE.)
2373 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE", &
2374 12 : r_val=hfx_fraction)
2375 12 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2376 2 : IF (admm_env%aux_exch_func_param) THEN
2377 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_A", &
2378 0 : r_val=admm_env%aux_x_param(1))
2379 : ! LibXC rescales the second parameter of the OPTX functional (see documentation there)
2380 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_B", &
2381 0 : r_val=admm_env%aux_x_param(2)/x_factor_c)
2382 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_GAMMA", &
2383 0 : r_val=admm_env%aux_x_param(3))
2384 : END IF
2385 : END IF
2386 :
2387 : ELSE
2388 : CALL section_vals_val_get(xc_fun_section, TRIM(name_x_func)//"%SCALE", &
2389 0 : r_val=scale_x)
2390 0 : scale_x = scale_x + hfx_fraction
2391 : CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE", &
2392 0 : r_val=scale_x)
2393 0 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2394 0 : CPASSERT(.NOT. admm_env%aux_exch_func_param)
2395 : END IF
2396 : END IF
2397 : #else
2398 : CALL cp_abort(__LOCATION__, "In order use a LibXC-based ADMM "// &
2399 : "exchange correction functionals, you have to compile and link against LibXC!")
2400 : #endif
2401 :
2402 : ELSE
2403 0 : CPABORT("Unknown exchange correction functional!")
2404 : END IF
2405 :
2406 : IF (debug_functional) THEN
2407 : iounit = cp_logger_get_default_io_unit(logger)
2408 : IF (iounit > 0) THEN
2409 : WRITE (iounit, "(A)") " ADMM Primary Basis Set Functional"
2410 : END IF
2411 : xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2412 : ifun = 0
2413 : funct_found = .FALSE.
2414 : DO
2415 : ifun = ifun + 1
2416 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2417 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2418 :
2419 : scale_x = -1000.0_dp
2420 : IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
2421 : CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
2422 : END IF
2423 : IF (xc_fun%section%name == "XWPBE") THEN
2424 : CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
2425 : IF (iounit > 0) THEN
2426 : WRITE (iounit, "(T5,A,T25,2F10.3)") TRIM(xc_fun%section%name), scale_x, hfx_fraction
2427 : END IF
2428 : ELSE
2429 : IF (iounit > 0) THEN
2430 : WRITE (iounit, "(T5,A,T25,F10.3)") TRIM(xc_fun%section%name), scale_x
2431 : END IF
2432 : END IF
2433 : END DO
2434 :
2435 : IF (iounit > 0) THEN
2436 : WRITE (iounit, "(A)") " Auxiliary Basis Set Functional"
2437 : END IF
2438 : xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL")
2439 : ifun = 0
2440 : funct_found = .FALSE.
2441 : DO
2442 : ifun = ifun + 1
2443 : xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2444 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2445 : scale_x = -1000.0_dp
2446 : IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
2447 : CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
2448 : END IF
2449 : IF (xc_fun%section%name == "XWPBE") THEN
2450 : CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
2451 : IF (iounit > 0) THEN
2452 : WRITE (iounit, "(T5,A,T25,2F10.3)") TRIM(xc_fun%section%name), scale_x, hfx_fraction
2453 : END IF
2454 : ELSE
2455 : IF (iounit > 0) THEN
2456 : WRITE (iounit, "(T5,A,T25,F10.3)") TRIM(xc_fun%section%name), scale_x
2457 : END IF
2458 : END IF
2459 : END DO
2460 : END IF
2461 :
2462 474 : END SUBROUTINE create_admm_xc_section
2463 :
2464 : ! **************************************************************************************************
2465 : !> \brief Add the hfx contributions to the Hamiltonian
2466 : !>
2467 : !> \param matrix_ks Kohn-Sham matrix (updated on exit)
2468 : !> \param rho_ao electron density expressed in terms of atomic orbitals
2469 : !> \param qs_env Quickstep environment
2470 : !> \param update_energy whether to update energy (default: yes)
2471 : !> \param recalc_integrals whether to recalculate integrals (default: value of HF%TREAT_LSD_IN_CORE)
2472 : !> \param external_hfx_sections ...
2473 : !> \param external_x_data ...
2474 : !> \param external_para_env ...
2475 : !> \note
2476 : !> Simplified version of subroutine hfx_ks_matrix()
2477 : ! **************************************************************************************************
2478 7125 : SUBROUTINE tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, &
2479 7125 : external_hfx_sections, external_x_data, external_para_env)
2480 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
2481 : TARGET :: matrix_ks, rho_ao
2482 : TYPE(qs_environment_type), POINTER :: qs_env
2483 : LOGICAL, INTENT(IN), OPTIONAL :: update_energy, recalc_integrals
2484 : TYPE(section_vals_type), OPTIONAL, POINTER :: external_hfx_sections
2485 : TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, TARGET :: external_x_data
2486 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: external_para_env
2487 :
2488 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddft_hfx_matrix'
2489 :
2490 : INTEGER :: handle, irep, ispin, mspin, n_rep_hf, &
2491 : nspins
2492 : LOGICAL :: distribute_fock_matrix, &
2493 : hfx_treat_lsd_in_core, &
2494 : my_update_energy, s_mstruct_changed
2495 : REAL(KIND=dp) :: eh1, ehfx
2496 7125 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp
2497 : TYPE(dft_control_type), POINTER :: dft_control
2498 7125 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
2499 : TYPE(mp_para_env_type), POINTER :: para_env
2500 : TYPE(qs_energy_type), POINTER :: energy
2501 : TYPE(section_vals_type), POINTER :: hfx_sections, input
2502 :
2503 7125 : CALL timeset(routineN, handle)
2504 :
2505 7125 : NULLIFY (dft_control, hfx_sections, input, para_env, matrix_ks_kp, rho_ao_kp)
2506 :
2507 : CALL get_qs_env(qs_env=qs_env, &
2508 : dft_control=dft_control, &
2509 : energy=energy, &
2510 : input=input, &
2511 : para_env=para_env, &
2512 : s_mstruct_changed=s_mstruct_changed, &
2513 7125 : x_data=x_data)
2514 :
2515 : ! This should probably be the HF section from the TDDFPT XC section!
2516 7125 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
2517 :
2518 7125 : IF (PRESENT(external_hfx_sections)) hfx_sections => external_hfx_sections
2519 7125 : IF (PRESENT(external_x_data)) x_data => external_x_data
2520 7125 : IF (PRESENT(external_para_env)) para_env => external_para_env
2521 :
2522 7125 : my_update_energy = .TRUE.
2523 7125 : IF (PRESENT(update_energy)) my_update_energy = update_energy
2524 :
2525 7125 : IF (PRESENT(recalc_integrals)) s_mstruct_changed = recalc_integrals
2526 :
2527 7125 : CPASSERT(dft_control%nimages == 1)
2528 7125 : nspins = dft_control%nspins
2529 :
2530 7125 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
2531 : CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
2532 7125 : i_rep_section=1)
2533 :
2534 7125 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
2535 7125 : distribute_fock_matrix = .TRUE.
2536 :
2537 7125 : mspin = 1
2538 7125 : IF (hfx_treat_lsd_in_core) mspin = nspins
2539 :
2540 7125 : matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
2541 7125 : rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
2542 :
2543 14110 : DO irep = 1, n_rep_hf
2544 : ! the real hfx calulation
2545 6985 : ehfx = 0.0_dp
2546 :
2547 14110 : IF (x_data(irep, 1)%do_hfx_ri) THEN
2548 : CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, &
2549 : rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
2550 308 : nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
2551 :
2552 : ELSE
2553 13354 : DO ispin = 1, mspin
2554 : CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
2555 6677 : s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
2556 13354 : ehfx = ehfx + eh1
2557 : END DO
2558 : END IF
2559 : END DO
2560 7125 : IF (my_update_energy) energy%ex = ehfx
2561 :
2562 7125 : CALL timestop(handle)
2563 7125 : END SUBROUTINE tddft_hfx_matrix
2564 :
2565 : END MODULE hfx_admm_utils
|