Plotting a Double Sum

Question

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.

Answer 1

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.