Proper implementation of scaling function

[cortner]

• PIbasis scaling should reduce to 1pbasis scaling
• 1pbasis scaling should "know"about the type of variables in the state
• make scaling for categorical basis

[davkovacs]

@cortner I just found this issue and was wondering what is going on with scaling in the new ACE? I think I am probably the only one currently using the new ACE for a project with IPfitting and it would be important for me to use scaling. But am getting errors. In the old ACE1 scaling works like this:

s = ACE.scaling(db.basis.BB[2], rlap_scal)

Now this doesn't work any more because the type of the first argument is ::ACEatoms.ACESitePotentialBasis I have figured out that replacing the above with

ACE.scaling(dB.basis.BB[2].models[1], rlap_scal)

is a valid argument for scaling, but I still get the following error:

MethodError: no method matching abs(::AtomicNumber)
Closest candidates are:
  abs(::Complex) at complex.jl:264
  abs(::ForwardDiff.Dual{T, V, N} where {V, N}) where T at /home/cdt1906/.julia/packages/ForwardDiff/PBzup/src/dual.jl:206
  abs(::T) where T<:GeometryBasics.OffsetInteger at /home/cdt1906/.julia/packages/GeometryBasics/WMp6v/src/offsetintegers.jl:58
  ...

Stacktrace:
  [1] (::ACE.var"#_absvaluep#270"{Float64})(x::AtomicNumber)
    @ ACE ~/.julia/dev/ACE/src/pibasis.jl:201
  [2] MappingRF
    @ ./reduce.jl:93 [inlined]
  [3] _foldl_impl(op::Base.MappingRF{ACE.var"#_absvaluep#270"{Float64}, Base.BottomRF{typeof(Base.add_sum)}}, init::Base._InitialValue, itr::NamedTuple{(:mu, :n, :l, :m), Tuple{AtomicNumber, Int64, Int64, Int64}})
    @ Base ./reduce.jl:58
  [4] foldl_impl(op::Base.MappingRF{ACE.var"#_absvaluep#270"{Float64}, Base.BottomRF{typeof(Base.add_sum)}}, nt::Base._InitialValue, itr::NamedTuple{(:mu, :n, :l, :m), Tuple{AtomicNumber, Int64, Int64, Int64}})
    @ Base ./reduce.jl:48
  [5] mapfoldl_impl(f::ACE.var"#_absvaluep#270"{Float64}, op::typeof(Base.add_sum), nt::Base._InitialValue, itr::NamedTuple{(:mu, :n, :l, :m), Tuple{AtomicNumber, Int64, Int64, Int64}})
    @ Base ./reduce.jl:44
  [6] mapfoldl(f::Function, op::Function, itr::NamedTuple{(:mu, :n, :l, :m), Tuple{AtomicNumber, Int64, Int64, Int64}}; init::Base._InitialValue)
    @ Base ./reduce.jl:160
  [7] mapfoldl
    @ ./reduce.jl:160 [inlined]
  [8] #mapreduce#218
    @ ./reduce.jl:287 [inlined]
  [9] mapreduce
    @ ./reduce.jl:287 [inlined]
 [10] #sum#221
    @ ./reduce.jl:501 [inlined]
 [11] sum(f::Function, a::NamedTuple{(:mu, :n, :l, :m), Tuple{AtomicNumber, Int64, Int64, Int64}})
    @ Base ./reduce.jl:501
 [12] scaling(pibasis::PIBasis{ACE.Product1pBasis{3, Tuple{Species1PBasis{2}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{StaticArrays.SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{StaticArrays.SVector{3, ComplexF64}}}}}}, ComplexF64}, typeof(real), Float64, ComplexF64}, p::Float64)
    @ ACE ~/.julia/dev/ACE/src/pibasis.jl:208
 [13] scaling(basis::SymmetricBasis{PIBasis{ACE.Product1pBasis{3, Tuple{Species1PBasis{2}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{StaticArrays.SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{StaticArrays.SVector{3, ComplexF64}}}}}}, ComplexF64}, typeof(real), Float64, ComplexF64}, ACE.Invariant{Float64}, ACE.O3{:l, :m}, typeof(real), ACE.Invariant{Float64}}, p::Float64)
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:322
 [14] top-level scope
    @ In[27]:1
 [15] eval
    @ ./boot.jl:360 [inlined]
 [16] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1094

