1*
2* $Id$
3*
4      subroutine integrate_kbpp_band_orb(version,kvec,
5     >                            nrho,drho,lmax,locp,
6     >                            wp,rho,f,cs,sn,
7     >                            nfft1,nfft2,nfft3,lmmax,
8     >                            G,borbs,
9     >                            ierr)
10      implicit none
11      integer          version
12      double precision kvec(3)
13      integer          nrho
14      double precision drho
15      integer          lmax
16      integer          locp
17      double precision wp(nrho,0:lmax)
18      double precision rho(nrho)
19      double precision f(nrho)
20      double precision cs(nrho)
21      double precision sn(nrho)
22
23      integer nfft1,nfft2,nfft3,lmmax
24      double precision G(nfft1,nfft2,nfft3,3)
25      double precision borbs(nfft1,nfft2,nfft3,lmmax)
26
27      integer ierr
28
29      integer np,taskid,MASTER
30      parameter (MASTER=0)
31
32#include "bafdecls.fh"
33
34*     *** local variables ****
35      integer lcount,task_count,nfft3d
36      integer k1,k2,k3,i,l
37      double precision pi,twopi,forpi
38      double precision p0,p1,p2,p3,p
39      double precision gx,gy,gz,a,q,d
40
41
42*     **** external functions ****
43      double precision dsum,simp
44      external         dsum,simp
45      logical value
46
47      call Parallel_np(np)
48      call Parallel_taskid(taskid)
49
50      nfft3d = (nfft1)*nfft2*nfft3
51      pi=4.0d0*datan(1.0d0)
52      twopi=2.0d0*pi
53      forpi=4.0d0*pi
54
55      IF(LMMAX.GT.16) THEN
56        WRITE(*,*)"non-local psp not generated: lmmax exceeds 16"
57        IERR=1
58        RETURN
59      ENDIF
60      IF((NRHO/2)*2.EQ.NRHO) THEN
61        WRITE(*,*)"non-local psp not generated: nrho is not odd"
62        IERR=2
63        RETURN
64      ENDIF
65
66      P0=DSQRT(FORPI)
67      P1=DSQRT(3.0d0*FORPI)
68      P2=DSQRT(15.0d0*FORPI)
69      P3=DSQRT(105.0d0*FORPI)
70
71
72*======================  Fourier transformation  ======================
73      call dcopy(lmmax*nfft3d,0.0d0,0,borbs,1)
74
75      task_count = -1
76      DO 700 k3=1,nfft3
77      DO 700 k2=1,nfft2
78      DO 700 k1=1,nfft1
79        task_count = task_count + 1
80        if (mod(task_count,np).ne.taskid) go to 700
81        gx=G(k1,k2,k3,1)+kvec(1)
82        gy=G(k1,k2,k3,2)+kvec(2)
83        gz=G(k1,k2,k3,3)+kvec(3)
84
85        Q=DSQRT(gx**2 + gy**2 + gz**2)
86
87        if (dabs(Q).gt.1.0d-9) then
88
89           gx=gx/Q
90           gy=gy/Q
91           gz=gz/Q
92           DO I=1,NRHO
93             CS(I)=DCOS(Q*RHO(I))
94             SN(I)=DSIN(Q*RHO(I))
95           END DO
96
97           lcount = lmmax+1
98           GO TO (500,400,300,200,100), LMAX+1
99
100
101*::::::::::::::::::::::::::::::  g-wave  ::::::::::::::::::::::::::::::
102  100      CONTINUE
103*::::::::::::::::::::::::::::::  f-wave  ::::::::::::::::::::::::::::::
104  200      CONTINUE
105           if (locp.ne.3) then
106              F(1)=0.0d0
107              do I=2,NRHO
108                A=SN(I)/(Q*RHO(I))
109                A=15.0d0*(A-CS(I))/(Q*RHO(I))**2 - 6*A + CS(I)
110                F(I)=A*WP(I,3)
111              end do
112              D=P3*SIMP(NRHO,F,DRHO)/Q
113
114*          **** fy(3x2-y2) component ****
115           lcount = lcount-1
116           borbs(k1,k2,k3,lcount)=D*GY*(3.0d0*(1.0d0-GZ*GZ)-4.0d0*GY*GY)
117     >                          /dsqrt(24.0d0)
118
119*          **** fxyz component ****
120           lcount = lcount-1
121           borbs(k1,k2,k3,lcount)=D*GX*GY*GZ
122
123*             **** fy(5z2-1) component ****
124           lcount = lcount-1
125           borbs(k1,k2,k3,lcount)=D*GY*(5.0d0*GZ*GZ-1.0d0)
126     >                          /dsqrt(40.0d0)
127
128*          **** fz(5z2-3) component ****
129           lcount = lcount-1
130           borbs(k1,k2,k3,lcount)=D*GZ*(5.0d0*GZ*GZ-3.0d0)
131     >                          /dsqrt(60.0d0)
132
133*          **** fx(5z2-1) component ****
134           lcount = lcount-1
135           borbs(k1,k2,k3,lcount)=D*GX*(5.0d0*GZ*GZ-1.0d0)
136     >                          /dsqrt(40.0d0)
137
138*          **** fz(x2-y2) component ****
139           lcount = lcount-1
140           borbs(k1,k2,k3,lcount)=D*GZ*(GX*GX - GY*GY)
141     >                          /2.0d0
142
143*          **** fx(x2-3y2) component ****
144           lcount = lcount-1
145           borbs(k1,k2,k3,lcount)=D*GX*(4.0d0*GX*GX-3.0d0*(1.0d0-GZ*GZ))
146     >                          /dsqrt(24.0d0)
147           end if
148
149
150*::::::::::::::::::::::::::::::  d-wave  ::::::::::::::::::::::::::::::
151  300      CONTINUE
152           if (locp.ne.2) then
153             F(1)=0.0d0
154             DO I=2,NRHO
155               A=3.0d0*(SN(I)/(Q*RHO(I))-CS(I))/(Q*RHO(I))-SN(I)
156               F(I)=A*WP(I,2)
157             END DO
158             D=P2*SIMP(NRHO,F,DRHO)/Q
159
160*         **** dxy component ****
161          lcount = lcount-1
162          borbs(k1,k2,k3,lcount)=D*GX*GY
163
164*         **** dyz component ****
165          lcount = lcount-1
166          borbs(k1,k2,k3,lcount)=D*GY*GZ
167
168*         **** d3z2-1 component ****
169          lcount = lcount-1
170          borbs(k1,k2,k3,lcount)=D*(3.0d0*GZ*GZ-1.0d0)
171     >                          /(2.0d0*dsqrt(3.0d0))
172
173*         **** dzx component ****
174          lcount = lcount-1
175          borbs(k1,k2,k3,lcount)=D*GZ*GX
176
177*         **** dx2-y2 component ****
178          lcount = lcount-1
179          borbs(k1,k2,k3,lcount)=D*(GX*GX-GY*GY)/(2.0d0)
180
181           end if
182
183*::::::::::::::::::::::::::::::  p-wave  ::::::::::::::::::::::::::::::
184  400      CONTINUE
185           if (locp.ne.1) then
186              F(1)=0.0d0
187              DO I=2,NRHO
188                F(I)=(SN(I)/(Q*RHO(I))-CS(I))*WP(I,1)
189              END DO
190              P=P1*SIMP(NRHO,F,DRHO)/Q
191
192*          **** py component ****
193           lcount = lcount-1
194           borbs(k1,k2,k3,lcount)=P*GY
195
196*          **** pz component ****
197           lcount = lcount-1
198           borbs(k1,k2,k3,lcount)=P*GZ
199
200*          **** px component ****
201           lcount = lcount-1
202           borbs(k1,k2,k3,lcount)=P*GX
203           end if
204
205*::::::::::::::::::::::::::::::  s-wave  :::::::::::::::::::::::::::::::
206  500      CONTINUE
207           if (locp.ne.0) then
208             DO I=1,NRHO
209               F(I)=SN(I)*WP(I,0)
210             END DO
211             lcount = lcount-1
212             borbs(k1,k2,k3,lcount)=P0*SIMP(NRHO,F,DRHO)/Q
213           end if
214
215  600      CONTINUE
216
217*:::::::::::::::::::::::::::::::  G+k=0  ::::::::::::::::::::::::::::::::
218      else
219
220         do l=1,lmmax
221           borbs(k1,k2,k3,l)=0.0d0
222         end do
223*        *** only j0 is non-zero at zero ****
224         if (locp.ne.0) then
225            DO  I=1,NRHO
226              F(I)=RHO(I)*WP(I,0)
227            END DO
228            borbs(k1,k2,k3,1)=P0*SIMP(NRHO,F,DRHO)
229         end if
230
231      end if
232
233
234  700 CONTINUE
235
236      call Parallel_Vector_SumAll(lmmax*nfft3d,borbs)
237
238      IERR=0
239      RETURN
240      END
241
242
243
244