TABLE OF CONTENTS


ABINIT/m_optic_tools [ Modules ]

[ Top ] [ Modules ]

NAME

 m_optic_tools

FUNCTION

  Helper functions used in the optic code

COPYRIGHT

 Copyright (C) 2002-2018 ABINIT group (SSharma,MVer,VRecoules,TD,YG)
 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 .

 COMMENTS

  Right now the routine sums over the k-points. In future linear tetrahedron method might be useful.
  Reference articles:
  1. S. Sharma, J. K. Dewhurst and C. Ambrosch-Draxl, Phys. Rev. B {\bf 67} 165332 2003 [[cite:Sharma2003]]
  2. J. L. P. Hughes and J. E. Sipe, Phys. Rev. B {\bf 53} 10 751 1996 [[cite:Hughes1996]]
  3. S. Sharma and C. Ambrosch-Draxl, Physica Scripta T 109 2004 [[cite:Sharma2004]]
  4. J. E. Sipe and Ed. Ghahramani, Phys. Rev. B {\bf 48} 11 705 1993 [[cite:Sipe1993]]

SOURCE

27 #if defined HAVE_CONFIG_H
28 #include "config.h"
29 #endif
30 
31 #include "abi_common.h"
32 
33 MODULE m_optic_tools
34 
35  use defs_basis
36  use m_errors
37  use m_abicore
38  use m_linalg_interfaces
39  use m_xmpi
40  use m_nctk
41 #ifdef HAVE_NETCDF
42  use netcdf
43 #endif
44 
45  use defs_datatypes,    only : ebands_t
46  use m_numeric_tools,   only : wrap2_pmhalf, c2r
47  use m_io_tools,        only : open_file
48 
49  implicit none
50 
51  private
52 
53  public :: sym2cart
54  public :: getwtk
55  public :: pmat2cart
56  public :: pmat_renorm
57  public :: linopt           ! Compute dielectric function for semiconductors
58  public :: nlinopt          ! Second harmonic generation susceptibility for semiconductors
59  public :: linelop          ! Linear electro-optic susceptibility for semiconductors
60  public :: nonlinopt
61 
62 CONTAINS  !===========================================================

