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