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_xcpz
8 : use m_juDFT
9 : !-----------------------------------------------------------------------
10 : ! Called in case of icorr=5 : spin-polarized exchange-correlation
11 : ! from Ceperley-Alder monte carlo results with parametrization
12 : ! of Perdew,Zunger,PRB 23, 5048 (1981)
13 :
14 : ! krla=1: Relativistic correction of exchange energy and potential
15 : ! related to Dirac kinetic energy, according to:
16 : ! A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
17 :
18 : ! be careful: calculation in rydberg!
19 :
20 : ! vxcpz calculates the XC-potential and
21 : ! excpz calculates the XC-energy
22 :
23 : ! based on a subroutine by S. Bluegel; r.pentcheva 22.01.96
24 : !-----------------------------------------------------------------------
25 :
26 : USE m_constants, ONLY : pi_const
27 : USE m_relcor
28 : IMPLICIT NONE
29 :
30 : REAL, PARAMETER, PRIVATE :: cvx = 1.221774115422 ! 2 * ( 3/(2*pi) )^(2/3)
31 : REAL, PARAMETER, PRIVATE :: d_15 = 1.e-15 , c76 = 7.0 / 6.0
32 : REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0
33 : REAL, PARAMETER, PRIVATE :: half = 0.5 , thrd = one/three
34 : REAL, PARAMETER, PRIVATE :: thrhalf = three * half , two = 2.0
35 : REAL, PARAMETER, PRIVATE :: c23 = two * thrd , c43 = four * thrd
36 : REAL, PARAMETER, PRIVATE :: ap = 0.0622 , af = 0.0311
37 : REAL, PARAMETER, PRIVATE :: bp = -0.0960 , bf = -0.0538
38 : REAL, PARAMETER, PRIVATE :: cp = 0.0040 , cf = 0.0014
39 : REAL, PARAMETER, PRIVATE :: dp = -0.0232 , df = -0.0096
40 : REAL, PARAMETER, PRIVATE :: gp = -0.2846 , gf = -0.1686
41 : REAL, PARAMETER, PRIVATE :: b1p = 1.0529 , b1f = 1.3981
42 : REAL, PARAMETER, PRIVATE :: b2p = 0.3334 , b2f = 0.2611
43 :
44 : CONTAINS
45 : !************************************************************************
46 127 : SUBROUTINE vxcpz( &
47 : iofile,krla,jspins, &
48 127 : mgrid,ngrid,rh, &
49 127 : vx,vxc)
50 : !************************************************************************
51 :
52 : ! .. Scalar Arguments ..
53 : INTEGER, INTENT (IN) :: jspins
54 : INTEGER, INTENT (IN) :: krla ! run mode parameters
55 : INTEGER, INTENT (IN) :: iofile ! file number for read and write
56 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
57 :
58 : ! .. Array Arguments ..
59 : REAL, INTENT (IN) :: rh(:,:) ! charge density
60 : REAL, INTENT (OUT) :: vx(:,:),vxc(:,:) ! x/xc potential
61 :
62 : ! .. Local Scalars ..
63 : REAL :: c_2, cbrt1, cbrt2, dfds, decds, dec1, dec2, vcf, vcp
64 :
65 : ! .. Local Arrays ..
66 127 : REAL, ALLOCATABLE :: psi(:) ! relativistic exchange potential corr.
67 :
68 : !-----s Intrinsic Functions
69 : INTRINSIC max
70 : REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
71 : REAL :: fothrd, thfpi, c_1, y1, y2, s, fs, rs
72 : REAL :: ecp, ecf
73 : INTEGER :: ir
74 :
75 127 : fothrd = c43
76 127 : thfpi = three / ( four * pi_const )
77 :
78 : !-----> evaluate relativistic corrections for exchange
79 :
80 381 : ALLOCATE ( psi(ngrid) )
81 : CALL relcor( &
82 : mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
83 127 : psi)
84 :
85 127 : IF ( jspins == 2 ) THEN ! spinpolarized calculation
86 :
87 : c_1 = one / ( two**fothrd - two )
88 306007 : DO ir = 1,ngrid
89 305980 : rh1 = max(d_15,rh(ir,1))
90 305980 : rh2 = max(d_15,rh(ir,jspins))
91 305980 : rho = rh1 + rh2
92 305980 : y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1 ! s = (rh2 - rh1) / rho
93 305980 : cbrt1 = (one-s) ** thrd
94 305980 : cbrt2 = (one+s) ** thrd
95 305980 : fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
96 305980 : dfds = c_1 * fothrd * (cbrt1 - cbrt2)
97 305980 : rs = ( thfpi/rho )**thrd
98 :
99 305980 : IF (rs >= one) THEN
100 42489 : ecp = fecl(rs,gp,b1p,b2p) ! correlation energy paramagnetic
101 42489 : ecf = fecl(rs,gf,b1f,b2f) ! correlation energy ferromagnetic
102 42489 : vcp = fvcl(ecp,rs,b1p,b2p) ! d(rho*ecp)/d(rho) = ecp - rs/3*d(ecp)/d(rs)
103 42489 : vcf = fvcl(ecf,rs,b1f,b2f) ! d(rho*ecf)/d(rho)
104 : ELSE
105 263491 : ecp = fecs(rs,ap,bp,cp,dp)
106 263491 : ecf = fecs(rs,af,bf,cf,df)
107 263491 : vcp = fvcs(rs,ap,bp,cp,dp)
108 263491 : vcf = fvcs(rs,af,bf,cf,df)
109 : ENDIF
110 305980 : decds = (ecf-ecp)*dfds ! = d(ec)/d(rho)
111 305980 : dec1 = vcp + (vcf-vcp)*fs ! = d(rho*ec)/d(rho)
112 305980 : dec2 = two/rho*decds ! = 2/rho*d(ec)/ds
113 :
114 305980 : c_2 = cvx / rs * psi(ir) ! exchange potential muxp=-cvx/rs= 4/3*ex
115 305980 : vxc(ir,1) =dec1+dec2*rh2 - c_2*cbrt1 ! muxp*(2x)**(1/3)
116 305980 : vxc(ir,jspins)=dec1-dec2*rh1 - c_2*cbrt2 ! calculate exchange correlation potential
117 : ! vc = ec +vcp + (vcf-vcp)f(s) - (ecf-ecp)df/ds*(s+/-1)
118 :
119 305980 : vx (ir,1) = - c_2*cbrt1
120 306007 : vx (ir,jspins)= - c_2*cbrt2
121 : ENDDO
122 :
123 100 : ELSEIF ( jspins == 1 ) THEN ! non-spinpolarized calculation
124 :
125 2964876 : DO ir = 1,ngrid
126 2964776 : rho = max(d_15,rh(ir,1))
127 2964776 : rs = ( thfpi/rho )**thrd
128 2964776 : IF (rs >= one) THEN
129 374448 : ecp = fecl(rs,gp,b1p,b2p)
130 374448 : vcp = fvcl(ecp,rs,b1p,b2p)
131 : ELSE
132 2590328 : ecp = fecs(rs,ap,bp,cp,dp)
133 2590328 : vcp = fvcs(rs,ap,bp,cp,dp)
134 : ENDIF
135 2964776 : vxc(ir,1) = vcp - cvx / rs * psi(ir)
136 :
137 2964876 : vx (ir,1) = cvx / rs * psi(ir)
138 : ENDDO
139 :
140 : ELSE
141 0 : WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
142 0 : CALL juDFT_error("vxcpz",calledby="xcpz")
143 : ENDIF
144 :
145 127 : DEALLOCATE (psi)
146 127 : RETURN
147 :
148 : END SUBROUTINE vxcpz
149 : !***********************************************************************
150 15 : SUBROUTINE excpz( &
151 : iofile,krla,jspins, &
152 15 : mgrid,ngrid,rh, &
153 15 : exc)
154 : !***********************************************************************
155 :
156 : ! .. Scalar Arguments ..
157 : INTEGER, INTENT (IN) :: jspins
158 : INTEGER, INTENT (IN) :: krla ! run mode parameters
159 : INTEGER, INTENT (IN) :: iofile ! file number for read and write
160 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
161 :
162 : ! .. Array Arguments ..
163 : REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
164 : REAL, INTENT (OUT) :: exc(mgrid) ! xc energy
165 :
166 : ! .. Local Scalars ..
167 : REAL :: ec, ex, cex
168 :
169 : ! .. Local Arrays ..
170 15 : REAL, ALLOCATABLE :: phi(:) ! relativistic exchange energy correct.
171 :
172 : !-----> Intrinsic Functions
173 : INTRINSIC max
174 :
175 : REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
176 : REAL :: fothrd, thfpi, c_1, y1, y2, s, fs, rs
177 : REAL :: ecp, ecf
178 : INTEGER :: ir
179 :
180 15 : fothrd = c43
181 15 : thfpi = three / ( four * pi_const )
182 15 : cex = cvx / c43
183 :
184 45 : ALLOCATE ( phi(ngrid) )
185 : CALL relcor( &
186 : mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
187 15 : phi)
188 :
189 15 : IF ( jspins == 2 ) THEN ! spinpolarized calculation
190 :
191 : c_1 = one / ( two**fothrd - two )
192 284383 : DO ir = 1,ngrid
193 284380 : rh1 = max(d_15,rh(ir,1))
194 284380 : rh2 = max(d_15,rh(ir,jspins))
195 284380 : rho = rh1 + rh2
196 284380 : rs = ( thfpi/rho )**thrd
197 284380 : y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1 ! s = (rh2 - rh1) / rho
198 284380 : fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
199 284380 : IF (rs >= one) THEN
200 38220 : ecp = fecl(rs,gp,b1p,b2p) ! correlation energy paramagnetic
201 38220 : ecf = fecl(rs,gf,b1f,b2f) ! correlation energy ferromagnetic
202 : ELSE
203 246160 : ecp = fecs(rs,ap,bp,cp,dp)
204 246160 : ecf = fecs(rs,af,bf,cf,df)
205 : ENDIF
206 :
207 284380 : ec = ecp + (ecf-ecp)*fs ! total correlation energy
208 284380 : ex = -cex/rs* (one + 0.2599210482*fs) ! ex = exp + (exf-exp)*f(s)
209 : ! exf-exp = (2**(1/3)-1) * exp
210 284383 : exc(ir) = ec + ex*phi(ir)
211 : ENDDO
212 :
213 12 : ELSEIF ( jspins == 1 ) THEN ! non-spinpolarized calculation
214 :
215 2866844 : DO ir = 1,ngrid
216 2866832 : rho = max(d_15,rh(ir,1))
217 2866832 : rs = ( thfpi/rho )**thrd
218 2866832 : IF (rs >= one) THEN
219 355232 : ecp = fecl(rs,gp,b1p,b2p)
220 : ELSE
221 2511600 : ecp = fecs(rs,ap,bp,cp,dp)
222 : ENDIF
223 2866832 : ex = -cex/rs
224 2866844 : exc(ir) = ecp + ex*phi(ir)
225 : ENDDO
226 :
227 : ELSE
228 0 : WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
229 0 : CALL juDFT_error("vxcpz",calledby="xcpz")
230 : ENDIF
231 :
232 15 : DEALLOCATE (phi)
233 15 : RETURN
234 :
235 : END SUBROUTINE excpz
236 :
237 : !--------------------------------------------------------------------
238 810389 : REAL FUNCTION fecl(r,g,b1,b2)
239 : REAL :: r,g,b1,b2
240 810389 : fecl = g / ( one + b1*sqrt(r) + b2*r )
241 0 : END FUNCTION fecl
242 459426 : REAL FUNCTION fvcl(ce,r,b1,b2)
243 : REAL :: ce,r,b1,b2
244 459426 : fvcl = ce* (one+c76*b1*sqrt(r)+c43*b2*r)/(one+b1*sqrt(r)+b2*r)
245 459426 : END FUNCTION fvcl
246 5611579 : REAL FUNCTION fecs(r,a,b,c,d)
247 : REAL :: r,a,b,c,d
248 : INTRINSIC alog
249 5611579 : fecs = a*alog(r) + b + c*r*alog(r) + d*r
250 0 : END FUNCTION fecs
251 3117310 : REAL FUNCTION fvcs(r,a,b,c,d)
252 : REAL :: r,a,b,c,d
253 : INTRINSIC alog
254 : fvcs = a*alog(r) + (b-a/three) + c23*c*r*alog(r) + &
255 3117310 : (two*d-c)*r/three
256 3117310 : END FUNCTION fvcs
257 : !--------------------------------------------------------------------
258 :
259 : END MODULE m_xcpz
|
