ACEsuit / ACE.jl

Parameterisation of Equivariant Properties of Particle Systems
66 stars 15 forks source link

Ms/equivarinat matrices #99

Closed MatthiasSachs closed 2 years ago

MatthiasSachs commented 2 years ago

All tests asserting the equivariance properties of EuclidianMatrix pass with tol = 1E-15.

The test println(@test(all(test_fio(basis; warntype = false)))) (see line 41 in test/test_EquivariantMatrix.jl) does NOT pass...

MatthiasSachs commented 2 years ago

test/test_EquivariantMatrix.jl not yet added to test collection. Thus, all tests are passing for now...

cortner commented 2 years ago

thank you - I'll take a look asap.

cortner commented 2 years ago

I'm looking into the read/write issue now.

cortner commented 2 years ago

Regarding the name, I think EquivariantMatrix is too general - it doesn't encode the specific symmetry. How about EuclideanMatrix? This would then also generalize nicely to EuclideanTensor if/when we get to that.

cortner commented 2 years ago

Another thing I've noticed is that you've based it off a fairly old commit. You'll need to rebase it onto the latest main please.

cortner commented 2 years ago

You'll see there that EuclideanVector has changed a bit and it would be good to model EuclideanMatrix on that. I hope that should require minimal work.

MatthiasSachs commented 2 years ago

You'll see there that EuclideanVector has changed a bit and it would be good to model EuclideanMatrix on that. I hope that should require minimal work.

This is mostly done. However, 1) I am not sure if the implementation of coco_init(::EuclideanMatrix{CT}) is correct? 2) Moreover, somehow there is an error thrown which relates to the conversion between SMatrix and EuclideanMatrix.

Would you mind taking a brief look at that (I think this might be much easier for you to resolve since you know the parts of the code well), or alternatively, we could meet at the beginning of next week?

cortner commented 2 years ago

Error message:

 [1] convert(T::Type, φ::EuclideanMatrix{Float64})
   @ ACE ~/gits/temp/ms.ACE.jl/src/properties.jl:28
 [2] push!(a::Vector{EuclideanMatrix{ComplexF64}}, item::EuclideanMatrix{Float64})
   @ Base ./array.jl:994
....

It happens in this line:

push!(vals, U[irow, icol])

The computed cc U[irow, icol] is of type item::EuclideanMatrix{Float64} and get's appended to a vector with eltype EuclideanMatrix{ComplexF64}. We need a convert implementation for this. THat's easy to do, but it concerns me. Why are we creating a complex-valued CC matrix but then fill it with real-valued CCs?

MatthiasSachs commented 2 years ago

Error message:

 [1] convert(T::Type, φ::EuclideanMatrix{Float64})
   @ ACE ~/gits/temp/ms.ACE.jl/src/properties.jl:28
 [2] push!(a::Vector{EuclideanMatrix{ComplexF64}}, item::EuclideanMatrix{Float64})
   @ Base ./array.jl:994
....

It happens in this line:

push!(vals, U[irow, icol])

The computed cc U[irow, icol] is of type item::EuclideanMatrix{Float64} and get's appended to a vector with eltype EuclideanMatrix{ComplexF64}. We need a convert implementation for this. THat's easy to do, but it concerns me. Why are we creating a complex-valued CC matrix but then fill it with real-valued CCs?

I fixed other parts of the code so that now there is no conversion required! It all seems to work now.

MatthiasSachs commented 2 years ago

sorry, just pushed a fix for a small bug that I noticed today in the morning.

MatthiasSachs commented 2 years ago

@cortner I think this PR is ready to merge. Do you approve?

cortner commented 2 years ago

I will take a look. We will also need to merge in the latest main branch, or rebase.

MatthiasSachs commented 2 years ago

@cortner: I just merged the upstream/main into this branch.

MatthiasSachs commented 2 years ago

Argghh... it's not obvious to me how to fix the newly appearing error after merging with upstream/main!

cortner commented 2 years ago

I wil try - no promises about time-scale

cortner commented 2 years ago

what error? I only get a warning?

MatthiasSachs commented 2 years ago

This one:

ERROR: LoadError: MethodError: namedtuple(::Tuple{}) is ambiguous. Candidates:
  namedtuple(names::Tuple{Vararg{Symbol, N}}) where N in NamedTupleTools at /Users/msachs2/.julia/packages/NamedTupleTools/GZxRn/src/NamedTupleTools.jl:196
  namedtuple(names::Tuple{Vararg{String, N}}) where N in NamedTupleTools at /Users/msachs2/.julia/packages/NamedTupleTools/GZxRn/src/NamedTupleTools.jl:198
Possible fix, define
  namedtuple(::Tuple{})
