1*
2* $Id$
3*
4
5
6      subroutine HGH_local(version,rlocal,
7     >                     zv,rloc,C1,C2,C3,C4,
8     >                     nfft1,nfft2,nfft3,
9     >                     G,
10     >                     vl)
11      implicit none
12      integer          version
13      double precision rlocal
14      double precision zv,rloc,C1,C2,C3,C4
15
16      integer nfft1,nfft2,nfft3
17      double precision G(nfft1/2+1,nfft2,nfft3,3)
18      double precision vl(nfft1/2+1,nfft2,nfft3)
19
20
21      integer np,taskid,MASTER
22      parameter (MASTER=0)
23
24*     *** local variables ****
25      integer task_count,nfft3d
26      integer k1,k2,k3
27      double precision pi,twopi,forpi
28      double precision b0,b1
29      double precision y,x,x2,x4,x6,a,a1,f,Q
30
31*     **** Error function parameters ****
32
33
34      call Parallel_np(np)
35      call Parallel_taskid(taskid)
36
37      nfft3d = (nfft1/2+1)*nfft2*nfft3
38      pi=4.0d0*datan(1.0d0)
39      twopi=2.0d0*pi
40      forpi=4.0d0*pi
41      b0 = -forpi*zv
42      b1 = dsqrt(8.0d0*pi*pi*pi)*rloc*rloc*rloc
43
44
45
46*======================  Fourier transformation  ======================
47      call dcopy(nfft3d,0.0d0,0,vl,1)
48      task_count = -1
49      DO 700 k3=1,nfft3
50      DO 700 k2=1,nfft2
51      DO 700 k1=1,(nfft1/2+1)
52        task_count = task_count + 1
53        if (mod(task_count,np).ne.taskid) go to 700
54
55        Q=DSQRT(G(k1,k2,k3,1)**2
56     >         +G(k1,k2,k3,2)**2
57     >         +G(k1,k2,k3,3)**2)
58
59        if ((k1.eq.1).and.(k2.eq.1).and.(k3.eq.1)) go to 700
60
61        x = Q*rloc
62        x2 = x*x
63        x4 = x2*x2
64        x6 = x4*x2
65        a = dexp(-0.5d0*x2)
66
67        f = b0*a/(Q*Q)
68     >    + b1*a
69     >        *( C1
70     >         + C2*(3.0d0-x2)
71     >         + C3*(15.0d0 - 10.0d0*x2 + x4)
72     >         + C4*(105.0d0 - 105.0d0*x2 + 21.0d0*x4 - x6))
73
74
75        if (version.eq.4) then
76          y = Q*rlocal
77          a1 = dexp(-0.25d0*y*y)
78          f = f - b0*a1/(Q*Q)
79        end if
80
81        vl(k1,k2,k3)= f
82
83  700 continue
84
85      call Parallel_Vector_SumAll(nfft3d,vl)
86
87      f = b1 *(C1+3.0d0*C2+15.0d0*C3+105.0d0*C4)-b0*0.5d0*rloc*rloc
88      if (version.eq.4) then
89        f = f + b0*(0.25d0*rlocal*rlocal)
90      end if
91      vl(1,1,1) = f
92
93
94
95      return
96      end
97
98
99
100