TABLE OF CONTENTS
ABINIT/exc_plot [ 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