"""
scattering states calculator
"""
from pathlib import Path
from nanotools.base import Base, Quantity
from nanotools.utils import to_quantity, ureg
from nanotools.totalenergy import TotalEnergy
from nanotools.twoprobe import TwoProbe, get_transport_dir
from nanotools.utils import dict_converter
from nanotools.jsonio import json_write
from matplotlib import pyplot as plt
import h5py
from nptyping import NDArray
import attr
from attr import field
import copy
import numpy as np
import os
import shutil
def converter_in_light(x):
if x is None:
return x
return np.array(x, dtype=np.int32)
[docs]
@attr.s
class PhotoCurrent(Base):
"""``Photocurrent`` class.
Attributes:
transport_axis
left_equal_right: boolean variable indicating whether the two lead structures are identical.
photon_energy (Quantity 1d array):
Input. Photon energies. Default unit is eV.
e.g. photon_energy = 1.55
photon_energy = [1.55, 1.56, 1.59]
elec_energy_resol (Quantity):
Input. Energy resolution for the energy grid on which photo-transmission and Green's
functions of electrons are computed. Default unit is eV.
in_light (integer array):
Input. in_light must have the same length as the number of atoms. 1 (0) indicates that
light is (not) shed on the corresponding atom. When in_light is None (default), all
atoms within the central device region are assumed to be under the light, and hence
subject to light induced excitation.
e.g.
in_light = [0 0 1 1 1 0 0]
photo_transmission (float 5d array):
computation result read out from nano_photo_out.h5, see _read_photo_transmission method.
5d array with each dimension corresponding to the electronic energy grid, photon energy,
k-points, spin, and lead.
theta (float scalar): input polarization angle, in radial unit.
phi (float scalar): input polarization angle, in radial unit.
polar_basis_vec (float 2d array):
Input. Two rows of vectors indicating the basis frame for light polarization.
e.g.
polar_basis_vec = [ [0,1,0] , [0,0,1] ]
polar_vector (complex 1d array):
input polarization vector of light.
e.g.
polar_vector = [1.0, 0.0, 1j]
If polar_vector is not specified by user, it will be calculated according to
a0 = cos(theta) * cos(phi) - 1j * sin(theta) * sin(phi)
a1 = sin(theta) * cos(phi) + 1j * cos(theta) * sin(phi)
polar_vector = a0 * e0 + a1 * e1
where e0, e1 are basis-vectors from polar_basis_vec
"""
center: TotalEnergy = attr.ib(
converter=lambda d: dict_converter(d, TotalEnergy),
validator=attr.validators.instance_of(TotalEnergy),
)
left: TotalEnergy = attr.ib(
default=None,
converter=attr.converters.optional(lambda d: dict_converter(d, TotalEnergy)),
validator=attr.validators.optional(attr.validators.instance_of(TotalEnergy)),
)
right: TotalEnergy = attr.ib(
default=None,
converter=attr.converters.optional(lambda d: dict_converter(d, TotalEnergy)),
validator=attr.validators.optional(attr.validators.instance_of(TotalEnergy)),
)
left_equal_right: bool = attr.ib(default=False) # are lead structures identical?
transport_axis: int = attr.ib(
default=-1, validator=attr.validators.in_([-1, 0, 1, 2])
)
photon_energy: Quantity = attr.ib(
default=np.array([0.0]) * ureg.eV,
converter=lambda x: to_quantity(x, "eV", shape=(-1)),
validator=attr.validators.instance_of(Quantity),
)
elec_energy_resol: Quantity = field(
default=0.025 * ureg.eV,
converter=attr.converters.optional(lambda x: to_quantity(x, "eV")),
validator=attr.validators.optional(
attr.validators.instance_of(Quantity),
),
)
in_light: np.ndarray = attr.ib(
default=None,
converter=attr.converters.optional(converter_in_light),
validator=attr.validators.optional(attr.validators.instance_of(np.ndarray)),
)
photo_transmission: NDArray = field(
default=None,
converter=attr.converters.optional(lambda x: np.array(x).ravel(order="F")),
validator=attr.validators.optional(attr.validators.instance_of(NDArray)),
)
polar_vector: NDArray = field(
default=None,
converter=attr.converters.optional(
lambda x: np.array(x).ravel(order="F").astype(np.cdouble)
),
validator=attr.validators.optional(attr.validators.instance_of(NDArray)),
)
polar_vector_real: NDArray = field(
default=None,
converter=attr.converters.optional(lambda x: np.array(x).ravel(order="F")),
validator=attr.validators.optional(attr.validators.instance_of(NDArray)),
)
polar_vector_imag: NDArray = field(
default=None,
converter=attr.converters.optional(lambda x: np.array(x).ravel(order="F")),
validator=attr.validators.optional(attr.validators.instance_of(NDArray)),
)
polar_basis_vec: NDArray = field(
default=None,
converter=attr.converters.optional(
lambda x: np.array(x).reshape((-1, 3), order="F")
),
validator=attr.validators.optional(attr.validators.instance_of(NDArray)),
)
theta: float = attr.ib(
default=0.0,
converter=float,
validator=attr.validators.instance_of(float),
)
phi: float = attr.ib(
default=0.0,
converter=float,
validator=attr.validators.instance_of(float),
)
classname: str = attr.ib()
@classname.default
def _classname_default_value(self):
return self.__class__.__name__
def __attrs_post_init__(self):
if self.transport_axis < 0:
self.transport_axis = get_transport_dir(self.center.system.cell.boundary)
if self.transport_axis not in [0, 1, 2]:
raise Exception("transport direction not specified.")
self.center = copy.deepcopy(self.center)
if self.left is None:
self.left = copy.deepcopy(self.center)
else:
self.left = copy.deepcopy(self.left)
if self.right is None:
self.left_equal_right = True
if self.left_equal_right:
print("Left and right lead structures are same.")
self.right = copy.deepcopy(self.left)
else:
self.right = copy.deepcopy(self.right)
# self.center.solver.restart.densityPath = "center_out.h5"
self._post_init()
self.polar_vector_real = np.real(self.polar_vector)
self.polar_vector_imag = np.imag(self.polar_vector)
def _post_init(self):
if self.polar_vector is not None:
return
if self.polar_basis_vec is None:
raise Exception("Must specify basis vectors [ [e0], [e1] ].")
e0 = self.polar_basis_vec[0, :]
e0 = e0 / np.linalg.norm(e0)
e1 = self.polar_basis_vec[1, :]
e1 = e1 / np.linalg.norm(e1)
if np.abs(np.dot(e0, e1)) > 1.0e-8:
raise Exception("The two basis vectors must be orthogonal.")
self.polar_vector = get_unit_vec_from_angles(e0, e1, self.theta, self.phi)
[docs]
@classmethod
def from_twoprobe(cls, twoprb, **kwargs):
"""initializes instance from TwoProbe object"""
if not isinstance(twoprb, TwoProbe):
raise Exception("Reading from not a TwoProbe object.")
calc = cls(
center=twoprb.center,
left=twoprb.left,
right=twoprb.right,
transport_axis=twoprb.transport_axis,
left_equal_right=twoprb.left_equal_right,
**kwargs,
)
return calc
def set_elec_energy_resol(self, resolution):
self.elec_energy_resol = to_quantity(resolution, "eV")
if self.center.system.pop.eta < (0.1 * self.elec_energy_resol):
print(
"Warning: eta for Green Function computation is way smaller than elec_energy_resol ..."
)
print("consider larger eta or adjusting elec_energy_resol down to eta.")
def set_eta(self, eta):
self.left.system.pop.set_eta(eta)
self.right.system.pop.set_eta(eta)
self.center.system.pop.set_eta(eta)
if self.center.system.pop.eta < (0.1 * self.elec_energy_resol):
print(
"Warning: eta for Green Function computation is way smaller than elec_energy_resol ..."
)
print("consider larger eta or adjusting elec_energy_resol down to eta.")
def set_photon_energy(self, energy):
self.photon_energy = to_quantity(energy, "eV", shape=(-1))
[docs]
def solve(self, photon_energy=None, input="nano_photo_in", output="nano_photo_out"):
"""Calculates transmission.
This method triggers the nanodcalplus_trsm executable.
Result is stored in self.dos.transmission
Args:
photon_energy (float / iterable): energy grid over which the transmission is to be calculated.
eg. energies=0.1
eg. energies=[0.1,0.2,0.3]
"""
self.polar_vector_real = np.real(self.polar_vector)
self.polar_vector_imag = np.imag(self.polar_vector)
if photon_energy is not None:
self.set_photon_energy(photon_energy)
if not (Path(self.center.solver.restart.densityPath).is_file()):
raise Exception("densityPath file not found for central region.")
if not (Path(self.left.solver.restart.densityPath).is_file()):
raise Exception("densityPath file not found for left lead.")
if not self.left_equal_right: # if different structure
if not (Path(self.right.solver.restart.densityPath).is_file()):
raise Exception("densityPath file not found for right lead.")
inputname = os.path.splitext(input)[0] + ".json"
self.left.solver.mpidist = self.center.solver.mpidist
self.right.solver.mpidist = self.center.solver.mpidist
self.write(inputname)
command, binname = self.center.solver.cmd.get_cmd("photocurr")
ret = command(inputname)
ret.check_returncode()
shutil.move("nano_photo_out.json", output + ".json")
self._update(output + ".json")
self.set_units("si")
self._read_photo_transmission()
def plot_photo_transmission(
self, photon=0, k=0, s=0, lead="left", show=True, save_to=None
):
self._read_photo_transmission()
fld = self.photo_transmission
photo_cutoff_expo = 7.0
bias = self.center.system.pop.bias.magnitude
sig = self.center.system.pop.sigma.magnitude
phot_e_max = np.amax(self.photon_energy.magnitude)
dE = self.elec_energy_resol.magnitude
Emin = min(-bias, 0) - photo_cutoff_expo * sig - phot_e_max
Emax = max(-bias, 0) + photo_cutoff_expo * sig + phot_e_max
egrid = np.array(range(np.rint((Emax - Emin) / dE).astype(int) + 1)) * dE + Emin
fig = plt.figure()
if lead == "left":
y_fld = fld[:, photon, k, s, 0]
else:
y_fld = fld[:, photon, k, s, 1]
plt.plot(egrid, y_fld)
plt.xlabel(f"Energy ({self.elec_energy_resol.u})")
plt.ylabel("Photo-Transmission")
plt.show()
fig.savefig("photo_transmission.png")
if save_to is not None:
dep = {"Energy": egrid, "PhotoTransmission": y_fld}
json_write(save_to, dep)
return fig
[docs]
def get_response(self):
"""Calculate the photo-response coefficient according to Nanodcal definition (named
"photocurrent" in Nanodcal output).
When calling this method, the photo-transmission calculation is assumed completed already.
You should find nano_photo_out.json and nano_photo_out.h5 in your work directory.
e.g.
cal = PhotoCurrent.read("nano_photo_out.json")
res = cal.get_response()
Returns:
a dimensionless 3d array, with each dimension corresponding to photon-energy, spin,
and lead.
"""
self._read_photo_transmission()
alpha = 0.0072973525693 # fine structure constant
# energy integration
tr = self.photo_transmission
tr = np.sum(tr, axis=0) # --> (photon, k, s, lead)
dE = self.elec_energy_resol.to("eV").magnitude
phe = self.photon_energy.to("eV").magnitude
# Nanodcal convention: multiply transmission with wght
# calculatephotocurrent.m
# decomphcurrent(jjphe,:,:,:,:,:) = decomphcurrent(jjphe,:,:,:,:,:)*(cPC.alpha/hbaromega(jjphe));
wght = alpha * (dE / phe)
tr = tr * wght[:, np.newaxis, np.newaxis, np.newaxis]
# k integration
wght = np.array(self.center.system.kpoint.kwght)
k4d = wght[np.newaxis, :, np.newaxis, np.newaxis]
tr = tr * k4d
tr = np.sum(tr, axis=1) # --> (photon, s, lead)
if self.center.system.hamiltonian.ispin == 1:
tr = 2.0 * tr # spin degenerate
return tr
[docs]
def get_photo_current(
self, flux=1.0 * ureg.watt / ureg.meter**2, mu=1.0, epsilon=1.0
):
"""Calculate photo-charge current.
When calling this method, the photo-transmission calculation is assumed completed already.
You should find nano_photo_out.json and nano_photo_out.h5 in your work directory.
e.g.
cal = PhotoCurrent.read("nano_photo_out.json")
c = cal.get_photo_current(flux=3.0 * ureg.watt / ureg.meter**2, mu=0.7, epsilon=1.5)
Args:
flux (power / area): energy flux of incident light
mu: relative magnetic permeability of material
epsilon: relative electric permittivity of material
Returns:
charge current, as 3d array, with each dimension corresponding to photon-energy,
spin, and lead. Its unit can be ampere, ampere/m, or ampere/m^2 depending on the
device cross-section.
"""
a0 = 5.29177210903e-11 # * ureg.meter # Bohr radius
respo = self.get_response()
flux_unit = ureg.watt / ureg.meter**2
pref = flux.to(flux_unit).magnitude / self.photon_energy.to("eV").magnitude
pref = pref * a0 * a0 * np.sqrt(mu / epsilon)
pref3d = pref[:, np.newaxis, np.newaxis]
respo = pref3d * respo
# compute cross section of device
self.set_units("si")
dims = np.array([0, 1, 2])
dims = dims[dims != self.transport_axis] # dimensions of cross section
bc = np.array(self.center.system.cell.boundary)
bc = bc[2 * dims]
v0 = self.center.system.cell.avec[dims[0], :].magnitude
v1 = self.center.system.cell.avec[dims[1], :].magnitude
if bc[0] == 0 and bc[1] == 0:
unit = "ampere / m ** 2"
crossec = np.linalg.norm(np.cross(v0, v1)) / 1e20
elif bc[0] == 0 and bc[1] != 0:
unit = "ampere / m"
crossec = np.linalg.norm(v0) / 1e10
elif bc[0] != 0 and bc[1] == 0:
unit = "ampere / m"
crossec = np.linalg.norm(v1) / 1e10
else:
unit = "ampere"
crossec = 1.0
return to_quantity(respo / crossec, unit)
[docs]
def get_pho_curr(self, flux=1.0 * ureg.watt / ureg.meter**2, mu=1.0, epsilon=1.0):
"""shadow implementation for get_photo_current, should be removed in future"""
self.set_units("si")
self._read_photo_transmission()
bias = self.center.system.pop.bias
dims = np.array([0, 1, 2])
dims = dims[dims != self.transport_axis] # dimensions of cross section
bc = np.array(self.center.system.cell.boundary)
bc = bc[2 * dims]
v0 = self.center.system.cell.avec[dims[0], :].magnitude
v1 = self.center.system.cell.avec[dims[1], :].magnitude
if bc[0] == 0 and bc[1] == 0:
unit = "ampere / m ** 2"
crossec = np.linalg.norm(np.cross(v0, v1)) / 1e20
elif bc[0] == 0 and bc[1] != 0:
unit = "ampere / m"
crossec = np.linalg.norm(v0) / 1e10
elif bc[0] != 0 and bc[1] == 0:
unit = "ampere / m"
crossec = np.linalg.norm(v1) / 1e10
else:
unit = "ampere"
crossec = 1.0
# energy integration
e2h = 3.87404585e-5 # e ** 2/h
dE = self.elec_energy_resol.to("eV").magnitude
tr = self.photo_transmission
tr = np.sum(tr, axis=0) * (dE * e2h / crossec) # (photon, k, s, lead)
# k integration
kwght = np.array(self.center.system.kpoint.kwght)
k4d = kwght[np.newaxis, :, np.newaxis, np.newaxis]
tr = tr * k4d
tr = np.sum(tr, axis=1) # (photon, s, lead)
flux_unit = ureg.watt / ureg.meter**2
pref = get_photo_current_prefactor(mu, epsilon)
pref = (
pref
* flux.to(flux_unit).magnitude
/ self.photon_energy.to("joule").magnitude ** 2
)
pref3d = pref[:, np.newaxis, np.newaxis]
tr = pref3d * tr
if self.center.system.hamiltonian.ispin == 1:
tr = 2.0 * tr # spin degenerate
return to_quantity(tr, unit)
def _read_photo_transmission(self, filename="nano_photo_out.h5"):
"""get transmission result from HDF5 file"""
f = h5py.File(filename, mode="r")
fld = f["photo_transmission"][0:]
fld = np.transpose(fld, [i for i in range(fld.ndim - 1, -1, -1)])
fld = np.asfortranarray(fld)
self.photo_transmission = fld
def get_unit_vec_from_angles(e0, e1, theta, phi):
a0 = np.cos(theta) * np.cos(phi) - 1j * np.sin(theta) * np.sin(phi)
a1 = np.sin(theta) * np.cos(phi) + 1j * np.cos(theta) * np.sin(phi)
v = a0 * e0 + a1 * e1
return v.astype(np.cdouble) / np.linalg.norm(v)
def get_photo_current_prefactor(mu=1.0, epsilon=1.0):
me = 9.1093837015e-31 # Electron mass
e = 1.60217662e-19 # Elementary charge
h = 6.62607004e-34 # * ureg.meter**2 * ureg.kg / ureg.second # Planck's constant
hbar = h / 2 / np.pi # Reduced Planck's constant
eps0 = 8.8541878128e-12 # Vacuum permittivity
a0 = 5.29177210903e-11 # * ureg.meter # Bohr radius
alpha = 0.0072973525693 # fine structure constant
return np.sqrt(mu / epsilon) * alpha * h * a0 * a0