TABLE OF CONTENTS


ABINIT/m_exc_spectra [ Modules ]

[ Top ] [ Modules ]

NAME

 m_exc_spectra

FUNCTION

  Routines to compute the macroscopic dielectric function in the Bethe-Salpeter code.

COPYRIGHT

 Copyright (C) 2009-2018 ABINIT and EXC groups (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi, Y. Gillet)
  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

SOURCE

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 MODULE m_exc_spectra
26 
27  use defs_basis
28  use defs_datatypes
29  use m_bs_defs
30  use m_abicore
31  use iso_c_binding
32  use m_xmpi
33  use m_errors
34  use m_nctk
35 #ifdef HAVE_NETCDF
36  use netcdf
37 #endif
38  use m_ebands
39  !use m_hdr,          only : hdr_free
40 
41  use defs_abitypes,     only : hdr_type
42  use m_io_tools,        only : open_file
43  use m_fstrings,        only : toupper, strcat, sjoin, int2char4
44  use m_numeric_tools,   only : simpson_int, simpson_cplx
45  use m_hide_blas,       only : xdotu,xdotc
46  use m_special_funcs,   only : dirac_delta
47  use m_crystal,         only : crystal_t
48  use m_crystal_io,      only : crystal_ncwrite
49  use m_bz_mesh,         only : kmesh_t
50  use m_eprenorms,       only : eprenorms_t, renorm_bst
51  use m_pawtab,          only : pawtab_type
52  use m_paw_hr,          only : pawhur_t
53  use m_wfd,             only : wfd_t
54  !use m_bse_io,          only : exc_amplitude
55  use m_wfd_optic,       only : calc_optical_mels
56 
57  implicit none
58 
59  private
60 
61  public :: build_spectra           ! Driver routine for the computation of optical spectra.
62  public :: exc_write_data          ! This routine drives the writing of the files produced by the Bethe-Salpeter code.
63  public :: exc_eps_rpa             ! Build epsilon within RPA and GW.
64  public :: mdfs_ncwrite            ! Writes the MDF.nc file with the final results.
65  public :: exc_write_tensor        ! Write of complex dielectric tensor
66  !public :: exc_eps_resonant       ! Build the macroscopic dielectric function with excitonic effects.

m_exc_spectra/build_spectra [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  build_spectra

FUNCTION

  Driver routine for the computation of optical spectra.

INPUTS

  usepaw=1 for PAW calculations, 0 otherwise.
  drude_plsmf=Drude plasma frequency.
  Bsp<excparam>=Data type gathering the paramenters used for the Bethe-Salpeter calculation.
    inclvkb=If different from 0, [Vnl,r] is included in the calculation of the matrix elements of the velocity operator.
  BS_files<excfiles>=filenames used in the Bethe-Salpeter part.
  Kmesh<kmesh_t>=the k-point sampling for the wave functions.
  Cryst<crystal_t>=Structure defining the crystalline structure.
  KS_BSt=The KS energies.
  QP_BSt=The QP energies.
  Psps <pseudopotential_type>=variables related to pseudopotentials.
  Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data
  Hur(Cryst%natom*usepaw)<pawhur_t>=Only for PAW and LDA+U, quantities used to evaluate the commutator [H_u,r].
  Wfd<wfd_t>=Handler for the wavefunctions.
    nsppol=Number of independent spin polarizations.
    nspinor=Number of spinorial components.
  comm=MPI communicator.

OUTPUT

  No output. The routine calls specialized routines where the computation and the output of the spectra is done.

PARENTS

      m_exc_diago

CHILDREN

      c_f_pointer

SOURCE

107 subroutine build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm,Epren)
108 
109 
110 !This section has been created automatically by the script Abilint (TD).
111 !Do not modify the following lines by hand.
112 #undef ABI_FUNC
113 #define ABI_FUNC 'build_spectra'
114 !End of the abilint section
115 
116  implicit none
117 
118 !Arguments ------------------------------------
119 !scalars
120  integer,intent(in) :: comm
121  real(dp),intent(in) :: drude_plsmf
122  type(excparam),intent(in) :: BSp
123  type(excfiles),intent(in) :: BS_files
124  type(pseudopotential_type),intent(in) :: Psps
125  type(kmesh_t),intent(in) :: Kmesh
126  type(crystal_t),intent(in) :: Cryst
127  type(ebands_t),intent(in) :: KS_BSt,QP_BSt
128  type(wfd_t),intent(inout) :: Wfd
129  type(eprenorms_t),optional,intent(in) :: Epren
130 !arrays
131  type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
132  type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
133 
134 !Local variables ------------------------------
135 !scalars
136  integer :: my_rank,master,iq,io,nsppol,lomo_min,max_band,ncid
137  integer :: itemp,ntemp
138  logical :: do_ep_renorm
139  real(dp) :: omegaev
140  complex(dpc) :: ks_avg,gw_avg,exc_avg
141  character(len=4) :: ts
142  character(len=fnlen) :: path,prefix
143  character(len=fnlen) :: filbseig, ost_fname
144  !character(len=500) :: msg
145  type(ebands_t) :: EPBSt, EP_QPBSt
146 !arrays
147  real(dp),allocatable :: dos_exc(:),dos_gw(:),dos_ks(:)
148  complex(dpc),allocatable :: eps_rpanlf(:,:),eps_gwnlf(:,:)
149  complex(dpc),allocatable :: eps_exc(:,:),opt_cvk(:,:,:,:,:)
150 
151 !************************************************************************
152 
153  my_rank = Wfd%my_rank
154  master  = Wfd%master
155  nsppol  = Wfd%nsppol
156 
157  do_ep_renorm = .False.
158  ntemp = 1
159  if (BSp%do_ep_renorm .and. PRESENT(Epren)) then
160    do_ep_renorm = .True.
161    ntemp = Epren%ntemp
162  end if
163 
164  !
165  ! =====================================================
166  ! === Calculate fcv(k)=<c k s|e^{-iqr}|v k s> in BZ ===
167  ! =====================================================
168  lomo_min=Bsp%lomo_min; max_band=Bsp%nbnds
169  ABI_MALLOC(opt_cvk,(lomo_min:max_band,lomo_min:max_band,BSp%nkbz,nsppol,BSp%nq))
170 
171  do iq=1,BSp%nq
172    call calc_optical_mels(Wfd,Kmesh,KS_BSt,Cryst,Psps,Pawtab,Hur,BSp%inclvkb,Bsp%lomo_spin,lomo_min,max_band,&
173 &                         BSp%nkbz,BSp%q(:,iq),opt_cvk(:,:,:,:,iq))
174  end do
175  !
176  ! ============================
177  ! ==== Make EPS EXCITONIC ====
178  ! ============================
179  if (my_rank==master) then ! Only master works.
180 
181    ABI_MALLOC(eps_exc,(BSp%nomega,BSp%nq))
182    ABI_MALLOC(dos_exc,(BSp%nomega))
183 
184    ABI_MALLOC(eps_rpanlf,(BSp%nomega,BSp%nq))
185    ABI_MALLOC(dos_ks,(BSp%nomega))
186 
187    ABI_MALLOC(eps_gwnlf ,(BSp%nomega,BSp%nq))
188    ABI_MALLOC(dos_gw,(BSp%nomega))
189 
190    do itemp = 1, ntemp
191 
192      call int2char4(itemp,ts)
193 
194      if(do_ep_renorm) then
195        prefix = TRIM("_T") // ts
196      else
197        prefix = ""
198      end if
199 
200      ost_fname = strcat(BS_files%out_basename,prefix,"_EXC_OST")
201 
202      !TODO for RPA
203      call ebands_copy(KS_BST, EPBSt)
204      call ebands_copy(QP_BST, EP_QPBSt)
205 
206      if (BS_files%in_eig /= BSE_NOFILE) then
207        filbseig = strcat(BS_files%in_eig,prefix)
208      else
209        filbseig = strcat(BS_files%out_eig,prefix)
210      end if
211 
212      if(do_ep_renorm) then
213        ! No scissor with KSBST
214        call renorm_bst(Epren, EPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.TRUE.)
215 
216        call renorm_bst(Epren, EP_QPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.FALSE.)
217      end if
218 
219 
220      if (BSp%use_coupling==0) then
221        call exc_eps_resonant(BSp,filbseig,ost_fname,lomo_min,max_band,BSp%nkbz,nsppol,opt_cvk,&
222 &        Cryst%ucvol,BSp%nomega,BSp%omega,eps_exc,dos_exc,elph_lifetime=do_ep_renorm)
223      else
224        call exc_eps_coupling(Bsp,BS_files,lomo_min,max_band,BSp%nkbz,nsppol,opt_cvk,&
225 &        Cryst%ucvol,BSp%nomega,BSp%omega,eps_exc,dos_exc)
226      end if
227      !
228      ! =======================================================
229      ! === Make EPS RPA and GW without local-field effects ===
230      ! =======================================================
231      call wrtout(std_out," Calculating RPA NLF and QP NLF epsilon","COLL")
232 
233      call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,Bsp%lomo_min,BSp%homo_spin,Kmesh,EPBSt,BSp%nq,nsppol,opt_cvk,&
234 &      Cryst%ucvol,BSp%broad,BSp%nomega,BSp%omega,eps_rpanlf,dos_ks)
235 
236      call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,Bsp%lomo_min,BSp%homo_spin,Kmesh,EP_QPBSt,BSp%nq,nsppol,opt_cvk,&
237 &      Cryst%ucvol,Bsp%broad,BSp%nomega,BSp%omega,eps_gwnlf,dos_gw)
238      !
239      ! =========================
240      ! === Write out Epsilon ===
241      ! =========================
242      !this is just for the automatic tests, It will be removed when fldiff
243      !will be able to compare two optical spectral
244      write(ab_out,*)" "
245      write(ab_out,*)"Macroscopic dielectric function:"
246      write(ab_out,*)"omega [eV] <KS_RPA_nlf>  <GW_RPA_nlf>  <BSE> "
247      do io=1,MIN(10,BSp%nomega)
248        omegaev = REAL(BSp%omega(io))*Ha_eV
249        ks_avg  = SUM( eps_rpanlf(io,:)) / Bsp%nq
250        gw_avg  = SUM( eps_gwnlf (io,:)) / Bsp%nq
251        exc_avg = SUM( eps_exc   (io,:)) / Bsp%nq
252        write(ab_out,'(7f9.4)')omegaev,ks_avg,gw_avg,exc_avg
253      end do
254      write(ab_out,*)" "
255 
256      !
257      ! Master node writes final results on file.
258      call exc_write_data(BSp,BS_files,"RPA_NLF_MDF",eps_rpanlf,prefix=prefix,dos=dos_ks)
259 
260      call exc_write_data(BSp,BS_files,"GW_NLF_MDF",eps_gwnlf,prefix=prefix,dos=dos_gw)
261 
262      call exc_write_data(BSp,BS_files,"EXC_MDF",eps_exc,prefix=prefix,dos=dos_exc)
263 
264      call wrtout(std_out," Checking Kramers Kronig on Excitonic Macroscopic Epsilon","COLL")
265      call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_exc(:,1))
266 
267      call wrtout(std_out," Checking Kramers Kronig on RPA NLF Macroscopic Epsilon","COLL")
268      call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_rpanlf(:,1))
269 
270      call wrtout(std_out," Checking Kramers Kronig on GW NLF Macroscopic Epsilon","COLL")
271      call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_gwnlf(:,1))
272 
273      call wrtout(std_out," Checking f-sum rule on Excitonic Macroscopic Epsilon","COLL")
274 
275      if (BSp%exchange_term>0) then
276        MSG_COMMENT(' f-sum rule should be checked without LF')
277      end if
278      call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_exc(:,1)),drude_plsmf)
279 
280      call wrtout(std_out," Checking f-sum rule on RPA NLF Macroscopic Epsilon","COLL")
281      call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_rpanlf(:,1)),drude_plsmf)
282 
283      call wrtout(std_out," Checking f-sum rule on GW NLF Macroscopic Epsilon","COLL")
284      call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_gwnlf(:,1)),drude_plsmf)
285 
286 #ifdef HAVE_NETCDF
287      path = strcat(BS_files%out_basename, strcat(prefix,"_MDF.nc"))
288      NCF_CHECK_MSG(nctk_open_create(ncid, path, xmpi_comm_self), sjoin("Creating MDF file:", path))
289      ! Write structure
290      NCF_CHECK(crystal_ncwrite(Cryst, ncid))
291      ! Write QP energies
292      NCF_CHECK(ebands_ncwrite(QP_BSt, ncid))
293      ! Write dielectric functions.
294      call mdfs_ncwrite(ncid, Bsp, eps_exc,eps_rpanlf,eps_gwnlf)
295      NCF_CHECK(nf90_close(ncid))
296 #else
297      ABI_UNUSED(ncid)
298 #endif
299 
300      !TODO
301      call ebands_free(EPBSt)
302      call ebands_free(EP_QPBSt)
303 
304    end do
305 
306    ABI_FREE(eps_rpanlf)
307    ABI_FREE(eps_gwnlf)
308    ABI_FREE(eps_exc)
309    ABI_FREE(dos_exc)
310    ABI_FREE(dos_ks)
311    ABI_FREE(dos_gw)
312  end if ! my_rank==master
313 
314  ABI_FREE(opt_cvk)
315 
316  call xmpi_barrier(comm)
317 
318 end subroutine build_spectra

