Custom generators of rotations in set_pair_coupling!

[kbarros]

To support custom spin representations on each site, the user must be able to specify the generators of rotation. For example, to work with explicit spin and orbital angular momenta, one should work in a tensor product space representation. But some models might also employ a direct sum of spin multiplets, and it's nice to give power users this level of control.

The proper, eventual solution will be a variant of Moment that can be passed to the system constructor. That might take a while, so this PR provides partial support, where the generators of rotation can be provided in the call to set_pair_coupling!. This new feature is not exposed in the public API.

What's correct in this PR:

• The symmetry consistency check on the pair coupling operator.
• Symmetry propagation of the coupling to other bonds.

What's missing and dangerous:

• No capability to specify a custom local magnetic moment dipole. In particular, this means that Zeeman coupling and dipole-dipole interactions will be wrong (don't try to use them). Also, spin-spin structure factor intensites will require the user to build a custom MeasureSpec object.
• Support only for set_pair_coupling!. Generators cannot yet be passed to set_onsite_coupling!.

Because of these limitations, there is an intentional choice to hide the new keyword argument gen from the public docs.

[kbarros]

Here's how this could look for two spin multiplets

bracket(A, B) = A*B - B*A

# First spin 1/2 doublet in direct sum space for one site
S1x = Hermitian(ComplexF64[0 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 0] / 2)
S1y = Hermitian(ComplexF64[0 0 0 0; 0 0 -im 0; 0 im 0 0; 0 0 0 0] / 2)
S1z = Hermitian(diagm(ComplexF64[0, 1, -1, 0]) / 2)
S1 = [S1x, S1y, S1z]
@assert bracket(S1x, S1y) ≈ im * S1z
@assert bracket(S1y, S1z) ≈ im * S1x
@assert bracket(S1z, S1x) ≈ im * S1y

# Second spin 1/2 doublet in direct sum space for one site
S2x = Hermitian(ComplexF64[0 0 0 1; 0 0 0 0; 0 0 0 0; 1 0 0 0] / 2)
S2y = Hermitian(ComplexF64[0 0 0 -im; 0 0 0 0; 0 0 0 0; im 0 0 0] / 2)
S2z = Hermitian(diagm(ComplexF64[1, 0, 0, -1]) / 2)
S2 = Hermitian.([S2x, S2y, S2z])

# Generator of rotation for both doublets on one site
gen = SVector{3}(Hermitian.(S1+S2))
@assert bracket(gen[1], gen[2]) ≈ im * gen[3]
@assert bracket(gen[2], gen[3]) ≈ im * gen[1]
@assert bracket(gen[3], gen[1]) ≈ im * gen[2]


# First doublet operator in product space of sites (i, j)
S1_i, S1_j = to_product_space(S1, S1)

# Second doublet operator in product space of sites (i, j)
S2_i, S2_j = to_product_space(S2, S2)

# Fake spin-3/2. We're actually using SU(4) to model the direct sum of two spin
# doublets. In principle, this requires specifying generator of rotation _and_
# the local magnetic dipole operator. With the Sunny#gen-spec branch, we support
# overriding the generator of rotation, but not yet the magnetic dipole.
sys = System(cryst, [1 => Moment(; s=3/2, g=2)], :SUN)

# "Usual" symmetry-allowed J also works for the bilinear coupling of the same
# doublet on two sites. It's crucial to specify the generator of rotation `gen`
# that acts on a direct sum of two doublets. Then the interactions will be
# correctly checked and propagated by spacegroup symmetries.
b = Bond(1, 1, [1, 0, 0])
J = [1 0 0; 0 3 0; 0 0 4]
set_pair_coupling!(sys, S1_i'*J*S1_j, b; gen) # OK
set_pair_coupling!(sys, S2_i'*J*S2_j, b; gen) # OK

# Conversely, symmetry-disallowed J will properly generate an error
J_bad = SA[0.0 1.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0]
set_pair_coupling!(sys, S1_i' * J_bad * S1_j, b; gen) # Error! Symmetry violating coupling

# Interestingly, this also fails the symmetry check. This failure is associated
# with symops that map bond to its reverse. It fails because doublets 1 and 2
# are distinct.
set_pair_coupling!(sys, S1_i'*J*S2_j, b; gen) # Error! Symmetry violating coupling

# For now, Zeeman coupling and dipole-dipole interactions will be WRONG. Also,
# to calculate structure factor, one will need a custom MeasureSpec with the
# correct dipole observable (probably you want to provide the dipole S1+S2 for
# each SU(4) site).