1*
2* $Id$
3*
4
5*     ****************************************
6*     *                                      *
7*     *              wgc_init                *
8*     *                                      *
9*     ****************************************
10      subroutine wgc_init(rho0)
11      implicit none
12      real*8 rho0
13
14#include "errquit.fh"
15#include "bafdecls.fh"
16#include "wgc.fh"
17
18
19*     **** local variables ****
20      real*8 one3rd
21      parameter (one3rd=1.0d0/3.0d0)
22      double precision toll
23      parameter (toll=1.0d-16)
24
25      integer npack0,nfft3d,G(3)
26      integer i,j,k
27      integer zero,qzero,pzero,taskid
28      integer nx,ny,nz
29      real*8  pi,gg,ss,aa1,aa2,tf,kf,lind,eta1,eta2,eta,vw
30      logical value
31      integer tmp1(2)
32
33*     **** external functions ****
34      integer  G_indx
35      external G_indx
36      real*8   control_wgc_alphabetalambda
37      external control_wgc_alphabetalambda
38
39
40      call nwpw_timing_start(7)
41      call Parallel2d_taskid_i(taskid)
42
43      call D3dB_nfft3d(1,nfft3d)
44      call Pack_npack(0,npack0)
45      G(1) = G_indx(1)
46      G(2) = G_indx(2)
47      G(3) = G_indx(3)
48
49      call D3dB_nx(1,nx)
50      call D3dB_ny(1,ny)
51      call D3dB_nz(1,nz)
52      wgc_scal1 = 1.0d0/(nx*ny*nz)
53
54      wgc_rho0    = rho0
55      wgc_tune1   = control_wgc_alphabetalambda(1)
56      wgc_tune2   = control_wgc_alphabetalambda(2)
57      wgc_lmd     = control_wgc_alphabetalambda(3)
58      wgc_eq_tune = (dabs(wgc_tune1-wgc_tune2).lt.1.0d-6)
59
60      !*******************************************************************************
61      !**** An extra scaling of 1/(N1*N2*N3) is done to the WGG_kernel(G) so that ****
62      !**** the rho(G) does not need to be scaled by 1/(N1*N2*N3)                 ****
63      !**** Might want to scal rho(G) properly to make the code more transparent  ****
64      !*******************************************************************************
65      ss = wgc_lmd*wgc_scal1/(2.0d0*wgc_tune1*wgc_tune2
66     >                       *rho0**(wgc_tune1+wgc_tune2-2.0d0))
67
68*     **** allocate vc memory ****
69      value = BA_alloc_get(mt_dbl,npack0,'wgc',wgc_hndl,wgc_indx)
70      if (.not.value) call errquit('wgc_init:out of heap',0,MA_ERR)
71
72      value = BA_push_get(mt_dbl,nfft3d,'tmp1',tmp1(2),tmp1(1))
73      if (.not.value) call errquit('wgc_init:out of stack',1,MA_ERR)
74
75
76*     ***** find the G==0 point in the lattice *****
77      i=0
78      j=0
79      k=0
80      call D3dB_ijktoindexp(1,i+1,j+1,k+1,zero,pzero)
81
82
83*     ***** form Vc = 4*pi/G**2  *****
84      pi = (4.0d0*datan(1.0d0))
85      kf = (3.0d0*pi*pi*wgc_rho0)**one3rd
86
87      do i = 1,nfft3d
88
89         gg  = ( dbl_mb(G(1)+i-1)*dbl_mb(G(1)+i-1)
90     >         + dbl_mb(G(2)+i-1)*dbl_mb(G(2)+i-1)
91     >         + dbl_mb(G(3)+i-1)*dbl_mb(G(3)+i-1) )
92
93         if (((pzero.eq.taskid) .and. (i.eq.zero)).or.
94     E    (abs(gg) .lt.toll))  then
95            dbl_mb(tmp1(1)+i-1) = 0.0d0
96         else
97            eta = dsqrt(gg)/(2.0d0*kf)
98            eta2 = eta*eta
99            tf = 3.0d0*rho0/(kf*kf)
100            vw = 4.0d0*rho0/gg
101            aa1 = (1.0d0-eta2)/(4.0d0*eta)
102            aa2 = dlog(dabs((1.0d0+eta)/(1.0d0-eta)))
103            lind = tf*(0.5d0 + aa1*aa2)
104            dbl_mb(tmp1(1)+i-1) = ss*(1.0d0/lind-1.0d0/vw-1.0d0/tf)
105         end if
106
107      end do
108      call Pack_t_pack(0,dbl_mb(tmp1(1)))
109      call Pack_t_Copy(0,dbl_mb(tmp1(1)),dbl_mb(wgc_indx))
110      call Pack_tt_dot(0,dbl_mb(wgc_indx),dbl_mb(wgc_indx),aa1)
111
112      value = BA_pop_stack(tmp1(2))
113      if (.not. value) call errquit('wgc_init:popping stack',3,MA_ERR)
114
115      call nwpw_timing_end(7)
116      return
117      end
118
119
120*     ****************************************
121*     *                                      *
122*     *              wgc_end                 *
123*     *                                      *
124*     ****************************************
125      subroutine wgc_end()
126      implicit none
127
128#include "bafdecls.fh"
129#include "errquit.fh"
130#include "wgc.fh"
131
132
133      if (.not.BA_free_heap(wgc_hndl))
134     >   call errquit('wgc_end:freeing heap',0,MA_ERR)
135
136      return
137      end
138
139
140*     ****************************************
141*     *                                      *
142*     *              wgc_rho                 *
143*     *                                      *
144*     ****************************************
145      real*8 function wgc_rho()
146      implicit none
147
148#include "wgc.fh"
149
150      wgc_rho = wgc_rho0
151      return
152      end
153
154*     ****************************************
155*     *                                      *
156*     *              wgc_alpha               *
157*     *                                      *
158*     ****************************************
159      real*8 function wgc_alpha()
160      implicit none
161
162#include "wgc.fh"
163
164      wgc_alpha = wgc_tune1
165      return
166      end
167
168*     ****************************************
169*     *                                      *
170*     *              wgc_beta                *
171*     *                                      *
172*     ****************************************
173      real*8 function wgc_beta()
174      implicit none
175
176#include "wgc.fh"
177
178      wgc_beta = wgc_tune2
179      return
180      end
181
182
183*     ****************************************
184*     *                                      *
185*     *              wgc_lambda              *
186*     *                                      *
187*     ****************************************
188      real*8 function wgc_lambda()
189      implicit none
190
191#include "wgc.fh"
192
193      wgc_lambda = wgc_lmd
194      return
195      end
196
197
198*     ****************************************
199*     *                                      *
200*     *              wgc_v                   *
201*     *                                      *
202*     ****************************************
203      subroutine wgc_v(ispin,dn,v_out)
204      implicit none
205      integer ispin
206      real*8  dn(*)
207      real*8  v_out(*)
208
209#include "bafdecls.fh"
210#include "errquit.fh"
211#include "wgc.fh"
212
213*     **** local variables ****
214      logical value
215      integer nfft3d,n2ft3d
216      integer rho1(2),tmp1(2)
217
218      call nwpw_timing_start(7)
219      call D3dB_nfft3d(1,nfft3d)
220      call D3dB_n2ft3d(1,n2ft3d)
221
222      value =      BA_push_get(mt_dbl,2*n2ft3d,'rho1',rho1(2),rho1(1))
223     >        .and.BA_push_get(mt_dbl,2*n2ft3d,'tmp1',tmp1(2),tmp1(1))
224      if (.not.value) call errquit('wgc_v:out of stack',0,MA_ERR)
225
226      !**** alpha==beta ****
227      if (wgc_eq_tune) then
228
229         !**** the rho^alpha(G) is not beinge scaled by 1/(N1*N2*N3) ****
230         call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dbl_mb(rho1(1)))
231         call D3dB_r_notZero_Ends(1,dbl_mb(rho1(1)))
232         call D3dB_r_Power1(1,wgc_tune1,dbl_mb(rho1(1)))
233         call D3dB_r_Zero_Ends(1,dbl_mb(rho1(1)))
234         !call D3dB_rc_fft3f(1,dbl_mb(rho1(1)))
235         call D3dB_rc_pfft3f(1,0,dbl_mb(rho1(1)))
236         call Pack_c_pack(0,dbl_mb(rho1(1)))
237
238         call Pack_tc_Mul(0,dbl_mb(wgc_indx),dbl_mb(rho1(1)),v_out)
239         call Pack_c_SMul1(0,2.0d0*wgc_tune1,v_out)
240         !call Pack_c_SMul1(0,wgc_tune1,v_out)
241         call Pack_c_unpack(0,v_out)
242         !call D3dB_cr_fft3b(1,v_out)
243         call D3dB_cr_pfft3b(1,0,v_out)
244
245         call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dbl_mb(rho1(1)))
246         call D3dB_r_notZero_Ends(1,dbl_mb(rho1(1)))
247         call D3dB_r_Power1(1,(wgc_tune1-1.0d0),dbl_mb(rho1(1)))
248         call D3dB_r_Zero_Ends(1,dbl_mb(rho1(1)))
249         call D3dB_rr_Multiply2(1,dbl_mb(rho1(1)),v_out)
250
251      !**** alpha!=beta ****
252      else
253         !**** the rho^alpha(G) is not being scaled by 1/(N1*N2*N3) ****
254         call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dbl_mb(rho1(1)))
255         call D3dB_r_notZero_Ends(1,dbl_mb(rho1(1)))
256         call D3dB_r_Power1(1,wgc_tune1,dbl_mb(rho1(1)))
257         call D3dB_r_Zero_Ends(1,dbl_mb(rho1(1)))
258         !call D3dB_rc_fft3f(1,dbl_mb(rho1(1)))
259         call D3dB_rc_pfft3f(1,0,dbl_mb(rho1(1)))
260         call Pack_c_pack(0,dbl_mb(rho1(1)))
261         call Pack_tc_Mul(0,dbl_mb(wgc_indx),
262     >                      dbl_mb(rho1(1)),
263     >                      v_out)
264         call Pack_c_SMul1(0,wgc_tune2,v_out)
265         call Pack_c_unpack(0,v_out)
266         !call D3dB_cr_fft3b(1,v_out)
267         call D3dB_cr_pfft3b(1,0,v_out)
268         call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dbl_mb(rho1(1)))
269         call D3dB_r_notZero_Ends(1,dbl_mb(rho1(1)))
270         call D3dB_r_Power1(1,(wgc_tune2-1.0d0),dbl_mb(rho1(1)))
271         call D3dB_rr_Multiply2(1,dbl_mb(rho1(1)),v_out)
272
273
274         !**** the rho^beta(G) is not being scaled by 1/(N1*N2*N3) ****
275         call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dbl_mb(rho1(1)))
276         call D3dB_r_notZero_Ends(1,dbl_mb(rho1(1)))
277         call D3dB_r_Power1(1,wgc_tune2,dbl_mb(rho1(1)))
278         call D3dB_r_Zero_Ends(1,dbl_mb(rho1(1)))
279         !call D3dB_rc_fft3f(1,dbl_mb(rho1(1)))
280         call D3dB_rc_pfft3f(1,0,dbl_mb(rho1(1)))
281         call Pack_c_pack(0,dbl_mb(rho1(1)))
282         call Pack_tc_Mul(0,dbl_mb(wgc_indx),
283     >                      dbl_mb(rho1(1)),
284     >                      dbl_mb(tmp1(1)))
285         call Pack_c_SMul1(0,wgc_tune1,dbl_mb(tmp1(1)))
286         call Pack_c_unpack(0,dbl_mb(tmp1(1)))
287         !call D3dB_cr_fft3b(1,dbl_mb(tmp1(1)))
288         call D3dB_cr_pfft3b(1,0,dbl_mb(tmp1(1)))
289         call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dbl_mb(rho1(1)))
290         call D3dB_r_notZero_Ends(1,dbl_mb(rho1(1)))
291         call D3dB_r_Power1(1,(wgc_tune1-1.0d0),dbl_mb(rho1(1)))
292         call D3dB_r_Zero_Ends(1,dbl_mb(rho1(1)))
293         call D3dB_rrr_MultiplyAdd(1,dbl_mb(rho1(1)),dbl_mb(tmp1(1)),
294     >                             v_out)
295      end if
296      call D3dB_r_Zero_Ends(1,v_out)
297
298      value =      BA_pop_stack(tmp1(2))
299     >        .and.BA_pop_stack(rho1(2))
300      if (.not.value) call errquit('wgc_v:popping stack',1,MA_ERR)
301
302      call nwpw_timing_end(7)
303
304      return
305      end
306
307
308*     ****************************************
309*     *                                      *
310*     *              wgc_e                   *
311*     *                                      *
312*     ****************************************
313      real*8 function wgc_e(ispin,dn)
314      implicit none
315      integer ispin
316      real*8  dn(*)
317
318#include "bafdecls.fh"
319#include "errquit.fh"
320#include "wgc.fh"
321
322
323*     **** local variables ****
324      logical value
325      integer npack0,nfft3d,n2ft3d
326      real*8 ec
327
328      integer tmp1(2),rho1(2),rho2(2)
329
330*     **** external functions ****
331      real*8   lattice_omega
332      external lattice_omega
333
334      call nwpw_timing_start(7)
335      call D3dB_nfft3d(1,nfft3d)
336      call D3dB_n2ft3d(1,n2ft3d)
337      call Pack_npack(0,npack0)
338
339      value =      BA_push_get(mt_dcpl,nfft3d,'rho1',rho1(2),rho1(1))
340     >        .and.BA_push_get(mt_dcpl,nfft3d,'rho2',rho2(2),rho2(1))
341     >        .and.BA_push_get(mt_dbl, npack0,'tmp1',tmp1(2),tmp1(1))
342      if (.not.value) call errquit('wgc_e:out of stack',0,MA_ERR)
343
344      !**** the rho^alpha(G) is not being scaled by 1/(N1*N2*N3) ****
345      call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dcpl_mb(rho1(1)))
346      call D3dB_r_notZero_Ends(1,dcpl_mb(rho1(1)))
347      call D3dB_r_Power1(1,wgc_tune1,dcpl_mb(rho1(1)))
348      call D3dB_r_Zero_Ends(1,dcpl_mb(rho1(1)))
349      !call D3dB_rc_fft3f(1,dcpl_mb(rho1(1)))
350      call D3dB_rc_pfft3f(1,0,dcpl_mb(rho1(1)))
351      call Pack_c_pack(0,dcpl_mb(rho1(1)))
352      if (wgc_eq_tune) then
353         call Pack_ct_Sqr(0,dcpl_mb(rho1(1)),dbl_mb(tmp1(1)))
354      else
355         call D3dB_rr_Sum(1,dn,dn(1+n2ft3d*(ispin-1)),dcpl_mb(rho2(1)))
356         call D3dB_r_notZero_Ends(1,dcpl_mb(rho2(1)))
357         call D3dB_r_Power1(1,wgc_tune2,dcpl_mb(rho2(1)))
358         call D3dB_r_Zero_Ends(1,dcpl_mb(rho2(1)))
359         !call D3dB_rc_fft3f(1,dcpl_mb(rho2(1)))
360         call D3dB_rc_pfft3f(1,0,dcpl_mb(rho2(1)))
361         call Pack_c_pack(0,dcpl_mb(rho2(1)))
362         call Pack_cct_conjgMul(0,dcpl_mb(rho1(1)),
363     >                            dcpl_mb(rho2(1)),
364     >                            dbl_mb(tmp1(1)))
365      end if
366      call Pack_tt_dot(0,dbl_mb(tmp1(1)),dbl_mb(wgc_indx),ec)
367
368      !**************************************************************************************************
369      !**** Note the energy is really omega*Sum(G) conj(rho^alpha(G))*rho^beta(G)|^2 * WGC_kernel(G) ****
370      !**** This funny form is because rho(G) is not scaled 1/(N1*N2*N3) and an                      ****
371      !**** extra scaling of 1/(N1*N2*N3) is done to the WGG_kernel(G)                               ****
372      !**** Might want to scal rho^alpha/beta(G) properly to make the code more transparent          ****
373      !**************************************************************************************************
374      ec = ec*lattice_omega()*wgc_scal1
375
376
377      value =      BA_pop_stack(tmp1(2))
378     >        .and.BA_pop_stack(rho2(2))
379     >        .and.BA_pop_stack(rho1(2))
380      if (.not.value) call errquit('wgc_e:popping stack',1,MA_ERR)
381
382      call nwpw_timing_end(7)
383
384      wgc_e = ec
385      return
386      end
387
388