1c $Id: xc_eval_fnl.F 27518 2015-09-17 09:34:57Z alogsdail $
2      subroutine xc_xn12(tol_rho, fac,lfac,nlfac, rho, delrho,
3     &                     amat, cmat, nq, ipol, ex,
4     &                     qwght, ldew, func)
5      implicit none
6************************************************************************
7*                                                                      *
8*  N12 evaluates the exchange part of the N12 and N12-SX               *
9*        functionals on the grid.                                      *
10*                                                                      *
11*    OUTPUT:                                                           *
12*      F      - Functional values                                      *
13*      D1F    - First derivatives with respect to RA, RB, GA, GB       *
14*                                                                      *
15*    INPUT:                                                            *
16*       ijzy - 1 N12                                                   *
17*       ijzy - 2 N12-SX                                                *
18*                                                                      *
19*      RA,B     - Spin densities                                       *
20*      D1RA,B   - Spin density gradients                               *
21*      NGrid    - number of grid points                                *
22*                                                                      *
23*  Analytic second derivatives are not available yet!!!                *
24*                                                                      *
25*  RP (09/12)                                                          *
26*                                                                      *
27************************************************************************
28#include "errquit.fh"
29#include "dft2drv.fh"
30c
31      double precision fac, Ex
32      integer nq, ipol
33      logical lfac, nlfac,ldew,   uselc
34      double precision func(*)  ! value of the functional [output]
35c
36c     Charge Density & Its Cube Root
37c
38      double precision rho(nq,ipol*(ipol+1)/2)
39c
40c     Charge Density Gradient
41c
42      double precision delrho(nq,3,ipol)
43c
44c     Quadrature Weights
45c
46      double precision qwght(nq)
47c
48c     Sampling Matrices for the XC Potential & Energy
49c
50      double precision Amat(nq,ipol), Cmat(nq,*)
51      double precision tol_rho
52      external  xc_xn12_os
53      call xc_os2cs(xc_xn12_os,
54     &     tol_rho, fac,lfac,nlfac, rho, delrho,
55     &                     amat, cmat, nq, ipol, ex,
56     &                     qwght, ldew, func)
57      return
58      end
59      subroutine xc_xn12_os(tol_rho, fac,lfac,nlfac, rho, delrho,
60     &                     amat, cmat, nq,  ex,
61     &                     qwght, ldew, func, fact_cs)
62      Implicit Real*8(A-H,O-Z)
63#include "errquit.fh"
64#include "dft2drv.fh"
65c
66c wrapper for n12 open-shell
67      double precision fac, Ex
68      double precision fact_cs
69      integer nq, ipol
70      logical lfac, nlfac,ldew,   uselc
71      double precision func(*)  ! value of the functional [output]
72c
73c     Charge Density & Its Cube Root
74c
75      double precision rho(nq)
76c
77c     Charge Density Gradient
78c
79      double precision delrho(nq,3)
80c
81c     Quadrature Weights
82c
83      double precision qwght(nq)
84c
85c     Sampling Matrices for the XC Potential & Energy
86c
87      double precision Amat(nq), Cmat(nq)
88      double precision tol_rho
89
90      call xc_xn12_0(tol_rho, fac,lfac,nlfac, rho, delrho,
91     &                     amat, cmat, nq,  ex,
92     &                     qwght, ldew, func, 1,fact_cs)
93      return
94      end
95      subroutine xc_xn12sx(tol_rho, fac,lfac,nlfac, rho, delrho,
96     &                     amat, cmat, nq, ipol, ex,
97     &                     qwght, ldew, func)
98      implicit none
99#include "errquit.fh"
100#include "dft2drv.fh"
101c
102      double precision fac, Ex
103      integer nq, ipol
104      logical lfac, nlfac,ldew,   uselc
105      double precision func(*)  ! value of the functional [output]
106      double precision rho(nq,ipol*(ipol+1)/2)
107      double precision delrho(nq,3,ipol)
108      double precision qwght(nq)
109      double precision Amat(nq,ipol), Cmat(nq,*)
110      double precision tol_rho
111      external  xc_xn12sx_os
112      call xc_os2cs(xc_xn12sx_os,
113     &     tol_rho, fac,lfac,nlfac, rho, delrho,
114     &                     amat, cmat, nq, ipol, ex,
115     &                     qwght, ldew, func)
116      return
117      end
118      subroutine xc_xn12sx_os(tol_rho, fac,lfac,nlfac, rho, delrho,
119     &                     amat, cmat, nq,  ex,
120     &                     qwght, ldew, func, fact_cs)
121      Implicit Real*8(A-H,O-Z)
122#include "errquit.fh"
123#include "dft2drv.fh"
124c
125c wrapper for n12 open-shell
126      double precision fac, Ex
127      double precision fact_cs
128      integer nq, ipol
129      logical lfac, nlfac,ldew,   uselc
130      double precision func(*)  ! value of the functional [output]
131c
132c     Charge Density & Its Cube Root
133c
134      double precision rho(nq)
135c
136c     Charge Density Gradient
137c
138      double precision delrho(nq,3)
139c
140c     Quadrature Weights
141c
142      double precision qwght(nq)
143c
144c     Sampling Matrices for the XC Potential & Energy
145c
146      double precision Amat(nq), Cmat(nq)
147      double precision tol_rho
148
149      call xc_xn12_0(tol_rho, fac,lfac,nlfac, rho, delrho,
150     &                     amat, cmat, nq,  ex,
151     &                     qwght, ldew, func, 2,fact_cs)
152      return
153      end
154      subroutine xc_xn12_0(tol_rho, fac,lfac,nlfac, rho, delrho,
155     &                     amat, cmat, nq,  ex,
156     &                     qwght, ldew, func, ijzy,fact_cs)
157      Implicit Real*8(A-H,O-Z)
158#include "errquit.fh"
159#include "dft2drv.fh"
160c
161      integer ijzy
162      double precision fac, Ex
163      double precision fact_cs
164      integer nq, ipol
165      logical lfac, nlfac,ldew,   uselc
166      double precision func(*)  ! value of the functional [output]
167c
168c     Charge Density & Its Cube Root
169c
170      double precision rho(nq)
171c
172c     Charge Density Gradient
173c
174      double precision delrho(nq,3)
175c
176c     Quadrature Weights
177c
178      double precision qwght(nq)
179c
180c     Sampling Matrices for the XC Potential & Energy
181c
182      double precision Amat(nq), Cmat(nq)
183      double precision tol_rho, pi
184c      Real*8 RA(*),RB(*),D1RA(NGrid,3),D1RB(NGrid,3),F(*),D1F(NGrid,*)
185      Save One,Two,Three,Four,Five,Six,Seven,Eight,Nine
186      Data One/1.0d0/,Two/2.0d0/,Three/3.0d0/,Four/4.0d0/,Five/5.0d0/,
187     $  Six/6.0d0/,Seven/7.0d0/,Eight/8.0d0/,Nine/9.0d0/
188C
189C
190      G    = 0.004d+0
191      ome  = 2.5d+0
192      If (ijzy.eq.1) then
193c
194c      N12
195c
196       CC00 =  1.00000D+00
197       CC01 =  5.07880D-01
198       CC02 =  1.68233D-01
199       CC03 =  1.28887D-01
200       CC10 =  8.60211D-02
201       CC11 = -1.71008D+01
202       CC12 =  6.50814D+01
203       CC13 = -7.01726D+01
204       CC20 = -3.90755D-01
205       CC21 =  5.13392D+01
206       CC22 = -1.66220D+02
207       CC23 =  1.42738D+02
208       CC30 =  4.03611D-01
209       CC31 = -3.44631D+01
210       CC32 =  7.61661D+01
211       CC33 = -2.41834D+00
212      Else If (ijzy.eq.2) then
213c
214c      N12-SX
215c
216       CC00 =  6.81116D-01
217       CC01 =  1.88858D+00
218       CC02 =  1.78590D+00
219       CC03 =  8.79456D-01
220       CC10 = -8.12270D-02
221       CC11 = -1.08723D+00
222       CC12 = -4.18682D+00
223       CC13 = -3.00000D+01
224       CC20 =  5.36236D-01
225       CC21 = -5.45678D+00
226       CC22 =  3.00000D+01
227       CC23 =  5.51105D+01
228       CC30 = -7.09913D-01
229       CC31 =  1.30001D+01
230       CC32 = -7.24877D+01
231       CC33 =  2.98363D+01
232      End If
233      F12   = Two * Six
234      F24   = Four * Six
235      F28   = Four * Seven
236      F2o3  = Two / Three
237      F3o2  = Three / Two
238      F1o3  = One / Three
239      F4o3  = Four / Three
240      F7o3  = Seven / Three
241      F8o3  = Eight / Three
242      F10o3 = F2o3 * Five
243      F28o9 = F28 / Nine
244      Pi    = ACos(-One)
245C
246      Ax = -F3o2*(F4o3*PI)**(-F1o3)
247C
248C     Alpha contributions ...
249C
250      Do 10 n = 1, nq
251        If(rho(n).gt.tol_rho) then
252          pX    = rho(n)
253          GamX2 =(delrho(n,1)*delrho(n,1) +
254     &              delrho(n,2)*delrho(n,2) +
255     &     delrho(n,3)*delrho(n,3))
256          S2    = GamX2*pX**(-F8o3)
257          U     = G*S2/(One+G*S2)
258          E     = Ax*pX**F4o3
259          FV    = ome*pX**F1o3/(One+ome*pX**F1o3)
260c
261       FN12 =CC00  +CC01      *U+CC02      *U**2+CC03      *U**3+
262     $    CC10*FV   +CC11*FV   *U+CC12*FV   *U**2+CC13*FV   *U**3+
263     $    CC20*FV**2+CC21*FV**2*U+CC22*FV**2*U**2+CC23*FV**2*U**3+
264     $    CC30*FV**3+CC31*FV**3*U+CC32*FV**3*U**2+CC33*FV**3*U**3
265c
266       ex = ex + E*FN12*qwght(n)*fac*fact_cs
267       if(ldew) func(n)=func(n) + E*FN12*fac*fact_cs
268c
269c First Derivative
270c
271          ER = F4o3*E/pX
272          S    = Sqrt(S2)
273          GamX = Sqrt(GamX2)
274          SR   = -F4o3*S/pX
275          SG   = S/GamX
276          US   = Two*G*S/((One+G*S*S)**2)
277c
278          dFVdR =  (ome/(Three*pX**F2o3))
279     $                *(One+ome*pX**F1o3)**(-Two)
280c
281        dFN12dU=CC01+2.0d0*CC02      *U+3.0d0*CC03      *U**2+
282     $     CC11*FV   +2.0d0*CC12*FV   *U+3.0d0*CC13*FV   *U**2+
283     $     CC21*FV**2+2.0d0*CC22*FV**2*U+3.0d0*CC23*FV**2*U**2+
284     $     CC31*FV**3+2.0d0*CC32*FV**3*U+3.0d0*CC33*FV**3*U**2
285        dFN12dV=
286     $     CC10+CC11   *U+CC12   *U**2+CC13   *U**3+
287     $     2.0d0*CC20*FV   +2.0d0*CC21*FV*U+
288     $     2.0d0*CC22*FV*U**2+2.0d0*CC23*FV*U**3+
289     $     3.0d0*CC30*FV**2+3.0d0*CC31*FV**2*U+
290     $     3.0d0*CC32*FV**2*U**2+3.0d0*CC33*FV**2*U**3
291c
292        dFN12dR  = dFN12dU*US*SR+dFN12dV*dFVdR
293        dFN12dG  = dFN12dU*US*SG
294c
295            amat(n) = amat(n) + (ER*FN12
296     $                      + E*dFN12dR)*fac
297            cmat(n) = cmat(n) + E*dFN12dG/(Two*GamX)*fac
298          endIf
299 10    Continue
300      Return
301      End
302      Subroutine xc_xn12_d2()
303      call errquit(' not coded ',0,0)
304      return
305      end
306