TABLE OF CONTENTS
ABINIT/m_exc_analyze [ 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-2024 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 .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_exc_analyze 24 25 use defs_basis 26 use m_abicore 27 use m_bs_defs 28 use m_xmpi 29 use m_errors 30 31 use defs_datatypes, only : pseudopotential_type, ebands_t 32 use m_io_tools, only : open_file 33 use m_numeric_tools, only : iseven, wrap2_zero_one 34 use m_bz_mesh, only : kmesh_t 35 use m_crystal, only : crystal_t 36 use m_wfd, only : wfdgw_t 37 use m_bse_io, only : exc_read_eigen 38 use m_pptools, only : printxsf 39 use m_pawrad, only : pawrad_type 40 use m_pawtab, only : pawtab_type,pawtab_get_lsize 41 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print 42 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 43 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 44 use m_paw_nhat, only : nhatgrid 45 46 implicit none 47 48 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
SOURCE
423 subroutine exc_den(BSp,BS_files,ngfft,nfftot,Kmesh,ktabr,Wfd) 424 425 !Arguments ------------------------------------ 426 !scalars 427 integer,intent(in) :: nfftot 428 type(excparam),intent(in) :: BSp 429 type(excfiles),intent(in) :: BS_files 430 type(kmesh_t),intent(in) :: Kmesh 431 type(wfdgw_t),intent(inout) :: Wfd 432 !arrays 433 integer,intent(in) :: ngfft(18) 434 integer,intent(in) :: ktabr(nfftot,BSp%nkbz) 435 436 !Local variables ------------------------------ 437 !scalars 438 integer :: nfft1,nfft2,nfft3,spin,ierr 439 integer :: it,itp,ik_ibz,ikp_bz,ik_bz,band,iv,ivp,ic,icp,ir,i1,i2,i3,iik 440 integer :: state,nstates,min_idx,eig_unt,den_unt,sden_unt,hsize,homo 441 real(dp) :: min_ene,cost 442 character(len=500) :: msg 443 character(len=fnlen) :: filbseig 444 !arrays 445 real(dp) :: n0(nfftot),rho_eh(nfftot),nexc(nfftot) 446 real(dp),allocatable :: exc_ene(:) 447 complex(dpc) :: rhor_h(nfftot),rhor_e(nfftot) 448 complex(gwpc),allocatable :: wfr(:,:,:), wfrk(:,:) 449 complex(dpc),allocatable :: exc_ene_cplx(:),exc_vec(:) 450 451 !************************************************************************ 452 453 ABI_CHECK(Wfd%nsppol==1,"nsppol==2 not coded") 454 455 if (BS_files%in_eig /= BSE_NOFILE) then 456 filbseig = BS_files%in_eig 457 else 458 filbseig = BS_files%out_eig 459 end if 460 461 call wrtout(std_out,' Calculating electron-hole, excited state density',"COLL") 462 call wrtout(std_out," Reading eigenstates from: "//TRIM(filbseig),"COLL") 463 ABI_ERROR("Not tested") 464 465 ! Prepare FFT tables to have u(r) on the ngfft_osc mesh. 466 !if ( ANY(ngfftf(1:3) /= Wfd%ngfft(1:3)) ) then 467 ! call wfd%change_ngfft(Cryst,Psps,ngfftf) 468 !end if 469 470 nfft1 = ngfft(1) 471 nfft2 = ngfft(2) 472 nfft3 = ngfft(3) 473 ABI_CHECK(nfftot==PRODUCT(ngfft(1:3)),"Mismatch in FFT size") 474 475 !allocate and load wavefunctions in real space 476 ABI_MALLOC_OR_DIE(wfr,(nfftot*Wfd%nspinor,BSp%nbnds,Wfd%nkibz), ierr) 477 478 spin=1 ! SPIN support is missing. 479 480 do ik_ibz=1,Wfd%nkibz 481 do band=1,BSp%nbnds 482 call wfd%get_ur(band,ik_ibz,spin,wfr(:,band,ik_ibz)) 483 end do 484 end do 485 486 if (open_file(filbseig,msg,newunit=eig_unt,form="unformatted", status="old", action="read") /= 0) then 487 ABI_ERROR(msg) 488 end if 489 490 read(eig_unt) hsize,nstates 491 492 if (nstates/=Bsp%nreh(1)) then 493 write(std_out,*) 'not resonant calculation' 494 close(eig_unt) 495 RETURN 496 end if 497 498 ABI_MALLOC(exc_ene_cplx,(nstates)) 499 ABI_MALLOC(exc_ene,(nstates)) 500 501 read(eig_unt) exc_ene_cplx(1:nstates) 502 503 ! Small imaginary part is always neglected. 504 exc_ene(:) = exc_ene_cplx(:) 505 ABI_FREE(exc_ene_cplx) 506 ! 507 ! Find the lowest non-negative eigenvalue. 508 min_ene = greatest_real 509 do state=1,nstates 510 if (exc_ene(state) < zero) cycle 511 if (exc_ene(state) < min_ene) then 512 min_ene = exc_ene(state) 513 min_idx = state 514 end if 515 end do 516 ABI_FREE(exc_ene) 517 518 write(std_out,*)' Lowest eigenvalue ', min_idx, min_ene*Ha_eV 519 ! 520 ! skip other eigenvectors 521 do state=1,min_idx-1 522 read(eig_unt) 523 end do 524 ! 525 ! read "lowest" eigenvector 526 ABI_MALLOC(exc_vec,(hsize)) 527 read(eig_unt) exc_vec 528 529 close(eig_unt) 530 531 ABI_MALLOC_OR_DIE(wfrk,(nfftot,BSp%nbnds), ierr) 532 533 !calculate ground state density 534 n0(:) = zero 535 spin = 1 536 homo = Bsp%homo_spin(spin) 537 do ik_bz = 1, BSp%nkbz 538 ik_ibz = Kmesh%tab(ik_bz) 539 iik = (3-Kmesh%tabi(ik_bz))/2 540 541 if (iik==1) then 542 do ir=1,nfftot 543 wfrk(ir,1:homo) = wfr(ktabr(ir,ik_bz),1:homo,ik_ibz) 544 end do 545 else 546 do ir=1,nfftot 547 wfrk(ir,1:homo) = conjg(wfr(ktabr(ir,ik_bz),1:homo,ik_ibz)) 548 end do 549 end if 550 551 do iv=1,homo 552 n0(:) = n0(:) + conjg(wfrk(:,band)) * wfrk(:,band) 553 end do 554 end do 555 ! 556 ! calculate electron and hole density 557 rhor_h=czero 558 rhor_e=czero 559 ABI_CHECK(BSp%nsppol==1,"nsppol=2 not coded") 560 561 do it=1,Bsp%nreh(1) 562 ik_bz = Bsp%Trans(it,1)%k 563 iv = Bsp%Trans(it,1)%v 564 ic = Bsp%Trans(it,1)%c 565 ik_ibz = Kmesh%tab(ik_bz) 566 iik = (3-Kmesh%tabi(ik_bz))/2 567 568 if (iik==1) then 569 do ir = 1, nfftot 570 wfrk(ir,:) = wfr(ktabr(ir,ik_bz),:,ik_ibz) 571 end do 572 else 573 do ir = 1, nfftot 574 wfrk(ir,:) = conjg(wfr(ktabr(ir,ik_bz),:,ik_ibz)) 575 end do 576 end if 577 ! 578 do itp=1,Bsp%nreh(1) 579 ikp_bz = Bsp%Trans(itp,1)%k 580 if (ikp_bz/=ik_bz) CYCLE 581 icp = Bsp%Trans(itp,1)%c 582 ivp = Bsp%Trans(itp,1)%v 583 if(icp==ic) then 584 rhor_h = rhor_h - conjg(exc_vec(it)) * exc_vec(itp) * wfrk(:,iv) * conjg(wfrk(:,ivp)) 585 end if 586 if(ivp==iv) then 587 rhor_e = rhor_e + conjg(exc_vec(it)) * exc_vec(itp) * conjg(wfrk(:,ic)) * wfrk(:,icp) 588 end if 589 end do 590 end do 591 592 ABI_FREE(exc_vec) 593 ! 594 ! calculate excited state density minus ground state density 595 ! n* - n0 = rhor_e + rhor_h 596 rho_eh(:) = rhor_e + rhor_h 597 ! 598 !calculate excited state density 599 ! n* = n0 + rhor_e + rhor_h 600 nexc(:) = n0(:) + rho_eh(:) 601 602 if (open_file('out.den',msg,newunit=den_unt) /= 0) then 603 ABI_ERROR(msg) 604 end if 605 606 ! here conversion to cartesian through a1, a2, a3 607 cost = zero 608 do i1=0,nfft1-1 609 do i2=0,nfft2-1 610 do i3=0,nfft3-1 611 ir=1 + i1 + i2*nfft1 + i3*nfft1*nfft2 612 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)) 613 cost = cost + n0(ir) 614 end do 615 end do 616 end do 617 618 write(std_out,*) 'density normalization constant ', cost 619 close(den_unt) 620 621 if (open_file('out.sden',msg,newunit=sden_unt) /= 0) then 622 ABI_ERROR(msg) 623 end if 624 625 ! we are looking for the plane between (100) (111) 626 ! so it's the place where (111) [ (100) x v ] = 0 (mixed product) 627 ! finally v2 - v3 = 0 628 do i2=0, nfft2-1 629 do i3=0, nfft3-1 630 if(i2==i3) then 631 write(std_out,*) i2, i3 632 write(sden_unt,*) (rho_eh(1+i1+i2*nfft1+i3*nfft1*nfft2),i1=0,nfft1-1) 633 end if 634 end do 635 end do 636 637 close(sden_unt) 638 639 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<wfdgw_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
SOURCE
86 subroutine exc_plot(Bsp,Bs_files,Wfd,Kmesh,Cryst,Psps,Pawtab,Pawrad,paw_add_onsite,spin_opt,which_fixed,eh_rcoord,nrcell,ngfftf) 87 88 !Arguments ------------------------------------ 89 !scalars 90 integer,intent(in) :: which_fixed,spin_opt 91 logical,intent(in) :: paw_add_onsite 92 type(excparam),intent(in) :: Bsp 93 type(excfiles),intent(in) :: BS_files 94 type(kmesh_t),intent(in) :: Kmesh 95 type(crystal_t),intent(in) :: Cryst 96 type(pseudopotential_type),intent(in) :: Psps 97 type(wfdgw_t),intent(inout) :: Wfd 98 !arrays 99 integer,intent(in) :: ngfftf(18),nrcell(3) 100 real(dp),intent(in) :: eh_rcoord(3) 101 type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 102 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*Wfd%usepaw) 103 104 !Local variables ------------------------------ 105 !scalars 106 integer,parameter :: cplex1=1 107 integer :: nsppol,usepaw,nspinor,comm,cond,val 108 integer :: optcut,optgr0,optgr1,optgr2,optrad 109 integer :: ik_bz,ierr,my_rank !ik_ibz, istwf_k, isym_k,itim_k, my_nbbp, npw_k, 110 integer :: spin,spin_start,spin_stop,reh 111 integer :: rt_idx,art_idx,ii,iatom,nr_tot 112 integer :: irc,ir1,ir2,ir3,wp1,wp2,wp3,wp_idx,eh_fft_idx,eh_rr,rr 113 integer :: hsize,xsf_unt,ncells,nvec 114 integer :: sc_natom,master 115 real(dp) :: ene_rt,k_dot_r12 116 complex(dpc) :: res_coeff,antres_coeff,eikr12 117 character(len=500) :: msg 118 character(len=fnlen) :: eig_fname,out_fname 119 !arrays 120 integer :: nrcl(3),bbox(3) 121 !integer :: bbp_distrb(Wfd%mband,Wfd%mband),got(Wfd%nproc) 122 integer,allocatable :: rcl2fft(:),sc_typat(:),l_size_atm(:),vec_idx(:) 123 real(dp),parameter :: origin0(3)=(/zero,zero,zero/) 124 real(dp) :: k_bz(3),eh_red(3),eh_shift(3),r12(3),sc_rprimd(3,3),scart_shift(3) 125 real(dp),allocatable :: rclred(:,:),exc_phi2(:),sc_xcart(:,:) !,sc_znucl(:) 126 complex(gwpc),allocatable :: exc_phi(:) 127 complex(gwpc),allocatable :: ur_v(:),ur_c(:) 128 complex(dpc),allocatable :: vec_list(:,:) 129 !logical :: bbp_mask(Wfd%mband,Wfd%mband) 130 type(pawcprj_type),allocatable :: Cp_v(:,:),Cp_c(:,:) 131 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 132 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 133 134 !************************************************************************ 135 136 ABI_WARNING("Exc plot is still under development") 137 138 call wrtout(std_out," Calculating excitonic wavefunctions in real space","COLL") 139 ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded") 140 ABI_CHECK(Wfd%usepaw==0,"PAW not coded") 141 ABI_CHECK(Wfd%nproc==1,"nproc>1 not supported") 142 143 ! If needed, prepare FFT tables to have u(r) on the ngfftf mesh. 144 if ( ANY(ngfftf(1:3) /= Wfd%ngfft(1:3)) ) then 145 call wfd%change_ngfft(Cryst,Psps,ngfftf) 146 end if 147 148 comm = Wfd%comm 149 my_rank = Wfd%my_rank 150 master = Wfd%master 151 152 nsppol = Wfd%nsppol 153 nspinor = Wfd%nspinor 154 usepaw = Wfd%usepaw 155 ! 156 ! TODO recheck this part. 157 ! Prepare tables needed for splining the onsite term on the FFT box. 158 if (Wfd%usepaw==1.and.paw_add_onsite) then 159 ABI_WARNING("Testing the calculation of AE PAW wavefunctions.") 160 ! Use a local pawfgrtab to make sure we use the correction in the paw spheres 161 ! the usual pawfgrtab uses r_shape which may not be the same as r_paw. 162 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 163 ABI_MALLOC(Pawfgrtab,(Cryst%natom)) 164 call pawfgrtab_init(Pawfgrtab,cplex1,l_size_atm,Wfd%nspden,Cryst%typat) 165 ABI_FREE(l_size_atm) 166 167 optcut=1 ! use rpaw to construct Pawfgrtab. 168 optgr0=0; optgr1=0; optgr2=0 ! dont need gY terms locally. 169 optrad=1 ! do store r-R. 170 171 call nhatgrid(Cryst%atindx1,Cryst%gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 172 & optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 173 !Pawfgrtab is ready to use 174 175 if (Wfd%pawprtvol>0) then 176 call pawfgrtab_print(Pawfgrtab,natom=Cryst%natom,unit=std_out,& 177 & prtvol=Wfd%pawprtvol,mode_paral="COLL") 178 end if 179 180 ABI_MALLOC(Paw_onsite,(Cryst%natom)) 181 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,Cryst%rprimd,Cryst%xcart,& 182 & Pawtab,Pawrad,Pawfgrtab) 183 end if 184 ! 185 ! Number of cells in the big box. Odd numbers are needed to place the e or the h in the middle. 186 do ii=1,3 187 nrcl(ii) = nrcell(ii) 188 if (iseven(nrcell(ii))) then 189 nrcl(ii) = nrcell(ii)+1 190 write(msg,'(2(a,i0))')" Enforcing odd number of cell replicas ",nrcell(ii)," --> ",nrcl(ii) 191 ABI_WARNING(msg) 192 end if 193 end do 194 195 ncells = PRODUCT(nrcl) 196 nr_tot = Wfd%nfftot * ncells ! Total number of points in the big box. 197 bbox = nrcl*ngfftf(1:3) 198 ! 199 ! rcl2fft: The image of the point in the small box. 200 ! rcl2red: The reduced coordinates of the point in the big box in terms of rprimd. 201 ABI_MALLOC(rcl2fft,(nr_tot)) 202 ABI_MALLOC(rclred,(3,nr_tot)) 203 204 irc = 0 205 do ir3=0,bbox(3)-1 ! Loop over the points in the big box. 206 do ir2=0,bbox(2)-1 207 do ir1=0,bbox(1)-1 208 irc = 1+irc 209 ! 210 wp1=MODULO(ir1,ngfftf(1)) ! The FFT index of the point wrapped into the original cell. 211 wp2=MODULO(ir2,ngfftf(2)) 212 wp3=MODULO(ir3,ngfftf(3)) 213 wp_idx = 1 + wp1 + wp2*ngfftf(1) + wp3*ngfftf(1)*ngfftf(2) 214 rcl2fft(irc) = wp_idx 215 rclred(1,irc) = DBLE(ir1)/ngfftf(1) ! Reduced coordinates in terms of the origina cell. 216 rclred(2,irc) = DBLE(ir2)/ngfftf(2) 217 rclred(3,irc) = DBLE(ir3)/ngfftf(3) 218 end do 219 end do 220 end do 221 ! 222 ! Wrap the point to [0,1[ where 1 is not included (tol12) 223 call wrap2_zero_one(eh_rcoord(1),eh_red(1),eh_shift(1)) 224 call wrap2_zero_one(eh_rcoord(2),eh_red(2),eh_shift(2)) 225 call wrap2_zero_one(eh_rcoord(3),eh_red(3),eh_shift(3)) 226 write(std_out,*)"Initial Position of (e|h) in reduced coordinates:",eh_red," shift: ",eh_shift 227 ! 228 ! Initial position on the FFT grid (the closest one) 229 eh_red(:)=NINT(eh_red(:)*ngfftf(1:3)) 230 ! 231 ! Not translate the (electron|hole) in the center of the big box. 232 do ii=1,3 233 if (nrcl(ii)/=1) eh_red(ii) = eh_red(ii) + (nrcl(ii)/2) * ngfftf(ii) 234 end do 235 write(std_out,*)"Position of (e|h) in the center of the big box in FFT coordinates:",eh_red(:) 236 237 eh_rr = 1 + eh_red(1) + eh_red(2)*bbox(1) + eh_red(3)*bbox(1)*bbox(2) 238 eh_fft_idx = rcl2fft(eh_rr) 239 240 write(std_out,*)" Reduced coordinates of (e|h) in the center of the big box ",rclred(:,eh_rr) 241 ! 242 ! Allocate wavefunctions. 243 ABI_MALLOC(ur_v,(Wfd%nfft*nspinor)) 244 ABI_MALLOC(ur_c,(Wfd%nfft*nspinor)) 245 246 if (usepaw==1) then 247 ABI_MALLOC(Cp_v,(Wfd%natom,nspinor)) 248 call pawcprj_alloc(Cp_v,0,Wfd%nlmn_atm) 249 ABI_MALLOC(Cp_c,(Wfd%natom,nspinor)) 250 call pawcprj_alloc(Cp_c,0,Wfd%nlmn_atm) 251 end if 252 ! 253 ! Read excitonic eigenvector. 254 hsize = SUM(BSp%nreh); if (BSp%use_coupling>0) hsize=2*hsize 255 nvec=1 256 257 ABI_MALLOC_OR_DIE(vec_list,(hsize,nvec), ierr) 258 259 ABI_MALLOC(vec_idx,(nvec)) 260 vec_idx = (/(ii, ii=1,nvec)/) 261 262 if (BS_files%in_eig /= BSE_NOFILE) then 263 eig_fname = BS_files%in_eig 264 else 265 eig_fname = BS_files%out_eig 266 end if 267 268 if (my_rank==master) then 269 call exc_read_eigen(eig_fname,hsize,nvec,vec_idx,vec_list,Bsp=Bsp) 270 end if 271 272 call xmpi_bcast(vec_list,master,comm,ierr) 273 274 ABI_FREE(vec_idx) 275 ! 276 ! Allocate the excitonic wavefunction on the big box. 277 ABI_MALLOC(exc_phi,(nr_tot)) 278 exc_phi=czero 279 ! 280 !got=0; 281 !bbp_mask=.FALSE.; bbp_mask(minb:maxb,minb:maxb)=.TRUE. 282 spin_start=1; spin_stop=Wfd%nsppol 283 if (spin_opt==1) spin_stop =1 284 if (spin_opt==2) spin_start=2 285 286 ! TODO 287 ! Distribute (b,b',k,s) states where k is in the *full* BZ. 288 !call wfd_distribute_bands(Wfd,ik_ibz,spin,my_nband,my_band_list,got,bmask) 289 290 do spin=spin_start,spin_stop 291 do reh=1,Bsp%nreh(spin) 292 ene_rt = Bsp%Trans(reh,spin)%en 293 ik_bz = Bsp%Trans(reh,spin)%k 294 val = Bsp%Trans(reh,spin)%v 295 cond = Bsp%Trans(reh,spin)%c 296 k_bz = Kmesh%bz(:,ik_bz) 297 ! 298 ! The resonant component. 299 rt_idx = reh + (spin-1)*Bsp%nreh(1) 300 res_coeff = vec_list(rt_idx,1) 301 ! 302 ! The anti-resonant component. 303 antres_coeff = czero 304 if (Bsp%use_coupling>0) then 305 art_idx = reh + (spin-1)*Bsp%nreh(1) + SUM(Bsp%nreh) 306 antres_coeff = vec_list(art_idx,1) 307 end if 308 309 call wfd%sym_ur(Cryst,Kmesh,val ,ik_bz,spin,ur_v) ! TODO recheck this one. 310 call wfd%sym_ur(Cryst,Kmesh,cond,ik_bz,spin,ur_c) 311 312 !call wfd_paw_get_aeur(Wfd,band,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae,ur_ae_onsite,ur_ps_onsite) 313 314 if (which_fixed==1) then ! electron 315 do rr=1,nr_tot 316 wp_idx = rcl2fft(rr) 317 r12 = eh_red - rclred(:,rr) 318 k_dot_r12 = two_pi * DOT_PRODUCT(k_bz,r12) 319 eikr12 = DCMPLX(COS(k_dot_r12), SIN(k_dot_r12)) 320 exc_phi(rr) = exc_phi(rr) + eikr12 * ur_v(wp_idx) * CONJG(ur_c(wp_idx)) 321 end do 322 else ! hole 323 ABI_ERROR("Not coded") 324 end if 325 ! 326 end do 327 end do 328 329 ABI_FREE(vec_list) 330 ABI_FREE(rcl2fft) 331 ABI_FREE(rclred) 332 ABI_FREE(ur_v) 333 ABI_FREE(ur_c) 334 ! 335 if (Wfd%usepaw==1) then 336 call pawcprj_free(Cp_v) 337 ABI_FREE(Cp_v) 338 call pawcprj_free(Cp_c) 339 ABI_FREE(Cp_c) 340 if (paw_add_onsite) then 341 call pawfgrtab_free(Pawfgrtab) 342 ABI_FREE(Pawfgrtab) 343 call paw_pwaves_lmn_free(Paw_onsite) 344 ABI_FREE(Paw_onsite) 345 end if 346 end if 347 ! 348 ! Collect results on master node. 349 call xmpi_sum_master(exc_phi,master,comm,ierr) 350 ! 351 ! The exciton density probability. 352 ABI_MALLOC(exc_phi2,(nr_tot)) 353 354 do rr=1,nr_tot 355 exc_phi2(rr) = ABS(exc_phi(rr))**2 356 end do 357 ! 358 ! Construct the supercell. 359 sc_natom = ncells*Cryst%natom 360 ABI_MALLOC(sc_xcart,(3,sc_natom)) 361 ABI_MALLOC(sc_typat,(sc_natom)) 362 363 ii=0 364 do ir3=0,nrcl(3)-1 ! Loop over the replicaed cells. 365 do ir2=0,nrcl(2)-1 366 do ir1=0,nrcl(1)-1 367 ! 368 sc_rprimd(:,1) = ir1 * Cryst%rprimd(:,1) ! Workspace. 369 sc_rprimd(:,2) = ir2 * Cryst%rprimd(:,2) 370 sc_rprimd(:,3) = ir3 * Cryst%rprimd(:,3) 371 scart_shift = SUM(sc_rprimd,DIM=2) 372 do iatom=1,Cryst%natom 373 ii=ii+1 374 sc_xcart(:,ii) = Cryst%xcart(:,iatom) + scart_shift 375 sc_typat(ii) = Cryst%typat(iatom) 376 end do 377 ! 378 end do 379 end do 380 end do 381 ! 382 ! Lattice vectors of the big box. 383 do ii=1,3 384 sc_rprimd(:,ii) = nrcl(ii) * Cryst%rprimd(:,ii) 385 end do 386 ! 387 ! Master writes the data. 388 if (my_rank==master) then 389 out_fname = TRIM(BS_files%out_basename)//"_EXC_WF" 390 if (open_file(out_fname,msg,newunit=xsf_unt,form='formatted') /= 0) then 391 ABI_ERROR(msg) 392 end if 393 394 call printxsf(bbox(1),bbox(2),bbox(3),exc_phi2,sc_rprimd,origin0,sc_natom,Cryst%ntypat,sc_typat,& 395 & sc_xcart,Cryst%znucl,xsf_unt,0) 396 397 close(xsf_unt) 398 end if 399 400 ABI_FREE(sc_xcart) 401 ABI_FREE(sc_typat) 402 403 ABI_FREE(exc_phi) 404 ABI_FREE(exc_phi2) 405 406 end subroutine exc_plot