Do you happen to know what is going on? And is there perhaps a 1 line fix that would allow me to keep using IPfitting and scaling to fit ACE2 potentials. I am using the dipole branch of IPfitting for this project.

[davkovacs]

@cortner if you had time could you please give me some hints how to fix this? Or whether it requires more substantial work?

[cortner]

sorry I missed this and thanks for the reminder. It should work. Can you post a short script so I can reproduce the error? (and which packages and package versions???)

[davkovacs]

Here is a little code that I would like to fix, but currently gives and error on the main branch of ACE, ACEatoms.

Basically the scaling is not implemented for some of the functions in ACE (and is not implemented for ACEatoms.ACESitePotentialBasis that should just call the not implemented functions in ACE)

using ACE, ACEatoms
using ACE: SymmetricBasis

function create_basis(maxdeg, maxorder, maxlevels, symmetry)
    """
    symmetry is either "invariant" or "vector"
    """
    weights = Dict(:l => 1.5, :n => 1.0)
    Bsel = SparseBasis(maxorder=maxorder, weight=weights, maxdegs=maxlevels, p=1.0);

    B1p = ACEatoms.ZμRnYlm_1pbasis(; species = [:O, :H], maxdeg=maxdeg,
                                   Bsel=Bsel, rin = 0.65, rcut = 4.5, r0 = 3.0,
                                    pin = 2)

    ACE.init1pspec!(B1p, Bsel)
    if symmetry == "vector"
        basis = SymmetricBasis(ACE.EuclideanVector(Float64), B1p, Bsel)
    elseif symmetry == "invariant"
        basis = SymmetricBasis(ACE.Invariant(), B1p, Bsel)
    else
        @warn "Unkown symmetry, using invariant or vector"
        return nothing
    end

    cH = randn(length(basis))
    cO = randn(length(basis))

    models = Dict(:H => ACE.LinearACEModel(basis, cH),
                  :O => ACE.LinearACEModel(basis, cO))
    Vsite = ACEatoms.ACESitePotential(models)

    Bsite = ACEatoms.basis(Vsite);
    println("Creating basis of length ", length(Bsite))
    return Bsite, Vsite
end

maxdeg = 8
maxorder = 2
maxlevels = Dict("default" => 8,
                 1 => 8,
                 2 => 6)
Bsite2, Vsite2 = create_basis(maxdeg, maxorder, maxlevels, "invariant");

s = ACE.scaling(Bsite2, 1.0)

Bsite2, Vsite2 = create_basis(maxdeg, maxorder, maxlevels, "vector");

s = ACE.scaling(Bsite2, 1.0)

[cortner]

I'm trying to reproduce with latest ACE and running into the weirdest bug here - I'll need to sort this out and then will come back to the scaling issue.

[cortner]

ok - so back to @davkovacs 's problem -- the issue is that in ACE.jl I assumed the species basis will use symbols, but it actually uses atomic number. I will put something mildly hacky together now, and we leave this issue open for the future.

[cortner]

@davkovacs --- I've now made ACEatoms and ACE interoperable on the implementation of the scaling functions. But I have not yet implemented the scaling function for a ACEatoms.ACESitePotentialBasis. This should be relatively straightforward to do now by putting together the scalings for the individual basis components e.g.

scal_H = ACE.scaling(Bsite2.models[1], 1)
scal_O = ACE.scaling(Bsite2.models[8], 1)

and then putting them together in the right fashion, whatever the ordering of the basis is? But let me know if it isn't clear and I can get back to it.

I've tagged ACEatoms 0.0.11 and ACE 0.12.32, which you need to use now.

[davkovacs]

I have created a pull request implementing this to ACEatoms.