1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
2#if !defined(NWAD_PRINT)
3C> \ingroup nwxc
4C> @{
5C>
6C> \file nwxc_c_op.F
7C> The OP correlation functional
8C>
9C> @}
10#endif
11#endif
12C>
13C> \ingroup nwxc_priv
14C> @{
15C>
16C> \brief Evaluate the OP correlation functional
17C>
18C> The OP correlation functional [1,2] is a functional designed to
19C> have few optimized parameters (only one in this case) and has
20C> a "progressive" form.
21C>
22C> ### References ###
23C>
24C> [1] T. Tsuneda, T. Suzumura, K. Hirao,
25C>     "A new one-parameter progressive Colle-Salvetti-type correlation
26C>     functional", J. Chem. Phys. <b>110</b>, 10664-10678 (1999), DOI:
27C>     <a href="https://doi.org/10.1063/1.479012">
28C>     10.1063/1.479012</a>.
29C>
30C> [2] T. Tsuneda, T. Suzumura, K. Hirao,
31C>     "A reexamination of exchange energy functionals",
32C>     J. Chem. Phys. <b>111</b>, 5656-5667 (1999), DOI:
33C>     <a href="https://doi.org/10.1063/1.479954">
34C>     10.1063/1.479954</a>.
35C>
36#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
37#if defined(NWAD_PRINT)
38      Subroutine nwxc_c_op_p(kop,param,tol_rho,ipol,nq,wght,rho,rgamma,
39     &                     func)
40#else
41      Subroutine nwxc_c_op(kop,param,tol_rho,ipol,nq,wght,rho,rgamma,
42     &                     func)
43#endif
44#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
45      Subroutine nwxc_c_op_d2(kop,param,tol_rho,ipol,nq,wght,rho,rgamma,
46     &                     func)
47#else
48      Subroutine nwxc_c_op_d3(kop,param,tol_rho,ipol,nq,wght,rho,rgamma,
49     &                     func)
50#endif
51c
52C$Id$
53c
54#include "nwad.fh"
55c
56      implicit none
57c
58#include "nwxc_param.fh"
59c
60c     Input and other parameters
61c
62      external kop             !< [Input] Subroutine to evaluate the
63      !< GGA exchange enhancement factor (see Eq.(14) in [1]).
64#if defined(NWAD_PRINT)
65#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
66      type(nwad_dble)::param(*)!< [Input] Parameters of functional
67      type(nwad_dble)::QABOP
68#else
69      double precision param(*)!< [Input] Parameters of functional
70      double precision QABOP
71#endif
72#else
73      double precision param(*)!< [Input] Parameters of functional
74      !< - param(1): \f$ q^{\alpha\beta}_{OP} \f$ see Eq.(27) in [1] and
75      !<   Table III in [2].
76      double precision QABOP
77#endif
78      double precision tol_rho !< [Input] The lower limit on the density
79      integer ipol             !< [Input] The number of spin channels
80      integer nq               !< [Input] The number of points
81      double precision wght    !< [Input] The weight of the functional
82c
83c     Charge Density
84c
85      type(nwad_dble)::rho(nq,*)    !< [Input] The density
86c
87c     Charge Density Gradient
88c
89      type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients
90c
91c     Sampling Matrices for the XC Potential
92c
93      type(nwad_dble)::func(nq)    !< [Output] The value of the functional
94c     double precision Amat(nq,*)   !< [Output] The derivative wrt rho
95c     double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
96c
97      double precision QAB88OP,QABPBOP
98c     Parameter (QAB88OP=2.3670D0,QABPBOP=2.3789D0)
99c
100c References:
101c    Tsuneda, Suzumura, Hirao, JCP 110, 10664 (1999)
102c    Tsuneda, Suzumura, Hirao, JCP 111, 5656 (1999)
103c
104c***************************************************************************
105c
106      integer n
107      type(nwad_dble)::rho13, rho43, gamma, x
108      type(nwad_dble)::kalpha,kbeta, rho13a, rho13b,rhoa,rhob
109c     type(nwad_dble)::banb, hbab, hbabx
110      type(nwad_dble)::banb, hbab
111      double precision dhdab,dhdabx,dkadra,dkbdrb,dkadxa,dkbdxb,
112     A     dbabdra,dbabdrb,dbabdga,dbabdgb,dkadga,dkbdgb,
113     A     dbabdka,dbabdkb
114c
115c     hbabx(x) = (1.5214d0*x + 0.5764d0)/
116c    /           (x**2.0d0*(x**2.0d0+1.1284d0*x+0.3183d0))
117c     dhdabx(x) = -(4.5642d0*x**4+5.7391d0*x**3+
118c    +     2.4355*x**2+0.3669d0*x)/
119c    /           ((x**4+1.1284d0*x**3+0.3183d0*x**2)**2)
120c
121c     if(whichf.eq.'be88') then
122c        QABOP=QAB88OP
123c     endif
124c     if(whichf.eq.'pb96') then
125c        QABOP=QABPBOP
126c     endif
127      QABOP = param(1)
128      if (ipol.eq.1) then
129c
130c        ======> SPIN-RESTRICTED <======
131c
132         do 10 n = 1, nq
133            if (rho(n,R_T).lt.tol_rho) goto 10
134c
135c           Spin alpha:
136c
137            rhoa=rho(n,R_T)*0.5d0
138            rho13a = (rhoa)**(1.d0/3.d0)
139            rho43 = rho13a**4.0d0
140            gamma = rgamma(n,G_TT)
141c           gamma = delrho(n,1,1)*delrho(n,1,1) +
142c    &              delrho(n,2,1)*delrho(n,2,1) +
143c    &              delrho(n,3,1)*delrho(n,3,1)
144            gamma = 0.25d0 * gamma
145            if (sqrt(gamma).gt.tol_rho)then
146               x = sqrt(gamma) / rho43
147               call kop(tol_rho,x,kalpha)
148c              dkadra = -(4d0/3d0)*x*dkadxa/rhoa
149c              dkadga = (dkadxa/rho43)*0.5d0/sqrt(gamma)
150            else
151               x=0d0
152               call kop(tol_rho,x,kalpha)
153c              dkadra = 0d0
154c              dkadga = 0d0
155            endif
156c
157c
158
159            banb = qabop * rho13a * kalpha *0.5d0
160
161            if(banb.ne.0.0d0) then
162c              dbabdra = banb*0.5d0*
163c    /              (1d0/(3d0*rhoa)+dkadra/kalpha)
164
165c              dbabdga = banb/kalpha*dkadga*0.5d0
166
167               hbab = hbabx(banb)
168c              dhdab = dhdabx(banb)
169            else
170c              dbabdra =0d0
171c              dbabdga =0d0
172               hbab = 0d0
173c              dhdab = 0d0
174            endif
175
176c           Ec = Ec - rhoa**2*hbab*qwght(n)*fac
177            func(n) = func(n) - rhoa**2.0d0*hbab*wght
178c           Amat(n,D1_RA) = Amat(n,D1_RA) -
179c    -           (rhoa*hbab + rhoa**2*dhdab*dbabdra)*wght
180
181c
182c           if (x.gt.tol_rho) then
183c               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) -
184c    -              rhoa**2*dhdab*dbabdga*wght
185c            endif
186c
187 10      continue
188c
189      else
190c
191c        ======> SPIN-UNRESTRICTED <======
192c
193         do 20 n = 1, nq
194            if (rho(n,R_A)+rho(n,R_B).lt.tol_rho) goto 20
195            if (rho(n,R_A).ge.tol_rho*0.5d0)  then
196c
197c           Spin alpha:
198c
199               rhoa=rho(n,R_A)
200               rho13a = rhoa**(1.d0/3.d0)
201               rho43 = rho13a*rhoa
202               gamma = rgamma(n,G_AA)
203c              gamma = delrho(n,1,1)*delrho(n,1,1) +
204c    &              delrho(n,2,1)*delrho(n,2,1) +
205c    &              delrho(n,3,1)*delrho(n,3,1)
206               if (sqrt(gamma).gt.tol_rho)then
207                  x = sqrt(gamma) / rho43
208                  call kop(tol_rho,x,kalpha)
209c                 dkadra = -(4d0/3d0)*x*dkadxa/rhoa
210c                 dkadga = dkadxa*0.5d0/(rho43*dsqrt(gamma))
211               else
212                  x = 0d0
213               endif
214            else
215               rhoa=0d0
216               rho13a=0d0
217               x = 0d0
218            endif
219            if(x.eq.0d0) then
220               call kop(tol_rho,x,kalpha)
221c              dkadra = 0d0
222c              dkadga = 0d0
223            endif
224c
225c           Spin beta:
226c
227            if (rho(n,R_B).ge.tol_rho*0.5d0) then
228c
229               rhob=rho(n,R_B)
230               rho13b = rhob**(1.d0/3.d0)
231               rho43 = rho13b*rhob
232               gamma = rgamma(n,G_BB)
233c              gamma = delrho(n,1,2)*delrho(n,1,2) +
234c    &              delrho(n,2,2)*delrho(n,2,2) +
235c    &              delrho(n,3,2)*delrho(n,3,2)
236               if (sqrt(gamma).gt.tol_rho)then
237                  x = sqrt(gamma) / rho43
238                  call kop(tol_rho,x, kbeta)
239c                 dkbdrb = -(4d0/3d0)*x*dkbdxb/rhob
240c                 dkbdgb = dkbdxb*0.5d0/(rho43*dsqrt(gamma))
241               else
242                  x = 0d0
243               endif
244            else
245               if(rho13a.eq.0.0d0) goto 20
246               rhob=0d0
247               rho13b=0d0
248               x=0d0
249            endif
250            if(x.eq.0d0) then
251               call kop(tol_rho,x, kbeta)
252c              dkbdrb = 0d0
253c              dkbdgb=  0d0
254            endif
255
256            banb = qabop*(rho13a*kalpha*rho13b*kbeta)/
257     /           (rho13a*kalpha+rho13b*kbeta)
258
259            if(banb.ne.0.0d0) then
260c              dbabdra = banb*kbeta*rho13b/
261c    /              (rho13a*kalpha+rho13b*kbeta)*
262c    /              (1d0/(3d0*rhoa)+dkadra/kalpha)
263c              dbabdrb = banb*kalpha*rho13a/
264c    /              (rho13a*kalpha+rho13b*kbeta)*
265c    /              (1d0/(3d0*rhob)+dkbdrb/kbeta)
266
267c              dbabdga = banb*rho13b*kbeta/
268c    /              ((rho13a*kalpha+rho13b*kbeta)*kalpha)*
269c    *              dkadga
270c              dbabdgb = banb*rho13a*kalpha/
271c    /              ((rho13a*kalpha+rho13b*kbeta)*kbeta)*
272c    *              dkbdgb
273
274               hbab = hbabx(banb)
275c              dhdab = dhdabx(banb)
276            else
277c              dbabdra =0d0
278c              dbabdrb =0d0
279c              dbabdga =0d0
280c              dbabdgb =0d0
281               hbab = 0d0
282c              dhdab = 0d0
283            endif
284
285c           Ec = Ec - rhoa*rhob*hbab*qwght(n)*fac
286            func(n) = func(n) - rhoa*rhob*hbab*wght
287c           Amat(n,D1_RA) = Amat(n,D1_RA) -
288c    -           (rhob*hbab + rhoa*rhob*dhdab*dbabdra)*wght
289c           Amat(n,D1_RB) = Amat(n,D1_RB) -
290c    -           (rhoa*hbab + rhoa*rhob*dhdab*dbabdrb)*wght
291c
292c
293c           if (x.gt.tol_rho) then
294c              Cmat(n,D1_GAA) = Cmat(n,D1_GAA) -
295c    -              rhoa*rhob*dhdab*dbabdga*wght
296c              Cmat(n,D1_GBB) = Cmat(n,D1_GBB) -
297c    -              rhoa*rhob*dhdab*dbabdgb*wght
298c           endif
299
300c
301c
302 20      continue
303c
304      endif
305c
306      return
307c
308      contains
309c
310c     The combination of statement functions and derived types with
311c     overloaded operators is not properly supported in the Fortran
312c     standard (apparently). Therefore the statement functions from
313c     the original subroutine had to be transformed into contained
314c     functions.
315c
316c     WARNING: This code is EXTREMELY EVIL! Although there appears to be
317c     only one function there are actually three with the same name,
318c     each one returning results of a different data type. The trick is
319c     that depending on the data type the the subroutine that contains
320c     these functions changes its name and hence these different
321c     functions of the same name do not lead to conflicts. The
322c     alternative would have been to add a forest of conditional
323c     compilation constructs (#ifdef's) to change the function names
324c     in the declarations and all places where they are used. That
325c     would have been extremely ugly, so we are between a rock and a
326c     hard place on this one.
327c
328      function hbabx(x) result(s)
329#include "nwad.fh"
330        implicit none
331        type(nwad_dble), intent(in) :: x
332        type(nwad_dble)             :: s
333        s = (1.5214d0*x + 0.5764d0)/
334     .      (x**2.0d0*(x**2.0d0+1.1284d0*x+0.3183d0))
335      end function
336c
337      end
338C>
339C> \brief The Becke88 exchange GGA enhancement factor
340C>
341#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
342#if defined(NWAD_PRINT)
343      subroutine nwxc_k_becke88_p(tol_rho,x, kalpha)
344#else
345      subroutine nwxc_k_becke88(tol_rho,x, kalpha)
346#endif
347#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
348      subroutine nwxc_k_becke88_d2(tol_rho,x, kalpha)
349#else
350      subroutine nwxc_k_becke88_d3(tol_rho,x, kalpha)
351#endif
352#include "nwad.fh"
353      implicit none
354c
355      double precision tol_rho
356      type(nwad_dble)::x,kalpha
357      double precision dkadxa
358c
359      double precision BETA, C
360      Parameter (BETA = 0.0042D0)
361      type(nwad_dble)::g,gdenom
362      double precision dgdenom,dg
363c     double precision arcsinh, darcsinh
364c     arcsinh(x)=log(x+dsqrt(1d0+x*x))
365c     darcsinh(x)=1d0/dsqrt(1d0+x*x)
366c
367c
368c     Uniform electron gas constant
369c
370      C =  3d0*(0.75d0/acos(-1d0))**(1d0/3d0)
371
372      if (x.gt.0d0)then
373         gdenom = 1d0 + 6d0*BETA*x*asinh(x)
374         g = 2d0*BETA*x*x / gdenom
375c        dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
376c        dg = g*(2d0/x-dgdenom/gdenom)
377
378         kalpha= C + g
379c        dkadxa = dg
380
381      else
382         kalpha= C
383c        dkadxa = 0d0
384      endif
385      return
386      end
387C>
388C> \brief The PBE96 exchange GGA enhancement factor
389C>
390#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
391#if defined(NWAD_PRINT)
392      subroutine nwxc_k_pbe96_p(tol_rho,x, kalpha)
393#else
394      subroutine nwxc_k_pbe96(tol_rho,x, kalpha)
395#endif
396#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
397      subroutine nwxc_k_pbe96_d2(tol_rho,x, kalpha)
398#else
399      subroutine nwxc_k_pbe96_d3(tol_rho,x, kalpha)
400#endif
401#include "nwad.fh"
402      implicit none
403c
404      double precision tol_rho
405      type(nwad_dble)::x,kalpha,deno
406      double precision dkadxa
407c
408      double precision pi,um, uk, umk
409      parameter(um=0.21951d0, uk=0.804d0, umk=um/uk)
410      double precision C
411      double precision forty8
412c
413c
414c     Uniform electron gas constant
415c
416      pi = acos(-1.d0)
417      C =  3d0*(0.75d0/pi)**(1d0/3d0)
418
419      if (x.gt.0d0)then
420         forty8=1d0/((48d0*pi*pi)**(2d0/3d0))
421         deno=1d0/(1d0+um*x*x*forty8/uk)
422         kalpha= C * (1d0 + uk - uk *deno)
423c        dkadxa = C * (2d0*um*x*deno*deno*
424c    *        forty8)
425
426      else
427         kalpha= C
428c        dkadxa = 0d0
429      endif
430      return
431      end
432C>
433C> \brief The Dirac exchange GGA enhancement factor
434C>
435C> Of course the Dirac functional is the exchange part for LDA
436C> so this subroutine just returns a constant.
437C>
438#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
439#if defined(NWAD_PRINT)
440      subroutine nwxc_k_dirac_p(tol_rho,x, kalpha)
441#else
442      subroutine nwxc_k_dirac(tol_rho,x, kalpha)
443#endif
444#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
445      subroutine nwxc_k_dirac_d2(tol_rho,x, kalpha)
446#else
447      subroutine nwxc_k_dirac_d3(tol_rho,x, kalpha)
448#endif
449#include "nwad.fh"
450      implicit none
451c
452      double precision tol_rho
453      type(nwad_dble)::x,kalpha
454      double precision dkadxa
455c
456      double precision pi
457      double precision C
458c
459c
460c     Uniform electron gas constant
461c
462      pi = acos(-1.d0)
463      C =  3d0*(0.75d0/pi)**(1d0/3d0)
464
465      kalpha= C
466c     dkadxa = 0d0
467
468      return
469      end
470#ifndef NWAD_PRINT
471#define NWAD_PRINT
472c
473c     Compile source again for the 2nd derivative case
474c
475#include "nwxc_c_op.F"
476#endif
477#ifndef SECOND_DERIV
478#define SECOND_DERIV
479c
480c     Compile source again for the 2nd derivative case
481c
482#include "nwxc_c_op.F"
483#endif
484#ifndef THIRD_DERIV
485#define THIRD_DERIV
486c
487c     Compile source again for the 3rd derivative case
488c
489#include "nwxc_c_op.F"
490#endif
491#undef NWAD_PRINT
492C>
493C> @}
494