TABLE OF CONTENTS


m_oscillators/calc_wfwfg [ Functions ]

[ Top ] [ 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

PARENTS

      calc_sigc_me,cchi0,cchi0q0,cohsex_me

CHILDREN

SOURCE

713 subroutine calc_wfwfg(ktabr_k,ktabi_k,spinrot,nr,nspinor,ngfft_gw,wfr_jb,wfr_kb,wfg2_jk)
714 
715 
716 !This section has been created automatically by the script Abilint (TD).
717 !Do not modify the following lines by hand.
718 #undef ABI_FUNC
719 #define ABI_FUNC 'calc_wfwfg'
720 !End of the abilint section
721 
722  implicit none
723 
724 !Arguments ------------------------------------
725 !scalars
726  integer,intent(in) :: ktabi_k,nr,nspinor
727 !arrays
728  integer,intent(in) :: ktabr_k(nr),ngfft_gw(18)
729  real(dp),intent(in) :: spinrot(4)
730  complex(gwpc),intent(in) :: wfr_jb(nr*nspinor),wfr_kb(nr*nspinor)
731  complex(gwpc),intent(out) :: wfg2_jk(nr*nspinor)
732 
733 !Local variables-------------------------------
734  integer,parameter :: ndat1=1
735  type(fftbox_plan3_t) :: plan
736 !arrays
737  complex(gwpc),allocatable :: wfr2_dpcplx(:),ujb_bz(:),ukb_bz(:)
738 
739 ! *************************************************************************
740 
741  ! There is no need to take into account phases arising from non-symmorphic
742  ! operations since the wavefunctions are evaluated at the same k-point.
743  ABI_MALLOC(wfr2_dpcplx, (nr * nspinor * ndat1))
744 
745  if (nspinor == 1) then
746    select case (ktabi_k)
747    case (1)
748      wfr2_dpcplx = GWPC_CONJG(wfr_jb(ktabr_k)) * wfr_kb(ktabr_k)
749    case (2)
750      ! Conjugate the product if time-reversal is used to reconstruct this k-point
751      wfr2_dpcplx = wfr_jb(ktabr_k) * GWPC_CONJG(wfr_kb(ktabr_k))
752    case default
753      MSG_ERROR(sjoin("Wrong ktabi_k:", itoa(ktabi_k)))
754    end select
755 
756  else if (nspinor == 2) then
757    ABI_MALLOC(ujb_bz, (nr * nspinor * ndat1))
758    ABI_MALLOC(ukb_bz, (nr * nspinor * ndat1))
759    ! Use wfr2_dpcplx as workspace array
760    call rotate_spinor(ktabi_k, ktabr_k, cone, spinrot, nr, nspinor, ndat1, wfr_jb, wfr2_dpcplx, ujb_bz)
761    call rotate_spinor(ktabi_k, ktabr_k, cone, spinrot, nr, nspinor, ndat1, wfr_kb, wfr2_dpcplx, ukb_bz)
762    wfr2_dpcplx = GWPC_CONJG(ujb_bz) * ukb_bz
763    ABI_FREE(ujb_bz)
764    ABI_FREE(ukb_bz)
765 
766  else
767    MSG_ERROR(sjoin("Wrong nspinor:", itoa(nspinor)))
768  end if
769 
770  ! Transform to Fourier space (result in wfg2_jk)
771  call fftbox_plan3_many(plan, nspinor, ngfft_gw(1:3), ngfft_gw(1:3), ngfft_gw(7), -1)
772  call fftbox_execute(plan, wfr2_dpcplx, wfg2_jk)
773  ABI_FREE(wfr2_dpcplx)
774 
775 end subroutine calc_wfwfg

m_oscillators/get_uug [ Functions ]

[ Top ] [ Functions ]

NAME

 get_uug

FUNCTION

 Calculate usug(G)= \intl \dr u1(r) exp(-i(q+G).r) u2(r) for ndat pair of wavefunctions

INPUTS

 npw=number of plane waves in the sphere.
 nr=number of FFT grid points
 ndat=Number of wavefunctions to transform.
 ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
 use_padfft=
             1 if matrix elements are calculated via zero-padded FFT.
             0 R-->G Transform in done on the full FFT box.
 igfftg0(npw)=index of G-G_o in the FFT array for each G in the sphere.
 u1(nr*ndat),u2(nr*ndat)=the two wavefunctions (periodic part)
 [trans1] = "C" if the complex conjugate of u1 should be used. Default is "N"
 [trans2] = "C" if the complex conjugate of u2 should be used. Default is "N"

OUTPUT

 usug(npw*ndat)=density of a pair of states, in reciprocal space

PARENTS

CHILDREN

SOURCE

253 subroutine get_uug(npw,nr,ndat,ngfft,use_padfft,igfftg0,gbound,u1,u2,usug,trans1,trans2)
254 
255 
256 !This section has been created automatically by the script Abilint (TD).
257 !Do not modify the following lines by hand.
258 #undef ABI_FUNC
259 #define ABI_FUNC 'get_uug'
260 !End of the abilint section
261 
262  implicit none
263 
264 !Arguments ------------------------------------
265 !scalars
266  integer,intent(in) :: npw,nr,use_padfft,ndat
267  character(len=*),optional,intent(in) :: trans1,trans2
268 !arrays
269  integer,intent(in) :: gbound(:,:) !gbound(2*mgfft+8,2)
270  integer,intent(in) :: igfftg0(npw),ngfft(18)
271  complex(gwpc),intent(in) :: u1(nr*ndat),u2(nr*ndat)
272  complex(gwpc),intent(out) :: usug(npw*ndat)
273 
274 !Local variables-------------------------------
275 !scalars
276  integer :: nx,ny,nz,ldx,ldy,ldz,mgfft
277  character(len=1) :: my_trans1,my_trans2
278  type(fftbox_plan3_t) :: plan
279 !arrays
280  complex(gwpc),allocatable :: u12prod(:)
281 
282 ! *************************************************************************
283 
284  ABI_MALLOC(u12prod,(nr*ndat))
285 
286  my_trans1 = "N"; if (PRESENT(trans1)) my_trans1 = toupper(trans1(1:1))
287  my_trans2 = "N"; if (PRESENT(trans2)) my_trans2 = toupper(trans2(1:1))
288 
289  if (my_trans1=="N" .and. my_trans2=="N") then
290    u12prod = u1 * u2
291  else if (my_trans1=="H" .and. my_trans2=="N") then
292    u12prod = GWPC_CONJG(u1) * u2
293  else if (my_trans1=="N" .and. my_trans2=="H") then
294    u12prod = u1 * GWPC_CONJG(u2)
295  else if (my_trans1=="H" .and. my_trans2=="H") then
296    u12prod = GWPC_CONJG(u1) * GWPC_CONJG(u2)
297  else
298    MSG_ERROR("Wrong combination of trans1, trans2")
299  end if
300 
301  ! Call zero-padded FFT routines if required.
302  if (use_padfft==1) then
303    nx =ngfft(1); ny =ngfft(2); nz =ngfft(3); mgfft = MAXVAL(ngfft(1:3))
304    ldx=nx; ldy=ny; ldz=nz
305    call fftpad(u12prod,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,-1,gbound)
306  else
307    call fftbox_plan3_many(plan,ndat,ngfft(1:3),ngfft(1:3),ngfft(7),-1)
308    call fftbox_execute(plan,u12prod)
309  end if
310 
311  ! From the FFT to the G-sphere.
312  call gw_box2gsph(nr,ndat,npw,igfftg0,u12prod,usug)
313 
314  ABI_FREE(u12prod)
315 
316 end subroutine get_uug

m_oscillators/gw_box2gsph [ Functions ]

[ Top ] [ 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.

PARENTS

      m_oscillators

CHILDREN

SOURCE

633 subroutine gw_box2gsph(nr,ndat,npw,igfftg0,iarrbox,oarrsph)
634 
635 
636 !This section has been created automatically by the script Abilint (TD).
637 !Do not modify the following lines by hand.
638 #undef ABI_FUNC
639 #define ABI_FUNC 'gw_box2gsph'
640 !End of the abilint section
641 
642  implicit none
643 
644 !Arguments ------------------------------------
645 !scalars
646  integer,intent(in) :: nr,ndat,npw
647 !arrays
648  integer,intent(in) :: igfftg0(npw)
649  complex(gwpc),intent(in) :: iarrbox(nr*ndat)
650  complex(gwpc),intent(out) :: oarrsph(npw*ndat)
651 
652 !Local variables-------------------------------
653 !scalars
654  integer :: ig,igfft,dat,pgsp,pfft
655 
656 ! *************************************************************************
657 
658  if (ndat==1) then
659    do ig=1,npw
660      igfft=igfftg0(ig)
661      if (igfft/=0) then
662        ! G-G0 belong to the FFT mesh.
663        oarrsph(ig) = iarrbox(igfft)
664      else
665        ! Set this component to zero.
666        oarrsph(ig) = czero_gw
667      end if
668    end do
669  else
670 !$OMP PARALLEL DO PRIVATE(pgsp,pfft,igfft)
671    do dat=1,ndat
672      pgsp = (dat-1)*npw
673      pfft = (dat-1)*nr
674      do ig=1,npw
675        igfft=igfftg0(ig)
676        if (igfft/=0) then
677          ! G-G0 belong to the FFT mesh.
678          oarrsph(ig+pgsp) = iarrbox(igfft+pfft)
679        else
680          ! Set this component to zero.
681          oarrsph(ig+pgsp) = czero_gw
682        end if
683      end do
684    end do
685  end if
686 
687 end subroutine gw_box2gsph

m_oscillators/rho_tw_g [ Functions ]

[ Top ] [ 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

PARENTS

      calc_sigc_me,calc_sigx_me,cchi0,cchi0q0,cchi0q0_intraband,cohsex_me
      exc_build_block,exc_build_ham,m_fft_prof,prep_calc_ucrpa

CHILDREN

SOURCE

102 subroutine rho_tw_g(nspinor,npwvec,nr,ndat,ngfft,map2sphere,use_padfft,igfftg0,gbound,&
103 & wfn1,i1,ktabr1,ktabp1,spinrot1,&
104 & wfn2,i2,ktabr2,ktabp2,spinrot2,&
105 & dim_rtwg,rhotwg) !& nhat12)
106 
107 
108 !This section has been created automatically by the script Abilint (TD).
109 !Do not modify the following lines by hand.
110 #undef ABI_FUNC
111 #define ABI_FUNC 'rho_tw_g'
112 !End of the abilint section
113 
114  implicit none
115 
116 !Arguments ------------------------------------
117 !scalars
118  integer,intent(in) :: i1,i2,npwvec,nr,nspinor,dim_rtwg,map2sphere,use_padfft,ndat
119  complex(dpc),intent(in) :: ktabp1,ktabp2
120 !arrays
121  integer,intent(in) :: gbound(:,:) !gbound(2*mgfft+8,2)
122  integer,intent(in) :: igfftg0(npwvec*map2sphere),ngfft(18)
123  integer,intent(in) :: ktabr1(nr),ktabr2(nr)
124  real(dp),intent(in) :: spinrot1(4),spinrot2(4)
125  complex(gwpc),intent(in) :: wfn1(nr*nspinor*ndat),wfn2(nr*nspinor*ndat)
126  complex(gwpc),intent(out) :: rhotwg(npwvec*dim_rtwg*ndat)
127 ! real(dp),optional,intent(in) :: nhat12(2,nr,nspinor**2*ndat)
128 
129 !Local variables-------------------------------
130 !scalars
131  integer :: ig,igfft,iab,spad1,spad2,spad0,nx,ny,nz,ldx,ldy,ldz,mgfft
132  type(fftbox_plan3_t) :: plan
133 !arrays
134  integer :: spinor_pad(2,4)
135  complex(gwpc),allocatable :: u12prod(:),cwavef1(:),cwavef2(:),cwork(:)
136 
137 ! *************************************************************************
138 
139  SELECT CASE (nspinor)
140  CASE (1)
141    ! Collinear case.
142    call ts_usug_kkp_bz(npwvec,nr,ndat,ngfft,map2sphere,use_padfft,igfftg0,gbound,&
143 &    wfn1,i1,ktabr1,ktabp1,&
144 &    wfn2,i2,ktabr2,ktabp2,rhotwg)
145 
146  CASE (2)
147    ! Spinorial case.
148    ABI_CHECK(ndat==1,"ndat != 1 not coded")
149    ABI_MALLOC(cwavef1, (nr * nspinor * ndat))
150    ABI_MALLOC(cwavef2, (nr * nspinor * ndat))
151    ABI_MALLOC(cwork, (nr * nspinor * ndat))
152 
153    call rotate_spinor(i1, ktabr1, ktabp1, spinrot1, nr, nspinor, ndat, wfn1, cwork, cwavef1)
154    call rotate_spinor(i2, ktabr2, ktabp2, spinrot2, nr, nspinor, ndat, wfn2, cwork, cwavef2)
155 
156    ABI_FREE(cwork)
157    ABI_MALLOC(u12prod, (nr))
158 
159    spinor_pad = reshape([0, 0, nr, nr, 0, nr, nr, 0], [2, 4])
160    rhotwg = czero_gw
161    do iab=1,2
162      spad1 = spinor_pad(1,iab); spad2=spinor_pad(2,iab)
163 
164      u12prod = GWPC_CONJG(cwavef1(spad1+1:spad1+nr)) * cwavef2(spad2+1:spad2+nr)
165      ! Add compensation charge.
166      !if (PRESENT(nhat12)) u12prod = u12prod + CMPLX(nhat12(1,:,iab),nhat12(2,:,iab))
167 
168      spad0 = (iab-1)*npwvec
169      SELECT CASE (map2sphere)
170      CASE (0)
171        ! Need results on the full FFT box thus cannot use zero-padded FFT.
172        call fftbox_plan3_many(plan,ndat,ngfft(1:3),ngfft(1:3),ngfft(7),-1)
173        call fftbox_execute(plan,u12prod)
174        if (dim_rtwg == 1) then
175          rhotwg(1:npwvec) = rhotwg(1:npwvec) + u12prod
176        else
177          rhotwg(spad0+1:spad0+npwvec) = u12prod
178        end if
179 
180      CASE (1)
181        ! Need results on the G-sphere. Call zero-padded FFT routines if required.
182        if (use_padfft == 1) then
183          nx =ngfft(1); ny =ngfft(2); nz =ngfft(3); mgfft = MAXVAL(ngfft(1:3))
184          ldx=nx; ldy=ny; ldz=nz
185          call fftpad(u12prod,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,-1,gbound)
186        else
187          call fftbox_plan3_many(plan,ndat,ngfft(1:3),ngfft(1:3),ngfft(7),-1)
188          call fftbox_execute(plan,u12prod)
189        end if
190 
191        ! Have to map FFT to G-sphere.
192        if (dim_rtwg == 1) then
193          do ig=1,npwvec
194            igfft = igfftg0(ig)
195            ! G-G0 belong to the FFT mesh.
196            if (igfft /= 0) rhotwg(ig) = rhotwg(ig) + u12prod(igfft)
197          end do
198        else
199          do ig=1,npwvec
200            igfft = igfftg0(ig)
201            ! G-G0 belong to the FFT mesh.
202            if (igfft /= 0) rhotwg(ig+spad0) = u12prod(igfft)
203          end do
204        end if
205 
206      CASE DEFAULT
207        MSG_BUG("Wrong map2sphere")
208      END SELECT
209    end do !iab
210 
211    ABI_FREE(u12prod)
212    ABI_FREE(cwavef1)
213    ABI_FREE(cwavef2)
214 
215  CASE DEFAULT
216    MSG_BUG('Wrong nspinor')
217  END SELECT
218 
219 end subroutine rho_tw_g

m_oscillators/rotate_spinor [ Functions ]

[ Top ] [ Functions ]

NAME

  rotate_spinor

FUNCTION

  Return a spinor in the full BZ from its symmetrical image in the IBZ.

INPUTS

OUTPUT

PARENTS

      m_oscillators

CHILDREN

SOURCE

898 subroutine rotate_spinor(itim_kbz, ktabr_kbz, ktabp_kbz, spinrot, nr, nspinor, ndat, ug_ibz, cwork, oug_bz)
899 
900 
901 !This section has been created automatically by the script Abilint (TD).
902 !Do not modify the following lines by hand.
903 #undef ABI_FUNC
904 #define ABI_FUNC 'rotate_spinor'
905 !End of the abilint section
906 
907  implicit none
908 
909 !Arguments ------------------------------------
910 !scalars
911  integer,intent(in) :: itim_kbz, nr, nspinor, ndat
912  complex(dpc),intent(in) :: ktabp_kbz
913 !arrays
914  integer,intent(in) :: ktabr_kbz(nr)
915  real(dp),intent(in) :: spinrot(4)
916  complex(gwpc),intent(in) :: ug_ibz(nr*nspinor*ndat)
917  complex(gwpc),intent(out) :: cwork(nr*nspinor*ndat), oug_bz(nr*nspinor*ndat)
918 
919 !Local variables ------------------------------
920 !scalars
921  integer :: ir,ir1,spad0,ispinor
922  complex(gwpc) :: u1a,u1b
923 !arrays
924  complex(dpc) :: spinrot_cmat1(2,2)
925 
926 !************************************************************************
927 
928  ABI_CHECK(ndat == 1, "ndat > 1 not coded")
929  ABI_CHECK(nspinor == 2, "nspinor should be 1")
930 
931  ! Apply Time-reversal if required.
932  ! \psi_{-k}^1 =  (\psi_k^2)^*
933  ! \psi_{-k}^2 = -(\psi_k^1)^*
934  if (itim_kbz == 1) then
935    cwork(:) = ug_ibz(:)
936  else if (itim_kbz == 2) then
937    cwork(1:nr) = GWPC_CONJG(ug_ibz(nr+1:2*nr))
938    cwork(nr+1:2*nr) = -GWPC_CONJG(ug_ibz(1:nr))
939  else
940    MSG_ERROR(sjoin('Wrong itim_kbz:', itoa(itim_kbz)))
941  end if
942 
943  do ispinor=1,nspinor
944    spad0 = (ispinor-1) * nr
945    do ir=1,nr
946      ir1 = ktabr_kbz(ir)
947      oug_bz(ir+spad0) = cwork(ir1+spad0) * ktabp_kbz
948    end do
949  end do
950 
951  ! Rotation in spinor space (same equations as in wfconv)
952  spinrot_cmat1 = spinrot_cmat(spinrot)
953  cwork = oug_bz
954  do ir=1,nr
955    u1a = cwork(ir); u1b = cwork(ir+nr)
956    oug_bz(ir)    = spinrot_cmat1(1, 1) * u1a + spinrot_cmat1(1, 2) * u1b
957    oug_bz(ir+nr) = spinrot_cmat1(2, 1) * u1a + spinrot_cmat1(2, 2) * u1b
958  end do
959 
960 end subroutine rotate_spinor

m_oscillators/sym_rhotwgq0 [ Functions ]

[ Top ] [ 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)

PARENTS

CHILDREN

SOURCE

814 function sym_rhotwgq0(itim_k,isym_k,dim_rtwg,npw,rhxtwg_in,Gsph) result(rhxtwg_sym)
815 
816 
817 !This section has been created automatically by the script Abilint (TD).
818 !Do not modify the following lines by hand.
819 #undef ABI_FUNC
820 #define ABI_FUNC 'sym_rhotwgq0'
821 !End of the abilint section
822 
823  implicit none
824 
825 !Arguments ------------------------------------
826 !scalars
827  integer,intent(in) :: npw,dim_rtwg,itim_k,isym_k
828  type(gsphere_t),intent(in) :: Gsph
829 !arrays
830  complex(gwpc),intent(in) :: rhxtwg_in(dim_rtwg*npw)
831  complex(gwpc) :: rhxtwg_sym(dim_rtwg*npw)
832 
833 !Local variables ------------------------------
834 !scalars
835  integer :: ig
836 
837 !************************************************************************
838 
839  ABI_CHECK(dim_rtwg == 1, "dim_rtwg/=1 not coded")
840 
841  SELECT CASE (isym_k)
842  CASE (1)
843    ! Fractional translation associated to E is assumed to be (zero,zero,zero).
844    SELECT CASE (itim_k)
845    CASE (1)
846      ! Identity, no time-reversal. No symmetrization is needed.
847      rhxtwg_sym(:) = rhxtwg_in(:)
848    CASE (2)
849      ! Identity + Time-reversal.
850      do ig=1,npw
851        rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = GWPC_CONJG(rhxtwg_in(ig))
852      end do
853    CASE DEFAULT
854      MSG_ERROR(sjoin("Wrong value of itim_k:", itoa(itim_k)))
855    END SELECT
856 
857  CASE DEFAULT
858    ! Rotate wavefunctions.
859    SELECT CASE (itim_k)
860    CASE (1)
861      ! no time-reversal, only rotation.
862      do ig=1,npw
863        rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = rhxtwg_in(ig) * Gsph%phmSGt(ig,isym_k)
864      end do
865    CASE (2)
866      ! time-reversal + spatial rotation.
867      do ig=1,npw
868        rhxtwg_sym( Gsph%rottb(ig,itim_k,isym_k) ) = GWPC_CONJG( rhxtwg_in(ig) * Gsph%phmSGt(ig,isym_k) )
869      end do
870    CASE DEFAULT
871      MSG_ERROR(sjoin("Wrong value of itim_k:", itoa(itim_k)))
872    END SELECT
873  END SELECT
874 
875 end function sym_rhotwgq0

m_oscillators/ts_usug_kkp_bz [ Functions ]

[ Top ] [ 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

PARENTS

      m_oscillators

CHILDREN

SOURCE

357 subroutine ts_usug_kkp_bz(npw,nr,ndat,ngfft,map2sphere,use_padfft,igfftg0,gbound,&
358 & u1,time1,ktabr1,ktabp1,&
359 & u2,time2,ktabr2,ktabp2,usug) !& nhat12)
360 
361 
362 !This section has been created automatically by the script Abilint (TD).
363 !Do not modify the following lines by hand.
364 #undef ABI_FUNC
365 #define ABI_FUNC 'ts_usug_kkp_bz'
366 !End of the abilint section
367 
368  implicit none
369 
370 !Arguments ------------------------------------
371 !scalars
372  integer,intent(in) :: time1,time2,npw,nr,map2sphere,use_padfft,ndat
373  complex(dpc),intent(in) :: ktabp1,ktabp2
374 !arrays
375  integer,intent(in) :: gbound(:,:) !gbound(2*mgfft+8,2)
376  integer,intent(in) :: igfftg0(npw*map2sphere),ngfft(18)
377  integer,intent(in) :: ktabr1(nr),ktabr2(nr)
378  complex(gwpc),intent(in) :: u1(nr*ndat),u2(nr*ndat)
379  complex(gwpc),intent(out) :: usug(npw*ndat)
380 
381 !Local variables-------------------------------
382 !scalars
383  integer :: nx,ny,nz,ldx,ldy,ldz,mgfft
384  type(fftbox_plan3_t) :: plan
385 !arrays
386  complex(gwpc),allocatable :: u12prod(:)
387 
388 ! *************************************************************************
389 
390  ! Form rho-twiddle(r)=u_1^*(r,b1,kbz1) u_2(r,b2,kbz2), to account for symmetries:
391  ! u(r,b,kbz)=e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t),b,kibz)
392  !           =e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t),b,kibz) for time-reversal
393  !
394  ABI_MALLOC(u12prod,(nr*ndat))
395  call usur_kkp_bz(nr,ndat,time1,ktabr1,ktabp1,u1,time2,ktabr2,ktabp2,u2,u12prod)
396 
397  ! Add compensation charge.
398  !if (PRESENT(nhat12)) u12prod = u1prod + CMPLX(nhat12(1,:,1),nhat12(2,:,1))
399 
400  SELECT CASE (map2sphere)
401  CASE (0)
402    ! Need results on the full FFT box thus cannot use zero-padded FFT.
403    call fftbox_plan3_many(plan,ndat,ngfft(1:3),ngfft(1:3),ngfft(7),-1)
404    call fftbox_execute(plan,u12prod)
405    call xcopy(nr*ndat,u12prod,1,usug,1)
406 
407  CASE (1)
408    ! Need results on the G-sphere. Call zero-padded FFT routines if required.
409    if (use_padfft==1) then
410      nx = ngfft(1); ny = ngfft(2); nz = ngfft(3); mgfft = MAXVAL(ngfft(1:3))
411      ldx=nx; ldy=ny; ldz=nz
412      call fftpad(u12prod,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat,mgfft,-1,gbound)
413    else
414      call fftbox_plan3_many(plan,ndat,ngfft(1:3),ngfft(1:3),ngfft(7),-1)
415      call fftbox_execute(plan,u12prod)
416    end if
417 
418    ! From the FFT to the G-sphere.
419    call gw_box2gsph(nr,ndat,npw,igfftg0,u12prod,usug)
420 
421  CASE DEFAULT
422    MSG_BUG("Wrong map2sphere")
423  END SELECT
424 
425  ABI_FREE(u12prod)
426 
427 end subroutine ts_usug_kkp_bz

