As part of an equation of state calculator that I am working on, I have to handle a virial expansion for the equation of state. As such I have to solve for the virial coefficent, $B(T)$, which is defined as

$$B(T)=-2\pi N_a\int_0^{\infty}(e^{\frac{-\phi}{\kappa_b T}}-1)R^2dR$$

Where $\phi$ is the radial potential of a given electron at a given radius away from a given atom. This function is fairly complecated, and so there are many approximations of this formula, several of which include a heavside step function. For example, if I assume that only hard sphere inteactions hold in the system

$$\phi = \infty\left(\Theta(R-0)-\Theta(R-\sigma)\right)$$

meaning that the potential is infinite for all points from radius $R=0$ to atomic radius $\sigma$, and zero for all other points. However, I am having trouble solving the integral that comes from these symbolically using Mathematica. In this case, the system should solve out as

$$\begin{align}B(T) &=-2\pi N_a\int_0^{\infty}(e^{\frac{-\infty\left(\Theta(R-0)-\Theta(R-\sigma)\right)}{\kappa_b T}}-1)R^2dR\\&= -2\pi N_a\left(\int_0^{\sigma}(e^{\frac{-\infty\left(1-0)\right)}{\kappa_b T}}-1)R^2dR + \int_{\sigma}^{\infty}(e^{\frac{-\infty\left(1-1)\right)}{\kappa_b T}}-1)R^2dR\right)\\&= -2\pi N_a\left(\int_0^{\sigma}(e^{-\infty}-1)R^2dR + \int_{\sigma}^{\infty}(e^{0}-1)R^2dR\right)\\&= -2\pi N_a\left(\int_0^{\sigma}(-1)R^2dR + \int_{\sigma}^{\infty}(0)R^2dR\right)\\&= -2\pi N_a\left(\int_0^{\sigma}(-1)R^2dR + 0\right)\\&= 2\pi N_a\int_0^{\sigma}R^2dR\\&= \frac{2\pi}{3} N_a\sigma^3 \end{align}$$

where $N_a$ is Avagadros Constant. What is important to note about this is that I have some 10 or so different sets of assumptions that I can make to define different versions of $\phi$, and that i would like to be able to solve these systems symbolically using Mathematica.

I have reduced what I have worked out so far to a near minimal example which you can see below.

$Assumptions := {Inequality[1000000, Greater, Temperature, GreaterEqual, 0], Radial\[LetterSpace]Position >= 0, Atomic\[LetterSpace]Radius >= 0, Second\[LetterSpace]Step >= 0}
Subscript[k, b] := Quantity[1, "BoltzmannConstant"]
Subscript[A, n] := Quantity[1, "AvogadroConstant"]
T := Quantity[Temperature, "Kelvins"]
R := Quantity[Radial\[LetterSpace]Position, "Picometers"]
\[Sigma] := Quantity[Atomic\[LetterSpace]Radius, "Picometers"]
Step[r_] := HeavisideTheta[QuantityMagnitude[UnitConvert[r, "Picometers"]]]
Potential\[LetterSpace]Gradient\[LetterSpace]Def = {\[Phi] -> Infinity*(Step[R + 0] - Step[R - \[Sigma]])}
Virial\[LetterSpace]B\[LetterSpace]Def = {B[T] -> -2*Pi*Subscript[A, n]*Integrate[(E^(-(\[Phi]/(Subscript[k, b]*T))) - 1)*R^2, {R, 0, Infinity}]}
Simplify[Virial\[LetterSpace]B\[LetterSpace]Def /. Potential\[LetterSpace]Gradient\[LetterSpace]Def]

This appears in the IDE as

Second Compiled Version

and when executed, this produces the error

Integrate: Missing or incompatible quantities encountered in integration limits {Radial_Position pm,0,Infinity}.

and fails to produce the desired result, instead producing a function somewhere along the first line of the derivation above.

What approaches can I use to solve this type of system symbolically using Mathematica?

Approaches I have tried

I have tried restricting the bounds of several of the variables using the $Assumptions tag, declaring the relevant values as being Constants, and redefining the $\phi$ function as using piecewise notation to see if any of these help with the integration, however these approaches have as of yet proved fruitless.

