Making `supercell()` work for arbitary matrices

[brainandforce]

The current approach to supercell() seems to work correctly when a diagonal matrix/vector is used for scaling (and certainly for a scalar), but it doesn't seem to work for an arbitrary matrix. In my test case (MgZn₂) using a matrix that converts the cell to an orthorhombic supercell (using the C-centered primitive > conventional matrix) causes some atomic positions to collide and be removed.

[brainandforce]

I think I've fixed this one with #128

[brainandforce]

Turns out I haven't. Here's an example for MgCu₂:

julia> MgCu2 = AtomList(
           RealBasis(3.5 * SMatrix{3,3}(0, 1, 1, 1, 0, 1, 1, 1, 0)),
           [
               AtomPosition(12, 0, 0, 0),
               AtomPosition(12, 1/4, 1/4, 1/4),
               AtomPosition(29, 5/8, 5/8, 5/8),
               AtomPosition(29, 1/8, 5/8, 5/8),
               AtomPosition(29, 5/8, 1/8, 5/8),
               AtomPosition(29, 5/8, 5/8, 1/8)
           ]
       )
AtomList{3}:
  Atomic positions:
    12  Mg  [  0.000000  0.000000  0.000000 ]
    12  Mg  [  0.250000  0.250000  0.250000 ]
    29  Cu  [  0.625000  0.625000  0.625000 ]
    29  Cu  [  0.125000  0.625000  0.625000 ]
    29  Cu  [  0.625000  0.125000  0.625000 ]
    29  Cu  [  0.625000  0.625000  0.125000 ]
  defined in terms of basis vectors:
    a: [  0.000000  3.500000  3.500000 ]   (4.949747)
    b: [  3.500000  0.000000  3.500000 ]   (4.949747)
    c: [  3.500000  3.500000  0.000000 ]   (4.949747)

julia> transform = [-1 1 1; 1 -1 1; 1 1 -1] * diagm([2,1,1])
3×3 Matrix{Int64}:
 -2   1   1
  2  -1   1
  2   1  -1