m_optic_tools/getwtk [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 getwtk

FUNCTION

 Routine called by the program optic
 Presumes kpts are the irreducible ones of a good uniform grid

INPUTS

  kpt(3,nkpt)=reduced coordinates of k points.
  nkpt = number of k points
  nsym=Number of symmetry operations.
  symrel(3,3,nsym)=symmetry operations

OUTPUT

  wtk(nkpt)=weight assigned to each k point.

PARENTS

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

162 subroutine getwtk(kpt,nkpt,nsym,symrel,wtk)
163 
164 
165 !This section has been created automatically by the script Abilint (TD).
166 !Do not modify the following lines by hand.
167 #undef ABI_FUNC
168 #define ABI_FUNC 'getwtk'
169 !End of the abilint section
170 
171  implicit none
172 
173 !Arguments -----------------------------------------------
174 ! in
175 ! out
176 !scalars
177  integer,intent(in) :: nkpt,nsym
178 !arrays
179  integer,intent(in) :: symrel(3,3,nsym)
180  real(dp),intent(in) :: kpt(3,nkpt)
181  real(dp),intent(out) :: wtk(nkpt)
182 
183 !Local variables -----------------------------------------
184 !scalars
185  integer :: ikpt,istar,isym,itim,new,nkpt_tot
186  real(dp) :: shift,timsign,tmp
187 !arrays
188  integer :: nstar(nkpt)
189  real(dp) :: dkpt(3),kptstar(3,2*nkpt*nsym),rsymrel(3,3,nsym),symkpt(3)
190  real(dp) :: tsymkpt(3)
191 
192 ! *************************************************************************
193 
194  do isym=1,nsym
195    rsymrel(:,:,isym) = dble(symrel(:,:,isym))
196  end do
197 
198 !for each kpt find star and accumulate nkpts
199  do ikpt=1,nkpt
200    write(std_out,*) ' getwtk : ikpt = ', ikpt
201    nstar(ikpt) = 0
202    kptstar(:,:) = zero
203    do isym=1,nsym
204 
205      call dgemv('N',3,3,one,rsymrel(:,:,isym),3,kpt(:,ikpt),1,zero,symkpt,1)
206 
207 !    is symkpt already in star?
208      do itim=0,1
209        timsign=one-itim*two
210        tsymkpt(:) = timsign*symkpt(:)
211        call wrap2_pmhalf(tsymkpt(1),tmp,shift) ;  tsymkpt(1) = tmp
212        call wrap2_pmhalf(tsymkpt(2),tmp,shift) ;  tsymkpt(2) = tmp
213        call wrap2_pmhalf(tsymkpt(3),tmp,shift) ;  tsymkpt(3) = tmp
214        new=1
215        do istar=1,nstar(ikpt)
216          dkpt(:) = abs(tsymkpt(:)-kptstar(:,istar))
217          if ( sum(dkpt) < 1.0d-6) then
218            new=0
219            exit
220          end if
221        end do
222        if (new==1) then
223          nstar(ikpt) = nstar(ikpt)+1
224          kptstar(:,nstar(ikpt)) = tsymkpt(:)
225        end if
226      end do
227 
228    end do
229 !  end do nsym
230 !  DEBUG
231 !  write(std_out,*) ' getwtk : nstar = ', nstar(ikpt)
232 !  write(std_out,*) ' getwtk : star = '
233 !  write(std_out,*)  kptstar(:,1:nstar(ikpt))
234 !  ENDDEBUG
235  end do
236 !end do nkpt
237 
238  nkpt_tot = sum(nstar)
239 !write(std_out,*) ' getwtk : nkpt_tot = ', nkpt_tot
240  do ikpt=1,nkpt
241    wtk(ikpt) = dble(nstar(ikpt))/dble(nkpt_tot)
242  end do
243 
244 end subroutine getwtk

m_optic_tools/linelop [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 linelop

FUNCTION

 Compute optical frequency dependent linear electro-optic susceptibility for semiconductors

INPUTS

  icomp=Sequential index associated to computed tensor components (used for netcdf output)
  itemp=Temperature index (used for netcdf output)
  nspin = number of spins(integer)
  omega = crystal volume in au (real)
  nkpt  = total number of kpoints (integer)
  wkpt(nkpt) = weights of kpoints (real)
  nsymcrys = number of crystal symmetry operations(integer)
  symcrys(3,3,nsymcrys) = symmetry operations in cartisian coordinates(real)
  nstval = total number of valence states(integer)
  evalv(nstval,nspin,nkpt) = eigen value for each band in Ha(real)
  occv(nstval,nspin,nkpt) = occupation number
  efermi = Fermi energy in Ha(real)
  pmat(nstval,nstval,nkpt,3,nspin) = momentum matrix elements in cartesian coordinates(complex)
  v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z
  nmesh = desired number of energy mesh points(integer)
  de = desired step in energy(real); nmesh*de=maximum energy for plotting
  sc = scissors shift in Ha(real)
  brod = broadening in Ha(real)
  tol = tolerance:how close to the singularity exact exact is calculated(real)
  fnam=root for filenames that will contain the output  :
   fnam1=trim(fnam)//'-ChiTotIm.out'
   fnam2=trim(fnam)//'-ChiTotRe.out'
   fnam3=trim(fnam)//'-ChiIm.out'
   fnam4=trim(fnam)//'-ChiRe.out'
   fnam5=trim(fnam)//'-ChiAbs.out'
  ncid=Netcdf id to save output data.

OUTPUT

  Calculates the second harmonic generation susceptibility on a desired energy mesh and
  for desired direction of polarisation. The output is in files named
  ChiEOTot.out : Im\chi_{v1v2v3}(\omega,\omega,0) and Re\chi_{v1v2v3}(\omega,\omega,0)
  ChiEOIm.out  : contributions to the Im\chi_{v1v2v3}(\omega,\omega,0) from various terms
  ChiEORe.out  : contributions to Re\chi_{v1v2v3}(\omega,\omega,-0) from various terms
  ChiEOAbs.out : abs\chi_{v1v2v3}(\omega,\omega,0). The headers in these files contain
  information about the calculation.

  NOTES:
    - The routine has been written using notations of Ref. 2
    - This routine does not symmetrize the tensor (up to now)
    - Sum over all the states and use occupation factors instead of looping only on resonant contributions

PARENTS

      optic

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

1626 subroutine linelop(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,evalv,occv,efermi, &
1627   pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,do_antiresonant,ncid,comm)
1628 
1629 
1630 !This section has been created automatically by the script Abilint (TD).
1631 !Do not modify the following lines by hand.
1632 #undef ABI_FUNC
1633 #define ABI_FUNC 'linelop'
1634 !End of the abilint section
1635 
1636  implicit none
1637 
1638 !Arguments ------------------------------------
1639 !no_abirules
1640 integer, intent(in) :: icomp,itemp,nspin, ncid
1641 real(dp), intent(in) :: omega
1642 integer, intent(in) :: nkpt
1643 real(dp), intent(in) :: wkpt(nkpt)
1644 integer, intent(in) :: nsymcrys
1645 real(dp), intent(in) :: symcrys(3,3,nsymcrys)
1646 integer, intent(in) :: nstval
1647 real(dp), intent(in) :: evalv(nstval,nkpt,nspin)
1648 real(dp), intent(in) :: occv(nstval,nkpt,nspin)
1649 real(dp), intent(in) :: efermi
1650 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin)
1651 integer, intent(in) :: v1
1652 integer, intent(in) :: v2
1653 integer, intent(in) :: v3
1654 integer, intent(in) :: nmesh
1655 integer, intent(in) :: comm
1656 real(dp), intent(in) :: de
1657 real(dp), intent(in) :: sc
1658 real(dp), intent(in) :: brod
1659 real(dp), intent(in) :: tol
1660 character(len=*), intent(in) :: fnam
1661 logical, intent(in) :: do_antiresonant
1662 
1663 !Local variables -------------------------
1664 !no_abirules
1665 ! present calculation related (user specific)
1666 integer :: iw
1667 integer :: i,j,k,lx,ly,lz
1668 integer :: isp,isym,ik
1669 integer :: ist1,istl,istn,istm
1670 real(dp) :: ha2ev
1671 real(dp) :: t1,t2,t3,tst
1672 real(dp) :: ene,totre,totabs,totim
1673 real(dp) :: el,en,em
1674 real(dp) :: emin,emax,my_emin,my_emax
1675 real(dp) :: const_esu,const_au,au2esu
1676 real(dp) :: wmn,wnm,wln,wnl,wml,wlm
1677 complex(dpc) :: idel,w,zi
1678 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5
1679 ! local allocatable arrays
1680 real(dp), allocatable :: s(:,:)
1681 real(dp), allocatable :: sym(:,:,:)
1682 integer :: start4(4),count4(4)
1683 ! DBYG
1684  integer :: istp
1685  real(dp) :: ep, wmp, wpn
1686  real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included !
1687  real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, fmn
1688  complex(dpc), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a}
1689  complex(dpc), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a}
1690  complex(dpc), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k)
1691  complex(dpc), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c
1692  complex(dpc), allocatable :: chi(:) ! \chi_{II}^{abc}(-\omega,\omega,0)
1693  complex(dpc), allocatable :: eta(:) ! \eta_{II}^{abc}(-\omega,\omega,0)
1694  complex(dpc), allocatable :: sigma(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0)
1695  complex(dpc), allocatable :: chi2tot(:)
1696  complex(dpc) :: num1, num2, den1, den2, term1, term2
1697  complex(dpc) :: chi1, chi1_1, chi1_2, chi2_1b, chi2_2b
1698  complex(dpc), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega)
1699  complex(dpc) :: eta1, eta2, eta2_1, eta2_2
1700  complex(dpc) :: sigma1, sigma1_1, sigma1_2, sigma2
1701  !Parallelism
1702  integer :: my_rank, nproc, master=0
1703  integer :: ierr
1704  integer :: my_k1, my_k2
1705  character(500) :: msg
1706  integer :: fout1,fout2,fout3,fout4,fout5
1707 
1708 ! *********************************************************************
1709 
1710  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
1711 
1712 !calculate the constant
1713  zi=(0._dp,1._dp)
1714  idel=zi*brod
1715 ! Disable symmetries for now
1716  const_au=-2._dp/(omega*dble(nsymcrys))
1717  au2esu=5.8300348177d-8
1718  const_esu=const_au*au2esu
1719  ha2ev=13.60569172*2._dp
1720 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1721 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
1722 !bohr: 5.2917ifc nlinopt.f907E-11
1723 !c: 2.99792458   velocity of sound
1724 !ry2ev: 13.60569172
1725 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2)
1726 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
1727 !mass comes from converting P_mn to r_mn
1728 !hbar^3 comes from converting all frequencies to energies in denominator
1729 !hbar^3 comes from operator for momentum (hbar/i nabla)
1730 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1731 !output file names
1732  fnam1=trim(fnam)//'-ChiEOTotIm.out'
1733  fnam2=trim(fnam)//'-ChiEOTotRe.out'
1734  fnam3=trim(fnam)//'-ChiEOIm.out'
1735  fnam4=trim(fnam)//'-ChiEORe.out'
1736  fnam5=trim(fnam)//'-ChiEOAbs.out'
1737 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
1738 ! fool proof:
1739 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
1740 !If there exists inversion symmetry exit with a mesg.
1741  tst=1.d-09
1742  do isym=1,nsymcrys
1743    t1=symcrys(1,1,isym)+1
1744    t2=symcrys(2,2,isym)+1
1745    t3=symcrys(3,3,isym)+1
1746 !  test if diagonal elements are -1
1747    if (abs(t1).lt.tst.and.abs(t2).lt.tst.and.abs(t3).lt.tst) then
1748 !    test if off-diagonal elements are zero
1749      if (abs(symcrys(1,2,isym)).lt.tst.and.abs(symcrys(1,3,isym)).lt.tst &
1750      .and.abs(symcrys(2,1,isym)).lt.tst.and.abs(symcrys(2,3,isym)).lt.tst.and.  &
1751      abs(symcrys(3,1,isym)).lt.tst.and.abs(symcrys(3,2,isym)).lt.tst) then
1752        write(std_out,*) '-----------------------------------------'
1753        write(std_out,*) '    the crystal has inversion symmetry   '
1754        write(std_out,*) '    the LEO susceptibility is zero       '
1755        write(std_out,*) '-----------------------------------------'
1756        MSG_ERROR("Aborting now")
1757      end if
1758    end if
1759  end do
1760 !check polarisation
1761  if (v1.le.0.or.v2.le.0.or.v3.le.0.or.v1.gt.3.or.v2.gt.3.or.v3.gt.3) then
1762    write(std_out,*) '---------------------------------------------'
1763    write(std_out,*) '    Error in linelop:                        '
1764    write(std_out,*) '    the polarisation directions incorrect    '
1765    write(std_out,*) '    1=x,  2=y  and 3=z                       '
1766    write(std_out,*) '---------------------------------------------'
1767    MSG_ERROR("Aborting now")
1768  end if
1769 !number of energy mesh points
1770  if (nmesh.le.0) then
1771    write(std_out,*) '---------------------------------------------'
1772    write(std_out,*) '    Error in linelop:                        '
1773    write(std_out,*) '    number of energy mesh points incorrect   '
1774    write(std_out,*) '    number has to be integer greater than 0  '
1775    write(std_out,*) '    nmesh*de = max energy for calculation    '
1776    write(std_out,*) '---------------------------------------------'
1777    MSG_ERROR("Aborting now")
1778  end if
1779 !step in energy
1780  if (de.le.0._dp) then
1781    write(std_out,*) '---------------------------------------------'
1782    write(std_out,*) '    Error in nlinopt:                        '
1783    write(std_out,*) '    energy step is incorrect                 '
1784    write(std_out,*) '    number has to real greater than 0.0      '
1785    write(std_out,*) '    nmesh*de = max energy for calculation    '
1786    write(std_out,*) '---------------------------------------------'
1787    MSG_ERROR("Aborting now")
1788  end if
1789 !broadening
1790  if (brod.gt.0.009) then
1791    write(std_out,*) '---------------------------------------------'
1792    write(std_out,*) '    ATTENTION: broadening is quite high      '
1793    write(std_out,*) '    ideally should be less than 0.005        '
1794    write(std_out,*) '---------------------------------------------'
1795  else if (brod.gt.0.015) then
1796    write(std_out,*) '----------------------------------------'
1797    write(std_out,*) '    ATTENTION: broadening is too high   '
1798    write(std_out,*) '    ideally should be less than 0.005   '
1799    write(std_out,*) '----------------------------------------'
1800  end if
1801 !tolerance
1802  if (tol.gt.0.006) then
1803    write(std_out,*) '----------------------------------------'
1804    write(std_out,*) '    ATTENTION: tolerance is too high    '
1805    write(std_out,*) '    ideally should be less than 0.004   '
1806    write(std_out,*) '----------------------------------------'
1807  end if
1808 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1809 !fool proof ends
1810 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1811 
1812 !
1813 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1814 !allocate local arrays
1815 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1816 
1817  ABI_MALLOC(enk,(nstval))
1818  ABI_MALLOC(delta,(nstval,nstval,3))
1819  ABI_MALLOC(rmnbc,(nstval,nstval,3,3))
1820  ABI_MALLOC(roverw,(nstval,nstval,3,3))
1821  ABI_MALLOC(rmna,(nstval,nstval,3))
1822  ABI_MALLOC(chi,(nmesh))
1823  ABI_MALLOC(eta,(nmesh))
1824  ABI_MALLOC(sigma,(nmesh))
1825  ABI_MALLOC(chi2,(nmesh))
1826  ABI_MALLOC(sym,(3,3,3))
1827  ABI_MALLOC(s,(3,3))
1828 
1829 !generate the symmetrizing tensor
1830  sym(:,:,:)=0._dp
1831  do isym=1,nsymcrys
1832    s(:,:)=symcrys(:,:,isym)
1833    do i=1,3
1834      do j=1,3
1835        do k=1,3
1836          sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
1837        end do
1838      end do
1839    end do
1840  end do
1841 
1842 !initialise
1843  delta(:,:,:)=0._dp
1844  rmnbc(:,:,:,:)=0._dp
1845  chi(:)=0._dp
1846  chi2(:) = 0._dp
1847  eta(:)=0._dp
1848  sigma(:)=0._dp
1849  my_emin=HUGE(0._dp)
1850  my_emax=-HUGE(0._dp)
1851 
1852  ! Split work
1853  call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr)
1854 
1855 ! loop over kpts
1856  do ik=my_k1,my_k2
1857    write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt
1858 !  loop over spins
1859    do isp=1,nspin
1860 !    Calculate the scissor corrected energies and the energy window
1861      do ist1=1,nstval
1862        en = evalv(ist1,ik,isp)
1863        my_emin=min(my_emin,en)
1864        my_emax=max(my_emax,en)
1865        if(en > efermi) then
1866          en = en + sc
1867        end if
1868        enk(ist1) = en
1869      end do
1870 
1871 !    calculate \Delta_nm and r_mn^a
1872      do istn=1,nstval
1873        en = enk(istn)
1874        do istm=1,nstval
1875          em = enk(istm)
1876          wmn = em - en
1877          delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp)
1878          if(abs(wmn) < tol) then
1879            rmna(istm,istn,1:3) = 0._dp
1880          else
1881            rmna(istm,istn,1:3)=-zi*pmat(istm,istn,ik,1:3,isp)/wmn
1882          end if
1883        end do
1884      end do
1885 !    calculate \r^b_mn;c
1886      do istm=1,nstval
1887        em = enk(istm)
1888        do istn=1,nstval
1889          en = enk(istn)
1890          wmn = em - en
1891          if(abs(wmn) > tol) then
1892            do ly = 1,3
1893              do lz = 1,3
1894                num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly))
1895                den1 = wmn
1896                term1 = num1/den1
1897                term2 = 0._dp
1898                do istp=1,nstval
1899                  ep = enk(istp)
1900                  wmp = em - ep
1901                  wpn = ep - en
1902                  num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly))
1903                  den2 = wmn
1904                  term2 = term2 + (num2/den2)
1905                end do
1906                rmnbc(istm,istn,ly,lz) = -term1-(zi*term2)
1907                roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz)
1908              end do
1909            end do
1910          end if
1911        end do
1912      end do
1913 
1914 !    initialise the factors
1915 !    start the calculation
1916      do istn=1,nstval
1917        en=enk(istn)
1918        if (do_antiresonant .and. en .ge. efermi) then
1919          cycle
1920        end if
1921        fn=occv(istn,ik,isp)
1922        do istm=1,nstval
1923          em=enk(istm)
1924          if (do_antiresonant .and. em .le. efermi) then
1925            cycle
1926          end if
1927          wmn=em-en
1928          wnm=-wmn
1929          fm = occv(istm,ik,isp)
1930          fnm = fn - fm
1931          fmn = fm - fn
1932          eta1 = 0._dp
1933          eta2_1 = 0._dp
1934          eta2_2 = 0._dp
1935          sigma1_1 = 0._dp
1936          sigma1_2 = 0._dp
1937          sigma2 = 0._dp
1938          if(abs(wmn) > tol) then
1939            do lx = 1,3
1940              do ly = 1,3
1941                do lz = 1,3
1942                  eta1 = eta1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*(roverw(istm,istn,lz,ly)))
1943                  eta2_1 = eta2_1 + sym(lx,ly,lz)*(fnm*(rmna(istn,istm,lx)*rmnbc(istm,istn,ly,lz)))
1944                  eta2_2 = eta2_2 + sym(lx,ly,lz)*(fnm*(rmnbc(istn,istm,lx,lz)*rmna(istm,istn,ly)))
1945                  sigma1_1 = sigma1_1 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,ly)*rmna(istm,istn,lz))/(wmn**2)
1946                  sigma1_2 = sigma1_2 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,lz)*rmna(istm,istn,ly))/(wmn**2)
1947                  sigma2 = sigma2 + sym(lx,ly,lz)*(fnm*rmnbc(istn,istm,lz,lx)*rmna(istm,istn,ly))/wmn
1948                end do
1949              end do
1950            end do
1951          end if
1952          chi1_1 = 0._dp
1953          chi1_2 = 0._dp
1954          chi2_1b = 0._dp
1955          chi2_2b = 0._dp
1956          chi2(:) = 0._dp
1957 !        Three band terms
1958          do istl=1,nstval
1959            el=enk(istl)
1960            fl = occv(istl,ik,isp)
1961            wlm = el-em
1962            wln = el-en
1963            wnl = en-el
1964            wml = em-el
1965            fnl = fn-fl
1966            fln = fl-fn
1967            fml = fm-fl
1968            do lx = 1,3
1969              do ly = 1,3
1970                do lz = 1,3
1971                  if(abs(wlm) > tol) then
1972                    chi1_1 = chi1_1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,lz)*rmna(istl,istn,ly))/(wlm)
1973                    chi2_1b = chi2_1b + sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*rmna(istl,istm,lz)*rmna(istm,istn,ly))/(wlm)
1974                  end if
1975                  if(abs(wln) > tol) then
1976                    chi1_2 = chi1_2 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,ly)*rmna(istl,istn,lz))/(wln)
1977                    chi2_2b = chi2_2b + sym(lx,ly,lz)*(fmn*rmna(istl,istm,lx)*rmna(istm,istn,ly)*rmna(istn,istl,lz))/(wnl)
1978                  end if
1979                end do
1980              end do
1981            end do
1982          end do
1983 
1984          sigma1 = 0.5_dp*(sigma1_1-sigma1_2)
1985          eta2 = 0.5_dp*(eta2_1-eta2_2)
1986          chi1 = chi1_1 + chi1_2
1987 !
1988 !        calculate over the desired energy mesh and sum over k-points
1989          do iw=1,nmesh
1990            w=(iw-1)*de+idel
1991            ! Better way to compute it
1992            chi(iw) = chi(iw) + 0.5_dp*wkpt(ik)*((chi1/(wmn-w)) + ((chi2_1b+chi2_2b)/(wmn-w)))*const_esu
1993            eta(iw) = eta(iw) + 0.5_dp*zi*wkpt(ik)*((eta1/(wmn-w)) + (eta2/((wmn-w)**2)))*const_esu
1994            sigma(iw) = sigma(iw) + 0.5_dp*zi*wkpt(ik)*((sigma1/(wmn-w))- (sigma2/(wmn-w)))*const_esu
1995          end do
1996        end do ! istn and istm
1997      end do
1998    end do ! spins
1999  end do ! k-points
2000 
2001  call xmpi_sum(chi,comm,ierr)
2002  call xmpi_sum(eta,comm,ierr)
2003  call xmpi_sum(sigma,comm,ierr)
2004  call xmpi_min(my_emin,emin,comm,ierr)
2005  call xmpi_max(my_emax,emax,comm,ierr)
2006 
2007  if (my_rank == master) then
2008 
2009    if (ncid /= nctk_noid) then
2010      start4 = [1, 1, icomp, itemp]
2011      count4 = [2, nmesh, 1, 1]
2012      ABI_MALLOC(chi2tot, (nmesh))
2013      chi2tot = chi + eta + sigma
2014 #ifdef HAVE_NETCDF
2015      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi"), c2r(chi), start=start4, count=count4))
2016      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_eta"), c2r(eta), start=start4, count=count4))
2017      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_sigma"), c2r(sigma), start=start4, count=count4))
2018      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi2tot"), c2r(chi2tot), start=start4, count=count4))
2019 #endif
2020      ABI_FREE(chi2tot)
2021    end if
2022 
2023   ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
2024    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
2025      MSG_ERROR(msg)
2026    end if
2027    if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
2028      MSG_ERROR(msg)
2029    end if
2030    if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
2031      MSG_ERROR(msg)
2032    end if
2033    if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
2034      MSG_ERROR(msg)
2035    end if
2036    if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
2037      MSG_ERROR(msg)
2038    end if
2039    ! write headers
2040    write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2041    write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
2042    write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
2043    write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
2044    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2045    write(fout1, '(a)' )' # Energy      Tot-Im Chi(-w,w,0)  Tot-Im Chi(-w,w,0)'
2046    write(fout1, '(a)' )' # eV          *10^-7 esu        *10^-12 m/V SI units '
2047    write(fout1, '(a)' )' # '
2048 
2049    write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2050    write(fout2, '(a,es16.6)') ' #tolerance:',tol
2051    write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2052    write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2053    write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2054    write(fout2, '(a)')' # Energy      Tot-Re Chi(-w,w,0)  Tot-Re Chi(-w,w,0)'
2055    write(fout2, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
2056    write(fout2, '(a)')' # '
2057 
2058    write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2059    write(fout3, '(a,es16.6)') ' #tolerance:',tol
2060    write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2061    write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2062    write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2063    write(fout3, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2064    write(fout3, '(a)')' # in esu'
2065    write(fout3, '(a)')' # '
2066 
2067    write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2068    write(fout4, '(a,es16.6)') ' #tolerance:',tol
2069    write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2070    write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2071    write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2072    write(fout4, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2073    write(fout4, '(a)')' # in esu'
2074    write(fout4, '(a)')' # '
2075 
2076    write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2077    write(fout5, '(a,es16.6)') ' #tolerance:',tol
2078    write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2079    write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2080    write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2081    write(fout5, '(a)')' # Energy(eV)  |TotChi(-w,w,0)|   |Tot Chi(-w,w,0)|'
2082    write(fout5, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
2083    write(fout5, '(a)')' # '
2084 
2085    totim=0._dp
2086    totre=0._dp
2087    totabs=0._dp
2088    do iw=2,nmesh
2089      ene=(iw-1)*de
2090      ene=ene*ha2ev
2091      totim=aimag(chi(iw)+eta(iw)+sigma(iw))/1.d-7
2092      write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2093      totim=0._dp
2094      totre=dble(chi(iw)+eta(iw)+sigma(iw))/1.d-7
2095      write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2096      totre=0._dp
2097      write(fout3,'(f15.6,3es15.6)') ene,aimag(chi(iw))/1.d-7,      &
2098      aimag(eta(iw))/1.d-7,aimag(sigma(iw))/1.d-7
2099      write(fout4,'(f15.6,3es15.6)') ene,dble(chi(iw))/1.d-7,       &
2100      dble(eta(iw))/1.d-7,dble(sigma(iw))/1.d-7
2101      totabs=abs(chi(iw)+eta(iw)+sigma(iw))/1.d-7
2102      write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2103      totabs=0._dp
2104    end do
2105 
2106    close(fout1)
2107    close(fout2)
2108    close(fout3)
2109    close(fout4)
2110    close(fout5)
2111    ! print information
2112    write(std_out,*) ' '
2113    write(std_out,*) 'information about calculation just performed:'
2114    write(std_out,*) ' '
2115    write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of LEO susceptibility'
2116    write(std_out,*) 'tolerance:',tol
2117    if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
2118    write(std_out,*) 'broadening:',brod,'Hartree'
2119    if (brod.gt.0.009) then
2120      write(std_out,*) ' '
2121      write(std_out,*) 'ATTENTION: broadening is quite high'
2122      write(std_out,*) ' '
2123    else if (brod.gt.0.015) then
2124      write(std_out,*) ' '
2125      write(std_out,*) 'ATTENTION: broadening is too high'
2126      write(std_out,*) ' '
2127    end if
2128    write(std_out,*) 'scissors shift:',sc,'Hartree'
2129    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree'
2130 
2131  end if
2132 
2133  ! deallocate local arrays
2134  ABI_FREE(enk)
2135  ABI_FREE(delta)
2136  ABI_FREE(rmnbc)
2137  ABI_FREE(roverw)
2138  ABI_FREE(rmna)
2139  ABI_FREE(chi)
2140  ABI_FREE(chi2)
2141  ABI_FREE(eta)
2142  ABI_FREE(sigma)
2143  ABI_FREE(s)
2144  ABI_FREE(sym)
2145 
2146 end subroutine linelop

m_optic_tools/linopt [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 linopt

FUNCTION

 Compute optical frequency dependent dielectric function for semiconductors

INPUTS

  icomp=Sequential index associated to computed tensor components (used for netcdf output)
  itemp=Temperature index (used for netcdf output)
  nspin=number of spins(integer)
  omega=crystal volume in au (real)
  nkpt=total number of kpoints (integer)
  wkpt(nkpt)=weights of kpoints (real)
  nsymcrys=number of crystal symmetry operations(integer)
  symcrys(3,3,nsymcrys)=symmetry operations in cartisian coordinates(real)
  nstval=total number of valence states(integer)
  occv(nstval,nkpt,nspin)=occupation number for each band(real)
  evalv(nstval,nkpt,nspin)=eigen value for each band in Ha(real)
  efermi=Fermi energy in Ha(real)
  pmat(nstval,nstval,nkpt,3,nspin)=momentum matrix elements in cartesian coordinates(complex)
  v1,v2=desired component of the dielectric function(integer) 1=x,2=y,3=z
  nmesh=desired number of energy mesh points(integer)
  de=desired step in energy(real); nmesh*de=maximum energy
  sc=scissors shift in Ha(real)
  brod=broadening in Ha(real)
  fnam=root for filename that will contain the output filename will be trim(fnam)//'-linopt.out'
  ncid=Netcdf id to save output data.

SIDE EFFECTS

  Dielectric function for semiconductors, on a desired energy mesh and for a desired
  direction of polarisation is written to file.
  The output is in a file named trim(fnam)//'-linopt.out' and contains
  Im(\epsilon_{v1v2}(\omega), Re(\epsilon_{v1v2}(\omega) and abs(\epsilon_{v1v2}(\omega).
  Comment:
  Right now the routine sums over the kpoints. In future linear tetrahedron method should be useful.

PARENTS

      optic

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

450 subroutine linopt(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,KSBSt,EPBSt,efermi,pmat, &
451   v1,v2,nmesh,de,sc,brod,fnam,ncid,comm)
452 
453 
454 !This section has been created automatically by the script Abilint (TD).
455 !Do not modify the following lines by hand.
456 #undef ABI_FUNC
457 #define ABI_FUNC 'linopt'
458 !End of the abilint section
459 
460  implicit none
461 
462 !Arguments ------------------------------------
463 !no_abirules
464 integer, intent(in) :: icomp,itemp,nspin,ncid
465 real(dp), intent(in) :: omega
466 integer, intent(in) :: nkpt
467 real(dp), intent(in) :: wkpt(nkpt)
468 integer, intent(in) :: nsymcrys
469 real(dp), intent(in) :: symcrys(3,3,nsymcrys)
470 integer, intent(in) :: nstval
471 type(ebands_t),intent(in) :: KSBSt,EPBSt
472 real(dp), intent(in) :: efermi
473 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin)
474 integer, intent(in) :: v1
475 integer, intent(in) :: v2
476 integer, intent(in) :: nmesh
477 real(dp), intent(in) :: de
478 real(dp), intent(in) :: sc
479 real(dp), intent(in) :: brod
480 character(len=*), intent(in) :: fnam
481 integer, intent(in) :: comm
482 
483 
484 
485 !Local variables -------------------------
486 !no_abirules
487 integer :: isp
488 integer :: i,j,isym,lx,ly,ik
489 integer :: ist1,ist2,iw
490 ! Parallelism
491 integer :: my_rank, nproc
492 integer,parameter :: master=0
493 integer :: my_k1, my_k2
494 #ifdef HAVE_NETCDF
495 integer :: ncerr
496 #endif
497 integer :: ierr
498 integer :: fout1
499 logical :: do_lifetime
500 complex(dpc) :: e1,e2,e12
501 complex(dpc) :: e1_ep,e2_ep,e12_ep
502 real(dp) :: deltav1v2
503 real(dp) :: ha2ev
504 real(dp) :: tmpabs
505 real(dp) :: renorm_factor,emin,emax
506 real(dp) :: ene
507 complex(dpc) :: b11,b12
508 complex(dpc) :: ieta,w
509 character(len=fnlen) :: fnam1
510 character(len=500) :: msg
511 ! local allocatable arrays
512 real(dp) :: s(3,3),sym(3,3)
513 complex(dpc), allocatable :: chi(:)
514 complex(dpc), allocatable :: eps(:)
515 
516 ! *********************************************************************
517 
518  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
519 
520  if (my_rank == master) then
521   !check polarisation
522    if (v1.le.0.or.v2.le.0.or.v1.gt.3.or.v2.gt.3) then
523      write(std_out,*) '---------------------------------------------'
524      write(std_out,*) '    Error in linopt:                         '
525      write(std_out,*) '    the polarisation directions incorrect    '
526      write(std_out,*) '    1=x and 2=y and 3=z                      '
527      write(std_out,*) '---------------------------------------------'
528      MSG_ERROR("Aborting now")
529    end if
530   !number of energy mesh points
531    if (nmesh.le.0) then
532      write(std_out,*) '---------------------------------------------'
533      write(std_out,*) '    Error in linopt:                         '
534      write(std_out,*) '    number of energy mesh points incorrect   '
535      write(std_out,*) '    number has to integer greater than 0     '
536      write(std_out,*) '    nmesh*de = max energy for calculation    '
537      write(std_out,*) '---------------------------------------------'
538      MSG_ERROR("Aborting now")
539    end if
540   !step in energy
541    if (de.le.0._dp) then
542      write(std_out,*) '---------------------------------------------'
543      write(std_out,*) '    Error in linopt:                         '
544      write(std_out,*) '    energy step is incorrect                 '
545      write(std_out,*) '    number has to real greater than 0.0      '
546      write(std_out,*) '    nmesh*de = max energy for calculation    '
547      write(std_out,*) '---------------------------------------------'
548      MSG_ERROR("Aborting now")
549    end if
550   !broadening
551    if (brod.gt.0.009) then
552      write(std_out,*) '---------------------------------------------'
553      write(std_out,*) '    ATTENTION: broadening is quite high      '
554      write(std_out,*) '    ideally should be less than 0.005        '
555      write(std_out,*) '---------------------------------------------'
556    else if (brod.gt.0.015) then
557      write(std_out,*) '----------------------------------------'
558      write(std_out,*) '    ATTENTION: broadening is too high   '
559      write(std_out,*) '    ideally should be less than 0.005   '
560      write(std_out,*) '----------------------------------------'
561    end if
562   !fermi energy
563    if(efermi<-1.0d4) then
564      write(std_out,*) '---------------------------------------------'
565      write(std_out,*) '    ATTENTION: Fermi energy seems extremely  '
566      write(std_out,*) '    low                                      '
567      write(std_out,*) '---------------------------------------------'
568    end if
569   !scissors operator
570    if (sc.lt.0._dp) then
571      write(std_out,*) '---------------------------------------------'
572      write(std_out,*) '    Error in linopt:                         '
573      write(std_out,*) '    scissors shift is incorrect              '
574      write(std_out,*) '    number has to real greater than 0.0      '
575      write(std_out,*) '---------------------------------------------'
576      MSG_ERROR("Aborting now")
577    end if
578 !fool proof end
579  end if
580 
581  ABI_CHECK(KSBSt%mband==nstval, "The number of bands in the BSt should be equal to nstval !")
582 
583  do_lifetime = allocated(EPBSt%lifetime)
584 ! TODO: activate this, and remove do_lifetime - always add it in even if 0.
585 ! if (.not. allocated(EPBSt%lifetime)) then
586 !   ABI_ALLOCATE(EPBSt%lifetime, (nstval, my_k2-my_k1+1, nspin))
587 !   EPBSt%lifetime = zero
588 ! end if
589 
590 !allocate local arrays
591  ABI_ALLOCATE(chi,(nmesh))
592  ABI_ALLOCATE(eps,(nmesh))
593  ieta=(0._dp,1._dp)*brod
594  renorm_factor=1._dp/(omega*dble(nsymcrys))
595  ha2ev=13.60569172*2._dp
596 !output file names
597  fnam1=trim(fnam)//'-linopt.out'
598 !construct symmetrisation tensor
599  sym(:,:)=0._dp
600  do isym=1,nsymcrys
601    s(:,:)=symcrys(:,:,isym)
602    do i=1,3
603      do j=1,3
604        sym(i,j)=sym(i,j)+s(i,v1)*s(j,v2)
605      end do
606    end do
607  end do
608 
609 !calculate the energy window
610  emin=0._dp
611  emax=0._dp
612  do ik=1,nkpt
613    do isp=1,nspin
614      do ist1=1,nstval
615        emin=min(emin,EPBSt%eig(ist1,ik,isp))
616        emax=max(emax,EPBSt%eig(ist1,ik,isp))
617      end do
618    end do
619  end do
620 
621  ! Split work
622  call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr)
623 
624 !start calculating linear optical response
625  chi(:)=0._dp
626 ! TODO: this loop should be outside the ik one, for speed and cache.
627  do isp=1,nspin
628    !do ik=1,nkpt
629    do ik=my_k1,my_k2
630      write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt
631      do ist1=1,nstval
632        e1=KSBSt%eig(ist1,ik,isp)
633        e1_ep=EPBSt%eig(ist1,ik,isp)
634 ! TODO: unless memory is a real issue, should set lifetimes to 0 and do this sum systematically
635 ! instead of putting an if statement in a loop! See above
636        if(do_lifetime) then
637          e1_ep = e1_ep + EPBSt%lifetime(ist1,ik,isp)*(0.0_dp,1.0_dp)
638        end if
639 !      if (e1.lt.efermi) then
640 !      do ist2=ist1,nstval
641        do ist2=1,nstval
642          e2=KSBSt%eig(ist2,ik,isp)
643          e2_ep=EPBSt%eig(ist2,ik,isp)
644          if(do_lifetime) then
645            e2_ep = e2_ep - EPBSt%lifetime(ist2,ik,isp)*(0.0_dp,1.0_dp)
646          end if
647 !        if (e2.gt.efermi) then
648          if (ist1.ne.ist2) then
649 !          scissors correction of momentum matrix
650            if(REAL(e1) > REAL(e2)) then
651              e12 = e1-e2+sc
652            else
653              e12 = e1-e2-sc
654            end if
655            if(REAL(e1_ep) > REAL(e2_ep)) then
656              e12_ep = e1_ep-e2_ep+sc
657            else
658              e12_ep = e1_ep-e2_ep-sc
659            end if
660 !          e12=e1-e2-sc
661            b11=0._dp
662 !          symmetrization of momentum matrix
663            do lx=1,3
664              do ly=1,3
665                b11=b11+(sym(lx,ly)*pmat(ist1,ist2,ik,lx,isp)* &
666                conjg(pmat(ist1,ist2,ik,ly,isp)))
667              end do
668            end do
669            b12=b11*renorm_factor*(1._dp/(e12**2))
670 !          calculate on the desired energy grid
671            do iw=2,nmesh
672              w=(iw-1)*de+ieta
673              chi(iw)=chi(iw)+(wkpt(ik)*(KSBSt%occ(ist1,ik,isp)-KSBSt%occ(ist2,ik,isp))* &
674              (b12/(-e12_ep-w)))
675            end do
676          end if
677        end do ! states
678 !      end if
679      end do
680    end do ! spin
681  end do ! k-points
682 
683  call xmpi_sum(chi,comm,ierr)
684 
685  ! calculate epsilon
686  eps(1) = zero
687  deltav1v2=zero; if (v1 == v2) deltav1v2=one
688  do iw=2,nmesh
689    ene=(iw-1)*de
690    ene=ene*ha2ev
691    eps(iw)=deltav1v2+4._dp*pi*chi(iw)
692  end do
693 
694  if (my_rank == master) then
695    !  open the output files
696    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
697      MSG_ERROR(msg)
698    end if
699    ! write the output
700    write(fout1, '(a,2i3,a)' )' #calculated the component:',v1,v2,'  of dielectric function'
701    write(std_out,*) 'calculated the component:',v1,v2,'  of dielectric function'
702    write(fout1, '(a,2es16.6)' ) ' #broadening:', real(ieta),aimag(ieta)
703    write(std_out,*) ' with broadening:',ieta
704    write(fout1, '(a,es16.6)' ) ' #scissors shift:',sc
705    write(std_out,*) 'and scissors shift:',sc
706    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
707    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
708    write(fout1,*)
709    write(fout1, '(a)' ) ' # Energy(eV)         Im(eps(w))'
710    do iw=2,nmesh
711      ene=(iw-1)*de
712      ene=ene*ha2ev
713      write(fout1, '(2es16.6)' ) ene,aimag(eps(iw))
714    end do
715    write(fout1,*)
716    write(fout1,*)
717    write(fout1, '(a)' ) ' # Energy(eV)         Re(eps(w))'
718    do iw=2,nmesh
719      ene=(iw-1)*de
720      ene=ene*ha2ev
721      write(fout1, '(2es16.6)' ) ene,dble(eps(iw))
722    end do
723    write(fout1,*)
724    write(fout1,*)
725    write(fout1, '(a)' )' # Energy(eV)         abs(eps(w))'
726    do iw=2,nmesh
727      ene=(iw-1)*de
728      ene=ene*ha2ev
729      write(fout1, '(2es16.6)' ) ene,abs(eps(iw))
730    end do
731    write(fout1,*)
732    write(fout1,*)
733    write(fout1, '(a)' )' # Energy(eV)         Im(refractive index(w)) aka kappa'
734    do iw=2,nmesh
735      ene=(iw-1)*de
736      ene=ene*ha2ev
737      write(fout1, '(2es16.6)' ) ene,sqrt(half*(abs(eps(iw)) - dble(eps(iw)) ))
738    end do
739    write(fout1,*)
740    write(fout1,*)
741    write(fout1, '(a)' )' # Energy(eV)         Re(refractive index(w)) aka n'
742    do iw=2,nmesh
743      ene=(iw-1)*de
744      ene=ene*ha2ev
745      write(fout1, '(2es16.6)' ) ene,sqrt(half*(abs(eps(iw)) + dble(eps(iw)) ))
746    end do
747    write(fout1,*)
748    write(fout1,*)
749    write(fout1, '(a)' )' # Energy(eV)         Reflectivity(w) from vacuum, at normal incidence'
750    do iw=2,nmesh
751      ene=(iw-1)*de
752      ene=ene*ha2ev
753      write(fout1, '(2es16.6)' ) ene, sqrt(half*(abs(eps(iw)) + dble(eps(iw)) ))
754    end do
755    write(fout1,*)
756    write(fout1,*)
757    write(fout1, '(a)' )' # Energy(eV)         absorption coeff (in m-1) = omega Im(eps) / c n(eps)'
758    do iw=2,nmesh
759      ene=(iw-1)*de
760      tmpabs=zero
761      if (abs(eps(iw)) + dble(eps(iw)) > zero) then
762        tmpabs = aimag(eps(iw))*ene / sqrt(half*( abs(eps(iw)) + dble(eps(iw)) )) / Sp_Lt / Bohr_meter
763      end if
764      write(fout1, '(2es16.6)' ) ha2ev*ene, tmpabs
765    end do
766 
767    ! close output file
768    close(fout1)
769 
770 #ifdef HAVE_NETCDF
771    if (ncid /= nctk_noid) then
772      ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_epsilon"), c2r(eps), start=[1, 1, icomp, itemp])
773      NCF_CHECK(ncerr)
774    end if
775 #endif
776  end if
777 
778  ABI_DEALLOCATE(chi)
779  ABI_FREE(eps)
780 
781 end subroutine linopt

m_optic_tools/nlinopt [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 nlinopt

FUNCTION

 Compute optical frequency dependent second harmonic generation susceptibility for semiconductors

INPUTS

  icomp=Sequential index associated to computed tensor components (used for netcdf output)
  itemp=Temperature index (used for netcdf output)
  nspin = number of spins(integer)
  omega = crystal volume in au (real)
  nkpt  = total number of kpoints (integer)
  wkpt(nkpt) = weights of kpoints (real)
  nsymcrys = number of crystal symmetry operations(integer)
  symcrys(3,3,nsymcrys) = symmetry operations in cartisian coordinates(real)
  nstval = total number of valence states(integer)
  evalv(nstval,nspin,nkpt) = eigen value for each band in Ha(real)
  efermi = Fermi energy in Ha(real)
  pmat(nstval,nstval,nkpt,3,nspin) = momentum matrix elements in cartesian coordinates(complex)
  v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z
  nmesh = desired number of energy mesh points(integer)
  de = desired step in energy(real); nmesh*de=maximum energy for plotting
  sc = scissors shift in Ha(real)
  brod = broadening in Ha(real)
  tol = tolerance:how close to the singularity exact exact is calculated(real)
  fnam=root for filenames that will contain the output  :
   fnam1=trim(fnam)//'-ChiTotIm.out'
   fnam2=trim(fnam)//'-ChiTotRe.out'
   fnam3=trim(fnam)//'-ChiIm.out'
   fnam4=trim(fnam)//'-ChiRe.out'
   fnam5=trim(fnam)//'-ChiAbs.out'
  ncid=Netcdf id to save output data.

OUTPUT

  Calculates the second harmonic generation susceptibility on a desired energy mesh and
  for desired direction of polarisation. The output is in files named
  ChiTot.out : Im\chi_{v1v2v3}(2\omega,\omega,-\omega) and Re\chi_{v1v2v3}(2\omega,\omega,-\omega)
  ChiIm.out  : contributions to the Im\chi_{v1v2v3}(2\omega,\omega,-\omega) from various terms
  ChiRe.out  : contributions to Re\chi_{v1v2v3}(2\omega,\omega,-\omega) from various terms
  ChiAbs.out : abs\chi_{v1v2v3}(2\omega,\omega,-\omega). The headers in these files contain
  information about the calculation.

PARENTS

      optic

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

 837 subroutine nlinopt(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,evalv,efermi, &
 838   pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,ncid,comm)
 839 
 840 
 841 !This section has been created automatically by the script Abilint (TD).
 842 !Do not modify the following lines by hand.
 843 #undef ABI_FUNC
 844 #define ABI_FUNC 'nlinopt'
 845 !End of the abilint section
 846 
 847  implicit none
 848 
 849 !Arguments ------------------------------------
 850 !no_abirules
 851 integer, intent(in) :: icomp,itemp,nspin, ncid
 852 real(dp), intent(in) :: omega
 853 integer, intent(in) :: nkpt
 854 real(dp), intent(in) :: wkpt(nkpt)
 855 integer, intent(in) :: nsymcrys
 856 real(dp), intent(in) :: symcrys(3,3,nsymcrys)
 857 integer, intent(in) :: nstval
 858 real(dp), intent(in) :: evalv(nstval,nkpt,nspin)
 859 real(dp), intent(in) :: efermi
 860 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin)
 861 integer, intent(in) :: v1
 862 integer, intent(in) :: v2
 863 integer, intent(in) :: v3
 864 integer, intent(in) :: nmesh
 865 integer, intent(in) :: comm
 866 real(dp), intent(in) :: de
 867 real(dp), intent(in) :: sc
 868 real(dp), intent(in) :: brod
 869 real(dp), intent(in) :: tol
 870 character(len=*), intent(in) :: fnam
 871 
 872 !Local variables -------------------------
 873 integer :: iw
 874 integer :: i,j,k,lx,ly,lz
 875 integer :: isp,isym,ik
 876 integer :: ist1,ist2,istl,istn,istm
 877 integer,parameter :: master=0
 878 integer :: my_rank, nproc
 879 integer :: my_k1, my_k2
 880 integer :: ierr
 881 integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7
 882 real(dp) :: f1,f2,f3
 883 real(dp) :: ha2ev
 884 real(dp) :: t1,t2,t3,tst
 885 real(dp) :: ene,totre,totabs,totim
 886 real(dp) :: e1,e2,el,en,em
 887 real(dp) :: emin,emax,my_emin,my_emax
 888 real(dp) :: const_esu,const_au,au2esu
 889 real(dp) :: wmn,wnm,wln,wnl,wml,wlm
 890 complex(dpc) :: idel,w,zi
 891 complex(dpc) :: mat2w,mat1w1,mat1w2,mat2w_tra,mat1w3_tra
 892 complex(dpc) :: b111,b121,b131,b112,b122,b132,b113,b123,b133
 893 complex(dpc) :: b241,b242,b243,b221,b222,b223,b211,b212,b213,b231
 894 complex(dpc) :: b311,b312,b313,b331
 895 complex(dpc) :: b24,b21_22,b11,b12_13,b31_32
 896 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7
 897 character(500) :: msg
 898 ! local allocatable arrays
 899 integer :: start4(4),count4(4)
 900 real(dp) :: s(3,3),sym(3,3,3)
 901 complex(dpc), allocatable :: px(:,:,:,:,:)
 902 complex(dpc), allocatable :: py(:,:,:,:,:)
 903 complex(dpc), allocatable :: pz(:,:,:,:,:)
 904 complex(dpc), allocatable :: delta(:,:,:)
 905 complex(dpc), allocatable :: inter2w(:)
 906 complex(dpc), allocatable :: inter1w(:)
 907 complex(dpc), allocatable :: intra2w(:)
 908 complex(dpc), allocatable :: intra1w(:)
 909 complex(dpc), allocatable :: intra1wS(:),chi2tot(:)
 910 
 911 ! *********************************************************************
 912 
 913  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
 914 
 915 !calculate the constant
 916  zi=(0._dp,1._dp)
 917  idel=zi*brod
 918  const_au=-2._dp/(omega*dble(nsymcrys))
 919  au2esu=5.8300348177d-8
 920  const_esu=const_au*au2esu
 921  ha2ev=13.60569172*2._dp
 922 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 923 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
 924 !bohr: 5.2917ifc nlinopt.f907E-11
 925 !c: 2.99792458   velocity of sound
 926 !ry2ev: 13.60569172
 927 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2)
 928 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
 929 !mass comes from converting P_mn to r_mn
 930 !hbar^3 comes from converting all frequencies to energies in denominator
 931 !hbar^3 comes from operator for momentum (hbar/i nabla)
 932 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 933 !output file names
 934  fnam1=trim(fnam)//'-ChiTotIm.out'
 935  fnam2=trim(fnam)//'-ChiTotRe.out'
 936  fnam3=trim(fnam)//'-ChiIm.out'
 937  fnam4=trim(fnam)//'-ChiRe.out'
 938  fnam5=trim(fnam)//'-ChiAbs.out'
 939  fnam6=trim(fnam)//'-ChiImDec.out'
 940  fnam7=trim(fnam)//'-ChiReDec.out'
 941 
 942  if(my_rank == master) then
 943   !If there exists inversion symmetry exit with a message.
 944    tst=1.d-09
 945    do isym=1,nsymcrys
 946      t1=symcrys(1,1,isym)+1
 947      t2=symcrys(2,2,isym)+1
 948      t3=symcrys(3,3,isym)+1
 949   !  test if diagonal elements are -1
 950      if (abs(t1).lt.tst.and.abs(t2).lt.tst.and.abs(t3).lt.tst) then
 951   !    test if off-diagonal elements are zero
 952        if (abs(symcrys(1,2,isym)).lt.tst.and.abs(symcrys(1,3,isym)).lt.tst &
 953        .and.abs(symcrys(2,1,isym)).lt.tst.and.abs(symcrys(2,3,isym)).lt.tst.and.  &
 954        abs(symcrys(3,1,isym)).lt.tst.and.abs(symcrys(3,2,isym)).lt.tst) then
 955          write(std_out,*) '-----------------------------------------'
 956          write(std_out,*) '    the crystal has inversion symmetry   '
 957          write(std_out,*) '    the SHG susceptibility is zero       '
 958          write(std_out,*) '-----------------------------------------'
 959          MSG_ERROR("Aborting now")
 960        end if
 961      end if
 962    end do
 963   !check polarisation
 964    if (v1.le.0.or.v2.le.0.or.v3.le.0.or.v1.gt.3.or.v2.gt.3.or.v3.gt.3) then
 965      write(std_out,*) '---------------------------------------------'
 966      write(std_out,*) '    Error in nlinopt:                        '
 967      write(std_out,*) '    the polarisation directions incorrect    '
 968      write(std_out,*) '    1=x,  2=y  and 3=z                       '
 969      write(std_out,*) '---------------------------------------------'
 970      MSG_ERROR("Aborting now")
 971    end if
 972   !number of energy mesh points
 973    if (nmesh.le.0) then
 974      write(std_out,*) '---------------------------------------------'
 975      write(std_out,*) '    Error in nlinopt:                        '
 976      write(std_out,*) '    number of energy mesh points incorrect   '
 977      write(std_out,*) '    number has to be integer greater than 0  '
 978      write(std_out,*) '    nmesh*de = max energy for calculation    '
 979      write(std_out,*) '---------------------------------------------'
 980      MSG_ERROR("Aborting now")
 981    end if
 982   !step in energy
 983    if (de.le.0._dp) then
 984      write(std_out,*) '---------------------------------------------'
 985      write(std_out,*) '    Error in nlinopt:                        '
 986      write(std_out,*) '    energy step is incorrect                 '
 987      write(std_out,*) '    number has to real greater than 0.0      '
 988      write(std_out,*) '    nmesh*de = max energy for calculation    '
 989      write(std_out,*) '---------------------------------------------'
 990      MSG_ERROR("Aborting now")
 991    end if
 992   !broadening
 993    if (brod.gt.0.009) then
 994      write(std_out,*) '---------------------------------------------'
 995      write(std_out,*) '    ATTENTION: broadening is quite high      '
 996      write(std_out,*) '    ideally should be less than 0.005        '
 997      write(std_out,*) '---------------------------------------------'
 998    else if (brod.gt.0.015) then
 999      write(std_out,*) '----------------------------------------'
1000      write(std_out,*) '    ATTENTION: broadening is too high   '
1001      write(std_out,*) '    ideally should be less than 0.005   '
1002      write(std_out,*) '----------------------------------------'
1003    end if
1004   !tolerance
1005    if (tol.gt.0.006) then
1006      write(std_out,*) '----------------------------------------'
1007      write(std_out,*) '    ATTENTION: tolerance is too high    '
1008      write(std_out,*) '    ideally should be less than 0.004   '
1009      write(std_out,*) '----------------------------------------'
1010    end if
1011  end if
1012 
1013  !allocate local arrays
1014  ABI_ALLOCATE(px,(nstval,nstval,3,3,3))
1015  ABI_ALLOCATE(py,(nstval,nstval,3,3,3))
1016  ABI_ALLOCATE(pz,(nstval,nstval,3,3,3))
1017  ABI_ALLOCATE(inter2w,(nmesh))
1018  ABI_ALLOCATE(inter1w,(nmesh))
1019  ABI_ALLOCATE(intra2w,(nmesh))
1020  ABI_ALLOCATE(intra1w,(nmesh))
1021  ABI_ALLOCATE(intra1wS,(nmesh))
1022  ABI_ALLOCATE(delta,(nstval,nstval,3))
1023 
1024 !generate the symmetrizing tensor
1025  sym(:,:,:)=0._dp
1026  do isym=1,nsymcrys
1027    s(:,:)=symcrys(:,:,isym)
1028    do i=1,3
1029      do j=1,3
1030        do k=1,3
1031          sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
1032        end do
1033      end do
1034    end do
1035  end do
1036  !DBYG
1037  ! Disable symmetries for now
1038  !sym(:,:,:) = 0._dp
1039  !sym(v1,v2,v3) = nsymcrys
1040  !ENDDBYG
1041 
1042  ! Split work
1043  call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr)
1044 
1045 !initialise
1046  inter2w(:)=0._dp
1047  inter1w(:)=0._dp
1048  intra2w(:)=0._dp
1049  intra1w(:)=0._dp
1050  intra1wS(:)=0._dp
1051  delta(:,:,:)=0._dp
1052 
1053  my_emin=HUGE(0._dp)
1054  my_emax=-HUGE(0._dp)
1055 ! loop over kpts
1056  do ik=my_k1,my_k2
1057    write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt
1058 ! loop over spins
1059    do isp=1,nspin
1060 !  loop over states
1061      do ist1=1,nstval
1062        e1=evalv(ist1,ik,isp)
1063        if (e1.lt.efermi) then   ! ist1 is a valence state
1064          do ist2=1,nstval
1065            e2=evalv(ist2,ik,isp)
1066            if (e2.gt.efermi) then ! ist2 is a conduction state
1067 !            symmetrize the momentum matrix elements
1068              do lx=1,3
1069                do ly=1,3
1070                  do lz=1,3
1071                    f1=sym(lx,ly,lz)+sym(lx,lz,ly)
1072                    f2=sym(ly,lx,lz)+sym(ly,lz,lx)
1073                    f3=sym(lz,lx,ly)+sym(lz,ly,lx)
1074                    px(ist1,ist2,lx,ly,lz)=f1*pmat(ist1,ist2,ik,lx,isp)
1075                    py(ist2,ist1,lx,ly,lz)=f2*pmat(ist2,ist1,ik,lx,isp)
1076                    pz(ist2,ist1,lx,ly,lz)=f3*pmat(ist2,ist1,ik,lx,isp)
1077                  end do
1078                end do
1079              end do
1080 !            end loop over states
1081            end if
1082          end do
1083        end if
1084      end do
1085 !    calculate the energy window and \Delta_nm
1086      do ist1=1,nstval
1087        my_emin=min(my_emin,evalv(ist1,ik,isp))
1088        my_emax=max(my_emax,evalv(ist1,ik,isp))
1089        do ist2=1,nstval
1090          delta(ist1,ist2,1:3)=pmat(ist1,ist1,ik,1:3,isp)-pmat(ist2,ist2,ik,1:3,isp)
1091        end do
1092      end do
1093 !    initialise the factors
1094 !    factors are named according to the Ref. article 2.
1095      b111=0._dp
1096      b121=0._dp
1097      b131=0._dp
1098      b112=0._dp
1099      b122=0._dp
1100      b132=0._dp
1101      b113=0._dp
1102      b123=0._dp
1103      b133=0._dp
1104      b211=0._dp
1105      b221=0._dp
1106      b212=0._dp
1107      b222=0._dp
1108      b213=0._dp
1109      b223=0._dp
1110      b231=0._dp
1111      b241=0._dp
1112      b242=0._dp
1113      b243=0._dp
1114      b311=0._dp
1115      b312=0._dp
1116      b313=0._dp
1117      b331=0._dp
1118 !    start the calculation
1119      do istn=1,nstval
1120        en=evalv(istn,ik,isp)
1121        if (en.lt.efermi) then    ! istn is a valence state
1122          do istm=1,nstval
1123            em=evalv(istm,ik,isp)
1124            if (em.gt.efermi) then   ! istm is a conduction state
1125              em = em + sc ! Should add the scissor to conduction energies
1126              wmn=em-en
1127              wnm=-wmn
1128 !            calculate the matrix elements for two band intraband term
1129              mat2w_tra=0._dp
1130              mat1w3_tra=0._dp
1131              do lx=1,3
1132                do ly=1,3
1133                  do lz=1,3
1134                    mat2w_tra=mat2w_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,lz,isp)    &
1135                    *delta(istm,istn,ly)
1136                    mat1w3_tra=mat1w3_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,ly,isp)  &
1137                    *delta(istm,istn,lz)
1138 !                  NOTE:: lx to ly m to n in pmat matrices respectively
1139 !                  Changes are made so that this (b3) term is according to paper
1140 !                  [[cite:Sipe1993]] (Ref. 4) rather than [[cite:Hughes1996]] (Ref 2) in which this term is incorrect
1141                  end do
1142                end do
1143              end do
1144              b331=mat1w3_tra/wnm
1145              b11=0._dp
1146              b12_13=0._dp
1147              b24=0._dp
1148              b31_32=0._dp
1149              b21_22=0._dp
1150 
1151              b231=8._dp*mat2w_tra/wmn
1152              b331=mat1w3_tra/(wnm)
1153 !            !!!!!!!!!!!!!!!!!!!
1154 !            istl < istn   !
1155 !            !!!!!!!!!!!!!!!!!!!
1156              do istl=1,istn-1           ! istl is a valence state below istn
1157                el=evalv(istl,ik,isp)
1158                wln=el-en                ! do not add sc to the valence bands!
1159                wml=em-el
1160                wnl=-wln
1161                wlm=-wml
1162 !              calculate the matrix elements for three band terms
1163                mat2w=0._dp
1164                mat1w1=0._dp
1165                mat1w2=0._dp
1166                do lx=1,3
1167                  do ly=1,3
1168                    do lz=1,3
1169 
1170                      mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp)   &
1171                      *pmat(istl,istn,ik,lz,isp))
1172 
1173                      mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1174                      *pmat(istn,istl,ik,ly,isp))
1175 
1176                      mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1177                      *pmat(istn,istl,ik,ly,isp))
1178                    end do
1179                  end do
1180                end do
1181                b111=mat2w*(1._dp/(wln+wlm))*(1._dp/wlm)
1182                b121=mat1w1*(1._dp/(wnm+wlm))*(1._dp/wlm)
1183                b131=mat1w2*(1._dp/wlm)
1184 !
1185                b221=0._dp
1186                b211=mat1w1/wml
1187                b241=-mat2w/wml
1188 !
1189                b311=mat1w2/wlm
1190                if (abs(wln).gt.tol) then
1191                  b111=b111/wln
1192                  b121=b121/wln
1193                  b131=b131/wln
1194                  b221=mat1w2/wln
1195                  b241=b241+(mat2w/wln)
1196                  b311=b311+(mat1w1/wln)
1197                else
1198                  b111=0._dp
1199                  b121=0._dp
1200                  b131=0._dp
1201                  b221=0._dp
1202                end if
1203                t1=wln-wnm
1204                if (abs(t1).gt.tol) then
1205                  b131=b131/t1
1206                else
1207                  b131=0._dp
1208                end if
1209                b11=b11-2._dp*b111
1210                b12_13=b12_13+b121+b131
1211                b21_22=b21_22-b211+b221
1212                b24=b24+2._dp*b241
1213                b31_32=b31_32+b311
1214 !              end loop over istl
1215              end do
1216 
1217 !            !!!!!!!!!!!!!!!!!!!!!!!!!!!
1218 !            istn < istl < istm    !
1219 !            !!!!!!!!!!!!!!!!!!!!!!!!!!!
1220              do istl=istn+1,istm-1
1221                el=evalv(istl,ik,isp)
1222 !              calculate the matrix elements for three band terms
1223                mat2w=0._dp
1224                mat1w1=0._dp
1225                mat1w2=0._dp
1226                do lx=1,3
1227                  do ly=1,3
1228                    do lz=1,3
1229 
1230                      mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp)   &
1231                      *pmat(istl,istn,ik,lz,isp))
1232 
1233                      mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1234                      *pmat(istn,istl,ik,ly,isp))
1235 
1236                      mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1237                      *pmat(istn,istl,ik,ly,isp))
1238                    end do
1239                  end do
1240                end do
1241                if (el.lt.efermi) then
1242                  wln=el-en
1243                  wnl=-wln
1244                  wml=em-el
1245                  wlm=-wml
1246                else
1247                  el=el+sc
1248                  wln=el-en
1249                  wnl=-wln
1250                  wml=em-el
1251                  wlm=-wml
1252                end if
1253 !
1254                b112=0._dp
1255                b122=mat1w1*(1._dp/(wnm+wlm))
1256                b132=mat1w2*(1._dp/(wnm+wnl))
1257                b242=0._dp
1258                b222=0._dp
1259                b212=0._dp
1260                b312=0._dp
1261                if (abs(wnl).gt.tol) then
1262                  b112=mat2w/wln
1263                  b122=b122/wnl
1264                  b132=b132/wnl
1265                  b242=mat2w/wln
1266                  b222=mat1w2/wln
1267                  b312=mat1w1/wln
1268                else
1269                  b122=0._dp
1270                  b132=0._dp
1271                end if
1272                if (abs(wlm).gt.tol) then
1273                  b112=b112/wml
1274                  b122=b122/wlm
1275                  b132=b132/wlm
1276                  b242=b242-(mat2w/wml)
1277                  b212=mat1w1/wml
1278                  b312=b312+(mat1w2/wlm)
1279                else
1280                  b112=0._dp
1281                  b122=0._dp
1282                  b132=0._dp
1283                  b212=0._dp
1284                end if
1285                t1=wlm-wnl
1286                if (abs(t1).gt.tol) then
1287                  b112=b112/t1
1288                else
1289                  b112=0._dp
1290                end if
1291                b11=b11+2._dp*b112
1292                b12_13=b12_13-b122+b132
1293                b24=b24+2._dp*b242
1294                b21_22=b21_22-b212+b222
1295                b31_32=b31_32+b312
1296 !              end loop over istl
1297              end do
1298 
1299 !            !!!!!!!!!!!!!!!!!!!!!
1300 !            istl > istm    !
1301 !            !!!!!!!!!!!!!!!!!!!!!
1302              do istl=istm+1,nstval
1303                el=evalv(istl,ik,isp)+sc
1304                wln=el-en
1305                wnl=-wln
1306                wml=em-el
1307                wlm=-wml
1308 !              calculate the matrix elements for three band terms
1309                mat2w=0._dp
1310                mat1w1=0._dp
1311                mat1w2=0._dp
1312                do lx=1,3
1313                  do ly=1,3
1314                    do lz=1,3
1315                      mat2w=mat2w+px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) &
1316                      *pmat(istl,istn,ik,lz,isp)
1317 
1318                      mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1319                      *pmat(istn,istl,ik,ly,isp))
1320 
1321                      mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1322                      *pmat(istn,istl,ik,ly,isp))
1323                    end do
1324                  end do
1325                end do
1326 !
1327                b113=mat2w*(1._dp/(wnl+wml))*(1._dp/wnl)
1328                b123=mat1w1*(1._dp/wnl)
1329                b133=mat1w2*(1._dp/wnl)*(1._dp/(wnl+wnm))
1330                b243=mat2w/wln
1331                b223=mat1w2/wln
1332                b213=0._dp
1333                b313=-1._dp*mat1w1/wnl
1334                if (abs(wml).gt.tol) then
1335                  b113=b113/wml
1336                  b123=b123/wml
1337                  b133=b133/wml
1338                  b243=b243-(mat2w/wml)
1339                  b213=mat1w1/wml
1340                  b313=b313+(mat1w2/wlm)
1341                else
1342                  b113=0._dp
1343                  b123=0._dp
1344                  b133=0._dp
1345                end if
1346                t1=wnm-wml
1347                if (abs(t1).gt.tol) then
1348                  b123=b123/t1
1349                else
1350                  b123=0._dp
1351                end if
1352                b11=b11+2._dp*b113
1353                b12_13=b12_13+b123-b133
1354                b21_22=b21_22-b213+b223
1355                b24=b24+2._dp*b243
1356                b31_32=b31_32+b313
1357 !              end loop over istl
1358              end do
1359 
1360              b11=b11*zi*(1._dp/wnm)*const_esu
1361              b12_13=b12_13*zi*(1._dp/wnm)*const_esu
1362              b24=(b24+b231)*zi*(1._dp/(wnm**3))*const_esu
1363              b21_22=(b21_22)*zi*(1._dp/(wnm**3))*const_esu
1364              b31_32=(b31_32-b331)*zi*(1._dp/(wmn**3))*const_esu*0.5_dp
1365 !            calculate over the desired energy mesh and sum over k-points
1366              do iw=1,nmesh
1367                w=(iw-1)*de+idel
1368                inter2w(iw)=inter2w(iw)+(wkpt(ik)*(b11/(wmn-2._dp*w))) ! Inter(2w) from chi
1369                inter1w(iw)=inter1w(iw)+(wkpt(ik)*(b12_13/(wmn-w))) ! Inter(1w) from chi
1370                intra2w(iw)=intra2w(iw)+(wkpt(ik)*(b24/(wmn-2._dp*w))) ! Intra(2w) from eta
1371                intra1w(iw)=intra1w(iw)+(wkpt(ik)*((b21_22)/(wmn-w))) ! Intra(1w) from eta
1372                intra1wS(iw)=intra1wS(iw)+(wkpt(ik)*((b31_32)/(wmn-w))) ! Intra(1w) from sigma
1373              end do
1374            end if
1375          end do ! istn and istm
1376        end if
1377      end do
1378    end do  ! spins
1379  end do ! k-points
1380 
1381  call xmpi_sum(inter2w,comm,ierr)
1382  call xmpi_sum(inter1w,comm,ierr)
1383  call xmpi_sum(intra2w,comm,ierr)
1384  call xmpi_sum(intra1w,comm,ierr)
1385  call xmpi_sum(intra1wS,comm,ierr)
1386  call xmpi_min(my_emin,emin,comm,ierr)
1387  call xmpi_max(my_emax,emax,comm,ierr)
1388 
1389  if (my_rank == master) then
1390    ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
1391 
1392    if (ncid /= nctk_noid) then
1393      start4 = [1, 1, icomp, itemp]
1394      count4 = [2, nmesh, 1, 1]
1395      ABI_MALLOC(chi2tot, (nmesh))
1396      chi2tot = inter2w + inter1w + intra2w + intra1w + intra1wS
1397 #ifdef HAVE_NETCDF
1398      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter2w"), c2r(inter2w), start=start4, count=count4))
1399      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter1w"), c2r(inter1w), start=start4, count=count4))
1400      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra2w"), c2r(intra2w), start=start4, count=count4))
1401      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1w"), c2r(intra1w), start=start4, count=count4))
1402      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1wS"), c2r(intra1wS), start=start4, count=count4))
1403      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_chi2tot"), c2r(chi2tot), start=start4, count=count4))
1404 #endif
1405      ABI_FREE(chi2tot)
1406    end if
1407 
1408    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
1409      MSG_ERROR(msg)
1410    end if
1411    if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
1412      MSG_ERROR(msg)
1413    end if
1414    if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
1415      MSG_ERROR(msg)
1416    end if
1417    if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
1418      MSG_ERROR(msg)
1419    end if
1420    if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
1421      MSG_ERROR(msg)
1422    end if
1423    if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then
1424      MSG_ERROR(msg)
1425    end if
1426    if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then
1427      MSG_ERROR(msg)
1428    end if
1429 !  write headers
1430    write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1431    write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
1432    write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
1433    write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
1434    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1435    write(fout1, '(a)' )' # Energy      Tot-Im Chi(-2w,w,w)  Tot-Im Chi(-2w,w,w)'
1436    write(fout1, '(a)' )' # eV          *10^-7 esu        *10^-12 m/V SI units '
1437    write(fout1, '(a)' )' # '
1438 
1439    write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1440    write(fout2, '(a,es16.6)') ' #tolerance:',tol
1441    write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1442    write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1443    write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1444    write(fout2, '(a)')' # Energy      Tot-Re Chi(-2w,w,w)  Tot-Re Chi(-2w,w,w)'
1445    write(fout2, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
1446    write(fout2, '(a)')' # '
1447 
1448    write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1449    write(fout3, '(a,es16.6)') ' #tolerance:',tol
1450    write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1451    write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1452    write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1453    write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
1454    write(fout3, '(a)')' # in esu'
1455    write(fout3, '(a)')' # '
1456 
1457    write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1458    write(fout4, '(a,es16.6)') ' #tolerance:',tol
1459    write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1460    write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1461    write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1462    write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
1463    write(fout4, '(a)')' # in esu'
1464    write(fout4, '(a)')' # '
1465 
1466    write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1467    write(fout5, '(a,es16.6)') ' #tolerance:',tol
1468    write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1469    write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1470    write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1471    write(fout5, '(a)')' # Energy(eV)  |TotChi(-2w,w,w)|   |Tot Chi(-2w,w,w)|'
1472    write(fout5, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
1473    write(fout5, '(a)')' # '
1474 
1475    write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1476    write(fout6, '(a,es16.6)') ' #tolerance:',tol
1477    write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1478    write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1479    write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1480    write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1481    write(fout6, '(a)')' # in esu'
1482    write(fout6, '(a)')' # '
1483 
1484    write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1485    write(fout7, '(a,es16.6)') ' #tolerance:',tol
1486    write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1487    write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1488    write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1489    write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1490    write(fout7, '(a)')' # in esu'
1491    write(fout7, '(a)')' # '
1492 
1493    totim=0._dp
1494    totre=0._dp
1495    totabs=0._dp
1496    do iw=2,nmesh
1497      ene=(iw-1)*de
1498      ene=ene*ha2ev
1499 
1500      totim=aimag(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1501      write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1502      totim=0._dp
1503 
1504      totre=dble(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1505      write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1506      totre=0._dp
1507 
1508      write(fout3,'(f15.6,4es15.6)') ene,aimag(inter2w(iw))/1.d-7,      &
1509      aimag(inter1w(iw))/1.d-7,aimag(intra2w(iw))/1.d-7, aimag(intra1w(iw)+intra1wS(iw))/1.d-7
1510 
1511      write(fout4,'(f15.6,4es15.6)') ene,dble(inter2w(iw))/1.d-7,       &
1512      dble(inter1w(iw))/1.d-7,dble(intra2w(iw))/1.d-7,dble(intra1w(iw)+intra1wS(iw))/1.d-7
1513 
1514      totabs=abs(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1515      write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1516      totabs=0._dp
1517 
1518      write(fout6,'(f15.6,4es15.6)') ene,aimag(inter2w(iw)+inter1w(iw))/1.d-7,      &
1519      aimag(intra2w(iw)+intra1w(iw))/1.d-7,aimag(intra1wS(iw))/1.d-7
1520 
1521      write(fout7,'(f15.6,4es15.6)') ene,dble(inter2w(iw)+inter1w(iw))/1.d-7,       &
1522      dble(intra2w(iw)+intra1w(iw))/1.d-7,dble(intra1wS(iw))/1.d-7
1523    end do
1524 
1525    close(fout1)
1526    close(fout2)
1527    close(fout3)
1528    close(fout4)
1529    close(fout5)
1530    close(fout6)
1531    close(fout7)
1532 !  print information
1533    write(std_out,*) ' '
1534    write(std_out,*) 'information about calculation just performed:'
1535    write(std_out,*) ' '
1536    write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of second order susceptibility'
1537    write(std_out,*) 'tolerance:',tol
1538    if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
1539    write(std_out,*) 'broadening:',brod,'Hartree'
1540    if (brod.gt.0.009) then
1541      write(std_out,*) ' '
1542      write(std_out,*) 'ATTENTION: broadening is quite high'
1543      write(std_out,*) ' '
1544    else if (brod.gt.0.015) then
1545      write(std_out,*) ' '
1546      write(std_out,*) 'ATTENTION: broadening is too high'
1547      write(std_out,*) ' '
1548    end if
1549    write(std_out,*) 'scissors shift:',sc,'Hartree'
1550    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree'
1551  end if
1552 
1553  ! deallocate local arrays
1554  ABI_DEALLOCATE(px)
1555  ABI_DEALLOCATE(py)
1556  ABI_DEALLOCATE(pz)
1557  ABI_DEALLOCATE(inter2w)
1558  ABI_DEALLOCATE(inter1w)
1559  ABI_DEALLOCATE(intra2w)
1560  ABI_DEALLOCATE(intra1w)
1561  ABI_DEALLOCATE(intra1wS)
1562  ABI_DEALLOCATE(delta)
1563 
1564 end subroutine nlinopt

m_optic_tools/nonlinopt [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 nonlinopt

FUNCTION

 Compute optical frequency dependent linear electro-optic susceptibility for semiconductors

INPUTS

  icomp=Sequential index associated to computed tensor components (used for netcdf output)
  itemp=Temperature index (used for netcdf output)
  nspin = number of spins(integer)
  omega = crystal volume in au (real)
  nkpt  = total number of kpoints (integer)
  wkpt(nkpt) = weights of kpoints (real)
  nsymcrys = number of crystal symmetry operations(integer)
  symcrys(3,3,nsymcrys) = symmetry operations in cartisian coordinates(real)
  nstval = total number of valence states(integer)
  evalv(nstval,nspin,nkpt) = eigen value for each band in Ha(real)
  occv(nstval,nspin,nkpt) = occupation number
  efermi = Fermi energy in Ha(real)
  pmat(nstval,nstval,nkpt,3,nspin) = momentum matrix elements in cartesian coordinates(complex)
  v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z
  nmesh = desired number of energy mesh points(integer)
  de = desired step in energy(real); nmesh*de=maximum energy for plotting
  sc = scissors shift in Ha(real)
  brod = broadening in Ha(real)
  tol = tolerance:how close to the singularity exact exact is calculated(real)
  fnam=root for filenames that will contain the output  :
   fnam1=trim(fnam)//'-ChiTotIm.out'
   fnam2=trim(fnam)//'-ChiTotRe.out'
   fnam3=trim(fnam)//'-ChiIm.out'
   fnam4=trim(fnam)//'-ChiRe.out'
   fnam5=trim(fnam)//'-ChiAbs.out'

OUTPUT

  Calculates the second harmonic generation susceptibility on a desired energy mesh and
  for desired direction of polarisation. The output is in files named
  ChiEOTot.out : Im\chi_{v1v2v3}(\omega,\omega,0) and Re\chi_{v1v2v3}(\omega,\omega,0)
  ChiEOIm.out  : contributions to the Im\chi_{v1v2v3}(\omega,\omega,0) from various terms
  ChiEORe.out  : contributions to Re\chi_{v1v2v3}(\omega,\omega,-0) from various terms
  ChiEOAbs.out : abs\chi_{v1v2v3}(\omega,\omega,0). The headers in these files contain
  information about the calculation.
  ncid=Netcdf id to save output data.

 COMMENTS
    - The routine has been written using notations of Ref. 2
    - This routine does not symmetrize the tensor (up to now)
    - Sum over all the states and use occupation factors instead of looping only on resonant contributions

PARENTS

      optic

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

2208 subroutine nonlinopt(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,evalv,occv,efermi, &
2209   pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,do_antiresonant,ncid,comm)
2210 
2211 
2212 !This section has been created automatically by the script Abilint (TD).
2213 !Do not modify the following lines by hand.
2214 #undef ABI_FUNC
2215 #define ABI_FUNC 'nonlinopt'
2216 !End of the abilint section
2217 
2218  implicit none
2219 
2220 !Arguments ------------------------------------
2221 !no_abirules
2222 integer, intent(in) :: icomp,itemp,nspin, ncid
2223 real(dp), intent(in) :: omega
2224 integer, intent(in) :: nkpt
2225 real(dp), intent(in) :: wkpt(nkpt)
2226 integer, intent(in) :: nsymcrys
2227 real(dp), intent(in) :: symcrys(3,3,nsymcrys)
2228 integer, intent(in) :: nstval
2229 real(dp), intent(in) :: evalv(nstval,nkpt,nspin)
2230 real(dp), intent(in) :: occv(nstval,nkpt,nspin)
2231 real(dp), intent(in) :: efermi
2232 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin)
2233 integer, intent(in) :: v1
2234 integer, intent(in) :: v2
2235 integer, intent(in) :: v3
2236 integer, intent(in) :: nmesh
2237 integer, intent(in) :: comm
2238 real(dp), intent(in) :: de
2239 real(dp), intent(in) :: sc
2240 real(dp), intent(in) :: brod
2241 real(dp), intent(in) :: tol
2242 character(len=*), intent(in) :: fnam
2243 logical, intent(in) :: do_antiresonant
2244 
2245 !Local variables -------------------------
2246 integer :: iw
2247 integer :: i,j,k,lx,ly,lz
2248 integer :: isp,isym,ik
2249 integer :: ist1,istl,istn,istm
2250 real(dp) :: ha2ev
2251 real(dp) :: t1,t2,t3,tst
2252 real(dp) :: ene,totre,totabs,totim
2253 real(dp) :: el,en,em
2254 real(dp) :: emin,emax, my_emin,my_emax
2255 real(dp) :: const_esu,const_au,au2esu
2256 real(dp) :: wmn,wnm,wln,wnl,wml,wlm
2257 complex(dpc) :: idel,w,zi
2258 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7
2259 ! local allocatable arrays
2260  integer :: start4(4),count4(4)
2261  real(dp) :: s(3,3),sym(3,3,3)
2262  integer :: istp
2263  real(dp) :: ep, wmp, wpn
2264  real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included !
2265  real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, flm
2266  complex(dpc), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a}
2267  complex(dpc), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a}
2268  complex(dpc), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k)
2269  complex(dpc), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c
2270  complex(dpc), allocatable :: chiw(:), chi2w(:) ! \chi_{II}^{abc}(-\omega,\omega,0)
2271  complex(dpc), allocatable :: etaw(:), eta2w(:) ! \eta_{II}^{abc}(-\omega,\omega,0)
2272  complex(dpc), allocatable :: sigmaw(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0)
2273  complex(dpc) :: num1, num2, den1, den2, term1, term2
2274  complex(dpc) :: chi1, chi2_1, chi2_2
2275  complex(dpc), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega)
2276  complex(dpc), allocatable :: eta1(:) ! Second term that depends on the frequency ! (omega)
2277  complex(dpc), allocatable :: chi2tot(:)
2278  complex(dpc) :: eta1_1, eta1_2, eta2_1, eta2_2
2279  complex(dpc) :: sigma2_1, sigma1
2280  complex(dpc), allocatable :: symrmn(:,:,:) ! (m,l,n) = 1/2*(rml^b rln^c+rml^c rln^b)
2281  complex(dpc) :: symrmnl(3,3), symrlmn(3,3), symrmln(3,3)
2282 !Parallelism
2283  integer :: my_rank, nproc
2284  integer,parameter :: master=0
2285  integer :: ierr
2286  integer :: my_k1, my_k2
2287  character(500) :: msg
2288  integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7
2289 
2290 ! *********************************************************************
2291 
2292  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
2293 
2294 !calculate the constant
2295  zi=(0._dp,1._dp)
2296  idel=zi*brod
2297  const_au=-2._dp/(omega*dble(nsymcrys))
2298  !const_au=-2._dp/(omega)
2299  au2esu=5.8300348177d-8
2300  const_esu=const_au*au2esu
2301  ha2ev=13.60569172*2._dp
2302 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2303 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
2304 !bohr: 5.2917ifc nlinopt.f907E-11
2305 !c: 2.99792458   velocity of sound
2306 !ry2ev: 13.60569172
2307 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2)
2308 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
2309 !mass comes from converting P_mn to r_mn
2310 !hbar^3 comes from converting all frequencies to energies in denominator
2311 !hbar^3 comes from operator for momentum (hbar/i nabla)
2312 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2313 !output file names
2314  fnam1=trim(fnam)//'-ChiSHGTotIm.out'
2315  fnam2=trim(fnam)//'-ChiSHGTotRe.out'
2316  fnam3=trim(fnam)//'-ChiSHGIm.out'
2317  fnam4=trim(fnam)//'-ChiSHGRe.out'
2318  fnam5=trim(fnam)//'-ChiSHGAbs.out'
2319  fnam6=trim(fnam)//'-ChiSHGImDec.out'
2320  fnam7=trim(fnam)//'-ChiSHGReDec.out'
2321 
2322 !If there exists inversion symmetry exit with a mesg.
2323  tst=1.d-09
2324  do isym=1,nsymcrys
2325    t1=symcrys(1,1,isym)+1
2326    t2=symcrys(2,2,isym)+1
2327    t3=symcrys(3,3,isym)+1
2328 !  test if diagonal elements are -1
2329    if (abs(t1).lt.tst.and.abs(t2).lt.tst.and.abs(t3).lt.tst) then
2330 !    test if off-diagonal elements are zero
2331      if (abs(symcrys(1,2,isym)).lt.tst.and.abs(symcrys(1,3,isym)).lt.tst &
2332      .and.abs(symcrys(2,1,isym)).lt.tst.and.abs(symcrys(2,3,isym)).lt.tst.and.  &
2333      abs(symcrys(3,1,isym)).lt.tst.and.abs(symcrys(3,2,isym)).lt.tst) then
2334        write(std_out,*) '-----------------------------------------'
2335        write(std_out,*) '    the crystal has inversion symmetry   '
2336        write(std_out,*) '    the SHG susceptibility is zero       '
2337        write(std_out,*) '-----------------------------------------'
2338        MSG_ERROR("Aborting now")
2339      end if
2340    end if
2341  end do
2342 
2343 !check polarisation
2344  if (v1.le.0.or.v2.le.0.or.v3.le.0.or.v1.gt.3.or.v2.gt.3.or.v3.gt.3) then
2345    write(std_out,*) '---------------------------------------------'
2346    write(std_out,*) '    Error in linelop:                        '
2347    write(std_out,*) '    the polarisation directions incorrect    '
2348    write(std_out,*) '    1=x,  2=y  and 3=z                       '
2349    write(std_out,*) '---------------------------------------------'
2350    MSG_ERROR("Aborting now")
2351  end if
2352 
2353 !number of energy mesh points
2354  if (nmesh.le.0) then
2355    write(std_out,*) '---------------------------------------------'
2356    write(std_out,*) '    Error in linelop:                        '
2357    write(std_out,*) '    number of energy mesh points incorrect   '
2358    write(std_out,*) '    number has to be integer greater than 0  '
2359    write(std_out,*) '    nmesh*de = max energy for calculation    '
2360    write(std_out,*) '---------------------------------------------'
2361    MSG_ERROR("Aborting now")
2362  end if
2363 
2364 !step in energy
2365  if (de.le.0._dp) then
2366    write(std_out,*) '---------------------------------------------'
2367    write(std_out,*) '    Error in nlinopt:                        '
2368    write(std_out,*) '    energy step is incorrect                 '
2369    write(std_out,*) '    number has to real greater than 0.0      '
2370    write(std_out,*) '    nmesh*de = max energy for calculation    '
2371    write(std_out,*) '---------------------------------------------'
2372    MSG_ERROR("Aborting now")
2373  end if
2374 
2375 !broadening
2376  if (brod.gt.0.009) then
2377    write(std_out,*) '---------------------------------------------'
2378    write(std_out,*) '    ATTENTION: broadening is quite high      '
2379    write(std_out,*) '    ideally should be less than 0.005        '
2380    write(std_out,*) '---------------------------------------------'
2381  else if (brod.gt.0.015) then
2382    write(std_out,*) '----------------------------------------'
2383    write(std_out,*) '    ATTENTION: broadening is too high   '
2384    write(std_out,*) '    ideally should be less than 0.005   '
2385    write(std_out,*) '----------------------------------------'
2386  end if
2387 
2388 !tolerance
2389  if (tol.gt.0.006) then
2390    write(std_out,*) '----------------------------------------'
2391    write(std_out,*) '    ATTENTION: tolerance is too high    '
2392    write(std_out,*) '    ideally should be less than 0.004   '
2393    write(std_out,*) '----------------------------------------'
2394  end if
2395 
2396  ! allocate local arrays
2397  ABI_MALLOC(enk,(nstval))
2398  ABI_MALLOC(delta,(nstval,nstval,3))
2399  ABI_MALLOC(rmnbc,(nstval,nstval,3,3))
2400  ABI_MALLOC(roverw,(nstval,nstval,3,3))
2401  ABI_MALLOC(rmna,(nstval,nstval,3))
2402  ABI_MALLOC(chiw,(nmesh))
2403  ABI_MALLOC(etaw,(nmesh))
2404  ABI_MALLOC(chi2w,(nmesh))
2405  ABI_MALLOC(eta2w,(nmesh))
2406  ABI_MALLOC(sigmaw,(nmesh))
2407  ABI_MALLOC(chi2,(nmesh))
2408  ABI_MALLOC(eta1,(nmesh))
2409  ABI_MALLOC(symrmn,(nstval,nstval,nstval))
2410 
2411 !generate the symmetrizing tensor
2412  sym(:,:,:)=0._dp
2413  do isym=1,nsymcrys
2414    s(:,:)=symcrys(:,:,isym)
2415    do i=1,3
2416      do j=1,3
2417        do k=1,3
2418          sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
2419        end do
2420      end do
2421    end do
2422  end do
2423 
2424 
2425 !initialise
2426  delta(:,:,:)=0._dp
2427  rmnbc(:,:,:,:)=0._dp
2428  chiw(:)=0._dp
2429  chi2w(:)=0._dp
2430  chi2(:) = 0._dp
2431  etaw(:)=0._dp
2432  eta2w(:)=0._dp
2433  sigmaw(:)=0._dp
2434  my_emin=HUGE(0._dp)
2435  my_emax=-HUGE(0._dp)
2436 
2437  ! Split work
2438  call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr)
2439 
2440 ! loop over kpts
2441  do ik=my_k1,my_k2
2442    write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt
2443    do isp=1,nspin
2444      ! Calculate the scissor corrected energies and the energy window
2445      do ist1=1,nstval
2446        en = evalv(ist1,ik,isp)
2447        my_emin=min(my_emin,en)
2448        my_emax=max(my_emax,en)
2449        if(en > efermi) then
2450          en = en + sc
2451        end if
2452        enk(ist1) = en
2453      end do
2454 
2455 !    calculate \Delta_nm and r_mn^a
2456      do istn=1,nstval
2457        en = enk(istn)
2458        do istm=1,nstval
2459          em = enk(istm)
2460          wmn = em - en
2461          delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp)
2462          if(abs(wmn) < tol) then
2463            rmna(istm,istn,1:3) = 0._dp
2464          else
2465            rmna(istm,istn,1:3)=pmat(istm,istn,ik,1:3,isp)/wmn
2466          end if
2467        end do
2468      end do
2469 !    calculate \r^b_mn;c
2470      do istm=1,nstval
2471        em = enk(istm)
2472        do istn=1,nstval
2473          en = enk(istn)
2474          wmn = em - en
2475          if (abs(wmn) < tol) then ! Degenerate energies
2476            rmnbc(istm,istn,:,:) = 0.0
2477            roverw(istm,istn,:,:) = 0.0
2478            cycle
2479          end if
2480          do ly = 1,3
2481            do lz = 1,3
2482              num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly))
2483              den1 = wmn
2484              term1 = num1/den1
2485              term2 = 0._dp
2486              do istp=1,nstval
2487                ep = enk(istp)
2488                wmp = em - ep
2489                wpn = ep - en
2490                num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly))
2491                den2 = wmn
2492                term2 = term2 + (num2/den2)
2493              end do
2494              rmnbc(istm,istn,ly,lz) = -term1-(zi*term2)
2495              roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz)
2496            end do
2497          end do
2498        end do
2499      end do
2500 
2501 !    initialise the factors
2502 !    start the calculation
2503      do istn=1,nstval
2504        en=enk(istn)
2505        fn=occv(istn,ik,isp)
2506        if(do_antiresonant .and. en .ge. efermi) then
2507          cycle
2508        end if
2509        do istm=1,nstval
2510          em=enk(istm)
2511          if (do_antiresonant .and. em .le. efermi) then
2512            cycle
2513          end if
2514          wmn=em-en
2515          wnm=-wmn
2516          fm = occv(istm,ik,isp)
2517          fnm = fn - fm
2518          if(abs(wmn) > tol) then
2519            chi1 = 0._dp
2520            chi2(:) = 0._dp
2521            chi2_1 = 0._dp
2522            chi2_2 = 0._dp
2523            eta1(:) = 0._dp
2524            eta1_1 = 0._dp
2525            eta1_2 = 0._dp
2526            eta2_1 = 0._dp
2527            eta2_2 = 0._dp
2528            sigma1 = 0._dp
2529            sigma2_1 = 0._dp
2530            ! Three band terms
2531            do istl=1,nstval
2532              el=enk(istl)
2533              fl = occv(istl,ik,isp)
2534              wlm = el-em
2535              wln = el-en
2536              wnl = -wln
2537              wml = em-el
2538              fnl = fn-fl
2539              fml = fm-fl
2540              flm = -fml
2541              fln = -fnl
2542              do ly = 1,3
2543                do lz = 1,3
2544                  symrmnl(ly,lz) = 0.5_dp*(rmna(istm,istn,ly)*rmna(istn,istl,lz)+rmna(istm,istn,lz)*rmna(istn,istl,ly))
2545                  symrlmn(ly,lz) = 0.5_dp*(rmna(istl,istm,ly)*rmna(istm,istn,lz)+rmna(istl,istm,lz)*rmna(istm,istn,ly))
2546                  symrmln(ly,lz) = 0.5_dp*(rmna(istm,istl,ly)*rmna(istl,istn,lz)+rmna(istm,istl,lz)*rmna(istl,istn,ly))
2547                end do
2548              end do
2549 
2550              do lx = 1,3
2551                do ly = 1,3
2552                  do lz = 1,3
2553                    sigma1 = sigma1 + sym(lx,ly,lz)*(wnl*rmna(istl,istm,lx)*symrmnl(ly,lz)-wlm*rmna(istn,istl,lx)*symrlmn(ly,lz))
2554                    eta2_2 = eta2_2 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*symrmln(ly,lz)*(wml-wln)
2555                    if(abs(wln-wml) > tol) then
2556                      chi1 = chi1 + sym(lx,ly,lz)*(rmna(istn,istm,lx)*symrmln(ly,lz))/(wln-wml)
2557                    end if
2558                    eta1_1 = eta1_1 + sym(lx,ly,lz)*wln*rmna(istn,istl,lx)*symrlmn(ly,lz)
2559                    eta1_2 = eta1_2 - sym(lx,ly,lz)*wml*rmna(istl,istm,lx)*symrmnl(ly,lz)
2560                    if(abs(wnl-wmn) > tol) then
2561                      chi2_1 = chi2_1 - sym(lx,ly,lz)*(fnm*rmna(istl,istm,lx)*symrmnl(ly,lz)/(wnl-wmn))
2562                    end if
2563                    if(abs(wmn-wlm) > tol) then
2564                      chi2_2 = chi2_2 - sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*symrlmn(ly,lz)/(wmn-wlm))
2565                    end if
2566                  end do
2567                end do
2568              end do
2569            end do
2570 
2571            ! Two band terms
2572            eta2_1 = 0.0_dp
2573            sigma2_1 = 0.0_dp
2574            do lx = 1,3
2575              do ly = 1,3
2576                do lz = 1,3
2577                  eta2_1 = eta2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp &
2578 &                    *(delta(istm,istn,ly)*rmna(istm,istn,lz)+delta(istm,istn,lz)*rmna(istm,istn,ly))
2579                  ! Correct version (Sipe 1993)
2580                  sigma2_1 = sigma2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp &
2581 &                    *(rmna(istm,istn,ly)*delta(istn,istm,lz)+rmna(istm,istn,lz)*delta(istn,istm,ly))
2582 
2583                  ! Incorrect version (Hughes 1996)
2584                  !sigma2_1 = fnm*delta(istn,istm,v1)*0.5_dp*(rmna(istm,istn,v2)*rmna(istn,istm,v3)+rmna(istm,istn,v3)*rmna(istn,istm,v2))
2585                end do
2586              end do
2587            end do
2588 !
2589 !          calculate over the desired energy mesh and sum over k-points
2590            do iw=1,nmesh
2591              w=(iw-1)*de+idel
2592              chi2w(iw) = chi2w(iw) + zi*wkpt(ik)*((2.0_dp*fnm*chi1/(wmn-2.0_dp*w)))*const_esu ! Inter(2w) from chi
2593              chiw(iw) = chiw(iw) + zi*wkpt(ik)*((chi2_1+chi2_2)/(wmn-w))*const_esu ! Inter(w) from chi
2594              eta2w(iw) = eta2w(iw) + zi*wkpt(ik)*(8.0_dp*(eta2_1/((wmn**2)*(wmn-2.0_dp*w))) &
2595 &                 + 2.0_dp*eta2_2/((wmn**2)*(wmn-2.0_dp*w)))*const_esu ! Intra(2w) from eta
2596              etaw(iw) = etaw(iw) + zi*wkpt(ik)*((eta1_1 + eta1_2)*fnm/((wmn**2)*(wmn-w)))*const_esu ! Intra(w) from eta
2597              sigmaw(iw) = sigmaw(iw) + 0.5_dp*zi*wkpt(ik)*(fnm*sigma1/((wmn**2)*(wmn-w)) &
2598 &                 + (sigma2_1/((wmn**2)*(wmn-w))))*const_esu ! Intra(1w) from sigma
2599            end do
2600          end if
2601        end do ! end loop over istn and istm
2602      end do
2603    end do ! spins
2604  end do ! k-points
2605 
2606  ! Collect info among the nodes
2607  call xmpi_min(my_emin,emin,comm,ierr)
2608  call xmpi_max(my_emax,emax,comm,ierr)
2609 
2610  call xmpi_sum(chiw,comm,ierr)
2611  call xmpi_sum(etaw,comm,ierr)
2612  call xmpi_sum(chi2w,comm,ierr)
2613  call xmpi_sum(eta2w,comm,ierr)
2614  call xmpi_sum(sigmaw,comm,ierr)
2615 
2616  ! Master writes the output
2617  if (my_rank == master) then
2618 
2619    if (ncid /= nctk_noid) then
2620      start4 = [1, 1, icomp, itemp]
2621      count4 = [2, nmesh, 1, 1]
2622      ABI_MALLOC(chi2tot, (nmesh))
2623      chi2tot = chiw + chi2w + etaw + eta2w + sigmaw
2624 #ifdef HAVE_NETCDF
2625      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2tot"), c2r(chi2tot), start=start4, count=count4))
2626      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chiw"), c2r(chiw), start=start4, count=count4))
2627      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_etaw"), c2r(etaw), start=start4, count=count4))
2628      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2w"), c2r(chi2w), start=start4, count=count4))
2629      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_eta2w"), c2r(eta2w), start=start4, count=count4))
2630      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_sigmaw"), c2r(sigmaw), start=start4, count=count4))
2631 #endif
2632      ABI_FREE(chi2tot)
2633    end if
2634 
2635    ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
2636    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
2637      MSG_ERROR(msg)
2638    end if
2639    if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
2640      MSG_ERROR(msg)
2641    end if
2642    if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
2643      MSG_ERROR(msg)
2644    end if
2645    if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
2646      MSG_ERROR(msg)
2647    end if
2648    if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
2649      MSG_ERROR(msg)
2650    end if
2651    if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then
2652      MSG_ERROR(msg)
2653    end if
2654    if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then
2655      MSG_ERROR(msg)
2656    end if
2657   !!write headers
2658    write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2659    write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
2660    write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
2661    write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
2662    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2663    write(fout1, '(a)' )' # Energy      Tot-Im Chi(-w,w,0)  Tot-Im Chi(-w,w,0)'
2664    write(fout1, '(a)' )' # eV          *10^-7 esu        *10^-12 m/V SI units '
2665    write(fout1, '(a)' )' # '
2666 
2667    write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2668    write(fout2, '(a,es16.6)') ' #tolerance:',tol
2669    write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2670    write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2671    write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2672    write(fout2, '(a)')' # Energy      Tot-Re Chi(-w,w,0)  Tot-Re Chi(-w,w,0)'
2673    write(fout2, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
2674    write(fout2, '(a)')' # '
2675 
2676    write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2677    write(fout3, '(a,es16.6)') ' #tolerance:',tol
2678    write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2679    write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2680    write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2681    write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
2682    write(fout3, '(a)')' # in esu'
2683    write(fout3, '(a)')' # '
2684 
2685    write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2686    write(fout4, '(a,es16.6)') ' #tolerance:',tol
2687    write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2688    write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2689    write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2690    write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
2691    write(fout4, '(a)')' # in esu'
2692    write(fout4, '(a)')' # '
2693 
2694    write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2695    write(fout5, '(a,es16.6)') ' #tolerance:',tol
2696    write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2697    write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2698    write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2699    write(fout5, '(a)')' # Energy(eV)  |TotChi(-w,w,0)|   |Tot Chi(-w,w,0)|'
2700    write(fout5, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
2701    write(fout5, '(a)')' # '
2702 
2703    write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2704    write(fout6, '(a,es16.6)') ' #tolerance:',tol
2705    write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2706    write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2707    write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2708    write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2709    write(fout6, '(a)')' # in esu'
2710    write(fout6, '(a)')' # '
2711 
2712    write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2713    write(fout7, '(a,es16.6)') ' #tolerance:',tol
2714    write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2715    write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2716    write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2717    write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2718    write(fout7, '(a)')' # in esu'
2719    write(fout7, '(a)')' # '
2720 
2721    totim=0._dp
2722    totre=0._dp
2723    totabs=0._dp
2724    do iw=2,nmesh
2725      ene=(iw-1)*de
2726      ene=ene*ha2ev
2727 
2728      totim=aimag(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7
2729      write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2730      totim=0._dp
2731 
2732      totre=dble(chiw(iw)+chi2w(iw)+eta2w(iw)+etaw(iw)+sigmaw(iw))/1.d-7
2733      write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2734      totre=0._dp
2735 
2736      write(fout3,'(f15.6,4es15.6)') ene,aimag(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7,     &
2737      aimag(eta2w(iw))/1.d-7,aimag(etaw(iw))/1.d-7+aimag(sigmaw(iw))/1.d-7
2738 
2739      write(fout4,'(f15.6,4es15.6)') ene,dble(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7,       &
2740      dble(eta2w(iw))/1.d-7,dble(etaw(iw))/1.d-7+dble(sigmaw(iw))/1.d-7
2741 
2742      totabs=abs(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7
2743      write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2744      totabs=0._dp
2745 
2746      write(fout6,'(f15.6,4es15.6)') ene,aimag(chi2w(iw)+chiw(iw))/1.d-7,      &
2747      aimag(eta2w(iw)+etaw(iw))/1.d-7,aimag(sigmaw(iw))/1.d-7
2748 
2749      write(fout7,'(f15.6,4es15.6)') ene,dble(chi2w(iw)+chiw(iw))/1.d-7,       &
2750      dble(eta2w(iw)+etaw(iw))/1.d-7,dble(sigmaw(iw))/1.d-7
2751    end do
2752 
2753    close(fout1)
2754    close(fout2)
2755    close(fout3)
2756    close(fout4)
2757    close(fout5)
2758    close(fout6)
2759    close(fout7)
2760 
2761    ! print information
2762    write(std_out,*) ' '
2763    write(std_out,*) 'information about calculation just performed:'
2764    write(std_out,*) ' '
2765    write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of second order susceptibility'
2766    write(std_out,*) 'tolerance:',tol
2767    if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
2768    write(std_out,*) 'broadening:',brod,'Hartree'
2769    if (brod.gt.0.009) then
2770      write(std_out,*) ' '
2771      write(std_out,*) 'ATTENTION: broadening is quite high'
2772      write(std_out,*) ' '
2773    else if (brod.gt.0.015) then
2774      write(std_out,*) ' '
2775      write(std_out,*) 'ATTENTION: broadening is too high'
2776      write(std_out,*) ' '
2777    end if
2778    write(std_out,*) 'scissors shift:',sc,'Hartree'
2779    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree'
2780 
2781  end if
2782 
2783  ! deallocate local arrays
2784  ABI_FREE(enk)
2785  ABI_FREE(delta)
2786  ABI_FREE(rmnbc)
2787  ABI_FREE(roverw)
2788  ABI_FREE(rmna)
2789  ABI_FREE(chiw)
2790  ABI_FREE(chi2w)
2791  ABI_FREE(chi2)
2792  ABI_FREE(etaw)
2793  ABI_FREE(eta1)
2794  ABI_FREE(symrmn)
2795  ABI_FREE(eta2w)
2796  ABI_FREE(sigmaw)
2797 
2798 end subroutine nonlinopt

m_optic_tools/pmat2cart [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 pmat2cart

FUNCTION

  turn momentum matrix elements to cartesian axes. To be used in optic calculation of linear and non-linear RPA dielectric matrices

INPUTS

  eigen11,eigen12,eigen13 = first order ddk eigen values = d eig_i,k / dk for 3 reduced directions
  mband=maximum number of bands
  nkpt = number of k-points
  nsppol=1 for unpolarized, 2 for spin-polarized
  rprimd(3,3)=dimensional real space primitive translations (bohr)

OUTPUT

  pmat(mband,mband,nkpt,3,nsppol) = matrix elements of momentum operator, in cartesian coordinates

PARENTS

      optic

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

274 subroutine pmat2cart(eigen11,eigen12,eigen13,mband,nkpt,nsppol,pmat,rprimd)
275 
276 
277 !This section has been created automatically by the script Abilint (TD).
278 !Do not modify the following lines by hand.
279 #undef ABI_FUNC
280 #define ABI_FUNC 'pmat2cart'
281 !End of the abilint section
282 
283  implicit none
284 
285 !Arguments -----------------------------------------------
286 !scalars
287  integer,intent(in) :: mband,nkpt,nsppol
288 !arrays
289  real(dp),intent(in) :: eigen11(2,mband,mband,nkpt,nsppol)
290  real(dp),intent(in) :: eigen12(2,mband,mband,nkpt,nsppol)
291  real(dp),intent(in) :: eigen13(2,mband,mband,nkpt,nsppol),rprimd(3,3)
292 !no_abirules
293  complex(dpc),intent(out) :: pmat(mband,mband,nkpt,3,nsppol)
294 
295 !Local variables -----------------------------------------
296 !scalars
297  integer :: iband1,iband2,ikpt,isppol
298 !arrays
299  real(dp) :: rprim(3,3)
300 
301 ! *************************************************************************
302 
303 !rescale the rprim
304  rprim(:,:) = rprimd(:,:) / two_pi
305 
306  do isppol=1,nsppol
307    do ikpt=1,nkpt
308      do iband1=1,mband
309        do iband2=1,mband
310          pmat(iband2,iband1,ikpt,:,isppol) =             &
311 &         rprim(:,1)*cmplx(eigen11(1,iband2,iband1,ikpt,isppol),eigen11(2,iband2,iband1,ikpt,isppol),kind=dp) &
312 &         +rprim(:,2)*cmplx(eigen12(1,iband2,iband1,ikpt,isppol),eigen12(2,iband2,iband1,ikpt,isppol),kind=dp) &
313 &         +rprim(:,3)*cmplx(eigen13(1,iband2,iband1,ikpt,isppol),eigen13(2,iband2,iband1,ikpt,isppol),kind=dp)
314        end do
315      end do
316    end do
317  end do
318 
319 end subroutine pmat2cart

m_optic_tools/pmat_renorm [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 pmat_renorm

FUNCTION

 Renormalize the momentum matrix elements according to the scissor shift which is imposed

INPUTS

  mband= number of bands
  nkpt = number of k-points
  nsppol=1 for unpolarized, 2 for spin-polarized
  efermi = Fermi level
  sc = scissor shift for conduction bands
  evalv = eigenvalues for ground state

OUTPUT

  pmat(mband,mband,nkpt,3,nsppol) = momentum matrix elements, renormalized by denominator change with scissor shift

PARENTS

      optic

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

350 subroutine pmat_renorm(efermi, evalv, mband, nkpt, nsppol, pmat, sc)
351 
352 
353 !This section has been created automatically by the script Abilint (TD).
354 !Do not modify the following lines by hand.
355 #undef ABI_FUNC
356 #define ABI_FUNC 'pmat_renorm'
357 !End of the abilint section
358 
359  implicit none
360 
361 !Arguments -----------------------------------------------
362 !scalars
363  integer, intent(in) :: nsppol
364  integer, intent(in) :: nkpt
365  integer, intent(in) :: mband
366  real(dp), intent(in) :: efermi
367  real(dp), intent(in) :: sc
368 !arrays
369  real(dp), intent(in) :: evalv(mband,nkpt,nsppol)
370  complex(dpc), intent(inout) :: pmat(mband,mband,nkpt,3,nsppol)
371 
372 !Local variables -----------------------------------------
373 !scalars
374  integer :: iband1,iband2,ikpt,isppol
375  real(dp) :: corec, e1, e2
376 
377 ! *************************************************************************
378 
379  if (abs(sc) < tol8) then
380    call wrtout(std_out,' No scissor shift to be applied. Returning to main optic routine.',"COLL")
381    return
382  end if
383 
384  do isppol=1,nsppol
385    do ikpt=1,nkpt
386      do iband1=1,mband ! valence states
387        e1 = evalv(iband1,ikpt,isppol)
388        if (e1 > efermi) cycle
389        do iband2=1,mband ! conduction states
390          e2 = evalv(iband2,ikpt,isppol)
391          if (e2 < efermi) cycle
392          corec = (e2+sc-e1)/(e2-e1)
393          pmat(iband2,iband1,ikpt,:,isppol) = corec * pmat(iband2,iband1,ikpt,:,isppol)
394          pmat(iband1,iband2,ikpt,:,isppol) = corec * pmat(iband1,iband2,ikpt,:,isppol)
395        end do
396      end do
397    end do
398  end do
399 
400 end subroutine pmat_renorm

m_optic_tools/sym2cart [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 sym2cart

FUNCTION

 Routine called by the program optic
 Convert to symmetry matrice in cartesian coordinates

INPUTS

      gprimd(3,3)=dimensional primitive translations for reciprocal space
      nsym=number of symmetries in group
      rprimd(3,3)=dimensional real space primitive translations (bohr)
      symrel(3,3,nsym)=symmetry matrices in terms of real space

OUTPUT

      symcart(3,3)=symmetry matrice in cartesian coordinates (reals)

PARENTS

      optic

CHILDREN

      xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum

SOURCE

 91 subroutine sym2cart(gprimd,nsym,rprimd,symrel,symcart)
 92 
 93 
 94 !This section has been created automatically by the script Abilint (TD).
 95 !Do not modify the following lines by hand.
 96 #undef ABI_FUNC
 97 #define ABI_FUNC 'sym2cart'
 98 !End of the abilint section
 99 
100  implicit none
101 
102 !Arguments -----------------------------------------------
103 ! in
104 ! out
105 !scalars
106  integer,intent(in) :: nsym
107 !arrays
108  integer,intent(in) :: symrel(3,3,nsym)
109  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
110  real(dp),intent(out) :: symcart(3,3,nsym)
111 
112 !Local variables-------------------------------
113 !scalars
114  integer :: isym
115 !arrays
116  real(dp) :: rsym(3,3),rsymcart(3,3),tmp(3,3)
117 
118 ! *************************************************************************
119 
120  do isym=1,nsym
121    rsym(:,:) = dble(symrel(:,:,isym))
122 !  write(std_out,*) 'rsym = ',rsym
123    call dgemm('N','N',3,3,3,one,rprimd,3,rsym,  3,zero,tmp,     3)
124    call dgemm('N','N',3,3,3,one,tmp,   3,gprimd,3,zero,rsymcart,3)
125 !  write(std_out,*) 'rsymcart = ',rsymcart
126    symcart(:,:,isym) = rsymcart(:,:)
127 ! purify symops in cartesian dp coordinates
128    where( abs(symcart(:,:,isym))<tol14)
129      symcart(:,:,isym) = zero
130    end where
131  end do
132 
133 end subroutine sym2cart