1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Contains some functions used in the context of adiabatic hybrid functionals 8!> \par History 9!> 01.2008 created [Manuel Guidon] 10!> \author Manuel Guidon 11! ************************************************************************************************** 12MODULE xc_adiabatic_methods 13 14 USE input_constants, ONLY: do_potential_coulomb 15 USE kinds, ONLY: dp 16 USE mathconstants, ONLY: oorootpi 17 USE qs_energy_types, ONLY: qs_energy_type 18 USE qs_environment_types, ONLY: get_qs_env,& 19 qs_environment_type 20 USE qs_mo_types, ONLY: get_mo_set,& 21 mo_set_p_type 22#include "./base/base_uses.f90" 23 24 IMPLICIT NONE 25 PRIVATE 26 27 PUBLIC :: rescale_MCY3_pade 28 29CONTAINS 30 31! ************************************************************************************************** 32!> \brief - Calculates rescaling factors for XC potentials and energy expression 33!> \param qs_env ... 34!> \param hf_energy Array of size 2 containing the two HF energies (Ex^{HF} and Ex^{HF,LR} 35!> \param energy QS energy type 36!> \param adiabatic_lambda , adiabatic_omega: Parameters for adiabatic connection 37!> \param adiabatic_omega ... 38!> \param scale_dEx1 scaling coefficient for xc-potentials to be calculated 39!> \param scale_ddW0 scaling coefficient for xc-potentials to be calculated 40!> \param scale_dDFA scaling coefficient for xc-potentials to be calculated 41!> \param scale_dEx2 scaling coefficient for xc-potentials to be calculated 42!> \param total_energy_xc will contain the full xc energy 43!> \par History 44!> 09.2007 created [Manuel Guidon] 45!> \author Manuel Guidon 46!> \note 47!> - Model for adiabatic connection: 48!> 49!> W_lambda = a + b*lambda/(1+c*lambda) 50!> Exc = a + b*(c-log(1+c)/c^2) 51!> a = Ex^{HF} 52!> b = -c1*2*omega/sqrt(PI)*nelectron 53!> c = -1/lambda - b/(Ex^{HF}-W_lambda^{BLYP} + c2*W_lambda^{B88,LR}-c3*W_lambda^{HF,LR} 54! ************************************************************************************************** 55 SUBROUTINE rescale_MCY3_pade(qs_env, hf_energy, energy, adiabatic_lambda, & 56 adiabatic_omega, scale_dEx1, scale_ddW0, scale_dDFA, & 57 scale_dEx2, total_energy_xc) 58 59 TYPE(qs_environment_type), POINTER :: qs_env 60 REAL(dp), INTENT(INOUT) :: hf_energy(*) 61 TYPE(qs_energy_type), POINTER :: energy 62 REAL(dp), INTENT(IN) :: adiabatic_lambda, adiabatic_omega 63 REAL(dp), INTENT(INOUT) :: scale_dEx1, scale_ddW0, scale_dDFA, & 64 scale_dEx2, total_energy_xc 65 66 INTEGER :: nelec_a, nelec_b, nelectron, nspins 67 LOGICAL :: do_swap_hf 68 REAL(dp) :: a, b, c, c1, da_dDFA, da_ddW0, da_dEx1, da_dEx2, db_dDFA, db_ddW0, db_dEx1, & 69 db_dEx2, dc_dDFA, dc_ddW0, dc_dEx1, dc_dEx2, dExc_da, dExc_db, dExc_dc, dfa_energy, & 70 swap_value 71 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 72 73 do_swap_hf = .FALSE. 74 !! Assume the first HF section is the Coulomb one 75 IF (qs_env%x_data(1, 1)%potential_parameter%potential_type /= do_potential_coulomb) do_swap_hf = .TRUE. 76 77 IF (do_swap_hf) THEN 78 swap_value = hf_energy(1) 79 hf_energy(1) = hf_energy(2) 80 hf_energy(2) = swap_value 81 END IF 82 83 c1 = 0.23163_dp 84 85 CALL get_qs_env(qs_env=qs_env, mos=mos) 86 CALL get_mo_set(mo_set=mos(1)%mo_set, nelectron=nelec_a) 87 nspins = SIZE(mos) 88 IF (nspins == 2) THEN 89 CALL get_mo_set(mo_set=mos(2)%mo_set, nelectron=nelec_b) 90 ELSE 91 nelec_b = 0 92 END IF 93 nelectron = nelec_a + nelec_b 94 dfa_energy = energy%exc + energy%exc1 95 a = hf_energy(1) 96 b = -c1*2.0_dp*adiabatic_omega*oorootpi*nelectron !-0.23163_dp*2.0_dp*0.2_dp*oorootpi*nelectron 97 c = -1.0_dp/adiabatic_lambda - b/(hf_energy(1) - dfa_energy - hf_energy(2)) 98 99 dExc_da = 1.0_dp 100 dExc_db = 1.0_dp/c - (LOG(ABS(1.0_dp + c))/(c*c)) 101 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) 102 103 da_dEx1 = 1.0_dp 104 da_ddW0 = 0.0_dp 105 da_dDFA = 0.0_dp 106 da_dEx2 = 0.0_dp 107 108 db_dEx1 = 0.0_dp 109 db_ddW0 = 1.0_dp 110 db_dDFA = 0.0_dp 111 db_dEx2 = 0.0_dp 112 113 dc_dEx1 = b/(hf_energy(1) - dfa_energy - hf_energy(2))**2 114 dc_ddW0 = -1.0_dp/(hf_energy(1) - dfa_energy - hf_energy(2)) 115 dc_dDFA = -dc_dEx1 116 dc_dEx2 = -dc_dEx1 117 118 scale_dEx1 = dExc_da*da_dEx1 + dExc_db*db_dEx1 + dExc_dc*dc_dEx1 119 scale_ddW0 = dExc_da*da_ddW0 + dExc_db*db_ddW0 + dExc_dc*dc_ddW0 120 scale_dDFA = dExc_da*da_dDFA + dExc_db*db_dDFA + dExc_dc*dc_dDFA 121 scale_dEx2 = dExc_da*da_dEx2 + dExc_db*db_dEx2 + dExc_dc*dc_dEx2 122 123 total_energy_xc = a + b/(c*c)*(c - LOG(ABS(1.0_dp + c))) 124 IF (do_swap_hf) THEN 125 swap_value = scale_dEx1 126 scale_dEx1 = scale_dEx2 127 scale_dEx2 = swap_value 128 END IF 129 END SUBROUTINE rescale_MCY3_pade 130 131END MODULE xc_adiabatic_methods 132