LCOV - code coverage report
Current view: top level - src - hfx_exx.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 94.0 % 301 283
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            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 to calculate EXX in RPA and energy correction methods
      10              : !> \par History
      11              : !>      07.2020 separated from mp2.F [F. Stein, code by Jan Wilhelm]
      12              : !>      06.2022 EXX contribution to the forces [A. Bussy]
      13              : !>      03.2023 Generalized for energy correction methods
      14              : !> \author Jan Wilhelm, Frederick Stein, Augustin Bussy, Fabian Belleflamme
      15              : ! **************************************************************************************************
      16              : MODULE hfx_exx
      17              :    USE admm_methods,                    ONLY: admm_projection_derivative
      18              :    USE admm_types,                      ONLY: admm_env_create,&
      19              :                                               admm_env_release,&
      20              :                                               admm_type,&
      21              :                                               get_admm_env
      22              :    USE cp_control_types,                ONLY: dft_control_type
      23              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      24              :                                               dbcsr_copy,&
      25              :                                               dbcsr_create,&
      26              :                                               dbcsr_p_type,&
      27              :                                               dbcsr_release,&
      28              :                                               dbcsr_scale,&
      29              :                                               dbcsr_set,&
      30              :                                               dbcsr_type
      31              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      32              :                                               copy_fm_to_dbcsr,&
      33              :                                               dbcsr_allocate_matrix_set,&
      34              :                                               dbcsr_deallocate_matrix_set
      35              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      36              :                                               cp_logger_get_default_unit_nr,&
      37              :                                               cp_logger_type
      38              :    USE hfx_admm_utils,                  ONLY: create_admm_xc_section,&
      39              :                                               tddft_hfx_matrix
      40              :    USE hfx_derivatives,                 ONLY: derivatives_four_center
      41              :    USE hfx_energy_potential,            ONLY: integrate_four_center
      42              :    USE hfx_ri,                          ONLY: hfx_ri_update_forces,&
      43              :                                               hfx_ri_update_ks
      44              :    USE hfx_types,                       ONLY: hfx_type
      45              :    USE input_constants,                 ONLY: do_admm_aux_exch_func_none
      46              :    USE input_section_types,             ONLY: section_vals_create,&
      47              :                                               section_vals_duplicate,&
      48              :                                               section_vals_get,&
      49              :                                               section_vals_get_subs_vals,&
      50              :                                               section_vals_release,&
      51              :                                               section_vals_set_subs_vals,&
      52              :                                               section_vals_type,&
      53              :                                               section_vals_val_get
      54              :    USE kinds,                           ONLY: dp
      55              :    USE machine,                         ONLY: m_walltime
      56              :    USE message_passing,                 ONLY: mp_para_env_type
      57              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      58              :    USE pw_env_types,                    ONLY: pw_env_get,&
      59              :                                               pw_env_type
      60              :    USE pw_methods,                      ONLY: pw_scale
      61              :    USE pw_pool_types,                   ONLY: pw_pool_type
      62              :    USE pw_types,                        ONLY: pw_r3d_rs_type
      63              :    USE qs_energy_types,                 ONLY: qs_energy_type
      64              :    USE qs_environment_types,            ONLY: get_qs_env,&
      65              :                                               qs_environment_type
      66              :    USE qs_integrate_potential,          ONLY: integrate_v_rspace
      67              :    USE qs_ks_reference,                 ONLY: ks_ref_potential
      68              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      69              :                                               qs_rho_type
      70              :    USE qs_vxc,                          ONLY: qs_vxc_create
      71              :    USE task_list_types,                 ONLY: task_list_type
      72              :    USE virial_types,                    ONLY: virial_type
      73              : 
      74              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
      75              : 
      76              : #include "./base/base_uses.f90"
      77              : 
      78              :    IMPLICIT NONE
      79              : 
      80              :    PRIVATE
      81              : 
      82              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_exx'
      83              : 
      84              :    PUBLIC :: calculate_exx, add_exx_to_rhs, calc_exx_admm_xc_contributions, exx_pre_hfx, exx_post_hfx
      85              : 
      86              : CONTAINS
      87              : 
      88              : ! **************************************************************************************************
      89              : !> \brief ...
      90              : !> \param qs_env ...
      91              : !> \param unit_nr ...
      92              : !> \param hfx_sections ...
      93              : !> \param x_data ...
      94              : !> \param do_gw ...
      95              : !> \param do_admm ...
      96              : !> \param calc_forces ...
      97              : !> \param reuse_hfx ...
      98              : !> \param do_im_time ...
      99              : !> \param E_ex_from_GW ...
     100              : !> \param E_admm_from_GW ...
     101              : !> \param t3 ...
     102              : ! **************************************************************************************************
     103          444 :    SUBROUTINE calculate_exx(qs_env, unit_nr, hfx_sections, x_data, &
     104              :                             do_gw, do_admm, calc_forces, reuse_hfx, do_im_time, &
     105              :                             E_ex_from_GW, E_admm_from_GW, t3)
     106              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     107              :       INTEGER, INTENT(IN)                                :: unit_nr
     108              :       TYPE(section_vals_type), POINTER                   :: hfx_sections
     109              :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     110              :       LOGICAL, INTENT(IN)                                :: do_gw, do_admm, calc_forces, reuse_hfx, &
     111              :                                                             do_im_time
     112              :       REAL(KIND=dp), INTENT(IN)                          :: E_ex_from_GW, E_admm_from_GW(2), t3
     113              : 
     114              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calculate_exx'
     115              : 
     116              :       INTEGER                                            :: handle, i, irep, ispin, mspin, n_rep_hf, &
     117              :                                                             nspins
     118              :       LOGICAL                                            :: calc_ints, hfx_treat_lsd_in_core, &
     119              :                                                             use_virial
     120              :       REAL(KIND=dp)                                      :: eh1, ehfx, t1, t2, tf1, tf2
     121          222 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_ks_aux_fit, rho_ao, &
     122          222 :                                                             rho_ao_aux_fit, rho_ao_resp
     123          222 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_2d, rho_ao_2d
     124              :       TYPE(dft_control_type), POINTER                    :: dft_control
     125              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     126              :       TYPE(qs_energy_type), POINTER                      :: energy
     127              :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit
     128              :       TYPE(section_vals_type), POINTER                   :: input
     129              :       TYPE(virial_type), POINTER                         :: virial
     130              : 
     131          222 :       CALL timeset(routineN, handle)
     132              : 
     133          222 :       t1 = m_walltime()
     134              : 
     135          222 :       NULLIFY (input, para_env, matrix_ks, matrix_ks_aux_fit, rho, rho_ao, virial, &
     136          222 :                dft_control, rho_aux_fit, rho_ao_aux_fit)
     137              : 
     138          222 :       CALL exx_pre_hfx(hfx_sections, x_data, reuse_hfx)
     139              : 
     140              :       CALL get_qs_env(qs_env=qs_env, &
     141              :                       input=input, &
     142              :                       para_env=para_env, &
     143              :                       energy=energy, &
     144              :                       rho=rho, &
     145              :                       matrix_ks=matrix_ks, &
     146              :                       virial=virial, &
     147          222 :                       dft_control=dft_control)
     148          222 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     149              : 
     150          222 :       IF (do_admm) THEN
     151              :          CALL get_admm_env(qs_env%admm_env, &
     152              :                            matrix_ks_aux_fit=matrix_ks_aux_fit, &
     153           42 :                            rho_aux_fit=rho_aux_fit)
     154           42 :          CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
     155              : 
     156           42 :          IF (qs_env%admm_env%do_gapw) THEN
     157            0 :             CPABORT("ADMM EXX only implmented with GPW")
     158              :          END IF
     159              :       END IF
     160              : 
     161          222 :       CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
     162              :       CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
     163          222 :                                 i_rep_section=1)
     164              : 
     165              :       ! put matrix_ks to zero
     166          474 :       DO i = 1, SIZE(matrix_ks)
     167          252 :          CALL dbcsr_set(matrix_ks(i)%matrix, 0.0_dp)
     168          474 :          IF (do_admm) THEN
     169           48 :             CALL dbcsr_set(matrix_ks_aux_fit(i)%matrix, 0.0_dp)
     170              :          END IF
     171              :       END DO
     172              : 
     173              :       ! take the exact exchange energy from GW or calculate it
     174          222 :       IF (do_gw) THEN
     175              : 
     176           56 :          IF (calc_forces) CPABORT("Not implemented")
     177              : 
     178           56 :          IF (qs_env%mp2_env%ri_g0w0%update_xc_energy) THEN
     179           24 :             CALL remove_exc_energy(energy)
     180           24 :             energy%total = energy%total + E_ex_from_GW
     181           24 :             energy%ex = E_ex_from_GW
     182           24 :             IF (do_admm) THEN
     183            4 :                energy%total = energy%total + E_admm_from_GW(1) + E_admm_from_GW(2)
     184            4 :                energy%exc = E_admm_from_GW(1)
     185            4 :                energy%exc_aux_fit = E_admm_from_GW(2)
     186              :             END IF
     187           24 :             t2 = m_walltime()
     188              : 
     189           24 :             IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Time=', t2 - t1 + t3
     190           24 :             IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX energy  =   ', energy%ex
     191           24 :             IF (do_admm .AND. unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') &
     192            2 :                'EXX ADMM XC correction  =   ', E_admm_from_GW(1) + E_admm_from_GW(2)
     193              :          END IF
     194              : 
     195              :       ELSE
     196              : 
     197          166 :          CALL remove_exc_energy(energy)
     198              : 
     199          166 :          nspins = dft_control%nspins
     200          166 :          mspin = 1
     201          166 :          IF (hfx_treat_lsd_in_core) mspin = nspins
     202              : 
     203          166 :          calc_ints = .TRUE.
     204          166 :          IF (reuse_hfx) calc_ints = .FALSE.
     205          166 :          IF (calc_forces .AND. do_im_time) calc_ints = .FALSE.
     206              : 
     207          166 :          ehfx = 0.0_dp
     208          166 :          IF (do_admm) THEN
     209           34 :             matrix_ks_2d(1:nspins, 1:1) => matrix_ks_aux_fit(1:nspins)
     210           34 :             rho_ao_2d(1:nspins, 1:1) => rho_ao_aux_fit(1:nspins)
     211              :          ELSE
     212          132 :             matrix_ks_2d(1:nspins, 1:1) => matrix_ks(1:nspins)
     213          132 :             rho_ao_2d(1:nspins, 1:1) => rho_ao(1:nspins)
     214              :          END IF
     215              : 
     216          332 :          DO irep = 1, n_rep_hf
     217              : 
     218          332 :             IF (x_data(irep, 1)%do_hfx_ri) THEN
     219              :                CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, matrix_ks_2d, ehfx, &
     220              :                                      rho_ao=rho_ao_2d, geometry_did_change=calc_ints, nspins=nspins, &
     221            8 :                                      hf_fraction=x_data(irep, 1)%general_parameter%fraction)
     222              :             ELSE
     223              : 
     224          316 :                DO ispin = 1, mspin
     225              : 
     226              :                   CALL integrate_four_center(qs_env, x_data, matrix_ks_2d, eh1, &
     227              :                                              rho_ao_2d, hfx_sections, para_env, &
     228          158 :                                              calc_ints, irep, .TRUE., ispin=ispin)
     229          316 :                   ehfx = ehfx + eh1
     230              :                END DO
     231              :             END IF
     232              :          END DO
     233              : 
     234              :          ! include the EXX contribution to the total energy
     235          166 :          energy%ex = ehfx
     236          166 :          energy%total = energy%total + energy%ex
     237              : 
     238          166 :          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     239          166 :          IF (use_virial) THEN
     240            0 :             virial%pv_calculate = .TRUE.
     241            0 :             virial%pv_fock_4c = 0.0_dp
     242              :          END IF
     243              : 
     244          166 :          IF (calc_forces) THEN
     245           58 :             tf1 = m_walltime()
     246              : 
     247              :             !Note: no need to remove xc forces: they are not even calculated in the first place
     248           58 :             NULLIFY (rho_ao_resp)
     249              : 
     250          116 :             DO irep = 1, n_rep_hf
     251              : 
     252          116 :                IF (x_data(irep, 1)%do_hfx_ri) THEN
     253              : 
     254              :                   CALL hfx_ri_update_forces(qs_env, x_data(irep, 1)%ri_data, nspins, &
     255              :                                             x_data(irep, 1)%general_parameter%fraction, &
     256              :                                             rho_ao=rho_ao_2d, rho_ao_resp=rho_ao_resp, &
     257            4 :                                             use_virial=use_virial)
     258              :                ELSE
     259              : 
     260              :                   CALL derivatives_four_center(qs_env, rho_ao_2d, rho_ao_resp, &
     261              :                                                hfx_sections, para_env, irep, &
     262           54 :                                                use_virial, external_x_data=x_data)
     263              : 
     264              :                END IF
     265              : 
     266              :             END DO !irep
     267              : 
     268           58 :             tf2 = m_walltime()
     269           58 :             IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Force Time=', tf2 - tf1
     270              :          END IF
     271              : 
     272          166 :          IF (use_virial) THEN
     273            0 :             virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
     274            0 :             virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
     275              :          END IF
     276              : 
     277              :          ! ADMM XC correction
     278          166 :          IF (do_admm) THEN
     279              : 
     280              :             CALL calc_exx_admm_xc_contributions(qs_env, matrix_ks, matrix_ks_aux_fit, x_data, &
     281              :                                                 energy%exc, energy%exc_aux_fit, calc_forces, &
     282           34 :                                                 use_virial)
     283              : 
     284              :             ! ADMM overlap forces
     285           34 :             IF (calc_forces) CALL admm_projection_derivative(qs_env, matrix_ks_aux_fit, rho_ao)
     286              : 
     287           34 :             energy%total = energy%total + energy%exc_aux_fit
     288           34 :             energy%total = energy%total + energy%exc
     289              : 
     290              :          END IF
     291              : 
     292          166 :          IF (use_virial) THEN
     293            0 :             virial%pv_calculate = .FALSE.
     294              :          END IF
     295              : 
     296          166 :          t2 = m_walltime()
     297              : 
     298          166 :          IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Time=', t2 - t1 + t3
     299          166 :          IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX energy  =   ', energy%ex
     300          166 :          IF (do_admm .AND. unit_nr > 0) THEN
     301           17 :             WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX ADMM XC correction  =   ', energy%exc + energy%exc_aux_fit
     302              :          END IF
     303              :       END IF
     304              : 
     305          222 :       CALL exx_post_hfx(qs_env, x_data, reuse_hfx)
     306              : 
     307          222 :       CALL timestop(handle)
     308              : 
     309          222 :    END SUBROUTINE calculate_exx
     310              : 
     311              : ! **************************************************************************************************
     312              : !> \brief Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian
     313              : !> \param rhs ...
     314              : !> \param qs_env ...
     315              : !> \param ext_hfx_section ...
     316              : !> \param x_data ...
     317              : !> \param recalc_integrals ...
     318              : !> \param do_admm ...
     319              : !> \param do_ec ...
     320              : !> \param do_exx ...
     321              : !> \param reuse_hfx ...
     322              : ! **************************************************************************************************
     323          230 :    SUBROUTINE add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, &
     324              :                              recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
     325              : 
     326              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN)       :: rhs
     327              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     328              :       TYPE(section_vals_type), POINTER                   :: ext_hfx_section
     329              :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     330              :       LOGICAL, INTENT(IN), OPTIONAL                      :: recalc_integrals, do_admm, do_ec, &
     331              :                                                             do_exx, reuse_hfx
     332              : 
     333              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'add_exx_to_rhs'
     334              : 
     335              :       INTEGER                                            :: handle, ispin, nao, nao_aux, nspins, &
     336              :                                                             unit_nr
     337              :       LOGICAL                                            :: calc_ints, do_hfx, my_do_ec, my_do_exx, &
     338              :                                                             my_recalc_integrals
     339              :       REAL(dp)                                           :: dummy_real1, dummy_real2
     340              :       TYPE(admm_type), POINTER                           :: admm_env
     341              :       TYPE(cp_logger_type), POINTER                      :: logger
     342          230 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: dbcsr_work, matrix_ks, matrix_s_aux, &
     343          230 :                                                             rho_ao, rho_ao_aux, work_admm
     344              :       TYPE(dbcsr_type)                                   :: dbcsr_tmp
     345              :       TYPE(dft_control_type), POINTER                    :: dft_control
     346              :       TYPE(pw_env_type), POINTER                         :: pw_env
     347              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     348              :       TYPE(pw_r3d_rs_type)                               :: vh_rspace
     349          230 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: vadmm_rspace, vtau_rspace, vxc_rspace
     350              :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit
     351              :       TYPE(section_vals_type), POINTER                   :: hfx_section
     352              :       TYPE(task_list_type), POINTER                      :: task_list_aux_fit
     353              : 
     354          230 :       NULLIFY (dbcsr_work, matrix_ks, rho, hfx_section, pw_env, vadmm_rspace, vtau_rspace, vxc_rspace, &
     355          230 :                auxbas_pw_pool, dft_control, rho_ao, rho_aux_fit, rho_ao_aux, work_admm, &
     356          230 :                matrix_s_aux, admm_env, task_list_aux_fit)
     357              : 
     358          460 :       logger => cp_get_default_logger()
     359          230 :       IF (logger%para_env%is_source()) THEN
     360          124 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     361              :       ELSE
     362              :          unit_nr = -1
     363              :       END IF
     364              : 
     365          230 :       CALL timeset(routineN, handle)
     366              : 
     367          230 :       my_recalc_integrals = .FALSE.
     368          230 :       IF (PRESENT(recalc_integrals)) my_recalc_integrals = recalc_integrals
     369          230 :       my_do_ec = .FALSE.
     370          230 :       IF (PRESENT(do_ec)) my_do_ec = do_ec
     371          230 :       my_do_exx = .TRUE.
     372          230 :       IF (PRESENT(do_exx)) my_do_exx = do_exx
     373              : 
     374              :       ! Strategy: we take the ks_matrix, remove the current xc contribution, and then add the RPA HF one
     375          230 :       CALL get_qs_env(qs_env, matrix_ks=matrix_ks, rho=rho, pw_env=pw_env, dft_control=dft_control)
     376          230 :       nspins = dft_control%nspins
     377              : 
     378              :       ! do_exx ; subtract XC and EX
     379              :       ! do_dcdft: subtract EX
     380              : 
     381          230 :       CALL dbcsr_allocate_matrix_set(dbcsr_work, nspins)
     382          468 :       DO ispin = 1, nspins
     383          238 :          ALLOCATE (dbcsr_work(ispin)%matrix)
     384          238 :          CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix)
     385          468 :          CALL dbcsr_set(dbcsr_work(ispin)%matrix, 0.0_dp)
     386              :       END DO
     387              : 
     388          230 :       IF (dft_control%do_admm) THEN
     389           80 :          CALL get_qs_env(qs_env, admm_env=admm_env)
     390              :          CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux, task_list_aux_fit=task_list_aux_fit, &
     391           80 :                            rho_aux_fit=rho_aux_fit)
     392           80 :          CALL dbcsr_allocate_matrix_set(work_admm, nspins)
     393          164 :          DO ispin = 1, nspins
     394           84 :             ALLOCATE (work_admm(ispin)%matrix)
     395           84 :             CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
     396           84 :             CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
     397          164 :             CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
     398              :          END DO
     399           80 :          CALL dbcsr_create(dbcsr_tmp, template=matrix_ks(1)%matrix)
     400           80 :          CALL dbcsr_copy(dbcsr_tmp, matrix_ks(1)%matrix)
     401              : 
     402           80 :          nao = admm_env%nao_orb
     403           80 :          nao_aux = admm_env%nao_aux_fit
     404           80 :          CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
     405              :       END IF
     406              : 
     407          230 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     408              : 
     409              :       ! Remove the standard XC + HFX contribution
     410          230 :       CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, dummy_real1, dummy_real2)
     411              : 
     412              :       ! Only remove standard XC and/or HFX
     413              :       ! (due to different requirements for RHS)
     414          230 :       IF (my_do_exx) THEN
     415              : 
     416           60 :          DO ispin = 1, nspins
     417           60 :             CALL dbcsr_copy(dbcsr_work(ispin)%matrix, matrix_ks(ispin)%matrix)
     418              :          END DO
     419              : 
     420           60 :          DO ispin = 1, nspins
     421           34 :             CALL pw_scale(vxc_rspace(ispin), -1.0_dp)
     422              :             CALL integrate_v_rspace(v_rspace=vxc_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, &
     423           34 :                                     calculate_forces=.FALSE.)
     424           60 :             IF (ASSOCIATED(vtau_rspace)) THEN
     425            0 :                CALL pw_scale(vtau_rspace(ispin), -1.0_dp)
     426              :                CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), hmat=dbcsr_work(ispin), qs_env=qs_env, &
     427            0 :                                        calculate_forces=.FALSE., compute_tau=.TRUE.)
     428              :             END IF
     429              :          END DO
     430              : 
     431              :       END IF ! do_exx
     432              : 
     433              :       ! Remove standard HFX
     434          230 :       IF (my_do_exx .OR. my_do_ec) THEN
     435              : 
     436          198 :          hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
     437          198 :          CALL section_vals_get(hfx_section, explicit=do_hfx)
     438          198 :          IF (do_hfx) THEN
     439          100 :             IF (dft_control%do_admm) THEN
     440              : 
     441              :                !note: factor -1.0 taken care of later, after HFX ADMM contribution is taken
     442           56 :                IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
     443           76 :                   DO ispin = 1, nspins
     444              :                      CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=work_admm(ispin), &
     445              :                                              qs_env=qs_env, calculate_forces=.FALSE., basis_type="AUX_FIT", &
     446           76 :                                              task_list_external=task_list_aux_fit)
     447              :                   END DO
     448              :                END IF
     449              : 
     450           56 :                CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., my_recalc_integrals)
     451              : 
     452          116 :                DO ispin = 1, nspins
     453           60 :                   CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
     454              :                   CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
     455           60 :                                      0.0_dp, admm_env%work_aux_orb)
     456              :                   CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
     457           60 :                                      0.0_dp, admm_env%work_orb_orb)
     458           60 :                   CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_tmp, keep_sparsity=.TRUE.)
     459          116 :                   CALL dbcsr_add(dbcsr_work(ispin)%matrix, dbcsr_tmp, 1.0_dp, -1.0_dp)
     460              :                END DO
     461              :             ELSE
     462           88 :                DO ispin = 1, nspins
     463           88 :                   CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp)
     464              :                END DO
     465           44 :                CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., my_recalc_integrals)
     466           88 :                DO ispin = 1, nspins
     467           88 :                   CALL dbcsr_scale(rho_ao(ispin)%matrix, -1.0_dp)
     468              :                END DO
     469              :             END IF
     470              :          END IF !do_hfx
     471              : 
     472              :       END IF ! do_exx
     473              : 
     474              :       ! Clean
     475          230 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     476          230 :       CALL auxbas_pw_pool%give_back_pw(vh_rspace)
     477          468 :       DO ispin = 1, nspins
     478          238 :          CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
     479          468 :          IF (ASSOCIATED(vtau_rspace)) THEN
     480           16 :             CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
     481              :          END IF
     482              :       END DO
     483          230 :       DEALLOCATE (vxc_rspace)
     484          230 :       IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
     485              : 
     486          230 :       IF (dft_control%do_admm) THEN
     487          164 :          DO ispin = 1, nspins
     488          164 :             IF (ASSOCIATED(vadmm_rspace)) THEN
     489           56 :                CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
     490              :             END IF
     491              :          END DO
     492           80 :          IF (ASSOCIATED(vadmm_rspace)) DEALLOCATE (vadmm_rspace)
     493              :       END IF
     494              : 
     495              :       !Add the HF contribution from RI_RPA/EC_ENV to the ks matrix
     496          230 :       CALL exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx)
     497              : 
     498          230 :       calc_ints = .TRUE.
     499          230 :       IF (reuse_hfx) calc_ints = .FALSE.
     500          230 :       IF (do_admm) THEN
     501           68 :          DO ispin = 1, nspins
     502           68 :             CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
     503              :          END DO
     504              : 
     505              :          CALL tddft_hfx_matrix(work_admm, rho_ao_aux, qs_env, .FALSE., calc_ints, ext_hfx_section, &
     506           32 :                                x_data)
     507              : 
     508              :          !ADMM XC correction
     509              :          CALL calc_exx_admm_xc_contributions(qs_env, dbcsr_work, work_admm, x_data, dummy_real1, &
     510           32 :                                              dummy_real2, .FALSE., .FALSE.)
     511              : 
     512           68 :          DO ispin = 1, nspins
     513           36 :             CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
     514              :             CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
     515           36 :                                0.0_dp, admm_env%work_aux_orb)
     516              :             CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
     517           36 :                                0.0_dp, admm_env%work_orb_orb)
     518           36 :             CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_tmp, keep_sparsity=.TRUE.)
     519           68 :             CALL dbcsr_add(dbcsr_work(ispin)%matrix, dbcsr_tmp, 1.0_dp, 1.0_dp)
     520              :          END DO
     521              : 
     522              :       ELSE
     523          198 :          CALL tddft_hfx_matrix(dbcsr_work, rho_ao, qs_env, .FALSE., calc_ints, ext_hfx_section, x_data)
     524              :       END IF
     525          230 :       CALL exx_post_hfx(qs_env, x_data, reuse_hfx)
     526              : 
     527              :       !Update the RHS
     528          468 :       DO ispin = 1, nspins
     529          468 :          CALL dbcsr_add(rhs(ispin)%matrix, dbcsr_work(ispin)%matrix, 1.0_dp, 1.0_dp)
     530              :       END DO
     531              : 
     532          230 :       CALL dbcsr_deallocate_matrix_set(dbcsr_work)
     533          230 :       IF (dft_control%do_admm) THEN
     534           80 :          CALL dbcsr_deallocate_matrix_set(work_admm)
     535           80 :          CALL dbcsr_release(dbcsr_tmp)
     536              :       END IF
     537              : 
     538          230 :       CALL timestop(handle)
     539              : 
     540          230 :    END SUBROUTINE add_exx_to_rhs
     541              : 
     542              : ! **************************************************************************************************
     543              : !> \brief get the ADMM XC section from the ri_rpa/ec_env type if available, create and store them otherwise
     544              : !> \param qs_env ...
     545              : !> \param x_data ...
     546              : !> \param xc_section_aux ...
     547              : !> \param xc_section_primary ...
     548              : ! **************************************************************************************************
     549           74 :    SUBROUTINE get_exx_admm_xc_sections(qs_env, x_data, xc_section_aux, xc_section_primary)
     550              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     551              :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     552              :       TYPE(section_vals_type), POINTER                   :: xc_section_aux, xc_section_primary
     553              : 
     554              :       INTEGER                                            :: natom
     555              :       TYPE(admm_type), POINTER                           :: qs_admm_env, tmp_admm_env
     556              :       TYPE(dft_control_type), POINTER                    :: dft_control
     557              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     558              :       TYPE(section_vals_type), POINTER                   :: xc_fun, xc_fun_empty, xc_section, &
     559              :                                                             xc_section_empty
     560              : 
     561           74 :       NULLIFY (qs_admm_env, tmp_admm_env, para_env, xc_section, xc_section_empty, xc_fun_empty, &
     562           74 :                xc_fun, dft_control)
     563              : 
     564           74 :       IF (ASSOCIATED(qs_env%mp2_env)) THEN
     565           26 :          IF (ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_aux) .AND. &
     566              :              ASSOCIATED(qs_env%mp2_env%ri_rpa%xc_section_primary)) THEN
     567           12 :             xc_section_aux => qs_env%mp2_env%ri_rpa%xc_section_aux
     568           12 :             xc_section_primary => qs_env%mp2_env%ri_rpa%xc_section_primary
     569              :          END IF
     570           48 :       ELSEIF (qs_env%energy_correction) THEN
     571           48 :          IF (ASSOCIATED(qs_env%ec_env%xc_section_aux) .AND. &
     572              :              ASSOCIATED(qs_env%ec_env%xc_section_primary)) THEN
     573           42 :             xc_section_aux => qs_env%ec_env%xc_section_aux
     574           42 :             xc_section_primary => qs_env%ec_env%xc_section_primary
     575              :          END IF
     576              :       END IF
     577              : 
     578           74 :       IF (.NOT. ASSOCIATED(xc_section_aux) .OR. .NOT. ASSOCIATED(xc_section_primary)) THEN
     579              : 
     580           20 :          CALL get_qs_env(qs_env, admm_env=qs_admm_env, natom=natom, para_env=para_env, dft_control=dft_control)
     581           20 :          CPASSERT(ASSOCIATED(qs_admm_env))
     582              : 
     583              :          !create XC section with XC_FUNCITONAL NONE (aka empty XC_FUNCTIONAL section)
     584           20 :          xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
     585           20 :          xc_fun => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     586           20 :          CALL section_vals_duplicate(xc_section, xc_section_empty)
     587           20 :          CALL section_vals_create(xc_fun_empty, xc_fun%section)
     588           20 :          CALL section_vals_set_subs_vals(xc_section_empty, "XC_FUNCTIONAL", xc_fun_empty)
     589              : 
     590              :          CALL admm_env_create(tmp_admm_env, dft_control%admm_control, qs_admm_env%mos_aux_fit, &
     591           20 :                               para_env, natom, qs_admm_env%nao_aux_fit)
     592              : 
     593              :          CALL create_admm_xc_section(x_data=x_data, xc_section=xc_section_empty, &
     594           20 :                                      admm_env=tmp_admm_env)
     595              : 
     596           20 :          CALL section_vals_duplicate(tmp_admm_env%xc_section_aux, xc_section_aux)
     597           20 :          CALL section_vals_duplicate(tmp_admm_env%xc_section_primary, xc_section_primary)
     598              : 
     599           20 :          IF (ASSOCIATED(qs_env%mp2_env)) THEN
     600           14 :             qs_env%mp2_env%ri_rpa%xc_section_aux => xc_section_aux
     601           14 :             qs_env%mp2_env%ri_rpa%xc_section_primary => xc_section_primary
     602              :          END IF
     603           20 :          IF (qs_env%energy_correction) THEN
     604            6 :             qs_env%ec_env%xc_section_aux => xc_section_aux
     605            6 :             qs_env%ec_env%xc_section_primary => xc_section_primary
     606              :          END IF
     607              : 
     608           20 :          CALL section_vals_release(xc_section_empty)
     609           20 :          CALL section_vals_release(xc_fun_empty)
     610           20 :          CALL admm_env_release(tmp_admm_env)
     611              : 
     612              :       END IF
     613              : 
     614           74 :    END SUBROUTINE get_exx_admm_xc_sections
     615              : 
     616              : ! **************************************************************************************************
     617              : !> \brief Calculate the RI_RPA%HF / EC_ENV%HF ADMM XC contributions to the KS matrices and the respective energies
     618              : !> \param qs_env ...
     619              : !> \param matrix_prim ...
     620              : !> \param matrix_aux ...
     621              : !> \param x_data ...
     622              : !> \param exc ...
     623              : !> \param exc_aux_fit ...
     624              : !> \param calc_forces ...
     625              : !> \param use_virial ...
     626              : ! **************************************************************************************************
     627          222 :    SUBROUTINE calc_exx_admm_xc_contributions(qs_env, matrix_prim, matrix_aux, x_data, exc, exc_aux_fit, &
     628              :                                              calc_forces, use_virial)
     629              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     630              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_prim, matrix_aux
     631              :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     632              :       REAL(dp), INTENT(INOUT)                            :: exc, exc_aux_fit
     633              :       LOGICAL, INTENT(IN)                                :: calc_forces, use_virial
     634              : 
     635              :       INTEGER                                            :: ispin, nspins
     636              :       REAL(dp), DIMENSION(3, 3)                          :: pv_loc
     637              :       TYPE(admm_type), POINTER                           :: admm_env
     638           74 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao, rho_ao_aux_fit
     639              :       TYPE(dft_control_type), POINTER                    :: dft_control
     640              :       TYPE(pw_env_type), POINTER                         :: pw_env
     641              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     642           74 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: v_dummy, v_rspace, v_rspace_aux_fit
     643              :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit
     644              :       TYPE(section_vals_type), POINTER                   :: xc_section_aux, xc_section_primary
     645              :       TYPE(virial_type), POINTER                         :: virial
     646              : 
     647           74 :       NULLIFY (xc_section_aux, xc_section_primary, rho, rho_aux_fit, v_dummy, v_rspace, v_rspace_aux_fit, &
     648           74 :                auxbas_pw_pool, pw_env, rho_ao, rho_ao_aux_fit, dft_control, admm_env)
     649              : 
     650           74 :       CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, rho=rho, admm_env=admm_env, virial=virial)
     651           74 :       CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
     652              : 
     653           74 :       nspins = dft_control%nspins
     654           74 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     655           74 :       CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
     656           74 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     657              : 
     658           74 :       CALL get_exx_admm_xc_sections(qs_env, x_data, xc_section_aux, xc_section_primary)
     659              : 
     660           74 :       IF (use_virial) virial%pv_xc = 0.0_dp
     661              :       CALL qs_vxc_create(qs_env%ks_env, rho_struct=rho_aux_fit, xc_section=xc_section_aux, &
     662           74 :                          vxc_rho=v_rspace_aux_fit, vxc_tau=v_dummy, exc=exc_aux_fit)
     663           74 :       IF (use_virial) THEN
     664            0 :          virial%pv_exc = virial%pv_exc - virial%pv_xc
     665            0 :          virial%pv_virial = virial%pv_virial - virial%pv_xc
     666              :       END IF
     667              : 
     668           74 :       IF (.NOT. dft_control%admm_control%aux_exch_func == do_admm_aux_exch_func_none) THEN
     669           74 :          IF (use_virial) pv_loc = virial%pv_virial
     670          158 :          DO ispin = 1, nspins
     671           84 :             CALL pw_scale(v_rspace_aux_fit(ispin), v_rspace_aux_fit(ispin)%pw_grid%dvol)
     672              :             CALL integrate_v_rspace(v_rspace=v_rspace_aux_fit(ispin), hmat=matrix_aux(ispin), &
     673              :                                     pmat=rho_ao_aux_fit(ispin), qs_env=qs_env, &
     674              :                                     basis_type="AUX_FIT", calculate_forces=calc_forces, &
     675          158 :                                     task_list_external=qs_env%admm_env%task_list_aux_fit)
     676              :          END DO
     677           74 :          IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
     678              :       END IF
     679              : 
     680           74 :       IF (ASSOCIATED(v_rspace_aux_fit)) THEN
     681          158 :          DO ispin = 1, nspins
     682          158 :             CALL auxbas_pw_pool%give_back_pw(v_rspace_aux_fit(ispin))
     683              :          END DO
     684           74 :          DEALLOCATE (v_rspace_aux_fit)
     685              :       END IF
     686           74 :       IF (ASSOCIATED(v_dummy)) THEN
     687            0 :          DO ispin = 1, nspins
     688            0 :             CALL auxbas_pw_pool%give_back_pw(v_dummy(ispin))
     689              :          END DO
     690            0 :          DEALLOCATE (v_dummy)
     691              :       END IF
     692              : 
     693           74 :       IF (use_virial) virial%pv_xc = 0.0_dp
     694              :       CALL qs_vxc_create(qs_env%ks_env, rho_struct=rho, xc_section=xc_section_primary, &
     695           74 :                          vxc_rho=v_rspace, vxc_tau=v_dummy, exc=exc)
     696           74 :       IF (use_virial) THEN
     697            0 :          virial%pv_exc = virial%pv_exc - virial%pv_xc
     698            0 :          virial%pv_virial = virial%pv_virial - virial%pv_xc
     699              :       END IF
     700              : 
     701           74 :       IF (.NOT. dft_control%admm_control%aux_exch_func == do_admm_aux_exch_func_none) THEN
     702           74 :          IF (use_virial) pv_loc = virial%pv_virial
     703          158 :          DO ispin = 1, nspins
     704           84 :             CALL pw_scale(v_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
     705              :             CALL integrate_v_rspace(v_rspace=v_rspace(ispin), hmat=matrix_prim(ispin), &
     706              :                                     pmat=rho_ao(ispin), qs_env=qs_env, &
     707          158 :                                     calculate_forces=calc_forces)
     708              :          END DO
     709           74 :          IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
     710              :       END IF
     711              : 
     712           74 :       IF (ASSOCIATED(v_rspace)) THEN
     713          158 :          DO ispin = 1, nspins
     714          158 :             CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
     715              :          END DO
     716           74 :          DEALLOCATE (v_rspace)
     717              :       END IF
     718           74 :       IF (ASSOCIATED(v_dummy)) THEN
     719            0 :          DO ispin = 1, nspins
     720            0 :             CALL auxbas_pw_pool%give_back_pw(v_dummy(ispin))
     721              :          END DO
     722            0 :          DEALLOCATE (v_dummy)
     723              :       END IF
     724              : 
     725           74 :    END SUBROUTINE calc_exx_admm_xc_contributions
     726              : 
     727              : ! **************************************************************************************************
     728              : !> \brief Prepare the external x_data for integration. Simply change the HFX fraction in case the
     729              : !>        qs_env%x_data is reused
     730              : !> \param ext_hfx_section ...
     731              : !> \param x_data ...
     732              : !> \param reuse_hfx ...
     733              : ! **************************************************************************************************
     734          620 :    SUBROUTINE exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx)
     735              :       TYPE(section_vals_type), POINTER                   :: ext_hfx_section
     736              :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     737              :       LOGICAL                                            :: reuse_hfx
     738              : 
     739              :       INTEGER                                            :: irep, n_rep_hf
     740              :       REAL(dp)                                           :: frac
     741              : 
     742          506 :       IF (.NOT. reuse_hfx) RETURN
     743              : 
     744          114 :       CALL section_vals_get(ext_hfx_section, n_repetition=n_rep_hf)
     745              : 
     746          228 :       DO irep = 1, n_rep_hf
     747          114 :          CALL section_vals_val_get(ext_hfx_section, "FRACTION", r_val=frac, i_rep_section=irep)
     748          342 :          x_data(irep, :)%general_parameter%fraction = frac
     749              :       END DO
     750              : 
     751              :    END SUBROUTINE exx_pre_hfx
     752              : 
     753              : ! **************************************************************************************************
     754              : !> \brief Revert back to the proper HFX fraction in case qs_env%x_data is reused
     755              : !> \param qs_env ...
     756              : !> \param x_data ...
     757              : !> \param reuse_hfx ...
     758              : ! **************************************************************************************************
     759          620 :    SUBROUTINE exx_post_hfx(qs_env, x_data, reuse_hfx)
     760              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     761              :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
     762              :       LOGICAL                                            :: reuse_hfx
     763              : 
     764              :       INTEGER                                            :: irep, n_rep_hf
     765              :       REAL(dp)                                           :: frac
     766              :       TYPE(section_vals_type), POINTER                   :: input, qs_hfx_section
     767              : 
     768          506 :       IF (.NOT. reuse_hfx) RETURN
     769              : 
     770          114 :       CALL get_qs_env(qs_env, input=input)
     771          114 :       qs_hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
     772          114 :       CALL section_vals_get(qs_hfx_section, n_repetition=n_rep_hf)
     773              : 
     774          228 :       DO irep = 1, n_rep_hf
     775          114 :          CALL section_vals_val_get(qs_hfx_section, "FRACTION", r_val=frac, i_rep_section=irep)
     776          342 :          x_data(irep, :)%general_parameter%fraction = frac
     777              :       END DO
     778              : 
     779              :    END SUBROUTINE exx_post_hfx
     780              : 
     781              : ! **************************************************************************************************
     782              : !> \brief ...
     783              : !> \param energy ...
     784              : ! **************************************************************************************************
     785          190 :    ELEMENTAL SUBROUTINE remove_exc_energy(energy)
     786              :       TYPE(qs_energy_type), INTENT(INOUT)                :: energy
     787              : 
     788              :       ! Remove the Exchange-correlation energy contributions from the total energy
     789              :       energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
     790          190 :                                      energy%exc_aux_fit + energy%exc1_aux_fit)
     791              : 
     792          190 :       energy%exc = 0.0_dp
     793          190 :       energy%exc1 = 0.0_dp
     794          190 :       energy%exc_aux_fit = 0.0_dp
     795          190 :       energy%exc1_aux_fit = 0.0_dp
     796          190 :       energy%ex = 0.0_dp
     797              : 
     798          190 :    END SUBROUTINE
     799              : 
     800              : END MODULE hfx_exx
     801              : 
        

Generated by: LCOV version 2.0-1