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-2018 ABINIT group (SSharma,MVer,VRecoules,TD,YG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt . COMMENTS Right now the routine sums over the k-points. In future linear tetrahedron method might be useful. Reference articles: 1. S. Sharma, J. K. Dewhurst and C. Ambrosch-Draxl, Phys. Rev. B {\bf 67} 165332 2003 2. J. L. P. Hughes and J. E. Sipe, Phys. Rev. B {\bf 53} 10 751 1996 3. S. Sharma and C. Ambrosch-Draxl, Physica Scripta T 109 2004 4. J. E. Sipe and Ed. Ghahramani, Phys. Rev. B {\bf 48} 11 705 1993
SOURCE
27 #if defined HAVE_CONFIG_H 28 #include "config.h" 29 #endif 30 31 #include "abi_common.h" 32 33 MODULE m_optic_tools 34 35 use defs_basis 36 use m_errors 37 use m_profiling_abi 38 use m_linalg_interfaces 39 use m_xmpi 40 use m_nctk 41 #ifdef HAVE_NETCDF 42 use netcdf 43 #endif 44 45 use defs_datatypes, only : ebands_t 46 use m_numeric_tools, only : wrap2_pmhalf, c2r 47 use m_io_tools, only : open_file 48 49 implicit none 50 51 private 52 53 public :: sym2cart 54 public :: getwtk 55 public :: pmat2cart 56 public :: pmat_renorm 57 public :: linopt ! Compute dielectric function for semiconductors 58 public :: nlinopt ! Second harmonic generation susceptibility for semiconductors 59 public :: linelop ! Linear electro-optic susceptibility for semiconductors 60 public :: nonlinopt 61 62 CONTAINS !===========================================================
m_optic_tools/getwtk [ Functions ]
[ Top ] [ m_optic_tools ] [ Functions ]
NAME
getwtk
FUNCTION
Routine called by the program optic Presumes kpts are the irreducible ones of a good uniform grid
INPUTS
kpt(3,nkpt)=reduced coordinates of k points. nkpt = number of k points nsym=Number of symmetry operations. symrel(3,3,nsym)=symmetry operations
OUTPUT
wtk(nkpt)=weight assigned to each k point.
PARENTS
CHILDREN
xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum
SOURCE
162 subroutine getwtk(kpt,nkpt,nsym,symrel,wtk) 163 164 165 !This section has been created automatically by the script Abilint (TD). 166 !Do not modify the following lines by hand. 167 #undef ABI_FUNC 168 #define ABI_FUNC 'getwtk' 169 !End of the abilint section 170 171 implicit none 172 173 !Arguments ----------------------------------------------- 174 ! in 175 ! out 176 !scalars 177 integer,intent(in) :: nkpt,nsym 178 !arrays 179 integer,intent(in) :: symrel(3,3,nsym) 180 real(dp),intent(in) :: kpt(3,nkpt) 181 real(dp),intent(out) :: wtk(nkpt) 182 183 !Local variables ----------------------------------------- 184 !scalars 185 integer :: ikpt,istar,isym,itim,new,nkpt_tot 186 real(dp) :: shift,timsign,tmp 187 !arrays 188 integer :: nstar(nkpt) 189 real(dp) :: dkpt(3),kptstar(3,2*nkpt*nsym),rsymrel(3,3,nsym),symkpt(3) 190 real(dp) :: tsymkpt(3) 191 192 ! ************************************************************************* 193 194 do isym=1,nsym 195 rsymrel(:,:,isym) = dble(symrel(:,:,isym)) 196 end do 197 198 !for each kpt find star and accumulate nkpts 199 do ikpt=1,nkpt 200 write(std_out,*) ' getwtk : ikpt = ', ikpt 201 nstar(ikpt) = 0 202 kptstar(:,:) = zero 203 do isym=1,nsym 204 205 call dgemv('N',3,3,one,rsymrel(:,:,isym),3,kpt(:,ikpt),1,zero,symkpt,1) 206 207 ! is symkpt already in star? 208 do itim=0,1 209 timsign=one-itim*two 210 tsymkpt(:) = timsign*symkpt(:) 211 call wrap2_pmhalf(tsymkpt(1),tmp,shift) ; tsymkpt(1) = tmp 212 call wrap2_pmhalf(tsymkpt(2),tmp,shift) ; tsymkpt(2) = tmp 213 call wrap2_pmhalf(tsymkpt(3),tmp,shift) ; tsymkpt(3) = tmp 214 new=1 215 do istar=1,nstar(ikpt) 216 dkpt(:) = abs(tsymkpt(:)-kptstar(:,istar)) 217 if ( sum(dkpt) < 1.0d-6) then 218 new=0 219 exit 220 end if 221 end do 222 if (new==1) then 223 nstar(ikpt) = nstar(ikpt)+1 224 kptstar(:,nstar(ikpt)) = tsymkpt(:) 225 end if 226 end do 227 228 end do 229 ! end do nsym 230 ! DEBUG 231 ! write(std_out,*) ' getwtk : nstar = ', nstar(ikpt) 232 ! write(std_out,*) ' getwtk : star = ' 233 ! write(std_out,*) kptstar(:,1:nstar(ikpt)) 234 ! ENDDEBUG 235 end do 236 !end do nkpt 237 238 nkpt_tot = sum(nstar) 239 !write(std_out,*) ' getwtk : nkpt_tot = ', nkpt_tot 240 do ikpt=1,nkpt 241 wtk(ikpt) = dble(nstar(ikpt))/dble(nkpt_tot) 242 end do 243 244 end subroutine getwtk
m_optic_tools/linelop [ Functions ]
[ Top ] [ m_optic_tools ] [ Functions ]
NAME
linelop
FUNCTION
Compute optical frequency dependent linear electro-optic susceptibility for semiconductors
INPUTS
icomp=Sequential index associated to computed tensor components (used for netcdf output) itemp=Temperature index (used for netcdf output) nspin = number of spins(integer) omega = crystal volume in au (real) nkpt = total number of kpoints (integer) wkpt(nkpt) = weights of kpoints (real) nsymcrys = number of crystal symmetry operations(integer) symcrys(3,3,nsymcrys) = symmetry operations in cartisian coordinates(real) nstval = total number of valence states(integer) evalv(nstval,nspin,nkpt) = eigen value for each band in Ha(real) occv(nstval,nspin,nkpt) = occupation number efermi = Fermi energy in Ha(real) pmat(nstval,nstval,nkpt,3,nspin) = momentum matrix elements in cartesian coordinates(complex) v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z nmesh = desired number of energy mesh points(integer) de = desired step in energy(real); nmesh*de=maximum energy for plotting sc = scissors shift in Ha(real) brod = broadening in Ha(real) tol = tolerance:how close to the singularity exact exact is calculated(real) fnam=root for filenames that will contain the output : fnam1=trim(fnam)//'-ChiTotIm.out' fnam2=trim(fnam)//'-ChiTotRe.out' fnam3=trim(fnam)//'-ChiIm.out' fnam4=trim(fnam)//'-ChiRe.out' fnam5=trim(fnam)//'-ChiAbs.out' ncid=Netcdf id to save output data.
OUTPUT
Calculates the second harmonic generation susceptibility on a desired energy mesh and for desired direction of polarisation. The output is in files named ChiEOTot.out : Im\chi_{v1v2v3}(\omega,\omega,0) and Re\chi_{v1v2v3}(\omega,\omega,0) ChiEOIm.out : contributions to the Im\chi_{v1v2v3}(\omega,\omega,0) from various terms ChiEORe.out : contributions to Re\chi_{v1v2v3}(\omega,\omega,-0) from various terms ChiEOAbs.out : abs\chi_{v1v2v3}(\omega,\omega,0). The headers in these files contain information about the calculation. NOTES: - The routine has been written using notations of Ref. 2 - This routine does not symmetrize the tensor (up to now) - Sum over all the states and use occupation factors instead of looping only on resonant contributions
PARENTS
optic
CHILDREN
xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum
SOURCE
1627 subroutine linelop(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,evalv,occv,efermi, & 1628 pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,do_antiresonant,ncid,comm) 1629 1630 1631 !This section has been created automatically by the script Abilint (TD). 1632 !Do not modify the following lines by hand. 1633 #undef ABI_FUNC 1634 #define ABI_FUNC 'linelop' 1635 !End of the abilint section 1636 1637 implicit none 1638 1639 !Arguments ------------------------------------ 1640 !no_abirules 1641 integer, intent(in) :: icomp,itemp,nspin, ncid 1642 real(dp), intent(in) :: omega 1643 integer, intent(in) :: nkpt 1644 real(dp), intent(in) :: wkpt(nkpt) 1645 integer, intent(in) :: nsymcrys 1646 real(dp), intent(in) :: symcrys(3,3,nsymcrys) 1647 integer, intent(in) :: nstval 1648 real(dp), intent(in) :: evalv(nstval,nkpt,nspin) 1649 real(dp), intent(in) :: occv(nstval,nkpt,nspin) 1650 real(dp), intent(in) :: efermi 1651 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin) 1652 integer, intent(in) :: v1 1653 integer, intent(in) :: v2 1654 integer, intent(in) :: v3 1655 integer, intent(in) :: nmesh 1656 integer, intent(in) :: comm 1657 real(dp), intent(in) :: de 1658 real(dp), intent(in) :: sc 1659 real(dp), intent(in) :: brod 1660 real(dp), intent(in) :: tol 1661 character(len=*), intent(in) :: fnam 1662 logical, intent(in) :: do_antiresonant 1663 1664 !Local variables ------------------------- 1665 !no_abirules 1666 ! present calculation related (user specific) 1667 integer :: iw 1668 integer :: i,j,k,lx,ly,lz 1669 integer :: isp,isym,ik 1670 integer :: ist1,istl,istn,istm 1671 real(dp) :: ha2ev 1672 real(dp) :: t1,t2,t3,tst 1673 real(dp) :: ene,totre,totabs,totim 1674 real(dp) :: el,en,em 1675 real(dp) :: emin,emax,my_emin,my_emax 1676 real(dp) :: const_esu,const_au,au2esu 1677 real(dp) :: wmn,wnm,wln,wnl,wml,wlm 1678 complex(dpc) :: idel,w,zi 1679 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5 1680 ! local allocatable arrays 1681 real(dp), allocatable :: s(:,:) 1682 real(dp), allocatable :: sym(:,:,:) 1683 integer :: start4(4),count4(4) 1684 ! DBYG 1685 integer :: istp 1686 real(dp) :: ep, wmp, wpn 1687 real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included ! 1688 real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, fmn 1689 complex(dpc), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a} 1690 complex(dpc), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a} 1691 complex(dpc), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k) 1692 complex(dpc), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c 1693 complex(dpc), allocatable :: chi(:) ! \chi_{II}^{abc}(-\omega,\omega,0) 1694 complex(dpc), allocatable :: eta(:) ! \eta_{II}^{abc}(-\omega,\omega,0) 1695 complex(dpc), allocatable :: sigma(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0) 1696 complex(dpc), allocatable :: chi2tot(:) 1697 complex(dpc) :: num1, num2, den1, den2, term1, term2 1698 complex(dpc) :: chi1, chi1_1, chi1_2, chi2_1b, chi2_2b 1699 complex(dpc), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega) 1700 complex(dpc) :: eta1, eta2, eta2_1, eta2_2 1701 complex(dpc) :: sigma1, sigma1_1, sigma1_2, sigma2 1702 !Parallelism 1703 integer :: my_rank, nproc, master=0 1704 integer :: ierr 1705 integer :: my_k1, my_k2 1706 character(500) :: msg 1707 integer :: fout1,fout2,fout3,fout4,fout5 1708 1709 ! ********************************************************************* 1710 1711 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 1712 1713 !calculate the constant 1714 zi=(0._dp,1._dp) 1715 idel=zi*brod 1716 ! Disable symmetries for now 1717 const_au=-2._dp/(omega*dble(nsymcrys)) 1718 au2esu=5.8300348177d-8 1719 const_esu=const_au*au2esu 1720 ha2ev=13.60569172*2._dp 1721 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1722 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev 1723 !bohr: 5.2917ifc nlinopt.f907E-11 1724 !c: 2.99792458 velocity of sound 1725 !ry2ev: 13.60569172 1726 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2) 1727 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3) 1728 !mass comes from converting P_mn to r_mn 1729 !hbar^3 comes from converting all frequencies to energies in denominator 1730 !hbar^3 comes from operator for momentum (hbar/i nabla) 1731 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1732 !output file names 1733 fnam1=trim(fnam)//'-ChiEOTotIm.out' 1734 fnam2=trim(fnam)//'-ChiEOTotRe.out' 1735 fnam3=trim(fnam)//'-ChiEOIm.out' 1736 fnam4=trim(fnam)//'-ChiEORe.out' 1737 fnam5=trim(fnam)//'-ChiEOAbs.out' 1738 !!!!!!!!!!!!!!!!!!!!!!!!!!!! 1739 ! fool proof: 1740 !!!!!!!!!!!!!!!!!!!!!!!!!!!! 1741 !If there exists inversion symmetry exit with a mesg. 1742 tst=1.d-09 1743 do isym=1,nsymcrys 1744 t1=symcrys(1,1,isym)+1 1745 t2=symcrys(2,2,isym)+1 1746 t3=symcrys(3,3,isym)+1 1747 ! test if diagonal elements are -1 1748 if (abs(t1).lt.tst.and.abs(t2).lt.tst.and.abs(t3).lt.tst) then 1749 ! test if off-diagonal elements are zero 1750 if (abs(symcrys(1,2,isym)).lt.tst.and.abs(symcrys(1,3,isym)).lt.tst & 1751 .and.abs(symcrys(2,1,isym)).lt.tst.and.abs(symcrys(2,3,isym)).lt.tst.and. & 1752 abs(symcrys(3,1,isym)).lt.tst.and.abs(symcrys(3,2,isym)).lt.tst) then 1753 write(std_out,*) '-----------------------------------------' 1754 write(std_out,*) ' the crystal has inversion symmetry ' 1755 write(std_out,*) ' the LEO susceptibility is zero ' 1756 write(std_out,*) '-----------------------------------------' 1757 MSG_ERROR("Aborting now") 1758 end if 1759 end if 1760 end do 1761 !check polarisation 1762 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 1763 write(std_out,*) '---------------------------------------------' 1764 write(std_out,*) ' Error in linelop: ' 1765 write(std_out,*) ' the polarisation directions incorrect ' 1766 write(std_out,*) ' 1=x, 2=y and 3=z ' 1767 write(std_out,*) '---------------------------------------------' 1768 MSG_ERROR("Aborting now") 1769 end if 1770 !number of energy mesh points 1771 if (nmesh.le.0) then 1772 write(std_out,*) '---------------------------------------------' 1773 write(std_out,*) ' Error in linelop: ' 1774 write(std_out,*) ' number of energy mesh points incorrect ' 1775 write(std_out,*) ' number has to be integer greater than 0 ' 1776 write(std_out,*) ' nmesh*de = max energy for calculation ' 1777 write(std_out,*) '---------------------------------------------' 1778 MSG_ERROR("Aborting now") 1779 end if 1780 !step in energy 1781 if (de.le.0._dp) then 1782 write(std_out,*) '---------------------------------------------' 1783 write(std_out,*) ' Error in nlinopt: ' 1784 write(std_out,*) ' energy step is incorrect ' 1785 write(std_out,*) ' number has to real greater than 0.0 ' 1786 write(std_out,*) ' nmesh*de = max energy for calculation ' 1787 write(std_out,*) '---------------------------------------------' 1788 MSG_ERROR("Aborting now") 1789 end if 1790 !broadening 1791 if (brod.gt.0.009) then 1792 write(std_out,*) '---------------------------------------------' 1793 write(std_out,*) ' ATTENTION: broadening is quite high ' 1794 write(std_out,*) ' ideally should be less than 0.005 ' 1795 write(std_out,*) '---------------------------------------------' 1796 else if (brod.gt.0.015) then 1797 write(std_out,*) '----------------------------------------' 1798 write(std_out,*) ' ATTENTION: broadening is too high ' 1799 write(std_out,*) ' ideally should be less than 0.005 ' 1800 write(std_out,*) '----------------------------------------' 1801 end if 1802 !tolerance 1803 if (tol.gt.0.006) then 1804 write(std_out,*) '----------------------------------------' 1805 write(std_out,*) ' ATTENTION: tolerance is too high ' 1806 write(std_out,*) ' ideally should be less than 0.004 ' 1807 write(std_out,*) '----------------------------------------' 1808 end if 1809 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1810 !fool proof ends 1811 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1812 1813 ! 1814 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1815 !allocate local arrays 1816 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1817 1818 ABI_MALLOC(enk,(nstval)) 1819 ABI_MALLOC(delta,(nstval,nstval,3)) 1820 ABI_MALLOC(rmnbc,(nstval,nstval,3,3)) 1821 ABI_MALLOC(roverw,(nstval,nstval,3,3)) 1822 ABI_MALLOC(rmna,(nstval,nstval,3)) 1823 ABI_MALLOC(chi,(nmesh)) 1824 ABI_MALLOC(eta,(nmesh)) 1825 ABI_MALLOC(sigma,(nmesh)) 1826 ABI_MALLOC(chi2,(nmesh)) 1827 ABI_MALLOC(sym,(3,3,3)) 1828 ABI_MALLOC(s,(3,3)) 1829 1830 !generate the symmetrizing tensor 1831 sym(:,:,:)=0._dp 1832 do isym=1,nsymcrys 1833 s(:,:)=symcrys(:,:,isym) 1834 do i=1,3 1835 do j=1,3 1836 do k=1,3 1837 sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3)) 1838 end do 1839 end do 1840 end do 1841 end do 1842 1843 !initialise 1844 delta(:,:,:)=0._dp 1845 rmnbc(:,:,:,:)=0._dp 1846 chi(:)=0._dp 1847 chi2(:) = 0._dp 1848 eta(:)=0._dp 1849 sigma(:)=0._dp 1850 my_emin=HUGE(0._dp) 1851 my_emax=-HUGE(0._dp) 1852 1853 ! Split work 1854 call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr) 1855 1856 ! loop over kpts 1857 do ik=my_k1,my_k2 1858 write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt 1859 ! loop over spins 1860 do isp=1,nspin 1861 ! Calculate the scissor corrected energies and the energy window 1862 do ist1=1,nstval 1863 en = evalv(ist1,ik,isp) 1864 my_emin=min(my_emin,en) 1865 my_emax=max(my_emax,en) 1866 if(en > efermi) then 1867 en = en + sc 1868 end if 1869 enk(ist1) = en 1870 end do 1871 1872 ! calculate \Delta_nm and r_mn^a 1873 do istn=1,nstval 1874 en = enk(istn) 1875 do istm=1,nstval 1876 em = enk(istm) 1877 wmn = em - en 1878 delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp) 1879 if(abs(wmn) < tol) then 1880 rmna(istm,istn,1:3) = 0._dp 1881 else 1882 rmna(istm,istn,1:3)=-zi*pmat(istm,istn,ik,1:3,isp)/wmn 1883 end if 1884 end do 1885 end do 1886 ! calculate \r^b_mn;c 1887 do istm=1,nstval 1888 em = enk(istm) 1889 do istn=1,nstval 1890 en = enk(istn) 1891 wmn = em - en 1892 if(abs(wmn) > tol) then 1893 do ly = 1,3 1894 do lz = 1,3 1895 num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly)) 1896 den1 = wmn 1897 term1 = num1/den1 1898 term2 = 0._dp 1899 do istp=1,nstval 1900 ep = enk(istp) 1901 wmp = em - ep 1902 wpn = ep - en 1903 num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly)) 1904 den2 = wmn 1905 term2 = term2 + (num2/den2) 1906 end do 1907 rmnbc(istm,istn,ly,lz) = -term1-(zi*term2) 1908 roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz) 1909 end do 1910 end do 1911 end if 1912 end do 1913 end do 1914 1915 ! initialise the factors 1916 ! start the calculation 1917 do istn=1,nstval 1918 en=enk(istn) 1919 if (do_antiresonant .and. en .ge. efermi) then 1920 cycle 1921 end if 1922 fn=occv(istn,ik,isp) 1923 do istm=1,nstval 1924 em=enk(istm) 1925 if (do_antiresonant .and. em .le. efermi) then 1926 cycle 1927 end if 1928 wmn=em-en 1929 wnm=-wmn 1930 fm = occv(istm,ik,isp) 1931 fnm = fn - fm 1932 fmn = fm - fn 1933 eta1 = 0._dp 1934 eta2_1 = 0._dp 1935 eta2_2 = 0._dp 1936 sigma1_1 = 0._dp 1937 sigma1_2 = 0._dp 1938 sigma2 = 0._dp 1939 if(abs(wmn) > tol) then 1940 do lx = 1,3 1941 do ly = 1,3 1942 do lz = 1,3 1943 eta1 = eta1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*(roverw(istm,istn,lz,ly))) 1944 eta2_1 = eta2_1 + sym(lx,ly,lz)*(fnm*(rmna(istn,istm,lx)*rmnbc(istm,istn,ly,lz))) 1945 eta2_2 = eta2_2 + sym(lx,ly,lz)*(fnm*(rmnbc(istn,istm,lx,lz)*rmna(istm,istn,ly))) 1946 sigma1_1 = sigma1_1 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,ly)*rmna(istm,istn,lz))/(wmn**2) 1947 sigma1_2 = sigma1_2 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,lz)*rmna(istm,istn,ly))/(wmn**2) 1948 sigma2 = sigma2 + sym(lx,ly,lz)*(fnm*rmnbc(istn,istm,lz,lx)*rmna(istm,istn,ly))/wmn 1949 end do 1950 end do 1951 end do 1952 end if 1953 chi1_1 = 0._dp 1954 chi1_2 = 0._dp 1955 chi2_1b = 0._dp 1956 chi2_2b = 0._dp 1957 chi2(:) = 0._dp 1958 ! Three band terms 1959 do istl=1,nstval 1960 el=enk(istl) 1961 fl = occv(istl,ik,isp) 1962 wlm = el-em 1963 wln = el-en 1964 wnl = en-el 1965 wml = em-el 1966 fnl = fn-fl 1967 fln = fl-fn 1968 fml = fm-fl 1969 do lx = 1,3 1970 do ly = 1,3 1971 do lz = 1,3 1972 if(abs(wlm) > tol) then 1973 chi1_1 = chi1_1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,lz)*rmna(istl,istn,ly))/(wlm) 1974 chi2_1b = chi2_1b + sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*rmna(istl,istm,lz)*rmna(istm,istn,ly))/(wlm) 1975 end if 1976 if(abs(wln) > tol) then 1977 chi1_2 = chi1_2 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,ly)*rmna(istl,istn,lz))/(wln) 1978 chi2_2b = chi2_2b + sym(lx,ly,lz)*(fmn*rmna(istl,istm,lx)*rmna(istm,istn,ly)*rmna(istn,istl,lz))/(wnl) 1979 end if 1980 end do 1981 end do 1982 end do 1983 end do 1984 1985 sigma1 = 0.5_dp*(sigma1_1-sigma1_2) 1986 eta2 = 0.5_dp*(eta2_1-eta2_2) 1987 chi1 = chi1_1 + chi1_2 1988 ! 1989 ! calculate over the desired energy mesh and sum over k-points 1990 do iw=1,nmesh 1991 w=(iw-1)*de+idel 1992 ! Better way to compute it 1993 chi(iw) = chi(iw) + 0.5_dp*wkpt(ik)*((chi1/(wmn-w)) + ((chi2_1b+chi2_2b)/(wmn-w)))*const_esu 1994 eta(iw) = eta(iw) + 0.5_dp*zi*wkpt(ik)*((eta1/(wmn-w)) + (eta2/((wmn-w)**2)))*const_esu 1995 sigma(iw) = sigma(iw) + 0.5_dp*zi*wkpt(ik)*((sigma1/(wmn-w))- (sigma2/(wmn-w)))*const_esu 1996 end do 1997 end do ! istn and istm 1998 end do 1999 end do ! spins 2000 end do ! k-points 2001 2002 call xmpi_sum(chi,comm,ierr) 2003 call xmpi_sum(eta,comm,ierr) 2004 call xmpi_sum(sigma,comm,ierr) 2005 call xmpi_min(my_emin,emin,comm,ierr) 2006 call xmpi_max(my_emax,emax,comm,ierr) 2007 2008 if (my_rank == master) then 2009 2010 if (ncid /= nctk_noid) then 2011 start4 = [1, 1, icomp, itemp] 2012 count4 = [2, nmesh, 1, 1] 2013 ABI_MALLOC(chi2tot, (nmesh)) 2014 chi2tot = chi + eta + sigma 2015 #ifdef HAVE_NETCDF 2016 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi"), c2r(chi), start=start4, count=count4)) 2017 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_eta"), c2r(eta), start=start4, count=count4)) 2018 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_sigma"), c2r(sigma), start=start4, count=count4)) 2019 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi2tot"), c2r(chi2tot), start=start4, count=count4)) 2020 #endif 2021 ABI_FREE(chi2tot) 2022 end if 2023 2024 ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000) 2025 if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then 2026 MSG_ERROR(msg) 2027 end if 2028 if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then 2029 MSG_ERROR(msg) 2030 end if 2031 if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then 2032 MSG_ERROR(msg) 2033 end if 2034 if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then 2035 MSG_ERROR(msg) 2036 end if 2037 if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then 2038 MSG_ERROR(msg) 2039 end if 2040 ! write headers 2041 write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3 2042 write(fout1, '(a,es16.6)' ) ' #tolerance:',tol 2043 write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha' 2044 write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha' 2045 write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2046 write(fout1, '(a)' )' # Energy Tot-Im Chi(-w,w,0) Tot-Im Chi(-w,w,0)' 2047 write(fout1, '(a)' )' # eV *10^-7 esu *10^-12 m/V SI units ' 2048 write(fout1, '(a)' )' # ' 2049 2050 write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3 2051 write(fout2, '(a,es16.6)') ' #tolerance:',tol 2052 write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2053 write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2054 write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2055 write(fout2, '(a)')' # Energy Tot-Re Chi(-w,w,0) Tot-Re Chi(-w,w,0)' 2056 write(fout2, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units ' 2057 write(fout2, '(a)')' # ' 2058 2059 write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2060 write(fout3, '(a,es16.6)') ' #tolerance:',tol 2061 write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2062 write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2063 write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2064 write(fout3, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)' 2065 write(fout3, '(a)')' # in esu' 2066 write(fout3, '(a)')' # ' 2067 2068 write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2069 write(fout4, '(a,es16.6)') ' #tolerance:',tol 2070 write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2071 write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2072 write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2073 write(fout4, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)' 2074 write(fout4, '(a)')' # in esu' 2075 write(fout4, '(a)')' # ' 2076 2077 write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2078 write(fout5, '(a,es16.6)') ' #tolerance:',tol 2079 write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2080 write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2081 write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2082 write(fout5, '(a)')' # Energy(eV) |TotChi(-w,w,0)| |Tot Chi(-w,w,0)|' 2083 write(fout5, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units ' 2084 write(fout5, '(a)')' # ' 2085 2086 totim=0._dp 2087 totre=0._dp 2088 totabs=0._dp 2089 do iw=2,nmesh 2090 ene=(iw-1)*de 2091 ene=ene*ha2ev 2092 totim=aimag(chi(iw)+eta(iw)+sigma(iw))/1.d-7 2093 write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 2094 totim=0._dp 2095 totre=dble(chi(iw)+eta(iw)+sigma(iw))/1.d-7 2096 write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 2097 totre=0._dp 2098 write(fout3,'(f15.6,3es15.6)') ene,aimag(chi(iw))/1.d-7, & 2099 aimag(eta(iw))/1.d-7,aimag(sigma(iw))/1.d-7 2100 write(fout4,'(f15.6,3es15.6)') ene,dble(chi(iw))/1.d-7, & 2101 dble(eta(iw))/1.d-7,dble(sigma(iw))/1.d-7 2102 totabs=abs(chi(iw)+eta(iw)+sigma(iw))/1.d-7 2103 write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 2104 totabs=0._dp 2105 end do 2106 2107 close(fout1) 2108 close(fout2) 2109 close(fout3) 2110 close(fout4) 2111 close(fout5) 2112 ! print information 2113 write(std_out,*) ' ' 2114 write(std_out,*) 'information about calculation just performed:' 2115 write(std_out,*) ' ' 2116 write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of LEO susceptibility' 2117 write(std_out,*) 'tolerance:',tol 2118 if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high' 2119 write(std_out,*) 'broadening:',brod,'Hartree' 2120 if (brod.gt.0.009) then 2121 write(std_out,*) ' ' 2122 write(std_out,*) 'ATTENTION: broadening is quite high' 2123 write(std_out,*) ' ' 2124 else if (brod.gt.0.015) then 2125 write(std_out,*) ' ' 2126 write(std_out,*) 'ATTENTION: broadening is too high' 2127 write(std_out,*) ' ' 2128 end if 2129 write(std_out,*) 'scissors shift:',sc,'Hartree' 2130 write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree' 2131 2132 end if 2133 2134 ! deallocate local arrays 2135 ABI_FREE(enk) 2136 ABI_FREE(delta) 2137 ABI_FREE(rmnbc) 2138 ABI_FREE(roverw) 2139 ABI_FREE(rmna) 2140 ABI_FREE(chi) 2141 ABI_FREE(chi2) 2142 ABI_FREE(eta) 2143 ABI_FREE(sigma) 2144 ABI_FREE(s) 2145 ABI_FREE(sym) 2146 2147 end subroutine linelop
m_optic_tools/linopt [ Functions ]
[ Top ] [ m_optic_tools ] [ Functions ]
NAME
linopt
FUNCTION
Compute optical frequency dependent dielectric function for semiconductors
INPUTS
icomp=Sequential index associated to computed tensor components (used for netcdf output) itemp=Temperature index (used for netcdf output) nspin=number of spins(integer) omega=crystal volume in au (real) nkpt=total number of kpoints (integer) wkpt(nkpt)=weights of kpoints (real) nsymcrys=number of crystal symmetry operations(integer) symcrys(3,3,nsymcrys)=symmetry operations in cartisian coordinates(real) nstval=total number of valence states(integer) occv(nstval,nkpt,nspin)=occupation number for each band(real) evalv(nstval,nkpt,nspin)=eigen value for each band in Ha(real) efermi=Fermi energy in Ha(real) pmat(nstval,nstval,nkpt,3,nspin)=momentum matrix elements in cartesian coordinates(complex) v1,v2=desired component of the dielectric function(integer) 1=x,2=y,3=z nmesh=desired number of energy mesh points(integer) de=desired step in energy(real); nmesh*de=maximum energy sc=scissors shift in Ha(real) brod=broadening in Ha(real) fnam=root for filename that will contain the output filename will be trim(fnam)//'-linopt.out' ncid=Netcdf id to save output data.
SIDE EFFECTS
Dielectric function for semiconductors, on a desired energy mesh and for a desired direction of polarisation is written to file. The output is in a file named trim(fnam)//'-linopt.out' and contains Im(\epsilon_{v1v2}(\omega), Re(\epsilon_{v1v2}(\omega) and abs(\epsilon_{v1v2}(\omega). Comment: Right now the routine sums over the kpoints. In future linear tetrahedron method should be useful.
PARENTS
optic
CHILDREN
xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum
SOURCE
451 subroutine linopt(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,KSBSt,EPBSt,efermi,pmat, & 452 v1,v2,nmesh,de,sc,brod,fnam,ncid,comm) 453 454 455 !This section has been created automatically by the script Abilint (TD). 456 !Do not modify the following lines by hand. 457 #undef ABI_FUNC 458 #define ABI_FUNC 'linopt' 459 !End of the abilint section 460 461 implicit none 462 463 !Arguments ------------------------------------ 464 !no_abirules 465 integer, intent(in) :: icomp,itemp,nspin,ncid 466 real(dp), intent(in) :: omega 467 integer, intent(in) :: nkpt 468 real(dp), intent(in) :: wkpt(nkpt) 469 integer, intent(in) :: nsymcrys 470 real(dp), intent(in) :: symcrys(3,3,nsymcrys) 471 integer, intent(in) :: nstval 472 type(ebands_t),intent(in) :: KSBSt,EPBSt 473 real(dp), intent(in) :: efermi 474 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin) 475 integer, intent(in) :: v1 476 integer, intent(in) :: v2 477 integer, intent(in) :: nmesh 478 real(dp), intent(in) :: de 479 real(dp), intent(in) :: sc 480 real(dp), intent(in) :: brod 481 character(len=*), intent(in) :: fnam 482 integer, intent(in) :: comm 483 484 485 486 !Local variables ------------------------- 487 !no_abirules 488 integer :: isp 489 integer :: i,j,isym,lx,ly,ik 490 integer :: ist1,ist2,iw 491 ! Parallelism 492 integer :: my_rank, nproc 493 integer,parameter :: master=0 494 integer :: my_k1, my_k2 495 #ifdef HAVE_NETCDF 496 integer :: ncerr 497 #endif 498 integer :: ierr 499 integer :: fout1 500 logical :: do_lifetime 501 complex(dpc) :: e1,e2,e12 502 complex(dpc) :: e1_ep,e2_ep,e12_ep 503 real(dp) :: deltav1v2 504 real(dp) :: ha2ev 505 real(dp) :: tmpabs 506 real(dp) :: renorm_factor,emin,emax 507 real(dp) :: ene 508 complex(dpc) :: b11,b12 509 complex(dpc) :: ieta,w 510 character(len=fnlen) :: fnam1 511 character(len=500) :: msg 512 ! local allocatable arrays 513 real(dp) :: s(3,3),sym(3,3) 514 complex(dpc), allocatable :: chi(:) 515 complex(dpc), allocatable :: eps(:) 516 517 ! ********************************************************************* 518 519 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 520 521 if (my_rank == master) then 522 !check polarisation 523 if (v1.le.0.or.v2.le.0.or.v1.gt.3.or.v2.gt.3) then 524 write(std_out,*) '---------------------------------------------' 525 write(std_out,*) ' Error in linopt: ' 526 write(std_out,*) ' the polarisation directions incorrect ' 527 write(std_out,*) ' 1=x and 2=y and 3=z ' 528 write(std_out,*) '---------------------------------------------' 529 MSG_ERROR("Aborting now") 530 end if 531 !number of energy mesh points 532 if (nmesh.le.0) then 533 write(std_out,*) '---------------------------------------------' 534 write(std_out,*) ' Error in linopt: ' 535 write(std_out,*) ' number of energy mesh points incorrect ' 536 write(std_out,*) ' number has to integer greater than 0 ' 537 write(std_out,*) ' nmesh*de = max energy for calculation ' 538 write(std_out,*) '---------------------------------------------' 539 MSG_ERROR("Aborting now") 540 end if 541 !step in energy 542 if (de.le.0._dp) then 543 write(std_out,*) '---------------------------------------------' 544 write(std_out,*) ' Error in linopt: ' 545 write(std_out,*) ' energy step is incorrect ' 546 write(std_out,*) ' number has to real greater than 0.0 ' 547 write(std_out,*) ' nmesh*de = max energy for calculation ' 548 write(std_out,*) '---------------------------------------------' 549 MSG_ERROR("Aborting now") 550 end if 551 !broadening 552 if (brod.gt.0.009) then 553 write(std_out,*) '---------------------------------------------' 554 write(std_out,*) ' ATTENTION: broadening is quite high ' 555 write(std_out,*) ' ideally should be less than 0.005 ' 556 write(std_out,*) '---------------------------------------------' 557 else if (brod.gt.0.015) then 558 write(std_out,*) '----------------------------------------' 559 write(std_out,*) ' ATTENTION: broadening is too high ' 560 write(std_out,*) ' ideally should be less than 0.005 ' 561 write(std_out,*) '----------------------------------------' 562 end if 563 !fermi energy 564 if(efermi<-1.0d4) then 565 write(std_out,*) '---------------------------------------------' 566 write(std_out,*) ' ATTENTION: Fermi energy seems extremely ' 567 write(std_out,*) ' low ' 568 write(std_out,*) '---------------------------------------------' 569 end if 570 !scissors operator 571 if (sc.lt.0._dp) then 572 write(std_out,*) '---------------------------------------------' 573 write(std_out,*) ' Error in linopt: ' 574 write(std_out,*) ' scissors shift is incorrect ' 575 write(std_out,*) ' number has to real greater than 0.0 ' 576 write(std_out,*) '---------------------------------------------' 577 MSG_ERROR("Aborting now") 578 end if 579 !fool proof end 580 end if 581 582 ABI_CHECK(KSBSt%mband==nstval, "The number of bands in the BSt should be equal to nstval !") 583 584 do_lifetime = allocated(EPBSt%lifetime) 585 ! TODO: activate this, and remove do_lifetime - always add it in even if 0. 586 ! if (.not. allocated(EPBSt%lifetime)) then 587 ! ABI_ALLOCATE(EPBSt%lifetime, (nstval, my_k2-my_k1+1, nspin)) 588 ! EPBSt%lifetime = zero 589 ! end if 590 591 !allocate local arrays 592 ABI_ALLOCATE(chi,(nmesh)) 593 ABI_ALLOCATE(eps,(nmesh)) 594 ieta=(0._dp,1._dp)*brod 595 renorm_factor=1._dp/(omega*dble(nsymcrys)) 596 ha2ev=13.60569172*2._dp 597 !output file names 598 fnam1=trim(fnam)//'-linopt.out' 599 !construct symmetrisation tensor 600 sym(:,:)=0._dp 601 do isym=1,nsymcrys 602 s(:,:)=symcrys(:,:,isym) 603 do i=1,3 604 do j=1,3 605 sym(i,j)=sym(i,j)+s(i,v1)*s(j,v2) 606 end do 607 end do 608 end do 609 610 !calculate the energy window 611 emin=0._dp 612 emax=0._dp 613 do ik=1,nkpt 614 do isp=1,nspin 615 do ist1=1,nstval 616 emin=min(emin,EPBSt%eig(ist1,ik,isp)) 617 emax=max(emax,EPBSt%eig(ist1,ik,isp)) 618 end do 619 end do 620 end do 621 622 ! Split work 623 call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr) 624 625 !start calculating linear optical response 626 chi(:)=0._dp 627 ! TODO: this loop should be outside the ik one, for speed and cache. 628 do isp=1,nspin 629 !do ik=1,nkpt 630 do ik=my_k1,my_k2 631 write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt 632 do ist1=1,nstval 633 e1=KSBSt%eig(ist1,ik,isp) 634 e1_ep=EPBSt%eig(ist1,ik,isp) 635 ! TODO: unless memory is a real issue, should set lifetimes to 0 and do this sum systematically 636 ! instead of putting an if statement in a loop! See above 637 if(do_lifetime) then 638 e1_ep = e1_ep + EPBSt%lifetime(ist1,ik,isp)*(0.0_dp,1.0_dp) 639 end if 640 ! if (e1.lt.efermi) then 641 ! do ist2=ist1,nstval 642 do ist2=1,nstval 643 e2=KSBSt%eig(ist2,ik,isp) 644 e2_ep=EPBSt%eig(ist2,ik,isp) 645 if(do_lifetime) then 646 e2_ep = e2_ep - EPBSt%lifetime(ist2,ik,isp)*(0.0_dp,1.0_dp) 647 end if 648 ! if (e2.gt.efermi) then 649 if (ist1.ne.ist2) then 650 ! scissors correction of momentum matrix 651 if(REAL(e1) > REAL(e2)) then 652 e12 = e1-e2+sc 653 else 654 e12 = e1-e2-sc 655 end if 656 if(REAL(e1_ep) > REAL(e2_ep)) then 657 e12_ep = e1_ep-e2_ep+sc 658 else 659 e12_ep = e1_ep-e2_ep-sc 660 end if 661 ! e12=e1-e2-sc 662 b11=0._dp 663 ! symmetrization of momentum matrix 664 do lx=1,3 665 do ly=1,3 666 b11=b11+(sym(lx,ly)*pmat(ist1,ist2,ik,lx,isp)* & 667 conjg(pmat(ist1,ist2,ik,ly,isp))) 668 end do 669 end do 670 b12=b11*renorm_factor*(1._dp/(e12**2)) 671 ! calculate on the desired energy grid 672 do iw=2,nmesh 673 w=(iw-1)*de+ieta 674 chi(iw)=chi(iw)+(wkpt(ik)*(KSBSt%occ(ist1,ik,isp)-KSBSt%occ(ist2,ik,isp))* & 675 (b12/(-e12_ep-w))) 676 end do 677 end if 678 end do ! states 679 ! end if 680 end do 681 end do ! spin 682 end do ! k-points 683 684 call xmpi_sum(chi,comm,ierr) 685 686 ! calculate epsilon 687 eps(1) = zero 688 deltav1v2=zero; if (v1 == v2) deltav1v2=one 689 do iw=2,nmesh 690 ene=(iw-1)*de 691 ene=ene*ha2ev 692 eps(iw)=deltav1v2+4._dp*pi*chi(iw) 693 end do 694 695 if (my_rank == master) then 696 ! open the output files 697 if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then 698 MSG_ERROR(msg) 699 end if 700 ! write the output 701 write(fout1, '(a,2i3,a)' )' #calculated the component:',v1,v2,' of dielectric function' 702 write(std_out,*) 'calculated the component:',v1,v2,' of dielectric function' 703 write(fout1, '(a,2es16.6)' ) ' #broadening:', real(ieta),aimag(ieta) 704 write(std_out,*) ' with broadening:',ieta 705 write(fout1, '(a,es16.6)' ) ' #scissors shift:',sc 706 write(std_out,*) 'and scissors shift:',sc 707 write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 708 write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 709 write(fout1,*) 710 write(fout1, '(a)' ) ' # Energy(eV) Im(eps(w))' 711 do iw=2,nmesh 712 ene=(iw-1)*de 713 ene=ene*ha2ev 714 write(fout1, '(2es16.6)' ) ene,aimag(eps(iw)) 715 end do 716 write(fout1,*) 717 write(fout1,*) 718 write(fout1, '(a)' ) ' # Energy(eV) Re(eps(w))' 719 do iw=2,nmesh 720 ene=(iw-1)*de 721 ene=ene*ha2ev 722 write(fout1, '(2es16.6)' ) ene,dble(eps(iw)) 723 end do 724 write(fout1,*) 725 write(fout1,*) 726 write(fout1, '(a)' )' # Energy(eV) abs(eps(w))' 727 do iw=2,nmesh 728 ene=(iw-1)*de 729 ene=ene*ha2ev 730 write(fout1, '(2es16.6)' ) ene,abs(eps(iw)) 731 end do 732 write(fout1,*) 733 write(fout1,*) 734 write(fout1, '(a)' )' # Energy(eV) Im(refractive index(w)) aka kappa' 735 do iw=2,nmesh 736 ene=(iw-1)*de 737 ene=ene*ha2ev 738 write(fout1, '(2es16.6)' ) ene,sqrt(half*(abs(eps(iw)) - dble(eps(iw)) )) 739 end do 740 write(fout1,*) 741 write(fout1,*) 742 write(fout1, '(a)' )' # Energy(eV) Re(refractive index(w)) aka n' 743 do iw=2,nmesh 744 ene=(iw-1)*de 745 ene=ene*ha2ev 746 write(fout1, '(2es16.6)' ) ene,sqrt(half*(abs(eps(iw)) + dble(eps(iw)) )) 747 end do 748 write(fout1,*) 749 write(fout1,*) 750 write(fout1, '(a)' )' # Energy(eV) Reflectivity(w) from vacuum, at normal incidence' 751 do iw=2,nmesh 752 ene=(iw-1)*de 753 ene=ene*ha2ev 754 write(fout1, '(2es16.6)' ) ene, sqrt(half*(abs(eps(iw)) + dble(eps(iw)) )) 755 end do 756 write(fout1,*) 757 write(fout1,*) 758 write(fout1, '(a)' )' # Energy(eV) absorption coeff (in m-1) = omega Im(eps) / c n(eps)' 759 do iw=2,nmesh 760 ene=(iw-1)*de 761 tmpabs=zero 762 if (abs(eps(iw)) + dble(eps(iw)) > zero) then 763 tmpabs = aimag(eps(iw))*ene / sqrt(half*( abs(eps(iw)) + dble(eps(iw)) )) / Sp_Lt / Bohr_meter 764 end if 765 write(fout1, '(2es16.6)' ) ha2ev*ene, tmpabs 766 end do 767 768 ! close output file 769 close(fout1) 770 771 #ifdef HAVE_NETCDF 772 if (ncid /= nctk_noid) then 773 ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_epsilon"), c2r(eps), start=[1, 1, icomp, itemp]) 774 NCF_CHECK(ncerr) 775 end if 776 #endif 777 end if 778 779 ABI_DEALLOCATE(chi) 780 ABI_FREE(eps) 781 782 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
838 subroutine nlinopt(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,evalv,efermi, & 839 pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,ncid,comm) 840 841 842 !This section has been created automatically by the script Abilint (TD). 843 !Do not modify the following lines by hand. 844 #undef ABI_FUNC 845 #define ABI_FUNC 'nlinopt' 846 !End of the abilint section 847 848 implicit none 849 850 !Arguments ------------------------------------ 851 !no_abirules 852 integer, intent(in) :: icomp,itemp,nspin, ncid 853 real(dp), intent(in) :: omega 854 integer, intent(in) :: nkpt 855 real(dp), intent(in) :: wkpt(nkpt) 856 integer, intent(in) :: nsymcrys 857 real(dp), intent(in) :: symcrys(3,3,nsymcrys) 858 integer, intent(in) :: nstval 859 real(dp), intent(in) :: evalv(nstval,nkpt,nspin) 860 real(dp), intent(in) :: efermi 861 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin) 862 integer, intent(in) :: v1 863 integer, intent(in) :: v2 864 integer, intent(in) :: v3 865 integer, intent(in) :: nmesh 866 integer, intent(in) :: comm 867 real(dp), intent(in) :: de 868 real(dp), intent(in) :: sc 869 real(dp), intent(in) :: brod 870 real(dp), intent(in) :: tol 871 character(len=*), intent(in) :: fnam 872 873 !Local variables ------------------------- 874 integer :: iw 875 integer :: i,j,k,lx,ly,lz 876 integer :: isp,isym,ik 877 integer :: ist1,ist2,istl,istn,istm 878 integer,parameter :: master=0 879 integer :: my_rank, nproc 880 integer :: my_k1, my_k2 881 integer :: ierr 882 integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7 883 real(dp) :: f1,f2,f3 884 real(dp) :: ha2ev 885 real(dp) :: t1,t2,t3,tst 886 real(dp) :: ene,totre,totabs,totim 887 real(dp) :: e1,e2,el,en,em 888 real(dp) :: emin,emax,my_emin,my_emax 889 real(dp) :: const_esu,const_au,au2esu 890 real(dp) :: wmn,wnm,wln,wnl,wml,wlm 891 complex(dpc) :: idel,w,zi 892 complex(dpc) :: mat2w,mat1w1,mat1w2,mat2w_tra,mat1w3_tra 893 complex(dpc) :: b111,b121,b131,b112,b122,b132,b113,b123,b133 894 complex(dpc) :: b241,b242,b243,b221,b222,b223,b211,b212,b213,b231 895 complex(dpc) :: b311,b312,b313,b331 896 complex(dpc) :: b24,b21_22,b11,b12_13,b31_32 897 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7 898 character(500) :: msg 899 ! local allocatable arrays 900 integer :: start4(4),count4(4) 901 real(dp) :: s(3,3),sym(3,3,3) 902 complex(dpc), allocatable :: px(:,:,:,:,:) 903 complex(dpc), allocatable :: py(:,:,:,:,:) 904 complex(dpc), allocatable :: pz(:,:,:,:,:) 905 complex(dpc), allocatable :: delta(:,:,:) 906 complex(dpc), allocatable :: inter2w(:) 907 complex(dpc), allocatable :: inter1w(:) 908 complex(dpc), allocatable :: intra2w(:) 909 complex(dpc), allocatable :: intra1w(:) 910 complex(dpc), allocatable :: intra1wS(:),chi2tot(:) 911 912 ! ********************************************************************* 913 914 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 915 916 !calculate the constant 917 zi=(0._dp,1._dp) 918 idel=zi*brod 919 const_au=-2._dp/(omega*dble(nsymcrys)) 920 au2esu=5.8300348177d-8 921 const_esu=const_au*au2esu 922 ha2ev=13.60569172*2._dp 923 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 924 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev 925 !bohr: 5.2917ifc nlinopt.f907E-11 926 !c: 2.99792458 velocity of sound 927 !ry2ev: 13.60569172 928 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2) 929 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3) 930 !mass comes from converting P_mn to r_mn 931 !hbar^3 comes from converting all frequencies to energies in denominator 932 !hbar^3 comes from operator for momentum (hbar/i nabla) 933 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 934 !output file names 935 fnam1=trim(fnam)//'-ChiTotIm.out' 936 fnam2=trim(fnam)//'-ChiTotRe.out' 937 fnam3=trim(fnam)//'-ChiIm.out' 938 fnam4=trim(fnam)//'-ChiRe.out' 939 fnam5=trim(fnam)//'-ChiAbs.out' 940 fnam6=trim(fnam)//'-ChiImDec.out' 941 fnam7=trim(fnam)//'-ChiReDec.out' 942 943 if(my_rank == master) then 944 !If there exists inversion symmetry exit with a message. 945 tst=1.d-09 946 do isym=1,nsymcrys 947 t1=symcrys(1,1,isym)+1 948 t2=symcrys(2,2,isym)+1 949 t3=symcrys(3,3,isym)+1 950 ! test if diagonal elements are -1 951 if (abs(t1).lt.tst.and.abs(t2).lt.tst.and.abs(t3).lt.tst) then 952 ! test if off-diagonal elements are zero 953 if (abs(symcrys(1,2,isym)).lt.tst.and.abs(symcrys(1,3,isym)).lt.tst & 954 .and.abs(symcrys(2,1,isym)).lt.tst.and.abs(symcrys(2,3,isym)).lt.tst.and. & 955 abs(symcrys(3,1,isym)).lt.tst.and.abs(symcrys(3,2,isym)).lt.tst) then 956 write(std_out,*) '-----------------------------------------' 957 write(std_out,*) ' the crystal has inversion symmetry ' 958 write(std_out,*) ' the SHG susceptibility is zero ' 959 write(std_out,*) '-----------------------------------------' 960 MSG_ERROR("Aborting now") 961 end if 962 end if 963 end do 964 !check polarisation 965 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 966 write(std_out,*) '---------------------------------------------' 967 write(std_out,*) ' Error in nlinopt: ' 968 write(std_out,*) ' the polarisation directions incorrect ' 969 write(std_out,*) ' 1=x, 2=y and 3=z ' 970 write(std_out,*) '---------------------------------------------' 971 MSG_ERROR("Aborting now") 972 end if 973 !number of energy mesh points 974 if (nmesh.le.0) then 975 write(std_out,*) '---------------------------------------------' 976 write(std_out,*) ' Error in nlinopt: ' 977 write(std_out,*) ' number of energy mesh points incorrect ' 978 write(std_out,*) ' number has to be integer greater than 0 ' 979 write(std_out,*) ' nmesh*de = max energy for calculation ' 980 write(std_out,*) '---------------------------------------------' 981 MSG_ERROR("Aborting now") 982 end if 983 !step in energy 984 if (de.le.0._dp) then 985 write(std_out,*) '---------------------------------------------' 986 write(std_out,*) ' Error in nlinopt: ' 987 write(std_out,*) ' energy step is incorrect ' 988 write(std_out,*) ' number has to real greater than 0.0 ' 989 write(std_out,*) ' nmesh*de = max energy for calculation ' 990 write(std_out,*) '---------------------------------------------' 991 MSG_ERROR("Aborting now") 992 end if 993 !broadening 994 if (brod.gt.0.009) then 995 write(std_out,*) '---------------------------------------------' 996 write(std_out,*) ' ATTENTION: broadening is quite high ' 997 write(std_out,*) ' ideally should be less than 0.005 ' 998 write(std_out,*) '---------------------------------------------' 999 else if (brod.gt.0.015) then 1000 write(std_out,*) '----------------------------------------' 1001 write(std_out,*) ' ATTENTION: broadening is too high ' 1002 write(std_out,*) ' ideally should be less than 0.005 ' 1003 write(std_out,*) '----------------------------------------' 1004 end if 1005 !tolerance 1006 if (tol.gt.0.006) then 1007 write(std_out,*) '----------------------------------------' 1008 write(std_out,*) ' ATTENTION: tolerance is too high ' 1009 write(std_out,*) ' ideally should be less than 0.004 ' 1010 write(std_out,*) '----------------------------------------' 1011 end if 1012 end if 1013 1014 !allocate local arrays 1015 ABI_ALLOCATE(px,(nstval,nstval,3,3,3)) 1016 ABI_ALLOCATE(py,(nstval,nstval,3,3,3)) 1017 ABI_ALLOCATE(pz,(nstval,nstval,3,3,3)) 1018 ABI_ALLOCATE(inter2w,(nmesh)) 1019 ABI_ALLOCATE(inter1w,(nmesh)) 1020 ABI_ALLOCATE(intra2w,(nmesh)) 1021 ABI_ALLOCATE(intra1w,(nmesh)) 1022 ABI_ALLOCATE(intra1wS,(nmesh)) 1023 ABI_ALLOCATE(delta,(nstval,nstval,3)) 1024 1025 !generate the symmetrizing tensor 1026 sym(:,:,:)=0._dp 1027 do isym=1,nsymcrys 1028 s(:,:)=symcrys(:,:,isym) 1029 do i=1,3 1030 do j=1,3 1031 do k=1,3 1032 sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3)) 1033 end do 1034 end do 1035 end do 1036 end do 1037 !DBYG 1038 ! Disable symmetries for now 1039 !sym(:,:,:) = 0._dp 1040 !sym(v1,v2,v3) = nsymcrys 1041 !ENDDBYG 1042 1043 ! Split work 1044 call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr) 1045 1046 !initialise 1047 inter2w(:)=0._dp 1048 inter1w(:)=0._dp 1049 intra2w(:)=0._dp 1050 intra1w(:)=0._dp 1051 intra1wS(:)=0._dp 1052 delta(:,:,:)=0._dp 1053 1054 my_emin=HUGE(0._dp) 1055 my_emax=-HUGE(0._dp) 1056 ! loop over kpts 1057 do ik=my_k1,my_k2 1058 write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt 1059 ! loop over spins 1060 do isp=1,nspin 1061 ! loop over states 1062 do ist1=1,nstval 1063 e1=evalv(ist1,ik,isp) 1064 if (e1.lt.efermi) then ! ist1 is a valence state 1065 do ist2=1,nstval 1066 e2=evalv(ist2,ik,isp) 1067 if (e2.gt.efermi) then ! ist2 is a conduction state 1068 ! symmetrize the momentum matrix elements 1069 do lx=1,3 1070 do ly=1,3 1071 do lz=1,3 1072 f1=sym(lx,ly,lz)+sym(lx,lz,ly) 1073 f2=sym(ly,lx,lz)+sym(ly,lz,lx) 1074 f3=sym(lz,lx,ly)+sym(lz,ly,lx) 1075 px(ist1,ist2,lx,ly,lz)=f1*pmat(ist1,ist2,ik,lx,isp) 1076 py(ist2,ist1,lx,ly,lz)=f2*pmat(ist2,ist1,ik,lx,isp) 1077 pz(ist2,ist1,lx,ly,lz)=f3*pmat(ist2,ist1,ik,lx,isp) 1078 end do 1079 end do 1080 end do 1081 ! end loop over states 1082 end if 1083 end do 1084 end if 1085 end do 1086 ! calculate the energy window and \Delta_nm 1087 do ist1=1,nstval 1088 my_emin=min(my_emin,evalv(ist1,ik,isp)) 1089 my_emax=max(my_emax,evalv(ist1,ik,isp)) 1090 do ist2=1,nstval 1091 delta(ist1,ist2,1:3)=pmat(ist1,ist1,ik,1:3,isp)-pmat(ist2,ist2,ik,1:3,isp) 1092 end do 1093 end do 1094 ! initialise the factors 1095 ! factors are named according to the Ref. article 2. 1096 b111=0._dp 1097 b121=0._dp 1098 b131=0._dp 1099 b112=0._dp 1100 b122=0._dp 1101 b132=0._dp 1102 b113=0._dp 1103 b123=0._dp 1104 b133=0._dp 1105 b211=0._dp 1106 b221=0._dp 1107 b212=0._dp 1108 b222=0._dp 1109 b213=0._dp 1110 b223=0._dp 1111 b231=0._dp 1112 b241=0._dp 1113 b242=0._dp 1114 b243=0._dp 1115 b311=0._dp 1116 b312=0._dp 1117 b313=0._dp 1118 b331=0._dp 1119 ! start the calculation 1120 do istn=1,nstval 1121 en=evalv(istn,ik,isp) 1122 if (en.lt.efermi) then ! istn is a valence state 1123 do istm=1,nstval 1124 em=evalv(istm,ik,isp) 1125 if (em.gt.efermi) then ! istm is a conduction state 1126 em = em + sc ! Should add the scissor to conduction energies 1127 wmn=em-en 1128 wnm=-wmn 1129 ! calculate the matrix elements for two band intraband term 1130 mat2w_tra=0._dp 1131 mat1w3_tra=0._dp 1132 do lx=1,3 1133 do ly=1,3 1134 do lz=1,3 1135 mat2w_tra=mat2w_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,lz,isp) & 1136 *delta(istm,istn,ly) 1137 mat1w3_tra=mat1w3_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,ly,isp) & 1138 *delta(istm,istn,lz) 1139 ! NOTE:: lx to ly m to n in pmat matrices respectively 1140 ! Changes are made so that this (b3) term is according to paper 1141 ! PRB48(Ref. 4) rather than PRB53(Ref 2) in which this term is incorrect 1142 end do 1143 end do 1144 end do 1145 b331=mat1w3_tra/wnm 1146 b11=0._dp 1147 b12_13=0._dp 1148 b24=0._dp 1149 b31_32=0._dp 1150 b21_22=0._dp 1151 1152 b231=8._dp*mat2w_tra/wmn 1153 b331=mat1w3_tra/(wnm) 1154 ! !!!!!!!!!!!!!!!!!!! 1155 ! istl < istn ! 1156 ! !!!!!!!!!!!!!!!!!!! 1157 do istl=1,istn-1 ! istl is a valence state below istn 1158 el=evalv(istl,ik,isp) 1159 wln=el-en ! do not add sc to the valence bands! 1160 wml=em-el 1161 wnl=-wln 1162 wlm=-wml 1163 ! calculate the matrix elements for three band terms 1164 mat2w=0._dp 1165 mat1w1=0._dp 1166 mat1w2=0._dp 1167 do lx=1,3 1168 do ly=1,3 1169 do lz=1,3 1170 1171 mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) & 1172 *pmat(istl,istn,ik,lz,isp)) 1173 1174 mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) & 1175 *pmat(istn,istl,ik,ly,isp)) 1176 1177 mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) & 1178 *pmat(istn,istl,ik,ly,isp)) 1179 end do 1180 end do 1181 end do 1182 b111=mat2w*(1._dp/(wln+wlm))*(1._dp/wlm) 1183 b121=mat1w1*(1._dp/(wnm+wlm))*(1._dp/wlm) 1184 b131=mat1w2*(1._dp/wlm) 1185 ! 1186 b221=0._dp 1187 b211=mat1w1/wml 1188 b241=-mat2w/wml 1189 ! 1190 b311=mat1w2/wlm 1191 if (abs(wln).gt.tol) then 1192 b111=b111/wln 1193 b121=b121/wln 1194 b131=b131/wln 1195 b221=mat1w2/wln 1196 b241=b241+(mat2w/wln) 1197 b311=b311+(mat1w1/wln) 1198 else 1199 b111=0._dp 1200 b121=0._dp 1201 b131=0._dp 1202 b221=0._dp 1203 end if 1204 t1=wln-wnm 1205 if (abs(t1).gt.tol) then 1206 b131=b131/t1 1207 else 1208 b131=0._dp 1209 end if 1210 b11=b11-2._dp*b111 1211 b12_13=b12_13+b121+b131 1212 b21_22=b21_22-b211+b221 1213 b24=b24+2._dp*b241 1214 b31_32=b31_32+b311 1215 ! end loop over istl 1216 end do 1217 1218 ! !!!!!!!!!!!!!!!!!!!!!!!!!!! 1219 ! istn < istl < istm ! 1220 ! !!!!!!!!!!!!!!!!!!!!!!!!!!! 1221 do istl=istn+1,istm-1 1222 el=evalv(istl,ik,isp) 1223 ! calculate the matrix elements for three band terms 1224 mat2w=0._dp 1225 mat1w1=0._dp 1226 mat1w2=0._dp 1227 do lx=1,3 1228 do ly=1,3 1229 do lz=1,3 1230 1231 mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) & 1232 *pmat(istl,istn,ik,lz,isp)) 1233 1234 mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) & 1235 *pmat(istn,istl,ik,ly,isp)) 1236 1237 mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) & 1238 *pmat(istn,istl,ik,ly,isp)) 1239 end do 1240 end do 1241 end do 1242 if (el.lt.efermi) then 1243 wln=el-en 1244 wnl=-wln 1245 wml=em-el 1246 wlm=-wml 1247 else 1248 el=el+sc 1249 wln=el-en 1250 wnl=-wln 1251 wml=em-el 1252 wlm=-wml 1253 end if 1254 ! 1255 b112=0._dp 1256 b122=mat1w1*(1._dp/(wnm+wlm)) 1257 b132=mat1w2*(1._dp/(wnm+wnl)) 1258 b242=0._dp 1259 b222=0._dp 1260 b212=0._dp 1261 b312=0._dp 1262 if (abs(wnl).gt.tol) then 1263 b112=mat2w/wln 1264 b122=b122/wnl 1265 b132=b132/wnl 1266 b242=mat2w/wln 1267 b222=mat1w2/wln 1268 b312=mat1w1/wln 1269 else 1270 b122=0._dp 1271 b132=0._dp 1272 end if 1273 if (abs(wlm).gt.tol) then 1274 b112=b112/wml 1275 b122=b122/wlm 1276 b132=b132/wlm 1277 b242=b242-(mat2w/wml) 1278 b212=mat1w1/wml 1279 b312=b312+(mat1w2/wlm) 1280 else 1281 b112=0._dp 1282 b122=0._dp 1283 b132=0._dp 1284 b212=0._dp 1285 end if 1286 t1=wlm-wnl 1287 if (abs(t1).gt.tol) then 1288 b112=b112/t1 1289 else 1290 b112=0._dp 1291 end if 1292 b11=b11+2._dp*b112 1293 b12_13=b12_13-b122+b132 1294 b24=b24+2._dp*b242 1295 b21_22=b21_22-b212+b222 1296 b31_32=b31_32+b312 1297 ! end loop over istl 1298 end do 1299 1300 ! !!!!!!!!!!!!!!!!!!!!! 1301 ! istl > istm ! 1302 ! !!!!!!!!!!!!!!!!!!!!! 1303 do istl=istm+1,nstval 1304 el=evalv(istl,ik,isp)+sc 1305 wln=el-en 1306 wnl=-wln 1307 wml=em-el 1308 wlm=-wml 1309 ! calculate the matrix elements for three band terms 1310 mat2w=0._dp 1311 mat1w1=0._dp 1312 mat1w2=0._dp 1313 do lx=1,3 1314 do ly=1,3 1315 do lz=1,3 1316 mat2w=mat2w+px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) & 1317 *pmat(istl,istn,ik,lz,isp) 1318 1319 mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) & 1320 *pmat(istn,istl,ik,ly,isp)) 1321 1322 mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) & 1323 *pmat(istn,istl,ik,ly,isp)) 1324 end do 1325 end do 1326 end do 1327 ! 1328 b113=mat2w*(1._dp/(wnl+wml))*(1._dp/wnl) 1329 b123=mat1w1*(1._dp/wnl) 1330 b133=mat1w2*(1._dp/wnl)*(1._dp/(wnl+wnm)) 1331 b243=mat2w/wln 1332 b223=mat1w2/wln 1333 b213=0._dp 1334 b313=-1._dp*mat1w1/wnl 1335 if (abs(wml).gt.tol) then 1336 b113=b113/wml 1337 b123=b123/wml 1338 b133=b133/wml 1339 b243=b243-(mat2w/wml) 1340 b213=mat1w1/wml 1341 b313=b313+(mat1w2/wlm) 1342 else 1343 b113=0._dp 1344 b123=0._dp 1345 b133=0._dp 1346 end if 1347 t1=wnm-wml 1348 if (abs(t1).gt.tol) then 1349 b123=b123/t1 1350 else 1351 b123=0._dp 1352 end if 1353 b11=b11+2._dp*b113 1354 b12_13=b12_13+b123-b133 1355 b21_22=b21_22-b213+b223 1356 b24=b24+2._dp*b243 1357 b31_32=b31_32+b313 1358 ! end loop over istl 1359 end do 1360 1361 b11=b11*zi*(1._dp/wnm)*const_esu 1362 b12_13=b12_13*zi*(1._dp/wnm)*const_esu 1363 b24=(b24+b231)*zi*(1._dp/(wnm**3))*const_esu 1364 b21_22=(b21_22)*zi*(1._dp/(wnm**3))*const_esu 1365 b31_32=(b31_32-b331)*zi*(1._dp/(wmn**3))*const_esu*0.5_dp 1366 ! calculate over the desired energy mesh and sum over k-points 1367 do iw=1,nmesh 1368 w=(iw-1)*de+idel 1369 inter2w(iw)=inter2w(iw)+(wkpt(ik)*(b11/(wmn-2._dp*w))) ! Inter(2w) from chi 1370 inter1w(iw)=inter1w(iw)+(wkpt(ik)*(b12_13/(wmn-w))) ! Inter(1w) from chi 1371 intra2w(iw)=intra2w(iw)+(wkpt(ik)*(b24/(wmn-2._dp*w))) ! Intra(2w) from eta 1372 intra1w(iw)=intra1w(iw)+(wkpt(ik)*((b21_22)/(wmn-w))) ! Intra(1w) from eta 1373 intra1wS(iw)=intra1wS(iw)+(wkpt(ik)*((b31_32)/(wmn-w))) ! Intra(1w) from sigma 1374 end do 1375 end if 1376 end do ! istn and istm 1377 end if 1378 end do 1379 end do ! spins 1380 end do ! k-points 1381 1382 call xmpi_sum(inter2w,comm,ierr) 1383 call xmpi_sum(inter1w,comm,ierr) 1384 call xmpi_sum(intra2w,comm,ierr) 1385 call xmpi_sum(intra1w,comm,ierr) 1386 call xmpi_sum(intra1wS,comm,ierr) 1387 call xmpi_min(my_emin,emin,comm,ierr) 1388 call xmpi_max(my_emax,emax,comm,ierr) 1389 1390 if (my_rank == master) then 1391 ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000) 1392 1393 if (ncid /= nctk_noid) then 1394 start4 = [1, 1, icomp, itemp] 1395 count4 = [2, nmesh, 1, 1] 1396 ABI_MALLOC(chi2tot, (nmesh)) 1397 chi2tot = inter2w + inter1w + intra2w + intra1w + intra1wS 1398 #ifdef HAVE_NETCDF 1399 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter2w"), c2r(inter2w), start=start4, count=count4)) 1400 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter1w"), c2r(inter1w), start=start4, count=count4)) 1401 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra2w"), c2r(intra2w), start=start4, count=count4)) 1402 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1w"), c2r(intra1w), start=start4, count=count4)) 1403 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1wS"), c2r(intra1wS), start=start4, count=count4)) 1404 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_chi2tot"), c2r(chi2tot), start=start4, count=count4)) 1405 #endif 1406 ABI_FREE(chi2tot) 1407 end if 1408 1409 if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then 1410 MSG_ERROR(msg) 1411 end if 1412 if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then 1413 MSG_ERROR(msg) 1414 end if 1415 if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then 1416 MSG_ERROR(msg) 1417 end if 1418 if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then 1419 MSG_ERROR(msg) 1420 end if 1421 if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then 1422 MSG_ERROR(msg) 1423 end if 1424 if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then 1425 MSG_ERROR(msg) 1426 end if 1427 if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then 1428 MSG_ERROR(msg) 1429 end if 1430 ! write headers 1431 write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3 1432 write(fout1, '(a,es16.6)' ) ' #tolerance:',tol 1433 write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha' 1434 write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha' 1435 write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 1436 write(fout1, '(a)' )' # Energy Tot-Im Chi(-2w,w,w) Tot-Im Chi(-2w,w,w)' 1437 write(fout1, '(a)' )' # eV *10^-7 esu *10^-12 m/V SI units ' 1438 write(fout1, '(a)' )' # ' 1439 1440 write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3 1441 write(fout2, '(a,es16.6)') ' #tolerance:',tol 1442 write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 1443 write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 1444 write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 1445 write(fout2, '(a)')' # Energy Tot-Re Chi(-2w,w,w) Tot-Re Chi(-2w,w,w)' 1446 write(fout2, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units ' 1447 write(fout2, '(a)')' # ' 1448 1449 write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3 1450 write(fout3, '(a,es16.6)') ' #tolerance:',tol 1451 write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 1452 write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 1453 write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 1454 write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)' 1455 write(fout3, '(a)')' # in esu' 1456 write(fout3, '(a)')' # ' 1457 1458 write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3 1459 write(fout4, '(a,es16.6)') ' #tolerance:',tol 1460 write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 1461 write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 1462 write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 1463 write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)' 1464 write(fout4, '(a)')' # in esu' 1465 write(fout4, '(a)')' # ' 1466 1467 write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3 1468 write(fout5, '(a,es16.6)') ' #tolerance:',tol 1469 write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 1470 write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 1471 write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 1472 write(fout5, '(a)')' # Energy(eV) |TotChi(-2w,w,w)| |Tot Chi(-2w,w,w)|' 1473 write(fout5, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units ' 1474 write(fout5, '(a)')' # ' 1475 1476 write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3 1477 write(fout6, '(a,es16.6)') ' #tolerance:',tol 1478 write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 1479 write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 1480 write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 1481 write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)' 1482 write(fout6, '(a)')' # in esu' 1483 write(fout6, '(a)')' # ' 1484 1485 write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3 1486 write(fout7, '(a,es16.6)') ' #tolerance:',tol 1487 write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 1488 write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 1489 write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 1490 write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)' 1491 write(fout7, '(a)')' # in esu' 1492 write(fout7, '(a)')' # ' 1493 1494 totim=0._dp 1495 totre=0._dp 1496 totabs=0._dp 1497 do iw=2,nmesh 1498 ene=(iw-1)*de 1499 ene=ene*ha2ev 1500 1501 totim=aimag(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7 1502 write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 1503 totim=0._dp 1504 1505 totre=dble(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7 1506 write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 1507 totre=0._dp 1508 1509 write(fout3,'(f15.6,4es15.6)') ene,aimag(inter2w(iw))/1.d-7, & 1510 aimag(inter1w(iw))/1.d-7,aimag(intra2w(iw))/1.d-7, aimag(intra1w(iw)+intra1wS(iw))/1.d-7 1511 1512 write(fout4,'(f15.6,4es15.6)') ene,dble(inter2w(iw))/1.d-7, & 1513 dble(inter1w(iw))/1.d-7,dble(intra2w(iw))/1.d-7,dble(intra1w(iw)+intra1wS(iw))/1.d-7 1514 1515 totabs=abs(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7 1516 write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 1517 totabs=0._dp 1518 1519 write(fout6,'(f15.6,4es15.6)') ene,aimag(inter2w(iw)+inter1w(iw))/1.d-7, & 1520 aimag(intra2w(iw)+intra1w(iw))/1.d-7,aimag(intra1wS(iw))/1.d-7 1521 1522 write(fout7,'(f15.6,4es15.6)') ene,dble(inter2w(iw)+inter1w(iw))/1.d-7, & 1523 dble(intra2w(iw)+intra1w(iw))/1.d-7,dble(intra1wS(iw))/1.d-7 1524 end do 1525 1526 close(fout1) 1527 close(fout2) 1528 close(fout3) 1529 close(fout4) 1530 close(fout5) 1531 close(fout6) 1532 close(fout7) 1533 ! print information 1534 write(std_out,*) ' ' 1535 write(std_out,*) 'information about calculation just performed:' 1536 write(std_out,*) ' ' 1537 write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of second order susceptibility' 1538 write(std_out,*) 'tolerance:',tol 1539 if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high' 1540 write(std_out,*) 'broadening:',brod,'Hartree' 1541 if (brod.gt.0.009) then 1542 write(std_out,*) ' ' 1543 write(std_out,*) 'ATTENTION: broadening is quite high' 1544 write(std_out,*) ' ' 1545 else if (brod.gt.0.015) then 1546 write(std_out,*) ' ' 1547 write(std_out,*) 'ATTENTION: broadening is too high' 1548 write(std_out,*) ' ' 1549 end if 1550 write(std_out,*) 'scissors shift:',sc,'Hartree' 1551 write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree' 1552 end if 1553 1554 ! deallocate local arrays 1555 ABI_DEALLOCATE(px) 1556 ABI_DEALLOCATE(py) 1557 ABI_DEALLOCATE(pz) 1558 ABI_DEALLOCATE(inter2w) 1559 ABI_DEALLOCATE(inter1w) 1560 ABI_DEALLOCATE(intra2w) 1561 ABI_DEALLOCATE(intra1w) 1562 ABI_DEALLOCATE(intra1wS) 1563 ABI_DEALLOCATE(delta) 1564 1565 end subroutine nlinopt
m_optic_tools/nonlinopt [ Functions ]
[ Top ] [ m_optic_tools ] [ Functions ]
NAME
nonlinopt
FUNCTION
Compute optical frequency dependent linear electro-optic susceptibility for semiconductors
INPUTS
icomp=Sequential index associated to computed tensor components (used for netcdf output) itemp=Temperature index (used for netcdf output) nspin = number of spins(integer) omega = crystal volume in au (real) nkpt = total number of kpoints (integer) wkpt(nkpt) = weights of kpoints (real) nsymcrys = number of crystal symmetry operations(integer) symcrys(3,3,nsymcrys) = symmetry operations in cartisian coordinates(real) nstval = total number of valence states(integer) evalv(nstval,nspin,nkpt) = eigen value for each band in Ha(real) occv(nstval,nspin,nkpt) = occupation number efermi = Fermi energy in Ha(real) pmat(nstval,nstval,nkpt,3,nspin) = momentum matrix elements in cartesian coordinates(complex) v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z nmesh = desired number of energy mesh points(integer) de = desired step in energy(real); nmesh*de=maximum energy for plotting sc = scissors shift in Ha(real) brod = broadening in Ha(real) tol = tolerance:how close to the singularity exact exact is calculated(real) fnam=root for filenames that will contain the output : fnam1=trim(fnam)//'-ChiTotIm.out' fnam2=trim(fnam)//'-ChiTotRe.out' fnam3=trim(fnam)//'-ChiIm.out' fnam4=trim(fnam)//'-ChiRe.out' fnam5=trim(fnam)//'-ChiAbs.out'
OUTPUT
Calculates the second harmonic generation susceptibility on a desired energy mesh and for desired direction of polarisation. The output is in files named ChiEOTot.out : Im\chi_{v1v2v3}(\omega,\omega,0) and Re\chi_{v1v2v3}(\omega,\omega,0) ChiEOIm.out : contributions to the Im\chi_{v1v2v3}(\omega,\omega,0) from various terms ChiEORe.out : contributions to Re\chi_{v1v2v3}(\omega,\omega,-0) from various terms ChiEOAbs.out : abs\chi_{v1v2v3}(\omega,\omega,0). The headers in these files contain information about the calculation. ncid=Netcdf id to save output data. COMMENTS - The routine has been written using notations of Ref. 2 - This routine does not symmetrize the tensor (up to now) - Sum over all the states and use occupation factors instead of looping only on resonant contributions
PARENTS
optic
CHILDREN
xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum
SOURCE
2209 subroutine nonlinopt(icomp,itemp,nspin,omega,nkpt,wkpt,nsymcrys,symcrys,nstval,evalv,occv,efermi, & 2210 pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,do_antiresonant,ncid,comm) 2211 2212 2213 !This section has been created automatically by the script Abilint (TD). 2214 !Do not modify the following lines by hand. 2215 #undef ABI_FUNC 2216 #define ABI_FUNC 'nonlinopt' 2217 !End of the abilint section 2218 2219 implicit none 2220 2221 !Arguments ------------------------------------ 2222 !no_abirules 2223 integer, intent(in) :: icomp,itemp,nspin, ncid 2224 real(dp), intent(in) :: omega 2225 integer, intent(in) :: nkpt 2226 real(dp), intent(in) :: wkpt(nkpt) 2227 integer, intent(in) :: nsymcrys 2228 real(dp), intent(in) :: symcrys(3,3,nsymcrys) 2229 integer, intent(in) :: nstval 2230 real(dp), intent(in) :: evalv(nstval,nkpt,nspin) 2231 real(dp), intent(in) :: occv(nstval,nkpt,nspin) 2232 real(dp), intent(in) :: efermi 2233 complex(dpc), intent(in) :: pmat(nstval,nstval,nkpt,3,nspin) 2234 integer, intent(in) :: v1 2235 integer, intent(in) :: v2 2236 integer, intent(in) :: v3 2237 integer, intent(in) :: nmesh 2238 integer, intent(in) :: comm 2239 real(dp), intent(in) :: de 2240 real(dp), intent(in) :: sc 2241 real(dp), intent(in) :: brod 2242 real(dp), intent(in) :: tol 2243 character(len=*), intent(in) :: fnam 2244 logical, intent(in) :: do_antiresonant 2245 2246 !Local variables ------------------------- 2247 integer :: iw 2248 integer :: i,j,k,lx,ly,lz 2249 integer :: isp,isym,ik 2250 integer :: ist1,istl,istn,istm 2251 real(dp) :: ha2ev 2252 real(dp) :: t1,t2,t3,tst 2253 real(dp) :: ene,totre,totabs,totim 2254 real(dp) :: el,en,em 2255 real(dp) :: emin,emax, my_emin,my_emax 2256 real(dp) :: const_esu,const_au,au2esu 2257 real(dp) :: wmn,wnm,wln,wnl,wml,wlm 2258 complex(dpc) :: idel,w,zi 2259 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7 2260 ! local allocatable arrays 2261 integer :: start4(4),count4(4) 2262 real(dp) :: s(3,3),sym(3,3,3) 2263 integer :: istp 2264 real(dp) :: ep, wmp, wpn 2265 real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included ! 2266 real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, flm 2267 complex(dpc), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a} 2268 complex(dpc), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a} 2269 complex(dpc), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k) 2270 complex(dpc), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c 2271 complex(dpc), allocatable :: chiw(:), chi2w(:) ! \chi_{II}^{abc}(-\omega,\omega,0) 2272 complex(dpc), allocatable :: etaw(:), eta2w(:) ! \eta_{II}^{abc}(-\omega,\omega,0) 2273 complex(dpc), allocatable :: sigmaw(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0) 2274 complex(dpc) :: num1, num2, den1, den2, term1, term2 2275 complex(dpc) :: chi1, chi2_1, chi2_2 2276 complex(dpc), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega) 2277 complex(dpc), allocatable :: eta1(:) ! Second term that depends on the frequency ! (omega) 2278 complex(dpc), allocatable :: chi2tot(:) 2279 complex(dpc) :: eta1_1, eta1_2, eta2_1, eta2_2 2280 complex(dpc) :: sigma2_1, sigma1 2281 complex(dpc), allocatable :: symrmn(:,:,:) ! (m,l,n) = 1/2*(rml^b rln^c+rml^c rln^b) 2282 complex(dpc) :: symrmnl(3,3), symrlmn(3,3), symrmln(3,3) 2283 !Parallelism 2284 integer :: my_rank, nproc 2285 integer,parameter :: master=0 2286 integer :: ierr 2287 integer :: my_k1, my_k2 2288 character(500) :: msg 2289 integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7 2290 2291 ! ********************************************************************* 2292 2293 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 2294 2295 !calculate the constant 2296 zi=(0._dp,1._dp) 2297 idel=zi*brod 2298 const_au=-2._dp/(omega*dble(nsymcrys)) 2299 !const_au=-2._dp/(omega) 2300 au2esu=5.8300348177d-8 2301 const_esu=const_au*au2esu 2302 ha2ev=13.60569172*2._dp 2303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2304 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev 2305 !bohr: 5.2917ifc nlinopt.f907E-11 2306 !c: 2.99792458 velocity of sound 2307 !ry2ev: 13.60569172 2308 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2) 2309 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3) 2310 !mass comes from converting P_mn to r_mn 2311 !hbar^3 comes from converting all frequencies to energies in denominator 2312 !hbar^3 comes from operator for momentum (hbar/i nabla) 2313 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2314 !output file names 2315 fnam1=trim(fnam)//'-ChiSHGTotIm.out' 2316 fnam2=trim(fnam)//'-ChiSHGTotRe.out' 2317 fnam3=trim(fnam)//'-ChiSHGIm.out' 2318 fnam4=trim(fnam)//'-ChiSHGRe.out' 2319 fnam5=trim(fnam)//'-ChiSHGAbs.out' 2320 fnam6=trim(fnam)//'-ChiSHGImDec.out' 2321 fnam7=trim(fnam)//'-ChiSHGReDec.out' 2322 2323 !If there exists inversion symmetry exit with a mesg. 2324 tst=1.d-09 2325 do isym=1,nsymcrys 2326 t1=symcrys(1,1,isym)+1 2327 t2=symcrys(2,2,isym)+1 2328 t3=symcrys(3,3,isym)+1 2329 ! test if diagonal elements are -1 2330 if (abs(t1).lt.tst.and.abs(t2).lt.tst.and.abs(t3).lt.tst) then 2331 ! test if off-diagonal elements are zero 2332 if (abs(symcrys(1,2,isym)).lt.tst.and.abs(symcrys(1,3,isym)).lt.tst & 2333 .and.abs(symcrys(2,1,isym)).lt.tst.and.abs(symcrys(2,3,isym)).lt.tst.and. & 2334 abs(symcrys(3,1,isym)).lt.tst.and.abs(symcrys(3,2,isym)).lt.tst) then 2335 write(std_out,*) '-----------------------------------------' 2336 write(std_out,*) ' the crystal has inversion symmetry ' 2337 write(std_out,*) ' the SHG susceptibility is zero ' 2338 write(std_out,*) '-----------------------------------------' 2339 MSG_ERROR("Aborting now") 2340 end if 2341 end if 2342 end do 2343 2344 !check polarisation 2345 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 2346 write(std_out,*) '---------------------------------------------' 2347 write(std_out,*) ' Error in linelop: ' 2348 write(std_out,*) ' the polarisation directions incorrect ' 2349 write(std_out,*) ' 1=x, 2=y and 3=z ' 2350 write(std_out,*) '---------------------------------------------' 2351 MSG_ERROR("Aborting now") 2352 end if 2353 2354 !number of energy mesh points 2355 if (nmesh.le.0) then 2356 write(std_out,*) '---------------------------------------------' 2357 write(std_out,*) ' Error in linelop: ' 2358 write(std_out,*) ' number of energy mesh points incorrect ' 2359 write(std_out,*) ' number has to be integer greater than 0 ' 2360 write(std_out,*) ' nmesh*de = max energy for calculation ' 2361 write(std_out,*) '---------------------------------------------' 2362 MSG_ERROR("Aborting now") 2363 end if 2364 2365 !step in energy 2366 if (de.le.0._dp) then 2367 write(std_out,*) '---------------------------------------------' 2368 write(std_out,*) ' Error in nlinopt: ' 2369 write(std_out,*) ' energy step is incorrect ' 2370 write(std_out,*) ' number has to real greater than 0.0 ' 2371 write(std_out,*) ' nmesh*de = max energy for calculation ' 2372 write(std_out,*) '---------------------------------------------' 2373 MSG_ERROR("Aborting now") 2374 end if 2375 2376 !broadening 2377 if (brod.gt.0.009) then 2378 write(std_out,*) '---------------------------------------------' 2379 write(std_out,*) ' ATTENTION: broadening is quite high ' 2380 write(std_out,*) ' ideally should be less than 0.005 ' 2381 write(std_out,*) '---------------------------------------------' 2382 else if (brod.gt.0.015) then 2383 write(std_out,*) '----------------------------------------' 2384 write(std_out,*) ' ATTENTION: broadening is too high ' 2385 write(std_out,*) ' ideally should be less than 0.005 ' 2386 write(std_out,*) '----------------------------------------' 2387 end if 2388 2389 !tolerance 2390 if (tol.gt.0.006) then 2391 write(std_out,*) '----------------------------------------' 2392 write(std_out,*) ' ATTENTION: tolerance is too high ' 2393 write(std_out,*) ' ideally should be less than 0.004 ' 2394 write(std_out,*) '----------------------------------------' 2395 end if 2396 2397 ! allocate local arrays 2398 ABI_MALLOC(enk,(nstval)) 2399 ABI_MALLOC(delta,(nstval,nstval,3)) 2400 ABI_MALLOC(rmnbc,(nstval,nstval,3,3)) 2401 ABI_MALLOC(roverw,(nstval,nstval,3,3)) 2402 ABI_MALLOC(rmna,(nstval,nstval,3)) 2403 ABI_MALLOC(chiw,(nmesh)) 2404 ABI_MALLOC(etaw,(nmesh)) 2405 ABI_MALLOC(chi2w,(nmesh)) 2406 ABI_MALLOC(eta2w,(nmesh)) 2407 ABI_MALLOC(sigmaw,(nmesh)) 2408 ABI_MALLOC(chi2,(nmesh)) 2409 ABI_MALLOC(eta1,(nmesh)) 2410 ABI_MALLOC(symrmn,(nstval,nstval,nstval)) 2411 2412 !generate the symmetrizing tensor 2413 sym(:,:,:)=0._dp 2414 do isym=1,nsymcrys 2415 s(:,:)=symcrys(:,:,isym) 2416 do i=1,3 2417 do j=1,3 2418 do k=1,3 2419 sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3)) 2420 end do 2421 end do 2422 end do 2423 end do 2424 2425 2426 !initialise 2427 delta(:,:,:)=0._dp 2428 rmnbc(:,:,:,:)=0._dp 2429 chiw(:)=0._dp 2430 chi2w(:)=0._dp 2431 chi2(:) = 0._dp 2432 etaw(:)=0._dp 2433 eta2w(:)=0._dp 2434 sigmaw(:)=0._dp 2435 my_emin=HUGE(0._dp) 2436 my_emax=-HUGE(0._dp) 2437 2438 ! Split work 2439 call xmpi_split_work(nkpt,comm,my_k1,my_k2,msg,ierr) 2440 2441 ! loop over kpts 2442 do ik=my_k1,my_k2 2443 write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt 2444 do isp=1,nspin 2445 ! Calculate the scissor corrected energies and the energy window 2446 do ist1=1,nstval 2447 en = evalv(ist1,ik,isp) 2448 my_emin=min(my_emin,en) 2449 my_emax=max(my_emax,en) 2450 if(en > efermi) then 2451 en = en + sc 2452 end if 2453 enk(ist1) = en 2454 end do 2455 2456 ! calculate \Delta_nm and r_mn^a 2457 do istn=1,nstval 2458 en = enk(istn) 2459 do istm=1,nstval 2460 em = enk(istm) 2461 wmn = em - en 2462 delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp) 2463 if(abs(wmn) < tol) then 2464 rmna(istm,istn,1:3) = 0._dp 2465 else 2466 rmna(istm,istn,1:3)=pmat(istm,istn,ik,1:3,isp)/wmn 2467 end if 2468 end do 2469 end do 2470 ! calculate \r^b_mn;c 2471 do istm=1,nstval 2472 em = enk(istm) 2473 do istn=1,nstval 2474 en = enk(istn) 2475 wmn = em - en 2476 if (abs(wmn) < tol) then ! Degenerate energies 2477 rmnbc(istm,istn,:,:) = 0.0 2478 roverw(istm,istn,:,:) = 0.0 2479 cycle 2480 end if 2481 do ly = 1,3 2482 do lz = 1,3 2483 num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly)) 2484 den1 = wmn 2485 term1 = num1/den1 2486 term2 = 0._dp 2487 do istp=1,nstval 2488 ep = enk(istp) 2489 wmp = em - ep 2490 wpn = ep - en 2491 num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly)) 2492 den2 = wmn 2493 term2 = term2 + (num2/den2) 2494 end do 2495 rmnbc(istm,istn,ly,lz) = -term1-(zi*term2) 2496 roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz) 2497 end do 2498 end do 2499 end do 2500 end do 2501 2502 ! initialise the factors 2503 ! start the calculation 2504 do istn=1,nstval 2505 en=enk(istn) 2506 fn=occv(istn,ik,isp) 2507 if(do_antiresonant .and. en .ge. efermi) then 2508 cycle 2509 end if 2510 do istm=1,nstval 2511 em=enk(istm) 2512 if (do_antiresonant .and. em .le. efermi) then 2513 cycle 2514 end if 2515 wmn=em-en 2516 wnm=-wmn 2517 fm = occv(istm,ik,isp) 2518 fnm = fn - fm 2519 if(abs(wmn) > tol) then 2520 chi1 = 0._dp 2521 chi2(:) = 0._dp 2522 chi2_1 = 0._dp 2523 chi2_2 = 0._dp 2524 eta1(:) = 0._dp 2525 eta1_1 = 0._dp 2526 eta1_2 = 0._dp 2527 eta2_1 = 0._dp 2528 eta2_2 = 0._dp 2529 sigma1 = 0._dp 2530 sigma2_1 = 0._dp 2531 ! Three band terms 2532 do istl=1,nstval 2533 el=enk(istl) 2534 fl = occv(istl,ik,isp) 2535 wlm = el-em 2536 wln = el-en 2537 wnl = -wln 2538 wml = em-el 2539 fnl = fn-fl 2540 fml = fm-fl 2541 flm = -fml 2542 fln = -fnl 2543 do ly = 1,3 2544 do lz = 1,3 2545 symrmnl(ly,lz) = 0.5_dp*(rmna(istm,istn,ly)*rmna(istn,istl,lz)+rmna(istm,istn,lz)*rmna(istn,istl,ly)) 2546 symrlmn(ly,lz) = 0.5_dp*(rmna(istl,istm,ly)*rmna(istm,istn,lz)+rmna(istl,istm,lz)*rmna(istm,istn,ly)) 2547 symrmln(ly,lz) = 0.5_dp*(rmna(istm,istl,ly)*rmna(istl,istn,lz)+rmna(istm,istl,lz)*rmna(istl,istn,ly)) 2548 end do 2549 end do 2550 2551 do lx = 1,3 2552 do ly = 1,3 2553 do lz = 1,3 2554 sigma1 = sigma1 + sym(lx,ly,lz)*(wnl*rmna(istl,istm,lx)*symrmnl(ly,lz)-wlm*rmna(istn,istl,lx)*symrlmn(ly,lz)) 2555 eta2_2 = eta2_2 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*symrmln(ly,lz)*(wml-wln) 2556 if(abs(wln-wml) > tol) then 2557 chi1 = chi1 + sym(lx,ly,lz)*(rmna(istn,istm,lx)*symrmln(ly,lz))/(wln-wml) 2558 end if 2559 eta1_1 = eta1_1 + sym(lx,ly,lz)*wln*rmna(istn,istl,lx)*symrlmn(ly,lz) 2560 eta1_2 = eta1_2 - sym(lx,ly,lz)*wml*rmna(istl,istm,lx)*symrmnl(ly,lz) 2561 if(abs(wnl-wmn) > tol) then 2562 chi2_1 = chi2_1 - sym(lx,ly,lz)*(fnm*rmna(istl,istm,lx)*symrmnl(ly,lz)/(wnl-wmn)) 2563 end if 2564 if(abs(wmn-wlm) > tol) then 2565 chi2_2 = chi2_2 - sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*symrlmn(ly,lz)/(wmn-wlm)) 2566 end if 2567 end do 2568 end do 2569 end do 2570 end do 2571 2572 ! Two band terms 2573 eta2_1 = 0.0_dp 2574 sigma2_1 = 0.0_dp 2575 do lx = 1,3 2576 do ly = 1,3 2577 do lz = 1,3 2578 eta2_1 = eta2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp & 2579 & *(delta(istm,istn,ly)*rmna(istm,istn,lz)+delta(istm,istn,lz)*rmna(istm,istn,ly)) 2580 ! Correct version (Sipe 1993) 2581 sigma2_1 = sigma2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp & 2582 & *(rmna(istm,istn,ly)*delta(istn,istm,lz)+rmna(istm,istn,lz)*delta(istn,istm,ly)) 2583 2584 ! Incorrect version (Hughes 1996) 2585 !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)) 2586 end do 2587 end do 2588 end do 2589 ! 2590 ! calculate over the desired energy mesh and sum over k-points 2591 do iw=1,nmesh 2592 w=(iw-1)*de+idel 2593 chi2w(iw) = chi2w(iw) + zi*wkpt(ik)*((2.0_dp*fnm*chi1/(wmn-2.0_dp*w)))*const_esu ! Inter(2w) from chi 2594 chiw(iw) = chiw(iw) + zi*wkpt(ik)*((chi2_1+chi2_2)/(wmn-w))*const_esu ! Inter(w) from chi 2595 eta2w(iw) = eta2w(iw) + zi*wkpt(ik)*(8.0_dp*(eta2_1/((wmn**2)*(wmn-2.0_dp*w))) & 2596 & + 2.0_dp*eta2_2/((wmn**2)*(wmn-2.0_dp*w)))*const_esu ! Intra(2w) from eta 2597 etaw(iw) = etaw(iw) + zi*wkpt(ik)*((eta1_1 + eta1_2)*fnm/((wmn**2)*(wmn-w)))*const_esu ! Intra(w) from eta 2598 sigmaw(iw) = sigmaw(iw) + 0.5_dp*zi*wkpt(ik)*(fnm*sigma1/((wmn**2)*(wmn-w)) & 2599 & + (sigma2_1/((wmn**2)*(wmn-w))))*const_esu ! Intra(1w) from sigma 2600 end do 2601 end if 2602 end do ! end loop over istn and istm 2603 end do 2604 end do ! spins 2605 end do ! k-points 2606 2607 ! Collect info among the nodes 2608 call xmpi_min(my_emin,emin,comm,ierr) 2609 call xmpi_max(my_emax,emax,comm,ierr) 2610 2611 call xmpi_sum(chiw,comm,ierr) 2612 call xmpi_sum(etaw,comm,ierr) 2613 call xmpi_sum(chi2w,comm,ierr) 2614 call xmpi_sum(eta2w,comm,ierr) 2615 call xmpi_sum(sigmaw,comm,ierr) 2616 2617 ! Master writes the output 2618 if (my_rank == master) then 2619 2620 if (ncid /= nctk_noid) then 2621 start4 = [1, 1, icomp, itemp] 2622 count4 = [2, nmesh, 1, 1] 2623 ABI_MALLOC(chi2tot, (nmesh)) 2624 chi2tot = chiw + chi2w + etaw + eta2w + sigmaw 2625 #ifdef HAVE_NETCDF 2626 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2tot"), c2r(chi2tot), start=start4, count=count4)) 2627 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chiw"), c2r(chiw), start=start4, count=count4)) 2628 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_etaw"), c2r(etaw), start=start4, count=count4)) 2629 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2w"), c2r(chi2w), start=start4, count=count4)) 2630 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_eta2w"), c2r(eta2w), start=start4, count=count4)) 2631 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_sigmaw"), c2r(sigmaw), start=start4, count=count4)) 2632 #endif 2633 ABI_FREE(chi2tot) 2634 end if 2635 2636 ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000) 2637 if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then 2638 MSG_ERROR(msg) 2639 end if 2640 if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then 2641 MSG_ERROR(msg) 2642 end if 2643 if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then 2644 MSG_ERROR(msg) 2645 end if 2646 if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then 2647 MSG_ERROR(msg) 2648 end if 2649 if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then 2650 MSG_ERROR(msg) 2651 end if 2652 if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then 2653 MSG_ERROR(msg) 2654 end if 2655 if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then 2656 MSG_ERROR(msg) 2657 end if 2658 !!write headers 2659 write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3 2660 write(fout1, '(a,es16.6)' ) ' #tolerance:',tol 2661 write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha' 2662 write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha' 2663 write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2664 write(fout1, '(a)' )' # Energy Tot-Im Chi(-w,w,0) Tot-Im Chi(-w,w,0)' 2665 write(fout1, '(a)' )' # eV *10^-7 esu *10^-12 m/V SI units ' 2666 write(fout1, '(a)' )' # ' 2667 2668 write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3 2669 write(fout2, '(a,es16.6)') ' #tolerance:',tol 2670 write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2671 write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2672 write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2673 write(fout2, '(a)')' # Energy Tot-Re Chi(-w,w,0) Tot-Re Chi(-w,w,0)' 2674 write(fout2, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units ' 2675 write(fout2, '(a)')' # ' 2676 2677 write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2678 write(fout3, '(a,es16.6)') ' #tolerance:',tol 2679 write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2680 write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2681 write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2682 write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)' 2683 write(fout3, '(a)')' # in esu' 2684 write(fout3, '(a)')' # ' 2685 2686 write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2687 write(fout4, '(a,es16.6)') ' #tolerance:',tol 2688 write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2689 write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2690 write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2691 write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)' 2692 write(fout4, '(a)')' # in esu' 2693 write(fout4, '(a)')' # ' 2694 2695 write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2696 write(fout5, '(a,es16.6)') ' #tolerance:',tol 2697 write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2698 write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2699 write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2700 write(fout5, '(a)')' # Energy(eV) |TotChi(-w,w,0)| |Tot Chi(-w,w,0)|' 2701 write(fout5, '(a)')' # eV *10^-7 esu *10^-12 m/V SI units ' 2702 write(fout5, '(a)')' # ' 2703 2704 write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2705 write(fout6, '(a,es16.6)') ' #tolerance:',tol 2706 write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2707 write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2708 write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2709 write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)' 2710 write(fout6, '(a)')' # in esu' 2711 write(fout6, '(a)')' # ' 2712 2713 write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3 2714 write(fout7, '(a,es16.6)') ' #tolerance:',tol 2715 write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha' 2716 write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha' 2717 write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha' 2718 write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)' 2719 write(fout7, '(a)')' # in esu' 2720 write(fout7, '(a)')' # ' 2721 2722 totim=0._dp 2723 totre=0._dp 2724 totabs=0._dp 2725 do iw=2,nmesh 2726 ene=(iw-1)*de 2727 ene=ene*ha2ev 2728 2729 totim=aimag(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7 2730 write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 2731 totim=0._dp 2732 2733 totre=dble(chiw(iw)+chi2w(iw)+eta2w(iw)+etaw(iw)+sigmaw(iw))/1.d-7 2734 write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 2735 totre=0._dp 2736 2737 write(fout3,'(f15.6,4es15.6)') ene,aimag(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7, & 2738 aimag(eta2w(iw))/1.d-7,aimag(etaw(iw))/1.d-7+aimag(sigmaw(iw))/1.d-7 2739 2740 write(fout4,'(f15.6,4es15.6)') ene,dble(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7, & 2741 dble(eta2w(iw))/1.d-7,dble(etaw(iw))/1.d-7+dble(sigmaw(iw))/1.d-7 2742 2743 totabs=abs(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7 2744 write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5) 2745 totabs=0._dp 2746 2747 write(fout6,'(f15.6,4es15.6)') ene,aimag(chi2w(iw)+chiw(iw))/1.d-7, & 2748 aimag(eta2w(iw)+etaw(iw))/1.d-7,aimag(sigmaw(iw))/1.d-7 2749 2750 write(fout7,'(f15.6,4es15.6)') ene,dble(chi2w(iw)+chiw(iw))/1.d-7, & 2751 dble(eta2w(iw)+etaw(iw))/1.d-7,dble(sigmaw(iw))/1.d-7 2752 end do 2753 2754 close(fout1) 2755 close(fout2) 2756 close(fout3) 2757 close(fout4) 2758 close(fout5) 2759 close(fout6) 2760 close(fout7) 2761 2762 ! print information 2763 write(std_out,*) ' ' 2764 write(std_out,*) 'information about calculation just performed:' 2765 write(std_out,*) ' ' 2766 write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of second order susceptibility' 2767 write(std_out,*) 'tolerance:',tol 2768 if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high' 2769 write(std_out,*) 'broadening:',brod,'Hartree' 2770 if (brod.gt.0.009) then 2771 write(std_out,*) ' ' 2772 write(std_out,*) 'ATTENTION: broadening is quite high' 2773 write(std_out,*) ' ' 2774 else if (brod.gt.0.015) then 2775 write(std_out,*) ' ' 2776 write(std_out,*) 'ATTENTION: broadening is too high' 2777 write(std_out,*) ' ' 2778 end if 2779 write(std_out,*) 'scissors shift:',sc,'Hartree' 2780 write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree' 2781 2782 end if 2783 2784 ! deallocate local arrays 2785 ABI_FREE(enk) 2786 ABI_FREE(delta) 2787 ABI_FREE(rmnbc) 2788 ABI_FREE(roverw) 2789 ABI_FREE(rmna) 2790 ABI_FREE(chiw) 2791 ABI_FREE(chi2w) 2792 ABI_FREE(chi2) 2793 ABI_FREE(etaw) 2794 ABI_FREE(eta1) 2795 ABI_FREE(symrmn) 2796 ABI_FREE(eta2w) 2797 ABI_FREE(sigmaw) 2798 2799 end subroutine nonlinopt
m_optic_tools/pmat2cart [ Functions ]
[ Top ] [ m_optic_tools ] [ Functions ]
NAME
pmat2cart
FUNCTION
turn momentum matrix elements to cartesian axes. To be used in optic calculation of linear and non-linear RPA dielectric matrices
INPUTS
eigen11,eigen12,eigen13 = first order ddk eigen values = d eig_i,k / dk for 3 reduced directions mband=maximum number of bands nkpt = number of k-points nsppol=1 for unpolarized, 2 for spin-polarized rprimd(3,3)=dimensional real space primitive translations (bohr)
OUTPUT
pmat(mband,mband,nkpt,3,nsppol) = matrix elements of momentum operator, in cartesian coordinates
PARENTS
optic
CHILDREN
xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum
SOURCE
274 subroutine pmat2cart(eigen11,eigen12,eigen13,mband,nkpt,nsppol,pmat,rprimd) 275 276 277 !This section has been created automatically by the script Abilint (TD). 278 !Do not modify the following lines by hand. 279 #undef ABI_FUNC 280 #define ABI_FUNC 'pmat2cart' 281 !End of the abilint section 282 283 implicit none 284 285 !Arguments ----------------------------------------------- 286 !scalars 287 integer,intent(in) :: mband,nkpt,nsppol 288 !arrays 289 real(dp),intent(in) :: eigen11(2,mband,mband,nkpt,nsppol) 290 real(dp),intent(in) :: eigen12(2,mband,mband,nkpt,nsppol) 291 real(dp),intent(in) :: eigen13(2,mband,mband,nkpt,nsppol),rprimd(3,3) 292 !no_abirules 293 complex(dpc),intent(out) :: pmat(mband,mband,nkpt,3,nsppol) 294 295 !Local variables ----------------------------------------- 296 !scalars 297 integer :: iband1,iband2,ikpt,isppol 298 !arrays 299 real(dp) :: rprim(3,3) 300 301 ! ************************************************************************* 302 303 !rescale the rprim 304 rprim(:,:) = rprimd(:,:) / two_pi 305 306 do isppol=1,nsppol 307 do ikpt=1,nkpt 308 do iband1=1,mband 309 do iband2=1,mband 310 pmat(iband2,iband1,ikpt,:,isppol) = & 311 & rprim(:,1)*cmplx(eigen11(1,iband2,iband1,ikpt,isppol),eigen11(2,iband2,iband1,ikpt,isppol),kind=dp) & 312 & +rprim(:,2)*cmplx(eigen12(1,iband2,iband1,ikpt,isppol),eigen12(2,iband2,iband1,ikpt,isppol),kind=dp) & 313 & +rprim(:,3)*cmplx(eigen13(1,iband2,iband1,ikpt,isppol),eigen13(2,iband2,iband1,ikpt,isppol),kind=dp) 314 end do 315 end do 316 end do 317 end do 318 319 end subroutine pmat2cart
m_optic_tools/pmat_renorm [ Functions ]
[ Top ] [ m_optic_tools ] [ Functions ]
NAME
pmat_renorm
FUNCTION
Renormalize the momentum matrix elements according to the scissor shift which is imposed
INPUTS
mband= number of bands nkpt = number of k-points nsppol=1 for unpolarized, 2 for spin-polarized efermi = Fermi level sc = scissor shift for conduction bands evalv = eigenvalues for ground state
OUTPUT
pmat(mband,mband,nkpt,3,nsppol) = momentum matrix elements, renormalized by denominator change with scissor shift
PARENTS
optic
CHILDREN
xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum
SOURCE
350 subroutine pmat_renorm(efermi, evalv, mband, nkpt, nsppol, pmat, sc) 351 352 353 !This section has been created automatically by the script Abilint (TD). 354 !Do not modify the following lines by hand. 355 #undef ABI_FUNC 356 #define ABI_FUNC 'pmat_renorm' 357 use interfaces_14_hidewrite 358 !End of the abilint section 359 360 implicit none 361 362 !Arguments ----------------------------------------------- 363 !scalars 364 integer, intent(in) :: nsppol 365 integer, intent(in) :: nkpt 366 integer, intent(in) :: mband 367 real(dp), intent(in) :: efermi 368 real(dp), intent(in) :: sc 369 !arrays 370 real(dp), intent(in) :: evalv(mband,nkpt,nsppol) 371 complex(dpc), intent(inout) :: pmat(mband,mband,nkpt,3,nsppol) 372 373 !Local variables ----------------------------------------- 374 !scalars 375 integer :: iband1,iband2,ikpt,isppol 376 real(dp) :: corec, e1, e2 377 378 ! ************************************************************************* 379 380 if (abs(sc) < tol8) then 381 call wrtout(std_out,' No scissor shift to be applied. Returning to main optic routine.',"COLL") 382 return 383 end if 384 385 do isppol=1,nsppol 386 do ikpt=1,nkpt 387 do iband1=1,mband ! valence states 388 e1 = evalv(iband1,ikpt,isppol) 389 if (e1 > efermi) cycle 390 do iband2=1,mband ! conduction states 391 e2 = evalv(iband2,ikpt,isppol) 392 if (e2 < efermi) cycle 393 corec = (e2+sc-e1)/(e2-e1) 394 pmat(iband2,iband1,ikpt,:,isppol) = corec * pmat(iband2,iband1,ikpt,:,isppol) 395 pmat(iband1,iband2,ikpt,:,isppol) = corec * pmat(iband1,iband2,ikpt,:,isppol) 396 end do 397 end do 398 end do 399 end do 400 401 end subroutine pmat_renorm
m_optic_tools/sym2cart [ Functions ]
[ Top ] [ m_optic_tools ] [ Functions ]
NAME
sym2cart
FUNCTION
Routine called by the program optic Convert to symmetry matrice in cartesian coordinates
INPUTS
gprimd(3,3)=dimensional primitive translations for reciprocal space nsym=number of symmetries in group rprimd(3,3)=dimensional real space primitive translations (bohr) symrel(3,3,nsym)=symmetry matrices in terms of real space
OUTPUT
symcart(3,3)=symmetry matrice in cartesian coordinates (reals)
PARENTS
optic
CHILDREN
xmpi_max,xmpi_min,xmpi_split_work,xmpi_sum
SOURCE
91 subroutine sym2cart(gprimd,nsym,rprimd,symrel,symcart) 92 93 94 !This section has been created automatically by the script Abilint (TD). 95 !Do not modify the following lines by hand. 96 #undef ABI_FUNC 97 #define ABI_FUNC 'sym2cart' 98 !End of the abilint section 99 100 implicit none 101 102 !Arguments ----------------------------------------------- 103 ! in 104 ! out 105 !scalars 106 integer,intent(in) :: nsym 107 !arrays 108 integer,intent(in) :: symrel(3,3,nsym) 109 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 110 real(dp),intent(out) :: symcart(3,3,nsym) 111 112 !Local variables------------------------------- 113 !scalars 114 integer :: isym 115 !arrays 116 real(dp) :: rsym(3,3),rsymcart(3,3),tmp(3,3) 117 118 ! ************************************************************************* 119 120 do isym=1,nsym 121 rsym(:,:) = dble(symrel(:,:,isym)) 122 ! write(std_out,*) 'rsym = ',rsym 123 call dgemm('N','N',3,3,3,one,rprimd,3,rsym, 3,zero,tmp, 3) 124 call dgemm('N','N',3,3,3,one,tmp, 3,gprimd,3,zero,rsymcart,3) 125 ! write(std_out,*) 'rsymcart = ',rsymcart 126 symcart(:,:,isym) = rsymcart(:,:) 127 ! purify symops in cartesian dp coordinates 128 where( abs(symcart(:,:,isym))<tol14) 129 symcart(:,:,isym) = zero 130 end where 131 end do 132 133 end subroutine sym2cart