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