I have problem making ND work. Define relevant needs and vectors

d2 = 1./2. {-Sqrt[3.], 1.};
d3 = {0., -1.};
a1 = Sqrt[3.] {1., 0.};
a2 = Sqrt[3.]/2. {1., Sqrt[3.]};
a3 = -a1 + a2;
k = {kx, ky};

Define matrices

H[kx_, ky_, \[Alpha]_, \[Beta]_, \[Gamma]_] = {{0., 
2*\[Gamma]*(Sin[k.a1] + Sin[k.a2] + Sin[k.a3]), \[Alpha]*(1 + 
E^(I*k.a1) + E^(I*k.a2)), 0.}, {2*\[Gamma]*(Sin[k.a1] + 
Sin[k.a2] + Sin[k.a3]), 0., 
0., \[Beta]*(1 + E^(I*k.a1) + E^(I*k.a2))}, {\[Alpha]*(1 + 
   E^(-I*k.a1) + E^(-I*k.a2)), 0., 0., 
2*\[Gamma]*(Sin[k.a1] + Sin[k.a2] + Sin[k.a3])}, {0., \[Beta]* 
(1 +E^(-I*k.a1) + E^(-I*k.a2)), 
2*\[Gamma]*(Sin[k.a1] + Sin[k.a2] + Sin[k.a3]), 0.}};
Htbg[kx_, ky_] = 
2*(H[kx, ky, 1., -0.2, 0.2] + 0.3456 IdentityMatrix[4]);

Solve for normalized eigenstates. Transpose makes the column elements the eigenvalues. Sorting is based on the first column (the eigenvalues).

\[Psi][kx_?NumericQ, ky_?NumericQ] := \[Psi][kx, ky] = 
 Orthogonalize[Sort[Transpose[Eigensystem[Htbg[kx, ky]]]] 
[[All, 2]],Conjugate[#1].#2 &];

This gives wildly fluctuating phase values in k space which is not suitable for taking derivatives. Impose smooth gauge by making the 1st element of each eigenvectors real. That is, by multiplying a phase.

\[Psi]sg[kx_, ky_] := Table[Exp[-I*
  ArcTan[Re[\[Psi][kx, ky][[i, 1]]], 
   Im[\[Psi][kx, ky][[i, 1]]]]] \[Psi][kx, ky][[i]], {i, 1, 

As example, consider the [[1,1]] element. This gives a smooth plot

Plot[Re[\[Psi]sg[x, 1.][[1, 1]]], {x, -1, 1}]

But numerical derivative fails

ND[Re[\[Psi]sg[kx, 1.][[1, 1]]] // Chop, kx, 0.]

I do not why this is happening and would appreciate if someone can help. Thank you.

P.S. I tried to composed a simpler example with stripped-down irrelevant details, but the resulting simple example worked just fine. I do not know why the above does not work.


1 Answer 1



[Psi]sg[kx_?NumericQ, ky_?NumericQ] := ...

will reduce the number of error messages you get, although not completely eliminate them.

As you mentioned, his works:

Plot[Re[\[Psi]sg[x, 1.][[1, 1]]], {x, -1, 1}]

If you now say:

ND[Re[\[Psi]sg[kx, 1.][[1, 1]]] // Chop, kx, 0.]

You get a bunch of error messages of the form:

Part::partd: Part specification \[Psi]sg[kx,1.][[1,1]] is longer than depth of object.

Part::partw: Part 3 of \[Psi][kx,1.] does not exist.

This is because MMA evaluates arguments before feeding them to functions. In your case it evaluates: "[Psi]sg[kx, 1.][[1, 1]]", feeding "[Psi]sg[kx, 1.]" a symbolic argument, what makes it to fail.

If you define "[Psi]sg[kx_?NumericQ, ky_?NumericQ]=.." this will not evaluate and "[Psi]sg[kx,1.][[1,1]]" gives you only one error message, that it can not take part "[[1,1]]". It nevertheless then feeds the argument to "Re" and finally "ND" and gives a result.

This will not happen inside "Plot", because "Plot" is special, it does not evaluate its argument before giving numerical values to the variables. Therefore, we may plot the derivative without problems provided you define "[Psi]sg[kx_?NumericQ, ky_?NumericQ]=..":

Plot[ND[Re[\[Psi]sg[kx, 1.][[1, 1]]] // Chop, kx, kx0], {kx0, -1, 1}]

![enter image description here

  • $\begingroup$ This answer is very helpful and I also now have a clearer understanding of how ND works. Thank you so much! $\endgroup$
    – khvillegas
    May 13, 2021 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.