1C> \ingroup nwpwxc
2C> @{
3C>
4C> \file nwpwxc_c_opt.F
5C> The OPT correlation functional
6C>
7C> @}
8C>
9C> \ingroup nwpwxc_priv
10C> @{
11C>
12C> \brief The OPT correlation functional
13C>
14C> The OPTimized correlation functional [1].
15C>
16C> [1] N.C. Handy, A.J. Cohen, "Dynamic correlation",
17C>     Mol. Phys. <b>99</b>, 607-615 (2001), DOI:
18C>     <A HREF="https://doi.org/10.1080/00268970010023435">
19C>     10.1080/00268970010023435</A>.
20c
21C$Id$
22c
23      subroutine nwpwxc_c_opt(tol_rho,ipol,nq,wght,rho,rgamma,func,
24     &                      Amat, Cmat)
25      implicit none
26c
27#include "nwpwxc_param.fh"
28c
29      double precision tol_rho !< [Input] The lower limit on the density
30      integer nq               !< [Input] The number of points
31      integer ipol             !< [Input] The number of spin channels
32      double precision wght    !< [Input] The weight of the functional
33c
34c     Charge Density
35c
36      double precision rho(nq,*) !< [Input] The density
37c
38c     Charge Density Gradient
39c
40      double precision rgamma(nq,*) !< [Input] The norm of the density gradients
41c
42c     Sampling Matrices for the XC Potential & Energy
43c
44      double precision func(nq)   !< [Output] The value of the functional
45      double precision Amat(nq,*) !< [Output] The derivative wrt rho
46      double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma
47c
48c References:
49c
50c    Handy NC, Cohen AJ, Mol Phys 99 (7); 607-615 2001
51c
52      integer l_rho,k_rho,l_delrho,k_delrho
53      double precision c1,c2
54!      parameter (c1=1.1015d0,c2=0.6625d0)
55      parameter (c1=1.d0,c2=0d0)
56c
57      integer iq   ! counter over grid points
58      integer ii   ! offset
59      integer iqt  ! upper limit
60      integer num  ! the number of grid points in current batch
61      integer maxp ! the maximum number of points in a batch
62      parameter (maxp = 16)
63      double precision funcl(maxp)
64      double precision rhol(maxp,2)
65      double precision rgammal(maxp,3)
66      double precision Amatl(maxp,2)
67      double precision Cmatl(maxp,3)
68      double precision fac
69c
70c***************************************************************************
71c
72c
73      if (ipol.eq.1) then
74c
75c        ======> SPIN-RESTRICTED <======
76c
77         do iq = 1, nq, maxp
78           iqt = min(nq,iq+maxp-1)
79           num = iqt-iq+1
80c
81c          c1*Ec[a,b]
82c
83           call dfill(maxp,0.0d0,funcl,1)
84           call dfill(maxp*2,0.0d0,rhol,1)
85           call dfill(maxp*3,0.0d0,rgammal,1)
86           call dfill(maxp*2,0.0d0,Amatl,1)
87           call dfill(maxp*3,0.0d0,Cmatl,1)
88           do ii = 0, num - 1
89              rhol(ii+1,R_T) = rho(iq+ii,R_T)
90              rgammal(ii+1,G_TT) = rgamma(iq+ii,G_TT)
91           enddo
92           fac = c1*wght
93           call nwpwxc_c_pw91lda(tol_rho,1,maxp,fac,rhol,funcl,Amatl)
94           call nwpwxc_c_p91(tol_rho,1,maxp,fac,rhol,rgammal,funcl,
95     +                     Amatl,Cmatl)
96           do ii = 0, num - 1
97              func(iq+ii) = func(iq+ii) + funcl(ii+1)
98              Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA)
99              Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA)
100     +                           + Cmatl(ii+1,D1_GAA)
101              Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB)
102     +                           + Cmatl(ii+1,D1_GAB)
103           enddo
104c
105c          {(c2-c1)*Ec[a,0]}(a) + {(c2-c1)*Ec[0,b]}(a)
106c
107           fac = (c2-c1)*wght
108           call dfill(maxp,0.0d0,funcl,1)
109           call dfill(maxp*2,0.0d0,Amatl,1)
110           call dfill(maxp*3,0.0d0,Cmatl,1)
111           do ii = 0, num - 1
112              rhol(ii+1,R_A) = 0.5d0*rhol(ii+1,R_T)
113              rgammal(ii+1,G_AA) = 0.25d0*rgammal(ii+1,G_TT)
114           enddo
115           call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl)
116           call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl,
117     +                     Amatl,Cmatl)
118           do ii = 0, num - 1
119              func(iq+ii) = func(iq+ii) + 2*funcl(ii+1)
120              Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA)
121     +                          + Amatl(ii+1,D1_RA) + Amatl(ii+1,D1_RB)
122              Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA)
123     +                           + (2.0d0)*Cmatl(ii+1,D1_GAA)
124              Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB)
125     +                           + (2.0d0)*Cmatl(ii+1,D1_GAB)
126           enddo
127         enddo
128c
129      else
130c
131c        ======> SPIN-UNRESTRICTED <======
132c
133         do iq = 1, nq, maxp
134           iqt = min(nq,iq+maxp-1)
135           num = iqt-iq+1
136c
137c          c1*Ec[a,b]
138c
139           call dfill(maxp,0.0d0,funcl,1)
140           call dfill(maxp*2,0.0d0,rhol,1)
141           call dfill(maxp*3,0.0d0,rgammal,1)
142           call dfill(maxp*2,0.0d0,Amatl,1)
143           call dfill(maxp*3,0.0d0,Cmatl,1)
144           do ii = 0, num - 1
145              rhol(ii+1,R_A) = rho(iq+ii,R_A)
146              rhol(ii+1,R_B) = rho(iq+ii,R_B)
147              rgammal(ii+1,G_AA) = rgamma(iq+ii,G_AA)
148              rgammal(ii+1,G_AB) = rgamma(iq+ii,G_AB)
149              rgammal(ii+1,G_BB) = rgamma(iq+ii,G_BB)
150           enddo
151           fac = c1*wght
152           call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl)
153           call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl,
154     +                     Amatl,Cmatl)
155           do ii = 0, num - 1
156              func(iq+ii) = func(iq+ii) + funcl(ii+1)
157              Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA)
158              Amat(iq+ii,D1_RB) = Amat(iq+ii,D1_RB) + Amatl(ii+1,D1_RB)
159              Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA)
160     +                           + Cmatl(ii+1,D1_GAA)
161              Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB)
162     +                           + Cmatl(ii+1,D1_GAB)
163              Cmat(iq+ii,D1_GBB) = Cmat(iq+ii,D1_GBB)
164     +                           + Cmatl(ii+1,D1_GBB)
165           enddo
166c
167c          (c2-c1)*Ec[a,0]
168c
169           fac = (c2-c1)*wght
170           call dfill(maxp  ,0.0d0,funcl,1)
171           call dfill(maxp*2,0.0d0,Amatl,1)
172           call dfill(maxp*3,0.0d0,Cmatl,1)
173           do ii = 0, num - 1
174              rhol(ii+1,R_A) = rho(iq+ii,R_A)
175              rhol(ii+1,R_B) = 0.0d0
176              rgammal(ii+1,G_AA) = rgamma(iq+ii,G_AA)
177              rgammal(ii+1,G_AB) = 0.0d0
178              rgammal(ii+1,G_BB) = 0.0d0
179           enddo
180           call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl)
181           call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl,
182     +                     Amatl,Cmatl)
183           do ii = 0, num - 1
184              func(iq+ii) = func(iq+ii) + funcl(ii+1)
185              Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA)
186              Amat(iq+ii,D1_RB) = Amat(iq+ii,D1_RB) + Amatl(ii+1,D1_RB)
187              Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA)
188     +                           + Cmatl(ii+1,D1_GAA)
189              Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB)
190     +                           + Cmatl(ii+1,D1_GAB)
191              Cmat(iq+ii,D1_GBB) = Cmat(iq+ii,D1_GBB)
192     +                           + Cmatl(ii+1,D1_GBB)
193           enddo
194c
195c          (c2-c1)*Ec[0,b]
196c
197           fac = (c2-c1)*wght
198           call dfill(maxp  ,0.0d0,funcl,1)
199           call dfill(maxp*2,0.0d0,Amatl,1)
200           call dfill(maxp*3,0.0d0,Cmatl,1)
201           do ii = 0, num - 1
202              rhol(ii+1,R_A) = 0.0d0
203              rhol(ii+1,R_B) = rho(iq+ii,R_B)
204              rgammal(ii+1,G_AA) = 0.0d0
205              rgammal(ii+1,G_AB) = 0.0d0
206              rgammal(ii+1,G_BB) = rgamma(iq+ii,G_BB)
207           enddo
208           call nwpwxc_c_pw91lda(tol_rho,2,maxp,fac,rhol,funcl,Amatl)
209           call nwpwxc_c_p91(tol_rho,2,maxp,fac,rhol,rgammal,funcl,
210     +                     Amatl,Cmatl)
211           do ii = 0, num - 1
212              func(iq+ii) = func(iq+ii) + funcl(ii+1)
213              Amat(iq+ii,D1_RA) = Amat(iq+ii,D1_RA) + Amatl(ii+1,D1_RA)
214              Amat(iq+ii,D1_RB) = Amat(iq+ii,D1_RB) + Amatl(ii+1,D1_RB)
215              Cmat(iq+ii,D1_GAA) = Cmat(iq+ii,D1_GAA)
216     +                           + Cmatl(ii+1,D1_GAA)
217              Cmat(iq+ii,D1_GAB) = Cmat(iq+ii,D1_GAB)
218     +                           + Cmatl(ii+1,D1_GAB)
219              Cmat(iq+ii,D1_GBB) = Cmat(iq+ii,D1_GBB)
220     +                           + Cmatl(ii+1,D1_GBB)
221           enddo
222         enddo
223c
224      endif
225c
226      return
227      end
228C>
229C> @}
230