I kinda figured you'd need more info. But I didn't want to spam everyone right away! Here is everything I have with annotations. In desperation I have also attached the file. Any tips at all would be appreciated.

P.S. I should say that the end goal of this will be to plot that value of T[n,a,d,kx,l,g] {k,-10,10} without it taking infinite time.

Remove["Global`*"];

$Assumptions = a > 0 && d > 0 && l > 0 && g > 0 && n \[Element] Integers;

Hamiltonian for region with SOI and Zeeman field:

Hsoz[qx_, l_, g_] = {{-g, 0, px + I py, -2 I l}, {0, g, 0, px + I py}, {px - I py, 0, -g, 
  0}, {2 I l, px - I py, 0, g}}

Evals[qx_, l_, g_] = Simplify[Eigenvalues[Hsoz[qx, l, g]]];

Ez1[qx_, l_, g_] = Evals[qx, l, g][[1]];
Ez2[qx_, l_, g_] = Evals[qx, l, g][[2]];
Ez3[qx_, l_, g_] = Evals[qx, l, g][[3]];
Ez4[qx_, l_, g_] = Evals[qx, l, g][[4]];

Evecs[qx_, l_, g_] = Simplify[Eigenvectors[Hsoz[qx, l, g]]];

Normalization below (I suspect there is a better way to do this)

n1[qx_, l_, g_] = g^2 + l^2 + Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];
n2[qx_, l_, g_] = g^2 + l^2 - Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];

Nz1[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez1[qx, l, g] + g)]^2)/
    Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2 + 
    Abs[Ez1[qx, l, g]*(n1[qx, l, g] - g^2) + 
       g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz2[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez2[qx, l, g] + g)]^2)/
    Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2 + 
    Abs[Ez2[qx, l, g]*(n1[qx, l, g] - g^2) + 
       g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz3[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez3[qx, l, g] + g)]^2)/
    Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2 + 
    Abs[Ez3[qx, l, g]*(n2[qx, l, g] - g^2) + 
       g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];
Nz4[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez4[qx, l, g] + g)]^2)/
    Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2 + 
    Abs[Ez4[qx, l, g]*(n2[qx, l, g] - g^2) + 
       g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];

wz1[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz1[qx, l, g]]) Evecs[qx, l, g][[1]]];
wz2[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz2[qx, l, g]]) Evecs[qx, l, g][[2]]];
wz3[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz3[qx, l, g]]) Evecs[qx, l, g][[3]]];
wz4[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz4[qx, l, g]]) Evecs[qx, l, g][[4]]];

Conservation of energy:

qz3x[kx_, l_, g_] = Sqrt[
  g^2 + kx^2 - 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];
qz4x[kx_, l_, g_] = Sqrt[
  g^2 + kx^2 + 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];

Generic wave equation constructed:

M1[x_, kx_, l_, g_] = 
  Simplify[wz3[qz3x[kx, l, g], l, g]*E^(I qz3x[kx, l, g] x)];
M2[x_, kx_, l_, g_] = 
  Simplify[wz4[qz4x[kx, l, g], l, g]*E^(I qz4x[kx, l, g] x)];
M3[x_, kx_, l_, g_] = 
  Simplify[wz3[-qz3x[kx, l, g], l, g]*E^(-I qz3x[kx, l, g] x)];
M4[x_, kx_, l_, g_] = 
  Simplify[wz4[-qz4x[kx, l, g], l, g]*E^(-I qz4x[kx, l, g] x)];

And in matrix form:

Mz[x_, kx_, l_, g_] = ({
    {M1[x, kx, l, g][[1]], M2[x, kx, l, g][[1]], M3[x, kx, l, g][[1]],
      M4[x, kx, l, g][[1]]},
    {M1[x, kx, l, g][[2]], M2[x, kx, l, g][[2]], M3[x, kx, l, g][[2]],
      M4[x, kx, l, g][[2]]},
    {M1[x, kx, l, g][[3]], M2[x, kx, l, g][[3]], M3[x, kx, l, g][[3]],
      M4[x, kx, l, g][[3]]},
    {M1[x, kx, l, g][[4]], M2[x, kx, l, g][[4]], M3[x, kx, l, g][[4]],
      M4[x, kx, l, g][[4]]}
   });

