1C> \ingroup nwpwxc
2C> @{
3C>
4C> \file nwpwxc_c_b97.F
5C> The B97 family of correlation functionals
6C>
7C> @}
8#include "nwpwxcP_xc_b97.h"
9C>
10C> \ingroup nwpwxc_priv
11C> @{
12C>
13C> \brief Evaluate the B97 family of correlation functionals
14C>
15C> This code evaluates correlation functionals from the
16C> B97 family of functionals [1,2].
17C>
18C> ### References ###
19C>
20C> [1] A.D. Becke, "Density-functional thermochemistry. V. Systematic
21C>     optimization of exchange-correlation functionals", J. Chem. Phys.
22C>     107 (1997) 8554-8560, DOI:
23C>     <a href="https://doi.org/10.1063/1.475007">
24C>     10.1063/1.475007</a>.
25C>
26C> [2] S. Grimme, "Semiempirical GGA-type density functional constructed
27C>     with a long-range dispersion correction", J. Comput. Chem. 27
28C>     (2006) 1787-1799, DOI:
29C>     <a href="https://doi.org/10.1002/jcc.20495">
30C>     10.1002/jcc.20495</a>.
31C>
32      Subroutine nwpwxc_c_b97(param,tol_rho,ipol,nq,wght,rho,rgamma,
33     &                      func,Amat,Cmat)
34c
35c     $Id$
36c
37      implicit none
38c
39#include "nwpwxc_param.fh"
40c
41c     Input and other parameters
42c
43      double precision param(*)!< [Input] Parameters of functional as
44                               !< defined in [1]:
45                               !< - param(1): \f$m\f$ of Eqs.(20).
46                               !< - param(2): \f$C_{C\sigma\sigma,0}\f$
47                               !< - param(3): \f$C_{C\alpha\beta,0}\f$
48                               !< - param(4): \f$C_{C\sigma\sigma,1}\f$
49                               !< - param(5): \f$C_{C\alpha\beta,1}\f$
50                               !< - param(6): \f$C_{C\sigma\sigma,2}\f$
51                               !< - param(7): \f$C_{C\alpha\beta,2}\f$
52                               !< - param(8): \f$C_{C\sigma\sigma,3}\f$
53                               !< - param(9): \f$C_{C\alpha\beta,3}\f$
54                               !< - param(10): \f$C_{C\sigma\sigma,4}\f$
55                               !< - param(11): \f$C_{C\alpha\beta,4}\f$
56      double precision tol_rho !< [Input] The lower limit on the density
57      integer ipol             !< [Input] The number of spin channels
58      integer nq               !< [Input] The number of points
59      double precision wght    !< [Input] The weight of the functional
60c
61c     Charge Density
62c
63      double precision rho(nq,*)    !< [Input] The density
64c
65c     Charge Density Gradient
66c
67      double precision rgamma(nq,*) !< [Input] The norm of the density gradients
68c
69c     Sampling Matrices for the XC Potential
70c
71      double precision func(nq)     !< [Output] The value of the functional
72      double precision Amat(nq,*)   !< [Output] The derivative wrt rho
73      double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
74c
75c     Local variables
76c
77      integer i
78      double precision rho_a(0:1)
79      double precision rho_b(0:1)
80      double precision FC(0:_FXC_NUMDERI)
81c
82c     Code
83c
84      if (ipol.eq.1) then
85        do i = 1, nq
86          rho_a(0) = rho(i,R_T)*0.5d0
87          rho_b(0) = rho_a(0)
88          rho_a(1) = rgamma(i,G_TT)*0.25d0
89          rho_b(1) = rho_a(1)
90          if (rho_a(0).gt.tol_rho) then
91            call nwpwxcp_c_b97(rho_a,rho_b,wght,tol_rho,FC,param)
92            func(i)        = func(i)        + FC(_FXC_E)*wght
93            Amat(i,D1_RA)  = Amat(i,D1_RA)  + FC(_FXC_RA)*wght
94            Cmat(i,D1_GAA) = Cmat(i,D1_GAA) + FC(_FXC_GAA)*wght
95          endif
96        enddo
97      else
98        do i = 1, nq
99          rho_a(0) = rho(i,R_A)
100          rho_b(0) = rho(i,R_B)
101          rho_a(1) = rgamma(i,G_AA)
102          rho_b(1) = rgamma(i,G_BB)
103          if (rho_a(0).gt.tol_rho.or.rho_b(0).gt.tol_rho) then
104            call nwpwxcp_c_b97(rho_a,rho_b,wght,tol_rho,FC,param)
105            func(i)        = func(i)        + FC(_FXC_E)*wght
106            Amat(i,D1_RA)  = Amat(i,D1_RA)  + FC(_FXC_RA)*wght
107            Amat(i,D1_RB)  = Amat(i,D1_RB)  + FC(_FXC_RB)*wght
108            Cmat(i,D1_GAA) = Cmat(i,D1_GAA) + FC(_FXC_GAA)*wght
109            Cmat(i,D1_GBB) = Cmat(i,D1_GBB) + FC(_FXC_GBB)*wght
110          endif
111        enddo
112      endif
113c
114      end
115c
116      Subroutine nwpwxc_c_b97_d2(param,tol_rho,ipol,nq,wght,rho,rgamma,
117     &                         func,Amat,Amat2,Cmat,Cmat2)
118c
119c     $Id$
120c
121      implicit none
122c
123#include "nwpwxc_param.fh"
124c
125c     Input and other parameters
126c
127      double precision param(*)!< [Input] Parameters of functional as
128                               !< defined in [1]:
129                               !< - param(1): \f$m\f$ of Eqs.(20).
130                               !< - param(2): \f$C_{C\sigma\sigma,0}\f$
131                               !< - param(3): \f$C_{C\alpha\beta,0}\f$
132                               !< - param(4): \f$C_{C\sigma\sigma,1}\f$
133                               !< - param(5): \f$C_{C\alpha\beta,1}\f$
134                               !< - param(6): \f$C_{C\sigma\sigma,2}\f$
135                               !< - param(7): \f$C_{C\alpha\beta,2}\f$
136                               !< - param(8): \f$C_{C\sigma\sigma,3}\f$
137                               !< - param(9): \f$C_{C\alpha\beta,3}\f$
138                               !< - param(10): \f$C_{C\sigma\sigma,4}\f$
139                               !< - param(11): \f$C_{C\alpha\beta,4}\f$
140      double precision tol_rho !< [Input] The lower limit on the density
141      integer ipol             !< [Input] The number of spin channels
142      integer nq               !< [Input] The number of points
143      double precision wght    !< [Input] The weight of the functional
144c
145c     Charge Density
146c
147      double precision rho(nq,*)    !< [Input] The density
148c
149c     Charge Density Gradient
150c
151      double precision rgamma(nq,*) !< [Input] The norm of the density gradients
152c
153c     Sampling Matrices for the XC Potential
154c
155      double precision func(nq)     !< [Output] The value of the functional
156      double precision Amat(nq,*)   !< [Output] The derivative wrt rho
157      double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
158      double precision Amat2(nq,*)  !< [Output] The 2nd derivative wrt rho
159      double precision Cmat2(nq,*)  !< [Output] The 2nd derivative wrt
160                                    !< rho and rgamma
161c
162c     Local variables
163c
164      integer i
165      double precision rho_a(0:1)
166      double precision rho_b(0:1)
167      double precision FC(0:_FXC_NUMDERI)
168c
169c     Code
170c
171      if (ipol.eq.1) then
172        do i = 1, nq
173          rho_a(0) = rho(i,R_T)*0.5d0
174          rho_b(0) = rho_a(0)
175          rho_a(1) = rgamma(i,G_TT)*0.25d0
176          rho_b(1) = rho_a(1)
177          if (rho_a(0).gt.tol_rho) then
178            call nwpwxcp_c_b97(rho_a,rho_b,wght,tol_rho,FC,param)
179            func(i)        = func(i)        + FC(_FXC_E)*wght
180            Amat(i,D1_RA)  = Amat(i,D1_RA)  + FC(_FXC_RA)*wght
181            Cmat(i,D1_GAA) = Cmat(i,D1_GAA) + FC(_FXC_GAA)*wght
182            Amat2(i,D2_RA_RA) = Amat2(i,D2_RA_RA) + FC(_FXC_RARA)*wght
183            Amat2(i,D2_RA_RB) = Amat2(i,D2_RA_RB) + FC(_FXC_RARB)*wght
184            Cmat2(i,D2_RA_GAA) = Cmat2(i,D2_RA_GAA)
185     &                         + FC(_FXC_RAGAA)*wght
186            Cmat2(i,D2_RA_GBB) = Cmat2(i,D2_RA_GBB)
187     &                         + FC(_FXC_RAGBB)*wght
188            Cmat2(i,D2_GAA_GAA) = Cmat2(i,D2_GAA_GAA)
189     &                          + FC(_FXC_GAAGAA)*wght
190            Cmat2(i,D2_GAA_GBB) = Cmat2(i,D2_GAA_GBB)
191     &                          + FC(_FXC_GAAGBB)*wght
192          endif
193        enddo
194      else
195        do i = 1, nq
196          rho_a(0) = rho(i,R_A)
197          rho_b(0) = rho(i,R_B)
198          rho_a(1) = rgamma(i,G_AA)
199          rho_b(1) = rgamma(i,G_BB)
200          if (rho_a(0).gt.tol_rho.or.rho_b(0).gt.tol_rho) then
201            call nwpwxcp_c_b97(rho_a,rho_b,wght,tol_rho,FC,param)
202            func(i)        = func(i)        + FC(_FXC_E)*wght
203            Amat(i,D1_RA)  = Amat(i,D1_RA)  + FC(_FXC_RA)*wght
204            Amat(i,D1_RB)  = Amat(i,D1_RB)  + FC(_FXC_RB)*wght
205            Cmat(i,D1_GAA) = Cmat(i,D1_GAA) + FC(_FXC_GAA)*wght
206            Cmat(i,D1_GBB) = Cmat(i,D1_GBB) + FC(_FXC_GBB)*wght
207            Amat2(i,D2_RA_RA) = Amat2(i,D2_RA_RA) + FC(_FXC_RARA)*wght
208            Amat2(i,D2_RA_RB) = Amat2(i,D2_RA_RB) + FC(_FXC_RARB)*wght
209            Amat2(i,D2_RB_RB) = Amat2(i,D2_RB_RB) + FC(_FXC_RBRB)*wght
210            Cmat2(i,D2_RA_GAA) = Cmat2(i,D2_RA_GAA)
211     &                         + FC(_FXC_RAGAA)*wght
212            Cmat2(i,D2_RA_GBB) = Cmat2(i,D2_RA_GBB)
213     &                         + FC(_FXC_RAGBB)*wght
214            Cmat2(i,D2_RB_GAA) = Cmat2(i,D2_RB_GAA)
215     &                         + FC(_FXC_RBGAA)*wght
216            Cmat2(i,D2_RB_GBB) = Cmat2(i,D2_RB_GBB)
217     &                         + FC(_FXC_RBGBB)*wght
218            Cmat2(i,D2_GAA_GAA) = Cmat2(i,D2_GAA_GAA)
219     &                          + FC(_FXC_GAAGAA)*wght
220            Cmat2(i,D2_GAA_GBB) = Cmat2(i,D2_GAA_GBB)
221     &                          + FC(_FXC_GAAGBB)*wght
222            Cmat2(i,D2_GBB_GBB) = Cmat2(i,D2_GBB_GBB)
223     &                          + FC(_FXC_GBBGBB)*wght
224          endif
225        enddo
226      endif
227c
228      end
229c
230      subroutine nwpwxcp_c_pwlda(ra,rb,FCLDA)
231      implicit none
232c
233      double precision ra
234      double precision rb
235      double precision FCLDA(0:_FCLDA_ELEMENTS)
236c
237      double precision ec
238      double precision rho(2)
239      double precision Amat(2)
240      double precision Amat2(3)
241c
242      rho(R_A)        = ra
243      rho(R_B)        = rb
244      ec              = 0.0d0
245      Amat(D1_RA)     = 0.0d0
246      Amat(D1_RB)     = 0.0d0
247      Amat2(D2_RA_RA) = 0.0d0
248      Amat2(D2_RA_RB) = 0.0d0
249      Amat2(D2_RB_RB) = 0.0d0
250c
251      call nwpwxc_c_pw91lda_d2(1.0d-20,2,1,1.0d0,rho,ec,Amat,Amat2)
252c
253      FCLDA(_FXC_E)    = ec
254      FCLDA(_FXC_RA)   = Amat(D1_RA)
255      FCLDA(_FXC_RB)   = Amat(D1_RB)
256      FCLDA(_FXC_RARA) = Amat2(D2_RA_RA)
257      FCLDA(_FXC_RARB) = Amat2(D2_RA_RB)
258      FCLDA(_FXC_RBRB) = Amat2(D2_RB_RB)
259c
260      end
261C>
262C> @}
263