1c
2c$Id: xc_xmvs15.F 21740 2012-01-11 00:25:15Z edo $
3c
4#include "dft2drv.fh"
5c     Made Very Simple (MVS) functional (Exchange part only)
6c           META GGA
7C         utilizes ingredients:
8c                              rho   -  density
9c                              delrho - gradient of density
10c                              tau - K.S kinetic energy density
11c                              tauW - von Weiszacker kinetic energy density
12c                              tauU - uniform-gas KE density
13c     References:
14c     J. Sun, J.P. Perdew, A. Ruzsinszky
15c     PNAS 2015, 112, 685-689
16c     DOI: 10.1073/pnas.1423145112
17
18      Subroutine xc_xmvs15(tol_rho, fac,  rho, delrho,
19     &                     Amat, Cmat, nq, ipol, Ex,
20     &                     qwght, ldew, func, tau,Mmat)
21      implicit none
22c
23      double precision fac, Ex
24      integer nq, ipol
25      logical ldew
26      double precision func(*)  ! value of the functional [output]
27c
28c     Charge Density & Its Cube Root
29c
30      double precision rho(nq,ipol*(ipol+1)/2)
31c
32c     Charge Density Gradient
33c
34      double precision delrho(nq,3,ipol)
35c
36c     Quadrature Weights
37c
38      double precision qwght(nq)
39c
40c     Sampling Matrices for the XC Potential & Energy
41c
42      double precision Amat(nq,ipol), Cmat(nq,*)
43c
44c     kinetic energy density   or  tau
45c
46      double precision tau(nq,ipol), Mmat(nq,*)
47      double precision tol_rho
48c
49      integer ispin,cmatpos
50c
51      if (ipol.eq.1 )then
52c
53c     SPIN-RESTRICTED
54c     Ex = Ex[n]
55c
56         call xc_xmvs15_cs(tol_rho, fac,  rho, delrho,
57     &                     Amat, Cmat, nq, Ex, 1d0,
58     &                     qwght, ldew, func, tau,Mmat)
59      else
60c
61c     SPIN-UNRESTRICTED
62c     Ex = Ex[2n_up]/2 + Ex[2n_down]/2
63
64         do ispin=1,2
65            if (ispin.eq.1) cmatpos=D1_GAA
66            if (ispin.eq.2) cmatpos=D1_GBB
67            call xc_xmvs15_cs(tol_rho, fac,
68     R           rho(1,ispin+1), delrho(1,1,ispin),
69     &           Amat(1,ispin), Cmat(1,cmatpos),
70     &           nq, Ex, 2d0,
71     &           qwght, ldew, func,
72     T           tau(1,ispin),Mmat(1,ispin))
73         enddo
74
75      endif
76      return
77      end
78      Subroutine xc_xmvs15_cs(tol_rho, fac,  rho, delrho,
79     &                     Amat, Cmat, nq, Ex, facttwo,
80     &                     qwght, ldew, func, tau,Mmat)
81      implicit none
82c
83      double precision fac, Ex
84      integer nq
85      logical ldew
86      double precision func(*)  ! value of the functional [output]
87c
88c     Charge Density & Its Cube Root
89c
90      double precision rho(*)
91c
92c     Charge Density Gradient
93c
94      double precision delrho(nq,3)
95c
96c     Quadrature Weights
97c
98      double precision qwght(nq)
99c
100c     Sampling Matrices for the XC Potential & Energy
101c
102      double precision Amat(nq), Cmat(nq)
103c
104c     kinetic energy density   or  tau
105c
106      double precision tau(nq,*), Mmat(nq)
107c
108      double precision facttwo ! 2 for o.s. 1 for c.s.
109c
110      integer n
111      double precision tol_rho, pi
112      double precision rrho, rho43, rho13, rho23, rho83
113      double precision tauN, tauW, tauU, gamma
114      double precision p, a, z, rz, g2, g, uge
115      double precision F83, F23, F53, F43, F13, F18
116      double precision afact2, Ax, rhoval, Pconst
117      double precision rK0, rB, rE1, rC1
118
119c     functional derivatives below FFFFFFFFFFFF
120
121      double precision derivr, derivg, derivt
122      double precision FaNum, BFaDen, FaDen, Fa
123      double precision dFaNumda, dFaDenda, dBFaDenda, dFada
124      double precision BFxDen, FxDen, FxNum
125      double precision dFxNumda, dFxDendp
126      double precision dadp, dadz, dpdg, dpdr, dzdr, dzdg, dzdt
127      double precision Fx, dFxdp, dFxdz, dFxdr, dFxdg, dFxdt
128
129c     functional derivatives above FFFFFFFFFFFF
130
131      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0)
132      parameter (F83=8.d0/3.d0, F23=2.d0/3.d0)
133      parameter (F18=1.d0/8.d0, F53=5.d0/3.d0)
134      parameter (rK0=0.174d0, rB=0.0233d0,
135     &           rE1=-1.66651d0, uge=10.d0/81.d0)
136c
137      pi=acos(-1d0)
138      Ax = (-0.75d0)*(3d0/pi)**F13
139      Pconst = (3.d0*pi**2)**F23
140      afact2=1d0/facttwo
141c
142      rC1=(20.d0*rK0/(27.d0*uge))**(4d0)-(1.d0+rE1)*(1.d0+rE1) ! gfortran error if this is defined as a parameter
143c
144      do n = 1, nq
145         if (rho(n).ge.tol_rho) then
146
147            rhoval=rho(n)*facttwo
148            rho43 = rhoval**F43  ! rho^4/3
149            rrho  = 1d0/rhoval   ! reciprocal of rho
150            rho13 = rho43*rrho
151            rho23 = rhoval**F23
152            rho83 = rhoval**F83
153
154            g2 = delrho(n,1)*delrho(n,1) +
155     &           delrho(n,2)*delrho(n,2) +
156     &           delrho(n,3)*delrho(n,3)
157
158            if (dsqrt(g2).gt.tol_rho)then
159               g2 = g2 *facttwo*facttwo
160               g = dsqrt(g2)
161            else
162               g  = tol_rho
163               g2 = tol_rho*tol_rho
164            endif
165
166            tauN = tau(n,1)*facttwo
167            tauW = F18*g2*rrho
168            tauU = 0.3d0*Pconst*rhoval**F53
169
170c
171c     Evaluate the Fx
172c
173            p=g2/(4d0*Pconst*rho83)
174            z=TauW/TauN
175            rz=TauN/TauW
176c
177c            a=(tauN-tauW)/tauU
178c
179            a = F53*p*(rz - 1d0)
180            if(a.lt.0d0)  a=0d0
181
182            FaNum = (1d0 - a)
183            BFaDen = (1d0 + rE1*a*a)**2d0 + rC1*a**4d0
184            FaDen = BFaDen**0.25d0
185            Fa = Fanum / FaDen
186
187            BFxDen = 1d0 + rB*p*p
188            FxDen = BFxDen**F18
189            FxNum = 1.d0 + rK0*Fa
190            Fx = FxNum / FxDen
191
192            Ex = Ex + Fx*Ax*rho43*qwght(n)*fac*afact2
193            if (ldew)  func(n)= func(n) + Fx*Ax*rho43*fac*afact2
194
195c     functional derivatives FFFFFFFFFFFFFFFFFFFFFFFFFFFF
196
197            dpdr = -F83*p*rrho
198            dpdg = p/g2
199
200            dzdr = -z*rrho
201            dzdg =  z/g2
202            dzdt = -z/TauN
203
204            dadp = a/p
205            dadz = F53*p*(-1d0*rz*rz)
206
207            dFaNumda = -1d0
208            dBFaDenda = 2d0*(1d0 + rE1*a*a)*2d0*rE1*a + 4d0*rC1*a**3d0
209            dFaDenda = 0.25d0*(FaDen/BFaDen)*dBFaDenda
210            dFada = (dFaNumda * FaDen - FaNum*dFaDenda)/(FaDen*FaDen)
211
212            dFxNumda = rK0 * dFada
213            dFxDendp = F18*(FxDen/BFxDen)*2d0*rB*p
214
215            dFxdp = dFxNumda*dadp/FxDen - FxNum*dFxDendp/(FxDen*FxDen)
216            dFxdz = dFxNumda*dadz/FxDen
217
218            dFxdr = dFxdz * dzdr + dFxdp * dpdr
219            dFxdg = dFxdz * dzdg + dFxdp * dpdg
220            dFxdt = dFxdz * dzdt
221
222            derivr = F43*Ax*rho13*Fx + Ax*rho43*dFxdr
223            derivg = Ax*rho43*dFxdg
224            derivt = Ax*rho43*dFxdt
225
226            Amat(n) = Amat(n) + derivr*fac
227c
228c     4x factor comes from gamma_aa = gamma_total/4
229c
230            Cmat(n)=  Cmat(n) + 2.0d0*derivg*fac
231            Mmat(n)=  Mmat(n) + 0.5d0*derivt*fac
232         endif
233      enddo
234      return
235      end
236
237      Subroutine xc_xmvs15_d2()
238      call errquit(' xmvs15: d2 not coded ',0,0)
239      return
240      end
241
242
243