1 Subroutine dftg_gridssw(grid_written, 2 , d_qwght,qwght, qxyz, xyz, Rij, rq, p, 3 , dzeta, d_p, 4 . ictr, nctrs_pruned, nq_orig,nq, 5 , lscreen) 6c 7C$Id$ 8c 9 implicit none 10#include "errquit.fh" 11 logical grid_written 12 integer nctrs_pruned ! [in] natoms after signf 13 integer nq ! [in] no. grid pts 14 integer nq_orig ! [in/out] no. grid pts after compression 15 logical lscreen ! [in] screen weights 16c 17 double precision qxyz(3,nq)! grid points [input] 18 double precision xyz(3,*) ! atom coordinates [input] 19c 20 double precision Rij(*) !interatomic distances [input] 21 integer ictr ! [in] ctr of quadr 22 double precision p(nctrs_pruned) 23 double precision rq(nq_orig,*) ! sum of p(n) [ output] 24 double precision qwght(nq_orig) ! weights [output] 25 double precision d_qwght(3,nq_orig,*) ! weight deriv [output] 26c 27 double precision d_p(3,nctrs_pruned), 28 & dzeta(3,nctrs_pruned) 29 integer i, j, ij 30 integer iind,jind, AA,NN 31 double precision mu, radwgh 32 double precision x, y, z, x0i, y0i, z0i 33 double precision sk 34 double precision toll,rag 35 parameter(toll=1d-13) 36 37 logical inotA,jnota 38c 39 integer n 40 double precision ass,distnc,distnn,wsum,dist2 41 double precision xa,ya,za 42 integer A,B,n1,n2 43 double precision asse 44#include "xc_erftab.fh" 45c 46c RE Stratmann, GE Scuseria, MJ Frisch, Chem Phys Lett 257, 213 (1996) 47c Evaluate Stratman space partitioning weight. Then, incorporate it 48c with weights from the single-center quadratures to form the total 49c multi-center quadrature weight. 50c The following 2 lines are to satisfy compiler warnings. 51c 52 NN = 1 53c call dfill(nq, 0.d0, zeta, 1) 54 call dfill(nq*3*nctrs_pruned, 0.d0, d_qwght, 1) 55 ass=ass_erf1 56 asse=ass-eps 57 do i = 1, nctrs_pruned 58 59 x0i = xyz(1,i) 60 y0i = xyz(2,i) 61 z0i = xyz(3,i) 62 63 do n = 1,nq 64 x = qxyz(1,n) - x0i 65 y = qxyz(2,n) - y0i 66 z = qxyz(3,n) - z0i 67 68 rq(n,i) = sqrt(x*x + y*y + z*z) 69 enddo 70 enddo 71 call a_dist(xyz, Rij, nctrs_pruned,.false.) 72 AA=ictr 73c 74c find nearest neighb 75c 76 distnn=1.d+10 77 x0i=xyz(1,AA) 78 y0i=xyz(2,AA) 79 z0i=xyz(3,AA) 80 do i=1,nctrs_pruned 81 if(i.ne.AA) then 82 distnc=(xyz(1,i)-x0i)*(xyz(1,i)-x0i)+ 83 + (xyz(2,i)-y0i)*(xyz(2,i)-y0i)+ 84 + (xyz(3,i)-z0i)*(xyz(3,i)-z0i) 85 if(distnc.lt.distnn) then 86 distnn=distnc 87 NN=i 88 endif 89 endif 90 enddo 91 92 93 radwgh=(1.d0-ass)*sqrt(distnn)*0.5d0 94 do n = 1,nq 95c 96c check if grid point is within sphere where w=1 97c 98 if(rq(n,AA).ge.radwgh+eps) then 99 n1=n 100 goto 31 101 endif 102 enddo 103c all inside 104 return 105 31 continue 106 107 dist2=asse*sqrt(distnn)*4d0 108 n2=nq 109 do n=nq,n1,-1 110 if ((rq(n,AA)-rq(n,nn)).gt.dist2) then 111 qwght(n)=0d0 112 else 113 n2=n 114 goto 32 115 endif 116 enddo 11732 continue 118#ifdef USE_OPENMP 119!$omp parallel do 120!$omp& default(shared) 121!$omp& private(n,p,d_p,dzeta) 122#endif 123 do n=n1,n2 124 call dftg_gridssw0(n,nn,AA,nctrs_pruned,nq,nq_orig, 125 l grid_written, 126 c toll,distnn, 127 d qxyz,xyz,d_p,dzeta,p,rij,rq, 128 q qwght,d_qwght) 129 enddo ! n loop 130#ifdef USE_OPENMP 131!$omp end parallel do 132#endif 133 return 134 end 135 subroutine dftg_gridssw0(n,nn,AA,nctrs_pruned,nq,nq_orig, 136 l grid_written, 137 c toll,distnn, 138 d qxyz,xyz,d_p,dzeta,p,rij,rq, 139 q qwght,d_qwght) 140 implicit none 141 integer n 142 integer AA,nn,nctrs_pruned,nq,nq_orig 143 logical grid_written 144 double precision toll 145 double precision distnn 146 double precision qxyz(3,nq)! grid points [input] 147 double precision xyz(3,*) ! atom coordinates [input] 148 double precision d_p(3,nctrs_pruned), 149 & dzeta(3,nctrs_pruned) 150 double precision Rij(*) !interatomic distances [input] 151 double precision p(nctrs_pruned) 152 double precision rq(nq_orig,*) ! sum of p(n) [ output] 153 double precision qwght(nq_orig) ! weights [output] 154 double precision d_qwght(3,nq_orig,*) ! weight deriv [output] 155c 156 double precision x 157#include "xc_erftab.fh" 158c 159 integer i,j,ij,iind,jind 160 integer A,B 161 logical inota,jnota 162 double precision zetan 163 double precision mu,mu1,dmu1,dmu2,dmu3,dmu1dmu 164 double precision dskdmu1,tmu,tmu1 165 double precision rag,sk,wsum 166 double precision xa,ya,za 167 double precision xi,yi,zi 168 double precision dmuba(3),dbpa(3),damuab(3),dbmuba(3) 169c 20 digits ln(10)=2.3025 170c 0.5d0 factor because of mu^2?? 171 double precision undovl 172 parameter(undovl=20d0*2.3025d0*0.5d0) 173 double precision asqrtpi,alpha_erf,asse,ass 174 asqrtpi=1d0/sqrt(4*datan(1d0)) 175 ass=ass_erf1 176 asse=ass-eps 177 alpha_erf=alpha_erf1 178 zetan=0d0 179 call dfill(3*nctrs_pruned, 0.d0, dzeta, 1) 180 call dfill(3*nctrs_pruned, 0.d0, d_p, 1) 181 182c 183c compute mu_AN 184c 185 mu=(rq(n,AA)-rq(n,nn))/sqrt(distnn) 186 if (mu.gt.asse) then 187 p(AA)=0d0 188 zetan=1d0 189 goto 1100 190 endif 191 192 call dfill(nctrs_pruned,1.d0,p,1) 193 do i = 2, nctrs_pruned 194 inota=i.ne.AA 195 rag=rq(n,i) 196 ij = (i*(i-1))/2 197 do j = 1, i-1 198 199 jnota=j.ne.AA 200c 201 ij=ij+1 202 mu = (rag - rq(n,j))*Rij(ij) 203 if (mu.ge.asse) then 204 p(i)=0.d0 205 206 elseif (mu.le.-asse) then 207 p(j)=0.d0 208 209 else 210 if(inota.and.jnota) then 211c 212c use interpolation for erfs 213c 214 sk=erf1c(mu) 215 if(mu.lt.0d0) sk=1d0-sk 216 else 217 sk=erf1(mu) 218 endif 219 p(i) = p(i)*sk 220 p(j) = p(j)*(1d0 - sk) 221 endif 222 enddo ! end loop over j 223 enddo ! end loop over i 224c 225c compute sum of partitioning weights for normalization 226c 227c 228 wsum=0.d0 229 do i = 1, nctrs_pruned 230 wsum=wsum+p(i) 231 enddo 232 if(abs(wsum).lt.toll) goto 300 233 zetan = 1d0/wsum 234 1100 continue 235 do A = 1, nctrs_pruned 236 if(abs(p(A)).gt.toll) then 237 iind=A 238 inota=A.ne.AA 239 xA = (qxyz(1,n) - xyz(1,A))/rq(n,A) 240 yA = (qxyz(2,n) - xyz(2,A))/rq(n,A) 241 zA = (qxyz(3,n) - xyz(3,A))/rq(n,A) 242c 243c derivation variable B 244c 245 do B = 1, nctrs_pruned 246 jnota=B.ne.AA 247 if (A.ne.B)then 248 jind = B 249c 250 if (A.ge.B)then 251 ij = (A*(A-1))/2 + B 252 else 253 ij = (B*(B-1))/2 + A 254 endif 255c 256c 257 mu = (rq(n,A) - rq(n,B))*Rij(ij) 258 if(abs(mu).lt.asse) then 259 mu1=mu/(1d0-mu*mu)*alpha_erf 260 sk=erf1c(mu) 261 if(mu.lt.0d0) sk=1d0-sk 262 if(mu1.gt.undovl) then 263 dskdmu1=0d0 264 tmu=0d0 265 else 266 dmu1 = Rij(ij)*(xyz(1,A)-xyz(1,B)) 267 dmu2 = Rij(ij)*(xyz(2,A)-xyz(2,B)) 268 dmu3 = Rij(ij)*(xyz(3,A)-xyz(3,B)) 269 dmu1dmu=((mu*mu+1d0)/(1d0-mu*mu)**2)*alpha_erf 270 dskdmu1=-exp(-mu1*mu1)*asqrtpi 271 tmu = dskdmu1*dmu1dmu 272 endif 273 else 274 tmu=0d0 275 sk=0d0 276 endif 277c 278 279 if(abs(sk).gt.toll.and.abs(tmu).gt.0d0) then 280c 281c compute D(B)mu(AB) 282c 283 xi = qxyz(1,n) - xyz(1,B) 284 yi = qxyz(2,n) - xyz(2,B) 285 zi = qxyz(3,n) - xyz(3,B) 286c 287c atomic size adjustment derivative 288c 289 dbmuba(1) = -(xi/rq(n,B) + mu*dmu1)*Rij(ij) 290 dbmuba(2) = -(yi/rq(n,B) + mu*dmu2)*Rij(ij) 291 dbmuba(3) = -(zi/rq(n,B) + mu*dmu3)*Rij(ij) 292c 293 tmu1=tmu*p(A)/sk 294c 295c term \Delta_B PA 296c 297 dBPA(1)= -tmu1*dbmuba(1) 298 dBPA(2)= -tmu1*dbmuba(2) 299 dBPA(3)= -tmu1*dbmuba(3) 300c 301 dzeta(1,B) = dzeta(1,B)+ dBPA(1) 302 dzeta(2,B) = dzeta(2,B)+ dBPA(2) 303 dzeta(3,B) = dzeta(3,B)+ dBPA(3) 304c 305 306 if (inota)then 307c 308c term \Delta_A PA (partial) 309c 310c 311c compute D(A)mu(AB) 312c 313 damuab(1) = -(xA+mu*dmu1)*Rij(ij) 314 damuab(2) = -(yA+mu*dmu2)*Rij(ij) 315 damuab(3) = -(zA+mu*dmu3)*Rij(ij) 316 dzeta(1,A) = dzeta(1,A)+tmu1*damuab(1) 317 dzeta(2,A) = dzeta(2,A)+tmu1*damuab(2) 318 dzeta(3,A) = dzeta(3,A)+tmu1*damuab(3) 319 else 320 d_p(1,B) = dBPA(1) 321 d_p(2,B) = dBPA(2) 322 d_p(3,B) = dBPA(3) 323 endif 324 endif 325 endif 326 enddo ! B loop 327 endif 328 enddo ! A loop 329 300 continue 330 if(.not.grid_written) then 331 if(abs(p(AA)).gt.toll)then 332 qwght(n) = (p(AA)*qwght(n))*zetan 333 else 334 qwght(n)=0d0 335 endif 336 endif 337c 338c compute \Delta_i W_ictr 339c 340c \Delta_B PA -\delta_B Z*PA/Z 341c 342 do B = 1, nctrs_pruned 343 if (B.ne.AA.and.abs(p(AA)).gt.toll) then 344 d_qwght(1,n,b)= (d_p(1,B)/p(AA) - 345 & dzeta(1,B)*zetan)*qwght(n) 346 d_qwght(2,n,b)= (d_p(2,B)/p(AA) - 347 & dzeta(2,B)*zetan)*qwght(n) 348 d_qwght(3,n,b)= (d_p(3,B)/p(AA) - 349 & dzeta(3,B)*zetan)*qwght(n) 350 endif 351 enddo 352 return 353 end 354 355