TABLE OF CONTENTS


ABINIT/m_optic_tools [ Modules ]

[ Top ] [ Modules ]

NAME

 m_optic_tools

FUNCTION

  Helper functions used in the optic code

COPYRIGHT

 Copyright (C) 2002-2022 ABINIT group (SSharma,MVer,VRecoules,TD,YG, NAP)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

 COMMENTS

  Right now the routine sums over the k-points. In future linear tetrahedron method might be useful.

  Reference articles:

      1. S. Sharma, J. K. Dewhurst and C. Ambrosch-Draxl, Phys. Rev. B {\bf 67} 165332 2003 [[cite:Sharma2003]]
      2. J. L. P. Hughes and J. E. Sipe, Phys. Rev. B {\bf 53} 10 751 1996 [[cite:Hughes1996]]
      3. S. Sharma and C. Ambrosch-Draxl, Physica Scripta T 109 2004 [[cite:Sharma2004]]
      4. J. E. Sipe and Ed. Ghahramani, Phys. Rev. B {\bf 48} 11 705 1993 [[cite:Sipe1993]]

SOURCE

28 #if defined HAVE_CONFIG_H
29 #include "config.h"
30 #endif
31 
32 #include "abi_common.h"
33 
34 module m_optic_tools
35 
36  use defs_basis
37  use m_errors
38  use m_abicore
39  use m_linalg_interfaces
40  use m_xmpi
41  use m_nctk
42 #ifdef HAVE_NETCDF
43  use netcdf
44 #endif
45 
46  use defs_datatypes,    only : ebands_t
47  use m_numeric_tools,   only : c2r
48  use m_io_tools,        only : open_file
49  use m_crystal,         only : crystal_t
50 
51  implicit none
52 
53  private
54 
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        ! nonlinear electro-optic susceptibility for semiconductors
61 
62 contains

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)
  nband_sum=Number of bands included in the sum. Must be <= mband
  pmat(mband,mband,nkpt,3,nsppol) = 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

SOURCE

