diff --git a/tp/calculate.py b/tp/calculate.py index 095722a..791e3a1 100644 --- a/tp/calculate.py +++ b/tp/calculate.py @@ -15,6 +15,8 @@ particle lifetime. mfp: particle mean free path. + fd_occupation: + fermion occupation. be_occupation: boson occupation. dfdde: @@ -191,6 +193,46 @@ def mfp(gamma, group_velocity, use_tprc=True): return mfp +def fd_occupation(energy, temperature, fermi_level, use_tprc=True): + """Calculates Fermi-Dirac occupation. + + Assumes singly-occupied bands: double as necessary. + + Arguments + --------- + + energy : array-like or float + energies (by default in eV). + temperature : array-like or float + temperature (by default in K). + fermi_level : array-like or float + fermi level (by default in eV). + + use_tprc : bool, optional + use custom unit conversions. Default: True. + + Returns + ------- + + array-like + occupations. + """ + + kb = physical_constants['Boltzmann constant in eV/K'][0] + + if use_tprc: + energy = to_tp('energy', energy) + temperature = to_tp('temperature', temperature) + fermi_level = to_tp('fermi_level', fermi_level) + + occupation = (np.exp(np.divide(np.add.outer(-fermi_level, energy), + kb * temperature[None, :, None, None])) + 1) ** -1 + + if use_tprc: + occupation = from_tp('occupation', occupation) + + return occupation + def be_occupation(frequency, temperature, use_tprc=True): """Calculates Bose-Einstein occupation. diff --git a/tp/cli/cli.py b/tp/cli/cli.py index 645df38..a550ae8 100644 --- a/tp/cli/cli.py +++ b/tp/cli/cli.py @@ -12,6 +12,8 @@ get amset get specific data from amset transport json or mesh h5. + occupation + get charge carrier occupation boltztrap get specific data from tp boltztrap hdf5. phono3py @@ -267,7 +269,77 @@ def get_amset(amset_file, quantity, dtype, doping, direction, temperature, spin, return -@get.command('boltztrap', no_args_is_help=True) +@get.command('occupation') +@inputs_function('amset_mesh_file', nargs=1) +@doping_type_option +@doping_option +@temperature_option +@click.option('--spin', + help='Spin direction.', + type=click.Choice(['up', 'down', 'avg'], case_sensitive=False), + default='avg', + show_default=True) +@click.option('-p', '--poscar', + help='POSCAR path. Required for --qpoint.', + type=click.Path(file_okay=True, dir_okay=False), + default='POSCAR', + show_default=True) + + +def get_occupation(amset_mesh_file, dtype, doping, temperature, spin, poscar): + """Prints AMSET occupation info at given conditions. + + Requires an AMSET mesh file. + """ + + from pymatgen.io.vasp.inputs import Poscar + from scipy.constants import physical_constants + kb = physical_constants['Boltzmann constant in eV/K'][0] + if 'energy' in tp.settings.conversions(): + kb *= tp.settings.conversions()['energy'] + if 'temperature' in tp.settings.conversions(): + kb /= tp.settings.conversions()['temperature'] + + data = tp.data.load.amset_mesh(amset_mesh_file, 'occupation energy vb_idx', + spin=spin, doping=dtype) + + v = Poscar.from_file(poscar).structure.lattice.volume + data = tp.data.utilities.resolve(data, 'occupation energy', + temperature=temperature, dtype=dtype, + doping=doping) + + avg_factor = 2 if spin in ['average', 'avg'] else 1 + cb = data['vb_idx'] + 1 + occ = np.average(data['occupation'], axis=1) + e = np.sum(occ[cb:]) + h = avg_factor * cb - np.sum(occ[:cb]) + e /= v * 1e-8 ** 3 + h /= v * 1e-8 ** 3 + c = tp.settings.conversions() + if 'doping' in c: + e *= c['doping'] + h *= c['doping'] + vbm = np.amax(data['energy'][cb-1]) + cbm = np.amin(data['energy'][cb]) + eg = cbm - vbm + + print('The bandgap is {:.3f} {}, and at {:d} {} 10kbT = {:.3f} {}.'.format( + eg, data['meta']['units']['energy'], + int(np.ceil(data['meta']['temperature'])), + data['meta']['units']['temperature'], + 10 * kb * int(np.ceil(data['meta']['temperature'])), + data['meta']['units']['energy'])) + print('Your chosen carrier concentration is {} {}.'.format( + data['meta']['doping'], data['meta']['units']['doping'])) + print('There are {:.3e} holes {} in the valence band and {:.3e} ' + 'electrons {} in the conduction band.'.format( + h, data['meta']['units']['doping'], + e, data['meta']['units']['doping'])) + + return + + +@get.command('boltztrap') @inputs_function('boltztrap_hdf5', nargs=1) @click.option('-q', '--quantity', help='Quantity to read.', @@ -735,8 +807,8 @@ def avg_rates(mesh_h5, mfp, total, x, crt, exclude, doping, direction, Requires AMSET mesh files. Plots scattering rates averaged over kpoints and weighted by the derivative of the Fermi-Dirac - distribution against tempearture or carrier concentration or both. - If one file is given, it will be used for both plots, or if more are + distribution against temperature or carrier concentration or both. + If one file is given, it will be used for both plots, or if two are specified the one with the most temperatures will be used for the temperature plot and the one with the most carrier concentrations will be used for the doping plot. x-limits only work for individual diff --git a/tp/data/load.py b/tp/data/load.py index 62db5c4..829ee6a 100644 --- a/tp/data/load.py +++ b/tp/data/load.py @@ -193,11 +193,12 @@ def amset_mesh(filename, quantities='all', doping='n', spin='avg'): channels, which are dealt with in the spin variable. Also accepts ibz_weights, the weights of the irreducible k-points, fd_weights, the weights of the energies wrt the - derivative of the Fermi-Dirac distribution, and - weighted_rates, scattering_rates weighted by fd_weights - and averaged over kpoints. "all" loads all quantities in - the file, which does not include ibz_weights, fd_weights or - weighted_rates. Loads dependent properties. Default: all. + derivative of the Fermi-Dirac distribution, weighted_rates, + scattering_rates weighted by fd_weights and averaged over + kpoints and occupation, the Fermi-Dirac occupation. + "all" loads all quantities in the file, which does not + include ibz_weights, fd_weights, weighted_rates or + occupation. Loads dependent properties. Default: all. doping : str, optional doing type (n or p). If there is more than one, defaults to @@ -205,7 +206,9 @@ def amset_mesh(filename, quantities='all', doping='n', spin='avg'): spin : str, optional spin. Accepts up, down or avg. If avg and there is only one - spin channel, selects that, else averages both. Default: avg. + spin channel, selects that, else averages both. If avg, + assumes non-spin separated bands and so multiplies occupation + by 2. Default: avg. Returns ------- @@ -224,14 +227,16 @@ def amset_mesh(filename, quantities='all', doping='n', spin='avg'): tnames = settings.to_tp() units = settings.units() dimensions = settings.dimensions() + dimensions['occupation'] = ['temperature', 'doping', 'kpoint', 'band'] if isinstance(quantities, str): quantities = quantities.split() quantities = [anames[q] if q in anames else q for q in quantities] # list of abbriviations and dependent quantites subs = {'weights': ['ibz_weights', 'fd_weights'], - 'weighted_rates': ['scattering_rates', 'weighted_rates'], + 'weighted_rates': ['scattering_rates', 'energies', 'fermi_levels', 'weighted_rates'], 'mean_free_path': ['scattering_rates', 'velocities', 'mean_free_path'], - 'weighted_mfp': ['scattering_rates', 'velocities', 'mean_free_path', 'weighted_mfp']} + 'weighted_mfp': ['scattering_rates', 'velocities', 'mean_free_path', 'energies', 'fermi_levels', 'weighted_mfp'], + 'occupation': ['energies', 'fermi_levels', 'occupation']} hasspin = ['energies', 'vb_index', 'scattering_rates', 'velocities'] l = len(quantities) @@ -256,6 +261,9 @@ def amset_mesh(filename, quantities='all', doping='n', spin='avg'): quantities.append(d) break + if 'temperatures' in quantities: + quantities.insert(0, quantities.pop(quantities.index('temperatures'))) + # load data def resolve_spin(data, q, spin): @@ -294,7 +302,7 @@ def resolve_spin(data, q, spin): if 'all' in quantities: q2 = quantities quantities = list(f.keys()) - for q in ['ibz_weights', 'fd_weights', 'weighted_rates']: + for q in ['ibz_weights', 'fd_weights', 'weighted_rates', 'occupation']: if q in q2: quantities.append(q) if 'doping' in quantities: @@ -310,30 +318,31 @@ def resolve_spin(data, q, spin): elif doping == 'p': di = np.where(d > 0)[0] + # spin2 so if avg selected doubles occupation if spin in ['avg', 'average']: if 'energies_up' in f and 'energies_down' in f: - spin = 'avg' + spin2 = 'avg' elif 'energies_up' in f: - spin = 'up' + spin2 = 'up' elif 'energies_down' in f: - spin = 'down' + spin2 = 'down' data = {'meta': {'doping_type': doping, 'electronic_source': 'amset', - 'spin': spin, + 'spin': spin2, 'units': {}, 'dimensions': {}}} for q in quantities: q2 = tnames[q] if q in tnames else q if q in hasspin: - data[q2] = resolve_spin(f, q, spin) + data[q2] = resolve_spin(f, q, spin2) elif q in f: data[q2] = f[q][()] elif q not in ['ibz_weights', 'fd_weights', 'weighted_rates', - 'mean_free_path', 'weighted_mfp']: + 'mean_free_path', 'weighted_mfp', 'occupation']: raise Exception('{} unrecognised. Quantity must be {}, ' 'ibz_weights, fd_weights, weighted_rates, ' - 'mean_free_path or weighted_mfp.' + 'mean_free_path or weighted_mfp or occupation.' ''.format(q, ', '.join(f))) if q == 'scattering_rates': data[q2] = [*list(data[q2]), list(np.sum(data[q2], axis=0))] @@ -341,7 +350,7 @@ def resolve_spin(data, q, spin): _, data['ibz_weights'] = np.unique(f['ir_to_full_kpoint_mapping'], return_counts=True) if q in ['fd_weights', 'weighted_rates', 'weighted_mfp']: - e = resolve_spin(f, 'energies', spin) + e = resolve_spin(f, 'energies', spin2) data['fd_weights'] = -tp.calculate.dfdde(e, f['fermi_levels'], f['temperatures'], f['doping'], @@ -366,6 +375,12 @@ def resolve_spin(data, q, spin): data[q2] = data['mean_free_path'] * \ data['normalised_weights'][None, :, :, :, :, None] data[q2] = data[q2].sum(axis=(3,4)) + if q == 'occupation': + data['occupation'] = tp.calculate.fd_occupation( + data['energy'], data['temperature'], + data['fermi_level']) + if spin in ['avg', 'average']: + data['occupation'] *= 2 for q2 in data: q = anames[q2] if q2 in anames else q2