The Gaussian function $\propto\exp((x-a)^2/b)$ with $b>0$ is one of the most common functions used in molecular modelling (e.g. Gaussian type orbitals).
What are some examples of applications of functions (in the literature) used other than the Gaussian in molecular modelling, and are there explanations as to why those are preferred over the Gaussian?
There are two considerations guiding the choice a basis for orbital expansion: 1. Compactness; 2. Efficiency of computations.
There are two common choices for basis functions (A) Gaussians and (B) Plane waves. Both of those allow for the most efficient way to evaluate the integrals needed to construct the Fock matrix (i.e. second derivative for kinetic energy, nuclear-electron attraction and electron-electron repulsion) - analytically.
(A) Gaussians is by far the most popular set for localized systems (e.g. molecules). The reason is because Gaussian functions are themselves localized and thus they can be used to compactly represent electron density localized around nuclei (i.e. you don't need too many Gaussian functions to do this). Plane waves resolve all of the simulation volume with equal accuracy, so a large part of the computational effort is wasted on (nearly) empty space. (B) Plane waves are much more popular for computations on periodic (condensed matter) systems (e.g. crystal structures). They enable even more efficient computations, since (unlike Gaussians) they are orthonormal and so there is no need to compute the overlap matrix S.
In practice, the sharp features (like the cusp near the nucleus) are still tough for the plane waves. This problem is circumvented by using the pseudopotentials, which smooth out the sharp features.
Another good reason to use either Gaussian or plane wave basis sets is that a lot of work has been done for you to develop really good Gaussian contractions and pseudopotentials. For a typical system there is no need to reinvent the wheel. That said, better basis sets are sought after and situationally (e.g. for a strange system) a different basis set could be much better.
I would like to expand on Roman Korol's answer a bit. He already lists GTO's and plane waves as they are the most common kind of basis functions. These are characteristic for the underlying models by which they are motivated. GTO's approximate the solutions to the hydrogen atom and are thus atom centered functions used for molecules. Plane waves on the other hand derive from the (nearly) free electron gas in a lattice and as such are suited for periodic systems.
There are, however, some lesser used types of basis functions. The ADF suite uses Slater type orbitals. While the idea of GTOs and STOs is of course similar, I would list them separately as their treatment in the code is quite different and, supposedly, they yield more accurate results. The reason GTOs are more commonly used is that for Gaussians you can exploit the Gaussian product theorem to evaluate the molecular integrals. It states (roughly) that the product of two Gaussians is another Gaussian centered between the two original ones. This can be used to get closed expressions for the integrals and calculating them becomes a matter of using recurrence relations. STOs are less convenient and the integrals have to be computed numerically.
The other type of basis functions I have seen are wavelets. These functions strike a balance between localization in real space (like GTOs) and in momentum space (like plane waves). They are used in the BigDFT program, but I am not very familiar with their properties.
Another slightly less commonly seen basis are the sinc functions, which are related to plane waves, but come at the problem from the perspective of a position, rather than momentum, space. They are delocalized functions, but sharply peaked at their center point and zero valued at the centers of other sinc functions, effectively partitioning space into a grid.
These have seen use in the Discrete Variable Representation method [1], which has the potential energy matrix elements being defined only at grid points and kinetic energy elements being defined by an infinite order finite difference for the second derivative.
Sinc functions have been mixed with localized functions like Gaussians [2] to try to combine the best parts, similar to the wavelets mentioned in Fuzzy's answer.
In the ONETEP program [3], periodic sinc functions are used as a basis to form nonorthogonal generalized Wannier functions (NGWFs), which are in turn optimized over a localization sphere during the calculation. The localization of these NGWFs allows for linear scaling DFT calculations, as it makes apparent the exponential decay of the density matrix.
The question was about "orbital basis sets" but explicitly mentions Gaussians; I guess the topic here is atomic orbital basis sets. In this case, the molecular orbitals are expanded as a linear combination of atomic orbitals (LCAO) as $ \psi_i({\bf r}) = \sum_{\alpha} C_{\alpha i} \chi_\alpha({\bf r})$; minimizing the Hartree-Fock / density functional energy with respect to the expansion coefficients $C_{\alpha i}$ typically leads to Roothaan-Hall type equations ${\bf FC}={\bf SCE}$ where ${\bf F}$ and ${\bf S}$ are the Fock and overlap matrices and ${\bf E}$ is a diagonal matrix of the energy eigenvalues.
Atomic-orbital basis sets have the form $ \chi_{\alpha}^{nlm}({\bf r})=R_{nl}(r)Y_l^m(\hat{\bf r})$ (non-relativistic case), and since few atomic orbitals already give qualitatively good results, you don't need a huge number of them, and you can diagonalize the Fock matrix exactly. You can also use atomic-orbital basis sets in the context of solid state calculations, in which case you actually use a periodic, symmetry adapted version $\chi_{\alpha}^{nlm}({\bf r};{\bf k})=\sum_{\bf g} e^{i{\bf k}\cdot{\bf g}} \chi_{\alpha}^{nlm}({\bf r}-{\bf R}_\alpha-{\bf g})$, where ${\bf R}_\mu$ are the coordinates of the center of the basis function, ${\bf g}$ sums over the lattice vectors, and ${\bf k}$ is the crystal momentum; but crystalline case is largely analogous to the molecular one.
The basis functions used in the expansion $ \psi_i({\bf r}) = \sum_{\alpha} C_{\alpha i} \chi_\alpha({\bf r})$ don't have to be atomic orbitals - using e.g. finite element basis functions or plane waves is perfectly fine as well - but then your basis set can become very large (10k, 100k, even 1M basis functions!), which makes solving the eigenvalue problem in the Roothaan-Hall equation prohibitively costly; this is why finite element and plane wave calculations typically use other ways to find the self-consistent solution.
Three kinds of radial functions are commonly used: Gaussian-type orbitals (GTOs) $R_{nl} = r^l \exp(-\zeta r^2)$, Slater-type orbitals (STOs) $R_{nl} = r^n \exp(-\zeta r)$, and numerical atomic orbitals (NAOs) $R_{nl} = u_{nl}(r)/r$. If you have few basis functions, then NAO $\gg$ STO $>$ GTO, since
However, if you use a large radial expansion, I don't think that the form of the individual basis functions matters that much, as large expansions accrue enough variational freedom to describe whatever electronic structure you throw at the problem. Traditionally one avoids large expansions due to pathological problems with overcompleteness; however, I have recently shown that a simple numerical trick can be used to overcome this issue in J. Chem. Phys. 151, 241102 (2019) and Phys. Rev. A 101, 032504 (2020).
GTOs have been the overwhelming favorite in calculations, since despite their drawbacks, they carry the huge benefit of analytic integrals evaluation: the integrals are both fast to calculate, and numerically exact. STOs and NAOs, on the other hand, require numerical quadrature. While quadrature can be made adaptively better, in many computations you also need to evaluate nuclear forces and Hessians, and especially the second derivatives for the Hessian can be tricky to compute accurately enough. The situation is, however, changing: while STOs have not become mainstream, there's been a lot of work with NAOs which have shown great accuracy and become usable even for coupled-cluster calculations, see the FHI-aims program.
For more details and discussion on solving the SCF equations, see e.g. our recent open access overview paper in Molecules 25, 1218 (2020).
You can also find a lengthier discussion on Gaussian-type orbitals, Slater-type orbitals and numerical atomic orbitals as well as other numerical approaches for electronic structure in my other recent open access review paper, Int. J. Quantum Chem. 119, e25968 (2019).
London Orbitals or Gauge-Including Atomic Orbitals (GIAOs)
These are used for computing magnetic properties of molecules.
In practical electronic structure calculations where a finite basis set is used, magnetic properties are not origin invariant.
In practice that means that if you computed, say, magnetizabilities or a circular dichroism spectra of a molecule centered at (0,0,0), and then re-did the calculation at, say 100 Angstroms away, (100,0,0), then you would get wildly different answers for the magnetic properties even if you have the same molecular geometry, wave function, energy, etc.
This is obviously not OK.
One solution is to make the orbitals themselves magnetic-field dependent by pre-multiplying the field-free atomic orbitals (commonly Gaussian, but not necessarily) with a magnetic-field dependent phase factor, or plane wave. For all practical purposes, this eliminates the origin dependence.
Mathematically, a London orbital or GIAO, $\chi^{\textrm{GIAO}}(\mathbf{r} - \mathbf{R})$, centered at $\mathbf{R}$ looks like a plane-wave/atomic-orbital hybrid:
$$\chi^{\textrm{GIAO}}(\mathbf{r} - \mathbf{R}) = \mathrm{exp}\left(\frac{i}{2}\left(\mathbf{R} \times \mathbf{B} \right) \cdot \left(\mathbf{r} - \mathbf{R}\right)\right)\chi^{\textrm{AO}}(\mathbf{r} - \mathbf{R})$$
Where the field-free atomic orbital $\chi^{\textrm{AO}}(\mathbf{r} - \mathbf{R})$ is multiplied by the magnetic field $\mathbf{B}$ dependent phase factor. (And $\mathbf{r}$ is the electronic coordinate vector.)
Implementations:
Lots of software, like Gaussian or DALTON, utilize GIAOs when computing molecular magnetic properties such as magnetizabilities, chiroptical properties, etc. Because these properties are defined in the limit as $\mathbf{B}\to \mathbf{0}$, you don't have to use a special London orbital basis set, as it's already taken into account by the theoretical method.
On the other hand, finite magnetic field calculations are much less common, but here are some programs:
A few papers:
