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