(git:15c1bfc)
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, &
767 reorder_rs_grid_ranks=.false., soft_valid=.false., &
768 basis_type=aux_basis_type, &
769 skip_load_balance_distributed=skip_load_balance_distributed, &
770 sab_orb_external=admm_env%sab_aux_fit)
771
772 !The aux_fit densities
773 CALL qs_rho_rebuild(admm_env%rho_aux_fit, qs_env=qs_env, admm=.true.)
774 CALL qs_rho_rebuild(admm_env%rho_aux_fit_buffer, qs_env=qs_env, admm=.true.)
775
776 CALL timestop(handle)
777
778 END SUBROUTINE admm_update_s_mstruct
779
780! **************************************************************************************************
781!> \brief Update the admm_gapw_env internals to the current qs_env (i.e. atomic positions)
782!> \param qs_env ...
783! **************************************************************************************************
784 SUBROUTINE update_admm_gapw(qs_env)
785
786 TYPE(qs_environment_type), POINTER :: qs_env
787
788 CHARACTER(len=*), PARAMETER :: routinen = 'update_admm_gapw'
789
790 INTEGER :: handle, ikind, nkind
791 LOGICAL :: paw_atom
792 LOGICAL, ALLOCATABLE, DIMENSION(:) :: aux_present, oce_present
793 REAL(dp) :: subcells
794 REAL(dp), ALLOCATABLE, DIMENSION(:) :: aux_radius, oce_radius
795 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
796 TYPE(admm_gapw_r3d_rs_type), POINTER :: admm_gapw_env
797 TYPE(admm_type), POINTER :: admm_env
798 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
799 TYPE(cell_type), POINTER :: cell
800 TYPE(dft_control_type), POINTER :: dft_control
801 TYPE(distribution_1d_type), POINTER :: distribution_1d
802 TYPE(distribution_2d_type), POINTER :: distribution_2d
803 TYPE(gto_basis_set_type), POINTER :: aux_fit_basis
804 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
805 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
806 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
807 POINTER :: sap_oce
808 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
809 TYPE(paw_proj_set_type), POINTER :: paw_proj
810 TYPE(qs_kind_type), DIMENSION(:), POINTER :: admm_kind_set, qs_kind_set
811 TYPE(qs_ks_env_type), POINTER :: ks_env
812
813 NULLIFY (ks_env, qs_kind_set, admm_kind_set, aux_fit_basis, cell, distribution_1d)
814 NULLIFY (distribution_2d, paw_proj, particle_set, molecule_set, admm_env, admm_gapw_env)
815 NULLIFY (dft_control, atomic_kind_set, sap_oce)
816
817 CALL timeset(routinen, handle)
818
819 CALL get_qs_env(qs_env, ks_env=ks_env, qs_kind_set=qs_kind_set, admm_env=admm_env, &
820 dft_control=dft_control)
821 admm_gapw_env => admm_env%admm_gapw_env
822 admm_kind_set => admm_gapw_env%admm_kind_set
823 nkind = SIZE(qs_kind_set)
824
825 !Update the task lisft for the AUX_FIT_SOFT basis
826 IF (ASSOCIATED(admm_gapw_env%task_list)) CALL deallocate_task_list(admm_gapw_env%task_list)
827 CALL allocate_task_list(admm_gapw_env%task_list)
828
829 !note: we set soft_valid to .FALSE. want to use AUX_FIT_SOFT and not the normal ORB SOFT basis
830 CALL generate_qs_task_list(ks_env, admm_gapw_env%task_list, reorder_rs_grid_ranks=.false., &
831 soft_valid=.false., basis_type="AUX_FIT_SOFT", &
832 skip_load_balance_distributed=dft_control%qs_control%skip_load_balance_distributed, &
833 sab_orb_external=admm_env%sab_aux_fit)
834
835 !Update the precomputed oce integrals
836 !a sap_oce neighbor list is required => build it here
837 ALLOCATE (aux_present(nkind), oce_present(nkind))
838 aux_present = .false.; oce_present = .false.
839 ALLOCATE (aux_radius(nkind), oce_radius(nkind))
840 aux_radius = 0.0_dp; oce_radius = 0.0_dp
841
842 DO ikind = 1, nkind
843 CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis, basis_type="AUX_FIT")
844 IF (ASSOCIATED(aux_fit_basis)) THEN
845 aux_present(ikind) = .true.
846 CALL get_gto_basis_set(aux_fit_basis, kind_radius=aux_radius(ikind))
847 END IF
848
849 !note: get oce info from admm_kind_set
850 CALL get_qs_kind(admm_kind_set(ikind), paw_atom=paw_atom, paw_proj_set=paw_proj)
851 IF (paw_atom) THEN
852 oce_present(ikind) = .true.
853 CALL get_paw_proj_set(paw_proj, rcprj=oce_radius(ikind))
854 END IF
855 END DO
856
857 ALLOCATE (pair_radius(nkind, nkind))
858 pair_radius = 0.0_dp
859 CALL pair_radius_setup(aux_present, oce_present, aux_radius, oce_radius, pair_radius)
860
861 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
862 distribution_2d=distribution_2d, local_particles=distribution_1d, &
863 particle_set=particle_set, molecule_set=molecule_set)
864 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
865
866 ALLOCATE (atom2d(nkind))
867 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
868 molecule_set, .false., particle_set)
869 CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, &
870 subcells=subcells, operator_type="ABBA", nlname="AUX_PAW-PRJ")
871 CALL atom2d_cleanup(atom2d)
872
873 !actually compute the oce matrices
874 CALL create_oce_set(admm_gapw_env%oce)
875 CALL allocate_oce_set(admm_gapw_env%oce, nkind)
876
877 !always compute the derivative, cheap anyways
878 CALL build_oce_matrices(admm_gapw_env%oce%intac, calculate_forces=.true., nder=1, &
879 qs_kind_set=admm_kind_set, particle_set=particle_set, &
880 sap_oce=sap_oce, eps_fit=dft_control%qs_control%gapw_control%eps_fit)
881
882 CALL release_neighbor_list_sets(sap_oce)
883
884 CALL timestop(handle)
885
886 END SUBROUTINE update_admm_gapw
887
888! **************************************************************************************************
889!> \brief Allocates the various ADMM KS matrices
890!> \param admm_env ...
891!> \param qs_env ...
892! **************************************************************************************************
893 SUBROUTINE admm_alloc_ks_matrices(admm_env, qs_env)
894
895 TYPE(admm_type), POINTER :: admm_env
896 TYPE(qs_environment_type), POINTER :: qs_env
897
898 CHARACTER(len=*), PARAMETER :: routinen = 'admm_alloc_ks_matrices'
899
900 INTEGER :: handle, ic, ispin
901 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit_dft_kp, &
902 matrix_ks_aux_fit_hfx_kp, &
903 matrix_ks_aux_fit_kp, &
904 matrix_s_aux_fit_kp
905 TYPE(dft_control_type), POINTER :: dft_control
906
907 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)
908
909 CALL timeset(routinen, handle)
910
911 CALL get_qs_env(qs_env, dft_control=dft_control)
912 CALL get_admm_env(admm_env, matrix_s_aux_fit_kp=matrix_s_aux_fit_kp)
913
914 CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit)
915 CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit_dft)
916 CALL kpoint_transitional_release(admm_env%matrix_ks_aux_fit_hfx)
917
918 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_kp, dft_control%nspins, dft_control%nimages)
919 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_dft_kp, dft_control%nspins, dft_control%nimages)
920 CALL dbcsr_allocate_matrix_set(matrix_ks_aux_fit_hfx_kp, dft_control%nspins, dft_control%nimages)
921
922 DO ispin = 1, dft_control%nspins
923 DO ic = 1, dft_control%nimages
924 ALLOCATE (matrix_ks_aux_fit_kp(ispin, ic)%matrix)
925 CALL dbcsr_create(matrix_ks_aux_fit_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, ic)%matrix, &
926 name="KOHN-SHAM_MATRIX for ADMM")
927 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
928 CALL dbcsr_set(matrix_ks_aux_fit_kp(ispin, ic)%matrix, 0.0_dp)
929
930 ALLOCATE (matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix)
931 CALL dbcsr_create(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, 1)%matrix, &
932 name="KOHN-SHAM_MATRIX for ADMM")
933 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
934 CALL dbcsr_set(matrix_ks_aux_fit_dft_kp(ispin, ic)%matrix, 0.0_dp)
935
936 ALLOCATE (matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix)
937 CALL dbcsr_create(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, template=matrix_s_aux_fit_kp(1, 1)%matrix, &
938 name="KOHN-SHAM_MATRIX for ADMM")
939 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, admm_env%sab_aux_fit)
940 CALL dbcsr_set(matrix_ks_aux_fit_hfx_kp(ispin, ic)%matrix, 0.0_dp)
941 END DO
942 END DO
943
944 CALL set_admm_env(admm_env, &
945 matrix_ks_aux_fit_kp=matrix_ks_aux_fit_kp, &
946 matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft_kp, &
947 matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx_kp)
948
949 CALL timestop(handle)
950
951 END SUBROUTINE admm_alloc_ks_matrices
952
953! **************************************************************************************************
954!> \brief Add the HFX K-point contribution to the real-space Hamiltonians
955!> \param qs_env ...
956!> \param matrix_ks ...
957!> \param energy ...
958!> \param calculate_forces ...
959! **************************************************************************************************
960 SUBROUTINE hfx_ks_matrix_kp(qs_env, matrix_ks, energy, calculate_forces)
961 TYPE(qs_environment_type), POINTER :: qs_env
962 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
963 TYPE(qs_energy_type), POINTER :: energy
964 LOGICAL, INTENT(in) :: calculate_forces
965
966 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ks_matrix_kp'
967
968 INTEGER :: handle, img, irep, ispin, n_rep_hf, &
969 nimages, nspins
970 LOGICAL :: do_adiabatic_rescaling, &
971 s_mstruct_changed, use_virial
972 REAL(dp) :: eh1, ehfx, eold
973 REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy
974 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit_im, matrix_ks_im
975 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks_aux_fit_hfx_kp, &
976 matrix_ks_aux_fit_kp, matrix_ks_orb, &
977 rho_ao_orb
978 TYPE(dft_control_type), POINTER :: dft_control
979 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
980 TYPE(mp_para_env_type), POINTER :: para_env
981 TYPE(pw_env_type), POINTER :: pw_env
982 TYPE(pw_poisson_type), POINTER :: poisson_env
983 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
984 TYPE(qs_rho_type), POINTER :: rho_orb
985 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
986 hfx_sections, input
987 TYPE(virial_type), POINTER :: virial
988
989 CALL timeset(routinen, handle)
990
991 NULLIFY (auxbas_pw_pool, dft_control, hfx_sections, input, &
992 para_env, poisson_env, pw_env, virial, matrix_ks_im, &
993 matrix_ks_orb, rho_ao_orb, matrix_h, matrix_ks_aux_fit_kp, &
994 matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx_kp)
995
996 CALL get_qs_env(qs_env=qs_env, &
997 dft_control=dft_control, &
998 input=input, &
999 matrix_h_kp=matrix_h, &
1000 para_env=para_env, &
1001 pw_env=pw_env, &
1002 virial=virial, &
1003 matrix_ks_im=matrix_ks_im, &
1004 s_mstruct_changed=s_mstruct_changed, &
1005 x_data=x_data)
1006
1007 ! No RTP
1008 IF (qs_env%run_rtp) cpabort("No RTP implementation with K-points HFX")
1009
1010 ! No adiabatic rescaling
1011 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
1012 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1013 IF (do_adiabatic_rescaling) cpabort("No adiabatic rescaling implementation with K-points HFX")
1014
1015 IF (dft_control%do_admm) THEN
1016 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit_kp, &
1017 matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, &
1018 matrix_ks_aux_fit_hfx_kp=matrix_ks_aux_fit_hfx_kp)
1019 END IF
1020
1021 nspins = dft_control%nspins
1022 nimages = dft_control%nimages
1023
1024 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1025 IF (use_virial .AND. calculate_forces) virial%pv_fock_4c = 0.0_dp
1026
1027 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1028 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1029
1030 ! *** Initialize the auxiliary ks matrix to zero if required
1031 IF (dft_control%do_admm) THEN
1032 DO ispin = 1, nspins
1033 DO img = 1, nimages
1034 CALL dbcsr_set(matrix_ks_aux_fit_kp(ispin, img)%matrix, 0.0_dp)
1035 END DO
1036 END DO
1037 END IF
1038 DO ispin = 1, nspins
1039 DO img = 1, nimages
1040 CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
1041 END DO
1042 END DO
1043
1044 ALLOCATE (hf_energy(n_rep_hf))
1045
1046 eold = 0.0_dp
1047
1048 DO irep = 1, n_rep_hf
1049
1050 ! fetch the correct matrices for normal HFX or ADMM
1051 IF (dft_control%do_admm) THEN
1052 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks_orb, rho_aux_fit=rho_orb)
1053 ELSE
1054 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_orb, rho=rho_orb)
1055 END IF
1056 CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb)
1057
1058 ! Finally the real hfx calulation
1059 ehfx = 0.0_dp
1060
1061 IF (.NOT. x_data(irep, 1)%do_hfx_ri) THEN
1062 cpabort("Only RI-HFX is implemented for K-points")
1063 END IF
1064
1065 CALL hfx_ri_update_ks_kp(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1066 rho_ao_orb, s_mstruct_changed, nspins, &
1067 x_data(irep, 1)%general_parameter%fraction)
1068
1069 IF (calculate_forces) THEN
1070 !Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force
1071 IF (dft_control%do_admm) THEN
1072 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.false.)
1073 END IF
1074
1075 CALL hfx_ri_update_forces_kp(qs_env, x_data(irep, 1)%ri_data, nspins, &
1076 x_data(irep, 1)%general_parameter%fraction, &
1077 rho_ao_orb, use_virial=use_virial)
1078
1079 IF (dft_control%do_admm) THEN
1080 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.true.)
1081 END IF
1082 END IF
1083
1084 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
1085 eh1 = ehfx - eold
1086 CALL pw_hfx(qs_env, eh1, hfx_sections, poisson_env, auxbas_pw_pool, irep)
1087 eold = ehfx
1088
1089 END DO
1090
1091 ! *** Set the total HFX energy
1092 energy%ex = ehfx
1093
1094 ! *** Add Core-Hamiltonian-Matrix ***
1095 DO ispin = 1, nspins
1096 DO img = 1, nimages
1097 CALL dbcsr_add(matrix_ks(ispin, img)%matrix, matrix_h(1, img)%matrix, &
1098 1.0_dp, 1.0_dp)
1099 END DO
1100 END DO
1101 IF (use_virial .AND. calculate_forces) THEN
1102 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1103 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1104 virial%pv_calculate = .false.
1105 END IF
1106
1107 !update the hfx aux_fit matrix
1108 IF (dft_control%do_admm) THEN
1109 DO ispin = 1, nspins
1110 DO img = 1, nimages
1111 CALL dbcsr_add(matrix_ks_aux_fit_hfx_kp(ispin, img)%matrix, matrix_ks_aux_fit_kp(ispin, img)%matrix, &
1112 0.0_dp, 1.0_dp)
1113 END DO
1114 END DO
1115 END IF
1116
1117 CALL timestop(handle)
1118
1119 END SUBROUTINE hfx_ks_matrix_kp
1120
1121! **************************************************************************************************
1122!> \brief Add the hfx contributions to the Hamiltonian
1123!>
1124!> \param qs_env ...
1125!> \param matrix_ks ...
1126!> \param rho ...
1127!> \param energy ...
1128!> \param calculate_forces ...
1129!> \param just_energy ...
1130!> \param v_rspace_new ...
1131!> \param v_tau_rspace ...
1132!> \par History
1133!> refactoring 03-2011 [MI]
1134! **************************************************************************************************
1135
1136 SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
1137 just_energy, v_rspace_new, v_tau_rspace)
1138
1139 TYPE(qs_environment_type), POINTER :: qs_env
1140 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
1141 TYPE(qs_rho_type), POINTER :: rho
1142 TYPE(qs_energy_type), POINTER :: energy
1143 LOGICAL, INTENT(in) :: calculate_forces, just_energy
1144 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
1145
1146 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ks_matrix'
1147
1148 INTEGER :: handle, img, irep, ispin, mspin, &
1149 n_rep_hf, nimages, ns, nspins
1150 LOGICAL :: distribute_fock_matrix, &
1151 do_adiabatic_rescaling, &
1152 hfx_treat_lsd_in_core, &
1153 s_mstruct_changed, use_virial
1154 REAL(dp) :: eh1, ehfx, ehfxrt, eold
1155 REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy
1156 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_1d, matrix_ks_aux_fit, &
1157 matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_im, matrix_ks_im, rho_ao_1d, rho_ao_resp
1158 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_h_im, matrix_ks_orb, &
1159 rho_ao_orb
1160 TYPE(dft_control_type), POINTER :: dft_control
1161 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1162 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1163 TYPE(mp_para_env_type), POINTER :: para_env
1164 TYPE(pw_env_type), POINTER :: pw_env
1165 TYPE(pw_poisson_type), POINTER :: poisson_env
1166 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1167 TYPE(qs_rho_type), POINTER :: rho_orb
1168 TYPE(rt_prop_type), POINTER :: rtp
1169 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
1170 hfx_sections, input
1171 TYPE(virial_type), POINTER :: virial
1172
1173 CALL timeset(routinen, handle)
1174
1175 NULLIFY (auxbas_pw_pool, dft_control, hfx_sections, input, &
1176 para_env, poisson_env, pw_env, virial, matrix_ks_im, &
1177 matrix_ks_orb, rho_ao_orb, matrix_h, matrix_h_im, matrix_ks_aux_fit, &
1178 matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx)
1179
1180 CALL get_qs_env(qs_env=qs_env, &
1181 dft_control=dft_control, &
1182 input=input, &
1183 matrix_h_kp=matrix_h, &
1184 matrix_h_im_kp=matrix_h_im, &
1185 para_env=para_env, &
1186 pw_env=pw_env, &
1187 virial=virial, &
1188 matrix_ks_im=matrix_ks_im, &
1189 s_mstruct_changed=s_mstruct_changed, &
1190 x_data=x_data)
1191
1192 IF (dft_control%do_admm) THEN
1193 CALL get_admm_env(qs_env%admm_env, mos_aux_fit=mo_array, matrix_ks_aux_fit=matrix_ks_aux_fit, &
1194 matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
1195 ELSE
1196 CALL get_qs_env(qs_env=qs_env, mos=mo_array)
1197 END IF
1198
1199 nspins = dft_control%nspins
1200 nimages = dft_control%nimages
1201
1202 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1203
1204 IF (use_virial .AND. calculate_forces) virial%pv_fock_4c = 0.0_dp
1205
1206 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1207 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1208 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
1209 i_rep_section=1)
1210 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
1211 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
1212
1213 ! *** Initialize the auxiliary ks matrix to zero if required
1214 IF (dft_control%do_admm) THEN
1215 DO ispin = 1, nspins
1216 CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp)
1217 END DO
1218 END IF
1219 DO ispin = 1, nspins
1220 DO img = 1, nimages
1221 CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
1222 END DO
1223 END DO
1224
1225 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
1226
1227 ALLOCATE (hf_energy(n_rep_hf))
1228
1229 eold = 0.0_dp
1230
1231 DO irep = 1, n_rep_hf
1232 ! Remember: Vhfx is added, energy is calclulated from total Vhfx,
1233 ! so energy of last iteration is correct
1234
1235 IF (do_adiabatic_rescaling .AND. hfx_treat_lsd_in_core) &
1236 cpabort("HFX_TREAT_LSD_IN_CORE not implemented for adiabatically rescaled hybrids")
1237 ! everything is calculated with adiabatic rescaling but the potential is not added in a first step
1238 distribute_fock_matrix = .NOT. do_adiabatic_rescaling
1239
1240 mspin = 1
1241 IF (hfx_treat_lsd_in_core) mspin = nspins
1242
1243 ! fetch the correct matrices for normal HFX or ADMM
1244 IF (dft_control%do_admm) THEN
1245 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit=matrix_ks_1d, rho_aux_fit=rho_orb)
1246 ns = SIZE(matrix_ks_1d)
1247 matrix_ks_orb(1:ns, 1:1) => matrix_ks_1d(1:ns)
1248 ELSE
1249 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_orb, rho=rho_orb)
1250 END IF
1251 CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb)
1252 ! Finally the real hfx calulation
1253 ehfx = 0.0_dp
1254
1255 IF (x_data(irep, 1)%do_hfx_ri) THEN
1256
1257 CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1258 mo_array, rho_ao_orb, &
1259 s_mstruct_changed, nspins, &
1260 x_data(irep, 1)%general_parameter%fraction)
1261 IF (dft_control%do_admm) THEN
1262 !for ADMMS, we need the exchange matrix k(d) for both spins
1263 DO ispin = 1, nspins
1264 CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
1265 name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1266 END DO
1267 END IF
1268
1269 ELSE
1270
1271 DO ispin = 1, mspin
1272 CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
1273 para_env, s_mstruct_changed, irep, distribute_fock_matrix, &
1274 ispin=ispin)
1275 ehfx = ehfx + eh1
1276 END DO
1277 END IF
1278
1279 IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
1280 !Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force
1281 IF (dft_control%do_admm) THEN
1282 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.false.)
1283 END IF
1284 NULLIFY (rho_ao_resp)
1285
1286 IF (x_data(irep, 1)%do_hfx_ri) THEN
1287
1288 CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
1289 x_data(irep, 1)%general_parameter%fraction, &
1290 rho_ao=rho_ao_orb, mos=mo_array, &
1291 rho_ao_resp=rho_ao_resp, &
1292 use_virial=use_virial)
1293
1294 ELSE
1295
1296 CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
1297 para_env, irep, use_virial)
1298
1299 END IF
1300
1301 !Scale auxiliary density matrix for ADMMP back with 1/gsi(ispin)
1302 IF (dft_control%do_admm) THEN
1303 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.true.)
1304 END IF
1305 END IF
1306
1307 !! If required, the calculation of the forces will be done later with adiabatic rescaling
1308 IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx
1309
1310 ! special case RTP/EMD we have a full complex density and HFX has a contribution from the imaginary part
1311 ehfxrt = 0.0_dp
1312 IF (qs_env%run_rtp) THEN
1313
1314 CALL get_qs_env(qs_env=qs_env, rtp=rtp)
1315 DO ispin = 1, nspins
1316 CALL dbcsr_set(matrix_ks_im(ispin)%matrix, 0.0_dp)
1317 END DO
1318 IF (dft_control%do_admm) THEN
1319 ! matrix_ks_orb => matrix_ks_aux_fit_im
1320 ns = SIZE(matrix_ks_aux_fit_im)
1321 matrix_ks_orb(1:ns, 1:1) => matrix_ks_aux_fit_im(1:ns)
1322 DO ispin = 1, nspins
1323 CALL dbcsr_set(matrix_ks_aux_fit_im(ispin)%matrix, 0.0_dp)
1324 END DO
1325 ELSE
1326 ! matrix_ks_orb => matrix_ks_im
1327 ns = SIZE(matrix_ks_im)
1328 matrix_ks_orb(1:ns, 1:1) => matrix_ks_im(1:ns)
1329 END IF
1330
1331 CALL qs_rho_get(rho_orb, rho_ao_im=rho_ao_1d)
1332 ns = SIZE(rho_ao_1d)
1333 rho_ao_orb(1:ns, 1:1) => rho_ao_1d(1:ns)
1334
1335 ehfxrt = 0.0_dp
1336
1337 IF (x_data(irep, 1)%do_hfx_ri) THEN
1338 CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_orb, ehfx, &
1339 mo_array, rho_ao_orb, &
1340 .false., nspins, &
1341 x_data(irep, 1)%general_parameter%fraction)
1342 IF (dft_control%do_admm) THEN
1343 !for ADMMS, we need the exchange matrix k(d) for both spins
1344 DO ispin = 1, nspins
1345 CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_orb(ispin, 1)%matrix, &
1346 name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1347 END DO
1348 END IF
1349
1350 ELSE
1351 DO ispin = 1, mspin
1352 CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
1353 para_env, .false., irep, distribute_fock_matrix, &
1354 ispin=ispin)
1355 ehfxrt = ehfxrt + eh1
1356 END DO
1357 END IF
1358
1359 IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN
1360 NULLIFY (rho_ao_resp)
1361
1362 IF (x_data(irep, 1)%do_hfx_ri) THEN
1363
1364 CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
1365 x_data(irep, 1)%general_parameter%fraction, &
1366 rho_ao=rho_ao_orb, mos=mo_array, &
1367 use_virial=use_virial)
1368
1369 ELSE
1370 CALL derivatives_four_center(qs_env, rho_ao_orb, rho_ao_resp, hfx_sections, &
1371 para_env, irep, use_virial)
1372 END IF
1373 END IF
1374
1375 !! If required, the calculation of the forces will be done later with adiabatic rescaling
1376 IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx + ehfxrt
1377
1378 IF (dft_control%rtp_control%velocity_gauge) THEN
1379 cpassert(ASSOCIATED(matrix_h_im))
1380 DO ispin = 1, nspins
1381 CALL dbcsr_add(matrix_ks_im(ispin)%matrix, matrix_h_im(1, 1)%matrix, &
1382 1.0_dp, 1.0_dp)
1383 END DO
1384 END IF
1385
1386 END IF
1387
1388 IF (.NOT. qs_env%run_rtp) THEN
1389 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1390 poisson_env=poisson_env)
1391 eh1 = ehfx - eold
1392 CALL pw_hfx(qs_env, eh1, hfx_sections, poisson_env, auxbas_pw_pool, irep)
1393 eold = ehfx
1394 END IF
1395
1396 END DO
1397
1398 ! *** Set the total HFX energy
1399 energy%ex = ehfx + ehfxrt
1400
1401 ! *** Add Core-Hamiltonian-Matrix ***
1402 DO ispin = 1, nspins
1403 DO img = 1, nimages
1404 CALL dbcsr_add(matrix_ks(ispin, img)%matrix, matrix_h(1, img)%matrix, &
1405 1.0_dp, 1.0_dp)
1406 END DO
1407 END DO
1408 IF (use_virial .AND. calculate_forces) THEN
1409 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1410 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1411 virial%pv_calculate = .false.
1412 END IF
1413
1414 !! If we perform adiabatic rescaling we are now able to rescale the xc-potential
1415 IF (do_adiabatic_rescaling) THEN
1416 CALL rescale_xc_potential(qs_env, matrix_ks, rho, energy, v_rspace_new, v_tau_rspace, &
1417 hf_energy, just_energy, calculate_forces, use_virial)
1418 END IF ! do_adiabatic_rescaling
1419
1420 !update the hfx aux_fit matrixIF (dft_control%do_admm) THEN
1421 IF (dft_control%do_admm) THEN
1422 DO ispin = 1, nspins
1423 CALL dbcsr_add(matrix_ks_aux_fit_hfx(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, &
1424 0.0_dp, 1.0_dp)
1425 END DO
1426 END IF
1427
1428 CALL timestop(handle)
1429
1430 END SUBROUTINE hfx_ks_matrix
1431
1432! **************************************************************************************************
1433!> \brief This routine modifies the xc section depending on the potential type
1434!> used for the HF exchange and the resulting correction term. Currently
1435!> three types of corrections are implemented:
1436!>
1437!> coulomb: Ex,hf = Ex,hf' + (PBEx-PBEx')
1438!> shortrange: Ex,hf = Ex,hf' + (XWPBEX-XWPBEX')
1439!> truncated: Ex,hf = Ex,hf' + ( (XWPBEX0-PBE_HOLE_TC_LR) -(XWPBEX0-PBE_HOLE_TC_LR)' )
1440!>
1441!> with ' denoting the auxiliary basis set and
1442!>
1443!> PBEx: PBE exchange functional
1444!> XWPBEX: PBE exchange hole for short-range potential (erfc(omega*r)/r)
1445!> XWPBEX0: PBE exchange hole for standard coulomb potential
1446!> PBE_HOLE_TC_LR: PBE exchange hole for longrange truncated coulomb potential
1447!>
1448!> Above explanation is correct for the deafult case. If a specific functional is requested
1449!> for the correction term (cfun), we get
1450!> Ex,hf = Ex,hf' + (cfun-cfun')
1451!> for all cases of operators.
1452!>
1453!> \param x_data ...
1454!> \param xc_section the original xc_section
1455!> \param admm_env the ADMM environment
1456!> \par History
1457!> 12.2009 created [Manuel Guidon]
1458!> 05.2021 simplify for case of no correction [JGH]
1459!> \author Manuel Guidon
1460! **************************************************************************************************
1461 SUBROUTINE create_admm_xc_section(x_data, xc_section, admm_env)
1462 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1463 TYPE(section_vals_type), POINTER :: xc_section
1464 TYPE(admm_type), POINTER :: admm_env
1465
1466 LOGICAL, PARAMETER :: debug_functional = .false.
1467#if defined (__LIBXC)
1468 REAL(kind=dp), PARAMETER :: x_factor_c = 0.930525736349100025_dp
1469#endif
1470
1471 CHARACTER(LEN=20) :: name_x_func
1472 INTEGER :: hfx_potential_type, ifun, iounit, nfun
1473 LOGICAL :: funct_found
1474 REAL(dp) :: cutoff_radius, hfx_fraction, omega, &
1475 scale_coulomb, scale_longrange, scale_x
1476 TYPE(cp_logger_type), POINTER :: logger
1477 TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_section
1478
1479 logger => cp_get_default_logger()
1480 NULLIFY (admm_env%xc_section_aux, admm_env%xc_section_primary)
1481
1482 !! ** Duplicate existing xc-section
1483 CALL section_vals_duplicate(xc_section, admm_env%xc_section_aux)
1484 CALL section_vals_duplicate(xc_section, admm_env%xc_section_primary)
1485 !** Now modify the auxiliary basis
1486 !** First remove all functionals
1487 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL")
1488
1489 !* Overwrite possible shortcut
1490 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1491 i_val=xc_funct_no_shortcut)
1492
1493 !** Get number of Functionals in the list
1494 ifun = 0
1495 nfun = 0
1496 DO
1497 ifun = ifun + 1
1498 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1499 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1500 nfun = nfun + 1
1501 END DO
1502
1503 ifun = 0
1504 DO ifun = 1, nfun
1505 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=1)
1506 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1507 CALL section_vals_remove_values(xc_fun)
1508 END DO
1509
1510 IF (ASSOCIATED(x_data)) THEN
1511 hfx_potential_type = x_data(1, 1)%potential_parameter%potential_type
1512 hfx_fraction = x_data(1, 1)%general_parameter%fraction
1513 ELSE
1514 cpwarn("ADMM requested without a DFT%XC%HF section. It will be ignored for the SCF.")
1515 admm_env%aux_exch_func = do_admm_aux_exch_func_none
1516 END IF
1517
1518 !in case of no admm exchange corr., no auxiliary exchange functional needed
1519 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
1520 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1521 i_val=xc_none)
1522 hfx_fraction = 0.0_dp
1523 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default) THEN
1524 ! default PBE Functional
1525 !! ** Add functionals evaluated with auxiliary basis
1526 SELECT CASE (hfx_potential_type)
1527 CASE (do_potential_coulomb)
1528 CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", &
1529 l_val=.true.)
1530 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1531 r_val=-hfx_fraction)
1532 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", &
1533 r_val=0.0_dp)
1534 CASE (do_potential_short)
1535 omega = x_data(1, 1)%potential_parameter%omega
1536 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1537 l_val=.true.)
1538 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1539 r_val=-hfx_fraction)
1540 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1541 r_val=0.0_dp)
1542 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1543 r_val=omega)
1544 CASE (do_potential_truncated)
1545 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1546 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1547 l_val=.true.)
1548 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1549 r_val=hfx_fraction)
1550 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1551 r_val=cutoff_radius)
1552 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1553 l_val=.true.)
1554 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1555 r_val=0.0_dp)
1556 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1557 r_val=-hfx_fraction)
1558 CASE (do_potential_long)
1559 omega = x_data(1, 1)%potential_parameter%omega
1560 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1561 l_val=.true.)
1562 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1563 r_val=hfx_fraction)
1564 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1565 r_val=-hfx_fraction)
1566 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1567 r_val=omega)
1568 CASE (do_potential_mix_cl)
1569 omega = x_data(1, 1)%potential_parameter%omega
1570 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1571 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1572 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1573 l_val=.true.)
1574 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1575 r_val=hfx_fraction*scale_longrange)
1576 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1577 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1578 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1579 r_val=omega)
1580 CASE (do_potential_mix_cl_trunc)
1581 omega = x_data(1, 1)%potential_parameter%omega
1582 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1583 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1584 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1585 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1586 l_val=.true.)
1587 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1588 r_val=hfx_fraction*(scale_longrange + scale_coulomb))
1589 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1590 r_val=cutoff_radius)
1591 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1592 l_val=.true.)
1593 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1594 r_val=hfx_fraction*scale_longrange)
1595 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1596 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1597 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1598 r_val=omega)
1599 CASE DEFAULT
1600 cpabort("Unknown potential operator!")
1601 END SELECT
1602
1603 !** Now modify the functionals for the primary basis
1604 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
1605 !* Overwrite possible shortcut
1606 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1607 i_val=xc_funct_no_shortcut)
1608
1609 SELECT CASE (hfx_potential_type)
1610 CASE (do_potential_coulomb)
1611 ifun = 0
1612 funct_found = .false.
1613 DO
1614 ifun = ifun + 1
1615 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1616 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1617 IF (xc_fun%section%name == "PBE") THEN
1618 funct_found = .true.
1619 END IF
1620 END DO
1621 IF (.NOT. funct_found) THEN
1622 CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", &
1623 l_val=.true.)
1624 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1625 r_val=hfx_fraction)
1626 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", &
1627 r_val=0.0_dp)
1628 ELSE
1629 CALL section_vals_val_get(xc_fun_section, "PBE%SCALE_X", &
1630 r_val=scale_x)
1631 scale_x = scale_x + hfx_fraction
1632 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", &
1633 r_val=scale_x)
1634 END IF
1635 CASE (do_potential_short)
1636 omega = x_data(1, 1)%potential_parameter%omega
1637 ifun = 0
1638 funct_found = .false.
1639 DO
1640 ifun = ifun + 1
1641 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1642 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1643 IF (xc_fun%section%name == "XWPBE") THEN
1644 funct_found = .true.
1645 END IF
1646 END DO
1647 IF (.NOT. funct_found) THEN
1648 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1649 l_val=.true.)
1650 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1651 r_val=hfx_fraction)
1652 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1653 r_val=0.0_dp)
1654 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1655 r_val=omega)
1656 ELSE
1657 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1658 r_val=scale_x)
1659 scale_x = scale_x + hfx_fraction
1660 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1661 r_val=scale_x)
1662 END IF
1663 CASE (do_potential_long)
1664 omega = x_data(1, 1)%potential_parameter%omega
1665 ifun = 0
1666 funct_found = .false.
1667 DO
1668 ifun = ifun + 1
1669 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1670 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1671 IF (xc_fun%section%name == "XWPBE") THEN
1672 funct_found = .true.
1673 END IF
1674 END DO
1675 IF (.NOT. funct_found) THEN
1676 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1677 l_val=.true.)
1678 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1679 r_val=-hfx_fraction)
1680 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1681 r_val=hfx_fraction)
1682 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1683 r_val=omega)
1684 ELSE
1685 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1686 r_val=scale_x)
1687 scale_x = scale_x - hfx_fraction
1688 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1689 r_val=scale_x)
1690 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1691 r_val=scale_x)
1692 scale_x = scale_x + hfx_fraction
1693 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1694 r_val=scale_x)
1695
1696 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1697 r_val=omega)
1698 END IF
1699 CASE (do_potential_truncated)
1700 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1701 ifun = 0
1702 funct_found = .false.
1703 DO
1704 ifun = ifun + 1
1705 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1706 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1707 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
1708 funct_found = .true.
1709 END IF
1710 END DO
1711 IF (.NOT. funct_found) THEN
1712 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1713 l_val=.true.)
1714 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1715 r_val=-hfx_fraction)
1716 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1717 r_val=cutoff_radius)
1718 ELSE
1719 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1720 r_val=scale_x)
1721 scale_x = scale_x - hfx_fraction
1722 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1723 r_val=scale_x)
1724 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1725 r_val=cutoff_radius)
1726 END IF
1727 ifun = 0
1728 funct_found = .false.
1729 DO
1730 ifun = ifun + 1
1731 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1732 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1733 IF (xc_fun%section%name == "XWPBE") THEN
1734 funct_found = .true.
1735 END IF
1736 END DO
1737 IF (.NOT. funct_found) THEN
1738 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1739 l_val=.true.)
1740 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1741 r_val=hfx_fraction)
1742 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1743 r_val=0.0_dp)
1744
1745 ELSE
1746 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1747 r_val=scale_x)
1748 scale_x = scale_x + hfx_fraction
1749 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1750 r_val=scale_x)
1751 END IF
1752 CASE (do_potential_mix_cl_trunc)
1753 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1754 omega = x_data(1, 1)%potential_parameter%omega
1755 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1756 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1757 ifun = 0
1758 funct_found = .false.
1759 DO
1760 ifun = ifun + 1
1761 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1762 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1763 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
1764 funct_found = .true.
1765 END IF
1766 END DO
1767 IF (.NOT. funct_found) THEN
1768 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1769 l_val=.true.)
1770 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1771 r_val=-hfx_fraction*(scale_coulomb + scale_longrange))
1772 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1773 r_val=cutoff_radius)
1774
1775 ELSE
1776 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1777 r_val=scale_x)
1778 scale_x = scale_x - hfx_fraction*(scale_coulomb + scale_longrange)
1779 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1780 r_val=scale_x)
1781 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1782 r_val=cutoff_radius)
1783 END IF
1784 ifun = 0
1785 funct_found = .false.
1786 DO
1787 ifun = ifun + 1
1788 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1789 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1790 IF (xc_fun%section%name == "XWPBE") THEN
1791 funct_found = .true.
1792 END IF
1793 END DO
1794 IF (.NOT. funct_found) THEN
1795 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1796 l_val=.true.)
1797 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1798 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
1799 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1800 r_val=-hfx_fraction*scale_longrange)
1801 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1802 r_val=omega)
1803
1804 ELSE
1805 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1806 r_val=scale_x)
1807 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
1808 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1809 r_val=scale_x)
1810 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1811 r_val=scale_x)
1812 scale_x = scale_x - hfx_fraction*scale_longrange
1813 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1814 r_val=scale_x)
1815
1816 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1817 r_val=omega)
1818 END IF
1819 CASE (do_potential_mix_cl)
1820 omega = x_data(1, 1)%potential_parameter%omega
1821 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1822 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1823 ifun = 0
1824 funct_found = .false.
1825 DO
1826 ifun = ifun + 1
1827 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1828 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1829 IF (xc_fun%section%name == "XWPBE") THEN
1830 funct_found = .true.
1831 END IF
1832 END DO
1833 IF (.NOT. funct_found) THEN
1834 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", &
1835 l_val=.true.)
1836 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1837 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
1838 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1839 r_val=-hfx_fraction*scale_longrange)
1840 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1841 r_val=omega)
1842
1843 ELSE
1844 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", &
1845 r_val=scale_x)
1846 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
1847 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", &
1848 r_val=scale_x)
1849
1850 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", &
1851 r_val=scale_x)
1852 scale_x = scale_x - hfx_fraction*scale_longrange
1853 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", &
1854 r_val=scale_x)
1855
1856 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", &
1857 r_val=omega)
1858 END IF
1859 END SELECT
1860 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default_libxc) THEN
1861 ! default PBE Functional
1862 !! ** Add functionals evaluated with auxiliary basis
1863#if defined (__LIBXC)
1864 SELECT CASE (hfx_potential_type)
1865 CASE (do_potential_coulomb)
1866 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1867 l_val=.true.)
1868 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1869 r_val=-hfx_fraction)
1870 CASE (do_potential_short)
1871 omega = x_data(1, 1)%potential_parameter%omega
1872 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1873 l_val=.true.)
1874 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1875 r_val=-hfx_fraction)
1876 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1877 r_val=omega)
1878 CASE (do_potential_truncated)
1879 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1880 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1881 l_val=.true.)
1882 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1883 r_val=hfx_fraction)
1884 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1885 r_val=cutoff_radius)
1886 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1887 l_val=.true.)
1888 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1889 r_val=-hfx_fraction)
1890 CASE (do_potential_long)
1891 omega = x_data(1, 1)%potential_parameter%omega
1892 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1893 l_val=.true.)
1894 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1895 r_val=hfx_fraction)
1896 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1897 r_val=omega)
1898 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1899 l_val=.true.)
1900 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1901 r_val=-hfx_fraction)
1902 CASE (do_potential_mix_cl)
1903 omega = x_data(1, 1)%potential_parameter%omega
1904 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1905 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1906 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1907 l_val=.true.)
1908 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1909 r_val=hfx_fraction*scale_longrange)
1910 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1911 r_val=omega)
1912 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1913 l_val=.true.)
1914 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1915 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1916 CASE (do_potential_mix_cl_trunc)
1917 omega = x_data(1, 1)%potential_parameter%omega
1918 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
1919 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
1920 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
1921 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
1922 l_val=.true.)
1923 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
1924 r_val=hfx_fraction*(scale_longrange + scale_coulomb))
1925 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
1926 r_val=cutoff_radius)
1927 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1928 l_val=.true.)
1929 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1930 r_val=hfx_fraction*scale_longrange)
1931 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1932 r_val=omega)
1933 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1934 l_val=.true.)
1935 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1936 r_val=-hfx_fraction*(scale_longrange + scale_coulomb))
1937 CASE DEFAULT
1938 cpabort("Unknown potential operator!")
1939 END SELECT
1940
1941 !** Now modify the functionals for the primary basis
1942 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
1943 !* Overwrite possible shortcut
1944 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
1945 i_val=xc_funct_no_shortcut)
1946
1947 SELECT CASE (hfx_potential_type)
1948 CASE (do_potential_coulomb)
1949 ifun = 0
1950 funct_found = .false.
1951 DO
1952 ifun = ifun + 1
1953 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1954 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1955 IF (xc_fun%section%name == "GGA_X_PBE") THEN
1956 funct_found = .true.
1957 END IF
1958 END DO
1959 IF (.NOT. funct_found) THEN
1960 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
1961 l_val=.true.)
1962 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1963 r_val=hfx_fraction)
1964 ELSE
1965 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
1966 r_val=scale_x)
1967 scale_x = scale_x + hfx_fraction
1968 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
1969 r_val=scale_x)
1970 END IF
1971 CASE (do_potential_short)
1972 omega = x_data(1, 1)%potential_parameter%omega
1973 ifun = 0
1974 funct_found = .false.
1975 DO
1976 ifun = ifun + 1
1977 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
1978 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
1979 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
1980 funct_found = .true.
1981 END IF
1982 END DO
1983 IF (.NOT. funct_found) THEN
1984 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
1985 l_val=.true.)
1986 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1987 r_val=hfx_fraction)
1988 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
1989 r_val=omega)
1990 ELSE
1991 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1992 r_val=scale_x)
1993 scale_x = scale_x + hfx_fraction
1994 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
1995 r_val=scale_x)
1996 END IF
1997 CASE (do_potential_long)
1998 omega = x_data(1, 1)%potential_parameter%omega
1999 ifun = 0
2000 funct_found = .false.
2001 DO
2002 ifun = ifun + 1
2003 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2004 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2005 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2006 funct_found = .true.
2007 END IF
2008 END DO
2009 IF (.NOT. funct_found) THEN
2010 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2011 l_val=.true.)
2012 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2013 r_val=-hfx_fraction)
2014 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2015 r_val=omega)
2016 ELSE
2017 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2018 r_val=scale_x)
2019 scale_x = scale_x - hfx_fraction
2020 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2021 r_val=scale_x)
2022
2023 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2024 r_val=omega)
2025 END IF
2026 ifun = 0
2027 funct_found = .false.
2028 DO
2029 ifun = ifun + 1
2030 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2031 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2032 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2033 funct_found = .true.
2034 END IF
2035 END DO
2036 IF (.NOT. funct_found) THEN
2037 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2038 l_val=.true.)
2039 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2040 r_val=hfx_fraction)
2041 ELSE
2042 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2043 r_val=scale_x)
2044 scale_x = scale_x + hfx_fraction
2045 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2046 r_val=scale_x)
2047 END IF
2048 CASE (do_potential_truncated)
2049 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
2050 ifun = 0
2051 funct_found = .false.
2052 DO
2053 ifun = ifun + 1
2054 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2055 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2056 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
2057 funct_found = .true.
2058 END IF
2059 END DO
2060 IF (.NOT. funct_found) THEN
2061 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
2062 l_val=.true.)
2063 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2064 r_val=-hfx_fraction)
2065 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2066 r_val=cutoff_radius)
2067
2068 ELSE
2069 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2070 r_val=scale_x)
2071 scale_x = scale_x - hfx_fraction
2072 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2073 r_val=scale_x)
2074 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2075 r_val=cutoff_radius)
2076 END IF
2077 ifun = 0
2078 funct_found = .false.
2079 DO
2080 ifun = ifun + 1
2081 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2082 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2083 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2084 funct_found = .true.
2085 END IF
2086 END DO
2087 IF (.NOT. funct_found) THEN
2088 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2089 l_val=.true.)
2090 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2091 r_val=hfx_fraction)
2092
2093 ELSE
2094 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2095 r_val=scale_x)
2096 scale_x = scale_x + hfx_fraction
2097 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2098 r_val=scale_x)
2099 END IF
2100 CASE (do_potential_mix_cl_trunc)
2101 cutoff_radius = x_data(1, 1)%potential_parameter%cutoff_radius
2102 omega = x_data(1, 1)%potential_parameter%omega
2103 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
2104 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
2105 ifun = 0
2106 funct_found = .false.
2107 DO
2108 ifun = ifun + 1
2109 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2110 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2111 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN
2112 funct_found = .true.
2113 END IF
2114 END DO
2115 IF (.NOT. funct_found) THEN
2116 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", &
2117 l_val=.true.)
2118 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2119 r_val=-hfx_fraction*(scale_coulomb + scale_longrange))
2120 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2121 r_val=cutoff_radius)
2122
2123 ELSE
2124 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2125 r_val=scale_x)
2126 scale_x = scale_x - hfx_fraction*(scale_coulomb + scale_longrange)
2127 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", &
2128 r_val=scale_x)
2129 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", &
2130 r_val=cutoff_radius)
2131 END IF
2132 ifun = 0
2133 funct_found = .false.
2134 DO
2135 ifun = ifun + 1
2136 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2137 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2138 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2139 funct_found = .true.
2140 END IF
2141 END DO
2142 IF (.NOT. funct_found) THEN
2143 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2144 l_val=.true.)
2145 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2146 r_val=-hfx_fraction*scale_longrange)
2147 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2148 r_val=omega)
2149
2150 ELSE
2151 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2152 r_val=scale_x)
2153 scale_x = scale_x - hfx_fraction*scale_longrange
2154 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2155 r_val=scale_x)
2156
2157 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2158 r_val=omega)
2159 END IF
2160 ifun = 0
2161 funct_found = .false.
2162 DO
2163 ifun = ifun + 1
2164 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2165 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2166 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2167 funct_found = .true.
2168 END IF
2169 END DO
2170 IF (.NOT. funct_found) THEN
2171 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2172 l_val=.true.)
2173 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2174 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
2175 ELSE
2176 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2177 r_val=scale_x)
2178 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
2179 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2180 r_val=scale_x)
2181 END IF
2182 CASE (do_potential_mix_cl)
2183 omega = x_data(1, 1)%potential_parameter%omega
2184 scale_coulomb = x_data(1, 1)%potential_parameter%scale_coulomb
2185 scale_longrange = x_data(1, 1)%potential_parameter%scale_longrange
2186 ifun = 0
2187 funct_found = .false.
2188 DO
2189 ifun = ifun + 1
2190 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2191 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2192 IF (xc_fun%section%name == "GGA_X_WPBEH") THEN
2193 funct_found = .true.
2194 END IF
2195 END DO
2196 IF (.NOT. funct_found) THEN
2197 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_SECTION_PARAMETERS_", &
2198 l_val=.true.)
2199 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2200 r_val=-hfx_fraction*scale_longrange)
2201 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2202 r_val=omega)
2203
2204 ELSE
2205 CALL section_vals_val_get(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2206 r_val=scale_x)
2207 scale_x = scale_x - hfx_fraction*scale_longrange
2208 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%SCALE", &
2209 r_val=scale_x)
2210
2211 CALL section_vals_val_set(xc_fun_section, "GGA_X_WPBEH%_OMEGA", &
2212 r_val=omega)
2213 END IF
2214 ifun = 0
2215 funct_found = .false.
2216 DO
2217 ifun = ifun + 1
2218 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2219 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2220 IF (xc_fun%section%name == "GGA_X_PBE") THEN
2221 funct_found = .true.
2222 END IF
2223 END DO
2224 IF (.NOT. funct_found) THEN
2225 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%_SECTION_PARAMETERS_", &
2226 l_val=.true.)
2227 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2228 r_val=hfx_fraction*(scale_coulomb + scale_longrange))
2229 ELSE
2230 CALL section_vals_val_get(xc_fun_section, "GGA_X_PBE%SCALE", &
2231 r_val=scale_x)
2232 scale_x = scale_x + hfx_fraction*(scale_coulomb + scale_longrange)
2233 CALL section_vals_val_set(xc_fun_section, "GGA_X_PBE%SCALE", &
2234 r_val=scale_x)
2235 END IF
2236 END SELECT
2237#else
2238 CALL cp_abort(__location__, "In order use a LibXC-based ADMM "// &
2239 "exchange correction functionals, you have to compile and link against LibXC!")
2240#endif
2241
2242 ! PBEX (always bare form), OPTX and Becke88 functional
2243 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex .OR. &
2244 admm_env%aux_exch_func == do_admm_aux_exch_func_opt .OR. &
2245 admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN
2246 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2247 name_x_func = 'PBE'
2248 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2249 name_x_func = 'OPTX'
2250 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN
2251 name_x_func = 'BECKE88'
2252 END IF
2253 !primary basis
2254 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2255 l_val=.true.)
2256 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2257 r_val=-hfx_fraction)
2258
2259 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2260 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_C", r_val=0.0_dp)
2261 END IF
2262
2263 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2264 IF (admm_env%aux_exch_func_param) THEN
2265 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A1", &
2266 r_val=admm_env%aux_x_param(1))
2267 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A2", &
2268 r_val=admm_env%aux_x_param(2))
2269 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%GAMMA", &
2270 r_val=admm_env%aux_x_param(3))
2271 END IF
2272 END IF
2273
2274 !** Now modify the functionals for the primary basis
2275 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2276 !* Overwrite possible L")
2277 !* Overwrite possible shortcut
2278 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
2279 i_val=xc_funct_no_shortcut)
2280
2281 ifun = 0
2282 funct_found = .false.
2283 DO
2284 ifun = ifun + 1
2285 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2286 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2287 IF (xc_fun%section%name == trim(name_x_func)) THEN
2288 funct_found = .true.
2289 END IF
2290 END DO
2291 IF (.NOT. funct_found) THEN
2292 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2293 l_val=.true.)
2294 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2295 r_val=hfx_fraction)
2296 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN
2297 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_C", &
2298 r_val=0.0_dp)
2299 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2300 IF (admm_env%aux_exch_func_param) THEN
2301 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A1", &
2302 r_val=admm_env%aux_x_param(1))
2303 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%A2", &
2304 r_val=admm_env%aux_x_param(2))
2305 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%GAMMA", &
2306 r_val=admm_env%aux_x_param(3))
2307 END IF
2308 END IF
2309
2310 ELSE
2311 CALL section_vals_val_get(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2312 r_val=scale_x)
2313 scale_x = scale_x + hfx_fraction
2314 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE_X", &
2315 r_val=scale_x)
2316 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN
2317 cpassert(.NOT. admm_env%aux_exch_func_param)
2318 END IF
2319 END IF
2320
2321 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc .OR. &
2322 admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc .OR. &
2323 admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc .OR. &
2324 admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
2325#if defined(__LIBXC)
2326 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex_libxc) THEN
2327 name_x_func = 'GGA_X_PBE'
2328 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2329 name_x_func = 'GGA_X_OPTX'
2330 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee_libxc) THEN
2331 name_x_func = 'GGA_X_B88'
2332 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_sx_libxc) THEN
2333 name_x_func = 'LDA_X'
2334 END IF
2335 !primary basis
2336 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2337 l_val=.true.)
2338 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE", &
2339 r_val=-hfx_fraction)
2340
2341 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2342 IF (admm_env%aux_exch_func_param) THEN
2343 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_A", &
2344 r_val=admm_env%aux_x_param(1))
2345 ! LibXC rescales the second parameter of the OPTX functional (see documentation there)
2346 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_B", &
2347 r_val=admm_env%aux_x_param(2)/x_factor_c)
2348 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_GAMMA", &
2349 r_val=admm_env%aux_x_param(3))
2350 END IF
2351 END IF
2352
2353 !** Now modify the functionals for the primary basis
2354 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2355 !* Overwrite possible L")
2356 !* Overwrite possible shortcut
2357 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", &
2358 i_val=xc_funct_no_shortcut)
2359
2360 ifun = 0
2361 funct_found = .false.
2362 DO
2363 ifun = ifun + 1
2364 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2365 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2366 IF (xc_fun%section%name == trim(name_x_func)) THEN
2367 funct_found = .true.
2368 END IF
2369 END DO
2370 IF (.NOT. funct_found) THEN
2371 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_SECTION_PARAMETERS_", &
2372 l_val=.true.)
2373 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE", &
2374 r_val=hfx_fraction)
2375 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2376 IF (admm_env%aux_exch_func_param) THEN
2377 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_A", &
2378 r_val=admm_env%aux_x_param(1))
2379 ! LibXC rescales the second parameter of the OPTX functional (see documentation there)
2380 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_B", &
2381 r_val=admm_env%aux_x_param(2)/x_factor_c)
2382 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%_GAMMA", &
2383 r_val=admm_env%aux_x_param(3))
2384 END IF
2385 END IF
2386
2387 ELSE
2388 CALL section_vals_val_get(xc_fun_section, trim(name_x_func)//"%SCALE", &
2389 r_val=scale_x)
2390 scale_x = scale_x + hfx_fraction
2391 CALL section_vals_val_set(xc_fun_section, trim(name_x_func)//"%SCALE", &
2392 r_val=scale_x)
2393 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt_libxc) THEN
2394 cpassert(.NOT. admm_env%aux_exch_func_param)
2395 END IF
2396 END IF
2397#else
2398 CALL cp_abort(__location__, "In order use a LibXC-based ADMM "// &
2399 "exchange correction functionals, you have to compile and link against LibXC!")
2400#endif
2401
2402 ELSE
2403 cpabort("Unknown exchange correction functional!")
2404 END IF
2405
2406 IF (debug_functional) THEN
2407 iounit = cp_logger_get_default_io_unit(logger)
2408 IF (iounit > 0) THEN
2409 WRITE (iounit, "(A)") " ADMM Primary Basis Set Functional"
2410 END IF
2411 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL")
2412 ifun = 0
2413 funct_found = .false.
2414 DO
2415 ifun = ifun + 1
2416 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2417 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2418
2419 scale_x = -1000.0_dp
2420 IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
2421 CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
2422 END IF
2423 IF (xc_fun%section%name == "XWPBE") THEN
2424 CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
2425 IF (iounit > 0) THEN
2426 WRITE (iounit, "(T5,A,T25,2F10.3)") trim(xc_fun%section%name), scale_x, hfx_fraction
2427 END IF
2428 ELSE
2429 IF (iounit > 0) THEN
2430 WRITE (iounit, "(T5,A,T25,F10.3)") trim(xc_fun%section%name), scale_x
2431 END IF
2432 END IF
2433 END DO
2434
2435 IF (iounit > 0) THEN
2436 WRITE (iounit, "(A)") " Auxiliary Basis Set Functional"
2437 END IF
2438 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL")
2439 ifun = 0
2440 funct_found = .false.
2441 DO
2442 ifun = ifun + 1
2443 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun)
2444 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
2445 scale_x = -1000.0_dp
2446 IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN
2447 CALL section_vals_val_get(xc_fun, "SCALE_X", r_val=scale_x)
2448 END IF
2449 IF (xc_fun%section%name == "XWPBE") THEN
2450 CALL section_vals_val_get(xc_fun, "SCALE_X0", r_val=hfx_fraction)
2451 IF (iounit > 0) THEN
2452 WRITE (iounit, "(T5,A,T25,2F10.3)") trim(xc_fun%section%name), scale_x, hfx_fraction
2453 END IF
2454 ELSE
2455 IF (iounit > 0) THEN
2456 WRITE (iounit, "(T5,A,T25,F10.3)") trim(xc_fun%section%name), scale_x
2457 END IF
2458 END IF
2459 END DO
2460 END IF
2461
2462 END SUBROUTINE create_admm_xc_section
2463
2464! **************************************************************************************************
2465!> \brief Add the hfx contributions to the Hamiltonian
2466!>
2467!> \param matrix_ks Kohn-Sham matrix (updated on exit)
2468!> \param rho_ao electron density expressed in terms of atomic orbitals
2469!> \param qs_env Quickstep environment
2470!> \param update_energy whether to update energy (default: yes)
2471!> \param recalc_integrals whether to recalculate integrals (default: value of HF%TREAT_LSD_IN_CORE)
2472!> \param external_hfx_sections ...
2473!> \param external_x_data ...
2474!> \param external_para_env ...
2475!> \note
2476!> Simplified version of subroutine hfx_ks_matrix()
2477! **************************************************************************************************
2478 SUBROUTINE tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, &
2479 external_hfx_sections, external_x_data, external_para_env)
2480 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
2481 TARGET :: matrix_ks, rho_ao
2482 TYPE(qs_environment_type), POINTER :: qs_env
2483 LOGICAL, INTENT(IN), OPTIONAL :: update_energy, recalc_integrals
2484 TYPE(section_vals_type), OPTIONAL, POINTER :: external_hfx_sections
2485 TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, TARGET :: external_x_data
2486 TYPE(mp_para_env_type), OPTIONAL, POINTER :: external_para_env
2487
2488 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddft_hfx_matrix'
2489
2490 INTEGER :: handle, irep, ispin, mspin, n_rep_hf, &
2491 nspins
2492 LOGICAL :: distribute_fock_matrix, &
2493 hfx_treat_lsd_in_core, &
2494 my_update_energy, s_mstruct_changed
2495 REAL(kind=dp) :: eh1, ehfx
2496 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp
2497 TYPE(dft_control_type), POINTER :: dft_control
2498 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
2499 TYPE(mp_para_env_type), POINTER :: para_env
2500 TYPE(qs_energy_type), POINTER :: energy
2501 TYPE(section_vals_type), POINTER :: hfx_sections, input
2502
2503 CALL timeset(routinen, handle)
2504
2505 NULLIFY (dft_control, hfx_sections, input, para_env, matrix_ks_kp, rho_ao_kp)
2506
2507 CALL get_qs_env(qs_env=qs_env, &
2508 dft_control=dft_control, &
2509 energy=energy, &
2510 input=input, &
2511 para_env=para_env, &
2512 s_mstruct_changed=s_mstruct_changed, &
2513 x_data=x_data)
2514
2515 ! This should probably be the HF section from the TDDFPT XC section!
2516 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
2517
2518 IF (PRESENT(external_hfx_sections)) hfx_sections => external_hfx_sections
2519 IF (PRESENT(external_x_data)) x_data => external_x_data
2520 IF (PRESENT(external_para_env)) para_env => external_para_env
2521
2522 my_update_energy = .true.
2523 IF (PRESENT(update_energy)) my_update_energy = update_energy
2524
2525 IF (PRESENT(recalc_integrals)) s_mstruct_changed = recalc_integrals
2526
2527 cpassert(dft_control%nimages == 1)
2528 nspins = dft_control%nspins
2529
2530 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
2531 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
2532 i_rep_section=1)
2533
2534 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
2535 distribute_fock_matrix = .true.
2536
2537 mspin = 1
2538 IF (hfx_treat_lsd_in_core) mspin = nspins
2539
2540 matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
2541 rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins)
2542
2543 DO irep = 1, n_rep_hf
2544 ! the real hfx calulation
2545 ehfx = 0.0_dp
2546
2547 IF (x_data(irep, 1)%do_hfx_ri) THEN
2548 CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_kp, ehfx, &
2549 rho_ao=rho_ao_kp, geometry_did_change=s_mstruct_changed, &
2550 nspins=nspins, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
2551
2552 ELSE
2553 DO ispin = 1, mspin
2554 CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
2555 s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
2556 ehfx = ehfx + eh1
2557 END DO
2558 END IF
2559 END DO
2560 IF (my_update_energy) energy%ex = ehfx
2561
2562 CALL timestop(handle)
2563 END SUBROUTINE tddft_hfx_matrix
2564
2565END 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)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Routines to calculate HFX energy and potential.
subroutine, public integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, geometry_did_change, irep, distribute_fock_matrix, ispin)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
Test routines for HFX caclulations using PW.
subroutine, public pw_hfx(qs_env, ehfx, hfx_section, poisson_env, auxbas_pw_pool, irep)
computes the Hartree-Fock energy brute force in a pw basis
RI-methods for HFX and K-points. \auhtor Augustin Bussy (01.2023)
Definition hfx_ri_kp.F:13
subroutine, public hfx_ri_update_forces_kp(qs_env, ri_data, nspins, hf_fraction, rho_ao, use_virial)
Update the K-points RI-HFX forces.
Definition hfx_ri_kp.F: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, reorder_rs_grid_ranks, skip_load_balance_distributed, soft_valid, basis_type, pw_env_external, sab_orb_external)
...
types for task lists
subroutine, public deallocate_task_list(task_list)
deallocates the components and the object itself
subroutine, public allocate_task_list(task_list)
allocates and initialised the components of the task_list_type
subroutine, public rescale_xc_potential(qs_env, ks_matrix, rho, energy, v_rspace_new, v_tau_rspace, hf_energy, just_energy, calculate_forces, use_virial)
A subtype of the admm_env that contains the extra data needed for an ADMM GAPW calculation.
Definition admm_types.F:85
stores some data used in wavefunction fitting
Definition admm_types.F:120
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
structure to store local (to a processor) ordered lists of integers.
distributes pairs on a 2d grid of processors
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F: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.