julia> MgCu2_sc = supercell(MgCu2, transform)
┌ Warning: Atoms 13 and 31 are very close or identical!
│ Atom 13: AtomPosition{3}:
│   12  Mg  [  0.250000  0.000000  0.500000 ]
│ Atom 31: AtomPosition{3}:
│   12  Mg  [  0.250000  0.000000  0.500000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 14 and 32 are very close or identical!
│ Atom 14: AtomPosition{3}:
│   12  Mg  [  0.375000  0.250000  0.750000 ]
│ Atom 32: AtomPosition{3}:
│   12  Mg  [  0.375000  0.250000  0.750000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 15 and 33 are very close or identical!
│ Atom 15: AtomPosition{3}:
│   29  Cu  [  0.562500  0.625000  0.125000 ]
│ Atom 33: AtomPosition{3}:
│   29  Cu  [  0.562500  0.625000  0.125000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 16 and 34 are very close or identical!
│ Atom 16: AtomPosition{3}:
│   29  Cu  [  0.562500  0.375000  0.875000 ]
│ Atom 34: AtomPosition{3}:
│   29  Cu  [  0.562500  0.375000  0.875000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 17 and 35 are very close or identical!
│ Atom 17: AtomPosition{3}:
│   29  Cu  [  0.437500  0.625000  0.875000 ]
│ Atom 35: AtomPosition{3}:
│   29  Cu  [  0.437500  0.625000  0.875000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 18 and 36 are very close or identical!
│ Atom 18: AtomPosition{3}:
│   29  Cu  [  0.437500  0.375000  0.125000 ]
│ Atom 36: AtomPosition{3}:
│   29  Cu  [  0.437500  0.375000  0.125000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 19 and 25 are very close or identical!
│ Atom 19: AtomPosition{3}:
│   12  Mg  [  0.250000  0.500000  0.000000 ]
│ Atom 25: AtomPosition{3}:
│   12  Mg  [  0.250000  0.500000  0.000000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 20 and 26 are very close or identical!
│ Atom 20: AtomPosition{3}:
│   12  Mg  [  0.375000  0.750000  0.250000 ]
│ Atom 26: AtomPosition{3}:
│   12  Mg  [  0.375000  0.750000  0.250000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 21 and 27 are very close or identical!
│ Atom 21: AtomPosition{3}:
│   29  Cu  [  0.562500  0.125000  0.625000 ]
│ Atom 27: AtomPosition{3}:
│   29  Cu  [  0.562500  0.125000  0.625000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 22 and 28 are very close or identical!
│ Atom 22: AtomPosition{3}:
│   29  Cu  [  0.562500  0.875000  0.375000 ]
│ Atom 28: AtomPosition{3}:
│   29  Cu  [  0.562500  0.875000  0.375000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 23 and 29 are very close or identical!
│ Atom 23: AtomPosition{3}:
│   29  Cu  [  0.437500  0.125000  0.375000 ]
│ Atom 29: AtomPosition{3}:
│   29  Cu  [  0.437500  0.125000  0.375000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
┌ Warning: Atoms 24 and 30 are very close or identical!
│ Atom 24: AtomPosition{3}:
│   29  Cu  [  0.437500  0.875000  0.625000 ]
│ Atom 30: AtomPosition{3}:
│   29  Cu  [  0.437500  0.875000  0.625000 ]
│ Distance (in supplied basis): 0.0
└ @ Xtal ~/git/Xtal.jl/src/atoms.jl:106
AtomList{3}:
  Atomic positions:
    12  Mg  [  0.000000  0.000000  0.000000 ]
    12  Mg  [  0.125000  0.250000  0.250000 ]
    29  Cu  [  0.312500  0.625000  0.625000 ]
    29  Cu  [  0.312500  0.375000  0.375000 ]
    29  Cu  [  0.187500  0.625000  0.375000 ]
    29  Cu  [  0.187500  0.375000  0.625000 ]
    12  Mg  [  0.000000  0.500000  0.500000 ]
    12  Mg  [  0.125000  0.750000  0.750000 ]
    29  Cu  [  0.312500  0.125000  0.125000 ]
    29  Cu  [  0.312500  0.875000  0.875000 ]
    29  Cu  [  0.187500  0.125000  0.875000 ]
    29  Cu  [  0.187500  0.875000  0.125000 ]
    12  Mg  [  0.250000  0.000000  0.500000 ]
    12  Mg  [  0.375000  0.250000  0.750000 ]
    29  Cu  [  0.562500  0.625000  0.125000 ]
    29  Cu  [  0.562500  0.375000  0.875000 ]
    29  Cu  [  0.437500  0.625000  0.875000 ]
    29  Cu  [  0.437500  0.375000  0.125000 ]
    12  Mg  [  0.250000  0.500000  0.000000 ]
    12  Mg  [  0.375000  0.750000  0.250000 ]
    29  Cu  [  0.562500  0.125000  0.625000 ]
    29  Cu  [  0.562500  0.875000  0.375000 ]
    29  Cu  [  0.437500  0.125000  0.375000 ]
    29  Cu  [  0.437500  0.875000  0.625000 ]
    12  Mg  [  0.500000  0.500000  0.500000 ]
    12  Mg  [  0.625000  0.750000  0.750000 ]
    29  Cu  [  0.812500  0.125000  0.125000 ]
    29  Cu  [  0.812500  0.875000  0.875000 ]
    29  Cu  [  0.687500  0.125000  0.875000 ]
    29  Cu  [  0.687500  0.875000  0.125000 ]
    12  Mg  [  0.500000  0.000000  0.000000 ]
    12  Mg  [  0.625000  0.250000  0.250000 ]
    29  Cu  [  0.812500  0.625000  0.625000 ]
    29  Cu  [  0.812500  0.375000  0.375000 ]
    29  Cu  [  0.687500  0.625000  0.375000 ]
    29  Cu  [  0.687500  0.375000  0.625000 ]
  defined in terms of basis vectors:
    a: [ 14.000000  0.000000  0.000000 ]   (14.000000)
    b: [  0.000000  7.000000  0.000000 ]   (7.000000)
    c: [  0.000000  0.000000  7.000000 ]   (7.000000)

It seems to be dropping a large number of atoms, leading to only a 6x increase in atom count rather than the expected 8x:

julia> det(transform)
8.0

julia> length(MgCu2)
6

julia> length(MgCu2_sc)
36

[brainandforce]

Of note: The matrix

julia> hcat([-1, 1, 1], [1, -1, 1], [1, 1, -1]) * diagm([2, 1, 1])
3×3 Matrix{Int64}:
 -2   1   1
  2  -1   1
  2   1  -1

shows this issue, but the matrices

julia> hcat([-1, 1, 1], [1, -1, 1], [1, 1, -1]) * diagm([1, 2, 1])
3×3 Matrix{Int64}:
 -1   2   1
  1  -2   1
  1   2  -1

and

julia> hcat([-1, 1, 1], [1, -1, 1], [1, 1, -1]) * diagm([1, 1, 2])
3×3 Matrix{Int64}:
 -1   1   2
  1  -1   2
  1   1  -2

do not.

[brainandforce]

I'm just noting this here: the solution to this might involve some kind of matrix decomposition that I'm not aware of yet. If we can decompose the problematic transform into a the product of the primitive>FCC matrix and the diagonal matrix, then we can execute the enlargement in two steps.

[brainandforce]

It might be best to use some variation of Bresenham's algorithm here to get all of the integer points contained within a supercell.

[brainandforce]

So I think I figured out the problem: The LU decomposition is distinctly different from the Hermite normal form decomposition I need.

[brainandforce]

Cross-referencing #131 here: this reference gives the algorithm for constructing a supercell: https://lan496.github.io/dsenum/supercell.html

The only issue is understanding what exactly it does. I get that the Smith normal form (SNF) is what provides the translations needed for constructing the supercell, but the issue is that the SNF requires that the factors exist in a canonical order, and that order in itself is not always sensible: if you are transforming a C-centered orthorhombic cell into its conventional setting, the factor of 2 will translate in the wrong direction if applied naively.

Therefore, the last part of the algorithm needs to permute the factors so that they copy the atoms along the correct dimension.