Branch data Line data Source code
1 : : !!****m* ABINIT/m_xciit
2 : : !! NAME
3 : : !! m_xciit
4 : : !!
5 : : !! FUNCTION
6 : : !! Exchange-correlation at finite temperature of an electron gas
7 : : !! Ichimaru S., Iyetomi H., Tanaka S., Phys. Rep. 149, 91-205 (1987) [[cite:Ichimaru1987]]
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2002-2026 ABINIT group (JFD,LK)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_xciit
24 : :
25 : : use defs_basis
26 : : use m_errors
27 : :
28 : : implicit none
29 : :
30 : : private
31 : : !!***
32 : :
33 : : public :: xciit
34 : : !!***
35 : :
36 : : contains
37 : : !!***
38 : :
39 : : !!****f* ABINIT/xciit
40 : : !! NAME
41 : : !! xciit
42 : : !!
43 : : !! FUNCTION
44 : : !! Exchange-correlation at finite temperature of an electron gas
45 : : !! Ichimaru S., Iyetomi H., Tanaka S., Phys. Rep. 149, 91-205 (1987) [[cite:Ichimaru1987]]
46 : : !!
47 : : !! INPUTS
48 : : !! temp= (electronic) temperature
49 : : !! npt=number of real space points
50 : : !! order=gives the maximal derivative of Exc computed.
51 : : !! rspts(npt)=Wigner-Seitz radii at each point
52 : : !!
53 : : !! OUTPUT
54 : : !! fxc(npt)=exchange-correlation free energy density (hartree)
55 : : !! tsxc(npt)=exchange-correlation entropy energy density (hartree)
56 : : !! vxc(npt)=exchange-correlation potential
57 : : !! --- optional output ---
58 : : !! [dvxc(npt)]=partial second derivatives of the xc energy
59 : : !!
60 : : !! SOURCE
61 : :
62 : 252 : subroutine xciit(fxc,tsxc,npt,order,rspts,temp,vxc, &
63 : : & dvxc)!Optional argument
64 : :
65 : : !Arguments ------------------------------------
66 : : !scalars
67 : : integer,intent(in) :: npt,order
68 : : real(dp),intent(in) :: temp
69 : : !arrays
70 : : real(dp),intent(in) :: rspts(npt)
71 : : real(dp),intent(out) :: fxc(npt),tsxc(npt),vxc(npt)
72 : : real(dp),intent(out),optional :: dvxc(npt)
73 : :
74 : : !Local variables-------------------------------
75 : : !scalars
76 : : integer :: ipt
77 : : real(dp) :: ef,deltavxc,Gamma,rs,rsm1,tt
78 : : character(len=500) :: msg
79 : :
80 : : ! *************************************************************************
81 : :
82 : : !Checks the values of order
83 [ - + ]: 252 : if(order<0.or.order>2)then
84 : 0 : write(msg, '(4a,i3,a)' ) ch10,&
85 : 0 : & 'With Ishimaru-Iyetomi-Tanka xc functional, the only',ch10,&
86 : 0 : & 'allowed values for order are 0, 1 or 2, while it is found to be ',order,'.'
87 : 0 : ABI_BUG(msg)
88 : : end if
89 : :
90 : : !Loop over grid points
91 [ + + ]: 972252 : do ipt=1,npt
92 : :
93 : 972000 : rs=rspts(ipt)
94 : 972000 : rsm1=one/rs
95 : : ! Step for the Vxc computation
96 : 972000 : deltavxc=0.01_dp
97 : : ! Compute ef
98 : 972000 : ef=0.5_dp*(9.0_dp*pi/4.0_dp)**(2.0_dp/3.0_dp)*rsm1**2
99 : : ! Compute temperature
100 : 972000 : tt=max(temp/ef,tol12)
101 : : ! Compute Gamma
102 : 972000 : Gamma=one/(tt*ef)/rs
103 : :
104 : : ! Exchange-correlation of Ichimaru functional
105 : 972000 : fxc(ipt)=fexsGamma(Gamma,tt)*rsm1
106 : : ! exc(ipt)=fxc(ipt) - tdexcsdtiit(rs,tt);
107 : 972000 : tsxc(ipt)=-tdexcsdtiit(rs,tt)
108 : : vxc(ipt)=(8.0_dp*(Fxc_iit(rs,tt,deltavxc)-Fxc_iit(rs,tt,-deltavxc)) &
109 : : & -(Fxc_iit(rs,tt,two*deltavxc)-Fxc_iit(rs,tt,-two*deltavxc)))/ &
110 : 972000 : & (12.0_dp*deltavxc*3.0_dp/(4.0_dp*pi)/rs**3)
111 [ - + ]: 972252 : if (order==2) then
112 : : dvxc(ipt)=(-30.0_dp*Fxc_iit(rs,tt,zero)+16.0_dp*(Fxc_iit(rs,tt,deltavxc)+Fxc_iit(rs,tt,-deltavxc)) &
113 : : & -(Fxc_iit(rs,tt,two*deltavxc)+Fxc_iit(rs,tt,two*deltavxc)))/ &
114 : 0 : & (12.0_dp*(deltavxc*3.0_dp/(4.0_dp*pi)/rs**3)**2)
115 : : end if
116 : : end do
117 : :
118 : : CONTAINS
119 : : !!***
120 : :
121 : : !!****f* ABINIT/fexsGamma
122 : : !!
123 : : !! NAME
124 : : !! fexsGamma
125 : : !!
126 : : !! FUNCTION
127 : : !! Free energy for the IIT finite temperautre XC functional
128 : : !!
129 : : !! INPUTS
130 : : !! Gamma= ?
131 : : !! t=temperature
132 : : !!
133 : : !! OUTPUT
134 : : !! fexsGamma=free energy
135 : : !!
136 : : !! SOURCE
137 : :
138 : 8748000 : function fexsGamma(Gamma,t)
139 : :
140 : : !Arguments ------------------------------------
141 : : real(dp) :: fexsGamma
142 : : real(dp),intent(in) :: Gamma,t
143 : : !Local variables-------------------------------
144 : : real(dp) :: lambda
145 : : real(dp) :: tanht,tanhst
146 : : real(dp) :: a,b,c,d,e
147 : : real(dp) :: bmcdse,amcse,sqrt4emd2
148 : :
149 : : ! *************************************************************************
150 : :
151 : 8748000 : lambda=(4.0_dp/(9.0_dp*pi))**(one/3.0_dp)
152 : 8748000 : tanht=tanh(one/t)
153 : 8748000 : tanhst=tanh(one/sqrt(t))
154 : :
155 : 8748000 : a=one/(pi*lambda)*(0.75_dp+3.04363_dp*t**2-0.09227_dp*t**3+1.7035_dp*t**4)/(one+8.31051_dp*t**2+5.1105_dp*t**4)*tanht
156 : 8748000 : b=(0.341308_dp+12.070873_dp*t**2+1.148889_dp*t**4)/(one+10.495346_dp*t**2+1.326623_dp*t**4)*sqrt(t)*tanhst
157 : 8748000 : e=(0.539409_dp+2.522206_dp*t**2+0.178484_dp*t**4)/(one+2.555501_dp*t**2+0.146319_dp*t**4)*t*tanht
158 : 8748000 : c=(0.872496_dp+0.025248_dp*exp(-1./t))*e
159 : 8748000 : d=(0.614925_dp+16.996055_dp*t**2+1.489056_dp*t**4)/(one+10.10935_dp*t**2+1.22184_dp*t**4)*sqrt(t)*tanhst
160 : :
161 : 8748000 : bmcdse=b-c*d/e
162 : 8748000 : amcse=a-c/e
163 : 8748000 : sqrt4emd2=sqrt(4.0_dp*e-d**2)
164 : :
165 : : fexsGamma=-one/Gamma*(c/e*Gamma+2.0_dp/e*bmcdse*sqrt(Gamma)+one/e*(amcse-d/e*bmcdse)*log(e*Gamma+d*sqrt(Gamma)+one)- &
166 : 8748000 : & 2.0_dp/(e*sqrt4emd2)*(d*amcse+(2.0_dp-d**2/e)*bmcdse)*(atan((2.0_dp*e*sqrt(Gamma)+d)/sqrt4emd2)-atan(d/sqrt4emd2)))
167 : :
168 : 8748000 : end function fexsGamma
169 : : !!***
170 : :
171 : : !!****f* ABINIT/Fxc_iit
172 : : !!
173 : : !! NAME
174 : : !! Fxc_iit
175 : : !!
176 : : !! FUNCTION
177 : : !! Auxiliary function for the IIT finite temperature XC functional
178 : : !!
179 : : !! INPUTS
180 : : !! deltavxc= ?
181 : : !! rs=Wigner-Seitz radius
182 : : !! t=temperature
183 : : !!
184 : : !! OUTPUT
185 : : !! Fxc_iit=auxiliary function
186 : : !!
187 : : !! SOURCE
188 : :
189 : 3888000 : function Fxc_iit(rs,t,deltavxc)
190 : :
191 : : !Arguments ------------------------------------
192 : : real(dp) :: Fxc_iit
193 : : real(dp),intent(in) :: rs,t,deltavxc
194 : : !Local variables-------------------------------
195 : : real(dp) :: newrs,newt,newGamma
196 : :
197 : : ! *************************************************************************
198 : :
199 : 3888000 : newrs=rs/(one+deltavxc)**(one/3.0_dp)
200 : 3888000 : newt=t/(one+deltavxc)**(2.0_dp/3.0_dp)
201 : 3888000 : newGamma=2.0_dp*(4.0_dp/(9.0_dp*pi))**(2.0_dp/3.0_dp)*newrs/newt
202 : 3888000 : Fxc_iit=3.0_dp/(4.0_dp*pi)*fexsGamma(newGamma,newt)/newrs**4
203 : :
204 : 3888000 : end function Fxc_iit
205 : : !!***
206 : :
207 : : !!****f* ABINIT/tdexcsdtiit
208 : : !!
209 : : !! NAME
210 : : !! tdexcsdtiit
211 : : !!
212 : : !! FUNCTION
213 : : !! Auxiliary function for the IIT finite temperature XC functional
214 : : !!
215 : : !! INPUTS
216 : : !! rs=Wigner-Seitz radius
217 : : !! t=temperature
218 : : !!
219 : : !! OUTPUT
220 : : !! tdexcsdtiit=auxiliary function
221 : : !!
222 : : !! SOURCE
223 : :
224 : 972000 : function tdexcsdtiit(rs,t)
225 : :
226 : : !Arguments ------------------------------------
227 : : real(dp) :: tdexcsdtiit
228 : : real(dp),intent(in) :: rs,t
229 : : !Local variables-------------------------------
230 : : real(dp) :: ef,Gamma
231 : : real(dp) :: deltat=1.0d-2
232 : :
233 : : ! *************************************************************************
234 : :
235 : 972000 : ef=half*(9.0_dp*pi/4.0_dp)**(2.0_dp/3.0_dp)/rs**2
236 : 972000 : Gamma=one/(t*ef)/rs
237 : : tdexcsdtiit=8.0_dp*(fexsGamma(Gamma/(one+deltat),(one+deltat)*t) &
238 : : & -fexsGamma(Gamma/(one-deltat),(one-deltat)*t)) &
239 : : & -(fexsGamma(Gamma/(one+2.0_dp*deltat),(one+2.0_dp*deltat)*t) &
240 : 972000 : & -fexsGamma(Gamma/(one-2.0_dp*deltat),(one-2.0_dp*deltat)*t))
241 : 972000 : tdexcsdtiit=t*tdexcsdtiit/(12.0_dp*deltat*t)/rs
242 : :
243 : 972000 : end function tdexcsdtiit
244 : : !!***
245 : :
246 : : end subroutine xciit
247 : : !!***
248 : :
249 : : end module m_xciit
250 : : !!***
|
