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