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