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_xcvwn
8 : use m_juDFT
9 : !-----------------------------------------------------------------------
10 : ! Called in case of icorr=4 : spin-polarized exchange-correlation
11 : ! potential from Ceperley-Alder monte carlo results with parametrization
12 : ! of Vosko,Wilk,Nusair, Can. J. Phys. 58, 1200 (1980)
13 : ! ( In case of non-spinpolarization the fit of the correlation
14 : ! energy and potential is essentially equivallent to the
15 : ! parametrization of Perdew,Zunger,PRB 23, 5048 (1981), but
16 : ! in case of spin polarization the magnetization dependent
17 : ! interpolation formula for correlation is better.
18 : ! Part of the subroutine was supplied by M. Manninen. )
19 :
20 : ! krla=1: Relativistic correction of exchange energy and potential
21 : ! related to Dirac kinetic energy, according to:
22 : ! A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
23 :
24 : ! be careful: calculation in rydberg!
25 :
26 : ! vxc calculates the XC-potential and
27 : ! exc calculates the XC-energy
28 :
29 : ! based on a subroutine by S. Bluegel; r.pentcheva 22.01.96
30 : !-----------------------------------------------------------------------
31 :
32 : USE m_constants, ONLY : pi_const
33 : USE m_relcor
34 : IMPLICIT NONE
35 :
36 : REAL, PARAMETER, PRIVATE :: cex = 0.91633058742 ! 3/2 * ( 3/(2*pi) )^(2/3)
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 :: thrd = one/three , two = 2.0
40 : REAL, PARAMETER, PRIVATE :: fothrd = four * thrd
41 : REAL, PARAMETER, PRIVATE :: ap = 0.0621814 , xp0 = -0.10498
42 : REAL, PARAMETER, PRIVATE :: af = 0.0310907 , xf0 = -0.32500
43 : REAL, PARAMETER, PRIVATE :: al =-0.03377373, xl0 = -0.0047584
44 : REAL, PARAMETER, PRIVATE :: bp = 3.72744 , cp = 12.9352
45 : REAL, PARAMETER, PRIVATE :: bf = 7.06042 , cf = 18.0578
46 : REAL, PARAMETER, PRIVATE :: bl = 1.13107 , cl = 13.0045
47 : REAL, PARAMETER, PRIVATE :: qp = 6.1519908 , fdd0 = 1.70992093
48 : REAL, PARAMETER, PRIVATE :: qf = 4.7309269 , ql = 7.123109
49 :
50 : CONTAINS
51 : !************************************************************************
52 855 : SUBROUTINE vxcvwn( &
53 : iofile,krla,jspins, &
54 855 : mgrid,ngrid,rh, &
55 855 : vx,vxc)
56 : !************************************************************************
57 :
58 : ! .. Scalar Arguments ..
59 : INTEGER, INTENT (IN) :: jspins
60 : INTEGER, INTENT (IN) :: krla ! run mode parameters
61 : INTEGER, INTENT (IN) :: iofile ! file number for read and write
62 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
63 :
64 : ! .. Array Arguments ..
65 : REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
66 : REAL, INTENT (OUT) :: vxc(mgrid,jspins) ! xc potential
67 : REAL, INTENT (OUT) :: vx (mgrid,jspins) ! x potential
68 :
69 : ! .. Local Scalars ..
70 : REAL :: s3, s4, c_2, cbrt1, cbrt2, dfds, decdrp, decdrf
71 : REAL :: dacdr, dbdr, dec1, dec2, cvx, mucp
72 : REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
73 : REAL :: x, y1, y2, s, thfpi, c_1, rs, beta, bs41, alc
74 : REAL :: xpx, xfx, xlx, xpx0, xfx0, xlx0, ecp, ecf, fs, ec
75 : INTEGER :: ir
76 :
77 : ! .. Local Arrays ..
78 855 : REAL, ALLOCATABLE :: psi(:) ! relativistic exchange potential corr.
79 :
80 855 : thfpi = three / ( four * pi_const )
81 :
82 : !-----> evaluate relativistic corrections for exchange
83 :
84 2565 : ALLOCATE ( psi(ngrid) )
85 : CALL relcor( &
86 : mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
87 855 : psi)
88 :
89 855 : cvx = fothrd * cex
90 855 : IF ( jspins == 2 ) THEN ! spinpolarized calculation
91 :
92 : c_1 = one / ( two**fothrd - two )
93 18801798 : DO ir = 1,ngrid
94 18801039 : rh1 = max(d_15,rh(ir,1))
95 18801039 : rh2 = max(d_15,rh(ir,jspins))
96 18801039 : rho = rh1 + rh2
97 18801039 : y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1 ! s = (rh2 - rh1) / rho
98 18801039 : s3 = s**3 ; s4 = s*s3 ! spin polarisation
99 18801039 : cbrt1 = (one-s) ** thrd
100 18801039 : cbrt2 = (one+s) ** thrd
101 18801039 : fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
102 18801039 : dfds = c_1 * fothrd * ( cbrt1-cbrt2 )
103 18801039 : rs = ( thfpi/rho )**thrd
104 18801039 : x = sqrt(rs)
105 18801039 : xpx = x*x + bp*x + cp
106 18801039 : xfx = x*x + bf*x + cf
107 18801039 : xlx = x*x + bl*x + cl
108 18801039 : xpx0 = xp0*xp0 + bp*xp0 + cp
109 18801039 : xfx0 = xf0*xf0 + bf*xf0 + cf
110 18801039 : xlx0 = xl0*xl0 + bl*xl0 + cl
111 18801039 : ecp = fec(ap,x,xpx,xp0,xpx0,bp,qp) ! paramagnetic correlation energy
112 18801039 : ecf = fec(af,x,xfx,xf0,xfx0,bf,qf) ! ferromagnetic correlation ener.
113 18801039 : alc = fec(al,x,xlx,xl0,xlx0,bl,ql) ! alpha_c
114 18801039 : beta = fdd0* (ecf-ecp)/alc - one ! beta = (f"(0)*(ecf-ecp))/alc -1
115 18801039 : bs41 = one + beta * s**4
116 :
117 18801039 : ec = ecp + alc * fs / fdd0 * bs41 ! total correlation energy
118 18801039 : decdrp = fdedr(rho,ap,x,xpx,xp0,bp,cp) ! d(ecp)/d(rho)
119 18801039 : decdrf = fdedr(rho,af,x,xfx,xf0,bf,cf) ! d(ecf)/d(rho)
120 18801039 : dacdr = fdedr(rho,al,x,xlx,xl0,bl,cl) ! d(alc)/d(rho)
121 :
122 18801039 : dbdr =fdd0* ((decdrf-decdrp)*alc- (ecf-ecp)*dacdr)/alc**2 ! d(beta)/d(rho)
123 18801039 : dec1 =ec + rho* (decdrp+ (dacdr*fs*bs41+alc*fs*dbdr*s4)/fdd0) ! mucp= d(rho*ec)/d(rho)
124 18801039 : dec2 =two*alc/ (fdd0*rho)* (dfds*bs41+four*fs*beta*s3) ! 2/rho*d(ec)/ds
125 :
126 18801039 : c_2 = cvx / rs * psi(ir) ! exchange potential muxp=-cvx/rs= 4/3*ex
127 18801039 : vxc(ir,1) =dec1+dec2*rh2 - c_2*cbrt1 ! muxp*(2x)**(1/3)
128 18801039 : vxc(ir,jspins)=dec1-dec2*rh1 - c_2*cbrt2 ! calculate exchange correlation potential
129 :
130 18801039 : vx (ir,1) = - c_2*cbrt1
131 759 : vx (ir,jspins)= - c_2*cbrt2
132 : ENDDO
133 :
134 96 : ELSEIF ( jspins == 1 ) THEN ! non-spinpolarized calculation
135 :
136 789808 : DO ir = 1,ngrid
137 789712 : rho = max(d_15,rh(ir,1))
138 789712 : rs = ( thfpi/rho )**thrd
139 789712 : x = sqrt(rs)
140 789712 : xpx = x*x + bp*x + cp
141 789712 : xpx0 = xp0*xp0 + bp*xp0 + cp
142 :
143 789712 : ecp = fec(ap,x,xpx,xp0,xpx0,bp,qp) ! correlation energy paramagn.
144 789712 : decdrp = fdedr(rho,ap,x,xpx,xp0,bp,cp) ! d(ecp)/d(rho)
145 789712 : mucp = ecp + rho*decdrp ! d(rho*ec)/d(rho)
146 :
147 789712 : vxc(ir,1) = mucp - cvx/rs*psi(ir) ! -1.221774/rs :exchange potential
148 :
149 96 : vx (ir,1) = - cvx/rs*psi(ir)
150 : ENDDO
151 :
152 : ELSE
153 0 : WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
154 0 : CALL juDFT_error("vxcvwn",calledby="xcvwn")
155 : ENDIF
156 :
157 855 : DEALLOCATE (psi)
158 855 : RETURN
159 :
160 : END SUBROUTINE vxcvwn
161 : !***********************************************************************
162 162 : SUBROUTINE excvwn( &
163 : iofile,krla,jspins, &
164 162 : mgrid,ngrid,rh, &
165 162 : exc)
166 : !***********************************************************************
167 :
168 : ! .. Scalar Arguments ..
169 : INTEGER, INTENT (IN) :: jspins
170 : INTEGER, INTENT (IN) :: krla ! run mode parameters
171 : INTEGER, INTENT (IN) :: iofile ! file number for read and write
172 : INTEGER, INTENT (IN) :: mgrid,ngrid ! mesh points
173 :
174 : ! .. Array Arguments ..
175 : REAL, INTENT (IN) :: rh(mgrid,jspins) ! charge density
176 : REAL, INTENT (OUT) :: exc(mgrid) ! xc energy
177 :
178 : ! .. Local Scalars ..
179 : REAL :: c_1,ex
180 : REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
181 : REAL :: x, y1, y2, s, thfpi, rs, beta, bs41, alc
182 : REAL :: xpx, xfx, xlx, xpx0, xfx0, xlx0, ecp, ecf, fs, ec
183 : INTEGER :: ir
184 :
185 : ! .. Local Arrays ..
186 162 : REAL, ALLOCATABLE :: phi(:) ! relativistic exchange energy correct.
187 :
188 162 : thfpi = three / ( four * pi_const )
189 :
190 486 : ALLOCATE ( phi(ngrid) )
191 : CALL relcor( &
192 : mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
193 162 : phi)
194 :
195 162 : IF ( jspins == 2 ) THEN ! spinpolarized calculation
196 :
197 : c_1 = one / ( two**fothrd - two )
198 18237056 : DO ir = 1,ngrid
199 18236902 : rh1 = max(d_15,rh(ir,1))
200 18236902 : rh2 = max(d_15,rh(ir,jspins))
201 18236902 : rho = rh1 + rh2
202 18236902 : y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1
203 18236902 : fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
204 18236902 : rs = ( thfpi/rho )**thrd
205 18236902 : x = sqrt(rs)
206 18236902 : xpx = x*x + bp*x + cp
207 18236902 : xfx = x*x + bf*x + cf
208 18236902 : xlx = x*x + bl*x + cl
209 18236902 : xpx0 = xp0*xp0 + bp*xp0 + cp
210 18236902 : xfx0 = xf0*xf0 + bf*xf0 + cf
211 18236902 : xlx0 = xl0*xl0 + bl*xl0 + cl
212 18236902 : ecp = fec(ap,x,xpx,xp0,xpx0,bp,qp) ! paramagnetic correlation energy
213 18236902 : ecf = fec(af,x,xfx,xf0,xfx0,bf,qf) ! ferromagnetic correlation ener.
214 18236902 : alc = fec(al,x,xlx,xl0,xlx0,bl,ql) ! alpha_c
215 18236902 : beta = fdd0* (ecf-ecp)/alc - one ! beta = (f"(0)*(ecf-ecp))/alc -1
216 18236902 : bs41 = one + beta * s**4
217 :
218 18236902 : ec = ecp + alc * fs / fdd0 * bs41 ! total correlation energy
219 18236902 : ex = -cex/rs* (one + 0.2599210482 * fs) ! ex = exp + (exf-exp)*f(s)
220 : ! exf - exp = (2**(1/3)-1) * exp
221 : ! = 0.25992105 * exp
222 18237056 : exc(ir) = ec + ex*phi(ir)
223 : ENDDO
224 :
225 8 : ELSEIF ( jspins == 1 ) THEN ! non-spinpolarized calculation
226 :
227 716592 : DO ir = 1,ngrid
228 716584 : rho = max(d_15,rh(ir,1))
229 716584 : rs = ( thfpi/rho )**thrd
230 716584 : x = sqrt(rs)
231 716584 : xpx = x*x + bp*x + cp
232 716584 : xpx0 = xp0*xp0 + bp*xp0 + cp
233 716584 : ecp = fec(ap,x,xpx,xp0,xpx0,bp,qp) ! correlation energy paramagn.
234 716584 : ex = -cex / rs ! exp: paramag. exchange energy
235 : ! (like in wigner formula)
236 716592 : exc(ir) = ecp + ex*phi(ir)
237 : ENDDO
238 :
239 : ELSE
240 0 : WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
241 0 : CALL juDFT_error("excvwn",calledby="xcvwn")
242 : ENDIF
243 :
244 162 : DEALLOCATE (phi)
245 162 : RETURN
246 :
247 : END SUBROUTINE excvwn
248 :
249 : !--------------------------------------------------------------------
250 112620119 : REAL FUNCTION fec(ai,z,ziz,zi0,ziz0,bi,qi)
251 : REAL :: ai,z,ziz,zi0,ziz0,bi,qi
252 : fec = ai* (alog(z*z/ziz)+ two*bi/qi * atan(qi/(two*z+bi)) - &
253 : bi*zi0/ziz0* ( alog((z-zi0)**2/ziz) + &
254 112620119 : two*(bi+ two*zi0)/qi*atan(qi/ (two*z+bi))) )
255 112620119 : END FUNCTION fec
256 19590751 : REAL FUNCTION fdedr(rho,ai,z,ziz,zi0,bi,ci)
257 : REAL :: rho,ai,z,ziz,zi0,bi,ci
258 19590751 : fdedr = -ai*z/(three*rho*ziz)*(ci/z-bi*zi0/(z-zi0))
259 0 : END FUNCTION fdedr
260 : !--------------------------------------------------------------------
261 :
262 : END MODULE m_xcvwn
|
