1#if defined(FUJITSU_VPP)
2!ocl scalar
3#endif
4#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
5      Subroutine xc_att_xc(rho,ipol,Ex,Amat,Cmat)
6#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
7      Subroutine xc_att_xc_d2(rho,ipol,Ex,Amat,Cmat,Amat2,Cmat2,Cmat3)
8#else
9      Subroutine xc_att_xc_d3(rho,ipol,Ex,Amat,Cmat,Amat2,Cmat2,Cmat3,
10     1                        Amat3,Cmat4,Cmat5,Cmat6)
11#endif
12c
13      implicit none
14c
15#include "stdio.fh"
16#include "case.fh"
17c
18      double precision rho, Ex, Amat, Cmat
19      integer ipol
20
21#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
22c     Second Derivatives of the Exchange Energy Functional
23      double precision Amat2, Cmat2, Cmat3
24#endif
25#ifdef THIRD_DERIV
26c     Third Derivatives of the Exchange Energy Functional
27      double precision Amat3, Cmat4, Cmat5, Cmat6
28#endif
29c
30c
31c References:
32c
33c
34c***************************************************************************
35c
36      double precision a, b, c, btmp,bfactor
37c
38      double precision a_first,a2_first,btmp_first, btmp1
39c
40      double precision sqrt_pi,t1,t2,t3,t4,t5,t6,t7
41      double precision alpha,beta, DERF
42      double precision f10, f01, b_first
43      double precision a2, a3, a4, a5, a6, a7, a8, a9, a10, a11
44      double precision ta, ta2, ta3, ta4, ta5, ta6, ta7, ta8, ta9,
45     1                 ta10
46      double precision f43, f23
47      double precision expf, erff
48
49      Parameter (sqrt_pi = 1.77245385090552d0)
50      Parameter (t7 = 2.666666666666666667d0)
51      Parameter (f43 = 4.0d0/3.0d0)
52      Parameter (f23 = 2.0d0/3.0d0)
53c
54#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
55c     Second Derivatives of the Exchange Energy Functional
56      double precision a_second, a2_second, f20
57      double precision b_second, btmp_second, t8
58      double precision a3_second
59      double precision f11, f02
60#endif
61#ifdef THIRD_DERIV
62      double precision a_third, a2_third, a3_third, a4_third
63      double precision f30, f21, f12, f03, f02a
64      double precision b_third, btmp_third
65      double precision t9
66#endif
67
68
69c calculate the a_sigma parameter
70
71c         write(luout,*) 'alpha',alpha
72c         write(luout,*) 'beta',beta
73c         write(luout,*) 'mu',mu
74c
75          if (ipol.eq.1) then
76             Ex = Ex/2d0
77             rho = rho/2d0
78          endif
79          if(ex.le.0) then
80             a = cam_omega*sqrt(-2d0*Ex)/(6d0*sqrt_pi*rho)
81          else
82             write(6,*) ' negative Ex ',Ex
83             call errquit(' xc_att_xc: negative Ex ',0,0)
84          endif
85          alpha = cam_alpha
86          beta = cam_beta
87c
88          f10 = Amat/(2d0*Ex) -1d0/rho
89          a_first = f10*a
90          f01 = Cmat/(2d0*Ex)
91          a2_first = f01*a
92#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
93          f20 = Amat2/(2d0*Ex) - Amat*Amat/(2d0*Ex*Ex)
94     &        + 1d0/(rho*rho)
95          f11 = Cmat2 - Amat*Cmat/Ex
96          f11 = f11/(2.0d0*Ex)
97c
98          f02 = Cmat3 - Cmat*Cmat/(2.0d0*Ex)
99          f02 = f02/(2.0d0*Ex)
100c
101          a_second = a*(f10*f10 + f20)
102
103c          a2_second = a*(f10*f01 + Cmat2/(2d0*Ex)
104c     &              - Amat*Cmat/(2d0*Ex*Ex))
105          a2_second = a*(f10*f01 + f11)
106
107c          a3_second = a*(Cmat3/(2d0*Ex) - Cmat*Cmat/(4d0*Ex*Ex))
108          a3_second = a*f02
109#endif
110#ifdef THIRD_DERIV
111c Amat3 = drdrdr
112c Cmat4 = drdrdg
113c Cmat5 = drdgdg
114c Cmat6 = dgdgdg
115c          f02a = Cmat3/(2.0d0*Ex) - Cmat*Cmat/(2.0d0*Ex*Ex)
116          f02a = Cmat3 - Cmat*Cmat/Ex
117          f02a = f02a/(2.0d0*Ex)
118c
119          f30 = Amat3/(2.0d0*Ex)
120     1        - 3.0d0*Amat2*Amat/(2.0d0*Ex*Ex)
121     2        + Amat*Amat*Amat/(Ex*Ex*Ex)
122     3        - 2.0d0/(rho*rho*rho)
123c
124          f21 = Cmat4/(2.0d0*Ex)
125     1        - Cmat2*Amat/(Ex*Ex)
126     2        - Amat2*Cmat/(2.0d0*Ex*Ex)
127     3        + Amat*Amat*Cmat/(Ex*Ex*Ex)
128c
129          f12 = Cmat5/(2.0d0*Ex)
130     1        - Cmat2*Cmat/(Ex*Ex)
131     2        - Amat*Cmat3/(2.0d0*Ex*Ex)
132     3        + Amat*Cmat*Cmat/(Ex*Ex*Ex)
133c
134          f03 = Cmat6/(2.0d0*Ex)
135     1        - Cmat3*Cmat/(Ex*Ex)
136     2        + Cmat*Cmat*Cmat/(2.0d0*Ex*Ex*Ex)
137c
138          a_third = a*( f10*f10*f10 + 3.0d0*f10*f20 + f30 )
139c
140          a2_third = a*( f10*f10*f01 + f20*f01 + 2.0d0*f10*f11 + f21 )
141c
142          a3_third = a*( f10*f01*f01 + 2.0d0*f11*f01 + f10*f02a + f12 )
143c
144          a4_third = a*( f01*f02 + f03 )
145#endif
146          a2 = a*a
147#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
148          a4 = a2*a2
149          a6 = a4*a2
150          a8 = a6*a2
151          a10 = a8*a2
152#endif
153c
154#ifdef THIRD_DERIV
155          a3 = a2*a
156          a5 = a4*a
157          a7 = a6*a
158          a9 = a8*a
159          a11 = a10*a
160#endif
161          ta = 2d0*a
162          ta2 = ta*ta
163          ta3 = ta2*ta
164          ta4 = ta3*ta
165          ta5 = ta4*ta
166          ta6 = ta5*ta
167          ta7 = ta6*ta
168          ta8 = ta7*ta
169          ta9 = ta8*ta
170          ta10 = ta9*ta
171c
172          if (a .lt. 1.0d-10) then
173            expf = 0.0d0
174            erff = 1.0d0
175          else
176            expf = exp(-1d0/(4d0*a2))
177            erff = DERF(1d0/(2d0*a))
178          endif
179c
180          if (a .lt. 0.14d0) then
181c             write(luout,*) 'a is small'
182c OLD CODE
183c              a = 2d0*a
184c              btmp = 1d0-(4d0/3d0)*sqrt_pi*a + 2d0*a*a
185c     &             - (2d0/3d0)*a*a*a*a
186c              btmp = 1d0-btmp
187c
188c              btmp_first = (4d0/3d0)*(-sqrt_pi + 3d0*a +
189c     &                   (2d0*exp(-1/(a*a)) - 2d0)*a*a*a)
190c              btmp_first = 2d0*btmp_first
191c              a = a /2d0
192c OLD CODE
193              btmp = 1.0d0 - f43*sqrt_pi*ta
194     1             + 2.0d0*ta2 - f23*ta4
195              btmp = 1.0d0 - btmp
196
197              btmp_first = f43*( -sqrt_pi + 3.0d0*ta +
198     &                           (2.0d0*expf - 2.0d0)*ta3 )
199              btmp_first = 2.0d0*btmp_first
200          else if (a .lt. 4.25d0) then
201c            write(luout,*) 'a is medium'
202c             stop
203             b = expf - 1d0
204c             c = 2d0*a*a*b + 0.5d0
205             c = 2d0*a2*b + 0.5d0
206c             btmp = (8d0/3d0)*a*(sqrt_pi*DERF(1/(2d0*a)) + 2d0*a*(b-c))
207             btmp = (8d0/3d0)*a*(sqrt_pi*erff + 2d0*a*(b-c))
208c OLD CODE
209c             t1 = 1/a
210c             t2 = a*a
211c             t3 = 1/t2
212cc             t4 = exp(-0.25d0*t3)
213c             t4 = expf
214c             t5 = t4 -1d0
215c             t6 = t4 -2d0*t2*t5 - 1.5d0
216c             btmp_first = -t7*a *
217c     &       (2*a*(t4/(2*a**3) - 4d0*a*t5 - t1*t4) + 2d0*t6 -t3*t4) -
218c     &         t7*(2*a*t6 + sqrt_pi*DERF(0.5d0*t1))
219c             btmp_first = -t7*a *
220c     &       (2*a*(t4/(2*a**3) - 4d0*a*t5 - t1*t4) + 2d0*t6 -t3*t4) -
221c     &         t7*(2*a*t6 + sqrt_pi*erff)
222c OLD CODE
223c Daniel (4-12-13) This is a simplified form of what was written above.
224c It is more efficient and likely more stable.
225             btmp_first = -2.0d0*t7*a *
226     &       ( -8.0d0*a2*b + expf - 3.0d0 ) - t7*sqrt_pi*erff
227          else
228c            write(luout,*) 'a is large'
229c            stop
230c OLD CODE
231c             a = 2d0*a
232c             btmp = 1d0 - 1d0/(9d0*a*a) + 1d0/(60d0*a**4d0) -
233c     &           1d0/(420d0*a**6d0) + 1d0/(3240d0*a**8d0) -
234c     &           1d0/(27720d0*a**10d0)
235c
236c             btmp_first = -1d0/(4.5d0*a**3) + 1d0/(15d0*a**5d0) -
237c     &                  1d0/(70d0*a**7d0) + 1d0/(405d0*a**9d0)
238c             btmp_first = btmp_first*2d0
239c             a = a /2d0
240c OLD CODE
241c
242             btmp = 1.0d0 - 1.0d0/(9.0d0*ta2) + 1.0d0/(60.0d0*ta4)
243     1            - 1.0d0/(420.0d0*ta6) + 1.0d0/(3240.0d0*ta8)
244     2            - 1.0d0/(27720.0d0*ta10)
245
246             btmp_first = -1.0d0/(4.5d0*ta3) + 1.0d0/(15.0d0*ta5)
247     1                  - 1.0d0/(70.0d0*ta7) + 1.0d0/(405.0d0*ta9)
248             btmp_first = btmp_first*2.0d0
249          end if
250#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
251          if (abs(a) .lt. 1d-40) then
252c             btmp_second = 16d0
253             btmp_second = 16.0d0
254c Daniel (4-5-13): Should this be 4.25d0? (This probably doesn't
255c matter since both exponentials and error functions vary slowly
256c as the argument gets large).
257          else if (a .ge. 5d0)  then
258c OLD CODE
259c             btmp_second = 1d0/(6d0*a**4d0) - 1d0/(48d0*a**6d0) +
260c     &                  1d0/(640d0*a**8d0) - 1d0/(11520d0*a**10d0)
261c OLD CODE
262c
263             btmp_second = 1.0d0/(6.0d0*a4) - 1.0d0/(48.0d0*a6)
264     1                   + 1.0d0/(640.0d0*a8) - 1.0d0/(11520.0d0*a10)
265
266          else
267c OLD CODE
268c             t1 = a*a
269c             t2 = 1d0/t1
270cc             t3 = exp(-0.25d0*t2)
271c             t3 = expf
272c             t4 = 1d0/(a*a*a)
273c             t5 = t3 - 1d0
274c             t6 = -t2*t3
275c             t8 = -t3/a + 0.5d0*t4*t3 - 4d0*a*t5
276cc
277c             btmp_second = -(8d0*a*(2d0*a*(t3/(4*a**6d0) -
278c     &       2d0*t3/(a**4d0) +t6 - 4d0*t5) -t3/(2*a**5d0) +
279c     &       4d0*t8 + 2d0*t4*t3)/3d0 + 16d0*(2d0*a*t8 +
280c     &       2d0*(t3 - 2d0*t1*t5-1.5d0) + t6)/3d0)
281c OLD CODE
282c Daniel (4-12-13): This is a simplified form of what was written above.
283c It is likely more stable and efficient.
284             btmp_second = 16.0d0 - 128.0d0*a2
285     &                   + (16.0d0 + 128.0d0*a2)*expf
286          end if
287#endif
288#ifdef THIRD_DERIV
289          if (abs(a) .lt. 1.0d-40) then
290             btmp_third = 0.0d0
291          else if (a .ge. 5.0d0) then
292             btmp_third = -2.0d0/(3.0d0*a5)
293     1                  + 1.0d0/(8.0d0*a7)
294     2                  - 1.0d0/(80.0d0*a9)
295     3                  + 1.0d0/(1152.0d0*a11)
296          else
297c OLD CODE
298c            t1 = 1d0/(a*a)
299c            t2 = expf
300c            t3 = t1*t1*t1
301c            t4 = t1*t1
302c            t5 = t3*a
303c            t6 = t2 - 1d0
304c            t8 = -t1*t2 - 2d0*t4*t2 + t3*t2/4d0 - 4d0*t6
305c            t9 = t4*a
306c
307c            btmp_third = -8d0*( a*( 2d0*a*( t2/(8*a**9d0)
308c     1                                    - 5d0*t2/(2d0*a**7d0)
309c     2                                    + 15d0*t5*t2/2d0 )
310c     3                            - t2/(4d0*a**8d0) + 6d0*t8
311c     4                            - 6d0*t4*t2 + 7d0*t3*t2/2d0 )/3d0
312c     5                        + ( 4d0*( -t2/a + t9*t2/2d0 - 4d0*a*t6 )
313c     6                          + 2d0*a*t8 + 2d0*t9*t2
314c     7                          - t5*t2/2d0 ) )
315c OLD CODE
316c
317c Daniel (4-12-13): This is a simplified form of what was written above.
318c It is likely more stable and efficient.
319             btmp_third = 8.0d0*( -32.0d0*a4
320     1                          + ( 1.0d0 + 8.0d0*a2
321     2                            + 32.0d0*a4 )*expf )/a3
322          endif
323#endif
324          bfactor = 1d0 - alpha - beta*btmp
325          b_first = beta*btmp_first
326#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
327          b_second = beta*btmp_second
328#endif
329#ifdef THIRD_DERIV
330          b_third = beta*btmp_third
331c
332          Amat3 = bfactor*Amat3
333     1          + 3.0d0*Amat2*b_first*a_first
334     2          + 3.0d0*Amat*( b_second*a_first*a_first
335     3                       + b_first*a_second )
336     4          + Ex*( b_third*a_first*a_first*a_first
337     5               + 3.0d0*b_second*a_first*a_second
338     6               + b_first*a_third )
339c
340          Cmat4 = bfactor*Cmat4
341     1          + 2.0d0*Cmat2*b_first*a_first
342     2          + Amat2*b_first*a2_first
343     3          + 2.0d0*Amat*( b_second*a_first*a2_first
344     4                       + b_first*a2_second )
345     5          + Cmat*( b_second*a_first*a_first
346     6                 + b_first*a_second )
347     7          + Ex*( b_third*a_first*a_first*a2_first
348     8               + b_second*( a2_first*a_second
349     9                          + 2.0d0*a_first*a2_second )
350     A               + b_first*a2_third )
351c
352          Cmat5 = bfactor*Cmat5
353     1          + 2.0d0*Cmat2*b_first*a2_first
354     2          + Amat*( b_second*a2_first*a2_first
355     3                 + b_first*a3_second )
356     4          + Cmat3*b_first*a_first
357     5          + 2.0d0*Cmat*( b_second*a_first*a2_first
358     6                       + b_first*a2_second )
359     7          + Ex*( b_third*a_first*a2_first*a2_first
360     8               + b_second*( a_first*a3_second
361     9                          + 2.0d0*a2_first*a2_second )
362     A               + b_first*a3_third )
363c
364          Cmat6 = bfactor*Cmat6
365     1          + 3.0d0*Cmat3*b_first*a2_first
366     2          + 3.0d0*Cmat*( b_second*a2_first*a2_first
367     3                       + b_first*a3_second )
368     4          + Ex*( b_third*a2_first*a2_first*a2_first
369     5               + 3.0d0*b_second*a2_first*a3_second
370     6               + b_first*a4_third )
371#endif
372c
373#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
374c          b_second = beta*btmp_second
375          Amat2 = bfactor*Amat2 + 2d0*Amat*b_first*a_first
376     &          + Ex*b_second*a_first*a_first
377     &          + Ex*b_first*a_second
378
379          Cmat2 = bfactor*Cmat2 + Amat*b_first*a2_first
380     &          + Cmat*b_first*a_first
381     &          + Ex*b_second*a_first*a2_first
382     &          + Ex*b_first*a2_second
383
384          Cmat3 = bfactor*Cmat3 + 2d0*Cmat*b_first*a2_first
385     &          + Ex*b_second*a2_first*a2_first
386     &          + Ex*b_first*a3_second
387#endif
388          Amat = bfactor*Amat + Ex*b_first*a_first
389          Cmat = bfactor*Cmat + Ex*b_first*a2_first
390          Ex = Ex*bfactor
391
392          if (ipol.eq.1) then
393             Ex = 2d0*Ex
394             rho = 2d0*rho
395          endif
396c
397      return
398      end
399#ifndef SECOND_DERIV
400#define SECOND_DERIV
401c
402c     Compile source again for the 2nd derivative case
403c
404#include "xc_att_xc.F"
405#endif
406#ifndef THIRD_DERIV
407#define THIRD_DERIV
408c
409c     Compile source again for the 3rd derivative case
410c
411#include "xc_att_xc.F"
412#endif
413c $Id$
414