Hamiltonian for region without SOI and Zeeman field:

L[x_, kx_] = 1./Sqrt[2.] ({{E^(I kx x), 0, -E^(-I kx x), 0},
{0, E^(I kx x), 0, -E^(-I kx x)},
{E^(I kx x), 0, E^(-I kx x), 0},
{0, E^(I kx x), 0, E^(-I kx x)}});

Generalising for potential barriers of width "d" every "a" units of distance

Sl[a_, kx_] = Simplify[L[2 a + d, kx].Inverse[L[a + d, kx]]];

Where it starts to really slow down

Smz[a_, d_, kx_, l_, g_] = Mz[a + d, kx, l, g].Inverse[Mz[a, kx, l, g]];

Sz[a_, d_, kx_, l_, g_] = Sl[a, kx].Smz[a, d, kx, l, g];

Impossible to compute beyond here. Even takes forever if I define every variable:

Anz[n_, a_, d_, kx_, l_, g_] := Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

Tferz[n_, a_, d_, kx_, l_, g_] := ({
   {1, 0, -Anz[n, a, d, kx, l, g][[1, 3]], -Anz[n, a, d, kx, l, g][[1,4]]},
   {0, 1, -Anz[n, a, d, kx, l, g][[2, 3]], -Anz[n, a, d, kx, l, g][[2,4]]},
   {0, 0, -Anz[n, a, d, kx, l, g][[3, 3]], -Anz[n, a, d, kx, l, g][[3,4]]},
   {0, 0, -Anz[n, a, d, kx, l, g][[4, 3]], -Anz[n, a, d, kx, l, g][[4,4]]}})

Auz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1,1]]}, {Anz[n, a, d, kx, l, g][[2,1]]}, {Anz[n, a, d, kx, l, g][[3,1]]}, {Anz[n, a, d, kx, l, g][[4,1]]}}

Tuz[n_, a_, d_, kx_, l_, g_] := Abs[Auz[n, a, d, kx, l, g][[1]]]^2 + Abs[Auz[n, a, d, kx, l, g][[2]]]^2

Ruz[n_, a_, d_, kx_, l_, g_] := 1 - Tuz[n, a, d, kx, l, g]

Adz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1, 2]]}, {Anz[n, a, d, kx, l, g][[2, 2]]}, {Anz[n, a, d, kx, l, g][[3, 2]]}, {Anz[n, a, d, kx, l, g][[4, 2]]}}

Tdz[n_, a_, d_, kx_, l_, g_] := Abs[Adz[n, a, d, kx, l, g][[1]]]^2 + Abs[Adz[n, a, d, kx, l, g][[2]]]^2
Rdz[n_, a_, d_, kx_, l_, g_] := 1 - Tdz[n, a, d, kx, l, g]

Tz[n_, a_, d_, kx_, l_, g_] := (Tuz[n, a, d, kx, l, g] + Tdz[n, a, d, kx, l, g])/2.
Rz[n_, a_, d_, kx_, l_, g_] := 1 - Tz[n, a, d, kx, l, g]

Attachments:

In the end I want to make plots of T[n,a,d,kx,l,g] {k,-10,10} for various values of the other variables (and eventually want to include an angular dependence but I'm a ways from that yet)

I have played around with := for the parts that take a long time to evaluate, but when I try to make a graph the kernel eventually times out and mathematica shuts down.

I don't quite know what you mean by "numerically" and "numerically directly"

Thank you for responding.

To add, and please put the numbers with a period in it: 3.0 (not 3) et cetera. So you will trigger the usage of machine numbers rather then exact arithmetic.

and the function T itself would also be useful ;)

DeleteCases does what is kind-a says: it delete certain cases, in this case 0.0 is removed from the list, as Tz blows up for k=0. No i meant that 'Tz' returns a vector (of 1 value) not just a value. is that correct?

I don't see a difference!

I sped it up a bit more (factor of 6-7). A single value now takes 15 milliseconds on my computer, which is good improvement compared to your 30 minutes (a factor 10^5)

