Python Model IV

Question

Model IV.

Method: Run the integral as a function of E, which outputs the current for each voltage value used. This is repeated for an array of v_values ​​values. The equation can be found below.

Although the limits of this equation range from -inf to inf , the limits should be limited, so (E + eV) ^ 2- \ Delta ^ 2> 0 and E ^ 2- \ Delta ^ 2> 0 to avoid poles. (\ Delta_1 = \ Delta_2). Therefore, there are currently two integrals with limits from -inf to -gap-e*v and gap to inf .

However, I keep returning math range error , although, in my opinion, I have eliminated the unpleasant values ​​of E using the restrictions mentioned above. Mouth bugs: http://pastie.org/private/o3ugxtxai8zbktyxtxuvg

Apologies for the vagueness of this matter. But can anyone see obvious errors or misuse of code?

My attempt:

 from scipy import integrate from numpy import * import scipy as sp import pylab as pl import numpy as np import math e = 1.60217646*10**(-19) r = 3000 gap = 400*10**(-6)*e g = (gap)**2 t = 0.02 k = 1.3806503*10**(-23) kt = k*t v_values = np.arange(0,0.001,0.0001) I=[] for v in v_values: val, err = integrate.quad(lambda E:(1/(e*r))*(abs(E)/np.sqrt(abs(E**2-g)))*(abs(E+e*v)/(np.sqrt(abs((E+e*v)**2-g))))*((1/(1+math.exp((E+e*v)/kt)))-(1/(1+math.exp(E/k*t)))),-inf,(-gap-e*v)*0.9) I.append(val) I = array(I) I2=[] for v in v_values: val2, err = integrate.quad(lambda E:(1/(e*r))*(abs(E)/np.sqrt(abs(E**2-g)))*(abs(E+e*v)/(np.sqrt(abs((E+e*v)**2-g))))*((1/(1+math.exp((E+e*v)/kt)))-(1/(1+math.exp(E/k*t)))),gap*0.9,inf) I2.append(val2) I2 = array(I2) I[np.isnan(I)] = 0 I[np.isnan(I2)] = 0 pl.plot(v_values,I,'-b',v_values,I2,'-b') pl.show() 

Answer 1

The best way to approach this problem in the end is to use a heavy function to prevent the E variable from exceeding the \Delta variable.

Answer 2

This question is best suited for the Computational Science site. However, there are a few points you should think about.

First, the integration range is the intersection of (-oo, -eV-gap) U (-eV+gap, +oo) and (-oo, -gap) U (gap, +oo) . Two cases are possible:

• if eV < 2*gap , then the permissible energy values ​​are in (-oo, -eV-gap) U (gap, +oo) ;
• if eV > 2*gap , then the permissible energy values ​​are in (-oo, -eV-gap) U (-eV+gap, -gap) U (gap, +oo) .

Secondly, you work in an area with very low temperatures. At t equal to 0.02 K, the denominator in the Boltzmann factor is 1.7 meV, and the energy gap is 400 meV. In this case, the exponent is huge for positive energies and soon goes beyond the double-precision floating-point numbers used by Python. Since this is the minimum positive energy possible, things will not improve at higher energies. At negative energies, the value will always be very close to zero. Note that at this temperature, the Fermi-Dirac distribution has a very sharp edge and resembles the reflected theta function. At E = gap you will have exp(E/kT) approximately 6.24E + 100. You can go out of range if E/kT > 709.78 or E > 3.06*gap .

Nevertheless, it makes no sense to go to such energies, since at this temperature the difference between the two Fermi functions very quickly becomes zero outside the interval [-eV, 0] , which completely falls inside the gap for a given temperature, when V < (2*gap)/e (0.8 mV). Therefore, one would expect that the current would be very close to zero when the bias voltage is less than 0.8 mV. When it is greater than 0.8 mV, then the main value of the integral would be obtained from the integrand in (-eV+gap, -gap) , although some nonzero value would come from the region near the singularity at the point E = gap , and some singularity at E = -eV-gap . You should not avoid features in DoS, otherwise you will not get the expected gaps (vertical lines) on curve I (V) (image taken from Wikipedia ):

Rather, you should get equivalent approximate expressions in the neighborhood of each feature and integrate them instead.

As you can see, there are many special cases for the value of the integrand, and you must consider all of them when calculating numerically. If you do not want to do this, you should probably refer to another mathematical package, such as Maple or Mathematica. They have much more complex numerical integration procedures and can directly process your formula.

Please note that this is not an attempt to answer your question, but rather a very long comment that does not fit in any comment field.

Answer 3

The reason for the math range error is that your exponentiality goes to infinity. Taking v = 0.0009 and E = 5.18e-23 , the expression exp((E + e*v) / kt) (I corrected the typo noted by Christ Liev in your Python expression) is exp(709.984..) , which goes beyond limits of the range that you can represent with double precision numbers (up to approx. 1E308).

Two additional notes:

• As others have already noted, you should probably re-scale your equation using a single system that delivers numbers in a smaller range. Perhaps atomic units are a possible choice, since it set e = 1 , but I did not try to convert your equation into it. (Probably your timestep will then become quite large, since in atomic units the time unit is about 1/40 fs).

• Usually, exponential notation is used for floating point numbers: e = 1.60217E-19 instead of e = 1.60217*10**(-19) .