m_exc_spectra/check_fsumrule [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  check_fsumrule

FUNCTION

   check f-sum rule
   \int_0^\infty d\omega \omega Im \epsilon_GG'(q,\omega) =
   = \frac{1}{2} \pi \omega_p^2  \frac{\rho(G-G')}{\rho(0)}
   versor(q+G) \dot versor(q+G')
   for q = G = G' = 0, it reads:
   \int_0^\infty d\omega \omega Im \epsilon_00(q=0,\omega) =
   = \pi \omega_p^2 / 2
   calculate only the second one
   calculate the integral to evaluate an omega_plasma^eff to compare with omega_plasma

INPUTS

  n=Number of frequencies.
  o(n)=Frequency mesh.
  e2(n)=imaginary part of epsilon_00
  omegaplasma=Drude plasma frequency.

OUTPUT

  Only checking.

PARENTS

      m_exc_spectra

CHILDREN

      wrtout

SOURCE

1551 subroutine check_fsumrule(n,o,e2,omegaplasma)
1552 
1553 
1554 !This section has been created automatically by the script Abilint (TD).
1555 !Do not modify the following lines by hand.
1556 #undef ABI_FUNC
1557 #define ABI_FUNC 'check_fsumrule'
1558 !End of the abilint section
1559 
1560  implicit none
1561 
1562 !Arguments ------------------------------------
1563 !scalars
1564  integer,intent(in) :: n
1565  real(dp),intent(in) :: omegaplasma
1566 !arrays
1567  real(dp),intent(in) :: o(n),e2(n)
1568 
1569 !Local variables ------------------------------
1570 !scalars
1571  integer :: ii,ip
1572  real(dp) :: omegap,domega,integral,omegaplasmaeff,fsumrule
1573  character(len=500) :: msg
1574 !arrays
1575  complex(dpc) :: intg(n)
1576 
1577 !************************************************************************
1578 
1579 ! calculate domega step and verify
1580  domega = (o(n) - o(1)) / (n-1)
1581 
1582  do ii=2,n
1583    if (domega-(o(ii)-o(ii-1)) > tol3) then
1584      MSG_WARNING("Frequency mesh not linear. Returning")
1585      return
1586    end if
1587  end do
1588 
1589  if (o(1) > 0.1/Ha_eV) then
1590    MSG_WARNING("First frequency is not zero. Returning")
1591    return
1592  end if
1593 
1594  if (e2(n) > 0.1) then
1595    write(msg,'(a,f12.6,3a,f12.6,2a)')&
1596 &   ' Im epsilon for omega= ',o(n)*Ha_eV,' eV ',ch10,&
1597 &   ' is not yet zero, epsilon_2= ',e2(n),ch10,&
1598 &   ' f-sum rule test could give wrong results.'
1599    MSG_WARNING(msg)
1600  end if
1601 
1602 ! integrate to obtain f-sum rule
1603  integral=zero
1604  do ip=1,n
1605    omegap=o(ip)
1606    integral = integral + omegap * e2(ip)
1607  end do
1608  integral = domega * integral
1609 
1610 !integrate with simpson to obtain f-sum rule
1611  do ip = 1, n
1612    omegap = o(ip)
1613    intg(ip) = omegap * e2(ip)
1614  end do
1615 
1616  integral = real(simpson_cplx(n,domega,intg))
1617  if(integral < 0) then
1618    MSG_ERROR("The integral of the imaginary of dielectric function is negative !!!")
1619  else
1620    omegaplasmaeff = sqrt(integral*two/pi)
1621  end if
1622 
1623  fsumrule = abs((omegaplasmaeff - omegaplasma)) / omegaplasma
1624 
1625 ! write data
1626  write(msg,'(3(a,f6.2,2a))')&
1627 &  " omega_plasma     = ",omegaplasma*Ha_eV,   " [eV]",ch10,&
1628 &  " omega_plasma^eff = ",omegaplasmaeff*Ha_eV," [eV]",ch10,&
1629 &  " the f-sum rule is verified within ",fsumrule*100,"%",ch10
1630  call wrtout(std_out,msg,"COLL")
1631 
1632 end subroutine check_fsumrule

m_exc_spectra/check_kramerskronig [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  check_kramerskronig

FUNCTION

   check Kramers Kronig
   \int_0^\infty d\omega' frac{\omega'}{\omega'^2 - \omega^2}
   Im \epsilon(\omega') = Re \epsilon(\omega)

INPUTS

  n=Number of frequency points.
  eps(n)=Dielectric function.
  o(n)=Frequency mesh.

OUTPUT

  Only checking.

PARENTS

      m_exc_spectra

CHILDREN

      wrtout

SOURCE

1410 subroutine check_kramerskronig(n,o,eps)
1411 
1412 
1413 !This section has been created automatically by the script Abilint (TD).
1414 !Do not modify the following lines by hand.
1415 #undef ABI_FUNC
1416 #define ABI_FUNC 'check_kramerskronig'
1417 !End of the abilint section
1418 
1419  implicit none
1420 
1421 !Arguments ------------------------------------
1422 !scalars
1423  integer,intent(in) :: n
1424 !arrays
1425  complex(dpc),intent(in) :: eps(n)
1426  real(dp),intent(in) :: o(n)
1427 
1428 !Local variables ------------------------------
1429 !scalars
1430  integer :: ii,ip
1431  real(dp) :: omega,omegap,domega,kk,kkrms,eav
1432  complex(dpc) c
1433  character(len=500) :: msg
1434 !arrays
1435  real(dp) :: e1kk(n)
1436  complex(dpc) :: intg(n)
1437 
1438 !************************************************************************
1439 ! init jmb
1440  e1kk=zero
1441  intg=(zero,zero)
1442 
1443 ! calculate domega step and verify all
1444  domega = (o(n) - o(1)) / (n-1)
1445 
1446  do ii=2,n
1447   if (domega-(o(ii)-o(ii-1)) > tol3) then
1448     MSG_WARNING("Frequency mesh not linear. Returning")
1449     return
1450   end if
1451  end do
1452 
1453  if(o(1) > 0.1/Ha_eV) then
1454    MSG_WARNING("First frequency is not zero. Returning")
1455    return
1456  end if
1457 
1458  if (aimag(eps(n)) > 0.1) then
1459    write(msg,'(a,f12.6,3a,f12.6,2a)')&
1460 &   ' Im epsilon for omega= ',o(n)*Ha_eV,'eV',ch10,&
1461 &   ' is not yet zero, epsilon_2= ',aimag(eps(n)),ch10,&
1462 &   ' Kramers Kronig test could give wrong results. '
1463    MSG_WARNING(msg)
1464  end if
1465 
1466 ! Fill array for kramers kronig.
1467  do ii=1,n
1468    omega=o(ii)
1469    c = (0.0,0.0)
1470    do ip=1,n
1471      if(ip == ii) cycle
1472      omegap = o(ip)
1473      c = c + omegap / (omegap**2-omega**2) * aimag(eps(ip))
1474    end do
1475    e1kk(ii) = one + two/pi * domega*real(c)
1476  end do
1477 
1478 !perform kramers kronig with simpson integration
1479  do ii=1,n
1480    omega=o(ii)
1481    do ip=1,n
1482      if (ip==ii) cycle
1483      omegap = o(ip)
1484      intg(ip) = omegap / (omegap**2 - omega**2) * aimag(eps(ip))
1485    end do
1486    c = simpson_cplx(n,domega,intg)
1487    e1kk(ii) = one + two/pi * real(c)
1488  end do
1489 
1490 !verify kramers kronig
1491  eav=zero; kk=zero; kkrms=zero
1492  do ii=1,n
1493    kk = kk + abs(real(eps(ii)) - e1kk(ii))
1494    kkrms = kkrms +(real(eps(ii)) - e1kk(ii))*(real(eps(ii)) - e1kk(ii))
1495    eav = eav + abs(real(eps(ii)))
1496  end do
1497 
1498  eav = eav/n
1499  kk = (kk/n)/eav
1500  kkrms = (kkrms/n) / (eav*eav)
1501 
1502  kk = abs(real(eps(1)) - e1kk(1)) / real(eps(1))
1503 
1504 ! write data
1505  write(msg,'(a,f7.2,a)')" The Kramers-Kronig is verified within ",100*kk,"%"
1506  call wrtout(std_out,msg,"COLL")
1507 
1508 ! write(std_out,'("# Kramers Kronig calculation of epsilon1")')
1509 ! write(std_out,'("# omega   epsilon1  epsilon1kk")')
1510 ! do ii=1,n
1511 !   write(std_out,'(f7.3,2e15.7)') o(ii)*Ha_eV, real(eps(ii)), e1kk(ii)
1512 ! end do
1513 
1514 end subroutine check_kramerskronig

m_exc_spectra/exc_eps_coupling [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  exc_eps_coupling

FUNCTION

  Make epsilon EXCITONIC with full COUPLING.

INPUTS

 Bsp
 nkbz=Number of points in the BZ
 lomo_min,max_band
 nomega=Number of frequencies
 omega(nomega)=frequency mesh.
 nsppol=Number of independent spin polarizations.
 ucvol=Unit cell volume.
 BS_files<excfiles>File names used in the Bethe-Salpeter code.
 opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol)=Matrix elements <b k|e^{-iqr}|b" k> for a given q in the full BZ.

OUTPUT

  eps_exc(nomega)=Macroscopic dielectric function with excitonic effects calculated including the COUPLING.
  dos_exc(nomega)=The DOS of the excitonic Hamiltonian

PARENTS

      m_exc_spectra

CHILDREN

      c_f_pointer

SOURCE

 931 subroutine exc_eps_coupling(Bsp,BS_files,lomo_min,max_band,nkbz,nsppol,opt_cvk,ucvol,nomega,omega,eps_exc,dos_exc)
 932 
 933 
 934 !This section has been created automatically by the script Abilint (TD).
 935 !Do not modify the following lines by hand.
 936 #undef ABI_FUNC
 937 #define ABI_FUNC 'exc_eps_coupling'
 938 !End of the abilint section
 939 
 940  implicit none
 941 
 942 !Arguments ------------------------------------
 943 !scalars
 944  integer,intent(in) :: lomo_min,max_band,nkbz,nomega,nsppol
 945  real(dp),intent(in) :: ucvol
 946  type(excfiles),intent(in) :: BS_files
 947  type(excparam),intent(in) :: BSp
 948 !arrays
 949  real(dp),intent(out) :: dos_exc(nomega)
 950  complex(dpc),intent(in) :: opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol,BSp%nq),omega(nomega)
 951  complex(dpc),intent(out) :: eps_exc(nomega,BSp%nq)
 952 
 953 !Local variables ------------------------------
 954 !scalars
 955  integer :: mi,it,ii,ib_v,ib_c,ik_bz,exc_size_read,nstates_read,eig_unt !,fform
 956  integer :: exc_size,iq,spin,tr_idx,tar_idx,nstates,iw,ll,ierr
 957  real(dp) :: fact,arg
 958  complex(dpc) :: eps,fam,famp
 959  character(len=500) :: msg,errmsg
 960  character(len=fnlen) :: filbseig
 961  logical :: do_lifetime
 962  !type(Hdr_type) :: tmp_Hdr
 963 !arrays
 964  complex(dpc),allocatable :: Ami(:),exc_ene(:),Sm1mi(:)
 965  complex(dpc),allocatable :: msfap(:,:),fa(:,:),fap(:,:)
 966 
 967 !************************************************************************
 968 
 969  call wrtout(std_out," Calculating absorption strength with full coupling","COLL")
 970 
 971  if (nsppol==2) then
 972    MSG_WARNING("nsppol==2 is still under development")
 973  end if
 974 
 975  ! Rank of the entire excitonic Hamiltonian including the coupling block.
 976  exc_size = 2*SUM(BSp%nreh); if (nsppol==2) exc_size = 2*(SUM(BSp%nreh) + BSp%nreh(2))
 977  nstates  = BSp%nstates
 978 
 979  ! TODO: four_pi comes from the bare Coulomb term hence the
 980  ! present implementation is not compatible with the cutoff technique.
 981  ! factor two is due to the occupation factors.
 982  fact=four_pi/(ucvol*nkbz); if (nsppol==1) fact=two*fact
 983 
 984  if (BS_files%in_eig /= BSE_NOFILE) then
 985    filbseig = BS_files%in_eig
 986  else
 987    filbseig = BS_files%out_eig
 988  end if
 989 
 990  call wrtout(std_out," Reading excitonic eigenstates from file: "//trim(filbseig),"COLL")
 991  if (open_file(filbseig,msg,newunit=eig_unt,form="unformatted", status="old", action="read") /= 0) then
 992    MSG_ERROR(msg)
 993  end if
 994 
 995  read(eig_unt, err=10, iomsg=errmsg) do_lifetime
 996 
 997  if (do_lifetime) then
 998    ABI_CHECK(.not. do_lifetime, "Finite lifetime with coupling is not supported yet !")
 999  end if
1000 
1001  read(eig_unt, err=10, iomsg=errmsg) exc_size_read, nstates_read
1002  ABI_CHECK(exc_size_read==exc_size,"wrong file")
1003  ABI_CHECK(nstates_read==nstates,"Partial diago not supported yet")
1004  !
1005  ! Read eigenvalues
1006  ABI_MALLOC(exc_ene,(nstates))
1007  read(eig_unt, err=10, iomsg=errmsg) exc_ene(:)
1008 
1009  ABI_MALLOC(fa,(nstates,BSp%nq))
1010  ABI_MALLOC(fap,(nstates,BSp%nq))
1011  ABI_STAT_MALLOC(Ami,(exc_size), ierr)
1012  ABI_CHECK(ierr==0, " out-of-memory Ami")
1013 
1014  do mi=1,nstates ! Loop on excitonic eigenvalues mi
1015    read(eig_unt, err=10, iomsg=errmsg) Ami(:)
1016 
1017    do iq=1,BSp%nq
1018      fam  = czero
1019      famp = czero
1020      do spin=1,nsppol
1021        do it=1,BSp%nreh(spin) ! Loop over transition t = (k,v,c)
1022          ik_bz = Bsp%Trans(it,spin)%k
1023          ib_v  = Bsp%Trans(it,spin)%v
1024          ib_c  = Bsp%Trans(it,spin)%c
1025          tr_idx  = it + (spin-1)*Bsp%nreh(1)
1026          if (nsppol==1) then
1027            tar_idx = it + Bsp%nreh(1)
1028          else
1029            if (spin==1) tar_idx = it + SUM(Bsp%nreh)
1030            if (spin==2) tar_idx = it + 2*Bsp%nreh(1)+Bsp%nreh(2)
1031          end if
1032 
1033          fam = fam + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * Ami(tr_idx) &
1034 &                  + CONJG(opt_cvk(ib_v,ib_c,ik_bz,spin,iq)) * Ami(tar_idx)
1035 
1036          famp = famp - opt_cvk(ib_c,ib_v,ik_bz,spin,iq) * CONJG(Ami(tr_idx)) &
1037 &                    + opt_cvk(ib_v,ib_c,ik_bz,spin,iq) * CONJG(Ami(tar_idx))
1038        end do
1039      end do
1040      ! Save results.
1041      fa (mi,iq) = fam
1042      fap(mi,iq) = famp
1043    end do
1044  end do ! mi
1045 
1046  ABI_FREE(Ami)
1047 
1048  ! Read O{-1} and sum over the eigenstates.
1049  ABI_MALLOC(msfap,(nstates,BSp%nq))
1050  ABI_MALLOC(Sm1mi,(nstates))
1051 
1052  do mi=1,nstates
1053    read(eig_unt, err=10, iomsg=errmsg) Sm1mi
1054    Sm1mi = DCONJG(Sm1mi) ! This gives the row since O^{-1} is Hermitian.
1055    do iq=1,BSp%nq
1056      msfap(mi,iq) = xdotu(exc_size,Sm1mi,1,fap(:,iq),1)
1057    end do
1058  end do
1059 
1060  ABI_FREE(Sm1mi)
1061 
1062  close(eig_unt, err=10, iomsg=errmsg)
1063  !
1064  ! === Calculate excitonic epsilon with coupling ===
1065  do iq=1,BSp%nq
1066    !
1067    do ii=1,nomega
1068      eps = czero
1069      do mi=1,nstates ! sum over all exciton eigenstates
1070        eps = eps - fa(mi,iq) * msfap(mi,iq) / (exc_ene(mi) - omega(ii))
1071      end do
1072      eps_exc(ii,iq) = one + fact * eps
1073    end do
1074    !
1075  end do
1076 
1077  ABI_FREE(fa)
1078  ABI_FREE(msfap)
1079  ABI_FREE(fap)
1080  !
1081  ! The excitonic DOS.
1082  dos_exc=zero
1083  do ll=1,nstates ! Sum over the calculate excitonic eigenstates.
1084    do iw=1,nomega
1085      arg = DBLE(omega(iw) - exc_ene(ll))
1086      dos_exc(iw) = dos_exc(iw) + dirac_delta(arg,Bsp%broad)
1087    end do
1088  end do
1089 
1090  ABI_FREE(exc_ene)
1091 
1092  return
1093 
1094 10 continue
1095  MSG_ERROR(errmsg)
1096 
1097 end subroutine exc_eps_coupling

m_exc_spectra/exc_eps_resonant [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  exc_eps_resonant

FUNCTION

  This routine builds the macroscopic dielectric function with excitonic effects.

INPUTS

 Bsp
 lomo_min,max_band
 nkbz=Number of points in the BZ
 nsppol=Number of independent polarizations.
 nomega=Number of frequencies
 omega(nomega)=frequency mesh (complex shift is already included)
 ucvol=Volume of the unit cell.
 opt_cvk(lomo_min:max_band,mib:max_band,nkbz,nsppol,Bsp%nq)=Matrix elements <b k|e^{-iqr}|b" k> for a given q in the full BZ.

OUTPUT

  eps_exc(nomega,Bsp%nq)=Macroscopic dielectric function with excitonic effects.
  dos_exc(nomega)=The DOS of the excitonic Hamiltonian

PARENTS

      m_exc_spectra

CHILDREN

      c_f_pointer

SOURCE

671 subroutine exc_eps_resonant(Bsp,filbseig,ost_fname,lomo_min,max_band,nkbz,nsppol,opt_cvk,&
672 &    ucvol,nomega,omega,eps_exc,dos_exc,elph_lifetime)
673 
674 
675 !This section has been created automatically by the script Abilint (TD).
676 !Do not modify the following lines by hand.
677 #undef ABI_FUNC
678 #define ABI_FUNC 'exc_eps_resonant'
679 !End of the abilint section
680 
681  implicit none
682 
683 !Arguments ------------------------------------
684 !scalars
685  integer,intent(in) :: lomo_min,max_band,nkbz,nomega,nsppol
686  real(dp),intent(in) :: ucvol
687  type(excparam),intent(in) :: BSp
688  character(len=fnlen),intent(in) :: filbseig,ost_fname
689  logical,optional,intent(in) :: elph_lifetime
690 !arrays
691  real(dp),intent(out) :: dos_exc(nomega)
692  complex(dpc),intent(in) :: opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol,BSp%nq),omega(nomega)
693  complex(dpc),intent(out) :: eps_exc(nomega,BSp%nq)
694 
695 !Local variables ------------------------------
696 !scalars
697  integer :: ll,it,iw,ib_v,ib_c,ik_bz,neig_read,eig_unt,exc_size,iq !fform,
698  integer :: spin,spad,hsize_read,nstates,ost_unt
699  logical :: do_ep_lifetime, file_do_lifetime
700  real(dp) :: fact,arg
701  complex(dpc) :: dotprod
702  character(len=500) :: msg,frm,errmsg
703  !type(Hdr_type) :: tmp_Hdr
704 !arrays
705  real(dp),allocatable :: exc_ene(:)
706  complex(dpc) :: ctemp(BSp%nq),dtemp(BSp%nq)
707  complex(dpc),allocatable :: ostrength(:,:),exc_ene_cplx(:),exc_state(:),exc_state2(:)
708 
709 !************************************************************************
710 
711  call wrtout(std_out," Calculating excitonic epsilon with antiresonant","COLL")
712 
713  if (nsppol==2) then
714    MSG_WARNING("nsppol==2 still under development")
715  end if
716 
717  exc_size = SUM(BSp%nreh)
718  nstates  = BSp%nstates
719 
720  do_ep_lifetime = .FALSE.
721  if(PRESENT(elph_lifetime)) then
722    do_ep_lifetime = elph_lifetime
723  end if
724 
725  if (ANY(Bsp%nreh/=Bsp%nreh(1))) then
726    write(msg,'(a,2(i0,1x))')"BSE does not support different number of transitions for the two spin channels. nreh: ",Bsp%nreh
727    MSG_WARNING(msg)
728  end if
729  !
730  ! TODO:
731  ! four_pi comes from the bare Coulomb term hence the
732  ! present implementation is not compatible with the cutoff technique.
733  fact=four_pi/(ucvol*nkbz); if (nsppol==1) fact=two*fact ! two to account for the occupation numbers.
734 
735  call wrtout(std_out," Reading excitonic eigenstates from file: "//TRIM(filbseig),"COLL")
736  if (open_file(filbseig,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then
737    MSG_ERROR(msg)
738  end if
739 
740  read(eig_unt, err=10, iomsg=errmsg) file_do_lifetime
741 
742  if(do_ep_lifetime .and. .not. file_do_lifetime) then
743   MSG_ERROR("Cannot do lifetime as the data is not present in the file !")
744  end if
745 
746  read(eig_unt, err=10, iomsg=errmsg) hsize_read,neig_read
747 
748  if (hsize_read /= exc_size) then
749    write(msg,'(2(a,i0))')" Wrong size of the Hamiltonian: read: ",hsize_read," expected= ",exc_size
750    MSG_ERROR(msg)
751  end if
752 
753  if (neig_read /= nstates) then
754    write(msg,'(2(a,i0))')" Wrong number of eigenstates: read: ",neig_read," expected= ",nstates
755    MSG_ERROR(msg)
756  end if
757  !
758  ! Read eigenvalues, ignore possibly small imaginary part.
759  ABI_MALLOC(exc_ene_cplx,(neig_read))
760  read(eig_unt, err=10, iomsg=errmsg) exc_ene_cplx
761 
762  ABI_MALLOC(exc_ene,(neig_read))
763  exc_ene = DBLE(exc_ene_cplx)
764  !ABI_FREE(exc_ene_cplx)
765  !
766  ! Calculate oscillator strength.
767  ABI_MALLOC(exc_state,(exc_size))
768  ABI_MALLOC(ostrength,(neig_read,BSp%nq))
769 
770  if (do_ep_lifetime) then
771    ABI_MALLOC(exc_state2,(exc_size))
772    do ll=1,neig_read ! Loop over excitonic eigenstates reported on file.
773      read(eig_unt, err=10, iomsg=errmsg) exc_state(:) ! Righteigenvector
774      read(eig_unt, err=10, iomsg=errmsg) exc_state2(:) ! Lefteigenvector
775 
776      ! Here assuming that eigenvectors are such as Xl_i' Xr_j = delta_ij
777      ! Otherwise, I need to invert the overlap matrix !
778 
779      ! Rescale the vectors so that they are "normalized" with respect to the other one !
780      dotprod = xdotc(exc_size,exc_state2(:),1,exc_state(:),1)
781      exc_state2(:) = exc_state2(:)/CONJG(dotprod)
782 
783      ctemp(:) = czero
784      dtemp(:) = czero
785      do spin=1,nsppol
786        spad=(spin-1)*BSp%nreh(1) ! Loop over spin channels.
787        do it=1,BSp%nreh(spin)    ! Loop over resonant transition t = (k,v,c,s)
788          ik_bz = Bsp%Trans(it,spin)%k
789          ib_v  = Bsp%Trans(it,spin)%v
790          ib_c  = Bsp%Trans(it,spin)%c
791          do iq=1,BSp%nq
792            ctemp(iq) = ctemp(iq) + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * exc_state(it+spad)
793            dtemp(iq) = dtemp(iq) + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * exc_state2(it+spad)
794          end do
795        end do ! it
796      end do
797      ostrength(ll,:) = ctemp(:)*CONJG(dtemp(:))
798    end do ! ll
799    ABI_FREE(exc_state2)
800  else
801    do ll=1,neig_read ! Loop over excitonic eigenstates reported on file.
802      read(eig_unt, err=10, iomsg=errmsg) exc_state(:)
803      if(file_do_lifetime) read(eig_unt, err=10, iomsg=errmsg)
804 
805      ctemp(:) = czero
806      do spin=1,nsppol
807        spad=(spin-1)*BSp%nreh(1) ! Loop over spin channels.
808        do it=1,BSp%nreh(spin)    ! Loop over resonant transition t = (k,v,c,s)
809          ik_bz = Bsp%Trans(it,spin)%k
810          ib_v  = Bsp%Trans(it,spin)%v
811          ib_c  = Bsp%Trans(it,spin)%c
812          do iq=1,BSp%nq
813            ctemp(iq) = ctemp(iq) + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * exc_state(it+spad)
814          end do
815        end do ! it
816      end do
817      ostrength(ll,:) = ctemp(:)*CONJG(ctemp(:))
818    end do ! ll
819  end if
820 
821 
822  close(eig_unt, err=10, iomsg=errmsg)
823  ABI_FREE(exc_state)
824 
825  if(do_ep_lifetime) then
826    eps_exc = one
827    do ll=1,neig_read ! Sum over all excitonic eigenstates read from file.
828      do iq=1,BSp%nq
829         do iw=1,nomega
830           eps_exc(iw,iq) = eps_exc(iw,iq) +  &
831   &         fact * ostrength(ll,iq) * (one/(exc_ene_cplx(ll) - omega(iw)) - one/(-DCONJG(exc_ene_cplx(ll)) - omega(iw)))
832         end do
833      end do !ll
834    end do !iw
835  else
836    eps_exc = one
837    do ll=1,neig_read ! Sum over all excitonic eigenstates read from file.
838      do iq=1,BSp%nq
839         do iw=1,nomega
840           eps_exc(iw,iq) = eps_exc(iw,iq) +  &
841   &         fact * ostrength(ll,iq) * (one/(exc_ene(ll) - omega(iw)) - one/(-exc_ene(ll) - omega(iw)))
842         end do
843      end do !ll
844    end do !iw
845  end if
846 
847  !
848  ! The excitonic DOS.
849  dos_exc=zero
850  do ll=1,neig_read ! Sum over the calculate excitonic eigenstates.
851    do iw=1,nomega
852      arg = ( DBLE(omega(iw)) - exc_ene(ll))
853      if(do_ep_lifetime) then
854        dos_exc(iw) = dos_exc(iw) + dirac_delta(arg,AIMAG(exc_ene_cplx(ll)))
855      else
856        dos_exc(iw) = dos_exc(iw) + dirac_delta(arg,Bsp%broad)
857      end if
858    end do
859  end do
860  !
861  ! Write the oscillator strengths to file.
862  if (open_file(ost_fname,msg,newunit=ost_unt,form="formatted",action="write") /= 0) then
863    MSG_ERROR(msg)
864  end if
865 
866  write(ost_unt,'("# Oscillator strengths of the excitonic states for the different q-polarizations.")')
867  !
868  ! Write the list of q-points.
869  write(ost_unt,*)"# List of q-points for the optical limit"
870  do iq=1,BSp%nq
871    write(ost_unt,'(a,3(f9.6,","),a)')'# q = ',BSp%q(:,iq),' [Reduced coords] '
872  end do
873 
874  write(ost_unt,*)"# E_lambda [eV]     ostrength(q=1) ostrength(q=2) .... "
875  write(frm,*)'(f8.4,',BSp%nq,'es12.4)'
876  do ll=1,neig_read
877    write(ost_unt,frm)exc_ene(ll)*Ha_eV,(ostrength(ll,iq), iq=1,BSp%nq)
878  end do
879 
880  close(ost_unt)
881 
882  ABI_FREE(ostrength)
883  ABI_FREE(exc_ene)
884  ABI_FREE(exc_ene_cplx)
885 
886  !call exc_amplitude(Bsp,filbseig,1,(/(ll,ll=1,10)/),"TEST_AMPLITUDE")
887  !call exc_amplitude(Bsp,filbseig,1,(/30/),"TEST_AMPLITUDE")
888 
889  return
890 
891  ! Handler IO-error
892 10 continue
893  MSG_ERROR(errmsg)
894 
895 end subroutine exc_eps_resonant

m_exc_spectra/exc_eps_rpa [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  exc_eps_rpa

FUNCTION

  Build epsilon within RPA and GW.

INPUTS

 nkbz=Number of points in the BZ
 nbnds=Number of bands
 lomo_spin(nsppol)
 lomo_min=Lowest occupied state
 homo=Number of occupied states.
 homo_spin(nsppol)
 nsppol=Number of independent spin polarizations.
 nomega=Number of frequencies
 omega(nomega)=Frequency mesh.
 ucvol=Unit cell volume.
 broad=Broadening used for the DOS.
 opt_cvk(nbnds,nbnds,nkbz)=Matrix elements <b k|e^{-iqr}|b" k> for a given q in the full BZ.

OUTPUT

  eps_rpa(nomega)=RPA spectrum without local-field effects.
  dos(nomega)=The DOS.

PARENTS

      m_exc_spectra,m_haydock

CHILDREN

      c_f_pointer

SOURCE

532 subroutine exc_eps_rpa(nbnds,lomo_spin,lomo_min,homo_spin,Kmesh,Bst,nq,nsppol,opt_cvk,ucvol,broad,nomega,omega,eps_rpa,dos)
533 
534 
535 !This section has been created automatically by the script Abilint (TD).
536 !Do not modify the following lines by hand.
537 #undef ABI_FUNC
538 #define ABI_FUNC 'exc_eps_rpa'
539 !End of the abilint section
540 
541  implicit none
542 
543 !Arguments ------------------------------------
544 !scalars
545  integer,intent(in) :: nbnds,lomo_min,nsppol,nomega,nq
546  real(dp),intent(in) :: ucvol,broad
547  type(kmesh_t),intent(in) :: Kmesh
548  type(ebands_t),intent(in) :: BSt
549 !arrays
550  integer,intent(in) :: lomo_spin(nsppol),homo_spin(nsppol)
551  real(dp),intent(out) :: dos(nomega)
552  complex(dpc),intent(in) :: omega(nomega)
553  complex(dpc),intent(in) :: opt_cvk(lomo_min:nbnds,lomo_min:nbnds,Kmesh%nbz,nsppol,nq)
554  complex(dpc),intent(out) :: eps_rpa(nomega,nq)
555 
556 !Local variables ------------------------------
557 !scalars
558  integer :: iw,ib_v,ib_c,ik_bz,ik_ibz,spin,iq
559  real(dp) :: fact,arg,ediff
560  real(dp) :: lifetime
561  complex(dpc) :: ctemp
562  logical :: do_lifetime
563 
564 !************************************************************************
565 
566  ! TODO: four_pi comes from the bare Coulomb term hence the
567  ! present implementation is not compatible with the cutoff technique.
568  fact=four_pi/(ucvol*Kmesh%nbz)
569  if (nsppol==1) fact=two*fact ! two accounts for the occupation factors.
570 
571  eps_rpa=czero; dos=zero
572 
573  do_lifetime = .FALSE.
574  do_lifetime = allocated(BSt%lifetime)
575 
576  !write(std_out,*)nsppol,Kmesh%nbz,lomo_min,homo,nbnds
577  !
578  ! Sum over all QP transitions.
579  do spin=1,nsppol
580    do ik_bz=1,Kmesh%nbz
581      ik_ibz = Kmesh%tab(ik_bz)
582      do ib_v=lomo_spin(spin),homo_spin(spin)
583        do ib_c=homo_spin(spin)+1,nbnds
584          !
585          ! TODO here energies are always assumed to be real.
586          ediff = BSt%eig(ib_c,ik_ibz,spin) - BSt%eig(ib_v,ik_ibz,spin)
587 
588          !
589          if(do_lifetime) then
590            lifetime = BSt%lifetime(ib_c,ik_ibz,spin) + BSt%lifetime(ib_v,ik_ibz,spin)
591            do iq=1,nq
592              ctemp = opt_cvk(ib_c,ib_v,ik_bz,spin,iq)
593              do iw=1,nomega
594                eps_rpa(iw,iq) = eps_rpa(iw,iq)  + ctemp * CONJG(ctemp) *&
595 &             (one/(ediff-j_dpc*lifetime-omega(iw)) + one/(ediff+j_dpc*lifetime+omega(iw)))
596              end do
597            end do
598            !
599            ! The JDOS at q=0
600            !if (ediff*Ha_eV < 0.3) then
601            !  write(std_out,*)"Small transition ",ik_ibz,ib_v,ib_c
602            !end if
603 
604            do iw=1,nomega
605              arg = DBLE(omega(iw)) - ediff
606              dos(iw) = dos(iw) + dirac_delta(arg,lifetime)
607            end do
608          else
609            do iq=1,nq
610              ctemp = opt_cvk(ib_c,ib_v,ik_bz,spin,iq)
611              do iw=1,nomega
612                eps_rpa(iw,iq) = eps_rpa(iw,iq)  + ctemp * CONJG(ctemp) *&
613 &             (one/(ediff-omega(iw)) + one/(ediff+omega(iw)))
614              end do
615            end do
616            !
617            ! The JDOS at q=0
618            !if (ediff*Ha_eV < 0.3) then
619            !  write(std_out,*)"Small transition ",ik_ibz,ib_v,ib_c
620            !end if
621 
622            do iw=1,nomega
623              arg = DBLE(omega(iw)) - ediff
624              dos(iw) = dos(iw) + dirac_delta(arg,broad)
625            end do
626          end if
627          !
628        end do !ib_c
629      end do !ib_v
630    end do !ik_bz
631  end do !spin
632 
633  dos = dos/Kmesh%nbz
634  eps_rpa = cone + fact*eps_rpa
635 
636 end subroutine exc_eps_rpa

m_exc_spectra/exc_write_data [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  exc_write_data

FUNCTION

  This routine drives the writing of the files produced by the Bethe-Salpeter code.

INPUTS

 BSp<excparam>=Bethe-Salpeter Parameters.
 what= "EXC_MDF"
       "RPA_NLF_MDF"
       "GW_NLF_MDF"
 [dos(nomega)]

OUTPUT

  Only writing.

SIDE EFFECTS

  eps(BSp%nomega,BSp%nq) = Macroscopic dielectric function to be written.

PARENTS

      m_exc_spectra,m_haydock

CHILDREN

      c_f_pointer

SOURCE

351 subroutine exc_write_data(BSp,BS_files,what,eps,prefix,dos)
352 
353 
354 !This section has been created automatically by the script Abilint (TD).
355 !Do not modify the following lines by hand.
356 #undef ABI_FUNC
357 #define ABI_FUNC 'exc_write_data'
358 !End of the abilint section
359 
360  implicit none
361 
362 !Arguments ------------------------------------
363 !scalars
364  character(len=*),intent(in) :: what
365  type(excparam),intent(in) :: BSp
366  type(excfiles),intent(in) :: BS_files
367  character(len=*),optional,intent(in) :: prefix
368 !arrays
369  real(dp),optional,intent(in) :: dos(BSp%nomega)
370  complex(dpc),intent(in) :: eps(BSp%nomega,BSp%nq)
371 
372 !Local variables ------------------------------
373 !scalars
374  integer :: io,iq,funt
375  real(dp) :: omegaev,step
376  !real(dp),parameter :: SMALL=5.0d-99
377 !arrays
378  real(dp) :: int_dos(BSp%nomega)
379  real(dp) :: tmp_eps(2,BSp%nq)
380  character(len=500) :: lf_type,block_type,wgg_type,frm,str_type,msg
381  character(len=fnlen) :: fname
382 
383 !************************************************************************
384 
385  if (PRESENT(prefix)) then
386    fname = strcat(BS_files%out_basename,prefix,'_',toupper(what))
387  else
388    fname = strcat(BS_files%out_basename,'_',toupper(what))
389  end if
390 
391  if (open_file(fname,msg,newunit=funt,form="formatted", action="write") /= 0) then
392    MSG_ERROR(msg)
393  end if
394 
395  select case (toupper(what))
396  case ("EXC_MDF")
397    call wrtout(ab_out," Writing EXC Macroscopic dielectric function to file: "//trim(fname),"COLL")
398 
399    write(funt,'("# Macroscopic dielectric function obtained with the BS equation.")')
400 
401    lf_type = 'WITHOUT LOCAL FIELD EFFECTS'
402    if (BSp%exchange_term>0) lf_type='LOCAL FIELD EFFECTS INCLUDED'
403    call bsp_calctype2str(Bsp,str_type)
404    write(funt,'("# ",a,"     " ,a)') TRIM(str_type), TRIM(lf_type)
405 
406    block_type = 'RESONANT-ONLY calculation'
407    if (BSp%use_coupling>0) block_type = 'RESONANT+COUPLING calculation'
408    write(funt,'("# ",a)') TRIM(block_type)
409 
410    if (BSp%use_coulomb_term) then
411      wgg_type = "Coulomb term constructed with full W(G1,G2)"
412      if ( BSp%use_diagonal_Wgg ) wgg_type = "Coulomb term constructed with diagonal approximation W(G1,G1)"
413      write(funt,'("# ",a)') TRIM(wgg_type)
414    end if
415 
416    write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
417 
418  case ("RPA_NLF_MDF")
419    call wrtout(ab_out," Writing KS-RPA macroscopic dielectric function without local fields to file: "//trim(fname),"COLL")
420    write(funt,'("# RPA macroscopic dielectric function without local fields")')
421 
422  case ("GW_NLF_MDF")
423    call wrtout(ab_out," Writing GW-RPA macroscopic dielectric function without local fields to file: "//trim(fname),"COLL")
424 
425    write(funt,'("# GW Macroscopic dielectric function without local field effects ")')
426    write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
427 
428  case default
429    MSG_ERROR("Unknown value for what: "//trim(what))
430  end select
431  !
432  ! Paramaters common to the different calculations.
433  if (BSp%algorithm /= BSE_ALGO_HAYDOCK) then
434    write(funt,'(a,i0)')"# nstates included in the diagonalization = ",BSp%nstates
435  end if
436 
437  if (BSp%algorithm == BSE_ALGO_HAYDOCK) then
438    write(funt,'(a,2f7.4)')'# Tolerance = ',BSp%haydock_tol
439  end if
440 
441  write(funt,'(a,i0)')"# npweps  = ",BSp%npweps
442  write(funt,'(a,i0)')"# npwwfn  = ",BSp%npwwfn
443  write(funt,'(a,i0)')"# nbands  = ",BSp%nbnds
444  write(funt,'(a,i0)')"# loband  = ",BSp%lomo_spin(1)
445  if (Bsp%nsppol==2) write(funt,'(a,i0)')"# loband(spin=2) = ",BSp%lomo_spin(2)
446  write(funt,'(a,i0)')"# nkibz   = ",BSp%nkibz
447  write(funt,'(a,i0)')"# nkbz    = ",BSp%nkbz
448  write(funt,'(a,f7.4,a)')'# Lorentzian broadening = ',BSp%broad*Ha_eV,' [eV]'
449  !
450  ! Write the list of q-points.
451  write(funt,'(a)')"# List of q-points for the optical limit:"
452  do iq=1,BSp%nq
453    write(funt,'(a,3(f9.6,","),a)')'# q = ',BSp%q(:,iq),' [Reduced coords] '
454  end do
455  !
456  ! Write spectra.
457  if (.not.PRESENT(dos)) then
458    write(funt,'(a)')"# omega [eV]    RE(eps(q=1)) IM(eps(q=1) RE(eps(q=2) ) ... "
459    !write(frm,*)'(f7.3,',2*BSp%nq,'es12.4)'
460    write(frm,*)'(f7.3,',2*BSp%nq,'(1x,f9.4))'
461    do io=1,BSp%nomega
462      omegaev = DBLE(BSp%omega(io))*Ha_eV
463      tmp_eps(1,:) = REAL (eps(io,:))
464      tmp_eps(2,:) = AIMAG(eps(io,:))
465      !where (ABS(tmp_eps) < SMALL) ! this to improve the portability of the automatic tests.
466      !  tmp_eps = zero
467      !end where
468      write(funt,frm) omegaev,(tmp_eps(:,iq), iq=1,BSp%nq)
469    end do
470 
471  else
472    write(funt,'(a)')"# omega [eV]    RE(eps(q=1)) IM(eps(q=1) RE(eps(q=2) ) ... DOS   IDOS"
473    step = DBLE(BSp%omega(2) - BSp%omega(1))
474    if ( ABS( step - DBLE((BSp%omega(BSp%nomega) - BSp%omega(BSp%nomega-1)))) > tol6 ) then
475      MSG_WARNING("Frequency mesh must be linear for using simpson_int")
476    end if
477    call simpson_int(Bsp%nomega,step,dos,int_dos)
478    !write(frm,*)'(f7.3,',2*BSp%nq,'es12.4,2es12.4)'
479    write(frm,*)'(f7.3,',2*BSp%nq,'(1x,f9.4,1x,f9.4,1x,f9.4))'
480    do io=1,BSp%nomega
481      omegaev = DBLE(BSp%omega(io))*Ha_eV
482      tmp_eps(1,:) = REAL (eps(io,:))
483      tmp_eps(2,:) = AIMAG(eps(io,:))
484      !where (ABS(tmp_eps) < SMALL) ! this to improve the portability of the automatic tests.
485      !  tmp_eps = zero
486      !end where
487      !write(funt,frm) omegaev,(eps(io,iq), iq=1,BSp%nq), dos(io), int_dos(io)
488      write(funt,frm) omegaev,(tmp_eps(:,iq), iq=1,BSp%nq), dos(io), int_dos(io)
489    end do
490  end if
491 
492  close(funt)
493 
494 end subroutine exc_write_data

m_exc_spectra/exc_write_tensor [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

  exc_write_tensor

FUNCTION

  This routine drives the writing of complex dielectric tensor

INPUTS

 BSp<excparam>=Bethe-Salpeter Parameters.
 what= "EXC_TSR_CART" or "EXC_TSR_RED"
       "RPA_NLF_TSR_CART" or "RPA_NLF_TSR_RED"
       "GW_NLF_TSR_CART" or "GW_NLF_TSR_RED"

OUTPUT

  Only writing.

SIDE EFFECTS

  tensor(BSp%nomega,6) = Complex dielectric tensor to be written

PARENTS

      m_haydock

CHILDREN

      c_f_pointer

SOURCE

1129 subroutine exc_write_tensor(BSp,BS_files,what,tensor)
1130 
1131 
1132 !This section has been created automatically by the script Abilint (TD).
1133 !Do not modify the following lines by hand.
1134 #undef ABI_FUNC
1135 #define ABI_FUNC 'exc_write_tensor'
1136 !End of the abilint section
1137 
1138  implicit none
1139 
1140 !Arguments ------------------------------------
1141 !scalars
1142  character(len=*),intent(in) :: what
1143  type(excparam),intent(in) :: BSp
1144  type(excfiles),intent(in) :: BS_files
1145 !arrays
1146  complex(dpc),intent(in) :: tensor(BSp%nomega,6)
1147 
1148 !Local variables ------------------------------
1149 !scalars
1150  integer :: io,iq,funt
1151  real(dp) :: omegaev
1152 !arrays
1153  character(len=500) :: lf_type,block_type,wgg_type,frm,str_type, msg
1154  character(len=fnlen) :: fname
1155 
1156 !************************************************************************
1157 
1158  fname = strcat(BS_files%out_basename,'_',toupper(what))
1159  if (open_file(fname,msg,newunit=funt,form="formatted", action="write") /= 0) then
1160    MSG_ERROR(msg)
1161  end if
1162 
1163  select case (toupper(what))
1164 
1165  case ("EXC_TSR_CART")
1166    write(funt,'("# Complex dielectric tensor (cart. coord.) obtained with the BS equation.")')
1167 
1168    lf_type = 'WITHOUT LOCAL FIELD EFFECTS'
1169    if (BSp%exchange_term>0) lf_type='LOCAL FIELD EFFECTS INCLUDED'
1170    call bsp_calctype2str(Bsp,str_type)
1171    write(funt,'("# ",a,"     " ,a)') TRIM(str_type), TRIM(lf_type)
1172 
1173    block_type = 'RESONANT-ONLY calculation'
1174    if (BSp%use_coupling>0) block_type = 'RESONANT+COUPLING calculation'
1175    write(funt,'("# ",a)') TRIM(block_type)
1176 
1177    if (BSp%use_coulomb_term) then
1178      wgg_type = "Coulomb term constructed with full W(G1,G2)"
1179      if ( BSp%use_diagonal_Wgg ) wgg_type = "Coulomb term constructed with diagonal approximation W(G1,G1)"
1180      write(funt,'("# ",a)') TRIM(wgg_type)
1181    end if
1182 
1183    write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1184 
1185  case ("EXC_TSR_RED")
1186    write(funt,'("# Complex dielectric tensor (red. coord.) obtained with the BS equation.")')
1187 
1188    lf_type = 'WITHOUT LOCAL FIELD EFFECTS'
1189    if (BSp%exchange_term>0) lf_type='LOCAL FIELD EFFECTS INCLUDED'
1190    call bsp_calctype2str(Bsp,str_type)
1191    write(funt,'("# ",a,"     " ,a)') TRIM(str_type), TRIM(lf_type)
1192 
1193    block_type = 'RESONANT-ONLY calculation'
1194    if (BSp%use_coupling>0) block_type = 'RESONANT+COUPLING calculation'
1195    write(funt,'("# ",a)') TRIM(block_type)
1196 
1197    if (BSp%use_coulomb_term) then
1198      wgg_type = "Coulomb term constructed with full W(G1,G2)"
1199      if ( BSp%use_diagonal_Wgg ) wgg_type = "Coulomb term constructed with diagonal approximation W(G1,G1)"
1200      write(funt,'("# ",a)') TRIM(wgg_type)
1201    end if
1202 
1203    write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1204 
1205  case ("RPA_NLF_TSR_CART")
1206    write(funt,'("# RPA complex dielectric tensor (cart. coord.) without local fields")')
1207 
1208  case ("RPA_NLF_TSR_RED")
1209    write(funt,'("# RPA complex dielectric tensor (red. coord.) without local fields")')
1210 
1211  case ("GW_NLF_TSR_CART")
1212    write(funt,'("# GW complex dielectric tensor (cart. coord.) without local field effects ")')
1213    write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1214 
1215  case ("GW_NLF_TSR_RED")
1216    write(funt,'("# GW complex dielectric tensor (red. coord.) without local field effects ")')
1217    write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1218 
1219  case default
1220    MSG_ERROR("Unknown value for what: "//TRIM(what))
1221  end select
1222  !
1223  ! Paramaters common to the different calculations.
1224  if (BSp%algorithm /= BSE_ALGO_HAYDOCK) then
1225    write(funt,'(a,i0)')"# nstates included in the diagonalization = ",BSp%nstates
1226  end if
1227 
1228  if (BSp%algorithm == BSE_ALGO_HAYDOCK) then
1229    write(funt,'(a,2f7.4)')'# Tolerance = ',BSp%haydock_tol
1230  end if
1231 
1232  write(funt,'(a,i0)')"# npweps  = ",BSp%npweps
1233  write(funt,'(a,i0)')"# npwwfn  = ",BSp%npwwfn
1234  write(funt,'(a,i0)')"# nbands  = ",BSp%nbnds
1235  write(funt,'(a,i0)')"# loband  = ",BSp%lomo_spin(1)
1236  if (Bsp%nsppol==2) write(funt,'(a,i0)')"# loband(spin=2) = ",BSp%lomo_spin(2)
1237  write(funt,'(a,i0)')"# nkibz   = ",BSp%nkibz
1238  write(funt,'(a,i0)')"# nkbz    = ",BSp%nkbz
1239  write(funt,'(a,f7.4,a)')'# Lorentzian broadening = ',BSp%broad*Ha_eV,' [eV]'
1240 
1241  !
1242  ! Write tensor.
1243  write(funt,'(3a)') "# omega [eV] RE(eps_11) IM(eps_11) RE(eps_22)", &
1244 &  "IM(eps_22) RE(eps_33) IM(eps_33) RE(eps_12) IM(eps_12)", &
1245 &  "RE(eps_13) IM(eps_13) RE(eps_23) IM(eps_23))"
1246  write(frm,*) '(f7.3,12es14.6)'
1247  do io=1,BSp%nomega
1248    omegaev = DBLE(BSp%omega(io))*Ha_eV
1249    write(funt,frm) omegaev,(tensor(io,iq), iq=1,6)
1250  end do
1251 
1252  close(funt)
1253 
1254 end subroutine exc_write_tensor

m_exc_spectra/mdfs_ncwrite [ Functions ]

[ Top ] [ m_exc_spectra ] [ Functions ]

NAME

 mdfs_ncwrite

FUNCTION

  Writes the MDF.nc file with the final results.

INPUTS

  ncid =NC file handle
  Bsp<excparam>=Data type gathering the paramenters used for the Bethe-Salpeter calculation.
  eps_exc = Excitonic MDF
  eps_rpanlf = KS-RPA MDF without local-field effects.
  eps_gwnlf = GW-RPA MDF without local-field effects.

OUTPUT

  Only writing.

PARENTS

      m_exc_spectra,m_haydock

CHILDREN

      c_f_pointer

SOURCE

1284 subroutine mdfs_ncwrite(ncid,Bsp,eps_exc,eps_rpanlf,eps_gwnlf)
1285 
1286 
1287 !This section has been created automatically by the script Abilint (TD).
1288 !Do not modify the following lines by hand.
1289 #undef ABI_FUNC
1290 #define ABI_FUNC 'mdfs_ncwrite'
1291 !End of the abilint section
1292 
1293  implicit none
1294 
1295 !Arguments ------------------------------------
1296 !scalars
1297  integer,intent(in) :: ncid
1298  type(excparam),intent(in) :: BSp
1299 !arrays
1300  complex(dpc),target,intent(in) :: eps_exc(BSp%nomega,BSp%nq)
1301  complex(dpc),target,intent(in) :: eps_rpanlf(BSp%nomega,BSp%nq)
1302  complex(dpc),target,intent(in) :: eps_gwnlf(BSp%nomega,BSp%nq)
1303 
1304 !Local variables-------------------------------
1305 !scalars
1306 #ifdef HAVE_NETCDF
1307  integer :: ncerr
1308  real(dp), ABI_CONTIGUOUS pointer :: rvals(:,:,:)
1309 
1310 ! *************************************************************************
1311  ! =========================
1312  ! === Write the dimensions
1313  ! =========================
1314 
1315  ncerr = nctk_defnwrite_ivars(ncid, [character(len=nctk_slen) :: &
1316 &  "mdf_version", "nsppol", "npwwfn", "npweps", "nkibz", "nkbz",&
1317 &  "nkibz_iterp", "nkbz_interp", "wtype", "interp_mode"],&
1318 &  [1, Bsp%nsppol, Bsp%npwwfn, Bsp%npweps, Bsp%nkibz, Bsp%nkbz, &
1319 &   Bsp%nkibz_interp, Bsp%nkbz_interp,Bsp%wtype, Bsp%interp_mode])
1320  NCF_CHECK(ncerr)
1321 
1322  ncerr = nctk_defnwrite_dpvars(ncid, [character(len=nctk_slen) :: &
1323 &  "ecutwfn", "ecuteps", "mbpt_sciss", "broad", "eps_inf"],&
1324 &  [Bsp%ecutwfn, Bsp%ecuteps, Bsp%mbpt_sciss, Bsp%broad, Bsp%eps_inf])
1325  NCF_CHECK(ncerr)
1326 
1327  ncerr = nctk_def_dims(ncid, [nctkdim_t("two", 2), nctkdim_t("three", 3), nctkdim_t("number_of_qpoints", Bsp%nq),&
1328    nctkdim_t("number_of_frequencies", Bsp%nomega), nctkdim_t("number_of_spins", bsp%nsppol)], defmode=.True.)
1329  NCF_CHECK(ncerr)
1330 
1331 ! Define variables.
1332 
1333 !arrays
1334  ncerr = nctk_def_arrays(ncid, [&
1335    nctkarr_t('qpoints', "dp", 'three, number_of_qpoints'),&
1336    nctkarr_t('wmesh', "dp", 'number_of_frequencies'),&
1337    nctkarr_t('nreh', "i", "number_of_spins"),&
1338    nctkarr_t('lomo_spin', "i", "number_of_spins"),&
1339    nctkarr_t('humo_spin', "i", "number_of_spins"),&
1340    nctkarr_t('exc_mdf', "dp", 'two, number_of_frequencies, number_of_qpoints'),&
1341    nctkarr_t('rpanlf_mdf', "dp", 'two, number_of_frequencies, number_of_qpoints'),&
1342    nctkarr_t('gwnlf_mdf', "dp", 'two, number_of_frequencies, number_of_qpoints')])
1343  NCF_CHECK(ncerr)
1344 
1345 ! Write data.
1346  NCF_CHECK(nctk_set_datamode(ncid))
1347  NCF_CHECK(nf90_put_var(ncid, vid('qpoints'), Bsp%q))
1348  NCF_CHECK(nf90_put_var(ncid, vid('nreh'), bsp%nreh))
1349  NCF_CHECK(nf90_put_var(ncid, vid('lomo_spin'), bsp%lomo_spin))
1350  NCF_CHECK(nf90_put_var(ncid, vid('humo_spin'), bsp%humo_spin))
1351 
1352  ! Write frequency in mesh in eV.
1353  NCF_CHECK(nf90_put_var(ncid, vid('wmesh'), REAL(Bsp%omega)*Ha_eV))
1354 
1355  call c_f_pointer(c_loc(eps_exc(1,1)), rvals, shape=[2, bsp%nomega, bsp%nq])
1356  NCF_CHECK(nf90_put_var(ncid, vid('exc_mdf'), rvals))
1357 
1358  call c_f_pointer(c_loc(eps_rpanlf(1,1)), rvals, shape=[2, bsp%nomega, bsp%nq])
1359  NCF_CHECK(nf90_put_var(ncid, vid('rpanlf_mdf'), rvals))
1360 
1361  call c_f_pointer(c_loc(eps_gwnlf(1,1)), rvals, shape=[2, bsp%nomega, bsp%nq])
1362  NCF_CHECK(nf90_put_var(ncid, vid("gwnlf_mdf"), rvals))
1363 
1364 #else
1365  MSG_ERROR("ETSF-IO support is not activated.")
1366 #endif
1367 
1368 contains
1369  integer function vid(vname)
1370 
1371 
1372 !This section has been created automatically by the script Abilint (TD).
1373 !Do not modify the following lines by hand.
1374 #undef ABI_FUNC
1375 #define ABI_FUNC 'vid'
1376 !End of the abilint section
1377 
1378    character(len=*),intent(in) :: vname
1379    vid = nctk_idname(ncid, vname)
1380  end function vid
1381 
1382 end subroutine mdfs_ncwrite