I have tried to convert the integral from $0$ to $\infty$ to the limit of the integral from $0$ to $x$ as $x\to\infty$. As of yet this has also not provided an appropriate answer.

  • $\begingroup$ I would introduce dimensionless quantities and consider Integrate[ r^2 (Exp[-\[Beta] v HeavisideTheta[ r - 0] HeavisideTheta[\[Sigma] - r]] - 1), {r, 0, Infinity}, Assumptions -> {\[Beta] > 0, v > 0, \[Sigma] > 0}]. You can then take the limit ` v -> Infinity` and recover your result. $\endgroup$ Oct 22, 2019 at 15:41
  • $\begingroup$ That seems like it should be a viable approach - let me tool around with it for a bit to see if I can get it to work with a symbolic approach, and I'll get back to you if I have any question. Thank you $\endgroup$ Oct 22, 2019 at 15:47
  • $\begingroup$ @b.gates.you.know.what I implemented this using Limit\[LetterSpace]of\[LetterSpace]Infinite\[LetterSpace]Integral = {Integrate[func_., {var_., 0, Infinity}] :> Limit[Integrate[func, {var, 0, int\[LetterSpace]var}], int\[LetterSpace]var -> Infinity]} and Simplify[Virial\[LetterSpace]B\[LetterSpace]Def /. Limit\[LetterSpace]of\[LetterSpace]Infinite\[LetterSpace]Integral /. Potential\[LetterSpace]Gradient\[LetterSpace]Def], and was able to get a good limit equation, but mathematica still is not handling the exponential heaviside step formulas in the integral. $\endgroup$ Oct 22, 2019 at 16:01
  • $\begingroup$ The new output takes the form of {B[Quantity[Temperature, "Kelvins"]] -> Limit[Integrate[(-1 + E^((-Infinity)*(-1 + HeavisideTheta[-Atomic\[LetterSpace]Radius + Radial\[LetterSpace]Position])* Quantity[-(1/Temperature), 1/("BoltzmannConstant"*"Kelvins")]))*Quantity[Radial\[LetterSpace]Position^2, "Picometers"^2], {Quantity[Radial\[LetterSpace]Position, "Picometers"], 0, int\[LetterSpace]var}], int\[LetterSpace]var -> Infinity]* Quantity[-2*Pi, "AvogadroConstant"]} and again produces the error mentioned above. $\endgroup$ Oct 22, 2019 at 16:04
  • $\begingroup$ @b.gates.you.know.what - as the above text is a bit hard to read, I have updated the picture above to include the above approach, and its output $\endgroup$ Oct 22, 2019 at 16:11

2 Answers 2


First of, do not use units in Mathematica code, this will only complicate things. I am sure you can keep track of units without having them spelled out within Mathematica.

Note that

$$ⅇ^{\frac{-\infty\left(\Theta(R-0)-\Theta(R-\sigma)\right)}{\kappa_b T}}=\Theta(-R)+\Theta(R-\sigma)$$

as you show in your calculation. So let us define an expression representing that

expϕ = HeavisideTheta[-r]+HeavisideTheta[r-s];

Then you can do your desired integral symbolically, first without taking the boundary values

indefiniteI = -2 π Subscript[Ν, a] Integrate[(expϕ - 1) r^2, r]

enter image description here

The boundary at zero is simple to evaluate:

Assuming[s > 0, indefiniteI /. r -> 0 // Simplify]


So we know that the result is going to be given by the boundary value at infinity. Trying to plug in r->Infinity right away confuses Mathematica, so instead we substitute in the fact that both step functions become equal to 1 as r becomes large

indefiniteI /. HeavisideTheta[_] -> 1

enter image description here

which is the desired result.

Alternatively, instead of directly replacing the step functions, you could expand the step function arguments to leading order in r, simplify, and then set r to infinity:

(indefiniteI /. HeavisideTheta[arg_] :> HeavisideTheta[ Series[arg, {r, Infinity, -Exponent[arg, r]}] // Normal] // Simplify)/.r->Infinity

enter image description here

which leads to the same result.

You could go through these steps 10 times, with updated expϕ expressions, or make a vector out of expϕ with the different cases as components and do the calculation in one go for example.

  • 1
    $\begingroup$ This did not quite answer the question in the way that I was hoping that it would, but it got me to the solution - I had missed that the $\phi$ definition went to an indeterminate value as written and forgotten to assume that $\sigma > 0$ . Thank you for your help $\endgroup$ Oct 23, 2019 at 17:51

From the input of both b.gates.you.know, and Kagaratsch, I was able to find that the issue lied in two places with my code - namely that the definition of $\phi$ went to an indeterminite value at $R>\sigma$, and that I was failing to assume that $\sigma >0$. From this, I can get the code down to

$Assumptions := {T > 0, Subscript[k, b] > 0, \[Sigma] > 0}
Potential\[LetterSpace]Gradient\[LetterSpace]Def := {\[Phi] -> Piecewise[{{Infinity, R < \[Sigma]}, {0, R >= \[Sigma]}}]}
Virial\[LetterSpace]B\[LetterSpace]Def := {B[T] -> -2*Pi*Subscript[A, n]*Integrate[(E^(-(\[Phi]/(Subscript[k, b]*T))) - 1)*R^2, {R, 0, Infinity}]}
Virial\[LetterSpace]B\[LetterSpace]Def /. Potential\[LetterSpace]Gradient\[LetterSpace]Def


Potential\[LetterSpace]Gradient\[LetterSpace]Def := {\[Phi] -> Piecewise[{{Infinity, R < \[Sigma]}, {0, R >= \[Sigma]}}]}
Virial\[LetterSpace]B\[LetterSpace]Def := {B[T] -> -2*Pi*Subscript[A, n]*Integrate[(E^(-(\[Phi]/(Subscript[k, b]*T))) - 1)*R^2, {R, 0, Infinity}]}
Assuming[{R > 0, Subscript[k, b] > 0, \[Sigma] > 0}, Virial\[LetterSpace]B\[LetterSpace]Def /. Potential\[LetterSpace]Gradient\[LetterSpace]Def]

Each of which renders the desired output of


while maintaining the symbolic nature of the solution

Thank you to both @B.Gates.You.Know and @Kagaratsch for helping me find this solution


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.