1*
2* $Id$
3*
4
5*     **********************************************
6*     *                                            *
7*     *   integrate_kbpp_band_stress_nonlocal_new  *
8*     *                                            *
9*     **********************************************
10
11      subroutine integrate_kbpp_band_stress_nonlocal_new(version,kvec,
12     >                            nrho,drho,lmax,locp,zv,
13     >                            vp,wp,rho,f,cs,sn,
14     >                            nfft1,nfft2,nfft3,lmmax,
15     >                            G,dvnl,
16     >                            nray,G_ray,dvnl_ray,
17     >                            ierr)
18      implicit none
19      integer          version
20      double precision kvec(3)
21      integer          nrho
22      double precision drho
23      integer          lmax
24      integer          locp
25      double precision zv
26      double precision vp(nrho,0:lmax)
27      double precision wp(nrho,0:lmax)
28      double precision rho(nrho)
29      double precision f(nrho)
30      double precision cs(nrho)
31      double precision sn(nrho)
32
33      integer nfft1,nfft2,nfft3,lmmax
34      double precision G(nfft1,nfft2,nfft3,3)
35      double precision dvnl(nfft1,nfft2,nfft3,3,lmmax)
36
37      integer nray
38      double precision G_ray(nray)
39      double precision dvnl_ray(nray,2,0:lmax,2)
40
41      integer ierr
42
43#include "errquit.fh"
44
45*     *** local variables ****
46      integer np,taskid,MASTER
47      parameter (MASTER=0)
48
49      integer lcount,task_count,nfft3d
50      integer k1,k2,k3,i,l,nx
51      double precision p,pp
52      double precision gx,gy,gz,q,d,dd
53      double precision duxdGx,duxdGy,duxdGz
54      double precision duydGx,duydGy,duydGz
55      double precision duzdGx,duzdGy,duzdGz
56      double precision sumx,sumy,sumz,dG
57      double precision T,dTdux,dTduy,dTduz
58
59
60
61*     **** external functions ****
62      double precision dsum,simp,nwpw_splint
63      external         dsum,simp,nwpw_splint
64
65      call Parallel_np(np)
66      call Parallel_taskid(taskid)
67
68      if(lmmax.gt.16) then
69         call errquit(
70     >    'integrate_kbpp_band_stress_nonlocal_new - lmax > f',0,
71     >       INPUT_ERR)
72      end if
73      if((nrho/2)*2.eq.nrho) then
74        call errquit(
75     >  'integrate_kbpp_band_stress_nonlocal_new - psp grid not odd',0,
76     >       INPUT_ERR)
77      end if
78
79      nfft3d = (nfft1)*nfft2*nfft3
80
81*======================  Fourier transformation  ======================
82      dG = G_ray(3)-G_ray(2)
83      call dcopy(3*lmmax*nfft3d,0.0d0,0,dvnl,1)
84      task_count = -1
85      do 700 k3=1,nfft3
86      do 700 k2=1,nfft2
87      do 700 k1=1,nfft1
88        task_count = task_count + 1
89        if (mod(task_count,np).ne.taskid) go to 700
90        gx=G(k1,k2,k3,1)+kvec(1)
91        gy=G(k1,k2,k3,2)+kvec(2)
92        gz=G(k1,k2,k3,3)+kvec(3)
93
94        q=dsqrt(gx**2 + gy**2 + gz**2)
95        nx = (q/dG) + 1.0d0
96
97        if (dabs(q).gt.1.0d-9) then
98
99           gx=gx/q
100           gy=gy/q
101           gz=gz/q
102           do i=1,nrho
103             cs(i)=dcos(q*rho(i))
104             sn(i)=dsin(q*rho(i))
105           end do
106
107*          **** calculate du_r/dG_s ****
108           duxdGx = 1.0d0/q -gx*gx/q
109           duxdGy = -gx*gy/q
110           duxdGz = -gx*gz/q
111
112           duydGx = -gy*gx/q
113           duydGy = 1.0d0/q - gy*gy/q
114           duydGz = -gy*gz/q
115
116           duzdGx = -gz*gx/q
117           duzdGy = -gz*gy/q
118           duzdGz = 1.0d0/q - gz*gz/q
119
120           lcount = lmmax+1
121           GO TO (500,400,300,200), LMAX+1
122
123*::::::::::::::::::::::::::::::  f-wave  ::::::::::::::::::::::::::::::
124  200      CONTINUE
125        if (locp.ne.3) then
126           D  = nwpw_splint(G_ray,dvnl_ray(1,1,3,1),
127     >                            dvnl_ray(1,1,3,2),nray,nx,Q)
128           DD = nwpw_splint(G_ray,dvnl_ray(1,2,3,1),
129     >                            dvnl_ray(1,2,3,2),nray,nx,Q)
130           lcount = lcount-1
131           T = gx*(4.0d0*gx*gx - 3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0)
132           dTdux = (12.0d0*gx*gx-3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0)
133           dTduy = 0.0d0
134           dTduz = 6.0d0*gx*gz/dsqrt(24.0d0)
135           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
136           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
137           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
138           dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
139           dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
140           dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
141
142           lcount = lcount-1
143           T = gy*(3.0d0*(1.0d0-gz*gz)-4.0d0*gy*gy)/dsqrt(24.0d0)
144           dTdux = 0.0d0
145           dTduy = (3.0d0*(1.0d0-gz*gz)-12.0d0*gy*gy)/dsqrt(24.0d0)
146           dTduz = -6.0d0*gy*gz/dsqrt(24.0d0)
147           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
148           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
149           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
150           dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
151           dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
152           dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
153
154           lcount = lcount-1
155           T =gz*(gx*gx - gy*gy)/2.0d0
156           dTdux =  gx*gz
157           dTduy = -gy*gz
158           dTduz = (gx*gx-gy*gy)/2.0d0
159           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
160           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
161           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
162           dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
163           dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
164           dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
165
166           lcount = lcount-1
167           T =gx*gy*gz
168           dTdux = gy*gz
169           dTduy = gx*gz
170           dTduz = gx*gy
171           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
172           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
173           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
174           dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
175           dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
176           dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
177
178           lcount = lcount-1
179           T =  gx*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
180           dTdux = (5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
181           dTduy = 0.0d0
182           dTduz = 10.0d0*gx*gz/dsqrt(40.0d0)
183           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
184           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
185           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
186           dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
187           dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
188           dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
189
190           lcount = lcount-1
191           T = gy*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
192           dTdux = 0.0d0
193           dTduy =(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
194           dTduz =10.0d0*gy*gz/dsqrt(40.0d0)
195           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
196           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
197           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
198           dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
199           dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
200           dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
201
202           lcount = lcount-1
203           T =gz*(5.0d0*gz*gz-3.0d0)/dsqrt(60.0d0)
204           dTdux = 0.0d0
205           dTduy = 0.0d0
206           dTduz =(15.0d0*gz*gz -3.0d0)/dsqrt(60.0d0)
207           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
208           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
209           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
210           dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
211           dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
212           dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
213        end if
214*::::::::::::::::::::::::::::::  d-wave  ::::::::::::::::::::::::::::::
215  300      CONTINUE
216        if (locp.ne.2) then
217          D  = nwpw_splint(G_ray,dvnl_ray(1,1,2,1),
218     >                           dvnl_ray(1,1,2,2),nray,nx,Q)
219          DD = nwpw_splint(G_ray,dvnl_ray(1,2,2,1),
220     >                           dvnl_ray(1,2,2,2),nray,nx,Q)
221          lcount = lcount-1
222          T = (3.0d0*gz*gz-1.0d0)/(2.0d0*dsqrt(3.0d0))
223          dTdux = 0.0d0
224          dTduy = 0.0d0
225          dTduz = 6.0d0*gz/(2.0d0*dsqrt(3.0d0))
226          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
227          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
228          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
229          dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
230          dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
231          dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
232
233          lcount = lcount-1
234          T = gx*gy
235          dTdux = gy
236          dTduy = gx
237          dTduz = 0.0d0
238          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
239          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
240          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
241          dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
242          dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
243          dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
244
245          lcount = lcount-1
246          T = gy*gz
247          dTdux = 0.0d0
248          dTduy = gz
249          dTduz = gy
250          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
251          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
252          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
253          dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
254          dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
255          dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
256
257          lcount = lcount-1
258          T = gz*gx
259          dTdux = gz
260          dTduy = 0.0d0
261          dTduz = gx
262          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
263          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
264          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
265          dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
266          dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
267          dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
268
269          lcount = lcount-1
270          T = (gx*gx-gy*gy)/2.0d0
271          dTdux = gx
272          dTduy = -gy
273          dTduz = 0.0d0
274          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
275          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
276          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
277          dvnl(k1,k2,k3,1,lcount)=DD*T*gx + D*sumx
278          dvnl(k1,k2,k3,2,lcount)=DD*T*gy + D*sumy
279          dvnl(k1,k2,k3,3,lcount)=DD*T*gz + D*sumz
280        end if
281*::::::::::::::::::::::::::::::  p-wave  ::::::::::::::::::::::::::::::
282  400      CONTINUE
283        if (locp.ne.1) then
284           P  = nwpw_splint(G_ray,dvnl_ray(1,1,1,1),
285     >                            dvnl_ray(1,1,1,2),nray,nx,Q)
286           PP = nwpw_splint(G_ray,dvnl_ray(1,2,1,1),
287     >                            dvnl_ray(1,2,1,2),nray,nx,Q)
288           lcount = lcount-1
289           T = gx
290           dTdux = 1.0d0
291           dTduy = 0.0d0
292           dTduz = 0.0d0
293           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
294           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
295           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
296           dvnl(k1,k2,k3,1,lcount)= PP*T*gx + P*sumx
297           dvnl(k1,k2,k3,2,lcount)= PP*T*gy + P*sumy
298           dvnl(k1,k2,k3,3,lcount)= PP*T*gz + P*sumz
299
300
301           lcount = lcount-1
302           T = gy
303           dTdux = 0.0d0
304           dTduy = 1.0d0
305           dTduz = 0.0d0
306           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
307           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
308           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
309           dvnl(k1,k2,k3,1,lcount)= PP*T*gx + P*sumx
310           dvnl(k1,k2,k3,2,lcount)= PP*T*gy + P*sumy
311           dvnl(k1,k2,k3,3,lcount)= PP*T*gz + P*sumz
312
313           lcount = lcount-1
314           T = gz
315           dTdux = 0.0d0
316           dTduy = 0.0d0
317           dTduz = 1.0d0
318           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
319           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
320           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
321           dvnl(k1,k2,k3,1,lcount)= PP*T*gx + P*sumx
322           dvnl(k1,k2,k3,2,lcount)= PP*T*gy + P*sumy
323           dvnl(k1,k2,k3,3,lcount)= PP*T*gz + P*sumz
324        end if
325*::::::::::::::::::::::::::::::  s-wave  :::::::::::::::::::::::::::::::
326  500      CONTINUE
327        if (locp.ne.0) then
328          P  = nwpw_splint(G_ray,dvnl_ray(1,1,0,1),
329     >                           dvnl_ray(1,1,0,2),nray,nx,Q)
330          lcount = lcount-1
331          dvnl(k1,k2,k3,1,lcount) = P *gx
332          dvnl(k1,k2,k3,2,lcount) = P *gy
333          dvnl(k1,k2,k3,3,lcount) = P *gz
334        end if
335
336
337  600      CONTINUE
338*:::::::::::::::::::::::::::::::  G+k=0  ::::::::::::::::::::::::::::::::
339      else
340
341        do l=1,lmmax
342          dvnl(1,1,1,1,l)=0.0d0
343          dvnl(1,1,1,2,l)=0.0d0
344          dvnl(1,1,1,3,l)=0.0d0
345        end do
346
347       end if
348
349  700 CONTINUE
350      call Parallel_Vector_Sumall(3*lmmax*nfft3d,dvnl)
351
352      ierr=0
353      RETURN
354      END
355
356
357
358