LCOV - code coverage report
Current view: top level - src - qs_vxc.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 88.7 % 293 260
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            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
      10              : !>
      11              : !>
      12              : !> \par History
      13              : !>     refactoring 03-2011 [MI]
      14              : !> \author MI
      15              : ! **************************************************************************************************
      16              : MODULE qs_vxc
      17              : 
      18              :    USE cell_types,                      ONLY: cell_type
      19              :    USE cp_control_types,                ONLY: dft_control_type
      20              :    USE input_constants,                 ONLY: sic_ad,&
      21              :                                               sic_eo,&
      22              :                                               sic_mauri_spz,&
      23              :                                               sic_mauri_us,&
      24              :                                               sic_none,&
      25              :                                               xc_none,&
      26              :                                               xc_vdw_fun_nonloc
      27              :    USE input_section_types,             ONLY: section_vals_type,&
      28              :                                               section_vals_val_get
      29              :    USE kinds,                           ONLY: dp
      30              :    USE message_passing,                 ONLY: mp_para_env_type
      31              :    USE pw_env_types,                    ONLY: pw_env_get,&
      32              :                                               pw_env_type
      33              :    USE pw_grids,                        ONLY: pw_grid_compare
      34              :    USE pw_methods,                      ONLY: pw_axpy,&
      35              :                                               pw_copy,&
      36              :                                               pw_multiply,&
      37              :                                               pw_scale,&
      38              :                                               pw_transfer,&
      39              :                                               pw_zero
      40              :    USE pw_pool_types,                   ONLY: pw_pool_type
      41              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      42              :                                               pw_r3d_rs_type
      43              :    USE qs_dispersion_nonloc,            ONLY: calculate_dispersion_nonloc
      44              :    USE qs_dispersion_types,             ONLY: qs_dispersion_type
      45              :    USE qs_ks_types,                     ONLY: get_ks_env,&
      46              :                                               qs_ks_env_type
      47              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      48              :                                               qs_rho_type
      49              :    USE virial_types,                    ONLY: virial_type
      50              :    USE xc,                              ONLY: calc_xc_density,&
      51              :                                               xc_exc_calc,&
      52              :                                               xc_vxc_pw_create1
      53              : #include "./base/base_uses.f90"
      54              : 
      55              :    IMPLICIT NONE
      56              : 
      57              :    PRIVATE
      58              : 
      59              :    ! *** Public subroutines ***
      60              :    PUBLIC :: qs_vxc_create, qs_xc_density
      61              : 
      62              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc'
      63              : 
      64              : CONTAINS
      65              : 
      66              : ! **************************************************************************************************
      67              : !> \brief calculates and allocates the xc potential, already reducing it to
      68              : !>      the dependence on rho and the one on tau
      69              : !> \param ks_env to get all the needed things
      70              : !> \param rho_struct density for which v_xc is calculated
      71              : !> \param xc_section ...
      72              : !> \param vxc_rho will contain the v_xc part that depend on rho
      73              : !>        (if one of the chosen xc functionals has it it is allocated and you
      74              : !>        are responsible for it)
      75              : !> \param vxc_tau will contain the kinetic tau part of v_xc
      76              : !>        (if one of the chosen xc functionals has it it is allocated and you
      77              : !>        are responsible for it)
      78              : !> \param exc ...
      79              : !> \param just_energy if true calculates just the energy, and does not
      80              : !>        allocate v_*_rspace
      81              : !> \param edisp ...
      82              : !> \param dispersion_env ...
      83              : !> \param adiabatic_rescale_factor ...
      84              : !> \param pw_env_external    external plane wave environment
      85              : !> \par History
      86              : !>      - 05.2002 modified to use the mp_allgather function each pe
      87              : !>        computes only part of the grid and this is broadcasted to all
      88              : !>        instead of summed.
      89              : !>        This scales significantly better (e.g. factor 3 on 12 cpus
      90              : !>        32 H2O) [Joost VdV]
      91              : !>      - moved to qs_ks_methods [fawzi]
      92              : !>      - sic alterations [Joost VandeVondele]
      93              : !> \author Fawzi Mohamed
      94              : ! **************************************************************************************************
      95       526636 :    SUBROUTINE qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, &
      96              :                             just_energy, edisp, dispersion_env, adiabatic_rescale_factor, &
      97              :                             pw_env_external)
      98              : 
      99              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     100              :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     101              :       TYPE(section_vals_type), POINTER                   :: xc_section
     102              :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: vxc_rho, vxc_tau
     103              :       REAL(KIND=dp), INTENT(out)                         :: exc
     104              :       LOGICAL, INTENT(in), OPTIONAL                      :: just_energy
     105              :       REAL(KIND=dp), INTENT(out), OPTIONAL               :: edisp
     106              :       TYPE(qs_dispersion_type), OPTIONAL, POINTER        :: dispersion_env
     107              :       REAL(KIND=dp), INTENT(in), OPTIONAL                :: adiabatic_rescale_factor
     108              :       TYPE(pw_env_type), OPTIONAL, POINTER               :: pw_env_external
     109              : 
     110              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_vxc_create'
     111              : 
     112              :       INTEGER                                            :: handle, ispin, mspin, myfun, &
     113              :                                                             nelec_spin(2), vdw
     114              :       LOGICAL :: compute_virial, do_adiabatic_rescaling, my_just_energy, rho_g_valid, &
     115              :          sic_scaling_b_zero, tau_r_valid, uf_grid, vdW_nl
     116              :       REAL(KIND=dp)                                      :: exc_m, factor, &
     117              :                                                             my_adiabatic_rescale_factor, &
     118              :                                                             my_scaling, nelec_s_inv
     119              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: virial_xc_tmp
     120              :       TYPE(cell_type), POINTER                           :: cell
     121              :       TYPE(dft_control_type), POINTER                    :: dft_control
     122              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     123       131659 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g, rho_m_gspace, rho_struct_g
     124              :       TYPE(pw_c1d_gs_type), POINTER                      :: rho_nlcc_g, tmp_g, tmp_g2
     125              :       TYPE(pw_env_type), POINTER                         :: pw_env
     126              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
     127       131659 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: my_vxc_rho, my_vxc_tau, rho_m_rspace, &
     128       131659 :                                                             rho_r, rho_struct_r, tau, tau_struct_r
     129              :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_nlcc, tmp_pw
     130              :       TYPE(virial_type), POINTER                         :: virial
     131              : 
     132       131659 :       CALL timeset(routineN, handle)
     133              : 
     134       131659 :       CPASSERT(.NOT. ASSOCIATED(vxc_rho))
     135       131659 :       CPASSERT(.NOT. ASSOCIATED(vxc_tau))
     136       131659 :       NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, my_vxc_rho, &
     137       131659 :                tmp_pw, tmp_g, tmp_g2, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, &
     138       131659 :                rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_struct_r, rho_struct_g, tau_struct_r)
     139              : 
     140       131659 :       exc = 0.0_dp
     141       131659 :       my_just_energy = .FALSE.
     142       131659 :       IF (PRESENT(just_energy)) my_just_energy = just_energy
     143       131659 :       my_adiabatic_rescale_factor = 1.0_dp
     144       131659 :       do_adiabatic_rescaling = .FALSE.
     145       131659 :       IF (PRESENT(adiabatic_rescale_factor)) THEN
     146           44 :          my_adiabatic_rescale_factor = adiabatic_rescale_factor
     147           44 :          do_adiabatic_rescaling = .TRUE.
     148              :       END IF
     149              : 
     150              :       CALL get_ks_env(ks_env, &
     151              :                       dft_control=dft_control, &
     152              :                       pw_env=pw_env, &
     153              :                       cell=cell, &
     154              :                       virial=virial, &
     155              :                       rho_nlcc=rho_nlcc, &
     156       131659 :                       rho_nlcc_g=rho_nlcc_g)
     157              : 
     158              :       CALL qs_rho_get(rho_struct, &
     159              :                       tau_r_valid=tau_r_valid, &
     160              :                       rho_g_valid=rho_g_valid, &
     161              :                       rho_r=rho_struct_r, &
     162              :                       rho_g=rho_struct_g, &
     163       131659 :                       tau_r=tau_struct_r)
     164              : 
     165       131659 :       compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
     166       131659 :       IF (compute_virial) THEN
     167        35958 :          virial%pv_xc = 0.0_dp
     168              :       END IF
     169              : 
     170              :       CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
     171       131659 :                                 i_val=myfun)
     172              :       CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
     173       131659 :                                 i_val=vdw)
     174              : 
     175       131659 :       vdW_nl = (vdw == xc_vdw_fun_nonloc)
     176              :       ! this combination has not been investigated
     177       131659 :       CPASSERT(.NOT. (do_adiabatic_rescaling .AND. vdW_nl))
     178              :       ! are the necessary inputs available
     179       131659 :       IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
     180              :          vdW_nl = .FALSE.
     181              :       END IF
     182       131659 :       IF (PRESENT(edisp)) edisp = 0.0_dp
     183              : 
     184       131659 :       IF (myfun /= xc_none .OR. vdW_nl) THEN
     185              : 
     186              :          ! test if the real space density is available
     187       118647 :          CPASSERT(ASSOCIATED(rho_struct))
     188       118647 :          IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) &
     189            0 :             CPABORT("nspins must be 1 or 2")
     190       118647 :          mspin = SIZE(rho_struct_r)
     191       118647 :          IF (dft_control%nspins == 2 .AND. mspin == 1) &
     192            0 :             CPABORT("Spin count mismatch")
     193              : 
     194              :          ! there are some options related to SIC here.
     195              :          ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD)
     196              :          ! SIC can             E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta))
     197              :          ! or compute          E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0)
     198              : 
     199              :          ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term
     200       118647 :          my_scaling = 1.0_dp
     201       118851 :          SELECT CASE (dft_control%sic_method_id)
     202              :          CASE (sic_none)
     203              :             ! all fine
     204              :          CASE (sic_mauri_spz, sic_ad)
     205              :             ! no idea yet what to do here in that case
     206          204 :             CPASSERT(.NOT. tau_r_valid)
     207              :          CASE (sic_mauri_us)
     208           92 :             my_scaling = 1.0_dp - dft_control%sic_scaling_b
     209              :             ! no idea yet what to do here in that case
     210           92 :             CPASSERT(.NOT. tau_r_valid)
     211              :          CASE (sic_eo)
     212              :             ! NOTHING TO BE DONE
     213              :          CASE DEFAULT
     214              :             ! this case has not yet been treated here
     215       118647 :             CPABORT("NYI")
     216              :          END SELECT
     217              : 
     218       118647 :          IF (dft_control%sic_scaling_b .EQ. 0.0_dp) THEN
     219              :             sic_scaling_b_zero = .TRUE.
     220              :          ELSE
     221       118547 :             sic_scaling_b_zero = .FALSE.
     222              :          END IF
     223              : 
     224       118647 :          IF (PRESENT(pw_env_external)) &
     225            0 :             pw_env => pw_env_external
     226       118647 :          CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
     227       118647 :          uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
     228              : 
     229       118647 :          IF (.NOT. uf_grid) THEN
     230       118633 :             rho_r => rho_struct_r
     231              : 
     232       118633 :             IF (tau_r_valid) THEN
     233         2734 :                tau => tau_struct_r
     234              :             END IF
     235              : 
     236              :             ! for gradient corrected functional the density in g space might
     237              :             ! be useful so if we have it, we pass it in
     238       118633 :             IF (rho_g_valid) THEN
     239       118613 :                rho_g => rho_struct_g
     240              :             END IF
     241              :          ELSE
     242           14 :             CPASSERT(rho_g_valid)
     243           56 :             ALLOCATE (rho_r(mspin))
     244           56 :             ALLOCATE (rho_g(mspin))
     245           28 :             DO ispin = 1, mspin
     246           14 :                CALL xc_pw_pool%create_pw(rho_g(ispin))
     247           28 :                CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
     248              :             END DO
     249           28 :             DO ispin = 1, mspin
     250           14 :                CALL xc_pw_pool%create_pw(rho_r(ispin))
     251           28 :                CALL pw_transfer(rho_g(ispin), rho_r(ispin))
     252              :             END DO
     253           14 :             IF (tau_r_valid) THEN
     254              :                ! tau with finer grids is not implemented (at least not correctly), which this asserts
     255            0 :                CPABORT("tau with finer grids not implemented")
     256              :             END IF
     257              :          END IF
     258              : 
     259              :          ! add the nlcc densities
     260       118647 :          IF (ASSOCIATED(rho_nlcc)) THEN
     261          354 :             factor = 1.0_dp
     262          708 :             DO ispin = 1, mspin
     263          354 :                CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
     264          708 :                CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
     265              :             END DO
     266              :          END IF
     267              : 
     268              :          !
     269              :          ! here the rho_r, rho_g, tau is what it should be
     270              :          ! we get back the right my_vxc_rho and my_vxc_tau as required
     271              :          !
     272       118647 :          IF (my_just_energy) THEN
     273              :             exc = xc_exc_calc(rho_r=rho_r, tau=tau, &
     274              :                               rho_g=rho_g, xc_section=xc_section, &
     275        10362 :                               pw_pool=xc_pw_pool)
     276              : 
     277              :          ELSE
     278              :             CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
     279              :                                    rho_g=rho_g, tau=tau, exc=exc, &
     280              :                                    xc_section=xc_section, &
     281              :                                    pw_pool=xc_pw_pool, &
     282              :                                    compute_virial=compute_virial, &
     283       108285 :                                    virial_xc=virial%pv_xc)
     284              :          END IF
     285              : 
     286              :          ! remove the nlcc densities (keep stuff in original state)
     287       118647 :          IF (ASSOCIATED(rho_nlcc)) THEN
     288          354 :             factor = -1.0_dp
     289          708 :             DO ispin = 1, mspin
     290          354 :                CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
     291          708 :                CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
     292              :             END DO
     293              :          END IF
     294              : 
     295              :          ! calclulate non-local vdW functional
     296              :          ! only if this XC_SECTION has it
     297              :          ! if yes, we use the dispersion_env from ks_env
     298              :          ! this is dangerous, as it assumes a special connection xc_section -> qs_env
     299       118647 :          IF (vdW_nl) THEN
     300          388 :             CALL get_ks_env(ks_env=ks_env, para_env=para_env)
     301              :             ! no SIC functionals allowed
     302          388 :             CPASSERT(dft_control%sic_method_id == sic_none)
     303              :             !
     304          388 :             CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
     305          388 :             IF (my_just_energy) THEN
     306              :                CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
     307            6 :                                                 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env)
     308              :             ELSE
     309              :                CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
     310          382 :                                                 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env, virial=virial)
     311              :             END IF
     312              :          END IF
     313              : 
     314              :          !! Apply rescaling to the potential if requested
     315       118647 :          IF (.NOT. my_just_energy) THEN
     316       108285 :             IF (do_adiabatic_rescaling) THEN
     317           24 :                IF (ASSOCIATED(my_vxc_rho)) THEN
     318           62 :                   DO ispin = 1, SIZE(my_vxc_rho)
     319           62 :                      CALL pw_scale(my_vxc_rho(ispin), my_adiabatic_rescale_factor)
     320              :                   END DO
     321              :                END IF
     322              :             END IF
     323              :          END IF
     324              : 
     325       118647 :          IF (my_scaling .NE. 1.0_dp) THEN
     326           92 :             exc = exc*my_scaling
     327           92 :             IF (ASSOCIATED(my_vxc_rho)) THEN
     328          180 :                DO ispin = 1, SIZE(my_vxc_rho)
     329          180 :                   CALL pw_scale(my_vxc_rho(ispin), my_scaling)
     330              :                END DO
     331              :             END IF
     332           92 :             IF (ASSOCIATED(my_vxc_tau)) THEN
     333            0 :                DO ispin = 1, SIZE(my_vxc_tau)
     334            0 :                   CALL pw_scale(my_vxc_tau(ispin), my_scaling)
     335              :                END DO
     336              :             END IF
     337              :          END IF
     338              : 
     339              :          ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer
     340              :          ! pw -> coeff
     341       118647 :          IF (ASSOCIATED(my_vxc_rho)) THEN
     342       108285 :             vxc_rho => my_vxc_rho
     343       108285 :             NULLIFY (my_vxc_rho)
     344              :          END IF
     345       118647 :          IF (ASSOCIATED(my_vxc_tau)) THEN
     346         2008 :             vxc_tau => my_vxc_tau
     347         2008 :             NULLIFY (my_vxc_tau)
     348              :          END IF
     349       118647 :          IF (uf_grid) THEN
     350           28 :             DO ispin = 1, SIZE(rho_r)
     351           28 :                CALL xc_pw_pool%give_back_pw(rho_r(ispin))
     352              :             END DO
     353           14 :             DEALLOCATE (rho_r)
     354           14 :             IF (ASSOCIATED(rho_g)) THEN
     355           28 :                DO ispin = 1, SIZE(rho_g)
     356           28 :                   CALL xc_pw_pool%give_back_pw(rho_g(ispin))
     357              :                END DO
     358           14 :                DEALLOCATE (rho_g)
     359              :             END IF
     360              :          END IF
     361              : 
     362              :          ! compute again the xc but now for Exc(m,o) and the opposite sign
     363       118647 :          IF (dft_control%sic_method_id .EQ. sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
     364          390 :             ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
     365           78 :             CALL xc_pw_pool%create_pw(rho_m_gspace(1))
     366           78 :             CALL xc_pw_pool%create_pw(rho_m_rspace(1))
     367           78 :             CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
     368           78 :             CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
     369           78 :             CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
     370           78 :             CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
     371              :             ! bit sad, these will be just zero...
     372           78 :             CALL xc_pw_pool%create_pw(rho_m_gspace(2))
     373           78 :             CALL xc_pw_pool%create_pw(rho_m_rspace(2))
     374           78 :             CALL pw_zero(rho_m_rspace(2))
     375           78 :             CALL pw_zero(rho_m_gspace(2))
     376              : 
     377           78 :             IF (my_just_energy) THEN
     378              :                exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
     379              :                                    rho_g=rho_m_gspace, xc_section=xc_section, &
     380           24 :                                    pw_pool=xc_pw_pool)
     381              :             ELSE
     382              :                ! virial untested
     383           54 :                CPASSERT(.NOT. compute_virial)
     384              :                CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
     385              :                                       rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
     386              :                                       xc_section=xc_section, &
     387              :                                       pw_pool=xc_pw_pool, &
     388              :                                       compute_virial=.FALSE., &
     389           54 :                                       virial_xc=virial_xc_tmp)
     390              :             END IF
     391              : 
     392           78 :             exc = exc - dft_control%sic_scaling_b*exc_m
     393              : 
     394              :             ! and take care of the potential only vxc_rho is taken into account
     395           78 :             IF (.NOT. my_just_energy) THEN
     396           54 :                CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
     397           54 :                CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
     398           54 :                CALL my_vxc_rho(1)%release()
     399           54 :                CALL my_vxc_rho(2)%release()
     400           54 :                DEALLOCATE (my_vxc_rho)
     401              :             END IF
     402              : 
     403          234 :             DO ispin = 1, 2
     404          156 :                CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
     405          234 :                CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
     406              :             END DO
     407           78 :             DEALLOCATE (rho_m_rspace)
     408           78 :             DEALLOCATE (rho_m_gspace)
     409              : 
     410              :          END IF
     411              : 
     412              :          ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
     413       118647 :          IF (dft_control%sic_method_id .EQ. sic_ad .AND. .NOT. sic_scaling_b_zero) THEN
     414              : 
     415              :             ! find out how many elecs we have
     416           26 :             CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
     417              : 
     418          130 :             ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
     419           78 :             DO ispin = 1, 2
     420           52 :                CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
     421           78 :                CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
     422              :             END DO
     423              : 
     424           78 :             DO ispin = 1, 2
     425           52 :                IF (nelec_spin(ispin) .GT. 0.0_dp) THEN
     426           52 :                   nelec_s_inv = 1.0_dp/nelec_spin(ispin)
     427              :                ELSE
     428              :                   ! does it matter if there are no electrons with this spin (H) ?
     429            0 :                   nelec_s_inv = 0.0_dp
     430              :                END IF
     431           52 :                CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
     432           52 :                CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
     433           52 :                CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
     434           52 :                CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
     435           52 :                CALL pw_zero(rho_m_rspace(2))
     436           52 :                CALL pw_zero(rho_m_gspace(2))
     437              : 
     438           52 :                IF (my_just_energy) THEN
     439              :                   exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
     440              :                                       rho_g=rho_m_gspace, xc_section=xc_section, &
     441           12 :                                       pw_pool=xc_pw_pool)
     442              :                ELSE
     443              :                   ! virial untested
     444           40 :                   CPASSERT(.NOT. compute_virial)
     445              :                   CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
     446              :                                          rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
     447              :                                          xc_section=xc_section, &
     448              :                                          pw_pool=xc_pw_pool, &
     449              :                                          compute_virial=.FALSE., &
     450           40 :                                          virial_xc=virial_xc_tmp)
     451              :                END IF
     452              : 
     453           52 :                exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m
     454              : 
     455              :                ! and take care of the potential only vxc_rho is taken into account
     456           78 :                IF (.NOT. my_just_energy) THEN
     457           40 :                   CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
     458           40 :                   CALL my_vxc_rho(1)%release()
     459           40 :                   CALL my_vxc_rho(2)%release()
     460           40 :                   DEALLOCATE (my_vxc_rho)
     461              :                END IF
     462              :             END DO
     463              : 
     464           78 :             DO ispin = 1, 2
     465           52 :                CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
     466           78 :                CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
     467              :             END DO
     468           26 :             DEALLOCATE (rho_m_rspace)
     469           26 :             DEALLOCATE (rho_m_gspace)
     470              : 
     471              :          END IF
     472              : 
     473              :          ! compute again the xc but now for Exc(n_down,n_down)
     474       118647 :          IF (dft_control%sic_method_id .EQ. sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
     475          276 :             ALLOCATE (rho_r(2))
     476           92 :             rho_r(1) = rho_struct_r(2)
     477           92 :             rho_r(2) = rho_struct_r(2)
     478           92 :             IF (rho_g_valid) THEN
     479          276 :                ALLOCATE (rho_g(2))
     480           92 :                rho_g(1) = rho_struct_g(2)
     481           92 :                rho_g(2) = rho_struct_g(2)
     482              :             END IF
     483              : 
     484           92 :             IF (my_just_energy) THEN
     485              :                exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, &
     486              :                                    rho_g=rho_g, xc_section=xc_section, &
     487           32 :                                    pw_pool=xc_pw_pool)
     488              :             ELSE
     489              :                ! virial untested
     490           60 :                CPASSERT(.NOT. compute_virial)
     491              :                CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
     492              :                                       rho_g=rho_g, tau=tau, exc=exc_m, &
     493              :                                       xc_section=xc_section, &
     494              :                                       pw_pool=xc_pw_pool, &
     495              :                                       compute_virial=.FALSE., &
     496           60 :                                       virial_xc=virial_xc_tmp)
     497              :             END IF
     498              : 
     499           92 :             exc = exc + dft_control%sic_scaling_b*exc_m
     500              : 
     501              :             ! and take care of the potential
     502           92 :             IF (.NOT. my_just_energy) THEN
     503              :                ! both go to minority spin
     504           60 :                CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
     505           60 :                CALL my_vxc_rho(1)%release()
     506           60 :                CALL my_vxc_rho(2)%release()
     507           60 :                DEALLOCATE (my_vxc_rho)
     508              :             END IF
     509           92 :             DEALLOCATE (rho_r, rho_g)
     510              : 
     511              :          END IF
     512              : 
     513              :          !
     514              :          ! cleanups
     515              :          !
     516       118647 :          IF (uf_grid .AND. (ASSOCIATED(vxc_rho) .OR. ASSOCIATED(vxc_tau))) THEN
     517              :             BLOCK
     518              :                TYPE(pw_r3d_rs_type) :: tmp_pw
     519              :                TYPE(pw_c1d_gs_type) :: tmp_g, tmp_g2
     520           14 :                CALL xc_pw_pool%create_pw(tmp_g)
     521           14 :                CALL auxbas_pw_pool%create_pw(tmp_g2)
     522           14 :                IF (ASSOCIATED(vxc_rho)) THEN
     523           28 :                   DO ispin = 1, SIZE(vxc_rho)
     524           14 :                      CALL auxbas_pw_pool%create_pw(tmp_pw)
     525           14 :                      CALL pw_transfer(vxc_rho(ispin), tmp_g)
     526           14 :                      CALL pw_transfer(tmp_g, tmp_g2)
     527           14 :                      CALL pw_transfer(tmp_g2, tmp_pw)
     528           14 :                      CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
     529           28 :                      vxc_rho(ispin) = tmp_pw
     530              :                   END DO
     531              :                END IF
     532           14 :                IF (ASSOCIATED(vxc_tau)) THEN
     533            0 :                   DO ispin = 1, SIZE(vxc_tau)
     534            0 :                      CALL auxbas_pw_pool%create_pw(tmp_pw)
     535            0 :                      CALL pw_transfer(vxc_tau(ispin), tmp_g)
     536            0 :                      CALL pw_transfer(tmp_g, tmp_g2)
     537            0 :                      CALL pw_transfer(tmp_g2, tmp_pw)
     538            0 :                      CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
     539            0 :                      vxc_tau(ispin) = tmp_pw
     540              :                   END DO
     541              :                END IF
     542           14 :                CALL auxbas_pw_pool%give_back_pw(tmp_g2)
     543           28 :                CALL xc_pw_pool%give_back_pw(tmp_g)
     544              :             END BLOCK
     545              :          END IF
     546       118647 :          IF (ASSOCIATED(tau) .AND. uf_grid) THEN
     547            0 :             DO ispin = 1, SIZE(tau)
     548            0 :                CALL xc_pw_pool%give_back_pw(tau(ispin))
     549              :             END DO
     550            0 :             DEALLOCATE (tau)
     551              :          END IF
     552              : 
     553              :       END IF
     554              : 
     555       131659 :       CALL timestop(handle)
     556              : 
     557       131659 :    END SUBROUTINE qs_vxc_create
     558              : 
     559              : ! **************************************************************************************************
     560              : !> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r)  or  E_xc(r)/rho(r)
     561              : !> \param ks_env to get all the needed things
     562              : !> \param rho_struct density
     563              : !> \param xc_section ...
     564              : !> \param dispersion_env ...
     565              : !> \param xc_ener will contain the xc energy density E_xc(r) - V_xc(r)*rho(r)
     566              : !> \param xc_den will contain the xc energy density E_xc(r)/rho(r)
     567              : !> \param vxc ...
     568              : !> \param vtau ...
     569              : !> \author JGH
     570              : ! **************************************************************************************************
     571          490 :    SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
     572           98 :                             xc_ener, xc_den, vxc, vtau)
     573              : 
     574              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     575              :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     576              :       TYPE(section_vals_type), POINTER                   :: xc_section
     577              :       TYPE(qs_dispersion_type), OPTIONAL, POINTER        :: dispersion_env
     578              :       TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL      :: xc_ener, xc_den
     579              :       TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL       :: vxc, vtau
     580              : 
     581              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_xc_density'
     582              : 
     583              :       INTEGER                                            :: handle, ispin, myfun, nspins, vdw
     584              :       LOGICAL                                            :: rho_g_valid, tau_r_valid, uf_grid, vdW_nl
     585              :       REAL(KIND=dp)                                      :: edisp, exc, factor, rho_cutoff
     586              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: vdum
     587              :       TYPE(cell_type), POINTER                           :: cell
     588              :       TYPE(dft_control_type), POINTER                    :: dft_control
     589              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     590           98 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g
     591              :       TYPE(pw_c1d_gs_type), POINTER                      :: rho_nlcc_g
     592              :       TYPE(pw_env_type), POINTER                         :: pw_env
     593              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
     594              :       TYPE(pw_r3d_rs_type)                               :: exc_r
     595           98 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r, tau_r, vxc_rho, vxc_tau
     596              :       TYPE(pw_r3d_rs_type), POINTER                      :: rho_nlcc
     597              : 
     598           98 :       CALL timeset(routineN, handle)
     599              : 
     600              :       CALL get_ks_env(ks_env, &
     601              :                       dft_control=dft_control, &
     602              :                       pw_env=pw_env, &
     603              :                       cell=cell, &
     604              :                       rho_nlcc=rho_nlcc, &
     605           98 :                       rho_nlcc_g=rho_nlcc_g)
     606              : 
     607              :       CALL qs_rho_get(rho_struct, &
     608              :                       tau_r_valid=tau_r_valid, &
     609              :                       rho_g_valid=rho_g_valid, &
     610              :                       rho_r=rho_r, &
     611              :                       rho_g=rho_g, &
     612           98 :                       tau_r=tau_r)
     613           98 :       nspins = dft_control%nspins
     614              : 
     615           98 :       CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
     616           98 :       CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
     617           98 :       vdW_nl = (vdw == xc_vdw_fun_nonloc)
     618           98 :       IF (PRESENT(xc_ener)) THEN
     619           32 :          IF (tau_r_valid) THEN
     620            0 :             CALL cp_warn(__LOCATION__, "Tau contribution will not be correctly handled")
     621              :          END IF
     622              :       END IF
     623           98 :       IF (vdW_nl) THEN
     624            0 :          CALL cp_warn(__LOCATION__, "vdW functional contribution will be ignored")
     625              :       END IF
     626              : 
     627           98 :       CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
     628           98 :       uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
     629           98 :       IF (uf_grid) THEN
     630            0 :          CALL cp_warn(__LOCATION__, "Fine grid option not possible with local energy")
     631            0 :          CPABORT("Fine Grid in Local Energy")
     632              :       END IF
     633              : 
     634           98 :       IF (PRESENT(xc_ener)) THEN
     635           32 :          CALL pw_zero(xc_ener)
     636              :       END IF
     637           98 :       IF (PRESENT(xc_den)) THEN
     638           66 :          CALL pw_zero(xc_den)
     639              :       END IF
     640           98 :       IF (PRESENT(vxc)) THEN
     641          138 :          DO ispin = 1, nspins
     642          138 :             CALL pw_zero(vxc(ispin))
     643              :          END DO
     644              :       END IF
     645           98 :       IF (PRESENT(vtau)) THEN
     646           40 :          DO ispin = 1, nspins
     647           40 :             CALL pw_zero(vtau(ispin))
     648              :          END DO
     649              :       END IF
     650              : 
     651           98 :       IF (myfun /= xc_none) THEN
     652              : 
     653           96 :          CPASSERT(ASSOCIATED(rho_struct))
     654           96 :          CPASSERT(dft_control%sic_method_id == sic_none)
     655              : 
     656              :          ! add the nlcc densities
     657           96 :          IF (ASSOCIATED(rho_nlcc)) THEN
     658            0 :             factor = 1.0_dp
     659            0 :             DO ispin = 1, nspins
     660            0 :                CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
     661            0 :                CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
     662              :             END DO
     663              :          END IF
     664           96 :          NULLIFY (vxc_rho, vxc_tau)
     665              :          CALL xc_vxc_pw_create1(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
     666              :                                 rho_g=rho_g, tau=tau_r, exc=exc, &
     667              :                                 xc_section=xc_section, &
     668              :                                 pw_pool=xc_pw_pool, &
     669              :                                 compute_virial=.FALSE., &
     670              :                                 virial_xc=vdum, &
     671           96 :                                 exc_r=exc_r)
     672              :          ! calclulate non-local vdW functional
     673              :          ! only if this XC_SECTION has it
     674              :          ! if yes, we use the dispersion_env from ks_env
     675              :          ! this is dangerous, as it assumes a special connection xc_section -> qs_env
     676           96 :          IF (vdW_nl) THEN
     677            0 :             CALL get_ks_env(ks_env=ks_env, para_env=para_env)
     678              :             ! no SIC functionals allowed
     679            0 :             CPASSERT(dft_control%sic_method_id == sic_none)
     680              :             !
     681            0 :             CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
     682              :             CALL calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
     683            0 :                                              .FALSE., vdw_pw_pool, xc_pw_pool, para_env)
     684              :          END IF
     685              : 
     686              :          ! remove the nlcc densities (keep stuff in original state)
     687           96 :          IF (ASSOCIATED(rho_nlcc)) THEN
     688            0 :             factor = -1.0_dp
     689            0 :             DO ispin = 1, dft_control%nspins
     690            0 :                CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
     691            0 :                CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
     692              :             END DO
     693              :          END IF
     694              :          !
     695           96 :          IF (PRESENT(xc_den)) THEN
     696           64 :             CALL pw_copy(exc_r, xc_den)
     697           64 :             rho_cutoff = 1.E-14_dp
     698           64 :             CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
     699              :          END IF
     700           96 :          IF (PRESENT(xc_ener)) THEN
     701           32 :             CALL pw_copy(exc_r, xc_ener)
     702           64 :             DO ispin = 1, nspins
     703           64 :                CALL pw_multiply(xc_ener, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
     704              :             END DO
     705              :          END IF
     706           96 :          IF (PRESENT(vxc)) THEN
     707          134 :             DO ispin = 1, nspins
     708          134 :                CALL pw_copy(vxc_rho(ispin), vxc(ispin))
     709              :             END DO
     710              :          END IF
     711           96 :          IF (PRESENT(vtau)) THEN
     712           40 :             DO ispin = 1, nspins
     713           40 :                CALL pw_copy(vxc_tau(ispin), vtau(ispin))
     714              :             END DO
     715              :          END IF
     716              :          ! remove arrays
     717           96 :          IF (ASSOCIATED(vxc_rho)) THEN
     718          198 :             DO ispin = 1, nspins
     719          198 :                CALL vxc_rho(ispin)%release()
     720              :             END DO
     721           96 :             DEALLOCATE (vxc_rho)
     722              :          END IF
     723           96 :          IF (ASSOCIATED(vxc_tau)) THEN
     724           40 :             DO ispin = 1, nspins
     725           40 :                CALL vxc_tau(ispin)%release()
     726              :             END DO
     727           20 :             DEALLOCATE (vxc_tau)
     728              :          END IF
     729           96 :          CALL exc_r%release()
     730              :          !
     731              :       END IF
     732              : 
     733           98 :       CALL timestop(handle)
     734              : 
     735           98 :    END SUBROUTINE qs_xc_density
     736              : 
     737              : END MODULE qs_vxc
        

Generated by: LCOV version 2.0-1