1*
2* $Id$
3*
4
5***  a = 0.2 am=0.8
6**** a*HF + (1-a)*Dirac - 0.72*Dirac + 0.72*Becke88 + 0.81*LYP + (1-0.81)*Vosko
7**** a*HF + (1-a-0.72)*Dirac + 0.72*Becke88 + 0.81*LYP + 0.19*Vosko
8**** a*HF + (am/0.8)*((0.28-0.2)*Dirac + 0.72*Becke88) + 0.81*LYP + 0.19*Vosko
9**** a*HF + (am/0.8)*(0.08*Dirac + 0.72*Becke88) + 0.81*LYP + 0.19*Vosko
10**** a*HF + (am)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko
11**** a*HF + (1-a)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko
12
13*     **************************************************
14*     *                                                *
15*     *           gen_B3LYP_BW_unrestricted             *
16*     *                                                *
17*     **************************************************
18
19      subroutine gen_B3LYP_BW_unrestricted(n2ft3d,
20     &                          dn_in,agr_in,
21     &                          x_parameter,c_parameter,
22     &                          xce,fn,fdn)
23
24      implicit none
25
26      integer   n2ft3d
27      real*8    dn_in(n2ft3d,2)
28      real*8    agr_in(n2ft3d,3)
29      real*8    xce(n2ft3d)
30      real*8    fn(n2ft3d,2)
31      real*8    fdn(n2ft3d,3)
32      real*8    x_parameter, c_parameter
33
34
35c*---- parameters given by vosko et al -----------------*
36      real*8 ap,af,x0p,x0f,bp,bf,cpt,cft
37      parameter (ap = 3.109070d-02, af = 1.554530d-02)
38      parameter (x0p=-1.049800d-01, x0f=-3.250000d-01)
39      parameter (bp = 3.727440d+00, bf = 7.060420d+00)
40      parameter (cpt = 1.293520d+01, cft = 1.805780d+01)
41c*------------------------------------------------------*
42
43*     constants calculated from vosko's parameters
44      real*8 xp,xf,qp,qf,xx0p,xx0f,fc1,fd1,crs
45      real*8 cp1,cp2,cp3,cp4,cp5,cp6,dp1,dp2,dp3,dp4,dp5,dp6,dp7
46      real*8 cf1,cf2,cf3,cf4,cf5,cf6,df1,df2,df3,df4,df5,df6,df7
47      parameter (xp  = -4.581653d-01,  xf  = -5.772521d-01)
48      parameter (qp  =  6.151991d+00,  qf  =  4.730927d+00)
49      parameter (xx0p=  1.255491d+01,  xx0f=  1.586879d+01)
50      parameter (cp1 =  3.109070d-02,  cf1 =  1.554530d-02)
51      parameter (cp2 =  9.690228d-04,  cf2 =  2.247860d-03)
52      parameter (cp3 =  1.049800d-01,  cf3 =  3.250000d-01)
53      parameter (cp4 =  3.878329d-02,  cf4 =  5.249122d-02)
54      parameter (cp5 =  3.075995d+00,  cf5 =  2.365463d+00)
55      parameter (cp6 =  1.863720d+00,  cf6 =  3.530210d+00)
56      parameter (dp1 =  6.218140d-02,  df1 =  3.109060d-02)
57      parameter (dp2 =  1.938045d-03,  df2 =  4.495720d-03)
58      parameter (dp3 =  1.049800d-01,  df3 =  3.250000d-01)
59      parameter (dp4 = -3.205972d-02,  df4 = -1.779316d-02)
60      parameter (dp5 = -1.192972d-01,  df5 = -1.241661d-01)
61      parameter (dp6 =  1.863720d+00,  df6 =  3.530210d+00)
62      parameter (dp7 =  9.461748d+00,  df7 =  5.595417d+00)
63      parameter (fc1 =  1.923661d+00, fd1  =  2.564881d+00)
64      parameter (crs =  7.876233d-01)
65      real*8 one6th
66      parameter (one6th=1.0d0/6.0d0)
67
68
69      real*8 ETA,ETA2,DNS_CUT
70      parameter (ETA = 0.0d0)
71      parameter (ETA2 = 1.0d-20)
72      parameter (DNS_CUT        =      1.0d-18)
73
74      real*8 tolrho,minagr
75      parameter (tolrho=2.0e-11,minagr=1.0e-12)
76      !parameter (tolrho=2.0e-8,minagr=1.0e-10)
77
78*     ***LYP parameters****************
79      real*8 a,b,c,d
80      parameter (a = 0.04918d0)
81      parameter (b = 0.132d0)
82      parameter (c = 0.2533d0)
83      parameter (d = 0.349d0)
84*     ***Becke88 parameters*************
85      real*8 beta
86      parameter (beta = 0.0042d0)
87
88      real*8 thrd,AA
89      parameter (thrd = 1.0d0/3.0d0)
90
91      integer j
92      real*8 twthrd,frthrd,fvthrd,snthrd,etthrd
93      real*8 pi, Cf,lda_c
94
95      real*8 nup,ndn,n,agrup,agrdn,agr
96      real*8 agr2, agrup2,agrdn2
97      real*8 n2,nup2,ndn2
98      real*8 gamma,gamma_u,gamma_uu
99      real*8 gamma_d,gamma_dd,gamma_du
100      real*8 F,F_u,F_uu
101      real*8 F_d,F_dd,F_du
102      real*8 n13d
103      real*8 enthrd
104      real*8 G,G_u,G_uu
105      real*8 G_d,G_dd,G_du
106      real*8 fc,fclda,H0,Hu,Hd
107      real*8 ce
108      real*8 Hu_u,Hu_d,Hd_d,Hd_u
109      real*8 fc_u,fc_d
110      real*8 H0_u,H0_d
111      real*8 fclda_u,fclda_d
112      real*8 fc_agr, fc_agrup,fc_agrdn
113
114      real*8 chiup,chidn
115      real*8 chiup2,chidn2,chiupSQ,chidnSQ
116      real*8 Kup,Kdn,F1up,F1dn,F2up,F2dn
117      real*8 xeup_pbe,fdnxup_pbe,fdagrxup_pbe
118      real*8 xedn_pbe,fdnxdn_pbe,fdagrxdn_pbe
119      real*8 xeup,xedn,xe
120      real*8 fdnxup,fdnxdn,fdagrxup,fdagrxdn
121      real*8  Q,Q1,P,P1,n_m,n2_m,n3_m,n4_m
122      real*8  n13d_m,n13d2_m,n13d3_m,d3
123      real*8  n_thrd,nupndn
124      real*8  n_mthrd,n_mfrthrd,n_mfvthrd
125      real*8   nup_etthrd,nup_fvthrd
126      real*8   ndn_etthrd,ndn_fvthrd
127      real*8 x,xi,ff,dff,xxp,dxxp,xxf,dxxf
128      real*8 ex_p,ux_p,ex_f,ux_f
129      real*8 ec_p,uc_p,ec_f,uc_f
130      real*8 xe_dirac,fnxup_dirac,fnxdn_dirac
131      real*8 ce_vosko,fncup_vosko,fncdn_vosko
132
133
134      twthrd = thrd*2.0d0
135      frthrd = thrd*4.0d0
136      fvthrd = thrd*5.0d0
137      snthrd = thrd*7.0d0
138      etthrd = thrd*8.0d0
139
140      pi = 4.0d0*datan(1.0d0)
141      Cf = dble(3.0d0*pi*pi)
142      Cf = dble(3.0d0*(Cf**(2.0d0/3.0d0))/10.0d0)
143
144      lda_c =(3.0d0/2.0d0)*(3.0d0/(4.0d0*pi))**(1.0d0/3.0d0)
145      AA = 2.0d0**thrd-1.0d0
146
147
148!$OMP DO
149      do j=1,n2ft3d
150       nup     = dn_in(j,1) + 0.5d0*ETA2
151       ndn     = dn_in(j,2) + 0.5d0*ETA2
152
153       agrup   = agr_in(j,1)
154       agrdn   = agr_in(j,2)
155       agrup2  = agrup*agrup
156       agrdn2  = agrdn*agrdn
157
158       n       = nup + ndn
159       agr     = agr_in(j,3)
160       agr2    = agr*agr
161
162       n2      = n*n
163       nup2    = nup*nup
164       ndn2    = ndn*ndn
165
166*      ***** LSDA terms ****
167       x = crs/n**one6th
168       xi = (nup-ndn)/n
169
170       ff =  ( (1.0d0+xi)**frthrd
171     >       + (1.0d0-xi)**frthrd - 2.0d0)/(2.0d0*AA)
172       dff = frthrd*(  (1.0d0+xi)**thrd
173     >               - (1.0d0-xi)**thrd) /(2.0d0*AA)
174
175*      **** dirac exchange ****
176       ex_p   = (xp/x**2)
177       ux_p   = frthrd*(xp/x**2)
178       ex_f   = (xf/x**2)
179       ux_f   = frthrd*(xf/x**2)
180       xe_dirac    = ex_p + ff*(ex_f-ex_p)
181       fnxup_dirac = ux_p + ff*(ux_f-ux_p) + (+1.0d0-xi)*dff*(ex_f-ex_p)
182       fnxdn_dirac = ux_p + ff*(ux_f-ux_p) + (-1.0d0-xi)*dff*(ex_f-ex_p)
183
184*      **** vosko correlation ****
185       xxp  = x**2 + bp*x + cpt
186       dxxp = 2.0d0*x + bp
187       ec_p = cp1*dlog(x**2/xxp) + cp2*dlog( (x+cp3)*(x+cp3)/xxp)
188     >      + cp4*datan(cp5/(x+cp6))
189       uc_p = ec_p
190     >      - one6th*x*(  dp1/x + dp2/(x+cp3) + dp4*dxxp/xxp
191     >                    + dp5/( (x+dp6)*(x+dp6)+dp7))
192
193       xxf  = x**2 + bf*x + cft
194       dxxf = 2.0d0*x + bf
195       ec_f = cf1*dlog(x**2/xxf) + cf2*dlog( (x+cf3)*(x+cf3)/xxf)
196     >     + cf4*datan(cf5/(x+cf6))
197       uc_f = ec_f
198     >      - one6th*x*(  df1/x + df2/(x+cf3) + df4*dxxf/xxf
199     >                    + df5/( (x+df6)*(x+df6)+df7))
200       ce_vosko    = ec_p + ff*(ec_f - ec_p)
201       fncup_vosko = uc_p + ff*(uc_f-uc_p) + (+1.0d0-xi)*dff*(ec_f-ec_p)
202       fncdn_vosko = uc_p + ff*(uc_f-uc_p) + (-1.0d0-xi)*dff*(ec_f-ec_p)
203
204       if ((dn_in(j,1)+dn_in(j,2)).lt.DNS_CUT) then
205          xe       = 0.0d0
206          fdnxup   = 0.0d0
207          fdnxdn   = 0.0d0
208          fdagrxup = 0.0d0
209          fdagrxdn = 0.0d0
210       else
211
212*      *******exchange part***************
213
214
215*      **************UP*******************
216          if (dn_in(j,1).lt.DNS_CUT) then
217             xeup     = 0.0d0
218             fdnxup   = 0.0d0
219             fdagrxup = 0.0d0
220          else
221             chiup = agrup/nup**(4.0d0/3.0d0)
222             chiup2 = chiup*chiup
223             chiupSQ = dsqrt(1.0d0+chiup2)
224
225             Kup = 6.0d0*beta*dlog(chiup+chiupSQ)
226             F1up = chiup2/(1.0d0 + chiup*Kup)
227             xeup = -nup**thrd*(lda_c + beta*F1up)
228             F2up = (2.0d0 + chiup*Kup
229     &              - 6.0d0*beta*chiup2/chiupSQ)
230     &              /(1.0d0+chiup*Kup)**2.0d0
231             fdnxup = -nup**(thrd)*(4.0d0/3.0d0)
232     &             *(lda_c+beta*(F1up-chiup2*F2up))
233             fdagrxup = -beta*chiup*F2up
234             if ((fdnxup-xeup).gt.0.0d0) then
235               call gen_PBE96_x_unrestricted(nup,agrup,
236     >                                xeup_pbe,fdnxup_pbe,fdagrxup_pbe)
237
238               call Becke_smalln_correction(nup,nup**thrd,1.0d0,
239     >                                      beta,lda_c,chiup,chiup2,
240     >                                      chiupSQ,Kup,F1up,F2up,
241     >                                xeup_pbe,fdnxup_pbe,fdagrxup_pbe,
242     >                                xeup,fdnxup,fdagrxup)
243
244*               call Becke_smalln_correction(nup,nup**thrd,agrup,
245*    >                               beta,lda_c,
246*    >                               chiup,chiup2,chiupSQ,Kup,F1up,F2up,
247*    >                               xeup,fdnxup,fdagrxup)
248             end if
249          end if
250*      ************END UP*****************
251
252
253*      *************DOWN******************
254          if (dn_in(j,2).lt.DNS_CUT) then
255             xedn     = 0.0d0
256             fdnxdn   = 0.0d0
257             fdagrxdn = 0.0d0
258          else
259             chidn = agrdn/ndn**(4.0d0/3.0d0)
260             chidn2 = chidn*chidn
261             chidnSQ = dsqrt(1.0d0+chidn2)
262
263             Kdn = 6.0d0*beta*dlog(chidn+chidnSQ)
264             F1dn = chidn2/(1.0d0 + chidn*Kdn)
265             xedn = -ndn**thrd*(lda_c + beta*F1dn)
266             F2dn = (2.0d0 + chidn*Kdn
267     &              - 6.0d0*beta*chidn2/chidnSQ)
268     &              /(1.0d0+chidn*Kdn)**2.0d0
269             fdnxdn = -ndn**(thrd)*(4.0d0/3.0d0)
270     &                *(lda_c+beta*(F1dn-chidn2*F2dn))
271             fdagrxdn = -beta*chidn*F2dn
272             if ((fdnxdn-xedn).gt.0.0d0) then
273               call gen_PBE96_x_unrestricted(ndn,agrdn,
274     >                                xedn_pbe,fdnxdn_pbe,fdagrxdn_pbe)
275
276               call Becke_smalln_correction(ndn,ndn**thrd,1.0d0,
277     >                                      beta,lda_c,chidn,chidn2,
278     >                                      chidnSQ,Kdn,F1dn,F2dn,
279     >                                xedn_pbe,fdnxdn_pbe,fdagrxdn_pbe,
280     >                                xedn,fdnxdn,fdagrxdn)
281
282*               call Becke_smalln_correction(ndn,ndn**thrd,agrdn,
283*    >                               beta,lda_c,
284*    >                               chidn,chidn2,chidnSQ,Kdn,F1dn,F2dn,
285*    >                               xedn,fdnxdn,fdagrxdn)
286             end if
287          end if
288*      ***********END DOWN****************
289
290          xe = (xeup*nup + xedn*ndn)/n
291
292*      *******end excange part************
293       end if
294
295
296*      *******correlation part************
297       n_thrd   = n**thrd
298       n_m     = 1.0d0/n
299       n2_m    = n_m*n_m
300       n3_m    = n2_m*n_m
301       n4_m    = n3_m*n_m
302       nupndn  = nup*ndn
303       n_mthrd = 1.0d0/n_thrd
304       n_mfrthrd = n_mthrd*n_m
305       n_mfvthrd = n_mfrthrd*n_mthrd
306
307       nup_etthrd = nup**etthrd
308       nup_fvthrd = nup**fvthrd
309       ndn_etthrd = ndn**etthrd
310       ndn_fvthrd = ndn**fvthrd
311
312
313
314
315       gamma   = (4.0d0*nupndn)*n2_m
316       gamma_u = 4.0d0*ndn*n2_m - 8.0d0*nupndn*n3_m
317       gamma_d = 4.0d0*nup*n2_m - 8.0d0*nupndn*n3_m
318
319       gamma_uu = -16.0d0*ndn*n3_m + 24.0d0*nupndn*n4_m
320       gamma_dd = -16.0d0*nup*n3_m + 24.0d0*nupndn*n4_m
321
322       gamma_du =  (6.0d0*gamma-4.0d0)*n2_m
323
324
325
326       d3         = d/3.0d0
327       n13d_m     = 1.0d0/(1.0d0 + d*n_mthrd)
328       n13d2_m    = n13d_m*n13d_m
329       n13d3_m    = n13d2_m*n13d_m
330       F        = gamma*n13d_m
331       F_u      = gamma_u*n13d_m + d3*gamma*n_mfrthrd*n13d2_m
332       F_d      = gamma_d*n13d_m + d3*gamma*n_mfrthrd*n13d2_m
333
334       F_uu   = gamma_uu*n13d_m + d3*(gamma_u+gamma_u)*n_mfrthrd*n13d2_m
335     &        - (4.0d0/9.0d0)*d*gamma*n_mfrthrd*n_m*n13d2_m
336     &        + (2.0d0/9.0d0)*d*d*gamma*n_mfrthrd*n_mfrthrd*n13d3_m
337
338       F_dd   = gamma_dd*n13d_m + d3*(gamma_d+gamma_d)*n_mfrthrd*n13d2_m
339     &        - (4.0d0/9.0d0)*d*gamma*n_mfrthrd*n_m*n13d2_m
340     &        + (2.0d0/9.0d0)*d*d*gamma*n_mfrthrd*n_mfrthrd*n13d3_m
341
342       F_du   = gamma_du*n13d_m + d3*(gamma_u+gamma_d)*n_mfrthrd*n13d2_m
343     &        - (4.0d0/9.0d0)*d*gamma*n_mfrthrd*n_m*n13d2_m
344     &        + (2.0d0/9.0d0)*d*d*gamma*n_mfrthrd*n_mfrthrd*n13d3_m
345
346
347       enthrd   = dexp(-c*n_mthrd)
348       Q        = enthrd*n_mfvthrd
349c       Q1       = (1.0d0/3.0d0)*c*n_mfrthrd  *enthrd
350c     &          - (5.0d0/3.0d0)*n_mfvthrd*n_m*enthrd
351       Q1 = (Q/3.0d0)*(c*n_mfrthrd - 5.0*n_m)
352       G        = F*Q
353
354       P  = (1.0d0/3.0d0)*c*n_mfrthrd - (5.0d0/3.0d0)*n_m
355       P1 = ((-4.0d0/9.0d0)*c*n_mfrthrd + (5.d0/3.0d0)*n_m)*n_m
356
357       G_u      = F_u*Q + G*P
358       G_d      = F_d*Q + G*P
359
360       G_uu     = F_uu*Q + F_u*Q1 + G_u*P + G*P1
361       G_dd     = F_dd*Q + F_d*Q1 + G_d*P + G*P1
362       G_du     = F_du*Q + F_d*Q1 + G_u*P + G*P1
363
364
365
366       fclda = -a*F*n - 2.0d0*a*b*G*Cf*(2.0d0**twthrd)
367     &         *(nup_etthrd + ndn_etthrd)
368
369       fclda_u = -a*F_u*n - a*F
370     &           - 2.0d0*a*b*G_u*Cf*(2.0d0**twthrd)
371     &           *(nup_etthrd + ndn_etthrd)
372     &           - 2.0d0*a*b*G*Cf*(2.0d0**twthrd)
373     &           *(8.0d0/3.0d0)*nup_fvthrd
374
375       fclda_d = -a*F_d*n - a*F
376     &           - 2.0d0*a*b*G_d*Cf*(2.0d0**twthrd)
377     &           *(nup_etthrd + ndn_etthrd)
378     &           - 2.0d0*a*b*G*Cf*(2.0d0**twthrd)
379     &           *(8.0d0/3.0d0)*ndn_fvthrd
380
381
382
383       H0 = (a*b/2.0d0)*(G
384     &    + (1.0d0/3.0d0)*(nup*G_d + ndn*G_u)
385     &    + (1.0d0/4.0d0)*(nup*G_u + ndn*G_d))
386       H0_u=(a*b/2.0d0)*(G_u
387     &     + (1.0d0/3.0d0)*(G_d + nup*G_du + ndn*G_uu)
388     &     + (1.0d0/4.0d0)*(G_u + nup*G_uu + ndn*G_du))
389       H0_d=(a*b/2.0d0)*(G_d
390     &     + (1.0d0/3.0d0)*(G_u + ndn*G_du + nup*G_dd)
391     &     + (1.0d0/4.0d0)*(G_d + ndn*G_dd + nup*G_du))
392
393
394       Hu = (a*b/18.0d0)*(G + (15.0d0/4.0d0)*nup*G_u
395     &                      -  (9.0d0/4.0d0)*ndn*G_d
396     &                      -        (3.0d0)*nup*G_d
397     &                      +  (3.0d0/2.0d0)*ndn*G_u)
398
399       Hu_u = (a*b/18.0d0)*(G_u + (15.0d0/4.0d0)*(G_u + nup*G_uu)
400     &                          -  (9.0d0/4.0d0)*(ndn*G_du)
401     &                          -        (3.0d0)*(G_d + nup*G_du)
402     &                          +  (3.0d0/2.0d0)*(ndn*G_uu))
403
404       Hu_d = (a*b/18.0d0)*(G_d + (15.0d0/4.0d0)*(nup*G_du)
405     &                          -  (9.0d0/4.0d0)*(G_d+ndn*G_dd)
406     &                          -        (3.0d0)*(nup*G_dd)
407     &                          +  (3.0d0/2.0d0)*(G_u + ndn*G_du))
408
409
410
411       Hd = (a*b/18.0d0)*(G + (15.0d0/4.0d0)*ndn*G_d
412     &                      -  (9.0d0/4.0d0)*nup*G_u
413     &                      -        (3.0d0)*ndn*G_u
414     &                      +  (3.0d0/2.0d0)*nup*G_d)
415
416       Hd_d = (a*b/18.0d0)*(G_d + (15.0d0/4.0d0)*(G_d + ndn*G_dd)
417     &                          -  (9.0d0/4.0d0)*(nup*G_du)
418     &                          -        (3.0d0)*(G_u + ndn*G_du)
419     &                          +  (3.0d0/2.0d0)*(nup*G_dd))
420
421       Hd_u = (a*b/18.0d0)*(G_u + (15.0d0/4.0d0)*(ndn*G_du)
422     &                          -  (9.0d0/4.0d0)*(G_u+nup*G_uu)
423     &                          -        (3.0d0)*(ndn*G_uu)
424     &                          +  (3.0d0/2.0d0)*(G_d + nup*G_du))
425
426
427       fc = fclda + H0*agr2 + Hu*agrup2 + Hd*agrdn2
428
429*    ***calculate derivatives w.r.t up and down density
430       fc_u = fclda_u + H0_u*agr2 + Hu_u*agrup2 + Hd_u*agrdn2
431       fc_d = fclda_d + H0_d*agr2 + Hu_d*agrup2 + Hd_d*agrdn2
432
433*    ***calculate derivatives w.r.t. up,down and total density gradients
434       fc_agr   = 2.0d0*H0*agr
435       fc_agrup = 2.0d0*Hu*agrup
436       fc_agrdn = 2.0d0*Hd*agrdn
437
438       ce = fc/n
439
440
441*      *******end correlation part********
442
443c       write(*,*) "n,nup,ndn,agrup,agrdn,:",j,n,nup,ndn,agrup,agrdn,agr
444c       write(*,*) "xe,ce         :",j,xe,ce
445c       write(*,*) "fdnxup,fdncup     :",j,fdnxup,fc_u
446c       write(*,*) "fdnxdn,fdncdn     :",j,fdnxdn,fc_d
447c       write(*,*) "fdagrxup,fdargrcup:",j,fdagrxup,
448c     >                                fc_agrup,fc_agr
449c       write(*,*) "fdagrxdn,fdargrcdn:",j,fdagrxdn,
450c     >                                fc_agrdn,fc_agr
451c
452c       write(*,*) "restricted fdagrc",fc_agr+0.5d0*(fc_agrdn+fc_agrup)
453c       write(*,*) "fc_lda,fdnc_lda:",fclda,fclda_u,fclda_d
454c       write(*,*) "Hu,Hd,H0:",Hu,Hd,H0
455c       write(*,*) "Ho:",H0+0.25d0*Hu+0.25d0*Hd
456c       write(*,*) "Ho_n:",0.5d0*(H0_u+H0_d+0.25d0*(Hu_u+Hu_d+Hd_d+Hd_u))
457c       write(*,*) "F,G:",F,G
458c       write(*,*) "F_u,G_u:",F_u,G_u
459c       write(*,*) "F_d,G_d:",F_d,G_d
460c       write(*,*) "G_uu,G_dd,G_du:",G_uu,G_dd,G_du
461c       write(*,*) "F_uu,F_dd,F_du:",F_uu,F_dd,F_du
462c       write(*,*) "Gamma's:",gamma,gamma_u,gamma_d,
463c     >             gamma_uu,gamma_dd,gamma_du
464c       write(*,*) "n13d",n13d
465c       write(*,*) "0.5/nup2,0.5/ndn2:",0.5d0/nup2,0.5d0/ndn2
466c       write(*,*) "fc:",fc
467c       write(*,*) "0.5d0*(F_uu+F_du)",0.5d0*(F_uu + F_du)
468c       write(*,*) "0.5d0*(G_uu+G_du)",0.5d0*(G_uu + G_du)
469c       write(*,*) "nup*G_u...",nup*G_u,ndn*G_d,nup*G_d,ndn*G_u
470c       write(*,*)
471
472
473**** a*HF + (1-a)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko
474**** x_parameter = (1-a)
475
476*      ***return blyp exchange correlation values***
477       xce(j)  = x_parameter*0.9d0*xe     + c_parameter*0.81d0*ce
478       fn(j,1) = x_parameter*0.9d0*fdnxup + c_parameter*0.81d0*fc_u
479       fn(j,2) = x_parameter*0.9d0*fdnxdn + c_parameter*0.81d0*fc_d
480
481       fdn(j,1) = x_parameter*0.9d0*fdagrxup + c_parameter*0.81*fc_agrup
482       fdn(j,2) = x_parameter*0.9d0*fdagrxdn + c_parameter*0.81*fc_agrdn
483       fdn(j,3) =                              c_parameter*0.81*fc_agr
484
485       xce(j) = xce(j) + x_parameter*0.10d0*xe_dirac
486     >                 + c_parameter*0.19d0*ce_vosko
487       fn(j,1)  = fn(j,1)  + x_parameter*0.10d0*fnxup_dirac
488     >                     + c_parameter*0.19d0*fncup_vosko
489       fn(j,2)  = fn(j,2)  + x_parameter*0.10d0*fnxdn_dirac
490     >                     + c_parameter*0.19d0*fncdn_vosko
491
492      end do
493!$OMP END DO
494      return
495      end
496
497*     **************************************************
498*     *                                                *
499*     *            gen_B3LYP_BW_restricted             *
500*     *                                                *
501*     **************************************************
502
503*     blyp restricted  calc.
504*
505*
506*
507
508*      subroutine gen_B3LYP_BW_restricted(n2ft3d,rho_in,agr_in,xce,xcp,fn,fdn)
509*      input:  n2ft3d                  grid
510*              rho_in                  density
511*              agr_in                  absolute gradient of density
512*              x_parameter:            scale parameter for exchange
513*              c_parameter:            scale parameter for correlation
514*      output: xce                     exchange correlation energy density
515*              fn                      d(n*exc)/dn
516*              fdn                     d(n*exc)/d(|grad n|)
517
518      subroutine gen_B3LYP_BW_restricted(n2ft3d,rho_in,agr_in,
519     &                               x_parameter,
520     &                               c_parameter,xce,fn,fdn)
521      implicit none
522
523      integer   n2ft3d
524      real*8    rho_in(n2ft3d)
525      real*8    agr_in(n2ft3d)
526      real*8    xce(n2ft3d)
527      real*8    fn(n2ft3d)
528      real*8    fdn(n2ft3d)
529      real*8    x_parameter, c_parameter
530
531c*---- parameters given by vosko et al -----------------*
532      real*8 ap,af,x0p,x0f,bp,bf,cpt,cft
533      parameter (ap = 3.109070d-02, af = 1.554530d-02)
534      parameter (x0p=-1.049800d-01, x0f=-3.250000d-01)
535      parameter (bp = 3.727440d+00, bf = 7.060420d+00)
536      parameter (cpt = 1.293520d+01, cft = 1.805780d+01)
537c*------------------------------------------------------*
538
539*     constants calculated from vosko's parameters
540      real*8 xp,xf,qp,qf,xx0p,xx0f,fc1,fd1,crs
541      real*8 cp1,cp2,cp3,cp4,cp5,cp6,dp1,dp2,dp3,dp4,dp5,dp6,dp7
542      real*8 cf1,cf2,cf3,cf4,cf5,cf6,df1,df2,df3,df4,df5,df6,df7
543      parameter (xp  = -4.581653d-01,  xf  = -5.772521d-01)
544      parameter (qp  =  6.151991d+00,  qf  =  4.730927d+00)
545      parameter (xx0p=  1.255491d+01,  xx0f=  1.586879d+01)
546      parameter (cp1 =  3.109070d-02,  cf1 =  1.554530d-02)
547      parameter (cp2 =  9.690228d-04,  cf2 =  2.247860d-03)
548      parameter (cp3 =  1.049800d-01,  cf3 =  3.250000d-01)
549      parameter (cp4 =  3.878329d-02,  cf4 =  5.249122d-02)
550      parameter (cp5 =  3.075995d+00,  cf5 =  2.365463d+00)
551      parameter (cp6 =  1.863720d+00,  cf6 =  3.530210d+00)
552      parameter (dp1 =  6.218140d-02,  df1 =  3.109060d-02)
553      parameter (dp2 =  1.938045d-03,  df2 =  4.495720d-03)
554      parameter (dp3 =  1.049800d-01,  df3 =  3.250000d-01)
555      parameter (dp4 = -3.205972d-02,  df4 = -1.779316d-02)
556      parameter (dp5 = -1.192972d-01,  df5 = -1.241661d-01)
557      parameter (dp6 =  1.863720d+00,  df6 =  3.530210d+00)
558      parameter (dp7 =  9.461748d+00,  df7 =  5.595417d+00)
559      parameter (fc1 =  1.923661d+00, fd1  =  2.564881d+00)
560      parameter (crs =  7.876233d-01)
561      real*8 for3rd,one6th
562      parameter (for3rd=4.0d0/3.0d0,one6th=1.0d0/6.0d0)
563
564
565
566*******local declarations***************************************
567      integer   i
568      real*8    Fc,Gc,C1,C2
569      real*8    n, n_thrd,n_fv,n_fr,n_tw
570      real*8    n_m,n_mthrd,n_mfrthrd,n_mfvthrd
571      real*8    agr,agr2, chi, chi2,chiSQ,sd
572      real*8    K
573      real*8    F1, F2
574      real*8    xe_pbe,fdnx_pbe,fdagrx_pbe
575      real*8    xe,fdnx,fdagrx,xe_dirac,fnx_dirac,x,xx
576      real*8    ce,fdnc, fdagrc,fdnc_lda,ce_vosko,fnc_vosko
577
578      real*8    fc_lda,Ho,Ho_n,Gc_n,Gc_nn,Fc_n,Fc_nn
579      real*8    P,P_n
580
581******* constants **********************************
582      real*8 pi,thrd,two_thrd
583      parameter (pi=3.14159265358979311599d0)
584      parameter (thrd=1.0d0/3.0d0)
585      parameter (two_thrd=1.25992104989487319066d0)  ! two_thrd = 2**thrd
586
587
588*******density cutoff parameters********************
589      real*8 DNS_CUT, ETA
590      parameter (DNS_CUT        =      1.0d-18)
591      parameter (ETA            =      1.0d-20)
592
593      real*8 tolrho,minagr
594      parameter (tolrho=2.0e-11,minagr=1.0e-12)
595****** Becke constants *****************************
596      real*8 lda_c,beta
597      parameter (beta = 0.0042d0)
598      parameter (lda_c = 0.93052573634910018540d0)   ! lda_c = (3/2)*(3/(4*pi))**(1/3)
599*******LYP correlation parameters a, b, c, d********
600      real*8 a,b,c,d,Cf
601      parameter (a              =      0.04918d0)
602      parameter (b              =      0.132d0)
603      parameter (c              =      0.2533d0)
604      parameter (d              =      0.349d0)
605      parameter (Cf =2.87123400018819108225d0)     ! Cf = (3/10)*(3*pi*pi)**(2/3)
606
607****** collated LYP parameters *********************
608      real*8 ho1,ho2,ho3,thrd_d,abCf,p2,p3,p4,p5,p6
609      parameter (ho1=(19.0d0/36.0d0)*a*b)
610      parameter (ho2=( 7.0d0/24.0d0)*a*b)
611      parameter (ho3=(59.0d0/72.0d0)*a*b)
612      parameter (thrd_d=thrd*d)
613      parameter (abCf=a*b*Cf)
614      parameter (p2=-5.0d0*thrd,p3=c*thrd)
615      parameter (p4=-4.0d0*thrd*thrd_d)
616      parameter (p5=5.0d0*thrd)
617      parameter (p6=-4.0d0*thrd*thrd*c)
618
619****************************************************
620
621
622!$OMP DO
623      do i=1,n2ft3d
624       n        = rho_in(i) + ETA
625       agr      = agr_in(i)
626       n_thrd 	= n**thrd
627
628*      ***** LDA terms ****
629       x = crs/n**one6th
630
631       !**** dirac exchange ****
632       xe_dirac   = (xp/x**2)
633       fnx_dirac = for3rd*(xp/x**2)
634
635       !**** vosko correlation ****
636       xx=1.0d0/(x*(x+bp)+cpt)
637       ce_vosko=cp1*dlog(xx*x**2)+cp2*dlog(xx*(x+cp3)**2)
638     >       +cp4*datan(cp5/(x+cp6))
639       fnc_vosko = ce_vosko-one6th*x*(
640     >           dp1/x+dp2/(x+dp3)+dp4*xx*(2.0d0*x+bp)
641     >          +dp5/((x+dp6)**2+dp7) )
642
643
644       if (rho_in(i).lt.DNS_CUT) then
645         xe     = 0.0d0
646         fdnx   = 0.0d0
647         fdagrx = 0.0d0
648       else
649
650******************************************************************
651*     *******calc. becke exchange energy density, fnx, fdnx*******
652*****************************************************************
653          sd       = 1.0d0/(n_thrd*n)
654          chi      = two_thrd*agr*sd
655          chi2     = chi*chi
656          chiSQ    = dsqrt(1.0d0+chi2)
657
658          K        = 6.0d0*beta*dlog(chi+chiSQ)
659          F1       = chi2/(1.0d0+chi*K)
660          xe       = -n_thrd*(lda_c+beta*F1)/two_thrd
661          F2       = (2.0d0 + chi*K-(chi2)*6.0d0*beta
662     &               /chiSQ)
663     &             /((1.0d0+chi*K)*(1.0d0+chi*K))
664          fdnx     = -(n_thrd/two_thrd)*dble(4.0d0/3.0d0)
665     &            *(lda_c+beta*(F1-chi2*F2))
666          fdagrx   = -beta*chi*F2
667          if ((fdnx-xe).gt.0.0d0) then
668              call gen_PBE96_x_restricted(n,agr,
669     >                                    xe_pbe,fdnx_pbe,fdagrx_pbe)
670              call Becke_smalln_correction(n,n**thrd,2.0d0,beta,
671     >                                   lda_c,chi,chi2,chiSQ,K,F1,F2,
672     >                                   xe_pbe,fdnx_pbe,fdagrx_pbe,
673     >                                   xe,fdnx,fdagrx)
674
675*             call Becke_smalln_correction(n,n_thrd,agr,beta,lda_c,
676*    >                                   chi,chi2,chiSQ,K,F1,F2,
677*    >                                   xe,fdnx,fdagrx)
678          end if
679       end if
680
681
682*******final result for restricted LYP*****************************
683       agr2 = agr*agr
684
685       n_m     = 1.0d0/n
686       n_mthrd = 1.0d0/n_thrd
687       n_mfrthrd = n_mthrd*n_m
688       n_mfvthrd = n_mfrthrd*n_mthrd
689
690       n_fv = n_thrd*n_thrd*n_thrd*n_thrd*n_thrd
691       n_fr = n_thrd*n_thrd*n_thrd*n_thrd
692       n_tw = n_thrd*n_thrd
693
694
695       Fc = (1.0d0/(1.0d0+d*n_mthrd))
696       Fc_n = thrd_d*n_mfrthrd*Fc*Fc
697       Fc_nn = thrd_d*(-4.0d0*thrd)*n_mfrthrd*n_m*Fc*Fc
698     >       + thrd_d*n_mfrthrd*2.0d0*Fc*Fc_n
699
700       Gc = Fc*dexp(-c*n_mthrd)*n_mfvthrd
701
702       P  = (thrd_d*Fc*n_mfrthrd - 5.0d0*thrd*n_m + c*thrd*n_mfrthrd)
703
704       P_n = thrd_d*Fc_n*n_mfrthrd
705     >     - thrd_d*Fc  *n_mfrthrd*(4.0d0*thrd*n_m)
706     >     + 5.0d0*thrd*n_m*n_m
707     >     - 4.0d0*thrd*thrd*c*n_mfrthrd*n_m
708
709c       P  = (thrd_d*Fc*n_mfrthrd + p2*n_m + p3*n_mfrthrd)
710c       P_n = thrd_d*Fc_n*n_mfrthrd
711c     >     + p4*Fc*n_mfrthrd*n_m
712c     >     + p5*n_m*n_m
713c     >     + p6*n_mfrthrd*n_m
714
715       Gc_n  = Gc*P
716       Gc_nn = Gc*P*P + Gc*P_n
717
718
719       fc_lda   = -a*Fc*n
720     >          - abCf*Gc*(n_fr*n_fr)
721       fdnc_lda = -a*Fc_n*n
722     >          - a*Fc
723     >          - abCf*Gc_n*n_fr*n_fr
724     >          - 8.0d0*thrd*abCf*Gc*n_fv
725
726c      Ho   = (19.0d0/36.0d0)*(a*b*Gc)   + (7.0d0/24.0d0)*a*b*Gc_n*n
727c      Ho_n = (59.0d0/72.0d0)*(a*b*Gc_n) + (7.0d0/24.0d0)*a*b*Gc_nn*n
728       Ho   = ho1*Gc   + ho2*Gc_n*n
729       Ho_n = ho3*Gc_n + ho2*Gc_nn*n
730
731       ce = (fc_lda + Ho*agr2)*n_m
732       fdnc = fdnc_lda + Ho_n*agr2
733
734       fdagrc = 2.0d0*Ho*agr
735
736c       write(*,*) "n,agr         :",i,n,agr
737c       write(*,*) "xe,ce         :",i,xe,ce
738c       write(*,*) "fdnx,fdnc     :",i,fdnx,fdnc
739c       write(*,*) "fdagrx,fdargrc:",i,fdagrx,fdagrc
740c       write(*,*) "fc_lda,fdnc_lda:",fc_lda,fdnc_lda
741c       write(*,*) "Ho :",Ho
742c       write(*,*) "Ho_n:",Ho_n
743c       write(*,*) "F,G:",Fc,Gc
744c       write(*,*) "F_n,G_n:",Fc_n,Gc_n
745c       write(*,*) "G_nn:",Gc_nn
746c       write(*,*) "F_nn:",Fc_nn
747c       write(*,*) "n13d:",1.0d0+d*n_mthrd
748c       write(*,*) "fc:",fc_lda + Ho*agr2
749c       write(*,*)
750
751**** a*HF + (1-a)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko
752**** x_parameter = (1-a)
753
754       xce(i) = x_parameter*0.9d0*xe     + c_parameter*0.81d0*ce
755       fn(i)  = x_parameter*0.9d0*fdnx   + c_parameter*0.81d0*fdnc
756       fdn(i) = x_parameter*0.9d0*fdagrx + c_parameter*0.81d0*fdagrc
757
758       xce(i) = xce(i) + x_parameter*0.10d0*xe_dirac
759     >                 + c_parameter*0.19d0*ce_vosko
760       fn(i)  = fn(i)  + x_parameter*0.10d0*fnx_dirac
761     >                 + c_parameter*0.19d0*fnc_vosko
762      end do
763!$OMP END DO
764      return
765      end
766