1c $Id$
2
3*    ************************************
4*    *                                  *
5*    *          nwpw_gga                *
6*    *                                  *
7*    ************************************
8      subroutine nwpw_gga(gga,ic,ispin,
9     >                   rho,agr,
10     >                   xce,fn,fdn,tmp)
11      implicit none
12      integer gga,ic,ispin
13      double precision rho(ic,ispin)
14      double precision agr(ic,*)
15
16      double precision xce(ic)
17      double precision fn(ic,ispin)
18      double precision fdn(ic,*)
19      double precision tmp(ic)
20
21#include "nwpwxc.fh"
22#include "util.fh"
23
24*     **** local variables ****
25      double precision grad(ispin+1)
26      double precision dtmp
27      double precision dncut
28      parameter(dncut = 1.0d-30)
29
30      integer i
31      logical use_nwpwxc
32
33
34      use_nwpwxc = nwpwxc_is_on()
35
36c     **** restricted calculations ****
37      if (ispin.eq.1) then
38
39!$OMP DO
40         do i=1,ic
41           tmp(i) = (rho(i,1)+rho(i,1))
42         end do
43!$OMP END DO
44
45cccccc USE_NWXC needs to be checked cccccccc
46         if (use_nwpwxc) then
47           !call D3db_rr_Sqr1(1,agr)
48!OMP DO
49           do i=1,ic
50              agr(i,1) = agr(i,1)**2
51           end do
52!OMP END DO
53!$OMP MASTER
54           call nwpwxc_eval_df(1,ic,tmp,agr,
55     >                       dtmp,xce,
56     >                       fn,fdn,dtmp)
57!$OMP END MASTER
58!$OMP BARRIER
59c
60c          Combine (df/d|grad a|) with (df/d(grad a|grad b))
61c
62           !call D3dB_r_SMul(1,0.5d0,fdn(1,2),fdn)
63c
64c          Calculate energy density from energy
65c
66!$OMP DO
67           do i = 1, ic
68              fdn(i,1) = fdn(i,1) + 0.5d0*fdn(i,2)
69              xce(i) = xce(i)/(tmp(i)+dncut)
70              agr(i,1) = dsqrt(agr(i,1))
71              fdn(i,1) = fdn(i,1)*agr(i,1)
72           enddo
73!$OMP END DO
74           !call D3dB_rr_Sqrt1(1,agr)
75           !call D3dB_rr_Mul2(1,agr,fdn)
76cccccc USE_NWXC needs to be checked cccccccc
77
78         else
79         if (gga.eq.10) then
80         call gen_PBE96_BW_restricted(ic,
81     >                                tmp,
82     >                                agr,
83     >                                1.0d0,1.0d0,
84     >                                xce,
85     >                                fn,
86     >                                fdn)
87         else if (gga.eq.11) then
88         call gen_BLYP_BW_restricted(ic,
89     >                                tmp,
90     >                                agr,
91     >                                1.0d0,1.0d0,
92     >                                xce,
93     >                                fn,
94     >                                fdn)
95
96         else if (gga.eq.12) then
97         call gen_revPBE_BW_restricted(ic,
98     >                                tmp,
99     >                                agr,
100     >                                1.0d0,1.0d0,
101     >                                xce,
102     >                                fn,
103     >                                fdn)
104         else
105          call errquit('bad gga',0,0)
106         end if
107         endif ! nwpwxc_is_on
108
109
110c     **** unrestricted calculations ****
111      else if (ispin.eq.2) then
112
113         if (use_nwpwxc) then
114!$OMP DO
115           do i = 1, ic
116             grad(1) = agr(i,1)**2
117             grad(2) = agr(i,3)**2
118             grad(3) = agr(i,2)**2
119             grad(2) = 0.5d0*(grad(2) - grad(1) - grad(3))
120             agr(i,1) = grad(1)
121             agr(i,2) = grad(2)
122             agr(i,3) = grad(3)
123           enddo
124!$OMP END DO
125!$OMP MASTER
126           call nwpwxc_eval_df(2,ic,rho,agr,
127     >                       dtmp,xce,
128     >                       fn,fdn,dtmp)
129!$OMP END MASTER
130!$OMP BARRIER
131!$OMP DO
132           do i = 1, ic
133             xce(i)   = xce(i)/(rho(i,1)+rho(i,2)+dncut)
134             dtmp     = fdn(i,2)
135             fdn(i,2) = fdn(i,3)
136             fdn(i,3) = dtmp
137             grad(1)  = dsqrt(agr(i,1))
138             grad(2)  = dsqrt(agr(i,3))
139             grad(3)  = dsqrt(2.0d0*agr(i,2)+agr(i,1)+agr(i,3))
140             agr(i,1) = grad(1)
141             agr(i,2) = grad(2)
142             agr(i,3) = grad(3)
143             fdn(i,1) = (2.0d0*fdn(i,1)-fdn(i,3))*agr(i,1)
144             fdn(i,2) = (2.0d0*fdn(i,2)-fdn(i,3))*agr(i,2)
145             fdn(i,3) = fdn(i,3)*agr(i,3)
146           enddo
147!$OMP END DO
148         else
149         if (gga.eq.10) then
150         call gen_PBE96_BW_unrestricted(ic,
151     >                                rho,
152     >                                agr,
153     >                                1.0d0,1.0d0,
154     >                                xce,
155     >                                fn,
156     >                                fdn)
157         else if (gga.eq.11) then
158         call gen_BLYP_BW_unrestricted(ic,
159     >                                rho,
160     >                                agr,
161     >                                1.0d0,1.0d0,
162     >                                xce,
163     >                                fn,
164     >                                fdn)
165
166         else if (gga.eq.12) then
167         call gen_revPBE_BW_unrestricted(ic,
168     >                                rho,
169     >                                agr,
170     >                                1.0d0,1.0d0,
171     >                                xce,
172     >                                fn,
173     >                                fdn)
174         else
175          call errquit('bad gga',0,0)
176         end if
177         end if ! nwpwxc_is_on
178
179      else
180        call errquit('bad spin',0,ispin)
181      end if
182      return
183      end
184