Line data Source code
1 : MODULE m_vxcwb91
2 : !.....-----------------------------------------------------------------
3 : !.....pw91 exchange-correlation potential in hartree.
4 : !.....------------------------------------------------------------------
5 : CONTAINS
6 0 : SUBROUTINE vxcwb91( &
7 0 : jspins,mirm,irmx,rh,agr,agru,agrd, &
8 0 : g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr, &
9 0 : vx,vxc, &
10 : idsprs,isprsv,sprsv)
11 :
12 : USE m_corl91
13 : USE m_corg91
14 : USE m_xch91
15 : USE m_constants, ONLY: pi_const
16 :
17 : IMPLICIT NONE
18 :
19 : ! .. Arguments ..
20 : INTEGER, INTENT (IN) :: jspins,irmx,mirm,isprsv,idsprs
21 : REAL, INTENT (IN) :: sprsv
22 :
23 : REAL, INTENT (IN) :: rh(mirm,jspins)
24 : REAL, INTENT (IN) :: agr(mirm),agru(mirm),agrd(mirm)
25 : REAL, INTENT (IN) :: g2r(mirm),g2ru(mirm),g2rd(mirm)
26 : REAL, INTENT (IN) :: gggr(mirm),gggru(mirm)
27 : REAL, INTENT (IN) :: gggrd(mirm),gzgr(mirm)
28 : REAL, INTENT (OUT) ::vx(mirm,jspins),vxc(mirm,jspins)
29 :
30 : ! .. local variables
31 : REAL :: pi,ro,zta,alf,alfc,c13,c23,c43,c53, &
32 : cedg,cedl,dbrod,dbrou,dsprs,ec,ecrs,eczta,fk,gz, &
33 : ro13,ro2,rod,rod3,rod43,rod53,rou,rou3,rou43,rou53, &
34 : rs,sd,sk,su,tc,td,tksg,tu,uc,ud,uu,vc, &
35 : vcgd,vcgu,vcld,vclu,vxgd,vxgu,vxld, &
36 : vxlu,wc,xced,xedg,xedgd,xedgu,xedl,xedld,xedlu, &
37 : rou13,rod13,xcptu,xcptd
38 :
39 : INTEGER :: i
40 : REAL, PARAMETER :: sml = 1.e-14
41 : REAL, PARAMETER :: huges = 1.e+6
42 : !.....-----------------------------------------------------------------
43 : !.....huges: if su(sdw).gt.huges, grad.expansion calc.non-excuted.
44 : !.....-----------------------------------------------------------------
45 :
46 0 : do 300 i = 1,irmx
47 :
48 0 : if(jspins == 1) then
49 0 : rou=rh(i,1)/2
50 0 : rou=max(rou,sml)
51 0 : rod=rou
52 : else
53 0 : rou=rh(i,1)
54 0 : rod=rh(i,jspins)
55 0 : rou=max(rou,sml)
56 0 : rod=max(rod,sml)
57 : endif
58 :
59 0 : zta=(rou-rod)/(rou+rod)
60 :
61 0 : if(zta > 1.e0-sml) zta = 1.e0 - sml
62 0 : if(zta < -1.e0-sml) zta = -1.e0 + sml
63 :
64 : !.....
65 : ! vxlu,vxld,vxgu,vxgd: exchange potential in ry.(local,grad),
66 : !c (up,dw).
67 : ! vclu,vcld,vcgu,vcgd: correl. potential in ry.(local,grad),
68 : !c (up,dw).
69 : ! all later in hartree.
70 : !.....
71 0 : vxlu = 0.0e0
72 0 : vclu = 0.0e0
73 0 : vxld = 0.0e0
74 0 : vcld = 0.0e0
75 0 : vxgu = 0.0e0
76 0 : vcgu = 0.0e0
77 0 : vxgd = 0.0e0
78 0 : vcgd = 0.0e0
79 :
80 : !.....
81 0 : ro = rh(i,1)
82 0 : IF (ro < sml) THEN
83 0 : ro = sml
84 : zta = 0.e0
85 : GO TO 200
86 : ENDIF
87 : !.....
88 0 : c13 = 1.e0/3.e0
89 0 : c23 = 2.e0/3.e0
90 0 : c43 = 4.e0/3.e0
91 0 : c53 = 5.e0/3.e0
92 : !..... alf=-3*(3/4*pai)**(1/3).
93 0 : alf = -1.861051473e0
94 : !.....
95 0 : ro2 = ro*ro
96 0 : ro13 = ro**c13
97 : !.....
98 0 : rou = ro* (1.e0+zta)/2.e0
99 0 : rou3 = rou**3
100 0 : rou13 = rou**c13
101 : rou43 = rou**c43
102 : !.....
103 0 : rod = ro - rou
104 0 : rod3 = rod**3
105 0 : rod13 = rod**c13
106 : rod43 = rod**c43
107 : !.....
108 : ! gr2=drr*drr
109 : ! gr2u=drru**2
110 : ! drrd=drr-drru
111 : ! gr2d=drrd**2
112 : ! ddrrd=ddrr-ddrru
113 : !.....
114 : !..... gz,gz2,gz3: for wang-perdew ssf.
115 0 : gz = ((1.e0+zta)**c23+ (1.e0-zta)**c23)/2.e0
116 : !.....
117 : !.....
118 0 : rs = 0.620350491e0/ro13
119 : !.....
120 : !..... exchange-potential, vxp,vxlu,vxld: v-exchange-(para,up,dw).
121 0 : vxlu = c43*alf*rou13
122 0 : vxld = c43*alf*rod13
123 :
124 : !.....
125 : !.... .gradient correction.
126 : !.....
127 0 : if(abs(agr(i)) < sml) go to 200
128 : !.....
129 :
130 : !.....
131 0 : dsprs = 1.e0
132 0 : if(idsprs == 1) dsprs = 1.e-19
133 0 : if(isprsv == 1) dsprs = dsprs*sprsv
134 :
135 : !.....
136 : ! agr,agru,agrd: abs(grad(rho)), for all, up, and down.
137 : !c gr2,gr2u,gr2d: grad(rho_all)**2, grad(rho_up)**2, grad(rho_d)**2.
138 : ! g2r,g2ru,g2rd: laplacian rho_all, _up and _down.
139 : ! gggru,-d: grad(rho)*grad(abs(grad(rho))) for all,up and down.
140 : ! grgru,-d: grad(rho_all)*grad(rhor_up) and for down.
141 :
142 : ! g2r=ddrr+2*drr/rv
143 : !.....
144 : rou53 = rou**c53
145 : !.....
146 : !..... edrru: d(abs(d(rou)/dr))/dr, edrrd for down.
147 : ! edrru=ddrru
148 : ! if(drru.lt.0.) edrru=-ddrru
149 : !.....
150 : ! agr,agbru,-d: abs(grad(rho)),for rou, rod.
151 : ! gggru,-d: grad(rho)*grad(abs(grad(rho))) for up and down.
152 : !..... su:at ro=2*rou. 1/(2(3*pai**2)**(1/3))*|grad(rou)|/rou**(4/3).
153 0 : su = 0.128278244e0*agru(i)/rou43
154 0 : if(su > huges) go to 200
155 : ! g2ru=ddrru+2*drru/rv
156 0 : tu = .016455307e0*g2ru(i)/rou53
157 0 : uu = 0.002110857e0*gggru(i)/rou3
158 :
159 0 : dbrou = rou*2
160 :
161 0 : call xch91(dbrou,su,uu,tu,xedlu,xedgu,vxlu,vxgu)
162 :
163 0 : vxgu = dsprs*vxgu
164 :
165 : !.....
166 : rod53 = rod**c53
167 : ! edrrd=ddrrd
168 : ! if(drrd.lt.0.) edrrd=-ddrrd
169 :
170 0 : sd = 0.128278244e0*agrd(i)/rod43
171 0 : if(sd > huges) go to 200
172 :
173 : ! g2rd=ddrrd+2*drrd/rv
174 :
175 0 : td = .016455307e0*g2rd(i)/rod53
176 0 : ud = 0.002110857e0*gggrd(i)/rod3
177 :
178 0 : dbrod = rod*2
179 :
180 0 : call xch91(dbrod,sd,ud,td,xedld,xedgd,vxld,vxgd)
181 :
182 0 : vxgd = dsprs*vxgd
183 :
184 : !.... cro: c(n) of (6),phys.rev..b33,8822('86). in ry.
185 : !.... dcdr: d(cro)/d(ro).
186 : !..... 0.001625816=1.745*f(=0.11)*cro(rs=0).
187 :
188 : ! pw91
189 :
190 0 : call corl91(rs,zta,ec,vclu,vcld,ecrs,eczta,alfc)
191 :
192 0 : vclu = vclu*2.e0
193 0 : vcld = vcld*2.e0
194 :
195 0 : fk = 1.91915829e0/rs
196 0 : sk = sqrt(4.e0*fk/pi_const)
197 0 : tksg = 2.e0*sk*gz
198 0 : tc = agr(i)/ (ro*tksg)
199 0 : uc = gggr(i)/ (ro2*tksg**3)
200 0 : vc = g2r(i)/ (ro*tksg**2)
201 0 : wc = gzgr(i)/ (ro*tksg**2)
202 :
203 : call corg91(fk,sk,gz,ec,ecrs,eczta,rs,zta,tc,uc,vc,wc,cedg, &
204 0 : vcgu,vcgd)
205 :
206 0 : vcgu = dsprs*vcgu*2.e0
207 0 : vcgd = dsprs*vcgd*2.e0
208 :
209 : 200 continue
210 :
211 0 : xcptu = vxlu + vclu + vxgu + vcgu
212 0 : xcptd = vxld + vcld + vxgd + vcgd
213 :
214 : ! change to hartree.
215 0 : vx(i,1) = vxlu + vxgu
216 0 : vx(i,jspins) = vxld + vxgd
217 :
218 0 : vxc(i,1) = xcptu
219 0 : vxc(i,jspins) = xcptd
220 :
221 0 : 300 END DO
222 :
223 0 : END SUBROUTINE vxcwb91
224 : END MODULE m_vxcwb91
|
