-
Notifications
You must be signed in to change notification settings - Fork 1
/
xc_adiabatic_methods.F
133 lines (114 loc) · 5.92 KB
/
xc_adiabatic_methods.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 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