# Plotting a Double Sum

I am attempting to plot a graph of a function that is the absolute value of the double sum of an exponential function, and I keep getting errors. Here is the code:

Ly = 10;
Lz = 12;
hbar = UnitConvert[Quantity[1, "\[HBar]"], "SIBase"];
h = UnitConvert[Quantity[1, "PlanckConstant"], "SIBase"];
mass = 9.11*10^(-31);
mstar = 0.1*mass;
q = 1.602*10^(-19);
u = 0.05;
V = 0.02;
k = 1.38*10^(-23);
T = 4;
Plot[((2*0.02*q^2)/h)*
Abs[Sum[(1 + E^((q*V - u + ((hbar^2)/(2*
mstar))*((π*m/Ly)^2 + (π*n/Lz)^2) - q*Vg)/(k*
T)))^(-1), {m, 1, Infinity}, {n, 1, Infinity}]], {Vg, 0, 7}]


I know it is a huge function. I have tried to use NSum, I have tried to use Sum, I have tried breaking up the function into pieces and assigning the pieces to variables to make it more manageable. But I keep getting errors:

NSum::nsnum: "Summand (or its derivative) 1/(1. +2.71828^((2.5*10^22 (-0.05+(<<1>>) <<1>>))/(1.38&))) is not numerical at point n = 46662."


In an attempt to isolate what might be going wrong, I tried to only sum over one of the indices by deleting the term with pi*n/Lz (and I did remember to delete

, {n, 1, Infinity}


as well). But I still got an error:

NSum::itraw: Raw object 1 cannot be used as an iterator. >>


This error was odd, since my iterator is specified as m, not a number.

I know from seeing the plot in textbooks and from entering my equation into Wolfram Alpha, given some value for Vg, that at the very least, a single sum WILL converge (W.A. couldn't do the double sum without running out of computational time), and that I should be able to plot this. Can anyone please help?

UPDATE:

Per the suggestions below, I have assigned numerical values to h and hbar rather than using UnitConvert, and I have switched from SI units to atomic units to eliminate the underflow problem. Now I am getting a new error when I try to plot the function:

SequenceLimit::seqlim: The general form of the sequence could not be determined, and the result may be incorrect.
Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[..., _SystemException].
SystemException["MemoryAllocationFailure"]


I tried to evaluate the sum at Vg = 2, using NSum since Sum was just giving me the symbolic form:

NSum[(1 +
E^((q*V -
u + ((hbar^2)/(2*mstar))*((\[Pi]*m/Ly)^2 + (\[Pi]*n/Lz)^2) -
q*2)/(k*T)))^(-1), {m, 1, Infinity}, {n, 1, Infinity}]


I got the same error about the "general form of the sequence could not be determined" as above, along with:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::ncvb: "NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in m near {m} = {29.6319}. NIntegrate obtained 1809.6488576380852 and 0.004501440304162219 for the integral and error estimates."
NSum::nsnum: "Summand (or its derivative) -((0.000690943\2.71828^(0.25 (-2.0011+5. (<<1>>)))\m)/(1. +2.71828^(0.25 (-<<19>>+5. <<1>>)))^2) is not numerical at point n = 16.."


Any additional help would be greatly appreciated.

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Feb 3 '16 at 9:10

TL;DR - Evaluate the function for just one value before trying to plot it. This will always save you hassle and headaches. And if your function involves a summation out to infinity, make sure it converges. Use atomic units when working on atomic-scale problems.

I personally stay away from the Units functionality altogether, and just enter the numbers in whatever unit system I'm working in (SI maybe, but more likely atomic units). To see why, let's evaluate the sum (but only keeping the first term) and substitute a numerical value for Vg

((2*0.02*q^2)/h)*
Abs[(1 + E^((q*V -
u + ((hbar^2)/(2*
mstar))*((π*m/Ly)^2 + (π*n/Lz)^2) - q*Vg)/(k*
T)))^(-1)] /. {Vg -> 1, m -> 14, n -> 1}


Right away you see an issue - you have units in the exponent, and units in the exponent always have to cancel. You can either assign units to all the other quantities, or just define h and hbar by their numeric values,

h = 6.626070*^-34;
hbar = h/(2 π);


Now that you have that done, try your calculation again,

which means that

The result is smaller than the value of \$MinNumber on this computer

So, now I reckon you ought to switch to atomic units since you are clearly working with electrons. This allows you to redefine

h = 2 π;
hbar = h/(2 π);
mass = 1.0;
mstar = 0.1*mass;
q = 1.0;
k = 1.0;


You still may need to convert some of the other constants, like Ly, Lz, V, etc.

• Fine TL;DR, +1. – user9660 Feb 3 '16 at 9:20
• Why thank you :-) – Jason B. Feb 3 '16 at 9:21
• Brilliant, thank you! I haven't converted it to atomic units yet, but I've already changed how I define h and hbar. I'm not getting a plot yet, but it is infinitely faster at giving me the empty x- and y- axes, and there are no error messages. I'll post an update once I try the unit conversion! – KFahyFoxtrot Feb 3 '16 at 9:29
• @KFahyFoxtrot - Definitely before you plot as a function of Vg, make sure that you get a number for just a single value of Vg`. Also, where do these equations come from, what are they trying to simulate? – Jason B. Feb 3 '16 at 9:30
• We are using the Landauer-Büttiker formalism to find the conductance of a 1D ballistic conduction channel as the applied gate voltage (Vg) varies. Quantum mechanics stuff. – KFahyFoxtrot Feb 3 '16 at 9:40