5
$\begingroup$

I have a matrix $M$ in Bloch basis (momentum space). Then I have to take Fourier transform of $M$ to go to real space. The size of this lattice in real space can be fixed to some value $N_{unit}$. Furthermore, I have to determine the eigenvectors of this real space matrix and study dynamics so it is very slow for bigger matrix sizes.

  1. Let us say $M$ is a square lattice with periodic boundaries can be reduced to a matrix hsquare (kx, ky) by taking Fourier transform, where $k_x$ and $k_y$ are the conjugate variables to $x$ and $y$.
  1. Then we take the inverse Fourier transform of hsquare (kx, ky) to go back to real space, i.e. matrix hsquareFin(x,y).
  2. Furthermore, I have to find the eigenvectors of hsquareFin(x,y). That I plot on top of underlying lattice which is square lattice here.
  3. Lastly, I am interested in the dynamics of how a pulse excited at a point in the real space will propagate according to the Schrödinger equation.

My MWE:

Nunit = 400; (* Number of unitcells *)

cond = {J1x -> 1/1, J2x -> 1, J1y -> 1/1};

hsquare[kx_, ky_] = (2 J1x Cos[kx] + 2 J1y Cos[ky] + x1) /. cond
     

square[kx_, ky_] = 2 Cos[kx] + 2 Cos[ky];


(*The Hamiltonian matrix for the Finite Case*)

DistributeDefinitions[Nunit, hsquare];
hsquareFinx[ky_] = 
  ParallelTable[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
       hsquare[kx, ky]/Sqrt[2 \[Pi]], kx, 
       x] /. {DiracDelta -> KroneckerDelta, x -> j - i, x1 -> i}) // 
   ArrayFlatten; 

DistributeDefinitions[hsquareFinx];
hsquareFin = 
  ParallelTable[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
       hsquareFinx[ky]/Sqrt[2 \[Pi]], ky, 
       y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) // 
   ArrayFlatten; 


graphsq[ky_] = 
  Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
       square[kx, ky]/Sqrt[2 \[Pi]], kx, 
       x] /. {DiracDelta -> KroneckerDelta, x -> j - i}) // 
   ArrayFlatten;

sqg = AdjacencyGraph[
   Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
        graphsq[ky]/Sqrt[2 \[Pi]], ky, 
        y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) // 
    ArrayFlatten, VertexLabels -> Automatic];  

The eigenvectors are calculated

graphcoordsq = GraphEmbedding[sqg]; (*Generating underlying lattice*)

