From a first principles bandstructure calculation I get an energy scalar field in three dimensions $E(x,y,z)$. It's now easy to plot a constant energy (contour)-surface for dedicated values $\epsilon=E(x,y,z)$ by using ListContourPlot3D. Lets assume the following example scalar field (e.g. the case of a graviatation potential field):

data = Table[1/Sqrt[x^2 + y^2 + z^2], {x, -1, 1, 0.1}, {y, -1, 1, 0.1}, {z, -1, 1,0.1}]
        /. ComplexInfinity -> 10^16 // Quiet;

I exchanged the infinite value at $x=y=z=0$ with a finite one to prevent any error messages coming from the visualisation routines. Now the iso-surfaces can be easily plotted for e.g. $\epsilon=2.0$

ListContourPlot3D[data, Contours -> {2.0}, DataRange -> {{-1, 1}, {-1, 1}, {-1, 1}}]

Isosurface for $\epsilon$=2.0

Coming back to the bandstructure calculation, the energy surfaces are much more complicated and cannot be described analytically for the general case as can be seen from the following plot of an energy surface for Silicon:

Isoenergy surface for Silicon

An important electronic property, the DOS (Density of States) is calculated through the integral $$\oint \limits_{E(x,y,z)=\epsilon}\frac{dS}{\lvert\nabla E(x,y,z)\rvert}.$$

and additional integrals of interest have the general form

$$\oint \limits_{E(x,y,z)=\epsilon}\frac{dS}{\lvert\nabla(\Delta E(x,y,z))\rvert}A(x,y,z)$$

with $A(x,y,z)$ as some additional scalar function (e.g. for effective mass tensor, $A(x,y,z)$ is the tensor product of the gradients of the energy).

Now for the question. How to most efficiently implement such integrals on numerically specified data (see above) over iso-surfaces in Mathematica?

EDIT: Trying out the proposal by Mastrok I run into convergence problems for e.g. $\epsilon$ = 2.0 of the integration (most probably due to the singularity in the coordinate center and the additional derivatives on the interpolation function when using Gauss' Theorem):

e = ListInterpolation[data, {{-1, 1}, {-1, 1}, {-1, 1}}];
dive = Simplify[Div[#/Norm[#]^2& @ Grad[e[x, y, z], {x, y, z}], {x, y, z}] /. 
  Derivative[1][Abs][x_] :> x/Abs[x] /. {x->#1, y->#2, z->#3} &];
NIntegrate[Piecewise[{{dive[x, y, z], e[x, y, z] < 2.0}}, 0] ,
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1}

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.
  • 1
    $\begingroup$ If I'm not mistaken, the integral is equal to $\displaystyle\frac{\mathrm d}{\mathrm d\epsilon}\iiint_{E(x,y,z)\le\epsilon}\mathrm dV$, is that right? It might be easier to compute that instead, but I'm not sure. $\endgroup$
    – user484
    Apr 23, 2015 at 6:07
  • $\begingroup$ Could you add to the Question a plot of a typical complicated energy surface over which you wish to integrate? $\endgroup$
    – bbgodfrey
    Apr 23, 2015 at 6:12
  • $\begingroup$ The literal implementation Integrate[1/Norm[Grad[e[x, y, z], {x, y, z}]], {x, y, z} ∈ ImplicitRegion[e[x, y, z] == ϵ, {x, y, z}]] seems to work when $E$ is given analytically, but not when it is an InterpolatingFunction constructed from the data. $\endgroup$
    – user484
    Apr 23, 2015 at 6:32
  • $\begingroup$ @bbgodfrey: Example added $\endgroup$
    – Rainer
    Apr 23, 2015 at 6:33
  • $\begingroup$ @Rahul, I did not check if the Gauss theorem yields an easier representation through a volume integral. To my view one would just exchange the problem of integrating over a complicated surface by the problem of integrating over a volume integral with bounds of similiar complexity. $\endgroup$
    – Rainer
    Apr 23, 2015 at 6:34

3 Answers 3


First note that,

$$ \oint_{E(\vec{x})=\epsilon}\frac{dS}{|\nabla E(\vec{x})|} =\oint_{E(\vec{x})=\epsilon}\frac{\nabla E(\vec{x})\cdot d\vec{S}}{|\nabla E(\vec{x})|^2} = \int_{E(\vec{x})\leq\epsilon} \nabla\cdot\left( \frac{\nabla E(\vec{x})}{|\nabla E(\vec{x})|^2}\right) dV$$

By now, the volume integral can be evaluated by

Integrate[f Boole[E[x,y,z]<epsilon],{x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax}]

where $f$ is the integrant and


forces the integrant to be 0 outside the region. Therefore, you need to have large enough {x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax} to cover your whole domain.

  • $\begingroup$ Note also that if you really want to do surface integral numerically, you can use $\delta$-function to impose the surface to be integrated. For the domain problem, since you are dealing with lattice (I guess), then it is much easier since you just need to set the domain to be one lattice size. $\endgroup$
    – mastrok
    Apr 23, 2015 at 6:38
  • $\begingroup$ thanks for your comprehensive answer. However when adopting your method to my example above, the numerical integration does not converge. I assume that the singularity in the coordinate center is causing this problem (see update in my question on this topic). $\endgroup$
    – Rainer
    Apr 26, 2015 at 8:36

By using the output of ListContourPlot3D and the new Mathematica 10.0 feature DiscretizeGraphics, one can nicely generate a meshed contour region which is suitable for NIntegrate. We can show this for the above example for energy contours from 1.0 to 2.0 with a step width of 0.1:

  e = ListInterpolation[data, {{-1, 1}, {-1, 1}, {-1, 1}}];
  f = Simplify[1/Norm[Grad[e[x, y, z], {x, y, z}]]];
  plot = ListContourPlot3D[data, Contours -> {ϵ}, 
   DataRange -> {{-1, 1}, {-1, 1}, {-1, 1}}];
  R = DiscretizeGraphics[Normal[plot /. (Lighting -> _) :> Lighting -> Automatic]];
  NIntegrate[f, {x, y, z} ∈ R], ϵ, 1.0, 2.0, 0.1}], ϵ]
, {ϵ, 1.0, 2.0, 0.1}], ϵ]
{12.5221, 8.55496, 6.03428, 4.38027, 3.25017, 2.46631, 1.90372, 1.49193, 1.18722, 0.952876, 0.777326}

