LCOV - code coverage report
Current view: top level - src - qs_ks_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:dc8eeee) Lines: 648 801 80.9 %
Date: 2025-05-15 08:34:30 Functions: 8 9 88.9 %

          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 routines that build the Kohn-Sham matrix (i.e calculate the coulomb
      10             : !>      and xc parts
      11             : !> \par History
      12             : !>      05.2002 moved from qs_scf (see there the history) [fawzi]
      13             : !>      JGH [30.08.02] multi-grid arrays independent from density and potential
      14             : !>      10.2002 introduced pools, uses updated rho as input,
      15             : !>              removed most temporary variables, renamed may vars,
      16             : !>              began conversion to LSD [fawzi]
      17             : !>      10.2004 moved calculate_w_matrix here [Joost VandeVondele]
      18             : !>              introduced energy derivative wrt MOs [Joost VandeVondele]
      19             : !> \author Fawzi Mohamed
      20             : ! **************************************************************************************************
      21             : 
      22             : MODULE qs_ks_utils
      23             :    USE admm_types,                      ONLY: admm_type,&
      24             :                                               get_admm_env
      25             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      26             :    USE cell_types,                      ONLY: cell_type
      27             :    USE cp_control_types,                ONLY: dft_control_type
      28             :    USE cp_dbcsr_api,                    ONLY: &
      29             :         dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_info, dbcsr_init_p, &
      30             :         dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_type
      31             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot,&
      32             :                                               dbcsr_scale_by_vector
      33             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      34             :                                               copy_fm_to_dbcsr,&
      35             :                                               cp_dbcsr_plus_fm_fm_t,&
      36             :                                               cp_dbcsr_sm_fm_multiply,&
      37             :                                               dbcsr_allocate_matrix_set,&
      38             :                                               dbcsr_deallocate_matrix_set
      39             :    USE cp_ddapc,                        ONLY: cp_ddapc_apply_CD
      40             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      41             :                                               cp_fm_struct_release,&
      42             :                                               cp_fm_struct_type
      43             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      44             :                                               cp_fm_get_info,&
      45             :                                               cp_fm_release,&
      46             :                                               cp_fm_set_all,&
      47             :                                               cp_fm_to_fm,&
      48             :                                               cp_fm_type
      49             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      50             :                                               cp_logger_type,&
      51             :                                               cp_to_string
      52             :    USE cp_output_handling,              ONLY: cp_p_file,&
      53             :                                               cp_print_key_finished_output,&
      54             :                                               cp_print_key_should_output,&
      55             :                                               cp_print_key_unit_nr
      56             :    USE hfx_admm_utils,                  ONLY: tddft_hfx_matrix
      57             :    USE hfx_derivatives,                 ONLY: derivatives_four_center
      58             :    USE hfx_types,                       ONLY: hfx_type
      59             :    USE input_constants,                 ONLY: &
      60             :         cdft_alpha_constraint, cdft_beta_constraint, cdft_charge_constraint, &
      61             :         cdft_magnetization_constraint, do_admm_aux_exch_func_none, do_ppl_grid, sic_ad, sic_eo, &
      62             :         sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none
      63             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      64             :                                               section_vals_type,&
      65             :                                               section_vals_val_get
      66             :    USE kahan_sum,                       ONLY: accurate_dot_product,&
      67             :                                               accurate_sum
      68             :    USE kinds,                           ONLY: default_string_length,&
      69             :                                               dp
      70             :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      71             :                                               kpoint_type
      72             :    USE lri_environment_methods,         ONLY: v_int_ppl_update
      73             :    USE lri_environment_types,           ONLY: lri_density_type,&
      74             :                                               lri_environment_type,&
      75             :                                               lri_kind_type
      76             :    USE lri_forces,                      ONLY: calculate_lri_forces,&
      77             :                                               calculate_ri_forces
      78             :    USE lri_ks_methods,                  ONLY: calculate_lri_ks_matrix,&
      79             :                                               calculate_ri_ks_matrix
      80             :    USE message_passing,                 ONLY: mp_para_env_type
      81             :    USE ps_implicit_types,               ONLY: MIXED_BC,&
      82             :                                               MIXED_PERIODIC_BC,&
      83             :                                               NEUMANN_BC,&
      84             :                                               PERIODIC_BC
      85             :    USE pw_env_types,                    ONLY: pw_env_get,&
      86             :                                               pw_env_type
      87             :    USE pw_methods,                      ONLY: pw_axpy,&
      88             :                                               pw_copy,&
      89             :                                               pw_integral_ab,&
      90             :                                               pw_integrate_function,&
      91             :                                               pw_scale,&
      92             :                                               pw_transfer,&
      93             :                                               pw_zero
      94             :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      95             :    USE pw_poisson_types,                ONLY: pw_poisson_implicit,&
      96             :                                               pw_poisson_type
      97             :    USE pw_pool_types,                   ONLY: pw_pool_type
      98             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      99             :                                               pw_r3d_rs_type
     100             :    USE qs_cdft_types,                   ONLY: cdft_control_type
     101             :    USE qs_charges_types,                ONLY: qs_charges_type
     102             :    USE qs_collocate_density,            ONLY: calculate_rho_elec
     103             :    USE qs_energy_types,                 ONLY: qs_energy_type
     104             :    USE qs_environment_types,            ONLY: get_qs_env,&
     105             :                                               qs_environment_type
     106             :    USE qs_force_types,                  ONLY: qs_force_type
     107             :    USE qs_integrate_potential,          ONLY: integrate_v_rspace,&
     108             :                                               integrate_v_rspace_diagonal,&
     109             :                                               integrate_v_rspace_one_center
     110             :    USE qs_kind_types,                   ONLY: get_qs_kind_set,&
     111             :                                               qs_kind_type
     112             :    USE qs_ks_qmmm_methods,              ONLY: qmmm_modify_hartree_pot
     113             :    USE qs_ks_types,                     ONLY: get_ks_env,&
     114             :                                               qs_ks_env_type
     115             :    USE qs_mo_types,                     ONLY: get_mo_set,&
     116             :                                               mo_set_type
     117             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
     118             :                                               qs_rho_type
     119             :    USE task_list_types,                 ONLY: task_list_type
     120             :    USE virial_types,                    ONLY: virial_type
     121             :    USE xc,                              ONLY: xc_exc_calc,&
     122             :                                               xc_vxc_pw_create1
     123             : #include "./base/base_uses.f90"
     124             : 
     125             :    IMPLICIT NONE
     126             : 
     127             :    PRIVATE
     128             : 
     129             :    LOGICAL, PARAMETER :: debug_this_module = .TRUE.
     130             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_utils'
     131             : 
     132             :    PUBLIC :: low_spin_roks, sic_explicit_orbitals, calc_v_sic_rspace, print_densities, &
     133             :              print_detailed_energy, compute_matrix_vxc, sum_up_and_integrate, &
     134             :              calculate_zmp_potential, get_embed_potential_energy
     135             : 
     136             : CONTAINS
     137             : 
     138             : ! **************************************************************************************************
     139             : !> \brief do ROKS calculations yielding low spin states
     140             : !> \param energy ...
     141             : !> \param qs_env ...
     142             : !> \param dft_control ...
     143             : !> \param do_hfx ...
     144             : !> \param just_energy ...
     145             : !> \param calculate_forces ...
     146             : !> \param auxbas_pw_pool ...
     147             : ! **************************************************************************************************
     148      102973 :    SUBROUTINE low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
     149             :                             calculate_forces, auxbas_pw_pool)
     150             : 
     151             :       TYPE(qs_energy_type), POINTER                      :: energy
     152             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     153             :       TYPE(dft_control_type), POINTER                    :: dft_control
     154             :       LOGICAL, INTENT(IN)                                :: do_hfx, just_energy, calculate_forces
     155             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     156             : 
     157             :       CHARACTER(*), PARAMETER                            :: routineN = 'low_spin_roks'
     158             : 
     159             :       INTEGER                                            :: handle, irep, ispin, iterm, k, k_alpha, &
     160             :                                                             k_beta, n_rep, Nelectron, Nspin, Nterms
     161      102973 :       INTEGER, DIMENSION(:), POINTER                     :: ivec
     162      102973 :       INTEGER, DIMENSION(:, :, :), POINTER               :: occupations
     163             :       LOGICAL                                            :: compute_virial, in_range, &
     164             :                                                             uniform_occupation
     165             :       REAL(KIND=dp)                                      :: ehfx, exc
     166             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_xc_tmp
     167      102973 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: energy_scaling, rvec, scaling
     168             :       TYPE(cell_type), POINTER                           :: cell
     169      102973 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_h, matrix_hfx, matrix_p, mdummy, &
     170      102973 :                                                             mo_derivs, rho_ao
     171      102973 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p2
     172             :       TYPE(dbcsr_type), POINTER                          :: dbcsr_deriv, fm_deriv, fm_scaled, &
     173             :                                                             mo_coeff
     174      102973 :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     175      102973 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     176             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     177      102973 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     178             :       TYPE(pw_env_type), POINTER                         :: pw_env
     179             :       TYPE(pw_pool_type), POINTER                        :: xc_pw_pool
     180             :       TYPE(pw_r3d_rs_type)                               :: work_v_rspace
     181      102973 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r, tau, vxc, vxc_tau
     182             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     183             :       TYPE(qs_rho_type), POINTER                         :: rho
     184             :       TYPE(section_vals_type), POINTER                   :: hfx_section, input, &
     185             :                                                             low_spin_roks_section, xc_section
     186             :       TYPE(virial_type), POINTER                         :: virial
     187             : 
     188      102611 :       IF (.NOT. dft_control%low_spin_roks) RETURN
     189             : 
     190         362 :       CALL timeset(routineN, handle)
     191             : 
     192         362 :       NULLIFY (ks_env, rho_ao)
     193             : 
     194             :       ! Test for not compatible options
     195         362 :       IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
     196           0 :          CALL cp_abort(__LOCATION__, "GAPW/GAPW_XC are not compatible with low spin ROKS method.")
     197             :       END IF
     198         362 :       IF (dft_control%do_admm) THEN
     199           0 :          CALL cp_abort(__LOCATION__, "ADMM not compatible with low spin ROKS method.")
     200             :       END IF
     201         362 :       IF (dft_control%do_admm) THEN
     202           0 :          IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
     203             :             CALL cp_abort(__LOCATION__, "ADMM with XC correction functional "// &
     204           0 :                           "not compatible with low spin ROKS method.")
     205             :          END IF
     206             :       END IF
     207         362 :       IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
     208             :           dft_control%qs_control%xtb) THEN
     209           0 :          CALL cp_abort(__LOCATION__, "SE/xTB/DFTB are not compatible with low spin ROKS method.")
     210             :       END IF
     211             : 
     212             :       CALL get_qs_env(qs_env, &
     213             :                       ks_env=ks_env, &
     214             :                       mo_derivs=mo_derivs, &
     215             :                       mos=mo_array, &
     216             :                       rho=rho, &
     217             :                       pw_env=pw_env, &
     218             :                       input=input, &
     219             :                       cell=cell, &
     220         362 :                       virial=virial)
     221             : 
     222         362 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     223             : 
     224         362 :       compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
     225         362 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     226         362 :       hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
     227             : 
     228             :       ! some assumptions need to be checked
     229             :       ! we have two spins
     230         362 :       CPASSERT(SIZE(mo_array, 1) == 2)
     231         362 :       Nspin = 2
     232             :       ! we want uniform occupations
     233         362 :       CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
     234         362 :       CPASSERT(uniform_occupation)
     235         362 :       CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
     236         362 :       CPASSERT(uniform_occupation)
     237         362 :       IF (do_hfx .AND. calculate_forces .AND. compute_virial) THEN
     238           0 :          CALL cp_abort(__LOCATION__, "ROKS virial with HFX not available.")
     239             :       END IF
     240             : 
     241         362 :       NULLIFY (dbcsr_deriv)
     242         362 :       CALL dbcsr_init_p(dbcsr_deriv)
     243         362 :       CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix)
     244         362 :       CALL dbcsr_set(dbcsr_deriv, 0.0_dp)
     245             : 
     246             :       ! basic info
     247         362 :       CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
     248         362 :       CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_alpha)
     249         362 :       CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff)
     250         362 :       CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_beta)
     251             : 
     252             :       ! read the input
     253         362 :       low_spin_roks_section => section_vals_get_subs_vals(input, "DFT%LOW_SPIN_ROKS")
     254             : 
     255         362 :       CALL section_vals_val_get(low_spin_roks_section, "ENERGY_SCALING", r_vals=rvec)
     256         362 :       Nterms = SIZE(rvec)
     257        1086 :       ALLOCATE (energy_scaling(Nterms))
     258        1810 :       energy_scaling = rvec !? just wondering, should this add up to 1, in which case we should cpp?
     259             : 
     260         362 :       CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", n_rep_val=n_rep)
     261         362 :       CPASSERT(n_rep == Nterms)
     262         362 :       CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec)
     263         362 :       Nelectron = SIZE(ivec)
     264         362 :       CPASSERT(Nelectron == k_alpha - k_beta)
     265        1448 :       ALLOCATE (occupations(2, Nelectron, Nterms))
     266        5430 :       occupations = 0
     267        1086 :       DO iterm = 1, Nterms
     268         724 :          CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec)
     269         724 :          CPASSERT(Nelectron == SIZE(ivec))
     270        4344 :          in_range = ALL(ivec >= 1) .AND. ALL(ivec <= 2)
     271         724 :          CPASSERT(in_range)
     272        2534 :          DO k = 1, Nelectron
     273        2172 :             occupations(ivec(k), k, iterm) = 1
     274             :          END DO
     275             :       END DO
     276             : 
     277             :       ! set up general data structures
     278             :       ! density matrices, kohn-sham matrices
     279             : 
     280         362 :       NULLIFY (matrix_p)
     281         362 :       CALL dbcsr_allocate_matrix_set(matrix_p, Nspin)
     282        1086 :       DO ispin = 1, Nspin
     283         724 :          ALLOCATE (matrix_p(ispin)%matrix)
     284             :          CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, &
     285         724 :                          name="density matrix low spin roks")
     286        1086 :          CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
     287             :       END DO
     288             : 
     289         362 :       NULLIFY (matrix_h)
     290         362 :       CALL dbcsr_allocate_matrix_set(matrix_h, Nspin)
     291        1086 :       DO ispin = 1, Nspin
     292         724 :          ALLOCATE (matrix_h(ispin)%matrix)
     293             :          CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, &
     294         724 :                          name="KS matrix low spin roks")
     295        1086 :          CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
     296             :       END DO
     297             : 
     298         362 :       IF (do_hfx) THEN
     299         220 :          NULLIFY (matrix_hfx)
     300         220 :          CALL dbcsr_allocate_matrix_set(matrix_hfx, Nspin)
     301         660 :          DO ispin = 1, Nspin
     302         440 :             ALLOCATE (matrix_hfx(ispin)%matrix)
     303             :             CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
     304         660 :                             name="HFX matrix low spin roks")
     305             :          END DO
     306             :       END IF
     307             : 
     308             :       ! grids in real and g space for rho and vxc
     309             :       ! tau functionals are not supported
     310         362 :       NULLIFY (tau, vxc_tau, vxc)
     311         362 :       CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
     312             : 
     313        1086 :       ALLOCATE (rho_r(Nspin))
     314        1086 :       ALLOCATE (rho_g(Nspin))
     315        1086 :       DO ispin = 1, Nspin
     316         724 :          CALL auxbas_pw_pool%create_pw(rho_r(ispin))
     317        1086 :          CALL auxbas_pw_pool%create_pw(rho_g(ispin))
     318             :       END DO
     319         362 :       CALL auxbas_pw_pool%create_pw(work_v_rspace)
     320             : 
     321             :       ! get mo matrices needed to construct the density matrices
     322             :       ! we will base all on the alpha spin matrix, obviously possible in ROKS
     323         362 :       CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
     324         362 :       NULLIFY (fm_scaled, fm_deriv)
     325         362 :       CALL dbcsr_init_p(fm_scaled)
     326         362 :       CALL dbcsr_init_p(fm_deriv)
     327         362 :       CALL dbcsr_copy(fm_scaled, mo_coeff)
     328         362 :       CALL dbcsr_copy(fm_deriv, mo_coeff)
     329             : 
     330        1086 :       ALLOCATE (scaling(k_alpha))
     331             : 
     332             :       ! for each term, add it with the given scaling factor to the energy, and compute the required derivatives
     333        1086 :       DO iterm = 1, Nterms
     334             : 
     335        2172 :          DO ispin = 1, Nspin
     336             :             ! compute the proper density matrices with the required occupations
     337        1448 :             CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
     338       11584 :             scaling = 1.0_dp
     339        4344 :             scaling(k_alpha - Nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
     340        1448 :             CALL dbcsr_copy(fm_scaled, mo_coeff)
     341        1448 :             CALL dbcsr_scale_by_vector(fm_scaled, scaling, side='right')
     342             :             CALL dbcsr_multiply('n', 't', 1.0_dp, mo_coeff, fm_scaled, &
     343        1448 :                                 0.0_dp, matrix_p(ispin)%matrix, retain_sparsity=.TRUE.)
     344             :             ! compute the densities on the grid
     345             :             CALL calculate_rho_elec(matrix_p=matrix_p(ispin)%matrix, &
     346             :                                     rho=rho_r(ispin), rho_gspace=rho_g(ispin), &
     347        2172 :                                     ks_env=ks_env)
     348             :          END DO
     349             : 
     350             :          ! compute the exchange energies / potential if needed
     351         724 :          IF (just_energy) THEN
     352             :             exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
     353          88 :                               pw_pool=xc_pw_pool)
     354             :          ELSE
     355         636 :             CPASSERT(.NOT. compute_virial)
     356             :             CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, &
     357             :                                    rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
     358         636 :                                    pw_pool=xc_pw_pool, compute_virial=.FALSE., virial_xc=virial_xc_tmp)
     359             :          END IF
     360             : 
     361         724 :          energy%exc = energy%exc + energy_scaling(iterm)*exc
     362             : 
     363         724 :          IF (do_hfx) THEN
     364             :             ! Add Hartree-Fock contribution
     365        1320 :             DO ispin = 1, Nspin
     366        1320 :                CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
     367             :             END DO
     368         440 :             ehfx = energy%ex
     369             :             CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, &
     370         440 :                                   recalc_integrals=.FALSE., update_energy=.TRUE.)
     371         440 :             energy%ex = ehfx + energy_scaling(iterm)*energy%ex
     372             :          END IF
     373             : 
     374             :          ! add the corresponding derivatives to the MO derivatives
     375        1086 :          IF (.NOT. just_energy) THEN
     376             :             ! get the potential in matrix form
     377        1908 :             DO ispin = 1, Nspin
     378        1272 :                CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
     379             :                ! use a work_v_rspace
     380        1272 :                CALL pw_axpy(vxc(ispin), work_v_rspace, energy_scaling(iterm)*vxc(ispin)%pw_grid%dvol, 0.0_dp)
     381             :                CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
     382        1272 :                                        qs_env=qs_env, calculate_forces=calculate_forces)
     383        1908 :                CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
     384             :             END DO
     385         636 :             DEALLOCATE (vxc)
     386             : 
     387         636 :             IF (do_hfx) THEN
     388             :                ! add HFX contribution
     389        1104 :                DO ispin = 1, Nspin
     390             :                   CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
     391        1104 :                                  1.0_dp, energy_scaling(iterm))
     392             :                END DO
     393         368 :                IF (calculate_forces) THEN
     394           8 :                   CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
     395           8 :                   IF (x_data(1, 1)%n_rep_hf /= 1) THEN
     396             :                      CALL cp_abort(__LOCATION__, "Multiple HFX section forces not compatible "// &
     397           0 :                                    "with low spin ROKS method.")
     398             :                   END IF
     399           8 :                   IF (x_data(1, 1)%do_hfx_ri) THEN
     400           0 :                      CALL cp_abort(__LOCATION__, "HFX_RI forces not compatible with low spin ROKS method.")
     401             :                   ELSE
     402           8 :                      irep = 1
     403           8 :                      NULLIFY (mdummy)
     404           8 :                      matrix_p2(1:Nspin, 1:1) => matrix_p(1:Nspin)
     405             :                      CALL derivatives_four_center(qs_env, matrix_p2, mdummy, hfx_section, para_env, &
     406             :                                                   irep, compute_virial, &
     407           8 :                                                   adiabatic_rescale_factor=energy_scaling(iterm))
     408             :                   END IF
     409             :                END IF
     410             : 
     411             :             END IF
     412             : 
     413             :             ! add this to the mo_derivs, again based on the alpha mo_coeff
     414        1908 :             DO ispin = 1, Nspin
     415             :                CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
     416        1272 :                                    0.0_dp, dbcsr_deriv, last_column=k_alpha)
     417             : 
     418       10176 :                scaling = 1.0_dp
     419        3816 :                scaling(k_alpha - Nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
     420        1272 :                CALL dbcsr_scale_by_vector(dbcsr_deriv, scaling, side='right')
     421        1908 :                CALL dbcsr_add(mo_derivs(1)%matrix, dbcsr_deriv, 1.0_dp, 1.0_dp)
     422             :             END DO
     423             : 
     424             :          END IF
     425             : 
     426             :       END DO
     427             : 
     428             :       ! release allocated memory
     429        1086 :       DO ispin = 1, Nspin
     430         724 :          CALL auxbas_pw_pool%give_back_pw(rho_r(ispin))
     431        1086 :          CALL auxbas_pw_pool%give_back_pw(rho_g(ispin))
     432             :       END DO
     433         362 :       DEALLOCATE (rho_r, rho_g)
     434         362 :       CALL dbcsr_deallocate_matrix_set(matrix_p)
     435         362 :       CALL dbcsr_deallocate_matrix_set(matrix_h)
     436         362 :       IF (do_hfx) THEN
     437         220 :          CALL dbcsr_deallocate_matrix_set(matrix_hfx)
     438             :       END IF
     439             : 
     440         362 :       CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
     441             : 
     442         362 :       CALL dbcsr_release_p(fm_deriv)
     443         362 :       CALL dbcsr_release_p(fm_scaled)
     444             : 
     445         362 :       DEALLOCATE (occupations)
     446         362 :       DEALLOCATE (energy_scaling)
     447         362 :       DEALLOCATE (scaling)
     448             : 
     449         362 :       CALL dbcsr_release_p(dbcsr_deriv)
     450             : 
     451         362 :       CALL timestop(handle)
     452             : 
     453      104783 :    END SUBROUTINE low_spin_roks
     454             : ! **************************************************************************************************
     455             : !> \brief do sic calculations on explicit orbitals
     456             : !> \param energy ...
     457             : !> \param qs_env ...
     458             : !> \param dft_control ...
     459             : !> \param poisson_env ...
     460             : !> \param just_energy ...
     461             : !> \param calculate_forces ...
     462             : !> \param auxbas_pw_pool ...
     463             : ! **************************************************************************************************
     464      102973 :    SUBROUTINE sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
     465             :                                     calculate_forces, auxbas_pw_pool)
     466             : 
     467             :       TYPE(qs_energy_type), POINTER                      :: energy
     468             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     469             :       TYPE(dft_control_type), POINTER                    :: dft_control
     470             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     471             :       LOGICAL, INTENT(IN)                                :: just_energy, calculate_forces
     472             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     473             : 
     474             :       CHARACTER(*), PARAMETER :: routineN = 'sic_explicit_orbitals'
     475             : 
     476             :       INTEGER                                            :: handle, i, Iorb, k_alpha, k_beta, Norb
     477      102973 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: sic_orbital_list
     478             :       LOGICAL                                            :: compute_virial, uniform_occupation
     479             :       REAL(KIND=dp)                                      :: ener, exc
     480             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_xc_tmp
     481             :       TYPE(cell_type), POINTER                           :: cell
     482             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     483             :       TYPE(cp_fm_type)                                   :: matrix_hv, matrix_v
     484      102973 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: mo_derivs_local
     485             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     486             :       TYPE(dbcsr_p_type)                                 :: orb_density_matrix_p, orb_h_p
     487      102973 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mo_derivs, rho_ao, tmp_dbcsr
     488             :       TYPE(dbcsr_type), POINTER                          :: orb_density_matrix, orb_h
     489      102973 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     490             :       TYPE(pw_c1d_gs_type)                               :: work_v_gspace
     491      102973 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     492             :       TYPE(pw_c1d_gs_type), TARGET                       :: orb_rho_g, tmp_g
     493             :       TYPE(pw_env_type), POINTER                         :: pw_env
     494             :       TYPE(pw_pool_type), POINTER                        :: xc_pw_pool
     495             :       TYPE(pw_r3d_rs_type)                               :: work_v_rspace
     496      102973 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r, tau, vxc, vxc_tau
     497             :       TYPE(pw_r3d_rs_type), TARGET                       :: orb_rho_r, tmp_r
     498             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     499             :       TYPE(qs_rho_type), POINTER                         :: rho
     500             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     501             :       TYPE(virial_type), POINTER                         :: virial
     502             : 
     503      102973 :       IF (dft_control%sic_method_id .NE. sic_eo) RETURN
     504             : 
     505          40 :       CALL timeset(routineN, handle)
     506             : 
     507          40 :       NULLIFY (tau, vxc_tau, mo_derivs, ks_env, rho_ao)
     508             : 
     509             :       ! generate the lists of orbitals that need sic treatment
     510             :       CALL get_qs_env(qs_env, &
     511             :                       ks_env=ks_env, &
     512             :                       mo_derivs=mo_derivs, &
     513             :                       mos=mo_array, &
     514             :                       rho=rho, &
     515             :                       pw_env=pw_env, &
     516             :                       input=input, &
     517             :                       cell=cell, &
     518          40 :                       virial=virial)
     519             : 
     520          40 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     521             : 
     522          40 :       compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
     523          40 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     524             : 
     525         120 :       DO i = 1, SIZE(mo_array) !fm->dbcsr
     526         120 :          IF (mo_array(i)%use_mo_coeff_b) THEN !fm->dbcsr
     527             :             CALL copy_dbcsr_to_fm(mo_array(i)%mo_coeff_b, &
     528          80 :                                   mo_array(i)%mo_coeff) !fm->dbcsr
     529             :          END IF !fm->dbcsr
     530             :       END DO !fm->dbcsr
     531             : 
     532          40 :       CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
     533             : 
     534             :       ! we have two spins
     535          40 :       CPASSERT(SIZE(mo_array, 1) == 2)
     536             :       ! we want uniform occupations
     537          40 :       CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
     538          40 :       CPASSERT(uniform_occupation)
     539          40 :       CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff, uniform_occupation=uniform_occupation)
     540          40 :       CPASSERT(uniform_occupation)
     541             : 
     542          40 :       NULLIFY (tmp_dbcsr)
     543          40 :       CALL dbcsr_allocate_matrix_set(tmp_dbcsr, SIZE(mo_derivs, 1))
     544         100 :       DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
     545             :          !
     546          60 :          NULLIFY (tmp_dbcsr(i)%matrix)
     547          60 :          CALL dbcsr_init_p(tmp_dbcsr(i)%matrix)
     548          60 :          CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix)
     549         100 :          CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp)
     550             :       END DO !fm->dbcsr
     551             : 
     552          40 :       k_alpha = 0; k_beta = 0
     553          60 :       SELECT CASE (dft_control%sic_list_id)
     554             :       CASE (sic_list_all)
     555             : 
     556          20 :          CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
     557          20 :          CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
     558             : 
     559          20 :          IF (SIZE(mo_array, 1) > 1) THEN
     560          20 :             CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
     561          20 :             CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
     562             :          END IF
     563             : 
     564          20 :          Norb = k_alpha + k_beta
     565          60 :          ALLOCATE (sic_orbital_list(3, Norb))
     566             : 
     567          80 :          iorb = 0
     568          80 :          DO i = 1, k_alpha
     569          60 :             iorb = iorb + 1
     570          60 :             sic_orbital_list(1, iorb) = 1
     571          60 :             sic_orbital_list(2, iorb) = i
     572          80 :             sic_orbital_list(3, iorb) = 1
     573             :          END DO
     574          60 :          DO i = 1, k_beta
     575          20 :             iorb = iorb + 1
     576          20 :             sic_orbital_list(1, iorb) = 2
     577          20 :             sic_orbital_list(2, iorb) = i
     578          40 :             IF (SIZE(mo_derivs, 1) == 1) THEN
     579           0 :                sic_orbital_list(3, iorb) = 1
     580             :             ELSE
     581          20 :                sic_orbital_list(3, iorb) = 2
     582             :             END IF
     583             :          END DO
     584             : 
     585             :       CASE (sic_list_unpaired)
     586             :          ! we have two spins
     587          20 :          CPASSERT(SIZE(mo_array, 1) == 2)
     588             :          ! we have them restricted
     589          20 :          CPASSERT(SIZE(mo_derivs, 1) == 1)
     590          20 :          CPASSERT(dft_control%restricted)
     591             : 
     592          20 :          CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
     593          20 :          CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
     594             : 
     595          20 :          CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
     596          20 :          CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
     597             : 
     598          20 :          Norb = k_alpha - k_beta
     599          60 :          ALLOCATE (sic_orbital_list(3, Norb))
     600             : 
     601          20 :          iorb = 0
     602         100 :          DO i = k_beta + 1, k_alpha
     603          40 :             iorb = iorb + 1
     604          40 :             sic_orbital_list(1, iorb) = 1
     605          40 :             sic_orbital_list(2, iorb) = i
     606             :             ! we are guaranteed to be restricted
     607          60 :             sic_orbital_list(3, iorb) = 1
     608             :          END DO
     609             : 
     610             :       CASE DEFAULT
     611          40 :          CPABORT("")
     612             :       END SELECT
     613             : 
     614             :       ! data needed for each of the orbs
     615          40 :       CALL auxbas_pw_pool%create_pw(orb_rho_r)
     616          40 :       CALL auxbas_pw_pool%create_pw(tmp_r)
     617          40 :       CALL auxbas_pw_pool%create_pw(orb_rho_g)
     618          40 :       CALL auxbas_pw_pool%create_pw(tmp_g)
     619          40 :       CALL auxbas_pw_pool%create_pw(work_v_gspace)
     620          40 :       CALL auxbas_pw_pool%create_pw(work_v_rspace)
     621             : 
     622          40 :       ALLOCATE (orb_density_matrix)
     623             :       CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, &
     624          40 :                       name="orb_density_matrix")
     625          40 :       CALL dbcsr_set(orb_density_matrix, 0.0_dp)
     626          40 :       orb_density_matrix_p%matrix => orb_density_matrix
     627             : 
     628          40 :       ALLOCATE (orb_h)
     629             :       CALL dbcsr_copy(orb_h, rho_ao(1)%matrix, &
     630          40 :                       name="orb_density_matrix")
     631          40 :       CALL dbcsr_set(orb_h, 0.0_dp)
     632          40 :       orb_h_p%matrix => orb_h
     633             : 
     634          40 :       CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
     635             : 
     636             :       CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=1, &
     637          40 :                                template_fmstruct=mo_coeff%matrix_struct)
     638          40 :       CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v")
     639          40 :       CALL cp_fm_create(matrix_hv, fm_struct_tmp, name="matrix_hv")
     640          40 :       CALL cp_fm_struct_release(fm_struct_tmp)
     641             : 
     642         200 :       ALLOCATE (mo_derivs_local(SIZE(mo_array, 1)))
     643         120 :       DO I = 1, SIZE(mo_array, 1)
     644          80 :          CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff)
     645         120 :          CALL cp_fm_create(mo_derivs_local(I), mo_coeff%matrix_struct)
     646             :       END DO
     647             : 
     648         120 :       ALLOCATE (rho_r(2))
     649          40 :       rho_r(1) = orb_rho_r
     650          40 :       rho_r(2) = tmp_r
     651          40 :       CALL pw_zero(tmp_r)
     652             : 
     653         120 :       ALLOCATE (rho_g(2))
     654          40 :       rho_g(1) = orb_rho_g
     655          40 :       rho_g(2) = tmp_g
     656          40 :       CALL pw_zero(tmp_g)
     657             : 
     658          40 :       NULLIFY (vxc)
     659             :       ! now apply to SIC correction to each selected orbital
     660         160 :       DO iorb = 1, Norb
     661             :          ! extract the proper orbital from the mo_coeff
     662         120 :          CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb)), mo_coeff=mo_coeff)
     663         120 :          CALL cp_fm_to_fm(mo_coeff, matrix_v, 1, sic_orbital_list(2, iorb), 1)
     664             : 
     665             :          ! construct the density matrix and the corresponding density
     666         120 :          CALL dbcsr_set(orb_density_matrix, 0.0_dp)
     667             :          CALL cp_dbcsr_plus_fm_fm_t(orb_density_matrix, matrix_v=matrix_v, ncol=1, &
     668         120 :                                     alpha=1.0_dp)
     669             : 
     670             :          CALL calculate_rho_elec(matrix_p=orb_density_matrix, &
     671             :                                  rho=orb_rho_r, rho_gspace=orb_rho_g, &
     672         120 :                                  ks_env=ks_env)
     673             : 
     674             :          ! compute the energy functional for this orbital and its derivative
     675             : 
     676         120 :          CALL pw_poisson_solve(poisson_env, orb_rho_g, ener, work_v_gspace)
     677             :          ! no PBC correction is done here, see "calc_v_sic_rspace" for SIC methods
     678             :          ! with PBC aware corrections
     679         120 :          energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener
     680         120 :          IF (.NOT. just_energy) THEN
     681          72 :             CALL pw_transfer(work_v_gspace, work_v_rspace)
     682          72 :             CALL pw_scale(work_v_rspace, -dft_control%sic_scaling_a*work_v_rspace%pw_grid%dvol)
     683          72 :             CALL dbcsr_set(orb_h, 0.0_dp)
     684             :          END IF
     685             : 
     686         120 :          IF (just_energy) THEN
     687             :             exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
     688          48 :                               pw_pool=xc_pw_pool)
     689             :          ELSE
     690          72 :             CPASSERT(.NOT. compute_virial)
     691             :             CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, &
     692             :                                    rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
     693          72 :                                    pw_pool=xc_pw_pool, compute_virial=compute_virial, virial_xc=virial_xc_tmp)
     694             :             ! add to the existing work_v_rspace
     695          72 :             CALL pw_axpy(vxc(1), work_v_rspace, -dft_control%sic_scaling_b*vxc(1)%pw_grid%dvol)
     696             :          END IF
     697         120 :          energy%exc = energy%exc - dft_control%sic_scaling_b*exc
     698             : 
     699         280 :          IF (.NOT. just_energy) THEN
     700             :             ! note, orb_h (which is being pointed to with orb_h_p) is zeroed above
     701             :             CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=orb_density_matrix_p, hmat=orb_h_p, &
     702          72 :                                     qs_env=qs_env, calculate_forces=calculate_forces)
     703             : 
     704             :             ! add this to the mo_derivs
     705          72 :             CALL cp_dbcsr_sm_fm_multiply(orb_h, matrix_v, matrix_hv, 1)
     706             :             ! silly trick, copy to an array of the right size and add to mo_derivs
     707          72 :             CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb)), 0.0_dp)
     708          72 :             CALL cp_fm_to_fm(matrix_hv, mo_derivs_local(sic_orbital_list(3, iorb)), 1, 1, sic_orbital_list(2, iorb))
     709             :             CALL copy_fm_to_dbcsr(mo_derivs_local(sic_orbital_list(3, iorb)), &
     710          72 :                                   tmp_dbcsr(sic_orbital_list(3, iorb))%matrix)
     711             :             CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, &
     712          72 :                            tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp)
     713             :             !
     714             :             ! need to deallocate vxc
     715          72 :             CALL xc_pw_pool%give_back_pw(vxc(1))
     716          72 :             CALL xc_pw_pool%give_back_pw(vxc(2))
     717          72 :             DEALLOCATE (vxc)
     718             : 
     719             :          END IF
     720             : 
     721             :       END DO
     722             : 
     723          40 :       CALL auxbas_pw_pool%give_back_pw(orb_rho_r)
     724          40 :       CALL auxbas_pw_pool%give_back_pw(tmp_r)
     725          40 :       CALL auxbas_pw_pool%give_back_pw(orb_rho_g)
     726          40 :       CALL auxbas_pw_pool%give_back_pw(tmp_g)
     727          40 :       CALL auxbas_pw_pool%give_back_pw(work_v_gspace)
     728          40 :       CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
     729             : 
     730          40 :       CALL dbcsr_deallocate_matrix(orb_density_matrix)
     731          40 :       CALL dbcsr_deallocate_matrix(orb_h)
     732          40 :       CALL cp_fm_release(matrix_v)
     733          40 :       CALL cp_fm_release(matrix_hv)
     734          40 :       CALL cp_fm_release(mo_derivs_local)
     735          40 :       DEALLOCATE (rho_r)
     736          40 :       DEALLOCATE (rho_g)
     737             : 
     738          40 :       CALL dbcsr_deallocate_matrix_set(tmp_dbcsr) !fm->dbcsr
     739             : 
     740          40 :       CALL timestop(handle)
     741             : 
     742      103133 :    END SUBROUTINE sic_explicit_orbitals
     743             : 
     744             : ! **************************************************************************************************
     745             : !> \brief do sic calculations on the spin density
     746             : !> \param v_sic_rspace ...
     747             : !> \param energy ...
     748             : !> \param qs_env ...
     749             : !> \param dft_control ...
     750             : !> \param rho ...
     751             : !> \param poisson_env ...
     752             : !> \param just_energy ...
     753             : !> \param calculate_forces ...
     754             : !> \param auxbas_pw_pool ...
     755             : ! **************************************************************************************************
     756      102973 :    SUBROUTINE calc_v_sic_rspace(v_sic_rspace, energy, &
     757             :                                 qs_env, dft_control, rho, poisson_env, just_energy, &
     758             :                                 calculate_forces, auxbas_pw_pool)
     759             : 
     760             :       TYPE(pw_r3d_rs_type), POINTER                      :: v_sic_rspace
     761             :       TYPE(qs_energy_type), POINTER                      :: energy
     762             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     763             :       TYPE(dft_control_type), POINTER                    :: dft_control
     764             :       TYPE(qs_rho_type), POINTER                         :: rho
     765             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     766             :       LOGICAL, INTENT(IN)                                :: just_energy, calculate_forces
     767             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     768             : 
     769             :       INTEGER                                            :: i, nelec, nelec_a, nelec_b, nforce
     770             :       REAL(kind=dp)                                      :: ener, full_scaling, scaling
     771      102973 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: store_forces
     772      102973 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     773             :       TYPE(pw_c1d_gs_type)                               :: work_rho, work_v
     774      102973 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     775      102973 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     776             : 
     777      102973 :       NULLIFY (mo_array, rho_g)
     778             : 
     779      102973 :       IF (dft_control%sic_method_id == sic_none) RETURN
     780         336 :       IF (dft_control%sic_method_id == sic_eo) RETURN
     781             : 
     782         296 :       IF (dft_control%qs_control%gapw) &
     783           0 :          CPABORT("sic and GAPW not yet compatible")
     784             : 
     785             :       ! OK, right now we like two spins to do sic, could be relaxed for AD
     786         296 :       CPASSERT(dft_control%nspins == 2)
     787             : 
     788         296 :       CALL auxbas_pw_pool%create_pw(work_rho)
     789         296 :       CALL auxbas_pw_pool%create_pw(work_v)
     790             : 
     791         296 :       CALL qs_rho_get(rho, rho_g=rho_g)
     792             : 
     793             :       ! Hartree sic corrections
     794         566 :       SELECT CASE (dft_control%sic_method_id)
     795             :       CASE (sic_mauri_us, sic_mauri_spz)
     796         270 :          CALL pw_copy(rho_g(1), work_rho)
     797         270 :          CALL pw_axpy(rho_g(2), work_rho, alpha=-1._dp)
     798         296 :          CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
     799             :       CASE (sic_ad)
     800             :          ! find out how many elecs we have
     801          26 :          CALL get_qs_env(qs_env, mos=mo_array)
     802          26 :          CALL get_mo_set(mo_set=mo_array(1), nelectron=nelec_a)
     803          26 :          CALL get_mo_set(mo_set=mo_array(2), nelectron=nelec_b)
     804          26 :          nelec = nelec_a + nelec_b
     805          26 :          CALL pw_copy(rho_g(1), work_rho)
     806          26 :          CALL pw_axpy(rho_g(2), work_rho)
     807          26 :          scaling = 1.0_dp/REAL(nelec, KIND=dp)
     808          26 :          CALL pw_scale(work_rho, scaling)
     809          26 :          CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
     810             :       CASE DEFAULT
     811         618 :          CPABORT("Unknown sic method id")
     812             :       END SELECT
     813             : 
     814             :       ! Correct for  DDAP charges (if any)
     815             :       ! storing whatever force might be there from previous decoupling
     816         296 :       IF (calculate_forces) THEN
     817          48 :          CALL get_qs_env(qs_env=qs_env, force=force)
     818          48 :          nforce = 0
     819         112 :          DO i = 1, SIZE(force)
     820         112 :             nforce = nforce + SIZE(force(i)%ch_pulay, 2)
     821             :          END DO
     822         144 :          ALLOCATE (store_forces(3, nforce))
     823         112 :          nforce = 0
     824         112 :          DO i = 1, SIZE(force)
     825         784 :             store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :)
     826         784 :             force(i)%ch_pulay(:, :) = 0.0_dp
     827         112 :             nforce = nforce + SIZE(force(i)%ch_pulay, 2)
     828             :          END DO
     829             :       END IF
     830             : 
     831             :       CALL cp_ddapc_apply_CD(qs_env, &
     832             :                              work_rho, &
     833             :                              ener, &
     834             :                              v_hartree_gspace=work_v, &
     835             :                              calculate_forces=calculate_forces, &
     836         296 :                              Itype_of_density="SPIN")
     837             : 
     838         566 :       SELECT CASE (dft_control%sic_method_id)
     839             :       CASE (sic_mauri_us, sic_mauri_spz)
     840         270 :          full_scaling = -dft_control%sic_scaling_a
     841             :       CASE (sic_ad)
     842          26 :          full_scaling = -dft_control%sic_scaling_a*nelec
     843             :       CASE DEFAULT
     844         296 :          CPABORT("Unknown sic method id")
     845             :       END SELECT
     846         296 :       energy%hartree = energy%hartree + full_scaling*ener
     847             : 
     848             :       ! add scaled forces, restoring the old
     849         296 :       IF (calculate_forces) THEN
     850          48 :          nforce = 0
     851         112 :          DO i = 1, SIZE(force)
     852             :             force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + &
     853         784 :                                       store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2))
     854         112 :             nforce = nforce + SIZE(force(i)%ch_pulay, 2)
     855             :          END DO
     856             :       END IF
     857             : 
     858         296 :       IF (.NOT. just_energy) THEN
     859         200 :          ALLOCATE (v_sic_rspace)
     860         200 :          CALL auxbas_pw_pool%create_pw(v_sic_rspace)
     861         200 :          CALL pw_transfer(work_v, v_sic_rspace)
     862             :          ! also take into account the scaling (in addition to the volume element)
     863             :          CALL pw_scale(v_sic_rspace, &
     864         200 :                        dft_control%sic_scaling_a*v_sic_rspace%pw_grid%dvol)
     865             :       END IF
     866             : 
     867         296 :       CALL auxbas_pw_pool%give_back_pw(work_rho)
     868         296 :       CALL auxbas_pw_pool%give_back_pw(work_v)
     869             : 
     870      103269 :    END SUBROUTINE calc_v_sic_rspace
     871             : 
     872             : ! **************************************************************************************************
     873             : !> \brief ...
     874             : !> \param qs_env ...
     875             : !> \param rho ...
     876             : ! **************************************************************************************************
     877      205902 :    SUBROUTINE print_densities(qs_env, rho)
     878             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     879             :       TYPE(qs_rho_type), POINTER                         :: rho
     880             : 
     881             :       INTEGER                                            :: img, ispin, n_electrons, output_unit
     882             :       REAL(dp)                                           :: tot1_h, tot1_s, tot_rho_r, trace, &
     883             :                                                             trace_tmp
     884      102951 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_r_arr
     885             :       TYPE(cell_type), POINTER                           :: cell
     886             :       TYPE(cp_logger_type), POINTER                      :: logger
     887      102951 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s, rho_ao
     888             :       TYPE(dft_control_type), POINTER                    :: dft_control
     889             :       TYPE(qs_charges_type), POINTER                     :: qs_charges
     890      102951 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     891             :       TYPE(section_vals_type), POINTER                   :: input, scf_section
     892             : 
     893      102951 :       NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
     894      102951 :                dft_control, tot_rho_r_arr, rho_ao)
     895             : 
     896      205902 :       logger => cp_get_default_logger()
     897             : 
     898             :       CALL get_qs_env(qs_env, &
     899             :                       qs_kind_set=qs_kind_set, &
     900             :                       cell=cell, qs_charges=qs_charges, &
     901             :                       input=input, &
     902             :                       matrix_s_kp=matrix_s, &
     903      102951 :                       dft_control=dft_control)
     904             : 
     905      102951 :       CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
     906             : 
     907      102951 :       scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
     908             :       output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
     909      102951 :                                          extension=".scfLog")
     910             : 
     911      102951 :       CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
     912      102951 :       n_electrons = n_electrons - dft_control%charge
     913      102951 :       tot_rho_r = accurate_sum(tot_rho_r_arr)
     914             : 
     915      102951 :       trace = 0
     916      102951 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN
     917        4146 :          DO ispin = 1, dft_control%nspins
     918        6752 :             DO img = 1, dft_control%nimages
     919        2606 :                CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
     920        5078 :                trace = trace + trace_tmp
     921             :             END DO
     922             :          END DO
     923             :       END IF
     924             : 
     925      102951 :       IF (output_unit > 0) THEN
     926         837 :          WRITE (UNIT=output_unit, FMT="(/,T3,A,T41,F20.10)") "Trace(PS):", trace
     927             :          WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
     928         837 :             "Electronic density on regular grids: ", &
     929         837 :             tot_rho_r, &
     930             :             tot_rho_r + &
     931         837 :             REAL(n_electrons, dp), &
     932         837 :             "Core density on regular grids:", &
     933         837 :             qs_charges%total_rho_core_rspace, &
     934        1674 :             qs_charges%total_rho_core_rspace - REAL(n_electrons + dft_control%charge, dp)
     935             :       END IF
     936      102951 :       IF (dft_control%qs_control%gapw) THEN
     937       13456 :          tot1_h = qs_charges%total_rho1_hard(1)
     938       13456 :          tot1_s = qs_charges%total_rho1_soft(1)
     939       16744 :          DO ispin = 2, dft_control%nspins
     940        3288 :             tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
     941       16744 :             tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
     942             :          END DO
     943       13456 :          IF (output_unit > 0) THEN
     944             :             WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
     945         399 :                "Hard and soft densities (Lebedev):", &
     946         798 :                tot1_h, tot1_s
     947             :             WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
     948         399 :                "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
     949         399 :                tot_rho_r + tot1_h - tot1_s, &
     950         399 :                "Total charge density (r-space):      ", &
     951             :                tot_rho_r + tot1_h - tot1_s &
     952         399 :                + qs_charges%total_rho_core_rspace, &
     953         399 :                "Total Rho_soft + Rho0_soft (g-space):", &
     954         798 :                qs_charges%total_rho_gspace
     955             :          END IF
     956             :          qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
     957       13456 :                                  qs_charges%total_rho_core_rspace
     958       89495 :       ELSE IF (dft_control%qs_control%gapw_xc) THEN
     959        2608 :          tot1_h = qs_charges%total_rho1_hard(1)
     960        2608 :          tot1_s = qs_charges%total_rho1_soft(1)
     961        2970 :          DO ispin = 2, dft_control%nspins
     962         362 :             tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
     963        2970 :             tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
     964             :          END DO
     965        2608 :          IF (output_unit > 0) THEN
     966             :             WRITE (UNIT=output_unit, FMT="(/,(T3,A,T41,2F20.10))") &
     967           0 :                "Hard and soft densities (Lebedev):", &
     968           0 :                tot1_h, tot1_s
     969             :             WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
     970           0 :                "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
     971           0 :                accurate_sum(tot_rho_r_arr) + tot1_h - tot1_s
     972             :          END IF
     973             :          qs_charges%background = tot_rho_r + &
     974        2608 :                                  qs_charges%total_rho_core_rspace
     975             :       ELSE
     976       86887 :          IF (output_unit > 0) THEN
     977             :             WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
     978         438 :                "Total charge density on r-space grids:     ", &
     979             :                tot_rho_r + &
     980         438 :                qs_charges%total_rho_core_rspace, &
     981         438 :                "Total charge density g-space grids:     ", &
     982         876 :                qs_charges%total_rho_gspace
     983             :          END IF
     984             :          qs_charges%background = tot_rho_r + &
     985       86887 :                                  qs_charges%total_rho_core_rspace
     986             :       END IF
     987      102951 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="()")
     988      102951 :       qs_charges%background = qs_charges%background/cell%deth
     989             : 
     990             :       CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
     991      102951 :                                         "PRINT%TOTAL_DENSITIES")
     992             : 
     993      102951 :    END SUBROUTINE print_densities
     994             : 
     995             : ! **************************************************************************************************
     996             : !> \brief Print detailed energies
     997             : !>
     998             : !> \param qs_env ...
     999             : !> \param dft_control ...
    1000             : !> \param input ...
    1001             : !> \param energy ...
    1002             : !> \param mulliken_order_p ...
    1003             : !> \par History
    1004             : !>    refactoring 04.03.2011 [MI]
    1005             : !> \author
    1006             : ! **************************************************************************************************
    1007      102951 :    SUBROUTINE print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
    1008             : 
    1009             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1010             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1011             :       TYPE(section_vals_type), POINTER                   :: input
    1012             :       TYPE(qs_energy_type), POINTER                      :: energy
    1013             :       REAL(KIND=dp), INTENT(IN)                          :: mulliken_order_p
    1014             : 
    1015             :       INTEGER                                            :: bc, n, output_unit, psolver
    1016             :       REAL(KIND=dp)                                      :: ddapc_order_p, implicit_ps_ehartree, &
    1017             :                                                             s2_order_p
    1018             :       TYPE(cp_logger_type), POINTER                      :: logger
    1019             :       TYPE(pw_env_type), POINTER                         :: pw_env
    1020             : 
    1021      102951 :       logger => cp_get_default_logger()
    1022             : 
    1023      102951 :       NULLIFY (pw_env)
    1024      102951 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    1025      102951 :       psolver = pw_env%poisson_env%parameters%solver
    1026             : 
    1027             :       output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", &
    1028      102951 :                                          extension=".scfLog")
    1029      102951 :       IF (output_unit > 0) THEN
    1030         490 :          IF (dft_control%do_admm) THEN
    1031             :             WRITE (UNIT=output_unit, FMT="((T3,A,T60,F20.10))") &
    1032           0 :                "Wfn fit exchange-correlation energy:            ", energy%exc_aux_fit
    1033           0 :             IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
    1034             :                WRITE (UNIT=output_unit, FMT="((T3,A,T60,F20.10))") &
    1035           0 :                   "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit
    1036             :             END IF
    1037             :          END IF
    1038         490 :          IF (dft_control%do_admm) THEN
    1039           0 :             IF (psolver .EQ. pw_poisson_implicit) THEN
    1040           0 :                implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
    1041           0 :                bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
    1042           0 :                SELECT CASE (bc)
    1043             :                CASE (MIXED_PERIODIC_BC, MIXED_BC)
    1044             :                   WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1045           0 :                      "Core Hamiltonian energy:                       ", energy%core, &
    1046           0 :                      "Hartree energy:                                ", implicit_ps_ehartree, &
    1047           0 :                      "Electric enthalpy:                             ", energy%hartree, &
    1048           0 :                      "Exchange-correlation energy:                   ", energy%exc + energy%exc_aux_fit
    1049             :                CASE (PERIODIC_BC, NEUMANN_BC)
    1050             :                   WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1051           0 :                      "Core Hamiltonian energy:                       ", energy%core, &
    1052           0 :                      "Hartree energy:                                ", energy%hartree, &
    1053           0 :                      "Exchange-correlation energy:                   ", energy%exc + energy%exc_aux_fit
    1054             :                END SELECT
    1055             :             ELSE
    1056             :                WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1057           0 :                   "Core Hamiltonian energy:                       ", energy%core, &
    1058           0 :                   "Hartree energy:                                ", energy%hartree, &
    1059           0 :                   "Exchange-correlation energy:                   ", energy%exc + energy%exc_aux_fit
    1060             :             END IF
    1061             :          ELSE
    1062             : !ZMP to print some variables at each step
    1063         490 :             IF (dft_control%apply_external_density) THEN
    1064             :                WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1065           0 :                   "DOING ZMP CALCULATION FROM EXTERNAL DENSITY    "
    1066             :                WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1067           0 :                   "Core Hamiltonian energy:                       ", energy%core, &
    1068           0 :                   "Hartree energy:                                ", energy%hartree
    1069         490 :             ELSE IF (dft_control%apply_external_vxc) THEN
    1070             :                WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1071           0 :                   "DOING ZMP READING EXTERNAL VXC                 "
    1072             :                WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1073           0 :                   "Core Hamiltonian energy:                       ", energy%core, &
    1074           0 :                   "Hartree energy:                                ", energy%hartree
    1075             :             ELSE
    1076         490 :                IF (psolver .EQ. pw_poisson_implicit) THEN
    1077           0 :                   implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
    1078           0 :                   bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
    1079           0 :                   SELECT CASE (bc)
    1080             :                   CASE (MIXED_PERIODIC_BC, MIXED_BC)
    1081             :                      WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1082           0 :                         "Core Hamiltonian energy:                       ", energy%core, &
    1083           0 :                         "Hartree energy:                                ", implicit_ps_ehartree, &
    1084           0 :                         "Electric enthalpy:                             ", energy%hartree, &
    1085           0 :                         "Exchange-correlation energy:                   ", energy%exc
    1086             :                   CASE (PERIODIC_BC, NEUMANN_BC)
    1087             :                      WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1088           0 :                         "Core Hamiltonian energy:                       ", energy%core, &
    1089           0 :                         "Hartree energy:                                ", energy%hartree, &
    1090           0 :                         "Exchange-correlation energy:                   ", energy%exc
    1091             :                   END SELECT
    1092             :                ELSE
    1093             :                   WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1094         490 :                      "Core Hamiltonian energy:                       ", energy%core, &
    1095         490 :                      "Hartree energy:                                ", energy%hartree, &
    1096         980 :                      "Exchange-correlation energy:                   ", energy%exc
    1097             :                END IF
    1098             :             END IF
    1099             :          END IF
    1100             : 
    1101         490 :          IF (dft_control%apply_external_density) THEN
    1102             :             WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1103           0 :                "Integral of the (density * v_xc):              ", energy%exc
    1104             :          END IF
    1105             : 
    1106         490 :          IF (energy%e_hartree /= 0.0_dp) &
    1107             :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1108         458 :             "Coulomb (electron-electron) energy:            ", energy%e_hartree
    1109         490 :          IF (energy%dispersion /= 0.0_dp) &
    1110             :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1111           0 :             "Dispersion energy:                             ", energy%dispersion
    1112         490 :          IF (energy%efield /= 0.0_dp) &
    1113             :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1114           0 :             "Electric field interaction energy:             ", energy%efield
    1115         490 :          IF (energy%gcp /= 0.0_dp) &
    1116             :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1117           0 :             "gCP energy:                                    ", energy%gcp
    1118         490 :          IF (dft_control%qs_control%gapw) THEN
    1119             :             WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1120          32 :                "GAPW| Exc from hard and soft atomic rho1:      ", energy%exc1 + energy%exc1_aux_fit, &
    1121          64 :                "GAPW| local Eh = 1 center integrals:           ", energy%hartree_1c
    1122             :          END IF
    1123         490 :          IF (dft_control%qs_control%gapw_xc) THEN
    1124             :             WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
    1125           0 :                "GAPW| Exc from hard and soft atomic rho1:      ", energy%exc1 + energy%exc1_aux_fit
    1126             :          END IF
    1127         490 :          IF (dft_control%dft_plus_u) THEN
    1128             :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1129           0 :                "DFT+U energy:", energy%dft_plus_u
    1130             :          END IF
    1131         490 :          IF (qs_env%qmmm) THEN
    1132             :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1133           0 :                "QM/MM Electrostatic energy:                    ", energy%qmmm_el
    1134           0 :             IF (qs_env%qmmm_env_qm%image_charge) THEN
    1135             :                WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1136           0 :                   "QM/MM image charge energy:                ", energy%image_charge
    1137             :             END IF
    1138             :          END IF
    1139         490 :          IF (dft_control%qs_control%mulliken_restraint) THEN
    1140             :             WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") &
    1141           0 :                "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken
    1142             :          END IF
    1143         490 :          IF (dft_control%qs_control%ddapc_restraint) THEN
    1144          40 :             DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
    1145             :                ddapc_order_p = &
    1146          20 :                   dft_control%qs_control%ddapc_restraint_control(n)%ddapc_order_p
    1147             :                WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") &
    1148          40 :                   "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n)
    1149             :             END DO
    1150             :          END IF
    1151         490 :          IF (dft_control%qs_control%s2_restraint) THEN
    1152           0 :             s2_order_p = dft_control%qs_control%s2_restraint_control%s2_order_p
    1153             :             WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") &
    1154           0 :                "S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint
    1155             :          END IF
    1156             : 
    1157             :       END IF ! output_unit
    1158             :       CALL cp_print_key_finished_output(output_unit, logger, input, &
    1159      102951 :                                         "DFT%SCF%PRINT%DETAILED_ENERGY")
    1160             : 
    1161      102951 :    END SUBROUTINE print_detailed_energy
    1162             : 
    1163             : ! **************************************************************************************************
    1164             : !> \brief compute matrix_vxc, defined via the potential created by qs_vxc_create
    1165             : !>        ignores things like tau functional, gapw, sic, ...
    1166             : !>         so only OK for GGA & GPW right now
    1167             : !> \param qs_env ...
    1168             : !> \param v_rspace ...
    1169             : !> \param matrix_vxc ...
    1170             : !> \par History
    1171             : !>    created 23.10.2012 [Joost VandeVondele]
    1172             : !> \author
    1173             : ! **************************************************************************************************
    1174           4 :    SUBROUTINE compute_matrix_vxc(qs_env, v_rspace, matrix_vxc)
    1175             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1176             :       TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN)     :: v_rspace
    1177             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_vxc
    1178             : 
    1179             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_matrix_vxc'
    1180             : 
    1181             :       INTEGER                                            :: handle, ispin
    1182             :       LOGICAL                                            :: gapw
    1183           4 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
    1184             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1185             : 
    1186           4 :       CALL timeset(routineN, handle)
    1187             : 
    1188             :       ! create the matrix using matrix_ks as a template
    1189           4 :       IF (ASSOCIATED(matrix_vxc)) THEN
    1190           0 :          CALL dbcsr_deallocate_matrix_set(matrix_vxc)
    1191             :       END IF
    1192           4 :       CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
    1193          18 :       ALLOCATE (matrix_vxc(SIZE(matrix_ks)))
    1194          10 :       DO ispin = 1, SIZE(matrix_ks)
    1195           6 :          NULLIFY (matrix_vxc(ispin)%matrix)
    1196           6 :          CALL dbcsr_init_p(matrix_vxc(ispin)%matrix)
    1197             :          CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
    1198           6 :                          name="Matrix VXC of spin "//cp_to_string(ispin))
    1199          10 :          CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
    1200             :       END DO
    1201             : 
    1202             :       ! and integrate
    1203           4 :       CALL get_qs_env(qs_env, dft_control=dft_control)
    1204           4 :       gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
    1205          10 :       DO ispin = 1, SIZE(matrix_ks)
    1206             :          CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
    1207             :                                  hmat=matrix_vxc(ispin), &
    1208             :                                  qs_env=qs_env, &
    1209             :                                  calculate_forces=.FALSE., &
    1210           6 :                                  gapw=gapw)
    1211             :          ! scale by the volume element... should really become part of integrate_v_rspace
    1212          10 :          CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
    1213             :       END DO
    1214             : 
    1215           4 :       CALL timestop(handle)
    1216             : 
    1217           4 :    END SUBROUTINE compute_matrix_vxc
    1218             : 
    1219             : ! **************************************************************************************************
    1220             : !> \brief Sum up all potentials defined  on the grid and integrate
    1221             : !>
    1222             : !> \param qs_env ...
    1223             : !> \param ks_matrix ...
    1224             : !> \param rho ...
    1225             : !> \param my_rho ...
    1226             : !> \param vppl_rspace ...
    1227             : !> \param v_rspace_new ...
    1228             : !> \param v_rspace_new_aux_fit ...
    1229             : !> \param v_tau_rspace ...
    1230             : !> \param v_tau_rspace_aux_fit ...
    1231             : !> \param v_sic_rspace ...
    1232             : !> \param v_spin_ddapc_rest_r ...
    1233             : !> \param v_sccs_rspace ...
    1234             : !> \param v_rspace_embed ...
    1235             : !> \param cdft_control ...
    1236             : !> \param calculate_forces ...
    1237             : !> \par History
    1238             : !>      - refactoring 04.03.2011 [MI]
    1239             : !>      - SCCS implementation (16.10.2013,MK)
    1240             : !> \author
    1241             : ! **************************************************************************************************
    1242       93637 :    SUBROUTINE sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, &
    1243             :                                    vppl_rspace, v_rspace_new, &
    1244             :                                    v_rspace_new_aux_fit, v_tau_rspace, &
    1245             :                                    v_tau_rspace_aux_fit, &
    1246             :                                    v_sic_rspace, v_spin_ddapc_rest_r, &
    1247             :                                    v_sccs_rspace, v_rspace_embed, cdft_control, &
    1248             :                                    calculate_forces)
    1249             : 
    1250             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1251             :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_matrix
    1252             :       TYPE(qs_rho_type), POINTER                         :: rho
    1253             :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: my_rho
    1254             :       TYPE(pw_r3d_rs_type), POINTER                      :: vppl_rspace
    1255             :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_rspace_new, v_rspace_new_aux_fit, &
    1256             :                                                             v_tau_rspace, v_tau_rspace_aux_fit
    1257             :       TYPE(pw_r3d_rs_type), POINTER                      :: v_sic_rspace, v_spin_ddapc_rest_r, &
    1258             :                                                             v_sccs_rspace
    1259             :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_rspace_embed
    1260             :       TYPE(cdft_control_type), POINTER                   :: cdft_control
    1261             :       LOGICAL, INTENT(in)                                :: calculate_forces
    1262             : 
    1263             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'sum_up_and_integrate'
    1264             : 
    1265             :       CHARACTER(LEN=default_string_length)               :: basis_type
    1266             :       INTEGER                                            :: handle, igroup, ikind, img, ispin, &
    1267             :                                                             nkind, nspins
    1268       93637 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    1269             :       LOGICAL                                            :: do_ppl, gapw, gapw_xc, lrigpw, rigpw
    1270             :       REAL(KIND=dp)                                      :: csign, dvol, fadm
    1271             :       TYPE(admm_type), POINTER                           :: admm_env
    1272       93637 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1273       93637 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ksmat, rho_ao, rho_ao_nokp, smat
    1274       93637 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_aux_fit, &
    1275       93637 :                                                             matrix_ks_aux_fit_dft, rho_ao_aux, &
    1276       93637 :                                                             rho_ao_kp
    1277             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1278             :       TYPE(kpoint_type), POINTER                         :: kpoints
    1279             :       TYPE(lri_density_type), POINTER                    :: lri_density
    1280             :       TYPE(lri_environment_type), POINTER                :: lri_env
    1281       93637 :       TYPE(lri_kind_type), DIMENSION(:), POINTER         :: lri_v_int
    1282             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1283             :       TYPE(pw_env_type), POINTER                         :: pw_env
    1284             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
    1285             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1286             :       TYPE(pw_r3d_rs_type), POINTER                      :: v_rspace, vee
    1287             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    1288             :       TYPE(qs_rho_type), POINTER                         :: rho_aux_fit
    1289             :       TYPE(task_list_type), POINTER                      :: task_list
    1290             : 
    1291       93637 :       CALL timeset(routineN, handle)
    1292             : 
    1293       93637 :       NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
    1294       93637 :                v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
    1295       93637 :                ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
    1296       93637 :                rho_ao_nokp, ks_env, admm_env, task_list)
    1297             : 
    1298             :       CALL get_qs_env(qs_env, &
    1299             :                       dft_control=dft_control, &
    1300             :                       pw_env=pw_env, &
    1301             :                       v_hartree_rspace=v_rspace, &
    1302       93637 :                       vee=vee)
    1303             : 
    1304       93637 :       CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
    1305       93637 :       CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
    1306       93637 :       gapw = dft_control%qs_control%gapw
    1307       93637 :       gapw_xc = dft_control%qs_control%gapw_xc
    1308       93637 :       do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
    1309             : 
    1310       93637 :       rigpw = dft_control%qs_control%rigpw
    1311       93637 :       lrigpw = dft_control%qs_control%lrigpw
    1312       93637 :       IF (lrigpw .OR. rigpw) THEN
    1313             :          CALL get_qs_env(qs_env, &
    1314             :                          lri_env=lri_env, &
    1315             :                          lri_density=lri_density, &
    1316         428 :                          atomic_kind_set=atomic_kind_set)
    1317             :       END IF
    1318             : 
    1319       93637 :       nspins = dft_control%nspins
    1320             : 
    1321             :       ! sum up potentials and integrate
    1322       93637 :       IF (ASSOCIATED(v_rspace_new)) THEN
    1323      185519 :          DO ispin = 1, nspins
    1324      100674 :             IF (gapw_xc) THEN
    1325             :                ! SIC not implemented (or at least not tested)
    1326        2796 :                CPASSERT(dft_control%sic_method_id == sic_none)
    1327             :                !Only the xc potential, because it has to be integrated with the soft basis
    1328        2796 :                CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
    1329             : 
    1330             :                ! add the xc  part due to v_rspace soft
    1331        2796 :                rho_ao => rho_ao_kp(ispin, :)
    1332        2796 :                ksmat => ks_matrix(ispin, :)
    1333             :                CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
    1334             :                                        pmat_kp=rho_ao, hmat_kp=ksmat, &
    1335             :                                        qs_env=qs_env, &
    1336             :                                        calculate_forces=calculate_forces, &
    1337        2796 :                                        gapw=gapw_xc)
    1338             : 
    1339             :                ! Now the Hartree potential to be integrated with the full basis
    1340        2796 :                CALL pw_copy(v_rspace, v_rspace_new(ispin))
    1341             :             ELSE
    1342             :                ! Add v_hartree + v_xc = v_rspace_new
    1343       97878 :                CALL pw_axpy(v_rspace, v_rspace_new(ispin), 1.0_dp, v_rspace_new(ispin)%pw_grid%dvol)
    1344             :             END IF ! gapw_xc
    1345      100674 :             IF (dft_control%qs_control%ddapc_explicit_potential) THEN
    1346         184 :                IF (dft_control%qs_control%ddapc_restraint_is_spin) THEN
    1347         184 :                   IF (ispin == 1) THEN
    1348          92 :                      CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
    1349             :                   ELSE
    1350          92 :                      CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), -1.0_dp)
    1351             :                   END IF
    1352             :                ELSE
    1353           0 :                   CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
    1354             :                END IF
    1355             :             END IF
    1356             :             ! CDFT constraint contribution
    1357      100674 :             IF (dft_control%qs_control%cdft) THEN
    1358       11336 :                DO igroup = 1, SIZE(cdft_control%group)
    1359        6644 :                   SELECT CASE (cdft_control%group(igroup)%constraint_type)
    1360             :                   CASE (cdft_charge_constraint)
    1361          16 :                      csign = 1.0_dp
    1362             :                   CASE (cdft_magnetization_constraint)
    1363          16 :                      IF (ispin == 1) THEN
    1364             :                         csign = 1.0_dp
    1365             :                      ELSE
    1366           8 :                         csign = -1.0_dp
    1367             :                      END IF
    1368             :                   CASE (cdft_alpha_constraint)
    1369        1944 :                      csign = 1.0_dp
    1370        1944 :                      IF (ispin == 2) CYCLE
    1371             :                   CASE (cdft_beta_constraint)
    1372        1944 :                      csign = 1.0_dp
    1373        1944 :                      IF (ispin == 1) CYCLE
    1374             :                   CASE DEFAULT
    1375        6644 :                      CPABORT("Unknown constraint type.")
    1376             :                   END SELECT
    1377             :                   CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_new(ispin), &
    1378       11336 :                                csign*cdft_control%strength(igroup))
    1379             :                END DO
    1380             :             END IF
    1381             :             ! functional derivative of the Hartree energy wrt the density in the presence of dielectric
    1382             :             ! (vhartree + v_eps); v_eps is nonzero only if the dielectric constant is defind as a function
    1383             :             ! of the charge density
    1384      100674 :             IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
    1385         440 :                dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
    1386         440 :                CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_new(ispin), dvol)
    1387             :             END IF
    1388             :             ! Add SCCS contribution
    1389      100674 :             IF (dft_control%do_sccs) THEN
    1390         132 :                CALL pw_axpy(v_sccs_rspace, v_rspace_new(ispin))
    1391             :             END IF
    1392             :             ! External electrostatic potential
    1393      100674 :             IF (dft_control%apply_external_potential) THEN
    1394             :                CALL qmmm_modify_hartree_pot(v_hartree=v_rspace_new(ispin), &
    1395         364 :                                             v_qmmm=vee, scale=-1.0_dp)
    1396             :             END IF
    1397      100674 :             IF (do_ppl) THEN
    1398          66 :                CPASSERT(.NOT. gapw)
    1399          66 :                CALL pw_axpy(vppl_rspace, v_rspace_new(ispin), vppl_rspace%pw_grid%dvol)
    1400             :             END IF
    1401             :             ! the electrostatic sic contribution
    1402      101034 :             SELECT CASE (dft_control%sic_method_id)
    1403             :             CASE (sic_none)
    1404             :                !
    1405             :             CASE (sic_mauri_us, sic_mauri_spz)
    1406         360 :                IF (ispin == 1) THEN
    1407         180 :                   CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
    1408             :                ELSE
    1409         180 :                   CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), 1.0_dp)
    1410             :                END IF
    1411             :             CASE (sic_ad)
    1412      100674 :                CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
    1413             :             CASE (sic_eo)
    1414             :                ! NOTHING TO BE DONE
    1415             :             END SELECT
    1416             :             ! DFT embedding
    1417      100674 :             IF (dft_control%apply_embed_pot) THEN
    1418         930 :                CALL pw_axpy(v_rspace_embed(ispin), v_rspace_new(ispin), v_rspace_embed(ispin)%pw_grid%dvol)
    1419         930 :                CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
    1420             :             END IF
    1421      100674 :             IF (lrigpw) THEN
    1422         440 :                lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
    1423         440 :                CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
    1424        1322 :                DO ikind = 1, nkind
    1425      287614 :                   lri_v_int(ikind)%v_int = 0.0_dp
    1426             :                END DO
    1427             :                CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
    1428         440 :                                                   lri_v_int, calculate_forces, "LRI_AUX")
    1429        1322 :                DO ikind = 1, nkind
    1430      573906 :                   CALL para_env%sum(lri_v_int(ikind)%v_int)
    1431             :                END DO
    1432         440 :                IF (lri_env%exact_1c_terms) THEN
    1433          36 :                   rho_ao => my_rho(ispin, :)
    1434          36 :                   ksmat => ks_matrix(ispin, :)
    1435             :                   CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), ksmat(1)%matrix, &
    1436             :                                                    rho_ao(1)%matrix, qs_env, &
    1437          36 :                                                    calculate_forces, "ORB")
    1438             :                END IF
    1439         440 :                IF (lri_env%ppl_ri) THEN
    1440           8 :                   CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
    1441             :                END IF
    1442      100234 :             ELSEIF (rigpw) THEN
    1443           0 :                lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
    1444           0 :                CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
    1445           0 :                DO ikind = 1, nkind
    1446           0 :                   lri_v_int(ikind)%v_int = 0.0_dp
    1447             :                END DO
    1448             :                CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
    1449           0 :                                                   lri_v_int, calculate_forces, "RI_HXC")
    1450           0 :                DO ikind = 1, nkind
    1451           0 :                   CALL para_env%sum(lri_v_int(ikind)%v_int)
    1452             :                END DO
    1453             :             ELSE
    1454      100234 :                rho_ao => my_rho(ispin, :)
    1455      100234 :                ksmat => ks_matrix(ispin, :)
    1456             :                CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
    1457             :                                        pmat_kp=rho_ao, hmat_kp=ksmat, &
    1458             :                                        qs_env=qs_env, &
    1459             :                                        calculate_forces=calculate_forces, &
    1460      100234 :                                        gapw=gapw)
    1461             :             END IF
    1462      185519 :             CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
    1463             :          END DO ! ispin
    1464             : 
    1465       85045 :          SELECT CASE (dft_control%sic_method_id)
    1466             :          CASE (sic_none)
    1467             :          CASE (sic_mauri_us, sic_mauri_spz, sic_ad)
    1468         200 :             CALL auxbas_pw_pool%give_back_pw(v_sic_rspace)
    1469       85045 :             DEALLOCATE (v_sic_rspace)
    1470             :          END SELECT
    1471       84845 :          DEALLOCATE (v_rspace_new)
    1472             : 
    1473             :       ELSE
    1474             :          ! not implemented (or at least not tested)
    1475        8792 :          CPASSERT(dft_control%sic_method_id == sic_none)
    1476        8792 :          CPASSERT(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
    1477       19652 :          DO ispin = 1, nspins
    1478             :             ! extra contribution attributed to the dependency of the dielectric constant to the charge density
    1479       10860 :             IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
    1480           0 :                dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
    1481           0 :                CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace, dvol)
    1482             :             END IF
    1483             :             ! Add SCCS contribution
    1484       10860 :             IF (dft_control%do_sccs) THEN
    1485           0 :                CALL pw_axpy(v_sccs_rspace, v_rspace)
    1486             :             END IF
    1487             :             ! DFT embedding
    1488       10860 :             IF (dft_control%apply_embed_pot) THEN
    1489          12 :                CALL pw_axpy(v_rspace_embed(ispin), v_rspace, v_rspace_embed(ispin)%pw_grid%dvol)
    1490          12 :                CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
    1491             :             END IF
    1492       19652 :             IF (lrigpw) THEN
    1493           0 :                lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
    1494           0 :                CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
    1495           0 :                DO ikind = 1, nkind
    1496           0 :                   lri_v_int(ikind)%v_int = 0.0_dp
    1497             :                END DO
    1498             :                CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
    1499           0 :                                                   lri_v_int, calculate_forces, "LRI_AUX")
    1500           0 :                DO ikind = 1, nkind
    1501           0 :                   CALL para_env%sum(lri_v_int(ikind)%v_int)
    1502             :                END DO
    1503           0 :                IF (lri_env%exact_1c_terms) THEN
    1504           0 :                   rho_ao => my_rho(ispin, :)
    1505           0 :                   ksmat => ks_matrix(ispin, :)
    1506             :                   CALL integrate_v_rspace_diagonal(v_rspace, ksmat(1)%matrix, &
    1507             :                                                    rho_ao(1)%matrix, qs_env, &
    1508           0 :                                                    calculate_forces, "ORB")
    1509             :                END IF
    1510           0 :                IF (lri_env%ppl_ri) THEN
    1511           0 :                   CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
    1512             :                END IF
    1513       10860 :             ELSEIF (rigpw) THEN
    1514           0 :                lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
    1515           0 :                CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
    1516           0 :                DO ikind = 1, nkind
    1517           0 :                   lri_v_int(ikind)%v_int = 0.0_dp
    1518             :                END DO
    1519             :                CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
    1520           0 :                                                   lri_v_int, calculate_forces, "RI_HXC")
    1521           0 :                DO ikind = 1, nkind
    1522           0 :                   CALL para_env%sum(lri_v_int(ikind)%v_int)
    1523             :                END DO
    1524             :             ELSE
    1525       10860 :                rho_ao => my_rho(ispin, :)
    1526       10860 :                ksmat => ks_matrix(ispin, :)
    1527             :                CALL integrate_v_rspace(v_rspace=v_rspace, &
    1528             :                                        pmat_kp=rho_ao, &
    1529             :                                        hmat_kp=ksmat, &
    1530             :                                        qs_env=qs_env, &
    1531             :                                        calculate_forces=calculate_forces, &
    1532       10860 :                                        gapw=gapw)
    1533             :             END IF
    1534             :          END DO
    1535             :       END IF ! ASSOCIATED(v_rspace_new)
    1536             : 
    1537             :       ! **** LRIGPW: KS matrix from integrated potential
    1538       93637 :       IF (lrigpw) THEN
    1539         428 :          CALL get_qs_env(qs_env, ks_env=ks_env)
    1540         428 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
    1541         428 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
    1542         868 :          DO ispin = 1, nspins
    1543         440 :             ksmat => ks_matrix(ispin, :)
    1544             :             CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, &
    1545         868 :                                          cell_to_index=cell_to_index)
    1546             :          END DO
    1547         428 :          IF (calculate_forces) THEN
    1548          22 :             CALL calculate_lri_forces(lri_env, lri_density, qs_env, rho_ao_kp, atomic_kind_set)
    1549             :          END IF
    1550       93209 :       ELSEIF (rigpw) THEN
    1551           0 :          CALL get_qs_env(qs_env, matrix_s=smat)
    1552           0 :          DO ispin = 1, nspins
    1553             :             CALL calculate_ri_ks_matrix(lri_env, lri_v_int, ks_matrix(ispin, 1)%matrix, &
    1554           0 :                                         smat(1)%matrix, atomic_kind_set, ispin)
    1555             :          END DO
    1556           0 :          IF (calculate_forces) THEN
    1557           0 :             rho_ao_nokp => rho_ao_kp(:, 1)
    1558           0 :             CALL calculate_ri_forces(lri_env, lri_density, qs_env, rho_ao_nokp, atomic_kind_set)
    1559             :          END IF
    1560             :       END IF
    1561             : 
    1562       93637 :       IF (ASSOCIATED(v_tau_rspace)) THEN
    1563        1840 :          IF (lrigpw .OR. rigpw) THEN
    1564           0 :             CPABORT("LRIGPW/RIGPW not implemented for meta-GGAs")
    1565             :          END IF
    1566        4066 :          DO ispin = 1, nspins
    1567        2226 :             CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
    1568             : 
    1569        2226 :             rho_ao => rho_ao_kp(ispin, :)
    1570        2226 :             ksmat => ks_matrix(ispin, :)
    1571             :             CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
    1572             :                                     pmat_kp=rho_ao, hmat_kp=ksmat, &
    1573             :                                     qs_env=qs_env, &
    1574             :                                     calculate_forces=calculate_forces, compute_tau=.TRUE., &
    1575        2226 :                                     gapw=gapw)
    1576        4066 :             CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
    1577             :          END DO
    1578        1840 :          DEALLOCATE (v_tau_rspace)
    1579             :       END IF
    1580             : 
    1581             :       ! Add contributions from ADMM if requested
    1582       93637 :       IF (dft_control%do_admm) THEN
    1583        9830 :          CALL get_qs_env(qs_env, admm_env=admm_env)
    1584             :          CALL get_admm_env(admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
    1585        9830 :                            matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
    1586        9830 :          CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
    1587        9830 :          IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN
    1588       15272 :             DO ispin = 1, nspins
    1589             :                ! Calculate the xc potential
    1590        8388 :                CALL pw_scale(v_rspace_new_aux_fit(ispin), v_rspace_new_aux_fit(ispin)%pw_grid%dvol)
    1591             : 
    1592             :                ! set matrix_ks_aux_fit_dft = matrix_ks_aux_fit(k_HF)
    1593       19476 :                DO img = 1, dft_control%nimages
    1594             :                   CALL dbcsr_copy(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit(ispin, img)%matrix, &
    1595       19476 :                                   name="DFT exch. part of matrix_ks_aux_fit")
    1596             :                END DO
    1597             : 
    1598             :                ! Add potential to ks_matrix aux_fit, skip integration if no DFT correction
    1599             : 
    1600        8388 :                IF (admm_env%aux_exch_func .NE. do_admm_aux_exch_func_none) THEN
    1601             : 
    1602             :                   !GPW by default. IF GAPW, then take relevant task list and basis
    1603        8388 :                   CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
    1604        8388 :                   basis_type = "AUX_FIT"
    1605        8388 :                   IF (admm_env%do_gapw) THEN
    1606        2154 :                      task_list => admm_env%admm_gapw_env%task_list
    1607        2154 :                      basis_type = "AUX_FIT_SOFT"
    1608             :                   END IF
    1609        8388 :                   fadm = 1.0_dp
    1610             :                   ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS,
    1611        8388 :                   IF (admm_env%do_admmp) THEN
    1612         222 :                      fadm = admm_env%gsi(ispin)**2
    1613        8166 :                   ELSE IF (admm_env%do_admms) THEN
    1614         384 :                      fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
    1615             :                   END IF
    1616             : 
    1617        8388 :                   rho_ao => rho_ao_aux(ispin, :)
    1618        8388 :                   ksmat => matrix_ks_aux_fit(ispin, :)
    1619             : 
    1620             :                   CALL integrate_v_rspace(v_rspace=v_rspace_new_aux_fit(ispin), &
    1621             :                                           pmat_kp=rho_ao, &
    1622             :                                           hmat_kp=ksmat, &
    1623             :                                           qs_env=qs_env, &
    1624             :                                           calculate_forces=calculate_forces, &
    1625             :                                           force_adm=fadm, &
    1626             :                                           gapw=.FALSE., & !even if actual GAPW calculation, want to use AUX_FIT_SOFT
    1627             :                                           basis_type=basis_type, &
    1628        8388 :                                           task_list_external=task_list)
    1629             :                END IF
    1630             : 
    1631             :                ! matrix_ks_aux_fit_dft(x_DFT)=matrix_ks_aux_fit_dft(old,k_HF)-matrix_ks_aux_fit(k_HF-x_DFT)
    1632       19476 :                DO img = 1, dft_control%nimages
    1633             :                   CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
    1634       19476 :                                  matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
    1635             :                END DO
    1636             : 
    1637       15272 :                CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
    1638             :             END DO
    1639        6884 :             DEALLOCATE (v_rspace_new_aux_fit)
    1640             :          END IF
    1641             :          ! Clean up v_tau_rspace_aux_fit, which is actually not needed
    1642        9830 :          IF (ASSOCIATED(v_tau_rspace_aux_fit)) THEN
    1643           0 :             DO ispin = 1, nspins
    1644           0 :                CALL auxbas_pw_pool%give_back_pw(v_tau_rspace_aux_fit(ispin))
    1645             :             END DO
    1646           0 :             DEALLOCATE (v_tau_rspace_aux_fit)
    1647             :          END IF
    1648             :       END IF
    1649             : 
    1650       93637 :       IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed)
    1651             : 
    1652       93637 :       CALL timestop(handle)
    1653             : 
    1654       93637 :    END SUBROUTINE sum_up_and_integrate
    1655             : 
    1656             : !**************************************************************************
    1657             : !> \brief Calculate the ZMP potential and energy as in Zhao, Morrison Parr
    1658             : !> PRA 50i, 2138 (1994)
    1659             : !> V_c^\lambda defined as int_rho-rho_0/r-r' or rho-rho_0 times a Lagrange
    1660             : !> multiplier, plus Fermi-Amaldi potential that should give the V_xc in the
    1661             : !> limit \lambda --> \infty
    1662             : !>
    1663             : !> \param qs_env ...
    1664             : !> \param v_rspace_new ...
    1665             : !> \param rho ...
    1666             : !> \param exc ...
    1667             : !> \author D. Varsano  [daniele.varsano@nano.cnr.it]
    1668             : ! **************************************************************************************************
    1669           0 :    SUBROUTINE calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
    1670             : 
    1671             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1672             :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_rspace_new
    1673             :       TYPE(qs_rho_type), POINTER                         :: rho
    1674             :       REAL(KIND=dp)                                      :: exc
    1675             : 
    1676             :       CHARACTER(*), PARAMETER :: routineN = 'calculate_zmp_potential'
    1677             : 
    1678             :       INTEGER                                            :: handle, my_val, nelectron, nspins
    1679             :       INTEGER, DIMENSION(2)                              :: nelectron_spin
    1680             :       LOGICAL                                            :: do_zmp_read, fermi_amaldi
    1681             :       REAL(KIND=dp)                                      :: lambda
    1682           0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_ext_r
    1683             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1684           0 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_ext_g, rho_g
    1685             :       TYPE(pw_env_type), POINTER                         :: pw_env
    1686             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
    1687             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1688             :       TYPE(pw_r3d_rs_type)                               :: v_xc_rspace
    1689           0 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
    1690             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    1691             :       TYPE(section_vals_type), POINTER                   :: ext_den_section, input
    1692             : 
    1693             : !, v_h_gspace, &
    1694             : 
    1695           0 :       CALL timeset(routineN, handle)
    1696           0 :       NULLIFY (auxbas_pw_pool)
    1697           0 :       NULLIFY (pw_env)
    1698           0 :       NULLIFY (poisson_env)
    1699           0 :       NULLIFY (v_rspace_new)
    1700           0 :       NULLIFY (dft_control)
    1701           0 :       NULLIFY (rho_r, rho_g, tot_rho_ext_r, rho_ext_g)
    1702             :       CALL get_qs_env(qs_env=qs_env, &
    1703             :                       pw_env=pw_env, &
    1704             :                       ks_env=ks_env, &
    1705             :                       rho=rho, &
    1706             :                       input=input, &
    1707             :                       nelectron_spin=nelectron_spin, &
    1708           0 :                       dft_control=dft_control)
    1709             :       CALL pw_env_get(pw_env=pw_env, &
    1710             :                       auxbas_pw_pool=auxbas_pw_pool, &
    1711           0 :                       poisson_env=poisson_env)
    1712           0 :       CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
    1713           0 :       nspins = 1
    1714           0 :       ALLOCATE (v_rspace_new(nspins))
    1715           0 :       CALL auxbas_pw_pool%create_pw(pw=v_rspace_new(1))
    1716           0 :       CALL auxbas_pw_pool%create_pw(pw=v_xc_rspace)
    1717             : 
    1718           0 :       CALL pw_zero(v_rspace_new(1))
    1719           0 :       do_zmp_read = dft_control%apply_external_vxc
    1720           0 :       IF (do_zmp_read) THEN
    1721           0 :          CALL pw_copy(qs_env%external_vxc, v_rspace_new(1))
    1722             :          exc = accurate_dot_product(v_rspace_new(1)%array, rho_r(1)%array)* &
    1723           0 :                v_rspace_new(1)%pw_grid%dvol
    1724             :       ELSE
    1725           0 :          BLOCK
    1726             :             REAL(KIND=dp)                                      :: factor
    1727             :             TYPE(pw_c1d_gs_type) :: rho_eff_gspace, v_xc_gspace
    1728           0 :             CALL auxbas_pw_pool%create_pw(pw=rho_eff_gspace)
    1729           0 :             CALL auxbas_pw_pool%create_pw(pw=v_xc_gspace)
    1730           0 :             CALL pw_zero(rho_eff_gspace)
    1731           0 :             CALL pw_zero(v_xc_gspace)
    1732           0 :             CALL pw_zero(v_xc_rspace)
    1733           0 :             factor = pw_integrate_function(rho_g(1))
    1734             :             CALL qs_rho_get(qs_env%rho_external, &
    1735             :                             rho_g=rho_ext_g, &
    1736           0 :                             tot_rho_r=tot_rho_ext_r)
    1737           0 :             factor = tot_rho_ext_r(1)/factor
    1738             : 
    1739           0 :             CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
    1740           0 :             CALL pw_axpy(rho_ext_g(1), rho_eff_gspace, alpha=-1.0_dp)
    1741           0 :             ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
    1742           0 :             CALL section_vals_val_get(ext_den_section, "LAMBDA", r_val=lambda)
    1743           0 :             CALL section_vals_val_get(ext_den_section, "ZMP_CONSTRAINT", i_val=my_val)
    1744           0 :             CALL section_vals_val_get(ext_den_section, "FERMI_AMALDI", l_val=fermi_amaldi)
    1745             : 
    1746           0 :             CALL pw_scale(rho_eff_gspace, a=lambda)
    1747           0 :             nelectron = nelectron_spin(1)
    1748           0 :             factor = -1.0_dp/nelectron
    1749           0 :             CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
    1750             : 
    1751           0 :             CALL pw_poisson_solve(poisson_env, rho_eff_gspace, vhartree=v_xc_gspace)
    1752           0 :             CALL pw_transfer(v_xc_gspace, v_rspace_new(1))
    1753           0 :             CALL pw_copy(v_rspace_new(1), v_xc_rspace)
    1754             : 
    1755           0 :             exc = 0.0_dp
    1756           0 :             exc = pw_integral_ab(v_rspace_new(1), rho_r(1))
    1757             : 
    1758             : !Note that this is not the xc energy but \int(\rho*v_xc)
    1759             : !Vxc---> v_rspace_new
    1760             : !Exc---> energy%exc
    1761           0 :             CALL auxbas_pw_pool%give_back_pw(rho_eff_gspace)
    1762           0 :             CALL auxbas_pw_pool%give_back_pw(v_xc_gspace)
    1763             :          END BLOCK
    1764             :       END IF
    1765             : 
    1766           0 :       CALL auxbas_pw_pool%give_back_pw(v_xc_rspace)
    1767             : 
    1768           0 :       CALL timestop(handle)
    1769             : 
    1770           0 :    END SUBROUTINE calculate_zmp_potential
    1771             : 
    1772             : ! **************************************************************************************************
    1773             : !> \brief ...
    1774             : !> \param qs_env ...
    1775             : !> \param rho ...
    1776             : !> \param v_rspace_embed ...
    1777             : !> \param dft_control ...
    1778             : !> \param embed_corr ...
    1779             : !> \param just_energy ...
    1780             : ! **************************************************************************************************
    1781         868 :    SUBROUTINE get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, &
    1782             :                                          just_energy)
    1783             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1784             :       TYPE(qs_rho_type), POINTER                         :: rho
    1785             :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_rspace_embed
    1786             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1787             :       REAL(KIND=dp)                                      :: embed_corr
    1788             :       LOGICAL                                            :: just_energy
    1789             : 
    1790             :       CHARACTER(*), PARAMETER :: routineN = 'get_embed_potential_energy'
    1791             : 
    1792             :       INTEGER                                            :: handle, ispin
    1793             :       REAL(KIND=dp)                                      :: embed_corr_local
    1794             :       TYPE(pw_env_type), POINTER                         :: pw_env
    1795             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1796         868 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r
    1797             : 
    1798         868 :       CALL timeset(routineN, handle)
    1799             : 
    1800         868 :       NULLIFY (auxbas_pw_pool)
    1801         868 :       NULLIFY (pw_env)
    1802         868 :       NULLIFY (rho_r)
    1803             :       CALL get_qs_env(qs_env=qs_env, &
    1804             :                       pw_env=pw_env, &
    1805         868 :                       rho=rho)
    1806             :       CALL pw_env_get(pw_env=pw_env, &
    1807         868 :                       auxbas_pw_pool=auxbas_pw_pool)
    1808         868 :       CALL qs_rho_get(rho, rho_r=rho_r)
    1809        3952 :       ALLOCATE (v_rspace_embed(dft_control%nspins))
    1810             : 
    1811         868 :       embed_corr = 0.0_dp
    1812             : 
    1813        2216 :       DO ispin = 1, dft_control%nspins
    1814        1348 :          CALL auxbas_pw_pool%create_pw(pw=v_rspace_embed(ispin))
    1815        1348 :          CALL pw_zero(v_rspace_embed(ispin))
    1816             : 
    1817        1348 :          CALL pw_copy(qs_env%embed_pot, v_rspace_embed(ispin))
    1818        1348 :          embed_corr_local = 0.0_dp
    1819             : 
    1820             :          ! Spin embedding potential in open-shell case
    1821        1348 :          IF (dft_control%nspins .EQ. 2) THEN
    1822         960 :             IF (ispin .EQ. 1) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), 1.0_dp)
    1823         960 :             IF (ispin .EQ. 2) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), -1.0_dp)
    1824             :          END IF
    1825             :          ! Integrate the density*potential
    1826        1348 :          embed_corr_local = pw_integral_ab(v_rspace_embed(ispin), rho_r(ispin))
    1827             : 
    1828        2216 :          embed_corr = embed_corr + embed_corr_local
    1829             : 
    1830             :       END DO
    1831             : 
    1832             :       ! If only energy requiested we delete the potential
    1833         868 :       IF (just_energy) THEN
    1834         692 :          DO ispin = 1, dft_control%nspins
    1835         692 :             CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
    1836             :          END DO
    1837         286 :          DEALLOCATE (v_rspace_embed)
    1838             :       END IF
    1839             : 
    1840         868 :       CALL timestop(handle)
    1841             : 
    1842         868 :    END SUBROUTINE get_embed_potential_energy
    1843             : 
    1844             : END MODULE qs_ks_utils

Generated by: LCOV version 1.15