SparseBasis basis selector bug?

[davkovacs]

I noticed that the maxlevels directory gets ignored in the SparseBasis Basis selector. I am using it for the covariant basis, and am finding that the code below returns the same set of basis functions for the two different maxlevels directories.

using ACE, ACEatoms
using ACE: SymmetricBasis, evaluate, SimpleSparseBasis

maxdeg = 6
maxorder = 3

maxlevels = Dict(1 => 6,
                 2 => 6,
                 3 => 0)
# maxlevels = Dict(1 => 6,
#                  2 => 6,
#                  3 => 6)

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.6, rcut = 8.0, r0 = 3.0)

ACE.init1pspec!(B1p, Bsel)
basis = SymmetricBasis(ACE.EuclideanVector(Float64), B1p, Bsel)

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);
@show length(Bsite);

I am on the dipole-2 branch of ACEatoms and on the 0.12.24 version of ACE. I tried using the most recent ACE main branch, but that has conflicts with ACEatoms, that is why I reverted to 0.12.24. (I think James created an issue for that in ACEatoms already. )

@cortner could you please verify if the above is expected behaviour (and give some explanation) or a bug indeed. @JPDarby have you observed this?

[JPDarby]

I just tried this with ACE 0.12.22 and the main branch of ACEatoms. I see the same length with both maxlevels if I run David's code (for both EuclideanVector AND Invariant) which seems suspicious... If I change the p-norm to 2 and make the differences more extreme I do see some difference (see below).

using ACE, ACEatoms
using ACE: SymmetricBasis, evaluate, SimpleSparseBasis

maxdeg = 6
maxorder = 3

maxlevels = Dict(1 => 3,
                 2 => 3,
                 3 => 4)
maxlevels2 = Dict(1 => 3,
                  2 => 5,
                  3 => 8)

weights = Dict(:l => 1.5, :n => 1.0)
Bsel = SparseBasis(maxorder=maxorder, weight=weights, maxdegs=maxlevels, p=2.0);
Bsel2 = SparseBasis(maxorder=maxorder, weight=weights, maxdegs=maxlevels2, p=2.0);

B1p = ACEatoms.ZμRnYlm_1pbasis(; species = [:O, :H], maxdeg=maxdeg,
                               Bsel=Bsel, rin = 0.6, rcut = 8.0, r0 = 3.0)

ACE.init1pspec!(B1p, Bsel)
basis = SymmetricBasis(ACE.EuclideanVector(Float64), B1p, Bsel)
basis2 = SymmetricBasis(ACE.EuclideanVector(Float64), B1p, Bsel2)

println("basis lengths are $(length(basis)) and $(length(basis2))")

basis = SymmetricBasis(ACE.Invariant(), B1p, Bsel)
basis2 = SymmetricBasis(ACE.Invariant(), B1p, Bsel2)

println("basis lengths are $(length(basis)) and $(length(basis2))")

[cortner]

Can we first get rid of the ACEatoms dependency please. Could one of you try to reproduce this on ACE 0.12.25, with just ACE (and no ACE atoms). Then try to reproduce it without the species as I don't see why it should be relevant.

[davkovacs]

@cortner Here is a minimal example still displaying the bug on the most recent main branch of ACE (v0.12.25) without any other dependency.

using ACE
using ACE: SymmetricBasis, SimpleSparseBasis

maxlevels1 = Dict(1 => 6,
                 2 => 6,
                 3 => 0)
maxlevels2 = Dict(1 => 6,
                 2 => 6,
                 3 => 6)

maxdeg = 6
maxorder = 3
weights = Dict(:l => 1.5, :n => 1.0)

Bsel1 = SparseBasis(maxorder=maxorder, weight=weights, maxdegs=maxlevels1, p=1.0);
B1p1 = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, maxL=maxdeg, Bsel=Bsel1, rin=0.6, rcut=8.0, r0=3.0)

Bsel2 = SparseBasis(maxorder=maxorder, weight=weights, maxdegs=maxlevels2, p=1.0);
B1p2 = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, maxL=maxdeg, Bsel=Bsel2, rin=0.6, rcut=8.0, r0=3.0)

ACE.init1pspec!(B1p1, Bsel1)
basis1 = SymmetricBasis(ACE.EuclideanVector(Float64), B1p1, Bsel1)

ACE.init1pspec!(B1p2, Bsel2)
basis2 = SymmetricBasis(ACE.EuclideanVector(Float64), B1p2, Bsel2)

@show length(basis1)
@show length(basis2);

[JPDarby]

I think the bug is in line 67 in pibasis.jl

admissible = bb -> (level(bb, Bsel, basis1p) <= maxlevel(Bsel, basis1p))

maxlevel(Bsel, basis1p)) returns the largest level in Bsel.maxlevels but it needs to depend on the correlation order of bb.

So in David's example maxlevel=6 for both cases whereas for my test it was 4 and 8, hence the difference

[cortner]

you are right, this should be

admissible = bb -> (level(bb, Bsel, basis1p) <= maxlevel(bb, Bsel, basis1p))

[cortner]

I'm fixing and testing it now.

[cortner]

Apparently to run this script we are forced to have a "default" key.

using ACE
using ACE: SymmetricBasis, SimpleSparseBasis

maxlevels1 = Dict("default" => 2, 1 => 6,
                 2 => 6,
                 3 => 0)
maxlevels2 = Dict("default" => 2, 1 => 6,
                 2 => 6,
                 3 => 6)

maxdeg = 6
maxorder = 3
weights = Dict(:l => 1.5, :n => 1.0)

Bsel1 = SparseBasis(maxorder=maxorder, weight=weights, maxdegs=maxlevels1, p=1.0);
B1p1 = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, maxL=maxdeg, Bsel=Bsel1, rin=0.6, rcut=8.0, r0=3.0)

Bsel2 = SparseBasis(maxorder=maxorder, weight=weights, maxdegs=maxlevels2, p=1.0);
B1p2 = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg, maxL=maxdeg, Bsel=Bsel2, rin=0.6, rcut=8.0, r0=3.0)

ACE.init1pspec!(B1p1, Bsel1)
basis1 = SymmetricBasis(ACE.EuclideanVector(Float64), B1p1, Bsel1)

ACE.init1pspec!(B1p2, Bsel2)
basis2 = SymmetricBasis(ACE.EuclideanVector(Float64), B1p2, Bsel2)

@show length(basis1)
@show length(basis2);

Is this the result you expect?

julia> @show length(basis1)
length(basis1) = 23
23

julia> @show length(basis2);
length(basis2) = 52

[cortner]

Word of warning! The ACE basis now has a constant term by default!!! Not for a Euclidean vector, because constants cannot be covariant. But e.g. for Invariants, covariant matrices, etc...

[JPDarby]

Does the Rn1pbasis also contain a constant (probably * cut-off function in some way)? And if so, is that basis function always the first one?

[JPDarby]

Is this the result you expect?

Yes I think so

[cortner]

Does the Rn1pbasis also contain a constant (probably * cut-off function in some way)? And if so, is that basis function always the first one?

Rn has no constant term as such, but it has 1 * fcut. This is not the constant term I mean above. I am referring to a zero-correlation. I.e. A_{} = 1 (empty).

See this test script on how to remove the zero-correlation.

[cortner]

I notice a branch co/level which appears to be dangling and hasn't been merged. Before I close this, I will investigate whether this fixes this and other problems.

[cortner]

now merged, and tagged as 0.12.27. (it's worth an update, I fixed some performance bugs over Christmas)