The Elk Code
 
Loading...
Searching...
No Matches
fxc_pwca.f90
Go to the documentation of this file.
1
2! Copyright (C) 2011 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
3! This file is distributed under the terms of the GNU Lesser General Public
4! License. See the file COPYING for license details.
5
6subroutine fxc_pwca(n,rhoup,rhodn,fxcuu,fxcud,fxcdd)
7implicit none
8! arguments
9integer, intent(in) :: n
10real(8), intent(in) :: rhoup(n)
11real(8), intent(in) :: rhodn(n)
12real(8), intent(out) :: fxcuu(n)
13real(8), intent(out) :: fxcud(n)
14real(8), intent(out) :: fxcdd(n)
15! local variables
16integer i
17real(8), parameter :: pi=3.1415926535897932385d0
18real(8), parameter :: thrd=1.d0/3.d0, thrd4=4.d0/3.d0
19real(8), parameter :: d2f0=1.709921d0
20real(8), parameter :: a(3)=[ 0.0310907d0, 0.01554535d0, 0.0168869d0 ]
21real(8), parameter :: a1(3)=[ 0.21370d0, 0.20548d0, 0.11125d0 ]
22real(8), parameter :: b1(3)=[ 7.5957d0, 14.1189d0, 10.357d0 ]
23real(8), parameter :: b2(3)=[ 3.5876d0, 6.1977d0, 3.6231d0 ]
24real(8), parameter :: b3(3)=[ 1.6382d0, 3.3662d0, 0.88026d0 ]
25real(8), parameter :: b4(3)=[ 0.49294d0, 0.62517d0, 0.49671d0 ]
26real(8) p1,p2,p3,rup,rdn,r,ri,ri2,ri3
27real(8) rs,rs2,rs12,rs32,rsi,rs12i,rs32i
28real(8) mz,z,z2,z3,z4,fz,dfz,d2fz
29real(8) drs,d2rs,dzu,d2zu,dzd,d2zd,d2zud
30real(8) ders,d2ers,dez,d2ez,d2ersz
31real(8) deu,d2eu,ded,d2ed,d2eud,ex
32real(8) ec0,dec0,d2ec0,ec1,dec1,d2ec1
33real(8) ac,dac,d2ac,a2,dt1,d2t1,dt2,d2t2
34real(8) t1,t2,t3,t4,t5,t6,t7,t8
35if (n < 1) then
36 write(*,*)
37 write(*,'("Error(fxc_pwca): n < 1 : ",I8)') n
38 write(*,*)
39 stop
40end if
41! prefactors
42t1=3.d0/(4.d0*pi)
43p1=t1**thrd
44p2=t1*(9.d0*pi/4.d0)**thrd
45p3=1.d0/(2.d0**thrd4-2.d0)
46do i=1,n
47 rup=rhoup(i); rdn=rhodn(i)
48! total density
49 r=rup+rdn
50 if ((rup < 0.d0).or.(rdn < 0.d0).or.(r < 1.d-20)) then
51 fxcuu(i)=0.d0
52 fxcud(i)=0.d0
53 fxcdd(i)=0.d0
54 cycle
55 end if
56 ri=1.d0/r
57 ri2=ri**2
58 ri3=ri2*ri
59 rs=p1*ri**thrd
60 rs2=rs**2
61 rs12=sqrt(rs)
62 rs32=rs12*rs
63 rsi=1.d0/rs
64 rs12i=1.d0/rs12
65 rs32i=1.d0/rs32
66 mz=rup-rdn
67 z=mz/r
68 z2=z**2
69 z3=z2*z
70 z4=z3*z
71! drs/drup = drs/drdn = drs/drho
72 drs=-thrd*rs*ri
73! d2rs/drup^2 = d^2rs/drn^2 = d^2rs/drho^2
74 d2rs=-thrd4*drs*ri
75! dz/drup, dz/drdn
76 t1=mz*ri2
77 dzu=ri-t1
78 dzd=-ri-t1
79! d^2z/drup^2, d^2z/drdn^2, d^2z/drup*drdn
80 t1=2.d0*mz*ri3
81 t2=2.d0*ri2
82 d2zu=t1-t2
83 d2zd=t1+t2
84 d2zud=t1
85!------------------!
86! exchange !
87!------------------!
88 t1=-p2*rsi/2.d0
89 t2=1.d0+z
90 t3=1.d0-z
91 t4=t2**thrd4
92 t5=t3**thrd4
93 t6=t4+t5
94! exchange energy density
95 ex=t1*t6
96! dex/drs
97 ders=-ex*rsi
98! d^2ex/drs^2
99 d2ers=-2.d0*ders*rsi
100! f(z)
101 fz=p3*(t6-2.d0)
102! dex/dz
103 t4=t4/t2
104 t5=t5/t3
105 t6=t4-t5
106 t7=thrd4*t6
107 dez=t1*t7
108! df/dz
109 dfz=p3*t7
110! d^2ex/dz^2
111 t4=t4/t2
112 t5=t5/t3
113 t6=t4+t5
114 t7=thrd4*thrd*t6
115 d2ez=t1*t7
116! d^2f/dz^2
117 d2fz=p3*t7
118! d^2ex/drs*dz
119 d2ersz=-dez*rsi
120! dex/drup, dex/drdn
121 t1=ders*drs
122 deu=t1+dez*dzu
123 ded=t1+dez*dzd
124! d^2ex/drup^2
125 t1=d2ers*drs
126 t2=d2ersz*drs
127 t3=ders*d2rs
128 t4=(t1+d2ersz*dzu)*drs+t3
129 t5=t2+d2ez*dzu
130 d2eu=t4+t5*dzu+dez*d2zu
131! d^2ex/drdn^2
132 d2ed=(t1+d2ersz*dzd)*drs+t3+(t2+d2ez*dzd)*dzd+dez*d2zd
133! d^2ex/drup*drdn
134 d2eud=t4+t5*dzd+dez*d2zud
135! calculate fxc
136 fxcuu(i)=2.d0*deu+r*d2eu
137 fxcud(i)=deu+ded+r*d2eud
138 fxcdd(i)=2.d0*ded+r*d2ed
139!---------------------!
140! correlation !
141!---------------------!
142! ec(rs,0)
143 a2=2.d0*a(1)
144 t1=a2*(b1(1)*rs12+b2(1)*rs+b3(1)*rs32+b4(1)*rs2)
145 dt1=a2*(0.5d0*b1(1)*rs12i+b2(1)+1.5d0*b3(1)*rs12+2.d0*b4(1)*rs)
146 d2t1=a2*(-0.25d0*b1(1)*rs32i+0.75d0*b3(1)*rs12i+2.d0*b4(1))
147 t3=1.d0/t1
148 t4=t3**2
149 t2=1.d0+t3
150 dt2=-dt1*t4
151 d2t2=t4*(2.d0*t3*dt1**2-d2t1)
152 t3=1.d0/t2
153 t4=1.d0+a1(1)*rs
154 t5=log(t2)
155 ec0=-a2*t4*t5
156 dec0=-a2*(a1(1)*t5+t4*t3*dt2)
157 d2ec0=-a2*(2.d0*a1(1)*t3*dt2+t4*t3*(d2t2-t3*dt2**2))
158! ec(rs,1)
159 a2=2.d0*a(2)
160 t1=a2*(b1(2)*rs12+b2(2)*rs+b3(2)*rs32+b4(2)*rs2)
161 dt1=a2*(0.5d0*b1(2)*rs12i+b2(2)+1.5d0*b3(2)*rs12+2.d0*b4(2)*rs)
162 d2t1=a2*(-0.25d0*b1(2)*rs32i+0.75d0*b3(2)*rs12i+2.d0*b4(2))
163 t3=1.d0/t1
164 t4=t3**2
165 t2=1.d0+t3
166 dt2=-dt1*t4
167 d2t2=t4*(2.d0*t3*dt1**2-d2t1)
168 t3=1.d0/t2
169 t4=1.d0+a1(2)*rs
170 t5=log(t2)
171 ec1=-a2*t4*t5
172 dec1=-a2*(a1(2)*t5+t4*t3*dt2)
173 d2ec1=-a2*(2.d0*a1(2)*t3*dt2+t4*t3*(d2t2-t3*dt2**2))
174! ac(rs)
175 a2=2.d0*a(3)
176 t1=a2*(b1(3)*rs12+b2(3)*rs+b3(3)*rs32+b4(3)*rs2)
177 dt1=a2*(0.5d0*b1(3)*rs12i+b2(3)+1.5d0*b3(3)*rs12+2.d0*b4(3)*rs)
178 d2t1=a2*(-0.25d0*b1(3)*rs32i+0.75d0*b3(3)*rs12i+2.d0*b4(3))
179 t3=1.d0/t1
180 t4=t3**2
181 t2=1.d0+t3
182 dt2=-dt1*t4
183 d2t2=t4*(2.d0*t3*dt1**2-d2t1)
184 t3=1.d0/t2
185 t4=1.d0+a1(3)*rs
186 t5=log(t2)
187 ac=a2*t4*t5
188 dac=a2*(a1(3)*t5+t4*t3*dt2)
189 d2ac=a2*(2.d0*a1(3)*t3*dt2+t4*t3*(d2t2-t3*dt2**2))
190! correlation energy density derivatives
191 t1=1.d0-z4
192 t2=(fz/d2f0)*t1
193 t3=ec1-ec0
194 t4=fz*z4
195! dec/drs
196 t5=dec1-dec0
197 ders=dec0+dac*t2+t5*t4
198! d^2ec/drs^2
199 t6=d2ec1-d2ec0
200 d2ers=d2ec0+d2ac*t2+t6*t4
201! dec/dz
202 t4=ac/d2f0
203 t6=4.d0*fz*z3
204 t7=dfz*t1-t6
205 t8=dfz*z4+t6
206 dez=t4*t7+t3*t8
207! d^2ec/drs*dz
208 d2ersz=(dac/d2f0)*t7+t5*t8
209! d^2ec/dz^2
210 t7=8.d0*dfz*z3
211 t8=12.d0*fz*z2
212 d2ez=t4*(d2fz*t1-t7-t8)+t3*(d2fz*z4+t7+t8)
213! dec/drup, dec/drdn
214 t1=ders*drs
215 deu=t1+dez*dzu
216 ded=t1+dez*dzd
217! d^2ec/drup^2
218 t1=d2ers*drs
219 t2=d2ersz*drs
220 t3=ders*d2rs
221 t4=(t1+d2ersz*dzu)*drs+t3
222 t5=t2+d2ez*dzu
223 d2eu=t4+t5*dzu+dez*d2zu
224! d^2ec/drdn^2
225 d2ed=(t1+d2ersz*dzd)*drs+t3+(t2+d2ez*dzd)*dzd+dez*d2zd
226! d^2ec/drup*drdn
227 d2eud=t4+t5*dzd+dez*d2zud
228! calculate fxc
229 fxcuu(i)=fxcuu(i)+2.d0*deu+r*d2eu
230 fxcud(i)=fxcud(i)+deu+ded+r*d2eud
231 fxcdd(i)=fxcdd(i)+2.d0*ded+r*d2ed
232end do
233end subroutine
234
fxc_pwca
subroutine fxc_pwca(n, rhoup, rhodn, fxcuu, fxcud, fxcdd)
Definition fxc_pwca.f90:7