TABLE OF CONTENTS


ABINIT/m_exc_analyze [ Modules ]

[ Top ] [ Modules ]

NAME

  m_exc_analyze

FUNCTION

  Postprocessing tools for BSE calculations.

COPYRIGHT

  Copyright (C) 1992-2009 EXC group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida)
  Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_exc_analyze
29 
30  use defs_basis
31  use defs_datatypes
32  use m_abicore
33  use m_bs_defs
34  use m_xmpi
35  use m_errors
36 
37  use m_io_tools,          only : open_file
38  use m_numeric_tools,     only : iseven, wrap2_zero_one
39  use m_bz_mesh,           only : kmesh_t, get_BZ_item
40  use m_crystal,           only : crystal_t
41  use m_wfd,               only : wfd_t, wfd_change_ngfft, wfd_get_cprj, wfd_sym_ur, wfd_get_ur
42  use m_bse_io,            only : exc_read_eigen
43  use m_pptools,           only : printxsf
44 
45  use m_pawrad,            only : pawrad_type
46  use m_pawtab,            only : pawtab_type,pawtab_get_lsize
47  use m_pawfgrtab,         only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
48  use m_pawcprj,           only : pawcprj_type, pawcprj_alloc, pawcprj_free
49  use m_paw_pwaves_lmn,    only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
50  use m_paw_nhat,          only : nhatgrid
51 
52  implicit none
53 
54  private

m_exc_analyze/exc_den [ Functions ]

[ Top ] [ m_exc_analyze ] [ Functions ]

NAME

  exc_den

FUNCTION

  This routines calculates the electron-hole excited state density.

INPUTS

  filbseig=Name of the file containing the excitonic eigenvectors and eigenvalues.

OUTPUT

PARENTS

      bethe_salpeter

CHILDREN

      wfd_get_ur,wrtout

SOURCE

