1*
2* $Id$
3*
4
5*     *********************************************
6*     *                                           *
7*     *               kbpp_calc_Gmax              *
8*     *                                           *
9*     *********************************************
10
11      double precision function kbpp_calc_Gmax(nfft1,nfft2,nfft3,unita)
12      implicit none
13      integer nfft1,nfft2,nfft3
14      double precision unita(3,3)
15
16*     *** local variables ****
17      integer k1,k2,k3
18      double precision qmax,gx,gy,gz,gg,ggmax
19      double precision unitg(3,3)
20
21      call get_unitg(unita,unitg)
22      k1 = nfft1/2
23      k2 = nfft2/2
24      k3 = nfft3/2
25      ggmax = 0.0d0
26      gx = unitg(1,1)*k1
27      gy = unitg(2,1)*k1
28      gz = unitg(3,1)*k1
29      gg = gx*gx + gy*gy + gz*gz
30      if (gg.gt.ggmax) ggmax = gg
31      gx = unitg(1,2)*k2
32      gy = unitg(2,2)*k2
33      gz = unitg(3,2)*k2
34      gg = gx*gx + gy*gy + gz*gz
35      if (gg.gt.ggmax) ggmax = gg
36      gx = unitg(1,1)*k3
37      gy = unitg(2,1)*k3
38      gz = unitg(3,1)*k3
39      gg = gx*gx + gy*gy + gz*gz
40      if (gg.gt.ggmax) ggmax = gg
41
42      qmax = dsqrt(ggmax)
43
44      qmax = qmax
45
46      kbpp_calc_Gmax = qmax
47      return
48      end
49
50*     *********************************************
51*     *                                           *
52*     *               kbpp_calc_dGmin             *
53*     *                                           *
54*     *********************************************
55
56      double precision function kbpp_calc_dGmin(unita)
57      implicit none
58      double precision unita(3,3)
59
60*     *** local variables ****
61      double precision gx,gy,gz,q,dqmin
62      double precision unitg(3,3)
63
64      call get_unitg(unita,unitg)
65
66      gx = unitg(1,1)
67      gy = unitg(2,1)
68      gz = unitg(3,1)
69      q = dsqrt(gx**2 + gy**2 + gz**2)
70      dqmin = q
71
72      gx = unitg(1,2)
73      gy = unitg(2,2)
74      gz = unitg(3,2)
75      q = dsqrt(gx**2 + gy**2 + gz**2)
76      if (q.lt.dqmin) dqmin = q
77
78      gx = unitg(1,3)
79      gy = unitg(2,3)
80      gz = unitg(3,3)
81      q = dsqrt(gx**2 + gy**2 + gz**2)
82      if (q.lt.dqmin) dqmin = q
83
84*     *** make fine dqmin ****
85      dqmin = 0.01d0*dqmin
86
87      kbpp_calc_dGmin = dqmin
88      return
89      end
90
91*     *********************************************
92*     *                                           *
93*     *               kbpp_calc_nray              *
94*     *                                           *
95*     *********************************************
96      integer function kbpp_calc_nray(nfft1,nfft2,nfft3,unita)
97      implicit none
98      integer nfft1,nfft2,nfft3
99      double precision unita(3,3)
100
101*     ***** local variables ****
102      double precision dG,Gmax
103      integer nray
104
105*     ***** external functions ****
106      real*8   kbpp_calc_dGmin,kbpp_calc_Gmax
107      external kbpp_calc_dGmin,kbpp_calc_Gmax
108
109      dG   = kbpp_calc_dGmin(unita)
110      Gmax = kbpp_calc_Gmax(nfft1,nfft2,nfft3,unita) + 2.0d0
111      nray = Gmax/dG + 1.0d0
112      if (nray.lt.10) nray = 10
113
114      kbpp_calc_nray = nray
115      return
116      end
117
118*     *********************************************
119*     *                                           *
120*     *               kbpp_generate_G_ray         *
121*     *                                           *
122*     *********************************************
123      subroutine kbpp_generate_G_ray(nfft1,nfft2,nfft3,unita,G_ray)
124      implicit none
125      integer nfft1,nfft2,nfft3
126      double precision unita(3,3)
127      double precision G_ray(*)
128
129*     **** local variables ***
130      integer i,nray
131      double precision dG
132*     ***** external functions ****
133      real*8   kbpp_calc_dGmin
134      integer  kbpp_calc_nray
135      external kbpp_calc_dGmin
136      external kbpp_calc_nray
137
138      dG   = kbpp_calc_dGmin(unita)
139      nray = kbpp_calc_nray(nfft1,nfft2,nfft3,unita)
140      do i=1,nray
141       G_ray(i) = dG*dble(i-1)
142      end do
143      return
144      end
145
146
147*     *********************************************
148*     *                                           *
149*     *               kbpp_filter_ray             *
150*     *                                           *
151*     *********************************************
152      subroutine kbpp_filter_ray(nray,G_ray,ecut,V_ray)
153      implicit none
154      integer nray
155      double precision G_ray(*)
156      double precision ecut
157      double precision V_ray(*)
158
159*     **** local variables ****
160      !real*8 ncut,eps
161      !parameter (eps=1.0d-12,ncut=15.0d0)
162      integer ncut
163      parameter (ncut=15)
164      integer i
165      double precision g,qmax,fac
166
167c      qmax = dsqrt(ecut+ecut)
168c     >      /(-log(1.0d0-(1.0d0-eps)**(1.0d0/ncut)))**(1.0d0/ncut)
169      qmax = dsqrt(ecut+ecut)
170      do i=1,nray
171         g = G_ray(i)
172         if (g.gt.(qmax-0.2d0)) then
173            fac = 1.0d0 - (1.0d0-exp(-(g/qmax)**ncut))**ncut
174            V_ray(i) = V_ray(i)*fac
175         end if
176      end do
177      return
178      end
179
180*     *********************************************
181*     *                                           *
182*     *               kbpp_filter                 *
183*     *                                           *
184*     *********************************************
185      subroutine kbpp_filter(nfft1,nfft2,nfft3,G,ecut,vl)
186      implicit none
187      integer nfft1,nfft2,nfft3
188      double precision G(nfft1/2+1,nfft2,nfft3,3)
189      double precision ecut
190      double precision vl(nfft1/2+1,nfft2,nfft3)
191
192*     **** local variables ****
193      !real*8 ncut,eps
194      !parameter (eps=1.0d-12,ncut=15.0d0)
195      integer ncut
196      parameter (ncut=15)
197      integer i,j,k
198      double precision q,qmax,fac
199
200c      qmax = dsqrt(ecut+ecut)
201c     >      /(-log(1.0d0-(1.0d0-eps)**(1.0d0/ncut)))**(1.0d0/ncut)
202      qmax = dsqrt(ecut+ecut)
203      do k=1,nfft3
204      do j=1,nfft2
205      do i=1,(nfft1/2+1)
206         q=DSQRT(G(i,j,k,1)**2
207     >          +G(i,j,k,2)**2
208     >          +G(i,j,k,3)**2)
209
210         if (q.gt.(qmax-0.2d0)) then
211            fac = 1.0d0 - (1.0d0-exp(-(q/qmax)**ncut))**ncut
212            vl(i,j,k) = vl(i,j,k)*fac
213         end if
214      end do
215      end do
216      end do
217      return
218      end
219