(git:4672866)
Loading...
Searching...
No Matches
hfx_admm_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Utilities for hfx and admm methods
10!>
11!>
12!> \par History
13!> refactoring 03-2011 [MI]
14!> Made GAPW compatible 12.2019 (A. Bussy)
15!> \author MI
16! **************************************************************************************************
21 USE admm_types, ONLY: admm_env_create,&
23 admm_type,&
31 USE cell_types, ONLY: cell_type,&
36 USE cp_dbcsr_api, ONLY: dbcsr_add,&
41 dbcsr_set,&
43 dbcsr_type_no_symmetry
51 USE cp_fm_types, ONLY: cp_fm_create,&
63 USE hfx_pw_methods, ONLY: pw_hfx
64 USE hfx_ri, ONLY: hfx_ri_update_forces,&
68 USE hfx_types, ONLY: hfx_type
69 USE input_constants, ONLY: &
86 USE kinds, ONLY: dp
90 USE kpoint_types, ONLY: get_kpoint_info,&
93 USE mathlib, ONLY: erfc_cutoff
99 USE pw_env_types, ONLY: pw_env_get,&
102 USE pw_pool_types, ONLY: pw_pool_type
103 USE pw_types, ONLY: pw_r3d_rs_type
109 USE qs_kind_types, ONLY: get_qs_kind,&
114 USE qs_ks_types, ONLY: qs_ks_env_type
116 USE qs_matrix_pools, ONLY: mpools_get
117 USE qs_mo_types, ONLY: allocate_mo_set,&
118 get_mo_set,&
130 USE qs_oce_types, ONLY: allocate_oce_set,&
135 USE qs_rho_types, ONLY: qs_rho_create,&
136 qs_rho_get,&
142 USE virial_types, ONLY: virial_type
144#include "./base/base_uses.f90"
145
146 IMPLICIT NONE
147
148 PRIVATE
149
150 ! *** Public subroutines ***
153
154 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_admm_utils'
155
156CONTAINS
157
158! **************************************************************************************************
159!> \brief ...
160!> \param qs_env ...
161!> \param calculate_forces ...
162!> \param ext_xc_section ...
163! **************************************************************************************************
164 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 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 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_aux_fit
185 TYPE(mp_para_env_type), POINTER :: para_env
186 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 CALL timeset(routinen, handle)
192
193 NULLIFY (admm_env, hfx_sections, mos, mos_aux_fit, para_env, virial, &
194 mo_coeff_aux_fit, xc_section, ks_env, dft_control, input, &
195 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 do_kpoints=do_kpoints)
208
209 calc_forces = .false.
210 IF (PRESENT(calculate_forces)) calc_forces = .true.
211
212 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
213 IF (PRESENT(ext_xc_section)) hfx_sections => section_vals_get_subs_vals(ext_xc_section, "HF")
214
215 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
216 IF (n_rep_hf > 1) &
217 cpabort("ADMM can handle only one HF section.")
218
219 IF (.NOT. ASSOCIATED(admm_env)) THEN
220 ! setup admm environment
221 CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
222 CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
223 CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
224 CALL set_qs_env(qs_env, admm_env=admm_env)
225 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
226 IF (PRESENT(ext_xc_section)) xc_section => ext_xc_section
227 CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
228 admm_env=admm_env)
229
230 ! Initialize the GAPW data types
231 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) &
232 CALL init_admm_gapw(qs_env)
233
234 ! ADMM neighbor lists and overlap matrices
235 CALL admm_init_hamiltonians(admm_env, qs_env, "AUX_FIT")
236
237 !The aux_fit task list and densities
238 ALLOCATE (admm_env%rho_aux_fit)
239 CALL qs_rho_create(admm_env%rho_aux_fit)
240 ALLOCATE (admm_env%rho_aux_fit_buffer)
241 CALL qs_rho_create(admm_env%rho_aux_fit_buffer)
242 CALL admm_update_s_mstruct(admm_env, qs_env, "AUX_FIT")
243 IF (admm_env%do_gapw) CALL update_admm_gapw(qs_env)
244
245 !The ADMM KS matrices
246 CALL admm_alloc_ks_matrices(admm_env, qs_env)
247
248 !The aux_fit MOs and derivatives
249 ALLOCATE (mos_aux_fit(dft_control%nspins))
250 DO ispin = 1, dft_control%nspins
251 CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, nelectron=nelectron, maxocc=maxocc)
252 CALL allocate_mo_set(mo_set=mos_aux_fit(ispin), &
253 nao=nao_aux_fit, &
254 nmo=nmo, &
255 nelectron=nelectron, &
256 n_el_f=real(nelectron, dp), &
257 maxocc=maxocc, &
258 flexible_electron_count=dft_control%relax_multiplicity)
259 END DO
260 admm_env%mos_aux_fit => mos_aux_fit
261
262 DO ispin = 1, dft_control%nspins
263 CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
264 CALL cp_fm_struct_create(aux_fit_fm_struct, context=blacs_env, para_env=para_env, &
265 nrow_global=nao_aux_fit, ncol_global=nmo)
266 CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, mo_coeff_b=mo_coeff_b)
267 IF (.NOT. ASSOCIATED(mo_coeff_aux_fit)) THEN
268 CALL init_mo_set(mos_aux_fit(ispin), fm_struct=aux_fit_fm_struct, &
269 name="qs_env%mo_aux_fit"//trim(adjustl(cp_to_string(ispin))))
270 END IF
271 CALL cp_fm_struct_release(aux_fit_fm_struct)
272
273 IF (.NOT. ASSOCIATED(mo_coeff_b)) THEN
274 CALL cp_fm_get_info(mos_aux_fit(ispin)%mo_coeff, ncol_global=nmo)
275 CALL dbcsr_init_p(mos_aux_fit(ispin)%mo_coeff_b)
276 CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit_kp)
277 CALL cp_dbcsr_m_by_n_from_row_template(mos_aux_fit(ispin)%mo_coeff_b, &
278 template=matrix_s_aux_fit_kp(1, 1)%matrix, &
279 n=nmo, sym=dbcsr_type_no_symmetry)
280 END IF
281 END DO
282
283 IF (qs_env%requires_mo_derivs) THEN
284 ALLOCATE (admm_env%mo_derivs_aux_fit(dft_control%nspins))
285 DO ispin = 1, dft_control%nspins
286 CALL get_mo_set(admm_env%mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit)
287 CALL cp_fm_create(admm_env%mo_derivs_aux_fit(ispin), mo_coeff_aux_fit%matrix_struct)
288 END DO
289 END IF
290
291 IF (do_kpoints) THEN
292 block
293 TYPE(kpoint_type), POINTER :: kpoints
294 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_aux_fit_kp
295 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools_aux_fit
296 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fm_struct
297 INTEGER :: ic, ik, ikk, is
298 INTEGER, PARAMETER :: nwork1 = 4
299 LOGICAL :: use_real_wfn
300
301 NULLIFY (ao_mo_fm_pools_aux_fit, mos_aux_fit_kp)
302
303 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
304 CALL get_kpoint_info(kpoints, use_real_wfn=use_real_wfn)
305
306 !Test combinations of input values. So far, only ADMM2 is availavle
307 IF (.NOT. admm_env%purification_method == do_admm_purify_none) &
308 cpabort("Only ADMM_PURIFICATION_METHOD NONE implemeted for ADMM K-points")
309 IF (.NOT. (dft_control%admm_control%method == do_admm_basis_projection &
310 .OR. dft_control%admm_control%method == do_admm_charge_constrained_projection)) &
311 cpabort("Only BASIS_PROJECTION and CHARGE_CONSTRAINED_PROJECTION implemented for KP")
312 IF (admm_env%do_admms .OR. admm_env%do_admmp .OR. admm_env%do_admmq) THEN
313 IF (use_real_wfn) cpabort("Only KP-HFX ADMM2 is implemented with REAL wavefunctions")
314 END IF
315
316 CALL kpoint_initialize_mos(kpoints, admm_env%mos_aux_fit, for_aux_fit=.true.)
317
318 CALL mpools_get(kpoints%mpools_aux_fit, ao_mo_fm_pools=ao_mo_fm_pools_aux_fit)
319 DO ik = 1, SIZE(kpoints%kp_aux_env)
320 mos_aux_fit_kp => kpoints%kp_aux_env(ik)%kpoint_env%mos
321 ikk = kpoints%kp_range(1) + ik - 1
322 DO ispin = 1, SIZE(mos_aux_fit_kp, 2)
323 DO ic = 1, SIZE(mos_aux_fit_kp, 1)
324 CALL get_mo_set(mos_aux_fit_kp(ic, ispin), mo_coeff=mo_coeff_aux_fit, mo_coeff_b=mo_coeff_b)
325
326 ! no sparse matrix representation of kpoint MO vectors
327 cpassert(.NOT. ASSOCIATED(mo_coeff_b))
328
329 IF (.NOT. ASSOCIATED(mo_coeff_aux_fit)) THEN
330 CALL init_mo_set(mos_aux_fit_kp(ic, ispin), &
331 fm_pool=ao_mo_fm_pools_aux_fit(ispin)%pool, &
332 name="kpoints_"//trim(adjustl(cp_to_string(ikk)))// &
333 "%mo_aux_fit"//trim(adjustl(cp_to_string(ispin))))
334 END IF
335 END DO
336 END DO
337 END DO
338
339 ALLOCATE (admm_env%scf_work_aux_fit(nwork1))
340
341 ! create an ao_ao parallel matrix structure
342 CALL cp_fm_struct_create(ao_ao_fm_struct, context=blacs_env, para_env=para_env, &
343 nrow_global=nao_aux_fit, &
344 ncol_global=nao_aux_fit)
345
346 DO is = 1, nwork1
347 CALL cp_fm_create(admm_env%scf_work_aux_fit(is), &
348 matrix_struct=ao_ao_fm_struct, &
349 name="SCF-WORK_MATRIX-AUX-"//trim(adjustl(cp_to_string(is))))
350 END DO
351 CALL cp_fm_struct_release(ao_ao_fm_struct)
352
353 ! Create and populate the internal ADMM overlap matrices at each KP
354 CALL kpoint_calc_admm_matrices(qs_env, calc_forces)
355
356 END block
357 END IF
358
359 ELSE IF (s_mstruct_changed) THEN
360 CALL admm_init_hamiltonians(admm_env, qs_env, "AUX_FIT")
361 CALL admm_update_s_mstruct(admm_env, qs_env, "AUX_FIT")
362 CALL admm_alloc_ks_matrices(admm_env, qs_env)
363 IF (admm_env%do_gapw) CALL update_admm_gapw(qs_env)
364 IF (do_kpoints) CALL kpoint_calc_admm_matrices(qs_env, calc_forces)
365 END IF
366
367 IF (admm_env%do_gapw .AND. dft_control%do_admm_dm) THEN
368 cpabort("GAPW ADMM not implemented for MCWEENY or NONE_DM purification.")
369 END IF
370
371 !ADMMS and ADMMP stress tensors only available for close-shell systesms, because virial cannot
372 !be scaled by gsi spin component wise
373 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
374 IF (use_virial .AND. admm_env%do_admms .AND. dft_control%nspins == 2) THEN
375 cpabort("ADMMS stress tensor is only available for closed-shell systems")
376 END IF
377 IF (use_virial .AND. admm_env%do_admmp .AND. dft_control%nspins == 2) THEN
378 cpabort("ADMMP stress tensor is only available for closed-shell systems")
379 END IF
380
381 IF (dft_control%do_admm_dm .AND. .NOT. ASSOCIATED(admm_env%admm_dm)) THEN
382 CALL admm_dm_create(admm_env%admm_dm, dft_control%admm_control, nspins=dft_control%nspins, natoms=natoms)
383 END IF
384
385 CALL timestop(handle)
386
387 END SUBROUTINE hfx_admm_init
388
389! **************************************************************************************************
390!> \brief Minimal setup routine for admm_env
391!> No forces
392!> No k-points
393!> No DFT correction terms
394!> \param qs_env ...
395!> \param mos ...
396!> \param admm_env ...
397!> \param admm_control ...
398!> \param basis_type ...
399! **************************************************************************************************
400 SUBROUTINE aux_admm_init(qs_env, mos, admm_env, admm_control, basis_type)
401
402 TYPE(qs_environment_type), POINTER :: qs_env
403 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
404 TYPE(admm_type), POINTER :: admm_env
405 TYPE(admm_control_type), POINTER :: admm_control
406 CHARACTER(LEN=*) :: basis_type
407
408 CHARACTER(LEN=*), PARAMETER :: routinen = 'aux_admm_init'
409
410 INTEGER :: handle, ispin, nao_aux_fit, natoms, &
411 nelectron, nmo
412 LOGICAL :: do_kpoints
413 REAL(dp) :: maxocc
414 TYPE(cp_blacs_env_type), POINTER :: blacs_env
415 TYPE(cp_fm_struct_type), POINTER :: aux_fit_fm_struct
416 TYPE(cp_fm_type), POINTER :: mo_coeff_aux_fit
417 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_aux_fit_kp
418 TYPE(dbcsr_type), POINTER :: mo_coeff_b
419 TYPE(dft_control_type), POINTER :: dft_control
420 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_aux_fit
421 TYPE(mp_para_env_type), POINTER :: para_env
422 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
423 TYPE(qs_ks_env_type), POINTER :: ks_env
424
425 CALL timeset(routinen, handle)
426
427 cpassert(.NOT. ASSOCIATED(admm_env))
428
429 CALL get_qs_env(qs_env, &
430 para_env=para_env, &
431 blacs_env=blacs_env, &
432 ks_env=ks_env, &
433 dft_control=dft_control, &
434 do_kpoints=do_kpoints)
435
436 cpassert(.NOT. do_kpoints)
437 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
438 cpabort("AUX ADMM not possible with GAPW")
439 END IF
440
441 ! setup admm environment
442 CALL get_qs_env(qs_env, natom=natoms, qs_kind_set=qs_kind_set)
443 CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type=basis_type)
444 !
445 CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
446 ! no XC correction used
447 NULLIFY (admm_env%xc_section_aux, admm_env%xc_section_primary)
448 ! ADMM neighbor lists and overlap matrices
449 CALL admm_init_hamiltonians(admm_env, qs_env, basis_type)
450 NULLIFY (admm_env%rho_aux_fit, admm_env%rho_aux_fit_buffer)
451 !The ADMM KS matrices
452 CALL admm_alloc_ks_matrices(admm_env, qs_env)
453 !The aux_fit MOs and derivatives
454 ALLOCATE (mos_aux_fit(dft_control%nspins))
455 DO ispin = 1, dft_control%nspins
456 CALL get_mo_set(mo_set=mos(ispin), nmo=nmo, nelectron=nelectron, maxocc=maxocc)
457 CALL allocate_mo_set(mo_set=mos_aux_fit(ispin), nao=nao_aux_fit, nmo=nmo, &
458 nelectron=nelectron, n_el_f=real(nelectron, dp), &
459 maxocc=maxocc, flexible_electron_count=0.0_dp)
460 END DO
461 admm_env%mos_aux_fit => mos_aux_fit
462
463 DO ispin = 1, dft_control%nspins
464 CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
465 CALL cp_fm_struct_create(aux_fit_fm_struct, context=blacs_env, para_env=para_env, &
466 nrow_global=nao_aux_fit, ncol_global=nmo)
467 CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, mo_coeff_b=mo_coeff_b)
468 IF (.NOT. ASSOCIATED(mo_coeff_aux_fit)) THEN
469 CALL init_mo_set(mos_aux_fit(ispin), fm_struct=aux_fit_fm_struct, &
470 name="mo_aux_fit"//trim(adjustl(cp_to_string(ispin))))
471 END IF
472 CALL cp_fm_struct_release(aux_fit_fm_struct)
473
474 IF (.NOT. ASSOCIATED(mo_coeff_b)) THEN
475 CALL cp_fm_get_info(mos_aux_fit(ispin)%mo_coeff, ncol_global=nmo)
476 CALL dbcsr_init_p(mos_aux_fit(ispin)%mo_coeff_b)
477 CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit_kp)
478 CALL cp_dbcsr_m_by_n_from_row_template(mos_aux_fit(ispin)%mo_coeff_b, &
479 template=matrix_s_aux_fit_kp(1, 1)%matrix, &
480 n=nmo, sym=dbcsr_type_no_symmetry)
481 END IF
482 END DO
483
484 CALL timestop(handle)
485
486 END SUBROUTINE aux_admm_init
487
488! **************************************************************************************************
489!> \brief Sets up the admm_gapw env
490!> \param qs_env ...
491! **************************************************************************************************
492 SUBROUTINE init_admm_gapw(qs_env)
493
494 TYPE(qs_environment_type), POINTER :: qs_env
495
496 INTEGER :: ikind, nkind
497 TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
498 TYPE(admm_type), POINTER :: admm_env
499 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
500 TYPE(dft_control_type), POINTER :: dft_control
501 TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, aux_fit_soft_basis, &
502 orb_basis, soft_basis
503 TYPE(mp_para_env_type), POINTER :: para_env
504 TYPE(qs_kind_type), DIMENSION(:), POINTER :: admm_kind_set, qs_kind_set
505 TYPE(section_vals_type), POINTER :: input
506
507 NULLIFY (admm_kind_set, aux_fit_basis, atomic_kind_set, aux_fit_soft_basis, &
508 dft_control, input, orb_basis, para_env, qs_kind_set, soft_basis)
509
510 CALL get_qs_env(qs_env, admm_env=admm_env, &
511 atomic_kind_set=atomic_kind_set, &
512 dft_control=dft_control, &
513 input=input, &
514 para_env=para_env, &
515 qs_kind_set=qs_kind_set)
516
517 admm_env%do_gapw = .true.
518 ALLOCATE (admm_env%admm_gapw_env)
519 admm_gapw_env => admm_env%admm_gapw_env
520 NULLIFY (admm_gapw_env%local_rho_set)
521 NULLIFY (admm_gapw_env%admm_kind_set)
522 NULLIFY (admm_gapw_env%task_list)
523
524 !Create a new kind set for the ADMM stuff (paw_proj soft AUX_FIT basis, etc)
525 nkind = SIZE(qs_kind_set)
526 ALLOCATE (admm_gapw_env%admm_kind_set(nkind))
527 admm_kind_set => admm_gapw_env%admm_kind_set
528
529 !In this new kind set, we want the AUX_FIT basis to be known as ORB, such that GAPW routines work
530 DO ikind = 1, nkind
531 !copying over simple data of interest from qs_kind_set
532 admm_kind_set(ikind)%name = qs_kind_set(ikind)%name
533 admm_kind_set(ikind)%element_symbol = qs_kind_set(ikind)%element_symbol
534 admm_kind_set(ikind)%natom = qs_kind_set(ikind)%natom
535 admm_kind_set(ikind)%hard_radius = qs_kind_set(ikind)%hard_radius
536 admm_kind_set(ikind)%max_rad_local = qs_kind_set(ikind)%max_rad_local
537 admm_kind_set(ikind)%gpw_type_forced = qs_kind_set(ikind)%gpw_type_forced
538 admm_kind_set(ikind)%ngrid_rad = qs_kind_set(ikind)%ngrid_rad
539 admm_kind_set(ikind)%ngrid_ang = qs_kind_set(ikind)%ngrid_ang
540
541 !copying potentials of interest from qs_kind_set
542 IF (ASSOCIATED(qs_kind_set(ikind)%all_potential)) THEN
543 CALL copy_potential(qs_kind_set(ikind)%all_potential, admm_kind_set(ikind)%all_potential)
544 END IF
545 IF (ASSOCIATED(qs_kind_set(ikind)%gth_potential)) THEN
546 CALL copy_potential(qs_kind_set(ikind)%gth_potential, admm_kind_set(ikind)%gth_potential)
547 END IF
548 IF (ASSOCIATED(qs_kind_set(ikind)%sgp_potential)) THEN
549 CALL copy_potential(qs_kind_set(ikind)%sgp_potential, admm_kind_set(ikind)%sgp_potential)
550 END IF
551
552 NULLIFY (orb_basis)
553 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis, basis_type="AUX_FIT")
554 CALL copy_gto_basis_set(aux_fit_basis, orb_basis)
555 CALL add_basis_set_to_container(admm_kind_set(ikind)%basis_sets, orb_basis, "ORB")
556 END DO
557
558 !Create the corresponding soft basis set (and projectors)
559 CALL init_gapw_basis_set(admm_kind_set, dft_control%qs_control, input, &
560 modify_qs_control=.false.)
561
562 !Make sure the basis and the projectors are well initialized
563 CALL init_interaction_radii(dft_control%qs_control, admm_kind_set)
564
565 !We also init the atomic grids and harmonics
566 CALL local_rho_set_create(admm_gapw_env%local_rho_set)
567 CALL init_rho_atom(admm_gapw_env%local_rho_set%rho_atom_set, &
568 atomic_kind_set, admm_kind_set, dft_control, para_env)
569
570 !Make sure that any NLCC potential is well initialized
571 CALL init_gapw_nlcc(admm_kind_set)
572
573 !Need to have access to the soft AUX_FIT basis from the qs_env => add it to the qs_kinds
574 DO ikind = 1, nkind
575 NULLIFY (aux_fit_soft_basis)
576 CALL get_qs_kind(admm_kind_set(ikind), basis_set=soft_basis, basis_type="ORB_SOFT")
577 CALL copy_gto_basis_set(soft_basis, aux_fit_soft_basis)
578 CALL add_basis_set_to_container(qs_kind_set(ikind)%basis_sets, aux_fit_soft_basis, "AUX_FIT_SOFT")
579 END DO
580
581 END SUBROUTINE init_admm_gapw
582
583! **************************************************************************************************
584!> \brief Builds the ADMM nmeighbor lists and overlap matrix on the model of qs_energies_init_hamiltonians()
585!> \param admm_env ...
586!> \param qs_env ...
587!> \param aux_basis_type ...
588! **************************************************************************************************
589 SUBROUTINE admm_init_hamiltonians(admm_env, qs_env, aux_basis_type)
590
591 TYPE(admm_type), POINTER :: admm_env
592 TYPE(qs_environment_type), POINTER :: qs_env
593 CHARACTER(len=*) :: aux_basis_type
594
595 CHARACTER(len=*), PARAMETER :: routinen = 'admm_init_hamiltonians'
596
597 INTEGER :: handle, hfx_pot, ikind, nkind
598 LOGICAL :: do_kpoints, mic, molecule_only
599 LOGICAL, ALLOCATABLE, DIMENSION(:) :: aux_fit_present, orb_present
600 REAL(dp) :: eps_schwarz, omega, pdist, roperator, &
601 subcells
602 REAL(dp), ALLOCATABLE, DIMENSION(:) :: aux_fit_radius, orb_radius
603 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
604 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
605 TYPE(cell_type), POINTER :: cell
606 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_aux_fit_kp, &
607 matrix_s_aux_fit_vs_orb_kp
608 TYPE(dft_control_type), POINTER :: dft_control
609 TYPE(distribution_1d_type), POINTER :: distribution_1d
610 TYPE(distribution_2d_type), POINTER :: distribution_2d
611 TYPE(gto_basis_set_type), POINTER :: aux_fit_basis_set, orb_basis_set
612 TYPE(kpoint_type), POINTER :: kpoints
613 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
614 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
615 TYPE(mp_para_env_type), POINTER :: para_env
616 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
617 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
618 TYPE(qs_ks_env_type), POINTER :: ks_env
619 TYPE(section_vals_type), POINTER :: hfx_sections, neighbor_list_section
620
621 NULLIFY (particle_set, cell, kpoints, distribution_1d, distribution_2d, molecule_set, &
622 atomic_kind_set, dft_control, neighbor_list_section, aux_fit_basis_set, orb_basis_set, &
623 ks_env, para_env, qs_kind_set, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb_kp)
624
625 CALL timeset(routinen, handle)
626
627 CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set, cell=cell, kpoints=kpoints, &
628 local_particles=distribution_1d, distribution_2d=distribution_2d, &
629 molecule_set=molecule_set, atomic_kind_set=atomic_kind_set, do_kpoints=do_kpoints, &
630 dft_control=dft_control, para_env=para_env, qs_kind_set=qs_kind_set)
631 ALLOCATE (orb_present(nkind), aux_fit_present(nkind))
632 ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind), pair_radius(nkind, nkind))
633 aux_fit_radius(:) = 0.0_dp
634
635 molecule_only = .false.
636 IF (dft_control%qs_control%do_kg) molecule_only = .true.
637 mic = molecule_only
638 IF (kpoints%nkp > 0) THEN
639 mic = .false.
640 ELSEIF (dft_control%qs_control%semi_empirical) THEN
641 mic = .true.
642 END IF
643
644 pdist = dft_control%qs_control%pairlist_radius
645
646 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
647 neighbor_list_section => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%NEIGHBOR_LISTS")
648
649 ALLOCATE (atom2d(nkind))
650 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
651 molecule_set, molecule_only, particle_set=particle_set)
652
653 DO ikind = 1, nkind
654 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
655 IF (ASSOCIATED(orb_basis_set)) THEN
656 orb_present(ikind) = .true.
657 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
658 ELSE
659 orb_present(ikind) = .false.
660 END IF
661
662 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type=aux_basis_type)
663 IF (ASSOCIATED(aux_fit_basis_set)) THEN
664 aux_fit_present(ikind) = .true.
665 CALL get_gto_basis_set(gto_basis_set=aux_fit_basis_set, kind_radius=aux_fit_radius(ikind))
666 ELSE
667 aux_fit_present(ikind) = .false.
668 END IF
669 END DO
670
671 IF (pdist < 0.0_dp) THEN
672 pdist = max(plane_distance(1, 0, 0, cell), &
673 plane_distance(0, 1, 0, cell), &
674 plane_distance(0, 0, 1, cell))
675 END IF
676
677 !In case of K-points, we need to add the HFX potential range to sab_aux_fit, because it is used
678 !to populate AUX density and KS matrices
679 roperator = 0.0_dp
680 IF (do_kpoints) THEN
681 hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
682 CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", i_val=hfx_pot)
683
684 SELECT CASE (hfx_pot)
685 CASE (do_potential_id)
686 roperator = 0.0_dp
687 CASE (do_potential_truncated)
688 CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
689 CASE (do_potential_mix_cl_trunc)
690 CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
691 CASE (do_potential_short)
692 CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%OMEGA", r_val=omega)
693 CALL section_vals_val_get(hfx_sections, "SCREENING%EPS_SCHWARZ", r_val=eps_schwarz)
694 CALL erfc_cutoff(eps_schwarz, omega, roperator)
695 CASE DEFAULT
696 cpabort("HFX potential not available for K-points (NYI)")
697 END SELECT
698 END IF
699
700 CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, pair_radius, pdist)
701 pair_radius = pair_radius + cutoff_screen_factor*roperator
702 CALL build_neighbor_lists(admm_env%sab_aux_fit, particle_set, atom2d, cell, pair_radius, &
703 mic=mic, molecular=molecule_only, subcells=subcells, nlname="sab_aux_fit")
704 CALL build_neighbor_lists(admm_env%sab_aux_fit_asymm, particle_set, atom2d, cell, pair_radius, &
705 mic=mic, symmetric=.false., molecular=molecule_only, subcells=subcells, &
706 nlname="sab_aux_fit_asymm")
707 CALL pair_radius_setup(aux_fit_present, orb_present, aux_fit_radius, orb_radius, pair_radius)
708 CALL build_neighbor_lists(admm_env%sab_aux_fit_vs_orb, particle_set, atom2d, cell, pair_radius, &
709 mic=mic, symmetric=.false., molecular=molecule_only, subcells=subcells, &
710 nlname="sab_aux_fit_vs_orb")
711
712 CALL write_neighbor_lists(admm_env%sab_aux_fit, particle_set, cell, para_env, neighbor_list_section, &
713 "/SAB_AUX_FIT", "sab_aux_fit", "AUX_FIT_ORBITAL AUX_FIT_ORBITAL")
714 CALL write_neighbor_lists(admm_env%sab_aux_fit_vs_orb, particle_set, cell, para_env, neighbor_list_section, &
715 "/SAB_AUX_FIT_VS_ORB", "sab_aux_fit_vs_orb", "ORBITAL AUX_FIT_ORBITAL")
716
717 CALL atom2d_cleanup(atom2d)
718
719 !The ADMM overlap matrices (initially in qs_core_hamiltonian.F)
720 CALL get_qs_env(qs_env, ks_env=ks_env)
721
722 CALL kpoint_transitional_release(admm_env%matrix_s_aux_fit)
723 CALL build_overlap_matrix(ks_env, matrixkp_s=matrix_s_aux_fit_kp, &
724 matrix_name="AUX_FIT_OVERLAP", &
725 basis_type_a=aux_basis_type, &
726 basis_type_b=aux_basis_type, &
727 sab_nl=admm_env%sab_aux_fit)
728 CALL set_2d_pointer(admm_env%matrix_s_aux_fit, matrix_s_aux_fit_kp)
729 CALL kpoint_transitional_release(admm_env%matrix_s_aux_fit_vs_orb)
730 CALL build_overlap_matrix(ks_env, matrixkp_s=matrix_s_aux_fit_vs_orb_kp, &
731 matrix_name="MIXED_OVERLAP", &
732 basis_type_a=aux_basis_type, &
733 basis_type_b="ORB", &
734 sab_nl=admm_env%sab_aux_fit_vs_orb)
735 CALL set_2d_pointer(admm_env%matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp)
736
737 CALL timestop(handle)
738
739 END SUBROUTINE admm_init_hamiltonians
740
741! **************************************************************************************************
742!> \brief Updates the ADMM task_list and density based on the model of qs_env_update_s_mstruct()
743!> \param admm_env ...
744!> \param qs_env ...
745!> \param aux_basis_type ...
746! **************************************************************************************************
747 SUBROUTINE admm_update_s_mstruct(admm_env, qs_env, aux_basis_type)
748
749 TYPE(admm_type), POINTER :: admm_env
750 TYPE(qs_environment_type), POINTER :: qs_env
751 CHARACTER(len=*) :: aux_basis_type
752
753 CHARACTER(len=*), PARAMETER :: routinen = 'admm_update_s_mstruct'
754
755 INTEGER :: handle
756 LOGICAL :: skip_load_balance_distributed
757 TYPE(dft_control_type), POINTER :: dft_control
758 TYPE(qs_ks_env_type), POINTER :: ks_env
759
760 NULLIFY (ks_env, dft_control)
761
762 CALL timeset(routinen, handle)
763
764 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control)
765
766 !The aux_fit task_list
767 skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
768 IF (ASSOCIATED(admm_env%task_list_aux_fit)) CALL deallocate_task_list(admm_env%task_list_aux_fit)
769 CALL allocate_task_list(admm_env%task_list_aux_fit)
770 CALL generate_qs_task_list(ks_env, admm_env%task_list_aux_fit, basis_type=aux_basis_type, &
771 reorder_rs_grid_ranks=.false., &
772 skip_load_balance_distributed=skip_load_balance_distributed, &
773 sab_orb_external=admm_env%sab_aux_fit)
774
775 !The aux_fit densities
776 CALL qs_rho_rebuild(admm_env%rho_aux_fit, qs_env=qs_env, admm=.true.)
777 CALL qs_rho_rebuild(admm_env%rho_aux_fit_buffer, qs_env=qs_env, admm=.true.)
778
779 CALL timestop(handle)
780
781 END SUBROUTINE admm_update_s_mstruct
782
783! **************************************************************************************************
784!> \brief Update the admm_gapw_env internals to the current qs_env (i.e. atomic positions)
785!> \param qs_env ...
786! **************************************************************************************************
787 SUBROUTINE update_admm_gapw(qs_env)
788
789 TYPE(qs_environment_type), POINTER :: qs_env
790
791 CHARACTER(len=*), PARAMETER :: routinen = 'update_admm_gapw'
792
793 INTEGER :: handle, ikind, nkind
794 LOGICAL :: paw_atom
795 LOGICAL, ALLOCATABLE, DIMENSION(:) :: aux_present, oce_present
796 REAL(dp) :: subcells
797 REAL(dp), ALLOCATABLE, DIMENSION(:) :: aux_radius, oce_radius
798 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
799 TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
800 TYPE(admm_type), POINTER :: admm_env
801 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
802 TYPE(cell_type), POINTER :: cell
803 TYPE(dft_control_type), POINTER :: dft_control
804 TYPE(distribution_1d_type), POINTER :: distribution_1d
805 TYPE(distribution_2d_type), POINTER :: distribution_2d
806 TYPE(gto_basis_set_type), POINTER :: aux_fit_basis
807 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
808 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
809 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
810 POINTER :: sap_oce
811 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
812 TYPE(paw_proj_set_type), POINTER :: paw_proj
813 TYPE(qs_kind_type), DIMENSION(:), POINTER :: admm_kind_set, qs_kind_set
814 TYPE(qs_ks_env_type), POINTER :: ks_env
815
816 NULLIFY (ks_env, qs_kind_set, admm_kind_set, aux_fit_basis, cell, distribution_1d)
817 NULLIFY (distribution_2d, paw_proj, particle_set, molecule_set, admm_env, admm_gapw_env)
818 NULLIFY (dft_control, atomic_kind_set, sap_oce)
819
820 CALL timeset(routinen, handle)
821
822 CALL get_qs_env(qs_env, ks_env=ks_env, qs_kind_set=qs_kind_set, admm_env=admm_env, &
823 dft_control=dft_control)
824 admm_gapw_env => admm_env%admm_gapw_env
825 admm_kind_set => admm_gapw_env%admm_kind_set
826 nkind = SIZE(qs_kind_set)
827
828 !Update the task lisft for the AUX_FIT_SOFT basis
829 IF (ASSOCIATED(admm_gapw_env%task_list)) CALL deallocate_task_list(admm_gapw_env%task_list)
830 CALL allocate_task_list(admm_gapw_env%task_list)
831
832 !note: we set soft_valid to .FALSE. want to use AUX_FIT_SOFT and not the normal ORB SOFT basis
833 CALL generate_qs_task_list(ks_env, admm_gapw_env%task_list, basis_type="AUX_FIT_SOFT", &
834 reorder_rs_grid_ranks=.false., &
835 skip_load_balance_distributed=dft_control%qs_control%skip_load_balance_distributed, &
836 sab_orb_external=admm_env%sab_aux_fit)
837
838 !Update the precomputed oce integrals
839 !a sap_oce neighbor list is required => build it here
840 ALLOCATE (aux_present(nkind), oce_present(nkind))
841 aux_present = .false.; oce_present = .false.
842 ALLOCATE (aux_radius(nkind), oce_radius(nkind))
843 aux_radius = 0.0_dp; oce_radius = 0.0_dp
844
845 DO ikind = 1, nkind
846 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis, basis_type="AUX_FIT")
847 IF (ASSOCIATED(aux_fit_basis)) THEN
848 aux_present(ikind) = .true.
849 CALL get_gto_basis_set(aux_fit_basis, kind_radius=aux_radius(ikind))
850 END IF
851
852 !note: get oce info from admm_kind_set
853 CALL get_qs_kind(admm_kind_set(ikind), paw_atom=paw_atom, paw_proj_set=paw_proj)
854 IF (paw_atom) THEN
855 oce_present(ikind) = .true.
856 CALL get_paw_proj_set(paw_proj, rcprj=oce_radius(ikind))
857 END IF
858 END DO
859
860 ALLOCATE (pair_radius(nkind, nkind))
861 pair_radius = 0.0_dp
862 CALL pair_radius_setup(aux_present, oce_present, aux_radius, oce_radius, pair_radius)
863
864 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
865 distribution_2d=distribution_2d, local_particles=distribution_1d, &
866 particle_set=particle_set, molecule_set=molecule_set)
867 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
868
869 ALLOCATE (atom2d(nkind))
870 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
871 molecule_set, .false., particle_set)
872 CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, &
873 subcells=subcells, operator_type="ABBA", nlname="AUX_PAW-PRJ")
874 CALL atom2d_cleanup(atom2d)
875
876 !actually compute the oce matrices
877 CALL create_oce_set(admm_gapw_env%oce)
878 CALL allocate_oce_set(admm_gapw_env%oce, nkind)
879
880 !always compute the derivative, cheap anyways
881 CALL build_oce_matrices(admm_gapw_env%oce%intac, calculate_forces=.true., nder=1, &
882 qs_kind_set=admm_kind_set, particle_set=particle_set, &
883 sap_oce=sap_oce, eps_fit=dft_control%qs_control%gapw_control%eps_fit)
884
885 CALL release_neighbor_list_sets(sap_oce)
886
887 CALL timestop(handle)
888
889 END SUBROUTINE update_admm_gapw
890
891! **************************************************************************************************
892!> \brief Allocates the various ADMM KS matrices
893!> \param admm_env ...
894!> \param qs_env ...
895! **************************************************************************************************
896 SUBROUTINE admm_alloc_ks_matrices(admm_env, qs_env)
897
898 TYPE(admm_type), POINTER :: admm_env
899 TYPE(qs_environment_type), POINTER :: qs_env
900
901 CHARACTER(len=*), PARAMETER :: routinen = 'admm_alloc_ks_matrices'
902
903 INTEGER :: handle, ic, ispin
904 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit_dft_kp, &
905 matrix_ks_aux_fit_hfx_kp, &
906 matrix_ks_aux_fit_kp, &
907 matrix_s_aux_fit_kp
908 TYPE(dft_control_type), POINTER :: dft_control
909
910 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)
911
912 CALL timeset(routinen, handle)
913
914 CALL get_qs_env(qs_env, dft_control=dft_control)
915 CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit_kp)
916
917 CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit)
918 CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit_dft)
919 CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit_hfx)
920
921 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_kp, dft_control%nspins, dft_control%nimages)
922 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_dft_kp, dft_control%nspins, dft_control%nimages)
923 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_hfx_kp, dft_control%nspins, dft_control%nimages)
924
925 DO ispin = 1, dft_control%nspins
926 DO ic = 1, dft_control%nimages
927 ALLOCATE (matrix_ks_aux_fit_kp(ispin, ic)%matrix)
928 CALL dbcsr_create(matrix_ks_aux_fit_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, ic)%matrix, &
929 name="KOHN-SHAM_MATRIX for ADMM")
930 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
931 CALL dbcsr_set(matrix_ks_aux_fit_kp(ispin, ic)%matrix, 0.0_dp)
932
933 ALLOCATE (matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix)
934 CALL dbcsr_create(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, 1)%matrix, &
935 name="KOHN-SHAM_MATRIX for ADMM")
936 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
937 CALL dbcsr_set(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, 0.0_dp)
938
939 ALLOCATE (matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix)
940 CALL dbcsr_create(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, 1)%matrix, &
941 name="KOHN-SHAM_MATRIX for ADMM")
942 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
943 CALL dbcsr_set(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, 0.0_dp)
944 END DO
945 END DO
946
947 CALL set_admm_env(admm_env, &
948 matrix_ks_aux_fit_kp=matrix_ks_aux_fit_kp, &
949 matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft_kp, &
950 matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx_kp)
951
952 CALL timestop(handle)
953
954 END SUBROUTINE admm_alloc_ks_matrices
955
956! **************************************************************************************************
957!> \brief Add the HFX K-point contribution to the real-space Hamiltonians
958!> \param qs_env ...
959!> \param matrix_ks ...
960!> \param energy ...
961!> \param calculate_forces ...
962! **************************************************************************************************
963 SUBROUTINE hfx_ks_matrix_kp(qs_env, matrix_ks, energy, calculate_forces)
964 TYPE(qs_environment_type), POINTER :: qs_env
965 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
966 TYPE(qs_energy_type), POINTER :: energy
967 LOGICAL, INTENT(in) :: calculate_forces
968
969 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ks_matrix_kp'
970
971 INTEGER :: handle, img, irep, ispin, n_rep_hf, &
972 nimages, nspins
973 LOGICAL :: do_adiabatic_rescaling, &
974 s_mstruct_changed, use_virial
975 REAL(dp) :: eh1, ehfx, eold
976 REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy
977 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit_im, matrix_ks_im
978 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks_aux_fit_hfx_kp, &
979 matrix_ks_aux_fit_kp, matrix_ks_orb, &
980 rho_ao_orb
981 TYPE(dft_control_type), POINTER :: dft_control
982 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
983 TYPE(mp_para_env_type), POINTER :: para_env
984 TYPE(pw_env_type), POINTER :: pw_env
985 TYPE(pw_poisson_type), POINTER :: poisson_env
986 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
987 TYPE(qs_rho_type), POINTER :: rho_orb
988 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
989 hfx_sections, input
990 TYPE(virial_type), POINTER :: virial
991
992 CALL timeset(routinen, handle)
993
994 NULLIFY (auxbas_pw_pool, dft_control, hfx_sections, input, &
995 para_env, poisson_env, pw_env, virial, matrix_ks_im, &
996 matrix_ks_orb, rho_ao_orb, matrix_h, matrix_ks_aux_fit_kp, &
997 matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx_kp)
998
999 CALL get_qs_env(qs_env=qs_env, &
1000 dft_control=dft_control, &
1001 input=input, &
1002 matrix_h_kp=matrix_h, &
1003 para_env=para_env, &
1004 pw_env=pw_env, &
1005 virial=virial, &
1006 matrix_ks_im=matrix_ks_im, &
1007 s_mstruct_changed=s_mstruct_changed, &
1008 x_data=x_data)
1009
1010 ! No RTP
1011 IF (qs_env%run_rtp) cpabort("No RTP implementation with K-points HFX")
1012
1013 ! No adiabatic rescaling
1014 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
1015 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1016 IF (do_adiabatic_rescaling) cpabort("No adiabatic rescaling implementation with K-points HFX")
1017
1018 IF (dft_control%do_admm) THEN
1019 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit_kp, &
1020 matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, &
1021 matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx_kp)
1022 END IF
1023
1024 nspins = dft_control%nspins
1025 nimages = dft_control%nimages
1026
1027 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1028 IF (use_virial .AND. calculate_forces) virial%pv_fock_4c = 0.0_dp
1029
1030 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1031 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1032
1033 ! *** Initialize the auxiliary ks matrix to zero if required
1034 IF (dft_control%do_admm) THEN
1035 DO ispin = 1, nspins
1036 DO img = 1, nimages
1037 CALL dbcsr_set(matrix_ks_aux_fit_kp(ispin, img)%matrix, 0.0_dp)
1038 END DO
1039 END DO
1040 END IF
1041 DO ispin = 1, nspins
1042 DO img = 1, nimages
1043 CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
1044 END DO
1045 END DO
1046
1047 ALLOCATE (hf_energy(n_rep_hf))
1048
1049 eold = 0.0_dp
1050
1051 DO irep = 1, n_rep_hf
1052
1053 ! fetch the correct matrices for normal HFX or ADMM
1054 IF (dft_control%do_admm) THEN
1055 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks_orb, rho_aux_fit=rho_orb)
1056 ELSE
1057 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_orb, rho=rho_orb)
1058 END IF
1059 CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb)
1060
1061 ! Finally the real hfx calulation
1062 ehfx = 0.0_dp
1063
1064 IF (.NOT. x_data(irep, 1)%do_hfx_ri) THEN
1065 cpabort("Only RI-HFX is implemented for K-points")
1066 END IF
1067
1068 CALL hfx_ri_update_ks_kp(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1069 rho_ao_orb, s_mstruct_changed, nspins, &
1070 x_data(irep, 1)%general_parameter%fraction)
1071
1072 IF (calculate_forces) THEN
1073 !Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force
1074 IF (dft_control%do_admm) THEN
1075 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.false.)
1076 END IF
1077
1078 CALL hfx_ri_update_forces_kp(qs_env, x_data(irep, 1)%ri_data, nspins, &
1079 x_data(irep, 1)%general_parameter%fraction, &
1080 rho_ao_orb, use_virial=use_virial)
1081
1082 IF (dft_control%do_admm) THEN
1083 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.true.)
1084 END IF
1085 END IF
1086
1087 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
1088 eh1 = ehfx - eold
1089 CALL pw_hfx(qs_env, eh1, hfx_sections, poisson_env, auxbas_pw_pool, irep)
1090 eold = ehfx
1091
1092 END DO
1093
1094 ! *** Set the total HFX energy
1095 energy%ex = ehfx
1096
1097 ! *** Add Core-Hamiltonian-Matrix ***
1098 DO ispin = 1, nspins
1099 DO img = 1, nimages
1100 CALL dbcsr_add(matrix_ks(ispin, img)%matrix, matrix_h(1, img)%matrix, &
1101 1.0_dp, 1.0_dp)
1102 END DO
1103 END DO
1104 IF (use_virial .AND. calculate_forces) THEN
1105 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1106 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1107 virial%pv_calculate = .false.
1108 END IF
1109
1110 !update the hfx aux_fit matrix
1111 IF (dft_control%do_admm) THEN
1112 DO ispin = 1, nspins
1113 DO img = 1, nimages
1114 CALL dbcsr_add(matrix_ks_aux_fit_hfx_kp(ispin, img)%matrix, matrix_ks_aux_fit_kp(ispin, img)%matrix, &
1115 0.0_dp, 1.0_dp)
1116 END DO
1117 END DO
1118 END IF
1119
1120 CALL timestop(handle)
1121
1122 END SUBROUTINE hfx_ks_matrix_kp
1123
1124! **************************************************************************************************
1125!> \brief Add the hfx contributions to the Hamiltonian
1126!>
1127!> \param qs_env ...
1128!> \param matrix_ks ...
1129!> \param rho ...
1130!> \param energy ...
1131!> \param calculate_forces ...
1132!> \param just_energy ...
1133!> \param v_rspace_new ...
1134!> \param v_tau_rspace ...
1135!> \param ext_xc_section ...
1136!> \par History
1137!> refactoring 03-2011 [MI]
1138! **************************************************************************************************
1139
1140 SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
1141 just_energy, v_rspace_new, v_tau_rspace, ext_xc_section)
1142
1143 TYPE(qs_environment_type), POINTER :: qs_env
1144 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
1145 TYPE(qs_rho_type), POINTER :: rho
1146 TYPE(qs_energy_type), POINTER :: energy
1147 LOGICAL, INTENT(in) :: calculate_forces, just_energy
1148 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
1149 TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
1150
1151 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ks_matrix'
1152
1153 INTEGER :: handle, img, irep, ispin, mspin, &
1154 n_rep_hf, nimages, ns, nspins
1155 LOGICAL :: distribute_fock_matrix, &
1156 do_adiabatic_rescaling, &
1157 hfx_treat_lsd_in_core, &
1158 s_mstruct_changed, use_virial
1159 REAL(dp) :: eh1, ehfx, ehfxrt, eold
1160 REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy
1161 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_1d, matrix_ks_aux_fit, &
1162 matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_im, matrix_ks_im, rho_ao_1d, rho_ao_resp
1163 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im, matrix_ks_orb, &
1164 rho_ao_orb
1165 TYPE(dft_control_type), POINTER :: dft_control
1166 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1167 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1168 TYPE(mp_para_env_type), POINTER :: para_env
1169 TYPE(pw_env_type), POINTER :: pw_env
1170 TYPE(pw_poisson_type), POINTER :: poisson_env
1171 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1172 TYPE(qs_rho_type), POINTER :: rho_orb
1173 TYPE(rt_prop_type), POINTER :: rtp
1174 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
1175 hfx_sections, input
1176 TYPE(virial_type), POINTER :: virial
1177
1178 CALL timeset(routinen, handle)
1179
1180 NULLIFY (auxbas_pw_pool, dft_control, hfx_sections, input, &
1181 para_env, poisson_env, pw_env, virial, matrix_ks_im, &
1182 matrix_ks_orb, rho_ao_orb, matrix_h, matrix_h_im, matrix_ks_aux_fit, &
1183 matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx)
1184
1185 CALL get_qs_env(qs_env=qs_env, &
1186 dft_control=dft_control, &
1187 input=input, &
1188 matrix_h_kp=matrix_h, &
1189 matrix_h_im_kp=matrix_h_im, &
1190 para_env=para_env, &
1191 pw_env=pw_env, &
1192 virial=virial, &
1193 matrix_ks_im=matrix_ks_im, &
1194 s_mstruct_changed=s_mstruct_changed, &
1195 x_data=x_data)
1196
1197 IF (dft_control%do_admm) THEN
1198 CALL get_admm_env(qs_env%admm_env, mos_aux_fit=mo_array, matrix_ks_aux_fit=matrix_ks_aux_fit, &
1199 matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
1200 ELSE
1201 CALL get_qs_env(qs_env=qs_env, mos=mo_array)
1202 END IF
1203
1204 nspins = dft_control%nspins
1205 nimages = dft_control%nimages
1206
1207 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1208
1209 IF (use_virial .AND. calculate_forces) virial%pv_fock_4c = 0.0_dp
1210
1211 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1212 IF (PRESENT(ext_xc_section)) hfx_sections => section_vals_get_subs_vals(ext_xc_section, "HF")
1213
1214 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1215 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
1216 i_rep_section=1)
1217 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
1218 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1219
1220 ! *** Initialize the auxiliary ks matrix to zero if required
1221 IF (dft_control%do_admm) THEN
1222 DO ispin = 1, nspins
1223 CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
1224 END DO
1225 END IF
1226 DO ispin = 1, nspins
1227 DO img = 1, nimages
1228 CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
1229 END DO
1230 END DO
1231
1232 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1233
1234 ALLOCATE (hf_energy(n_rep_hf))
1235
1236 eold = 0.0_dp
1237
1238 DO irep = 1, n_rep_hf
1239 ! Remember: Vhfx is added, energy is calclulated from total Vhfx,
1240 ! so energy of last iteration is correct
1241
1242 IF (do_adiabatic_rescaling .AND. hfx_treat_lsd_in_core) &
1243 cpabort("HFX_TREAT_LSD_IN_CORE not implemented for adiabatically rescaled hybrids")
1244 ! everything is calculated with adiabatic rescaling but the potential is not added in a first step
1245 distribute_fock_matrix = .NOT. do_adiabatic_rescaling
1246
1247 mspin = 1
1248 IF (hfx_treat_lsd_in_core) mspin = nspins
1249
1250 ! fetch the correct matrices for normal HFX or ADMM
1251 IF (dft_control%do_admm) THEN
1252 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_1d, rho_aux_fit=rho_orb)
1253 ns = SIZE(matrix_ks_1d)
1254 matrix_ks_orb(1:ns, 1:1) => matrix_ks_1d(1:ns)
1255 ELSE
1256 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_orb, rho=rho_orb)
1257 END IF
1258 CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb)
1259 ! Finally the real hfx calulation
1260 ehfx = 0.0_dp
1261
1262 IF (x_data(irep, 1)%do_hfx_ri) THEN
1263
1264 CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1265 mo_array, rho_ao_orb, &
1266 s_mstruct_changed, nspins, &
1267 x_data(irep, 1)%general_parameter%fraction)
1268 IF (dft_control%do_admm) THEN
1269 !for ADMMS, we need the exchange matrix k(d) for both spins
1270 DO ispin = 1, nspins
1271 CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
1272 name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1273 END DO
1274 END IF
1275
1276 ELSE
1277
1278 DO ispin = 1, mspin
1279 CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
1280 para_env, s_mstruct_changed, irep, distribute_fock_matrix, &
1281 ispin=ispin)
1282 ehfx = ehfx + eh1
1283 END DO
1284 END IF
1285
1286 IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
1287 !Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force
1288 IF (dft_control%do_admm) THEN
1289 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.false.)
1290 END IF
1291 NULLIFY (rho_ao_resp)
1292
1293 IF (x_data(irep, 1)%do_hfx_ri) THEN
1294
1295 CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
1296 x_data(irep, 1)%general_parameter%fraction, &
1297 rho_ao=rho_ao_orb, mos=mo_array, &
1298 rho_ao_resp=rho_ao_resp, &
1299 use_virial=use_virial)
1300
1301 ELSE
1302
1303 CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
1304 para_env, irep, use_virial)
1305
1306 END IF
1307
1308 !Scale auxiliary density matrix for ADMMP back with 1/gsi(ispin)
1309 IF (dft_control%do_admm) THEN
1310 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.true.)
1311 END IF
1312 END IF
1313
1314 !! If required, the calculation of the forces will be done later with adiabatic rescaling
1315 IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx
1316
1317 ! special case RTP/EMD we have a full complex density and HFX has a contribution from the imaginary part
1318 ehfxrt = 0.0_dp
1319 IF (qs_env%run_rtp) THEN
1320
1321 CALL get_qs_env(qs_env=qs_env, rtp=rtp)
1322 DO ispin = 1, nspins
1323 CALL dbcsr_set(matrix_ks_im(ispin)%matrix, 0.0_dp)
1324 END DO
1325 IF (dft_control%do_admm) THEN
1326 ! matrix_ks_orb => matrix_ks_aux_fit_im
1327 ns = SIZE(matrix_ks_aux_fit_im)
1328 matrix_ks_orb(1:ns, 1:1) => matrix_ks_aux_fit_im(1:ns)
1329 DO ispin = 1, nspins
1330 CALL dbcsr_set(matrix_ks_aux_fit_im(ispin)%matrix, 0.0_dp)
1331 END DO
1332 ELSE
1333 ! matrix_ks_orb => matrix_ks_im
1334 ns = SIZE(matrix_ks_im)
1335 matrix_ks_orb(1:ns, 1:1) => matrix_ks_im(1:ns)
1336 END IF
1337
1338 CALL qs_rho_get(rho_orb, rho_ao_im=rho_ao_1d)
1339 ns = SIZE(rho_ao_1d)
1340 rho_ao_orb(1:ns, 1:1) => rho_ao_1d(1:ns)
1341
1342 ehfxrt = 0.0_dp
1343
1344 IF (x_data(irep, 1)%do_hfx_ri) THEN
1345 CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1346 mo_array, rho_ao_orb, &
1347 .false., nspins, &
1348 x_data(irep, 1)%general_parameter%fraction)
1349 IF (dft_control%do_admm) THEN
1350 !for ADMMS, we need the exchange matrix k(d) for both spins
1351 DO ispin = 1, nspins
1352 CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
1353 name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1354 END DO
1355 END IF
1356
1357 ELSE
1358 DO ispin = 1, mspin
1359 CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
1360 para_env, .false., irep, distribute_fock_matrix, &
1361 ispin=ispin)
1362 ehfxrt = ehfxrt + eh1
1363 END DO
1364 END IF
1365
1366 IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
1367 NULLIFY (rho_ao_resp)
1368
1369 IF (x_data(irep, 1)%do_hfx_ri) THEN
1370
1371 CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
1372 x_data(irep, 1)%general_parameter%fraction, &
1373 rho_ao=rho_ao_orb, mos=mo_array, &
1374 use_virial=use_virial)
1375
1376 ELSE
1377 CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
1378 para_env, irep, use_virial)
1379 END IF
1380 END IF
1381
1382 !! If required, the calculation of the forces will be done later with adiabatic rescaling
1383 IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx + ehfxrt
1384
1385 IF (dft_control%rtp_control%velocity_gauge) THEN
1386 cpassert(ASSOCIATED(matrix_h_im))
1387 DO ispin = 1, nspins
1388 CALL dbcsr_add(matrix_ks_im(ispin)%matrix, matrix_h_im(1, 1)%matrix, &
1389 1.0_dp, 1.0_dp)
1390 END DO
1391 END IF
1392
1393 END IF
1394
1395 IF (.NOT. qs_env%run_rtp) THEN
1396 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1397 poisson_env=poisson_env)
1398 eh1 = ehfx - eold
1399 CALL pw_hfx(qs_env, eh1, hfx_sections, poisson_env, auxbas_pw_pool, irep)
1400 eold = ehfx
1401 END IF
1402
1403 END DO
1404
1405 ! *** Set the total HFX energy
1406 energy%ex = ehfx + ehfxrt
1407
1408 ! *** Add Core-Hamiltonian-Matrix ***
1409 DO ispin = 1, nspins
1410 DO img = 1, nimages
1411 CALL dbcsr_add(matrix_ks(ispin, img)%matrix, matrix_h(1, img)%matrix, &
1412 1.0_dp, 1.0_dp)
1413 END DO
1414 END DO
1415 IF (use_virial .AND. calculate_forces) THEN
1416 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1417 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1418 virial%pv_calculate = .false.
1419 END IF
1420
1421 !! If we perform adiabatic rescaling we are now able to rescale the xc-potential
1422 IF (do_adiabatic_rescaling) THEN
1423 CALL rescale_xc_potential(qs_env, matrix_ks, rho, energy, v_rspace_new, v_tau_rspace, &
1424 hf_energy, just_energy, calculate_forces, use_virial)
1425 END IF ! do_adiabatic_rescaling
1426
1427 !update the hfx aux_fit matrixIF (dft_control%do_admm) THEN
1428 IF (dft_control%do_admm) THEN
1429 DO ispin = 1, nspins
1430 CALL dbcsr_add(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, &
1431 0.0_dp, 1.0_dp)
1432 END DO
1433 END IF
1434
1435 CALL timestop(handle)
1436
1437 END SUBROUTINE hfx_ks_matrix
1438
1439! **************************************************************************************************
1440!> \brief This routine modifies the xc section depending on the potential type
1441!> used for the HF exchange and the resulting correction term. Currently
1442!> three types of corrections are implemented:
1443!>
1444!> coulomb: Ex,hf = Ex,hf' + (PBEx-PBEx')
1445!> shortrange: Ex,hf = Ex,hf' + (XWPBEX-XWPBEX')
1446!> truncated: Ex,hf = Ex,hf' + ( (XWPBEX0-PBE_HOLE_TC_LR) -(XWPBEX0-PBE_HOLE_TC_LR)' )
1447!>
1448!> with ' denoting the auxiliary basis set and
1449!>
1450!> PBEx: PBE exchange functional
1451!> XWPBEX: PBE exchange hole for short-range potential (erfc(omega*r)/r)
1452!> XWPBEX0: PBE exchange hole for standard coulomb potential
1453!> PBE_HOLE_TC_LR: PBE exchange hole for longrange truncated coulomb potential
1454!>
1455!> Above explanation is correct for the deafult case. If a specific functional is requested
1456!> for the correction term (cfun), we get
1457!> Ex,hf = Ex,hf' + (cfun-cfun')
1458!> for all cases of operators.
1459!>
1460!> \param x_data ...
1461!> \param xc_section the original xc_section
1462!> \param admm_env the ADMM environment
1463!> \par History
1464!> 12.2009 created [Manuel Guidon]
1465!> 05.2021 simplify for case of no correction [JGH]
1466!> \author Manuel Guidon
1467! **************************************************************************************************
1468 SUBROUTINE create_admm_xc_section(x_data, xc_section, admm_env)
1469 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1470 TYPE(section_vals_type), POINTER :: xc_section
1471 TYPE(admm_type), POINTER :: admm_env
1472
1473 LOGICAL, PARAMETER :: debug_functional = .false.
1474#if defined (__LIBXC)
1475 REAL(kind=dp), PARAMETER :: x_factor_c = 0.930525736349100025_dp
1476#endif
1477
1478 CHARACTER(LEN=20) :: name_x_func
1479 INTEGER :: hfx_potential_type, ifun, iounit, nfun
1480 LOGICAL :: funct_found
1481 REAL(dp) :: cutoff_radius, hfx_fraction, omega, &
1482 scale_coulomb, scale_longrange, scale_x
1483 TYPE(cp_logger_type), POINTER :: logger
1484 TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_section
1485
1486 logger => cp_get_default_logger()
1487 NULLIFY (admm_env%xc_section_aux, admm_env%xc_section_primary)
1488
1489 !! ** Duplicate existing xc-section
1490 CALL section_vals_duplicate(xc_section, admm_env%xc_section_aux)
1491 CALL section_vals_duplicate(xc_section, admm_env%xc_section_primary)
1492 !** Now modify the auxiliary basis
1493 !** First remove all functionals
1494 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL")
1495
1496 !* Overwrite possible shortcut
1497 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1498 i_val=xc_funct_no_shortcut)
1499
1500 !** Get number of Functionals in the list
1501 ifun = 0
1502 nfun = 0
1503 DO
1504 ifun = ifun + 1
1505 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1506 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1507 nfun = nfun + 1
1508 END DO
1509
1510 ifun = 0
1511 DO ifun = 1, nfun
1512 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=1)
1513 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1514 CALL section_vals_remove_values(xc_fun)
1515 END DO
1516
1517 IF (ASSOCIATED(x_data)) THEN
1518 hfx_potential_type = x_data(1, 1)%potential_parameter%potential_type
1519 hfx_fraction = x_data(1, 1)%general_parameter%fraction
1520 ELSE
1521 cpwarn("ADMM requested without a DFT%XC%HF section. It will be ignored for the SCF.")
1522 admm_env%aux_exch_func = do_admm_aux_exch_func_none
1523 END IF
1524
1525 !in case of no admm exchange corr., no auxiliary exchange functional needed
1526 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
1527 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1528 i_val=xc_none)
1529 hfx_fraction = 0.0_dp
1530 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default) THEN
1531 ! default PBE Functional
1532 !! ** Add functionals evaluated with auxiliary basis
1533 SELECT CASE (hfx_potential_type)
1534 CASE (do_potential_coulomb)
1535 CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", &
1536 l_val=.true.)
1537 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1538 r_val=-hfx_fraction)
1539 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", &
1540 r_val=0.0_dp)
1541 CASE (do_potential_short)
1542 omega = x_data(1, 1)%potential_parameter%omega
1543 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1544 l_val=.true.)
1545 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1546 r_val=-hfx_fraction)
1547 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1548 r_val=0.0_dp)
1549 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1550 r_val=omega)
1551 CASE (do_potential_truncated)
1552 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1553 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1554 l_val=.true.)
1555 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1556 r_val=hfx_fraction)
1557 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1558 r_val=cutoff_radius)
1559 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1560 l_val=.true.)
1561 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1562 r_val=0.0_dp)
1563 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1564 r_val=-hfx_fraction)
1565 CASE (do_potential_long)
1566 omega = x_data(1, 1)%potential_parameter%omega
1567 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1568 l_val=.true.)
1569 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1570 r_val=hfx_fraction)
1571 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1572 r_val=-hfx_fraction)
1573 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1574 r_val=omega)
1575 CASE (do_potential_mix_cl)
1576 omega = x_data(1, 1)%potential_parameter%omega
1577 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1578 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1579 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1580 l_val=.true.)
1581 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1582 r_val=hfx_fraction*scale_longrange)
1583 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1584 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1585 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1586 r_val=omega)
1587 CASE (do_potential_mix_cl_trunc)
1588 omega = x_data(1, 1)%potential_parameter%omega
1589 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1590 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1591 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1592 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1593 l_val=.true.)
1594 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1595 r_val=hfx_fraction*(scale_longrange + scale_coulomb))
1596 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1597 r_val=cutoff_radius)
1598 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1599 l_val=.true.)
1600 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1601 r_val=hfx_fraction*scale_longrange)
1602 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1603 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1604 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1605 r_val=omega)
1606 CASE DEFAULT
1607 cpabort("Unknown potential operator!")
1608 END SELECT
1609
1610 !** Now modify the functionals for the primary basis
1611 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
1612 !* Overwrite possible shortcut
1613 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1614 i_val=xc_funct_no_shortcut)
1615
1616 SELECT CASE (hfx_potential_type)
1617 CASE (do_potential_coulomb)
1618 ifun = 0
1619 funct_found = .false.
1620 DO
1621 ifun = ifun + 1
1622 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1623 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1624 IF (xc_fun%section%name == "PBE") THEN
1625 funct_found = .true.
1626 END IF
1627 END DO
1628 IF (.NOT. funct_found) THEN
1629 CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", &
1630 l_val=.true.)
1631 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1632 r_val=hfx_fraction)
1633 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", &
1634 r_val=0.0_dp)
1635 ELSE
1636 CALL section_vals_val_get(xc_fun_section, "PBE%SCALE_X", &
1637 r_val=scale_x)
1638 scale_x = scale_x + hfx_fraction
1639 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1640 r_val=scale_x)
1641 END IF
1642 CASE (do_potential_short)
1643 omega = x_data(1, 1)%potential_parameter%omega
1644 ifun = 0
1645 funct_found = .false.
1646 DO
1647 ifun = ifun + 1
1648 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1649 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1650 IF (xc_fun%section%name == "XWPBE") THEN
1651 funct_found = .true.
1652 END IF
1653 END DO
1654 IF (.NOT. funct_found) THEN
1655 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1656 l_val=.true.)
1657 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1658 r_val=hfx_fraction)
1659 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1660 r_val=0.0_dp)
1661 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1662 r_val=omega)
1663 ELSE
1664 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1665 r_val=scale_x)
1666 scale_x = scale_x + hfx_fraction
1667 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1668 r_val=scale_x)
1669 END IF
1670 CASE (do_potential_long)
1671 omega = x_data(1, 1)%potential_parameter%omega
1672 ifun = 0
1673 funct_found = .false.
1674 DO
1675 ifun = ifun + 1
1676 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1677 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1678 IF (xc_fun%section%name == "XWPBE") THEN
1679 funct_found = .true.
1680 END IF
1681 END DO
1682 IF (.NOT. funct_found) THEN
1683 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1684 l_val=.true.)
1685 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1686 r_val=-hfx_fraction)
1687 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1688 r_val=hfx_fraction)
1689 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1690 r_val=omega)
1691 ELSE
1692 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1693 r_val=scale_x)
1694 scale_x = scale_x - hfx_fraction
1695 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1696 r_val=scale_x)
1697 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1698 r_val=scale_x)
1699 scale_x = scale_x + hfx_fraction
1700 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1701 r_val=scale_x)
1702
1703 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1704 r_val=omega)
1705 END IF
1706 CASE (do_potential_truncated)
1707 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1708 ifun = 0
1709 funct_found = .false.
1710 DO
1711 ifun = ifun + 1
1712 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1713 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1714 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
1715 funct_found = .true.
1716 END IF
1717 END DO
1718 IF (.NOT. funct_found) THEN
1719 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1720 l_val=.true.)
1721 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1722 r_val=-hfx_fraction)
1723 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1724 r_val=cutoff_radius)
1725 ELSE
1726 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1727 r_val=scale_x)
1728 scale_x = scale_x - hfx_fraction
1729 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1730 r_val=scale_x)
1731 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1732 r_val=cutoff_radius)
1733 END IF
1734 ifun = 0
1735 funct_found = .false.
1736 DO
1737 ifun = ifun + 1
1738 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1739 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1740 IF (xc_fun%section%name == "XWPBE") THEN
1741 funct_found = .true.
1742 END IF
1743 END DO
1744 IF (.NOT. funct_found) THEN
1745 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1746 l_val=.true.)
1747 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1748 r_val=hfx_fraction)
1749 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1750 r_val=0.0_dp)
1751
1752 ELSE
1753 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1754 r_val=scale_x)
1755 scale_x = scale_x + hfx_fraction
1756 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1757 r_val=scale_x)
1758 END IF
1759 CASE (do_potential_mix_cl_trunc)
1760 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1761 omega = x_data(1, 1)%potential_parameter%omega
1762 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1763 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1764 ifun = 0
1765 funct_found = .false.
1766 DO
1767 ifun = ifun + 1
1768 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1769 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1770 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
1771 funct_found = .true.
1772 END IF
1773 END DO
1774 IF (.NOT. funct_found) THEN
1775 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1776 l_val=.true.)
1777 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1778 r_val=-hfx_fraction*(scale_coulomb + scale_longrange))
1779 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1780 r_val=cutoff_radius)
1781
1782 ELSE
1783 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1784 r_val=scale_x)
1785 scale_x = scale_x - hfx_fraction*(scale_coulomb + scale_longrange)
1786 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1787 r_val=scale_x)
1788 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1789 r_val=cutoff_radius)
1790 END IF
1791 ifun = 0
1792 funct_found = .false.
1793 DO
1794 ifun = ifun + 1
1795 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1796 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1797 IF (xc_fun%section%name == "XWPBE") THEN
1798 funct_found = .true.
1799 END IF
1800 END DO
1801 IF (.NOT. funct_found) THEN
1802 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1803 l_val=.true.)
1804 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1805 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
1806 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1807 r_val=-hfx_fraction*scale_longrange)
1808 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1809 r_val=omega)
1810
1811 ELSE
1812 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1813 r_val=scale_x)
1814 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
1815 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1816 r_val=scale_x)
1817 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1818 r_val=scale_x)
1819 scale_x = scale_x - hfx_fraction*scale_longrange
1820 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1821 r_val=scale_x)
1822
1823 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1824 r_val=omega)
1825 END IF
1826 CASE (do_potential_mix_cl)
1827 omega = x_data(1, 1)%potential_parameter%omega
1828 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1829 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1830 ifun = 0
1831 funct_found = .false.
1832 DO
1833 ifun = ifun + 1
1834 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1835 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1836 IF (xc_fun%section%name == "XWPBE") THEN
1837 funct_found = .true.
1838 END IF
1839 END DO
1840 IF (.NOT. funct_found) THEN
1841 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1842 l_val=.true.)
1843 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1844 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
1845 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1846 r_val=-hfx_fraction*scale_longrange)
1847 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1848 r_val=omega)
1849
1850 ELSE
1851 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1852 r_val=scale_x)
1853 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
1854 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1855 r_val=scale_x)
1856
1857 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1858 r_val=scale_x)
1859 scale_x = scale_x - hfx_fraction*scale_longrange
1860 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1861 r_val=scale_x)
1862
1863 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1864 r_val=omega)
1865 END IF
1866 END SELECT
1867 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default_libxc) THEN
1868 ! default PBE Functional
1869 !! ** Add functionals evaluated with auxiliary basis
1870#if defined (__LIBXC)
1871 SELECT CASE (hfx_potential_type)
1872 CASE (do_potential_coulomb)
1873 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1874 l_val=.true.)
1875 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1876 r_val=-hfx_fraction)
1877 CASE (do_potential_short)
1878 omega = x_data(1, 1)%potential_parameter%omega
1879 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1880 l_val=.true.)
1881 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1882 r_val=-hfx_fraction)
1883 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1884 r_val=omega)
1885 CASE (do_potential_truncated)
1886 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1887 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1888 l_val=.true.)
1889 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1890 r_val=hfx_fraction)
1891 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1892 r_val=cutoff_radius)
1893 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1894 l_val=.true.)
1895 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1896 r_val=-hfx_fraction)
1897 CASE (do_potential_long)
1898 omega = x_data(1, 1)%potential_parameter%omega
1899 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1900 l_val=.true.)
1901 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1902 r_val=hfx_fraction)
1903 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1904 r_val=omega)
1905 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1906 l_val=.true.)
1907 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1908 r_val=-hfx_fraction)
1909 CASE (do_potential_mix_cl)
1910 omega = x_data(1, 1)%potential_parameter%omega
1911 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1912 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1913 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1914 l_val=.true.)
1915 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1916 r_val=hfx_fraction*scale_longrange)
1917 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1918 r_val=omega)
1919 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1920 l_val=.true.)
1921 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1922 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1923 CASE (do_potential_mix_cl_trunc)
1924 omega = x_data(1, 1)%potential_parameter%omega
1925 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1926 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1927 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1928 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1929 l_val=.true.)
1930 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1931 r_val=hfx_fraction*(scale_longrange + scale_coulomb))
1932 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1933 r_val=cutoff_radius)
1934 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1935 l_val=.true.)
1936 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1937 r_val=hfx_fraction*scale_longrange)
1938 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1939 r_val=omega)
1940 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1941 l_val=.true.)
1942 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1943 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1944 CASE DEFAULT
1945 cpabort("Unknown potential operator!")
1946 END SELECT
1947
1948 !** Now modify the functionals for the primary basis
1949 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
1950 !* Overwrite possible shortcut
1951 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1952 i_val=xc_funct_no_shortcut)
1953
1954 SELECT CASE (hfx_potential_type)
1955 CASE (do_potential_coulomb)
1956 ifun = 0
1957 funct_found = .false.
1958 DO
1959 ifun = ifun + 1
1960 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1961 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1962 IF (xc_fun%section%name == "GGA_X_PBE") THEN
1963 funct_found = .true.
1964 END IF
1965 END DO
1966 IF (.NOT. funct_found) THEN
1967 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1968 l_val=.true.)
1969 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1970 r_val=hfx_fraction)
1971 ELSE
1972 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
1973 r_val=scale_x)
1974 scale_x = scale_x + hfx_fraction
1975 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1976 r_val=scale_x)
1977 END IF
1978 CASE (do_potential_short)
1979 omega = x_data(1, 1)%potential_parameter%omega
1980 ifun = 0
1981 funct_found = .false.
1982 DO
1983 ifun = ifun + 1
1984 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1985 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1986 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
1987 funct_found = .true.
1988 END IF
1989 END DO
1990 IF (.NOT. funct_found) THEN
1991 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1992 l_val=.true.)
1993 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1994 r_val=hfx_fraction)
1995 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1996 r_val=omega)
1997 ELSE
1998 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1999 r_val=scale_x)
2000 scale_x = scale_x + hfx_fraction
2001 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2002 r_val=scale_x)
2003 END IF
2004 CASE (do_potential_long)
2005 omega = x_data(1, 1)%potential_parameter%omega
2006 ifun = 0
2007 funct_found = .false.
2008 DO
2009 ifun = ifun + 1
2010 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2011 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2012 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2013 funct_found = .true.
2014 END IF
2015 END DO
2016 IF (.NOT. funct_found) THEN
2017 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2018 l_val=.true.)
2019 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2020 r_val=-hfx_fraction)
2021 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2022 r_val=omega)
2023 ELSE
2024 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2025 r_val=scale_x)
2026 scale_x = scale_x - hfx_fraction
2027 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2028 r_val=scale_x)
2029
2030 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2031 r_val=omega)
2032 END IF
2033 ifun = 0
2034 funct_found = .false.
2035 DO
2036 ifun = ifun + 1
2037 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2038 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2039 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2040 funct_found = .true.
2041 END IF
2042 END DO
2043 IF (.NOT. funct_found) THEN
2044 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2045 l_val=.true.)
2046 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2047 r_val=hfx_fraction)
2048 ELSE
2049 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2050 r_val=scale_x)
2051 scale_x = scale_x + hfx_fraction
2052 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2053 r_val=scale_x)
2054 END IF
2055 CASE (do_potential_truncated)
2056 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
2057 ifun = 0
2058 funct_found = .false.
2059 DO
2060 ifun = ifun + 1
2061 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2062 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2063 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
2064 funct_found = .true.
2065 END IF
2066 END DO
2067 IF (.NOT. funct_found) THEN
2068 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
2069 l_val=.true.)
2070 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2071 r_val=-hfx_fraction)
2072 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2073 r_val=cutoff_radius)
2074
2075 ELSE
2076 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2077 r_val=scale_x)
2078 scale_x = scale_x - hfx_fraction
2079 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2080 r_val=scale_x)
2081 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2082 r_val=cutoff_radius)
2083 END IF
2084 ifun = 0
2085 funct_found = .false.
2086 DO
2087 ifun = ifun + 1
2088 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2089 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2090 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2091 funct_found = .true.
2092 END IF
2093 END DO
2094 IF (.NOT. funct_found) THEN
2095 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2096 l_val=.true.)
2097 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2098 r_val=hfx_fraction)
2099
2100 ELSE
2101 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2102 r_val=scale_x)
2103 scale_x = scale_x + hfx_fraction
2104 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2105 r_val=scale_x)
2106 END IF
2107 CASE (do_potential_mix_cl_trunc)
2108 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
2109 omega = x_data(1, 1)%potential_parameter%omega
2110 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
2111 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
2112 ifun = 0
2113 funct_found = .false.
2114 DO
2115 ifun = ifun + 1
2116 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2117 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2118 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
2119 funct_found = .true.
2120 END IF
2121 END DO
2122 IF (.NOT. funct_found) THEN
2123 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
2124 l_val=.true.)
2125 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2126 r_val=-hfx_fraction*(scale_coulomb + scale_longrange))
2127 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2128 r_val=cutoff_radius)
2129
2130 ELSE
2131 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2132 r_val=scale_x)
2133 scale_x = scale_x - hfx_fraction*(scale_coulomb + scale_longrange)
2134 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2135 r_val=scale_x)
2136 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2137 r_val=cutoff_radius)
2138 END IF
2139 ifun = 0
2140 funct_found = .false.
2141 DO
2142 ifun = ifun + 1
2143 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2144 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2145 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2146 funct_found = .true.
2147 END IF
2148 END DO
2149 IF (.NOT. funct_found) THEN
2150 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2151 l_val=.true.)
2152 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2153 r_val=-hfx_fraction*scale_longrange)
2154 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2155 r_val=omega)
2156
2157 ELSE
2158 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2159 r_val=scale_x)
2160 scale_x = scale_x - hfx_fraction*scale_longrange
2161 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2162 r_val=scale_x)
2163
2164 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2165 r_val=omega)
2166 END IF
2167 ifun = 0
2168 funct_found = .false.
2169 DO
2170 ifun = ifun + 1
2171 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2172 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2173 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2174 funct_found = .true.
2175 END IF
2176 END DO
2177 IF (.NOT. funct_found) THEN
2178 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2179 l_val=.true.)
2180 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2181 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
2182 ELSE
2183 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2184 r_val=scale_x)
2185 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
2186 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2187 r_val=scale_x)
2188 END IF
2189 CASE (do_potential_mix_cl)
2190 omega = x_data(1, 1)%potential_parameter%omega
2191 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
2192 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
2193 ifun = 0
2194 funct_found = .false.
2195 DO
2196 ifun = ifun + 1
2197 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2198 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2199 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2200 funct_found = .true.
2201 END IF
2202 END DO
2203 IF (.NOT. funct_found) THEN
2204 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2205 l_val=.true.)
2206 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2207 r_val=-hfx_fraction*scale_longrange)
2208 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2209 r_val=omega)
2210
2211 ELSE
2212 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2213 r_val=scale_x)
2214 scale_x = scale_x - hfx_fraction*scale_longrange
2215 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2216 r_val=scale_x)
2217
2218 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2219 r_val=omega)
2220 END IF
2221 ifun = 0
2222 funct_found = .false.
2223 DO
2224 ifun = ifun + 1
2225 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2226 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2227 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2228 funct_found = .true.
2229 END IF
2230 END DO
2231 IF (.NOT. funct_found) THEN
2232 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2233 l_val=.true.)
2234 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2235 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
2236 ELSE
2237 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2238 r_val=scale_x)
2239 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
2240 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2241 r_val=scale_x)
2242 END IF
2243 END SELECT
2244#else
2245 CALL cp_abort(__location__, "In order use a LibXC-based ADMM "// &
2246 "exchange correction functionals, you have to compile and link against LibXC!")
2247#endif
2248
2249 ! PBEX (always bare form), OPTX and Becke88 functional
2250 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex .OR. &
2251 admm_env%aux_exch_func == do_admm_aux_exch_func_opt .OR. &
2252 admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN
2253 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2254 name_x_func = 'PBE'
2255 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2256 name_x_func = 'OPTX'
2257 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN
2258 name_x_func = 'BECKE88'
2259 END IF
2260 !primary basis
2261 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2262 l_val=.true.)
2263 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2264 r_val=-hfx_fraction)
2265
2266 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2267 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_C", r_val=0.0_dp)
2268 END IF
2269
2270 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2271 IF (admm_env%aux_exch_func_param) THEN
2272 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A1", &
2273 r_val=admm_env%aux_x_param(1))
2274 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A2", &
2275 r_val=admm_env%aux_x_param(2))
2276 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%GAMMA", &
2277 r_val=admm_env%aux_x_param(3))
2278 END IF
2279 END IF
2280
2281 !** Now modify the functionals for the primary basis
2282 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2283 !* Overwrite possible L")
2284 !* Overwrite possible shortcut
2285 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
2286 i_val=xc_funct_no_shortcut)
2287
2288 ifun = 0
2289 funct_found = .false.
2290 DO
2291 ifun = ifun + 1
2292 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2293 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2294 IF (xc_fun%section%name == trim(name_x_func)) THEN
2295 funct_found = .true.
2296 END IF
2297 END DO
2298 IF (.NOT. funct_found) THEN
2299 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2300 l_val=.true.)
2301 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2302 r_val=hfx_fraction)
2303 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2304 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_C", &
2305 r_val=0.0_dp)
2306 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2307 IF (admm_env%aux_exch_func_param) THEN
2308 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A1", &
2309 r_val=admm_env%aux_x_param(1))
2310 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A2", &
2311 r_val=admm_env%aux_x_param(2))
2312 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%GAMMA", &
2313 r_val=admm_env%aux_x_param(3))
2314 END IF
2315 END IF
2316
2317 ELSE
2318 CALL section_vals_val_get(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2319 r_val=scale_x)
2320 scale_x = scale_x + hfx_fraction
2321 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2322 r_val=scale_x)
2323 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2324 cpassert(.NOT. admm_env%aux_exch_func_param)
2325 END IF
2326 END IF
2327
2328 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc .OR. &
2329 admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc .OR. &
2330 admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc .OR. &
2331 admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
2332#if defined(__LIBXC)
2333 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc) THEN
2334 name_x_func = 'GGA_X_PBE'
2335 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2336 name_x_func = 'GGA_X_OPTX'
2337 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
2338 name_x_func = 'GGA_X_B88'
2339 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc) THEN
2340 name_x_func = 'LDA_X'
2341 END IF
2342 !primary basis
2343 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2344 l_val=.true.)
2345 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE", &
2346 r_val=-hfx_fraction)
2347
2348 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2349 IF (admm_env%aux_exch_func_param) THEN
2350 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_A", &
2351 r_val=admm_env%aux_x_param(1))
2352 ! LibXC rescales the second parameter of the OPTX functional (see documentation there)
2353 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_B", &
2354 r_val=admm_env%aux_x_param(2)/x_factor_c)
2355 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_GAMMA", &
2356 r_val=admm_env%aux_x_param(3))
2357 END IF
2358 END IF
2359
2360 !** Now modify the functionals for the primary basis
2361 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2362 !* Overwrite possible L")
2363 !* Overwrite possible shortcut
2364 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
2365 i_val=xc_funct_no_shortcut)
2366
2367 ifun = 0
2368 funct_found = .false.
2369 DO
2370 ifun = ifun + 1
2371 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2372 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2373 IF (xc_fun%section%name == trim(name_x_func)) THEN
2374 funct_found = .true.
2375 END IF
2376 END DO
2377 IF (.NOT. funct_found) THEN
2378 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2379 l_val=.true.)
2380 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE", &
2381 r_val=hfx_fraction)
2382 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2383 IF (admm_env%aux_exch_func_param) THEN
2384 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_A", &
2385 r_val=admm_env%aux_x_param(1))
2386 ! LibXC rescales the second parameter of the OPTX functional (see documentation there)
2387 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_B", &
2388 r_val=admm_env%aux_x_param(2)/x_factor_c)
2389 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_GAMMA", &
2390 r_val=admm_env%aux_x_param(3))
2391 END IF
2392 END IF
2393
2394 ELSE
2395 CALL section_vals_val_get(xc_fun_section, trim(name_x_func)//"%SCALE", &
2396 r_val=scale_x)
2397 scale_x = scale_x + hfx_fraction
2398 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE", &
2399 r_val=scale_x)
2400 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2401 cpassert(.NOT. admm_env%aux_exch_func_param)
2402 END IF
2403 END IF
2404#else
2405 CALL cp_abort(__location__, "In order use a LibXC-based ADMM "// &
2406 "exchange correction functionals, you have to compile and link against LibXC!")
2407#endif
2408
2409 ELSE
2410 cpabort("Unknown exchange correction functional!")
2411 END IF
2412
2413 IF (debug_functional) THEN
2414 iounit = cp_logger_get_default_io_unit(logger)
2415 IF (iounit > 0) THEN
2416 WRITE (iounit, "(A)") " ADMM Primary Basis Set Functional"
2417 END IF
2418 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2419 ifun = 0
2420 funct_found = .false.
2421 DO
2422 ifun = ifun + 1
2423 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2424 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2425
2426 scale_x = -1000.0_dp
2427 IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
2428 CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
2429 END IF
2430 IF (xc_fun%section%name == "XWPBE") THEN
2431 CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
2432 IF (iounit > 0) THEN
2433 WRITE (iounit, "(T5,A,T25,2F10.3)") trim(xc_fun%section%name), scale_x, hfx_fraction
2434 END IF
2435 ELSE
2436 IF (iounit > 0) THEN
2437 WRITE (iounit, "(T5,A,T25,F10.3)") trim(xc_fun%section%name), scale_x
2438 END IF
2439 END IF
2440 END DO
2441
2442 IF (iounit > 0) THEN
2443 WRITE (iounit, "(A)") " Auxiliary Basis Set Functional"
2444 END IF
2445 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL")
2446 ifun = 0
2447 funct_found = .false.
2448 DO
2449 ifun = ifun + 1
2450 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2451 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2452 scale_x = -1000.0_dp
2453 IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
2454 CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
2455 END IF
2456 IF (xc_fun%section%name == "XWPBE") THEN
2457 CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
2458 IF (iounit > 0) THEN
2459 WRITE (iounit, "(T5,A,T25,2F10.3)") trim(xc_fun%section%name), scale_x, hfx_fraction
2460 END IF
2461 ELSE
2462 IF (iounit > 0) THEN
2463 WRITE (iounit, "(T5,A,T25,F10.3)") trim(xc_fun%section%name), scale_x
2464 END IF
2465 END IF
2466 END DO
2467 END IF
2468
2469 END SUBROUTINE create_admm_xc_section
2470
2471! **************************************************************************************************
2472!> \brief Add the hfx contributions to the Hamiltonian
2473!>
2474!> \param matrix_ks Kohn-Sham matrix (updated on exit)
2475!> \param rho_ao electron density expressed in terms of atomic orbitals
2476!> \param qs_env Quickstep environment
2477!> \param update_energy whether to update energy (default: yes)
2478!> \param recalc_integrals whether to recalculate integrals (default: value of HF%TREAT_LSD_IN_CORE)
2479!> \param external_hfx_sections ...
2480!> \param external_x_data ...
2481!> \param external_para_env ...
2482!> \note
2483!> Simplified version of subroutine hfx_ks_matrix()
2484! **************************************************************************************************
2485 SUBROUTINE tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, &
2486 external_hfx_sections, external_x_data, external_para_env)
2487 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
2488 TARGET :: matrix_ks, rho_ao
2489 TYPE(qs_environment_type), POINTER :: qs_env
2490 LOGICAL, INTENT(IN), OPTIONAL :: update_energy, recalc_integrals
2491 TYPE(section_vals_type), OPTIONAL, POINTER :: external_hfx_sections
2492 TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, TARGET :: external_x_data
2493 TYPE(mp_para_env_type), OPTIONAL, POINTER :: external_para_env
2494
2495 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddft_hfx_matrix'
2496
2497 INTEGER :: handle, irep, ispin, mspin, n_rep_hf, &
2498 nspins
2499 LOGICAL :: distribute_fock_matrix, &
2500 hfx_treat_lsd_in_core, &
2501 my_update_energy, s_mstruct_changed
2502 REAL(kind=dp) :: eh1, ehfx
2503 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp
2504 TYPE(dft_control_type), POINTER :: dft_control
2505 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
2506 TYPE(mp_para_env_type), POINTER :: para_env
2507 TYPE(qs_energy_type), POINTER :: energy
2508 TYPE(section_vals_type), POINTER :: hfx_sections, input
2509
2510 CALL timeset(routinen, handle)
2511
2512 NULLIFY (dft_control, hfx_sections, input, para_env, matrix_ks_kp, rho_ao_kp)
2513
2514 CALL get_qs_env(qs_env=qs_env, &
2515 dft_control=dft_control, &
2516 energy=energy, &
2517 input=input, &
2518 para_env=para_env, &
2519 s_mstruct_changed=s_mstruct_changed, &
2520 x_data=x_data)
2521
2522 ! This should probably be the HF section from the TDDFPT XC section!
2523 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
2524
2525 IF (PRESENT(external_hfx_sections)) hfx_sections => external_hfx_sections
2526 IF (PRESENT(external_x_data)) x_data => external_x_data
2527 IF (PRESENT(external_para_env)) para_env => external_para_env
2528
2529 my_update_energy = .true.
2530 IF (PRESENT(update_energy)) my_update_energy = update_energy
2531
2532 IF (PRESENT(recalc_integrals)) s_mstruct_changed = recalc_integrals
2533
2534 cpassert(dft_control%nimages == 1)
2535 nspins = dft_control%nspins
2536
2537 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
2538 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
2539 i_rep_section=1)
2540
2541 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
2542 distribute_fock_matrix = .true.
2543
2544 mspin = 1
2545 IF (hfx_treat_lsd_in_core) mspin = nspins
2546
2547 matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
2548 rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
2549
2550 DO irep = 1, n_rep_hf
2551 ! the real hfx calulation
2552 ehfx = 0.0_dp
2553
2554 IF (x_data(irep, 1)%do_hfx_ri) THEN
2555 CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, &
2556 rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
2557 nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
2558
2559 ELSE
2560 DO ispin = 1, mspin
2561 CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
2562 s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
2563 ehfx = ehfx + eh1
2564 END DO
2565 END IF
2566 END DO
2567 IF (my_update_energy) energy%ex = ehfx
2568
2569 CALL timestop(handle)
2570 END SUBROUTINE tddft_hfx_matrix
2571
2572END MODULE hfx_admm_utils
Types and set/get functions for auxiliary density matrix methods.
subroutine, public admm_dm_create(admm_dm, admm_control, nspins, natoms)
Create a new admm_dm type.
Contains ADMM methods which require molecular orbitals.
subroutine, public scale_dm(qs_env, rho_ao_orb, scale_back)
Scale density matrix by gsi(ispin), is needed for force scaling in ADMMP.
subroutine, public kpoint_calc_admm_matrices(qs_env, calculate_forces)
Fill the ADMM overlp and basis change matrices in the KP env based on the real-space array.
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:593
subroutine, public set_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Set routine for the ADMM env.
Definition admm_types.F:677
subroutine, public admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit, blacs_env_ext)
creates ADMM environment, initializes the basic types
Definition admm_types.F:220
Define the atomic kind types and their sub types.
subroutine, public add_basis_set_to_container(container, basis_set, basis_set_type)
...
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
subroutine, public copy_gto_basis_set(basis_set_in, basis_set_out)
...
Handles all functions related to the CELL.
Definition cell_types.F:15
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition cell_types.F:257
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_m_by_n_from_row_template(matrix, template, n, sym)
Utility function to create dbcsr matrix, m x n matrix (n arbitrary) with the same processor grid and ...
pool for for elements that are retained and released
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition of the atomic potential types.
Utilities for hfx and admm methods.
subroutine, public hfx_admm_init(qs_env, calculate_forces, ext_xc_section)
...
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data, external_para_env)
Add the hfx contributions to the Hamiltonian.
subroutine, public aux_admm_init(qs_env, mos, admm_env, admm_control, basis_type)
Minimal setup routine for admm_env No forces No k-points No DFT correction terms.
subroutine, public hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, just_energy, v_rspace_new, v_tau_rspace, ext_xc_section)
Add the hfx contributions to the Hamiltonian.
subroutine, public create_admm_xc_section(x_data, xc_section, admm_env)
This routine modifies the xc section depending on the potential type used for the HF exchange and the...
subroutine, public hfx_ks_matrix_kp(qs_env, matrix_ks, energy, calculate_forces)
Add the HFX K-point contribution to the real-space Hamiltonians.
Routines to calculate derivatives with respect to basis function origin.
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data, nspins)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Routines to calculate HFX energy and potential.
subroutine, public integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, geometry_did_change, irep, distribute_fock_matrix, ispin, nspins)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
Test routines for HFX caclulations using PW.
subroutine, public pw_hfx(qs_env, ehfx, hfx_section, poisson_env, auxbas_pw_pool, irep)
computes the Hartree-Fock energy brute force in a pw basis
RI-methods for HFX and K-points. \auhtor Augustin Bussy (01.2023)
Definition hfx_ri_kp.F:13
subroutine, public hfx_ri_update_forces_kp(qs_env, ri_data, nspins, hf_fraction, rho_ao, use_virial)
Update the K-points RI-HFX forces.
Definition hfx_ri_kp.F:861
subroutine, public hfx_ri_update_ks_kp(qs_env, ri_data, ks_matrix, ehfx, rho_ao, geometry_did_change, nspins, hf_fraction)
Update the KS matrices for each real-space image.
Definition hfx_ri_kp.F:455
RI-methods for HFX.
Definition hfx_ri.F:12
subroutine, public hfx_ri_update_ks(qs_env, ri_data, ks_matrix, ehfx, mos, rho_ao, geometry_did_change, nspins, hf_fraction)
...
Definition hfx_ri.F:1041
subroutine, public hfx_ri_update_forces(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, mos, use_virial, resp_only, rescale_factor)
the general routine that calls the relevant force code
Definition hfx_ri.F:3044
Types and set/get functions for HFX.
Definition hfx_types.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_admm_aux_exch_func_opt_libxc
integer, parameter, public do_admm_purify_none
integer, parameter, public xc_funct_no_shortcut
integer, parameter, public do_admm_aux_exch_func_sx_libxc
integer, parameter, public do_admm_aux_exch_func_bee
integer, parameter, public do_potential_mix_cl
integer, parameter, public do_admm_basis_projection
integer, parameter, public do_admm_aux_exch_func_default_libxc
integer, parameter, public do_admm_aux_exch_func_opt
integer, parameter, public do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_bee_libxc
integer, parameter, public do_admm_aux_exch_func_pbex_libxc
integer, parameter, public do_admm_aux_exch_func_default
integer, parameter, public do_potential_truncated
integer, parameter, public do_admm_charge_constrained_projection
integer, parameter, public do_potential_id
integer, parameter, public do_potential_coulomb
integer, parameter, public do_potential_short
integer, parameter, public do_potential_mix_cl_trunc
integer, parameter, public xc_none
integer, parameter, public do_potential_long
integer, parameter, public do_admm_aux_exch_func_pbex
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
type(section_vals_type) function, pointer, public section_vals_get_subs_vals2(section_vals, i_section, i_rep_section)
returns the values of the n-th non default subsection (null if no such section exists (not so many no...
subroutine, public section_vals_remove_values(section_vals)
removes the values of a repetition of the section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_duplicate(section_vals_in, section_vals_out, i_rep_start, i_rep_end)
creates a deep copy from section_vals_in to section_vals_out
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Routines needed for kpoint calculation.
subroutine, public kpoint_initialize_mos(kpoint, mos, added_mos, for_aux_fit)
Initialize a set of MOs and density matrix for each kpoint (kpoint group)
Datatype to translate between k-points (2d) and gamma-point (1d) code.
subroutine, public kpoint_transitional_release(this)
Release the matrix set, using the right pointer.
subroutine, public set_2d_pointer(this, ptr_2d)
Assigns a 2D pointer.
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
real(kind=dp), parameter, public cutoff_screen_factor
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
Definition mathlib.F:1693
Interface to the message passing library MPI.
Define the data structure for the molecule information.
Define the data structure for the particle information.
subroutine, public get_paw_proj_set(paw_proj_set, csprj, chprj, first_prj, first_prjs, last_prj, local_oce_sphi_h, local_oce_sphi_s, maxl, ncgauprj, nsgauprj, nsatbas, nsotot, nprj, o2nindex, n2oindex, rcprj, rzetprj, zisomin, zetprj)
Get informations about a paw projectors set.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Calculate the interaction radii for the operator matrix calculation.
subroutine, public init_interaction_radii(qs_control, qs_kind_set)
Initialize all the atomic kind radii for a given threshold value.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
subroutine, public init_gapw_nlcc(qs_kind_set)
...
subroutine, public init_gapw_basis_set(qs_kind_set, qs_control, force_env_section, modify_qs_control)
...
subroutine, public local_rho_set_create(local_rho_set)
...
wrapper for the pools of matrixes
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Generate the atomic neighbor lists.
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep)
Build simple pair neighbor lists.
subroutine, public write_neighbor_lists(ab, particle_set, cell, para_env, neighbor_list_section, nl_type, middle_name, nlname)
Write a set of neighbor lists to the output unit.
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
Routines for the construction of the coefficients for the expansion of the atomic densities rho1_hard...
subroutine, public build_oce_matrices(intac, calculate_forces, nder, qs_kind_set, particle_set, sap_oce, eps_fit)
Set up the sparse matrix for the coefficients of one center expansions This routine uses the same log...
subroutine, public allocate_oce_set(oce_set, nkind)
Allocate and initialize the matrix set of oce coefficients.
subroutine, public create_oce_set(oce_set)
...
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:120
subroutine, public init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external)
rebuilds rho (if necessary allocating and initializing it)
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
subroutine, public qs_rho_create(rho)
Allocates a new instance of rho.
Types and set_get for real time propagation depending on runtype and diagonalization method different...
generate the tasks lists used by collocate and integrate routines
subroutine, public generate_qs_task_list(ks_env, task_list, basis_type, reorder_rs_grid_ranks, skip_load_balance_distributed, pw_env_external, sab_orb_external)
...
types for task lists
subroutine, public deallocate_task_list(task_list)
deallocates the components and the object itself
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type
subroutine, public rescale_xc_potential(qs_env, ks_matrix, rho, energy, v_rspace_new, v_tau_rspace, hf_energy, just_energy, calculate_forces, use_virial)
A subtype of the admm_env that contains the extra data needed for an ADMM GAPW calculation.
Definition admm_types.F:85
stores some data used in wavefunction fitting
Definition admm_types.F:120
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
structure to store local (to a processor) ordered lists of integers.
distributes pairs on a 2d grid of processors
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:510
Contains information about kpoints.
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.