Line data Source code
1 : MODULE m_corg91
2 : !.....-----------------------------------------------------------------
3 : ! gga91 correlation
4 : !.....-----------------------------------------------------------------
5 : CONTAINS
6 0 : SUBROUTINE corg91(fk,sk,gz,ec,ecrs,eczta,rs,zta,t,uu,vv,ww,h, &
7 : dvcup,dvcdn)
8 : !.....-----------------------------------------------------------------
9 : ! input
10 : ! rs: seitz radius
11 : ! zta: relative spin polarization
12 : ! t: abs(grad d)/(d*2.*ks*gz)
13 : ! uu: (grad d)*grad(abs(grad d))/(d**2 * (2*ks*gz)**3)
14 : ! vv: (laplacian d)/(d * (2*ks*gz)**2)
15 : ! ww: (grad d)*(gradzta)/(d * (2*ks*gz)**2
16 : ! output
17 : ! h: nonlocal part of correlation energy per electron
18 : ! dvcup,-dn: nonlocal parts of correlation potentials.
19 :
20 : ! with ks=sqrt(4*kf/pai), gz=[(1+zta)**(2/3)+(1-zta)**(2/3)]/2, &
21 : ! kf=cbrt(3*pai**2*d).
22 : !.....-----------------------------------------------------------------
23 : !.....-----------------------------------------------------------------
24 : ! .. previously untyped names ..
25 : IMPLICIT NONE
26 :
27 : REAL :: dvcdn,dvcup,ec,ecrs,eczta,fk,gz,h,rs,sk,t,uu,vv,ww,zta,a4, &
28 : alf,argmx,b,b2,b2fac,bec,bet,bg,c1,c2,c3,c4,c5,c6,cc,cc0, &
29 : ccrs,coeff,comm,cx,delt,fact0,fact1,fact2,fact3,fact4,fact5, &
30 : gm,gz3,gz4,h0,h0b,h0bt,h0rs,h0rst,h0t,h0tt,h0z,h0zt,h1,h1rs, &
31 : h1rst,h1t,h1tt,h1z,h1zt,hrs,hrst,ht,htt,hz,hzt,pon,pref,q4, &
32 : q5,q6,q7,q8,q9,r0,r1,r2,r3,r4,rs2,rs3,rsthrd,t2,t4,t6,thrd2, &
33 : thrdm,xnu,r1t2,expsm
34 : ! ..
35 : DATA xnu,cc0,cx,alf/15.75592e0,0.00423500,-0.001667212e0,0.0900/
36 : DATA c1,c2,c3,c4/0.00256800,0.02326600,7.389e-6,8.723e0/
37 : DATA c5,c6,a4/0.472e0,7.389e-2,100.00/
38 : DATA thrdm,thrd2/-0.333333333333e0,0.666666666667e0/
39 : !.....-----------------------------------------------------------------
40 : ! expsm: argument of exponential-smallest.
41 0 : expsm=minexponent(expsm)/1.5
42 0 : argmx = 174.0
43 0 : bet = xnu*cc0
44 0 : delt = 2.e0*alf/bet
45 0 : gz3 = gz**3
46 0 : gz4 = gz3*gz
47 0 : pon = -delt*ec/ (gz3*bet)
48 0 : IF (pon > argmx) THEN
49 : b = 0.
50 : ELSE
51 0 : b = delt/ (exp(pon)-1.00)
52 : ENDIF
53 0 : b2 = b*b
54 0 : t2 = t*t
55 0 : t4 = t2*t2
56 0 : t6 = t4*t2
57 0 : rs2 = rs*rs
58 0 : rs3 = rs2*rs
59 0 : q4 = 1.00 + b*t2
60 0 : q5 = 1.00 + b*t2 + b2*t4
61 0 : q6 = c1 + c2*rs + c3*rs2
62 0 : q7 = 1.00 + c4*rs + c5*rs2 + c6*rs3
63 0 : cc = -cx + q6/q7
64 0 : r0 = (sk/fk)**2
65 0 : r1 = a4*r0*gz4
66 0 : coeff = cc - cc0 - 3.e0*cx/7.00
67 0 : r2 = xnu*coeff*gz3
68 : ! tagu
69 0 : r1t2=max(-r1*t2,expsm)
70 0 : r3 = exp(r1t2)
71 : ! tagu
72 0 : h0 = gz3* (bet/delt)*log(1.00+delt*q4*t2/q5)
73 0 : h1 = r3*r2*t2
74 0 : h = h0 + h1
75 : ! local correlation option:
76 : ! h=0.0e0
77 :
78 : ! energy done. now the potential:
79 :
80 0 : ccrs = (c2+2.*c3*rs)/q7 - q6* (c4+2.*c5*rs+3.*c6*rs2)/q7**2
81 0 : rsthrd = rs/3.e0
82 0 : r4 = rsthrd*ccrs/coeff
83 0 : gm = ((1.00+zta)**thrdm- (1.00-zta)**thrdm)/3.00
84 0 : IF (pon > argmx) THEN
85 : b2fac = 0.
86 : ELSE
87 0 : b2fac = b2* (delt/b+1.00)
88 : ENDIF
89 0 : bg = -3.e0*ec*b2fac/ (bet*gz4)
90 0 : bec = b2fac/ (bet*gz3)
91 0 : q8 = q5*q5 + delt*q4*q5*t2
92 0 : q9 = 1.00 + 2.00*b*t2
93 0 : h0b = -bet*gz3*b*t6* (2.e0+b*t2)/q8
94 0 : h0rs = -rsthrd*h0b*bec*ecrs
95 0 : fact0 = 2.e0*delt - 6.00*b
96 0 : fact1 = q5*q9 + q4*q9*q9
97 0 : h0bt = 2.e0*bet*gz3*t4* ((q4*q5*fact0-delt*fact1)/q8)/q8
98 0 : h0rst = rsthrd*t2*h0bt*bec*ecrs
99 0 : h0z = 3.e0*gm*h0/gz + h0b* (bg*gm+bec*eczta)
100 0 : h0t = 2.*bet*gz3*q9/q8
101 0 : h0zt = 3.e0*gm*h0t/gz + h0bt* (bg*gm+bec*eczta)
102 0 : fact2 = q4*q5 + b*t2* (q4*q9+q5)
103 0 : fact3 = 2.e0*b*q5*q9 + delt*fact2
104 0 : h0tt = 4.e0*bet*gz3*t* (2.e0*b/q8- (q9* (fact3/q8))/q8)
105 0 : h1rs = r3*r2*t2* (-r4+r1*t2/3.e0)
106 0 : fact4 = 2.e0 - r1*t2
107 0 : h1rst = r3*r2*t2* (2.e0*r4* (1.00-r1*t2)-thrd2*r1*t2*fact4)
108 0 : h1z = gm*r3*r2*t2* (3.e0-4.e0*r1*t2)/gz
109 0 : h1t = 2.e0*r3*r2* (1.00-r1*t2)
110 0 : h1zt = 2.e0*gm*r3*r2* (3.e0-11.00*r1*t2+4.e0*r1*r1*t4)/gz
111 0 : h1tt = 4.e0*r3*r2*r1*t* (-2.e0+r1*t2)
112 0 : hrs = h0rs + h1rs
113 0 : hrst = h0rst + h1rst
114 0 : ht = h0t + h1t
115 0 : htt = h0tt + h1tt
116 0 : hz = h0z + h1z
117 0 : hzt = h0zt + h1zt
118 0 : comm = h + hrs + hrst + t2*ht/6.00 + 7.00*t2*t*htt/6.00
119 0 : pref = hz - gm*t2*ht/gz
120 0 : fact5 = gm* (2.e0*ht+t*htt)/gz
121 0 : comm = comm - pref*zta - uu*htt - vv*ht - ww* (hzt-fact5)
122 0 : dvcup = comm + pref
123 0 : dvcdn = comm - pref
124 :
125 : ! local correlation option:
126 : ! dvcup=0.0e0
127 : ! dvcdn=0.0e0
128 :
129 0 : END SUBROUTINE corg91
130 : END MODULE m_corg91
|
