How do biaxial potentials like the Gay-Berne potential work?

I am trying to understand and implement the Gay-Berne potential, but I do not quite understand how one goes about constructing them. I am following this paper to parse through the idea.

I understand how Lennard-Jones/Coulombic potentials work in MD engines: every particle has some $$\epsilon, \sigma, q$$ associated with it, and these potentials are purely a function of the distance $$r=|\mathbf{r}_{ij}|$$ between two particles. It does not depend upon the orientation of the two particles. $$U = U(r)$$

In the above paper, a particle $$i$$ at position $$\mathbf{r}_i$$ and orientation $$\omega$$ with respect to the laboratory frame is represented by a Gaussian function $$G_i (\omega _i, \mathbf{r})$$ which is $$G_i(\omega _i, \mathbf{r}) = \exp [-\frac{1}{2}(\mathbf{r}-\mathbf{r}_i)^T M_i ^T S_i ^{-2} M_i (\mathbf{r}-\mathbf{r}_i)]$$ where S is a diagonal 'shape' matrix with elements $$\sigma _x, \sigma _y, \sigma _z$$, the axes of the ellipsoid representing the molecule measures with respect to a unit length $$\sigma _0$$, and $$M_i$$ is the rotation matrix transforming from laboratory to a molecular frame.

What I don't understand is, why do we need a Gaussian function to denote a particle? A 3-dimensional ellipse can be represented with the shape matrix S, and its orientation which can be given by the matrix M. Where does $$\omega$$ figure into the above definition? If anybody has any recommendations on how to go about understanding the implementation of this idea, I would greatly appreciate it.

To start, they mention just after equation 6 of your linked paper that $$\omega=(\alpha, \beta,\gamma)$$, the Euler angles that define the orientation. This set of Euler angles is used to define the rotation matrix $$M$$ in a succinct way.