The result is very close to the analytical solution. The differences are coming from the low sampling rate of the original data generation:

anasol = Table[
  Simplify[1/Norm[Grad[1/Sqrt[x^2 + y^2 + z^2], {x, y, z}]]], 
   {x, y, z} ∈ 
    ImplicitRegion[x^2 + y^2 + z^2 == (1/ϵ)^2, {x, y, z}]],
  {ϵ, 1.0, 2.0, 0.1}]
{12.5664, 8.583, 6.06017, 4.39984, 3.27113, 2.48225, 1.91748, 1.50458, 1.19707, 0.964263, 0.785398}

I've been try to perform calculations of exactly this type, and trying the different methods suggested. The energy surface I need is just as convoluted as the silicon Fermi surface shown in the OP and highly computationally intensive, and the scalar function I am trying to calculate the integral of over the surface (the function called $A(x,y,z)$ in the OP) is also very computationally intensive. Finding a mesh for the surface, and the gradient were too much, and there isn't enough control provided in the available routines. So in my case none of the methods above have been able to work with anything like a reasonable convergence in an acceptable computation time. The best method for speed that I've found is kind of a brute-force method, and it produces the entire spectrum in a smaller timescale than the surface integral method for a single energy (both in my case, and for the OP minimal code data). As suggested in earlier replies, one can turn the surface integral into a volume integral. Unlike the volume integral version suggested by Mastrok, which still requires finding the surface and gradient (which was far too computationally demanding for me), a much simpler volume integral is produced by $$ \int_{E(x,y,z)=\epsilon}dS \frac{A(x,y,z)}{|\nabla E(x,y,z)|}= \int A(x,y,z)\delta(E(x,y,z)-\epsilon) dV $$ where $\delta(E(x,y,z)-\epsilon)$ is a Dirac delta function. The key here is to allow the delta function to have a finite width in energy, and remember that it has unit area so that it has a height equal to the reciprocal of the width. Then you can just evaluate A and E on a grid, and then bin and sum the A values according to the E's (treating x,y,z as dummy parameters).

r = 0.1; (* resolution in x,y,z *)
data = 
 Table[1/Sqrt[x^2 + y^2 + z^2], {x, -1, 1, r}, {y, -1, 1, r}, {z, -1, 
     1, r}] /. ComplexInfinity -> 10^16 // Quiet; 
dataA = 
 Table[1, {x, -1, 1, r}, {y, -1, 1, r}, {z, -1, 1, r}];
ew = 0.1; (* energy width of the delta function *)
surfint = HistogramList[
   WeightedData[Flatten[data], Flatten[dataA]], {1 - ew/2, 2 + ew/2, 
surfint = surfint[[2]] r^3/ew; (* take care of the resolution of the grid and the delta function width *)
ListPlot[{surfint, anasol}, DataRange -> {1.0, 2.0}]

The above code produces the following (where anasol, shown in yellow, is the exact result from another earlier answer). Note that in this example dataA is just a grid of 1's so that the result could be compared with the exact solution, but it could be any function of x,y,z. The result is a little more noisy than other answers, but far faster, and so the noise can be improved by reducing the resolution and playing with bin widths. enter image description here


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.