diff --git a/.github/workflows/UnitTest.yml b/.github/workflows/UnitTest.yml
index 5119a31..868e072 100644
--- a/.github/workflows/UnitTest.yml
+++ b/.github/workflows/UnitTest.yml
@@ -7,8 +7,8 @@ on:
     branches:
       - main
   pull_request:
-  schedule:
-    - cron: '20 00 1 * *'
+    branches:
+      - main
 
 jobs:
   test:
@@ -17,19 +17,19 @@ jobs:
       fail-fast: false
       matrix:
         julia-version:
-          - '1.6' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
-          - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
-          - 'nightly'
+          - "1.9" # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - "1" # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
+          - "nightly"
         os:
           - ubuntu-latest
         arch:
           - x64
         include:
           - os: windows-latest
-            julia-version: '1'
+            julia-version: "1"
             arch: x64
           - os: macOS-latest
-            julia-version: '1'
+            julia-version: "1"
             arch: x64
     env:
       PYTHON: ""
diff --git a/Project.toml b/Project.toml
index 7814280..9bfe1b3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,11 +1,10 @@
 name = "E3NN"
 uuid = "1c50a8ea-cbe2-4d3e-83e0-d59f5e8851b3"
-authors = [
-    "Deeptendu Santra <deeptendu.santra@protonmail.com> and contributors",
-]
+authors = ["Deeptendu Santra <deeptendu.santra@protonmail.com> and contributors"]
 version = "0.1.0"
 
 [deps]
+BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLUtils = "f1d291b0-491e-4a28-83b9-f70985020b54"
@@ -17,6 +16,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [compat]
+BlockDiagonals = "0.2.0"
 CoordinateTransformations = "0.6"
 Quaternions = "0.7"
 Rotations = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index dc5c943..dfa65cd 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,2 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-e3nn = "1c50a8ea-cbe2-4d3e-83e0-d59f5e8851b3"
diff --git a/docs/make.jl b/docs/make.jl
index 5230522..ccc5532 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,20 +1,22 @@
-using Documenter, e3nn
+using Documenter, E3NN
 
 mathengine = MathJax3()
 prettyurls = get(ENV, "CI", nothing) == "true"
 
 makedocs(
-    modules = [e3nn],
+    modules = [E3NN],
     doctest = false,
     clean = true,
-    sitename = "e3nn.jl",
+    sitename = "E3NN.jl",
     format = Documenter.HTML(;
         mathengine,
         prettyurls,
         assets = ["assets/favicon.ico"],
         collapselevel = 3
     ),
-    pages = ["Home" => "index.md", "API Reference" => ["o3" => "api/o3.md"]]
+    pages = ["Home" => "index.md", "Irreps" => "irreps.md",
+        "Examples" => "examples.md",
+        "API Reference" => ["Irreps" => "api/irreps.md"]]
 )
 
 deploydocs(repo = "github.com/Dsantra92/e3nn.jl.git")
diff --git a/docs/src/api/irreps.md b/docs/src/api/irreps.md
new file mode 100644
index 0000000..7f19aac
--- /dev/null
+++ b/docs/src/api/irreps.md
@@ -0,0 +1,7 @@
+## Irreducible Representation
+
+```@autodocs
+Modules = [E3NN.O3 ]
+Pages   = ["O3/irreps.jl"]
+Private = false
+```
diff --git a/docs/src/api/o3.md b/docs/src/api/o3.md
index fdcdc7f..283b938 100644
--- a/docs/src/api/o3.md
+++ b/docs/src/api/o3.md
@@ -1,9 +1,9 @@
 ```@meta
-CurrentModule = e3nn
+CurrentModule = E3NN
 ```
 
 ```@autodocs
-Modules = [o3]
-Pages   = ["o3/o3.jl", "o3/irreps.jl"]
+Modules = [O3]
+Pages   = ["O3/O3.jl"]
 Private = false
 ```
diff --git a/docs/src/assets/irreps.png b/docs/src/assets/irreps.png
new file mode 100644
index 0000000..ecf9a92
Binary files /dev/null and b/docs/src/assets/irreps.png differ
diff --git a/docs/src/examples.md b/docs/src/examples.md
new file mode 100644
index 0000000..df635b4
--- /dev/null
+++ b/docs/src/examples.md
@@ -0,0 +1 @@
+# Examples
diff --git a/docs/src/index.md b/docs/src/index.md
index e25d663..bd13650 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,3 +1,25 @@
-# e3nn
+# E3nn: E(3) Equivariant Neural Networks
 
-Julia implementation of e3nn python library.
+[E3nn.jl](https://github.com/Dsantra92/e3nn.jl) provides a pure julia implementation of the [e3nn framework](https://e3nn.org/).
+It aims to provide fast and extendable APIs for working with $\mathbb{E}^3$ Equivariant Neural Networks and other performing related operations.
+It is built on top of [GraphNeuralNetworks.jl]() and [Flux.jl]().
+The APIs are generally consistent with e3nn's [PyTorch](https://github.com/e3nn/e3nn/) and [JAX](https://github.com/e3nn/e3nn-jax/) libraries but there might be certain differences.
+
+## Installation
+
+The package isn't available through the Julia Package Manager (yet).
+If you want to use it, you can add it directly from GitHub.
+
+```julia-repl
+julia> using Pkg
+
+julia> Pkg.add("https://github.com/Dsantra92/e3nn.jl.git")
+```
+
+!!! note
+The package is still in under early stages of development and the APIs might change in the future.
+
+## Getting started
+
+We highly recommend starting with [Introduction to Irreps]() to gain some understanding of the primary components of the network.
+If you are looking for a bit more hands on experience, have a look at the [examples]().
diff --git a/docs/src/irreps.md b/docs/src/irreps.md
new file mode 100644
index 0000000..ad06126
--- /dev/null
+++ b/docs/src/irreps.md
@@ -0,0 +1,182 @@
+# Irreducible Representations of O(3)
+
+A group representation $(D, V)$ describes the action of a group $G$ on a vector space $V$. Formally, this is written as:
+
+```math
+D: G \to \text{GL}(V)
+```
+
+where $\text{GL}(V)$ denotes the general linear group of $V$, i.e., the group of invertible linear transformations on $V$.
+
+A representation of a group is a mapping(a literal matrix) from the group elements to the general linear group of matrices.
+The irreducible representations, or _irreps_, are the "minimal" non-trivial representations of a group.
+They are fundamental in the sense that any finite-dimensional representation of a group can be decomposed into a direct sum of irreps.
+
+The orthogonal group $O(3)$ is the group of $3×3$ orthogonal matrices.
+Geometrically, $O(3)$ represents all rotations and reflections in three-dimensional space that preserve distances and angles.
+Any physical quantity, under the action of $O(3)$, transforms with a representation of $O(3)$
+
+The special orthogonal group $SO(3)$ is a subgroup of $O(3)$. $SO(3)$ represents all rotations in three-dimensional space, excluding reflections.
+$O(3)$ can be therefore be decomposed as a direct product of $SO(3)$ and the inversion group ${e, I}$:
+
+$O(3) \cong SO(3) \times \{e, I\}$
+
+where $e$ is the identity element of the group(here $3 \times 3$ identity matrix) and $I$ represents inversion(reflection) ($-1 \times$ identity matrix).
+
+So, any element $g \in O(3)$ can be uniquely written as $g = ri$, where $r \in SO(3)$ and $i \in {e, I}$.
+
+The irreps of $O(3)$ are the product of the irreps of $SO(3)$ and inversion.
+The irreps of the $SO(3)$ group are called [Wigner matrices](https://en.wikipedia.org/wiki/Wigner_D-matrix) $D^L$.
+The Inversion group has 2 irreps:
+
+- The trivial representation $\sigma_+$, where inversion does nothing (multiplies by 1)
+- The sign representation $\sigma_-$, where inversion flips the sign (multiplies by -1).
+
+```math
+\sigma_p(g) = \left \{ \begin{array}{l} 1 \text{ if } g = e \\ p \text{ if } g = I \end{array} \right..
+```
+
+Therefore, the irreps of an element $g \in O(3)$, where $g = ri$ is given by, $\sigma_p(i)D^L(r)$. In our library, we represent it using [`Irrep`](@ref).
+
+And, [`Irreps`](@ref) represents the direct sum of irreps of $O(3)$:
+
+```math
+g = r i \mapsto \bigoplus_{j=1}^n m_j \sigma_{p_j}(i) D^{L_j}(r)
+```
+
+where $n \in \mathbb{N}$, $m_j \in \mathbb{N}$, $p_j \in \{+1, -1\}$, and $L_j \in \mathbb{N}_0$ for $j = 1, \ldots, n$.
+
+![image](./assets/irreps.png)
+
+## Irrep
+
+The [`Irrep`](@ref) struct stores the degree $l$ and the parity $p$ of the representation $\sigma_p(i)D^l(r)$.
+
+For a small recap:
+
+- The variable $l$ represents the angular momentum quantum number or degree of the irrep.
+- For $l$ = 0 , the representation is scalar, for $l$ = 1, the representation is a vector. For higher $l$, the representations are tensors of increasing rank, with dimension (2$l$ + 1).
+- The parity $p$ indicates whether the irrep is invariant under reflection through the origin.
+- A positive parity $( p = 1 )$ means the object behaves symmetrically under inversion, while negative parity $( p = -1 )$ implies antisymmetry, flipping sign upon inversion.
+  Together, $l$ and $p$ fully characterize how a function or field transforms under both rotation and inversion.
+
+```jldoctest
+julia> using E3NN.O3
+
+julia> Irrep("1e")
+1e
+
+julia> Irrep(1, -1)
+1o
+
+```
+
+Get the dimension of the representaion.
+
+```jldoctest
+julia> dim(Irrep("1e"))
+3
+```
+
+## Irreps
+
+The Irreps struct represents a direct sum of irreducible representations.
+It is a collection of Irrep objects, each with an associated multiplicity.
+This is the primary way to describe the structure of feature vectors in an equivariant neural network.
+
+### Construction
+
+You can create an Irreps object in several ways:From a string, which is the most common method:
+
+```jldoctest
+julia> Irreps("3x0e + 2x1o")
+3x0e+2x1o
+```
+
+From a single Irrep, which will have a multiplicity of 1:
+
+```jldoctest
+julia> Irreps(Irrep("2e"))
+1x2e
+```
+
+From a list (Vector or Tuple) of Irrep or MulIrrep objects:
+
+```jldoctest
+julia> Irreps([Irrep("0e"), Irrep("1o")])
+1x0e+1x1o
+
+julia> Irreps([(2, Irrep("1e")), (3, Irrep("2o"))])
+2x1e+3x2o
+```
+
+An empty Irreps represents the zero-dimensional vector space.
+
+```jldoctest
+julia> Irreps()
+
+julia> sprint(show, Irreps())
+""
+```
+
+Properties and OperationsYou can inspect the properties of an Irreps object with several utility functions.The total dimension of the vector space is given by dim:
+
+```jldoctest
+julia> dim(Irreps("3x0e + 2x1o"))
+9
+```
+
+The total number of irreps (sum of multiplicities) is given by num\*irreps:
+
+```jldoctest
+julia> num_irreps(Irreps("3x0e + 2x1o"))
+5
+```
+
+The maximum `l` value is given by lmax:
+
+```jldoctest
+julia> lmax(Irreps("3x0e + 2x1o"))
+1
+```
+
+You can perform a direct sum using the + operator. The result is automatically regrouped into a canonical form.
+
+```jldoctest
+julia> Irreps("1x1o") + Irreps("2x0e") + Irreps("1x1o")
+2x0e+2x1o
+```
+
+You can scale the multiplicities using the \* operator.
+
+```jldoctest
+julia> 3 * Irreps("1x0e + 2x1o")
+3x0e+6x1o
+```
+
+### Manipulation and Filtering
+
+The library provides powerful functions to manipulate Irreps.
+`regroup` sorts the irreps and combines identical ones by summing their multiplicities.
+This is useful for creating a canonical representation.
+
+```jldoctest
+julia> regroup(Irreps("2x1o + 1x0e + 3x1o"))
+1x0e+5x1o
+```
+
+`filter` allows you to select a subset of the irreps based on various criteria.
+
+```jldoctest
+julia> irreps = Irreps("3x0e + 2x1o + 4x2e")
+3x0e+2x1o+4x2e
+
+julia> filter(irreps, lmax=1)
+3x0e+2x1o
+
+julia> filter(irreps, keep=Irreps("0e + 2e"))
+3x0e+4x2e
+
+julia> filter(irreps, drop=mulir -> mulir.mul > 2)
+2x1o
+```
diff --git a/src/E3NN.jl b/src/E3NN.jl
new file mode 100644
index 0000000..9050f59
--- /dev/null
+++ b/src/E3NN.jl
@@ -0,0 +1,8 @@
+module E3NN
+
+# O3
+include("O3/O3.jl")
+using .O3
+export Irrep
+
+end
diff --git a/src/O3/O3.jl b/src/O3/O3.jl
new file mode 100644
index 0000000..30af468
--- /dev/null
+++ b/src/O3/O3.jl
@@ -0,0 +1,8 @@
+module O3
+
+include("irreps.jl")
+export Irrep, Irreps
+export dim, isscalar, lmax, ls, count, mul_gcd, num_irreps,
+       remove_zero_multiplicities, simplify, sort, unify, regroup
+
+end
diff --git a/src/O3/irreps.jl b/src/O3/irreps.jl
new file mode 100644
index 0000000..51ecd98
--- /dev/null
+++ b/src/O3/irreps.jl
@@ -0,0 +1,367 @@
+import Base.Iterators
+
+"""
+Irrreducible representations of ``O(3)``.
+
+This struct does not contain any data; it is a structure that describes the representation.
+It is typically used as an argument in other parts of the library to define the input and output representations of
+functions.
+
+# Fields:
+- `l::Int`: non-negative integer, the degree of the representation, ``l = 0, 1, \\dots``
+- `p::Int`: the parity of the representation, either 1 (even) or -1 (odd)
+
+The degree ``l`` corresponds to the angular momentum quantum number and defines the dimension of the irrep, with ``l = 0`` representing a scalar, ``l = 1`` a vector, and higher ``l`` values corresponding to tensor-like objects.
+"""
+struct Irrep{T <: Integer}
+    l::T
+    p::T
+
+    @doc"""
+         Irrep(l::Integer, p::Integer)
+
+     Instantiate a new [`O3.Irrep`](@ref) object.
+
+     # Examples:
+     Create a scalar representation (``l=0``) of even parity:
+     ```jldoctest
+     julia> Irrep(0, 1)
+     0e
+     ```
+
+     Create a pseudotensor representation (``l=2``) of odd parity:
+     ```jldoctest
+     julia> Irrep(2, -1)
+     2o
+     ```
+     """
+    function Irrep(l::Integer, p::Integer)
+        (l >= 0) ||
+            throw(ArgumentError("l must be zero or positive integer, got $l"))
+        (p in (-1, 1)) ||
+            throw(ArgumentError("parity(p) must be one of (-1, 1) got $p"))
+        promoted_l, promoted_p = promote(l, p)
+        return new{typeof(promoted_l)}(promoted_l, promoted_p)
+    end
+end
+
+"""
+    Irrep(ir::T) where {T <: AbstractString}
+
+Instantiate a new [`O3.Irrep`](@ref) object from it's string representation.
+
+The string representation should be of the form `l` followed by "e" or "o" for even or odd parity, respectively.
+There can also be a "y" at the end, which is used to represent the parity of the spherical harmonics.
+
+- "e" for even parity, translates to `p`=1
+- "o" for odd parity, translates to `p`=-1
+- "y" for the parity of the spherical harmonics, translates to `p`=``(-1)^l``
+
+# Examples
+Create a vector representation (``l=1``) of the parity of the spherical harmonics (``-1^l``) gives odd parity):
+```jldoctest
+julia> Irrep("1y")
+1o
+```
+"""
+function Irrep(ir::AbstractString)::Irrep{Int}
+    name = strip(ir)
+    try
+        l = parse(Int, name[1:(end - 1)])
+        (l >= 0) ||
+            throw(ArgumentError("l must be zero or positive integer, got $l"))
+        p = Dict('e' => 1, 'o' => -1, 'y' => (-1)^l)[name[end]]
+        return Irrep(l, p)
+    catch
+        throw(ArgumentError("Cannot convert string $name to an Irrep"))
+    end
+end
+
+"""
+    Irrep(ir::Tuple{Integer, Integer})
+
+Instantiate a new [`O3.Irrep`](@ref) object from a tuple of (`l`, `p`).
+"""
+Irrep(ir::Tuple{Integer, Integer}) = Irrep(ir...)
+
+Irrep(ir::Irrep) = ir
+
+Base.show(io::IO, ir::Irrep) = print(io, "$(ir.l)$(ir.p == 1 ? 'e' : 'o')")
+
+function Base.iterate(ir::Irrep, state::Int = 1)
+    if state == 1
+        return (ir.l, 2);
+    elseif state == 2
+        return (ir.p, 3);
+    else
+        return nothing;
+    end
+end
+
+# Core collection traits
+# Just a pragmatic approach
+# Not sure if all of it is necessary
+Base.length(::Irrep) = 2
+Base.eltype(::Type{Irrep{T}}) where {T} = T
+Base.first(ir::Irrep) = ir.l
+Base.last(ir::Irrep) = ir.p
+
+function Base.isless(a::Irrep, b::Irrep)
+    key_a = (a.l, -a.p * (-1)^a.l)
+    key_b = (b.l, -b.p * (-1)^b.l)
+    return isless(key_a, key_b)
+end
+
+"""
+    dim(x::Irreps)
+
+Returns the dimension representation of `Irreps`, ``2l+1``
+"""
+dim(x::Irrep) = 2 * x.l + 1
+
+"""
+    isscalar(x::Irrep)
+
+Check if the irrep is a scalar representation.
+Equivalent to `x.l == 0 && x.p == 1`.
+"""
+isscalar(x::Irrep) = (x.l == 0) && (x.p == 1)
+
+function iterator(::Type{Irrep}; lmax::Union{Int, Nothing} = nothing)
+    return Channel{Irrep}() do ch
+        for l in 0:typemax(Int)
+            p1 = (-1)^l
+            put!(ch, Irrep(l, p1))
+            put!(ch, Irrep(l, -p1))
+            if !isnothing(lmax) && l == lmax
+                break
+            end
+        end
+    end
+end
+
+function (T::Type{Irrep})(; lmax::Union{Int, Nothing} = nothing)
+    iterator(T; lmax = lmax)
+end
+
+struct MulIrrep{T}
+    mul::Int # might reconsider using T here, but not sure if it's worth the constraint
+    ir::Irrep{T}
+end
+
+Base.show(io::IO, mul_ir::MulIrrep) = print(io, "$(mul_ir.mul)x$(mul_ir.ir)")
+
+Base.isless(a::MulIrrep, b::MulIrrep) = isless((a.ir, a.mul), (b.ir, b.mul))
+
+function Base.iterate(mul_ir::MulIrrep, state::Int = 1)
+    if state == 1
+        return (mul_ir.mul, 2);
+    elseif state == 2
+        return (mul_ir.ir, 3);
+    else
+        return nothing;
+    end
+end
+
+dim(mul_ir::MulIrrep) = mul_ir.mul * dim(mul_ir.ir)
+
+"""
+    Irreps
+
+Direct sum of irreducible representations of ``O(3)``.
+
+This struct does not contain any data, it is a structure that describes the representation.
+
+# Examples
+```jldoctest
+# Create a representation of 100 l=0 of even parity and 50 pseudo-vectors.
+julia> x = Irreps([(100, (0, 1)), (50, (1, 1))])
+100x0e+50x1e
+
+julia> dim(x)
+250
+
+# Create a representation of 100 l=0 of even parity and 50 pseudo-vectors.
+julia> Irreps("100x0e + 50x1e")
+100x0e+50x1e
+
+julia> Irreps("100x0e + 50x1e + 0x2e")
+100x0e+50x1e+0x2e
+
+julia> Irreps("100x0e + 50x1e + 0x2e") |> lmax
+1
+
+julia> Irrep("2e") in irs("0e + 2e")
+true
+
+# Empty Irreps
+julia> Irreps(), Irreps("")
+(, )
+```
+"""
+struct Irreps{T}
+    _irreps::Tuple{Vararg{MulIrrep{T}}}
+end
+
+Irreps(ir::Irrep{T}) where {T} = Irreps((MulIrrep(1, ir),))
+Irreps(mul_ir::MulIrrep{T}) where {T} = Irreps((mul_ir,))
+Irreps(mul_ir::Tuple{Int, Irrep{T}}) where {T} = Irreps((MulIrrep(mulir...),))
+Irreps() = Irreps{Int}(tuple()) # Default to Int for empty constructor
+Irreps(irs::Irreps{T}) where {T} = irs
+
+function Irreps(irs::AbstractString)
+    _irs = strip(irs)
+    isempty(_irs) && return Irreps()
+    mulirreps = Vector{MulIrrep{Int}}()
+    if _irs != ""
+        for mul_ir in split(_irs, "+")
+            if occursin("x", mul_ir)
+                mul, ir = split(mul_ir, "x")
+                mul = parse(Int, mul)
+                (mul >= 0) ||
+                    throw(ArgumentError("mul should be greater than 0 got $mul"))
+                irrep = Irrep(ir)
+            else
+                mul = 1
+                irrep = Irrep(mul_ir)
+            end
+            push!(mulirreps, MulIrrep(mul, irrep))
+        end
+    end
+    return Irreps(tuple(mulirreps...))
+end
+
+function Irreps(irs::Union{Vector, Tuple})
+    isempty(irs) && return Irreps()
+
+    # potentially mixed types
+    initial_list = map(irs) do x
+        if typeof(x) <: AbstractString
+            # fix case for mixture of irrep and irrpes in a vector
+            # if occursin("+", mul_irrep)
+            #     irreps = Irreps(mul_irrep)
+            # end
+            MulIrrep(1, Irrep(x))
+        elseif x isa MulIrrep
+            x
+        elseif x isa Tuple{Int, Irrep}
+            MulIrrep(x...)
+        elseif x isa Irrep
+            MulIrrep(1, x)
+        elseif length(x) == 2
+            mul, ir = x
+            MulIrrep(mul, Irrep(ir))
+        else
+            throw(ArgumentError("Invalid type in list for Irreps constructor: $(typeof(x))"))
+        end
+    end
+
+    # common promoted type for the Irrep's T
+    promoted_T = mapreduce(x -> eltype(x.ir), promote_type, initial_list)
+
+    final_list = map(initial_list) do mul_ir
+        new_ir = Irrep(
+            convert(promoted_T, mul_ir.ir.l), convert(promoted_T, mul_ir.ir.p))
+        MulIrrep(mul_ir.mul, new_ir)
+    end
+
+    return Irreps(tuple(final_list...))
+end
+
+function Base.show(io::IO, irs::Irreps)
+    isempty(irs._irreps) && return print(io, "")
+    print(io, join(string.(irs._irreps), "+"))
+end
+
+Base.iterate(irs::Irreps, args...) = iterate(irs._irreps, args...)
+Base.getindex(irs::Irreps, i) = getindex(irs._irreps, i)
+Base.length(irs::Irreps) = length(irs._irreps)
+Base.eltype(irs::Irreps) = eltype(irs._irreps)
+function Base.in(ir::Irrep, irs::Irreps)
+    any(mul_ir -> mul_ir.ir == ir, irs._irreps)
+end
+
+function Base.:+(a::Irreps, b::Irreps)
+    return Irreps((a._irreps..., b._irreps...))
+end
+
+function Base.:*(n::Integer, irs::Irreps)
+    @assert n >= 0 "Multiplicity must be non-negative."
+    return Irreps(tuple([MulIrrep(n * mul_ir.mul, mul_ir.ir) for mul_ir in irs]...))
+end
+
+Base.:*(irs::Irreps, n::Integer) = n * irs # Commutativity
+
+function Base.:*(n::Integer, ir::Irrep)
+    @assert n >= 0 "Multiplicity must be non-negative."
+    return Irreps((MulIrrep(n, ir),))
+end
+
+Base.:+(ir1::Irrep, ir2::Irrep) = Irreps(ir1) + Irreps(ir2)
+
+function Base.div(irs::Irreps, n::Integer)
+    Irreps(tuple([MulIrrep(div(mul_ir.mul, mul_ir.ir), mul_ir.ir)
+                  for mul_ir in irs]...))
+end
+function Base.fld(irs::Irreps, n::Integer)
+    Irreps(tuple([MulIrrep(fld(mul_ir.mul, n), mul_ir.ir) for mul_ir in irs]...))
+end
+
+dim(irs::Irreps) = sum(dim, irs._irreps; init = 0)
+num_irreps(irs::Irreps) = sum(x -> x.mul, irs._irreps; init = 0)
+ls(irs::Irreps) = [mul_ir.ir.l for mul_ir in irs for _ in 1:mul_ir.mul]
+function lmax(irs::Irreps)
+    isempty(irs._irreps) ? -1 : maximum(x -> x.ir.l, irs._irreps)
+end
+function count(ir::Irrep, irs::Irreps)
+    sum(mul_ir.mul for mul_ir in irs if mul_ir.ir == ir; init = 0)
+end
+function isscalar(irs::Irreps)
+    all(mul_ir -> mul_ir.ir.l == 0 && mul_ir.ir.p == 1, irs._irreps)
+end
+function mul_gcd(irs::Irreps)
+    isempty(irs._irreps) ? 0 : gcd(x -> x.mul, irs._irreps)
+end
+
+function Base.repeat(irs::Irreps, n::Integer)
+    Irreps(tuple(repeat(collect(irs._irreps), n)...))
+end
+
+function remove_zero_multiplicities(irs::Irreps)
+    Irreps(tuple([mul_ir for mul_ir in irs if mul_ir.mul > 0]...))
+end
+
+function unify(irs::Irreps)
+    out = Vector{MulIrrep}{eltype(irs)}()
+    for mul_ir in irs
+        if !isempty(out) && last(out).ir == mul_ir.ir
+            out[end] = MulIrrep(last(out).mul + mul_ir.mul, mul_ir.ir)
+        else
+            push!(out, mul_ir)
+        end
+    end
+    return Irreps(tuple(out...))
+end
+
+simplify(irs::Irreps) = unify(remove_zero_multiplicities(irs))
+
+function regroup(irs::Irreps)
+    counts = Dict{Irrep, Int}()
+    for mul_ir in irs
+        counts[mul_ir.ir] = get(counts, mul_ir.ir, 0) + mul_ir.mul
+    end
+    sorted_irreps = sort(collect(keys(counts)))
+    return Irreps(tuple([MulIrrep(counts[ir], ir)
+                         for ir in sorted_irreps if counts[ir] > 0]...))
+end
+
+function Base.sort(irs::Irreps)
+    v = collect(irs._irreps)
+    p = sortperm(v)
+    inv_p = invperm(p)
+    return Irreps(tuple(v[p]...))
+end
+
+function Base.filter(f::Function, irs::Irreps)
+    return Irreps(tuple([x for x in irs if f(x)]...))
+end
diff --git a/src/O3/irreps1.jl b/src/O3/irreps1.jl
new file mode 100644
index 0000000..ca80c09
--- /dev/null
+++ b/src/O3/irreps1.jl
@@ -0,0 +1,138 @@
+
+"""
+    Irreps
+
+Direct sum of irreducible representations of ``O(3)``.
+
+This struct does not contain any data, it is a structure that describes the representation.
+
+# Examples
+```jldoctest
+# Create a representation of 100 l=0 of even parity and 50 pseudo-vectors.
+julia> x = Irreps([(100, (0, 1)), (50, (1, 1))])
+100x0e+50x1e
+
+julia> dim(x)
+250
+
+# Create a representation of 100 l=0 of even parity and 50 pseudo-vectors.
+julia> Irreps("100x0e + 50x1e")
+100x0e+50x1e
+
+julia> Irreps("100x0e + 50x1e + 0x2e")
+100x0e+50x1e+0x2e
+
+julia> Irreps("100x0e + 50x1e + 0x2e") |> lmax
+1
+
+julia> Irrep("2e") in Irreps("0e + 2e")
+true
+
+# Empty Irreps
+julia> Irreps(), Irreps("")
+(, )
+```
+"""
+struct Irreps
+    irreps::Tuple{Vararg{MulIrrep}}
+end
+
+Irreps(ir::Irrep{T}) where {T <: Integer} = Irreps((MulIrrep(1, ir),))
+Irreps(mulir::MulIrrep{T}) where {T <: Integer} = Irreps((mulir,))
+function Irreps(mulir::Tuple{Int, Irrep{T}}) where {T <: Integer}
+    Irreps((MulIrrep(mulir...),))
+end
+Irreps() = Irreps{Int}(tuple())
+
+function Base.show(io::IO, irreps::Irreps)
+    isempty(irreps.irreps) && return print(io, "")
+    print(io, join(string.(irreps.irreps), " + "))
+end
+
+Base.iterate(irreps::Irreps, args...) = iterate(irreps.irreps, args...)
+Base.getindex(irreps::Irreps, i) = getindex(irreps.irreps, i)
+Base.length(irreps::Irreps) = length(irreps.irreps)
+Base.in(ir::Irrep, irreps::Irreps) = any(mulir -> mulir.ir == ir, irreps.irreps)
+
+function Base.:+(a::Irreps, b::Irreps)
+    return regroup(Irreps((a.irreps..., b.irreps...)))
+end
+
+Base.:+(a::Irrep, b::Irrep) = Irreps(a) + Irreps(b)
+
+function Base.:*(n::Integer, irreps::Irreps)
+    @assert n >= 0 "Multiplicity must be non-negative."
+    return Irreps(tuple([(n * mul, ir) for (mul, ir) in irreps]...))
+end
+
+function Base.:*(n::Integer, ir::Irrep)
+    @assert n >= 0 "Multiplicity must be non-negative."
+    return Irreps(((n, ir),))
+end
+
+Base.:*(ir::Irrep, n::Integer) = n * ir
+
+function Base.div(irreps::Irreps, n::Integer)
+    Irreps(tuple([(div(mul, n), ir) for (mul, ir) in irreps]...))
+end
+function Base.repeat(irreps::Irreps, n::Integer)
+    Irreps(tuple(repeat(collect(irreps.irreps), n)...))
+end
+
+function Base.show(io::IO, irreps::Irreps)
+    print(io, join(["$(mul)x$(ir)" for (mul, ir) in irreps.irreps], " + "))
+end
+
+dim(irreps::Irreps) = sum(mul * dim(ir) for (mul, ir) in irreps; init = 0)
+num_irreps(irreps::Irreps) = sum(first, irreps.irreps; init = 0)
+ls(irreps::Irreps) = [ir.l for (mul, ir) in irreps for _ in 1:mul]
+function lmax(irreps::Irreps)
+    isempty(irreps.irreps) ? -1 : maximum(ir.l for (_, ir) in irreps)
+end
+function count(ir::Irrep, irreps::Irreps)
+    sum(mul for (mul, irr) in irreps if irr == ir; init = 0)
+end
+function mul_gcd(irreps::Irreps)
+    isempty(irreps.irreps) ? 0 : gcd(first(x) for x in irreps)
+end
+isscalar(irreps::Irreps) = all(ir.l == 0 && ir.p == 1 for (_, ir) in irreps)
+
+function remove_zero_multiplicities(irreps::Irreps)
+    Irreps(tuple([(mul, ir) for (mul, ir) in irreps if mul > 0]...))
+end
+
+function unify(irreps::Irreps)
+    out = Vector{Tuple{Int, Irrep}}()
+    for (mul, ir) in irreps
+        if !isempty(out) && last(out)[2] == ir
+            out[end] = (last(out)[1] + mul, ir)
+        else
+            push!(out, (mul, ir))
+        end
+    end
+    return Irreps(tuple(out...))
+end
+
+simplify(irreps::Irreps) = unify(remove_zero_multiplicities(irreps))
+
+function regroup(irreps::Irreps)
+    counts = Dict{Irrep, Int}()
+    for (mul, ir) in irreps
+        counts[ir] = get(counts, ir, 0) + mul
+    end
+    sorted_irreps = sort(collect(keys(counts)))
+    return Irreps(tuple([(counts[ir], ir)
+                         for ir in sorted_irreps if counts[ir] > 0]...))
+end
+
+function _sort_with_p_info(irreps::Irreps)
+    v = collect(irreps.irreps)
+    p = sortperm(v, by = x -> x[2])
+    inv_p = invperm(p)
+    sorted_irreps = Irreps(tuple(v[p]...))
+    return (irreps = sorted_irreps, p = p, inv = inv_p)
+end
+
+function Base.sort(irreps::Irreps)
+    return _sort_with_p_info(irreps)[irreps]
+end
diff --git a/src/o3/rotations.jl b/src/O3/rotations.jl
similarity index 85%
rename from src/o3/rotations.jl
rename to src/O3/rotations.jl
index 71b656c..afaad00 100644
--- a/src/o3/rotations.jl
+++ b/src/O3/rotations.jl
@@ -28,11 +28,14 @@ function (::Type{Q})(; α::Real, β::Real, γ::Real) where {Q <: Quaternion}
     return QuatRotation(RotYXY(promote(α, β, γ)...)).q |> Q
 end
 
-function euler_angles(x::T) where {T <: Union{QuatRotation, RotMatrix3, AngleAxis}}
+function euler_angles(x::T) where {T <:
+                                   Union{QuatRotation, RotMatrix3, AngleAxis}}
     Rotations.params(RotYXY(x))
 end
 
-euler_angles(q::Q) where {Q <: Quaternion} = euler_angles(QuatRotation(q, false))
+function euler_angles(q::Q) where {Q <: Quaternion}
+    euler_angles(QuatRotation(q, false))
+end
 
 function (::Type{R})(q::Quaternion) where {R <: RotMatrix3}
     QuatRotation(q) |> R
@@ -44,7 +47,7 @@ end
 
 function CartesianToSphericalAngles(x::SVector{3, T}) where {T <: Real}
 
-    # done in e3nn to remove NaNs
+    # done in E3NN to remove NaNs
     # need to check for Julia
     normalize!(x, 2)
     clamp!(x, -1, 1)
diff --git a/src/o3/s2grid.jl b/src/O3/s2grid.jl
similarity index 97%
rename from src/o3/s2grid.jl
rename to src/O3/s2grid.jl
index 28ccabf..df4e206 100644
--- a/src/o3/s2grid.jl
+++ b/src/O3/s2grid.jl
@@ -47,7 +47,8 @@ function Base.show(io::IO, s::SphericalSignal)
 end
 
 function *(s::SphericalSignal, scalar::Number)
-    SphericalSignal(s.grid_values * scalar, s.quadrature, p_val = s.p_val, p_arg = s.p_arg)
+    SphericalSignal(
+        s.grid_values * scalar, s.quadrature, p_val = s.p_val, p_arg = s.p_arg)
 end
 *(scalar::Number, s::SphericalSignal) = s * scalar
 /(s::SphericalSignal, scalar::Number) = s * (1 / scalar)
diff --git a/src/o3/spherical_harmonics.jl b/src/O3/spherical_harmonics.jl
similarity index 86%
rename from src/o3/spherical_harmonics.jl
rename to src/O3/spherical_harmonics.jl
index c4a0bcb..a744c50 100644
--- a/src/o3/spherical_harmonics.jl
+++ b/src/O3/spherical_harmonics.jl
@@ -8,7 +8,7 @@ struct SphericalHarmonics
     _is_range_lmax::Bool
 
     function SphericalHarmonics(
-            irreps_out::Union{AbstractVector{Int}, o3.Irreps},
+            irreps_out::Union{AbstractVector{Int}, O3.Irreps},
             normalize::Bool,
             normalization::String,
             irreps_in
@@ -18,15 +18,15 @@ struct SphericalHarmonics
         if irreps_out isa Irreps && isnothing(irreps_in)
             for (mul, (l, p)) in irreps_out
                 if l % 2 == 1 & p == 1
-                    irreps_in = o3.Irreps("1e")
+                    irreps_in = O3.Irreps("1e")
                 end
             end
         end
         if isnothing(irreps_in)
-            irreps_in = o3.Irreps("1o")
+            irreps_in = O3.Irreps("1o")
         end
-        irreps_in = o3.Irreps(irreps_in)
-        if !(irreps_in in (o3.Irreps("1x1o"), o3.Irreps("1x1e")))
+        irreps_in = O3.Irreps(irreps_in)
+        if !(irreps_in in (O3.Irreps("1x1o"), O3.Irreps("1x1e")))
             throw(
                 ArgumentError(
                 "irreps_in for SphericalHarmonics must be either a vector (`1x1o`) or a pseudovector (`1x1e`), not `$(irreps_in)`",
@@ -34,7 +34,7 @@ struct SphericalHarmonics
             )
         end
         input_p = irreps_in[1].irrep.p
-        if irreps_out isa o3.Irreps
+        if irreps_out isa O3.Irreps
             ls = []
             for (mul, (l, p)) in irreps_out
                 if p != input_p^l
@@ -49,15 +49,16 @@ struct SphericalHarmonics
         else
             ls = Vector(irreps_out)
         end
-        irreps_out = o3.Irreps([(1, (l, input_p^l)) for l in ls])
-        irreps_out = o3.simplify(irreps_out)
+        irreps_out = O3.Irreps([(1, (l, input_p^l)) for l in ls])
+        irreps_out = O3.simplify(irreps_out)
         lmax = maximum(ls)
         _lmax = 11
         if lmax > _lmax
             throw(ErrorException("spherical_harmonics maximum l implemented is $_lmax."))
         end
 
-        return new(normalize, normalization, ls, lmax, ls == collect(range(0, maximum(ls))))
+        return new(normalize, normalization, ls, lmax,
+            ls == collect(range(0, maximum(ls))))
     end
 end
 
@@ -67,8 +68,9 @@ function SphericalHarmonics(
         normalization::String,
         irreps_in::Any
 )
-    irreps_out_rep = o3.Irreps(irreps_out)
-    return SphericalHarmonics(irreps_out_rep, normalize, normalization, irreps_in)
+    irreps_out_rep = O3.Irreps(irreps_out)
+    return SphericalHarmonics(
+        irreps_out_rep, normalize, normalization, irreps_in)
 end
 
 function (sh::SphericalHarmonics)(x::AbstractArray)
@@ -88,7 +90,8 @@ function (sh::SphericalHarmonics)(x::AbstractArray)
     if sh.normalization == "integral"
         sh_matrix /= sqrt(4 * pi)
     elseif sh.normalization == "norm"
-        sh_matrix /= cat([sqrt(2 * l + 1) * ones(2 * l + 1) for l in sh._ls_list])
+        sh_matrix /= cat([sqrt(2 * l + 1) * ones(2 * l + 1)
+                          for l in sh._ls_list])
     end
 
     return sh_matrix
@@ -123,7 +126,8 @@ function _spherical_harmonics(
 
     if lmax == 2
         return stack(
-            [sh_0_0, sh_1_0, sh_1_1, sh_1_2, sh_2_0, sh_2_1, sh_2_2, sh_2_3, sh_2_4],
+            [sh_0_0, sh_1_0, sh_1_1, sh_1_2, sh_2_0,
+                sh_2_1, sh_2_2, sh_2_3, sh_2_4],
             dims = ndims(sh_0_0) + 1
         )
     end
@@ -653,7 +657,8 @@ function _spherical_harmonics(
         (1 / 48) *
         sqrt(2) *
         (
-            sqrt(1122) * sh_7_10 * y - sqrt(102) * sh_7_11 * z - sqrt(102) * sh_7_3 * x -
+            sqrt(1122) * sh_7_10 * y - sqrt(102) * sh_7_11 * z -
+            sqrt(102) * sh_7_3 * x -
             sqrt(561) * sh_7_5 * x + sqrt(561) * sh_7_9 * z
         )
     )
@@ -666,14 +671,16 @@ function _spherical_harmonics(
         (1 / 80) *
         sqrt(2) *
         (
-            -sqrt(85) * sh_7_1 * x + sqrt(2210) * sh_7_11 * z + sqrt(2210) * sh_7_12 * y -
+            -sqrt(85) * sh_7_1 * x + sqrt(2210) * sh_7_11 * z +
+            sqrt(2210) * sh_7_12 * y -
             sqrt(85) * sh_7_13 * z - sqrt(2210) * sh_7_3 * x
         )
     )
     sh_8_14 = (
         -1 / 240 * sqrt(510) * sh_7_0 * x +
         (1 / 240) * sqrt(46410) * sh_7_12 * z +
-        (1 / 60) * sqrt(1785) * sh_7_13 * y - 1 / 240 * sqrt(510) * sh_7_14 * z -
+        (1 / 60) * sqrt(1785) * sh_7_13 * y -
+        1 / 240 * sqrt(510) * sh_7_14 * z -
         1 / 240 * sqrt(46410) * sh_7_2 * x
     )
     sh_8_15 = (
@@ -772,7 +779,8 @@ function _spherical_harmonics(
     end
 
     sh_9_0 = (1 / 6) * sqrt(38) * (sh_8_0 * z + sh_8_16 * x)
-    sh_9_1 = (1 / 9) * sqrt(19) * (sh_8_0 * y + 2 * sh_8_1 * z + 2 * sh_8_15 * x)
+    sh_9_1 = (1 / 9) * sqrt(19) *
+             (sh_8_0 * y + 2 * sh_8_1 * z + 2 * sh_8_15 * x)
     sh_9_2 = (
         -1 / 306 * sqrt(646) * sh_8_0 * z +
         (4 / 153) * sqrt(646) * sh_8_1 * y +
@@ -788,73 +796,88 @@ function _spherical_harmonics(
         (1 / 306) * sqrt(67830) * sh_8_3 * z
     )
     sh_9_4 = (
-        (1 / 306) * sqrt(58786) * sh_8_12 * x + (1 / 153) * sqrt(969) * sh_8_14 * x -
+        (1 / 306) * sqrt(58786) * sh_8_12 * x +
+        (1 / 153) * sqrt(969) * sh_8_14 * x -
         1 / 153 * sqrt(969) * sh_8_2 * z +
         (2 / 153) * sqrt(4522) * sh_8_3 * y +
         (1 / 306) * sqrt(58786) * sh_8_4 * z
     )
     sh_9_5 = (
-        (1 / 153) * sqrt(12597) * sh_8_11 * x + (1 / 153) * sqrt(1615) * sh_8_13 * x -
+        (1 / 153) * sqrt(12597) * sh_8_11 * x +
+        (1 / 153) * sqrt(1615) * sh_8_13 * x -
         1 / 153 * sqrt(1615) * sh_8_3 * z +
         (1 / 153) * sqrt(20995) * sh_8_4 * y +
         (1 / 153) * sqrt(12597) * sh_8_5 * z
     )
     sh_9_6 = (
-        (1 / 153) * sqrt(10659) * sh_8_10 * x + (1 / 306) * sqrt(9690) * sh_8_12 * x -
+        (1 / 153) * sqrt(10659) * sh_8_10 * x +
+        (1 / 306) * sqrt(9690) * sh_8_12 * x -
         1 / 306 * sqrt(9690) * sh_8_4 * z +
         (2 / 51) * sqrt(646) * sh_8_5 * y +
         (1 / 153) * sqrt(10659) * sh_8_6 * z
     )
     sh_9_7 = (
-        (1 / 306) * sqrt(13566) * sh_8_11 * x - 1 / 306 * sqrt(13566) * sh_8_5 * z +
+        (1 / 306) * sqrt(13566) * sh_8_11 * x -
+        1 / 306 * sqrt(13566) * sh_8_5 * z +
         (1 / 153) * sqrt(24871) * sh_8_6 * y +
         (1 / 306) * sqrt(35530) * sh_8_7 * z +
         (1 / 306) * sqrt(35530) * sh_8_9 * x
     )
     sh_9_8 = (
-        (1 / 153) * sqrt(4522) * sh_8_10 * x - 1 / 153 * sqrt(4522) * sh_8_6 * z +
+        (1 / 153) * sqrt(4522) * sh_8_10 * x -
+        1 / 153 * sqrt(4522) * sh_8_6 * z +
         (4 / 153) * sqrt(1615) * sh_8_7 * y +
         (1 / 51) * sqrt(1615) * sh_8_8 * x
     )
-    sh_9_9 = (1 / 51) * sqrt(323) * (-2 * sh_8_7 * x + 3 * sh_8_8 * y - 2 * sh_8_9 * z)
+    sh_9_9 = (1 / 51) * sqrt(323) *
+             (-2 * sh_8_7 * x + 3 * sh_8_8 * y - 2 * sh_8_9 * z)
     sh_9_10 = (
-        -1 / 153 * sqrt(4522) * sh_8_10 * z - 1 / 153 * sqrt(4522) * sh_8_6 * x +
+        -1 / 153 * sqrt(4522) * sh_8_10 * z -
+        1 / 153 * sqrt(4522) * sh_8_6 * x +
         (1 / 51) * sqrt(1615) * sh_8_8 * z +
         (4 / 153) * sqrt(1615) * sh_8_9 * y
     )
     sh_9_11 = (
-        (1 / 153) * sqrt(24871) * sh_8_10 * y - 1 / 306 * sqrt(13566) * sh_8_11 * z -
-        1 / 306 * sqrt(13566) * sh_8_5 * x - 1 / 306 * sqrt(35530) * sh_8_7 * x +
+        (1 / 153) * sqrt(24871) * sh_8_10 * y -
+        1 / 306 * sqrt(13566) * sh_8_11 * z -
+        1 / 306 * sqrt(13566) * sh_8_5 * x -
+        1 / 306 * sqrt(35530) * sh_8_7 * x +
         (1 / 306) * sqrt(35530) * sh_8_9 * z
     )
     sh_9_12 = (
-        (1 / 153) * sqrt(10659) * sh_8_10 * z + (2 / 51) * sqrt(646) * sh_8_11 * y -
+        (1 / 153) * sqrt(10659) * sh_8_10 * z +
+        (2 / 51) * sqrt(646) * sh_8_11 * y -
         1 / 306 * sqrt(9690) * sh_8_12 * z - 1 / 306 * sqrt(9690) * sh_8_4 * x -
         1 / 153 * sqrt(10659) * sh_8_6 * x
     )
     sh_9_13 = (
-        (1 / 153) * sqrt(12597) * sh_8_11 * z + (1 / 153) * sqrt(20995) * sh_8_12 * y -
+        (1 / 153) * sqrt(12597) * sh_8_11 * z +
+        (1 / 153) * sqrt(20995) * sh_8_12 * y -
         1 / 153 * sqrt(1615) * sh_8_13 * z - 1 / 153 * sqrt(1615) * sh_8_3 * x -
         1 / 153 * sqrt(12597) * sh_8_5 * x
     )
     sh_9_14 = (
-        (1 / 306) * sqrt(58786) * sh_8_12 * z + (2 / 153) * sqrt(4522) * sh_8_13 * y -
+        (1 / 306) * sqrt(58786) * sh_8_12 * z +
+        (2 / 153) * sqrt(4522) * sh_8_13 * y -
         1 / 153 * sqrt(969) * sh_8_14 * z - 1 / 153 * sqrt(969) * sh_8_2 * x -
         1 / 306 * sqrt(58786) * sh_8_4 * x
     )
     sh_9_15 = (
         -1 / 306 * sqrt(1938) * sh_8_1 * x +
         (1 / 306) * sqrt(67830) * sh_8_13 * z +
-        (1 / 51) * sqrt(1615) * sh_8_14 * y - 1 / 306 * sqrt(1938) * sh_8_15 * z -
+        (1 / 51) * sqrt(1615) * sh_8_14 * y -
+        1 / 306 * sqrt(1938) * sh_8_15 * z -
         1 / 306 * sqrt(67830) * sh_8_3 * x
     )
     sh_9_16 = (
         -1 / 306 * sqrt(646) * sh_8_0 * x +
         (2 / 153) * sqrt(4845) * sh_8_14 * z +
-        (4 / 153) * sqrt(646) * sh_8_15 * y - 1 / 306 * sqrt(646) * sh_8_16 * z -
+        (4 / 153) * sqrt(646) * sh_8_15 * y -
+        1 / 306 * sqrt(646) * sh_8_16 * z -
         2 / 153 * sqrt(4845) * sh_8_2 * x
     )
-    sh_9_17 = (1 / 9) * sqrt(19) * (-2 * sh_8_1 * x + 2 * sh_8_15 * z + sh_8_16 * y)
+    sh_9_17 = (1 / 9) * sqrt(19) *
+              (-2 * sh_8_1 * x + 2 * sh_8_15 * z + sh_8_16 * y)
     sh_9_18 = (1 / 6) * sqrt(38) * (-sh_8_0 * x + sh_8_16 * z)
     if lmax == 9
         return stack(
@@ -985,31 +1008,36 @@ function _spherical_harmonics(
         (1 / 95) * sqrt(6783) * sh_9_3 * z
     )
     sh_10_4 = (
-        (3 / 95) * sqrt(665) * sh_9_14 * x + (3 / 190) * sqrt(133) * sh_9_16 * x -
+        (3 / 95) * sqrt(665) * sh_9_14 * x +
+        (3 / 190) * sqrt(133) * sh_9_16 * x -
         3 / 190 * sqrt(133) * sh_9_2 * z +
         (4 / 95) * sqrt(399) * sh_9_3 * y +
         (3 / 95) * sqrt(665) * sh_9_4 * z
     )
     sh_10_5 = (
-        (21 / 380) * sqrt(190) * sh_9_13 * x + (1 / 190) * sqrt(1995) * sh_9_15 * x -
+        (21 / 380) * sqrt(190) * sh_9_13 * x +
+        (1 / 190) * sqrt(1995) * sh_9_15 * x -
         1 / 190 * sqrt(1995) * sh_9_3 * z +
         (3 / 38) * sqrt(133) * sh_9_4 * y +
         (21 / 380) * sqrt(190) * sh_9_5 * z
     )
     sh_10_6 = (
-        (7 / 380) * sqrt(1482) * sh_9_12 * x + (3 / 380) * sqrt(1330) * sh_9_14 * x -
+        (7 / 380) * sqrt(1482) * sh_9_12 * x +
+        (3 / 380) * sqrt(1330) * sh_9_14 * x -
         3 / 380 * sqrt(1330) * sh_9_4 * z +
         (21 / 95) * sqrt(19) * sh_9_5 * y +
         (7 / 380) * sqrt(1482) * sh_9_6 * z
     )
     sh_10_7 = (
-        (3 / 190) * sqrt(1729) * sh_9_11 * x + (21 / 380) * sqrt(38) * sh_9_13 * x -
+        (3 / 190) * sqrt(1729) * sh_9_11 * x +
+        (21 / 380) * sqrt(38) * sh_9_13 * x -
         21 / 380 * sqrt(38) * sh_9_5 * z +
         (7 / 190) * sqrt(741) * sh_9_6 * y +
         (3 / 190) * sqrt(1729) * sh_9_7 * z
     )
     sh_10_8 = (
-        (3 / 190) * sqrt(1463) * sh_9_10 * x + (7 / 190) * sqrt(114) * sh_9_12 * x -
+        (3 / 190) * sqrt(1463) * sh_9_10 * x +
+        (7 / 190) * sqrt(114) * sh_9_12 * x -
         7 / 190 * sqrt(114) * sh_9_6 * z +
         (6 / 95) * sqrt(266) * sh_9_7 * y +
         (3 / 190) * sqrt(1463) * sh_9_8 * z
@@ -1020,42 +1048,50 @@ function _spherical_harmonics(
         (1 / 190) * sqrt(21945) * sh_9_9 * x
     )
     sh_10_10 = (
-        -3 / 190 * sqrt(1995) * sh_9_10 * z - 3 / 190 * sqrt(1995) * sh_9_8 * x +
+        -3 / 190 * sqrt(1995) * sh_9_10 * z -
+        3 / 190 * sqrt(1995) * sh_9_8 * x +
         (1 / 19) * sqrt(399) * sh_9_9 * y
     )
     sh_10_11 = (
-        (3 / 190) * sqrt(4389) * sh_9_10 * y - 3 / 190 * sqrt(798) * sh_9_11 * z -
+        (3 / 190) * sqrt(4389) * sh_9_10 * y -
+        3 / 190 * sqrt(798) * sh_9_11 * z -
         3 / 190 * sqrt(798) * sh_9_7 * x + (1 / 190) * sqrt(21945) * sh_9_9 * z
     )
     sh_10_12 = (
-        (3 / 190) * sqrt(1463) * sh_9_10 * z + (6 / 95) * sqrt(266) * sh_9_11 * y -
+        (3 / 190) * sqrt(1463) * sh_9_10 * z +
+        (6 / 95) * sqrt(266) * sh_9_11 * y -
         7 / 190 * sqrt(114) * sh_9_12 * z - 7 / 190 * sqrt(114) * sh_9_6 * x -
         3 / 190 * sqrt(1463) * sh_9_8 * x
     )
     sh_10_13 = (
-        (3 / 190) * sqrt(1729) * sh_9_11 * z + (7 / 190) * sqrt(741) * sh_9_12 * y -
+        (3 / 190) * sqrt(1729) * sh_9_11 * z +
+        (7 / 190) * sqrt(741) * sh_9_12 * y -
         21 / 380 * sqrt(38) * sh_9_13 * z - 21 / 380 * sqrt(38) * sh_9_5 * x -
         3 / 190 * sqrt(1729) * sh_9_7 * x
     )
     sh_10_14 = (
-        (7 / 380) * sqrt(1482) * sh_9_12 * z + (21 / 95) * sqrt(19) * sh_9_13 * y -
+        (7 / 380) * sqrt(1482) * sh_9_12 * z +
+        (21 / 95) * sqrt(19) * sh_9_13 * y -
         3 / 380 * sqrt(1330) * sh_9_14 * z - 3 / 380 * sqrt(1330) * sh_9_4 * x -
         7 / 380 * sqrt(1482) * sh_9_6 * x
     )
     sh_10_15 = (
-        (21 / 380) * sqrt(190) * sh_9_13 * z + (3 / 38) * sqrt(133) * sh_9_14 * y -
+        (21 / 380) * sqrt(190) * sh_9_13 * z +
+        (3 / 38) * sqrt(133) * sh_9_14 * y -
         1 / 190 * sqrt(1995) * sh_9_15 * z - 1 / 190 * sqrt(1995) * sh_9_3 * x -
         21 / 380 * sqrt(190) * sh_9_5 * x
     )
     sh_10_16 = (
-        (3 / 95) * sqrt(665) * sh_9_14 * z + (4 / 95) * sqrt(399) * sh_9_15 * y -
+        (3 / 95) * sqrt(665) * sh_9_14 * z +
+        (4 / 95) * sqrt(399) * sh_9_15 * y -
         3 / 190 * sqrt(133) * sh_9_16 * z - 3 / 190 * sqrt(133) * sh_9_2 * x -
         3 / 95 * sqrt(665) * sh_9_4 * x
     )
     sh_10_17 = (
         -3 / 380 * sqrt(266) * sh_9_1 * x +
         (1 / 95) * sqrt(6783) * sh_9_15 * z +
-        (3 / 190) * sqrt(2261) * sh_9_16 * y - 3 / 380 * sqrt(266) * sh_9_17 * z -
+        (3 / 190) * sqrt(2261) * sh_9_16 * y -
+        3 / 380 * sqrt(266) * sh_9_17 * z -
         1 / 95 * sqrt(6783) * sh_9_3 * x
     )
     sh_10_18 = (
@@ -1220,93 +1256,114 @@ function _spherical_harmonics(
         (1 / 154) * sqrt(18354) * sh_10_3 * z
     )
     sh_11_4 = (
-        (1 / 154) * sqrt(16422) * sh_10_16 * x + (1 / 77) * sqrt(161) * sh_10_18 * x -
+        (1 / 154) * sqrt(16422) * sh_10_16 * x +
+        (1 / 77) * sqrt(161) * sh_10_18 * x -
         1 / 77 * sqrt(161) * sh_10_2 * z +
         (2 / 77) * sqrt(966) * sh_10_3 * y +
         (1 / 154) * sqrt(16422) * sh_10_4 * z
     )
     sh_11_5 = (
-        (2 / 231) * sqrt(8211) * sh_10_15 * x + (1 / 231) * sqrt(2415) * sh_10_17 * x -
+        (2 / 231) * sqrt(8211) * sh_10_15 * x +
+        (1 / 231) * sqrt(2415) * sh_10_17 * x -
         1 / 231 * sqrt(2415) * sh_10_3 * z +
         (1 / 231) * sqrt(41055) * sh_10_4 * y +
         (2 / 231) * sqrt(8211) * sh_10_5 * z
     )
     sh_11_6 = (
-        (2 / 77) * sqrt(805) * sh_10_14 * x + (1 / 154) * sqrt(1610) * sh_10_16 * x -
+        (2 / 77) * sqrt(805) * sh_10_14 * x +
+        (1 / 154) * sqrt(1610) * sh_10_16 * x -
         1 / 154 * sqrt(1610) * sh_10_4 * z +
         (4 / 77) * sqrt(322) * sh_10_5 * y +
         (2 / 77) * sqrt(805) * sh_10_6 * z
     )
     sh_11_7 = (
-        (1 / 22) * sqrt(230) * sh_10_13 * x + (1 / 22) * sqrt(46) * sh_10_15 * x -
+        (1 / 22) * sqrt(230) * sh_10_13 * x +
+        (1 / 22) * sqrt(46) * sh_10_15 * x -
         1 / 22 * sqrt(46) * sh_10_5 * z +
         (1 / 11) * sqrt(115) * sh_10_6 * y +
         (1 / 22) * sqrt(230) * sh_10_7 * z
     )
     sh_11_8 = (
-        (1 / 66) * sqrt(1794) * sh_10_12 * x + (1 / 33) * sqrt(138) * sh_10_14 * x -
+        (1 / 66) * sqrt(1794) * sh_10_12 * x +
+        (1 / 33) * sqrt(138) * sh_10_14 * x -
         1 / 33 * sqrt(138) * sh_10_6 * z +
         (4 / 33) * sqrt(69) * sh_10_7 * y +
         (1 / 66) * sqrt(1794) * sh_10_8 * z
     )
     sh_11_9 = (
-        (1 / 77) * sqrt(2093) * sh_10_11 * x + (1 / 77) * sqrt(966) * sh_10_13 * x -
+        (1 / 77) * sqrt(2093) * sh_10_11 * x +
+        (1 / 77) * sqrt(966) * sh_10_13 * x -
         1 / 77 * sqrt(966) * sh_10_7 * z +
         (1 / 77) * sqrt(6279) * sh_10_8 * y +
         (1 / 77) * sqrt(2093) * sh_10_9 * z
     )
     sh_11_10 = (
-        (1 / 77) * sqrt(3542) * sh_10_10 * x + (1 / 154) * sqrt(4830) * sh_10_12 * x -
-        1 / 154 * sqrt(4830) * sh_10_8 * z + (2 / 77) * sqrt(1610) * sh_10_9 * y
+        (1 / 77) * sqrt(3542) * sh_10_10 * x +
+        (1 / 154) * sqrt(4830) * sh_10_12 * x -
+        1 / 154 * sqrt(4830) * sh_10_8 * z +
+        (2 / 77) * sqrt(1610) * sh_10_9 * y
     )
     sh_11_11 = (
-        (1 / 21) * sqrt(483) * sh_10_10 * y - 1 / 231 * sqrt(26565) * sh_10_11 * z -
+        (1 / 21) * sqrt(483) * sh_10_10 * y -
+        1 / 231 * sqrt(26565) * sh_10_11 * z -
         1 / 231 * sqrt(26565) * sh_10_9 * x
     )
     sh_11_12 = (
-        (1 / 77) * sqrt(3542) * sh_10_10 * z + (2 / 77) * sqrt(1610) * sh_10_11 * y -
-        1 / 154 * sqrt(4830) * sh_10_12 * z - 1 / 154 * sqrt(4830) * sh_10_8 * x
+        (1 / 77) * sqrt(3542) * sh_10_10 * z +
+        (2 / 77) * sqrt(1610) * sh_10_11 * y -
+        1 / 154 * sqrt(4830) * sh_10_12 * z -
+        1 / 154 * sqrt(4830) * sh_10_8 * x
     )
     sh_11_13 = (
-        (1 / 77) * sqrt(2093) * sh_10_11 * z + (1 / 77) * sqrt(6279) * sh_10_12 * y -
+        (1 / 77) * sqrt(2093) * sh_10_11 * z +
+        (1 / 77) * sqrt(6279) * sh_10_12 * y -
         1 / 77 * sqrt(966) * sh_10_13 * z - 1 / 77 * sqrt(966) * sh_10_7 * x -
         1 / 77 * sqrt(2093) * sh_10_9 * x
     )
     sh_11_14 = (
-        (1 / 66) * sqrt(1794) * sh_10_12 * z + (4 / 33) * sqrt(69) * sh_10_13 * y -
+        (1 / 66) * sqrt(1794) * sh_10_12 * z +
+        (4 / 33) * sqrt(69) * sh_10_13 * y -
         1 / 33 * sqrt(138) * sh_10_14 * z - 1 / 33 * sqrt(138) * sh_10_6 * x -
         1 / 66 * sqrt(1794) * sh_10_8 * x
     )
     sh_11_15 = (
-        (1 / 22) * sqrt(230) * sh_10_13 * z + (1 / 11) * sqrt(115) * sh_10_14 * y -
+        (1 / 22) * sqrt(230) * sh_10_13 * z +
+        (1 / 11) * sqrt(115) * sh_10_14 * y -
         1 / 22 * sqrt(46) * sh_10_15 * z - 1 / 22 * sqrt(46) * sh_10_5 * x -
         1 / 22 * sqrt(230) * sh_10_7 * x
     )
     sh_11_16 = (
-        (2 / 77) * sqrt(805) * sh_10_14 * z + (4 / 77) * sqrt(322) * sh_10_15 * y -
-        1 / 154 * sqrt(1610) * sh_10_16 * z - 1 / 154 * sqrt(1610) * sh_10_4 * x -
+        (2 / 77) * sqrt(805) * sh_10_14 * z +
+        (4 / 77) * sqrt(322) * sh_10_15 * y -
+        1 / 154 * sqrt(1610) * sh_10_16 * z -
+        1 / 154 * sqrt(1610) * sh_10_4 * x -
         2 / 77 * sqrt(805) * sh_10_6 * x
     )
     sh_11_17 = (
-        (2 / 231) * sqrt(8211) * sh_10_15 * z + (1 / 231) * sqrt(41055) * sh_10_16 * y -
-        1 / 231 * sqrt(2415) * sh_10_17 * z - 1 / 231 * sqrt(2415) * sh_10_3 * x -
+        (2 / 231) * sqrt(8211) * sh_10_15 * z +
+        (1 / 231) * sqrt(41055) * sh_10_16 * y -
+        1 / 231 * sqrt(2415) * sh_10_17 * z -
+        1 / 231 * sqrt(2415) * sh_10_3 * x -
         2 / 231 * sqrt(8211) * sh_10_5 * x
     )
     sh_11_18 = (
-        (1 / 154) * sqrt(16422) * sh_10_16 * z + (2 / 77) * sqrt(966) * sh_10_17 * y -
+        (1 / 154) * sqrt(16422) * sh_10_16 * z +
+        (2 / 77) * sqrt(966) * sh_10_17 * y -
         1 / 77 * sqrt(161) * sh_10_18 * z - 1 / 77 * sqrt(161) * sh_10_2 * x -
         1 / 154 * sqrt(16422) * sh_10_4 * x
     )
     sh_11_19 = (
         -1 / 154 * sqrt(322) * sh_10_1 * x +
         (1 / 154) * sqrt(18354) * sh_10_17 * z +
-        (1 / 77) * sqrt(3059) * sh_10_18 * y - 1 / 154 * sqrt(322) * sh_10_19 * z -
+        (1 / 77) * sqrt(3059) * sh_10_18 * y -
+        1 / 154 * sqrt(322) * sh_10_19 * z -
         1 / 154 * sqrt(18354) * sh_10_3 * x
     )
     sh_11_20 = (
         -1 / 462 * sqrt(966) * sh_10_0 * x +
         (1 / 231) * sqrt(45885) * sh_10_18 * z +
-        (2 / 231) * sqrt(4830) * sh_10_19 * y - 1 / 231 * sqrt(45885) * sh_10_2 * x -
+        (2 / 231) * sqrt(4830) * sh_10_19 * y -
+        1 / 231 * sqrt(45885) * sh_10_2 * x -
         1 / 462 * sqrt(966) * sh_10_20 * z
     )
     sh_11_21 = (
@@ -1486,13 +1543,15 @@ function _spherical_harmonics(
         (5 / 276) * sqrt(2415) * sh_11_3 * z
     )
     sh_12_4 = (
-        (5 / 276) * sqrt(2185) * sh_11_18 * x - 5 / 276 * sqrt(69) * sh_11_2 * z +
+        (5 / 276) * sqrt(2185) * sh_11_18 * x -
+        5 / 276 * sqrt(69) * sh_11_2 * z +
         (5 / 276) * sqrt(69) * sh_11_20 * x +
         (5 / 69) * sqrt(115) * sh_11_3 * y +
         (5 / 276) * sqrt(2185) * sh_11_4 * z
     )
     sh_12_5 = (
-        (5 / 184) * sqrt(874) * sh_11_17 * x + (5 / 276) * sqrt(115) * sh_11_19 * x -
+        (5 / 184) * sqrt(874) * sh_11_17 * x +
+        (5 / 276) * sqrt(115) * sh_11_19 * x -
         5 / 276 * sqrt(115) * sh_11_3 * z +
         (5 / 276) * sqrt(2185) * sh_11_4 * y +
         (5 / 184) * sqrt(874) * sh_11_5 * z
@@ -1501,25 +1560,29 @@ function _spherical_harmonics(
         (5 / 552) *
         sqrt(3) *
         (
-            sqrt(2346) * sh_11_16 * x + sqrt(230) * sh_11_18 * x - sqrt(230) * sh_11_4 * z +
+            sqrt(2346) * sh_11_16 * x + sqrt(230) * sh_11_18 * x -
+            sqrt(230) * sh_11_4 * z +
             12 * sqrt(23) * sh_11_5 * y +
             sqrt(2346) * sh_11_6 * z
         )
     )
     sh_12_7 = (
-        (5 / 138) * sqrt(391) * sh_11_15 * x + (5 / 552) * sqrt(966) * sh_11_17 * x -
+        (5 / 138) * sqrt(391) * sh_11_15 * x +
+        (5 / 552) * sqrt(966) * sh_11_17 * x -
         5 / 552 * sqrt(966) * sh_11_5 * z +
         (5 / 276) * sqrt(2737) * sh_11_6 * y +
         (5 / 138) * sqrt(391) * sh_11_7 * z
     )
     sh_12_8 = (
-        (5 / 138) * sqrt(345) * sh_11_14 * x + (5 / 276) * sqrt(322) * sh_11_16 * x -
+        (5 / 138) * sqrt(345) * sh_11_14 * x +
+        (5 / 276) * sqrt(322) * sh_11_16 * x -
         5 / 276 * sqrt(322) * sh_11_6 * z +
         (10 / 69) * sqrt(46) * sh_11_7 * y +
         (5 / 138) * sqrt(345) * sh_11_8 * z
     )
     sh_12_9 = (
-        (5 / 552) * sqrt(4830) * sh_11_13 * x + (5 / 92) * sqrt(46) * sh_11_15 * x -
+        (5 / 552) * sqrt(4830) * sh_11_13 * x +
+        (5 / 92) * sqrt(46) * sh_11_15 * x -
         5 / 92 * sqrt(46) * sh_11_7 * z +
         (5 / 92) * sqrt(345) * sh_11_8 * y +
         (5 / 552) * sqrt(4830) * sh_11_9 * z
@@ -1527,40 +1590,49 @@ function _spherical_harmonics(
     sh_12_10 = (
         (5 / 552) * sqrt(4186) * sh_11_10 * z +
         (5 / 552) * sqrt(4186) * sh_11_12 * x +
-        (5 / 184) * sqrt(230) * sh_11_14 * x - 5 / 184 * sqrt(230) * sh_11_8 * z +
+        (5 / 184) * sqrt(230) * sh_11_14 * x -
+        5 / 184 * sqrt(230) * sh_11_8 * z +
         (5 / 138) * sqrt(805) * sh_11_9 * y
     )
     sh_12_11 = (
         (5 / 276) * sqrt(3289) * sh_11_10 * y +
         (5 / 276) * sqrt(1794) * sh_11_11 * x +
-        (5 / 552) * sqrt(2530) * sh_11_13 * x - 5 / 552 * sqrt(2530) * sh_11_9 * z
+        (5 / 552) * sqrt(2530) * sh_11_13 * x -
+        5 / 552 * sqrt(2530) * sh_11_9 * z
     )
     sh_12_12 = (
-        -5 / 276 * sqrt(1518) * sh_11_10 * x + (5 / 23) * sqrt(23) * sh_11_11 * y -
+        -5 / 276 * sqrt(1518) * sh_11_10 * x +
+        (5 / 23) * sqrt(23) * sh_11_11 * y -
         5 / 276 * sqrt(1518) * sh_11_12 * z
     )
     sh_12_13 = (
-        (5 / 276) * sqrt(1794) * sh_11_11 * z + (5 / 276) * sqrt(3289) * sh_11_12 * y -
-        5 / 552 * sqrt(2530) * sh_11_13 * z - 5 / 552 * sqrt(2530) * sh_11_9 * x
+        (5 / 276) * sqrt(1794) * sh_11_11 * z +
+        (5 / 276) * sqrt(3289) * sh_11_12 * y -
+        5 / 552 * sqrt(2530) * sh_11_13 * z -
+        5 / 552 * sqrt(2530) * sh_11_9 * x
     )
     sh_12_14 = (
         -5 / 552 * sqrt(4186) * sh_11_10 * x +
         (5 / 552) * sqrt(4186) * sh_11_12 * z +
-        (5 / 138) * sqrt(805) * sh_11_13 * y - 5 / 184 * sqrt(230) * sh_11_14 * z -
+        (5 / 138) * sqrt(805) * sh_11_13 * y -
+        5 / 184 * sqrt(230) * sh_11_14 * z -
         5 / 184 * sqrt(230) * sh_11_8 * x
     )
     sh_12_15 = (
-        (5 / 552) * sqrt(4830) * sh_11_13 * z + (5 / 92) * sqrt(345) * sh_11_14 * y -
+        (5 / 552) * sqrt(4830) * sh_11_13 * z +
+        (5 / 92) * sqrt(345) * sh_11_14 * y -
         5 / 92 * sqrt(46) * sh_11_15 * z - 5 / 92 * sqrt(46) * sh_11_7 * x -
         5 / 552 * sqrt(4830) * sh_11_9 * x
     )
     sh_12_16 = (
-        (5 / 138) * sqrt(345) * sh_11_14 * z + (10 / 69) * sqrt(46) * sh_11_15 * y -
+        (5 / 138) * sqrt(345) * sh_11_14 * z +
+        (10 / 69) * sqrt(46) * sh_11_15 * y -
         5 / 276 * sqrt(322) * sh_11_16 * z - 5 / 276 * sqrt(322) * sh_11_6 * x -
         5 / 138 * sqrt(345) * sh_11_8 * x
     )
     sh_12_17 = (
-        (5 / 138) * sqrt(391) * sh_11_15 * z + (5 / 276) * sqrt(2737) * sh_11_16 * y -
+        (5 / 138) * sqrt(391) * sh_11_15 * z +
+        (5 / 276) * sqrt(2737) * sh_11_16 * y -
         5 / 552 * sqrt(966) * sh_11_17 * z - 5 / 552 * sqrt(966) * sh_11_5 * x -
         5 / 138 * sqrt(391) * sh_11_7 * x
     )
@@ -1574,25 +1646,30 @@ function _spherical_harmonics(
         )
     )
     sh_12_19 = (
-        (5 / 184) * sqrt(874) * sh_11_17 * z + (5 / 276) * sqrt(2185) * sh_11_18 * y -
+        (5 / 184) * sqrt(874) * sh_11_17 * z +
+        (5 / 276) * sqrt(2185) * sh_11_18 * y -
         5 / 276 * sqrt(115) * sh_11_19 * z - 5 / 276 * sqrt(115) * sh_11_3 * x -
         5 / 184 * sqrt(874) * sh_11_5 * x
     )
     sh_12_20 = (
-        (5 / 276) * sqrt(2185) * sh_11_18 * z + (5 / 69) * sqrt(115) * sh_11_19 * y -
+        (5 / 276) * sqrt(2185) * sh_11_18 * z +
+        (5 / 69) * sqrt(115) * sh_11_19 * y -
         5 / 276 * sqrt(69) * sh_11_2 * x - 5 / 276 * sqrt(69) * sh_11_20 * z -
         5 / 276 * sqrt(2185) * sh_11_4 * x
     )
     sh_12_21 = (
         -5 / 552 * sqrt(138) * sh_11_1 * x +
         (5 / 276) * sqrt(2415) * sh_11_19 * z +
-        (5 / 92) * sqrt(161) * sh_11_20 * y - 5 / 552 * sqrt(138) * sh_11_21 * z -
+        (5 / 92) * sqrt(161) * sh_11_20 * y -
+        5 / 552 * sqrt(138) * sh_11_21 * z -
         5 / 276 * sqrt(2415) * sh_11_3 * x
     )
     sh_12_22 = (
-        -5 / 552 * sqrt(46) * sh_11_0 * x - 5 / 552 * sqrt(10626) * sh_11_2 * x +
+        -5 / 552 * sqrt(46) * sh_11_0 * x -
+        5 / 552 * sqrt(10626) * sh_11_2 * x +
         (5 / 552) * sqrt(10626) * sh_11_20 * z +
-        (5 / 138) * sqrt(253) * sh_11_21 * y - 5 / 552 * sqrt(46) * sh_11_22 * z
+        (5 / 138) * sqrt(253) * sh_11_21 * y -
+        5 / 552 * sqrt(46) * sh_11_22 * z
     )
     sh_12_23 = -5 / 24 * sqrt(22) * sh_11_1 * x +
                (5 / 24) * sqrt(22) * sh_11_21 * z +
diff --git a/src/o3/wigner.jl b/src/O3/wigner.jl
similarity index 89%
rename from src/o3/wigner.jl
rename to src/O3/wigner.jl
index b3027ae..970695a 100644
--- a/src/o3/wigner.jl
+++ b/src/O3/wigner.jl
@@ -25,8 +25,10 @@ function su2_clebsch_gordan(j1::T, j2::T, j3::T) where {T <: Real}
         for m1 in [x / 2 for x in (-Int(2 * j1)):2:Int(2 * j1)]
             for m2 in [x / 2 for x in (-Int(2 * j2)):2:Int(2 * j2)]
                 if abs(m1 + m2) <= j3
-                    val = _su2_cg(promote(j1, m1), promote(j2, m2), promote(j3, m1 + m2))
-                    mat[Int(j1 + m1 + 1), Int(j2 + m2 + 1), Int(j3 + m1 + m2 + 1)] = val
+                    val = _su2_cg(
+                        promote(j1, m1), promote(j2, m2), promote(j3, m1 + m2))
+                    mat[Int(j1 + m1 + 1), Int(j2 + m2 + 1),
+                        Int(j3 + m1 + m2 + 1)] = val
                 end
             end
         end
@@ -76,7 +78,8 @@ function clebsch_gordan(l1::Int, l2::Int, l3::Int)
     Q3 = change_basis_real_to_complex(l3)
     result = zeros(ComplexF64, size(Q1, 2), size(Q2, 2), size(Q3, 2))
     Q3_conj_T = conj(transpose(Q3))
-    @tensor result[j, l, m] := Q1[i, j] * Q2[k, l] * Q3_conj_T[m, n] * C[i, k, n]
+    @tensor result[
+        j, l, m] := Q1[i, j] * Q2[k, l] * Q3_conj_T[m, n] * C[i, k, n]
     @assert all(abs.(imag.(C)) .< 1e-5)
     return real.(result)
 end
@@ -121,4 +124,6 @@ function wigner_D(l::Int, angles::SVector{3, T}) where {T <: Real}
     return wigner_D(l, angles[1], angles[2], angles[3])
 end
 
-Broadcast.broadcast(::typeof(wigner_D), l, α, β, γ) = broadcast(wigner_D, Ref(l), α, β, γ)
+function Broadcast.broadcast(::typeof(wigner_D), l, α, β, γ)
+    broadcast(wigner_D, Ref(l), α, β, γ)
+end
diff --git a/src/e3nn.jl b/src/e3nn.jl
deleted file mode 100644
index f25c53a..0000000
--- a/src/e3nn.jl
+++ /dev/null
@@ -1,6 +0,0 @@
-module e3nn
-
-# o3
-include("o3/o3.jl")
-
-end
diff --git a/src/o3/irreps.jl b/src/o3/irreps.jl
deleted file mode 100644
index a436953..0000000
--- a/src/o3/irreps.jl
+++ /dev/null
@@ -1,252 +0,0 @@
-using Random: GLOBAL_RNG, AbstractRNG
-
-struct Irrep
-    l::Int
-    p::Int
-
-    function Irrep(l::Int, p::Int)
-        (l >= 0) || throw(ArgumentError("l must be zero or positive integer, got $l"))
-        (p in (-1, 1)) || throw(ArgumentError("parity(p) must be one of (-1, 1) got $p"))
-        return new(l, p)
-    end
-end
-
-function Irrep(l::T) where {T <: AbstractString}
-    name = strip(l)
-    try
-        l = parse(Int, name[1:(end - 1)])
-        (l >= 0) || throw(ArgumentError("l must be zero or positive integer, got $l"))
-        p = Dict('e' => 1, 'o' => -1, 'y' => (-1)^l)[name[end]]
-        return Irrep(l, p)
-    catch
-        throw(ArgumentError("Cannot convert string $name to an Irrep"))
-    end
-end
-
-function Irrep(l::Tuple)
-    @assert length(l) == 2
-    l, p = l
-    return Irrep(l, p)
-end
-
-Irrep(l::Irrep) = l
-
-dim(x::Irrep) = 2 * x.l + 1
-
-isscalar(x::Irrep) = (x.l == 0) && (x.p == 1)
-
-Base.iterate(x::Irrep, args...) = iterate((x.l, x.p), args...)
-
-function Base.:*(x1::Irrep, x2::Irrep)
-    p = x1.p * x2.p
-    lmin = abs(x1.l - x2.l)
-    lmax = x1.l + x2.l
-    return (Irrep(l, p) for l in lmin:lmax)
-end
-
-function Base.:*(i::Int, x::Irrep)
-    return Irreps([(i, x)])
-end
-
-function Base.:+(x1::Irrep, x2::Irrep)
-    return Irreps(x1) + Irreps(x2)
-end
-
-Base.isless(x1::Irrep, x2::Irrep) = Base.isless((x1.l, x1.p), (x2.l, x2.p))
-
-function Base.show(io::IO, x::Irrep)
-    p = Dict(+1 => "e", -1 => "o")[x.p]
-    s = "$(x.l)$p"
-    print(io, s)
-end
-
-function D_from_angles(x::Irrep, α, β, γ, k)
-    throw(error("Not Implemnted yet"))
-end
-
-function D_from_quaternion(x::Irrep, q, k)
-    throw(error("Not Implemnted yet"))
-end
-
-function D_from_matrix(x::Irrep, R)
-    throw(error("Not Implemnted yet"))
-end
-
-function D_from_matrix(x::Irrep, axis, angle)
-    throw(error("Not Implemnted yet"))
-end
-
-struct MulIrrep
-    mul::Int
-    irrep::Irrep
-end
-
-function Base.show(io::IO, mx::MulIrrep)
-    s = "$(mx.mul)x$(mx.irrep)"
-    print(io, s)
-end
-
-dim(mx::MulIrrep) = mx.mul * dim(mx.irrep)
-
-Base.convert(::Type{Tuple}, x::MulIrrep) = (x.mul, x.irrep)
-struct Irreps
-    irreps::Vector{MulIrrep}
-end
-
-Irreps(irrep::Irrep) = Irreps([MulIrrep(1, irrep)])
-Irreps(irreps::Irreps) = irreps
-
-function Irreps(irreps::T) where {T <: AbstractString}
-    mulirreps = MulIrrep[]
-    if strip(irreps) != ""
-        for mul_irrep in split(irreps, "+")
-            if occursin("x", mul_irrep)
-                mul, irrep = split(mul_irrep, "x")
-                mul = parse(Int, mul)
-                (mul >= 0) || throw(ArgumentError("mul should be greater than 0 got $mul"))
-                irrep = Irrep(irrep)
-            else
-                mul = 1
-                irrep = Irrep(mul_irrep)
-            end
-            push!(mulirreps, MulIrrep(mul, irrep))
-        end
-    end
-    return Irreps(mulirreps)
-end
-
-# fallback
-function Irreps(irreps::AbstractVector)
-    mul = nothing
-    irrep = nothing
-    out = MulIrrep[]
-    for mul_irrep in irreps
-        if typeof(mul_irrep) <: AbstractString
-            # fix case for mixture of irrep and irrpes in a vector
-            # if occursin("+", mul_irrep)
-            #     irreps = Irreps(mul_irrep)
-            # end
-            mul = 1
-            irrep = Irrep(mul_irrep)
-        elseif mul_irrep isa Irrep
-            mul = 1
-            irrep = mul_irrep
-        elseif mul_irrep isa MulIrrep
-            mul, irrep = mul_irrep
-        elseif length(mul_irrep) == 2
-            mul, irrep = mul_irrep
-            irrep = Irrep(irrep)
-        else
-            throw(ArgumentError("Unable to interpret $(mul_irrep) as an irrep."))
-        end
-        push!(out, MulIrrep(mul, irrep))
-    end
-    return Irreps(out)
-end
-
-function Base.show(io::IO, xs::Irreps)
-    s = join(["$(mul_irrep)" for mul_irrep in xs.irreps], "+")
-    print(io, s)
-end
-
-Base.:(==)(x1::Irreps, x2::Irreps) = x1.irreps == x2.irreps
-Base.:+(x1::Irreps, x2::Irreps) = Irreps(vcat(x1.irreps, x2.irreps))
-Base.:*(i::Int, xs::Irreps) = Irreps(repeat(xs.irreps, i))
-Base.:*(xs::Irreps, i::Int) = Base.:*(i::Int, xs::Irreps)
-
-Base.length(xs::Irreps) = length(xs.irreps)
-
-Base.getindex(xs::Irreps, i::Int) = xs.irreps[i]
-Base.getindex(xs::Irreps, idx::AbstractRange) = xs.irreps[idx] |> Irreps
-Base.firstindex(xs::Irreps) = Base.firstindex(xs.irreps)
-Base.lastindex(xs::Irreps) = Base.lastindex(xs.irreps)
-
-Base.in(x::Irrep, xs::Irreps) = x ∈ [mx.irrep for mx in xs.irreps]
-function Base.count(x::Irrep, xs::Irreps)
-    sum([mx.mul for mx in xs.irreps if mx.irrep == x], init = 0)
-end
-function Base.iterate(xs::Irreps, state = 1)
-    state > length(xs) ? nothing :
-    ((xs[state].mul, (xs[state].irrep.l, xs[state].irrep.p)), state + 1)
-end
-
-"""
-Representation of spherical harmonics.
-"""
-function spherical_harmonics(lmax::Int, p::Int = -1)::Irreps
-    return Irreps([(1, (l, p^l)) for l in 1:lmax])
-end
-
-function Base.randn(
-        xs::Irreps,
-        dims,
-        normalization::String,
-        rng::AbstractRNG,
-        ::Type{T}
-) where {T <: Number}
-    di = dims[end]
-    lsize = dims[:di]
-    rsize = dims[di + 1 :]
-
-    if normalization == "component"
-        return randn(rng, T, lsize..., dim(xs), rsize...)
-        # implement the norm
-    end
-end
-
-"""
-Simplify the representaions.
-"""
-function simplify(xs::Irreps)::Irreps
-    out = []
-    for (mul, irrep) in xs
-        if length(out) != 0 && out[end][2] == irrep
-            out[end] = (out[end][1] + mul, irrep)
-        elseif mul > 0
-            push!(out, (mul, irrep))
-        end
-    end
-    return Irreps(out)
-end
-
-"""
-Remove any irreps with multiplicities of zero.
-"""
-remove_zero_multiplicities(xs::Irreps) = [(mul, irreps) for (mul, irreps) in xs if mul > 0] |>
-                                         Irreps
-
-function Base.sort(xs::Irreps)::Irreps
-    out = [(mx.irrep, i, mx.mul) for (i, mx) in enumerate(xs)]
-    out = Base.sort(out)
-    sorted_irreps = Irreps([(mul, irrep) for (irrep, _, mul) in out])
-    return sorted_irreps
-end
-
-dim(xs::Irreps) = sum([mx.mul * dim(mx.irrep) for mx in xs], init = 0)
-
-num_irreps(xs::Irreps) = sum([mx.mul for mx in xs], init = 0)
-
-ls(xs::Irreps) = [mx.irrep.l for mx in xs for _ in 1:(mx.mul)]
-
-function lmax(xs::Irreps)::Int
-    if length(xs) == 0
-        throw(ArgumentError("Cannot get lmax of empty Irreps"))
-    end
-    return maximum(ls(xs))
-end
-
-function D_from_angles(xs::Irreps, α, β, γ, k)
-    throw(error("Not Implemnted yet"))
-end
-
-function D_from_quaternion(xs::Irreps, q, k)
-    throw(error("Not Implemnted yet"))
-end
-
-function D_from_matrix(xs::Irreps, R)
-    throw(error("Not Implemnted yet"))
-end
-
-function D_from_axis_angle(xs::Irreps, axis, angle)
-    throw(error("Not Implemnted yet"))
-end
diff --git a/src/o3/o3.jl b/src/o3/o3.jl
deleted file mode 100644
index d800d94..0000000
--- a/src/o3/o3.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-module o3
-
-include("irreps.jl")
-export Irrep, Irreps, dim, isscalar, spherical_harmonics, simplify
-export remove_zero_multiplicities, num_irreps, ls, lmax
-# not implemented yet
-export D_from_angles, D_from_quaternion, D_from_axis_angle, D_from_matrix
-
-include("wigner.jl")
-export clebsch_gordan, so3_generators, su2_generators, wigner_D
-
-include("rotations.jl")
-using .rot
-export euler_angles, CartesianToSphericalAngles, SphercialAnglesToCartesian
-
-include("spherical_harmonics.jl")
-export spherical_harmonics, SphericalHarmonics
-
-include("s2grid.jl")
-using .S2Grid
-
-end
diff --git a/test/o3/irreps.jl b/test/O3/irreps.jl
similarity index 83%
rename from test/o3/irreps.jl
rename to test/O3/irreps.jl
index 4e81636..7e1da3d 100644
--- a/test/o3/irreps.jl
+++ b/test/O3/irreps.jl
@@ -1,4 +1,4 @@
-using e3nn.o3
+using E3NN.O3
 using Test
 
 @testset "Irreps" begin
@@ -37,14 +37,16 @@ using Test
 
     @testset "arithmetic" begin
         @test 3 * Irrep("6o") == Irreps("3x6o")
-        products = [ir for ir in Irrep("1o") * Irrep("2e")]
-        @test products == [Irrep("1o"), Irrep("2o"), Irrep("3o")]
+        # We are considering product of irrep as tensor_product
+        # products = [ir for ir in Irrep("1o") * Irrep("2e")]
+        # @test products == [Irrep("1o"), Irrep("2o"), Irrep("3o")]
 
         @test Irrep("4o") + Irrep("7e") == Irreps("4o + 7e")
-        @test Irreps("1o + 4o") + Irreps("1o + 7e") == Irreps("1o + 4o + 1o + 7e")
+        @test Irreps("1o + 4o") + Irreps("1o + 7e") ==
+              Irreps("1o + 4o + 1o + 7e")
 
-        @test 2 * Irreps("2x2e + 4x1o") == Irreps("2x2e + 4x1o + 2x2e + 4x1o")
-        @test Irreps("2x2e + 4x1o") * 2 == Irreps("2x2e + 4x1o + 2x2e + 4x1o")
+        @test 2 * Irreps("2x2e + 4x1o") == Irreps("4x2e + 8x1o")
+        @test Irreps("2x2e + 4x1o") * 2 == Irreps("4x2e + 8x1o")
     end
 
     @testset "empty" begin
diff --git a/test/o3/wigner.jl b/test/o3/wigner.jl
deleted file mode 100644
index 9d04ce1..0000000
--- a/test/o3/wigner.jl
+++ /dev/null
@@ -1,13 +0,0 @@
-using e3nn.o3
-using Test
-using Rotations
-using MLUtils: batch
-
-@testset "WignerD" begin
-    @testset "basic" begin
-        R = rand(RotYXY, 10)
-        angles = R .|> Rotations.params
-        D = wigner_D.(1, angles)
-        @test batch(R - D) .|> abs |> maximum < 1e-8
-    end
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index eafee2d..62bb449 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,6 @@
-using e3nn
+using E3NN
 using Test
 
-@testset "e3nn.jl" begin
-    include("o3/irreps.jl")
+@testset "E3NN.jl" begin
+    include("O3/irreps.jl")
 end