Why is ListPointPlot3D extremely slow when plotting many points with per-point Style colors, and how can this be optimized?

I am generating a 3D band structure from a slab Hamiltonian and color-coding each eigenvalue by its IPR (defined from the corresponding eignvector, e.g., inverse participation ratio). The numerical diagonalization and IPR calculation are fast, but ListPointPlot3D becomes extremely slow and almost unresponsive when the number of k-points increases.

A minimal version of my code is below.

1. Define the matrix:

hamSlabXZ[{kx_, kz_}, t3_] := 
 Module[{ny = 5, H0, HL, HR, h00, h01, h10, hL00, hL01, hR00, hR10, 
   HfullCell},
  h00 = ( {
     {0, ta, ta},
     {ta, 0, ta},
     {ta, ta, 0}
    } );
  h01 = ( {
     {0, 0, 0},
     {tz, 0, 0},
     {tb, tz, 0}
    } ); h10 = ( {
     {0, tz, tb},
     {0, 0, tz},
     {0, 0, 0}
    } );
  H0 = KroneckerProduct[IdentityMatrix[ny], h00]
    + KroneckerProduct[DiagonalMatrix[Array[1 &, ny - 1], 1], h01]
    + KroneckerProduct[DiagonalMatrix[Array[1 &, ny - 1], -1], h10];
  hL00 = ( {
      {0, tb, tz},
      {0, 0, 0},
      {0, tz, 0}
     } )*Exp[-I*kx];
  hL01 = ( {
      {0, tz, 0},
      {0, 0, 0},
      {tz, tb, 0}
     } )*Exp[-I*kx];
  HL = KroneckerProduct[IdentityMatrix[ny], hL00]
    + KroneckerProduct[DiagonalMatrix[Array[1 &, ny - 1], 1], hL01];
  hR00 = ( {
      {0, 0, 0},
      {tb, 0, tz},
      {tz, 0, 0}
     } )*Exp[I*kx];
  hR10 = ( {
      {0, 0, tz},
      {tz, 0, tb},
      {0, 0, 0}
     } )*Exp[I*kx];
  HR = KroneckerProduct[IdentityMatrix[ny], hR00]
    + KroneckerProduct[DiagonalMatrix[Array[1 &, ny - 1], -1], hR10];
  HfullCell = H0 + HL + HR;
  
  HfullCell /. {ta -> -1, tb -> -2.5, tz -> 2*t3*Cos[kz]}
  ]

2. Get the eigenvalues in (kx,kz) plane, together with the corresponding IPR by specifying Style, and storing them with data2DAll

surfband2D[kx_, kz_] := 
 Module[{ham, evals, evecs, ord, evalsSort, evecsSort, iprList, 
   iprNormal, colorlist},
  ham = N@hamSlabXZ[{kx, kz}, -0.02];
  {evals, evecs} = Eigensystem[ham];
  ord = Ordering[evals];
  evalsSort = evals[[ord]];
  evecsSort = evecs[[ord]];
  iprList = Map[Total[(Abs[#]^4)] &, evecsSort];
  iprNormal = iprList/Max[iprList];
  colorlist = 
   Map[Blend[{Lighter[Blue, 2/3], Lighter[Gray, 1/3], Orange, 
       Red}, #] &, iprNormal
    ];
  MapThread[Style, {evalsSort, colorlist}]
  ]
data2DAll = 
  ParallelTable[
   surfband2D[kkx, kkz], {kkx, -\[Pi], \[Pi], 
    0.05 \[Pi]}, {kkz, -\[Pi], \[Pi], 0.05 \[Pi]}];

3. ListPointPlot it:

ListPointPlot3D[Transpose[data2DAll[[All, All, All]], {3, 2, 1}], 
 AxesLabel -> {kx, kz, Ene}, 
 BoxRatios -> {1, 1, 1.5}]

The result may look like the picture blow, but is getting significantly slow with the increase of (kx,kz) mesh grids.

Questions

1. Is the use of per-point Style inside ListPointPlot3D the main cause of the slowdown?

2. What is the recommended way to efficiently color a large number of 3D points by an additional scalar quantity (such as IPR)?

3. A similar issue is found in here, where the top answer adopts a solution by replacing ListPlot with Graphics, just like the answer made by @Craig Carter below this question. So, a further question arises: Is that possible to plot the entire surface via ListPlot3D rather than play with point-dots?

Any insight into the internal rendering behavior of ListPointPlot3D or suggested alternative approaches would be appreciated.