TABLE OF CONTENTS
- ABINIT/m_rta
- m_rta/compute_rta
- m_rta/compute_rta_mobility
- m_rta/ibte_calc_tensors
- m_rta/ibte_driver
- m_rta/print_rta_txt_files
- m_rta/rta_driver
- m_rta/rta_free
- m_rta/rta_ncwrite
- m_rta/rta_new
- m_rta/rta_t
- m_rta/write_tensor
ABINIT/m_rta [ Modules ]
NAME
m_rta
FUNCTION
This module provides objects and procedures to compute transport properties by solving the linearized Boltzmann transport equation (BTE) in the relaxation-time approximation (RTA). or within the Iterative Bolztmann Equation (IBTE). The RTA has two different flavors: Self-energy Relaxation Time Approximation (SERTA) in which the back-scattering term is completely ignored and the Momentum-Relaxation Time Approximation (MRTA) in which backscattering is partly accounter for by multiplying the e-ph self-energy integrand function by the efficiency factor alpha that depends on the incoming/outgoing electron group velocity. The implementation assumes e-ph scattering although additional scattering mechanims (e.g. ionized impurities) can be easily included once an appropriate model is added to the ab-initio e-ph scattering rates.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (HM, MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
Dimensional analysis for conductivity (sigma), mobility (mu). [sigma] = Siemens/m with S = Ampere/Volt = Ohm^-1 [mu] = S L^2 Q
SOURCE
30 #if defined HAVE_CONFIG_H 31 #include "config.h" 32 #endif 33 34 #include "abi_common.h" 35 36 module m_rta 37 38 use defs_basis 39 use m_abicore 40 use m_xmpi 41 use m_errors 42 use m_copy 43 use m_ebands 44 use m_nctk 45 use m_wfk 46 use m_ephtk 47 use m_sigmaph 48 use m_dtset 49 use m_dtfil 50 use m_krank 51 #ifdef HAVE_NETCDF 52 use netcdf 53 #endif 54 55 use defs_datatypes, only : pseudopotential_type, ebands_t 56 use m_io_tools, only : open_file 57 use m_time, only : cwtime, cwtime_report 58 use m_crystal, only : crystal_t 59 use m_numeric_tools, only : bisect, simpson_int, safe_div, arth 60 use m_fstrings, only : strcat, sjoin, itoa, ltoa, stoa, ftoa, yesno 61 use m_kpts, only : kpts_timrev_from_kptopt, kpts_map 62 use m_occ, only : occ_fd, occ_dfde 63 use m_pawtab, only : pawtab_type 64 use m_ddk, only : ddkstore_t 65 66 implicit none 67 68 private
m_rta/compute_rta [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
compute_rta
FUNCTION
INPUTS
cryst<crystal_t>=Crystalline structure dtset<dataset_type>=All input variables for this dataset. dtfil<datafiles_type>=variables related to files. comm=MPI communicator.
SOURCE
672 subroutine compute_rta(self, cryst, dtset, dtfil, comm) 673 674 !Arguments ------------------------------------ 675 integer,intent(in) :: comm 676 class(rta_t),intent(inout) :: self 677 type(dataset_type),intent(in) :: dtset 678 type(datafiles_type),intent(in) :: dtfil 679 type(crystal_t),intent(in) :: cryst 680 681 !Local variables ------------------------------ 682 integer,parameter :: nvecs0 = 0, master = 0 683 integer :: nsppol, nkibz, ib, ik_ibz, iw, spin, ii, jj, itemp, irta, itens, iscal, cnt 684 integer :: ntens, edos_intmeth, ifermi, iel, nvals, my_rank 685 #ifdef HAVE_NETCDF 686 integer :: ncid 687 #endif 688 !character(len=500) :: msg 689 character(len=fnlen) :: path 690 real(dp) :: emin, emax, edos_broad, edos_step, max_occ, kT, Tkelv, linewidth, fact0, cpu, wall, gflops 691 !arrays 692 integer :: unts(2) 693 real(dp) :: vr(3), dummy_vecs(1,1,1,1,1), work_33(3,3), S_33(3,3), mat33(3,3) 694 real(dp),allocatable :: vv_tens(:,:,:,:,:,:,:), out_valsdos(:,:,:,:), dummy_dosvecs(:,:,:,:,:) 695 real(dp),allocatable :: out_tensdos(:,:,:,:,:,:), tau_vals(:,:,:,:,:), l0inv_33nw(:,:,:) 696 697 !************************************************************************ 698 699 call cwtime(cpu, wall, gflops, "start") 700 my_rank = xmpi_comm_rank(comm) 701 unts = [std_out, ab_out] 702 703 ! Basic dimensions 704 nsppol = self%ebands%nsppol; nkibz = self%ebands%nkpt 705 706 ! Allocate v x v tensors with and without the lifetimes. Eq 8 of [[cite:Madsen2018]] 707 ! The total number of tensorial entries is ntens and accounts for nrta 708 ! Remember that we haven't computed all the k-points in the IBZ hence we can have zero linewidths 709 ! or very small values when the states are at the band edge so we use safe_dif to avoid SIGFPE. 710 ! Also, note how we store only the states in the energy window. 711 712 nvals = self%ntemp * self%nrta 713 ABI_CALLOC(tau_vals, (self%ntemp, self%nrta, self%bmin:self%bmax, nkibz, nsppol)) 714 715 ntens = (1 + self%ntemp) * self%nrta 716 ABI_CALLOC(vv_tens, (3, 3, 1 + self%ntemp, self%nrta, self%bmin:self%bmax, nkibz, nsppol)) 717 718 cnt = 0 719 do spin=1,nsppol 720 do ik_ibz=1,nkibz 721 !cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism. 722 do ib=self%bmin,self%bmax 723 724 vr(:) = self%vbks(:, ib, ik_ibz, spin) 725 ! Store outer product (v_bks x v_bks) in vv_tens. This part does not depend on T and irta. 726 do ii=1,3 727 do jj=1,3 728 vv_tens(ii, jj, 1, 1:self%nrta, ib, ik_ibz, spin) = vr(ii) * vr(jj) 729 end do 730 end do 731 732 ! Multiply by the lifetime (SERTA and MRTA) 733 do irta=1,self%nrta 734 do itemp=1,self%ntemp 735 linewidth = self%linewidths(itemp, ib, ik_ibz, spin, irta) 736 call safe_div(vv_tens(:,:, 1, irta, ib, ik_ibz, spin), two * linewidth, zero, & 737 vv_tens(:,:, 1 + itemp, irta, ib, ik_ibz, spin)) 738 call safe_div(one, two * linewidth, zero, tau_vals(itemp, irta, ib, ik_ibz, spin)) 739 end do 740 end do 741 742 end do 743 end do 744 end do 745 746 !call xmpi_sum(vv_tens, comm, ierr) 747 !call xmpi_sum(tau_vals, comm, ierr) 748 call cwtime_report(" compute_rta_loop1", cpu, wall, gflops) 749 750 ! Compute DOS and VV_DOS and VV_TAU_DOS 751 ! Define integration method and mesh step. 752 edos_intmeth = 2; if (dtset%prtdos /= 0) edos_intmeth = dtset%prtdos 753 edos_step = dtset%dosdeltae 754 if (edos_step == zero) edos_step = 0.001 755 !if (edos_step == zero) edos_step = ten / Ha_meV 756 edos_broad = dtset%tsmear 757 758 ! Set default energy range for DOS 759 ! If sigma_erange is set, get emin and emax from this variable 760 ! MG: TODO This value should be read from SIGEPH 761 ! Recheck metals 762 if (self%assume_gap) then 763 emin = huge(one); emax = -huge(one) 764 do spin=1,self%ebands%nsppol 765 if (dtset%sigma_erange(1) >= zero) emin = min(emin, self%gaps%vb_max(spin) + tol2 * eV_Ha - dtset%sigma_erange(1)) 766 if (dtset%sigma_erange(2) >= zero) emax = max(emax, self%gaps%cb_min(spin) - tol2 * eV_Ha + dtset%sigma_erange(2)) 767 end do 768 ABI_CHECK(emin /= huge(one), "Cannot initialize emin") 769 ABI_CHECK(emax /= -huge(one), "Cannot initialize emax") 770 else 771 emin = minval(self%eminmax_spin(1, :)); emin = emin - tol1 * abs(emin) 772 emax = maxval(self%eminmax_spin(2, :)); emax = emax + tol1 * abs(emax) 773 end if 774 775 ! Compute DOS, vv_dos and vvtau_DOS (v x v tau) 776 ! 777 ! out_valsdos: (nw, 2, nvals, nsppol) array with DOS for scalar quantities if nvals > 0 778 ! out_tensdos: (nw, 2, 3, 3, ntens, nsppol) array with DOS weighted by tensorial terms if ntens > 0 779 ! 780 ! Vectors and tensors are in Cartesian coordinates. 781 ! Note how we compute the DOS only between [emin, emax] to save time and memory 782 ! this implies that IDOS and edos%ifermi are ill-defined 783 784 self%edos = ebands_get_edos_matrix_elements(self%ebands, cryst, self%bsize, & 785 nvals, tau_vals, nvecs0, dummy_vecs, ntens, vv_tens, & 786 edos_intmeth, edos_step, edos_broad, & 787 out_valsdos, dummy_dosvecs, out_tensdos, comm, & 788 brange=[self%bmin, self%bmax], erange=[emin, emax]) 789 790 if (my_rank == master) then 791 call self%edos%print(unit=std_out, header="Computation of DOS, VV_DOS and VVTAU_DOS") 792 call self%edos%print(unit=ab_out, header="Computation of DOS, VV_DOS and VVTAU_DOS") 793 end if 794 795 call cwtime_report(" compute_rta_edos", cpu, wall, gflops) 796 797 ! Unpack data stored in out_tensdos with shape (nw, 2, 3, 3, ntens, nsppol) 798 self%nw = self%edos%nw 799 ABI_MALLOC(self%tau_dos, (self%nw, self%ntemp, nsppol, self%nrta)) 800 ! TODO: Exchange dims? 801 ABI_MALLOC(self%vv_dos, (self%nw, 3, 3, nsppol)) 802 ABI_MALLOC(self%vvtau_dos, (self%nw, 3, 3, self%ntemp, nsppol, self%nrta)) 803 804 do irta=1,self%nrta 805 do spin=1,nsppol 806 do itemp=1,self%ntemp+1 807 808 itens = itemp + (irta - 1) * (self%ntemp + 1) 809 if (itemp == 1) then 810 self%vv_dos(:,:,:,spin) = out_tensdos(:, 1, :, :, itens, spin) 811 else 812 self%vvtau_dos(:,:,:, itemp-1, spin, irta) = out_tensdos(:, 1, :, :, itens, spin) 813 end if 814 815 end do 816 end do 817 end do 818 819 ! Transfer data for tau(e) 820 do irta=1,self%nrta 821 do spin=1,nsppol 822 do itemp=1,self%ntemp 823 iscal = itemp + (irta - 1) * self%ntemp 824 self%tau_dos(:, itemp, spin, irta) = out_valsdos(:, 1, iscal, spin) 825 end do 826 end do 827 end do 828 829 ! Free memory 830 ABI_SFREE(out_tensdos) 831 ABI_SFREE(tau_vals) 832 ABI_SFREE(out_valsdos) 833 ABI_SFREE(dummy_dosvecs) 834 ABI_SFREE(vv_tens) 835 836 ! Compute Onsager coefficients. Eq 9 of [[cite:Madsen2018]] 837 ! See also Eqs 41, page 11 of https://arxiv.org/pdf/1402.6979.pdf 838 ! 839 ! L^\alpha(\mu, T) = \int de \sigma(e, T) (e - mu)^\alpha (-df/de) 840 ! 841 ! with \sigma(e, T) stored in vvtau_dos 842 843 ABI_MALLOC(self%l0, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 844 ABI_MALLOC(self%l1, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 845 ABI_MALLOC(self%l2, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 846 847 call onsager(0, self%l0) 848 call onsager(1, self%l1) 849 call onsager(2, self%l2) 850 851 call cwtime_report(" compute_rta_onsanger", cpu, wall, gflops) 852 853 ! Compute transport tensors, Eqs 12-15 of [[cite:Madsen2018]] and convert to SI units. 854 ABI_CALLOC(self%sigma, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 855 ABI_CALLOC(self%seebeck, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 856 ABI_CALLOC(self%kappa, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 857 ABI_CALLOC(self%pi, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 858 ABI_CALLOC(self%zte, (3, 3, self%nw, self%nsppol, self%ntemp, self%nrta)) 859 860 ! Sigma = L0 861 fact0 = (siemens_SI / Bohr_meter / cryst%ucvol) 862 self%sigma = fact0 * self%l0 863 864 ! Used to stored L0^-1 865 ABI_MALLOC(l0inv_33nw, (3, 3, self%nw)) 866 867 do irta=1,self%nrta 868 do spin=1,nsppol 869 do itemp=1,self%ntemp 870 871 TKelv = self%kTmesh(itemp) / kb_HaK; if (TKelv < one) Tkelv = one 872 873 ! S = -1/T L0^-1 L1 = -1/T sigma L1 874 do iw=1,self%nw 875 call inv33(self%l0(:, :, iw, spin, itemp, irta), work_33) 876 l0inv_33nw(:,:,iw) = work_33 877 self%seebeck(:,:,iw,spin,itemp,irta) = - (volt_SI / TKelv) * matmul(work_33, self%l1(:,:,iw,spin,itemp,irta)) 878 end do 879 880 ! kappa = 1/T [L2 - L1 L0^-1 L1] 881 ! HM: Check why do we need minus sign here to get consistent results with Boltztrap! 882 ! MG: Likely because of a different definition of kappa. 883 do iw=1,self%nw 884 work_33 = self%l1(:, :, iw, spin, itemp, irta) 885 work_33 = self%l2(:, :, iw, spin, itemp, irta) - matmul(work_33, matmul(l0inv_33nw(:, :, iw), work_33)) 886 !self%kappa(:,:,iw,spin, itemp,spin,irta) = - (volt_SI**2 * fact0 / TKelv) * work_33 887 self%kappa(:,:,iw,spin,itemp,irta) = + (volt_SI**2 * fact0 / TKelv) * work_33 888 end do 889 890 ! Peltier pi = -L1 L0^-1 891 do iw=1,self%nw 892 work_33 = self%l1(:, :, iw, spin, itemp, irta) 893 self%pi(:,:,iw,spin,itemp,irta) = - volt_SI * matmul(work_33, l0inv_33nw(:, :, iw)) 894 end do 895 896 ! ZT: S^T sigma S k^-1 T (tensor form with k=k_electronic only): 897 do iw=1,self%nw 898 S_33 = self%seebeck(:,:,iw,spin,itemp,irta) 899 S_33 = matmul(matmul(transpose(S_33), self%sigma(:,:,iw,spin,itemp,irta)), S_33) 900 call inv33(self%kappa(:,:,iw,spin,itemp,irta), work_33) 901 self%zte(:,:,iw,spin,itemp,irta) = matmul(S_33, work_33) * TKelv 902 end do 903 904 end do 905 end do 906 end do 907 908 ABI_FREE(l0inv_33nw) 909 910 ! Compute the index of the Fermi level and handle possible out of range condition. 911 ifermi = bisect(self%edos%mesh, self%ebands%fermie) 912 if (ifermi == 0 .or. ifermi == self%nw) then 913 ABI_ERROR("Bisection could not find the index of the Fermi level in edos%mesh!") 914 end if 915 916 max_occ = two / (self%nspinor * self%nsppol) 917 918 ! Conductivity 919 ABI_MALLOC(self%conductivity, (3, 3, self%ntemp, self%nsppol, self%nrta)) 920 do spin=1,self%nsppol 921 do itemp=1,self%ntemp 922 do irta=1,self%nrta 923 do jj=1,3 924 do ii=1,3 925 self%conductivity(ii,jj,itemp,spin,irta) = self%sigma(ii, jj, ifermi, spin, itemp, irta) * 0.01 !m^-1 to cm^-1 926 end do 927 end do 928 end do 929 end do ! itemp 930 end do ! spin 931 932 ABI_MALLOC(self%resistivity, (3, 3, self%ntemp, self%nrta)) 933 do irta=1,self%nrta 934 do itemp=1,self%ntemp 935 work_33 = sum(self%conductivity(:,:,itemp,:,irta), dim=3) 936 call inv33(work_33, mat33); mat33 = 1e+6_dp * mat33 937 self%resistivity(:, :, itemp, irta) = mat33 938 end do 939 end do 940 941 ! Mobility 942 ABI_MALLOC(self%n, (self%nw, self%ntemp, 2)) 943 ABI_MALLOC(self%mobility, (3, 3, self%nw, self%ntemp, 2, self%nsppol, self%nrta)) 944 945 do spin=1,self%nsppol 946 do itemp=1,self%ntemp 947 ! Compute carrier density 948 kT = self%kTmesh(itemp) 949 950 ! MG TODO: I think that here we should use mu_e instead of ifermi. 951 ! Compute carrier density of electrons (ifermi:self%nw) 952 do iw=1,self%nw ! doping 953 self%n(iw,itemp,1) = carriers(self%edos%mesh, self%edos%dos(:,spin) * max_occ, ifermi, self%nw, & 954 kT, self%edos%mesh(iw)) / cryst%ucvol / Bohr_meter**3 955 end do 956 957 ! Compute carrier density of holes (1:ifermi) 958 do iw=1,self%nw ! doping 959 self%n(iw,itemp,2) = carriers(self%edos%mesh, self%edos%dos(:,spin) * max_occ, 1, ifermi, & 960 kT, self%edos%mesh(iw)) / cryst%ucvol / Bohr_meter**3 961 end do 962 963 self%n(:,itemp,2) = self%n(self%nw,itemp,2) - self%n(:,itemp,2) 964 965 ! Compute mobility 966 do irta=1,self%nrta 967 do iel=1,2 968 do iw=1,self%nw 969 do jj=1,3 970 do ii=1,3 971 call safe_div(self%sigma(ii, jj, iw, spin, itemp, irta) * 100**2, & 972 e_Cb * self%n(iw, itemp, iel), & 973 zero, self%mobility(ii, jj, iw, itemp, iel, spin, irta)) 974 end do 975 end do 976 end do 977 end do 978 end do 979 end do ! itemp 980 end do ! spin 981 982 ! Compute RTA mobility 983 call self%compute_rta_mobility(cryst, comm) 984 985 if (my_rank == master) then 986 ! Print RTA results to stdout and other external txt files (for the test suite) 987 call self%print_rta_txt_files(cryst, dtset, dtfil) 988 989 ! Creates the netcdf file used to store the results of the calculation. 990 #ifdef HAVE_NETCDF 991 path = strcat(dtfil%filnam_ds(4), "_RTA.nc") 992 call wrtout(unts, ch10//sjoin("- Writing RTA transport results to:", path)) 993 NCF_CHECK(nctk_open_create(ncid, path , xmpi_comm_self)) 994 call self%rta_ncwrite(cryst, dtset, ncid) 995 NCF_CHECK(nf90_close(ncid)) 996 #endif 997 end if 998 999 call cwtime_report(" compute_rta", cpu, wall, gflops) 1000 1001 contains 1002 1003 real(dp) function carriers(wmesh, dos, istart, istop, kT, mu) 1004 1005 !Arguments ------------------------------------------- 1006 real(dp),intent(in) :: kT, mu 1007 real(dp),intent(in) :: wmesh(self%nw), dos(self%nw) 1008 integer,intent(in) :: istart, istop 1009 1010 !Local variables ------------------------------------- 1011 integer :: iw 1012 real(dp) :: kernel(self%nw), integral(self%nw) 1013 1014 kernel = zero 1015 do iw=istart,istop 1016 kernel(iw) = dos(iw) * occ_fd(wmesh(iw), kT, mu) 1017 end do 1018 call simpson_int(self%nw, edos_step, kernel, integral) 1019 carriers = integral(self%nw) 1020 1021 end function carriers 1022 1023 ! Compute L^\alpha(\mu, T) = \int de \sigma(e, T) (e - mu)^\alpha (-df/de) 1024 subroutine onsager(order, lorder) 1025 1026 !Arguments ------------------------------------------- 1027 integer,intent(in) :: order 1028 real(dp),intent(out) :: lorder(3, 3, self%nw, self%nsppol, self%ntemp, self%nrta) 1029 1030 !Local variables ------------------------------------- 1031 integer :: spin, iw, imu, irta 1032 real(dp) :: mu, ee, kT 1033 real(dp) :: kernel(self%nw,3,3,self%nsppol), integral(self%nw) 1034 1035 ! Get spin degeneracy 1036 max_occ = two / (self%nspinor * self%nsppol) 1037 1038 do irta=1,self%nrta 1039 do itemp=1,self%ntemp 1040 kT = self%kTmesh(itemp) 1041 ! Loop over chemical potentials mu 1042 do imu=1,self%nw 1043 mu = self%edos%mesh(imu) 1044 1045 ! Build integrand for given mu 1046 do iw=1,self%nw 1047 ee = self%edos%mesh(iw) 1048 if (order > 0) then 1049 kernel(iw,:,:,:) = - max_occ * self%vvtau_dos(iw,:,:,itemp,:,irta) * (ee - mu)** order * occ_dfde(ee, kT, mu) 1050 else 1051 kernel(iw,:,:,:) = - max_occ * self%vvtau_dos(iw,:,:,itemp,:,irta) * occ_dfde(ee, kT, mu) 1052 end if 1053 end do 1054 1055 ! Integrate with simpson_int 1056 do spin=1,self%nsppol 1057 do jj=1,3 1058 do ii=1,3 1059 call simpson_int(self%nw, edos_step, kernel(:,ii,jj, spin), integral) 1060 lorder(ii, jj, imu, spin, itemp, irta) = integral(self%nw) 1061 end do 1062 end do 1063 end do 1064 1065 end do ! imu 1066 end do ! itemp 1067 end do ! irta 1068 1069 end subroutine onsager 1070 1071 end subroutine compute_rta
m_rta/compute_rta_mobility [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
compute_rta_mobility
FUNCTION
INPUTS
cryst<crystal_t>=Crystalline structure comm=MPI communicator.
SOURCE
1088 subroutine compute_rta_mobility(self, cryst, comm) 1089 1090 !Arguments ------------------------------------ 1091 class(rta_t),intent(inout) :: self 1092 type(crystal_t),intent(in) :: cryst 1093 integer,intent(in) :: comm 1094 1095 !Local variables ------------------------------ 1096 integer :: nsppol, nkibz, ib, ik_ibz, spin, ii, jj, itemp, ieh, cnt, nprocs, irta, time_opt 1097 real(dp) :: eig_nk, mu_e, linewidth, fact, fact0, max_occ, kT, wtk, cpu, wall, gflops 1098 real(dp) :: vr(3), vv_tens(3,3), vv_tenslw(3,3) !, tmp_tens(3,3) 1099 1100 !************************************************************************ 1101 1102 call cwtime(cpu, wall, gflops, "start") 1103 1104 nprocs = xmpi_comm_size(comm) 1105 nkibz = self%ebands%nkpt; nsppol = self%ebands%nsppol 1106 1107 time_opt = 0 ! This to preserve the previous behaviour in which TR was not used. 1108 !time_opt = -1 ! This to preserve the previous behaviour in which TR was not used. 1109 1110 ABI_CALLOC(self%mobility_mu, (3, 3, 2, nsppol, self%ntemp, self%nrta)) 1111 ABI_CALLOC(self%conductivity_mu, (3, 3, 2, nsppol, self%ntemp, self%nrta)) 1112 1113 ! Compute (e/h) carriers per unit cell at the different temperatures. 1114 ABI_CALLOC(self%n_ehst, (2, self%nsppol, self%ntemp)) 1115 call ebands_get_carriers(self%ebands, self%ntemp, self%kTmesh, self%transport_mu_e, self%n_ehst) 1116 1117 ! Compute conductivity_mu i.e. results in which lifetimes have been computed in a consistent way 1118 ! with the same the Fermi level. In all the other cases, indeed, we assume that tau does not depend on ef. 1119 ! 1120 ! sigma_RTA = -S e^2 / (N_k Omega) sum_\nk (v_\nk \otimes v_\nk) \tau_\nk (df^0/de_\nk) 1121 ! 1122 ! with S the spin degeneracy factor. 1123 ! 1124 ! TODO: Implement other tensors. Compare these results with the ones obtained with spectral sigma 1125 ! In principle, they should be the same, in practice the integration of sigma requires enough resolution 1126 ! around the band edge. 1127 !print *, "in RTA max velocities", maxval(abs(self%vbks)) 1128 cnt = 0 1129 do spin=1,nsppol 1130 do ik_ibz=1,nkibz 1131 !cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism. 1132 wtk = self%ebands%wtk(ik_ibz) 1133 1134 do ib=self%bmin,self%bmax 1135 eig_nk = self%ebands%eig(ib, ik_ibz, spin) 1136 1137 ! Store outer product in vv_tens 1138 vr(:) = self%vbks(:, ib, ik_ibz, spin) 1139 ! Don't remove this if: it makes the loop a bit faster and, most importantly, 1140 ! it prevents intel from miscompiling the code. 1141 if (all(abs(vr) == zero)) cycle 1142 1143 do ii=1,3 1144 do jj=1,3 1145 vv_tens(ii, jj) = vr(ii) * vr(jj) 1146 end do 1147 end do 1148 1149 ! Symmetrize tensor. 1150 !print *, "intens", vv_tens 1151 vv_tens = cryst%symmetrize_cart_tens33(vv_tens, time_opt) 1152 !print *, "out_tens", vv_tens 1153 1154 ! Multiply by the lifetime (SERTA or MRTA) 1155 do irta=1,self%nrta 1156 do itemp=1,self%ntemp 1157 kT = self%kTmesh(itemp) 1158 mu_e = self%transport_mu_e(itemp) 1159 ieh = 2; if (eig_nk >= mu_e) ieh = 1 1160 linewidth = self%linewidths(itemp, ib, ik_ibz, spin, irta) 1161 !print *, linewidth, wtk, occ_dfde(eig_nk, kT, mu_e), "tens", vv_tens 1162 call safe_div( - wtk * vv_tens * occ_dfde(eig_nk, kT, mu_e), two * linewidth, zero, vv_tenslw) 1163 self%conductivity_mu(:, :, ieh, spin, itemp, irta) = self%conductivity_mu(:, :, ieh, spin, itemp, irta) & 1164 + vv_tenslw(:, :) 1165 end do 1166 end do 1167 end do 1168 1169 end do ! ik_ibz 1170 end do ! spin 1171 1172 !call xmpi_sum(self%conductivity_mu, comm, ierr) 1173 1174 ! Get units conversion factor including spin degeneracy. 1175 max_occ = two / (self%nspinor * self%nsppol) 1176 fact0 = max_occ * (siemens_SI / Bohr_meter / cryst%ucvol) / 100 1177 self%conductivity_mu = fact0 * self%conductivity_mu ! siemens cm^-1 1178 1179 ! Scale by the carrier concentration 1180 fact = 100**3 / e_Cb 1181 do irta=1,self%nrta 1182 do spin=1,nsppol 1183 do itemp=1,self%ntemp 1184 do ieh=1,2 ! e/h 1185 call safe_div(fact * self%conductivity_mu(:,:,ieh,spin,itemp, irta), & 1186 self%n_ehst(ieh, spin, itemp) / cryst%ucvol / Bohr_meter**3, zero, & 1187 self%mobility_mu(:,:,ieh,spin,itemp,irta)) 1188 end do 1189 end do 1190 end do 1191 end do 1192 1193 call cwtime_report(" compute_rta_mobility", cpu, wall, gflops) 1194 1195 end subroutine compute_rta_mobility
m_rta/ibte_calc_tensors [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
ibte_calc_tensors
FUNCTION
INPUTS
cryst<crystal_t>=Crystalline structure comm=MPI communicator.
SOURCE
2141 subroutine ibte_calc_tensors(self, cryst, itemp, kT, mu_e, fk, onsager, sigma_eh, mob_eh, fsum_eh, comm) 2142 2143 !Arguments ------------------------------------ 2144 class(rta_t),intent(inout) :: self 2145 type(crystal_t),intent(in) :: cryst 2146 integer,intent(in) :: itemp 2147 real(dp),intent(in) :: kT, mu_e 2148 real(dp),intent(in) :: fk(3, self%ebands%nkpt, self%bmin:self%bmax, self%nsppol) 2149 real(dp),intent(out) :: sigma_eh(3,3,2,self%nsppol), mob_eh(3,3,2,self%nsppol) 2150 real(dp),intent(out) :: fsum_eh(3,2,self%nsppol), onsager(3,3,3,self%nsppol) 2151 integer,intent(in) :: comm 2152 2153 !Local variables ------------------------------ 2154 !scalars 2155 integer :: nsppol, nkibz, ib, ik_ibz, spin, ii, jj, ieh, cnt, nprocs, ia, time_opt 2156 real(dp) :: eig_nk, fact, fact0, max_occ, wtk, emu_alpha 2157 !arrays 2158 real(dp) :: vr(3), vv_tens(3,3) 2159 2160 !************************************************************************ 2161 2162 ABI_UNUSED(kt) 2163 2164 nprocs = xmpi_comm_size(comm) 2165 2166 ! Copy important dimensions 2167 nkibz = self%ebands%nkpt; nsppol = self%ebands%nsppol 2168 time_opt = 0 ! This to preserve the previous behaviour in which TR was not used. 2169 2170 ! sigma_IBTE = (-S e^ / omega sum_\nk) (v_\nk \otimes F_\nk) 2171 ! with S the spin degeneracy factor. 2172 sigma_eh = zero; fsum_eh = zero; onsager = zero 2173 2174 ! Compute mobility_mu i.e. results in which lifetimes have been computed in a consistent way 2175 ! with the same the Fermi level. In all the other cases, indeed, we assume that tau does not depend on ef. 2176 ! 2177 ! TODO: Implement other tensors. Compare these results with the ones obtained with spectral sigma 2178 ! In principle, they should be the same, in practice the integration of sigma requires enough resolution 2179 ! around the band edge. 2180 cnt = 0 2181 do spin=1,nsppol 2182 do ik_ibz=1,nkibz 2183 !cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism. 2184 wtk = self%ebands%wtk(ik_ibz) 2185 2186 do ib=self%bmin,self%bmax 2187 eig_nk = self%ebands%eig(ib, ik_ibz, spin) 2188 2189 ! Compute outer product in vv_tens and symmetrize tensor. 2190 vr(:) = self%vbks(:, ib, ik_ibz, spin) 2191 do ia=1,3 2192 if (ia == 1) then 2193 emu_alpha = one 2194 else 2195 emu_alpha = (eig_nk - mu_e) ** (ia - 1) 2196 end if 2197 2198 do ii=1,3 2199 do jj=1,3 2200 vv_tens(ii, jj) = vr(ii) * fk(jj, ik_ibz, ib, spin) * emu_alpha 2201 end do 2202 end do 2203 vv_tens = cryst%symmetrize_cart_tens33(vv_tens, time_opt) 2204 2205 if (ia == 1) then 2206 ieh = 2; if (eig_nk >= mu_e) ieh = 1 2207 sigma_eh(:,:,ieh,spin) = sigma_eh(:,:,ieh,spin) - wtk * vv_tens 2208 fsum_eh(:,ieh,spin) = fsum_eh(:,ieh,spin) + wtk * cryst%symmetrize_cart_vec3(fk(:, ik_ibz, ib, spin), time_opt) 2209 end if 2210 onsager(:,:,ia,spin) = onsager(:,:,ia,spin) - wtk * vv_tens 2211 end do ! ia 2212 2213 end do ! ib 2214 end do ! ik_ibz 2215 end do ! spin 2216 2217 !call xmpi_sum(sigma_eh, comm, ierr) 2218 !call xmpi_sum(onsager, comm, ierr) 2219 ! Get units conversion factor including spin degeneracy. 2220 max_occ = two / (self%nspinor * self%nsppol) 2221 fact0 = max_occ * (siemens_SI / Bohr_meter / cryst%ucvol) / 100 2222 fact = 100**3 / e_Cb 2223 2224 sigma_eh = fact0 * sigma_eh ! siemens cm^-1 2225 fsum_eh = fsum_eh / cryst%ucvol 2226 2227 ! Scale by the carrier concentration. 2228 do spin=1,nsppol 2229 do ieh=1,2 2230 call safe_div(sigma_eh(:,:,ieh,spin) * fact, & 2231 self%n_ehst(ieh, spin, itemp) / cryst%ucvol / Bohr_meter**3, zero, mob_eh(:,:,ieh,spin)) 2232 end do 2233 end do 2234 2235 end subroutine ibte_calc_tensors
m_rta/ibte_driver [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
ibte_driver
FUNCTION
Driver to compute transport properties within the IBTE.
INPUTS
dtfil<datafiles_type>=variables related to files. ngfftc(18)=Coarse FFT meshe dtset<dataset_type>=All input variables for this dataset. ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...) cryst<crystal_t>=Crystalline structure pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. comm=MPI communicator.
SOURCE
1592 subroutine ibte_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm) 1593 1594 !Arguments ------------------------------------ 1595 !scalars 1596 integer, intent(in) :: comm 1597 type(datafiles_type),intent(in) :: dtfil 1598 type(dataset_type),intent(in) :: dtset 1599 type(crystal_t),intent(in) :: cryst 1600 type(ebands_t),intent(in) :: ebands 1601 type(pseudopotential_type),intent(in) :: psps 1602 !arrays 1603 integer,intent(in) :: ngfftc(18) 1604 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 1605 1606 !Local variables ------------------------------ 1607 integer,parameter :: master = 0 1608 integer :: spin, ikcalc, nkcalc, nbsum, nbcalc, itemp, iter, ierr, bsize 1609 integer :: nkibz, nsppol, band_k, ik_ibz, bmin, bmax, band_sum, ntemp, ii, jj, iq_sum, btype, nsp 1610 integer :: ikq_ibz, isym_kq, trev_kq, cnt, tag, nprocs, receiver, my_rank, isym, itime, isym_lgk 1611 #ifdef HAVE_NETCDF 1612 integer :: ncid, grp_ncid, ncerr 1613 #endif 1614 real(dp) :: kT, mu_e, e_nk, dfde_nk, tau_nk, lw_nk, max_adiff, cpu, wall, gflops, btype_fact, abs_tol, rtmp 1615 logical :: send_data 1616 character(len=500) :: msg 1617 character(len=fnlen) :: path 1618 type(rta_t) :: ibte 1619 !arrays 1620 integer :: unts(2), dims(4) 1621 logical,allocatable :: converged(:) 1622 real(dp) :: vec3(3), sym_vec(3), mat33(3,3), f_kq(3), work33(3,3) 1623 real(dp) :: fsum_eh(3,2,ebands%nsppol), max_adiff_spin(ebands%nsppol) 1624 real(dp) :: onsager(3,3,3,ebands%nsppol) 1625 real(dp),pointer :: sig_p(:,:,:,:), mob_p(:,:,:,:) 1626 real(dp),target,allocatable :: ibte_sigma(:,:,:,:,:), ibte_mob(:,:,:,:,:), ibte_rho(:,:,:) 1627 real(dp),allocatable :: grp_srate(:,:,:,:), fkn_in(:,:,:,:), fkn_out(:,:,:,:), fkn_serta(:,:,:,:), taukn_serta(:,:,:,:) 1628 character(len=2) :: components(3) 1629 1630 type :: scatk_t 1631 1632 integer :: rank = xmpi_undefined_rank 1633 1634 integer :: nq_ibzk_eff 1635 ! Number of effective q-points in the IBZ(k) 1636 1637 integer :: lgk_nsym 1638 ! Number of symmetry operations in the little group of k. 1639 1640 integer,allocatable :: lgk_sym2glob(:,:) 1641 ! lgk_sym2glob(2, lgk_nsym) 1642 ! Mapping isym_lg --> [isym, itime] 1643 ! where isym is the index of the operation in the global array **crystal%symrec** 1644 ! and itim is 2 if time-reversal T must be included else 1. Depends on ikcalc 1645 1646 integer,allocatable :: kq_symtab(:,:) 1647 ! kq_symtab(6, nq_ibzk_eff) 1648 1649 real(dp),allocatable :: vals(:,:,:,:) 1650 ! (nq_ibzk_eff, nbsum, nbcalc, ntemp) 1651 end type scatk_t 1652 1653 type(scatk_t),target,allocatable :: sr(:,:) 1654 type(scatk_t),pointer :: sr_p 1655 1656 ! ************************************************************************* 1657 1658 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 1659 unts = [std_out, ab_out] 1660 1661 call wrtout(unts, ch10//' Entering IBTE driver.') 1662 call wrtout(unts, sjoin("- Reading SERTA lifetimes and e-ph scattering operator from:", & 1663 dtfil%filsigephin), newlines=1, do_flush=.True.) 1664 1665 ! Initialize IBTE object 1666 ibte = rta_new(dtset, dtfil, ngfftc, cryst, ebands, pawtab, psps, comm) 1667 1668 ! Compute RTA transport quantities 1669 call ibte%compute_rta(cryst, dtset, dtfil, comm) 1670 1671 nkcalc = ibte%nkcalc 1672 nkibz = ibte%ebands%nkpt; nsppol = ibte%nsppol; ntemp = ibte%ntemp 1673 bmin = ibte%bmin; bmax = ibte%bmax 1674 !call wrtout(std_out, sjoin(" nkcalc", itoa(nkcalc), "bmin:", itoa(bmin), "bmax:", itoa(bmax))) 1675 1676 !call ibte%read_scattering() 1677 ! Loops and memory are distributed over k-points and collinear spins 1678 ABI_MALLOC(sr, (nkcalc, nsppol)) 1679 cnt = 0 1680 do spin=1,nsppol 1681 do ikcalc=1,nkcalc 1682 cnt = cnt + 1 1683 sr(ikcalc, spin)%rank = mod(cnt, nprocs) 1684 end do 1685 end do 1686 1687 call cwtime(cpu, wall, gflops, "start") 1688 #ifdef HAVE_NETCDF 1689 ! Master reads and sends data to the rank treating (ikcalc, spin). 1690 if (my_rank == master) then 1691 NCF_CHECK(nctk_open_read(ncid, dtfil%filsigephin, xmpi_comm_self)) 1692 end if 1693 1694 do spin=1,nsppol 1695 do ikcalc=1,nkcalc 1696 sr_p => sr(ikcalc, spin) 1697 receiver = sr_p%rank 1698 send_data = master /= receiver 1699 if (.not. any(my_rank == [master, receiver])) cycle 1700 !call wrtout(std_out, sjoin(" Sending data from my_rank:", itoa(my_rank), " to:", itoa(receiver))) 1701 1702 if (my_rank == master) then 1703 ! Get ncid of group used to store scattering rate for this k-point. 1704 ncerr = nf90_inq_ncid(ncid, strcat("srate_k", itoa(ikcalc), "_s", itoa(spin)), grp_ncid) 1705 if (ncerr /= NF90_NOERR) then 1706 ABI_ERROR("Cannot find collision terms in SIGEPH file. Rerun eph_task -4 step with ibte_prep 1.") 1707 end if 1708 NCF_CHECK(nctk_get_dim(grp_ncid, "nq_ibzk_eff", sr_p%nq_ibzk_eff)) 1709 NCF_CHECK(nctk_get_dim(grp_ncid, "nbsum", nbsum)) 1710 NCF_CHECK(nctk_get_dim(grp_ncid, "nbcalc", nbcalc)) 1711 NCF_CHECK(nctk_get_dim(grp_ncid, "lgk_nsym", sr_p%lgk_nsym)) 1712 dims = [sr_p%nq_ibzk_eff, nbsum, nbcalc, sr_p%lgk_nsym] 1713 end if 1714 1715 if (send_data) then 1716 tag = size(dims) 1717 if (my_rank == master) call xmpi_send(dims, receiver, tag, comm, ierr) 1718 if (my_rank == receiver) then 1719 call xmpi_recv(dims, master, tag, comm, ierr) 1720 sr_p%nq_ibzk_eff = dims(1); nbsum = dims(2); nbcalc = dims(3); sr_p%lgk_nsym = dims(4) 1721 end if 1722 end if 1723 1724 ! Note that the size along the (n, m) axis does not depend on the kcalc index. 1725 ABI_CALLOC(sr_p%vals, (sr_p%nq_ibzk_eff, bmin:bmax, bmin:bmax, ntemp)) 1726 ABI_MALLOC(sr_p%kq_symtab, (6, sr_p%nq_ibzk_eff)) 1727 ABI_MALLOC(sr_p%lgk_sym2glob, (2, sr_p%lgk_nsym)) 1728 1729 if (my_rank == master) then 1730 NCF_CHECK(nf90_get_var(grp_ncid, nctk_idname(grp_ncid, "kq_symtab"), sr_p%kq_symtab)) 1731 NCF_CHECK(nf90_get_var(grp_ncid, nctk_idname(grp_ncid, "lgk_sym2glob"), sr_p%lgk_sym2glob)) 1732 1733 ! Note that on file, we have: 1734 ! 1735 ! nctkarr_t("srate", "dp", "nq_ibzk_eff, nbsum, nbcalc, ntemp") 1736 ! 1737 ! but in terms of n, m indices we have that the: 1738 ! 1739 ! n index: bstart_ks bstop_ks 1740 ! m index: bsum_start up to bsum_stop to account for phonon emission/absorption. 1741 ! 1742 ! so we have to insert the values in bmin:bmax slice. 1743 ! TODO: Recheck this part. 1744 ABI_MALLOC(grp_srate, (sr_p%nq_ibzk_eff, nbsum, nbcalc, ntemp)) 1745 NCF_CHECK(nf90_get_var(grp_ncid, nctk_idname(grp_ncid, "srate"), grp_srate)) 1746 ii = ibte%bstart_ks(ikcalc, spin) 1747 jj = ibte%bstop_ks(ikcalc, spin) 1748 sr_p%vals(:, bmin:bmax, ii:jj, :) = grp_srate(:, 1:bmax-bmin+1, 1:nbcalc, :) 1749 ABI_SFREE(grp_srate) 1750 end if 1751 1752 if (send_data) then 1753 tag = size(sr_p%vals) 1754 if (my_rank == master) then 1755 tag = tag + 1; call xmpi_send(sr_p%vals, receiver, tag, comm, ierr) 1756 tag = tag + 1; call xmpi_send(sr_p%kq_symtab, receiver, tag, comm, ierr) 1757 tag = tag + 1; call xmpi_send(sr_p%lgk_sym2glob, receiver, tag, comm, ierr) 1758 end if 1759 if (my_rank == receiver) then 1760 tag = tag + 1; call xmpi_recv(sr_p%vals, master, tag, comm, ierr) 1761 tag = tag + 1; call xmpi_recv(sr_p%kq_symtab, master, tag, comm, ierr) 1762 tag = tag + 1; call xmpi_recv(sr_p%lgk_sym2glob, master, tag, comm, ierr) 1763 end if 1764 end if 1765 1766 if (send_data .and. my_rank /= receiver) call free_sr_ks(ikcalc, spin) 1767 1768 end do ! spin 1769 end do ! ikcalc 1770 1771 if (my_rank == master) then 1772 NCF_CHECK(nf90_close(ncid)) 1773 end if 1774 #endif 1775 call cwtime_report(" sigeph IO", cpu, wall, gflops) 1776 1777 ! Solve the linearized BTE with B = 0. 1778 ! 1779 ! F_\nk = e df/de_\nk v_\nk \tau^0 + \tau^0 \sum_{mq} Srate_{nk,mq} F_{m,k+q} 1780 ! 1781 ! where F is a vector in Cartesian coordinates and tau^0 is the SERTA relaxation time. 1782 ! 1783 ! Take advantage of the following symmetry properties: 1784 ! 1785 ! 1. F_k = F_Sk. 1786 ! 2. F_{-k} = -F_k if TR symmetry. 1787 ! 3. The q-space integration is reduced to the IBZ(k) using the symmetries of the little group of k. 1788 1789 !call ibte%solve_ibte(solver_type=1) 1790 1791 bsize = bmax - bmin + 1 1792 ABI_CALLOC(fkn_in, (3, nkibz, bmin:bmax, nsppol)) 1793 ABI_CALLOC(fkn_out, (3, nkibz, bmin:bmax, nsppol)) 1794 ABI_CALLOC(fkn_serta, (3, nkibz, bmin:bmax, nsppol)) 1795 ABI_CALLOC(taukn_serta, (3, nkibz, bmin:bmax, nsppol)) 1796 ABI_MALLOC(ibte_sigma, (3, 3, 2, nsppol, ntemp)) 1797 ABI_MALLOC(ibte_mob, (3, 3, 2, nsppol, ntemp)) 1798 ABI_MALLOC(converged, (ntemp)) 1799 1800 abs_tol = dtset%ibte_abs_tol 1801 1802 ! If the fermi level is inside the gap, F_k is gonna be very small 1803 ! hence once should use a much smaller tolerance to converge. 1804 ! According to numerical tests, a reasonable value of abs_tol can be estimated from the free carrier density using: 1805 ! 1806 ! 1e-20 * e_density (in cm**-3) 1807 ! 1808 ! These are the numerical values used to derive the fit: 1809 1810 ! max_adiff = np.array([9e-12, 5.6e-6, 2.7e-6, 8.2e-60, 2.9e-46, 6.9e-6, 2.9e-8, 9.5E+00, 7.3E-04, 9.2E-04, 8.4E-04]) 1811 ! e_density = np.array([0.97e8, 0.86e13, 0.26e14, 0.1e-40, 0.89e-26, 0.45e14, 0.28e12, 0.10E+19, 0.12E+16, 0.10E+15, 0.90E+14]) 1812 1813 if (abs_tol <= zero) then 1814 rtmp = minval(ibte%n_ehst, mask=ibte%n_ehst > zero) * ibte%nsppol 1815 abs_tol = 1e-20 * rtmp / cryst%ucvol / Bohr_cm**3 1816 call wrtout(std_out, " Input ibte_abs_tol <= zero ==> computing abs tolerance from carrier density") 1817 call wrtout(std_out, " using: abs_tol = 1e-20 * e_density (in cm**-3)") 1818 call wrtout(std_out, sjoin(" abs_tol:", ftoa(abs_tol), " from carrier_density:", ftoa(rtmp / cryst%ucvol / Bohr_cm**3))) 1819 end if 1820 1821 if (my_rank == master) then 1822 path = strcat(dtfil%filnam_ds(4), "_RTA.nc") 1823 call wrtout(unts, ch10//sjoin("- Writing IBTE transport results to:", path)) 1824 NCF_CHECK(nctk_open_modify(ncid, path , xmpi_comm_self)) 1825 1826 ncerr = nctk_def_dims(ncid, [ & 1827 nctkdim_t("nkibz", nkibz), nctkdim_t("bsize", bsize)], defmode=.True.) 1828 NCF_CHECK(ncerr) 1829 1830 ncerr = nctk_def_arrays(ncid, [ & 1831 nctkarr_t('fkn_out_sigma', "dp", "three, nkibz, bsize, nsppol, ntemp"), & 1832 nctkarr_t('ibte_sigma', "dp", "three, three, two, nsppol, ntemp"), & 1833 nctkarr_t('ibte_mob', "dp", "three, three, two, nsppol, ntemp"), & 1834 nctkarr_t('ibte_rho', "dp", "three, three, ntemp") & 1835 ], defmode=.True.) 1836 NCF_CHECK(ncerr) 1837 1838 NCF_CHECK(nctk_set_datamode(ncid)) 1839 end if 1840 1841 cnt = 0 1842 btype = 1 1843 do itemp=1,ntemp 1844 !do itemp=ntemp, 1, -1 1845 cnt = cnt + 1 1846 kT = ibte%kTmesh(itemp) 1847 mu_e = ibte%eph_mu_e(itemp) 1848 sig_p => ibte_sigma(:,:,:,:,itemp) 1849 mob_p => ibte_mob(:,:,:,:,itemp) 1850 1851 ! Precompute tau_serta and fkn_serta for this T: f^'_nk v_\nk * \tau^0 1852 do spin=1,nsppol 1853 do ikcalc=1,nkcalc 1854 !ik_ibz = ibte%kcalc2ibz(ikcalc, 1) 1855 ik_ibz = ibte%kcalc2ebands(1, ikcalc) 1856 do band_k=ibte%bstart_ks(ikcalc, spin), ibte%bstop_ks(ikcalc, spin) 1857 lw_nk = ibte%linewidths(itemp, band_k, ik_ibz, spin, 1) ! 1 --> SERTA linewidths. 1858 call safe_div(one, two * lw_nk, zero, tau_nk) 1859 taukn_serta(:, ik_ibz, band_k, spin) = tau_nk 1860 e_nk = ebands%eig(band_k, ik_ibz, spin) 1861 dfde_nk = occ_dfde(e_nk, kT, mu_e) 1862 btype_fact = one 1863 if (btype == 2) btype_fact = (e_nk - mu_e) ! / (Kt / kb_HaK) 1864 if (btype == 3) btype_fact = (e_nk - mu_e) ** 2 ! / (Kt / kb_HaK) 1865 fkn_serta(:, ik_ibz, band_k, spin) = tau_nk * dfde_nk * btype_fact * ibte%vbks(:, band_k, ik_ibz, spin) 1866 end do 1867 end do 1868 end do 1869 1870 call wrtout(std_out, sjoin(" Begin IBTE looop for itemp:", itoa(itemp), ", KT:", ftoa(kT / kb_HaK), "[K]"), & 1871 pre_newlines=1, newlines=1) 1872 1873 if (ibte%assume_gap) then 1874 write(msg, "(a5,1x,a9,*(1x, a16))")" ITER", "max_adiff", "mobility_e+h", "sum_k(df_k)" 1875 else 1876 write(msg, "(a5,1x,a9,*(1x, a16))")" ITER", "max_adiff", "conductivity", "sum_k(df_k)" 1877 end if 1878 call wrtout(std_out, msg) 1879 1880 ! iter = 0 --> Compute SERTA transport tensors just for initial reference. 1881 call ibte_calc_tensors(ibte, cryst, itemp, kT, mu_e, fkn_serta, onsager, sig_p, mob_p, fsum_eh, comm) 1882 1883 ! Print mobility for semiconductors, conductivity for metals. 1884 if (ibte%assume_gap) then 1885 do spin=1,nsppol 1886 mat33 = sum(mob_p(:,:,:,spin), dim=3) 1887 write(msg, "(i5,1x,es9.1, *(1x, f16.2))") & 1888 0, zero, mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin)) 1889 end do 1890 else 1891 do spin=1,nsppol 1892 mat33 = sum(sig_p(:,:,:,spin), dim=3) 1893 write(msg, "(i5,1x,es9.1, *(1x, es16.6))") & 1894 0, zero, mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin)) 1895 end do 1896 end if 1897 call wrtout(std_out, msg) 1898 1899 fkn_in = fkn_serta 1900 ! TODO: B-field 1901 ! Initialize fkn_in either from SERTA or from previous T. 1902 !if (cnt == 1) fkn_in = fkn_serta 1903 !if (cnt > 1 ) fkn_in = fkn_out 1904 fkn_out = zero 1905 1906 ! Begin iterative solver. 1907 iter_loop: do iter=1,dtset%ibte_niter 1908 1909 ! Loop over the nk index in F_nk. 1910 do spin=1,nsppol 1911 do ikcalc=1,nkcalc 1912 sr_p => sr(ikcalc, spin) 1913 if (sr_p%rank /= my_rank) cycle ! MPI parallelism 1914 ik_ibz = ibte%kcalc2ebands(1, ikcalc) 1915 do band_k=ibte%bstart_ks(ikcalc, spin), ibte%bstop_ks(ikcalc, spin) 1916 1917 ! Summing over the q-points in the effective IBZ(k) and the m band index. 1918 ! Results stored in vec3. Integration weights are already included. 1919 vec3 = zero 1920 do band_sum=ibte%bmin, ibte%bmax 1921 do iq_sum=1, sr_p%nq_ibzk_eff 1922 ikq_ibz = sr_p%kq_symtab(1, iq_sum); isym_kq = sr_p%kq_symtab(2, iq_sum) 1923 trev_kq = sr_p%kq_symtab(6, iq_sum) !; g0_kq = sr_p%kq_symtab(3:5, iq_sum) 1924 ! Build F_{m,k+q} in the effective IBZ(k) from fkn_in using symmetries (need k+q --> IBZ map) 1925 ! Use transpose(R) because we are using the tables for the wavefunctions 1926 ! In this case listkk has been called with symrec and use_symrec=False 1927 ! so q_bz = S^T q_ibz where S is the isym_kq symmetry 1928 ! vkq = matmul(transpose(cryst%symrel_cart(:,:,isym_kq)), vkq) 1929 mat33 = transpose(cryst%symrel_cart(:,:,isym_kq)) 1930 f_kq = matmul(mat33, fkn_in(:, ikq_ibz, band_sum, spin)) 1931 if (trev_kq == 1) f_kq = -f_kq 1932 vec3 = vec3 + sr_p%vals(iq_sum, band_sum, band_k, itemp) * f_kq(:) 1933 end do ! iq_sum 1934 end do ! band_k 1935 1936 ! Symmetrize intermediate results using the operations of the litte group of k. 1937 sym_vec = zero 1938 do isym_lgk=1,sr_p%lgk_nsym 1939 isym = sr_p%lgk_sym2glob(1, isym_lgk) 1940 itime = sr_p%lgk_sym2glob(2, isym_lgk) 1941 mat33 = transpose(cryst%symrel_cart(:,:,isym)) 1942 !if(itime == 1) mat33 = -mat33 ! FIXME: here there's a different convention for TR used in m_lgroup 1943 if (itime == 2) mat33 = -mat33 1944 sym_vec = sym_vec + matmul(mat33, vec3) 1945 end do 1946 sym_vec = taukn_serta(:, ik_ibz, band_k, spin) * sym_vec / sr_p%lgk_nsym 1947 fkn_out(:, ik_ibz, band_k, spin) = fkn_serta(:, ik_ibz, band_k, spin) + sym_vec 1948 1949 end do ! band_k 1950 end do ! ikcalc 1951 end do ! spin 1952 1953 call xmpi_sum(fkn_out, comm, ierr) 1954 if (my_rank == master) then 1955 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "fkn_out_sigma"), fkn_out, start=[1,1,1,1,itemp])) 1956 end if 1957 1958 do spin=1,nsppol 1959 max_adiff_spin(spin) = maxval(abs(fkn_out(:,:,:,spin) - fkn_in(:,:,:,spin))) 1960 end do 1961 max_adiff = maxval(max_adiff_spin) 1962 1963 ! Compute transport tensors from fkn_out 1964 call ibte_calc_tensors(ibte, cryst, itemp, kT, mu_e, fkn_out, onsager, sig_p, mob_p, fsum_eh, comm) 1965 1966 ! Print mobility for semiconductors or conductivity for metals. 1967 if (ibte%assume_gap) then 1968 do spin=1,nsppol 1969 mat33 = sum(mob_p(:,:,:,spin), dim=3) 1970 write(msg, "(i5,1x,es9.1,*(1x, f16.2))") & 1971 iter, max_adiff_spin(spin), mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin)) 1972 end do 1973 else 1974 do spin=1,nsppol 1975 mat33 = sum(sig_p(:,:,:,spin), dim=3) 1976 write(msg, "(i5,1x,es9.1,*(1x, es16.6))") & 1977 iter, max_adiff_spin(spin), mat33(1,1), mat33(2,2), mat33(3,3), sum(fsum_eh(:,:,spin)) 1978 end do 1979 end if 1980 call wrtout(std_out, msg) 1981 1982 ! Check for convergence by testing max_k |F_k^i - F_k^{i-1}|. 1983 converged(itemp) = max_adiff < abs_tol 1984 if (converged(itemp)) then 1985 call wrtout(std_out, sjoin(" IBTE solver converged after:", itoa(iter), & 1986 "iterations within ibte_abs_tol:", ftoa(abs_tol)), pre_newlines=1) 1987 exit iter_loop 1988 else 1989 ! Linear mixing of fkn_in and fkn_out. 1990 fkn_in = (one - dtset%ibte_alpha_mix) * fkn_in + dtset%ibte_alpha_mix * fkn_out 1991 fkn_out = zero 1992 end if 1993 1994 end do iter_loop 1995 1996 if (.not. converged(itemp)) then 1997 msg = sjoin("Not converged after:", itoa(dtset%ibte_niter), "max iterations") 1998 call wrtout(ab_out, msg, pre_newlines=1, newlines=1) 1999 ABI_WARNING(msg) 2000 end if 2001 end do ! itemp 2002 2003 ABI_MALLOC(ibte_rho, (3, 3, ntemp)) 2004 do itemp=1,ntemp 2005 work33 = sum(ibte_sigma(:,:,:,1,itemp), dim=3) 2006 if (ibte%nsppol == 2) work33 = work33 + sum(ibte_sigma(:,:,:,2,itemp), dim=3) 2007 call inv33(work33, mat33) 2008 ibte_rho(:, :, itemp) = 1e+6_dp * mat33 2009 end do 2010 2011 if (my_rank == master) then 2012 ! Write final results to main output. 2013 components = ["xx", "yy", "zz"] 2014 if (ibte%assume_gap) then 2015 ! SemiConductor 2016 do ii=1,3 2017 call wrtout(unts, sjoin(" Cartesian component of IBTE mobility tensor:", components(ii))) 2018 write(msg, "(a16,2(a32),a16)") 'Temperature [K]', 'e/h density [cm^-3]', 'e/h mobility [cm^2/Vs]', "Converged" 2019 call wrtout(unts, msg) 2020 2021 do spin=1,ibte%nsppol 2022 if (ibte%nsppol == 2) call wrtout(unts, sjoin(" For spin:", stoa(spin)), newlines=1) 2023 2024 do itemp=1,ibte%ntemp 2025 write(msg,"(f16.2,2e16.2,2f16.2,a16)") & 2026 ibte%kTmesh(itemp) / kb_HaK, & 2027 ibte%n_ehst(1, spin, itemp) / cryst%ucvol / Bohr_cm **3, & 2028 ibte%n_ehst(2, spin, itemp) / cryst%ucvol / Bohr_cm **3, & 2029 ibte_mob(ii, ii, 1, spin, itemp), ibte_mob(ii, ii, 2, spin, itemp), & 2030 yesno(converged(itemp)) 2031 call wrtout(unts, msg) 2032 end do ! itemp 2033 end do ! spin 2034 call wrtout(unts, ch10) 2035 end do ! ii 2036 2037 else 2038 ! Metals. Print conductivity (spin resolved) and resistivity (no spin resolved) 2039 do ii=1,2 2040 if (ii == 1) msg = " Conductivity [Siemens cm^-1] using IBTE" 2041 if (ii == 2) msg = " Resistivity [micro-Ohm cm] using IBTE" 2042 call wrtout(unts, msg) 2043 2044 nsp = ibte%nsppol; if (ii == 2) nsp = 1 2045 do spin=1,nsp 2046 if (ibte%nsppol == 2) call wrtout(unts, sjoin(" For spin:", stoa(spin)), newlines=1) 2047 write(msg, "(5a16)") 'Temperature (K)', 'xx', 'yy', 'zz', "Converged" 2048 call wrtout(unts, msg) 2049 do itemp=1,ibte%ntemp 2050 if (ii == 1) then 2051 mat33 = sum(ibte_sigma(:,:,:,spin,itemp), dim=3) 2052 else 2053 mat33 = ibte_rho(:,:,itemp) 2054 end if 2055 write(msg,"(f16.2,3e16.6,a16)") & 2056 ibte%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3), yesno(converged(itemp)) 2057 call wrtout(unts, msg) 2058 end do ! itemp 2059 end do ! spin 2060 call wrtout(unts, ch10) 2061 end do ! ii 2062 end if 2063 2064 !pre = "_IBTE" 2065 !call ibte%write_tensor(dtset, irta, "sigma", ibte%sigma(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SIGMA")) 2066 !call ibte%write_tensor(dtset, irta, "seebeck", ibte%seebeck(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SBK")) 2067 !call ibte%write_tensor(dtset, irta, "kappa", ibte%kappa(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_KAPPA")) 2068 !call ibte%write_tensor(dtset, irta, "zte", ibte%zte(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_ZTE")) 2069 !call ibte%write_tensor(dtset, irta, "pi", ibte%pi(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_PI")) 2070 2071 ! Print IBTE results to stdout and other external txt files (for the test suite) 2072 !call ibte%print_rta_txt_files(cryst, dtset, dtfil) 2073 ! Creates the netcdf file used to store the results of the calculation. 2074 #ifdef HAVE_NETCDF 2075 !path = strcat(dtfil%filnam_ds(4), "_RTA.nc") 2076 !call wrtout(unts, ch10//sjoin("- Writing IBTE transport results to:", path)) 2077 !NCF_CHECK(nctk_open_modify(ncid, path , xmpi_comm_self)) 2078 2079 !ncerr = nctk_def_arrays(ncid, [ & 2080 ! nctkarr_t('ibte_sigma', "dp", "three, three, two, nsppol, ntemp"), & 2081 ! nctkarr_t('ibte_mob', "dp", "three, three, two, nsppol, ntemp"), & 2082 ! nctkarr_t('ibte_rho', "dp", "three, three, ntemp") & 2083 !], defmode=.True.) 2084 !NCF_CHECK(ncerr) 2085 2086 ! Write data. 2087 NCF_CHECK(nctk_set_datamode(ncid)) 2088 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_sigma"), ibte_sigma)) 2089 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_mob"), ibte_mob)) 2090 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ibte_rho"), ibte_rho)) 2091 NCF_CHECK(nf90_close(ncid)) 2092 #endif 2093 2094 end if ! master 2095 2096 ! Free memory 2097 ABI_FREE(fkn_serta) 2098 ABI_FREE(taukn_serta) 2099 ABI_FREE(fkn_in) 2100 ABI_FREE(fkn_out) 2101 ABI_FREE(ibte_sigma) 2102 ABI_FREE(ibte_mob) 2103 ABI_FREE(ibte_rho) 2104 ABI_FREE(converged) 2105 2106 do spin=1,nsppol 2107 do ikcalc=1,nkcalc 2108 call free_sr_ks(ikcalc, spin) 2109 end do 2110 end do 2111 ABI_FREE(sr) 2112 2113 call ibte%free() 2114 2115 contains 2116 2117 subroutine free_sr_ks(ikc, isp) 2118 integer,intent(in) :: ikc, isp 2119 ABI_SFREE(sr(ikc, isp)%vals) 2120 ABI_SFREE(sr(ikc, isp)%lgk_sym2glob) 2121 ABI_SFREE(sr(ikc, isp)%kq_symtab) 2122 end subroutine free_sr_ks 2123 2124 end subroutine ibte_driver
m_rta/print_rta_txt_files [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
print_rta_txt_files
FUNCTION
INPUTS
cryst<crystal_t>=Crystalline structure dtset<dataset_type>=All input variables for this dataset. dtfil<datafiles_type>=variables related to files.
SOURCE
1342 subroutine print_rta_txt_files(self, cryst, dtset, dtfil) 1343 1344 !Arguments -------------------------------------- 1345 class(rta_t),intent(in) :: self 1346 type(crystal_t),intent(in) :: cryst 1347 type(dataset_type),intent(in) :: dtset 1348 type(datafiles_type),intent(in) :: dtfil 1349 1350 !Local variables -------------------------------- 1351 integer :: itemp, spin, irta, ii, nsp 1352 character(len=500) :: msg, pre, rta_type 1353 integer :: unts(2) 1354 character(len=2) :: components(3) 1355 real(dp) :: mat33(3,3) ! work33(3,3), 1356 1357 !************************************************************************ 1358 1359 unts = [std_out, ab_out] 1360 call wrtout(unts, ch10//' Transport (RTA) calculation results:', newlines=1) 1361 components = ["xx", "yy", "zz"] 1362 1363 do irta=1,self%nrta 1364 if (irta == 1) rta_type = "SERTA" 1365 if (irta == 2) rta_type = "MRTA" 1366 1367 if (self%assume_gap) then 1368 ! SemiConductor 1369 do ii=1,3 1370 call wrtout(unts, sjoin(" Cartesian component of", rta_type, "mobility tensor:", components(ii))) 1371 write(msg, "(a16,a32,a32)") 'Temperature [K]', 'e/h density [cm^-3]', 'e/h mobility [cm^2/Vs]' 1372 call wrtout(unts, msg) 1373 1374 do spin=1,self%nsppol 1375 if (self%nsppol == 2) call wrtout(unts, sjoin(" For spin:", stoa(spin)), newlines=1) 1376 1377 do itemp=1,self%ntemp 1378 write(msg,"(f16.2,2e16.2,2f16.2)") & 1379 self%kTmesh(itemp) / kb_HaK, & 1380 self%n_ehst(1, spin, itemp) / cryst%ucvol / Bohr_cm**3, & 1381 self%n_ehst(2, spin, itemp) / cryst%ucvol / Bohr_cm**3, & 1382 self%mobility_mu(ii, ii, 1, spin, itemp, irta), & 1383 self%mobility_mu(ii, ii, 2, spin, itemp, irta) 1384 call wrtout(unts, msg) 1385 end do ! itemp 1386 end do ! spin 1387 call wrtout(unts, ch10) 1388 end do ! ii 1389 1390 else 1391 ! Metals. Print conductivity (spin resolved) and resistivity (no spin resolved) 1392 do ii=1,2 1393 if (ii == 1) msg = sjoin(" Conductivity [Siemens cm^-1] using ", rta_type, "approximation") 1394 if (ii == 2) msg = sjoin(" Resistivity [micro-Ohm cm] using ", rta_type, "approximation") 1395 call wrtout(unts, msg) 1396 1397 nsp = self%nsppol; if (ii == 2) nsp = 1 1398 do spin=1,nsp 1399 if (nsp == 2) call wrtout(unts, sjoin(" For spin:", stoa(spin)), newlines=1) 1400 write(msg, "(4a16)") 'Temperature (K)', 'xx', 'yy', 'zz' 1401 call wrtout(unts, msg) 1402 do itemp=1,self%ntemp 1403 if (ii == 1) then 1404 mat33 = self%conductivity(:,:,itemp,spin,irta) 1405 else 1406 mat33 = self%resistivity(:,:,itemp,irta) 1407 end if 1408 write(msg,"(f16.2,3e16.6)") self%kTmesh(itemp) / kb_HaK, mat33(1,1), mat33(2,2), mat33(3,3) 1409 call wrtout(unts, msg) 1410 end do !itemp 1411 end do !spin 1412 call wrtout(unts, ch10) 1413 end do 1414 end if 1415 1416 end do ! irta 1417 1418 do irta=1,self%nrta 1419 select case (irta) 1420 case (1) 1421 pre = "_SERTA" 1422 case (2) 1423 pre = "_MRTA" 1424 case default 1425 ABI_ERROR(sjoin("Don't know how to handle irta:", itoa(irta))) 1426 end select 1427 call self%write_tensor(dtset, irta, "sigma", self%sigma(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SIGMA")) 1428 call self%write_tensor(dtset, irta, "seebeck", self%seebeck(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_SBK")) 1429 call self%write_tensor(dtset, irta, "kappa", self%kappa(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_KAPPA")) 1430 call self%write_tensor(dtset, irta, "zte", self%zte(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_ZTE")) 1431 call self%write_tensor(dtset, irta, "pi", self%pi(:,:,:,:,:,irta), strcat(dtfil%filnam_ds(4), pre, "_PI")) 1432 end do 1433 1434 end subroutine print_rta_txt_files
m_rta/rta_driver [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
rta_driver
FUNCTION
Driver to compute transport properties within the RTA.
INPUTS
dtfil<datafiles_type>=variables related to files. ngfftc(18)=Coarse FFT mesh. dtset<dataset_type>=All input variables for this dataset. ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...) cryst<crystal_t>=Crystalline structure pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. comm=MPI communicator.
SOURCE
284 subroutine rta_driver(dtfil, ngfftc, dtset, ebands, cryst, pawtab, psps, comm) 285 286 !Arguments ------------------------------------ 287 !scalars 288 integer, intent(in) :: comm 289 type(datafiles_type),intent(in) :: dtfil 290 type(dataset_type),intent(in) :: dtset 291 type(crystal_t),intent(in) :: cryst 292 type(ebands_t),intent(in) :: ebands 293 type(pseudopotential_type),intent(in) :: psps 294 !arrays 295 integer,intent(in) :: ngfftc(18) 296 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 297 298 !Local variables ------------------------------ 299 type(rta_t) :: rta 300 !arrays 301 integer :: unts(2) 302 303 ! ************************************************************************* 304 305 unts = [std_out, ab_out] 306 call wrtout(unts, ch10//' Entering transport RTA computation driver.') 307 call wrtout(unts, sjoin("- Reading carrier lifetimes from:", dtfil%filsigephin), newlines=1, do_flush=.True.) 308 309 ! Initialize RTA object 310 rta = rta_new(dtset, dtfil, ngfftc, cryst, ebands, pawtab, psps, comm) 311 312 ! Compute RTA transport quantities 313 call rta%compute_rta(cryst, dtset, dtfil, comm) 314 315 call rta%free() 316 317 end subroutine rta_driver
m_rta/rta_free [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
rta_free
FUNCTION
Free dynamic memory.
INPUTS
SOURCE
1529 subroutine rta_free(self) 1530 1531 !Arguments -------------------------------------- 1532 class(rta_t),intent(inout) :: self 1533 1534 ABI_SFREE(self%n) 1535 ABI_SFREE(self%vv_dos) 1536 ABI_SFREE(self%vvtau_dos) 1537 ABI_SFREE(self%tau_dos) 1538 ABI_SFREE(self%bstart_ks) 1539 ABI_SFREE(self%bstop_ks) 1540 ABI_SFREE(self%nbcalc_ks) 1541 !ABI_SFREE(self%kcalc2ibz) 1542 ABI_SFREE(self%kcalc2ebands) 1543 ABI_SFREE(self%kTmesh) 1544 ABI_SFREE(self%eminmax_spin) 1545 ABI_SFREE(self%eph_mu_e) 1546 ABI_SFREE(self%transport_mu_e) 1547 ABI_SFREE(self%vbks) 1548 ABI_SFREE(self%linewidths) 1549 ABI_SFREE(self%l0) 1550 ABI_SFREE(self%l1) 1551 ABI_SFREE(self%l2) 1552 ABI_SFREE(self%sigma) 1553 ABI_SFREE(self%mobility) 1554 ABI_SFREE(self%conductivity) 1555 ABI_SFREE(self%resistivity) 1556 ABI_SFREE(self%seebeck) 1557 ABI_SFREE(self%kappa) 1558 ABI_SFREE(self%zte) 1559 ABI_SFREE(self%pi) 1560 ABI_SFREE(self%mobility_mu) 1561 ABI_SFREE(self%conductivity_mu) 1562 ABI_SFREE(self%n_ehst) 1563 1564 call ebands_free(self%ebands) 1565 call self%gaps%free() 1566 call self%edos%free() 1567 1568 end subroutine rta_free
m_rta/rta_ncwrite [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
rta_ncwrite
FUNCTION
INPUTS
cryst<crystal_t>=Crystalline structure dtset<dataset_type>=All input variables for this dataset. ncid=Netcdf file handle.
SOURCE
1213 subroutine rta_ncwrite(self, cryst, dtset, ncid) 1214 1215 !Arguments -------------------------------------- 1216 class(rta_t),intent(in) :: self 1217 type(crystal_t),intent(in) :: cryst 1218 type(dataset_type),intent(in) :: dtset 1219 integer,intent(in) :: ncid 1220 1221 !Local variables -------------------------------- 1222 integer :: ncerr, ii 1223 real(dp) :: cpu, wall, gflops 1224 real(dp) :: work(dtset%nsppol) 1225 1226 !************************************************************************ 1227 1228 call cwtime(cpu, wall, gflops, "start") 1229 1230 #ifdef HAVE_NETCDF 1231 ! Write to netcdf file 1232 NCF_CHECK(cryst%ncwrite(ncid)) 1233 NCF_CHECK(ebands_ncwrite(self%ebands, ncid)) 1234 NCF_CHECK(self%edos%ncwrite(ncid)) 1235 1236 !nctk_copy from sigeph? 1237 ! nctkarr_t("eph_ngqpt_fine", "int", "three"), & 1238 1239 ncerr = nctk_def_dims(ncid, [ & 1240 nctkdim_t("ntemp", self%ntemp), nctkdim_t("nrta", self%nrta), nctkdim_t("nsppol", self%nsppol)], defmode=.True.) 1241 NCF_CHECK(ncerr) 1242 1243 ncerr = nctk_def_arrays(ncid, [ & 1244 nctkarr_t('transport_ngkpt', "int", "three"), & 1245 nctkarr_t('sigma_erange', "dp", "two"), & 1246 nctkarr_t('kTmesh', "dp", "ntemp"), & 1247 nctkarr_t('transport_mu_e', "dp", "ntemp"), & 1248 nctkarr_t('n_ehst', "dp", "two, nsppol, ntemp"), & 1249 nctkarr_t('eph_mu_e', "dp", "ntemp"), & 1250 nctkarr_t('vb_max', "dp", "nsppol"), & 1251 nctkarr_t('cb_min', "dp", "nsppol"), & 1252 nctkarr_t('vv_dos', "dp", "edos_nw, three, three, nsppol"), & 1253 nctkarr_t('vvtau_dos', "dp", "edos_nw, three, three, ntemp, nsppol, nrta"), & 1254 nctkarr_t('tau_dos', "dp", "edos_nw, ntemp, nsppol, nrta"), & 1255 nctkarr_t('L0', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1256 nctkarr_t('L1', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1257 nctkarr_t('L2', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1258 nctkarr_t('sigma', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1259 nctkarr_t('kappa', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1260 nctkarr_t('zte', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1261 nctkarr_t('seebeck', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1262 nctkarr_t('pi', "dp", "three, three, edos_nw, nsppol, ntemp, nrta"), & 1263 nctkarr_t('mobility', "dp", "three, three, edos_nw, ntemp, two, nsppol, nrta"), & 1264 nctkarr_t('conductivity', "dp", "three, three, ntemp, nsppol, nrta"), & 1265 nctkarr_t('resistivity', "dp", "three, three, ntemp, nrta"), & 1266 nctkarr_t('N', "dp", "edos_nw, ntemp, two"), & 1267 !nctkarr_t('conductivity_mu',"dp", "three, three, two, nsppol, ntemp, nrta")], & 1268 nctkarr_t('mobility_mu', "dp", "three, three, two, nsppol, ntemp, nrta")], & 1269 defmode=.True.) 1270 NCF_CHECK(ncerr) 1271 1272 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "assume_gap"]) 1273 NCF_CHECK(ncerr) 1274 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: & 1275 "eph_extrael", "eph_fermie", "transport_extrael", "transport_fermie"]) 1276 NCF_CHECK(ncerr) 1277 1278 ! Write data. 1279 ii = 0; if (self%assume_gap) ii = 1 1280 ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: "assume_gap"], [ii], datamode=.True.) 1281 1282 ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: & 1283 "eph_extrael", "eph_fermie", "transport_extrael", "transport_fermie"], & 1284 [self%eph_extrael, self%eph_fermie, self%transport_extrael, self%transport_fermie]) 1285 NCF_CHECK(ncerr) 1286 1287 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "transport_ngkpt"), dtset%transport_ngkpt)) 1288 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "sigma_erange"), dtset%sigma_erange)) 1289 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kTmesh"), self%kTmesh)) 1290 if (self%assume_gap) then 1291 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vb_max"), self%gaps%vb_max)) 1292 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "cb_min"), self%gaps%cb_min)) 1293 else 1294 ! Set vbm and cbm to fermie if metal. 1295 work(:) = self%ebands%fermie 1296 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vb_max"), work)) 1297 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "cb_min"), work)) 1298 end if 1299 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_mu_e"), self%eph_mu_e)) 1300 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "transport_mu_e"), self%transport_mu_e)) 1301 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "n_ehst"), self%n_ehst)) 1302 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vv_dos"), self%vv_dos)) 1303 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vvtau_dos"), self%vvtau_dos)) 1304 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "tau_dos"), self%tau_dos)) 1305 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "L0"), self%l0)) 1306 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "L1"), self%l1)) 1307 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "L2"), self%l2)) 1308 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "sigma"), self%sigma)) 1309 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kappa"), self%kappa)) 1310 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "zte"), self%zte)) 1311 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "seebeck"), self%seebeck)) 1312 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "pi"), self%pi)) 1313 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "N"), self%n)) 1314 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "mobility"), self%mobility)) 1315 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "conductivity"), self%conductivity)) 1316 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "resistivity"), self%resistivity)) 1317 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "mobility_mu"), self%mobility_mu)) 1318 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "conductivity_mu"), self%conductivity_mu)) 1319 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "resistivity_mu"), self%resistivity_mu)) 1320 #endif 1321 1322 call cwtime_report(" rta_ncwrite", cpu, wall, gflops) 1323 1324 end subroutine rta_ncwrite
m_rta/rta_new [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
rta_new
FUNCTION
Build object to compute RTA transport quantities.
INPUTS
dtset<dataset_type>=All input variables for this dataset. dtfil<datafiles_type>=variables related to files. cryst<crystal_t>=Crystalline structure ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...) comm=MPI communicator.
SOURCE
338 type(rta_t) function rta_new(dtset, dtfil, ngfftc, cryst, ebands, pawtab, psps, comm) result (new) 339 340 !Arguments ------------------------------------- 341 integer, intent(in) :: comm 342 type(dataset_type),intent(in) :: dtset 343 type(datafiles_type),intent(in) :: dtfil 344 type(crystal_t),intent(in) :: cryst 345 type(ebands_t),intent(in) :: ebands 346 type(pseudopotential_type),intent(in) :: psps 347 !arrays 348 integer,intent(in) :: ngfftc(18) 349 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 350 351 !Local variables ------------------------------ 352 integer,parameter :: sppoldbl1 = 1, master = 0 353 integer :: ierr, spin, nprocs, my_rank, timrev, ik_ibz, ib, irta, itemp, ndat, nsppol, idat, mband, ikpt 354 real(dp) :: cpu, wall, gflops 355 character(len=500) :: msg 356 character(len=fnlen) :: wfk_fname_dense 357 type(ebands_t) :: tmp_ebands, ebands_dense 358 type(klinterp_t) :: klinterp 359 type(ddkstore_t) :: ds 360 type(sigmaph_t) :: sigmaph 361 type(krank_t) :: krank 362 !arrays 363 integer :: kptrlatt(3,3), unts(2), sigma_ngkpt(3) 364 integer,allocatable :: indkk(:,:) 365 real(dp) :: extrael_fermie(2), sigma_erange(2) 366 real(dp),allocatable :: values_bksd(:,:,:,:), vals_bsd(:,:,:), tmp_array4(:,:,:,:), tmp_array5(:,:,:,:,:) 367 368 !************************************************************************ 369 370 call cwtime(cpu, wall, gflops, "start") 371 unts = [std_out, ab_out] 372 373 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 374 375 ! Use sigma_erange to understand if we are dealing with a metal or a semiconductor. 376 new%assume_gap = (.not. all(dtset%sigma_erange < zero) .or. dtset%gw_qprange /= 0) 377 378 ! Read data from SIGEPH file. 379 sigmaph = sigmaph_read(dtfil%filsigephin, dtset, xmpi_comm_self, msg, ierr, keep_open=.true., & 380 extrael_fermie=extrael_fermie, sigma_ngkpt=sigma_ngkpt, sigma_erange=sigma_erange) 381 ABI_CHECK(ierr == 0, msg) 382 383 !if (any(sigma_erange /= zero)) then 384 ! ABI_CHECK(all(dtset%sigma_erange /= zero), "sigma_erange is required in input with a value compatible with SIGEPH.nc") 385 ! ! Make sure that the two values are consistent 386 ! ! Cannot switch to metallic case if SigmaPH file was produced assuming gapped-system 387 ! if (.not. (all(dtset%sigma_erange < zero) .eqv. all(sigma_erange < zero))) then 388 ! ABI_ERROR("The values of sigma_erange from input and SIGEPH are not compatible") 389 ! end if 390 !end if 391 392 ! How many RTA approximations have we computed in sigmaph? (SERTA, MRTA ...?) 393 new%nrta = 2; if (sigmaph%mrta == 0) new%nrta = 1 394 395 ! Copy important arrays from sigmaph file. 396 ! Allocate temperature arrays (use same values as the ones used in the SIGEPH calculation). 397 new%ntemp = sigmaph%ntemp 398 call alloc_copy(sigmaph%kTmesh, new%kTmesh) 399 400 new%nkcalc = sigmaph%nkcalc 401 call alloc_copy(sigmaph%bstart_ks, new%bstart_ks) 402 call alloc_copy(sigmaph%bstop_ks, new%bstop_ks) 403 call alloc_copy(sigmaph%nbcalc_ks, new%nbcalc_ks) 404 !call alloc_copy(sigmaph%kcalc2ibz, new%kcalc2ibz) 405 406 new%bmin = minval(sigmaph%bstart_ks); new%bmax = maxval(sigmaph%bstop_ks) 407 !new%bmin = 1; new%bmax = ebands%mband ! This for debugging purposes, results should not change 408 new%bsize = new%bmax - new%bmin + 1 409 410 new%nsppol = ebands%nsppol; new%nspinor = ebands%nspinor 411 nsppol = new%nsppol 412 413 ABI_MALLOC(new%eminmax_spin, (2, nsppol)) 414 new%eminmax_spin = ebands_get_minmax(ebands, "eig") 415 416 if (new%assume_gap) then 417 ! Get gaps 418 new%gaps = ebands_get_gaps(ebands, ierr) 419 if (ierr /= 0) then 420 do spin=1, nsppol 421 ABI_WARNING(trim(new%gaps%errmsg_spin(spin))) 422 new%gaps%vb_max(spin) = ebands%fermie - 1 * eV_Ha 423 new%gaps%cb_min(spin) = ebands%fermie + 1 * eV_Ha 424 end do 425 !ABI_ERROR("ebands_get_gaps returned non-zero exit status. See above warning messages...") 426 ABI_WARNING("ebands_get_gaps returned non-zero exit status. See above warning messages...") 427 end if 428 if (my_rank == master) call new%gaps%print(unit=std_out); call new%gaps%print(unit=ab_out) 429 end if 430 431 ! ================================================= 432 ! Read lifetimes and new%ebands from SIGEPH.nc file 433 ! ================================================= 434 ! After this point we have: 435 ! 436 ! vbks(3, bmin:bmax, nkpt, nsppol) 437 ! linewidths(self%ntemp, bmin:bmax, nkpt, nsppol, 2) 438 ! 439 if (any(dtset%sigma_ngkpt /= 0)) then 440 ! If integrals are computed with the sigma_ngkpt k-mesh, we need to downsample ebands. 441 !call wrtout(unts, sjoin(" SIGMAPH file used sigma_ngkpt:", ltoa(sigma_ngkpt))) 442 call wrtout(unts, sjoin(" Computing integrals with downsampled sigma_ngkpt:", ltoa(dtset%sigma_ngkpt))) 443 kptrlatt = 0 444 kptrlatt(1,1) = dtset%sigma_ngkpt(1); kptrlatt(2,2) = dtset%sigma_ngkpt(2); kptrlatt(3,3) = dtset%sigma_ngkpt(3) 445 446 tmp_ebands = ebands_downsample(ebands, cryst, kptrlatt, dtset%sigma_nshiftk, dtset%sigma_shiftk) 447 new%ebands = sigmaph%get_ebands(cryst, tmp_ebands, [new%bmin, new%bmax], & 448 new%kcalc2ebands, new%linewidths, new%vbks, xmpi_comm_self) 449 call ebands_free(tmp_ebands) 450 else 451 !call wrtout(unts, sjoin(" Computing integrals with SIGEPH k-mesh:", ebands_kmesh2str(ebands)) 452 new%ebands = sigmaph%get_ebands(cryst, ebands, [new%bmin, new%bmax], & 453 new%kcalc2ebands, new%linewidths, new%vbks, xmpi_comm_self) 454 kptrlatt = new%ebands%kptrlatt 455 end if 456 457 !print *, "linewidth_serta", maxval(abs(new%linewidths(:,:,:,:,1))) 458 !print *, "linewidth_mrta", maxval(abs(new%linewidths(:,:,:,:,2))) 459 !print *, "max velocities", maxval(abs(new%vbks)) 460 461 if ( & 462 dtset%useria == 888 .and. & 463 (dtset%getwfkfine /= 0 .or. dtset%irdwfkfine /= 0 .or. dtset%getwfkfine_filepath /= ABI_NOFILE)) then 464 465 ! In principle only getwfkfine_filepath is used here 466 wfk_fname_dense = trim(dtfil%fnameabi_wfkfine) 467 ABI_CHECK(nctk_try_fort_or_ncfile(wfk_fname_dense, msg) == 0, msg) 468 469 call wrtout(unts, " EPH double grid interpolation: will read energies from: "//trim(wfk_fname_dense), newlines=1) 470 mband = new%ebands%mband 471 472 !ebands_dense = wfk_read_ebands(wfk_fname_dense, comm) 473 474 tmp_ebands = wfk_read_ebands(wfk_fname_dense, comm) 475 ebands_dense = ebands_chop(tmp_ebands, 1, mband) 476 call ebands_free(tmp_ebands) 477 if (my_rank == master) then 478 write(std_out, *)" Using kptrlatt: ", ebands_dense%kptrlatt 479 write(std_out, *)" shiftk: ", ebands_dense%shiftk 480 end if 481 ABI_CHECK_IEQ(mband, ebands_dense%mband, "Inconsistent number of bands for the fine and dense grid:") 482 483 ! Compute v_{nk} on the dense grid in Cartesian coordinates. 484 ! vdiago(3, bmin:bmax, nkpt, nsppol) 485 ! NB: We select bands in [bmin:bmax] but all k-points in the IBZ are computed! 486 ds%only_diago = .True.; ds%bmin = new%bmin; ds%bmax = new%bmax; ds%mode = "cart" 487 call ds%compute_ddk(wfk_fname_dense, "", dtset, psps, pawtab, ngfftc, comm) 488 489 ! Transfer data to new%vbks 490 ABI_MOVE_ALLOC(ds%vdiago, new%vbks) 491 call ds%free() 492 !print *, "vbks:", new%vbks 493 494 ! Linear interpolation in k-space of the linewidths from input SIGEPH to the dense IBZ provided by fine WFK file. 495 ! First of all transfer linewidths to values_bksd to prepare call to klinterp_new. 496 ndat = new%ntemp * new%nrta 497 ABI_MALLOC(values_bksd, (new%bmin:new%bmax, new%ebands%nkpt, nsppol, ndat)) 498 499 do irta=1,new%nrta 500 do spin=1,nsppol 501 do ik_ibz=1,new%ebands%nkpt 502 do ib=new%bmin,new%bmax 503 do itemp=1,new%ntemp 504 idat = itemp + new%ntemp * (irta - 1) 505 values_bksd(ib, ik_ibz, spin, idat) = new%linewidths(itemp, ib, ik_ibz, spin, irta) ! - t(e) 506 end do 507 end do 508 end do 509 end do 510 end do 511 512 ! Build linear interpolator for linewidths (use only bsize bands) 513 klinterp = klinterp_new(cryst, new%ebands%kptrlatt, new%ebands%nshiftk, new%ebands%shiftk, new%ebands%kptopt, & 514 new%ebands%kptns, new%bsize, new%ebands%nkpt, nsppol, ndat, values_bksd, comm) 515 ABI_FREE(values_bksd) 516 517 ! HERE we re-malloc new%ebands and %linewidths on the fine k-mesh. 518 ! The call must be executed here, once klinterp has been built. 519 ! After this point we can use new%ebands to allocate stuff. 520 call ebands_move_alloc(ebands_dense, new%ebands) 521 522 ! Unlinke the ebands stored in SIGEPH, the eigens read from WFK_FINE have not been 523 ! shifted with the scissors operator or updated according to extrael_fermie so do it now. 524 call ephtk_update_ebands(dtset, new%ebands, "GS energies read from WFK_FINE") 525 526 ! And now interpolate linewidths on the fine k-mesh 527 ! Note: k-points that close to the edge of the pocket may get zero linewidths 528 ! One may fix the problem by using lw(e). 529 ABI_REMALLOC(new%linewidths, (new%ntemp, new%bmin:new%bmax, new%ebands%nkpt, nsppol, new%nrta)) 530 ABI_MALLOC(vals_bsd, (new%bmin:new%bmax, nsppol, ndat)) 531 532 ierr = 0 533 do ik_ibz=1,new%ebands%nkpt 534 535 call klinterp%eval_bsd(new%ebands%kptns(:, ik_ibz), vals_bsd) 536 !vals_bsd = vals_bsd + lw(e) 537 538 if (any(vals_bsd < zero)) then 539 ierr = ierr + 1 540 where (vals_bsd < zero) vals_bsd = zero 541 end if 542 543 ! Transfer data. 544 do spin=1,nsppol 545 do irta=1,new%nrta 546 do itemp=1,new%ntemp 547 idat = itemp + new%ntemp * (irta - 1) 548 do ib=new%bmin,new%bmax 549 new%linewidths(itemp, ib, ik_ibz, spin, irta) = vals_bsd(ib, spin, idat) 550 end do 551 end do 552 end do 553 end do 554 555 end do ! ik_ibz 556 557 if (ierr /= 0) then 558 ! This should never happen for linear interpolation. 559 ABI_WARNING(sjoin("Linear interpolation produced:", itoa(ierr), " k-points with negative linewidths")) 560 end if 561 562 ABI_FREE(vals_bsd) 563 call klinterp%free() 564 end if 565 566 ! FIXME: I think transport_ngkpt is buggy, wrong ne(T), weird zeros if MRTA ... 567 ! Do we really need this option? Can't we replace it with sigma_ngkpt and eph_task 7? 568 569 if (any(dtset%transport_ngkpt /= 0)) then 570 ! Perform further downsampling (usefull for debugging purposes) 571 call wrtout(unts, " Downsampling the k-mesh before computing transport:") 572 call wrtout(unts, sjoin(" Using transport_ngkpt: ", ltoa(dtset%transport_ngkpt))) 573 kptrlatt = 0 574 kptrlatt(1, 1) = dtset%transport_ngkpt(1) 575 kptrlatt(2, 2) = dtset%transport_ngkpt(2) 576 kptrlatt(3, 3) = dtset%transport_ngkpt(3) 577 tmp_ebands = ebands_downsample(new%ebands, cryst, kptrlatt, 1, [zero, zero, zero]) 578 579 ! Map the points of the downsampled bands to dense ebands 580 timrev = kpts_timrev_from_kptopt(ebands%kptopt) 581 ABI_MALLOC(indkk, (6, tmp_ebands%nkpt)) 582 583 krank = krank_from_kptrlatt(new%ebands%nkpt, new%ebands%kptns, new%ebands%kptrlatt, compute_invrank=.False.) 584 585 if (kpts_map("symrec", timrev, cryst, krank, tmp_ebands%nkpt, tmp_ebands%kptns, indkk) /= 0) then 586 write(msg, '(3a)' ) & 587 "Error while downsampling ebands in the transport driver",ch10, & 588 "The k-point could not be generated from a symmetrical one." 589 ABI_ERROR(msg) 590 end if 591 592 call krank%free() 593 594 ! Downsampling linewidths and velocities. 595 ABI_MOVE_ALLOC(new%linewidths, tmp_array5) 596 ABI_REMALLOC(new%linewidths, (new%ntemp, new%bmin:new%bmax, tmp_ebands%nkpt, nsppol, new%nrta)) 597 do ikpt=1,tmp_ebands%nkpt 598 new%linewidths(:,:,ikpt,:,:) = tmp_array5(:,:,indkk(1, ikpt),:,:) 599 end do 600 ABI_FREE(tmp_array5) 601 602 ABI_MOVE_ALLOC(new%vbks, tmp_array4) 603 ABI_REMALLOC(new%vbks, (3, new%bmin:new%bmax, tmp_ebands%nkpt, nsppol)) 604 do ikpt=1,tmp_ebands%nkpt 605 new%vbks(:,:,ikpt,:) = tmp_array4(:,:,indkk(1, ikpt),:) 606 end do 607 ABI_FREE(tmp_array4) 608 609 !print *, "after downsampling linewidths" 610 !print *, "linewidth_serta", maxval(abs(new%linewidths(:,:,:,:,1))) 611 !print *, "linewidth_mrta", maxval(abs(new%linewidths(:,:,:,:,2))) 612 613 ABI_FREE(indkk) 614 call ebands_move_alloc(tmp_ebands, new%ebands) 615 end if 616 617 ! Same doping case as in sigmaph file. 618 ABI_MALLOC(new%eph_mu_e, (new%ntemp)) 619 ABI_MALLOC(new%transport_mu_e, (new%ntemp)) 620 621 new%eph_extrael = extrael_fermie(1) 622 new%eph_fermie = extrael_fermie(2) 623 new%transport_fermie = dtset%eph_fermie 624 new%transport_extrael = dtset%eph_extrael 625 new%eph_mu_e = sigmaph%mu_e 626 new%transport_mu_e = sigmaph%mu_e 627 628 if (new%transport_fermie /= zero) new%transport_mu_e = new%transport_fermie 629 630 if (new%transport_fermie == zero .and. new%transport_extrael /= new%eph_extrael) then 631 632 if (new%transport_extrael /= new%eph_extrael) then 633 write(msg,'(2(a,e18.8),3a)') & 634 ' extrael from SIGEPH: ',new%transport_extrael, ' and input file: ',new%eph_extrael, "differ", ch10, & 635 ' Will recompute the chemical potential' 636 call wrtout(std_out, msg) 637 end if 638 639 ! Compute Fermi level for different T values. 640 call ebands_get_muT_with_fd(ebands, new%ntemp, new%kTmesh, dtset%spinmagntarget, dtset%prtvol, new%transport_mu_e, comm) 641 end if 642 643 ! TODO: Implement possible change of sigma_erange, useful for convergence studies 644 ! 1) Run sigmaph with relatively large sigma_erange. 645 ! 2) Decrease energy window in the trasport part to analyze the behaviour of transport tensors. 646 647 ! sigmaph is not needed anymore. Free it. 648 sigmaph%ncid = nctk_noid 649 call sigmaph%free() 650 651 call cwtime_report(" rta_new", cpu, wall, gflops) 652 653 end function rta_new
m_rta/rta_t [ Types ]
NAME
rta_t
FUNCTION
Container for transport quantities in the RTA
SOURCE
86 type,public :: rta_t 87 88 integer :: nsppol 89 ! Number of independent spin polarizations. 90 91 integer :: nspinor 92 ! Number of spinor components. 93 94 integer :: nkcalc 95 ! Number of computed k-points i.e. k-points inside the sigma_erange energy window. 96 97 integer :: ntemp 98 ! Number of temperatures. 99 100 integer :: nw 101 ! Number of energies (chemical potentials) at which transport quantities are computed 102 ! Same number of energies used in DOS. 103 104 integer :: bmin, bmax, bsize 105 ! Only bands between bmin and bmax are considered in the integrals 106 ! as we don't compute linewidths for all states. 107 ! bmin = minval(%bstart_ks); bmax = maxval(%bstop_ks) 108 ! bisze = bmax - bmin + 1 109 110 integer :: nrta 111 ! Number of relaxation-time approximations used (1 for SERTA, 2 for MRTA) 112 113 real(dp) :: eph_extrael 114 ! Extra electrons per unit cell used to compute SERTA lifetimes in sigmaph. 115 116 real(dp) :: eph_fermie 117 ! Fermi level specified in the input file when computing the SIGEPH file. 118 119 real(dp) :: transport_extrael 120 ! Extra electrons per unit cell specified in the input file when computing the SIGEPH file. 121 122 real(dp) :: transport_fermie 123 ! Fermi level specified in the input file when computing the SIGEPH file. 124 125 logical :: assume_gap 126 ! True if we are dealing with a semiconductor. 127 ! This parameter is initialized from the value of sigma_erange provided by the user. 128 129 integer,allocatable :: bstart_ks(:,:) 130 ! bstart_ks(nkcalc, nsppol) 131 ! Initial KS band index included in self-energy matrix elements for each k-point in kcalc. 132 ! Depends on spin because all degenerate states should be included when symmetries are used. 133 134 integer,allocatable :: bstop_ks(:,:) 135 ! bstop_ks(nkcalc, nsppol) 136 137 integer,allocatable :: nbcalc_ks(:,:) 138 ! nbcalc_ks(nkcalc, nsppol) 139 ! Number of bands included in self-energy matrix elements for each k-point in kcalc. 140 ! Depends on spin because all degenerate states should be included when symmetries are used. 141 142 !integer,allocatable :: kcalc2ibz(:,:) 143 !kcalc2ibz(nkcalc, 6)) 144 ! Mapping ikcalc --> IBZ as reported by listkk. 145 146 integer,allocatable :: kcalc2ebands(:,:) 147 ! Mapping ikcalc --> ebands IBZ 148 ! Note that this array is not necessarily equation to kcalc2ibz computed in sigmaph 149 ! because we may have used sigma_nkpt to downsample the initial nkpt mesh. 150 ! This array is computed in get_ebands and is equal to kcalc2ibz if sigma_nkpt == ngkpt 151 152 real(dp),allocatable :: kTmesh(:) 153 ! (%ntemp) 154 ! List of k * T temperatures at which to compute the transport 155 156 real(dp),allocatable :: eph_mu_e(:) 157 ! (%ntemp) 158 ! Chemical potential at this carrier concentration and temperature from sigeph (lifetime) 159 160 real(dp),allocatable :: transport_mu_e(:) 161 ! (%ntemp) 162 ! Chemical potential at this carrier concentration and temperature 163 164 real(dp),allocatable :: eminmax_spin(:,:) 165 ! (2, %nsppol)) 166 ! min/Max energy of the original ebands object 167 168 real(dp),allocatable :: linewidths(:,:,:,:,:) 169 ! (ntemp, bmin:bmax, nkpt, nsppol, nrta) 170 ! Linewidth in the IBZ computed in the SERTA/MRTA. 171 ! Non-zero only for the kcalc k-points. 172 173 real(dp),allocatable :: vbks(:,:,:,:) 174 ! (3, bmin:bmax, nkpt, nsppol)) 175 ! band velocity in Cartesian coordinates in the IBZ 176 ! Non-zero only for the kcalc k-points. 177 178 type(gaps_t) :: gaps 179 ! gaps of original ebands object. Only if assume_gap 180 181 type(ebands_t) :: ebands 182 ! bandstructure object used to compute the transport properties 183 ! Allocate using only the relevant bands for transport 184 ! including valence states to allow to compute different doping 185 186 type(edos_t) :: edos 187 ! electronic density of states 188 ! edos%mesh is the mesh used for vv_dos, vvtau_dos and tau_dos 189 ! (%nw) 190 191 real(dp),allocatable :: tau_dos(:,:,:,:) 192 ! tau(e) (isotropic average for tau_nk for SERTA and MRTA. 193 ! (nw, ntemp, nsppol, nrta) 194 195 real(dp),allocatable :: vv_dos(:,:,:,:) 196 ! (v x v) DOS 197 ! (nw, 3, 3, nsppol) 198 199 real(dp),allocatable :: vvtau_dos(:,:,:,:,:,:) 200 ! (v x v * tau) DOS 201 ! (nw, 3, 3, ntemp, nsppol, nrta) 202 203 real(dp),allocatable :: n_ehst(:,:,:) 204 ! (2, %nsppol, %ntemp) 205 ! Number of electrons (e) and holes (h) per unit cell 206 ! The first dimension is for electrons/holes. 207 ! If nsppol == 2, the second dimension is the number of e/h for spin else the total number of e/h summed over spins. 208 209 real(dp),allocatable :: l0(:,:,:,:,:,:), l1(:,:,:,:,:,:), l2(:,:,:,:,:,:) 210 ! (3, 3, nw, nsppol, ntemp, nrta) 211 ! Onsager coeficients in Cartesian coordinates 212 213 real(dp),allocatable :: sigma(:,:,:,:,:,:) 214 real(dp),allocatable :: seebeck(:,:,:,:,:,:) 215 real(dp),allocatable :: kappa(:,:,:,:,:,:) 216 real(dp),allocatable :: pi(:,:,:,:,:,:) 217 real(dp),allocatable :: zte(:,:,:,:,:,:) 218 ! (3, 3, nw, nsppol, ntemp, nrta) 219 ! Transport coefficients in Cartesian coordinates 220 221 real(dp),allocatable :: mobility(:,:,:,:,:,:,:) 222 ! Mobility 223 ! (3, 3, nw, %ntemp, 2, %nsppol, nrta) 224 ! 5-th index is for e-h 225 226 real(dp),allocatable :: conductivity(:,:,:,:,:) 227 ! Conductivity at the Fermi level 228 ! (3, 3, %ntemp, %nsppol, nrta) 229 230 real(dp),allocatable :: resistivity(:,:,:,:) 231 ! (3, 3, %ntemp, nrta) 232 233 real(dp),allocatable :: n(:,:,:) 234 ! (nw, ntemp, 2) carrier density for e/h (n/cm^3) 235 236 real(dp),allocatable :: mobility_mu(:,:,:,:,:,:) 237 ! (3, 3, 2, nsppol, ntemp, nrta) 238 ! mobility for electrons and holes (third dimension) at transport_mu_e(ntemp) 239 ! Third dimension is for electron/hole 240 241 real(dp),allocatable :: conductivity_mu(:,:,:,:,:,:) 242 ! (3, 3, 2, nsppol, ntemp, nrta) 243 ! Conductivity in Siemens * cm-1 244 ! computed by summing over k-points rather that by performing an energy integration). 245 246 contains 247 248 procedure :: compute_rta 249 procedure :: compute_rta_mobility 250 procedure :: print_rta_txt_files 251 procedure :: write_tensor 252 procedure :: free => rta_free 253 procedure :: rta_ncwrite 254 255 end type rta_t
m_rta/write_tensor [ Functions ]
[ Top ] [ m_rta ] [ Functions ]
NAME
FUNCTION
INPUTS
SOURCE
1448 subroutine write_tensor(self, dtset, irta, header, values, path) 1449 1450 !Arguments -------------------------------------- 1451 class(rta_t),intent(in) :: self 1452 type(dataset_type),intent(in) :: dtset 1453 integer,intent(in) :: irta 1454 character(len=*),intent(in) :: header 1455 real(dp),intent(in) :: values(:,:,:,:,:) 1456 character(len=*),intent(in) :: path 1457 1458 !Local variables -------------------------------- 1459 integer :: itemp, iw, ount 1460 character(len=500) :: msg, rta_type 1461 real(dp),allocatable :: tmp_values(:,:,:,:,:) 1462 1463 !************************************************************************ 1464 1465 if (open_file(trim(path), msg, newunit=ount, form="formatted", action="write", status='unknown') /= 0) then 1466 ABI_ERROR(msg) 1467 end if 1468 1469 if (irta == 1) rta_type = "RTA type: Self-energy relaxation time approximation (SERTA)" 1470 if (irta == 2) rta_type = "RTA type: Momentum relaxation time approximation (MRTA)" 1471 1472 ! write header 1473 write(ount, "(2a)")"# ", trim(header) 1474 write(ount, "(2a)")"# ", trim(rta_type) 1475 ! TODO: Units ? 1476 write(ount, "(a, 3(i0, 1x))")"#", dtset%transport_ngkpt 1477 write(ount, "(a)")"#" 1478 1479 ! This to improve portability of the unit tests. 1480 call alloc_copy(values, tmp_values) 1481 where (abs(values) > tol30) 1482 tmp_values = values 1483 else where 1484 tmp_values = zero 1485 end where 1486 1487 ! (nw, 3, 3, nsppol, ntemp) 1488 if (self%nsppol == 1) then 1489 do itemp=1, self%ntemp 1490 write(ount, "(/, a, 1x, f16.2)")"# T = ", self%kTmesh(itemp) / kb_HaK 1491 write(ount, "(a)")"# Energy [Ha], (xx, yx, zx, xy, yy, zy, xz, yz, zz) Cartesian components of tensor." 1492 do iw=1,self%nw 1493 write(ount, "(10(es16.6))")self%edos%mesh(iw), tmp_values(:, :, iw, 1, itemp) 1494 end do 1495 end do 1496 write(ount, "(a)")"" 1497 else 1498 do itemp=1, self%ntemp 1499 write(ount, "(/, a, 1x, f16.2)")"# T = ", self%kTmesh(itemp) / kb_HaK 1500 write(ount, "(a)") & 1501 "# Energy [Ha], (xx, yx, zx, xy, yy, zy, xz, yz, zz) Cartesian components of tensor for spin up followed by spin down." 1502 do iw=1,self%nw 1503 write(ount, "(19(es16.6))")self%edos%mesh(iw), tmp_values(:, :, iw, 1, itemp), tmp_values(:, :, iw, 2, itemp) 1504 end do 1505 end do 1506 write(ount, "(a)")"" 1507 end if 1508 1509 close(ount) 1510 1511 ABI_FREE(tmp_values) 1512 1513 end subroutine write_tensor