[ase-users] I-V curve from ase.transport

[Vincent Pohl]

Dear Mikkel,

thank you very much for your quick response!

Concerning the units of the method "get_current": How can I be in units of
2e/h (G_0)? Shouldn't it be in units of current, i.e. 2e^2/h*V?

I have updated the "get_current()" function as you suggested, also adding
the possibility of calculating the current for spin-paired and
spin-unpaired systems. Did I do this correctly? I have further updated the
method "fermidistribution" and will include that in the final commit.

With best regards,
Vincent

-------------------------------------

def get_current(bias, E, T_E, T = 0., spinpol=False):
    '''Returns the current as a function of the
    bias voltage.

    **Parameters:**
    bias : {float, (M,) ndarray}, units: V
      Specifies the bias voltage.
    E : {(N,) ndarray}, units: eV
      Contains energy grid of the transmission function.
    T_E {(N,) ndarray}, units: unitless
      Contains the transmission function.
    T : {float}, units: K, optional
      Specifies the temperature.
    spinpol: {bool}, optional
      Specifies wheter the current should be
      calculated assuming degenerate spins

    **Returns:**
    I : {float, (M,) ndarray}, units: 2e/h*eV
      Contains the electric current.
    '''
    if not isinstance(bias, float):
        bias = bias[np.newaxis]
        C = E[:, np.newaxis]
        T_E = T_E[:, np.newaxis]

    fl = fermidistribution(E - bias, units.kB * T)
    fr = fermidistribution(E + bias, units.kB * T)

    if spinpol:
        return 0.5*np.trapz((fl - fr) * T_E, x=E, axis=0)
    else:
        return np.trapz((fl - fr) * T_E, x=E, axis=0)

2016-11-07 10:20 GMT+01:00 Mikkel Strange <mikst at fysik.dtu.dk>:

> Dear Vincent,
>
> 1) the method "get_transmission" calculates the transmission (unit-less).
> In the zero-bias and zero temperature limit, the conductance is G = G_0 *
> T(E_F), where E_F is the Fermi-level.
>
> 2) At a quick glance your method "get_current" looks ok, but perhaps the
> current should just be returned
> in units of 2e/h?, i.e. without the factor "2. * units._e /
> units._hplanck".
> Also, one could perhaps try to use/extend the Fermi-distribution  all
> ready in tools.py?
>
> Best wishes,
> Mikkel
> ------------------------------
> *From:* ase-users-bounces at listserv.fysik.dtu.dk [
> ase-users-bounces at listserv.fysik.dtu.dk] on behalf of Vincent Pohl via
> ase-users [ase-users at listserv.fysik.dtu.dk]
> *Sent:* Friday, November 04, 2016 4:10 PM
> *To:* ase-users at listserv.fysik.dtu.dk
> *Subject:* [ase-users] I-V curve from ase.transport
>
> Dear ASE users,
>
> I want to compute I-V curves for an electron transport though a
> nanojunction. I have used the transport module of ASE to compute the
> transmission function. Afterwards, I wrote a function to evaluate the
> current as a function of the bias (See below). I am very new in the field,
> and I am not quiet sure, if I have understood everything correctly:
>
> 1) Am I right in assuming that the function "T = calc.get_transmission()"
> computes the quantum conductance (G(E) in units of G_0 = 2e^2/h) and not
> the transmission function which is unitless? If I am right, isn't the name
> of the function quiet misleading?
>
> 2) Concerning my function to compute the I-V curve (get_current, see
> below), can anyone please check its correctness? I am not quiet sure, if I
> have used the correct formula, and if I have handled the ase.units
> correctly.
>
> If everything is correct, shall I add the function to the ASE transport
> module and commit the changes. I think such a function might be very useful.
>
> Thank you very much in advance for your help.
>
> With best regards,
> Vincent
>
> ----------------------------------
>
> from __future__ import division
> from ase.transport.calculators import TransportCalculator
> from scipy import integrate
> import numpy as np
> from ase import units
> import pylab as plt
>
> h = np.array((0,)).reshape((1,1))
> h1 = np.array((0, -1, -1, 0)).reshape(2,2)
> energies = np.arange(-3, 3, 0.1)
> tcalc = TransportCalculator(h=h, h1=h1, energies=energies)
> G_E = tcalc.get_transmission() # Quantum conductance in units of G_0=2e^2/h
> T_E = 2 * G_E                  # Transmission function (Unitless)
>
> def get_current(bias, E, T_E, T = 0):
>     '''Returns the current as a function of the
>     bias voltage.
>
>     **Parameters:**
>     bias : {float, (M,) ndarray}, units: V
>       Specifies the bias voltage.
>     E : {(N,) ndarray}, units: eV
>       Contains energy grid of the transmission function.
>     T_E {(N,) ndarray}, units: unitless
>       Contains the transmission function.
>     T : {float}, units: K, optional
>       Specifies the temperature.
>
>     **Returns:**
>     I : {float, (M,) ndarray}, units: A
>       Contains the electric current.
>     '''
>     if not isinstance(bias, float):
>         bias = bias[np.newaxis]
>         E = E[:, np.newaxis]
>         T_E = T_E[:, np.newaxis]
>     assert T >= 0., 'Negative temperature given!'
>     if T==0:
>         fl = (E - bias / 2. <= 0).astype(int)
>         fr = (E + bias / 2. <= 0).astype(int)
>     else:
>         beta = 1. / (units.kB * T)
>         fl = 1. / (1. + np.exp(beta * (E - bias / 2.)))
>         fr = 1. / (1. + np.exp(beta * (E + bias / 2.)))
>     E = units._e * E
>     I = 2. * units._e / units._hplanck * integrate.simps((fl - fr) * T_E,
> x=E, axis=0)
>     return I
>
> bias = np.arange(0, 2, .1)
> plt.plot(bias, get_current(bias, energies, T_E, T = 0.))
> plt.xlabel('U [V]')
> plt.ylabel('I [A]')
> plt.show()
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20161108/45903b2f/attachment.html>