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