Stacktrace:
  [1] select(nt::NamedTuple{(:be, :n, :l, :m), Tuple{Symbol, Int64, Int64, Int64}}, ks::Tuple{})
    @ NamedTupleTools ~/.julia/packages/NamedTupleTools/GZxRn/src/NamedTupleTools.jl:270
  [2] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [3] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [4] getindex
    @ ./broadcast.jl:575 [inlined]
  [5] copy
    @ ./broadcast.jl:922 [inlined]
  [6] materialize
    @ ./broadcast.jl:883 [inlined]
  [7] _sparsify_component!(basis1p::ACE.EllipsoidBondEnvelope{Float64}, keep::Vector{NamedTuple{(:be, :n, :l, :m), Tuple{Symbol, Int64, Int64, Int64}}})
    @ ACE ~/.julia/dev/ACE/src/product_1pbasis.jl:297
  [8] sparsify!(basis1p::ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, keep::Vector{NamedTuple})
    @ ACE ~/.julia/dev/ACE/src/product_1pbasis.jl:277
  [9] clean_1pbasis!(basis::PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64})
    @ ACE ~/.julia/dev/ACE/src/pibasis.jl:221
 [10] clean_pibasis!(basis::ACE.SymmetricBasis{PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64}, ACE.Invariant{Float64}, ACE.O3{:l, :m}, typeof(real), ACE.Invariant{Float64}}; atol::Float64)
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:250
 [11] clean_pibasis!(basis::ACE.SymmetricBasis{PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64}, ACE.Invariant{Float64}, ACE.O3{:l, :m}, typeof(real), ACE.Invariant{Float64}})
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:243
 [12] ACE.SymmetricBasis(φ::ACE.Invariant{Float64}, symgrp::ACE.O3{:l, :m}, pibasis::PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64}, _real::Function)
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:176
 [13] #SymmetricBasis#313
    @ ~/.julia/dev/ACE/src/symmbasis.jl:103 [inlined]
 [14] ACE.SymmetricBasis(φ::ACE.Invariant{Float64}, basis1p::ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, symgrp::ACE.O3{:l, :m}, Bsel::ACE.CategorySparseBasis; isreal::Bool, kwargs::Base.Iterators.Pairs{Symbol, ACE.Utils.var"#5#7", Tuple{Symbol}, NamedTuple{(:filterfun,), Tuple{ACE.Utils.var"#5#7"}}})
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:96
 [15] #SymmetricBasis#310
    @ ~/.julia/dev/ACE/src/symmbasis.jl:85 [inlined]
 [16] SymmetricBond_basis(ϕ::ACE.Invariant{Float64}, env::ACE.EllipsoidBondEnvelope{Float64}, Bsel::SparseBasis; RnYlm::Nothing, bondsymmetry::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ACE.Utils ~/.julia/dev/ACE/src/utils.jl:97
 [17] SymmetricBond_basis(ϕ::ACE.Invariant{Float64}, env::ACE.EllipsoidBondEnvelope{Float64}, Bsel::SparseBasis)
    @ ACE.Utils ~/.julia/dev/ACE/src/utils.jl:78
 [18] top-level scope
    @ ~/.julia/dev/ACE/test/test_bondbasisselectors.jl:38
in expression starting at /Users/msachs2/.julia/dev/ACE/test/test_bondbasisselectors.jl:36

This is the error I get on my own computer. But I presume that the error encountered in 4b82aec is the same.

cortner commented 2 years ago

ok, I can reproduce

cortner commented 2 years ago

Ok - the problem is here:

env = ACE.EllipsoidBondEnvelope(r0cut, rcut, zcut;floppy=false, λ= .5)
property = Invariant() ... 
basis = ACE.Utils.SymmetricBond_basis(property, env, Bsel; )

you are trying to create a symmetric basis but it isn't a basis, so the basis spec contains only empty tuples. That's a weird edge case and I'm not even sure I want to allow this... Is this every going to occur outside of a test?

cortner commented 2 years ago

See here : once we sparsity (in this case just "clean up") ...

function _sparsify_component!(basis1p, keep)
   # get rid of all info we don't need 
   syms = symbols(basis1p)
   keep1 = unique( select.(keep, Ref(syms)) )
   .... 

I can't "select" a subset of the basis specs because there is no spec. It is at odds with the entire design of the basis framework.

MatthiasSachs commented 2 years ago

Ok - the problem is here:

env = ACE.EllipsoidBondEnvelope(r0cut, rcut, zcut;floppy=false, λ= .5)
property = Invariant() ... 
basis = ACE.Utils.SymmetricBond_basis(property, env, Bsel; )

you are trying to create a symmetric basis but it isn't a basis, so the basis spec contains only empty tuples. That's a weird edge case and I'm not even sure I want to allow this... Is this every going to occur outside of a test?

I can see in the error message that the spec is empty but I don't yet understand why this is the case? I create a special basis selector within the call of SymmetricBond_basis and then call a constructor for SymmetricBasis using this basis selector. Why is what I am instantiating not a basis? Is the bond envelope env and the fact that we don't associate a symbol with that causing problems here? So far this is the standard way for me to create bond bases...

cortner commented 2 years ago

It’s just a multiplier by itself - it isn’t multiplying a basis

cortner commented 2 years ago

Any thoughts on this? Maybe there is a way to just change the tests? Or do you need me to figure out a way around this? I'M just a bit overwhelmed right now - sorry for the delay.