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