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