TABLE OF CONTENTS


ABINIT/exc_plot [ Functions ]

[ Top ] [ Functions ]

NAME

  exc_plot

FUNCTION

  Plots the excitonic wavefunction in real space.

COPYRIGHT

 Copyright (C) 2009-2018 ABINIT group (MG)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

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