TABLE OF CONTENTS
- ABINIT/m_frohlich
- m_frohlich/frohlich_calc_polaronmass
- m_frohlich/frohlich_calc_zpr
- m_frohlich/frohlich_free_el
- m_frohlich/frohlich_free_ph
- m_frohlich/frohlich_init_el
- m_frohlich/frohlich_init_ph
- m_frohlich/frohlich_t
- m_frohlich/frohlichmodel_polaronmass
- m_frohlich/frohlichmodel_zpr
- m_frohlich/lk_hamiltonian
- m_frohlich/mat3inv
ABINIT/m_frohlich [ Modules ]
NAME
m_frohlich
FUNCTION
Description
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (VV, XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_frohlich 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_crystal 28 use m_ebands 29 use m_efmas_defs 30 use m_ifc 31 use m_dtset 32 33 use m_fstrings, only : sjoin, itoa 34 use m_gaussian_quadrature, only : cgqf 35 36 implicit none 37 38 private
m_frohlich/frohlich_calc_polaronmass [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlich_calc_polaronmass
FUNCTION
Description
INPUTS
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
OUTPUT
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
SOURCE
657 subroutine frohlich_calc_polaronmass(self) 658 659 !Arguments ------------------------------------ 660 !scalars 661 class(frohlich_t), intent(inout) :: self 662 !arrays 663 664 !Local variables------------------------------- 665 !scalars 666 integer :: ii 667 integer :: ikdir, nkdir 668 integer :: ik, nkgrid 669 integer :: lwork, info 670 integer :: nu, iband 671 integer :: ixi, nxi 672 integer :: iqdir 673 integer :: sigma 674 real(dp) :: deltak, minefmas 675 real(dp) :: xi, qlen 676 real(dp) :: tmp 677 !arrays 678 real(dp) :: id33(3, 3) 679 real(dp) :: kpt(3), unit_kdir(3, 3) 680 real(dp) :: k_plus_q(3) 681 real(dp), allocatable :: eigenvec(:,:), eigenval(:) 682 real(dp), allocatable :: invefmas(:,:) 683 real(dp), allocatable :: lutt_eigenvec(:,:,:,:), lutt_eigenval(:,:,:) 684 real(dp), allocatable :: lk_h(:,:) 685 real(dp), allocatable :: phfreq(:), investar(:) 686 real(dp), allocatable :: intsum(:,:,:,:) 687 real(dp), allocatable :: zpr(:,:,:), zpr_ddk(:,:) 688 real(dp), allocatable :: intsuminv(:,:) 689 complex(dpc), allocatable :: work(:) 690 691 ! ************************************************************************* 692 693 ! TODO: check if the electronic and phonon parts are initialized 694 ABI_MALLOC(self%invpolmas, (self%ndeg, 3)) 695 696 if (self%kind == 2) then 697 ! Polaron effective masses in the triply-degenerate case 698 ! (Eq. (86) of [Guster2021]) 699 700 ! Initialization of the diagonalization routine 701 ! TODO: this part is common for many routines -> initialize once and store 702 ! necessary variables as local state constants in the datatype 703 ABI_MALLOC(eigenvec, (self%ndeg, self%ndeg)) 704 ABI_MALLOC(eigenval, (self%ndeg)) 705 ABI_MALLOC(work, (3*self%ndeg - 2)) 706 lwork = -1 707 call dsyev('V', 'U', self%ndeg, eigenvec(:,:), self%ndeg, eigenval(:), & 708 work, lwork, info) 709 lwork = int(work(1)) 710 ABI_FREE(work) 711 ABI_MALLOC(work, (lwork)) 712 713 ! Inverse effectiv mass tensor for 100, 110 and 111 directions 714 ! used to obtain the material-dependent characteristic wavevector length 715 nkdir = 3 716 unit_kdir(:, 1) = (/1, 0, 0/) 717 unit_kdir(:, 2) = (/1, 1, 0/)/sqrt(2.0) 718 unit_kdir(:, 3) = (/1, 1, 1/)/sqrt(3.0) 719 ABI_MALLOC(invefmas, (self%ndeg, nkdir)) 720 ABI_MALLOC(self%invefmas, (self%ndeg, nkdir)) 721 722 ! 100: 2A, 2B (two-fold) 723 invefmas(1, 1) = two*self%band_params(1) 724 invefmas(2, 1) = two*self%band_params(2) 725 invefmas(3, 1) = two*self%band_params(2) 726 ! 110: A + B + C, 2B, A + B - C 727 invefmas(1, 2) = self%band_params(1) + self%band_params(2) + & 728 self%band_params(3) 729 invefmas(2, 2) = min(two*self%band_params(2), self%band_params(1) + & 730 self%band_params(2) - self%band_params(3)) 731 invefmas(3, 2) = max(two*self%band_params(2), self%band_params(1) + & 732 self%band_params(2) - self%band_params(3)) 733 ! 111: 2/3*(A + 2B + 2C), 2/3*(A + 2B - C) (two-fold) 734 invefmas(1, 3) = two*(self%band_params(1) + two*self%band_params(2) + & 735 two*self%band_params(3)) / three 736 invefmas(2, 3) = two*(self%band_params(1) + two*self%band_params(2) - & 737 self%band_params(3)) / three 738 invefmas(3, 3) = invefmas(2, 3) 739 740 self%invefmas(:,:) = invefmas(:,:) 741 742 ! Setting up the parameter to obtain the second derivative of the ZPR with 743 ! the finite differences (FD) method (Eqs. (80), (86) of [Guster2021]) 744 ! FD: number of points 745 nkgrid = 3 746 ! FD: material-dependent length scale <--> lowest optical frequency 747 iqdir = 1 748 minefmas = one / maxval(abs(invefmas(:,:))) 749 deltak = sqrt(two*minefmas*self%phfreq_qdir(4, iqdir) / 1000.0) 750 751 ! FD: Luttinger eigenvalues and eigenvectors in the symmetry inequivalent 752 ! cubic directions to be used in the ZPR calculations 753 ABI_MALLOC(lutt_eigenvec, (self%ndeg, self%ndeg, nkgrid, nkdir)) 754 ABI_MALLOC(lutt_eigenval, (self%ndeg, nkgrid, nkdir)) 755 756 lutt_eigenval(:,:,:) = zero 757 do ikdir=1,nkdir 758 do ik=1,nkgrid 759 kpt(:) = (ik - one)*deltak*unit_kdir(:, ikdir) 760 call lk_hamiltonian(self%band_params(:), kpt(:), & 761 lutt_eigenvec(:,:, ik, ikdir)) 762 call dsyev('V', 'U', self%ndeg, lutt_eigenvec(:,:, ik, ikdir), & 763 self%ndeg, lutt_eigenval(:, ik, ikdir), work(:), lwork, info) 764 enddo 765 enddo 766 767 ! Cubic system is assumed: phonon frequencies and permitivity at Gamma 768 ! do not depend on the q-vector direction 769 iqdir = 1 770 ABI_MALLOC(phfreq, (3*self%natom)) 771 ABI_MALLOC(investar, (3*self%natom)) 772 phfreq(:) = self%phfreq_qdir(:, iqdir) 773 investar(:) = self%investar(:, iqdir) 774 775 ! Effective mass sign 776 sigma = 1 777 if (self%band_params(1) < 0) sigma = -1 778 779 ! 3x3 Identity matrix 780 id33(:,:) = zero 781 do ii=1,3 782 id33(ii, ii) = one 783 enddo 784 785 ! FD: main loop 786 ! nqidr = 2*ntheta**2 (ntheta = nxi) 787 nxi = sqrt(one*self%nqdir/2) 788 ABI_MALLOC(intsum, (self%ndeg, self%ndeg, nkgrid, nkdir)) 789 ABI_MALLOC(zpr, (nkgrid, self%ndeg, nkdir)) 790 ABI_MALLOC(self%zpr_k, (nkgrid, self%ndeg, nkdir)) 791 ABI_MALLOC(lk_h, (self%ndeg, self%ndeg)) 792 ABI_MALLOC(intsuminv, (self%ndeg, self%ndeg)) 793 794 ! Values of ZPR (Eq. (86) of [Guster2021]) to be used in FD 795 zpr(:,:,:) = zero 796 do nu=4,3*self%natom 797 798 ! Summation over the infrared-active phonon modes 799 if (self%isiractive(nu)) then 800 801 do ikdir=1,nkdir 802 do iband=1,self%ndeg 803 do ik=1,nkgrid 804 kpt(:) = (ik - one)*deltak*unit_kdir(:, ikdir) 805 806 do ixi=0,nxi 807 xi = ixi*pi / (two*nxi) 808 if (ixi == nxi) xi = xi - tol8 809 810 ! q-vector length = (omega/A)^{1/2}*tan(xi) 811 qlen = sqrt(phfreq(nu) / abs(self%band_params(1)))*tan(xi) 812 do iqdir=1,self%nqdir 813 k_plus_q(:) = kpt(:) + qlen*self%unit_qdir(:, iqdir) 814 call lk_hamiltonian(self%band_params(:), k_plus_q(:), & 815 lk_h(:,:)) 816 817 intsum(:,:, ik, ikdir) = & 818 abs(lutt_eigenval(iband, ik, ikdir)*id33(:,:)) - (sigma * & 819 lk_h(:,:) + phfreq(nu)*id33(:,:)) 820 821 call mat3inv(intsum(:,:, ik, ikdir), intsuminv(:,:)) 822 823 tmp = dot_product(lutt_eigenvec(:, iband, 2, ikdir), & 824 matmul(intsuminv(:,:), lutt_eigenvec(:, iband, 2, ikdir))) 825 826 zpr(ik, iband, ikdir) = zpr(ik, iband, ikdir) + & 827 investar(nu)*phfreq(nu)*tmp*sqrt(phfreq(nu) / & 828 abs(self%band_params(1))) / cos(xi)**2 * & 829 self%weights_qdir(iqdir) 830 enddo 831 enddo 832 enddo 833 enddo 834 enddo 835 836 endif 837 enddo 838 zpr(:,:,:) = half*quarter*piinv*sigma/nxi * zpr(:,:,:) 839 self%zpr_k(:,:,:) = zpr(:,:,:) 840 841 ! FD: actual finite differences to obtain the second derivative of ZPR 842 ! (Eq. (80) of [Guster2021]) 843 ABI_MALLOC(zpr_ddk, (self%ndeg, nkdir)) 844 845 zpr_ddk(:,:) = zero 846 do ikdir=1,nkdir 847 do iband=1,self%ndeg 848 ! Copied from the previous implementation 849 ! Does not look like a usual FD expression... 850 zpr_ddk(iband, ikdir) = four/three *(zpr(2, iband, ikdir) - & 851 zpr(1, iband, ikdir) - (zpr(3, iband, ikdir) - & 852 zpr(1, iband, ikdir)) / 16.0_dp) * two/deltak**2 853 enddo 854 enddo 855 856 do ikdir=1,nkdir 857 do iband=1,self%ndeg 858 self%invpolmas(iband, ikdir) = & 859 two/deltak**2 * lutt_eigenval(iband, 2, ikdir) + zpr_ddk(iband, ikdir) 860 enddo 861 enddo 862 863 ABI_FREE(eigenvec) 864 ABI_FREE(eigenval) 865 ABI_FREE(work) 866 ABI_FREE(invefmas) 867 ABI_FREE(lutt_eigenvec) 868 ABI_FREE(lutt_eigenval) 869 ABI_FREE(phfreq) 870 ABI_FREE(investar) 871 ABI_FREE(intsum) 872 ABI_FREE(zpr) 873 ABI_FREE(lk_h) 874 ABI_FREE(zpr_ddk) 875 ABI_FREE(intsuminv) 876 877 endif 878 879 end subroutine frohlich_calc_polaronmass
m_frohlich/frohlich_calc_zpr [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlich_calc_zpr
FUNCTION
Description
INPUTS
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
OUTPUT
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
SOURCE
898 subroutine frohlich_calc_zpr(self) 899 900 !Arguments ------------------------------------ 901 !scalars 902 class(frohlich_t), intent(inout) :: self 903 !arrays 904 905 !Local variables------------------------------- 906 !scalars 907 integer :: iqdir, nu 908 integer :: iband, jband 909 integer :: lwork, info 910 logical :: sign_warn 911 !arrays 912 real(dp), allocatable :: efmas_qdir(:,:) 913 real(dp), allocatable :: invefmas_avg(:) 914 real(dp), allocatable :: eigenval(:), rwork(:) 915 complex(dpc), allocatable :: eigenvec(:,:), work(:) 916 complex(dpc), allocatable :: f3d(:,:) 917 logical, allocatable :: efmas_pos(:) 918 919 ! ************************************************************************* 920 921 ! TODO: check if the electronic and phonon parts are initialized 922 923 ! Initialization of the diagonalization routine for degenerate case 924 ABI_MALLOC(eigenval, (self%ndeg)) 925 if (self%ndeg > 1) then 926 ABI_MALLOC(eigenvec, (self%ndeg, self%ndeg)) 927 lwork = -1 928 ABI_MALLOC(rwork, (3*self%ndeg - 2)) 929 ABI_MALLOC(work, (1)) 930 call zheev('V', 'U', self%ndeg, eigenvec(:,:), self%ndeg, eigenval(:), & 931 work(:), lwork, rwork(:), info) 932 lwork = int(work(1)) 933 ABI_FREE(work) 934 ABI_MALLOC(work, (lwork)) 935 endif 936 937 ABI_MALLOC(f3d, (self%ndeg, self%ndeg)) 938 ABI_MALLOC(invefmas_avg, (self%ndeg)) 939 ABI_MALLOC(efmas_pos, (self%ndeg)) 940 ABI_MALLOC(efmas_qdir, (self%ndeg, self%nqdir)) 941 ABI_MALLOC(self%sqrt_efmas_avg, (self%ndeg)) 942 ABI_MALLOC(self%saddle_warn, (self%ndeg)) 943 ABI_MALLOC(self%zpr_band, (self%ndeg)) 944 945 ! Effective mass tensor in q-vector directions 946 do iqdir=1,self%nqdir 947 948 ! Band curvature tensor (inverse effective mass) 949 do iband=1,self%ndeg 950 do jband=1,self%ndeg 951 f3d(iband, jband) = dot_product(self%unit_qdir(:, iqdir), matmul( & 952 self%eig2_diag_cart(:,:, iband, jband), self%unit_qdir(:, iqdir))) 953 enddo 954 enddo 955 956 ! Band curvature tensor diagonalization 957 if (self%ndeg == 1) then 958 eigenval(1) = f3d(1, 1) 959 else 960 eigenvec(:,:) = f3d(:,:) 961 eigenval(:) = zero 962 work(:) = zero 963 rwork(:) = zero 964 call zheev('V', 'U', self%ndeg, eigenvec(:,:), self%ndeg, eigenval(:), & 965 work(:), lwork, rwork(:), info) 966 ABI_CHECK(info == 0, sjoin("zheev returned info: ", itoa(info))) 967 endif 968 969 efmas_qdir(:, iqdir) = one/eigenval(:) 970 enddo 971 972 ! Integration over the q-sphere: square root of the effective mass average and 973 ! ZPR for each band (Eq. (17) of [deMelo2023]) 974 self%saddle_warn(:) = .false. 975 self%sqrt_efmas_avg(:) = zero 976 self%zpr_band(:) = zero 977 invefmas_avg(:) = zero 978 do iqdir=1,self%nqdir 979 ! square root of the effective mass average 980 self%sqrt_efmas_avg(:) = self%sqrt_efmas_avg(:) + & 981 self%weights_qdir(iqdir)*sqrt(abs(efmas_qdir(:, iqdir))) 982 983 ! ZPR 984 do nu=4,3*self%natom 985 self%zpr_band(:) = self%zpr_band(:) + (self%weights_qdir(iqdir) * & 986 self%investar(nu, iqdir)*sqrt(self%phfreq_qdir(nu, iqdir) * & 987 abs(efmas_qdir(:, iqdir)))) 988 enddo 989 990 ! Inverse effective mass average: used to obtain the sign of the ZPR 991 invefmas_avg(:) = invefmas_avg(:) + & 992 self%weights_qdir(iqdir)/efmas_qdir(:, iqdir) 993 994 ! Check for saddle-points 995 if (iqdir == 1) efmas_pos(:) = (efmas_qdir(:, iqdir) > 0) 996 do iband=1,self%ndeg 997 if (efmas_pos(iband) .neqv. (efmas_qdir(iband, iqdir) > 0)) then 998 self%saddle_warn(iband) = .true. 999 endif 1000 enddo 1001 enddo 1002 self%sqrt_efmas_avg(:) = quarter*piinv*self%sqrt_efmas_avg(:) 1003 self%zpr_band(:) = sqrthalf*quarter*piinv*self%zpr_band(:) 1004 invefmas_avg(:) = quarter*piinv*invefmas_avg(:) 1005 1006 ! Check for the sign problem and caclulate total ZPR and polaron formation 1007 ! energy in the weak-coupling regime 1008 self%sign_warn = .false. 1009 do iband=1,self%ndeg 1010 if (self%saddle_warn(iband) .or. & 1011 (efmas_pos(iband) .neqv. efmas_pos(1))) then 1012 self%sign_warn = .true. 1013 else 1014 self%sqrt_efmas_avg(iband) = & 1015 sign(self%sqrt_efmas_avg(iband), invefmas_avg(iband)) 1016 1017 self%zpr_band(iband) = sign(self%zpr_band(iband), invefmas_avg(iband)) 1018 endif 1019 enddo 1020 1021 if (.not. sign_warn) then 1022 self%zpr = -sum(self%zpr_band(:)) / self%ndeg 1023 self%sqrt_efmas_tot = sum(self%sqrt_efmas_avg(:)) / self%ndeg 1024 self%enpol_wc = -self%sqrt_efmas_tot*sum(self%dielavg(:)) 1025 endif 1026 1027 ABI_FREE(eigenval) 1028 ABI_FREE(f3d) 1029 ABI_FREE(invefmas_avg) 1030 ABI_FREE(efmas_pos) 1031 ABI_FREE(efmas_qdir) 1032 if (self%ndeg > 1) then 1033 ABI_FREE(eigenvec) 1034 ABI_FREE(rwork) 1035 ABI_FREE(work) 1036 endif 1037 1038 end subroutine frohlich_calc_zpr
m_frohlich/frohlich_free_el [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlich_free_el
FUNCTION
Deallocate dynamic memory related to the electronic subspace and other related quantities
INPUTS
SOURCE
573 subroutine frohlich_free_el(self) 574 575 !Arguments ------------------------------------ 576 !scalars 577 class(frohlich_t),intent(inout) :: self 578 579 !Local variables------------------------------- 580 581 ! ************************************************************************* 582 583 self%isinitel = .false. 584 585 ! real 586 ABI_SFREE(self%sqrt_efmas_avg) 587 ABI_SFREE(self%zpr_band) 588 ABI_SFREE(self%invpolmas) 589 ABI_SFREE(self%zpr_k) 590 ABI_SFREE(self%invefmas) 591 592 ! complex 593 ABI_SFREE(self%eig2_diag_cart) 594 595 ! logical 596 ABI_SFREE(self%saddle_warn) 597 598 end subroutine frohlich_free_el
m_frohlich/frohlich_free_ph [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlich_free_ph
FUNCTION
Deallocate dynamic memory related to the vibrational subspace
INPUTS
SOURCE
613 subroutine frohlich_free_ph(self) 614 615 !Arguments ------------------------------------ 616 !scalars 617 class(frohlich_t),intent(inout) :: self 618 619 !Local variables------------------------------- 620 621 ! ************************************************************************* 622 623 self%isinitph = .false. 624 625 ! real 626 ABI_SFREE(self%unit_qdir) 627 ABI_SFREE(self%weights_qdir) 628 ABI_SFREE(self%dielt_qdir) 629 ABI_SFREE(self%phfreq_qdir) 630 ABI_SFREE(self%polarity_qdir) 631 ABI_SFREE(self%proj_polarity_qdir) 632 ABI_SFREE(self%investar) 633 ABI_SFREE(self%dielavg) 634 635 ! logical 636 ABI_SFREE(self%isiractive) 637 638 end subroutine frohlich_free_ph
m_frohlich/frohlich_init_el [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlich_init_el
FUNCTION
Description
INPUTS
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model cryst<crystal_t> = Structure defining the unit cell %rprimd(3, 3) = real space primitvie vectors kpt(3) = k-point characterizing the electronic subspace ndeg = Number of degenerate bands eig2_diag(3, 3, ndeg, ndeg) = Band curvature double tensor
OUTPUT
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
SOURCE
1062 subroutine frohlich_init_el(self, cryst, kpt, ndeg, eig2_diag) 1063 1064 !Arguments ------------------------------------ 1065 !scalars 1066 class(frohlich_t), intent(inout) :: self 1067 type(crystal_t), intent(in) :: cryst 1068 integer :: ndeg 1069 !arrays 1070 real(dp), intent(in) :: kpt(3) 1071 complex(dpc), intent(in) :: eig2_diag(3, 3, ndeg, ndeg) 1072 1073 !Local variables------------------------------- 1074 !scalars 1075 integer :: lwork, info 1076 integer :: iband, jband 1077 integer :: idir, ipar 1078 !arrays 1079 real(dp) :: unit_kdir(3, 3) 1080 real(dp) :: eigenval(ndeg), lutt_eigenval(ndeg, ndeg) 1081 real(dp), allocatable :: rwork(:) 1082 complex(dpc) :: eigenvec(ndeg, ndeg), lutt_eigenvec(ndeg, ndeg) 1083 complex(dpc), allocatable :: work(:) 1084 logical :: lutt_found(3) 1085 1086 ! ************************************************************************* 1087 1088 self%kpt(:) = kpt(:) 1089 self%ndeg = ndeg 1090 1091 ! Initialize the band curvatutre double tensor in Cartesian coordiantes 1092 ABI_MALLOC(self%eig2_diag_cart, (3, 3, ndeg, ndeg)) 1093 1094 do iband=1,ndeg 1095 do jband=1,ndeg 1096 self%eig2_diag_cart(:,:, iband, jband) = one/two_pi**2 * & 1097 matmul(matmul(cryst%rprimd(:,:),eig2_diag(:,:, iband, jband)), & 1098 transpose(cryst%rprimd(:,:))) 1099 enddo 1100 enddo 1101 1102 ! Determine the type of the Fr\"ohlich model 1103 ! TODO: check if a material actually has cubic symmetry for kind = 3? 1104 if (ndeg == 1) then 1105 ! Standard or Anisotropic Fr\"ohlich model 1106 self%kind = 1 1107 1108 ! Assuming the single band is parabolic in the three Cartesian directions 1109 unit_kdir(:, 1) = (/1, 0, 0/) 1110 unit_kdir(:, 2) = (/0, 1, 0/) 1111 unit_kdir(:, 3) = (/0, 0, 1/) 1112 1113 do idir=1,3 1114 self%band_params(idir) = dot_product(unit_kdir(:, idir), & 1115 matmul(self%eig2_diag_cart(:,:, 1, 1), unit_kdir(:, idir))) 1116 enddo 1117 1118 else if (ndeg == 3) then 1119 ! Cubic generalized Fr\"ohlich model with triply degenerate bands 1120 self%kind = 2 1121 1122 ! Symmety inequivalent cubic directions to obtain the Luttinger parameters 1123 unit_kdir(:, 1) = (/1, 0, 0/) 1124 unit_kdir(:, 2) = (/1, 1, 0/)/sqrt(2.0) 1125 unit_kdir(:, 3) = (/1, 1, 1/)/sqrt(3.0) 1126 1127 ! Initialize the diagonalization routine 1128 lwork = -1 1129 ABI_MALLOC(rwork, (3*ndeg - 2)) 1130 ABI_MALLOC(work, (1)) 1131 call zheev('V', 'U', ndeg, eigenvec(:,:), ndeg, eigenval, work(:), lwork, & 1132 rwork(:), info) 1133 lwork=int(work(1)) 1134 ABI_FREE(work) 1135 ABI_MALLOC(work, (lwork)) 1136 1137 lutt_eigenval(:, :) = zero 1138 ! Inverse effective mass tensor in the symmetry inequivalent directions 1139 do idir=1,3 1140 1141 do iband=1,ndeg 1142 do jband=1,ndeg 1143 lutt_eigenvec(iband, jband) = dot_product(unit_kdir(:, idir), & 1144 matmul(self%eig2_diag_cart(:, :, iband, jband), unit_kdir(:, idir))) 1145 enddo 1146 enddo 1147 1148 work(:) = zero 1149 rwork(:) = zero 1150 call zheev('V', 'U', ndeg, lutt_eigenvec(:,:), ndeg, & 1151 lutt_eigenval(idir, :), work(:), lwork, rwork(:), info) 1152 ABI_CHECK(info == 0, sjoin("zheev returned info: ", itoa(info))) 1153 enddo 1154 1155 ABI_FREE(work) 1156 ABI_FREE(rwork) 1157 1158 ! TODO: this looks ugly 1159 ! is there any better wa to analyze and set the Luttinger parameters? 1160 1161 ! Check degeneracies in the (100) direction, get A and B luttinger params 1162 ! Inverse effective masses are: 2*A and 2*B (twofold) 1163 if (abs(lutt_eigenval(1, 2) - lutt_eigenval(1, 3)) < tol5) then 1164 self%band_params(1)=half*lutt_eigenval(1, 1) 1165 self%band_params(2)=half*half*(lutt_eigenval(1, 2) + lutt_eigenval(1, 3)) 1166 else if (abs(lutt_eigenval(1, 2) - lutt_eigenval(1, 1)) < tol5) then 1167 self%band_params(1) = half*half*(lutt_eigenval(1, 2) + lutt_eigenval(1, 1)) 1168 self%band_params(2) = half*lutt_eigenval(1, 3) 1169 else 1170 self%lutt_warn(1) = .true. 1171 endif 1172 1173 ! Check degeneracies in the (111) direction, get C luttinger parameter 1174 ! Inverse effective masses are 2/3*(A + 2B - C) (twofold) and 1175 ! 2/3*(A + 2B + 2C) 1176 if (abs(lutt_eigenval(3, 2) - lutt_eigenval(3, 3)) < tol5) then 1177 self%band_params(3) = self%band_params(1) + 2*self%band_params(2) - & 1178 onehalf*half*(lutt_eigenval(3, 2) + lutt_eigenval(3, 3)) 1179 else if (abs(lutt_eigenval(3, 2) - lutt_eigenval(3, 1)) < tol5) then 1180 self%band_params(3) = self%band_params(1) + 2*self%band_params(2) - & 1181 onehalf*half*(lutt_eigenval(3, 2) + lutt_eigenval(3, 1)) 1182 else 1183 self%lutt_warn(2) = .true. 1184 endif 1185 1186 ! Verifty that the (110) direction inverse effective masses are coherent 1187 ! with the Luttinger parameters: 2*B, A + B + C, A + B - C 1188 lutt_found(:) = .false. 1189 do ipar=1,ndeg 1190 if (abs(lutt_eigenval(2, ipar) - 2*self%band_params(2)) < tol4) then 1191 lutt_found(1) = .true. 1192 else if (abs(lutt_eigenval(2, ipar) - (self%band_params(1) + & 1193 self%band_params(2) - self%band_params(3))) < tol4) then 1194 lutt_found(2) = .true. 1195 else if (abs(lutt_eigenval(2, ipar) - (self%band_params(1) + & 1196 self%band_params(2) + self%band_params(3))) < tol4) then 1197 lutt_found(3) = .true. 1198 endif 1199 enddo 1200 1201 if (.not. (all(lutt_found))) self%lutt_warn(3) = .true. 1202 1203 else 1204 ! Generalized Fr\"ohlich model 1205 self%kind = 0 1206 endif 1207 1208 ! The electronic subspace has been initialized 1209 self%isinitel = .true. 1210 1211 end subroutine frohlich_init_el
m_frohlich/frohlich_init_ph [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlich_init_ph
FUNCTION
Description
INPUTS
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model cryst<crystal_t> = Structure defining the unit cell %natom = number of atoms in the unit cell %ucvol = real space unit cell volume efmas_ntheta = number of points required for spherical integration used to obtain the effective mass tensor [Laflamme2016] ifc<ifc_type> = Dynamical matrix and interatomic force constants
OUTPUT
self<frohlich_t> = Datatype gathering information on the Fr\"ohlich model
SOURCE
1236 subroutine frohlich_init_ph(self, cryst, efmas_ntheta, ifc) 1237 1238 !Arguments ------------------------------------ 1239 !scalars 1240 class(frohlich_t), intent(inout) :: self 1241 type(crystal_t), intent(in) :: cryst 1242 type(ifc_type), intent(in) :: ifc 1243 integer, intent(in) :: efmas_ntheta 1244 !arrays 1245 1246 !Local variables------------------------------- 1247 !scalars 1248 integer :: ntheta, nphi, nqdir 1249 integer :: iphi, itheta, iqdir, nu 1250 ! integer :: ikpt, ideg, iband, jband 1251 ! integer :: ndeg 1252 real(dp) :: weight 1253 real(dp) :: weight_phi, phi_radians 1254 real(dp) :: costheta, sintheta, cosphi, sinphi 1255 !arrays 1256 ! real(dp) :: kpt(3) 1257 real(dp), allocatable :: gq_points_theta(:), gq_weights_theta(:) 1258 real(dp), allocatable :: gq_points_cosphi(:), gq_points_sinphi(:) 1259 1260 ! ************************************************************************* 1261 1262 self%natom = cryst%natom 1263 self%ucvol = cryst%ucvol 1264 self%gmet(:,:) = cryst%gmet(:,:) 1265 1266 ! Initialization of integrals 1267 ! This part allocates and initializes arrays unit_qdir(3,nqdir) and 1268 ! weights_qdir(nqdir) used for spherical integration over q-points directions 1269 ntheta = efmas_ntheta 1270 nphi = 2*ntheta 1271 nqdir = nphi*ntheta 1272 1273 ABI_MALLOC(gq_points_theta, (ntheta)) 1274 ABI_MALLOC(gq_weights_theta, (ntheta)) 1275 ABI_MALLOC(gq_points_cosphi, (nphi)) 1276 ABI_MALLOC(gq_points_sinphi, (nphi)) 1277 1278 self%nqdir = nqdir 1279 ABI_MALLOC(self%unit_qdir, (3, nqdir)) 1280 ABI_MALLOC(self%weights_qdir, (nqdir)) 1281 1282 call cgqf(ntheta, 1, zero, zero, zero, pi, gq_points_theta, gq_weights_theta) 1283 weight_phi = two*pi/real(nphi, dp) 1284 1285 do iphi=1,nphi 1286 phi_radians = weight_phi * (iphi-1) 1287 gq_points_cosphi(iphi) = cos(phi_radians) 1288 gq_points_sinphi(iphi) = sin(phi_radians) 1289 enddo 1290 1291 nqdir = 0 1292 do itheta=1,ntheta 1293 costheta = cos(gq_points_theta(itheta)) 1294 sintheta = sin(gq_points_theta(itheta)) 1295 weight = gq_weights_theta(itheta)*weight_phi*sintheta 1296 1297 do iphi=1,nphi 1298 cosphi = gq_points_cosphi(iphi) 1299 sinphi = gq_points_sinphi(iphi) 1300 nqdir = nqdir + 1 1301 1302 self%unit_qdir(1, nqdir) = sintheta*cosphi 1303 self%unit_qdir(2, nqdir) = sintheta*sinphi 1304 self%unit_qdir(3, nqdir) = costheta 1305 self%weights_qdir(nqdir) = weight 1306 enddo 1307 enddo 1308 1309 ABI_FREE(gq_points_theta) 1310 ABI_FREE(gq_weights_theta) 1311 ABI_FREE(gq_points_cosphi) 1312 ABI_FREE(gq_points_sinphi) 1313 1314 ! Initialization of the generalized Fr\"ohlich parameters 1315 1316 ABI_MALLOC(self%dielt_qdir, (nqdir)) 1317 ABI_MALLOC(self%phfreq_qdir, (3*self%natom, nqdir)) 1318 ABI_MALLOC(self%polarity_qdir, (3, 3*self%natom, nqdir)) 1319 ABI_MALLOC(self%proj_polarity_qdir, (3*self%natom, nqdir)) 1320 ABI_MALLOC(self%investar, (3*self%natom, nqdir)) 1321 ABI_MALLOC(self%dielavg, (3*self%natom)) 1322 ABI_MALLOC(self%isiractive, (3*self%natom)) 1323 1324 ! Phonon frequencies and mode polarity vectors for each q-vector direction 1325 call ifc%calcnwrite_nana_terms(cryst, nqdir, self%unit_qdir(:,:), & 1326 phfrq2l=self%phfreq_qdir, polarity2l=self%polarity_qdir) 1327 1328 ! High-frequency dielectric constant for each q-vector direction 1329 do iqdir=1,nqdir 1330 self%dielt_qdir(iqdir) = dot_product(self%unit_qdir(:, iqdir), & 1331 matmul(ifc%dielt(:,:), self%unit_qdir(:, iqdir))) 1332 enddo 1333 1334 self%investar(:,:) = zero 1335 self%dielavg(:) = zero 1336 ! Prjections of mode polarities on unit q-vectors, 1337 ! inverse effective dielectric constants for each mode and q-vector directions 1338 ! and dielectric average over q-vectors (Eqs. (22), (26) of [deMelo2023]) 1339 do iqdir=1,nqdir 1340 do nu=4,3*self%natom 1341 self%proj_polarity_qdir(nu, iqdir) = & 1342 dot_product(self%unit_qdir(:, iqdir), self%polarity_qdir(:, nu, iqdir)) 1343 1344 self%investar(nu, iqdir) = (self%proj_polarity_qdir(nu, iqdir) / & 1345 self%dielt_qdir(iqdir) / self%phfreq_qdir(nu, iqdir))**2 * & 1346 four*pi/self%ucvol 1347 1348 self%dielavg(nu) = self%dielavg(nu) + (self%weights_qdir(iqdir) * & 1349 self%investar(nu, iqdir)*sqrt(self%phfreq_qdir(nu, iqdir))) 1350 enddo 1351 enddo 1352 self%dielavg(:) = sqrthalf*quarter*piinv * self%dielavg(:) 1353 1354 self%isiractive(:) = .false. 1355 ! Check for the infrared-active phonon modes: dielectric average is non-zero 1356 do nu=4,3*self%natom 1357 if (self%dielavg(nu) > tol10) self%isiractive(nu) = .true. 1358 enddo 1359 1360 ! Rough correction for the ZPR in the infrared limit of the Fr\"ohlich model 1361 self%zpr_gamma = zero 1362 do nu=4,3*self%natom 1363 do iqdir=1,self%nqdir 1364 self%zpr_gamma = self%zpr_gamma + & 1365 self%investar(nu, iqdir)*self%weights_qdir(iqdir) 1366 enddo 1367 enddo 1368 self%zpr_gamma = quarter*piinv*self%zpr_gamma 1369 self%zpr_gamma = two*(three*quarter*piinv/self%ucvol)**third * self%zpr_gamma 1370 1371 ! The phonon subspace has been initialized 1372 self%isinitph = .true. 1373 1374 end subroutine frohlich_init_ph
m_frohlich/frohlich_t [ Types ]
[ Top ] [ m_frohlich ] [ Types ]
NAME
frohlich_t
FUNCTION
Description
SOURCE
52 type,public :: frohlich_t 53 54 integer :: kind = 0 55 ! Type of the Fr\"ohlich model 56 ! used to access the calculations of various properties 57 ! 0 -> ndeg != 1 && ndeg != 3, generalized Fr\"ohlich model 58 ! 1 -> ndeg = 1, standard (possibly anisotropic) Fr\"ohlich model 59 ! 2 -> ndeg = 3, cubic generalized Fr\"ohlich model with 3-fold degeneracy 60 61 ! Geometry ------------------------------------ 62 63 real(dp) :: ucvol 64 ! Real space unit cell volume 65 66 real(dp) :: gmet(3,3) 67 ! Reciprocal space metric 68 69 70 ! Electroncic subspace ----------------------- 71 72 logical :: isinitel = .false. 73 ! Flag indicating that the electronic subspace has been initialized 74 75 real(dp) :: kpt(3) 76 ! k-point characterizing the electronic subspace, i.e. the k-point at which 77 ! the effective mass tensor is obtained (usually, CBM or VBM) 78 79 integer :: ndeg 80 ! Number of degenerate bands taken into account 81 82 complex(dpc), allocatable :: eig2_diag_cart(:,:,:,:) 83 ! Band curvature double tensor in Cartesian coordinates 84 ! (3, 3, ndeg, ndeg) 85 86 real(dp) :: band_params(3) 87 ! Parameters describing electronic bands 88 ! The meaning of this arrayd depends on the value of the kind variable 89 ! kind = 0 -> undefined 90 ! kind = 1 -> inverse effetive masses along the 100, 010 and 001 directions 91 ! kind = 2 -> Luttinger-Kohn parameters A, B, C 92 93 logical :: lutt_warn(3) = .false. 94 95 real(dp), allocatable :: sqrt_efmas_avg(:) 96 ! Square root of effective mass averaged over q-sphere for each band 97 ! (ndeg) 98 99 logical, allocatable :: saddle_warn(:) 100 ! Signals if the k-point characterizing the electronic subspace is a 101 ! saddle-point for each of the degenerate bands 102 ! (ndeg) 103 104 real(dp) :: sqrt_efmas_tot 105 ! Total square root of effective mass average 106 107 real(dp), allocatable :: invefmas(:,:) 108 ! Inverse electronic effective masses along 109 ! kind = 1 -> 100, 010, 001 directions 110 ! kind = 2 -> 100, 110, 111 directions 111 ! (ndeg, 3) 112 113 114 ! Vibrational (phonon) subspace --------------- 115 116 logical :: isinitph = .false. 117 ! Flag indicating that the phonon subspace has been initialized 118 119 integer :: natom 120 ! Number of atoms in the cell 121 122 integer :: nqdir 123 ! Number of points for spherical integration used to compute ZPR and other 124 ! quantities 125 126 real(dp), allocatable :: unit_qdir(:,:) 127 ! Unit q-vectors representing reciprocal space directions used to compute 128 ! ZPR and other quantities 129 ! (3, nqdir) 130 131 real(dp), allocatable :: weights_qdir(:) 132 ! Gaussian quadrature weights used for spherical intragraion over q-vectors 133 ! (nqdir) 134 135 real(dp), allocatable :: dielt_qdir(:) 136 ! High-frequency dielectric constant for each q-vector direction 137 ! (nqdir) 138 139 real(dp), allocatable :: phfreq_qdir(:,:) 140 ! Phonon frequencies for each mode and q-vector direction 141 ! (3*natom, nqdir) 142 143 real(dp), allocatable :: polarity_qdir(:,:,:) 144 ! Mode polarity vectors for each mode and q-vector direction 145 ! (3, 3*natom, nqdir) 146 147 real(dp), allocatable :: proj_polarity_qdir(:,:) 148 ! Projections of mode polarity vectors for each mode and q-vector direction 149 ! (3*natom, nqdir) 150 151 real(dp), allocatable :: investar(:,:) 152 ! Inverse effective dielectric constant for each mode and q-vector 153 ! (3*natom, nqdir) 154 155 real(dp), allocatable :: dielavg(:) 156 ! Dielectric average over q-vectors for each mode (Eq. (26) of [deMelo2023]) 157 ! (3*natom) 158 159 logical, allocatable :: isiractive(:) 160 ! Flags to detect the infrared-active phonon modes 161 ! (3*natom) 162 163 real(dp) :: dielt_eff 164 ! Effective dielectric constant in the strong-coupling regime 165 166 real(dp) :: phfreq_eff 167 ! Effective LO phonon frequency in the strong-coupling regime 168 169 170 ! Weak-coupling parameters -------------------- 171 172 real(dp), allocatable :: zpr_band(:) 173 ! Zero-point renormalization energy for each band (Eq. (17) of [deMelo2023]) 174 ! (ndeg) 175 176 logical :: sign_warn = .false. 177 ! Sginals an error if a saddle-point is encountered or bands contribute to 178 ! the ZPR with different signs 179 180 real(dp) :: zpr_gamma 181 ! Correction to the ZPR taking into account the infrared divergence of the 182 ! electron-phonon coupling in the Fr\"ohlich model 183 184 real(dp) :: zpr 185 ! Total ZPR for these bands and k-point 186 187 real(dp) :: enpol_wc 188 ! Fr\"ohlich polaron formation energy in the weak-coupling regime 189 ! (Eq. (26) of [deMelo2023]) 190 191 real(dp), allocatable :: zpr_k(:,:,:) 192 ! Direction dependent ZPR (Eq. (86) of [Guster2021]) 193 ! (3, ndeg, 3) 194 195 real(dp), allocatable :: invpolmas(:,:) 196 ! Inverse Fr\"ohlich polaron effective masses in the weak-coupling regime 197 ! (Sec. III A, B of [Guster2021]) 198 ! kind = 1 -> 100, 010, 001 directions 199 ! kind = 2 -> 100, 110, 111 directions 200 ! (ndeg, 3) 201 202 203 contains 204 205 procedure :: init_ph => frohlich_init_ph 206 ! Initialization of a vibrational (phonon) subspace parameters 207 208 procedure :: init_el => frohlich_init_el 209 ! Initialization of an electronic subspace parameters 210 211 procedure :: free_ph => frohlich_free_ph 212 ! Free memory allocated to the vibrational subspace 213 214 procedure :: free_el => frohlich_free_el 215 ! Free memory allocated to the electronic subspace 216 ! and other related quantities 217 218 procedure :: calc_zpr => frohlich_calc_zpr 219 ! Calculate the zero-point renormalization energy corresponding to the 220 ! weak-coupling treatment of the Fr\"ohlich model 221 222 procedure :: calc_polaronmass => frohlich_calc_polaronmass 223 ! Calculate the polaron effective mass corresponding to the weak-coupling 224 ! treatment of the Fr\"ohlich model; available for kind = 1 or 2 225 226 ! procedure :: ncwrite => frohlich_ncwrite 227 ! Write main dimensions and header on a netcdf file 228 229 230 end type frohlich_t
m_frohlich/frohlichmodel_polaronmass [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlichmodel_polaronmass
FUNCTION
Main routine to compute the polaron effective masses of the generalized Fr\"ohlich model and other related quantities
INPUTS
cryst<crystal_t>=Structure defining the unit cell dtset<dataset_type>=All input variables for this dataset. efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band degeneracy at each k point efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= double tensor datastructure efmasval(:,:)%eig2_diag band curvature double tensor ifc<ifc_type>=contains the dynamical matrix and the IFCs.
NOTES
This routine has to be merged with the frohlichmodel_zpr routine, and their text output needs to be refined. For now, they are being kept to be compatible with the legacy unit tests.
SOURCE
264 subroutine frohlichmodel_polaronmass(frohlich, cryst, dtset, efmasdeg, & 265 efmasval, ifc) 266 267 !Arguments ------------------------------------ 268 !scalars 269 type(frohlich_t),intent(inout) :: frohlich 270 type(crystal_t),intent(in) :: cryst 271 type(dataset_type),intent(in) :: dtset 272 type(ifc_type),intent(in) :: ifc 273 !arrays 274 type(efmasdeg_type), intent(in) :: efmasdeg(:) 275 type(efmasval_type), intent(in) :: efmasval(:,:) 276 277 !Local variables------------------------------- 278 !scalar 279 integer :: nu, ikpt, ideg, ndeg 280 ! integer :: iband 281 integer :: iqdir 282 !arrays 283 real(dp) :: kpt(3) 284 285 ! ************************************************************************* 286 287 ! Initalize phonon and dielectric subspace 288 call frohlich%init_ph(cryst, dtset%efmas_ntheta, ifc) 289 290 ! For each k-point (with possible degeneracy), initialize an inistance of the 291 ! generalize Fr\"ohlich model and calculate related quantities 292 do ikpt=1,dtset%nkpt 293 kpt(:) = dtset%kptns(:, ikpt) 294 295 do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2) 296 ndeg = efmasdeg(ikpt)%degs_bounds(2, ideg) - & 297 efmasdeg(ikpt)%degs_bounds(1, ideg) + 1 298 299 call frohlich%init_el(cryst, kpt, ndeg, efmasval(ideg, ikpt)%eig2_diag) 300 call frohlich%calc_zpr() 301 call frohlich%calc_polaronmass() 302 303 ! Luttinger parameters if applicable 304 if (frohlich%ndeg == 3) then 305 306 if (.not. (any(frohlich%saddle_warn))) then 307 if (any(frohlich%lutt_warn)) then 308 309 write(ab_out, '(2a)') ch10, & 310 ' Luttinger parameters could not be determined:' 311 312 if (frohlich%lutt_warn(1)) then 313 write(ab_out, '(a)') ' Predicted degeneracies for & 314 °_dim = 3 are not met for (100) direction.' 315 endif 316 317 if (frohlich%lutt_warn(2)) then 318 write(ab_out, '(a)') ' Predicted degeneracies for & 319 °_dim = 3 are not met for (111) direction.' 320 endif 321 322 if (frohlich%lutt_warn(3)) then 323 write(ab_out, '(a)') ' Predicted degeneracies for & 324 °_dim = 3 are not met for (110) direction.' 325 endif 326 327 write(ab_out, '(a)') ch10 328 else 329 write(ab_out, '(a,3f14.6)') & 330 ' Luttinger parameters (A, B, C) (a.u.): ', & 331 frohlich%band_params(:) 332 endif 333 endif 334 335 ! Print inverse electronic effective masses in the output 336 write(ab_out, '(a)') repeat('-', 80) 337 write(ab_out, '(a)') ' Polaron properties from the generalized & 338 &Froehlich model' 339 write(ab_out, '(a)') repeat('-', 80) 340 write(ab_out, '(a)') ' Polar modes' 341 write(ab_out, '(a)') ' ## Frequency(meV) Epsilon*' 342 343 ! For a cubic material, dielectric tensor and phonon frequencies at Gamma 344 ! do not depend on the q-vector direction 345 iqdir = 1 346 do nu = 4,3*cryst%natom 347 if (frohlich%isiractive(nu)) then 348 write(ab_out,'(2x,i3,5x,f15.6,5x,f15.6)') nu, & 349 frohlich%phfreq_qdir(nu, iqdir)*Ha_eV*1000.0_dp, & 350 one/frohlich%investar(nu, iqdir) 351 endif 352 enddo 353 354 write(ab_out, '(a)') ' ' 355 write(ab_out, '(a,f10.2)') ' ZPR (meV): ', & 356 frohlich%zpr_k(1, 1, 1)*Ha_eV*1000.0_dp 357 write(ab_out, '(a)') ' ' 358 write(ab_out, '(a)') ' Electronic effective mass (a.u.) & 359 &along 3 directions' 360 write(ab_out, '(a, 3f15.6)')' Direction 100: ', & 361 one/frohlich%invefmas(:, 1) 362 write(ab_out, '(a, 3f15.6)')' Direction 110: ', & 363 one/frohlich%invefmas(:, 2) 364 write(ab_out, '(a, 3f15.6)')' Direction 111: ', & 365 one/frohlich%invefmas(:, 3) 366 367 ! Print inverse polaron effective masses in the output 368 write(ab_out, '(a)') ' ' 369 write(ab_out, '(a)') ' Polaron effective mass (a.u.) along 3 directions' 370 write(ab_out, '(a, 3f15.6)') ' Direction 100: ', & 371 one/frohlich%invpolmas(:, 1) 372 write(ab_out, '(a, 3f15.6)') ' Direction 110: ', & 373 one/frohlich%invpolmas(:, 2) 374 write(ab_out, '(a, 3f15.6)') ' Direction 111: ', & 375 one/frohlich%invpolmas(:, 3) 376 write(ab_out, '(a)')' ' 377 write(ab_out, '(a)')' Sum rule of inverse polaron masses check-up & 378 &(for convergence purposes):' 379 write(ab_out,'(a, 3f15.6)')' Direction 100: ', & 380 sum(frohlich%invpolmas(:, 1)) 381 write(ab_out,'(a, 3f15.6)')' Direction 110: ', & 382 sum(frohlich%invpolmas(:, 2)) 383 write(ab_out,'(a, 3f15.6)')' Direction 111: ', & 384 sum(frohlich%invpolmas(:, 3)) 385 endif 386 387 call frohlich%free_el() 388 enddo 389 enddo 390 391 call frohlich%free_ph() 392 393 end subroutine frohlichmodel_polaronmass
m_frohlich/frohlichmodel_zpr [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
frohlichmodel_zpr
FUNCTION
Main routine to compute the ZPR of the generalized Fr\"ohlich model and other related quantities
INPUTS
cryst<crystal_t>=Structure defining the unit cell dtset<dataset_type>=All input variables for this dataset. efmasdeg(nkpt_rbz) <type(efmasdeg_type)>= information about the band degeneracy at each k point efmasval(mband,nkpt_rbz) <type(efmasdeg_type)>= double tensor datastructure efmasval(:,:)%eig2_diag band curvature double tensor ifc<ifc_type>=contains the dynamical matrix and the IFCs.
SOURCE
416 subroutine frohlichmodel_zpr(frohlich, cryst, dtset, efmasdeg, efmasval, ifc) 417 418 !Arguments ------------------------------------ 419 !scalars 420 type(frohlich_t),intent(inout) :: frohlich 421 type(crystal_t),intent(in) :: cryst 422 type(dataset_type),intent(in) :: dtset 423 type(ifc_type),intent(in) :: ifc 424 !arrays 425 type(efmasdeg_type), intent(in) :: efmasdeg(:) 426 type(efmasval_type), intent(in) :: efmasval(:,:) 427 428 !Local variables------------------------------- 429 !scalar 430 integer :: nu, ikpt, ideg, ndeg, iband 431 !arrays 432 real(dp) :: kpt(3) 433 434 ! ************************************************************************* 435 436 ! Initalize phonon and dielectric subspace 437 call frohlich%init_ph(cryst, dtset%efmas_ntheta, ifc) 438 439 ! Dielectric average 440 write(ab_out, '(a)') repeat('-', 80) 441 write(ab_out, '(a)') ' Dielectric average (EQ. 25 Melo2022)' 442 write(ab_out, '(a)') repeat('-', 80) 443 write(ab_out, '(a)') & 444 ' Mode <1/epsilon*SQRT(w_LO/2)> Cumulative sum' 445 446 do nu=4,3*cryst%natom 447 write(ab_out, '(i5,f28.12,f28.12)') nu, & 448 frohlich%dielavg(nu), sum(frohlich%dielavg(1:nu)) 449 enddo 450 write(ab_out, '(a)') repeat('-', 80) 451 452 ! Infrared ZPR correction (does not depend on the electronic subspace) 453 write(ab_out, '(6a,f14.6,a,f14.6,a)') ch10, & 454 ' Rough correction to the ZPR, to take into account the missing q=0 piece & 455 &using Frohlich model:', ch10, & 456 ' (+ for occupied states, - for unoccupied states) * zpr_q0_fact / & 457 &(Nqpt_full_bz)**(1/3) ', ch10, & 458 ' where Nqpt_full_bz=number of q wavevectors in full BZ, and zpr_q0_fact=',& 459 frohlich%zpr_gamma, ' Ha=', frohlich%zpr_gamma*Ha_eV, ' eV' 460 461 ! For each k-point (with possible degeneracy), initialize an inistance of the 462 ! generalize Fr\"ohlich model and calculate ZPR and other related quantities 463 do ikpt=1,dtset%nkpt 464 kpt(:) = dtset%kptns(:, ikpt) 465 466 do ideg=efmasdeg(ikpt)%deg_range(1),efmasdeg(ikpt)%deg_range(2) 467 ndeg = efmasdeg(ikpt)%degs_bounds(2, ideg) - & 468 efmasdeg(ikpt)%degs_bounds(1, ideg) + 1 469 470 call frohlich%init_el(cryst, kpt, ndeg, efmasval(ideg, ikpt)%eig2_diag) 471 call frohlich%calc_zpr() 472 473 ! Print ZPR results for each k-point 474 if (ndeg == 1) then 475 write(ab_out, '(2a,3(f6.3,a),i5)') ch10, & 476 ' - At k-point (', kpt(1), ',', kpt(2), ',', kpt(3), '), band ', & 477 efmasdeg(ikpt)%degs_bounds(1,ideg) 478 else 479 write(ab_out, '(2a,3(f6.3,a),i5,a,i5)') ch10, & 480 ' - At k-point (', kpt(1), ',', kpt(2), ',', kpt(3), '), bands ', & 481 efmasdeg(ikpt)%degs_bounds(1,ideg), ' through ', & 482 efmasdeg(ikpt)%degs_bounds(2,ideg) 483 endif 484 485 ! Luttinger parameters if applicable 486 if (frohlich%kind == 2) then 487 488 if (.not. (any(frohlich%saddle_warn))) then 489 if (any(frohlich%lutt_warn)) then 490 491 write(ab_out, '(2a)') ch10, & 492 ' Luttinger parameters could not be determined:' 493 494 if (frohlich%lutt_warn(1)) then 495 write(ab_out, '(a)') ' Predicted degeneracies for & 496 °_dim = 3 are not met for (100) direction.' 497 endif 498 499 if (frohlich%lutt_warn(2)) then 500 write(ab_out, '(a)') ' Predicted degeneracies for & 501 °_dim = 3 are not met for (111) direction.' 502 endif 503 504 if (frohlich%lutt_warn(3)) then 505 write(ab_out, '(a)') ' Predicted degeneracies for & 506 °_dim = 3 are not met for (110) direction.' 507 endif 508 509 write(ab_out, '(a)') ch10 510 else 511 write(ab_out, '(a,3f14.6)') & 512 ' Luttinger parameters (A, B, C) [at. units]: ', & 513 frohlich%band_params(:) 514 endif 515 endif 516 endif 517 518 ! Effective mass average and ZPR 519 do iband=1,ndeg 520 if (frohlich%saddle_warn(iband)) then 521 write(ab_out, '(a,i5,a)') ' Band ', & 522 efmasdeg(ikpt)%degs_bounds(1, ideg) + iband - 1, ' SADDLE POINT - & 523 &Frohlich effective mass and ZPR cannot be defined. ' 524 else 525 write(ab_out, '(a,i5,a,f14.10)') ' Band ', & 526 efmasdeg(ikpt)%degs_bounds(1, ideg) + iband - 1, ' Angular average & 527 &effective mass for Frohlich model (<m**0.5>)**2= ', & 528 sign(frohlich%sqrt_efmas_avg(iband)**2, & 529 frohlich%sqrt_efmas_avg(iband)) 530 endif 531 enddo 532 533 if (.not. frohlich%sign_warn) then 534 write(ab_out, '(2a)') & 535 ' Angular and band average effective mass and ZPR for Frohlich model.' 536 537 write(ab_out, '(a,es16.6)') ' Value of (<<m**0.5>>)**2 = ', & 538 (sum(abs(frohlich%sqrt_efmas_avg(:))) / ndeg)**2 539 540 write(ab_out, '(a,es16.6)') ' Absolute Value of <<m**0.5>> = ', & 541 (sum(abs(frohlich%sqrt_efmas_avg(:))) / ndeg) 542 543 write(ab_out, '(a,es16.6,a,es16.6,a)') & 544 ' ZPR from Frohlich model = ', frohlich%zpr, ' Ha=', & 545 frohlich%zpr*Ha_eV,' eV' 546 else 547 write(ab_out, '(a)') ' Angular and band average effective mass for & 548 &Frohlich model cannot be defined because of a sign problem.' 549 endif 550 551 call frohlich%free_el() 552 enddo 553 enddo 554 555 call frohlich%free_ph() 556 557 end subroutine frohlichmodel_zpr
m_frohlich/lk_hamiltonian [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
lk_hamiltonian
FUNCTION
Construct a Luttinger-Kohn Hamiltonian at a given k-point for a set of the Luttinger-Kohn parameters A, B, C
INPUTS
band_params(3) = Luttinger-Kohn parameters A, B, C kpt(3) = wavevector at which the Hamiltonian is obtained (dimensional)
OUTPUT
h_lk(3, 3) = Luttingher-Kohn Hamiltonian, H(k)
SOURCE
1452 subroutine lk_hamiltonian(band_params, kpt, h_lk) 1453 1454 !Arguments ------------------------------------ 1455 !arrays 1456 real(dp), intent(in) :: band_params(3) 1457 real(dp), intent(in) :: kpt(3) 1458 real(dp), intent(out) :: h_lk(3, 3) 1459 1460 !Local variables------------------------------- 1461 1462 ! ************************************************************************* 1463 1464 h_lk(1, 1) = band_params(1)*kpt(1)**2 + band_params(2)*(kpt(2)**2 + kpt(3)**2) 1465 h_lk(2, 2) = band_params(1)*kpt(2)**2 + band_params(2)*(kpt(1)**2 + kpt(3)**2) 1466 h_lk(3, 3) = band_params(1)*kpt(3)**2 + band_params(2)*(kpt(1)**2 + kpt(2)**2) 1467 h_lk(1, 2) = band_params(3)*kpt(1)*kpt(2) 1468 h_lk(1, 3) = band_params(3)*kpt(1)*kpt(3) 1469 h_lk(2, 3) = band_params(3)*kpt(2)*kpt(3) 1470 ! Symmetric matrix 1471 h_lk(2, 1) = h_lk(1, 2) 1472 h_lk(3, 1) = h_lk(1, 3) 1473 h_lk(3, 2) = h_lk(2, 3) 1474 1475 end subroutine lk_hamiltonian
m_frohlich/mat3inv [ Functions ]
[ Top ] [ m_frohlich ] [ Functions ]
NAME
mat3inv
FUNCTION
Inverts a 3x3 matrix of real elements
INPUTS
mm(3, 3) = real matrix to be inverted
OUTPUT
mit(3, 3) = inverse of the input matrix
SOURCE
1393 subroutine mat3inv(mm, mit) 1394 1395 !Arguments ------------------------------------ 1396 !arrays 1397 real(dp), intent(in) :: mm(3, 3) 1398 real(dp), intent(out) :: mit(3, 3) 1399 1400 !Local variables------------------------------- 1401 !scalars 1402 real(dp) :: dd 1403 character(len=500) :: msg 1404 !arrays 1405 real(dp) :: tt(3, 3) 1406 1407 ! ************************************************************************* 1408 1409 ! Minors and determinant 1410 tt(1, 1) = mm(2, 2)*mm(3, 3) - mm(3, 2)*mm(2, 3) 1411 tt(1, 2) = mm(1, 3)*mm(3, 2) - mm(1, 2)*mm(3, 3) 1412 tt(1, 3) = mm(1, 2)*mm(2, 3) - mm(1, 3)*mm(2, 2) 1413 tt(2, 1) = mm(2, 3)*mm(3, 1) - mm(2, 1)*mm(3, 3) 1414 tt(2, 2) = mm(1, 1)*mm(3, 3) - mm(3, 1)*mm(1, 3) 1415 tt(2, 3) = mm(1, 3)*mm(2, 1) - mm(1, 1)*mm(2, 3) 1416 tt(3, 1) = mm(2, 1)*mm(3, 2) - mm(2, 2)*mm(3, 1) 1417 tt(3, 2) = mm(1, 2)*mm(3, 1) - mm(1, 1)*mm(3, 2) 1418 tt(3, 3) = mm(1, 1)*mm(2, 2) - mm(2, 1)*mm(1, 2) 1419 1420 dd = mm(1, 1)*tt(1, 1) + mm(2, 1)*tt(2, 1) + mm(3, 1)*tt(3, 1) 1421 1422 ! Make sure the matrix is not singular 1423 if (dd /=0) then 1424 mit(:,:) = tt(:,:) / dd 1425 else 1426 write(msg, '(2a,2x,9(i0,1x),a)') 'Attempting to invert real array',ch10, & 1427 mm,' ==> determinant is zero.' 1428 ABI_ERROR(msg) 1429 endif 1430 1431 end subroutine mat3inv