Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_xcbh
8 : use m_juDFT
9 : !-----------------------------------------------------------------------
10 : ! Called in case of icorr=2,3 : spin-polarized exchange-correlation
11 : ! potential of U. von Barth and L. Hedin, J.Phys.C5,1629 (1972)
12 : ! icorr = 2: parametrization of Moruzzi,Janak,Williams
13 : ! icorr = 3: parametrization of von Barth and Hedin
14 :
15 : ! krla=1: Relativistic correction of exchange energy and potential
16 : ! related to Dirac kinetic energy, according to:
17 : ! A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
18 :
19 : ! be careful: calculation in rydberg!
20 :
21 : ! vxcbh calculates the XC-potential and
22 : ! excbh calculates the XC-energy
23 :
24 : ! based on a subroutine by S. Bluegel; r.pentcheva 22.01.96
25 : !-----------------------------------------------------------------------
26 :
27 : USE m_constants, ONLY : pi_const
28 : USE m_relcor
29 : IMPLICIT NONE
30 :
31 : REAL, PARAMETER, PRIVATE :: ff = 3.847322101863 ! 1 / ( 2^(1/3) - 1 )
32 : REAL, PARAMETER, PRIVATE :: cvx = 1.221774115422 ! 2 * ( 3/(2*pi) )^(2/3)
33 : REAL, PARAMETER, PRIVATE :: cpmjw = 0.045 , cfmjw = 0.0225
34 : REAL, PARAMETER, PRIVATE :: cpvbh = 0.0504 , cfvbh = 0.0254
35 : REAL, PARAMETER, PRIVATE :: rpmjw = 21.0 , rfmjw = 52.916684096
36 : REAL, PARAMETER, PRIVATE :: rpvbh = 30.0 , rfvbh = 75.0
37 : REAL, PARAMETER, PRIVATE :: d_15 = 1.e-15
38 : REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0
39 : REAL, PARAMETER, PRIVATE :: half = 0.5 , thrd = one/three
40 : REAL, PARAMETER, PRIVATE :: hfthrd = 0.79370052705 ! 2^(-1/3)
41 : REAL, PARAMETER, PRIVATE :: thrhalf = three * half
42 : REAL, PARAMETER, PRIVATE :: fothrd = four * thrd , two = 2.0
43 :
44 : CONTAINS
45 : !************************************************************************
46 0 : SUBROUTINE vxcbh &
47 : (iofile,xcpot,jspins, &
48 0 : mgrid,ngrid,rh, &
49 0 : vx,vxc)
50 : !************************************************************************
51 : USE m_types_xcpot_data
52 :
53 : ! .. Scalar Arguments ..
54 : INTEGER, INTENT (IN) :: jspins
55 : TYPE(t_xcpot_data), INTENT (IN) :: xcpot ! run mode parameters
56 : INTEGER, INTENT (IN) :: iofile ! file number for read and write
57 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
58 :
59 : ! .. Array Arguments ..
60 : REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
61 : REAL, INTENT (OUT) :: vxc(mgrid,jspins) ! xc potential
62 : REAL, INTENT (OUT) :: vx (mgrid,jspins) ! x potential
63 :
64 : ! .. Local Scalars ..
65 : REAL :: txthrd,tythrd,muxp,mucp,mucf,ecfmp,tauc,mucnm
66 : REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
67 : REAL :: x, y, cp, cf, rp, rf, rs, ecprs, ecfrs
68 : INTEGER :: ir
69 :
70 : ! .. Local Arrays ..
71 0 : REAL, ALLOCATABLE :: psi(:) ! relativistic exchange potential corr.
72 :
73 : !---- Intrinsic Functions
74 : INTRINSIC alog,max
75 :
76 : !-----> evaluate relativistic corrections for exchange
77 :
78 0 : ALLOCATE ( psi(ngrid) )
79 : CALL relcor( &
80 : mgrid,ngrid,jspins,xcpot%krla, .TRUE. ,rh, &
81 0 : psi)
82 :
83 : !-----> select exchange correlation potential
84 :
85 0 : IF (xcpot%is_mjw) THEN
86 : cp = cpmjw ; cf = cfmjw
87 : rp = rpmjw ; rf = rfmjw
88 0 : ELSEIF (xcpot%is_bh) THEN
89 : cp = cpvbh ; cf = cfvbh
90 : rp = rpvbh ; rf = rfvbh
91 : ELSE
92 0 : WRITE (iofile,FMT=2000)
93 0 : CALL juDFT_error("BUG:vxcbh",calledby="xcbh")
94 : END IF
95 : 2000 FORMAT (13x,'set key for exchange-correlation potential')
96 :
97 : !-----> calculate exchange correlation potential
98 :
99 0 : IF ( jspins == 2) THEN ! spinpolarized calculation
100 :
101 0 : DO ir = 1,ngrid ! loop over realspace gridpoints
102 0 : rh1 = max(d_15,rh(ir,1))
103 0 : rh2 = max(d_15,rh(ir,jspins))
104 0 : rho = rh1 + rh2
105 0 : x = rh1/rho
106 0 : y = rh2/rho
107 0 : txthrd = (2*x)**thrd
108 0 : tythrd = (2*y)**thrd
109 0 : rs= (four*pi_const*rho/three)**thrd
110 0 : rs = 1/rs
111 :
112 0 : ecprs = -cp*fc(rs/rp) ! calculate correlation energy
113 0 : ecfrs = -cf*fc(rs/rf) ! p : paramagnetic, f : ferromagnetic
114 : ! x : exchange, c : correlation
115 0 : muxp = -psi(ir)* (cvx/rs) ! paramagnetic exchange potential
116 : ! (psi contains rel. corr.)
117 0 : mucp = -cp*alog(one+rp/rs) ! calculate correlation potential
118 0 : mucf = -cf*alog(one+rf/rs)
119 0 : ecfmp = fothrd * (ecfrs-ecprs)
120 0 : tauc = mucf - mucp - ecfmp
121 0 : mucnm = mucp + tauc*fex(x) - ff*ecfmp
122 :
123 0 : vxc(ir,1) = mucnm + (muxp+ff*ecfmp)*txthrd ! collect correlation
124 0 : vxc(ir,jspins) = mucnm + (muxp+ff*ecfmp)*tythrd ! and exchange parts
125 :
126 0 : vx (ir,1) = muxp*txthrd
127 0 : vx (ir,jspins) = muxp*tythrd
128 : ENDDO
129 :
130 0 : ELSEIF ( jspins == 1 ) THEN ! non - spinpolarized calculation
131 :
132 0 : DO ir = 1, ngrid ! loop over realspace gridpoints
133 0 : rh1 = max(d_15,rh(ir,1))
134 0 : rs = (four*pi_const*rh1/three)**thrd
135 0 : rs = 1/rs
136 0 : muxp = -psi(ir) * (cvx/rs) ! paramagnetic exchange potential
137 : ! (psi contains rel. corr.)
138 0 : mucp = -cp* alog(one+rp/rs) ! calculate correlation potential
139 0 : vxc(ir,1) = mucp + muxp ! collect correlation & exchange part
140 :
141 0 : vx (ir,1) = muxp
142 : ENDDO
143 :
144 : ELSE
145 0 : WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
146 0 : CALL juDFT_error("vxcbh",calledby="xcbh")
147 : ENDIF
148 :
149 0 : DEALLOCATE (psi)
150 0 : RETURN
151 :
152 : END SUBROUTINE vxcbh
153 : !***********************************************************************
154 0 : SUBROUTINE excbh &
155 : (iofile,xcpot,jspins, &
156 0 : mgrid,ngrid,rh, &
157 0 : exc)
158 : !***********************************************************************
159 : USE m_types_xcpot_data
160 :
161 : ! .. Scalar Arguments ..
162 : INTEGER, INTENT (IN) :: jspins
163 : TYPE(t_xcpot_data), INTENT (IN) :: xcpot ! run mode parameters
164 : INTEGER, INTENT (IN) :: iofile ! file number for read and write
165 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
166 :
167 : ! .. Array Arguments ..
168 : REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
169 : REAL, INTENT (OUT) :: exc(mgrid) ! xc energy
170 :
171 : ! .. Local Scalars ..
172 : REAL :: thfpi,thrquart,exprs,exfrs,excprs,excfrs
173 : REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
174 : REAL :: x, y, cp, cf, rp, rf, rs, ecprs, ecfrs
175 : INTEGER :: ir
176 :
177 : ! .. Local Arrays ..
178 0 : REAL, ALLOCATABLE :: phi(:) ! relativistic exchange energy correct.
179 :
180 : !-----> Intrinsic Functions
181 : INTRINSIC alog,max
182 :
183 0 : thrquart = 0.75
184 0 : thfpi = thrquart/pi_const
185 :
186 0 : ALLOCATE ( phi(ngrid) )
187 : CALL relcor( &
188 : mgrid,ngrid,jspins,xcpot%krla, .FALSE. ,rh, &
189 0 : phi)
190 :
191 : !-----> select exchange correlation potential
192 :
193 0 : IF (xcpot%is_mjw) THEN
194 : cp = cpmjw ; cf = cfmjw
195 : rp = rpmjw ; rf = rfmjw
196 0 : ELSEIF (xcpot%is_bh) THEN
197 : cp = cpvbh ; cf = cfvbh
198 : rp = rpvbh ; rf = rfvbh
199 : ELSE
200 0 : WRITE (iofile,FMT=2000)
201 0 : CALL juDFT_error("excbh",calledby="xcbh")
202 : END IF
203 : 2000 FORMAT (13x,'set key for exchange-correlation potential')
204 :
205 0 : IF ( jspins == 2) THEN ! spinpolarized calculation
206 :
207 0 : DO ir = 1,ngrid ! loop over realspace gridpoints
208 0 : rh1 = max(d_15,rh(ir,1))
209 0 : rh2 = max(d_15,rh(ir,jspins))
210 0 : rho = rh1 + rh2
211 0 : x = rh1/rho
212 0 : rs= (thfpi/rho)**thrd
213 :
214 0 : exprs = -phi(ir)*thrquart*cvx/rs ! first exchange energy
215 0 : exfrs = (2.0**thrd)*exprs ! phi contains rel. corr.
216 :
217 0 : ecprs = -cp*fc(rs/rp) ! calculate correlation energy
218 0 : ecfrs = -cf*fc(rs/rf) ! p: paramagnetic, f: ferromagn.
219 :
220 0 : excprs = exprs + ecprs ! now add correlation energy
221 0 : excfrs = exfrs + ecfrs
222 :
223 0 : exc(ir) = excprs + (excfrs-excprs)*fex(x) ! collect all terms
224 : ENDDO
225 :
226 0 : ELSEIF ( jspins == 1 ) THEN ! non - spinpolarized calculation
227 :
228 0 : DO ir = 1,ngrid ! loop over realspace gridpoints
229 0 : rh1 = max(d_15,rh(ir,1))
230 0 : rs = (thfpi/rh1)**thrd
231 0 : exprs = -phi(ir)*thrquart*cvx/rs ! exchange energy ; phi contains
232 : ! relativistic correctionS
233 0 : ecprs = -cp*fc(rs/rp) ! calculate correlation energy
234 0 : exc(ir) = exprs + ecprs ! add correlation energy
235 : ENDDO
236 :
237 : ELSE
238 0 : WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
239 0 : CALL juDFT_error("excbh",calledby="xcbh")
240 : ENDIF
241 :
242 0 : DEALLOCATE (phi)
243 0 : RETURN
244 :
245 : END SUBROUTINE excbh
246 : !--------------------------------------------------------------------
247 0 : REAL FUNCTION fc(x)
248 : REAL :: x
249 : fc = (one+(x)*(x)*(x))*alog(one+one/(x)) &
250 0 : + half*(x) - (x)*(x) - thrd
251 0 : END FUNCTION fc
252 0 : REAL FUNCTION fex(x)
253 : REAL :: x
254 0 : fex = ff/hfthrd*((x)**fothrd +(1-(x))**fothrd - hfthrd)
255 0 : END FUNCTION fex
256 : !--------------------------------------------------------------------
257 :
258 : END MODULE m_xcbh
|
