Line data Source code
1 : MODULE m_easypbe
2 : !----------------------------------------------------------------------
3 : ! easypbe---exchpbe
4 : ! --corpbe --- pbecor2
5 : !----------------------------------------------------------------------
6 : ! easypbe is a driver for the pbe subroutines, using simple inputs
7 : ! k. burke, may 14, 1996.
8 : !----------------------------------------------------------------------
9 : CONTAINS
10 147738515 : SUBROUTINE easypbe(xcpot, &
11 : up,agrup,delgrup,uplap, &
12 : dn,agrdn,delgrdn,dnlap, &
13 : agr,delgr,lcor,lpot, &
14 : exlsd,vxuplsd,vxdnlsd, &
15 : eclsd,vcuplsd,vcdnlsd, &
16 : expbe,vxuppbe,vxdnpbe, &
17 : ecpbe,vcuppbe,vcdnpbe, &
18 : vxupsr,vxdnsr)
19 :
20 : USE m_exchpbe
21 : USE m_corpbe
22 : USE m_types_xcpot_data
23 : IMPLICIT NONE
24 :
25 : ! .. Arguments ..
26 : TYPE(t_xcpot_data),INTENT(IN)::xcpot
27 : INTEGER, INTENT (IN) :: lcor ! flag to do correlation(=0=>don't)
28 : INTEGER, INTENT (IN) :: lpot ! flag to do potential (=0=>don't)
29 : REAL, INTENT (IN) :: up,dn ! density (spin up & down)
30 : REAL, INTENT (IN) :: agrup,agrdn ! |grad up|, |grad down|
31 : REAL, INTENT (IN) :: delgrup,delgrdn ! delgrup=(grad up).(grad |grad up|) (in pw91: gggru)
32 : REAL, INTENT (IN) :: uplap,dnlap ! grad^2 up=laplacian of up (in pw91: g2ru)
33 : REAL, INTENT (IN) :: agr,delgr ! agr=|grad rho|,
34 : ! delgr=(grad rho).(grad |grad rho|) (in pw91: gggr)
35 : REAL, INTENT (OUT) :: vxuplsd,vxdnlsd ! up/down lsd exchange potential
36 : REAL, INTENT (OUT) :: vcuplsd,vcdnlsd ! up/down lsd correlation potential
37 : REAL, INTENT (OUT) :: exlsd,eclsd ! lsd exchange / correlation energy density
38 : REAL, INTENT (OUT) :: vxuppbe,vxdnpbe ! as above, but pbe quantities
39 : REAL, INTENT (OUT) :: vcuppbe,vcdnpbe
40 : REAL, INTENT (OUT) :: expbe,ecpbe ! note that : exlsd=int d^3r rho(r) exlsd(r)
41 : REAL, INTENT (OUT) :: vxupsr,vxdnsr
42 :
43 : ! .. local variables ..
44 : ! for exchange:
45 : REAL :: fk ! local fermi wavevector for 2*up=(3 pi^2 (2up))^(1/3)
46 : REAL :: s ! dimensionless density gradient=|grad rho|/ (2*fk*rho)_(rho=2*up)
47 : REAL :: u ! delgrad/(rho^2*(2*fk)**3)_(rho=2*up)
48 : REAL :: v ! laplacian/(rho*(2*fk)**2)_(rho=2*up)
49 : ! for correlation:
50 : REAL :: zet ! (up-dn)/rho
51 : REAL :: g ! phi(zeta)
52 : REAL :: rs ! (3/(4pi*rho))^(1/3)=local seitz radius=alpha/fk
53 : REAL :: sk ! ks=thomas-fermi screening wavevector=sqrt(4fk/pi)
54 : REAL :: twoksg ! 2*ks*phi
55 : REAL :: t ! correlation dimensionless gradient=|grad rho|/(2*ks*phi*rho)
56 : REAL :: uu ! delgrad/(rho^2*twoksg^3)
57 : REAL :: rholap ! laplacian
58 : REAL :: vv ! laplacian/(rho*twoksg^2)
59 : REAL :: ww ! (|grad up|^2-|grad dn|^2-zet*|grad rho|^2)/(rho*twoksg)^2
60 : REAL :: ec ! lsd correlation energy
61 : REAL :: vcup ! lsd up correlation potential
62 : REAL :: vcdn ! lsd down correlation potential
63 : REAL :: h ! gradient correction to correlation energy
64 : REAL :: dvcup ! gradient correction to up correlation potential
65 : REAL :: dvcdn ! gradient correction to down correlation potential
66 :
67 : REAL :: rdum,eta,exdnlsd,exdnpbe,exuplsd,exuppbe,rho,rho2
68 : ! ..
69 : REAL, PARAMETER :: thrd = 1.e0/3.e0
70 : REAL, PARAMETER :: thrd2 = 2.e0*thrd
71 : REAL, PARAMETER :: pi = 3.1415926535897932384626433832795e0
72 : REAL, PARAMETER :: pi32 = 29.608813203268075856503472999628e0 ! 3 pi**2
73 : REAL, PARAMETER :: alpha=1.91915829267751300662482032624669e0 ! (9pi/4)**thrd
74 :
75 147738515 : rho2 = 2.e0*up
76 :
77 : !----------------------------------------------------------------------
78 : ! pbe exchange
79 : ! use ex[up,dn]=0.5*(ex[2*up]+ex[2*dn]) (i.e., exact spin-scaling)
80 : ! do up exchange
81 : !----------------------------------------------------------------------
82 :
83 147738515 : IF (rho2 > 1e-18) THEN
84 147738515 : fk = (pi32*rho2)**thrd
85 147738515 : s = 2.e0*agrup/ (2.e0*fk*rho2)
86 147738515 : u = 4.e0*delgrup/ (rho2*rho2* (2.e0*fk)**3)
87 147738515 : v = 2.e0*uplap/ (rho2* (2.e0*fk)**2)
88 :
89 147738515 : CALL exchpbe(xcpot,rho2,s,u,v,0,lpot,exuplsd,vxuplsd,rdum)
90 147738515 : CALL exchpbe(xcpot,rho2,s,u,v,1,lpot,exuppbe,vxuppbe,vxupsr)
91 :
92 : ELSE
93 :
94 0 : exuplsd = 0.e0
95 0 : vxuplsd = 0.e0
96 0 : exuppbe = 0.e0
97 0 : vxuppbe = 0.e0
98 0 : vxupsr = 0.e0
99 :
100 : ENDIF
101 :
102 : ! repeat for down
103 147738515 : rho2 = 2.e0*dn
104 :
105 147738515 : IF (rho2 > 1e-18) THEN
106 :
107 147738515 : fk = (pi32*rho2)**thrd
108 147738515 : s = 2.e0*agrdn/ (2.e0*fk*rho2)
109 147738515 : u = 4.e0*delgrdn/ (rho2*rho2* (2.e0*fk)**3)
110 147738515 : v = 2.e0*dnlap/ (rho2* (2.e0*fk)**2)
111 :
112 147738515 : CALL exchpbe(xcpot,rho2,s,u,v,0,lpot,exdnlsd,vxdnlsd,rdum)
113 147738515 : CALL exchpbe(xcpot,rho2,s,u,v,1,lpot,exdnpbe,vxdnpbe,vxdnsr)
114 :
115 : ELSE
116 :
117 0 : exdnlsd = 0.e0
118 0 : vxdnlsd = 0.e0
119 0 : exdnpbe = 0.e0
120 0 : vxdnpbe = 0.e0
121 0 : vxdnsr = 0.e0
122 :
123 : ENDIF
124 :
125 : ! construct total density and contribution to ex
126 147738515 : rho = up + dn
127 147738515 : exlsd = (exuplsd*up+exdnlsd*dn)/rho
128 147738515 : expbe = (exuppbe*up+exdnpbe*dn)/rho
129 :
130 147738515 : IF (lcor == 0) RETURN
131 : !----------------------------------------------------------------------
132 : ! now do correlation
133 : !----------------------------------------------------------------------
134 :
135 147738515 : IF (rho < 1.e-18) RETURN
136 :
137 147738515 : zet = (up-dn)/rho
138 :
139 : ! 999+
140 : ! eta: eta should not be smaller than 1.e-16.
141 : ! otherwise will cause floating invalid.
142 : ! if bigger, the last digit may differ
143 : ! from the run by aix without this zet-guard.
144 147738515 : eta = 1.e-16
145 147738515 : zet = min(zet,1.0-eta)
146 147738515 : zet = max(zet,-1.0+eta)
147 : ! 999-
148 :
149 147738515 : g = ((1.e0+zet)**thrd2+ (1.e0-zet)**thrd2)/2.e0
150 147738515 : fk = (pi32*rho)**thrd
151 147738515 : rs = alpha/fk
152 147738515 : sk = sqrt(4.e0*fk/pi)
153 147738515 : twoksg = 2.e0*sk*g
154 147738515 : t = agr/ (twoksg*rho)
155 147738515 : uu = delgr/ (rho*rho*twoksg**3)
156 147738515 : rholap = uplap + dnlap
157 147738515 : vv = rholap/ (rho*twoksg**2)
158 147738515 : ww = (agrup**2-agrdn**2-zet*agr**2)/ (rho*rho*twoksg**2)
159 :
160 : CALL corpbe( &
161 : xcpot%is_PBEs,rs,zet,t,uu,vv,ww,1,lpot, &
162 147738515 : ec,vcup,vcdn,h,dvcup,dvcdn)
163 :
164 147738515 : eclsd = ec
165 147738515 : ecpbe = ec + h
166 147738515 : vcuplsd = vcup
167 147738515 : vcdnlsd = vcdn
168 147738515 : vcuppbe = vcup + dvcup
169 147738515 : vcdnpbe = vcdn + dvcdn
170 :
171 : END SUBROUTINE easypbe
172 : END MODULE m_easypbe
|
