1*
2* $Id$
3*
4
5*     **************************************
6*     *                                    *
7*     *  integrate_kbpp_e_band_stress_local  *
8*     *                                    *
9*     **************************************
10
11      subroutine integrate_kbpp_e_band_stress_local(version,
12     >                            nrho,drho,lmax,locp,nmax,
13     >                            n_extra,n_expansion,zv,
14     >                            vp,wp,rho,f,cs,sn,
15     >                            nfft1,nfft2,nfft3,nprj,
16     >                            G,dvl,
17     >                            semicore,rho_sc_r,rho_sc_k,
18     >                            ierr)
19      implicit none
20      integer          version
21      integer          nrho
22      double precision drho
23      integer          lmax
24      integer          locp
25      integer          nmax
26      integer          n_extra,n_expansion(0:lmax)
27      double precision zv
28      double precision vp(nrho,0:lmax)
29      double precision wp(nrho,0:(lmax+n_extra))
30      double precision rho(nrho)
31      double precision f(nrho)
32      double precision cs(nrho)
33      double precision sn(nrho)
34
35      integer nfft1,nfft2,nfft3,nprj
36      double precision G(nfft1,nfft2,nfft3,3)
37      double precision dvl(nfft1,nfft2,nfft3)
38
39      logical semicore
40      double precision rho_sc_r(nrho,2)
41      double precision rho_sc_k(nfft1,nfft2,nfft3,4)
42
43      integer ierr
44
45      integer np,taskid,MASTER
46      parameter (MASTER=0)
47
48*     *** local variables ****
49      integer task_count,nfft3d
50      integer k1,k2,k3,i,l,n,nb
51      double precision pi,twopi,forpi
52      double precision q
53      !integer indx(5,0:3)
54
55*     **** external functions ****
56      double precision simp
57      external         simp
58
59      call Parallel_np(np)
60      call Parallel_taskid(taskid)
61
62      nfft3d = (nfft1)*nfft2*nfft3
63      pi=4.0d0*datan(1.0d0)
64      twopi=2.0d0*pi
65      forpi=4.0d0*pi
66
67      if ((nrho/2)*2.eq.nrho) then
68        ierr=2
69        return
70      end if
71
72
73*::::::::::::::::::  Define non-local pseudopotential  ::::::::::::::::
74      do l=0,lmax
75        if (l.ne.locp) then
76          do I=1,nrho
77            vp(i,l)=vp(i,l)-vp(i,locp)
78          end do
79        end if
80      end do
81
82*======================  Fourier transformation  ======================
83      call dcopy(nfft3d,0.0d0,0,dvl,1)
84      call dcopy(4*nfft3d,0.0d0,0,rho_sc_k,1)
85      task_count = -1
86      DO 700 k3=1,nfft3
87      DO 700 k2=1,nfft2
88      DO 700 k1=1,nfft1
89        task_count = task_count + 1
90        if (mod(task_count,np).ne.taskid) go to 700
91
92        q=dsqrt(G(k1,k2,k3,1)**2
93     >         +G(k1,k2,k3,2)**2
94     >         +G(k1,k2,k3,3)**2)
95
96        if ((k1.eq.1).and.(k2.eq.1).and.(k3.eq.1)) go to 700
97
98        DO I=1,nrho
99          cs(I)=dcos(q*rho(I))
100          sn(I)=dsin(q*rho(I))
101        END DO
102
103*::::::::::::::::::::::::::::::  local  :::::::::::::::::::::::::::::::
104  600   CONTINUE
105
106        do  i=1,nrho
107          f(i)=rho(i)*vp(i,locp)*(rho(i)*cs(i)-sn(i)/q)
108        end do
109        dvl(k1,k2,k3)= simp(nrho,f,drho)*forpi/q
110     >   + zv*forpi/(q*q)*(2.0d0*cs(nrho)/q + rho(nrho)*sn(nrho))
111
112
113
114*::::::::::::::::::::: semicore density :::::::::::::::::::::::::::::::
115        if (semicore) then
116
117           do  i=1,nrho
118             f(i)=rho(i)*dsqrt(rho_sc_r(i,1))*(rho(i)*cs(i)-sn(i)/q)
119           end do
120           rho_sc_k(k1,k2,k3,1)= simp(nrho,f,drho)*forpi/q
121        end if
122
123  700 CONTINUE
124
125      call Parallel_Vector_SumAll(4*nfft3d,rho_sc_k)
126      call Parallel_Vector_SumAll(nfft3d,dvl)
127*:::::::::::::::::::::::::::::::  G=0  ::::::::::::::::::::::::::::::::
128
129      dvl(1,1,1)= 0.0d0
130
131      ierr=0
132      return
133      end
134
135