Please test every subpart critically, and make sure everything is as it should be. I suspect there is at least one more error somewhere as it returns {value} rather than value, but without knowing what to expect, it is hard to guess what you actually want to do...

To further speed it up, you want to write the M-function as a single function that return 4 values, so you can turn wz* in a single function, and als Nz and Ez*, so the eigenvalues/eigenvectors are only calculated once, not 4x (8? 16?) times. Furthermore you could use the function Eigensystem to get the eigenvalues and eigenvectors at the same time, much faster (2x speed up probably).

good luck!

Attachments:

Once you have some sensible answers from all your sub-functions and you need a big speed up, just continue this thread. I will be notified and others will also see it. I guess it can done in a millisecond or so, and when parallelised even faster of course...

Yes, please use Normalize as it also deals with the case {0,0,0} which (when not programmed properly) blows up if you just divide by the length of the vector...

I am far from an expert, but I know that as a rookie one can be far from optimal. On the other hand, it could also be that the problem you are trying to do is just too hard (it certainly looks that way). I will look at your notebooks some more. I was hoping to see something besides code. I am not very familiar with condensed matter.

Best,

OL.

Typical values I'm trying now

n=5. (this really should be exactly 5 but if this runs faster 5.0 will do)

a=0.23

d=0.12

kx= 1. (this is what I want to plot over so 1. is fine for now)

l=0.01

g=0.09

kx! Yes. I caught that mistake as I was sending you the notebook but I didn't change it.

I suppose it should be a 1D vector cos that is for 20 values of k (?)

What is the DeleteCases command?

Other optimizations would be great if you have time. I am such a novice at this. I'm literally transposing a handwritten work-through of this problem.

Ooh. That will plot Tz.

Unfortunately I just noticed an error in the notebook I sent you

Anz[n_, a_, d_, kx_, l_, g_] := 
 Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

should read

Anz[n_, a_, d_, kx_, l_, g_] := 
 Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

(I had tried to make a simplification that didn't work.)

Changing that the Table[] still evaluates very quickly. And the graph of Tz comes out in 60 seconds.

It is, however not between 1 and 0 which is troubling. I may have more errors than I thought.

Have a look at the following functions:

AbsArg
AngleVector
ReIm
ArcTan[x, y]

The first gives the magnitude and the 'angle' (argument) of a complex number at once. The second can construct a vector by a radius-angle pair The third gives you real and Imaginary part of a number as a pair. And perhaps use Arctan with two arguments; so it takes care of the 4 quadrants (the signs of x and y et cetera).

The problem is in your eigensystem function:

First symbolic finding the general eigenvalues and eigenvectors may involve certain assumptions, not valid for all cases! Take for example the simple case of this polynomial:

In the general case 3x^2/x is 3x but not for x=0 ! So 'underneath' the symbolic calculation there are always some assumptions (unavoidable in some sense). I'm not sure why the general solution of the eigensystem 'works' for your input but not the non-general case!

I'm not quite sure, but I guess generally it has to do with cases there the numbers become '0' and you have to do some division. So is it possible for your elements in your matrix to be zero? Just check what you expect, you can put some values in there and try. I put 3 random values for qx qy and l in there...

Also you can check the symbolic result, is that what you want?

<<Progress Update>>

Now I have successfully expanded this for two barriers (or a square potential well) and am attempting to expand fro arbitrarily many barriers. Getting errors I don't understand though. Also still can't delay evaluation without getting the wrong result.

Attachments:

Delayed evaluation works for me, but only if you delay everything but the function that rely on Eigensystem commands...

Thanks for replying. I don't know how far down the thread you have read.

It's the basic square potential well problem. My tweaks are that I want a generic wave including spin, arbitrarily many identical barriers uniformly spaced ("a" apart and "d" wide) in graphene, with spin orbit interaction and a zeeman field under the potential.

My approach is to first try to replicate some simplified systems (only spin orbit, only one barrier, normal incidence etc) and then to generalise entirely so I can have any number of barriers at any incidence angle and with both SOI and zeeman.

I understand the maths of it (for the most part) but I am having huge trouble translating it to mathematica language cos I am far from fluent. I am making all of the rookie mistakes and it's costing me computing time, memory, and accuracy!