xc_adiabatic_methods.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Contains some functions used in the context of adiabatic hybrid functionals
!> \par History
!>      01.2008 created [Manuel Guidon]
!> \author Manuel Guidon
! **************************************************************************************************
MODULE xc_adiabatic_methods
 
   USE input_constants,                 ONLY: do_potential_coulomb
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: oorootpi
   USE qs_energy_types,                 ONLY: qs_energy_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_mo_types,                     ONLY: get_mo_set,&
                                              mo_set_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   PUBLIC :: rescale_mcy3_pade
 
CONTAINS
 
! **************************************************************************************************
!> \brief - Calculates rescaling factors for XC potentials and energy expression
!> \param qs_env ...
!> \param hf_energy Array of size 2 containing the two HF energies (Ex^{HF} and Ex^{HF,LR}
!> \param energy QS energy type
!> \param adiabatic_lambda , adiabatic_omega: Parameters for adiabatic connection
!> \param adiabatic_omega ...
!> \param scale_dEx1 scaling coefficient for xc-potentials to be calculated
!> \param scale_ddW0 scaling coefficient for xc-potentials to be calculated
!> \param scale_dDFA scaling coefficient for xc-potentials to be calculated
!> \param scale_dEx2 scaling coefficient for xc-potentials to be calculated
!> \param total_energy_xc will contain the full xc energy
!> \par History
!>      09.2007 created [Manuel Guidon]
!> \author Manuel Guidon
!> \note
!>      - Model for adiabatic connection:
!>
!>         W_lambda = a + b*lambda/(1+c*lambda)
!>         Exc = a + b*(c-log(1+c)/c^2)
!>         a = Ex^{HF}
!>         b = -c1*2*omega/sqrt(PI)*nelectron
!>         c = -1/lambda - b/(Ex^{HF}-W_lambda^{BLYP} + c2*W_lambda^{B88,LR}-c3*W_lambda^{HF,LR}
! **************************************************************************************************
   SUBROUTINE rescale_mcy3_pade(qs_env, hf_energy, energy, adiabatic_lambda, &
                                adiabatic_omega, scale_dEx1, scale_ddW0, scale_dDFA, &
                                scale_dEx2, total_energy_xc)
 
      TYPE(qs_environment_type), POINTER                 :: qs_env
      REAL(dp), INTENT(INOUT)                            :: hf_energy(*)
      TYPE(qs_energy_type), POINTER                      :: energy
      REAL(dp), INTENT(IN)                               :: adiabatic_lambda, adiabatic_omega
      REAL(dp), INTENT(INOUT)                            :: scale_dex1, scale_ddw0, scale_ddfa, &
                                                            scale_dex2, total_energy_xc
 
      INTEGER                                            :: nelec_a, nelec_b, nelectron, nspins
      LOGICAL                                            :: do_swap_hf
      REAL(dp) :: a, b, c, c1, da_ddfa, da_ddw0, da_dex1, da_dex2, db_ddfa, db_ddw0, db_dex1, &
         db_dex2, dc_ddfa, dc_ddw0, dc_dex1, dc_dex2, dexc_da, dexc_db, dexc_dc, dfa_energy, &
         swap_value
      TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
 
      do_swap_hf = .false.
      !! Assume the first HF section is the Coulomb one
      IF (qs_env%x_data(1, 1)%potential_parameter%potential_type /= do_potential_coulomb) do_swap_hf = .true.
 
      IF (do_swap_hf) THEN
         swap_value = hf_energy(1)
         hf_energy(1) = hf_energy(2)
         hf_energy(2) = swap_value
      END IF
 
      c1 = 0.23163_dp
 
      CALL get_qs_env(qs_env=qs_env, mos=mos)
      CALL get_mo_set(mo_set=mos(1), nelectron=nelec_a)
      nspins = SIZE(mos)
      IF (nspins == 2) THEN
         CALL get_mo_set(mo_set=mos(2), nelectron=nelec_b)
      ELSE
         nelec_b = 0
      END IF
      nelectron = nelec_a + nelec_b
      dfa_energy = energy%exc + energy%exc1
      a = hf_energy(1)
      b = -c1*2.0_dp*adiabatic_omega*oorootpi*nelectron !-0.23163_dp*2.0_dp*0.2_dp*oorootpi*nelectron
      c = -1.0_dp/adiabatic_lambda - b/(hf_energy(1) - dfa_energy - hf_energy(2))
 
      dexc_da = 1.0_dp
      dexc_db = 1.0_dp/c - (log(abs(1.0_dp + c))/(c*c))
      dexc_dc = -b/(c*c*c*(1.0_dp + c))*(2.0_dp*c + c*c - 2.0_dp*log(abs(1.0_dp + c)) - 2.0_dp*log(abs(1.0_dp + c))*c)
 
      da_dex1 = 1.0_dp
      da_ddw0 = 0.0_dp
      da_ddfa = 0.0_dp
      da_dex2 = 0.0_dp
 
      db_dex1 = 0.0_dp
      db_ddw0 = 1.0_dp
      db_ddfa = 0.0_dp
      db_dex2 = 0.0_dp
 
      dc_dex1 = b/(hf_energy(1) - dfa_energy - hf_energy(2))**2
      dc_ddw0 = -1.0_dp/(hf_energy(1) - dfa_energy - hf_energy(2))
      dc_ddfa = -dc_dex1
      dc_dex2 = -dc_dex1
 
      scale_dex1 = dexc_da*da_dex1 + dexc_db*db_dex1 + dexc_dc*dc_dex1
      scale_ddw0 = dexc_da*da_ddw0 + dexc_db*db_ddw0 + dexc_dc*dc_ddw0
      scale_ddfa = dexc_da*da_ddfa + dexc_db*db_ddfa + dexc_dc*dc_ddfa
      scale_dex2 = dexc_da*da_dex2 + dexc_db*db_dex2 + dexc_dc*dc_dex2
 
      total_energy_xc = a + b/(c*c)*(c - log(abs(1.0_dp + c)))
      IF (do_swap_hf) THEN
         swap_value = scale_dex1
         scale_dex1 = scale_dex2
         scale_dex2 = swap_value
      END IF
   END SUBROUTINE rescale_mcy3_pade
 
END MODULE xc_adiabatic_methods