1330 subroutine linelop(icomp, itemp, nband_sum, cryst, ks_ebands, &
1331                    pmat,v1,v2,v3,nmesh,de,sc,brod,tol,fnam,do_antiresonant,ncid,comm)
1332 
1333 !Arguments ------------------------------------
1334  integer, intent(in) :: icomp, itemp, nband_sum, ncid
1335  type(crystal_t),intent(in) :: cryst
1336  type(ebands_t),intent(in) :: ks_ebands
1337  complex(dpc), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
1338  integer, intent(in) :: v1, v2, v3
1339  integer, intent(in) :: nmesh
1340  integer, intent(in) :: comm
1341  real(dp), intent(in) :: de
1342  real(dp), intent(in) :: sc
1343  real(dp), intent(in) :: brod
1344  real(dp), intent(in) :: tol
1345  character(len=*), intent(in) :: fnam
1346  logical, intent(in) :: do_antiresonant
1347 
1348 !Local variables -------------------------
1349  integer,parameter :: master = 0
1350  integer :: iw
1351  integer :: i,j,k,lx,ly,lz
1352  integer :: isp,isym,ik
1353  integer :: ist1,istl,istn,istm, mband
1354  real(dp) :: ha2ev
1355  real(dp) :: ene,totre,totabs,totim
1356  real(dp) :: el,en,em
1357  real(dp) :: emin,emax,my_emin,my_emax
1358  real(dp) :: const_esu,const_au,au2esu
1359  real(dp) :: wmn,wnm,wln,wnl,wml,wlm
1360  complex(dpc) :: idel,w,zi
1361  character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5
1362 ! local allocatable arrays
1363  real(dp), allocatable :: s(:,:), sym(:,:,:)
1364  integer :: start4(4),count4(4)
1365  integer :: istp
1366  real(dp) :: ep, wmp, wpn
1367  real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included !
1368  real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, fmn
1369  complex(dpc), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a}
1370  complex(dpc), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a}
1371  complex(dpc), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k)
1372  complex(dpc), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c
1373  complex(dpc), allocatable :: chi(:) ! \chi_{II}^{abc}(-\omega,\omega,0)
1374  complex(dpc), allocatable :: eta(:) ! \eta_{II}^{abc}(-\omega,\omega,0)
1375  complex(dpc), allocatable :: sigma(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0)
1376  complex(dpc), allocatable :: chi2tot(:)
1377  complex(dpc) :: num1, num2, den1, den2, term1, term2
1378  complex(dpc) :: chi1, chi1_1, chi1_2, chi2_1b, chi2_2b
1379  complex(dpc), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega)
1380  complex(dpc) :: eta1, eta2, eta2_1, eta2_2
1381  complex(dpc) :: sigma1, sigma1_1, sigma1_2, sigma2
1382  !Parallelism
1383  integer :: my_rank, nproc, ierr, my_k1, my_k2
1384  integer :: fout1,fout2,fout3,fout4,fout5
1385  character(500) :: msg
1386 
1387 ! *********************************************************************
1388 
1389  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
1390 
1391 !calculate the constant
1392  zi=(0._dp,1._dp)
1393  idel=zi*brod
1394 ! Disable symmetries for now
1395  const_au=-2._dp/(cryst%ucvol*dble(cryst%nsym))
1396  au2esu=5.8300348177d-8
1397  const_esu=const_au*au2esu
1398  ha2ev=13.60569172*2._dp
1399 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1400 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
1401 !bohr: 5.2917ifc nlinopt.f907E-11
1402 !c: 2.99792458   velocity of sound
1403 !ry2ev: 13.60569172
1404 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2)
1405 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
1406 !mass comes from converting P_mn to r_mn
1407 !hbar^3 comes from converting all frequencies to energies in denominator
1408 !hbar^3 comes from operator for momentum (hbar/i nabla)
1409 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1410 !output file names
1411  fnam1=trim(fnam)//'-ChiEOTotIm.out'
1412  fnam2=trim(fnam)//'-ChiEOTotRe.out'
1413  fnam3=trim(fnam)//'-ChiEOIm.out'
1414  fnam4=trim(fnam)//'-ChiEORe.out'
1415  fnam5=trim(fnam)//'-ChiEOAbs.out'
1416 
1417  ! If there exists inversion symmetry exit with a mesg.
1418  if (cryst%idx_spatial_inversion() /= 0) then
1419    write(std_out,*) '-----------------------------------------'
1420    write(std_out,*) '    the crystal has inversion symmetry   '
1421    write(std_out,*) '    the LEO susceptibility is zero       '
1422    write(std_out,*) '-----------------------------------------'
1423    ABI_ERROR("Aborting now")
1424  end if
1425 
1426  ! check polarisation
1427  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
1428    write(std_out,*) '---------------------------------------------'
1429    write(std_out,*) '    Error in linelop:                        '
1430    write(std_out,*) '    the polarisation directions incorrect    '
1431    write(std_out,*) '    1=x,  2=y  and 3=z                       '
1432    write(std_out,*) '---------------------------------------------'
1433    ABI_ERROR("Aborting now")
1434  end if
1435 
1436  ! number of energy mesh points
1437  if (nmesh.le.0) then
1438    write(std_out,*) '---------------------------------------------'
1439    write(std_out,*) '    Error in linelop:                        '
1440    write(std_out,*) '    number of energy mesh points incorrect   '
1441    write(std_out,*) '    number has to be integer greater than 0  '
1442    write(std_out,*) '    nmesh*de = max energy for calculation    '
1443    write(std_out,*) '---------------------------------------------'
1444    ABI_ERROR("Aborting now")
1445  end if
1446 
1447  ! step in energy
1448  if (de.le.zero) then
1449    write(std_out,*) '---------------------------------------------'
1450    write(std_out,*) '    Error in linelop:                        '
1451    write(std_out,*) '    energy step is incorrect                 '
1452    write(std_out,*) '    number has to real greater than 0.0      '
1453    write(std_out,*) '    nmesh*de = max energy for calculation    '
1454    write(std_out,*) '---------------------------------------------'
1455    ABI_ERROR("Aborting now")
1456  end if
1457 
1458  ! broadening
1459  if (brod.gt.0.009) then
1460    write(std_out,*) '---------------------------------------------'
1461    write(std_out,*) '    ATTENTION: broadening is quite high      '
1462    write(std_out,*) '    ideally should be less than 0.005        '
1463    write(std_out,*) '---------------------------------------------'
1464  else if (brod.gt.0.015) then
1465    write(std_out,*) '----------------------------------------'
1466    write(std_out,*) '    ATTENTION: broadening is too high   '
1467    write(std_out,*) '    ideally should be less than 0.005   '
1468    write(std_out,*) '----------------------------------------'
1469  end if
1470 
1471  ! tolerance
1472  if (tol.gt.0.006) then
1473    write(std_out,*) '----------------------------------------'
1474    write(std_out,*) '    ATTENTION: tolerance is too high    '
1475    write(std_out,*) '    ideally should be less than 0.004   '
1476    write(std_out,*) '----------------------------------------'
1477  end if
1478 
1479  mband = ks_ebands%mband
1480  ABI_MALLOC(enk, (mband))
1481  ABI_MALLOC(delta, (mband, mband, 3))
1482  ABI_MALLOC(rmnbc,(mband,mband, 3, 3))
1483  ABI_MALLOC(roverw,(mband, mband, 3, 3))
1484  ABI_MALLOC(rmna, (mband, mband, 3))
1485  ABI_MALLOC(chi, (nmesh))
1486  ABI_MALLOC(eta, (nmesh))
1487  ABI_MALLOC(sigma, (nmesh))
1488  ABI_MALLOC(chi2, (nmesh))
1489  ABI_MALLOC(sym, (3, 3, 3))
1490  ABI_MALLOC(s, (3, 3))
1491 
1492  ABI_CHECK(nband_sum <= mband, "nband_sum <= mband")
1493 
1494  ! generate the symmetrizing tensor
1495  sym(:,:,:)=zero
1496  do isym=1,cryst%nsym
1497    s(:,:)=cryst%symrel_cart(:,:,isym)
1498    do i=1,3
1499      do j=1,3
1500        do k=1,3
1501          sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
1502        end do
1503      end do
1504    end do
1505  end do
1506 
1507  ! initialise
1508  delta(:,:,:)=zero
1509  rmnbc(:,:,:,:)=zero
1510  chi(:)=zero
1511  chi2(:) = zero
1512  eta(:)=zero
1513  sigma(:)=zero
1514  my_emin=HUGE(zero)
1515  my_emax=-HUGE(zero)
1516 
1517  ! Split work
1518  call xmpi_split_work(ks_ebands%nkpt,comm,my_k1,my_k2)
1519 
1520  ! loop over kpts
1521  do ik=my_k1,my_k2
1522    write(std_out,*) "P-",my_rank,": ",ik,'of',ks_ebands%nkpt
1523    do isp=1,ks_ebands%nsppol
1524      ! Calculate the scissor corrected energies and the energy window
1525      do ist1=1,nband_sum
1526        en = ks_ebands%eig(ist1,ik,isp)
1527        my_emin=min(my_emin,en)
1528        my_emax=max(my_emax,en)
1529        if(en > ks_ebands%fermie) then
1530          en = en + sc
1531        end if
1532        enk(ist1) = en
1533      end do
1534 
1535      ! calculate \Delta_nm and r_mn^a
1536      do istn=1,nband_sum
1537        en = enk(istn)
1538        do istm=1,nband_sum
1539          em = enk(istm)
1540          wmn = em - en
1541          delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp)
1542          if(abs(wmn) < tol) then
1543            rmna(istm,istn,1:3) = zero
1544          else
1545            rmna(istm,istn,1:3)=-zi*pmat(istm,istn,ik,1:3,isp)/wmn
1546          end if
1547        end do
1548      end do
1549 
1550      ! calculate \r^b_mn;c
1551      do istm=1,nband_sum
1552        em = enk(istm)
1553        do istn=1,nband_sum
1554          en = enk(istn)
1555          wmn = em - en
1556          if(abs(wmn) > tol) then
1557            do ly = 1,3
1558              do lz = 1,3
1559                num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly))
1560                den1 = wmn
1561                term1 = num1/den1
1562                term2 = zero
1563                do istp=1,nband_sum
1564                  ep = enk(istp)
1565                  wmp = em - ep
1566                  wpn = ep - en
1567                  num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly))
1568                  den2 = wmn
1569                  term2 = term2 + (num2/den2)
1570                end do
1571                rmnbc(istm,istn,ly,lz) = -term1-(zi*term2)
1572                roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz)
1573              end do
1574            end do
1575          end if
1576        end do
1577      end do
1578 
1579      ! initialise the factors
1580      ! start the calculation
1581      do istn=1,nband_sum
1582        en=enk(istn)
1583        if (do_antiresonant .and. en .ge. ks_ebands%fermie) then
1584          cycle
1585        end if
1586        fn=ks_ebands%occ(istn,ik,isp)
1587        do istm=1,nband_sum
1588          em=enk(istm)
1589          if (do_antiresonant .and. em .le. ks_ebands%fermie) then
1590            cycle
1591          end if
1592          wmn=em-en
1593          wnm=-wmn
1594          fm = ks_ebands%occ(istm,ik,isp)
1595          fnm = fn - fm
1596          fmn = fm - fn
1597          eta1 = zero
1598          eta2_1 = zero
1599          eta2_2 = zero
1600          sigma1_1 = zero
1601          sigma1_2 = zero
1602          sigma2 = zero
1603          if(abs(wmn) > tol) then
1604            do lx = 1,3
1605              do ly = 1,3
1606                do lz = 1,3
1607                  eta1 = eta1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*(roverw(istm,istn,lz,ly)))
1608                  eta2_1 = eta2_1 + sym(lx,ly,lz)*(fnm*(rmna(istn,istm,lx)*rmnbc(istm,istn,ly,lz)))
1609                  eta2_2 = eta2_2 + sym(lx,ly,lz)*(fnm*(rmnbc(istn,istm,lx,lz)*rmna(istm,istn,ly)))
1610                  sigma1_1 = sigma1_1 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,ly)*rmna(istm,istn,lz))/(wmn**2)
1611                  sigma1_2 = sigma1_2 + sym(lx,ly,lz)*(fnm*delta(istn,istm,lx)*rmna(istn,istm,lz)*rmna(istm,istn,ly))/(wmn**2)
1612                  sigma2 = sigma2 + sym(lx,ly,lz)*(fnm*rmnbc(istn,istm,lz,lx)*rmna(istm,istn,ly))/wmn
1613                end do
1614              end do
1615            end do
1616          end if
1617          chi1_1 = zero
1618          chi1_2 = zero
1619          chi2_1b = zero
1620          chi2_2b = zero
1621          chi2(:) = zero
1622          ! Three band terms
1623          do istl=1,nband_sum
1624            el=enk(istl)
1625            fl = ks_ebands%occ(istl,ik,isp)
1626            wlm = el-em
1627            wln = el-en
1628            wnl = en-el
1629            wml = em-el
1630            fnl = fn-fl
1631            fln = fl-fn
1632            fml = fm-fl
1633            do lx = 1,3
1634              do ly = 1,3
1635                do lz = 1,3
1636                  if(abs(wlm) > tol) then
1637                    chi1_1 = chi1_1 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,lz)*rmna(istl,istn,ly))/(wlm)
1638                    chi2_1b = chi2_1b + sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*rmna(istl,istm,lz)*rmna(istm,istn,ly))/(wlm)
1639                  end if
1640                  if(abs(wln) > tol) then
1641                    chi1_2 = chi1_2 + sym(lx,ly,lz)*(fnm*rmna(istn,istm,lx)*rmna(istm,istl,ly)*rmna(istl,istn,lz))/(wln)
1642                    chi2_2b = chi2_2b + sym(lx,ly,lz)*(fmn*rmna(istl,istm,lx)*rmna(istm,istn,ly)*rmna(istn,istl,lz))/(wnl)
1643                  end if
1644                end do
1645              end do
1646            end do
1647          end do
1648 
1649          sigma1 = 0.5_dp*(sigma1_1-sigma1_2)
1650          eta2 = 0.5_dp*(eta2_1-eta2_2)
1651          chi1 = chi1_1 + chi1_2
1652 
1653          !  calculate over the desired energy mesh and sum over k-points
1654          do iw=1,nmesh
1655            w=(iw-1)*de+idel
1656            ! Better way to compute it
1657            chi(iw) = chi(iw) + 0.5_dp*ks_ebands%wtk(ik)*((chi1/(wmn-w)) + ((chi2_1b+chi2_2b)/(wmn-w)))*const_esu
1658            eta(iw) = eta(iw) + 0.5_dp*zi*ks_ebands%wtk(ik)*((eta1/(wmn-w)) + (eta2/((wmn-w)**2)))*const_esu
1659            sigma(iw) = sigma(iw) + 0.5_dp*zi*ks_ebands%wtk(ik)*((sigma1/(wmn-w))- (sigma2/(wmn-w)))*const_esu
1660          end do
1661        end do ! istn and istm
1662      end do
1663    end do ! spins
1664  end do ! k-points
1665 
1666  call xmpi_sum(chi,comm,ierr)
1667  call xmpi_sum(eta,comm,ierr)
1668  call xmpi_sum(sigma,comm,ierr)
1669  call xmpi_min(my_emin,emin,comm,ierr)
1670  call xmpi_max(my_emax,emax,comm,ierr)
1671 
1672  if (my_rank == master) then
1673 
1674    if (ncid /= nctk_noid) then
1675      start4 = [1, 1, icomp, itemp]
1676      count4 = [2, nmesh, 1, 1]
1677      ABI_MALLOC(chi2tot, (nmesh))
1678      chi2tot = chi + eta + sigma
1679 #ifdef HAVE_NETCDF
1680      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi"), c2r(chi), start=start4, count=count4))
1681      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_eta"), c2r(eta), start=start4, count=count4))
1682      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_sigma"), c2r(sigma), start=start4, count=count4))
1683      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo_chi2tot"), c2r(chi2tot), start=start4, count=count4))
1684 #endif
1685      ABI_FREE(chi2tot)
1686    end if
1687 
1688   ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
1689    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
1690      ABI_ERROR(msg)
1691    end if
1692    if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
1693      ABI_ERROR(msg)
1694    end if
1695    if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
1696      ABI_ERROR(msg)
1697    end if
1698    if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
1699      ABI_ERROR(msg)
1700    end if
1701    if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
1702      ABI_ERROR(msg)
1703    end if
1704    ! write headers
1705    write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1706    write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
1707    write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
1708    write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
1709    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1710    write(fout1, '(a)' )' # Energy      Tot-Im Chi(-w,w,0)  Tot-Im Chi(-w,w,0)'
1711    write(fout1, '(a)' )' # eV          *10^-7 esu        *10^-12 m/V SI units '
1712    write(fout1, '(a)' )' # '
1713 
1714    write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1715    write(fout2, '(a,es16.6)') ' #tolerance:',tol
1716    write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1717    write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1718    write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1719    write(fout2, '(a)')' # Energy      Tot-Re Chi(-w,w,0)  Tot-Re Chi(-w,w,0)'
1720    write(fout2, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
1721    write(fout2, '(a)')' # '
1722 
1723    write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1724    write(fout3, '(a,es16.6)') ' #tolerance:',tol
1725    write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1726    write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1727    write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1728    write(fout3, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1729    write(fout3, '(a)')' # in esu'
1730    write(fout3, '(a)')' # '
1731 
1732    write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1733    write(fout4, '(a,es16.6)') ' #tolerance:',tol
1734    write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1735    write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1736    write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1737    write(fout4, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1738    write(fout4, '(a)')' # in esu'
1739    write(fout4, '(a)')' # '
1740 
1741    write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1742    write(fout5, '(a,es16.6)') ' #tolerance:',tol
1743    write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1744    write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1745    write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1746    write(fout5, '(a)')' # Energy(eV)  |TotChi(-w,w,0)|   |Tot Chi(-w,w,0)|'
1747    write(fout5, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
1748    write(fout5, '(a)')' # '
1749 
1750    totim=zero
1751    totre=zero
1752    totabs=zero
1753    do iw=2,nmesh
1754      ene=(iw-1)*de
1755      ene=ene*ha2ev
1756      totim=aimag(chi(iw)+eta(iw)+sigma(iw))/1.d-7
1757      write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1758      totim=zero
1759      totre=dble(chi(iw)+eta(iw)+sigma(iw))/1.d-7
1760      write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1761      totre=zero
1762      write(fout3,'(f15.6,3es15.6)') ene,aimag(chi(iw))/1.d-7,      &
1763      aimag(eta(iw))/1.d-7,aimag(sigma(iw))/1.d-7
1764      write(fout4,'(f15.6,3es15.6)') ene,dble(chi(iw))/1.d-7,       &
1765      dble(eta(iw))/1.d-7,dble(sigma(iw))/1.d-7
1766      totabs=abs(chi(iw)+eta(iw)+sigma(iw))/1.d-7
1767      write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1768      totabs=zero
1769    end do
1770 
1771    close(fout1)
1772    close(fout2)
1773    close(fout3)
1774    close(fout4)
1775    close(fout5)
1776    ! print information
1777    write(std_out,*) ' '
1778    write(std_out,*) 'information about calculation just performed:'
1779    write(std_out,*) ' '
1780    write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of LEO susceptibility'
1781    write(std_out,*) 'tolerance:',tol
1782    if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
1783    write(std_out,*) 'broadening:',brod,'Hartree'
1784    if (brod.gt.0.009) then
1785      write(std_out,*) ' '
1786      write(std_out,*) 'ATTENTION: broadening is quite high'
1787      write(std_out,*) ' '
1788    else if (brod.gt.0.015) then
1789      write(std_out,*) ' '
1790      write(std_out,*) 'ATTENTION: broadening is too high'
1791      write(std_out,*) ' '
1792    end if
1793    write(std_out,*) 'scissors shift:',sc,'Hartree'
1794    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree'
1795 
1796  end if
1797 
1798  ! deallocate local arrays
1799  ABI_FREE(enk)
1800  ABI_FREE(delta)
1801  ABI_FREE(rmnbc)
1802  ABI_FREE(roverw)
1803  ABI_FREE(rmna)
1804  ABI_FREE(chi)
1805  ABI_FREE(chi2)
1806  ABI_FREE(eta)
1807  ABI_FREE(sigma)
1808  ABI_FREE(s)
1809  ABI_FREE(sym)
1810 
1811 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)
  nband_sum=Number of bands included in the sum. Must be <= mband
  pmat(mband,mband,nkpt,3,nsppol)=momentum matrix elements in cartesian coordinates(complex)
  v1,v2=desired component of the dielectric function(integer) 1=x,2=y,3=z
  nmesh=desired number of energy mesh points(integer)
  de=desired step in energy(real); nmesh*de=maximum energy
  sc=scissors shift in Ha(real)
  brod=broadening in Ha(real)
  fnam=root for filename that will contain the output filename will be trim(fnam)//'-linopt.out'
  ncid=Netcdf id to save output data.
  prtlincompmatrixelements=if set to 1, the matrix elements are dumped in the _OPTIC.nc file for post processing.

SIDE EFFECTS

  Dielectric function for semiconductors, on a desired energy mesh and for a desired
  direction of polarisation is written to file.
  The output is in a file named trim(fnam)//'-linopt.out' and contains
  Im(\epsilon_{v1v2}(\omega), Re(\epsilon_{v1v2}(\omega) and abs(\epsilon_{v1v2}(\omega).

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

  Comment:
    Right now the routine sums over the kpoints. In future linear tetrahedron method should be useful.

SOURCE

230 subroutine linopt(icomp, itemp, nband_sum, cryst, ks_ebands, EPBSt, pmat, &
231   v1, v2, nmesh, de, sc, brod, fnam, ncid, prtlincompmatrixelements, comm)
232 
233 !Arguments ------------------------------------
234 integer, intent(in) :: icomp,itemp,nband_sum, ncid
235 type(crystal_t), intent(in) :: cryst
236 type(ebands_t),intent(in) :: ks_ebands,EPBSt
237 complex(dpc), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
238 integer, intent(in) :: v1, v2, nmesh
239 real(dp), intent(in) :: de, sc, brod
240 character(len=*), intent(in) :: fnam
241 integer, intent(in) :: comm
242 integer, intent(in) :: prtlincompmatrixelements
243 
244 !Local variables -------------------------
245 integer,parameter :: master=0
246 integer :: isp,i,j,isym,lx,ly,ik,ist1,ist2,iw,nkpt
247 integer :: my_rank, nproc, my_k1, my_k2, ierr, fout1, mband, nsppol
248 #ifdef HAVE_NETCDF
249 integer :: ncerr
250 #endif
251 logical :: do_linewidth
252 real(dp) :: deltav1v2, ha2ev, tmpabs, renorm_factor,emin,emax
253 real(dp) :: ene,abs_eps,re_eps
254 complex(dpc) :: e1,e2,e12, e1_ep,e2_ep,e12_ep, b11,b12, ieta, w
255 character(len=fnlen) :: fnam1
256 character(len=500) :: msg
257 ! allocatable arrays
258 real(dp) :: s(3,3),sym(3,3)
259 real(dp), allocatable :: im_refract(:),re_refract(:)
260 complex(dpc), allocatable :: chi(:,:), matrix_elements(:,:,:,:), renorm_eigs(:,:,:), eps(:)
261 
262 ! *********************************************************************
263 
264  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
265  nkpt = ks_ebands%nkpt
266  nsppol = ks_ebands%nsppol
267  mband = ks_ebands%mband
268  ABI_CHECK(nband_sum <= mband, "nband_sum <= mband")
269 
270  if (my_rank == master) then
271    ! check polarisation
272    if (v1.le.0.or.v2.le.0.or.v1.gt.3.or.v2.gt.3) then
273      write(std_out,*) '---------------------------------------------'
274      write(std_out,*) '    Error in linopt:                         '
275      write(std_out,*) '    the polarisation directions incorrect    '
276      write(std_out,*) '    1=x and 2=y and 3=z                      '
277      write(std_out,*) '---------------------------------------------'
278      ABI_ERROR("Aborting now")
279    end if
280    ! number of energy mesh points
281    if (nmesh.le.0) then
282      write(std_out,*) '---------------------------------------------'
283      write(std_out,*) '    Error in linopt:                         '
284      write(std_out,*) '    number of energy mesh points incorrect   '
285      write(std_out,*) '    number has to integer greater than 0     '
286      write(std_out,*) '    nmesh*de = max energy for calculation    '
287      write(std_out,*) '---------------------------------------------'
288      ABI_ERROR("Aborting now")
289    end if
290    ! step in energy
291    if (de.le.zero) then
292      write(std_out,*) '---------------------------------------------'
293      write(std_out,*) '    Error in linopt:                         '
294      write(std_out,*) '    energy step is incorrect                 '
295      write(std_out,*) '    number has to real greater than 0.0      '
296      write(std_out,*) '    nmesh*de = max energy for calculation    '
297      write(std_out,*) '---------------------------------------------'
298      ABI_ERROR("Aborting now")
299    end if
300    ! broadening
301    if (brod.gt.0.009) then
302      write(std_out,*) '---------------------------------------------'
303      write(std_out,*) '    ATTENTION: broadening is quite high      '
304      write(std_out,*) '    ideally should be less than 0.005        '
305      write(std_out,*) '---------------------------------------------'
306    else if (brod.gt.0.015) then
307      write(std_out,*) '----------------------------------------'
308      write(std_out,*) '    ATTENTION: broadening is too high   '
309      write(std_out,*) '    ideally should be less than 0.005   '
310      write(std_out,*) '----------------------------------------'
311    end if
312    ! fermi energy
313    if(ks_ebands%fermie<-1.0d4) then
314      write(std_out,*) '---------------------------------------------'
315      write(std_out,*) '    ATTENTION: Fermi energy seems extremely  '
316      write(std_out,*) '    low                                      '
317      write(std_out,*) '---------------------------------------------'
318    end if
319    ! scissors operator
320    if (sc.lt.zero) then
321      write(std_out,*) '---------------------------------------------'
322      write(std_out,*) '    Error in linopt:                         '
323      write(std_out,*) '    scissors shift is incorrect              '
324      write(std_out,*) '    number has to be greater than 0.0      '
325      write(std_out,*) '---------------------------------------------'
326      ABI_ERROR("Aborting now")
327    end if
328  end if
329 
330  do_linewidth = allocated(EPBSt%linewidth)
331 ! TODO: activate this, and remove do_linewidth - always add it in even if 0.
332 ! if (.not. allocated(EPBSt%linewidth)) then
333 !   ABI_MALLOC(EPBSt%linewidth, (1, mband, my_k2-my_k1+1, nsppol))
334 !   EPBSt%linewidth = zero
335 ! end if
336 
337  ABI_MALLOC(chi, (nmesh, nsppol))
338  ABI_MALLOC(eps, (nmesh))
339  ABI_MALLOC(im_refract, (nmesh))
340  ABI_MALLOC(re_refract, (nmesh))
341  ieta=(zero, 1._dp)*brod
342  renorm_factor=1._dp/(cryst%ucvol*dble(cryst%nsym))
343  ha2ev=13.60569172*2._dp
344 
345  ! output file names
346  fnam1=trim(fnam)//'-linopt.out'
347 
348  ! construct symmetrisation tensor
349  sym = zero
350  do isym=1,cryst%nsym
351    s(:,:)=cryst%symrel_cart(:,:,isym)
352    do i=1,3
353      do j=1,3
354        sym(i,j)=sym(i,j)+s(i,v1)*s(j,v2)
355      end do
356    end do
357  end do
358 
359  ! calculate the energy window
360  emin=zero
361  emax=zero
362  do ik=1,nkpt
363    do isp=1,nsppol
364      do ist1=1,nband_sum
365        emin=min(emin,EPBSt%eig(ist1,ik,isp))
366        emax=max(emax,EPBSt%eig(ist1,ik,isp))
367      end do
368    end do
369  end do
370 
371  ! Split work
372  call xmpi_split_work(nkpt,comm,my_k1,my_k2)
373  ! if we print matrix elements, allocate full arrays for each process
374  ! this is not optimized memory-wise since we could just allocate what is needed
375  ! however we would need to write all data using mpi-io.
376  if (prtlincompmatrixelements == 1) then
377    ABI_CALLOC(matrix_elements, (mband, mband, nkpt, nsppol))
378    ABI_CALLOC(renorm_eigs, (mband, nkpt, nsppol))
379  endif
380 
381  ! start calculating linear optical response
382  chi(:,:)=zero
383  do isp=1,nsppol
384    do ik=my_k1,my_k2
385      write(std_out,*) "P-",my_rank,": ",ik,'of',nkpt
386      do ist1=1,nband_sum
387        e1=ks_ebands%eig(ist1,ik,isp)
388        e1_ep=EPBSt%eig(ist1,ik,isp)
389        ! TODO: unless memory is a real issue, should set lifetimes to 0 and do this sum systematically
390        ! instead of putting an if statement in a loop! See above
391        if(do_linewidth) then
392          e1_ep = e1_ep + EPBSt%linewidth(1,ist1,ik,isp)*(0.0_dp,1.0_dp)
393        end if
394        do ist2=1,nband_sum
395          e2=ks_ebands%eig(ist2,ik,isp)
396          e2_ep=EPBSt%eig(ist2,ik,isp)
397          if(do_linewidth) then
398            e2_ep = e2_ep - EPBSt%linewidth(1,ist2,ik,isp)*(0.0_dp,1.0_dp)
399          end if
400          if (ist1.ne.ist2) then
401            ! scissors correction of momentum matrix
402            if(REAL(e1) > REAL(e2)) then
403              e12 = e1-e2+sc
404            else
405              e12 = e1-e2-sc
406            end if
407            if(REAL(e1_ep) > REAL(e2_ep)) then
408              e12_ep = e1_ep-e2_ep+sc
409            else
410              e12_ep = e1_ep-e2_ep-sc
411            end if
412            ! e12=e1-e2-sc
413            b11=zero
414            ! symmetrization of momentum matrix
415            do lx=1,3
416              do ly=1,3
417                b11=b11+(sym(lx,ly)*pmat(ist1,ist2,ik,lx,isp)* &
418                conjg(pmat(ist1,ist2,ik,ly,isp)))
419              end do
420            end do
421            b12=b11*renorm_factor*(1._dp/(e12**2))
422            ! store data for printing if necessary
423            if (prtlincompmatrixelements == 1) then
424              matrix_elements(ist1,ist2,ik,isp) = b12
425              renorm_eigs(ist1,ik,isp) = e1_ep
426              renorm_eigs(ist2,ik,isp) = e2_ep
427            endif
428            ! calculate on the desired energy grid
429            do iw=2,nmesh
430              w=(iw-1)*de+ieta
431              chi(iw,isp)=chi(iw,isp)+(ks_ebands%wtk(ik)*(ks_ebands%occ(ist1,ik,isp)-ks_ebands%occ(ist2,ik,isp))* &
432              (b12/(-e12_ep-w)))
433            end do ! frequencies
434          end if
435        end do  ! states 2
436      end do  ! states 1
437    end do ! k points
438  end do ! spin
439 
440  call xmpi_sum(chi,comm,ierr)
441  if (prtlincompmatrixelements == 1) then
442    ! gather all data to main process in order to write them using a single process
443    ! in the netcdf file. This could be avoided by doing mpiio.
444    call xmpi_sum(matrix_elements,comm,ierr)
445    call xmpi_sum(renorm_eigs,comm,ierr)
446  endif
447 
448  ! calculate epsilon
449  eps(1) = zero
450  deltav1v2=zero; if (v1 == v2) deltav1v2=one
451  do iw=2,nmesh
452    eps(iw)=deltav1v2+4._dp*pi*sum(chi(iw,:))
453  end do
454 
455  if (my_rank == master) then
456    !  open the output files
457    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
458      ABI_ERROR(msg)
459    end if
460    ! write output
461    write(fout1, '(a,2i3,a)' )' #calculated the component:',v1,v2,'  of dielectric function'
462    write(std_out,*) 'calculated the component:',v1,v2,'  of dielectric function'
463    write(fout1, '(a,2es16.6)' ) ' #broadening:', real(ieta),aimag(ieta)
464    write(std_out,*) ' with broadening:',ieta
465    write(fout1, '(a,es16.6)' ) ' #scissors shift:',sc
466    write(std_out,*) 'and scissors shift:',sc
467    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
468    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
469    write(fout1,*)
470    if(nsppol==1)write(fout1, '(a)' ) ' # Energy(eV)         Im(eps(w))'
471    if(nsppol==2)write(fout1, '(a)' ) ' # Energy(eV)         Im(eps(w))         Spin up       Spin down '
472    do iw=2,nmesh
473      ene=(iw-1)*de*ha2ev
474      if(nsppol==1)write(fout1, '(2es16.6)' ) ene,aimag(eps(iw))
475      if(nsppol==2)write(fout1, '(4es16.6)' ) ene,aimag(eps(iw)),4._dp*pi*aimag(chi(iw,1)),4._dp*pi*aimag(chi(iw,2))
476    end do
477    write(fout1,*)
478    write(fout1,*)
479    if(nsppol==1)write(fout1, '(a)' ) ' # Energy(eV)         Re(eps(w))'
480    if(nsppol==2)write(fout1, '(a)' ) ' # Energy(eV)         Re(eps(w))         Spin up       Spin down    +delta(diag) '
481    do iw=2,nmesh
482      ene=(iw-1)*de*ha2ev
483      if(nsppol==1)write(fout1, '(2es16.6)' ) ene,dble(eps(iw))
484      if(nsppol==2)write(fout1, '(5es16.6)' ) ene,dble(eps(iw)),4._dp*pi*dble(chi(iw,1)),4._dp*pi*dble(chi(iw,2)),deltav1v2
485    end do
486    write(fout1,*)
487    write(fout1,*)
488    write(fout1, '(a)' )' # Energy(eV)         abs(eps(w))'
489    do iw=2,nmesh
490      ene=(iw-1)*de*ha2ev
491      abs_eps=abs(eps(iw))
492      re_eps=dble(eps(iw))
493      write(fout1, '(2es16.6)' ) ene,abs_eps
494      re_refract(iw)=sqrt(half*(abs_eps+re_eps))
495      im_refract(iw)=sqrt(half*(abs_eps-re_eps))
496    end do
497    write(fout1,*)
498    write(fout1,*)
499    write(fout1, '(a)' )' # Energy(eV)         Im(refractive index(w)) aka kappa'
500    do iw=2,nmesh
501      ene=(iw-1)*de*ha2ev
502      write(fout1, '(2es16.6)' ) ene,im_refract(iw)
503    end do
504    write(fout1,*)
505    write(fout1,*)
506    write(fout1, '(a)' )' # Energy(eV)         Re(refractive index(w)) aka n'
507    do iw=2,nmesh
508      ene=(iw-1)*de*ha2ev
509      write(fout1, '(2es16.6)' ) ene,re_refract(iw)
510    end do
511    write(fout1,*)
512    write(fout1,*)
513    write(fout1, '(a)' )' # Energy(eV)         Reflectivity(w) from vacuum, at normal incidence'
514    do iw=2,nmesh
515      ene=(iw-1)*de*ha2ev
516      write(fout1, '(2es16.6)' ) ene, ((re_refract(iw)-one)**2+im_refract(iw)**2)/((re_refract(iw)+one)**2+im_refract(iw)**2)
517    end do
518    write(fout1,*)
519    write(fout1,*)
520    write(fout1, '(a)' )' # Energy(eV)         absorption coeff (in 10^6 m-1) = omega Im(eps) / c n(eps)'
521    do iw=2,nmesh
522      ene=(iw-1)*de
523      tmpabs=zero
524      if ( re_refract(iw) > tol10 ) then
525        tmpabs = aimag(eps(iw))*ene / re_refract(iw) / Sp_Lt / Bohr_meter * 1.0d-6
526      end if
527      write(fout1, '(2es16.6)' ) ha2ev*ene, tmpabs
528    end do
529 
530    ! close output file
531    close(fout1)
532 
533 #ifdef HAVE_NETCDF
534    if (ncid /= nctk_noid) then
535      ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_epsilon"), c2r(eps), start=[1, 1, icomp, itemp])
536      NCF_CHECK(ncerr)
537    end if
538    if (prtlincompmatrixelements == 1) then
539      ! write matrix elements and other quantities used to build the chi tensor.
540      write(std_out, '(a)') 'Writing linopt matrix elements in _OPTIC.nc file.'
541      ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_matrix_elements"), c2r(matrix_elements),&
542                           start=[1, 1, 1, 1, 1, icomp, itemp])
543      NCF_CHECK(ncerr)
544      ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_renorm_eigs"), c2r(renorm_eigs), start=[1, 1, 1, 1])
545      NCF_CHECK(ncerr)
546 
547      ! write occupations and kpt weights
548      ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_occupations"), ks_ebands%occ, start=[1, 1, 1])
549      NCF_CHECK(ncerr)
550      ncerr = nf90_put_var(ncid, nctk_idname(ncid, "linopt_wkpts"), ks_ebands%wtk, start=[1])
551      NCF_CHECK(ncerr)
552      write(std_out, '(a)') 'Writing linopt matrix elements done.'
553    endif
554 
555 #endif
556  end if ! rank == master
557 
558  ABI_FREE(chi)
559  ABI_FREE(eps)
560  ABI_FREE(im_refract)
561  ABI_FREE(re_refract)
562 
563  ABI_SFREE(matrix_elements)
564  ABI_SFREE(renorm_eigs)
565 
566 end subroutine linopt

m_optic_tools/nlinopt [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 nlinopt

FUNCTION

 Compute 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)
  nband_sum=Number of bands included in the sum. Must be <= mband
  fermie = Fermi energy in Ha(real)
  pmat(mband,mband,nkpt,3,nsppol) = 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.
  See eqs. (A4)-A(11) of S. Sharma et al Phys. Rev. B 67, 165332 (2003)

SOURCE

 610 subroutine nlinopt(icomp, itemp, nband_sum, cryst, ks_ebands, pmat, &
 611                    v1, v2, v3, nmesh, de, sc, brod, tol, fnam, ncid, comm)
 612 
 613 !Arguments ------------------------------------
 614 integer, intent(in) :: icomp, itemp, nband_sum, ncid
 615 type(crystal_t),intent(in) :: cryst
 616 type(ebands_t),intent(in) :: ks_ebands
 617 complex(dpc), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
 618 integer, intent(in) :: v1, v2, v3, nmesh, comm
 619 real(dp), intent(in) :: de, sc, brod, tol
 620 character(len=*), intent(in) :: fnam
 621 
 622 !Local variables -------------------------
 623 integer,parameter :: master=0
 624 integer :: iw, mband,i,j,k,lx,ly,lz
 625 integer :: isp,isym,ik,ist1,ist2,istl,istn,istm
 626 integer :: my_rank, nproc, my_k1, my_k2, ierr
 627 integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7
 628 real(dp) :: f1,f2,f3, ha2ev
 629 real(dp) :: ene,totre,totabs,totim
 630 real(dp) :: e1,e2,el,en,em,emin,emax,my_emin,my_emax
 631 real(dp) :: const_esu,const_au,au2esu,wmn,wnm,wln,wnl,wml,wlm, t1
 632 complex(dpc) :: idel,w,zi
 633 complex(dpc) :: mat2w,mat1w1,mat1w2,mat2w_tra,mat1w3_tra
 634 complex(dpc) :: b111,b121,b131,b112,b122,b132,b113,b123,b133
 635 complex(dpc) :: b241,b242,b243,b221,b222,b223,b211,b212,b213,b231
 636 complex(dpc) :: b311,b312,b313,b331
 637 complex(dpc) :: b24,b21_22,b11,b12_13,b31_32
 638 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7
 639 character(500) :: msg
 640 ! local allocatable arrays
 641 integer :: start4(4),count4(4)
 642 real(dp) :: s(3,3),sym(3,3,3)
 643 complex(dpc), allocatable :: px(:,:,:,:,:), py(:,:,:,:,:), pz(:,:,:,:,:)
 644 complex(dpc), allocatable :: delta(:,:,:), inter2w(:), inter1w(:)
 645 complex(dpc), allocatable :: intra2w(:), intra1w(:), intra1wS(:),chi2tot(:)
 646 
 647 ! *********************************************************************
 648 
 649  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
 650  mband = ks_ebands%mband
 651 
 652 !calculate the constant
 653  zi=(0._dp,1._dp)
 654  idel=zi*brod
 655  const_au=-2._dp/(cryst%ucvol*dble(cryst%nsym))
 656  au2esu=5.8300348177d-8
 657  const_esu=const_au*au2esu
 658  ha2ev=13.60569172*2._dp
 659 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 660 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
 661 !bohr: 5.2917ifc nlinopt.f907E-11
 662 !c: 2.99792458   velocity of sound
 663 !ry2ev: 13.60569172
 664 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2)
 665 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
 666 !mass comes from converting P_mn to r_mn
 667 !hbar^3 comes from converting all frequencies to energies in denominator
 668 !hbar^3 comes from operator for momentum (hbar/i nabla)
 669 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 670 !output file names
 671  fnam1=trim(fnam)//'-ChiTotIm.out'
 672  fnam2=trim(fnam)//'-ChiTotRe.out'
 673  fnam3=trim(fnam)//'-ChiIm.out'
 674  fnam4=trim(fnam)//'-ChiRe.out'
 675  fnam5=trim(fnam)//'-ChiAbs.out'
 676  fnam6=trim(fnam)//'-ChiImDec.out'
 677  fnam7=trim(fnam)//'-ChiReDec.out'
 678 
 679  if(my_rank == master) then
 680    ! If there exists inversion symmetry exit with a message.
 681    if (cryst%idx_spatial_inversion() /= 0) then
 682      write(std_out,*) '-------------------------------------------'
 683      write(std_out,*) '    The crystal has inversion symmetry     '
 684      write(std_out,*) '    The SHG susceptibility is zero         '
 685      write(std_out,*) '    Action : set num_nonlin_comp to zero   '
 686      write(std_out,*) '-------------------------------------------'
 687      ABI_ERROR("Aborting now")
 688    end if
 689    ! check polarisation
 690    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
 691      write(std_out,*) '---------------------------------------------'
 692      write(std_out,*) '    Error in nlinopt:                        '
 693      write(std_out,*) '    Incorrect polarisation directions        '
 694      write(std_out,*) '    1=x,  2=y  and 3=z                       '
 695      write(std_out,*) '    Action : check your input file,          '
 696      write(std_out,*) '    use only 1, 2 or 3 to define directions  '
 697      write(std_out,*) '---------------------------------------------'
 698      ABI_ERROR("Aborting now")
 699    end if
 700    !number of energy mesh points
 701    if (nmesh.le.0) then
 702      write(std_out,*) '---------------------------------------------'
 703      write(std_out,*) '    Error in nlinopt:                        '
 704      write(std_out,*) '    number of energy mesh points incorrect   '
 705      write(std_out,*) '    number has to be integer greater than 0  '
 706      write(std_out,*) '    nmesh*de = max energy for calculation    '
 707      write(std_out,*) '---------------------------------------------'
 708      ABI_ERROR("Aborting now")
 709    end if
 710    !step in energy
 711    if (de.le.zero) then
 712      write(std_out,*) '---------------------------------------------'
 713      write(std_out,*) '    Error in nlinopt:                        '
 714      write(std_out,*) '    energy step is incorrect                 '
 715      write(std_out,*) '    number has to real greater than 0.0      '
 716      write(std_out,*) '    nmesh*de = max energy for calculation    '
 717      write(std_out,*) '---------------------------------------------'
 718      ABI_ERROR("Aborting now")
 719    end if
 720    !broadening
 721    if (brod.gt.0.009) then
 722      write(std_out,*) '---------------------------------------------'
 723      write(std_out,*) '    WARNING : broadening is quite high       '
 724      write(std_out,*) '    ideally should be less than 0.005        '
 725      write(std_out,*) '---------------------------------------------'
 726    else if (brod.gt.0.015) then
 727      write(std_out,*) '----------------------------------------'
 728      write(std_out,*) '    WARNING : broadening is too high    '
 729      write(std_out,*) '    ideally should be less than 0.005   '
 730      write(std_out,*) '----------------------------------------'
 731    end if
 732    !tolerance
 733    if (tol.gt.0.006) then
 734      write(std_out,*) '----------------------------------------'
 735      write(std_out,*) '    WARNING : tolerance is too high     '
 736      write(std_out,*) '    ideally should be less than 0.004   '
 737      write(std_out,*) '----------------------------------------'
 738    end if
 739  end if
 740 
 741  !allocate local arrays
 742  ABI_MALLOC(px, (mband, mband, 3, 3, 3))
 743  ABI_MALLOC(py, (mband, mband, 3, 3, 3))
 744  ABI_MALLOC(pz,(mband,mband, 3, 3, 3))
 745  ABI_MALLOC(inter2w, (nmesh))
 746  ABI_MALLOC(inter1w, (nmesh))
 747  ABI_MALLOC(intra2w, (nmesh))
 748  ABI_MALLOC(intra1w, (nmesh))
 749  ABI_MALLOC(intra1wS, (nmesh))
 750  ABI_MALLOC(delta, (mband, mband, 3))
 751 
 752  !generate the symmetrizing tensor
 753  sym = zero
 754  do isym=1,cryst%nsym
 755    s(:,:)=cryst%symrel_cart(:,:,isym)
 756    do i=1,3
 757      do j=1,3
 758        do k=1,3
 759          sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
 760        end do
 761      end do
 762    end do
 763  end do
 764  ! Disable symmetries for now
 765  !sym(:,:,:) = zero
 766  !sym(v1,v2,v3) = nsym
 767 
 768  ! Split work
 769  call xmpi_split_work(ks_ebands%nkpt, comm, my_k1, my_k2)
 770 
 771  ! initialise
 772  inter2w(:)=zero
 773  inter1w(:)=zero
 774  intra2w(:)=zero
 775  intra1w(:)=zero
 776  intra1wS(:)=zero
 777  delta(:,:,:)=zero
 778 
 779  my_emin=HUGE(zero)
 780  my_emax=-HUGE(zero)
 781 
 782  ! loop over kpts
 783  do ik=my_k1,my_k2
 784    write(std_out,*) "P-",my_rank,": ",ik,'of',ks_ebands%nkpt
 785    ! loop over spins
 786    do isp=1,ks_ebands%nsppol
 787      !  loop over states
 788      do ist1=1,nband_sum
 789        e1 = ks_ebands%eig(ist1,ik,isp)
 790        if (e1.lt.ks_ebands%fermie) then   ! ist1 is a valence state
 791          do ist2=1,nband_sum
 792            e2 = ks_ebands%eig(ist2,ik,isp)
 793            if (e2.gt.ks_ebands%fermie) then ! ist2 is a conduction state
 794              ! symmetrize the momentum matrix elements
 795              do lx=1,3
 796                do ly=1,3
 797                  do lz=1,3
 798                    f1=sym(lx,ly,lz)+sym(lx,lz,ly)
 799                    f2=sym(ly,lx,lz)+sym(ly,lz,lx)
 800                    f3=sym(lz,lx,ly)+sym(lz,ly,lx)
 801                    px(ist1,ist2,lx,ly,lz)=f1*pmat(ist1,ist2,ik,lx,isp)
 802                    py(ist2,ist1,lx,ly,lz)=f2*pmat(ist2,ist1,ik,lx,isp)
 803                    pz(ist2,ist1,lx,ly,lz)=f3*pmat(ist2,ist1,ik,lx,isp)
 804                  end do
 805                end do
 806              end do ! end loop over states
 807            end if
 808          end do
 809        end if
 810      end do
 811 
 812      ! calculate the energy window and \Delta_nm
 813      do ist1=1,nband_sum
 814        my_emin=min(my_emin, ks_ebands%eig(ist1,ik,isp))
 815        my_emax=max(my_emax, ks_ebands%eig(ist1,ik,isp))
 816        do ist2=1,nband_sum
 817          delta(ist1,ist2,1:3)=pmat(ist1,ist1,ik,1:3,isp)-pmat(ist2,ist2,ik,1:3,isp)
 818        end do
 819      end do
 820      ! initialise the factors
 821      ! factors are named according to the Ref. article 2.
 822      b111=zero
 823      b121=zero
 824      b131=zero
 825      b112=zero
 826      b122=zero
 827      b132=zero
 828      b113=zero
 829      b123=zero
 830      b133=zero
 831      b211=zero
 832      b221=zero
 833      b212=zero
 834      b222=zero
 835      b213=zero
 836      b223=zero
 837      b231=zero
 838      b241=zero
 839      b242=zero
 840      b243=zero
 841      b311=zero
 842      b312=zero
 843      b313=zero
 844      b331=zero
 845      ! start the calculation
 846      do istn=1,nband_sum
 847        en=ks_ebands%eig(istn,ik,isp)
 848        if (en.lt.ks_ebands%fermie) then  ! istn is a valence state
 849          do istm=1,nband_sum
 850            em=ks_ebands%eig(istm,ik,isp)
 851            if (em.gt.ks_ebands%fermie) then  ! istm is a conduction state
 852              em = em + sc ! Should add the scissor to conduction energies
 853              wmn=em-en
 854              wnm=-wmn
 855              ! calculate the matrix elements for two band intraband term
 856              mat2w_tra=zero
 857              mat1w3_tra=zero
 858              do lx=1,3
 859                do ly=1,3
 860                  do lz=1,3
 861 !                  See Sharma03, A10, last term of first line, corrected ! Should be Delta^b_mn r c_mn, see A2.
 862                    mat2w_tra=mat2w_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,lz,isp)    &
 863                    *delta(istm,istn,ly)
 864 !                  See Sharma03, A11, last term.
 865                    mat1w3_tra=mat1w3_tra+px(istn,istm,lx,ly,lz)*pmat(istm,istn,ik,ly,isp)  &
 866                    *delta(istm,istn,lz)
 867                    ! NOTE:: lx to ly m to n in pmat matrices respectively
 868                    ! Changes are made so that this (b3) term is according to paper
 869                    ! [[cite:Sipe1993]] (Ref. 4) rather than [[cite:Hughes1996]] (Ref 2) in which this term is incorrect
 870                  end do
 871                end do
 872              end do
 873              b331=mat1w3_tra/wnm
 874              b231=8._dp*mat2w_tra/wmn
 875 
 876              b11=zero
 877              b12_13=zero
 878              b24=zero
 879              b31_32=zero
 880              b21_22=zero
 881 
 882              ! istl < istn
 883              do istl=1,istn-1  ! istl is a valence state below istn
 884                el=ks_ebands%eig(istl,ik,isp)
 885                wln=el-en  ! do not add sc to the valence bands!
 886                wml=em-el
 887                wnl=-wln
 888                wlm=-wml
 889                ! calculate the matrix elements for three band terms
 890                mat2w=zero
 891                mat1w1=zero
 892                mat1w2=zero
 893                do lx=1,3
 894                  do ly=1,3
 895                    do lz=1,3
 896 
 897                      mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp)   &
 898                      *pmat(istl,istn,ik,lz,isp))
 899 
 900                      mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
 901                      *pmat(istn,istl,ik,ly,isp))
 902 
 903                      mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
 904                      *pmat(istn,istl,ik,ly,isp))
 905                    end do
 906                  end do
 907                end do
 908                b111=mat2w*(1._dp/(wln+wlm))*(1._dp/wlm)
 909                b121=mat1w1*(1._dp/(wnm+wlm))*(1._dp/wlm)
 910                b131=mat1w2*(1._dp/wlm)
 911                b221=zero
 912                b211=mat1w1/wml
 913                b241=-mat2w/wml
 914                b311=mat1w2/wlm
 915 
 916                if (abs(wln).gt.tol) then
 917                  b111=b111/wln
 918                  b121=b121/wln
 919                  b131=b131/wln
 920                  b221=mat1w2/wln
 921                  b241=b241+(mat2w/wln)
 922                  b311=b311+(mat1w1/wln)
 923                else
 924                  b111=zero
 925                  b121=zero
 926                  b131=zero
 927                  b221=zero
 928                end if
 929                t1=wln-wnm
 930                if (abs(t1).gt.tol) then
 931                  b131=b131/t1
 932                else
 933                  b131=zero
 934                end if
 935                b11=b11-2._dp*b111
 936                b12_13=b12_13+b121+b131
 937                b21_22=b21_22-b211+b221
 938                b24=b24+2._dp*b241
 939                b31_32=b31_32+b311
 940              end do ! istl
 941 
 942              ! istn < istl < istm
 943              do istl=istn+1,istm-1
 944                el=ks_ebands%eig(istl,ik,isp)
 945                ! calculate the matrix elements for three band terms
 946                mat2w=zero
 947                mat1w1=zero
 948                mat1w2=zero
 949                do lx=1,3
 950                  do ly=1,3
 951                    do lz=1,3
 952 
 953                      mat2w=mat2w+(px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp)   &
 954                      *pmat(istl,istn,ik,lz,isp))
 955 
 956                      mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
 957                      *pmat(istn,istl,ik,ly,isp))
 958 
 959                      mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
 960                      *pmat(istn,istl,ik,ly,isp))
 961                    end do
 962                  end do
 963                end do
 964                if (el.lt.ks_ebands%fermie) then
 965                  wln=el-en
 966                  wnl=-wln
 967                  wml=em-el
 968                  wlm=-wml
 969                else
 970                  el=el+sc
 971                  wln=el-en
 972                  wnl=-wln
 973                  wml=em-el
 974                  wlm=-wml
 975                end if
 976 !
 977                b112=zero
 978                b122=mat1w1*(1._dp/(wnm+wlm))
 979                b132=mat1w2*(1._dp/(wnm+wnl))
 980                b242=zero
 981                b222=zero
 982                b212=zero
 983                b312=zero
 984                if (abs(wnl).gt.tol) then
 985                  b112=mat2w/wln
 986                  b122=b122/wnl
 987                  b132=b132/wnl
 988                  b242=mat2w/wln
 989                  b222=mat1w2/wln
 990                  b312=mat1w1/wln
 991                else
 992                  b122=zero
 993                  b132=zero
 994                end if
 995                if (abs(wlm).gt.tol) then
 996                  b112=b112/wml
 997                  b122=b122/wlm
 998                  b132=b132/wlm
 999                  b242=b242-(mat2w/wml)
