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