1*
2* $Id$
3*
4
5
6      subroutine integrate_kbppv3(version,rlocal,
7     >                            nrho,drho,lmax,locp,zv,
8     >                            vp,wp,rho,f,cs,sn,
9     >                            nfft1,nfft2,nfft3,lmmax,
10     >                            G,vl,vnl,
11     >                            n_prj,l_prj,m_prj,b_prj,vnlnrm,
12     >                            semicore,rho_sc_r,rho_sc_k,
13     >                            ierr)
14      implicit none
15      integer          version
16      double precision rlocal
17      integer          nrho
18      double precision drho
19      integer          lmax
20      integer          locp
21      double precision zv
22      double precision vp(nrho,0:lmax)
23      double precision wp(nrho,0:lmax)
24      double precision rho(nrho)
25      double precision f(nrho)
26      double precision cs(nrho)
27      double precision sn(nrho)
28
29      integer nfft1,nfft2,nfft3,lmmax
30      double precision G(nfft1/2+1,nfft2,nfft3,3)
31      double precision vl(nfft1/2+1,nfft2,nfft3)
32      double precision vnl(nfft1/2+1,nfft2,nfft3,lmmax)
33      integer          n_prj(lmmax),l_prj(lmmax),m_prj(lmmax)
34      integer          b_prj(lmmax)
35      double precision vnlnrm(0:lmax)
36
37      logical semicore
38      double precision rho_sc_r(nrho,2)
39      double precision rho_sc_k(nfft1/2+1,nfft2,nfft3,4)
40
41      integer ierr
42
43      integer np,taskid,MASTER
44      parameter (MASTER=0)
45
46*     *** local variables ****
47      logical fast_erf
48      integer lcount,task_count,nfft3d
49      integer k1,k2,k3,i,l
50      double precision pi,twopi,forpi
51      double precision p0,p1,p2,p3,p
52      double precision gx,gy,gz,a,q,d
53
54*     **** Error function parameters ****
55      real*8 yerf,xerf
56      real*8 c1,c2,c3,c4,c5,c6
57      parameter (c1=0.07052307840d0,c2=0.04228201230d0)
58      parameter (c3=0.00927052720d0)
59      parameter (c4=0.00015201430d0,c5=0.00027656720d0)
60      parameter (c6=0.00004306380d0)
61
62*     **** external functions ****
63      logical          control_fast_erf
64      external         control_fast_erf
65      double precision dsum,simp,util_erf
66      external         dsum,simp,util_erf
67
68      call Parallel_np(np)
69      call Parallel_taskid(taskid)
70
71      fast_erf = control_fast_erf()
72
73      nfft3d = (nfft1/2+1)*nfft2*nfft3
74      pi=4.0d0*datan(1.0d0)
75      twopi=2.0d0*pi
76      forpi=4.0d0*pi
77
78      IF(LMMAX.GT.16) THEN
79        IERR=1
80        RETURN
81      ENDIF
82      IF((NRHO/2)*2.EQ.NRHO) THEN
83        IERR=2
84        RETURN
85      ENDIF
86
87      P0=DSQRT(FORPI)
88      P1=DSQRT(3.0d0*FORPI)
89      P2=DSQRT(15.0d0*FORPI)
90      P3=DSQRT(105.0d0*FORPI)
91
92*::::::::::::::::::  Define non-local pseudopotential  ::::::::::::::::
93      do l=0,lmax
94        if (l.ne.locp) then
95          do I=1,nrho
96            vp(i,l)=vp(i,l)-vp(i,locp)
97          end do
98        end if
99      end do
100
101*:::::::::::::::::::::  Normarization constants  ::::::::::::::::::::::
102      lcount = 0
103      do l=0,lmax
104        if (l.ne.locp) then
105          do i=1,nrho
106            f(I)=vp(I,L)*wp(I,L)**2
107          end do
108          a=simp(nrho,f,drho)
109          vnlnrm(l) = (1.0d0/a)
110        else
111          vnlnrm(l) = 0.0d0
112        end if
113      end do
114
115*======================  Fourier transformation  ======================
116      call dcopy(nfft3d,0.0d0,0,vl,1)
117      call dcopy(lmmax*nfft3d,0.0d0,0,vnl,1)
118      call dcopy(4*nfft3d,0.0d0,0,rho_sc_k,1)
119      task_count = -1
120      DO 700 k3=1,nfft3
121      DO 700 k2=1,nfft2
122      DO 700 k1=1,(nfft1/2+1)
123        task_count = task_count + 1
124        if (mod(task_count,np).ne.taskid) go to 700
125
126        Q=DSQRT(G(k1,k2,k3,1)**2
127     >         +G(k1,k2,k3,2)**2
128     >         +G(k1,k2,k3,3)**2)
129
130        if ((k1.eq.1).and.(k2.eq.1).and.(k3.eq.1)) go to 700
131
132
133        GX=G(k1,k2,k3,1)/Q
134        GY=G(k1,k2,k3,2)/Q
135        GZ=G(k1,k2,k3,3)/Q
136        DO I=1,NRHO
137          CS(I)=DCOS(Q*RHO(I))
138          SN(I)=DSIN(Q*RHO(I))
139        END DO
140
141        lcount = lmmax+1
142        GO TO (500,400,300,200), LMAX+1
143
144
145*::::::::::::::::::::::::::::::  f-wave  ::::::::::::::::::::::::::::::
146  200   CONTINUE
147        if (locp.ne.3) then
148           F(1)=0.0d0
149           do I=2,NRHO
150             A=SN(I)/(Q*RHO(I))
151             A=15.0d0*(A-CS(I))/(Q*RHO(I))**2 - 6*A + CS(I)
152             F(I)=A*WP(I,3)*VP(I,3)
153           end do
154           D=P3*SIMP(NRHO,F,DRHO)/Q
155
156           lcount = lcount-1
157           vnl(k1,k2,k3,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY)
158     >                          /dsqrt(24.0d0)
159           lcount = lcount-1
160           vnl(k1,k2,k3,lcount)=D*GX*GY*GZ
161           lcount = lcount-1
162           vnl(k1,k2,k3,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0)
163     >                          /dsqrt(40.0d0)
164           lcount = lcount-1
165           vnl(k1,k2,k3,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0)
166     >                          /dsqrt(60.0d0)
167           lcount = lcount-1
168           vnl(k1,k2,k3,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0)
169     >                          /dsqrt(40.0d0)
170           lcount = lcount-1
171           vnl(k1,k2,k3,lcount)=D*GZ*(GX*GX - GY*GY)
172     >                          /2.0d0
173           lcount = lcount-1
174           vnl(k1,k2,k3,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ))
175     >                          /dsqrt(24.0d0)
176c           lcount = lcount-1
177c           vnl(k1,k2,k3,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ))
178c     >                          /dsqrt(24.0d0)
179c           lcount = lcount-1
180c           vnl(k1,k2,k3,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY)
181c     >                          /dsqrt(24.0d0)
182c           lcount = lcount-1
183c           vnl(k1,k2,k3,lcount)=D*GZ*(GX*GX - GY*GY)
184c     >                          /2.0d0
185c           lcount = lcount-1
186c           vnl(k1,k2,k3,lcount)=D*GX*GY*GZ
187c           lcount = lcount-1
188c           vnl(k1,k2,k3,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0)
189c     >                          /dsqrt(40.0d0)
190c           lcount = lcount-1
191c           vnl(k1,k2,k3,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0)
192c     >                          /dsqrt(40.0d0)
193c           lcount = lcount-1
194c           vnl(k1,k2,k3,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0)
195c     >                          /dsqrt(60.0d0)
196        end if
197
198
199
200*::::::::::::::::::::::::::::::  d-wave  ::::::::::::::::::::::::::::::
201  300   CONTINUE
202        if (locp.ne.2) then
203          F(1)=0.0d0
204          DO I=2,NRHO
205            A=3.0d0*(SN(I)/(Q*RHO(I))-CS(I))/(Q*RHO(I))-SN(I)
206            F(I)=A*WP(I,2)*VP(I,2)
207          END DO
208          D=P2*SIMP(NRHO,F,DRHO)/Q
209
210          lcount = lcount-1
211          vnl(k1,k2,k3,lcount)=D*GX*GY
212          lcount = lcount-1
213          vnl(k1,k2,k3,lcount)=D*GY*GZ
214          lcount = lcount-1
215          vnl(k1,k2,k3,lcount)=D*(3.0d0*GZ*GZ-1.0d0)
216     >                          /(2.0d0*dsqrt(3.0d0))
217          lcount = lcount-1
218          vnl(k1,k2,k3,lcount)=D*GZ*GX
219          lcount = lcount-1
220          vnl(k1,k2,k3,lcount)=D*(GX*GX-GY*GY)/(2.0d0)
221
222c          lcount = lcount-1
223c          vnl(k1,k2,k3,lcount)=D*(3.0d0*GZ*GZ-1.0d0)
224c     >                          /(2.0d0*dsqrt(3.0d0))
225c          lcount = lcount-1
226c          vnl(k1,k2,k3,lcount)=D*GX*GY
227c          lcount = lcount-1
228c          vnl(k1,k2,k3,lcount)=D*GY*GZ
229c          lcount = lcount-1
230c          vnl(k1,k2,k3,lcount)=D*GZ*GX
231c          lcount = lcount-1
232c          vnl(k1,k2,k3,lcount)=D*(GX*GX-GY*GY)/(2.0d0)
233        end if
234
235*::::::::::::::::::::::::::::::  p-wave  ::::::::::::::::::::::::::::::
236  400   CONTINUE
237        if (locp.ne.1) then
238           F(1)=0.0d0
239           DO I=2,NRHO
240             F(I)=(SN(I)/(Q*RHO(I))-CS(I))*WP(I,1)*VP(I,1)
241           END DO
242           P=P1*SIMP(NRHO,F,DRHO)/Q
243           lcount = lcount-1
244           vnl(k1,k2,k3,lcount)=P*GY
245           lcount = lcount-1
246           vnl(k1,k2,k3,lcount)=P*GZ
247           lcount = lcount-1
248           vnl(k1,k2,k3,lcount)=P*GX
249c           lcount = lcount-1
250c           vnl(k1,k2,k3,lcount)=P*GX
251c           lcount = lcount-1
252c           vnl(k1,k2,k3,lcount)=P*GY
253c           lcount = lcount-1
254c           vnl(k1,k2,k3,lcount)=P*GZ
255        end if
256
257*::::::::::::::::::::::::::::::  s-wave  :::::::::::::::::::::::::::::::
258  500   CONTINUE
259        if (locp.ne.0) then
260          DO I=1,NRHO
261            F(I)=SN(I)*WP(I,0)*VP(I,0)
262          END DO
263          lcount = lcount-1
264          vnl(k1,k2,k3,lcount)=P0*SIMP(NRHO,F,DRHO)/Q
265        end if
266
267*::::::::::::::::::::::::::::::  local  :::::::::::::::::::::::::::::::
268  600   CONTINUE
269
270
271        if (version.eq.3) then
272        DO  I=1,NRHO
273          F(I)=RHO(I)*VP(I,locp)*SN(I)
274        END DO
275        vl(k1,k2,k3)=SIMP(NRHO,F,DRHO)*FORPI/Q-ZV*FORPI*CS(NRHO)/(Q*Q)
276        end if
277
278        if (version.eq.4) then
279        if (fast_erf) then
280           do I=1,NRHO
281             xerf=RHO(I)/rlocal
282             yerf = (1.0d0
283     >            + xerf*(c1 + xerf*(c2
284     >            + xerf*(c3 + xerf*(c4
285     >            + xerf*(c5 + xerf*c6))))))**4
286             yerf = (1.0d0 - 1.0d0/yerf**4)
287             F(I)=(RHO(I)*VP(I,locp)+ZV*yerf)*SN(I)
288           end do
289        else
290           do I=1,NRHO
291             xerf=RHO(I)/rlocal
292             yerf = util_erf(xerf)
293             F(I)=(RHO(I)*VP(I,locp)+ZV*yerf)*SN(I)
294           end do
295        end if
296        vl(k1,k2,k3)=SIMP(NRHO,F,DRHO)*FORPI/Q
297        end if
298
299
300*::::::::::::::::::::: semicore density :::::::::::::::::::::::::::::::
301        if (semicore) then
302           do i=1,nrho
303              f(i) = rho(i)*dsqrt(rho_sc_r(i,1))*sn(i)
304           end do
305           rho_sc_k(k1,k2,k3,1) = SIMP(nrho,f,drho)*forpi/Q
306
307           do i=1,nrho
308             f(i)=(sn(i)/(Q*rho(i))-cs(i))*rho_sc_r(i,2)*rho(i)
309           end do
310           P = SIMP(nrho,f,drho)*forpi/Q
311           rho_sc_k(k1,k2,k3,2)=P*GX
312           rho_sc_k(k1,k2,k3,3)=P*GY
313           rho_sc_k(k1,k2,k3,4)=P*GZ
314
315        end if
316
317  700 CONTINUE
318
319      call Parallel_Vector_SumAll(4*nfft3d,rho_sc_k)
320      call Parallel_Vector_SumAll(nfft3d,vl)
321      call Parallel_Vector_Sumall(lmmax*nfft3d,vnl)
322*:::::::::::::::::::::::::::::::  G=0  ::::::::::::::::::::::::::::::::
323
324      if (version.eq.3) then
325      DO I=1,NRHO
326        F(I)=VP(I,locp)*RHO(I)**2
327      END DO
328      vl(1,1,1)=FORPI*SIMP(NRHO,F,DRHO)+TWOPI*ZV*RHO(NRHO)**2
329      end if
330
331      if (version.eq.4) then
332      if (fast_erf) then
333         do I=1,NRHO
334           xerf=RHO(I)/rlocal
335           yerf = (1.0d0
336     >          + xerf*(c1 + xerf*(c2
337     >          + xerf*(c3 + xerf*(c4
338     >          + xerf*(c5 + xerf*c6))))))**4
339           yerf = (1.0d0 - 1.0d0/yerf**4)
340           F(I)=(VP(I,locp)*RHO(I)+ZV*yerf)*RHO(I)
341         end do
342      else
343         do I=1,NRHO
344           xerf=RHO(I)/rlocal
345           yerf = util_erf(xerf)
346           F(I)=(VP(I,locp)*RHO(I)+ZV*yerf)*RHO(I)
347         end do
348      end if
349      vl(1,1,1)=FORPI*SIMP(NRHO,F,DRHO)
350      end if
351
352*     **** semicore density ****
353      if (semicore) then
354         do i=1,nrho
355            f(i) = dsqrt(rho_sc_r(i,1))*rho(i)**2
356         end do
357         rho_sc_k(1,1,1,1) = forpi*SIMP(nrho,f,drho)
358         rho_sc_k(1,1,1,2) = 0.0d0
359         rho_sc_k(1,1,1,3) = 0.0d0
360         rho_sc_k(1,1,1,4) = 0.0d0
361      end if
362
363      do l=1,lmmax
364        vnl(1,1,1,l)=0.0d0
365      end do
366*     *** only j0 is non-zero at zero ****
367      if (locp.ne.0) then
368         DO  I=1,NRHO
369           F(I)=RHO(I)*WP(I,0)*VP(I,0)
370         END DO
371         vnl(1,1,1,1)=P0*SIMP(NRHO,F,DRHO)
372      end if
373
374
375*     ********************************
376*     **** define n_prj and l_prj ****
377*     ********************************
378      lcount = lmmax+1
379      GO TO (950,940,930,920), lmax+1
380
381        !::::::  f-wave  :::::::
382  920   CONTINUE
383        if (locp.ne.3) then
384          lcount = lcount-1
385          n_prj(lcount) = 1
386          l_prj(lcount) = 3
387          m_prj(lcount) = -3
388          b_prj(lcount) = 4
389
390          lcount = lcount-1
391          n_prj(lcount) = 1
392          l_prj(lcount) = 3
393          m_prj(lcount) = -2
394          b_prj(lcount) = 4
395
396          lcount = lcount-1
397          n_prj(lcount) = 1
398          l_prj(lcount) = 3
399          m_prj(lcount) = -1
400          b_prj(lcount) = 4
401
402          lcount = lcount-1
403          n_prj(lcount) = 1
404          l_prj(lcount) = 3
405          m_prj(lcount) = 0
406          b_prj(lcount) = 4
407
408          lcount = lcount-1
409          n_prj(lcount) = 1
410          l_prj(lcount) = 3
411          m_prj(lcount) = 1
412          b_prj(lcount) = 4
413
414          lcount = lcount-1
415          n_prj(lcount) = 1
416          l_prj(lcount) = 3
417          m_prj(lcount) = 2
418          b_prj(lcount) = 4
419
420          lcount = lcount-1
421          n_prj(lcount) = 1
422          l_prj(lcount) = 3
423          m_prj(lcount) = 3
424          b_prj(lcount) = 4
425        end if
426
427
428        !::::  d-wave  ::::
429  930   CONTINUE
430        if (locp.ne.2) then
431          lcount = lcount-1
432          n_prj(lcount) = 1
433          l_prj(lcount) = 2
434          m_prj(lcount) = -2
435          b_prj(lcount) = 3
436
437          lcount = lcount-1
438          n_prj(lcount) = 1
439          l_prj(lcount) = 2
440          m_prj(lcount) = -1
441          b_prj(lcount) = 3
442
443          lcount = lcount-1
444          n_prj(lcount) = 1
445          l_prj(lcount) = 2
446          m_prj(lcount) = 0
447          b_prj(lcount) = 3
448
449          lcount = lcount-1
450          n_prj(lcount) = 1
451          l_prj(lcount) = 2
452          m_prj(lcount) = 1
453          b_prj(lcount) = 3
454
455          lcount = lcount-1
456          n_prj(lcount) = 1
457          l_prj(lcount) = 2
458          m_prj(lcount) = 2
459          b_prj(lcount) = 3
460        end if
461
462
463        !::::  p-wave  ::::
464  940   CONTINUE
465        if (locp.ne.1) then
466          lcount = lcount-1
467          n_prj(lcount) = 1
468          l_prj(lcount) = 1
469          m_prj(lcount) = -1
470          b_prj(lcount) = 2
471
472          lcount = lcount-1
473          n_prj(lcount) = 1
474          l_prj(lcount) = 1
475          m_prj(lcount) = 0
476          b_prj(lcount) = 2
477
478          lcount = lcount-1
479          n_prj(lcount) = 1
480          l_prj(lcount) = 1
481          m_prj(lcount) = 1
482          b_prj(lcount) = 2
483        end if
484
485
486        !::::  s-wave  ::::
487  950   CONTINUE
488        if (locp.ne.0) then
489          lcount = lcount-1
490          n_prj(lcount) = 1
491          l_prj(lcount) = 0
492          m_prj(lcount) = 0
493          b_prj(lcount) = 1
494        end if
495
496
497      IERR=0
498      RETURN
499      END
500
501
502
503