1*
2* $Id$
3*
4      subroutine ewald_strfact_init()
5      implicit none
6#include "bafdecls.fh"
7#include "errquit.fh"
8#include "ewald_strfac.fh"
9cccccccc arguments
10cccccccc locals
11      integer nion
12      logical value
13cccccccc externals
14      integer  ion_nion,ewald_npack
15      external ion_nion,ewald_npack
16      integer  ewald_grid_nx,ewald_grid_ny,ewald_grid_nz
17      external ewald_grid_nx,ewald_grid_ny,ewald_grid_nz
18cccccccc exec
19      nion  = ion_nion()
20
21      npack = ewald_npack()
22      nx    = ewald_grid_nx()
23      ny    = ewald_grid_ny()
24      nz    = ewald_grid_nz()
25      notzero_npack = (npack.gt.0)
26
27      value=BA_alloc_get(mt_dcpl,(nx*nion),'ewx1',ewx1(2),ewx1(1))
28      value=value.and.
29     >      BA_alloc_get(mt_dcpl,(ny*nion),'ewx2',ewx2(2),ewx2(1))
30      value=value.and.
31     >      BA_alloc_get(mt_dcpl,(nz*nion),'ewx3',ewx3(2),ewx3(1))
32
33      if (notzero_npack) then
34         value=BA_alloc_get(mt_int,npack,'i_indxe',i_indx(2),i_indx(1))
35         value=value.and.
36     >         BA_alloc_get(mt_int,npack,'j_indxe',j_indx(2),j_indx(1))
37         value=value.and.
38     >         BA_alloc_get(mt_int,npack,'k_indxe',k_indx(2),k_indx(1))
39      end if
40      if (.not.value)
41     >  call errquit("EWALD_STRUCTFACTOR_INIT: OUT OF HEAP MEMORY",0,
42     >               MA_ERR)
43      return
44      end
45cccccccccccccccccccccccc
46      subroutine ewald_strfact_end()
47      implicit none
48#include "bafdecls.fh"
49#include "errquit.fh"
50#include "ewald_strfac.fh"
51      logical value
52cccccccc exec
53      value=          BA_free_heap(ewx1(2))
54      value=value.and.BA_free_heap(ewx2(2))
55      value=value.and.BA_free_heap(ewx3(2))
56
57      if (notzero_npack) then
58         value=value.and.BA_free_heap(i_indx(2))
59         value=value.and.BA_free_heap(j_indx(2))
60         value=value.and.BA_free_heap(k_indx(2))
61      end if
62      if (.not.value) then
63        call errquit("ewald_strfact_end: error free heap",0,MA_ERR)
64      end if
65      return
66      end
67ccccccccccccccccccccccccccccc
68      integer function ewald_strfact_i_indx()
69      implicit none
70#include "ewald_strfac.fh"
71      ewald_strfact_i_indx = i_indx(1)
72      return
73      end
74      integer function ewald_strfact_j_indx()
75      implicit none
76#include "ewald_strfac.fh"
77      ewald_strfact_j_indx = j_indx(1)
78      return
79      end
80      integer function ewald_strfact_k_indx()
81      implicit none
82#include "ewald_strfac.fh"
83      ewald_strfact_k_indx = k_indx(1)
84      return
85      end
86ccccccccccccccccccccccccccccc
87      subroutine ewald_phafac()
88      implicit none
89#include "bafdecls.fh"
90#include "ewald_strfac.fh"
91cccccccc locals
92      integer i,k
93      integer nxh,nyh,nzh
94      integer nion
95      complex*16 cw1,cw2,cw3
96      real*8 sw1,sw2,sw3,pi
97ccccccccc externals
98      integer ion_nion,ewald_grid_nx
99      integer ewald_grid_ny,ewald_grid_nz
100      real*8 lattice_unitg,ion_rion
101      external ion_nion
102      external ewald_grid_nz
103      external ewald_grid_ny,ewald_grid_nx
104      external lattice_unitg,ion_rion
105
106ccccccccc openmp
107      integer tid,nthr
108      integer Parallel_threadid, Parallel_nthreads
109      external Parallel_threadid, Parallel_nthreads
110
111ccccccc start of code
112      call nwpw_timing_start(8)
113
114      tid = Parallel_threadid()
115      nthr = Parallel_nthreads()
116
117      pi=4.0d0*datan(1.0d0)
118      nion=ion_nion()
119
120      nxh=nx/2
121      nyh=ny/2
122      nzh=nz/2
123c     Each thread processes independend ions (if USE_OPENMP is defined)
124      do i=tid+1,nion,nthr
125        sw1=lattice_unitg(1,1)*ion_rion(1,i)
126     >     +lattice_unitg(2,1)*ion_rion(2,i)
127     >     +lattice_unitg(3,1)*ion_rion(3,i)+pi
128
129        sw2=lattice_unitg(1,2)*ion_rion(1,i)
130     >     +lattice_unitg(2,2)*ion_rion(2,i)
131     >     +lattice_unitg(3,2)*ion_rion(3,i)+pi
132
133        sw3=lattice_unitg(1,3)*ion_rion(1,i)
134     >     +lattice_unitg(2,3)*ion_rion(2,i)
135     >     +lattice_unitg(3,3)*ion_rion(3,i)+pi
136
137        cw1=dcmplx(dcos(sw1),-dsin(sw1))
138        cw2=dcmplx(dcos(sw2),-dsin(sw2))
139        cw3=dcmplx(dcos(sw3),-dsin(sw3))
140        dcpl_mb(ewx1(1)+(i-1)*nx)=dcmplx(1.0d0,0.0d0)
141        dcpl_mb(ewx2(1)+(i-1)*ny)=dcmplx(1.0d0,0.0d0)
142        dcpl_mb(ewx3(1)+(i-1)*nz)=dcmplx(1.0d0,0.0d0)
143
144        do k=1,nxh
145          dcpl_mb(ewx1(1)+k+(i-1)*nx)
146     >         =dcpl_mb(ewx1(1)+k-1+(i-1)*nx)*cw1
147          dcpl_mb(ewx1(1)+nx-k+(i-1)*nx)
148     >         =dconjg(dcpl_mb(ewx1(1)+k+(i-1)*nx))
149        end do
150
151        do k=1,nyh
152          dcpl_mb(ewx2(1)+k+(i-1)*ny)
153     >         = dcpl_mb(ewx2(1)+k-1+(i-1)*ny)*cw2
154          dcpl_mb(ewx2(1)+ny-k+(i-1)*ny)
155     >         =dconjg(dcpl_mb(ewx2(1)+k+(i-1)*ny))
156        end do
157
158        do k=1,nzh
159          dcpl_mb(ewx3(1)+k+(i-1)*nz)
160     >         = dcpl_mb(ewx3(1)+k-1+(i-1)*nz)*cw3
161          dcpl_mb(ewx3(1)+nz-k+(i-1)*nz)
162     >         =dconjg(dcpl_mb(ewx3(1)+k+(i-1)*nz))
163        end do
164        dcpl_mb(ewx1(1)+nxh+(i-1)*nx)=dcmplx(0.0d0, 0.0d0)
165        dcpl_mb(ewx2(1)+nyh+(i-1)*ny)=dcmplx(0.0d0, 0.0d0)
166        dcpl_mb(ewx3(1)+nzh+(i-1)*nz)=dcmplx(0.0d0, 0.0d0)
167      end do
168      call nwpw_timing_end(8)
169
170      return
171      end
172ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
173      subroutine ewald_strfac(ii,exi)
174      implicit none
175      integer ii
176      complex*16 exi(*)
177#include "bafdecls.fh"
178#include "ewald_strfac.fh"
179ccccccc locals
180
181      call nwpw_timing_start(8)
182      if (notzero_npack) then
183         call ewald_strfac_sub(npack,
184     >                     int_mb(i_indx(1)),
185     >                     int_mb(j_indx(1)),
186     >                     int_mb(k_indx(1)),
187     >                     dcpl_mb(ewx1(1)+(ii-1)*nx),
188     >                     dcpl_mb(ewx2(1)+(ii-1)*ny),
189     >                     dcpl_mb(ewx3(1)+(ii-1)*nz),
190     >                     exi)
191      end if
192      call nwpw_timing_end(8)
193      return
194      end
195cccccccccccccccccccccccccccccccccccccccccccccccccccccc
196      subroutine ewald_strfac_sub(npack,
197     >                            indxi,indxj,indxk,
198     >                            ex1,ex2,ex3,
199     >                            exi)
200      implicit none
201      integer npack,indxi(*),indxj(*),indxk(*)
202      complex*16 ex1(*),ex2(*),ex3(*),exi(*)
203      integer indx
204!$OMP DO
205      do indx=1,npack
206        exi(indx) = ex1(indxi(indx))*ex2(indxj(indx))*ex3(indxk(indx))
207      end do
208!$OMP END DO
209      return
210      end
211cccccccccccccccccccccccccccccccccccccccccccccccc
212
213
214
215
216
217
218
219