egvecsq = Eigenvectors[hsquareFin // N];

Manipulate[
 Show[ListDensityPlot[
   Table[
    Append[Chop@graphcoordsq[[i]], Chop@(Abs[egvecsq[[n, i]]]^2)] // 
     N, {i, Length[egvecsq[[n]]]}], PlotLegends -> Automatic, 
   PlotLabel -> "Probability density"], sqg], {n, 1, Length[egvecsq], 
  1}]

Is there a way to speed up the process? Because I also have to study the time dynamics of the eigenstates later.

$\endgroup$
15
  • 3
    $\begingroup$ It seems as if your matrixes have many 0 entries. I have not waited for the computations to finish. But I think the resulting matrices will be banded with a very simple structure. You should seek to exploit that, e.g. by using SparseArray. $\endgroup$ Commented Jan 21, 2022 at 18:52
  • $\begingroup$ @Shamina It could be better to explain you problem as mathematical problem as well. $\endgroup$ Commented Jan 22, 2022 at 5:16
  • 3
    $\begingroup$ Yes. The Table in the function hsquareFin should not loop over all {i,j} pairs. Try to loop over those pairs that lead to a nonzero entry. Afterwards, you can use the SparseArray constructor to construct the sparse matrix from the nonzero positions and the nonzero values. $\endgroup$ Commented Jan 24, 2022 at 17:35
  • 1
    $\begingroup$ @CATrevillian Thanks a lot! Let me have a look into that. I don’t think I’ve fully got it but I am going to give it a try. $\endgroup$
    – Shamina
    Commented Jan 28, 2022 at 21:38
  • 2
    $\begingroup$ @Shamina Do you understand that we have to compute matrix $Nunit^2\times Nunit^2$, therefore for Nunit=400 it gonna be $160000\times 160000$. What is the reason to compute so big matrix? For Nunit<=100 there is algorithm faster then your code. $\endgroup$ Commented Jan 29, 2022 at 4:37

1 Answer 1

5
+50
$\begingroup$

This code shows how fast SparseArray compare to ParallelTable for moderate $N_{unit}\le 10^2$. Note, that during computation code generates $N_{unit}^2\times N_{unit}^2$ matrixes hsquareFin1=hsquareFin and solutions egvecsq1=egvecsq. For $N_{unit}>100$ the computation time is limited by time of Eigenvectors evaluation, and therefore we can't decrease computation time for Nunit=400 as it required. To compare computation time we can organize Module as follows

Clear["Global`*"]

cond = {J1x -> 1, J2x -> 1, J1y -> 1};

hsquare[kx_, ky_] = (2 J1x Cos[kx] + 2 J1y Cos[ky] + x1) /. cond;


square[kx_, ky_] = 2 Cos[kx] + 2 Cos[ky];


fcomp[nn_] := 
 Module[{Nunit = nn}, 
  t1 = AbsoluteTiming[DistributeDefinitions[Nunit, hsquare];
     hsquareFinx[ky_] = 
      ParallelTable[#, {i, Nunit}, {j, 
           Nunit}] &@(InverseFourierTransform[
           hsquare[kx, ky]/Sqrt[2 \[Pi]], kx, 
           x] /. {DiracDelta -> KroneckerDelta, x -> j - i, 
           x1 -> i}) // ArrayFlatten;
     
     DistributeDefinitions[hsquareFinx];
     hsquareFin = 
      ParallelTable[#, {i, Nunit}, {j, 
           Nunit}] &@(InverseFourierTransform[
           hsquareFinx[ky]/Sqrt[2 \[Pi]], ky, 
           y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) // 
       ArrayFlatten;
     graphsq[ky_] = 
      Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
           square[kx, ky]/Sqrt[2 \[Pi]], kx, 
           x] /. {DiracDelta -> KroneckerDelta, x -> j - i}) // 
       ArrayFlatten;
     
     sqg = 
      AdjacencyGraph[
       Table[#, {i, Nunit}, {j, Nunit}] &@(InverseFourierTransform[
            graphsq[ky]/Sqrt[2 \[Pi]], ky, 
            y] /. {DiracDelta -> KroneckerDelta, y -> j - i}) // 
        ArrayFlatten, VertexLabels -> Automatic];
     (*The eigenvectors are calculated*)
     
     graphcoordsq = 
      GraphEmbedding[sqg];(*Generating underlying lattice*)
     egvecsq = Eigenvectors[hsquareFin // N // Normal];][[1]];
  
  t2 = AbsoluteTiming[
     s1 = 
      With[{n = Nunit}, 
       SparseArray[{Band[{1, 1}] -> 
          Flatten[Table[Table[i, {i, n}], {n}]], 
         Band[{1, n + 1}] -> 1, Band[{n + 1, 1}] -> 1, 
         Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}], 
         Band[{2, 1}, {n^2, n^2}] -> 
          Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]];
     s2 = 
      With[{n = Nunit}, 
       SparseArray[{Band[{1, n + 1}] -> 1, Band[{n + 1, 1}] -> 1, 
         Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}], 
         Band[{2, 1}, {n^2, n^2}] -> 
          Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]]; 
     sqg1 = 
      AdjacencyGraph[s2 // ArrayFlatten, VertexLabels -> Automatic]; 
     graphcoordsq1 = 
      GraphEmbedding[sqg1];(*Generating underlying lattice*)
     hsquareFin1 = (s1 // ArrayFlatten); 
     egvecsq1 = Eigenvectors[hsquareFin1 // N // Normal];][[1]]; {nn, 
   t1, t2}]

Here t1 is code proposed by Shamina, and t2 is my code with hand made array. First we call

fcomp[100]

Out[]= {100, 115.673, 82.1941} 

We see that computation times t1, t2 not so differ due to large Eigenvectors evaluation time. We also can compare matrixes

{hsquareFin1 == hsquareFin, sqg1 == sqg, egvecsq1 == egvecsq}

Out[]= {True, True, True}

Now we can prepare list for range 5<=Nunit<=40 and plot it separately and together with point Nunit=100

lst = Table[fcomp[i], {i, 5, 40}]

Out[]= {{5, 0.164432, 0.0043509}, {6, 0.205826, 0.0050914}, {7, 
  0.220041, 0.0064627}, {8, 0.253057, 0.0075893}, {9, 0.295401, 
  0.0123119}, {10, 0.319457, 0.0125374}, {11, 0.364707, 
  0.0142791}, {12, 0.384723, 0.016998}, {13, 0.429124, 
  0.0211146}, {14, 0.493359, 0.0247619}, {15, 0.513203, 
  0.0269879}, {16, 0.593006, 0.0354736}, {17, 0.648195, 
  0.0353599}, {18, 0.685608, 0.0465485}, {19, 0.743417, 
  0.0658679}, {20, 0.844621, 0.0619918}, {21, 0.890146, 
  0.0794076}, {22, 0.962821, 0.0800831}, {23, 1.04226, 0.108473}, {24,
   1.15235, 0.111343}, {25, 1.2167, 0.129602}, {26, 1.32637, 
  0.158706}, {27, 1.41704, 0.162905}, {28, 1.50065, 0.182941}, {29, 
  1.67359, 0.218453}, {30, 1.73249, 0.219222}, {31, 1.92085, 
  0.252828}, {32, 2.05668, 0.287723}, {33, 2.26161, 0.313731}, {34, 
  2.35782, 0.356584}, {35, 2.57131, 0.425878}, {36, 2.74442, 
  0.454498}, {37, 2.99826, 0.524451}, {38, 3.26045, 0.578484}, {39, 
  3.5149, 0.66792}, {40, 3.75821, 0.66004}}

Figure 1

To visualize result we can use this code

fm[nn_] := 
 Module[{Nunit = nn}, 
  s1 = With[{n = Nunit}, 
    SparseArray[{Band[{1, 1}] -> 
       Flatten[Table[Table[i, {i, n}], {n}]], Band[{1, n + 1}] -> 1, 
      Band[{n + 1, 1}] -> 1, 
      Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}], 
      Band[{2, 1}, {n^2, n^2}] -> 
       Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]];
  s2 = With[{n = Nunit}, 
    SparseArray[{Band[{1, n + 1}] -> 1, Band[{n + 1, 1}] -> 1, 
      Band[{1, 2}, {n^2, n^2}] -> Table[If[i == n, 0, 1], {i, n}], 
      Band[{2, 1}, {n^2, n^2}] -> 
       Table[If[i == n, 0, 1], {i, n}]}, {n^2, n^2}]]; 
  sqg1 = AdjacencyGraph[s2 // ArrayFlatten, 
    VertexLabels -> Automatic]; graphcoordsq1 = GraphEmbedding[sqg1]; 
  hsquareFin1 = (s1 // ArrayFlatten); 
  egvecsq1 = Eigenvectors[hsquareFin1 // N // Normal]; {egvecsq1, 
   graphcoordsq1}]

With[{nn = 12}, fm[nn]; 
 Manipulate[
  Show[ListDensityPlot[
    Table[Append[Chop@graphcoordsq1[[i]], 
       Chop@(Abs[egvecsq1[[j, i]]]^2)] // N, {i, 
      Length[egvecsq1[[j]]]}], PlotLegends -> Automatic, 
    Frame -> False, PlotRange -> {{0, nn - 1/2}, {0, nn - .5}, All}, 
    ColorFunction -> Hue], sqg1], {j, 1, Length[egvecsq1], 1}, 
  FrameLabel -> {"", "", "Probability density"}, ImageMargins -> 0, 
  FrameMargins -> Small, Frame -> False]]

Figure 2

$\endgroup$

Your Answer

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

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