1*
2* $Id$
3*
4
5*     *******************************
6*     *			  	    *
7*     *	        cloak_init          *
8*     *				    *
9*     *******************************
10      subroutine cloak_init()
11      implicit none
12
13#include "bafdecls.fh"
14#include "errquit.fh"
15#include "cloak_common.fh"
16
17      logical value
18      integer nfft3d
19
20*     **** allocate masker memory ****
21      call C3dB_nfft3d(1,nfft3d)
22
23      value = BA_alloc_get(mt_log,nfft3d,
24     >                     'masker',masker(2),masker(1))
25      if (.not. value) call errquit('out of heap memory',0, MA_ERR)
26
27      return
28      end
29
30*     *******************************
31*     *				    *
32*     *	        cloak_end           *
33*     *				    *
34*     *******************************
35      subroutine cloak_end()
36      implicit none
37#include "errquit.fh"
38
39#include "bafdecls.fh"
40#include "cloak_common.fh"
41
42      logical value
43
44      value = BA_free_heap(masker(2))
45      if (.not. value) call errquit('out of heap memory',0, MA_ERR)
46
47      return
48      end
49
50*     *******************************
51*     *				    *
52*     *	        cloak_set           *
53*     *				    *
54*     *******************************
55      subroutine cloak_set(kvector,ecut)
56      implicit none
57      real*8 kvector(3),ecut
58
59#include "bafdecls.fh"
60#include "cloak_common.fh"
61
62*     **** local variables ****
63      integer nfft3d
64      integer i,j,k,p,q,index
65      integer k1,k2,k3
66      integer nx,ny,nz
67      integer nxh,nyh,nzh
68      real*8  ggcut,g1,g2,g3,gg
69      integer taskid
70
71*     **** external functions ***
72      real*8   lattice_unitg_frozen
73      external lattice_unitg_frozen
74
75      call Parallel3d_taskid_i(taskid)
76      call C3dB_nfft3d(1,nfft3d)
77      call C3dB_nx(1,nx)
78      call C3dB_ny(1,ny)
79      call C3dB_nz(1,nz)
80      nxh = nx/2
81      nyh = ny/2
82      nzh = nz/2
83
84*     **** set all masker on ****
85      do i=1,nfft3d
86            log_mb(masker(1)+i-1) = .true.
87      end do
88
89*     **** get fermi sphere cut-off ****
90      nwave = 0
91      ggcut = 2.0d0*ecut
92
93
94*     **** undo masker in sphere defined by ggcut ****
95      do k3 = -nzh+1, nzh-1
96         do k2 = -nyh+1, nyh-1
97            do k1 = -nxh+1, nxh-1
98               g1 = k1*lattice_unitg_frozen(1,1)
99     >            + k2*lattice_unitg_frozen(1,2)
100     >            + k3*lattice_unitg_frozen(1,3)
101     >            + kvector(1)
102               g2 = k1*lattice_unitg_frozen(2,1)
103     >            + k2*lattice_unitg_frozen(2,2)
104     >            + k3*lattice_unitg_frozen(2,3)
105     >            + kvector(2)
106               g3 = k1*lattice_unitg_frozen(3,1)
107     >            + k2*lattice_unitg_frozen(3,2)
108     >            + k3*lattice_unitg_frozen(3,3)
109     >            + kvector(3)
110               i=k1
111               j=k2
112               k=k3
113               if (i .lt. 0) i = i + nx
114               if (j .lt. 0) j = j + ny
115               if (k .lt. 0) k = k + nz
116
117                !call C3dB_ktoqp(1,k+1,q,p)
118                call C3dB_ijktoindexp(1,i+1,j+1,k+1,index,p)
119                if (p .eq. taskid) then
120                  gg = g1*g1 + g2*g2 + g3*g3
121                  if (gg.lt.ggcut) then
122c                    index = (q-1)*(nx)*ny
123c    >                     + j*(nx)
124c    >                     + i+1
125                     log_mb(masker(1)+index-1) = .false.
126                     nwave = nwave + 1
127                  end if
128               end if
129            end do
130         end do
131      end do
132      nwave_all = nwave
133      call C3dB_ISumAll(nwave_all)
134
135      return
136      end
137
138*     *******************************
139*     *		 		    *
140*     *	        cloak_C             *
141*     *				    *
142*     *******************************
143      subroutine cloak_C(A)
144      implicit none
145      complex*16 A(*)
146
147#include "bafdecls.fh"
148#include "cloak_common.fh"
149
150*     **** local variables ****
151      integer i,nfft3d
152
153      call nwpw_timing_start(9)
154
155      call C3dB_nfft3d(1,nfft3d)
156      do i=1,nfft3d
157         if (log_mb(masker(1)+i-1)) A(i) = dcmplx(0.0d0,0.0d0)
158      end do
159
160      call nwpw_timing_end(9)
161      return
162      end
163
164*     *******************************
165*     *				    *
166*     *	        cloak_R             *
167*     *				    *
168*     *******************************
169      subroutine cloak_R(A)
170      implicit none
171      real*8    A(*)
172
173#include "bafdecls.fh"
174#include "cloak_common.fh"
175
176*     **** local variables ****
177      integer i,nfft3d
178
179      call C3dB_nfft3d(1,nfft3d)
180      do i=1,nfft3d
181         if (log_mb(masker(1)+i-1)) A(i) = 0.0d0
182      end do
183
184      return
185      end
186
187
188*     *******************************
189*     *			  	    *
190*     *	        cloak_nwave         *
191*     *				    *
192*     *******************************
193      integer function cloak_nwave()
194      implicit none
195
196#include "cloak_common.fh"
197
198      cloak_nwave = nwave
199      return
200      end
201
202*     *******************************
203*     *				    *
204*     *	        cloak_nwave_all     *
205*     *				    *
206*     *******************************
207      integer function cloak_nwave_all()
208      implicit none
209
210#include "cloak_common.fh"
211
212      cloak_nwave_all = nwave_all
213      return
214      end
215
216*     *******************************
217*     *				    *
218*     *	        cloak_masker        *
219*     *				    *
220*     *******************************
221      logical function cloak_masker(i)
222      implicit none
223      integer i
224
225#include "bafdecls.fh"
226#include "cloak_common.fh"
227
228      cloak_masker = log_mb(masker(1)+i-1)
229      return
230      end
231
232