TABLE OF CONTENTS
- m_oscillators/calc_wfwfg
- m_oscillators/gw_box2gsph
- m_oscillators/rho_tw_g
- m_oscillators/rotate_spinor
- m_oscillators/sym_rhotwgq0
- m_oscillators/ts_usug_kkp_bz
- m_oscillators/usur_kkp_bz
m_oscillators/calc_wfwfg [ Functions ]
NAME
calc_wfwfg
FUNCTION
Calculate the Fourier transform of the product u_{bk}^*(r) u_{b"k}(r)
Return values on the FFT box.
INPUTS
nspinor=number of spinorial components. spinrot(4)=components of the spinor rotation matrix
OUTPUT
SOURCE
555 subroutine calc_wfwfg(ktabr_k, ktabi_k, spinrot, nr, nspinor, ngfft_gw, wfr_jb, wfr_kb, wfg2_jk) 556 557 !Arguments ------------------------------------ 558 !scalars 559 integer,intent(in) :: ktabi_k,nr,nspinor 560 !arrays 561 integer,intent(in) :: ktabr_k(nr),ngfft_gw(18) 562 real(dp),intent(in) :: spinrot(4) 563 complex(gwpc),intent(in) :: wfr_jb(nr*nspinor),wfr_kb(nr*nspinor) 564 complex(gwpc),intent(out) :: wfg2_jk(nr*nspinor) 565 566 !Local variables------------------------------- 567 integer,parameter :: ndat1 = 1, fftcache0 = 0, gpu_option_0 = 0 568 type(fftbox_plan3_t) :: plan 569 !arrays 570 complex(gwpc),allocatable :: wfr2_dpcplx(:),ujb_bz(:),ukb_bz(:) 571 572 ! ************************************************************************* 573 574 ! There is no need to take into account phases arising from non-symmorphic 575 ! operations since the wavefunctions are evaluated at the same k-point. 576 ABI_MALLOC(wfr2_dpcplx, (nr * nspinor * ndat1)) 577 578 if (nspinor == 1) then 579 select case (ktabi_k) 580 case (1) 581 wfr2_dpcplx = GWPC_CONJG(wfr_jb(ktabr_k)) * wfr_kb(ktabr_k) 582 case (2) 583 ! Conjugate the product if time-reversal is used to reconstruct this k-point 584 wfr2_dpcplx = wfr_jb(ktabr_k) * GWPC_CONJG(wfr_kb(ktabr_k)) 585 case default 586 ABI_ERROR(sjoin("Wrong ktabi_k:", itoa(ktabi_k))) 587 end select 588 589 else if (nspinor == 2) then 590 ABI_MALLOC(ujb_bz, (nr * nspinor * ndat1)) 591 ABI_MALLOC(ukb_bz, (nr * nspinor * ndat1)) 592 ! Use wfr2_dpcplx as workspace array 593 call rotate_spinor(ktabi_k, ktabr_k, cone, spinrot, nr, nspinor, ndat1, wfr_jb, wfr2_dpcplx, ujb_bz) 594 call rotate_spinor(ktabi_k, ktabr_k, cone, spinrot, nr, nspinor, ndat1, wfr_kb, wfr2_dpcplx, ukb_bz) 595 wfr2_dpcplx = GWPC_CONJG(ujb_bz) * ukb_bz 596 ABI_FREE(ujb_bz) 597 ABI_FREE(ukb_bz) 598 599 else 600 ABI_ERROR(sjoin("Wrong nspinor:", itoa(nspinor))) 601 end if 602 603 ! Transform to Fourier space (result in wfg2_jk) 604 call plan%init(nspinor, ngfft_gw(1:3), ngfft_gw(1:3), ngfft_gw(7), fftcache0, gpu_option_0) 605 call plan%execute(wfr2_dpcplx, wfg2_jk, -1) 606 call plan%free() 607 ABI_FREE(wfr2_dpcplx) 608 609 end subroutine calc_wfwfg
m_oscillators/gw_box2gsph [ Functions ]
NAME
gw_box2gsph
FUNCTION
Trasnfer data from the FFT box to the G-sphere.
INPUTS
nr=number of FFT grid points ndat=Number of wavefunctions to transform. npw=number of plane waves in the sphere igfftg0(npw)=index of G-G_o in the FFT array for each G in the sphere. iarrbox(nr*ndat)=Input array on the FFT mesh
OUTPUT
oarrsph(npw*ndat)=output array on the sphere.
SOURCE
489 subroutine gw_box2gsph(nr, ndat, npw, igfftg0, iarrbox, oarrsph) 490 491 !Arguments ------------------------------------ 492 !scalars 493 integer,intent(in) :: nr,ndat,npw 494 !arrays 495 integer,intent(in) :: igfftg0(npw) 496 complex(gwpc),intent(in) :: iarrbox(nr*ndat) 497 complex(gwpc),intent(out) :: oarrsph(npw*ndat) 498 499 !Local variables------------------------------- 500 !scalars 501 integer :: ig,igfft,dat,pgsp,pfft 502 503 ! ************************************************************************* 504 505 if (ndat==1) then 506 do ig=1,npw 507 igfft=igfftg0(ig) 508 if (igfft/=0) then 509 ! G-G0 belongs to the FFT mesh. 510 oarrsph(ig) = iarrbox(igfft) 511 else 512 ! Set this component to zero. 513 oarrsph(ig) = czero_gw 514 end if 515 end do 516 else 517 !$OMP PARALLEL DO PRIVATE(pgsp,pfft,igfft) 518 do dat=1,ndat 519 pgsp = (dat-1)*npw 520 pfft = (dat-1)*nr 521 do ig=1,npw 522 igfft=igfftg0(ig) 523 if (igfft/=0) then 524 ! G-G0 belongs to the FFT mesh. 525 oarrsph(ig+pgsp) = iarrbox(igfft+pfft) 526 else 527 ! Set this component to zero. 528 oarrsph(ig+pgsp) = czero_gw 529 end if 530 end do 531 end do 532 end if 533 534 end subroutine gw_box2gsph
m_oscillators/rho_tw_g [ Functions ]
NAME
rho_tw_g
FUNCTION
Calculate rhotwg(G) = <wfn1| exp(-i(q+G).r) |wfn2>
INPUTS
dim_rtwg=Define the size of the output array rhotwg
=== for nspinor==1 ===
dim_rtwg=1
=== for nspinor==2 ===
dim_rtwg=1 if the sum of the matrix elements is wanted.
dim_rtwg=2 if <up|up>, <dwn|dwn> matrix elements are required
map2sphere= 1 to retrieve Fourier components indexed according to igfftg0.
0 to retrieve Fourier components indexed according to the FFT box.
NOTE: If map2sphere==0 npwvec must be equal to nr
use_padfft= Only compatible with map2sphere 1.
1 if matrix elements are calculated via zero-padded FFT.
0 R-->G Transform in done on the full FFT box.
igfftg0(npwvec*map2sphere)=index of G-G_o in the FFT array for each G in the sphere.
i1=1 if kbz1 = Sk1, 2 if kbz1 = -Sk_1 (k_1 is in the IBZ)
i2=1 if kbz2 = Sk2, 2 if kbz2 = -Sk_2 (k_2 is in the IBZ)
ktabr1(nr),ktabr2(nr)= tables R^-1(r-t) for the two k-points
ktabp1,ktabp2 = phase factors for non-simmorphic symmetries e^{-i 2\pi kbz.\tau}
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
npwvec=number of plane waves (in the sphere if map2sphere==1, in the FFT box if map2sphere==1)
nr=number of FFT grid points
ndat=Number of wavefunctions to transform.
nspinor=number of spinorial components.
spinrot1(4),spinrot2(4)=components of the spinor rotation matrix. See getspinrot
wfn1(nr*nspinor*ndat),wfn2(nr*nspinor*ndat)=the two wavefunctions (periodic part)
[nhat12(2,nr,nspinor**2*ndat)]=Compensation charge in real space to be added to \Psi_1^*\Psi_2 -- Only for PAW.
OUTPUT
rhotwg(npwvec)=density of a pair of states, in reciprocal space
SOURCE
90 subroutine rho_tw_g(nspinor, npwvec, nr, ndat, ngfft, map2sphere, use_padfft, igfftg0, gbound, & 91 wfn1, i1, ktabr1, ktabp1, spinrot1, & 92 wfn2, i2, ktabr2, ktabp2, spinrot2, & 93 dim_rtwg, rhotwg) !& nhat12) 94 95 !Arguments ------------------------------------ 96 !scalars 97 integer,intent(in) :: i1,i2,npwvec,nr,nspinor,dim_rtwg,map2sphere,use_padfft,ndat 98 complex(dpc),intent(in) :: ktabp1, ktabp2 99 !arrays 100 integer,intent(in) :: gbound(:,:) !gbound(2*mgfft+8,2) 101 integer,intent(in) :: igfftg0(npwvec*map2sphere),ngfft(18) 102 integer,intent(in) :: ktabr1(nr),ktabr2(nr) 103 real(dp),intent(in) :: spinrot1(4),spinrot2(4) 104 complex(gwpc),intent(in) :: wfn1(nr*nspinor*ndat),wfn2(nr*nspinor*ndat) 105 complex(gwpc),intent(out) :: rhotwg(npwvec*dim_rtwg*ndat) 106 ! real(dp),optional,intent(in) :: nhat12(2,nr,nspinor**2*ndat) 107 108 !Local variables------------------------------- 109 !scalars 110 integer :: fftcache0 = 0, gpu_option_0 = 0 111 integer :: ig,igfft,iab,spad1,spad2,spad0,nx,ny,nz,ldx,ldy,ldz,mgfft 112 type(fftbox_plan3_t) :: plan 113 !arrays 114 integer :: spinor_pad(2,4) 115 complex(gwpc),allocatable :: u12prod(:),cwavef1(:),cwavef2(:),cwork(:) 116 117 ! ************************************************************************* 118 119 SELECT CASE (nspinor) 120 CASE (1) 121 ! Collinear case. 122 call ts_usug_kkp_bz(npwvec,nr,ndat,ngfft,map2sphere,use_padfft,igfftg0,gbound,& 123 wfn1,i1,ktabr1,ktabp1,& 124 wfn2,i2,ktabr2,ktabp2,rhotwg) 125 126 CASE (2) 127 ! Spinorial case. 128 ABI_CHECK(ndat==1,"ndat != 1 not coded") 129 ABI_MALLOC(cwavef1, (nr * nspinor * ndat)) 130 ABI_MALLOC(cwavef2, (nr * nspinor * ndat)) 131 ABI_MALLOC(cwork, (nr * nspinor * ndat)) 132 133 call rotate_spinor(i1, ktabr1, ktabp1, spinrot1, nr, nspinor, ndat, wfn1, cwork, cwavef1) 134 call rotate_spinor(i2, ktabr2, ktabp2, spinrot2, nr, nspinor, ndat, wfn2, cwork, cwavef2) 135 136 ABI_FREE(cwork) 137 ABI_MALLOC(u12prod, (nr)) 138 139 spinor_pad = reshape([0, 0, nr, nr, 0, nr, nr, 0], [2, 4]) 140 rhotwg = czero_gw 141 do iab=1,2 142 spad1 = spinor_pad(1,iab); spad2=spinor_pad(2,iab) 143 144 u12prod = GWPC_CONJG(cwavef1(spad1+1:spad1+nr)) * cwavef2(spad2+1:spad2+nr) 145 ! Add compensation charge. 146 !if (PRESENT(nhat12)) u12prod = u12prod + CMPLX(nhat12(1,:,iab),nhat12(2,:,iab)) 147 148 spad0 = (iab-1)*npwvec 149 SELECT CASE (map2sphere) 150 CASE (0) 151 ! Need results on the full FFT box thus cannot use zero-padded FFT. 152 call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0) 153 call plan%execute(u12prod, -1) 154 call plan%free() 155 if (dim_rtwg == 1) then 156 rhotwg(1:npwvec) = rhotwg(1:npwvec) + u12prod 157 else 158 rhotwg(spad0+1:spad0+npwvec) = u12prod 159 end if 160 161 CASE (1) 162 ! Need results on the G-sphere. Call zero-padded FFT routines if required. 163 if (use_padfft == 1) then 164 nx = ngfft(1); ny = ngfft(2); nz = ngfft(3); mgfft = maxval(ngfft(1:3)) 165 ldx = nx; ldy = ny; ldz = nz 166 call fftpad(u12prod, ngfft, nx, ny, nz, ldx, ldy, ldz, ndat, mgfft, -1, gbound) 167 else 168 call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0) 169 call plan%execute(u12prod, -1) 170 call plan%free() 171 end if 172 173 ! Have to map FFT to G-sphere. 174 if (dim_rtwg == 1) then 175 do ig=1,npwvec 176 igfft = igfftg0(ig) 177 ! G-G0 belong to the FFT mesh. 178 if (igfft /= 0) rhotwg(ig) = rhotwg(ig) + u12prod(igfft) 179 end do 180 else 181 do ig=1,npwvec 182 igfft = igfftg0(ig) 183 ! G-G0 belong to the FFT mesh. 184 if (igfft /= 0) rhotwg(ig+spad0) = u12prod(igfft) 185 end do 186 end if 187 188 CASE DEFAULT 189 ABI_BUG("Wrong map2sphere") 190 END SELECT 191 end do !iab 192 193 ABI_FREE(u12prod) 194 ABI_FREE(cwavef1) 195 ABI_FREE(cwavef2) 196 197 CASE DEFAULT 198 ABI_BUG('Wrong nspinor') 199 END SELECT 200 201 end subroutine rho_tw_g
m_oscillators/rotate_spinor [ Functions ]
NAME
rotate_spinor
FUNCTION
Return a spinor in the full BZ from its symmetrical image in the IBZ.
INPUTS
OUTPUT
SOURCE
714 subroutine rotate_spinor(itim_kbz, ktabr_kbz, ktabp_kbz, spinrot, nr, nspinor, ndat, ug_ibz, cwork, oug_bz) 715 716 !Arguments ------------------------------------ 717 !scalars 718 integer,intent(in) :: itim_kbz, nr, nspinor, ndat 719 complex(dpc),intent(in) :: ktabp_kbz 720 !arrays 721 integer,intent(in) :: ktabr_kbz(nr) 722 real(dp),intent(in) :: spinrot(4) 723 complex(gwpc),intent(in) :: ug_ibz(nr*nspinor*ndat) 724 complex(gwpc),intent(out) :: cwork(nr*nspinor*ndat), oug_bz(nr*nspinor*ndat) 725 726 !Local variables ------------------------------ 727 !scalars 728 integer :: ir,ir1,spad0,ispinor 729 complex(gwpc) :: u1a,u1b 730 !arrays 731 complex(dpc) :: spinrot_cmat1(2,2) 732 733 !************************************************************************ 734 735 ABI_CHECK(ndat == 1, "ndat > 1 not coded") 736 ABI_CHECK(nspinor == 2, "nspinor should be 1") 737 738 ! Apply Time-reversal if required. 739 ! \psi_{-k}^1 = (\psi_k^2)^* 740 ! \psi_{-k}^2 = -(\psi_k^1)^* 741 if (itim_kbz == 1) then 742 cwork(:) = ug_ibz(:) 743 else if (itim_kbz == 2) then 744 cwork(1:nr) = GWPC_CONJG(ug_ibz(nr+1:2*nr)) 745 cwork(nr+1:2*nr) = -GWPC_CONJG(ug_ibz(1:nr)) 746 else 747 ABI_ERROR(sjoin('Wrong itim_kbz:', itoa(itim_kbz))) 748 end if 749 750 do ispinor=1,nspinor 751 spad0 = (ispinor-1) * nr 752 do ir=1,nr 753 ir1 = ktabr_kbz(ir) 754 oug_bz(ir+spad0) = cwork(ir1+spad0) * ktabp_kbz 755 end do 756 end do 757 758 ! Rotation in spinor space (same equations as in wfconv) 759 spinrot_cmat1 = spinrot_cmat(spinrot) 760 cwork = oug_bz 761 do ir=1,nr 762 u1a = cwork(ir); u1b = cwork(ir+nr) 763 oug_bz(ir) = spinrot_cmat1(1, 1) * u1a + spinrot_cmat1(1, 2) * u1b 764 oug_bz(ir+nr) = spinrot_cmat1(2, 1) * u1a + spinrot_cmat1(2, 2) * u1b 765 end do 766 767 end subroutine rotate_spinor
m_oscillators/sym_rhotwgq0 [ Functions ]
NAME
sym_rhotwgq0
FUNCTION
Symmetrization of the oscillator matrix elements <k-q,b1|exp(-i(q+G).r)|k,b2> in the special case of q=0. The matrix elements in the full BZ is obtained from the matrix elements in the IBZ by rotating the wavefunctions and taking into account time reversal symmetry. strictly speaking the symmetrization can be performed only for non-degenerate states.
INPUTS
Gsph<gsphere_t>=Info on the G-sphere used to describe wavefunctions and W (the largest one is actually stored). npw=Number of G-vectors dim_rtwg=Number of spin-spin combinations, 1 for collinear spin, 4 is nspinor==2 (TODO NOT CODED) itim_k=2 if time reversal is used to reconstruct the k in the BZ, 1 otherwise. isym_k=The index of the symmetry symrec rotains k_IBZ onto k_BZ. rhxtwg_in(dim_rtwg*npw)=The input matrix elements in the IBZ.
OUTPUT
rhxtwg_sym(dim_rtwg*npw)=The symmetrized matrix elements in the BZ.
NOTES
Let M_{G}(k,q) =<k-q,b1|exp(-i(q+G).r)|k,b2>
At q ==0, supposing non-degenerate bands, one obtains:
1) M_{ SG}( Sk) = e^{-iSG.t} M_{G} (k)
2) M_{-SG}(-Sk) = e^{+iSG.t} M_{G}^* (k)
SOURCE
644 function sym_rhotwgq0(itim_k, isym_k, dim_rtwg, npw, rhxtwg_in, Gsph) result(rhxtwg_sym) 645 646 !Arguments ------------------------------------ 647 !scalars 648 integer,intent(in) :: npw,dim_rtwg,itim_k,isym_k 649 type(gsphere_t),intent(in) :: Gsph 650 !arrays 651 complex(gwpc),intent(in) :: rhxtwg_in(dim_rtwg*npw) 652 complex(gwpc) :: rhxtwg_sym(dim_rtwg*npw) 653 654 !Local variables ------------------------------ 655 !scalars 656 integer :: ig 657 658 !************************************************************************ 659 660 ABI_CHECK(dim_rtwg == 1, "dim_rtwg/=1 not coded") 661 662 SELECT CASE (isym_k) 663 CASE (1) 664 ! Fractional translation associated to E is assumed to be (zero,zero,zero). 665 SELECT CASE (itim_k) 666 CASE (1) 667 ! Identity, no time-reversal. No symmetrization is needed. 668 rhxtwg_sym(:) = rhxtwg_in(:) 669 CASE (2) 670 ! Identity + Time-reversal. 671 do ig=1,npw 672 rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = GWPC_CONJG(rhxtwg_in(ig)) 673 end do 674 CASE DEFAULT 675 ABI_ERROR(sjoin("Wrong value of itim_k:", itoa(itim_k))) 676 END SELECT 677 678 CASE DEFAULT 679 ! Rotate wavefunctions. 680 SELECT CASE (itim_k) 681 CASE (1) 682 ! no time-reversal, only rotation. 683 do ig=1,npw 684 rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = rhxtwg_in(ig) * Gsph%phmSGt(ig,isym_k) 685 end do 686 CASE (2) 687 ! time-reversal + spatial rotation. 688 do ig=1,npw 689 rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = GWPC_CONJG( rhxtwg_in(ig) * Gsph%phmSGt(ig,isym_k) ) 690 end do 691 CASE DEFAULT 692 ABI_ERROR(sjoin("Wrong value of itim_k:", itoa(itim_k))) 693 END SELECT 694 END SELECT 695 696 end function sym_rhotwgq0
m_oscillators/ts_usug_kkp_bz [ Functions ]
NAME
ts_usug_kkp_bz
FUNCTION
Calculate usug(G) = <u1|exp(-i(q+G).r)|u2> for ndat pair of wavefunctions TODO: The routine is thread-safe hence it can be called within an OMP parallel region.
INPUTS
map2sphere= 1 to retrieve Fourier components indexed according to igfftg0.
0 to retrieve Fourier components indexed according to the FFT box.
NOTE: If map2sphere==0 npw must be equal to nr
use_padfft= Only compatible with map2sphere 1.
1 if matrix elements are calculated via zero-padded FFT.
0 R-->G Transform in done on the full FFT box.
igfftg0(npw*map2sphere)=index of G-G_o in the FFT array for each G in the sphere.
time1=1 if kbz1 = Sk1, 2 if kbz1 = -Sk_1 (k_1 is in the IBZ)
time2=1 if kbz2 = Sk2, 2 if kbz2 = -Sk_2 (k_2 is in the IBZ)
ktabr1(nr),ktabr2(nr)= tables R^-1(r-t) for the two k-points
ktabp1,ktabp2 = phase factors for non-simmorphic symmetries e^{-i 2\pi kbz.\tau}
ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
npw=number of plane waves (in the sphere if map2sphere==1, in the FFT box if map2sphere==1)
nr=number of FFT grid points
ndat=Number of wavefunctions to transform.
u1(nr*ndat),u2(nr*ndat)=the two wavefunctions (periodic part)
OUTPUT
usug(npw*ndat)=density of a pair of states, in reciprocal space
SOURCE
237 subroutine ts_usug_kkp_bz(npw, nr, ndat, ngfft, map2sphere, use_padfft, igfftg0, gbound, & 238 u1, time1, ktabr1, ktabp1, & 239 u2, time2, ktabr2, ktabp2, usug) !& nhat12) 240 241 !Arguments ------------------------------------ 242 !scalars 243 integer,intent(in) :: time1,time2,npw,nr,map2sphere,use_padfft,ndat 244 complex(dpc),intent(in) :: ktabp1,ktabp2 245 !arrays 246 integer,intent(in) :: gbound(:,:) !gbound(2*mgfft+8,2) 247 integer,intent(in) :: igfftg0(npw*map2sphere),ngfft(18) 248 integer,intent(in) :: ktabr1(nr),ktabr2(nr) 249 complex(gwpc),intent(in) :: u1(nr*ndat),u2(nr*ndat) 250 complex(gwpc),intent(out) :: usug(npw*ndat) 251 252 !Local variables------------------------------- 253 !scalars 254 integer :: fftcache0 = 0, gpu_option_0 = 0 255 integer :: nx,ny,nz,ldx,ldy,ldz,mgfft 256 type(fftbox_plan3_t) :: plan 257 !arrays 258 complex(gwpc),allocatable :: u12prod(:) 259 260 ! ************************************************************************* 261 262 ! Form rho-twiddle(r) = u_1^*(r,b1,kbz1) u_2(r,b2,kbz2), to account for symmetries: 263 ! 264 ! u(r,b,kbz) = e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t), b, kibz) 265 ! = e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t), b, kibz) for time-reversal symmetry. 266 ! 267 ABI_MALLOC(u12prod,(nr*ndat)) 268 call usur_kkp_bz(nr,ndat,time1,ktabr1,ktabp1,u1,time2,ktabr2,ktabp2,u2,u12prod) 269 270 ! Add compensation charge. 271 !if (PRESENT(nhat12)) u12prod = u1prod + CMPLX(nhat12(1,:,1),nhat12(2,:,1)) 272 273 SELECT CASE (map2sphere) 274 CASE (0) 275 ! Need results on the full FFT box thus cannot use zero-padded FFT. 276 call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0) 277 call plan%execute(u12prod, -1) 278 call plan%free() 279 call xcopy(nr*ndat,u12prod,1,usug,1) 280 281 CASE (1) 282 ! Need results on the G-sphere. Call zero-padded FFT routines if required. 283 if (use_padfft==1) then 284 nx = ngfft(1); ny = ngfft(2); nz = ngfft(3); mgfft = MAXVAL(ngfft(1:3)) 285 ldx=nx; ldy=ny; ldz=nz 286 call fftpad(u12prod,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,-1,gbound) 287 else 288 call plan%init(ndat, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0) 289 call plan%execute(u12prod, -1) 290 call plan%free() 291 end if 292 293 ! From the FFT to the G-sphere. 294 call gw_box2gsph(nr,ndat,npw,igfftg0,u12prod,usug) 295 296 CASE DEFAULT 297 ABI_BUG("Wrong map2sphere") 298 END SELECT 299 300 ABI_FREE(u12prod) 301 302 end subroutine ts_usug_kkp_bz
m_oscillators/usur_kkp_bz [ Functions ]
NAME
usur_kkp_bz
FUNCTION
Calculate u1_kbz^*(r) u2_kbz(r) in real space from the symmetric images in the IBZ. Does not support spinor wavefunctions.
INPUTS
nr=number of FFT grid points ndat=Number of wavefunctions to transform. u1(nr*ndat),u2(nr*ndat)=the two wavefunctions in the IBZ (periodic part) time1=1 if kbz1 = Sk1, 2 if kbz1 = -Sk_1 (k_1 is in the IBZ) time2=1 if kbz2 = Sk2, 2 if kbz2 = -Sk_2 (k_2 is in the IBZ) ktabr1(nr),ktabr2(nr)= tables R^-1(r-t) for the two k-points ktabp1,ktabp2 = phase factors for non-simmorphic symmetries e^{-i 2\pi kbz.\tau}
OUTPUT
u12prod(nr*dat) = u1_kbz^*(r) u2_kbz(r) for the ndat pairs.
SOURCE
329 subroutine usur_kkp_bz(nr, ndat, time1, ktabr1, ktabp1, u1, time2, ktabr2, ktabp2, u2, u12prod) 330 331 !Arguments ------------------------------------ 332 !scalars 333 integer,intent(in) :: nr,ndat,time1,time2 334 complex(dpc),intent(in) :: ktabp1,ktabp2 335 !arrays 336 integer,intent(in) :: ktabr1(nr),ktabr2(nr) 337 complex(gwpc),intent(in) :: u1(nr*ndat),u2(nr*ndat) 338 complex(gwpc),intent(out) :: u12prod(nr*ndat) 339 340 !Local variables------------------------------- 341 !scalars 342 integer :: ir,dat,padat 343 complex(gwpc) :: my_ktabp1,my_ktabp2 344 !arrays 345 complex(gwpc),allocatable :: u1_bz(:),u2_bz(:) 346 347 ! ************************************************************************* 348 349 ! Form rho-twiddle(r)=u_1^*(r,b1,kbz1) u_2(r,b2,kbz2), to account for symmetries: 350 ! u(r,b,kbz)=e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t),b,kibz) 351 ! =e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t),b,kibz) for time-reversal 352 ! 353 ABI_MALLOC(u1_bz,(nr*ndat)) 354 ABI_MALLOC(u2_bz,(nr*ndat)) 355 356 my_ktabp1 = ktabp1 357 my_ktabp2 = ktabp2 358 359 if (ndat==1) then 360 do ir=1,nr 361 u1_bz(ir) = u1(ktabr1(ir))*my_ktabp1 362 end do 363 do ir=1,nr 364 u2_bz(ir) = u2(ktabr2(ir))*my_ktabp2 365 end do 366 else 367 !$OMP PARALLEL PRIVATE(padat) 368 !$OMP DO 369 do dat=1,ndat 370 padat = (dat-1)*nr 371 do ir=1,nr 372 u1_bz(ir+padat) = u1(ktabr1(ir)+padat)*my_ktabp1 373 end do 374 end do 375 !$OMP END DO NOWAIT 376 !$OMP DO 377 do dat=1,ndat 378 padat = (dat-1)*nr 379 do ir=1,nr 380 u2_bz(ir+padat) = u2(ktabr2(ir)+padat)*my_ktabp2 381 end do 382 end do 383 !$OMP END DO NOWAIT 384 !$OMP END PARALLEL 385 end if 386 387 ! Treat time-reversal. 388 SELECT CASE (time1) 389 CASE (1) 390 if (ndat==1) then 391 if (time2==1) then 392 do ir=1,nr 393 u12prod(ir) = GWPC_CONJG(u1_bz(ir)) * u2_bz(ir) 394 end do 395 else if (time2==2) then 396 do ir=1,nr 397 u12prod(ir) = GWPC_CONJG(u1_bz(ir)) * GWPC_CONJG(u2_bz(ir)) 398 end do 399 else 400 ABI_ERROR("Wrong time2") 401 end if 402 else 403 if (time2==1) then 404 !$OMP PARALLEL DO PRIVATE(padat) 405 do dat=1,ndat 406 padat = (dat-1)*nr 407 do ir=1,nr 408 u12prod(ir+padat) = GWPC_CONJG(u1_bz(ir+padat)) * u2_bz(ir+padat) 409 end do 410 end do 411 else if (time2==2) then 412 !$OMP PARALLEL DO PRIVATE(padat) 413 do dat=1,ndat 414 padat = (dat-1)*nr 415 do ir=1,nr 416 u12prod(ir+padat) = GWPC_CONJG(u1_bz(ir+padat)) * GWPC_CONJG(u2_bz(ir+padat)) 417 end do 418 end do 419 else 420 ABI_ERROR("Wrong time2") 421 end if 422 end if 423 424 CASE (2) 425 if (ndat==1) then 426 if (time2==1) then 427 do ir=1,nr 428 u12prod(ir) = u1_bz(ir) * u2_bz(ir) 429 end do 430 else if (time2==2) then 431 do ir=1,nr 432 u12prod(ir) = u1_bz(ir) * GWPC_CONJG(u2_bz(ir)) 433 end do 434 else 435 ABI_ERROR("Wrong time2") 436 end if 437 else 438 if (time2==1) then 439 !$OMP PARALLEL DO PRIVATE(padat) 440 do dat=1,ndat 441 padat = (dat-1)*nr 442 do ir=1,nr 443 u12prod(ir+padat) = u1_bz(ir+padat) * u2_bz(ir+padat) 444 end do 445 end do 446 else if (time2==2) then 447 !$OMP PARALLEL DO PRIVATE(padat) 448 do dat=1,ndat 449 padat = (dat-1)*nr 450 do ir=1,nr 451 u12prod(ir+padat) = u1_bz(ir+padat) * GWPC_CONJG(u2_bz(ir+padat)) 452 end do 453 end do 454 else 455 ABI_ERROR("Wrong time2") 456 end if 457 end if 458 CASE DEFAULT 459 ABI_ERROR("Wrong time1") 460 END SELECT 461 462 ABI_FREE(u1_bz) 463 ABI_FREE(u2_bz) 464 465 end subroutine usur_kkp_bz