SOAP and LODE radial integrals

On this page, we describe the exact mathematical expression that are implemented in the radial integral and the splined radial integral classes i.e. Splined radial integrals. Note that this page assumes knowledge of spherical expansion & friends and currently serves as a reference page for the developers to support the implementation.

Preliminaries

In this subsection, we briefly provide all the preliminary knowledge that is needed to understand what the radial integral class is doing. The actual explanation for what is computed in the radial integral class can be found in the next subsection (1.2). The spherical expansion coefficients anlm|ρi are completely determined by specifying two ingredients:

  • the atomic density function g(r) as implemented in Atomic density, often chosen to be a Gaussian or Delta function, that defined the type of density under consideration. For a given central atom i in the system, the total density function ρi(r) around is then defined as ρi(r)=jg(rrij).

  • the radial basis functions Rnl(r) as implementated Radial basis functions, on which the density ρi is projected. To be more precise, the actual basis functions are of the form

    Bnlm(r)=Rnl(r)Ylm(r^),

    where Ylm(r^) are the real spherical harmonics evaluated at the point r^, i.e. at the spherical angles (θ,ϕ) that determine the orientation of the unit vector r^=r/r.

The spherical expansion coefficient nlm|ρi (we omit the atom type index a for simplicity) is then defined as

nlm|ρi=d3rBnlm(r)ρi(r)=d3rRnl(r)Ylm(r^)ρi(r).

In practice, the atom centered density ρi is a superposition of the neighbor contributions, namely ρi(r)=jg(rrij). Due to linearity of integration, evaluating the integral can then be simplified to

nlm|ρi=d3rRnl(r)Ylm(r^)ρi(r)=d3rRnl(r)Ylm(r^)(jg(rrij))=jd3rRnl(r)Ylm(r^)g(rrij)=jnlm|g;rij.

Thus, instead of having to compute integrals for arbitrary densities ρi, we have reduced our problem to the evaluation of integrals of the form

nlm|g;rij=d3rRnl(r)Ylm(r^)g(rrij),

which are completely specified by

  • the density function g(r)

  • the radial basis Rnl(r)

  • the position of the neighbor atom rij relative to the center atom

The radial integral class

In the previous subsection, we have explained how the computation of the spherical expansion coefficients can be reduced to integrals of the form

nlm|g;rij=d3rRnl(r)Ylm(r^)g(rrij).

If the atomic density is spherically symmetric, i.e. if g(r)=g(r) this integral can always be written in the following form:

nlm|g;rij=Ylm(r^ij)Inl(rij).

The key point is that the dependence on the vectorial position rij is split into a factor that contains information about the orientation of this vector, namely Ylm(r^ij), which is just the spherical harmonic evaluated at r^ij, and a remaining part that captures the dependence on the distance of atom j from the central atom i, namely Inl(rij), which we shall call the radial integral. The radial integral class computes and outputs this radial part Inl(rij). Since the angular part is just the usual spherical harmonic, this is the part that also depends on the choice of atomic density g(r), as well as the radial basis Rnl(r). In the following, for users only interested in a specific type of density, we provide the explicit expressions of Inl(r) for the Delta and Gaussian densities, followed by the general expression.

Delta Densities

Here, we consider the especially simple special case where the atomic density function g(r)=δ(r). Then:

nlm|g;rij=d3rRnl(r)Ylm(r^)g(rrij)=d3rRnl(r)Ylm(r^)δ(rrij)=Rnl(r)Ylm(r^ij)=Bnlm(rij).

Thus, in this particularly simple case, the radial integral is simply the radial basis function evaluated at the pair distance rij, and we see that the integrals have indeed the form presented above.

Gaussian Densities

Here, we consider another popular use case, where the atomic density function is a Gaussian. In featomic, we use the convention

g(r)=1(πσ2)3/4er22σ2.

The prefactor was chosen such that the “L2-norm” of the Gaussian

g2=d3r|g(r)|2=1,

but does not affect the following calculations in any way. With these conventions, it can be shown that the integral has the desired form

nlm|g;rij=d3rRnl(r)Ylm(r^)g(rrij)=Ylm(r^ij)Inl(rij)

with

Inl(rij)=1(πσ2)3/44πerij22σ20drr2Rnl(r)er22σ2il(rrijσ2),

where il is a modified spherical Bessel function. The first factor, of course, is just the normalization factor of the Gaussian density. See the next two subsections for a derivation of this formula.

Derivation of the General Case

We now derive an explicit formula for radial integral that works for any density. Let g(r) be a generic spherically symmetric density function. Our goal will be to show that

nlm|g;rij=Ylm(r^ij)[2π0drr2Rnl(r)11d(cosθ)Pl(cosθ)g(r2+rij22rrijcosθ)]

and thus we have the desired form nlm|g;rij=Ylm(r^ij)Inl(rij) with

Inl(rij)=2π0drr2Rnl(r)11duPl(u)g(r2+rij22rriju),

where Pl(x) is the l-th Legendre polynomial.

Derivation of the explicit radial integral for Gaussian densities

Denoting by θ(r,rij) the angle between a generic position vector r and the vector rij, we can write

g(rrij)=1(πσ2)3/4e(rrij)22σ2=1(πσ2)3/4e(rij)22σ2e(r22rrij)2σ2,

where the first factor no longer depends on the integration variable r.

Analytical Expressions for the GTO Basis

While the above integrals are hard to compute in general, the GTO basis is one of the few sets of basis functions for which many of the integrals can be evaluated analytically. This is also useful to test the correctness of more numerical implementations.

The primitive basis functions are defined as

Rnl(r)=Rn(r)=rner22σn2

In this form, the basis functions are not yet orthonormal, which requires an extra linear transformation. Since this transformation can also be applied after computing the integrals, we simply evaluate the radial integral with respect to these primitive basis functions.

Real Space Integral for Gaussian Densities

We now evaluate

Inl(rij)=1(πσ2)3/44πerij22σ20drr2Rnl(r)er22σ2il(rrijσ2)=1(πσ2)3/44πerij22σ20drr2rner22σn2er22σ2il(rrijσ2),

the result of which can be conveniently expressed using a=12σ2, bn=12σn2, neff=n+l+32 and leff=l+32 as

Inl(rij)=1(πσ2)3/4π32Γ(neff)Γ(leff)(arij)l(a+b)neffM(neff,leff,a2rij2a2+b2),

where M(a,b,z) is the confluent hypergeometric function (hyp1f1).