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