m_oscillators/usur_kkp_bz [ Functions ]

[ Top ] [ 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 iin 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.

PARENTS

      m_oscillators

CHILDREN

SOURCE

459 subroutine usur_kkp_bz(nr,ndat,time1,ktabr1,ktabp1,u1,time2,ktabr2,ktabp2,u2,u12prod)
460 
461 
462 !This section has been created automatically by the script Abilint (TD).
463 !Do not modify the following lines by hand.
464 #undef ABI_FUNC
465 #define ABI_FUNC 'usur_kkp_bz'
466 !End of the abilint section
467 
468  implicit none
469 
470 !Arguments ------------------------------------
471 !scalars
472  integer,intent(in) :: nr,ndat,time1,time2
473  complex(dpc),intent(in) :: ktabp1,ktabp2
474 !arrays
475  integer,intent(in) :: ktabr1(nr),ktabr2(nr)
476  complex(gwpc),intent(in) :: u1(nr*ndat),u2(nr*ndat)
477  complex(gwpc),intent(out) :: u12prod(nr*ndat)
478 
479 !Local variables-------------------------------
480 !scalars
481  integer :: ir,dat,padat
482  complex(gwpc) :: my_ktabp1,my_ktabp2
483 !arrays
484  complex(gwpc),allocatable :: u1_bz(:),u2_bz(:)
485 
486 ! *************************************************************************
487 
488  ! Form rho-twiddle(r)=u_1^*(r,b1,kbz1) u_2(r,b2,kbz2), to account for symmetries:
489  ! u(r,b,kbz)=e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t),b,kibz)
490  !           =e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t),b,kibz) for time-reversal
491  !
492  ABI_MALLOC(u1_bz,(nr*ndat))
493  ABI_MALLOC(u2_bz,(nr*ndat))
494 
495  my_ktabp1 = ktabp1
496  my_ktabp2 = ktabp2
497 
498  if (ndat==1) then
499    do ir=1,nr
500      u1_bz(ir) = u1(ktabr1(ir))*my_ktabp1
501    end do
502    do ir=1,nr
503      u2_bz(ir) = u2(ktabr2(ir))*my_ktabp2
504    end do
505  else
506 !$OMP PARALLEL PRIVATE(padat)
507 !$OMP DO
508    do dat=1,ndat
509      padat = (dat-1)*nr
510      do ir=1,nr
511        u1_bz(ir+padat) = u1(ktabr1(ir)+padat)*my_ktabp1
512      end do
513    end do
514 !$OMP END DO NOWAIT
515 !$OMP DO
516    do dat=1,ndat
517      padat = (dat-1)*nr
518      do ir=1,nr
519        u2_bz(ir+padat) = u2(ktabr2(ir)+padat)*my_ktabp2
520      end do
521    end do
522 !$OMP END DO NOWAIT
523 !$OMP END PARALLEL
524  end if
525 
526  ! Treat time-reversal.
527  SELECT CASE (time1)
528  CASE (1)
529    if (ndat==1) then
530      if (time2==1) then
531        do ir=1,nr
532          u12prod(ir) = GWPC_CONJG(u1_bz(ir)) * u2_bz(ir)
533        end do
534      else if (time2==2) then
535        do ir=1,nr
536          u12prod(ir) = GWPC_CONJG(u1_bz(ir)) * GWPC_CONJG(u2_bz(ir))
537        end do
538      else
539        MSG_ERROR("Wrong time2")
540      end if
541    else
542      if (time2==1) then
543 !$OMP PARALLEL DO PRIVATE(padat)
544        do dat=1,ndat
545          padat = (dat-1)*nr
546          do ir=1,nr
547            u12prod(ir+padat) = GWPC_CONJG(u1_bz(ir+padat)) * u2_bz(ir+padat)
548          end do
549        end do
550      else if (time2==2) then
551 !$OMP PARALLEL DO PRIVATE(padat)
552        do dat=1,ndat
553          padat = (dat-1)*nr
554          do ir=1,nr
555            u12prod(ir+padat) = GWPC_CONJG(u1_bz(ir+padat)) * GWPC_CONJG(u2_bz(ir+padat))
556          end do
557        end do
558      else
559        MSG_ERROR("Wrong time2")
560      end if
561    end if
562 
563  CASE (2)
564    if (ndat==1) then
565      if (time2==1) then
566        do ir=1,nr
567          u12prod(ir) = u1_bz(ir) * u2_bz(ir)
568        end do
569      else if (time2==2) then
570        do ir=1,nr
571          u12prod(ir) = u1_bz(ir) * GWPC_CONJG(u2_bz(ir))
572        end do
573      else
574        MSG_ERROR("Wrong time2")
575      end if
576    else
577      if (time2==1) then
578 !$OMP PARALLEL DO PRIVATE(padat)
579        do dat=1,ndat
580          padat = (dat-1)*nr
581          do ir=1,nr
582            u12prod(ir+padat) = u1_bz(ir+padat) * u2_bz(ir+padat)
583          end do
584        end do
585      else if (time2==2) then
586 !$OMP PARALLEL DO PRIVATE(padat)
587        do dat=1,ndat
588          padat = (dat-1)*nr
589          do ir=1,nr
590            u12prod(ir+padat) = u1_bz(ir+padat) * GWPC_CONJG(u2_bz(ir+padat))
591          end do
592        end do
593      else
594        MSG_ERROR("Wrong time2")
595      end if
596    end if
597  CASE DEFAULT
598    MSG_ERROR("Wrong time1")
599  END SELECT
600 
601  ABI_FREE(u1_bz)
602  ABI_FREE(u2_bz)
603 
604 end subroutine usur_kkp_bz