LCOV - code coverage report
Current view: top level - src - qs_ks_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 80.9 % 801 648
Test Date: 2025-07-25 12:55:17 Functions: 88.9 % 9 8

            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       103397 :    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       103397 :       INTEGER, DIMENSION(:), POINTER                     :: ivec
     162       103397 :       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       103397 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: energy_scaling, rvec, scaling
     168              :       TYPE(cell_type), POINTER                           :: cell
     169       103397 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_h, matrix_hfx, matrix_p, mdummy, &
     170       103397 :                                                             mo_derivs, rho_ao
     171       103397 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p2
     172              :       TYPE(dbcsr_type), POINTER                          :: dbcsr_deriv, fm_deriv, fm_scaled, &
     173              :                                                             mo_coeff
     174       103397 :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     175       103397 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     176              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     177       103397 :       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       103397 :       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       103035 :       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       105207 :    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       103397 :    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       103397 :       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       103397 :       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       103397 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mo_derivs, rho_ao, tmp_dbcsr
     488              :       TYPE(dbcsr_type), POINTER                          :: orb_density_matrix, orb_h
     489       103397 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     490              :       TYPE(pw_c1d_gs_type)                               :: work_v_gspace
     491       103397 :       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       103397 :       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       103397 :       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       103557 :    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       103397 :    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       103397 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: store_forces
     772       103397 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     773              :       TYPE(pw_c1d_gs_type)                               :: work_rho, work_v
     774       103397 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     775       103397 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     776              : 
     777       103397 :       NULLIFY (mo_array, rho_g)
     778              : 
     779       103397 :       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       103693 :    END SUBROUTINE calc_v_sic_rspace
     871              : 
     872              : ! **************************************************************************************************
     873              : !> \brief ...
     874              : !> \param qs_env ...
     875              : !> \param rho ...
     876              : ! **************************************************************************************************
     877       206750 :    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       103375 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_r_arr
     885              :       TYPE(cell_type), POINTER                           :: cell
     886              :       TYPE(cp_logger_type), POINTER                      :: logger
     887       103375 :       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       103375 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     891              :       TYPE(section_vals_type), POINTER                   :: input, scf_section
     892              : 
     893       103375 :       NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
     894       103375 :                dft_control, tot_rho_r_arr, rho_ao)
     895              : 
     896       206750 :       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       103375 :                       dft_control=dft_control)
     904              : 
     905       103375 :       CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
     906              : 
     907       103375 :       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       103375 :                                          extension=".scfLog")
     910              : 
     911       103375 :       CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
     912       103375 :       n_electrons = n_electrons - dft_control%charge
     913       103375 :       tot_rho_r = accurate_sum(tot_rho_r_arr)
     914              : 
     915       103375 :       trace = 0
     916       103375 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN
     917         4130 :          DO ispin = 1, dft_control%nspins
     918         6728 :             DO img = 1, dft_control%nimages
     919         2598 :                CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
     920         5062 :                trace = trace + trace_tmp
     921              :             END DO
     922              :          END DO
     923              :       END IF
     924              : 
     925       103375 :       IF (output_unit > 0) THEN
     926          833 :          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          833 :             "Electronic density on regular grids: ", &
     929          833 :             tot_rho_r, &
     930              :             tot_rho_r + &
     931          833 :             REAL(n_electrons, dp), &
     932          833 :             "Core density on regular grids:", &
     933          833 :             qs_charges%total_rho_core_rspace, &
     934         1666 :             qs_charges%total_rho_core_rspace - REAL(n_electrons + dft_control%charge, dp)
     935              :       END IF
     936       103375 :       IF (dft_control%qs_control%gapw) THEN
     937        13718 :          tot1_h = qs_charges%total_rho1_hard(1)
     938        13718 :          tot1_s = qs_charges%total_rho1_soft(1)
     939        17048 :          DO ispin = 2, dft_control%nspins
     940         3330 :             tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
     941        17048 :             tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
     942              :          END DO
     943        13718 :          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        13718 :                                  qs_charges%total_rho_core_rspace
     958        89657 :       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        87049 :          IF (output_unit > 0) THEN
     977              :             WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
     978          434 :                "Total charge density on r-space grids:     ", &
     979              :                tot_rho_r + &
     980          434 :                qs_charges%total_rho_core_rspace, &
     981          434 :                "Total charge density g-space grids:     ", &
     982          868 :                qs_charges%total_rho_gspace
     983              :          END IF
     984              :          qs_charges%background = tot_rho_r + &
     985        87049 :                                  qs_charges%total_rho_core_rspace
     986              :       END IF
     987       103375 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="()")
     988       103375 :       qs_charges%background = qs_charges%background/cell%deth
     989              : 
     990              :       CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
     991       103375 :                                         "PRINT%TOTAL_DENSITIES")
     992              : 
     993       103375 :    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       103375 :    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       103375 :       logger => cp_get_default_logger()
    1022              : 
    1023       103375 :       NULLIFY (pw_env)
    1024       103375 :       CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
    1025       103375 :       psolver = pw_env%poisson_env%parameters%solver
    1026              : 
    1027              :       output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", &
    1028       103375 :                                          extension=".scfLog")
    1029       103375 :       IF (output_unit > 0) THEN
    1030          486 :          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          486 :          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          486 :             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          486 :             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          486 :                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          486 :                      "Core Hamiltonian energy:                       ", energy%core, &
    1095          486 :                      "Hartree energy:                                ", energy%hartree, &
    1096          972 :                      "Exchange-correlation energy:                   ", energy%exc
    1097              :                END IF
    1098              :             END IF
    1099              :          END IF
    1100              : 
    1101          486 :          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          486 :          IF (energy%e_hartree /= 0.0_dp) &
    1107              :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1108          454 :             "Coulomb (electron-electron) energy:            ", energy%e_hartree
    1109          486 :          IF (energy%dispersion /= 0.0_dp) &
    1110              :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1111            0 :             "Dispersion energy:                             ", energy%dispersion
    1112          486 :          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          486 :          IF (energy%gcp /= 0.0_dp) &
    1116              :             WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
    1117            0 :             "gCP energy:                                    ", energy%gcp
    1118          486 :          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          486 :          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          486 :          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          486 :          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          486 :          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          486 :          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          486 :          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       103375 :                                         "DFT%SCF%PRINT%DETAILED_ENERGY")
    1160              : 
    1161       103375 :    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            6 :    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            6 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks
    1184              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1185              : 
    1186            6 :       CALL timeset(routineN, handle)
    1187              : 
    1188              :       ! create the matrix using matrix_ks as a template
    1189            6 :       IF (ASSOCIATED(matrix_vxc)) THEN
    1190            0 :          CALL dbcsr_deallocate_matrix_set(matrix_vxc)
    1191              :       END IF
    1192            6 :       CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
    1193           28 :       ALLOCATE (matrix_vxc(SIZE(matrix_ks)))
    1194           16 :       DO ispin = 1, SIZE(matrix_ks)
    1195           10 :          NULLIFY (matrix_vxc(ispin)%matrix)
    1196           10 :          CALL dbcsr_init_p(matrix_vxc(ispin)%matrix)
    1197              :          CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
    1198           10 :                          name="Matrix VXC of spin "//cp_to_string(ispin))
    1199           16 :          CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
    1200              :       END DO
    1201              : 
    1202              :       ! and integrate
    1203            6 :       CALL get_qs_env(qs_env, dft_control=dft_control)
    1204            6 :       gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
    1205           16 :       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           10 :                                  gapw=gapw)
    1211              :          ! scale by the volume element... should really become part of integrate_v_rspace
    1212           16 :          CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
    1213              :       END DO
    1214              : 
    1215            6 :       CALL timestop(handle)
    1216              : 
    1217            6 :    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        93997 :    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        93997 :       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        93997 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1273        93997 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ksmat, rho_ao, rho_ao_nokp, smat
    1274        93997 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_aux_fit, &
    1275        93997 :                                                             matrix_ks_aux_fit_dft, rho_ao_aux, &
    1276        93997 :                                                             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        93997 :       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        93997 :       CALL timeset(routineN, handle)
    1292              : 
    1293        93997 :       NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
    1294        93997 :                v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
    1295        93997 :                ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
    1296        93997 :                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        93997 :                       vee=vee)
    1303              : 
    1304        93997 :       CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
    1305        93997 :       CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
    1306        93997 :       gapw = dft_control%qs_control%gapw
    1307        93997 :       gapw_xc = dft_control%qs_control%gapw_xc
    1308        93997 :       do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
    1309              : 
    1310        93997 :       rigpw = dft_control%qs_control%rigpw
    1311        93997 :       lrigpw = dft_control%qs_control%lrigpw
    1312        93997 :       IF (lrigpw .OR. rigpw) THEN
    1313              :          CALL get_qs_env(qs_env, &
    1314              :                          lri_env=lri_env, &
    1315              :                          lri_density=lri_density, &
    1316          420 :                          atomic_kind_set=atomic_kind_set)
    1317              :       END IF
    1318              : 
    1319        93997 :       nspins = dft_control%nspins
    1320              : 
    1321              :       ! sum up potentials and integrate
    1322        93997 :       IF (ASSOCIATED(v_rspace_new)) THEN
    1323       185861 :          DO ispin = 1, nspins
    1324       100810 :             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        98014 :                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       100810 :             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       100810 :             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       100810 :             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       100810 :             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       100810 :             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       100810 :             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       101170 :             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       100810 :                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       100810 :             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       100810 :             IF (lrigpw) THEN
    1422          432 :                lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
    1423          432 :                CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
    1424         1298 :                DO ikind = 1, nkind
    1425       280382 :                   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          432 :                                                   lri_v_int, calculate_forces, "LRI_AUX")
    1429         1298 :                DO ikind = 1, nkind
    1430       559466 :                   CALL para_env%sum(lri_v_int(ikind)%v_int)
    1431              :                END DO
    1432          432 :                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          432 :                IF (lri_env%ppl_ri) THEN
    1440            8 :                   CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
    1441              :                END IF
    1442       100378 :             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       100378 :                rho_ao => my_rho(ispin, :)
    1455       100378 :                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       100378 :                                        gapw=gapw)
    1461              :             END IF
    1462       185861 :             CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
    1463              :          END DO ! ispin
    1464              : 
    1465        85251 :          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        85251 :             DEALLOCATE (v_sic_rspace)
    1470              :          END SELECT
    1471        85051 :          DEALLOCATE (v_rspace_new)
    1472              : 
    1473              :       ELSE
    1474              :          ! not implemented (or at least not tested)
    1475         8946 :          CPASSERT(dft_control%sic_method_id == sic_none)
    1476         8946 :          CPASSERT(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
    1477        20006 :          DO ispin = 1, nspins
    1478              :             ! extra contribution attributed to the dependency of the dielectric constant to the charge density
    1479        11060 :             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        11060 :             IF (dft_control%do_sccs) THEN
    1485            0 :                CALL pw_axpy(v_sccs_rspace, v_rspace)
    1486              :             END IF
    1487              :             ! DFT embedding
    1488        11060 :             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        20006 :             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        11060 :             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        11060 :                rho_ao => my_rho(ispin, :)
    1526        11060 :                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        11060 :                                        gapw=gapw)
    1533              :             END IF
    1534              :          END DO
    1535              :       END IF ! ASSOCIATED(v_rspace_new)
    1536              : 
    1537              :       ! **** LRIGPW: KS matrix from integrated potential
    1538        93997 :       IF (lrigpw) THEN
    1539          420 :          CALL get_qs_env(qs_env, ks_env=ks_env)
    1540          420 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
    1541          420 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
    1542          852 :          DO ispin = 1, nspins
    1543          432 :             ksmat => ks_matrix(ispin, :)
    1544              :             CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, &
    1545          852 :                                          cell_to_index=cell_to_index)
    1546              :          END DO
    1547          420 :          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        93577 :       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        93997 :       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        93997 :       IF (dft_control%do_admm) THEN
    1583         9910 :          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         9910 :                            matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
    1586         9910 :          CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
    1587         9910 :          IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN
    1588        15432 :             DO ispin = 1, nspins
    1589              :                ! Calculate the xc potential
    1590         8468 :                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        21556 :                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        21556 :                                   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         8468 :                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         8468 :                   CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
    1604         8468 :                   basis_type = "AUX_FIT"
    1605         8468 :                   IF (admm_env%do_gapw) THEN
    1606         2234 :                      task_list => admm_env%admm_gapw_env%task_list
    1607         2234 :                      basis_type = "AUX_FIT_SOFT"
    1608              :                   END IF
    1609         8468 :                   fadm = 1.0_dp
    1610              :                   ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS,
    1611         8468 :                   IF (admm_env%do_admmp) THEN
    1612          222 :                      fadm = admm_env%gsi(ispin)**2
    1613         8246 :                   ELSE IF (admm_env%do_admms) THEN
    1614          432 :                      fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
    1615              :                   END IF
    1616              : 
    1617         8468 :                   rho_ao => rho_ao_aux(ispin, :)
    1618         8468 :                   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         8468 :                                           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        21556 :                DO img = 1, dft_control%nimages
    1633              :                   CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
    1634        21556 :                                  matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
    1635              :                END DO
    1636              : 
    1637        15432 :                CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
    1638              :             END DO
    1639         6964 :             DEALLOCATE (v_rspace_new_aux_fit)
    1640              :          END IF
    1641              :          ! Clean up v_tau_rspace_aux_fit, which is actually not needed
    1642         9910 :          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        93997 :       IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed)
    1651              : 
    1652        93997 :       CALL timestop(handle)
    1653              : 
    1654        93997 :    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 2.0-1