1 2* ****************************** 3* * * 4* * dipole_Efield_init * 5* * * 6* ****************************** 7 subroutine dipole_Efield_init(ispin,ne,efield0,efield0_center) 8 implicit none 9 integer ispin,ne(2) 10 real*8 efield0(3),efield0_center(3) 11 12#include "bafdecls.fh" 13#include "errquit.fh" 14 15 !**** dipole_efield common block **** 16 real*8 efield(3),efield_center(3) 17 real*8 dx(2),dy(2),dz(2) 18 real*8 nx,ny,nz,pcharge 19 real*8 dipole(3),Edipole,Pdipole 20 real*8 bv(3,6),wts(6) 21 integer rank 22 integer Tc(2),Ts(2),rgrid(2) 23 integer ispin0,ne0(2) 24 common /dipole_efield/ efield,efield_center, 25 > dx,dy,dz, 26 > nx,ny,nz,pcharge,dipole, 27 > Edipole,Pdipole, 28 > bv,wts,rank, 29 > Tc,Ts,rgrid, 30 > ispin0,ne0 31 32c *** local variables *** 33 logical value 34 integer n2ft3d 35 36 ispin0 = ispin 37 ne0(1) = ne(1) 38 ne0(2) = ne(2) 39 40 efield(1) = efield0(1) 41 efield(2) = efield0(2) 42 efield(3) = efield0(3) 43 44 efield_center(1) = efield0_center(1) 45 efield_center(2) = efield0_center(2) 46 efield_center(3) = efield0_center(3) 47 48 call dipole_Efield_set_ion() 49 call dipole_Efield_set_bv() 50 51c **** allocate from heap **** 52 call D3dB_n2ft3d(1,n2ft3d) 53 value = BA_alloc_get(mt_dbl,n2ft3d,'Tc',Tc(2),Tc(1)) 54 value = value.and.BA_alloc_get(mt_dbl,n2ft3d,'Ts',Ts(2),Ts(1)) 55 value = value.and. 56 > BA_alloc_get(mt_dbl,(3*n2ft3d),'rgrid',rgrid(2),rgrid(1)) 57 if (.not.value) 58 > call errquit('dipole_efield_init:out of memory',0,MA_ERR) 59 60c **** generate rgrid **** 61 call lattice_r_grid(dbl_mb(rgrid(1))) 62 63 return 64 end 65 66* ****************************** 67* * * 68* * dipole_Efield_end * 69* * * 70* ****************************** 71 subroutine dipole_Efield_end() 72 implicit none 73 74#include "bafdecls.fh" 75#include "errquit.fh" 76 77 !**** dipole_efield common block **** 78 real*8 efield(3),efield_center(3) 79 real*8 dx(2),dy(2),dz(2) 80 real*8 nx,ny,nz,pcharge 81 real*8 dipole(3),Edipole,Pdipole 82 real*8 bv(3,6),wts(6) 83 integer rank 84 integer Tc(2),Ts(2),rgrid(2) 85 integer ispin0,ne0(2) 86 common /dipole_efield/ efield,efield_center, 87 > dx,dy,dz, 88 > nx,ny,nz,pcharge,dipole, 89 > Edipole,Pdipole, 90 > bv,wts,rank, 91 > Tc,Ts,rgrid, 92 > ispin0,ne0 93 94c *** local variables *** 95 logical value 96 97c **** deallocate from heap **** 98 value = BA_free_heap(Tc(2)) 99 value = value.and.BA_free_heap(Ts(2)) 100 value = value.and.BA_free_heap(rgrid(2)) 101 if (.not.value) 102 > call errquit('dipole_efield_end:freeing memory',0,MA_ERR) 103 104 return 105 end 106 107 108* ****************************** 109* * * 110* * dipole_Efield_set_ion * 111* * * 112* ****************************** 113 subroutine dipole_Efield_set_ion() 114 implicit none 115 116 !**** dipole_efield common block **** 117 real*8 efield(3),efield_center(3) 118 real*8 dx(2),dy(2),dz(2) 119 real*8 nx,ny,nz,pcharge 120 real*8 dipole(3),Edipole,Pdipole 121 real*8 bv(3,6),wts(6) 122 integer rank 123 integer Tc(2),Ts(2),rgrid(2) 124 integer ispin0,ne0(2) 125 common /dipole_efield/ efield,efield_center, 126 > dx,dy,dz, 127 > nx,ny,nz,pcharge,dipole, 128 > Edipole,Pdipole, 129 > bv,wts,rank, 130 > Tc,Ts,rgrid, 131 > ispin0,ne0 132 133* *** local variables *** 134 integer i,ia 135 136* *** external functions *** 137 integer ion_katm,ion_nion 138 external ion_katm,ion_nion 139 real*8 ion_rion,psp_zv 140 external ion_rion,psp_zv 141 142 nx=0.0d0 143 ny=0.0d0 144 nz=0.0d0 145 pcharge = 0.0d0 146 do i=1,ion_nion() 147 ia=ion_katm(i) 148 nx=nx+psp_zv(ia)*ion_rion(1,i) 149 ny=ny+psp_zv(ia)*ion_rion(2,i) 150 nz=nz+psp_zv(ia)*ion_rion(3,i) 151 pcharge = pcharge + psp_zv(ia) 152 end do 153 154 return 155 end 156 157 158* ****************************** 159* * * 160* * dipole_Efield_set_bv * 161* * * 162* ****************************** 163 subroutine dipole_Efield_set_bv() 164 implicit none 165 166 !**** dipole_efield common block **** 167 real*8 efield(3),efield_center(3) 168 real*8 dx(2),dy(2),dz(2) 169 real*8 nx,ny,nz,pcharge 170 real*8 dipole(3),Edipole,Pdipole 171 real*8 bv(3,6),wts(6) 172 integer rank 173 integer Tc(2),Ts(2),rgrid(2) 174 integer ispin0,ne0(2) 175 common /dipole_efield/ efield,efield_center, 176 > dx,dy,dz, 177 > nx,ny,nz,pcharge,dipole, 178 > Edipole,Pdipole, 179 > bv,wts,rank, 180 > Tc,Ts,rgrid, 181 > ispin0,ne0 182 183* *** local variables *** 184 integer tmp_len 185 parameter (tmp_len=140) 186 187 integer i,j,info 188 real*8 wrk(6,6),bmat(3,3),ixmat(3,6) 189 190 real*8 xx,yy,zz,tmp1(tmp_len) 191 real*8 tx,ty,tz,alpha,scal 192 193* *** external functions *** 194 real*8 lattice_unitg,lattice_omega 195 external lattice_unitg,lattice_omega 196 197c *** Silvestrelli G1 *** 198 ixmat(1,1)=1.0d0 199 ixmat(2,1)=0.0d0 200 ixmat(3,1)=0.0d0 201 202c *** Silvestrelli G4 *** 203 ixmat(1,2)=1.0d0 204 ixmat(2,2)=1.0d0 205 ixmat(3,2)=0.0d0 206 207c *** Silvestrelli G5 *** 208 ixmat(1,3)=1.0d0 209 ixmat(2,3)=0.0d0 210 ixmat(3,3)=1.0d0 211 212c *** Silvestrelli G2 *** 213 ixmat(1,4)=0.0d0 214 ixmat(2,4)=1.0d0 215 ixmat(3,4)=0.0d0 216 217c *** Silvestrelli G6 *** 218 ixmat(1,5)=0.0d0 219 ixmat(2,5)=1.0d0 220 ixmat(3,5)=1.0d0 221 222c *** Silvestrelli G3 *** 223 ixmat(1,6)=0.0d0 224 ixmat(2,6)=0.0d0 225 ixmat(3,6)=1.0d0 226 227 do i=1,3 228 bmat(i,1)=lattice_unitg(1,i) 229 bmat(i,2)=lattice_unitg(2,i) 230 bmat(i,3)=lattice_unitg(3,i) 231 end do 232 233 do i=1,6 234 xx=0.0d0 235 yy=0.0d0 236 zz=0.0d0 237 do j=1,3 238 xx=xx+bmat(j,1)*ixmat(j,i) 239 yy=yy+bmat(j,2)*ixmat(j,i) 240 zz=zz+bmat(j,3)*ixmat(j,i) 241 end do 242 bv(1,i)=xx 243 bv(2,i)=yy 244 bv(3,i)=zz 245 end do 246 247 do i=1,6 248 wrk(1,i)=bv(1,i)*bv(1,i) 249 wrk(2,i)=bv(1,i)*bv(2,i) 250 wrk(3,i)=bv(1,i)*bv(3,i) 251 wrk(4,i)=bv(2,i)*bv(2,i) 252 wrk(5,i)=bv(2,i)*bv(3,i) 253 wrk(6,i)=bv(3,i)*bv(3,i) 254 wts(i)=0.0d0 255 end do 256 257* *** scal=(2*pi/L)**2 *** 258 scal = lattice_omega()**(1.0d0/3.0d0) 259 scal = 8.0d0*datan(1.0d0)/scal 260 scal = scal*scal 261 wts(1)=scal 262 wts(4)=scal 263 wts(6)=scal 264 call DGELS('N',6,6,1,wrk,6,wts,6,tmp1,tmp_len,info) 265 if (info.ne.0) then 266 write(*,*)"Illegal argument in call to dgels" 267 call flush(6) 268 end if 269 rank=0 270 do i=1,6 271 if (dabs(wts(i)).gt.1.e-6) then 272 rank=rank+1 273 wrk(1,rank)=bv(1,i) 274 wrk(2,rank)=bv(2,i) 275 wrk(3,rank)=bv(3,i) 276 wrk(4,rank)=wts(i) 277 end if 278 end do 279 do i=1,rank 280 bv(1,i)=wrk(1,i) 281 bv(2,i)=wrk(2,i) 282 bv(3,i)=wrk(3,i) 283 wts(i)=wrk(4,i) 284 end do 285 286 return 287 end 288 289 290 291* ****************************** 292* * * 293* * dipole_Efield_print * 294* * * 295* ****************************** 296 subroutine dipole_Efield_print(iunit) 297 implicit none 298 integer iunit 299 300 !**** dipole_efield common block **** 301 real*8 efield(3),efield_center(3) 302 real*8 dx(2),dy(2),dz(2) 303 real*8 nx,ny,nz,pcharge 304 real*8 dipole(3),Edipole,Pdipole 305 real*8 bv(3,6),wts(6) 306 integer rank 307 integer Tc(2),Ts(2),rgrid(2) 308 integer ispin0,ne0(2) 309 common /dipole_efield/ efield,efield_center, 310 > dx,dy,dz, 311 > nx,ny,nz,pcharge,dipole, 312 > Edipole,Pdipole, 313 > bv,wts,rank, 314 > Tc,Ts,rgrid, 315 > ispin0,ne0 316 317* *** local variables *** 318 real*8 autoDebye 319 parameter (autoDebye=2.5416d0) 320 321 real*8 xx 322 xx = dsqrt(dipole(1)**2 + dipole(2)**2 + dipole(3)**2) 323 324 write(iunit,1771) 325 write(iunit,1772) 'spin up ', 326 > dx(1)/dble(ne0(1)), 327 > dy(1)/dble(ne0(1)), 328 > dz(1)/dble(ne0(1)) 329 if (ne0(ispin0).ne.0) 330 > write(iunit,1772) 'spin down ', 331 > dx(ispin0)/dble(ne0(ispin0)), 332 > dy(ispin0)/dble(ne0(ispin0)), 333 > dz(ispin0)/dble(ne0(ispin0)) 334 write(iunit,1772) 'electronic', 335 > (dx(1)+dx(ispin0))/dble(ne0(1)+ne0(ispin0)), 336 > (dy(1)+dy(ispin0))/dble(ne0(1)+ne0(ispin0)), 337 > (dz(1)+dz(ispin0))/dble(ne0(1)+ne0(ispin0)) 338 write(iunit,1772) 'ionic ', 339 > nx/pcharge, 340 > ny/pcharge, 341 > nz/pcharge 342 write(iunit,1778) 343 write(iunit,1774) dipole 344 write(iunit,1775) xx,xx*autoDebye 345 346 return 347*::::::::::::::::::::::::::: format ::::::::::::::::::::::::::::::::: 348 1771 FORMAT(//'== Center of Charge =='/) 349 1772 FORMAT(A10,' (',F10.4,',',F10.4,',',F10.4,' )') 350 1773 FORMAT(//'== Wannier Crystal Dipole =='/) 351 1774 FORMAT('mu = (',F10.4,',',F10.4,',',F10.4,' ) au') 352 1775 FORMAT('|mu| = ',F10.4,' au, ',F10.4,' Debye') 353 1776 FORMAT(/"ELECTRONIC DIPOLES") 354 1777 FORMAT("DX =",F11.5," DY= ",F11.5," DZ= ",F11.5) 355 1778 FORMAT(//'== Crystal Dipole (Resta) =='/) 356 1780 FORMAT("NUCLEAR DIPOLES") 357 1785 FORMAT("TOTAL DIPOLES") 358 end 359 360* ****************************** 361* * * 362* * dipole_Efield_gen_TcTs * 363* * * 364* ****************************** 365 subroutine dipole_Efield_set_TcTs(i,ms,neq,n2ft3d,psi_r,psi_r2,W) 366 implicit none 367 integer i,ms,neq(2),n2ft3d 368 real*8 psi_r(*),psi_r2(*) 369 complex*16 W(*) 370 371#include "bafdecls.fh" 372#include "errquit.fh" 373 374 !**** dipole_efield common block **** 375 real*8 efield(3),efield_center(3) 376 real*8 dx(2),dy(2),dz(2) 377 real*8 nx,ny,nz,pcharge 378 real*8 dipole(3),Edipole,Pdipole 379 real*8 bv(3,6),wts(6) 380 integer rank 381 integer Tc(2),Ts(2),rgrid(2) 382 integer ispin0,ne0(2) 383 common /dipole_efield/ efield,efield_center, 384 > dx,dy,dz, 385 > nx,ny,nz,pcharge,dipole, 386 > Edipole,Pdipole, 387 > bv,wts,rank, 388 > Tc,Ts,rgrid, 389 > ispin0,ne0 390 391* **** local variables **** 392 logical value 393 integer j,k,n1,n2,n3 394 real*8 br,scal1 395 integer Wc(2),Ws(2) 396 397* **** external functions **** 398 logical Dneall_m_push_get,Dneall_m_pop_stack 399 external Dneall_m_push_get,Dneall_m_pop_stack 400 401 402 call D3dB_nx(1,n1) 403 call D3dB_ny(1,n2) 404 call D3dB_nz(1,n3) 405 scal1 = 1.0d0/dble(n1*n2*n3) 406 407c **** allocate memory from stack **** 408 value = Dneall_m_push_get(ms,Wc) 409 value = value.and.Dneall_m_push_get(ms,Ws) 410 if (.not. value) 411 > call errquit('dipole_Efield_set_TcTs:out of stack memory',0,0) 412 413 414 !*** generate Tc and Ts *** 415!$OMP DO 416 do k=1,n2ft3d 417 br = bv(1,i)*dbl_mb(rgrid(1)+(k-1)*3) 418 > + bv(2,i)*dbl_mb(rgrid(1)+(k-1)*3 + 1) 419 > + bv(3,i)*dbl_mb(rgrid(1)+(k-1)*3 + 2) 420 dbl_mb(Tc(1)+k-1) = cos(br) 421 dbl_mb(Ts(1)+k-1) = -sin(br) 422 end do 423!$OMP END DO 424 425* **** generate W = <psi_r(i)|exp(-i b*r)|psi_r(j)> **** 426 do j=1,neq(ms) 427 call D3dB_rr_Mul(1,dbl_mb(Tc(1)), 428 > psi_r(1+(j-1+(ms-1)*neq(1))*n2ft3d), 429 > psi_r2(1+(j-1+(ms-1)*neq(1))*n2ft3d)) 430 end do 431 call Dneall_ggm_sym_Multiply(ms,psi_r,psi_r2,n2ft3d, 432 > dbl_mb(Wc(1))) 433 434 call Dneall_m_scal(ms,scal1,dbl_mb(Wc(1))) 435 436 do j=1,neq(ms) 437 call D3dB_rr_Mul(1,dbl_mb(Ts(1)), 438 > psi_r(1+(j-1+(ms-1)*neq(1))*n2ft3d), 439 > psi_r2(1+(j-1+(ms-1)*neq(1))*n2ft3d)) 440 end do 441 call Dneall_ggm_sym_Multiply(ms,psi_r,psi_r2,n2ft3d, 442 > dbl_mb(Ws(1))) 443 444 call Dneall_m_scal(ms,scal1,dbl_mb(Ws(1))) 445 call Dneall_mmtow_Cmplx(ms,dbl_mb(Wc(1)),dbl_mb(Ws(1)),W) 446 447 448* **** pop memory *** 449 value = Dneall_m_pop_stack(Ws) 450 value = value.and.Dneall_m_pop_stack(Wc) 451 if (.not. value) 452 > call errquit('dipole_Efield_set_TcTs:popping stack memory',1,0) 453 454 return 455 end 456 457 458c subroutine dipole_Efield_gen_TcTs(Z,n2ft3d,ZTcTs) 459c implicit none 460c complex*16 Z 461c integer n2ft3d 462c real*8 Tc(*),Ts(*) 463c real*8 ZTcTs(*) 464c integer i 465c do i=1,n2ft3d 466c ZTcTs(i) = dble(Z)*Tc(i) + dimag(Z)*Ts(i) 467c end do 468c return 469c end 470 471* ****************************** 472* * * 473* * dipole_Efield_add_dipole * 474* * * 475* ****************************** 476 subroutine dipole_Efield_add_dipole(ms,r,Z) 477 implicit none 478 integer ms,r 479 complex*16 Z(*) 480 481 !**** dipole_efield common block **** 482 real*8 efield(3),efield_center(3) 483 real*8 dx(2),dy(2),dz(2) 484 real*8 nx,ny,nz,pcharge 485 real*8 dipole(3),Edipole,Pdipole 486 real*8 bv(3,6),wts(6) 487 integer rank 488 integer Tc(2),Ts(2),rgrid(2) 489 integer ispin0,ne0(2) 490 common /dipole_efield/ efield,efield_center, 491 > dx,dy,dz, 492 > nx,ny,nz,pcharge,dipole, 493 > Edipole,Pdipole, 494 > bv,wts,rank, 495 > Tc,Ts,rgrid, 496 > ispin0,ne0 497 498* **** local variables *** 499 integer i 500 real*8 scal 501 complex*16 arg 502 503* **** external functions **** 504 real*8 lattice_omega 505 external lattice_omega 506 507 scal = lattice_omega()**(1.0d0/3.0d0) 508 scal = 8.0d0*datan(1.0d0)/scal 509 scal = scal*scal 510 511!$OMP MASTER 512 !*** really just want complex eigenvalues of X here *** 513 do i=1,ne0(ms) 514 arg = -wts(r)*datan2(dimag(Z(i)),dble(Z(i))) 515 516 dx(ms) = dx(ms) + bv(1,r)*arg/scal 517 dy(ms) = dy(ms) + bv(2,r)*arg/scal 518 dz(ms) = dz(ms) + bv(3,r)*arg/scal 519 end do !* i * 520!$OMP END MASTER 521 522 return 523 end 524 525* ****************************** 526* * * 527* * dipole_Efield_gen_Cij * 528* * * 529* ****************************** 530 subroutine dipole_Efield_gen_Cij(N,Z,U,V,ZENOM,Creal,Cimg) 531 implicit none 532 integer N 533 complex*16 Z(N),U(N,N),V(N,N),ZENOM(N) 534 real*8 Creal(N,N),Cimg(N,N) 535 536* **** local variables **** 537 integer i,j,k 538 complex*16 zzz 539 540!$OMP DO 541 do k=1,N 542 ZENOM(k) = dcmplx(0.0d0,0.0d0) 543 do i=1,N 544 ZENOM(k) = ZENOM(k) + dconjg(U(i,k))*V(i,k) 545 end do 546 !write(*,*) "k,ZENOM=",k,ZENOM(k) 547 end do 548!$OMP END DO 549 550!$OMP DO 551 do i=1,N 552 do j=1,N 553 zzz = dcmplx(0.0d0,0.0d0) 554 do k=1,N 555c zzz = zzz 556c > + (dcmplx(-dimag(Z(k)),dble(Z(k))) 557c > /(dble(Z(k))**2 + dimag(Z(k))**2)) 558c > * (dconjg(U(i,k))*V(j,k))/ZENOM(k) 559 560 zzz = zzz 561 > + (dcmplx(dimag(Z(k)),dble(Z(k))) 562 > /(dble(Z(k))**2 + dimag(Z(k))**2)) 563 > * (dconjg(U(i,k))*V(j,k))/ZENOM(k) 564 565 end do 566 Creal(i,j) = dble(zzz) 567 Cimg(i,j) = dimag(zzz) 568 end do 569 end do 570!$OMP END DO 571 572 return 573 end 574 575 576* ****************************** 577* * * 578* * dipole_Efield_gen_hpsi_r * 579* * * 580* ****************************** 581 subroutine dipole_Efield_gen_hpsi_r(ms,r,neq,n2ft3d, 582 > psi_r,psi_r2, 583 > Creal,Cimg,hpsi_r) 584 implicit none 585 integer ms,r,neq(2),n2ft3d 586 real*8 psi_r(*),psi_r2(*) 587 real*8 Creal(*),Cimg(*) 588 real*8 hpsi_r(*) 589 590#include "bafdecls.fh" 591#include "errquit.fh" 592 593 !**** dipole_efield common block **** 594 real*8 efield(3),efield_center(3) 595 real*8 dx(2),dy(2),dz(2) 596 real*8 nx,ny,nz,pcharge 597 real*8 dipole(3),Edipole,Pdipole 598 real*8 bv(3,6),wts(6) 599 integer rank 600 integer Tc(2),Ts(2),rgrid(2) 601 integer ispin0,ne0(2) 602 common /dipole_efield/ efield,efield_center, 603 > dx,dy,dz, 604 > nx,ny,nz,pcharge,dipole, 605 > Edipole,Pdipole, 606 > bv,wts,rank, 607 > Tc,Ts,rgrid, 608 > ispin0,ne0 609 610* **** local variables **** 611 integer j,k 612 real*8 alpha,scal,br 613 614* *** external variables **** 615 real*8 lattice_omega 616 external lattice_omega 617 618 619 scal = lattice_omega()**(1.0d0/3.0d0) 620 scal = 8.0d0*datan(1.0d0)/scal 621 scal = scal*scal 622 623 alpha = wts(r)*( efield(1)*bv(1,r) 624 > + efield(2)*bv(2,r) 625 > + efield(3)*bv(3,r)) 626 > / scal 627 628 !write(*,*) "alpha=",alpha,Creal(1),Cimg(1) 629 630 631 !*** generate Tc and Ts *** 632!$OMP DO 633 do k=1,n2ft3d 634 br = bv(1,r)*dbl_mb(rgrid(1)+(k-1)*3) 635 > + bv(2,r)*dbl_mb(rgrid(1)+(k-1)*3 + 1) 636 > + bv(3,r)*dbl_mb(rgrid(1)+(k-1)*3 + 2) 637 dbl_mb(Tc(1)+k-1) = cos(br) 638 dbl_mb(Ts(1)+k-1) = -sin(br) 639 end do 640!$OMP END DO 641 642 643* **** generate W = <psi_r(i)|exp(-i b*r)|psi_r(j)> **** 644 do j=1,neq(ms) 645 call D3dB_rr_Mul(1,dbl_mb(Tc(1)), 646 > psi_r(1+(j-1+(ms-1)*neq(1))*n2ft3d), 647 > psi_r2(1+(j-1+(ms-1)*neq(1))*n2ft3d)) 648 end do 649 call Dneall_gmg_Multiply(ms,psi_r2,n2ft3d, 650 > Creal,-alpha, 651 > hpsi_r,1.0d0) 652 653 do j=1,neq(ms) 654 call D3dB_rr_Mul(1,dbl_mb(Ts(1)), 655 > psi_r(1+(j-1+(ms-1)*neq(1))*n2ft3d), 656 > psi_r2(1+(j-1+(ms-1)*neq(1))*n2ft3d)) 657 end do 658 call Dneall_gmg_Multiply(ms,psi_r2,n2ft3d, 659 > Cimg,alpha, 660 > hpsi_r,1.0d0) 661 662 663 return 664 end 665 666 667* ****************************** 668* * * 669* * dipole_Efield_Vnl * 670* * * 671* ****************************** 672 subroutine dipole_Efield_Vnl(ispin,neq,n2ft3d,psi_r,hpsi_r,edpol) 673 implicit none 674 integer ispin,neq(2),n2ft3d 675 real*8 psi_r(n2ft3d,*) 676 real*8 hpsi_r(n2ft3d,*) 677 real*8 edpol 678 679#include "bafdecls.fh" 680#include "errquit.fh" 681#include "util.fh" 682#include "stdio.fh" 683 684 !**** dipole_efield common block **** 685 real*8 efield(3),efield_center(3) 686 real*8 dx(2),dy(2),dz(2) 687 real*8 nx,ny,nz,pcharge 688 real*8 dipole(3),Edipole,Pdipole 689 real*8 bv(3,6),wts(6) 690 integer rank 691 integer Tc(2),Ts(2),rgrid(2) 692 integer ispin0,ne0(2) 693 common /dipole_efield/ efield,efield_center, 694 > dx,dy,dz, 695 > nx,ny,nz,pcharge,dipole, 696 > Edipole,Pdipole, 697 > bv,wts,rank, 698 > Tc,Ts,rgrid, 699 > ispin0,ne0 700 701* **** local variables **** 702 logical value 703 integer i,j,k,r,ms 704 integer n1,n2,n3 705 706 real*8 b(3) 707 real*8 xx,yy,zz,dv 708 real*8 tx,ty,tz,scal,scal1 709 complex*16 arg,wx 710 711 integer X(2,6),Xeig(2,6) 712 integer U(2,6),V(2,6),UV(2) 713 integer Creal(2),Cimg(2) 714 integer psi_r2(2) 715 716* **** external functions **** 717 logical Dneall_w_push_get,Dneall_w_pop_stack,control_print 718 external Dneall_w_push_get,Dneall_w_pop_stack,control_print 719 logical Dneall_m_push_get,Dneall_m_pop_stack 720 external Dneall_m_push_get,Dneall_m_pop_stack 721 real*8 lattice_omega,Dneall_m_trace 722 external lattice_omega,Dneall_m_trace 723 724 725 call D3dB_nx(1,n1) 726 call D3dB_ny(1,n2) 727 call D3dB_nz(1,n3) 728 scal1 = 1.0d0/dble(n1*n2*n3) 729 730* *** scal=(2*pi/L)**2 *** 731 scal = lattice_omega()**(1.0d0/3.0d0) 732 scal = 8.0d0*datan(1.0d0)/scal 733 scal = scal*scal 734 735* ***** allocate X,Y,Z **** 736 value = .true. 737 do j=1,6 738 value = value.and.Dneall_w_push_get(1,X(1,j)) 739 value = value.and. 740 > BA_push_get(mt_dcpl,ne0(1),'Xeig',Xeig(2,j),Xeig(1,j)) 741 value = value.and.Dneall_w_push_get(1,U(1,j)) 742 value = value.and.Dneall_w_push_get(1,V(1,j)) 743 end do 744 value = value.and.BA_push_get(mt_dcpl,ne0(1),'UV',UV(2),UV(1)) 745 value = value.and.Dneall_m_push_get(0,Creal) 746 value = value.and.Dneall_m_push_get(1,Cimg) 747 value = value.and. 748 > BA_push_get(mt_dbl,(neq(1)+neq(2))*n2ft3d, 749 > 'Hpsi_r',psi_r2(2),psi_r2(1)) 750 if (.not.value) 751 > call errquit('dipole_efield_Vnl:out of stack',0,MA_ERR) 752 753 754!$OMP MASTER 755 dx(1)=0.0d0 756 dy(1)=0.0d0 757 dz(1)=0.0d0 758 dx(2)=0.0d0 759 dy(2)=0.0d0 760 dz(2)=0.0d0 761!$OMP END MASTER 762 763 764 do ms=1,ispin 765 do r=1,rank 766 call dipole_Efield_set_TcTs(r,ms,neq,n2ft3d, 767 > psi_r,dbl_mb(psi_r2(1)),dcpl_mb(X(1,r))) 768 769 call Dneall_w_eigenvaluesvectors(ms,dcpl_mb(X(1,r)), 770 > dcpl_mb(Xeig(1,r)), 771 > dcpl_mb(U(1,r)), 772 > dcpl_mb(V(1,r))) 773 774 call dipole_Efield_add_dipole(ms,r,dcpl_mb(Xeig(1,r))) 775 call dipole_Efield_gen_Cij(ne0(ms), 776 > dcpl_mb(Xeig(1,r)), 777 > dcpl_mb(U(1,r)), 778 > dcpl_mb(V(1,r)), 779 > dcpl_mb(UV(1)), 780 > dbl_mb(Creal(1)),dbl_mb(Cimg(1))) 781 call dipole_Efield_gen_hpsi_r(ms,r,neq,n2ft3d, 782 > psi_r, 783 > dbl_mb(psi_r2(1)), 784 > dbl_mb(Creal(1)),dbl_mb(Cimg(1)), 785 > hpsi_r) 786 end do !* r * 787 end do !* ms * 788 789 call Dneall_ggm_sym_Multiply(0,psi_r,hpsi_r,n2ft3d, 790 > dbl_mb(Creal(1))) 791!$OMP MASTER 792 Pdipole = scal1*Dneall_m_trace(0,dbl_mb(Creal(1))) 793 794 795cccccccccccccccccccccccccccccccccccccccccccccc 796c Molecular dipoles from Resta's theory! 797ccccccccccccccccccccccccccccccccccccccccccccc 798 tx=nx-dx(1)-dx(ispin) 799 ty=ny-dy(1)-dy(ispin) 800 tz=nz-dz(1)-dz(ispin) 801 xx = dsqrt(tx*tx + ty*ty + tz*tz) 802 dipole(1) = tx 803 dipole(2) = ty 804 dipole(3) = tz 805!$OMP END MASTER 806!$OMP BARRIER 807 808 !write(*,*) "DipOlE=",dipole 809 edpol = dipole(1) * efield(1) 810 > + dipole(2) * efield(2) 811 > + dipole(3) * efield(3) 812 813 Edipole = edpol 814 815 816c **** pop stack **** 817 value = BA_pop_stack(psi_r2(2)) 818 value = value.and.Dneall_m_pop_stack(Cimg) 819 value = value.and.Dneall_m_pop_stack(Creal) 820 value = value.and.BA_pop_stack(UV(2)) 821 do j=6,1,-1 822 value = value.and.Dneall_w_pop_stack(V(1,j)) 823 value = value.and.Dneall_w_pop_stack(U(1,j)) 824 value = value.and.BA_pop_stack(Xeig(2,j)) 825 value = value.and.Dneall_w_pop_stack(X(1,j)) 826 end do 827 if (.not.value) 828 > call errquit('dipole_efield_vnl:pop stack',2,MA_ERR) 829 830 831 !write(*,*) "ispin,r,Pdipole,Edipole = ", ispin,r,Pdipole,Edipole 832 !write(*,*) 833 834 return 835 end 836 837 838* ****************************** 839* * * 840* * dipole_Efield_e * 841* * * 842* ****************************** 843 real*8 function dipole_Efield_e() 844 implicit none 845 846 !**** dipole_efield common block **** 847 real*8 efield(3),efield_center(3) 848 real*8 dx(2),dy(2),dz(2) 849 real*8 nx,ny,nz,pcharge 850 real*8 dipole(3),Edipole,Pdipole 851 real*8 bv(3,6),wts(6) 852 integer rank 853 integer Tc(2),Ts(2),rgrid(2) 854 integer ispin0,ne0(2) 855 common /dipole_efield/ efield,efield_center, 856 > dx,dy,dz, 857 > nx,ny,nz,pcharge,dipole, 858 > Edipole,Pdipole, 859 > bv,wts,rank, 860 > Tc,Ts,rgrid, 861 > ispin0,ne0 862 863 dipole_efield_e = Edipole 864 return 865 end 866 867 868* ****************************** 869* * * 870* * dipole_Efield_p * 871* * * 872* ****************************** 873 real*8 function dipole_Efield_p() 874 implicit none 875 876 !**** dipole_efield common block **** 877 real*8 efield(3),efield_center(3) 878 real*8 dx(2),dy(2),dz(2) 879 real*8 nx,ny,nz,pcharge 880 real*8 dipole(3),Edipole,Pdipole 881 real*8 bv(3,6),wts(6) 882 integer rank 883 integer Tc(2),Ts(2),rgrid(2) 884 integer ispin0,ne0(2) 885 common /dipole_efield/ efield,efield_center, 886 > dx,dy,dz, 887 > nx,ny,nz,pcharge,dipole, 888 > Edipole,Pdipole, 889 > bv,wts,rank, 890 > Tc,Ts,rgrid, 891 > ispin0,ne0 892 893 dipole_efield_p = Pdipole 894 return 895 end 896 897 898 899 900