453 subroutine exc_den(BSp,BS_files,ngfft,nfftot,Kmesh,ktabr,Wfd)
454 
455 
456 !This section has been created automatically by the script Abilint (TD).
457 !Do not modify the following lines by hand.
458 #undef ABI_FUNC
459 #define ABI_FUNC 'exc_den'
460 !End of the abilint section
461 
462  implicit none
463 
464 !Arguments ------------------------------------
465 !scalars
466  integer,intent(in) :: nfftot
467  type(excparam),intent(in) :: BSp
468  type(excfiles),intent(in) :: BS_files
469  type(kmesh_t),intent(in) :: Kmesh
470  type(wfd_t),intent(inout) :: Wfd
471 !arrays
472  integer,intent(in) :: ngfft(18)
473  integer,intent(in) :: ktabr(nfftot,BSp%nkbz)
474 
475 !Local variables ------------------------------
476 !scalars
477  integer :: nfft1,nfft2,nfft3,spin,ierr
478  integer :: it,itp,ik_ibz,ikp_bz,ik_bz,band,iv,ivp,ic,icp,ir,i1,i2,i3,iik
479  integer :: state,nstates,min_idx,eig_unt,den_unt,sden_unt,hsize,homo
480  real(dp) :: min_ene,cost
481  character(len=500) :: msg
482  character(len=fnlen) :: filbseig
483 !arrays
484  real(dp) :: n0(nfftot),rho_eh(nfftot),nexc(nfftot)
485  real(dp),allocatable :: exc_ene(:)
486  complex(dpc) :: rhor_h(nfftot),rhor_e(nfftot)
487  complex(gwpc),allocatable :: wfr(:,:,:), wfrk(:,:)
488  complex(dpc),allocatable :: exc_ene_cplx(:),exc_vec(:)
489 
490 !************************************************************************
491 
492  ABI_CHECK(Wfd%nsppol==1,"nsppol==2 not coded")
493 
494  if (BS_files%in_eig /= BSE_NOFILE) then
495    filbseig = BS_files%in_eig
496  else
497    filbseig = BS_files%out_eig
498  end if
499 
500  call wrtout(std_out,' Calculating electron-hole, excited state density',"COLL")
501  call wrtout(std_out," Reading eigenstates from: "//TRIM(filbseig),"COLL")
502  MSG_ERROR("Not tested")
503 
504  ! Prepare FFT tables to have u(r) on the ngfft_osc mesh.
505  !if ( ANY(ngfftf(1:3) /= Wfd%ngfft(1:3)) ) then
506  !  call wfd_change_ngfft(Wfd,Cryst,Psps,ngfftf)
507  !end if
508 
509  nfft1 = ngfft(1)
510  nfft2 = ngfft(2)
511  nfft3 = ngfft(3)
512  ABI_CHECK(nfftot==PRODUCT(ngfft(1:3)),"Mismatch in FFT size")
513 
514 !allocate and load wavefunctions in real space
515  ABI_STAT_MALLOC(wfr,(nfftot*Wfd%nspinor,BSp%nbnds,Wfd%nkibz), ierr)
516  ABI_CHECK(ierr==0, "out of memory: exc_den, wfr")
517 
518  spin=1 ! SPIN support is missing.
519 
520  do ik_ibz=1,Wfd%nkibz
521    do band=1,BSp%nbnds
522      call wfd_get_ur(Wfd,band,ik_ibz,spin,wfr(:,band,ik_ibz))
523    end do
524  end do
525 
526  if (open_file(filbseig,msg,newunit=eig_unt,form="unformatted", status="old", action="read") /= 0) then
527    MSG_ERROR(msg)
528  end if
529 
530  read(eig_unt) hsize,nstates
531 
532  if (nstates/=Bsp%nreh(1)) then
533    write(std_out,*) 'not resonant calculation'
534    close(eig_unt)
535    RETURN
536  end if
537 
538  ABI_MALLOC(exc_ene_cplx,(nstates))
539  ABI_MALLOC(exc_ene,(nstates))
540 
541  read(eig_unt) exc_ene_cplx(1:nstates)
542 
543  ! Small imaginary part is always neglected.
544  exc_ene(:) = exc_ene_cplx(:)
545  ABI_FREE(exc_ene_cplx)
546  !
547  ! Find the lowest non-negative eigenvalue.
548  min_ene = greatest_real
549  do state=1,nstates
550    if (exc_ene(state) < zero) cycle
551    if (exc_ene(state) < min_ene) then
552      min_ene = exc_ene(state)
553      min_idx = state
554    end if
555  end do
556  ABI_FREE(exc_ene)
557 
558  write(std_out,*)' Lowest eigenvalue ', min_idx, min_ene*Ha_eV
559  !
560  ! skip other eigenvectors
561  do state=1,min_idx-1
562    read(eig_unt)
563  end do
564  !
565  ! read "lowest" eigenvector
566  ABI_MALLOC(exc_vec,(hsize))
567  read(eig_unt) exc_vec
568 
569  close(eig_unt)
570 
571  ABI_STAT_MALLOC(wfrk,(nfftot,BSp%nbnds), ierr)
572  ABI_CHECK(ierr==0, 'out of memory: exc_den, wfrk')
573 
574 !calculate ground state density
575  n0(:) = zero
576  spin = 1
577  homo = Bsp%homo_spin(spin)
578  do ik_bz = 1, BSp%nkbz
579    ik_ibz = Kmesh%tab(ik_bz)
580    iik = (3-Kmesh%tabi(ik_bz))/2
581 
582    if (iik==1) then
583      do ir=1,nfftot
584        wfrk(ir,1:homo) = wfr(ktabr(ir,ik_bz),1:homo,ik_ibz)
585      end do
586    else
587      do ir=1,nfftot
588        wfrk(ir,1:homo) = conjg(wfr(ktabr(ir,ik_bz),1:homo,ik_ibz))
589      end do
590    end if
591 
592    do iv=1,homo
593      n0(:) = n0(:) + conjg(wfrk(:,band)) * wfrk(:,band)
594    end do
595  end do
596  !
597  ! calculate electron and hole density
598  rhor_h=czero
599  rhor_e=czero
600  ABI_CHECK(BSp%nsppol==1,"nsppol=2 not coded")
601 
602  do it=1,Bsp%nreh(1)
603    ik_bz = Bsp%Trans(it,1)%k
604    iv    = Bsp%Trans(it,1)%v
605    ic    = Bsp%Trans(it,1)%c
606    ik_ibz = Kmesh%tab(ik_bz)
607    iik = (3-Kmesh%tabi(ik_bz))/2
608 
609    if (iik==1) then
610      do ir = 1, nfftot
611        wfrk(ir,:) = wfr(ktabr(ir,ik_bz),:,ik_ibz)
612      end do
613    else
614      do ir = 1, nfftot
615        wfrk(ir,:) = conjg(wfr(ktabr(ir,ik_bz),:,ik_ibz))
616      end do
617    end if
618    !
619    do itp=1,Bsp%nreh(1)
620      ikp_bz = Bsp%Trans(itp,1)%k
621      if (ikp_bz/=ik_bz) CYCLE
622      icp = Bsp%Trans(itp,1)%c
623      ivp = Bsp%Trans(itp,1)%v
624      if(icp==ic) then
625        rhor_h = rhor_h - conjg(exc_vec(it)) * exc_vec(itp) * wfrk(:,iv) * conjg(wfrk(:,ivp))
626      end if
627      if(ivp==iv) then
628        rhor_e = rhor_e + conjg(exc_vec(it)) * exc_vec(itp) * conjg(wfrk(:,ic)) * wfrk(:,icp)
629      end if
630    end do
631  end do
632 
633  ABI_FREE(exc_vec)
634  !
635  ! calculate excited state density minus ground state density
636  ! n* - n0 = rhor_e + rhor_h
637  rho_eh(:) = rhor_e + rhor_h
638  !
639  !calculate excited state density
640  ! n* = n0 + rhor_e + rhor_h
641  nexc(:) = n0(:) + rho_eh(:)
642 
643  if (open_file('out.den',msg,newunit=den_unt) /= 0) then
644    MSG_ERROR(msg)
645  end if
646 
647  ! here conversion to cartesian through a1, a2, a3
648  cost = zero
649  do i1=0,nfft1-1
650    do i2=0,nfft2-1
651      do i3=0,nfft3-1
652        ir=1 + i1 + i2*nfft1 + i3*nfft1*nfft2
653        write(den_unt,'(3i3,2x,5e11.4)')i1,i2,i3,n0(ir),nexc(ir),rho_eh(ir),real(rhor_e(ir)),real(rhor_h(ir))
654        cost = cost + n0(ir)
655      end do
656    end do
657  end do
658 
659  write(std_out,*) 'density normalization constant ', cost
660  close(den_unt)
661 
662  if (open_file('out.sden',msg,newunit=sden_unt) /= 0) then
663    MSG_ERROR(msg)
664  end if
665 
666  ! we are looking for the plane between (100) (111)
667  ! so it's the place where (111) [ (100) x v ] = 0 (mixed product)
668  ! finally v2 - v3 = 0
669  do i2=0, nfft2-1
670    do i3=0, nfft3-1
671      if(i2==i3) then
672        write(std_out,*) i2, i3
673        write(sden_unt,*) (rho_eh(1+i1+i2*nfft1+i3*nfft1*nfft2),i1=0,nfft1-1)
674      end if
675    end do
676  end do
677 
678  close(sden_unt)
679 
680 end subroutine exc_den

m_exc_analyze/m_exc_analyze [ Functions ]

[ Top ] [ m_exc_analyze ] [ Functions ]

NAME

  exc_plot

FUNCTION

  Plots the excitonic wavefunction in real space.

INPUTS

 Bsp<excparam>=Container storing the input variables of the run.
 BS_files<excfiles>=filenames used in the Bethe-Salpeter part.
 Kmesh<kmesh_t>=Info on the k-point sampling for wave functions.
 Cryst<crystal_t>=Structure defining the crystalline structure.
 KS_Bst<ebands_t>
 Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data
 Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data.
 Psps <pseudopotential_type>=variables related to pseudopotentials.
 Wfd<wfd_t>=Handler for the wavefunctions.
 ngfftf(18)=Info on the dense FFT mesh used for plotting the excitonic wavefunctions.
 nrcell(3)=Number of cell replicas (the code will enforce odd number so that the e or the
  h can be centered in the bix box.
 eh_rcoord(3)=Reduced coordinated of the (electron|hole)
 which_fixed= 1 to fix the electron at eh_red, 2 for the hole.
 spin_opt=1 for resolving the up-component, 2 for the down component, 3 if spin resolution is not wanted.
 paw_add_onsite=.TRUE. if the onsite contribution is taken into account.

OUTPUT

PARENTS

CHILDREN

      exc_read_eigen,nhatgrid,paw_pwaves_lmn_free,paw_pwaves_lmn_init
      pawcprj_alloc,pawcprj_free,pawfgrtab_free,pawfgrtab_init
      pawfgrtab_print,pawtab_get_lsize,printxsf,wfd_change_ngfft,wfd_sym_ur
      wrap2_zero_one,wrtout,xmpi_bcast,xmpi_sum_master

SOURCE

100 subroutine exc_plot(Bsp,Bs_files,Wfd,Kmesh,Cryst,Psps,Pawtab,Pawrad,paw_add_onsite,spin_opt,which_fixed,eh_rcoord,nrcell,ngfftf)
101 
102 
103 !This section has been created automatically by the script Abilint (TD).
104 !Do not modify the following lines by hand.
105 #undef ABI_FUNC
106 #define ABI_FUNC 'exc_plot'
107 !End of the abilint section
108 
109  implicit none
110 
111 !Arguments ------------------------------------
112 !scalars
113  integer,intent(in) :: which_fixed,spin_opt
114  logical,intent(in) :: paw_add_onsite
115  type(excparam),intent(in) :: Bsp
116  type(excfiles),intent(in) :: BS_files
117  type(kmesh_t),intent(in) :: Kmesh
118  type(crystal_t),intent(in) :: Cryst
119  type(pseudopotential_type),intent(in) :: Psps
120  type(wfd_t),intent(inout) :: Wfd
121 !arrays
122  integer,intent(in) :: ngfftf(18),nrcell(3)
123  real(dp),intent(in) :: eh_rcoord(3)
124  type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
125  type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*Wfd%usepaw)
126 
127 !Local variables ------------------------------
128 !scalars
129  integer,parameter :: cplex1=1
130  integer :: nsppol,usepaw,nspinor,comm,cond,val
131  integer :: optcut,optgr0,optgr1,optgr2,optrad
132  integer :: ik_bz,ierr,my_rank !ik_ibz, istwf_k, isym_k,itim_k, my_nbbp, npw_k,
133  integer :: spin,spin_start,spin_stop,reh
134  integer :: rt_idx,art_idx,ii,iatom,nr_tot
135  integer :: irc,ir1,ir2,ir3,wp1,wp2,wp3,wp_idx,eh_fft_idx,eh_rr,rr
136  integer :: hsize,xsf_unt,ncells,nvec
137  integer :: sc_natom,master
138  real(dp) :: ene_rt,k_dot_r12
139  complex(dpc) :: res_coeff,antres_coeff,eikr12
140  character(len=500) :: msg
141  character(len=fnlen) :: eig_fname,out_fname
142 !arrays
143  integer :: nrcl(3),bbox(3)
144  !integer :: bbp_distrb(Wfd%mband,Wfd%mband),got(Wfd%nproc)
145  integer,allocatable :: rcl2fft(:),sc_typat(:),l_size_atm(:),vec_idx(:)
146  real(dp),parameter :: origin0(3)=(/zero,zero,zero/)
147  real(dp) :: k_bz(3),eh_red(3),eh_shift(3),r12(3),sc_rprimd(3,3),scart_shift(3)
148  real(dp),allocatable :: rclred(:,:),exc_phi2(:),sc_xcart(:,:) !,sc_znucl(:)
149  complex(gwpc),allocatable :: exc_phi(:)
150  complex(gwpc),allocatable :: ur_v(:),ur_c(:)
151  complex(dpc),allocatable :: vec_list(:,:)
152  !logical :: bbp_mask(Wfd%mband,Wfd%mband)
153  type(pawcprj_type),allocatable :: Cp_v(:,:),Cp_c(:,:)
154  type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
155  type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
156 
157 !************************************************************************
158 
159  MSG_WARNING("Exc plot is still under development")
160 
161  call wrtout(std_out," Calculating excitonic wavefunctions in real space","COLL")
162  ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded")
163  ABI_CHECK(Wfd%usepaw==0,"PAW not coded")
164  ABI_CHECK(Wfd%nproc==1,"nproc>1 not supported")
165 
166  ! If needed, prepare FFT tables to have u(r) on the ngfftf mesh.
167  if ( ANY(ngfftf(1:3) /= Wfd%ngfft(1:3)) ) then
168    call wfd_change_ngfft(Wfd,Cryst,Psps,ngfftf)
169  end if
170 
171  comm    = Wfd%comm
172  my_rank = Wfd%my_rank
173  master  = Wfd%master
174 
175  nsppol  = Wfd%nsppol
176  nspinor = Wfd%nspinor
177  usepaw  = Wfd%usepaw
178  !
179  ! TODO recheck this part.
180  ! Prepare tables needed for splining the onsite term on the FFT box.
181  if (Wfd%usepaw==1.and.paw_add_onsite) then
182    MSG_WARNING("Testing the calculation of AE PAW wavefunctions.")
183    ! Use a local pawfgrtab to make sure we use the correction in the paw spheres
184    ! the usual pawfgrtab uses r_shape which may not be the same as r_paw.
185    call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
186    ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom))
187    call pawfgrtab_init(Pawfgrtab,cplex1,l_size_atm,Wfd%nspden,Cryst%typat)
188    ABI_FREE(l_size_atm)
189 
190    optcut=1                     ! use rpaw to construct Pawfgrtab.
191    optgr0=0; optgr1=0; optgr2=0 ! dont need gY terms locally.
192    optrad=1                     ! do store r-R.
193 
194    call nhatgrid(Cryst%atindx1,Cryst%gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,&
195 &    optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
196    !Pawfgrtab is ready to use
197 
198    if (Wfd%pawprtvol>0) then
199      call pawfgrtab_print(Pawfgrtab,natom=Cryst%natom,unit=std_out,&
200 &                         prtvol=Wfd%pawprtvol,mode_paral="COLL")
201    end if
202 
203    ABI_DT_MALLOC(Paw_onsite,(Cryst%natom))
204    call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,Cryst%rprimd,Cryst%xcart,&
205 &                           Pawtab,Pawrad,Pawfgrtab)
206  end if
207  !
208  ! Number of cells in the big box. Odd numbers are needed to place the e or the h in the middle.
209  do ii=1,3
210    nrcl(ii) = nrcell(ii)
211    if (iseven(nrcell(ii))) then
212      nrcl(ii) = nrcell(ii)+1
213      write(msg,'(2(a,i0))')" Enforcing odd number of cell replicas ",nrcell(ii)," --> ",nrcl(ii)
214      MSG_WARNING(msg)
215    end if
216  end do
217 
218  ncells = PRODUCT(nrcl)
219  nr_tot = Wfd%nfftot * ncells  ! Total number of points in the big box.
220  bbox = nrcl*ngfftf(1:3)
221  !
222  ! rcl2fft: The image of the point in the small box.
223  ! rcl2red: The reduced coordinates of the point in the big box in terms of rprimd.
224  ABI_MALLOC(rcl2fft,(nr_tot))
225  ABI_MALLOC(rclred,(3,nr_tot))
226 
227  irc = 0
228  do ir3=0,bbox(3)-1 ! Loop over the points in the big box.
229    do ir2=0,bbox(2)-1
230      do ir1=0,bbox(1)-1
231        irc = 1+irc
232        !
233        wp1=MODULO(ir1,ngfftf(1)) ! The FFT index of the point wrapped into the original cell.
234        wp2=MODULO(ir2,ngfftf(2))
235        wp3=MODULO(ir3,ngfftf(3))
236        wp_idx = 1 + wp1 + wp2*ngfftf(1) + wp3*ngfftf(1)*ngfftf(2)
237        rcl2fft(irc)  = wp_idx
238        rclred(1,irc) = DBLE(ir1)/ngfftf(1) ! Reduced coordinates in terms of the origina cell.
239        rclred(2,irc) = DBLE(ir2)/ngfftf(2)
240        rclred(3,irc) = DBLE(ir3)/ngfftf(3)
241      end do
242    end do
243  end do
244  !
245  ! Wrap the point to [0,1[ where 1 is not included (tol12)
246  call wrap2_zero_one(eh_rcoord(1),eh_red(1),eh_shift(1))
247  call wrap2_zero_one(eh_rcoord(2),eh_red(2),eh_shift(2))
248  call wrap2_zero_one(eh_rcoord(3),eh_red(3),eh_shift(3))
249  write(std_out,*)"Initial Position of (e|h) in reduced coordinates:",eh_red," shift: ",eh_shift
250  !
251  ! Initial position on the FFT grid (the closest one)
252  eh_red(:)=NINT(eh_red(:)*ngfftf(1:3))
253  !
254  ! Not translate the (electron|hole) in the center of the big box.
255  do ii=1,3
256    if (nrcl(ii)/=1) eh_red(ii) = eh_red(ii) + (nrcl(ii)/2) * ngfftf(ii)
257  end do
258  write(std_out,*)"Position of (e|h) in the center of the big box in FFT coordinates:",eh_red(:)
259 
260  eh_rr = 1 + eh_red(1) + eh_red(2)*bbox(1) + eh_red(3)*bbox(1)*bbox(2)
261  eh_fft_idx = rcl2fft(eh_rr)
262 
263  write(std_out,*)" Reduced coordinates of (e|h) in the center of the big box ",rclred(:,eh_rr)
264  !
265  ! Allocate wavefunctions.
266  ABI_MALLOC(ur_v,(Wfd%nfft*nspinor))
267  ABI_MALLOC(ur_c,(Wfd%nfft*nspinor))
268 
269  if (usepaw==1) then
270    ABI_DT_MALLOC(Cp_v,(Wfd%natom,nspinor))
271    call pawcprj_alloc(Cp_v,0,Wfd%nlmn_atm)
272    ABI_DT_MALLOC(Cp_c,(Wfd%natom,nspinor))
273    call pawcprj_alloc(Cp_c,0,Wfd%nlmn_atm)
274  end if
275  !
276  ! Read excitonic eigenvector.
277  hsize = SUM(BSp%nreh); if (BSp%use_coupling>0) hsize=2*hsize
278  nvec=1
279 
280  ABI_STAT_MALLOC(vec_list,(hsize,nvec), ierr)
281  ABI_CHECK(ierr==0, "out of memory in vec_list")
282 
283  ABI_MALLOC(vec_idx,(nvec))
284  vec_idx = (/(ii, ii=1,nvec)/)
285 
286  if (BS_files%in_eig /= BSE_NOFILE) then
287    eig_fname = BS_files%in_eig
288  else
289    eig_fname = BS_files%out_eig
290  end if
291 
292  if (my_rank==master) then
293    call exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,Bsp=Bsp)
294  end if
295 
296  call xmpi_bcast(vec_list,master,comm,ierr)
297 
298  ABI_FREE(vec_idx)
299  !
300  ! Allocate the excitonic wavefunction on the big box.
301  ABI_MALLOC(exc_phi,(nr_tot))
302  exc_phi=czero
303  !
304  !got=0;
305  !bbp_mask=.FALSE.; bbp_mask(minb:maxb,minb:maxb)=.TRUE.
306  spin_start=1; spin_stop=Wfd%nsppol
307  if (spin_opt==1) spin_stop =1
308  if (spin_opt==2) spin_start=2
309 
310  ! TODO
311  ! Distribute (b,b',k,s) states where k is in the *full* BZ.
312  !call wfd_distribute_bands(Wfd,ik_ibz,spin,my_nband,my_band_list,got,bmask)
313 
314  do spin=spin_start,spin_stop
315    do reh=1,Bsp%nreh(spin)
316      ene_rt = Bsp%Trans(reh,spin)%en
317      ik_bz  = Bsp%Trans(reh,spin)%k
318      val    = Bsp%Trans(reh,spin)%v
319      cond   = Bsp%Trans(reh,spin)%c
320      k_bz   = Kmesh%bz(:,ik_bz)
321      !
322      ! The resonant component.
323      rt_idx = reh + (spin-1)*Bsp%nreh(1)
324      res_coeff = vec_list(rt_idx,1)
325      !
326      ! The anti-resonant component.
327      antres_coeff = czero
328      if (Bsp%use_coupling>0) then
329        art_idx  = reh + (spin-1)*Bsp%nreh(1) + SUM(Bsp%nreh)
330        antres_coeff = vec_list(art_idx,1)
331      end if
332 
333      call wfd_sym_ur(Wfd,Cryst,Kmesh,val ,ik_bz,spin,ur_v) ! TODO recheck this one.
334      call wfd_sym_ur(Wfd,Cryst,Kmesh,cond,ik_bz,spin,ur_c)
335 
336      !call wfd_paw_get_aeur(Wfd,band,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae,ur_ae_onsite,ur_ps_onsite)
337 
338      if (which_fixed==1) then ! electron
339        do rr=1,nr_tot
340          wp_idx = rcl2fft(rr)
341          r12 = eh_red - rclred(:,rr)
342          k_dot_r12 = two_pi * DOT_PRODUCT(k_bz,r12)
343          eikr12 = DCMPLX(COS(k_dot_r12), SIN(k_dot_r12))
344          exc_phi(rr) = exc_phi(rr) + eikr12 * ur_v(wp_idx) * CONJG(ur_c(wp_idx))
345        end do
346      else ! hole
347        MSG_ERROR("Not coded")
348      end if
349      !
350   end do
351  end do
352 
353  ABI_FREE(vec_list)
354  ABI_FREE(rcl2fft)
355  ABI_FREE(rclred)
356  ABI_FREE(ur_v)
357  ABI_FREE(ur_c)
358  !
359  if (Wfd%usepaw==1) then
360    call pawcprj_free(Cp_v)
361    ABI_DT_FREE(Cp_v)
362    call pawcprj_free(Cp_c)
363    ABI_DT_FREE(Cp_c)
364    if (paw_add_onsite) then
365      call pawfgrtab_free(Pawfgrtab)
366      ABI_DT_FREE(Pawfgrtab)
367      call paw_pwaves_lmn_free(Paw_onsite)
368      ABI_DT_FREE(Paw_onsite)
369    end if
370  end if
371  !
372  ! Collect results on master node.
373  call xmpi_sum_master(exc_phi,master,comm,ierr)
374  !
375  ! The exciton density probability.
376  ABI_MALLOC(exc_phi2,(nr_tot))
377 
378  do rr=1,nr_tot
379    exc_phi2(rr) = ABS(exc_phi(rr))**2
380  end do
381  !
382  ! Construct the supercell.
383  sc_natom = ncells*Cryst%natom
384  ABI_MALLOC(sc_xcart,(3,sc_natom))
385  ABI_MALLOC(sc_typat,(sc_natom))
386 
387  ii=0
388  do ir3=0,nrcl(3)-1  ! Loop over the replicaed cells.
389    do ir2=0,nrcl(2)-1
390      do ir1=0,nrcl(1)-1
391       !
392       sc_rprimd(:,1) = ir1 * Cryst%rprimd(:,1) ! Workspace.
393       sc_rprimd(:,2) = ir2 * Cryst%rprimd(:,2)
394       sc_rprimd(:,3) = ir3 * Cryst%rprimd(:,3)
395       scart_shift = SUM(sc_rprimd,DIM=2)
396       do iatom=1,Cryst%natom
397         ii=ii+1
398         sc_xcart(:,ii) = Cryst%xcart(:,iatom) + scart_shift
399         sc_typat(ii)   = Cryst%typat(iatom)
400       end do
401       !
402     end do
403    end do
404  end do
405  !
406  ! Lattice vectors of the big box.
407  do ii=1,3
408    sc_rprimd(:,ii) = nrcl(ii) * Cryst%rprimd(:,ii)
409  end do
410  !
411  ! Master writes the data.
412  if (my_rank==master) then
413    out_fname = TRIM(BS_files%out_basename)//"_EXC_WF"
414    if (open_file(out_fname,msg,newunit=xsf_unt,form='formatted') /= 0) then
415      MSG_ERROR(msg)
416    end if
417 
418    call printxsf(bbox(1),bbox(2),bbox(3),exc_phi2,sc_rprimd,origin0,sc_natom,Cryst%ntypat,sc_typat,&
419 &    sc_xcart,Cryst%znucl,xsf_unt,0)
420 
421    close(xsf_unt)
422  end if
423 
424  ABI_FREE(sc_xcart)
425  ABI_FREE(sc_typat)
426 
427  ABI_FREE(exc_phi)
428  ABI_FREE(exc_phi2)
429 
430 end subroutine exc_plot