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