TABLE OF CONTENTS
- ABINIT/m_optic_tools
- m_optic_tools/getwtk
- m_optic_tools/linelop
- m_optic_tools/linopt
- m_optic_tools/nlinopt
- m_optic_tools/nonlinopt
- m_optic_tools/pmat2cart
- m_optic_tools/pmat_renorm
- m_optic_tools/sym2cart
ABINIT/m_optic_tools [ 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