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