lscat

#! /usr/bin/env python3

from __future__ import print_function, division, absolute_import
import sys
import logging

from   scattering.fullscattering import FullScattering

'''
lscat: user script to create more elaborate scattering data
'''

def main_tempmode():
  ''' load LRTC energy file '''
  scatobj = FullScattering('tb-60-60-60-irr.hdf5')
  spins, nkp, nbands = scatobj.getDependencies()

  ''' define temperature grid in Kelvin - logarithmic steps or linear steps'''
  nt = 100 # steps
  temps = scatobj.defineTemperatures(tmin = 10, \
                                     tmax = 300, \
                                     nt   = nt, \
                                     tlog = True)

  ''' create scattering numpy array '''
  gamma = np.ones((nt, spins, nkp, nbands), dtype=np.float64)
  gamma *= 1e-7

  ''' load the array into the instantiated object '''
  scatobj.defineScatteringRates(gamma)

  ''' create output file '''
  scatobj.createOutput('scattering_file.hdf5')


def main_mumode():
  ''' load LRTC energy file '''
  scatobj = FullScattering('tb-60-60-60-irr.hdf5')
  spins, nkp, nbands = scatobj.getDependencies()

  ''' define chemical potential grid in eV
      if muabs: take the values as they are
      if not muabs: interpret them as difference to mu_DFT '''
  nt = 100 # steps
  mus = scatobj.defineChemicalPotentials(temp = 300, \
                                         mumin = -3, \
                                         mumax = +3, \
                                         nmu   = nt, \
                                         muabs = False)

  ''' create scattering numpy array unit: eV '''
  gamma = np.ones((nt, spins, nkp, nbands), dtype=np.float64)
  gamma *= 1e-7

  ''' create quasi particle weight numpy array: range (0,1] '''
  qpweight = np.ones_like(gamma, dtype=np.float64)
  qpweight *= 0.5 # Z = 0.5 throughout

  ''' load the array into the instantiated object '''
  scatobj.defineScatteringRates(scattering = gamma, \
                                qpweight   = qpweight, \
                                bandshift  = None)

  ''' create output file '''
  scatobj.createOutput('scattering_file.hdf5')

def main_momentum_dependence():
  scatobj = FullScattering('tb-60-60-60-irr.hdf5')
  spins, nkp, nbands = scatobj.getDependencies()
  nt = 100 # steps
  temps = scatobj.defineTemperatures(tmin = 10, \
                                     tmax = 300, \
                                     nt   = nt, \
                                     tlog = True)

  ''' get momentum grid '''
  kgrid = scatobj.getMomentumGrid() # shape: nkp, 3

  ''' create functional dependency for all temperatures, spins and bands '''
  gamma = np.zeros((nt, spins, nkp, nbands), dtype=np.float64)
  gammak = np.sum(np.sin(kgrid*2*np.pi)**2,axis=1)
  gamma[...] = gammak[None,None,:,None]
  gamma /= 3000.
  gamma[gamma<1e-10] = 1e-10 # introduce scattering cut-off

  scatobj.defineScatteringRates(gamma)
  scatobj.createOutput('scattering_file.hdf5')

def main_energy_direct_dependence():
  scatobj = FullScattering('tb-60-60-60-irr.hdf5')
  spins, nkp, nbands = scatobj.getDependencies()
  nt = 100 # steps
  temps = scatobj.defineTemperatures(tmin = 10, \
                                     tmax = 300, \
                                     nt   = nt, \
                                     tlog = True)

  ''' get energy grid and dft mu '''
  energies = scatobj.getEnergies()  # shape: spins, nkp, bands
  mudft    = scatobj.mudft

  ''' define gamma(spin,nkp,band) ... temperature independent
      if we have a direct dependence gamma(k) ~ energy(k)
      we can simply use the direct dependence '''
  gamma = np.zeros((nt, spins, nkp, nbands), dtype=np.float64)
  gamma[0,...]  = (energies-mudft)**2. / 10000.
  gamma[1:,...] = gamma[0,...]

  ''' truncate below 1e-10 eV '''
  gamma[gamma<1e-10] = 1e-10

  scatobj.defineScatteringRates(gamma)
  scatobj.createOutput('scattering_file.hdf5')

def main_energy_generic_dependence():
  scatobj = FullScattering('tb-60-60-60-irr.hdf5')
  spins, nkp, nbands = scatobj.getDependencies()
  nt = 100 # steps
  temps = scatobj.defineTemperatures(tmin = 10, \
                                     tmax = 300, \
                                     nt   = nt, \
                                     tlog = True)

  ''' get energy grid and dft mu '''
  energies = scatobj.getEnergies()  # shape: spins, nkp, bands
  mudft    = scatobj.mudft

  ''' define energy range '''
  enemin = energies.min()
  enemax = energies.max()
  print('min energy: {}\nmax energy: {}'.format(enemin,enemax))
  eneaxis = np.linspace(enemin,enemax,5001)

  ''' define generic gamma(energy)
      this method is necessary when the dependence is not direct
      but indirect via a generic distribution Gamma(epsilon) '''
  gammaeps = (eneaxis-mudft)**2 / 10000.

  try:
    import matplotlib.pyplot
    plt.plot(eneaxis,gammaeps)
    plt.show()
  except:
    pass

  ''' get the mapping eps(spin,nkp,band) <-> axis element '''
  where = np.argmin(np.abs(energies[None,...] - eneaxis[:,None,None,None]), axis=0)

  ''' define gamma(spin,nkp,band) ... temperature independent '''
  gamma = np.zeros((nt, spins, nkp, nbands), dtype=np.float64)
  gamma[0,...]  = gammaeps[where]
  gamma[1:,...] = gamma[0,...]

  ''' truncate below 1e-10 eV '''
  gamma[gamma<1e-10] = 1e-10

  scatobj.defineScatteringRates(gamma)
  scatobj.createOutput('scattering_file.hdf5')