1000                  b212=mat1w1/wml
1001                  b312=b312+(mat1w2/wlm)
1002                else
1003                  b112=zero
1004                  b122=zero
1005                  b132=zero
1006                  b212=zero
1007                end if
1008                t1=wlm-wnl
1009                if (abs(t1).gt.tol) then
1010                  b112=b112/t1
1011                else
1012                  b112=zero
1013                end if
1014                b11=b11+2._dp*b112
1015                b12_13=b12_13-b122+b132
1016                b24=b24+2._dp*b242
1017                b21_22=b21_22-b212+b222
1018                b31_32=b31_32+b312
1019              end do ! istl
1020 
1021              ! istl > istm    !
1022              do istl=istm+1,nband_sum
1023                el=ks_ebands%eig(istl,ik,isp)+sc
1024                wln=el-en
1025                wnl=-wln
1026                wml=em-el
1027                wlm=-wml
1028                ! calculate the matrix elements for three band terms
1029                mat2w=zero
1030                mat1w1=zero
1031                mat1w2=zero
1032                do lx=1,3
1033                  do ly=1,3
1034                    do lz=1,3
1035                      mat2w=mat2w+px(istn,istm,lx,ly,lz)*pmat(istm,istl,ik,ly,isp) &
1036                      *pmat(istl,istn,ik,lz,isp)
1037 
1038                      mat1w1=mat1w1+(py(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1039                      *pmat(istn,istl,ik,ly,isp))
1040 
1041                      mat1w2=mat1w2+(pz(istm,istn,lx,ly,lz)*pmat(istl,istm,ik,lz,isp) &
1042                      *pmat(istn,istl,ik,ly,isp))
1043                    end do
1044                  end do
1045                end do
1046 
1047                b113=mat2w*(1._dp/(wnl+wml))*(1._dp/wnl)
1048                b123=mat1w1*(1._dp/wnl)
1049                b133=mat1w2*(1._dp/wnl)*(1._dp/(wnl+wnm))
1050                b243=mat2w/wln
1051                b223=mat1w2/wln
1052                b213=zero
1053                b313=-1._dp*mat1w1/wnl
1054                if (abs(wml).gt.tol) then
1055                  b113=b113/wml
1056                  b123=b123/wml
1057                  b133=b133/wml
1058                  b243=b243-(mat2w/wml)
1059                  b213=mat1w1/wml
1060                  b313=b313+(mat1w2/wlm)
1061                else
1062                  b113=zero
1063                  b123=zero
1064                  b133=zero
1065                end if
1066                t1=wnm-wml
1067                if (abs(t1).gt.tol) then
1068                  b123=b123/t1
1069                else
1070                  b123=zero
1071                end if
1072                b11=b11+2._dp*b113
1073                b12_13=b12_13+b123-b133
1074                b21_22=b21_22-b213+b223
1075                b24=b24+2._dp*b243
1076                b31_32=b31_32+b313
1077              end do ! istl
1078 
1079              b11=b11*zi*(1._dp/wnm)*const_esu
1080              b12_13=b12_13*zi*(1._dp/wnm)*const_esu
1081              b24=(b24+b231)*zi*(1._dp/(wnm**3))*const_esu
1082              b21_22=(b21_22)*zi*(1._dp/(wnm**3))*const_esu
1083              b31_32=(b31_32-b331)*zi*(1._dp/(wmn**3))*const_esu*0.5_dp
1084              ! calculate over the desired energy mesh and sum over k-points
1085              do iw=1,nmesh
1086                w=(iw-1)*de+idel
1087                inter2w(iw)=inter2w(iw)+(ks_ebands%wtk(ik)*(b11/(wmn-2._dp*w))) ! Inter(2w) from chi
1088                inter1w(iw)=inter1w(iw)+(ks_ebands%wtk(ik)*(b12_13/(wmn-w))) ! Inter(1w) from chi
1089                intra2w(iw)=intra2w(iw)+(ks_ebands%wtk(ik)*(b24/(wmn-2._dp*w))) ! Intra(2w) from eta
1090                intra1w(iw)=intra1w(iw)+(ks_ebands%wtk(ik)*((b21_22)/(wmn-w))) ! Intra(1w) from eta
1091                intra1wS(iw)=intra1wS(iw)+(ks_ebands%wtk(ik)*((b31_32)/(wmn-w))) ! Intra(1w) from sigma
1092              end do
1093            end if
1094          end do ! istm
1095        end if
1096      end do ! istn
1097    end do  ! spins
1098  end do ! k-points
1099 
1100  call xmpi_sum(inter2w,comm,ierr)
1101  call xmpi_sum(inter1w,comm,ierr)
1102  call xmpi_sum(intra2w,comm,ierr)
1103  call xmpi_sum(intra1w,comm,ierr)
1104  call xmpi_sum(intra1wS,comm,ierr)
1105  call xmpi_min(my_emin,emin,comm,ierr)
1106  call xmpi_max(my_emax,emax,comm,ierr)
1107 
1108  if (my_rank == master) then
1109    ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
1110 
1111    if (ncid /= nctk_noid) then
1112      start4 = [1, 1, icomp, itemp]
1113      count4 = [2, nmesh, 1, 1]
1114      ABI_MALLOC(chi2tot, (nmesh))
1115      chi2tot = inter2w + inter1w + intra2w + intra1w + intra1wS
1116 #ifdef HAVE_NETCDF
1117      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter2w"), c2r(inter2w), start=start4, count=count4))
1118      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_inter1w"), c2r(inter1w), start=start4, count=count4))
1119      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra2w"), c2r(intra2w), start=start4, count=count4))
1120      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1w"), c2r(intra1w), start=start4, count=count4))
1121      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_intra1wS"), c2r(intra1wS), start=start4, count=count4))
1122      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "shg_chi2tot"), c2r(chi2tot), start=start4, count=count4))
1123 #endif
1124      ABI_FREE(chi2tot)
1125    end if
1126 
1127    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
1128      ABI_ERROR(msg)
1129    end if
1130    if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
1131      ABI_ERROR(msg)
1132    end if
1133    if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
1134      ABI_ERROR(msg)
1135    end if
1136    if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
1137      ABI_ERROR(msg)
1138    end if
1139    if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
1140      ABI_ERROR(msg)
1141    end if
1142    if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then
1143      ABI_ERROR(msg)
1144    end if
1145    if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then
1146      ABI_ERROR(msg)
1147    end if
1148    ! write headers
1149    write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1150    write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
1151    write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
1152    write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
1153    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1154    write(fout1, '(a)' )' # Energy      Tot-Im Chi(-2w,w,w)  Tot-Im Chi(-2w,w,w)'
1155    write(fout1, '(a)' )' # eV          *10^-7 esu        *10^-12 m/V SI units '
1156    write(fout1, '(a)' )' # '
1157 
1158    write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
1159    write(fout2, '(a,es16.6)') ' #tolerance:',tol
1160    write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1161    write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1162    write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1163    write(fout2, '(a)')' # Energy      Tot-Re Chi(-2w,w,w)  Tot-Re Chi(-2w,w,w)'
1164    write(fout2, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
1165    write(fout2, '(a)')' # '
1166 
1167    write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1168    write(fout3, '(a,es16.6)') ' #tolerance:',tol
1169    write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1170    write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1171    write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1172    write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
1173    write(fout3, '(a)')' # in esu'
1174    write(fout3, '(a)')' # '
1175 
1176    write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1177    write(fout4, '(a,es16.6)') ' #tolerance:',tol
1178    write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1179    write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1180    write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1181    write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
1182    write(fout4, '(a)')' # in esu'
1183    write(fout4, '(a)')' # '
1184 
1185    write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1186    write(fout5, '(a,es16.6)') ' #tolerance:',tol
1187    write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1188    write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1189    write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1190    write(fout5, '(a)')' # Energy(eV)  |TotChi(-2w,w,w)|   |Tot Chi(-2w,w,w)|'
1191    write(fout5, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
1192    write(fout5, '(a)')' # '
1193 
1194    write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1195    write(fout6, '(a,es16.6)') ' #tolerance:',tol
1196    write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1197    write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1198    write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1199    write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1200    write(fout6, '(a)')' # in esu'
1201    write(fout6, '(a)')' # '
1202 
1203    write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3
1204    write(fout7, '(a,es16.6)') ' #tolerance:',tol
1205    write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
1206    write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
1207    write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
1208    write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
1209    write(fout7, '(a)')' # in esu'
1210    write(fout7, '(a)')' # '
1211 
1212    totim=zero
1213    totre=zero
1214    totabs=zero
1215    do iw=2,nmesh
1216      ene=(iw-1)*de
1217      ene=ene*ha2ev
1218 
1219      totim=aimag(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1220      write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1221      totim=zero
1222 
1223      totre=dble(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1224      write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1225      totre=zero
1226 
1227      write(fout3,'(f15.6,4es15.6)') ene,aimag(inter2w(iw))/1.d-7,      &
1228      aimag(inter1w(iw))/1.d-7,aimag(intra2w(iw))/1.d-7, aimag(intra1w(iw)+intra1wS(iw))/1.d-7
1229 
1230      write(fout4,'(f15.6,4es15.6)') ene,dble(inter2w(iw))/1.d-7,       &
1231      dble(inter1w(iw))/1.d-7,dble(intra2w(iw))/1.d-7,dble(intra1w(iw)+intra1wS(iw))/1.d-7
1232 
1233      totabs=abs(inter2w(iw)+inter1w(iw)+intra2w(iw)+intra1w(iw)+intra1wS(iw))/1.d-7
1234      write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
1235      totabs=zero
1236 
1237      write(fout6,'(f15.6,4es15.6)') ene,aimag(inter2w(iw)+inter1w(iw))/1.d-7,      &
1238      aimag(intra2w(iw)+intra1w(iw))/1.d-7,aimag(intra1wS(iw))/1.d-7
1239 
1240      write(fout7,'(f15.6,4es15.6)') ene,dble(inter2w(iw)+inter1w(iw))/1.d-7,       &
1241      dble(intra2w(iw)+intra1w(iw))/1.d-7,dble(intra1wS(iw))/1.d-7
1242    end do
1243 
1244    close(fout1)
1245    close(fout2)
1246    close(fout3)
1247    close(fout4)
1248    close(fout5)
1249    close(fout6)
1250    close(fout7)
1251    ! print information
1252    write(std_out,*) ' '
1253    write(std_out,*) 'information about calculation just performed:'
1254    write(std_out,*) ' '
1255    write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of second order susceptibility'
1256    write(std_out,*) 'tolerance:',tol
1257    if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
1258    write(std_out,*) 'broadening:',brod,'Hartree'
1259    if (brod.gt.0.009) then
1260      write(std_out,*) ' '
1261      write(std_out,*) 'ATTENTION: broadening is quite high'
1262      write(std_out,*) ' '
1263    else if (brod.gt.0.015) then
1264      write(std_out,*) ' '
1265      write(std_out,*) 'ATTENTION: broadening is too high'
1266      write(std_out,*) ' '
1267    end if
1268    write(std_out,*) 'scissors shift:',sc,'Hartree'
1269    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree'
1270  end if
1271 
1272  ! deallocate local arrays
1273  ABI_FREE(px)
1274  ABI_FREE(py)
1275  ABI_FREE(pz)
1276  ABI_FREE(inter2w)
1277  ABI_FREE(inter1w)
1278  ABI_FREE(intra2w)
1279  ABI_FREE(intra1w)
1280  ABI_FREE(intra1wS)
1281  ABI_FREE(delta)
1282 
1283 end subroutine nlinopt

m_optic_tools/nonlinopt [ Functions ]

[ Top ] [ m_optic_tools ] [ Functions ]

NAME

 nonlinopt

FUNCTION

 Compute the frequency dependent nonlinear electro-optic susceptibility for semiconductors

INPUTS

  icomp=Sequential index associated to computed tensor components (used for netcdf output)
  itemp=Temperature index (used for netcdf output)
  nband_sum=Number of bands included in the sum. Must be <= mband
  pmat(mband,mband,nkpt,3,nsppol) = momentum matrix elements in cartesian coordinates(complex)
  v1,v2,v3 = desired component of the dielectric function(integer) 1=x,2=y,3=z
  nmesh = desired number of energy mesh points(integer)
  de = desired step in energy(real); nmesh*de=maximum energy for plotting
  sc = scissors shift in Ha(real)
  brod = broadening in Ha(real)
  tol = tolerance:how close to the singularity exact exact is calculated(real)
  fnam=root for filenames that will contain the output  :
   fnam1=trim(fnam)//'-ChiTotIm.out'
   fnam2=trim(fnam)//'-ChiTotRe.out'
   fnam3=trim(fnam)//'-ChiIm.out'
   fnam4=trim(fnam)//'-ChiRe.out'
   fnam5=trim(fnam)//'-ChiAbs.out'

OUTPUT

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

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

SOURCE

1858 subroutine nonlinopt(icomp, itemp, nband_sum, cryst, ks_ebands, &
1859                       pmat, v1, v2, v3, nmesh, de, sc, brod, tol, fnam, do_antiresonant, ncid, comm)
1860 
1861 !Arguments ------------------------------------
1862 integer, intent(in) :: icomp, itemp, nband_sum, ncid
1863 type(crystal_t),intent(in) :: cryst
1864 type(ebands_t),intent(in) :: ks_ebands
1865 complex(dpc), intent(in) :: pmat(ks_ebands%mband, ks_ebands%mband, ks_ebands%nkpt, 3, ks_ebands%nsppol)
1866 integer, intent(in) :: v1, v2, v3
1867 integer, intent(in) :: nmesh
1868 integer, intent(in) :: comm
1869 real(dp), intent(in) :: de, sc, brod, tol
1870 character(len=*), intent(in) :: fnam
1871 logical, intent(in) :: do_antiresonant
1872 
1873 !Local variables -------------------------
1874 integer :: iw,i,j,k,lx,ly,lz,mband
1875 integer :: isp,isym,ik,ist1,istl,istn,istm
1876 real(dp) :: ha2ev
1877 real(dp) :: ene,totre,totabs,totim
1878 real(dp) :: el,en,em
1879 real(dp) :: emin,emax, my_emin,my_emax
1880 real(dp) :: const_esu,const_au,au2esu
1881 real(dp) :: wmn,wnm,wln,wnl,wml,wlm !, t1
1882 complex(dpc) :: idel,w,zi
1883 character(len=fnlen) :: fnam1,fnam2,fnam3,fnam4,fnam5,fnam6,fnam7
1884 ! local allocatable arrays
1885  integer :: start4(4),count4(4)
1886  real(dp) :: s(3,3),sym(3,3,3)
1887  integer :: istp
1888  real(dp) :: ep, wmp, wpn, wtk
1889  real(dp), allocatable :: enk(:) ! (n) = \omega_n(k), with scissor included !
1890  real(dp) :: fn, fm, fl, fnm, fnl, fml, fln, flm
1891  complex(dpc), allocatable :: delta(:,:,:) ! (m,n,a) = \Delta_{mn}^{a}
1892  complex(dpc), allocatable :: rmna(:,:,:) ! (m,n,a) = r_{mn}^{a}
1893  complex(dpc), allocatable :: rmnbc(:,:,:,:) ! (m,n,b,c) = r^b_{mn;c}(k)
1894  complex(dpc), allocatable :: roverw(:,:,:,:) ! (m,n,b,c) = [r^b_{mn}(k)/w_{mn(k)];c
1895  complex(dpc), allocatable :: chiw(:), chi2w(:) ! \chi_{II}^{abc}(-\omega,\omega,0)
1896  complex(dpc), allocatable :: etaw(:), eta2w(:) ! \eta_{II}^{abc}(-\omega,\omega,0)
1897  complex(dpc), allocatable :: sigmaw(:) ! \frac{i}{\omega} \sigma_{II}^{abc}(-\omega,\omega,0)
1898  complex(dpc) :: num1, num2, den1, den2, term1, term2
1899  complex(dpc) :: chi1, chi2_1, chi2_2
1900  complex(dpc), allocatable :: chi2(:) ! Second term that depends on the frequency ! (omega)
1901  complex(dpc), allocatable :: eta1(:) ! Second term that depends on the frequency ! (omega)
1902  complex(dpc), allocatable :: chi2tot(:)
1903  complex(dpc) :: eta1_1, eta1_2, eta2_1, eta2_2
1904  complex(dpc) :: sigma2_1, sigma1
1905  complex(dpc), allocatable :: symrmn(:,:,:) ! (m,l,n) = 1/2*(rml^b rln^c+rml^c rln^b)
1906  complex(dpc) :: symrmnl(3,3), symrlmn(3,3), symrmln(3,3)
1907 !Parallelism
1908  integer :: my_rank, nproc
1909  integer,parameter :: master = 0
1910  integer :: ierr
1911  integer :: my_k1, my_k2
1912  character(500) :: msg
1913  integer :: fout1,fout2,fout3,fout4,fout5,fout6,fout7
1914 
1915 ! *********************************************************************
1916 
1917  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
1918 
1919 !calculate the constant
1920  zi=(0._dp,1._dp)
1921  idel=zi*brod
1922  const_au=-2._dp/(cryst%ucvol*dble(cryst%nsym))
1923  !const_au=-2._dp/(cryst%ucvol)
1924  au2esu=5.8300348177d-8
1925  const_esu=const_au*au2esu
1926  ha2ev=13.60569172*2._dp
1927 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1928 !5.8300348177d-8 : au2esu : bohr*c*10^4/4pi*2*ry2ev
1929 !bohr: 5.2917ifc nlinopt.f907E-11
1930 !c: 2.99792458   velocity of sound
1931 !ry2ev: 13.60569172
1932 !au2esu=(5.29177E-11*2.99792458*1.0E4)/(13.60569172*2)
1933 !this const includes (e^3*hbar^3*hbar^3)/(vol*hbar^5*m_e^3)
1934 !mass comes from converting P_mn to r_mn
1935 !hbar^3 comes from converting all frequencies to energies in denominator
1936 !hbar^3 comes from operator for momentum (hbar/i nabla)
1937 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1938 !output file names
1939  fnam1=trim(fnam)//'-ChiSHGTotIm.out'
1940  fnam2=trim(fnam)//'-ChiSHGTotRe.out'
1941  fnam3=trim(fnam)//'-ChiSHGIm.out'
1942  fnam4=trim(fnam)//'-ChiSHGRe.out'
1943  fnam5=trim(fnam)//'-ChiSHGAbs.out'
1944  fnam6=trim(fnam)//'-ChiSHGImDec.out'
1945  fnam7=trim(fnam)//'-ChiSHGReDec.out'
1946 
1947  ! If there exists inversion symmetry exit with a mesg.
1948  if (cryst%idx_spatial_inversion() /= 0) then
1949    write(std_out,*) '-----------------------------------------'
1950    write(std_out,*) '    the crystal has inversion symmetry   '
1951    write(std_out,*) '    the nl electro-optical susceptibility'
1952    write(std_out,*) '    is zero                              '
1953    write(std_out,*) '-----------------------------------------'
1954    ABI_ERROR("Aborting now")
1955  end if
1956 
1957  ! check polarisation
1958  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
1959    write(std_out,*) '---------------------------------------------'
1960    write(std_out,*) '    Error in nonlinopt:                        '
1961    write(std_out,*) '    the polarisation directions incorrect    '
1962    write(std_out,*) '    1=x,  2=y  and 3=z                       '
1963    write(std_out,*) '---------------------------------------------'
1964    ABI_ERROR("Aborting now")
1965  end if
1966 
1967  ! number of energy mesh points
1968  if (nmesh.le.0) then
1969    write(std_out,*) '---------------------------------------------'
1970    write(std_out,*) '    Error in nonlinopt:                        '
1971    write(std_out,*) '    number of energy mesh points incorrect   '
1972    write(std_out,*) '    number has to be integer greater than 0  '
1973    write(std_out,*) '    nmesh*de = max energy for calculation    '
1974    write(std_out,*) '---------------------------------------------'
1975    ABI_ERROR("Aborting now")
1976  end if
1977 
1978  ! step in energy
1979  if (de.le.zero) then
1980    write(std_out,*) '---------------------------------------------'
1981    write(std_out,*) '    Error in nonlinopt:                        '
1982    write(std_out,*) '    energy step is incorrect                 '
1983    write(std_out,*) '    number has to real greater than 0.0      '
1984    write(std_out,*) '    nmesh*de = max energy for calculation    '
1985    write(std_out,*) '---------------------------------------------'
1986    ABI_ERROR("Aborting now")
1987  end if
1988 
1989  ! broadening
1990  if (brod.gt.0.009) then
1991    write(std_out,*) '---------------------------------------------'
1992    write(std_out,*) '    ATTENTION: broadening is quite high      '
1993    write(std_out,*) '    ideally should be less than 0.005        '
1994    write(std_out,*) '---------------------------------------------'
1995  else if (brod.gt.0.015) then
1996    write(std_out,*) '----------------------------------------'
1997    write(std_out,*) '    ATTENTION: broadening is too high   '
1998    write(std_out,*) '    ideally should be less than 0.005   '
1999    write(std_out,*) '----------------------------------------'
2000  end if
2001 
2002  ! tolerance
2003  if (tol.gt.0.006) then
2004    write(std_out,*) '----------------------------------------'
2005    write(std_out,*) '    ATTENTION: tolerance is too high    '
2006    write(std_out,*) '    ideally should be less than 0.004   '
2007    write(std_out,*) '----------------------------------------'
2008  end if
2009 
2010  ! allocate local arrays
2011  mband = ks_ebands%mband
2012  ABI_CHECK(nband_sum <= mband, "nband_sum <= mband")
2013  ABI_MALLOC(enk, (mband))
2014  ABI_MALLOC(delta, (mband, mband, 3))
2015  ABI_MALLOC(rmnbc, (mband, mband, 3, 3))
2016  ABI_MALLOC(roverw, (mband, mband, 3, 3))
2017  ABI_MALLOC(rmna, (mband, mband, 3))
2018  ABI_MALLOC(chiw, (nmesh))
2019  ABI_MALLOC(etaw, (nmesh))
2020  ABI_MALLOC(chi2w, (nmesh))
2021  ABI_MALLOC(eta2w, (nmesh))
2022  ABI_MALLOC(sigmaw, (nmesh))
2023  ABI_MALLOC(chi2, (nmesh))
2024  ABI_MALLOC(eta1, (nmesh))
2025  ABI_MALLOC(symrmn, (mband, mband, mband))
2026 
2027  ! generate the symmetrizing tensor
2028  sym = zero
2029  do isym=1,cryst%nsym
2030    s(:,:)=cryst%symrel_cart(:,:,isym)
2031    do i=1,3
2032      do j=1,3
2033        do k=1,3
2034          sym(i,j,k)=sym(i,j,k)+(s(i,v1)*s(j,v2)*s(k,v3))
2035        end do
2036      end do
2037    end do
2038  end do
2039 
2040  ! initialise
2041  delta = zero
2042  rmnbc = zero
2043  chiw = zero
2044  chi2w = zero
2045  chi2 = zero
2046  etaw = zero
2047  eta2w = zero
2048  sigmaw = zero
2049  my_emin=HUGE(one)
2050  my_emax=-HUGE(one)
2051 
2052  ! Split work
2053  call xmpi_split_work(ks_ebands%nkpt, comm, my_k1, my_k2)
2054 
2055 ! loop over kpts
2056  do ik=my_k1,my_k2
2057    write(std_out,*) "P-",my_rank,": ",ik,'of ', ks_ebands%nkpt
2058    do isp=1,ks_ebands%nsppol
2059      ! Calculate the scissor corrected energies and the energy window
2060      do ist1=1,nband_sum
2061        en = ks_ebands%eig(ist1,ik,isp)
2062        my_emin=min(my_emin,en)
2063        my_emax=max(my_emax,en)
2064        if(en > ks_ebands%fermie) then
2065          en = en + sc
2066        end if
2067        enk(ist1) = en
2068      end do
2069 
2070      ! calculate \Delta_nm and r_mn^a
2071      do istn=1,nband_sum
2072        en = enk(istn)
2073        do istm=1,nband_sum
2074          em = enk(istm)
2075          wmn = em - en
2076          delta(istn,istm,1:3)=pmat(istn,istn,ik,1:3,isp)-pmat(istm,istm,ik,1:3,isp)
2077          if(abs(wmn) < tol) then
2078            rmna(istm,istn,1:3) = zero
2079          else
2080            rmna(istm,istn,1:3)=pmat(istm,istn,ik,1:3,isp)/wmn
2081          end if
2082        end do
2083      end do
2084 
2085      ! calculate \r^b_mn;c
2086      do istm=1,nband_sum
2087        em = enk(istm)
2088        do istn=1,nband_sum
2089          en = enk(istn)
2090          wmn = em - en
2091          if (abs(wmn) < tol) then ! Degenerate energies
2092            rmnbc(istm,istn,:,:) = zero
2093            roverw(istm,istn,:,:) = zero
2094            cycle
2095          end if
2096          do ly = 1,3
2097            do lz = 1,3
2098              num1 = (rmna(istm,istn,ly)*delta(istm,istn,lz))+(rmna(istm,istn,lz)*delta(istm,istn,ly))
2099              den1 = wmn
2100              term1 = num1/den1
2101              term2 = zero
2102              do istp=1,nband_sum
2103                ep = enk(istp)
2104                wmp = em - ep
2105                wpn = ep - en
2106                num2 = (wmp*rmna(istm,istp,ly)*rmna(istp,istn,lz))-(wpn*rmna(istm,istp,lz)*rmna(istp,istn,ly))
2107                den2 = wmn
2108                term2 = term2 + (num2/den2)
2109              end do
2110              rmnbc(istm,istn,ly,lz) = -term1-(zi*term2)
2111              roverw(istm,istn,ly,lz) = (rmnbc(istm,istn,ly,lz)/wmn) - (rmna(istm,istn,ly)/(wmn**2))*delta(istm,istn,lz)
2112            end do
2113          end do
2114        end do
2115      end do
2116 
2117      ! initialise the factors
2118      ! start the calculation
2119      do istn=1,nband_sum
2120        en=enk(istn)
2121        fn=ks_ebands%occ(istn,ik,isp)
2122        if(do_antiresonant .and. en .ge. ks_ebands%fermie) then
2123          cycle
2124        end if
2125        do istm=1,nband_sum
2126          em=enk(istm)
2127          if (do_antiresonant .and. em .le. ks_ebands%fermie) then
2128            cycle
2129          end if
2130          wmn=em-en
2131          wnm=-wmn
2132          fm = ks_ebands%occ(istm,ik,isp)
2133          fnm = fn - fm
2134          if(abs(wmn) > tol) then
2135            chi1 = zero
2136            chi2(:) = zero
2137            chi2_1 = zero
2138            chi2_2 = zero
2139            eta1(:) = zero
2140            eta1_1 = zero
2141            eta1_2 = zero
2142            eta2_1 = zero
2143            eta2_2 = zero
2144            sigma1 = zero
2145            sigma2_1 = zero
2146            ! Three band terms
2147            do istl=1,nband_sum
2148              el=enk(istl)
2149              fl = ks_ebands%occ(istl,ik,isp)
2150              wlm = el-em
2151              wln = el-en
2152              wnl = -wln
2153              wml = em-el
2154              fnl = fn-fl
2155              fml = fm-fl
2156              flm = -fml
2157              fln = -fnl
2158              do ly = 1,3
2159                do lz = 1,3
2160                  symrmnl(ly,lz) = 0.5_dp*(rmna(istm,istn,ly)*rmna(istn,istl,lz)+rmna(istm,istn,lz)*rmna(istn,istl,ly))
2161                  symrlmn(ly,lz) = 0.5_dp*(rmna(istl,istm,ly)*rmna(istm,istn,lz)+rmna(istl,istm,lz)*rmna(istm,istn,ly))
2162                  symrmln(ly,lz) = 0.5_dp*(rmna(istm,istl,ly)*rmna(istl,istn,lz)+rmna(istm,istl,lz)*rmna(istl,istn,ly))
2163                end do
2164              end do
2165 
2166              do lx = 1,3
2167                do ly = 1,3
2168                  do lz = 1,3
2169                    sigma1 = sigma1 + sym(lx,ly,lz)*(wnl*rmna(istl,istm,lx)*symrmnl(ly,lz)-wlm*rmna(istn,istl,lx)*symrlmn(ly,lz))
2170                    eta2_2 = eta2_2 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*symrmln(ly,lz)*(wml-wln)
2171                    if(abs(wln-wml) > tol) then
2172                      chi1 = chi1 + sym(lx,ly,lz)*(rmna(istn,istm,lx)*symrmln(ly,lz))/(wln-wml)
2173                    end if
2174                    eta1_1 = eta1_1 + sym(lx,ly,lz)*wln*rmna(istn,istl,lx)*symrlmn(ly,lz)
2175                    eta1_2 = eta1_2 - sym(lx,ly,lz)*wml*rmna(istl,istm,lx)*symrmnl(ly,lz)
2176                    if(abs(wnl-wmn) > tol) then
2177                      chi2_1 = chi2_1 - sym(lx,ly,lz)*(fnm*rmna(istl,istm,lx)*symrmnl(ly,lz)/(wnl-wmn))
2178                    end if
2179                    if(abs(wmn-wlm) > tol) then
2180                      chi2_2 = chi2_2 - sym(lx,ly,lz)*(fnm*rmna(istn,istl,lx)*symrlmn(ly,lz)/(wmn-wlm))
2181                    end if
2182                  end do
2183                end do
2184              end do
2185            end do
2186 
2187            ! Two band terms
2188            eta2_1 = zero
2189            sigma2_1 = zero
2190            do lx = 1,3
2191              do ly = 1,3
2192                do lz = 1,3
2193                  eta2_1 = eta2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp &
2194                     *(delta(istm,istn,ly)*rmna(istm,istn,lz)+delta(istm,istn,lz)*rmna(istm,istn,ly))
2195                  ! Correct version (Sipe 1993)
2196                  sigma2_1 = sigma2_1 + sym(lx,ly,lz)*fnm*rmna(istn,istm,lx)*0.5_dp &
2197                     *(rmna(istm,istn,ly)*delta(istn,istm,lz)+rmna(istm,istn,lz)*delta(istn,istm,ly))
2198 
2199                  ! Incorrect version (Hughes 1996)
2200                  !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))
2201                end do
2202              end do
2203            end do
2204 
2205            ! calculate over the desired energy mesh and sum over k-points
2206            wtk = ks_ebands%wtk(ik)
2207            do iw=1,nmesh
2208              w=(iw-1)*de+idel
2209              chi2w(iw) = chi2w(iw) + zi*wtk*((2.0_dp*fnm*chi1/(wmn-2.0_dp*w)))*const_esu ! Inter(2w) from chi
2210              chiw(iw) = chiw(iw) + zi*wtk*((chi2_1+chi2_2)/(wmn-w))*const_esu ! Inter(w) from chi
2211              eta2w(iw) = eta2w(iw) + zi*wtk*(8.0_dp*(eta2_1/((wmn**2)*(wmn-2.0_dp*w))) &
2212                  + 2.0_dp*eta2_2/((wmn**2)*(wmn-2.0_dp*w)))*const_esu ! Intra(2w) from eta
2213              etaw(iw) = etaw(iw) + zi*wtk*((eta1_1 + eta1_2)*fnm/((wmn**2)*(wmn-w)))*const_esu ! Intra(w) from eta
2214              sigmaw(iw) = sigmaw(iw) + 0.5_dp*zi*wtk*(fnm*sigma1/((wmn**2)*(wmn-w)) &
2215                  + (sigma2_1/((wmn**2)*(wmn-w))))*const_esu ! Intra(1w) from sigma
2216            end do
2217          end if
2218        end do ! end loop over istn and istm
2219      end do
2220    end do ! spins
2221  end do ! k-points
2222 
2223  ! Collect info among the nodes
2224  call xmpi_min(my_emin,emin,comm,ierr)
2225  call xmpi_max(my_emax,emax,comm,ierr)
2226 
2227  call xmpi_sum(chiw,comm,ierr)
2228  call xmpi_sum(etaw,comm,ierr)
2229  call xmpi_sum(chi2w,comm,ierr)
2230  call xmpi_sum(eta2w,comm,ierr)
2231  call xmpi_sum(sigmaw,comm,ierr)
2232 
2233  ! Master writes the output
2234  if (my_rank == master) then
2235 
2236    if (ncid /= nctk_noid) then
2237      start4 = [1, 1, icomp, itemp]
2238      count4 = [2, nmesh, 1, 1]
2239      ABI_MALLOC(chi2tot, (nmesh))
2240      chi2tot = chiw + chi2w + etaw + eta2w + sigmaw
2241 #ifdef HAVE_NETCDF
2242      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2tot"), c2r(chi2tot), start=start4, count=count4))
2243      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chiw"), c2r(chiw), start=start4, count=count4))
2244      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_etaw"), c2r(etaw), start=start4, count=count4))
2245      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_chi2w"), c2r(chi2w), start=start4, count=count4))
2246      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_eta2w"), c2r(eta2w), start=start4, count=count4))
2247      NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "leo2_sigmaw"), c2r(sigmaw), start=start4, count=count4))
2248 #endif
2249      ABI_FREE(chi2tot)
2250    end if
2251 
2252    ! write output in SI units and esu (esu to SI(m/v)=(value_esu)*(4xpi)/30000)
2253    if (open_file(fnam1,msg,newunit=fout1,action='WRITE',form='FORMATTED') /= 0) then
2254      ABI_ERROR(msg)
2255    end if
2256    if (open_file(fnam2,msg,newunit=fout2,action='WRITE',form='FORMATTED') /= 0) then
2257      ABI_ERROR(msg)
2258    end if
2259    if (open_file(fnam3,msg,newunit=fout3,action='WRITE',form='FORMATTED') /= 0) then
2260      ABI_ERROR(msg)
2261    end if
2262    if (open_file(fnam4,msg,newunit=fout4,action='WRITE',form='FORMATTED') /= 0) then
2263      ABI_ERROR(msg)
2264    end if
2265    if (open_file(fnam5,msg,newunit=fout5,action='WRITE',form='FORMATTED') /= 0) then
2266      ABI_ERROR(msg)
2267    end if
2268    if (open_file(fnam6,msg,newunit=fout6,action='WRITE',form='FORMATTED') /= 0) then
2269      ABI_ERROR(msg)
2270    end if
2271    if (open_file(fnam7,msg,newunit=fout7,action='WRITE',form='FORMATTED') /= 0) then
2272      ABI_ERROR(msg)
2273    end if
2274 
2275    ! write headers
2276    write(fout1, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2277    write(fout1, '(a,es16.6)' ) ' #tolerance:',tol
2278    write(fout1, '(a,es16.6,a)' ) ' #broadening:',brod,'Ha'
2279    write(fout1, '(a,es16.6,a)' ) ' #scissors shift:',sc,'Ha'
2280    write(fout1, '(a,es16.6,a,es16.6,a)' ) ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2281    write(fout1, '(a)' )' # Energy      Tot-Im Chi(-w,w,0)  Tot-Im Chi(-w,w,0)'
2282    write(fout1, '(a)' )' # eV          *10^-7 esu        *10^-12 m/V SI units '
2283    write(fout1, '(a)' )' # '
2284 
2285    write(fout2, '(a,3i3)' ) ' #calculated the component:',v1,v2,v3
2286    write(fout2, '(a,es16.6)') ' #tolerance:',tol
2287    write(fout2, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2288    write(fout2, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2289    write(fout2, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2290    write(fout2, '(a)')' # Energy      Tot-Re Chi(-w,w,0)  Tot-Re Chi(-w,w,0)'
2291    write(fout2, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
2292    write(fout2, '(a)')' # '
2293 
2294    write(fout3, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2295    write(fout3, '(a,es16.6)') ' #tolerance:',tol
2296    write(fout3, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2297    write(fout3, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2298    write(fout3, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2299    write(fout3, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
2300    write(fout3, '(a)')' # in esu'
2301    write(fout3, '(a)')' # '
2302 
2303    write(fout4, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2304    write(fout4, '(a,es16.6)') ' #tolerance:',tol
2305    write(fout4, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2306    write(fout4, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2307    write(fout4, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2308    write(fout4, '(a)')' # Energy(eV) Inter(2w) inter(1w) intra(2w) intra(1w)'
2309    write(fout4, '(a)')' # in esu'
2310    write(fout4, '(a)')' # '
2311 
2312    write(fout5, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2313    write(fout5, '(a,es16.6)') ' #tolerance:',tol
2314    write(fout5, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2315    write(fout5, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2316    write(fout5, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2317    write(fout5, '(a)')' # Energy(eV)  |TotChi(-w,w,0)|   |Tot Chi(-w,w,0)|'
2318    write(fout5, '(a)')' # eV          *10^-7 esu        *10^-12 m/V SI units '
2319    write(fout5, '(a)')' # '
2320 
2321    write(fout6, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2322    write(fout6, '(a,es16.6)') ' #tolerance:',tol
2323    write(fout6, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2324    write(fout6, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2325    write(fout6, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2326    write(fout6, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2327    write(fout6, '(a)')' # in esu'
2328    write(fout6, '(a)')' # '
2329 
2330    write(fout7, '(a,3i3)') ' #calculated the component:',v1,v2,v3
2331    write(fout7, '(a,es16.6)') ' #tolerance:',tol
2332    write(fout7, '(a,es16.6,a)') ' #broadening:',brod,'Ha'
2333    write(fout7, '(a,es16.6,a)') ' #scissors shift:',sc,'Ha'
2334    write(fout7, '(a,es16.6,a,es16.6,a)') ' #energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Ha'
2335    write(fout7, '(a)')' # Energy(eV) Chi(w) Eta(w) Sigma(w)'
2336    write(fout7, '(a)')' # in esu'
2337    write(fout7, '(a)')' # '
2338 
2339    totim=zero
2340    totre=zero
2341    totabs=zero
2342    do iw=2,nmesh
2343      ene=(iw-1)*de
2344      ene=ene*ha2ev
2345 
2346      totim=aimag(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7
2347      write(fout1,'(f15.6,2es15.6)') ene,totim,totim*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2348      totim=zero
2349 
2350      totre=dble(chiw(iw)+chi2w(iw)+eta2w(iw)+etaw(iw)+sigmaw(iw))/1.d-7
2351      write(fout2,'(f15.6,2es15.6)') ene,totre,totre*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2352      totre=zero
2353 
2354      write(fout3,'(f15.6,4es15.6)') ene,aimag(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7,     &
2355      aimag(eta2w(iw))/1.d-7,aimag(etaw(iw))/1.d-7+aimag(sigmaw(iw))/1.d-7
2356 
2357      write(fout4,'(f15.6,4es15.6)') ene,dble(chi2w(iw))/1.d-7,aimag(chiw(iw))/1.d-7,       &
2358      dble(eta2w(iw))/1.d-7,dble(etaw(iw))/1.d-7+dble(sigmaw(iw))/1.d-7
2359 
2360      totabs=abs(chiw(iw)+chi2w(iw)+etaw(iw)+eta2w(iw)+sigmaw(iw))/1.d-7
2361      write(fout5,'(f15.6,2es15.6)') ene,totabs,totabs*4._dp*pi*(1._dp/30000._dp)*(1._dp/1.d-5)
2362      totabs=zero
2363 
2364      write(fout6,'(f15.6,4es15.6)') ene,aimag(chi2w(iw)+chiw(iw))/1.d-7,      &
2365      aimag(eta2w(iw)+etaw(iw))/1.d-7,aimag(sigmaw(iw))/1.d-7
2366 
2367      write(fout7,'(f15.6,4es15.6)') ene,dble(chi2w(iw)+chiw(iw))/1.d-7,       &
2368      dble(eta2w(iw)+etaw(iw))/1.d-7,dble(sigmaw(iw))/1.d-7
2369    end do
2370 
2371    close(fout1)
2372    close(fout2)
2373    close(fout3)
2374    close(fout4)
2375    close(fout5)
2376    close(fout6)
2377    close(fout7)
2378 
2379    ! print information
2380    write(std_out,*) ' '
2381    write(std_out,*) 'information about calculation just performed:'
2382    write(std_out,*) ' '
2383    write(std_out,*) 'calculated the component:',v1,v2,v3 ,'of the nonlinear electro-optical susceptibility'
2384    write(std_out,*) 'tolerance:',tol
2385    if (tol.gt.0.008) write(std_out,*) 'ATTENTION: tolerance is too high'
2386    write(std_out,*) 'broadening:',brod,'Hartree'
2387    if (brod.gt.0.009) then
2388      write(std_out,*) ' '
2389      write(std_out,*) 'ATTENTION: broadening is quite high'
2390      write(std_out,*) ' '
2391    else if (brod.gt.0.015) then
2392      write(std_out,*) ' '
2393      write(std_out,*) 'ATTENTION: broadening is too high'
2394      write(std_out,*) ' '
2395    end if
2396    write(std_out,*) 'scissors shift:',sc,'Hartree'
2397    write(std_out,*) 'energy window:',(emax-emin)*ha2ev,'eV',(emax-emin),'Hartree'
2398 
2399  end if
2400 
2401  ! deallocate local arrays
2402  ABI_FREE(enk)
2403  ABI_FREE(delta)
2404  ABI_FREE(rmnbc)
2405  ABI_FREE(roverw)
2406  ABI_FREE(rmna)
2407  ABI_FREE(chiw)
2408  ABI_FREE(chi2w)
2409  ABI_FREE(chi2)
2410  ABI_FREE(etaw)
2411  ABI_FREE(eta1)
2412  ABI_FREE(symrmn)
2413  ABI_FREE(eta2w)
2414  ABI_FREE(sigmaw)
2415 
2416 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

SOURCE

 87 subroutine pmat2cart(eigen11, eigen12, eigen13, mband, nkpt, nsppol, pmat, rprimd)
 88 
 89 !Arguments -----------------------------------------------
 90 !scalars
 91  integer,intent(in) :: mband,nkpt,nsppol
 92 !arrays
 93  real(dp),intent(in) :: eigen11(2,mband,mband,nkpt,nsppol)
 94  real(dp),intent(in) :: eigen12(2,mband,mband,nkpt,nsppol)
 95  real(dp),intent(in) :: eigen13(2,mband,mband,nkpt,nsppol),rprimd(3,3)
 96 !no_abirules
 97  complex(dpc),intent(out) :: pmat(mband,mband,nkpt,3,nsppol)
 98 
 99 !Local variables -----------------------------------------
100 !scalars
101  integer :: iband1,iband2,ikpt,isppol
102 !arrays
103  real(dp) :: rprim(3,3)
104 
105 ! *************************************************************************
106 
107  !rescale the rprim
108  rprim(:,:) = rprimd(:,:) / two_pi
109 
110  do isppol=1,nsppol
111    do ikpt=1,nkpt
112      do iband1=1,mband
113        do iband2=1,mband
114          pmat(iband2,iband1,ikpt,:,isppol) =  &
115           rprim(:,1)*cmplx(eigen11(1,iband2,iband1,ikpt,isppol),eigen11(2,iband2,iband1,ikpt,isppol),kind=dp) &
116          +rprim(:,2)*cmplx(eigen12(1,iband2,iband1,ikpt,isppol),eigen12(2,iband2,iband1,ikpt,isppol),kind=dp) &
117          +rprim(:,3)*cmplx(eigen13(1,iband2,iband1,ikpt,isppol),eigen13(2,iband2,iband1,ikpt,isppol),kind=dp)
118        end do
119      end do
120    end do
121  end do
122 
123 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
  fermie = Fermi level
  sc = scissor shift for conduction bands
  eig = ground state eigenvalues

OUTPUT

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

SOURCE

148 subroutine pmat_renorm(fermie, eig, mband, nkpt, nsppol, pmat, sc)
149 
150 !Arguments -----------------------------------------------
151 !scalars
152  integer, intent(in) :: nsppol
153  integer, intent(in) :: nkpt
154  integer, intent(in) :: mband
155  real(dp), intent(in) :: fermie
156  real(dp), intent(in) :: sc
157 !arrays
158  real(dp), intent(in) :: eig(mband,nkpt,nsppol)
159  complex(dpc), intent(inout) :: pmat(mband,mband,nkpt,3,nsppol)
160 
161 !Local variables -----------------------------------------
162 !scalars
163  integer :: iband1,iband2,ikpt,isppol
164  real(dp) :: corec, e1, e2
165 
166 ! *************************************************************************
167 
168  if (abs(sc) < tol8) then
169    call wrtout(std_out,' No scissor shift to be applied. Returning to main optic routine.',"COLL")
170    return
171  end if
172 
173  do isppol=1,nsppol
174    do ikpt=1,nkpt
175      do iband1=1,mband ! valence states
176        e1 = eig(iband1,ikpt,isppol)
177        if (e1 > fermie) cycle
178        do iband2=1,mband ! conduction states
179          e2 = eig(iband2,ikpt,isppol)
180          if (e2 < fermie) cycle
181          corec = (e2+sc-e1)/(e2-e1)
182          pmat(iband2,iband1,ikpt,:,isppol) = corec * pmat(iband2,iband1,ikpt,:,isppol)
183          pmat(iband1,iband2,ikpt,:,isppol) = corec * pmat(iband1,iband2,ikpt,:,isppol)
184        end do
185      end do
186    end do
187  end do
188 
189 end subroutine pmat_renorm