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-2021 ABINIT group (SSharma,MVer,VRecoules,TD,YG, NAP)
 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

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

150 subroutine getwtk(kpt,nkpt,nsym,symrel,wtk)
151 
152 !Arguments -----------------------------------------------
153 !scalars
154  integer,intent(in) :: nkpt,nsym
155 !arrays
156  integer,intent(in) :: symrel(3,3,nsym)
157  real(dp),intent(in) :: kpt(3,nkpt)
158  real(dp),intent(out) :: wtk(nkpt)
159 
160 !Local variables -----------------------------------------
161 !scalars
162  integer :: ikpt,istar,isym,itim,new,nkpt_tot
163  real(dp) :: shift,timsign,tmp
164 !arrays
165  integer :: nstar(nkpt)
166  real(dp) :: dkpt(3),kptstar(3,2*nkpt*nsym),rsymrel(3,3,nsym),symkpt(3)
167  real(dp) :: tsymkpt(3)
168 
169 ! *************************************************************************
170 
171  do isym=1,nsym
172    rsymrel(:,:,isym) = dble(symrel(:,:,isym))
173  end do
174 
175 !for each kpt find star and accumulate nkpts
176  do ikpt=1,nkpt
177    write(std_out,*) ' getwtk : ikpt = ', ikpt
178    nstar(ikpt) = 0
179    kptstar(:,:) = zero
180    do isym=1,nsym
181 
182      call dgemv('N',3,3,one,rsymrel(:,:,isym),3,kpt(:,ikpt),1,zero,symkpt,1)
183 
184 !    is symkpt already in star?
185      do itim=0,1
186        timsign=one-itim*two
187        tsymkpt(:) = timsign*symkpt(:)
188        call wrap2_pmhalf(tsymkpt(1),tmp,shift) ;  tsymkpt(1) = tmp
189        call wrap2_pmhalf(tsymkpt(2),tmp,shift) ;  tsymkpt(2) = tmp
190        call wrap2_pmhalf(tsymkpt(3),tmp,shift) ;  tsymkpt(3) = tmp
191        new=1
192        do istar=1,nstar(ikpt)
193          dkpt(:) = abs(tsymkpt(:)-kptstar(:,istar))
194          if ( sum(dkpt) < 1.0d-6) then
195            new=0
196            exit
197          end if
198        end do
199        if (new==1) then
200          nstar(ikpt) = nstar(ikpt)+1
201          kptstar(:,nstar(ikpt)) = tsymkpt(:)
202        end if
203      end do
204 
205    end do
206 !  end do nsym
207 !  DEBUG
208 !  write(std_out,*) ' getwtk : nstar = ', nstar(ikpt)
209 !  write(std_out,*) ' getwtk : star = '
210 !  write(std_out,*)  kptstar(:,1:nstar(ikpt))
211 !  ENDDEBUG
212  end do
213 !end do nkpt
214 
215  nkpt_tot = sum(nstar)
216 !write(std_out,*) ' getwtk : nkpt_tot = ', nkpt_tot
217  do ikpt=1,nkpt
218    wtk(ikpt) = dble(nstar(ikpt))/dble(nkpt_tot)
219  end do
220 
221 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

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

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).

  If 'prtlincompmatrixelements' is set to 1, the matrix elements and other quantities used to build
  the chi tensor are stored in the _OPTIC.nc file as well. This includes the matrix elements,
  the occupations, the renormalized but unshifted eigenvalues and the kpts weights.

  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

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

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

