1*     ************************************************
2*     *                                              *
3*     *                 nwpw_m06_x                   *
4*     *                                              *
5*     ************************************************
6      subroutine nwpw_m06_x(pi,thrd,fvthrd,etthrd,
7     >                      a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,
8     >                      C1,C2,Cx,P23,
9     >                      n,agr,tau,
10     >                      xe,dfdnx,dfdagrx,dfdtaux)
11      implicit none
12*     ***** input *****
13      real*8 pi,thrd,fvthrd,etthrd
14      real*8 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11
15      real*8 C1,C2,Cx,P23
16      real*8 n,agr,tau
17*     ***** output *****
18      real*8 xe,dfdnx,dfdagrx,dfdtaux
19*     ***** local declarations *****
20      real*8 n_13,n_23,n_53,n_83,agr2
21      real*8 tauU,t,dt_dn,dt_dtau
22      real*8 w,dw_dt,w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11
23      real*8 fw,dfw_dw,dfw_dn,dfw_dtau
24      real*8 x,dx_dn,dx_dagr
25      real*8 enum,eden,eden2,etmp,detmp_dx
26      real*8 Fpbe,dFpbe_dn,dFpbe_dagr
27
28      n_13  = n**thrd
29      n_23  = n_13*n_13
30      n_53  = n_23*n
31      n_83  = n_53*n
32      agr2  = agr*agr
33
34      tauU    = P23*n_53
35      t       = tauU/tau
36      dt_dn   = fvthrd*t/n
37      dt_dtau = -t/tau
38
39      w     = (t - 1.0d0)/(t + 1.0d0)
40      dw_dt = 2.0d0/((1.0d0 + t)**2.0d0)
41
42      w1  = w
43      w2  = w1*w
44      w3  = w2*w
45      w4  = w3*w
46      w5  = w4*w
47      w6  = w5*w
48      w7  = w6*w
49      w8  = w7*w
50      w9  = w8*w
51      w10 = w9*w
52      w11 = w10*w
53
54      fw       = a0    + a1*w1 + a2*w2 + a3*w3 + a4*w4   + a5*w5
55     &         + a6*w6 + a7*w7 + a8*w8 + a9*w9 + a10*w10 + a11*w11
56      dfw_dw   = a1            + 2.0d0*a2*w1     + 3.0d0*a3*w2
57     &         + 4.0d0*a4*w3   + 5.0d0*a5*w4     + 6.0d0*a6*w5
58     &         + 7.0d0*a7*w6   + 8.0d0*a8*w7     + 9.0d0*a9*w8
59     &         + 10.0d0*a10*w9 + 11.0d0*a11*w10
60      dfw_dn   = dfw_dw*dw_dt*dt_dn
61      dfw_dtau = dfw_dw*dw_dt*dt_dtau
62
63      x       =  agr2/n_83
64      dx_dn   = -etthrd*x/n
65      dx_dagr =  2.0d0*x/agr
66
67      enum     =  C1*x
68      eden     =  1.0d0 + C2*x
69      eden2    =  eden*eden
70      etmp     =  -enum/eden
71      detmp_dx =  -(C1*eden - C2*enum)/eden2
72
73      Fpbe       = n_13*(Cx + etmp)
74      dFpbe_dn   = n_13*detmp_dx*dx_dn   + thrd*Fpbe/n
75      dFpbe_dagr = n_13*detmp_dx*dx_dagr
76
77      xe      = Fpbe*fw
78      dfdnx   = n*(dFpbe_dn*fw + Fpbe*dfw_dn) + xe
79      dfdagrx = n*dFpbe_dagr*fw
80      dfdtaux = n*Fpbe*dfw_dtau
81
82      return
83      end
84*     ************************************************
85*     *                                              *
86*     *                nwpw_m06_c_ss                 *
87*     *                                              *
88*     ************************************************
89      subroutine nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4,
90     >                         n,xs,dxs_dn,dxs_dagr,zs,dzs_dn,dzs_dtau,
91     >                         ess0c,dess0c_dn,
92     >                         cess,dfss_dn,dfss_dagr,dfss_dtau)
93      implicit none
94*     ***** input *****
95      real*8 cgss,css0,css1,css2,css3,css4
96      real*8 n,xs,dxs_dn,dxs_dagr,zs,dzs_dn,dzs_dtau
97*     ***** output *****
98      real*8 ess0c,dess0c_dn
99      real*8 cess,dfss_dn,dfss_dagr,dfss_dtau
100*     ***** local declarations *****
101      real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau
102      real*8 U,dU_dx,U2,U3,U4,gss,dgss_dU
103      real*8 rs,drs_dn
104      real*8 pwc,dpwc_drs,dummy,ness0c
105*     ***** constants *****
106      real*8 pi,thrd
107      parameter (pi     = 3.14159265358979311599d0)
108      parameter (thrd   = 1.0d0/3.0d0)
109*     ***** density cutoff parameters *****
110      real*8 tol
111      parameter (tol = 1.0d-18)
112*     ***** M06 constants *****
113      real*8 cf
114      parameter (cf = 9.115599720d0)
115
116      rs = (0.75d0/(pi*n))**thrd
117      drs_dn = -thrd*rs/n
118
119      call gen_PW91_c_rz(tol,rs,1.0d0,pwc,dpwc_drs,dummy)
120
121      ess0c = pwc
122      ness0c = n*ess0c
123      dess0c_dn = dpwc_drs*drs_dn
124
125      U = cgss*xs/(1.0d0 + cgss*xs)
126      dU_dx = cgss/((1.0d0 + cgss*xs)**2.0d0)
127
128      U2 = U*U
129      U3 = U2*U
130      U4 = U3*U
131
132      gss = css0 + css1*U + css2*U2 + css3*U3 + css4*U4
133      dgss_dU = css1 + 2.0d0*css2*U + 3.0d0*css3*U2 + 4.0d0*css4*U3
134
135      D = 1.0d0 - 0.25d0*xs/(zs+cf)
136      if (D .lt. 0.0d0) then
137        D       = 0.0d0
138        dD_dn   = 0.0d0
139        dD_dagr = 0.0d0
140        dD_dtau = 0.0d0
141      else
142        dD_dx   = -1.0d0/(4.0d0*(zs + cf))
143        dD_dz   =  xs/(4.0d0*(zs + cf)*(zs + cf))
144        dD_dn   =  dD_dx*dxs_dn + dD_dz*dzs_dn
145        dD_dagr =  dD_dx*dxs_dagr
146        dD_dtau =  dD_dz*dzs_dtau
147      end if
148
149      cess = ess0c*gss*D
150      dfss_dn = cess + n*dess0c_dn*gss*D
151     &        + ness0c*dgss_dU*dU_dx*dxs_dn*D + ness0c*gss*dD_dn
152      dfss_dagr = ness0c*(dgss_dU*dU_dx*dxs_dagr*D + gss*dD_dagr)
153      dfss_dtau = ness0c*gss*dD_dtau
154
155      return
156      end
157
158*     ************************************************
159*     *                                              *
160*     *              nwpw_m06_c_restricted           *
161*     *                                              *
162*     ************************************************
163      subroutine nwpw_m06_c_restricted(cgss,css0,css1,css2,css3,css4,
164     >                             cgopp,copp0,copp1,copp2,copp3,copp4,
165     >                             n,nup,xup,dxup_dnup,dxup_dagrup,
166     >                             zup,dzup_dnup,dzup_dtauup,rs,drs_dn,
167     >                             ce,fnc,fdnc,fdtauc)
168      implicit none
169*     ***** input *****
170      real*8 cgss,css0,css1,css2,css3,css4
171      real*8 cgopp,copp0,copp1,copp2,copp3,copp4
172      real*8 n,nup,xup,dxup_dnup,dxup_dagrup
173      real*8 zup,dzup_dnup,dzup_dtauup
174      real*8 rs,drs_dn
175*     ***** output *****
176      real*8 ce,fnc,fdnc,fdtauc
177*     ***** local declarations *****
178      real*8 euu0c,deuu0c_dn
179      real*8 euuc,dfuuc_dn,dfuuc_dagr,dfuuc_dtau
180      real*8 eud0c,deud0c_drs,dummy,eud1c,fud1c,dfud1c_dn
181      real*8 x,U,dU_dx,U2,U3,U4,gopp,dgopp_dU
182      real*8 eudc,dfudc_dn,dfudc_dagr,dfudc_dtau
183*     ***** density cutoff parameters *****
184      real*8 tol
185      parameter (tol = 1.0d-18)
186
187*     ***** same spin (alpha-alpha or beta-beta) *****
188      call nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4,
189     >                   nup,xup,dxup_dnup,dxup_dagrup,
190     >                   zup,dzup_dnup,dzup_dtauup,
191     >                   euu0c,deuu0c_dn,
192     >                   euuc,dfuuc_dn,dfuuc_dagr,dfuuc_dtau)
193
194*     ***** opposite spin (alpha-beta) *****
195      call gen_PW91_c_rz(tol,rs,0.0d0,eud0c,deud0c_drs,dummy)
196      eud1c = eud0c - euu0c
197      fud1c = n*eud1c
198      dfud1c_dn  = eud0c + n*deud0c_drs*drs_dn
199     &           - euu0c - nup*deuu0c_dn
200
201      x   = 2.0d0*xup
202      U = cgopp*x/(1.0d0 + cgopp*x)
203      dU_dx = cgopp/((1.0d0 + cgopp*x)**2.0d0)
204
205      U2 = U*U
206      U3 = U2*U
207      U4 = U3*U
208
209      gopp = copp0 + copp1*U + copp2*U2 + copp3*U3 + copp4*U4
210      dgopp_dU = copp1          + 2.0d0*copp2*U
211     &         + 3.0d0*copp3*U2 + 4.0d0*copp4*U3
212
213      eudc = eud1c*gopp
214      dfudc_dn = dfud1c_dn*gopp + fud1c*dgopp_dU*dU_dx*dxup_dnup
215      dfudc_dagr = fud1c*dgopp_dU*dU_dx*dxup_dagrup
216      dfudc_dtau = 0.0d0
217
218      ce     = eudc        + euuc
219      fnc    = dfudc_dn    + dfuuc_dn
220      fdnc   = dfudc_dagr  + dfuuc_dagr
221      fdtauc =               dfuuc_dtau
222
223      return
224      end
225
226*     ************************************************
227*     *                                              *
228*     *            nwpw_m06_c_unrestricted           *
229*     *                                              *
230*     ************************************************
231      subroutine nwpw_m06_c_unrestricted(cgss,css0,css1,css2,css3,css4,
232     >                             cgopp,copp0,copp1,copp2,copp3,copp4,
233     >                             n,nup,agrup,tauup,ndn,agrdn,taudn,
234     >              xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup,
235     >              xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn,
236     >                   rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn,
237     >                   ce,fnupc,fndnc,fdnupc,fdndnc,fdtauupc,fdtaudnc)
238      implicit none
239*     ***** input *****
240      real*8 cgss,css0,css1,css2,css3,css4
241      real*8 cgopp,copp0,copp1,copp2,copp3,copp4
242      real*8 n,nup,agrup,tauup,ndn,agrdn,taudn
243      real*8 xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup
244      real*8 xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn
245      real*8 rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn
246*     ***** output *****
247      real*8 ce,fnupc,fndnc,fdnupc,fdndnc,fdtauupc,fdtaudnc
248*     ***** local declarations *****
249      real*8 euu0c,deuu0c_dnup
250      real*8 euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup
251      real*8 edd0c,dedd0c_dndn
252      real*8 eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn
253      real*8 eud0c,deud0c_drs,deud0c_dzeta
254      real*8 x,U,dU_dx,U2,U3,U4,gopp,dgopp_dU
255      real*8 eudc,eud1c,fud1c,dfud1c_dnup,dfud1c_dndn
256      real*8 dfudc_dnup,dfudc_dagrup,dfudc_dtauup
257      real*8 dfudc_dndn,dfudc_dagrdn,dfudc_dtaudn
258*     ***** density cutoff parameters *****
259      real*8 tol
260      parameter (tol = 1.0d-18)
261
262*     ***** same spin *****
263*     ***** alpha-alpha *****
264      call nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4,
265     >                   nup,xup,dxup_dnup,dxup_dagrup,
266     >                   zup,dzup_dnup,dzup_dtauup,
267     >                   euu0c,deuu0c_dnup,
268     >                   euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup)
269*     ***** beta-beta *****
270      call nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4,
271     >                   ndn,xdn,dxdn_dndn,dxdn_dagrdn,
272     >                   zdn,dzdn_dndn,dzdn_dtaudn,
273     >                   edd0c,dedd0c_dndn,
274     >                   eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn)
275
276*     ***** opposite spin (alpha-beta) *****
277      call gen_PW91_c_rz(tol,rs,zeta,eud0c,deud0c_drs,deud0c_dzeta)
278      eud1c = eud0c - (euu0c*nup + edd0c*ndn)/n
279      fud1c = n*eud1c
280
281      x   = xup + xdn
282      U = cgopp*x/(1.0d0 + cgopp*x)
283      dU_dx = cgopp/((1.0d0 + cgopp*x)**2.0d0)
284
285      U2 = U*U
286      U3 = U2*U
287      U4 = U3*U
288
289      gopp = copp0 + copp1*U + copp2*U2 + copp3*U3 + copp4*U4
290      dgopp_dU = copp1          + 2.0d0*copp2*U
291     &         + 3.0d0*copp3*U2 + 4.0d0*copp4*U3
292
293      eudc = eud1c*gopp
294      dfud1c_dnup = eud0c + n*(deud0c_drs*drs_dn
295     &            + deud0c_dzeta*dzeta_dnup) - euu0c
296     &            - nup*deuu0c_dnup
297      dfud1c_dndn = eud0c + n*(deud0c_drs*drs_dn
298     &            + deud0c_dzeta*dzeta_dndn) - edd0c
299     &            - ndn*dedd0c_dndn
300
301      dfudc_dnup = dfud1c_dnup*gopp + fud1c*dgopp_dU*dU_dx*dxup_dnup
302      dfudc_dndn = dfud1c_dndn*gopp + fud1c*dgopp_dU*dU_dx*dxdn_dndn
303      dfudc_dagrup = fud1c*dgopp_dU*dU_dx*dxup_dagrup
304      dfudc_dagrdn = fud1c*dgopp_dU*dU_dx*dxdn_dagrdn
305      dfudc_dtauup = 0.0d0
306      dfudc_dtaudn = 0.0d0
307
308      ce       = eudc          + (euuc*nup  + eddc*ndn)/n
309      fnupc    = dfudc_dnup    + dfuuc_dnup
310      fndnc    = dfudc_dndn    + dfddc_dndn
311      fdnupc   = dfudc_dagrup  + dfuuc_dagrup
312      fdndnc   = dfudc_dagrdn  + dfddc_dagrdn
313      fdtauupc =                 dfuuc_dtauup
314      fdtaudnc =                 dfddc_dtaudn
315
316      return
317      end
318