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