Merge branch 'joss-paper' of https://github.com/SMTG-UCL/ThermoPlotter …

…into joss-paper

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.
 

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

tp/data/load.py

@@ -193,19 +193,22 @@ 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
  n, else this is ignored.
 
  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,38 +318,39 @@ 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))]
  if q in ['ibz_weights', 'weighted_rates', 'weighted_mfp']:
  _, 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