LCOV - code coverage report
Current view: top level - src - hfx_admm_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:3158929) Lines: 902 1055 85.5 %
Date: 2025-07-22 22:41:15 Functions: 11 11 100.0 %

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

Generated by: LCOV version 1.15