ACEsuit / ACE.jl

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

Bug / instability of inner cutoff #118

Open davkovacs opened 2 years ago

davkovacs commented 2 years ago

I am trying to fit ACE models using ACEatoms, and when I plot the dimer curves of the potentials I get the 1-body and 2-body parts behave as expected, but the ACE component of the basis does not go to 0 at the inner cutoff, but it explodes: image Code used to define the basis:

function create_basis(maxdeg, maxorder, maxlevels, symmetry, rin, rcut)
    """
    Args:
    maxdeg::Int 
    maxorder::Int: maximum correlation order
    maxlevels::Dict maximum degree for each cor order
    symmetry::string either \"invariant\" or \"vector\"
    rin::Float Inner cutoff
    rcut::Float Outer cutoff

    Returns:
    Bsite
    Vsite
    """
    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 = rin, rcut = rcut, r0 = 3.0)

    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
davkovacs commented 2 years ago

Turns out the problem was that I didn't pass pin=2 to the basis construction.

I think we should have a better behavior than this if someone specifies an inner cutoff, but doesn't specify pin.

cortner commented 2 years ago

There is always an rin even if you don’t choose it because it specifies the domain of orthogonality for the polynomials

maybe it could be better documented but if don’t see how to change the behaviour

davkovacs commented 2 years ago

I did choose rin, and that makes perfect sense that you always have it. The problem was the not specifying pin=2 lead to 10^11 or -10^11, so completelycrazy behaviour of ACE where it shouldn’t have been defined.

cortner commented 2 years ago

So we could consider making pin = nothing the current pin = 0 and pin=0 becomes a jump to zero. But this isn’t going to solve anything. You still have to specify pin to make sure you get the desired behaviour.

Maybe this could be resolved by REQUIRING pin once rin is specified.

davkovacs commented 2 years ago

Yes, requiring pin might be the best solution. It shouldn’t happen to anyone ever again, that you have rin Specified, but the ACE is nonzero outside it.