Implemented gen-biquad interaction in dipole mode

[Hao-Phys]

The implementation passes the check for the scalar biquadratic interactions.
However, a unit test for general interactions still needs to be included. I cannot find a simple general biquadratic interaction compatible with a given symmetry.

[Lazersmoke]

Putting this code here so people have it available:

The following applies to Sunny.diamond_crystal():

Bond(1, 1, [0, 0, 1])
Distance 1, coordination 6
Connects [0, 0, 0] to [0, 0, 1]
Allowed exchange matrix:[A ⋅  ⋅ -I  B ⋅  ⋅ ⋅
                         ⋅ A  ⋅  B -I ⋅  ⋅ ⋅
                         ⋅ ⋅  C  ⋅  ⋅ ⋅  D J
                         I B  ⋅  H  ⋅ ⋅  ⋅ ⋅
                         B I  ⋅  ⋅  H ⋅  ⋅ ⋅
                         ⋅ ⋅  ⋅  ⋅  ⋅ E  ⋅ ⋅
                         ⋅ ⋅  D  ⋅  ⋅ ⋅  F K
                         ⋅ ⋅ -J  ⋅  ⋅ ⋅ -K G]
function physical_basis_su3()
    S = 1
    dipole_operators = spin_matrices(; N=2S+1)
    Sx, Sy, Sz = dipole_operators

    # we choose a particular basis of the quadrupole operators listed in Appendix B of *Phys. Rev. B 104, 104409*
    quadrupole_operators = Vector{Matrix{ComplexF64}}(undef, 5)
    quadrupole_operators[1] = -(Sx * Sz + Sz * Sx)
    quadrupole_operators[2] = -(Sy * Sz + Sz * Sy)
    quadrupole_operators[3] = Sx * Sx - Sy * Sy
    quadrupole_operators[4] = Sx * Sy + Sy * Sx
    quadrupole_operators[5] = √3 * Sz * Sz - 1/√3 * S * (S+1) * I

    [dipole_operators; quadrupole_operators]
end

but note that the 3x5 off-diagonal block of the exchange matrix is disallowed by time-reversal symmetry due to having odd number of spin operators. via Kip, this is the code needed to implement this:

S = spin_irrep_label(sys, 1)
Sx, Sy, Sz = spin_matrices(S)
Q = [-(Sx * Sz + Sz * Sx),
     -(Sy * Sz + Sz * Sy),
     Sx * Sx - Sy * Sy,
     Sx * Sy + Sy * Sx,
     √3 * Sz * Sz - 1/√3 * S * (S+1) * I]
Qi, Qj = to_product_space(Q, Q)
biquad = [
    1.4 0 0 0 0
    0 1.4 0 0 0
    0 0 1.1 0 0
    0 0 0 1.2 +1.5
    0 0 0 -1.5 1.3
]
bond = Bond(1, 1, [0, 0, 1])
set_pair_coupling!(sys, 0.01(Qi'*biquad*Qj), bond)```

[kbarros]

I made a few tweaks to enhance readability. One notable thing is that it can be 10x slower to turn on quadrupoles (e.g., for the powder_averaging.jl tutorial), so I have it disabled unless a nonzero biquad coupling actually appears. In the future, we can precompute the rotated 5x5 coupling matrix which will speed things up a lot.