1
2*
3* $Id$
4*
5
6      subroutine coulomb_init()
7      implicit none
8#include "errquit.fh"
9
10#include "bafdecls.fh"
11
12*     **** common block used for coulomb.f ****
13c     real*8 vc(nfft3d)
14      integer vc_indx,vc_hndl
15      common / vc_block / vc_indx,vc_hndl
16
17*     **** local variables ****
18      integer npack0,nfft3d,G(3)
19      integer i,j,k
20      integer zero,qzero,pzero,taskid
21      integer nx,ny,nxh,nyh
22      real*8  fourpi,gg
23      logical value
24      integer tmp1(2)
25
26*     **** external functions ****
27*     real*8 G(nfft3d,3)
28      integer  G_indx
29      external G_indx
30      double precision toll
31      parameter (toll=1d-16)
32
33      call nwpw_timing_start(7)
34      call D3dB_nfft3d(1,nfft3d)
35      call Pack_npack(0,npack0)
36      G(1) = G_indx(1)
37      G(2) = G_indx(2)
38      G(3) = G_indx(3)
39
40*     **** allocate vc memory ****
41      value = BA_alloc_get(mt_dbl,npack0,'vc',vc_hndl,vc_indx)
42      if (.not. value) call errquit('out of heap memory',0, MA_ERR)
43
44      value = BA_push_get(mt_dbl,nfft3d,'tmp1',tmp1(2),tmp1(1))
45      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
46
47      call Parallel2d_taskid_i(taskid)
48      call D3dB_nx(1,nx)
49      call D3dB_ny(1,ny)
50      nxh=nx/2
51      nyh=ny/2
52
53*     ***** find the G==0 point in the lattice *****
54      i=0
55      j=0
56      k=0
57c     call D3dB_ktoqp(1,k+1,qzero,pzero)
58c     zero = (qzero-1)*(nxh+1)*ny
59c    >     + j*(nxh+1)
60c    >     + i+1
61      call D3dB_ijktoindexp(1,i+1,j+1,k+1,zero,pzero)
62
63*     ***** form Vc = 4*pi/G**2  *****
64      fourpi = 4.0d0*(4.0d0*datan(1.0d0))
65      do i = 1,nfft3d
66
67         gg  = ( dbl_mb(G(1)+i-1)*dbl_mb(G(1)+i-1)
68     >         + dbl_mb(G(2)+i-1)*dbl_mb(G(2)+i-1)
69     >         + dbl_mb(G(3)+i-1)*dbl_mb(G(3)+i-1) )
70
71         if (((pzero.eq.taskid) .and. (i.eq.zero)).or.
72     E    (abs(gg) .lt.toll))  then
73            dbl_mb(tmp1(1)+i-1) = 0.0d0
74         else
75            dbl_mb(tmp1(1)+i-1) = fourpi/gg
76         end if
77
78      end do
79      call Pack_t_pack(0,dbl_mb(tmp1(1)))
80      call Pack_t_Copy(0,dbl_mb(tmp1(1)),dbl_mb(vc_indx))
81      value = BA_pop_stack(tmp1(2))
82
83      call nwpw_timing_end(7)
84
85
86      return
87      end
88
89      subroutine coulomb_end()
90      implicit none
91#include "bafdecls.fh"
92
93*     **** common block used for coulomb.f ****
94      integer vc_indx,vc_hndl
95      common / vc_block / vc_indx,vc_hndl
96      logical value
97
98      value = BA_free_heap(vc_hndl)
99      return
100      end
101
102
103      subroutine coulomb_v(dng,vc_out)
104      implicit none
105      complex*16 dng(*)
106      complex*16 vc_out(*)
107
108#include "bafdecls.fh"
109
110
111*     **** common block used for coulomb.f ****
112      integer vc_indx,vc_hndl
113      common / vc_block / vc_indx,vc_hndl
114
115      call nwpw_timing_start(7)
116      call Pack_tc_Mul(0,dbl_mb(vc_indx),dng,vc_out)
117      call nwpw_timing_end(7)
118
119      return
120      end
121
122
123      real*8 function coulomb_e(dng)
124      implicit none
125      complex*16 dng(*)
126
127#include "bafdecls.fh"
128#include "errquit.fh"
129
130*     **** common block used for coulomb.f ****
131*     real*8 vc(nfft3d)
132*     common / vc_block / vc
133      integer vc_indx,vc_hndl
134      common / vc_block / vc_indx,vc_hndl
135
136
137*     **** local variables ****
138      integer npack0
139      real*8 ec
140
141c     real*8  tmp1(*)
142      integer tmp1(2)
143      logical value
144
145*     **** external functions ****
146      real*8   lattice_omega
147      external lattice_omega
148
149      call nwpw_timing_start(7)
150      call Pack_npack(0,npack0)
151      value = BA_push_get(mt_dbl,npack0,'tmp1',tmp1(2),tmp1(1))
152      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
153
154      call Pack_ct_Sqr(0,dng,dbl_mb(tmp1(1)))
155      call Pack_ttp_dot(0,dbl_mb(tmp1(1)),dbl_mb(vc_indx),ec)
156
157      ec = 0.5d0*ec*lattice_omega()
158
159      value = BA_pop_stack(tmp1(2))
160      call nwpw_timing_end(7)
161
162      coulomb_e = ec
163      return
164      end
165
166
167
168
169
170      subroutine coulomb_euv(dng,euv)
171*
172* $Id$
173*
174      implicit none
175#include "errquit.fh"
176      complex*16 dng(*)
177      real*8 euv(3,3)
178
179#include "bafdecls.fh"
180
181*     **** common block used for coulomb.f ****
182      integer vc_indx,vc_hndl
183      common / vc_block / vc_indx,vc_hndl
184
185
186*     **** local variables ****
187      integer npack0,nfft3d,G(2,3)
188      integer i,j
189      integer u,v,s
190      logical value
191
192      real*8 pi,fourpi,scal,ss,sum
193      real*8 hm(3,3),Bus(3,3),ecoul
194      integer tmp1(2),tmp2(2)
195
196*     **** external functions ****
197c     real*8 G(nfft3d,3)
198      integer  G_indx
199      external G_indx
200
201      real*8   lattice_unitg,lattice_omega,coulomb_e
202      external lattice_unitg,lattice_omega,coulomb_e
203
204
205
206      pi     = 4.0d0*datan(1.0d0)
207      fourpi = 4.0d0*pi
208      scal   = 1.0d0/(2.0d0*pi)
209
210*     *** define hm ****
211      do j=1,3
212      do i=1,3
213         hm(i,j) = scal*lattice_unitg(i,j)
214      end do
215      end do
216
217
218      call D3dB_nfft3d(1,nfft3d)
219      call Pack_npack(0,npack0)
220
221      value = BA_push_get(mt_dbl,nfft3d,
222     >                     'G1',G(2,1),G(1,1))
223      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
224      value = BA_push_get(mt_dbl,nfft3d,
225     >                     'G2',G(2,2),G(1,2))
226      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
227      value = BA_push_get(mt_dbl,nfft3d,
228     >                     'G3',G(2,3),G(1,3))
229      if (.not. value) call errquit('out of stack  memory',0, MA_ERR)
230
231      value = BA_push_get(mt_dbl,npack0,'tmp1',tmp1(2),tmp1(1))
232      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
233
234      value = BA_push_get(mt_dbl,npack0,'tmp2',tmp2(2),tmp2(1))
235      if (.not. value) call errquit('out of stack memory',0, MA_ERR)
236
237      call dcopy(nfft3d,dbl_mb(G_indx(1)),1,dbl_mb(G(1,1)),1)
238      call dcopy(nfft3d,dbl_mb(G_indx(2)),1,dbl_mb(G(1,2)),1)
239      call dcopy(nfft3d,dbl_mb(G_indx(3)),1,dbl_mb(G(1,3)),1)
240      call Pack_t_pack(0,dbl_mb(G(1,1)))
241      call Pack_t_pack(0,dbl_mb(G(1,2)))
242      call Pack_t_pack(0,dbl_mb(G(1,3)))
243
244*     **** tmp2(G) = (n(G)**2) * (4*pi/G**2)**2  ****
245      call Pack_ct_Sqr(0,dng,dbl_mb(tmp1(1)))
246      call Pack_tt_Mul(0,dbl_mb(tmp1(1)),dbl_mb(vc_indx),
247     >                                   dbl_mb(tmp2(1)))
248c      call Pack_tt_Mul(0,dbl_mb(tmp2(1)),dbl_mb(vc_indx),
249c     >                                   dbl_mb(tmp2(1)))
250      call Pack_tt_Mul2(0,dbl_mb(vc_indx),dbl_mb(tmp2(1)))
251
252
253*     **** Bus = Sum(G) (omega/4*pi)*tmp2(G)*Gu*Gs ****
254      call dcopy(9,0.0d0,0,Bus,1)
255      ss     = lattice_omega()/fourpi
256      do u=1,3
257      do s=u,3
258        call Pack_tt_Mul(0,dbl_mb(G(1,u)),
259     >                     dbl_mb(G(1,s)),
260     >                     dbl_mb(tmp1(1)))
261        call Pack_tt_dot(0,dbl_mb(tmp1(1)),dbl_mb(tmp2(1)),sum)
262
263        Bus(u,s) = ss*sum
264      end do
265      end do
266      do u=1,3
267      do s=u+1,3
268         Bus(s,u) = Bus(u,s)
269      end do
270      end do
271
272      ecoul = coulomb_e(dng)
273      do v=1,3
274      do u=1,3
275         euv(u,v) = -ecoul*hm(u,v)
276         do s=1,3
277            euv(u,v) = euv(u,v) + Bus(u,s)*hm(s,v)
278         end do
279      end do
280      end do
281
282      value = BA_pop_stack(tmp2(2))
283      value = BA_pop_stack(tmp1(2))
284      value = BA_pop_stack(G(2,3))
285      value = BA_pop_stack(G(2,2))
286      value = BA_pop_stack(G(2,1))
287
288      return
289      end
290
291
292*     **********************************************
293*     *                                            *
294*     *              coulomb_efg                   *
295*     *                                            *
296*     **********************************************
297      subroutine coulomb_efg(dng,efg_smoothr,efg_smoothi)
298
299      implicit none
300
301      complex*16 dng(*)
302      real*8 efg_smoothr(3,3,*),efg_smoothi(3,3,*) ! real and complex parts
303
304#include "bafdecls.fh"
305#include "errquit.fh"
306
307*     **** common block used for coulomb.f ****
308      integer vc_indx,vc_hndl
309      common / vc_block / vc_indx,vc_hndl
310
311*     **** local variables ****
312      logical value
313      integer npack0,nfft3d
314      integer ii,k
315      integer u,v
316      integer Gtmp(2),G(3),tmp1(2),tmp2(2),exi(2)
317      real*8 w,gg
318      integer mu,nu,nion
319      real*8 rdng,cdng,gvec(3),termgg
320      real*8 pi,two_pi,four_pi,sqrt_pi
321      real*8 phase,cosgr,singr
322      real*8 r1,r2,r3
323      complex*16 zsum
324
325*     **** external functions ****
326      integer  G_indx,ion_nion
327      external G_indx,ion_nion
328      real*8   lattice_omega,ion_rion
329      external lattice_omega,ion_rion
330
331      pi = 4.0d0*datan(1.0d0)
332      sqrt_pi  = dsqrt(pi)
333      two_pi   = 2.0d0*pi
334      four_pi  = 4.0d0*pi
335
336      call D3dB_nfft3d(1,nfft3d)
337      call Pack_npack(0,npack0)
338
339      value = BA_push_get(mt_dbl,3*nfft3d,'Gtmp',Gtmp(2),Gtmp(1))
340      G(1) = Gtmp(1)
341      G(2) = Gtmp(1)+nfft3d
342      G(3) = Gtmp(2)+nfft3d
343      value = value.and.
344     >        BA_push_get(mt_dbl,npack0,'tmp1',tmp1(2),tmp1(1))
345      value = value.and.
346     >        BA_push_get(mt_dbl,npack0,'tmp2',tmp2(2),tmp2(1))
347      value = value.and.
348     >        BA_push_get(mt_dcpl,npack0,'exi',exi(2),exi(1))
349      if (.not. value)
350     > call errquit('coulomb_efg:out of stack memory',0, MA_ERR)
351
352      call dcopy(nfft3d,dbl_mb(G_indx(1)),1,dbl_mb(G(1)),1)
353      call dcopy(nfft3d,dbl_mb(G_indx(2)),1,dbl_mb(G(2)),1)
354      call dcopy(nfft3d,dbl_mb(G_indx(3)),1,dbl_mb(G(3)),1)
355      call Pack_t_pack(0,dbl_mb(G(1)))
356      call Pack_t_pack(0,dbl_mb(G(2)))
357      call Pack_t_pack(0,dbl_mb(G(3)))
358
359*     **** initialize ***
360      nion = ion_nion()
361      do mu = 1,3
362        do nu = 1,3
363          do ii=1,nion
364             efg_smoothr(mu,nu,ii) = 0.0d0
365             efg_smoothi(mu,nu,ii) = 0.0d0
366          end do ! ii
367        end do ! nu
368      end do ! mu
369
370      !*** compute (-1/3)gg ***
371       do k=1,npack0
372          dbl_mb(tmp1(1)+k-1) = -( dbl_mb(G(1)+k-1)**2
373     >                           + dbl_mb(G(2)+k-1)**2
374     >                           + dbl_mb(G(3)+k-1)**2)/3.0d0
375       end do ! k
376c
377       do ii=1,nion
378         call strfac_pack(0,ii,dcpl_mb(exi(1)))
379         do mu=1,3
380            call Pack_ttt_dzaxpy(0,dbl_mb(G(mu)),
381     >                             dbl_mb(G(mu)),
382     >                             dbl_mb(tmp1(1)),
383     >                             dbl_mb(tmp2(1)))
384            call Pack_tt_Mul2(0,dbl_mb(vc_indx),dbl_mb(tmp2(1)))
385
386            zsum = dcmplx(0.0d0,0.0d0)
387            do k=1,npack0
388               zsum = zsum
389     >              - 2.0d0*dng(k)
390     >                     *dcpl_mb(exi(1)+k-1)
391     >                     *dbl_mb(tmp2(1)+k-1)
392            end do
393            efg_smoothr(mu,mu,ii) = dble(zsum)
394            efg_smoothi(mu,mu,ii) = dimag(zsum)
395
396            do nu=mu+1,3
397               call Pack_tt_Mul(0,dbl_mb(G(mu)),
398     >                            dbl_mb(G(nu)),
399     >                            dbl_mb(tmp2(1)))
400               call Pack_tt_Mul2(0,dbl_mb(vc_indx),dbl_mb(tmp2(1)))
401
402               zsum = dcmplx(0.0d0,0.0d0)
403               do k=1,npack0
404                  zsum = zsum - 2.0d0*dng(k)
405     >                               *dcpl_mb(exi(1)+k-1)
406     >                               *dbl_mb(tmp2(1)+k-1)
407               end do
408               efg_smoothr(mu,nu,ii) = dble(zsum)
409               efg_smoothr(nu,mu,ii) = dble(zsum)
410               efg_smoothi(mu,nu,ii) = dimag(zsum)
411               efg_smoothi(nu,mu,ii) = dimag(zsum)
412
413            end do  ! nu
414         end do ! mu
415      end do ! ii
416
417c      do k=1,npack0
418c        gvec(1) = dbl_mb(G(1)+k-1)
419c        gvec(2) = dbl_mb(G(2)+k-1)
420c        gvec(3) = dbl_mb(G(3)+k-1)
421c        gg = gvec(1)*gvec(1) + gvec(2)*gvec(2) + gvec(3)*gvec(3)
422c        if (gg .gt. 0.d0) then  ! eliminate g=o term
423c          do ii=1,nion
424c            r1 = ion_rion(1,ii)
425c            r2 = ion_rion(2,ii)
426c            r3 = ion_rion(3,ii)
427c            phase = (gvec(1)*r1+gvec(2)*r2+gvec(3)*r3)
428c            cosgr = cos(phase)  ! cos(G.R)
429c            singr = sin(phase)  ! sin(G.R)
430c            do mu=1,3
431c             do nu=1,3
432c                termgg = gvec(mu)*gvec(nu)
433c                if (mu == nu) termgg = termgg - gg/3.d0
434c                rdng = real(dng(k))
435c                cdng = aimag(dng(k))
436c!
437c                efg_smoothr(mu,nu,ii) = efg_smoothr(mu,nu,ii) -   ! real part
438c     &            2.0d0*four_pi*termgg*(rdng*cosgr-cdng*singr)/gg
439c                efg_smoothi(mu,nu,ii) = efg_smoothi(mu,nu,ii) -   ! imag part
440c     &            2.0d0*four_pi*termgg*(rdng*singr+cdng*cosgr)/gg
441c             end do ! nu
442c            end do ! mu
443c          end do ! ii
444c        end if ! gg ne 0
445c      end do ! k
446
447      call D3dB_Vector_SumAll(9*nion,efg_smoothr)
448      call D3dB_Vector_SumAll(9*nion,efg_smoothi)
449!
450      value =           BA_pop_stack(exi(2))
451      value = value.and.BA_pop_stack(tmp2(2))
452      value = value.and.BA_pop_stack(tmp1(2))
453      value = value.and. BA_pop_stack(Gtmp(2))
454      if (.not. value) call errquit('coulomb_efg:popstack',1,MA_ERR)
455!
456      return
457      end
458