1#include "dft2drv.fh"
2c    Bc95 correlation functional
3c           META GGA
4C         utilizes ingredients:
5c                              rho   -  density
6c                              delrho - gradient of density
7c                              tau (tauN)- K.S kinetic energy density
8
9      Subroutine xc_bc95(tol_rho, cfac, lcfac, nlcfac, rho, delrho,
10     &                     nq, ipol, Ec, qwght, ldew, func,
11     &                     tau, Amat, Cmat, Mmat,ijmswitch)
12
13
14c
15c$Id$
16c
17c  Reference
18c    Becke, A. D. J. Chem. Phys. 1996, 104, 1040.
19c
20      implicit none
21c
22c
23c
24c     Input and other parameters
25c
26      integer ipol, nq
27
28      double precision cfac
29      logical lcfac, nlcfac
30
31      logical lfac, nlfac
32      double precision fac
33      double precision tol_rho
34c
35      logical ldew
36      double precision func(*)
37c
38c     Threshold parameters
39c
40      double precision DTol,F1, F2, F3, F4,COpp
41      Data F1/1.0d0/,F2/2.0d0/,
42     & F3/3.0d0/,F4/4.0d0/
43c
44c     Correlation energy
45c
46      double precision Ec
47c
48c     Charge Density
49c
50      double precision rho(nq,ipol*(ipol+1)/2)
51c
52c     Charge Density Gradient
53c
54      double precision delrho(nq,3,ipol), gammaval, gam12
55
56c
57c     Kinetic Energy Density
58c
59      double precision tau(nq,ipol), tauN
60
61c
62c     Quadrature Weights
63c
64      double precision qwght(nq)
65c
66c     Sampling Matrices for the XC Potential
67c
68      double precision Amat(nq,ipol), Cmat(nq,*)
69      double precision Mmat(nq,*)
70
71      integer n, ijmswitch
72
73c    call to the bc95css subroutine
74      double precision PA,GAA,TA,FA,FPA,FGA,FTA,EUA,EUEGA,ChiA,EUPA
75     &,ChiAP,ChiAG
76      double precision PB,GBB,TB,FB,FPB,FGB,FTB,EUB,EUEGB,ChiB,EUPB
77     &,ChiBP,ChiBG
78c
79      double precision  sop
80      double precision Pi, F6, F43, Pi34, F13,
81     &RS,RSP,Zeta,dZdA,dZdB,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ
82      double precision P, EUEG,Denom, DenPA, DenPB, DenGA, DenGB
83      double precision EUEGPA,EUEGPB
84
85
86c
87c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
88c
89      DTol=tol_rho
90      sop=1.0d0
91      if (ijmswitch.eq.1) then
92C     Parameters for PW6B95 Correlation
93        COpp=0.00262d0
94      elseif (ijmswitch.eq.2) then
95C     Parameters for PWB6K Correlation
96        COpp=0.00353d0
97      else
98C     Parameters for B95 (data statement is unsafe as values may survive
99C     between calls)
100        COpp=0.0031d0
101      endif
102      Pi = F4*ATan(F1)
103      F6=6.0d0
104      F43 = F4 / F3
105      Pi34 = F3 / (F4*Pi)
106      F13 = F1 / F3
107
108      do 20 n = 1, nq
109       if (rho(n,1).lt.DTol) goto 20
110       if (ipol.eq.1) then
111c
112c    get the density, gradient, and tau for the alpha spin from the total
113c
114         PA = rho(n,1)/F2
115         GAA = (    delrho(n,1,1)*delrho(n,1,1) +
116     &                 delrho(n,2,1)*delrho(n,2,1) +
117     &                 delrho(n,3,1)*delrho(n,3,1))/4.0d0
118c  In the bc95css subroutine, we use 2*TA as the tau, so we do not divide
119c  the tau by 2 here
120
121         TA = tau(n,1)
122!         TA=0.0005d0
123
124         Call bc95ss(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
125     &                ChiA,EUPA,ChiAP,ChiAG,ijmswitch)
126         PB = PA
127         GBB = GAA
128         TB = TA
129         FB = FA
130         FPB = FPA
131         FGB = FGA
132         FTB = FTA
133         EUB = EUA
134         ChiB = ChiA
135         EUPB = EUPA
136         ChiBP = ChiAP
137         ChiBG = ChiAG
138
139         Ec = Ec + 2.0d0*FA*qwght(n)            !factor of 2 account for both spin
140         if(ldew) func(n)=func(n)+ 2.0d0*FA
141         Amat(n,1)=Amat(n,1)+ FPA
142         Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + FGA
143         Mmat(n,1)=  Mmat(n,1) + FTA
144#if 0
145      write (0,'(A,3F20.6)') " Amat Cmat Mmat",FPA,FGA,FTA
146#endif
147
148
149c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
150      else  ! ipol=2
151c
152c        ======> SPIN-UNRESTRICTED <======
153c
154c
155c       alpha
156c
157
158         PA = rho(n,2)
159         if (PA.le.DTol) go to 25
160         GAA =   delrho(n,1,1)*delrho(n,1,1) +
161     &           delrho(n,2,1)*delrho(n,2,1) +
162     &          delrho(n,3,1)*delrho(n,3,1)
163c
164c  In the bc95css subroutine, we use 2*TA as the tau
165c
166         TA = tau(n,1)*2.0d0
167
168         Call bc95ss(PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
169     &                ChiA,EUPA,ChiAP,ChiAG,ijmswitch)
170         Ec = Ec + FA*qwght(n)
171         if(ldew) func(n)=func(n)+ FA
172         Amat(n,1)=Amat(n,1)+ FPA
173         Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + FGA
174c      2*0.5=1.0 for Mmat
175         Mmat(n,1)=  Mmat(n,1) + FTA
176#if 0
177      write (0,'(A,3F20.6)') "AAmat Cmat Mmat",FPA,FGA,FTA
178#endif
179c
180c  In the bc95css subroutine, we use 2*TA as the tau,
181c
182c
183c       Beta
184c
185 25       continue
186         PB = rho(n,3)
187         if(PB.le.DTol) go to 30
188         GBB =   delrho(n,1,2)*delrho(n,1,2) +
189     &           delrho(n,2,2)*delrho(n,2,2) +
190     &          delrho(n,3,2)*delrho(n,3,2)
191
192         TB = tau(n,2)*2.0d0
193
194         Call bc95ss(PB,GBB,TB,FB,FPB,FGB,FTB,EUB,
195     &                ChiB,EUPB,ChiBP,ChiBG,ijmswitch)
196         Ec = Ec + FB*qwght(n)
197         if(ldew) func(n)=func(n)+ FB
198         Amat(n,2)= Amat(n,2)+ FPB
199         Cmat(n,D1_GBB)=  Cmat(n,D1_GBB) + FGB
200         Mmat(n,2)=  Mmat(n,2) + FTB
201#if 0
202      write (0,'(A,3F20.6)') "BAmat Cmat Mmat",FPB,FGB,FTB
203#endif
204      endif
205 30   continue
206      P = rho(n,1)
207      If(PA.gt.DTol.and.PB.gt.DTol) then
208          RS = (Pi34/P) ** F13
209          RSP = -RS/(F3*P)
210          Zeta = (PA-PB)/P
211          dZdA = (F1-Zeta)/P
212          dZdB = (-F1-Zeta)/P
213          Call lsdac(dtol,
214     D         RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,
215     $      d2LdZZ)
216          EUEG = P*PotLC - EUA - EUB
217          Denom = F1 + COpp*(ChiA+ChiB)
218          Ec = Ec + sop*EUEG*qwght(n)/Denom
219          if(ldew) func(n)=func(n)+ sop*EUEG/Denom
220          DenPA = COpp*ChiAP
221          DenPB = COpp*ChiBP
222          DenGA = COpp*ChiAG
223          DenGB = COpp*ChiBG
224          EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA
225          EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB
226          if (ipol.eq.1) then
227            Amat(n,1) = Amat(n,1) +
228     &               sop*(EUEGPA/Denom - EUEG*DenPA/Denom**2)
229            Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) - sop*(EUEG*DenGA/Denom**2)
230          else
231            Amat(n,1) = Amat(n,1) +
232     &                 sop*(EUEGPA/Denom - EUEG*DenPA/Denom**2)
233            Amat(n,2) = Amat(n,2) +
234     &                 sop*(EUEGPB/Denom - EUEG*DenPB/Denom**2)
235            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) - sop*EUEG*DenGA/Denom**2
236            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) - sop*EUEG*DenGB/Denom**2
237          endif
238      endIf
239c      write (*,*) "Amat(n,1),Cmat(n,1),Mmat(n,1)",Amat(n,1),Cmat(n,1)
240c     & ,Mmat(n,1)
241c      stop
24220    continue
243      end
244
245      Subroutine xc_bc95_d2()
246      call errquit(' bc95: d2 not coded ',0,0)
247      return
248      end
249
250
251
252
253      Subroutine bc95ss(PX,GX,TX,F,FP,FG,FT,EUEG,Chi,EUEGP,
254     &                   ChiP,ChiG,ijmswitch)
255      Implicit none
256C
257C     Compute the same-spin part of the bc95 correlation functional for one grid
258C     point and one spin-case.
259C
260C
261      integer ijmswitch
262      double precision PX, GX, TX, F, FP, FG, FT, DTol
263      double precision EUEG, Chi, EUEGP, ChiP, ChiG,Css
264      double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11
265      double precision Pi, Pi34, F13, F23, F43, F53, F83, F113
266      double precision RS, FDUEG, D,  RSP,DUEG, Denom, PotLC
267      double precision E, DenomG, DenomP, DUEGP, DP, DG, DT
268      double precision d2LdSS,d2LdSZ,d2LdZZ,dLdS,dLdZ
269
270
271
272      Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/,
273     $  F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/
274C
275      if (ijmswitch.eq.1) then
276C     Parameters for PW6B95 Correlation
277       Css=0.03668d0
278      elseif (ijmswitch.eq.2) then
279C     Parameters for PWB6K Correlation
280       Css=0.04120d0
281      else
282C     Parameters for B95 (data statement is unsafe as values may survive
283C     between calls)
284       Css=0.038d0
285      endif
286      DTol =1.0d-6
287      If(PX.le.DTol) then
288        EUEG = Zero
289        Chi = Zero
290        EUEGP = Zero
291        ChiP = Zero
292        ChiG = Zero
293        PX = Zero
294        GX = Zero
295        TX = Zero
296        F  = Zero
297        FP = Zero
298        FG = Zero
299        FT = Zero
300      else
301        Pi = F4*ATan(F1)
302        Pi34 = F3 / (F4*Pi)
303        F13 = F1 / F3
304        F23 = F2 / F3
305        F43 = F2 * F23
306        F53 = F5 / F3
307        F83 = F8 / F3
308        F113 = F11 / F3
309        FDUEG = (F3/F5)*(F6*Pi*Pi)**F23
310        RS = (Pi34/PX) ** F13
311        Call lsdac(dtol,
312     D       RS,F1,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
313        EUEG = PX*PotLC
314        D = TX - Pt25*GX/PX
315        DUEG = FDUEG*PX**F53
316        Chi = GX/PX**F83
317        Denom = F1 + Css*Chi
318        E = D*EUEG/(DUEG*Denom*Denom)
319c        write (*,*) "ijmswitch, Css, E= ",ijmswitch, Css, E
320c        stop
321        F = E
322c
323        RSP = -RS/(F3*Px)
324        ChiG = F1/PX**F83
325        ChiP = -F83*Chi/PX
326        DenomG = Css*ChiG
327        DenomP = Css*ChiP
328        DUEGP = F53*DUEG/PX
329        DP = Pt25*GX/PX**2
330        DG = -Pt25/PX
331        DT = F1
332        EUEGP = PotLC + PX*dLdS*RSP
333        FP = DP*EUEG/(DUEG*Denom*Denom) +
334     $      D*EUEGP/(DUEG*Denom*Denom)
335     $      - D*EUEG*DUEGP/(DUEG*Denom)**2 -
336     $      F2*D*EUEG*DenomP/(DUEG*Denom*Denom*Denom)
337        FG =DG*EUEG/(DUEG*Denom*Denom) -
338     $      F2*D*EUEG*DenomG/(DUEG*Denom*Denom*Denom)
339        FT =DT*EUEG/(DUEG*Denom*Denom)
340       Endif
341       Return
342       End
343
344
345