1*
2* $Id$
3*
4
5*     **********************************
6*     *                                *
7*     *        integrate_d_stress      *
8*     *                                *
9*     **********************************
10
11      subroutine integrate_d_stress(version,rlocal,
12     >                            nrho,drho,lmax,locp,zv,
13     >                            vp,wp,rho,f,cs,sn,
14     >                            nfft3d,lmmax,
15     >                            G,dvl,dvnl,
16     >                            semicore,rho_sc_r,rho_sc_k,
17     >                            ierr)
18      implicit none
19      integer          version
20      double precision rlocal
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 nfft3d,lmmax
34      double precision G(nfft3d,3)
35      double precision dvl(nfft3d)
36      double precision dvnl(nfft3d,3,lmmax)
37
38      logical semicore
39      double precision rho_sc_r(nrho,2)
40      double precision rho_sc_k(nfft3d,4)
41
42      integer ierr
43
44#include "errquit.fh"
45
46      integer np,taskid,MASTER
47      integer np_i,np_j,taskid_i,taskid_j,countj
48      parameter (MASTER=0)
49
50*     *** local variables ****
51      integer lcount
52      integer k1,k2,k3,i,l,pzero,zero
53      double precision pi,twopi,forpi
54      double precision p0,p1,p2,p3,p,pp
55      double precision gx,gy,gz,a,q,d,dd
56      double precision duxdGx,duxdGy,duxdGz
57      double precision duydGx,duydGy,duydGz
58      double precision duzdGx,duzdGy,duzdGz
59      double precision sumx,sumy,sumz
60      double precision T,dTdux,dTduy,dTduz
61
62*     **** external functions ****
63      double precision dsum,simp
64      external         dsum,simp
65      double precision tollz
66      parameter(tollz=1d-16)
67
68      if (version.ne.3) then
69         call errquit('integrate_stress - unit cell is aperiodic',0,
70     &       INPUT_ERR)
71      end if
72      call Parallel2d_np_i(np_i)
73      call Parallel2d_np_j(np_j)
74      call Parallel2d_taskid_i(taskid_i)
75      call Parallel2d_taskid_j(taskid_j)
76
77      pi=4.0d0*datan(1.0d0)
78      twopi=2.0d0*pi
79      forpi=4.0d0*pi
80
81      IF(LMMAX.GT.16) THEN
82        IERR=1
83        RETURN
84      ENDIF
85      IF((nrho/2)*2.Eq.nrho) THEN
86        IERR=2
87        RETURN
88      ENDIF
89
90      P0=dsqrt(forpi)
91      P1=dsqrt(3.0d0*forpi)
92      P2=dsqrt(15.0d0*forpi)
93      P3=dsqrt(105.0d0*forpi)
94
95*::::::::::::::::::  Define non-local pseudopotential  ::::::::::::::::
96      do l=0,lmax
97        if (l.ne.locp) then
98          do I=1,nrho
99            vp(i,l)=vp(i,l)-vp(i,locp)
100          end do
101        end if
102      end do
103
104*======================  Fourier transformation  ======================
105      call dcopy(nfft3d,0.0d0,0,dvl,1)
106      call dcopy(3*lmmax*nfft3d,0.0d0,0,dvnl,1)
107      call dcopy(4*nfft3d,0.0d0,0,rho_sc_k,1)
108
109*     ***** find the G==0 point in the lattice *****
110      call D3dB_ijktoindexp(1,1,1,1,zero,pzero)
111
112      countj = -1
113      DO 700 k1=1,nfft3d
114
115        countj = mod(countj+1,np_j)
116        if (countj.ne.taskid_j) go to 700
117        if ((pzero.eq.taskid_i).and.(k1.eq.zero)) go to 700
118
119        q=DSqRT(G(k1,1)**2
120     >         +G(k1,2)**2
121     >         +G(k1,3)**2)
122        if(abs(q).lt.tollz) goto 700
123
124        gx=G(k1,1)/q
125        gy=G(k1,2)/q
126        gz=G(k1,3)/q
127        DO i=1,nrho
128          cs(i)=DCOS(q*rho(i))
129          sn(i)=DSIN(q*rho(i))
130        END DO
131
132*       **** calculate du_r/dG_s ****
133        duxdGx = 1.0d0/q -gx*gx/q
134        duxdGy = -gx*gy/q
135        duxdGz = -gx*gz/q
136
137        duydGx = -gy*gx/q
138        duydGy = 1.0d0/q - gy*gy/q
139        duydGz = -gy*gz/q
140
141        duzdGx = -gz*gx/q
142        duzdGy = -gz*gy/q
143        duzdGz = 1.0d0/q - gz*gz/q
144
145        lcount = lmmax+1
146        GO TO (500,400,300,200), LMAX+1
147
148
149*::::::::::::::::::::::::::::::  f-wave  ::::::::::::::::::::::::::::::
150  200   CONTINUE
151        if (locp.ne.3) then
152           F(1)=0.0d0
153           do i=2,nrho
154             A=sn(i)/(q*rho(i))
155             A=15.0d0*(A-cs(i))/(q*rho(i))**2 - 6*A + cs(i)
156             f(i)=A*wp(i,3)*vp(i,3)
157           end do
158           D=P3*SIMP(nrho,F,drho)/q
159
160           F(1)=0.0d0
161           do i=2,nrho
162             A= -60.0d0*sn(i)/(rho(i)**3 * q**5)
163     >        +  60.0d0*cs(i)/(rho(i)**2 * q**4)
164     >        +  27.0d0*sn(i)/(rho(i)    * q**3)
165     >        -   7.0d0*cs(i)/(q**2)
166     >        -   rho(i)*sn(i)/q
167             f(i)=A*wp(i,3)*vp(i,3)
168           end do
169           DD=P3*SIMP(nrho,F,drho)
170
171           lcount = lcount-1
172           T = gy*(3.0d0*(1.0d0-gz*gz)-4.0d0*gy*gy)/dsqrt(24.0d0)
173           dTdux = 0.0d0
174           dTduy = (3.0d0*(1.0d0-gz*gz)-12.0d0*gy*gy)/dsqrt(24.0d0)
175           dTduz = -6.0d0*gy*gz/dsqrt(24.0d0)
176           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
177           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
178           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
179           dvnl(k1,1,lcount)=DD*T*gx + D*sumx
180           dvnl(k1,2,lcount)=DD*T*gy + D*sumy
181           dvnl(k1,3,lcount)=DD*T*gz + D*sumz
182
183           lcount = lcount-1
184           T =gx*gy*gz
185           dTdux = gy*gz
186           dTduy = gx*gz
187           dTduz = gx*gy
188           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
189           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
190           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
191           dvnl(k1,1,lcount)=DD*T*gx + D*sumx
192           dvnl(k1,2,lcount)=DD*T*gy + D*sumy
193           dvnl(k1,3,lcount)=DD*T*gz + D*sumz
194
195           lcount = lcount-1
196           T = gy*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
197           dTdux = 0.0d0
198           dTduy =(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
199           dTduz =10.0d0*gy*gz/dsqrt(40.0d0)
200           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
201           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
202           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
203           dvnl(k1,1,lcount)=DD*T*gx + D*sumx
204           dvnl(k1,2,lcount)=DD*T*gy + D*sumy
205           dvnl(k1,3,lcount)=DD*T*gz + D*sumz
206
207           lcount = lcount-1
208           T =gz*(5.0d0*gz*gz-3.0d0)/dsqrt(60.0d0)
209           dTdux = 0.0d0
210           dTduy = 0.0d0
211           dTduz =(15.0d0*gz*gz -3.0d0)/dsqrt(60.0d0)
212           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
213           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
214           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
215           dvnl(k1,1,lcount)=DD*T*gx + D*sumx
216           dvnl(k1,2,lcount)=DD*T*gy + D*sumy
217           dvnl(k1,3,lcount)=DD*T*gz + D*sumz
218
219           lcount = lcount-1
220           T =  gx*(5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
221           dTdux = (5.0d0*gz*gz-1.0d0)/dsqrt(40.0d0)
222           dTduy = 0.0d0
223           dTduz = 10.0d0*gx*gz/dsqrt(40.0d0)
224           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
225           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
226           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
227           dvnl(k1,1,lcount)=DD*T*gx + D*sumx
228           dvnl(k1,2,lcount)=DD*T*gy + D*sumy
229           dvnl(k1,3,lcount)=DD*T*gz + D*sumz
230
231           lcount = lcount-1
232           T =gz*(gx*gx - gy*gy)/2.0d0
233           dTdux =  gx*gz
234           dTduy = -gy*gz
235           dTduz = (gx*gx-gy*gy)/2.0d0
236           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
237           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
238           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
239           dvnl(k1,1,lcount)=DD*T*gx + D*sumx
240           dvnl(k1,2,lcount)=DD*T*gy + D*sumy
241           dvnl(k1,3,lcount)=DD*T*gz + D*sumz
242
243           lcount = lcount-1
244           T = gx*(4.0d0*gx*gx - 3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0)
245           dTdux = (12.0d0*gx*gx-3.0d0*(1.0d0-gz*gz))/dsqrt(24.0d0)
246           dTduy = 0.0d0
247           dTduz = 6.0d0*gx*gz/dsqrt(24.0d0)
248           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
249           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
250           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
251           dvnl(k1,1,lcount)=DD*T*gx + D*sumx
252           dvnl(k1,2,lcount)=DD*T*gy + D*sumy
253           dvnl(k1,3,lcount)=DD*T*gz + D*sumz
254        end if
255
256
257
258*::::::::::::::::::::::::::::::  d-wave  ::::::::::::::::::::::::::::::
259  300   CONTINUE
260        if (locp.ne.2) then
261          F(1)=0.0d0
262          DO i=2,nrho
263            A=3.0d0*(sn(i)/(q*rho(i))-cs(i))/(q*rho(i))-sn(i)
264            f(i)=A*wp(i,2)*vp(i,2)
265          END DO
266          D=P2*SIMP(nrho,F,drho)/q
267
268          F(1)=0.0d0
269          DO i=2,nrho
270            A= -9.0d0*sn(i)/(rho(i)**2 * q**4)
271     >       +  9.0d0*cs(i)/(rho(i)    * q**3)
272     >       +  4.0d0*sn(i)/(q**2)
273     >       -  rho(i)*cs(i)/q
274            f(i)=A*wp(i,2)*vp(i,2)
275          END DO
276          DD=P2*SIMP(nrho,F,drho)
277
278          lcount = lcount-1
279          T = gx*gy
280          dTdux = gy
281          dTduy = gx
282          dTduz = 0.0d0
283          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
284          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
285          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
286          dvnl(k1,1,lcount)=DD*T*gx + D*sumx
287          dvnl(k1,2,lcount)=DD*T*gy + D*sumy
288          dvnl(k1,3,lcount)=DD*T*gz + D*sumz
289
290          lcount = lcount-1
291          T = gy*gz
292          dTdux = 0.0d0
293          dTduy = gz
294          dTduz = gy
295          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
296          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
297          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
298          dvnl(k1,1,lcount)=DD*T*gx + D*sumx
299          dvnl(k1,2,lcount)=DD*T*gy + D*sumy
300          dvnl(k1,3,lcount)=DD*T*gz + D*sumz
301
302          lcount = lcount-1
303          T = (3.0d0*gz*gz-1.0d0)/(2.0d0*dsqrt(3.0d0))
304          dTdux = 0.0d0
305          dTduy = 0.0d0
306          dTduz = 6.0d0*gz/(2.0d0*dsqrt(3.0d0))
307          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
308          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
309          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
310          dvnl(k1,1,lcount)=DD*T*gx + D*sumx
311          dvnl(k1,2,lcount)=DD*T*gy + D*sumy
312          dvnl(k1,3,lcount)=DD*T*gz + D*sumz
313
314          lcount = lcount-1
315          T = gz*gx
316          dTdux = gz
317          dTduy = 0.0d0
318          dTduz = gx
319          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
320          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
321          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
322          dvnl(k1,1,lcount)=DD*T*gx + D*sumx
323          dvnl(k1,2,lcount)=DD*T*gy + D*sumy
324          dvnl(k1,3,lcount)=DD*T*gz + D*sumz
325
326          lcount = lcount-1
327          T = (gx*gx-gy*gy)/2.0d0
328          dTdux = gx
329          dTduy = -gy
330          dTduz = 0.0d0
331          sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
332          sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
333          sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
334          dvnl(k1,1,lcount)=DD*T*gx + D*sumx
335          dvnl(k1,2,lcount)=DD*T*gy + D*sumy
336          dvnl(k1,3,lcount)=DD*T*gz + D*sumz
337        end if
338
339*::::::::::::::::::::::::::::::  p-wave  ::::::::::::::::::::::::::::::
340  400   CONTINUE
341        if (locp.ne.1) then
342           F(1)=0.0d0
343           DO i=2,nrho
344             f(i)=(sn(i)/(q*rho(i)) - cs(i)) * wp(i,1)*vp(i,1)
345           END DO
346           P=P1*SIMP(nrho,F,drho)/q
347
348           F(1)=0.0d0
349           DO i=2,nrho
350             f(i)=wp(i,1)*vp(i,1)* ( -2.0d0*sn(i)/(rho(i) * q**3)
351     >                              + 2.0d0*cs(i)/(q**2)
352     >                              + rho(i)*sn(i)/q)
353           END DO
354           PP=P1*SIMP(nrho,F,drho)
355
356           lcount = lcount-1
357           T = gy
358           dTdux = 0.0d0
359           dTduy = 1.0d0
360           dTduz = 0.0d0
361           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
362           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
363           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
364           dvnl(k1,1,lcount)= PP*T*gx + P*sumx
365           dvnl(k1,2,lcount)= PP*T*gy + P*sumy
366           dvnl(k1,3,lcount)= PP*T*gz + P*sumz
367
368           lcount = lcount-1
369           T = gz
370           dTdux = 0.0d0
371           dTduy = 0.0d0
372           dTduz = 1.0d0
373           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
374           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
375           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
376           dvnl(k1,1,lcount)= PP*T*gx + P*sumx
377           dvnl(k1,2,lcount)= PP*T*gy + P*sumy
378           dvnl(k1,3,lcount)= PP*T*gz + P*sumz
379
380           lcount = lcount-1
381           T = gx
382           dTdux = 1.0d0
383           dTduy = 0.0d0
384           dTduz = 0.0d0
385           sumx = dTdux*duxdGx + dTduy*duydGx + dTduz*duzdGx
386           sumy = dTdux*duxdGy + dTduy*duydGy + dTduz*duzdGy
387           sumz = dTdux*duxdGz + dTduy*duydGz + dTduz*duzdGz
388           dvnl(k1,1,lcount)= PP*T*gx + P*sumx
389           dvnl(k1,2,lcount)= PP*T*gy + P*sumy
390           dvnl(k1,3,lcount)= PP*T*gz + P*sumz
391
392        end if
393
394*::::::::::::::::::::::::::::::  s-wave  :::::::::::::::::::::::::::::::
395  500   CONTINUE
396        if (locp.ne.0) then
397          DO i=1,nrho
398            f(i)=wp(i,0)*vp(i,0) * ( -sn(i)/(q**2)
399     >                              + rho(i)*cs(i)/q)
400          END DO
401          P = P0*SIMP(nrho,F,drho)
402          lcount = lcount-1
403          dvnl(k1,1,lcount) = P *gx
404          dvnl(k1,2,lcount) = P *gy
405          dvnl(k1,3,lcount) = P *gz
406        end if
407
408*::::::::::::::::::::::::::::::  local  :::::::::::::::::::::::::::::::
409  600   CONTINUE
410
411        do  i=1,nrho
412          f(i)=rho(i)*vp(i,locp)*(rho(i)*cs(i)-sn(i)/q)
413        end do
414        dvl(k1)= SIMP(nrho,f,drho)*forpi/q
415     >   + zv*forpi/(q*q)*(2.0d0*cs(nrho)/q + rho(nrho)*sn(nrho))
416
417
418
419*::::::::::::::::::::: semicore density :::::::::::::::::::::::::::::::
420        if (semicore) then
421
422           do  i=1,nrho
423             f(i)=rho(i)*dsqrt(rho_sc_r(i,1))*(rho(i)*cs(i)-sn(i)/q)
424           end do
425           rho_sc_k(k1,1)= SIMP(nrho,f,drho)*forpi/q
426        end if
427
428  700 CONTINUE
429      call D1dB_Vector_SumAll(nfft3d,dvl)
430      call D1dB_Vector_Sumall(3*lmmax*nfft3d,dvnl)
431      call D1dB_Vector_SumAll(nfft3d,rho_sc_k)
432
433
434
435*:::::::::::::::::::::::::::::::  G=0  ::::::::::::::::::::::::::::::::
436      if (pzero.eq.taskid_i) then
437      dvl(zero)= 0.0d0
438      do l=1,lmmax
439        dvnl(zero,1,l)=0.0d0
440        dvnl(zero,2,l)=0.0d0
441        dvnl(zero,3,l)=0.0d0
442      end do
443      end if
444
445
446      IERR=0
447      RETURN
448      END
449
450
451
452