1      subroutine x_rks_s(r,f,dfdra)
2      implicit none
3c
4c     This subroutine evaluates the spin polarised exchange functional
5c     in the Local Density Approximation [1], and the corresponding
6c     potential. Often this functional is referred to as the Dirac
7c     functional [2] or Slater functional.
8c
9c     [1] F. Bloch, Zeitschrift fuer Physik, Vol. 57 (1929) 545.
10c
11c     [2] P.A.M. Dirac, Proceedings of the Cambridge Philosophical
12c         Society, Vol. 26 (1930) 376.
13c
14c     Parameters:
15c
16c     r     the total electron density
17c     f     On return the functional value
18c     dfdra On return the derivative of f with respect to alpha electron
19c           density
20c
21c
22c     Ax = -3/4*(6/pi)**(1/3)
23c     Bx = -(6/pi)**(1/3)
24c     C  = (1/2)**(1/3)
25      real*8 Ax, Bx, C
26      parameter(Ax    = -0.930525736349100025002010218071667d0)
27      parameter(Bx    = -1.24070098179880003333601362409556d0)
28      parameter(C     =  0.793700525984099737375852819636154d0)
29c
30      real*8 third
31      parameter(third =  0.333333333333333333333333333333333d0)
32c
33      real*8 r
34      real*8 f, dfdra
35c
36      real*8 ra13
37c
38      ra13  = C*r**third
39      f     = Ax*r*ra13
40      dfdra = Bx*ra13
41c
42      end
43c-----------------------------------------------------------------------
44      subroutine x_uks_s(ra,rb,f,dfdra,dfdrb)
45      implicit none
46c
47c     This subroutine evaluates the spin polarised exchange functional
48c     in the Local Density Approximation [1], and the corresponding
49c     potential. Often this functional is referred to as the Dirac
50c     functional [2] or Slater functional.
51c
52c     [1] F. Bloch, Zeitschrift fuer Physik, Vol. 57 (1929) 545.
53c
54c     [2] P.A.M. Dirac, Proceedings of the Cambridge Philosophical
55c         Society, Vol. 26 (1930) 376.
56c
57c     Parameters:
58c
59c     ra    the alpha electron density
60c     rb    the beta  electron density
61c     f     On return the functional value
62c     dfdra On return the derivative of f with respect to ra
63c     dfdrb On return the derivative of f with respect to rb
64c
65c
66c     Ax = -3/4*(6/pi)**(1/3)
67c     Bx = -(6/pi)**(1/3)
68      real*8 Ax, Bx
69      parameter(Ax    = -0.930525736349100025002010218071667d0)
70      parameter(Bx    = -1.24070098179880003333601362409556d0)
71c
72      real*8 third
73      parameter(third =  0.333333333333333333333333333333333d0)
74c
75      real*8 ra, rb
76      real*8 f, dfdra, dfdrb
77c
78      real*8 ra13, rb13
79c
80      ra13  = ra**third
81      rb13  = rb**third
82      f     = Ax*(ra*ra13+rb*rb13)
83      dfdra = Bx*ra13
84      dfdrb = Bx*rb13
85c
86      end
87      subroutine c_rks_vwn5(r,f,dfdra)
88      implicit none
89c
90c     This subroutine evaluates the Vosko, Wilk and Nusair correlation
91c     functional number 5 [1] for the closed shell case, with the
92c     parametrisation as given in table 5.
93c
94c     The original code was obtained from Dr. Phillip Young,
95c     with corrections from Dr. Paul Sherwood.
96c
97c     [1] S.H. Vosko, L. Wilk, and M. Nusair
98c         "Accurate spin-dependent electron liquid correlation energies
99c          for local spin density calculations: a critical analysis",
100c         Can.J.Phys, Vol. 58 (1980) 1200-1211.
101c
102c     Parameters:
103c
104c     r      the total electron density
105c     f      On return the functional value
106c     dfdra  On return the derivative of f with respect to the alpha
107c            electron density
108c
109      real*8 r
110      real*8 f, dfdra
111c
112      real*8 t4,t5,t6,t7
113      real*8 a2,b2,c2,d2
114      real*8 P1,P2,P3,P4
115      real*8 alpha_rho13
116      real*8 srho,srho13
117      real*8 iv,iv2,inv,i1,i2,i3
118      real*8 pp1,pp2
119c
120      real*8 n13, n43, n16, one, two, four
121      parameter(n13 = 0.333333333333333333333333333333d0)
122      parameter(n43 = 1.333333333333333333333333333333d0)
123      parameter(n16 = 0.166666666666666666666666666666d0)
124      parameter(one = 1.0d0)
125      parameter(two = 2.0d0)
126      parameter(four= 4.0d0)
127C
128C VWN interpolation parameters
129C
130C paramagnetic
131      a2 =  0.0621814d0
132      b2 =  3.72744d0
133      c2 = 12.9352d0
134      d2 = -0.10498d0
135C
136C t4 = (1/(4/3)*pi)**(1/3)
137      t4=0.620350490899399531d0
138C
139C t5 = 0.5/(2**(1/3)-1)
140      t5 = 1.92366105093153617d0
141C
142C t6 = 2/(3*(2**(1/3)-1))
143      t6 = 2.56488140124204822d0
144C
145C t7 = 2.25*(2**(1/3)-1)
146      t7 = 0.584822362263464735d0
147C
148C Paramagnetic interpolation constants
149C
150      P1 = 6.15199081975907980d0
151      P2 = a2 *0.5d0
152      P3 = 0.193804554230887489d-02 *0.5d0
153      P4 = 0.775665897562260176d-01 *0.5d0
154C
155C closed shell case
156      srho        = r
157      srho13      = srho**n13
158      alpha_rho13 = (0.5d0**n13)*srho
159      iv2         = T4/srho13
160      iv          = sqrt(iv2)
161C
162C paramagnetic
163      inv = 1.0d0/(iv2+b2*iv+c2)
164      i1  = log(iv2*inv)
165      i2  = log((iv-d2)*(iv-d2)*inv)
166c corrected b1->b2 ps Apr98
167      i3  = atan(P1/(2.0d0*iv+b2))
168      pp1 = P2*i1 + P3*i2 + P4*i3
169      pp2 = a2*(1.0d0/iv-iv*inv*(1.0d0+b2/(iv-d2)))
170C
171      f     = pp1*srho
172      dfdra = pp1 - n16*iv*pp2
173c
174      end
175c-----------------------------------------------------------------------
176      subroutine c_uks_vwn5(ra,rb,f,dfdra,dfdrb)
177      implicit none
178c
179c     This subroutine evaluates the Vosko, Wilk and Nusair correlation
180c     functional number 5 [1], with the parametrisation as given in
181c     table 5.
182c
183c     The original code was obtained from Dr. Phillip Young,
184c     with corrections from Dr. Paul Sherwood.
185c
186c     [1] S.H. Vosko, L. Wilk, and M. Nusair
187c         "Accurate spin-dependent electron liquid correlation energies
188c          for local spin density calculations: a critical analysis",
189c         Can.J.Phys, Vol. 58 (1980) 1200-1211.
190c
191c     Parameters:
192c
193c     ra     the alpha-electron density
194c     rb     the beta-electron density
195c     f      On return the functional value
196c     dfdra  On return the derivative of f with respect to ra
197c     dfdrb  On return the derivative of f with respect to rb
198c
199      real*8 ra, rb
200      real*8 f, dfdra, dfdrb
201c
202      real*8 t1,t2,t4,t5,t6,t7
203      real*8 a1,b1,c1,d1
204      real*8 a2,b2,c2,d2
205      real*8 a3,b3,c3,d3
206      real*8 S1,S2,S3,S4
207      real*8 P1,P2,P3,P4
208      real*8 F1,F2,F3,F4
209      real*8 inter1,inter2
210      real*8 alpha_rho13,beta_rho13
211      real*8 srho,srho13
212      real*8 iv,iv2,inv,i1,i2,i3
213      real*8 vwn1,vwn2
214      real*8 ss1,ss2,pp1,pp2,ff1,ff2
215      real*8 zeta,zeta3,zeta4
216      real*8 tau,dtau,v
217c
218      real*8 n13, n43, n16, one, two, four, tn13, tn43
219c     tn13 = 2**(1/3)
220c     tn43 = 2**(4/3)
221      parameter(n13 = 0.333333333333333333333333333333d0)
222      parameter(n43 = 1.333333333333333333333333333333d0)
223      parameter(n16 = 0.166666666666666666666666666667d0)
224      parameter(tn13= 1.25992104989487316476721060727823d0)
225      parameter(tn43= 2.51984209978974632953442121455646d0)
226      parameter(one = 1.0d0)
227      parameter(two = 2.0d0)
228      parameter(four= 4.0d0)
229C
230C VWN interpolation parameters
231C
232C spin stiffness
233      a1 = -0.337737278807791058d-01
234      b1 =  1.13107d0
235      c1 = 13.0045d0
236      d1 = -0.00475840d0
237C paramagnetic
238      a2 =  0.0621814d0
239      b2 =  3.72744d0
240      c2 = 12.9352d0
241      d2 = -0.10498d0
242C ferromagnetic
243c try cadpac/nwchem value (.5*a2)
244      a3 = 0.0310907d0
245      b3 =  7.06042d0
246      c3 = 18.0578d0
247      d3 = -0.32500d0
248C
249C t4 = (1/(4/3)*pi)**(1/3)
250      t4=0.620350490899399531d0
251C
252C t5 = 0.5/(2**(1/3)-1)
253      t5 = 1.92366105093153617d0
254C
255C t6 = 2/(3*(2**(1/3)-1))
256      t6 = 2.56488140124204822d0
257C
258C t7 = 2.25*(2**(1/3)-1)
259      t7 = 0.584822362263464735d0
260C
261C Spin stiffness interpolation constants
262C
263      S1 = 7.12310891781811772d0
264      S2 = a1 *0.5d0
265      S3 = -0.139834647015288626d-04 *0.5d0
266      S4 = -0.107301836977671539d-01 *0.5d0
267C
268C Paramagnetic interpolation constants
269C
270      P1 = 6.15199081975907980d0
271      P2 = a2 *0.5d0
272      P3 = 0.193804554230887489d-02 *0.5d0
273      P4 = 0.775665897562260176d-01 *0.5d0
274C
275C Ferromagnetic interpolation constants
276C
277      F1 = 4.73092690956011364d0
278      F2 = a3 *0.5d0
279c
280c      F3 = -0.244185082989490298d-02 *0.5d0
281c      F4 = -0.570212323620622186d-01 *0.5d0
282c
283c  try nwchem values
284c
285      F3 = 0.224786709554261133d-02
286      F4 = 0.524913931697809227d-01
287C
288C Interpolation intervals
289C
290      inter1 =  1.0d0-1.0d-10
291      inter2 = -1.0d0+1.0d-10
292C
293C open shell case
294      alpha_rho13 = ra**n13
295      beta_rho13  = rb**n13
296      srho        = ra+rb
297      srho13      = srho**n13
298      iv2         = T4/srho13
299      iv          = sqrt(iv2)
300C
301C spin-stiffness
302      inv = 1.0d0/(iv2+b1*iv+c1)
303      i1  = log(iv2*inv)
304      i2  = log((iv-d1)*(iv-d1)*inv)
305      i3  = atan(S1/(2.0d0*iv+b1))
306      ss1 = S2*i1 + S3*i2 + S4*i3
307      ss2 = a1*(1.0d0/iv-iv*inv*(1.0d0+b1/(iv-d1)))
308C
309C paramagnetic
310      inv = 1.0d0/(iv2+b2*iv+c2)
311      i1  = log(iv2*inv)
312      i2  = log((iv-d2)*(iv-d2)*inv)
313c corrected b1->b2 ps Apr98
314      i3  = atan(P1/(2.0d0*iv+b2))
315      pp1 = P2*i1 + P3*i2 + P4*i3
316      pp2 = a2*(1.0d0/iv-iv*inv*(1.0d0+b2/(iv-d2)))
317C
318C ferromagnetic
319      inv = 1.0d0/(iv2+b3*iv+c3)
320      i1  = log(iv2*inv)
321      i2  = log((iv-d3)*(iv-d3)*inv)
322      i3  = atan(F1/(2.0d0*iv+b3))
323      ff1 = F2*i1 + F3*i2 + F4*i3
324      ff2 = a3*(1.0d0/iv-iv*inv*(1.0d0+b3/(iv-d3)))
325C
326C polarisation function
327c
328      zeta  = (ra-rb)/srho
329      zeta3 = zeta*zeta*zeta
330      zeta4 = zeta3*zeta
331      if(zeta.gt.inter1) then
332         vwn1 = (tn43-two)*t5
333         vwn2 = (tn13)*t6
334      elseif(zeta.lt.inter2) then
335         vwn1 = (tn43-two)*t5
336         vwn2 = -(tn13)*t6
337      else
338         vwn1 = ((one+zeta)**n43+(one-zeta)**n43-two)*t5
339         vwn2 = ((one+zeta)**n13-(one-zeta)**n13)*t6
340      endif
341      ss1  = ss1*t7
342      ss2  = ss2*t7
343      tau  = ff1-pp1-ss1
344      dtau = ff2-pp2-ss2
345c
346      v = pp1+vwn1*(ss1+tau*zeta4)
347      f = v*srho
348c
349      t1 = v - n16*iv*(pp2+vwn1*(ss2+dtau*zeta4))
350      t2 = vwn2*(ss1+tau*zeta4)+vwn1*four*tau*zeta3
351      dfdra = t1+t2*(one-zeta)
352      dfdrb = t1-t2*(one+zeta)
353c
354      end
355