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