def main_banddependence():
  scatobj = FullScattering('tb-60-60-60-irr.hdf5')
  spins, nkp, nbands = scatobj.getDependencies()
  nt = 100 # steps
  temps = scatobj.defineTemperatures(tmin = 10, \
                                     tmax = 300, \
                                     nt   = nt, \
                                     tlog = True)

  ''' scattering that is k-independent
      nkp <- 1  ... removes k-dependence (saves storage)
      first three bands different to the rest '''
  nkp = 1
  gamma = np.ones((nt, spins, nkp, nbands), dtype=np.float64)
  gamma[:,:,:,:3] = 1e-5
  gamma[:,:,:,3:] = 3e-5

  scatobj.defineScatteringRates(gamma)
  scatobj.createOutput('scattering_file.hdf5')

def main_tempddependence():
  scatobj = FullScattering('tb-60-60-60-irr.hdf5')
  spins, nkp, nbands = scatobj.getDependencies()
  nt = 100 # steps
  temps = scatobj.defineTemperatures(tmin = 10, \
                                     tmax = 300, \
                                     nt   = nt, \
                                     tlog = True)

  ''' define temperature dependence
      remove all momentum and band dependence '''
  nkp = nbands = 1
  gamma = np.zeros((nt, spins, nkp, nbands), dtype=np.float64)
  gammaT = 1e-8 + 1e-10 * temps
  gamma[...] = gammaT[:,None,None,None]

  scatobj.defineScatteringRates(gamma)
  scatobj.createOutput('scattering_file.hdf5')

def setLogger():
  ''' some logging '''
  class log_formatter(logging.Formatter):
    def format(self, record):
      if record.levelname == 'INFO':
        return record.msg
      else:
        return '{}: {}: {}'.format(record.filename, record.levelname, record.msg)
  logger = logging.getLogger()
  logger.setLevel(logging.DEBUG) # root logger
  console = logging.StreamHandler()
  console.setFormatter(log_formatter())
  console.setLevel(logging.INFO) # console logger print INFO, WARNING, ERROR and CRITICAL
  logger.addHandler(console)
  return logger

def main():
  template = '''#! /usr/bin/env python

from __future__ import print_function, division, absolute_import
import numpy as np
from  scattering.fullscattering import FullScattering

\'\'\' define object via lrtc energy file \'\'\'
scatobj = FullScattering('tb-60-60-60-irr.hdf5')

\'\'\' get dependencies \'\'\'
spins, nkp, nbands = scatobj.getDependencies()
kgrid              = scatobj.getMomentumGrid() # [ nkp, 3 ]
energies           = scatobj.getEnergies()     # [ spins, nkp, nbands ]
mudft              = scatobj.mudft             # DFT chemical potential

\'\'\' define temperature grid:
    tmin: minimum temperature [K]
    tmax: maximum temperature [K]
    nt  : number of temperature steps
    optional tlog: True  - logarithmic steps
                   False - linear steps
                   (default: False)
\'\'\'
nt = 100 # steps
scatobj.defineTemperatures(tmin = 10, \\
                           tmax = 300, \\
                           nt   = nt, \\
                           tlog = True)

\'\'\' or chemical potential grid:
    temp:  temperature [K]
    mumin: minimum chemical potential [eV]
    mumax: maximum chemical potential [eV]
    nmu  : number of chemical potential steps
    optional mlog : True  - logarithmic steps
                    False - linear steps
                    (default: False)
    optional muabs: True  - provided mumin/mumax values are absolute values
                    False - provided mumin/mumax values are relative to DFT chemical potential
                    (default: False)
\'\'\'
# scatobj.defineChemicalPotentials(temp = 300, \\
#                                  mumin = -3, \\
#                                  mumax = +3, \\
#                                  nmu   = nt, \\
#                                  mlog  = False, \\
#                                  muabs = False)

\'\'\' define Gamma, Z and energy shift dependence \'\'\'
gamma = np.ones((nt, spins, nkp, nbands), dtype=np.float64)
gamma *= 1e-3
qpweight = np.ones_like(gamma, dtype=np.float64)
bandshift = np.zeros_like(gamma, dtype=np.float64)

\'\'\' create output.
    note that qpweight and bandshift are optional parameters
    if qpweight is not provided it will be set to a default of Z=1 for all points
\'\'\'
scatobj.defineScatteringRates(gamma)
# scatobj.defineScatteringRates(gamma, qpweight=qpweight, bandshift=bandshift)
scatobj.createOutput('scattering_file.hdf5')'''
  with open('lscat_template','w') as scat:
    scat.write(template)
    scat.flush()
  print('\nSaved lscat_template into your currently active directory.')

if __name__ == '__main__':
  main()