1* 2* $Id$ 3* 4 integer function c_ewald_ncut() 5 implicit none 6 7 8* **** common block for c_ewald.f **** 9 integer ncut 10 real*8 rcut,cewald,alpha 11 integer vg(2),rcell(2) 12 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 13 14 c_ewald_ncut = ncut 15 return 16 end 17 18 real*8 function c_ewald_rcut() 19 implicit none 20 21* **** common block for c_ewald.f **** 22 integer ncut 23 real*8 rcut,cewald,alpha 24 integer vg(2),rcell(2) 25 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 26 27 c_ewald_rcut = rcut 28 return 29 end 30 31 integer function c_ewald_nshl3d() 32 implicit none 33 34 35* **** common block for c_ewald.f **** 36 integer ncut 37 real*8 rcut,cewald,alpha 38 integer vg(2),rcell(2) 39 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 40 41 c_ewald_nshl3d = (2*ncut+1)**3 42 return 43 end 44 45 46 real*8 function c_ewald_mandelung() 47 implicit none 48 49 50* **** common block for ewald.f **** 51 integer ncut 52 real*8 rcut,cewald,alpha 53 integer vg(2),rcell(2) 54 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 55 56 c_ewald_mandelung = alpha 57 return 58 end 59 60 61 subroutine c_ewald_end() 62 implicit none 63 64#include "bafdecls.fh" 65 66* **** common block for ewald.f **** 67 integer ncut 68 real*8 rcut,cewald,alpha 69 integer vg(2),rcell(2) 70 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 71 72 logical value 73 74 value = BA_free_heap(vg(2)) 75 value = BA_free_heap(rcell(2)) 76 77 return 78 end 79 80 81 subroutine c_ewald_init() 82 implicit none 83 84#include "bafdecls.fh" 85#include "errquit.fh" 86 87 88* **** common block for ewald.f **** 89 integer ncut 90 real*8 rcut,cewald,alpha 91 integer vg(2),rcell(2) 92 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 93 94* **** local variables **** 95 integer nfft3d,G(3) 96 integer nshl3d 97 integer i,j,k,l 98 real*8 pi,fourpi,gg,w 99 real*8 rs 100 real*8 zz,z 101 integer taskid,pzero,qzero,zero 102 integer nx,ny,nxh 103 logical value 104 real*8 kv(3),ecut 105 106* **** external functions **** 107 integer control_ncut 108 real*8 control_rcut,control_ecut 109 integer ion_nion,ion_katm,c_G_indx 110 real*8 lattice_omega,lattice_unita,cpsp_zv 111 112 external control_ncut 113 external control_rcut,control_ecut 114 external ion_nion,ion_katm,c_G_indx 115 external lattice_omega,lattice_unita,cpsp_zv 116 117 118* **** allocate vg memory **** 119 call C3dB_nfft3d(1,nfft3d) 120 value = BA_alloc_get(mt_dbl,nfft3d,'vg',vg(2),vg(1)) 121 if (.not. value) 122 > call errquit('c_ewald_init:out of heap memory',0,MA_ERR) 123 124 G(1) = c_G_indx(1) 125 G(2) = c_G_indx(2) 126 G(3) = c_G_indx(3) 127 128* **** get constants **** 129 pi = 4.0d0*datan(1.0d0) 130 fourpi = 4.0d0*pi 131 132 call Parallel_taskid(taskid) 133 call C3dB_nx(1,nx) 134 call C3dB_ny(1,ny) 135 nxh=nx/2 136 137* ***** find the G==0 index ****** 138 i=0 139 j=0 140 k=0 141c call C3dB_ktoqp(1,k+1,qzero,pzero) 142c zero = (qzero-1)*(nx)*ny 143c > + j*(nx) 144c > + i+1 145 call C3dB_ijktoindexp(1,i+1,j+1,k+1,zero,pzero) 146 147 148 149* ***** initialize common block and find w ***** 150 ncut = control_ncut() 151 rcut = control_rcut() 152 if (ncut.le.0) ncut=1 153 if (rcut.le.0.0d0) then 154 rs = lattice_unita(1,1)**2 155 > + lattice_unita(2,1)**2 156 > + lattice_unita(3,1)**2 157 rs = dsqrt(rs) 158 rcut=rs/pi 159 160 rs = lattice_unita(1,2)**2 161 > + lattice_unita(2,2)**2 162 > + lattice_unita(3,2)**2 163 rs = dsqrt(rs) 164 w=rs/pi 165 if (w.lt.rcut) rcut = w 166 167 rs = lattice_unita(1,3)**2 168 > + lattice_unita(2,3)**2 169 > + lattice_unita(3,3)**2 170 rs = dsqrt(rs) 171 w=rs/pi 172 if (w.lt.rcut) rcut = w 173 end if 174 175 w = 0.25d0*rcut*rcut 176 177 178* ***** initialize Vg ***** 179 do i=1,nfft3d 180 gg = ( dbl_mb(G(1)+i-1)*dbl_mb(G(1)+i-1) 181 > + dbl_mb(G(2)+i-1)*dbl_mb(G(2)+i-1) 182 > + dbl_mb(G(3)+i-1)*dbl_mb(G(3)+i-1) ) 183 184 if ((pzero.eq.taskid) .and. (i.eq.zero)) then 185 dbl_mb(vg(1)+i-1) = 0.0d0 186 else 187 dbl_mb(vg(1)+i-1) = (fourpi/gg)*exp(-w*gg) 188 end if 189 end do 190 191 kv(1) = 0.0d0 192 kv(2) = 0.0d0 193 kv(3) = 0.0d0 194 ecut = control_ecut() 195 call cloak_set(kv,ecut) 196 call cloak_R(dbl_mb(vg(1))) 197 198 199* **** set the Mandelung constant **** 200 call mandelung_set(alpha) 201 202 203* **** ewald summation **** 204 rs = (3.0d0*lattice_omega()/fourpi)**(1.0d0/3.0d0) 205 zz = 0.0d0 206 z = 0.0d0 207 do i=1,ion_nion() 208 zz = zz + cpsp_zv(ion_katm(i))**2 209 z = z + cpsp_zv(ion_katm(i)) 210 end do 211 call C3dB_r_dsum(1,dbl_mb(vg(1)),cewald) 212 cewald = -0.5d0*zz*(alpha/rs + cewald/lattice_omega()) 213 > -0.5d0*(z*z-zz)*rcut*rcut*pi/lattice_omega() 214 215* **** allocate rcell memory **** 216 nshl3d=(2*ncut+1)**3 217 value = BA_alloc_get(mt_dbl,(3*nshl3d),'rcell',rcell(2), 218 > rcell(1)) 219 if (.not. value) 220 > call errquit('c_ewald_init:out of heap memory',0,MA_ERR) 221 222 223* **** get lattice vectors in real space **** 224 l=0 225 do k=-ncut,ncut 226 do j=-ncut,ncut 227 do i=-ncut,ncut 228 l = l+1 229 dbl_mb(rcell(1)+ (l-1) ) 230 > = i*lattice_unita(1,1) 231 > + j*lattice_unita(1,2) 232 > + k*lattice_unita(1,3) 233 dbl_mb(rcell(1)+(l-1)+nshl3d) 234 > = i*lattice_unita(2,1) 235 > + j*lattice_unita(2,2) 236 > + k*lattice_unita(2,3) 237 dbl_mb(rcell(1)+(l-1)+2*nshl3d) 238 > = i*lattice_unita(3,1) 239 > + j*lattice_unita(3,2) 240 > + k*lattice_unita(3,3) 241 242 end do 243 end do 244 end do 245 246 247 return 248 end 249 250* *********************************** 251* * * 252* * c_ewald_e * 253* * * 254* *********************************** 255 real*8 function c_ewald_e() 256 implicit none 257 258#include "bafdecls.fh" 259#include "errquit.fh" 260 261* **** common block for ewald.f **** 262 integer ncut 263 real*8 rcut,cewald,alpha 264 integer vg(2),rcell(2) 265 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 266 267 268* **** local variables **** 269 integer dutask,taskid,np 270 integer i,j,ii,l 271 real*8 w,dx,dy,dz,x,y,z,r,zz 272 real*8 yerfc 273 real*8 energy,etmp 274 275* **** temporary workspace variables **** 276c complex*16 exi(nfft3d) 277c complex*16 s(nfft3d) 278c real*8 tmp3(nfft3d*2) 279 integer nfft3d,nshl3d 280 integer exi(2),s(2),tmp3(2),ft(2) 281 logical value 282 283* **** external functions **** 284 integer ion_nion,ion_katm,c_ewald_nshl3d 285 real*8 lattice_omega,cpsp_zv,dsum,ion_rion,util_erfc 286 external ion_nion,ion_katm,c_ewald_nshl3d 287 external lattice_omega,cpsp_zv,dsum,ion_rion,util_erfc 288 289 call Parallel_np(np) 290 call Parallel_taskid(taskid) 291 292* **** allocate temp workspace **** 293 call C3dB_nfft3d(1,nfft3d) 294 nshl3d = c_ewald_nshl3d() 295 value = BA_push_get(mt_dcpl,nfft3d,'exi',exi(2),exi(1)) 296 value = value.and. 297 > BA_push_get(mt_dcpl,nfft3d,'s',s(2),s(1)) 298 value = value.and. 299 > BA_push_get(mt_dbl, nfft3d,'tmp3',tmp3(2),tmp3(1)) 300 value = value.and. 301 > BA_push_get(mt_dbl, (3*nshl3d),'ft',ft(2),ft(1)) 302 if (.not. value) 303 > call errquit('c_ewald_e:out of stack memory',0,MA_ERR) 304 305* **** get the structure factor **** 306 call dcopy((2*nfft3d),0.0d0,0,dcpl_mb(s(1)),1) 307 do ii=1,ion_nion() 308 call cstrfac(ii,dcpl_mb(exi(1))) 309 call C3dB_cc_daxpy(1,cpsp_zv(ion_katm(ii)), 310 > dcpl_mb(exi(1)), 311 > dcpl_mb(s(1))) 312 313 end do 314 315 316* **** calculate the ewald energy **** 317 call C3dB_cr_Sqr(1,dcpl_mb(s(1)),dbl_mb(tmp3(1))) 318 call C3dB_rr_dot(1,dbl_mb(tmp3(1)),dbl_mb(vg(1)),energy) 319 energy = 0.5d0*energy/lattice_omega() + cewald 320 321* *** need to make parallel **** 322 dutask = 0 323 etmp = 0.0d0 324 do i=1,ion_nion()-1 325 do j=i+1,ion_nion() 326 if (dutask.eq.taskid) then 327 dx = ion_rion(1,i) - ion_rion(1,j) 328 dy = ion_rion(2,i) - ion_rion(2,j) 329 dz = ion_rion(3,i) - ion_rion(3,j) 330 zz = cpsp_zv(ion_katm(i)) * cpsp_zv(ion_katm(j)) 331 do l=1,nshl3d 332 x = dbl_mb(rcell(1)+(l-1)) + dx 333 y = dbl_mb(rcell(1)+(l-1)+nshl3d) + dy 334 z = dbl_mb(rcell(1)+(l-1)+2*nshl3d) + dz 335 r = dsqrt(x*x+y*y+z*z) 336 w = r/rcut 337 338c erfc=1.0d0/(1.0d0+w*(b1+w*(b2+w*(b3 339c > +w*(b4+w*(b5+w*b6))))))**4 340c dbl_mb(ft(1)+(l-1))=zz*erfc**4/r 341 yerfc = util_erfc(w) 342 dbl_mb(ft(1)+(l-1))=zz*yerfc/r 343 end do 344 etmp = etmp + dsum(nshl3d,dbl_mb(ft(1)),1) 345 end if 346 dutask = mod(dutask+1,np) 347 end do 348 end do 349 if (np.gt.1) call C3dB_SumAll(etmp) 350 energy = energy + etmp 351 352 353* **** deallocate temp workspace **** 354 value = BA_pop_stack(ft(2)) 355 value = value.and.BA_pop_stack(tmp3(2)) 356 value = value.and.BA_pop_stack(s(2)) 357 value = value.and.BA_pop_stack(exi(2)) 358 if (.not. value) 359 > call errquit('c_ewald_e:error popping stack',0,MA_ERR) 360 361 c_ewald_e = energy 362 return 363 end 364 365 366* *********************************** 367* * * 368* * c_ewald_f * 369* * * 370* *********************************** 371 372 subroutine c_ewald_f(fion) 373 implicit none 374 real*8 fion(3,*) 375 376#include "bafdecls.fh" 377#include "errquit.fh" 378 379* **** common block for ewald.f **** 380 integer ncut 381 real*8 rcut,cewald,alpha 382 integer vg(2),rcell(2) 383 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 384 385 386* **** expansion coefficient of the error function **** 387 real*8 cerfc 388 parameter (cerfc=1.128379167d0) 389 390 391* **** local variables **** 392 integer dutask,taskid,np 393 integer i,j,l,ii 394 real*8 w,dx,dy,dz,x,y,z,r,zz 395 real*8 yerfc 396 real*8 sum,scal2,f 397 real*8 sw1,sw2,sw3 398 399* **** temporary workspace variables **** 400c complex*16 exi(nfft3d) 401c complex*16 s(nfft3d) 402c real*8 tmp3(nfft3d*2) 403 integer nfft3d,nshl3d,nion 404 integer exi(2),s(2),tmp3(2),ft(2) 405 integer fx(2),fy(2),fz(2) 406 logical value 407 408* **** external functions **** 409 integer ion_nion,ion_katm,c_G_indx,c_ewald_nshl3d 410 external ion_nion,ion_katm,c_G_indx,c_ewald_nshl3d 411 real*8 lattice_omega,cpsp_zv,dsum,ion_rion,util_erfc 412 external lattice_omega,cpsp_zv,dsum,ion_rion,util_erfc 413 414 call Parallel_np(np) 415 call Parallel_taskid(taskid) 416 nion = ion_nion() 417 418* **** allocate temp workspace **** 419 call C3dB_nfft3d(1,nfft3d) 420 nshl3d = c_ewald_nshl3d() 421 value = BA_push_get(mt_dcpl,nfft3d,'exi',exi(2),exi(1)) 422 value = value.and. 423 > BA_push_get(mt_dcpl,nfft3d,'s',s(2),s(1)) 424 value = value.and. 425 > BA_push_get(mt_dbl, nfft3d,'tmp3',tmp3(2),tmp3(1)) 426 value = value.and. 427 > BA_push_get(mt_dbl, (3*nshl3d),'ft',ft(2),ft(1)) 428 value = value.and. 429 > BA_push_get(mt_dbl, (nion),'fx',fx(2),fx(1)) 430 value = value.and. 431 > BA_push_get(mt_dbl, (nion),'fy',fy(2),fy(1)) 432 value = value.and. 433 > BA_push_get(mt_dbl, (nion),'fz',fz(2),fz(1)) 434 if (.not. value) 435 > call errquit('c_ewald_f:out of stack memory',0,MA_ERR) 436 437 438 scal2 = 1.0d0/lattice_omega() 439 call dcopy(nion,0.0d0,0,dbl_mb(fx(1)),1) 440 call dcopy(nion,0.0d0,0,dbl_mb(fy(1)),1) 441 call dcopy(nion,0.0d0,0,dbl_mb(fz(1)),1) 442 443* **** get the structure factor **** 444 call dcopy((2*nfft3d),0.0d0,0,dcpl_mb(s(1)),1) 445 do ii=1,nion 446 call cstrfac(ii,dcpl_mb(exi(1))) 447 call C3dB_cc_daxpy(1,cpsp_zv(ion_katm(ii)),dcpl_mb(exi(1)), 448 > dcpl_mb(s(1))) 449 end do 450 451 do ii=1,nion 452 call cstrfac(ii,dcpl_mb(exi(1))) 453 454 do i=1,nfft3d 455 dbl_mb(tmp3(1)+i-1) 456 > = ( dble(dcpl_mb(exi(1)+i-1)) 457 > *dimag(dcpl_mb(s(1)+i-1)) 458 > - dimag(dcpl_mb(exi(1)+i-1)) 459 > *dble(dcpl_mb(s(1)+i-1)) 460 > )*dbl_mb(vg(1)+i-1) 461 end do 462 463 call C3dB_rr_idot(1,dbl_mb(c_G_indx(1)),dbl_mb(tmp3(1)),sum) 464* fion(1,ii) = fion(1,ii) + sum*cpsp_zv(ion_katm(ii))*scal2 465 dbl_mb(fx(1)+ii-1) = dbl_mb(fx(1)+ii-1) 466 > + sum*cpsp_zv(ion_katm(ii))*scal2 467 468 call C3dB_rr_idot(1,dbl_mb(c_G_indx(2)),dbl_mb(tmp3(1)),sum) 469* fion(2,ii) = fion(2,ii) + sum*cpsp_zv(ion_katm(ii))*scal2 470 dbl_mb(fy(1)+ii-1) = dbl_mb(fy(1)+ii-1) 471 > + sum*cpsp_zv(ion_katm(ii))*scal2 472 473 call C3dB_rr_idot(1,dbl_mb(c_G_indx(3)),dbl_mb(tmp3(1)),sum) 474* fion(3,ii) = fion(3,ii) + sum*cpsp_zv(ion_katm(ii))*scal2 475 dbl_mb(fz(1)+ii-1) = dbl_mb(fz(1)+ii-1) 476 > + sum*cpsp_zv(ion_katm(ii))*scal2 477 end do 478 479 480 481 dutask=0 482 do i=1,nion-1 483 do j=i+1,nion 484 if (dutask.eq.taskid) then 485 dx = ion_rion(1,i) - ion_rion(1,j) 486 dy = ion_rion(2,i) - ion_rion(2,j) 487 dz = ion_rion(3,i) - ion_rion(3,j) 488 zz = cpsp_zv(ion_katm(i)) * cpsp_zv(ion_katm(j)) 489 do l=1,nshl3d 490 x = dbl_mb(rcell(1)+(l-1)) + dx 491 y = dbl_mb(rcell(1)+(l-1)+ nshl3d) + dy 492 z = dbl_mb(rcell(1)+(l-1)+2*nshl3d) + dz 493 r = dsqrt(x*x+y*y+z*z) 494 w = r/rcut 495 496c erfc=(1.0d0+w*(b1+w*(b2+w*(b3 497c > +w*(b4+w*(b5+w*b6))))))**4 498c erfc = 1.0d0/erfc**4 499 yerfc = util_erfc(w) 500 f = zz*(yerfc+cerfc*w*dexp(-w*w))/r**3 501 dbl_mb(ft(1)+(l-1)) =x*f 502 dbl_mb(ft(1)+(l-1)+nshl3d) =y*f 503 dbl_mb(ft(1)+(l-1)+2*nshl3d)=z*f 504 end do 505 sw1 = dsum(nshl3d,dbl_mb(ft(1)),1) 506 sw2 = dsum(nshl3d,dbl_mb(ft(1)+ nshl3d),1) 507 sw3 = dsum(nshl3d,dbl_mb(ft(1)+2*nshl3d),1) 508 509* fion(1,i) = fion(1,i) + sw1 510* fion(2,i) = fion(2,i) + sw2 511* fion(3,i) = fion(3,i) + sw3 512* fion(1,j) = fion(1,j) - sw1 513* fion(2,j) = fion(2,j) - sw2 514* fion(3,j) = fion(3,j) - sw3 515 516 dbl_mb(fx(1)+i-1) = dbl_mb(fx(1)+i-1) + sw1 517 dbl_mb(fy(1)+i-1) = dbl_mb(fy(1)+i-1) + sw2 518 dbl_mb(fz(1)+i-1) = dbl_mb(fz(1)+i-1) + sw3 519 520 dbl_mb(fx(1)+j-1) = dbl_mb(fx(1)+j-1) - sw1 521 dbl_mb(fy(1)+j-1) = dbl_mb(fy(1)+j-1) - sw2 522 dbl_mb(fz(1)+j-1) = dbl_mb(fz(1)+j-1) - sw3 523 524 end if 525 dutask = mod((dutask+1),np) 526 end do 527 end do 528 529 if (np.gt.1) call C3dB_Vector_SumAll(nion,dbl_mb(fx(1))) 530 if (np.gt.1) call C3dB_Vector_SumAll(nion,dbl_mb(fy(1))) 531 if (np.gt.1) call C3dB_Vector_SumAll(nion,dbl_mb(fz(1))) 532 do i=1,nion 533 fion(1,i) = fion(1,i) + dbl_mb(fx(1)+i-1) 534 fion(2,i) = fion(2,i) + dbl_mb(fy(1)+i-1) 535 fion(3,i) = fion(3,i) + dbl_mb(fz(1)+i-1) 536 end do 537 538* **** deallocate temp workspace **** 539 value = BA_pop_stack(fz(2)) 540 value = value.and.BA_pop_stack(fy(2)) 541 value = value.and.BA_pop_stack(fx(2)) 542 value = value.and.BA_pop_stack(ft(2)) 543 value = value.and.BA_pop_stack(tmp3(2)) 544 value = value.and.BA_pop_stack(s(2)) 545 value = value.and.BA_pop_stack(exi(2)) 546 if (.not. value) 547 > call errquit('c_ewald_f:error popping stack memory',0,MA_ERR) 548 549 return 550 end 551 552 553* *********************************** 554* * * 555* * c_ewald_stress * 556* * * 557* *********************************** 558 559 subroutine c_ewald_stress(stress) 560 implicit none 561 real*8 stress(3,3) 562 563#include "bafdecls.fh" 564#include "errquit.fh" 565 566 integer N 567 parameter (N=40) 568 569* **** common block for c_ewald.f **** 570 integer ncut 571 real*8 rcut,cewald,alpha 572 integer vg(2),rcell(2) 573 common / c_ewald_block / vg,rcell,cewald,alpha,rcut,ncut 574 575 576* **** common block used for coulomb.f **** 577 integer vc_indx,vc_hndl 578 common / c_vc_block / vc_indx,vc_hndl 579 580 581* **** expansion coefficient of the error function **** 582 real*8 cerfc 583 parameter (cerfc=1.128379167d0) 584c real*8 cerfc,b1,b2,b3,b4,b5,b6 585c parameter (b1=0.0705230784d0,b2=0.0422820123d0,b3=0.0092705272d0) 586c parameter (b4=0.0001520143d0,b5=0.0002765672d0,b6=0.0000430638d0) 587 588* **** local variables **** 589 logical value 590 integer npack0,nfft3d 591 integer i,ii,j,l 592 integer n1,n2,n3 593 integer u,v,s 594 real*8 pi,fourpi,scal 595 real*8 zz,z 596 real*8 Cus(3,3),hm(3,3),energy,sum,ss,rs 597 real*8 ea,ax,ay,az,epsilon 598 real*8 dx,dy,dz,w 599 real*8 unita(3,3),unitg(3,3) 600 601 integer G(2,3),H(2),F(2),tmp1(2),tmp2(2),exi(2),strf(2) 602 integer nshl3d 603 604* **** external functions **** 605 integer ion_katm,ion_nion,c_G_indx,c_ewald_nshl3d 606 real*8 cpsp_zv,lattice_unitg,lattice_unita,lattice_omega 607 real*8 util_erfc,ion_rion 608 external ion_katm,ion_nion,c_G_indx,c_ewald_nshl3d 609 external cpsp_zv,lattice_unitg,lattice_unita,lattice_omega 610 external util_erfc,ion_rion 611 612 pi = 4.0d0*datan(1.0d0) 613 fourpi = 4.0d0*pi 614 scal = 1.0d0/(2.0d0*pi) 615 616* *** define hm,unita,unitg **** 617 do v=1,3 618 do u=1,3 619 hm(u,v) = scal*lattice_unitg(u,v) 620 unitg(u,v) = lattice_unitg(u,v) 621 unita(u,v) = lattice_unita(u,v) 622 end do 623 end do 624 625 626 call C3dB_nfft3d(1,nfft3d) 627 call Cram_npack(0,npack0) 628 629 630 zz = 0.0d0 631 z = 0.0d0 632 do i=1,ion_nion() 633 zz = zz + cpsp_zv(ion_katm(i))**2 634 z = z + cpsp_zv(ion_katm(i)) 635 end do 636 637* **** Miscellaneous contributions - stress from cewald term **** 638 do v=1,3 639 do u=1,3 640 stress(u,v) = 0.5d0*z*z*pi*rcut*rcut/lattice_omega() 641 > *hm(u,v) 642 end do 643 end do 644 645 646* **** G-space contributions **** 647 648* **** get the structure factor **** 649 value = BA_push_get(mt_dbl,npack0,'H',H(2),H(1)) 650 value = value.and.BA_push_get(mt_dcpl,nfft3d,'exi',exi(2),exi(1)) 651 value = value.and. 652 > BA_push_get(mt_dcpl,npack0,'strf',strf(2),strf(1)) 653 if (.not. value) 654 > call errquit('c_ewald_stress:out of stack memory',0, MA_ERR) 655 656 call dcopy((2*npack0),0.0d0,0,dcpl_mb(strf(1)),1) 657 do ii=1,ion_nion() 658 call cstrfac(ii,dcpl_mb(exi(1))) 659 call Cram_c_pack(0,dcpl_mb(exi(1))) 660 call Cram_cc_daxpy(0,cpsp_zv(ion_katm(ii)),dcpl_mb(exi(1)), 661 > dcpl_mb(strf(1))) 662 end do 663 call Cram_cr_Sqr(0,dcpl_mb(strf(1)),dbl_mb(H(1))) 664 value = BA_pop_stack(strf(2)) 665 value = value.and.BA_pop_stack(exi(2)) 666 if (.not. value) 667 > call errquit('c_ewald_stress:error popping stack',0,MA_ERR) 668 669* **** calculate the ewald energy **** 670 value = BA_push_get(mt_dbl,nfft3d,'F',F(2),F(1)) 671 if (.not. value) 672 > call errquit('c_ewald_stress:out of stack memory',0,MA_ERR) 673 call dcopy(nfft3d,dbl_mb(vg(1)), 1,dbl_mb(F(1)), 1) 674 call Cram_r_Pack(0,dbl_mb(F(1))) 675 676 call Cram_rr_dot(0,dbl_mb(F(1)),dbl_mb(H(1)),energy) 677 energy = -0.5d0*energy/lattice_omega() 678 679 680 do v=1,3 681 do u=1,3 682 stress(u,v) = stress(u,v) + energy*hm(u,v) 683 end do 684 end do 685 686* **** tmp2(G) = F(G)*H(G)/G**2 + F(G)*H(G)*rcut*rcut/4 **** 687 value = BA_push_get(mt_dbl,npack0,'tmp1',tmp1(2),tmp1(1)) 688 value = value.and. 689 > BA_push_get(mt_dbl,npack0,'tmp2',tmp2(2),tmp2(1)) 690 if (.not. value) 691 > call errquit('c_ewald_stress:out of stack memory',0,MA_ERR) 692 693 call Cram_rr_Mul(0,dbl_mb(F(1)), 694 > dbl_mb(H(1)), 695 > dbl_mb(tmp1(1))) 696 ss = 0.25d0*rcut*rcut 697 call Cram_r_SMul(0,ss,dbl_mb(tmp1(1)), 698 > dbl_mb(tmp2(1))) 699 ss = 1.0d0/fourpi 700c call Cram_r_SMul(0,ss,dbl_mb(tmp1(1)), 701c > dbl_mb(tmp1(1))) 702 call Cram_r_SMul1(0,ss,dbl_mb(tmp1(1))) 703c call Cram_rr_Mul(0,dbl_mb(tmp1(1)), 704c > dbl_mb(vc_indx), 705c > dbl_mb(tmp1(1))) 706 call Cram_rr_Mul2(0,dbl_mb(vc_indx), 707 > dbl_mb(tmp1(1))) 708 709c call Cram_rr_Sum(0,dbl_mb(tmp1(1)), 710c > dbl_mb(tmp2(1)), 711c > dbl_mb(tmp2(1))) 712 call Cram_rr_Sum2(0,dbl_mb(tmp1(1)), 713 > dbl_mb(tmp2(1))) 714 715 716* **** calculate Cus **** 717 value = BA_push_get(mt_dbl,nfft3d, 718 > 'G1',G(2,1),G(1,1)) 719 value = value.and.BA_push_get(mt_dbl,nfft3d, 720 > 'G2',G(2,2),G(1,2)) 721 value = value.and.BA_push_get(mt_dbl,nfft3d, 722 > 'G3',G(2,3),G(1,3)) 723 if (.not. value) 724 > call errquit('c_ewald_stress:out of stack memory',0,MA_ERR) 725 call dcopy(nfft3d,dbl_mb(c_G_indx(1)),1,dbl_mb(G(1,1)),1) 726 call dcopy(nfft3d,dbl_mb(c_G_indx(2)),1,dbl_mb(G(1,2)),1) 727 call dcopy(nfft3d,dbl_mb(c_G_indx(3)),1,dbl_mb(G(1,3)),1) 728 call Cram_r_pack(0,dbl_mb(G(1,1))) 729 call Cram_r_pack(0,dbl_mb(G(1,2))) 730 call Cram_r_pack(0,dbl_mb(G(1,3))) 731 732 call dcopy(9,0.0d0,0,Cus,1) 733c ss = -1.0d0/lattice_omega() 734 ss = 1.0d0/lattice_omega() 735 do u=1,3 736 do s=u,3 737 call Cram_rr_Mul(0,dbl_mb(G(1,u)), 738 > dbl_mb(G(1,s)), 739 > dbl_mb(tmp1(1))) 740 call Cram_rr_dot(0,dbl_mb(tmp1(1)),dbl_mb(tmp2(1)),sum) 741 Cus(u,s) = ss*sum 742 end do 743 end do 744 do u=1,3 745 do s=u+1,3 746 Cus(s,u) = Cus(u,s) 747 end do 748 end do 749 do v=1,3 750 do u=1,3 751 do s=1,3 752 stress(u,v) = stress(u,v) + Cus(u,s)*hm(s,v) 753 end do 754 end do 755 end do 756 757 value = BA_pop_stack(G(2,3)) 758 value = value.and.BA_pop_stack(G(2,2)) 759 value = value.and.BA_pop_stack(G(2,1)) 760 value = value.and.BA_pop_stack(tmp2(2)) 761 value = value.and.BA_pop_stack(tmp1(2)) 762 value = value.and.BA_pop_stack(F(2)) 763 value = value.and.BA_pop_stack(H(2)) 764 if (.not. value) 765 > call errquit('c_ewald_stress:error popping stack',0,MA_ERR) 766 767 768* **** R-space contributions **** 769 770 771* **** calculate alpha1 - stress from cewald term***** 772 call dcopy(9,0.0d0,0,Cus,1) 773 rs = (3.0d0*lattice_omega()/(4.0d0*pi))**(1.0d0/3.0d0) 774 epsilon = 1.0d0/rcut 775 sum = 0.0d0 776 do n1=(-N+1),(N-1) 777 do n2=(-N+1),(N-1) 778 do n3=(-N+1),(N-1) 779 if (.not.((n1.eq.0).and.(n2.eq.0).and.(n3.eq.0))) then 780 ax = n1*unita(1,1) 781 > + n2*unita(1,2) 782 > + n3*unita(1,3) 783 784 ay = n1*unita(2,1) 785 > + n2*unita(2,2) 786 > + n3*unita(2,3) 787 788 az = n1*unita(3,1) 789 > + n2*unita(3,2) 790 > + n3*unita(3,3) 791 792 ea = dsqrt(ax*ax + ay*ay + az*az) 793 w = ea*epsilon 794 795 ss = util_erfc(w)/ea 796 > + 2.0d0*epsilon/dsqrt(pi)*dexp(-w*w) 797 ss = -(0.5d0*zz)*ss/(ea*ea) 798 Cus(1,1) = Cus(1,1) + ss * ax*ax 799 Cus(1,2) = Cus(1,2) + ss * ax*ay 800 Cus(1,3) = Cus(1,3) + ss * ax*az 801 802 Cus(2,1) = Cus(2,1) + ss * ay*ax 803 Cus(2,2) = Cus(2,2) + ss * ay*ay 804 Cus(2,3) = Cus(2,3) + ss * ay*az 805 806 Cus(3,1) = Cus(3,1) + ss * az*ax 807 Cus(3,2) = Cus(3,2) + ss * az*ay 808 Cus(3,3) = Cus(3,3) + ss * az*az 809 810 end if 811 end do 812 end do 813 end do 814 815c do u=1,3 816c do s=u+1,3 817c Cus(s,u) = Cus(u,s) 818c end do 819c end do 820 821 do v=1,3 822 do u=1,3 823 do s=1,3 824 stress(u,v) = stress(u,v) + Cus(u,s)*hm(s,v) 825 end do 826 end do 827 end do 828 829* *** need to make parallel **** 830* **** calculate erfc contribution ***** 831 nshl3d = c_ewald_nshl3d() 832 call dcopy(9,0.0d0,0,Cus,1) 833 epsilon = 1.0d0/rcut 834 do i=1,ion_nion()-1 835 do j=i+1,ion_nion() 836 dx = ion_rion(1,i) - ion_rion(1,j) 837 dy = ion_rion(2,i) - ion_rion(2,j) 838 dz = ion_rion(3,i) - ion_rion(3,j) 839 zz = cpsp_zv(ion_katm(i)) * cpsp_zv(ion_katm(j)) 840 do l=1,nshl3d 841 ax = dbl_mb(rcell(1)+(l-1)) + dx 842 ay = dbl_mb(rcell(1)+(l-1)+nshl3d) + dy 843 az = dbl_mb(rcell(1)+(l-1)+2*nshl3d) + dz 844 ea = dsqrt(ax*ax+ay*ay+az*az) 845 w = ea*epsilon 846 847 ss = -util_erfc(w)/ea 848 > - 2.0d0*epsilon/dsqrt(pi)*exp(-w*w) 849 ss = ss/(ea*ea) 850 Cus(1,1) = Cus(1,1) + ss * ax*ax * zz 851 Cus(1,2) = Cus(1,2) + ss * ax*ay * zz 852 Cus(1,3) = Cus(1,3) + ss * ax*az * zz 853 Cus(2,2) = Cus(2,2) + ss * ay*ay * zz 854 Cus(2,3) = Cus(2,3) + ss * ay*az * zz 855 Cus(3,3) = Cus(3,3) + ss * az*az * zz 856 857 end do 858 end do 859 end do 860 do u=1,3 861 do s=u+1,3 862 Cus(s,u) = Cus(u,s) 863 end do 864 end do 865 866 do v=1,3 867 do u=1,3 868 do s=1,3 869 stress(u,v) = stress(u,v) + Cus(u,s)*hm(s,v) 870 end do 871 end do 872 end do 873 874 return 875 end 876