m_optic_tools/nonlinopt [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 nonlinopt

FUNCTION

 Compute the frequency dependent nonlinear 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 nonlinear electro-optical 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

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

252 subroutine pmat2cart(eigen11,eigen12,eigen13,mband,nkpt,nsppol,pmat,rprimd)
253 
254 !Arguments -----------------------------------------------
255 !scalars
256  integer,intent(in) :: mband,nkpt,nsppol
257 !arrays
258  real(dp),intent(in) :: eigen11(2,mband,mband,nkpt,nsppol)
259  real(dp),intent(in) :: eigen12(2,mband,mband,nkpt,nsppol)
260  real(dp),intent(in) :: eigen13(2,mband,mband,nkpt,nsppol),rprimd(3,3)
261 !no_abirules
262  complex(dpc),intent(out) :: pmat(mband,mband,nkpt,3,nsppol)
263 
264 !Local variables -----------------------------------------
265 !scalars
266  integer :: iband1,iband2,ikpt,isppol
267 !arrays
268  real(dp) :: rprim(3,3)
269 
270 ! *************************************************************************
271 
272 !rescale the rprim
273  rprim(:,:) = rprimd(:,:) / two_pi
274 
275  do isppol=1,nsppol
276    do ikpt=1,nkpt
277      do iband1=1,mband
278        do iband2=1,mband
279          pmat(iband2,iband1,ikpt,:,isppol) =             &
280 &         rprim(:,1)*cmplx(eigen11(1,iband2,iband1,ikpt,isppol),eigen11(2,iband2,iband1,ikpt,isppol),kind=dp) &
281 &         +rprim(:,2)*cmplx(eigen12(1,iband2,iband1,ikpt,isppol),eigen12(2,iband2,iband1,ikpt,isppol),kind=dp) &
282 &         +rprim(:,3)*cmplx(eigen13(1,iband2,iband1,ikpt,isppol),eigen13(2,iband2,iband1,ikpt,isppol),kind=dp)
283        end do
284      end do
285    end do
286  end do
287 
288 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

319 subroutine pmat_renorm(efermi, evalv, mband, nkpt, nsppol, pmat, sc)
320 
321 !Arguments -----------------------------------------------
322 !scalars
323  integer, intent(in) :: nsppol
324  integer, intent(in) :: nkpt
325  integer, intent(in) :: mband
326  real(dp), intent(in) :: efermi
327  real(dp), intent(in) :: sc
328 !arrays
329  real(dp), intent(in) :: evalv(mband,nkpt,nsppol)
330  complex(dpc), intent(inout) :: pmat(mband,mband,nkpt,3,nsppol)
331 
332 !Local variables -----------------------------------------
333 !scalars
334  integer :: iband1,iband2,ikpt,isppol
335  real(dp) :: corec, e1, e2
336 
337 ! *************************************************************************
338 
339  if (abs(sc) < tol8) then
340    call wrtout(std_out,' No scissor shift to be applied. Returning to main optic routine.',"COLL")
341    return
342  end if
343 
344  do isppol=1,nsppol
345    do ikpt=1,nkpt
346      do iband1=1,mband ! valence states
347        e1 = evalv(iband1,ikpt,isppol)
348        if (e1 > efermi) cycle
349        do iband2=1,mband ! conduction states
350          e2 = evalv(iband2,ikpt,isppol)
351          if (e2 < efermi) cycle
352          corec = (e2+sc-e1)/(e2-e1)
353          pmat(iband2,iband1,ikpt,:,isppol) = corec * pmat(iband2,iband1,ikpt,:,isppol)
354          pmat(iband1,iband2,ikpt,:,isppol) = corec * pmat(iband1,iband2,ikpt,:,isppol)
355        end do
356      end do
357    end do
358  end do
359 
360 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

 90 subroutine sym2cart(gprimd,nsym,rprimd,symrel,symcart)
 91 
 92 !Arguments -----------------------------------------------
 93 !scalars
 94  integer,intent(in) :: nsym
 95 !arrays
 96  integer,intent(in) :: symrel(3,3,nsym)
 97  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
 98  real(dp),intent(out) :: symcart(3,3,nsym)
 99 
100 !Local variables-------------------------------
101 !scalars
102  integer :: isym
103 !arrays
104  real(dp) :: rsym(3,3),rsymcart(3,3),tmp(3,3)
105 
106 ! *************************************************************************
107 
108  do isym=1,nsym
109    rsym(:,:) = dble(symrel(:,:,isym))
110 !  write(std_out,*) 'rsym = ',rsym
111    call dgemm('N','N',3,3,3,one,rprimd,3,rsym,  3,zero,tmp,     3)
112    call dgemm('N','N',3,3,3,one,tmp,   3,gprimd,3,zero,rsymcart,3)
113 !  write(std_out,*) 'rsymcart = ',rsymcart
114    symcart(:,:,isym) = rsymcart(:,:)
115 ! purify symops in cartesian dp coordinates
116    where( abs(symcart(:,:,isym))<tol14)
117      symcart(:,:,isym) = zero
118    end where
119  end do
120 
121 end subroutine sym2cart