Changing Cell

[jaemolihm]

Hi, thanks for this nice package!

I would like to ask what do you think about making following changes to the Cell struct.

1. Allow magmom to be nothing (nonmagnetic), Vector{T} (collinear), and Vector{MVector{3,M}} (noncollinear).
2. Use Vector{MVector{3,P}} instead of MMatrix{3,N,P} for position. The latter causes type instability when constructing Cell. (This was also suggested in Fix type-instability in lattice-field of Cell #84)

[singularitti]

1. Hi @jaemolihm, I have never used magmom. Does the spglib C-API support non-collinear argument? Could you give me an example?
2. Could you tell me how MMatrix{3,N,P} makes it type unstable? An example will be helpful. Thank you!

[jaemolihm]

1. Yes it does. The following is an example, where setting magmoms reduces the symmetry of the system.

import spglib
lattice = [[1,0,0],[0,1,0],[0,0,1]]
positions = [[-0.1, -0.1, -0.1], [0.1, 0.1, 0.1]]
atoms = [1, 1]
magmoms = [[0., 0., 1.], [0., 0., -1.]]
>>> spglib.get_symmetry((lattice, positions, atoms))["rotations"].shape
(12, 3, 3)
>>> spglib.get_symmetry((lattice, positions, atoms, magmoms))["rotations"].shape
(2, 3, 3)

The python function is
https://github.com/spglib/spglib/blob/develop/python/spglib/spglib.py#L163-L300
It calls the C function spgms_get_symmetry_with_site_tensors:
https://github.com/spglib/spglib/blob/6bc78f9f1caf346fdd9e35e9eb1873f841d71303/src/spglib.c#L465-L472

2. It's because N, the number of atoms, is part of the type and it is not known at compile time.
   Example:

using StaticArrays

struct MyCell1{N,P}
    positions::MMatrix{3,N,P}
end

function MyCell1(positions_::AbstractVector)
    positions = reduce(hcat, positions_)
    N, P = size(positions, 2), eltype(positions)
    return MyCell1{N,P}(positions)
end

struct MyCell2{P}
    positions::Vector{MVector{3, P}}
end

function MyCell2(positions_::AbstractVector)
    positions = MVector{3}.(positions_)
    P = eltype(eltype(positions))
    return MyCell2{P}(positions)
end

pos = [[0., 0., 0.], [0.5, 0.5, 0.5]]

@code_warntype MyCell1(pos)
@code_warntype MyCell2(pos)

Example output:

julia> @code_warntype MyCell1(pos)
MethodInstance for MyCell1(::Vector{Vector{Float64}})
  from MyCell1(positions_::AbstractVector) in Main at REPL[3]:1
Arguments
  #self#::Type{MyCell1}
  positions_::Vector{Vector{Float64}}
Locals
  P::Type{Float64}
  N::Int64
  positions::Matrix{Float64}
Body::MyCell1{_A, Float64} where _A
1 ─      (positions = Main.reduce(Main.hcat, positions_))
│   %2 = Main.size(positions, 2)::Int64
│   %3 = Main.eltype(positions)::Core.Const(Float64)
│        (N = %2)
│        (P = %3)
│   %6 = Core.apply_type(Main.MyCell1, N, P::Core.Const(Float64))::Type{MyCell1{_A, Float64}} where _A
│   %7 = (%6)(positions)::MyCell1{_A, Float64} where _A
└──      return %7


julia> @code_warntype MyCell2(pos)
MethodInstance for MyCell2(::Vector{Vector{Float64}})
  from MyCell2(positions_::AbstractVector) in Main at REPL[5]:1
Arguments
  #self#::Type{MyCell2}
  positions_::Vector{Vector{Float64}}
Locals
  P::Type{Float64}
  positions::Vector{MVector{3, Float64}}
Body::MyCell2{Float64}
1 ─ %1 = Core.apply_type(Main.MVector, 3)::Core.Const(MVector{3})
│   %2 = Base.broadcasted(%1, positions_)::Core.PartialStruct(Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, Type{MVector{3}}, Tuple{Vector{Vector{Float64}}}}, Any[Core.Const(MVector{3}), Tuple{Vector{Vector{Float64}}}, Core.Const(nothing)])
│        (positions = Base.materialize(%2))
│   %4 = Main.eltype(positions)::Core.Const(MVector{3, Float64})
│        (P = Main.eltype(%4))
│   %6 = Core.apply_type(Main.MyCell2, P::Core.Const(Float64))::Core.Const(MyCell2{Float64})
│   %7 = (%6)(positions)::MyCell2{Float64}
└──      return %7

[singularitti]

Hi @jaemolihm, I removed the constraint on N in Cell in version v0.6.0. For the first request, I still need a little time.

[singularitti]

1 should be fixed in recent updates. I still have to test, though.

[singularitti]

Hi @jaemolihm, the magnetic code is now available on the main branch's head. Please try.