Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Proper implementation of scaling function #68

Open
cortner opened this issue Oct 25, 2021 · 8 comments
Open

Proper implementation of scaling function #68

cortner opened this issue Oct 25, 2021 · 8 comments

Comments

@cortner
Copy link
Member

cortner commented Oct 25, 2021

  • 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
Copy link

davkovacs commented Feb 8, 2022

@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
Copy link

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

@cortner
Copy link
Member Author

cortner commented Feb 22, 2022

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
Copy link

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
Copy link
Member Author

cortner commented Feb 22, 2022

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
Copy link
Member Author

cortner commented Feb 22, 2022

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
Copy link
Member Author

cortner commented Feb 22, 2022

@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
Copy link

I have created a pull request implementing this to ACEatoms.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants