Species dependent transforms for v0.8.x

[cortner]

cf #77

@davkovacs @LarsSchaaf you can follow along here if useful.

[cortner]

@davkovacs @LarsSchaaf this is ready to be tested. I don't really want to merge it yet though until you tell me it behaves in the way you hope/expect. Please look through test/polynomials/test_multitrans.jl to see how to use the functionality and ask questions here?

[davkovacs]

Thank @cortner I had a look, just to clarify what the parameters mean:

(:Fe, :C) => PolyTransform(2, d)

means that the Fe - C distance (irrespective of whether it is an Fe or a C center will be transformed based on the PolyTransform(2, d) which is (1 + d) / (1+x)^2. Here we should set d to be the minimum distance as in the example? I would intuitively set it to the typical bond length rather than the minimum of the training set.

[cortner]

d is not the minimum bond-length, it is the nearest-neighbour bond-length. Around the NN distance is where most of the action happens, where you need the best resolution. That's how you should choose d. But you should set it to anything you like and then see what happens. I just made some choice for the sake of having an example to test.

irrespective of whether it is an Fe or a C center

it will be symmetric unless you also specify (:C, :Fe) => ... see last line in the example.

[davkovacs]

Makes sense, thank you! I think this will be a big help when there is a dataset with large variation in the sizes of atoms.

[cortner]

So merge and tag as is, or wait until you have some tests?

[davkovacs]

I am quite busy at the moment and don’t have a system where I would immediately need it (no metals and hydrogens etc. ) Maybe Lars could try it, and i am also happy to try it next week, I can probably come up with some tests.

[cortner]

@LarsSchaaf This should now provide what you need (in principle!) - I'd be grateful if you can test it. I don't have the time to write very sophisticated tests covering all corner cases, so just report back any problems you encounter.

[cortner]

Usage:

transforms = Dict(
     (:Fe, :C) => PolyTransform(2, (rnn(:Fe)+rnn(:C)) / 2), 
     (:C, :Al) => PolyTransform(2, (rnn(:Al)+rnn(:C)) / 2), 
     (:Fe, :Al) => PolyTransform(2, (rnn(:Al)+rnn(:Fe)) / 2), 
     (:Fe, :Fe) => PolyTransform(2, rnn(:Fe)), 
     (:Al, :Al) => PolyTransform(2, rnn(:Al)), 
     (:C, :C) => PolyTransform(2, rnn(:C)), 
     (:Al, :Fe) => PolyTransform(2, (rnn(:Al)+rnn(:Fe)) / 2 + 1)  
   )

# define cutoffs as (rin, rcut) pairs. If (S1, S2) and (S2, S1) are both 
# specified then the cutoff is non-symmetric, If only one is specified, then 
# it will be symmetric 
cutoffs = Dict(
   (:Fe, :C) => (1.5, 5.0), 
   (:C, :Al) => (0.7, 6.0), 
   (:Fe, :Al) => (2.2, 4.5), 
   (:Fe, :Fe) => (2.0, 5.0), 
   (:Al, :Al) => (2.0, 5.0), 
   (:C, :C) => (1.5, 5.2), 
   (:Al, :Fe) => (1.5, 5.0)  )

trans = multitransform(transforms, cutoffs=cutoffs)
Rn = transformed_jacobi(maxdeg, trans; pcut = 2)

Note that the calling convention for transformed_jacobi has now changed - but only if you construct it with a MultiTransform - since that transform already contains the cutoff information.

I look forward to seeing your results!

[WillBaldwin0]

I believe that the @assert statements in this section of code are left over from an earlier version of the multitransform:

function transform(t::MultiTransform{NZ, TT}, r::Number) where {NZ, TT}
   @assert (TT <: AffineT) "transform(::MultiTransfrom, r) is only defined if rin, rcut are specified during construction"
   x = sum( transform(_t, r)  for _t in t.transforms ) / NZ^2 
   @assert (abs(abs(x) - 1) <= 1e-7) "transform(::MultiTransfrom, r) is only defined for r = rin, rcut"
   return x  
end

should this be removed or modified before merging?

[davkovacs]

@WillBaldwin0 have you ended up implementing the exponential distance transform?

[cortner]

@WillBaldwin0 -- thanks again for flagging this. I'll try to look into it today.

Re exponential distance transform - this is implemented in ACE. Or did you just mean whether you've tried to use it?

[cortner]

Also ... I'm still very fond of the AgnesiTransform - anybody want to take that on?

[WillBaldwin0]

I will probably be using the Agnesi transform. Today we were looking at an exponential scaling which would scale with r_in and r_out, so that the free parameter would be transferable between different pairs of atoms with different inner and outer cutoffs.

I will have a look at the exponential transform in ACE too.

[cortner]

Better to extend the transforms in ACE than create a new one. If you file an issue with the functional form you want, I can probably put it together in 15-30 minutes.

[LarsSchaaf]

I was wondering if @davkovacs @WillBaldwin0, you managed to get the pair potential to work? The transform works perfectly for the rip_basis.

[LarsSchaaf]

Concerning Preliminary Results, I'm not finding any advantage in using the altered transforms currently. In the plot below the RMSEs remain lowest for the species independent setting (Constant). This is without the PairPotential, so these RMSE's need to be taken with a pinch of salt. I'm also dealing with a very small data set.

[LarsSchaaf]

In the above, I compare setting the cutoff and transform species independent (constant), setting both based on the covalent bond radii of the involved species (r_in = 0.9*cbr, r0 = cbr) and a mixture of both (r_in set as previously, r0 set to a constant).

[cortner]

Have you tried keeping the outer cutoff the same for all, and only adjusting the inner cutoff? TBH, I am not sure that I would expect better accuracy with the species-dependent transforms; the reason is that at long distances the space is compressed a lot anyhow and in transformed coordinates increasing cutoff makes almost no difference.

But I could imagine that the inner cutoff can make a difference, not with accuracy but with extrapolation quality.

[cortner]

@WillBaldwin0 - I now think this method should be removed entirely and have done so. If you run into problems let me know.

[WillBaldwin0]

Hi @cortner, I'm trying to get a species dependent pair potential working. I've set up a multitransform based RPIbasis by doing

trans = multitransform(transforms, cutoffs=cutoffs)
Pr = transformed_jacobi(get_maxn(Deg, maxdeg, species), trans; pcut = 2)
P1 = BasicPSH1pBasis(Pr; species = species, D = Deg)
Bsite = RPIBasis(P1, N, Deg, maxdeg)

and this works fine for fitting. It also works fine when I combine it with a non-species dependent pair potential:

Bpair = pair_basis(
    species=species,
    r0=1.0,
    maxdeg=deg_pair,
    rcut=r_cut_pair,
    rin=0.0,
    pin=0,
)

Bpair = PolyPairBasis(pair_rbasis, species)
B = JuLIP.MLIPs.IPSuperBasis([Bpair, Bsite]);

However trying to make a species dependent pair potential seems to have some problems. Building a basis as follows:

pair_trans = multitransform(transforms, cutoffs=cutoffs_pair)
pair_rbasis = transformed_jacobi(deg_pair, pair_trans; pcut = 2, pin = 0)
Bpair = PolyPairBasis(pair_rbasis, species)
B = JuLIP.MLIPs.IPSuperBasis([Bpair, Bsite]);

I find that the basis can be created, but when trying to assemble the LsqDB I get the error

MethodError: no method matching transform_d(::ACE.Transforms.MultiTransform{3, ACE.Transforms.AffineT{Float64, AgnesiTransform{Float64, Int64}}}, ::Float64)
Closest candidates are:
  transform_d(::ACE.Transforms.MultiTransform, ::Number, ::AtomicNumber, ::AtomicNumber) 

Is this a bug or should a species dependent pair potential be specified differently?

(p.s. to get this to run this far I also had to remove the assert statements that I mentioned in an earlier post)

[cortner]

this is a bug - I just forgot about it, thank you for reporting it. looking at it now.

[cortner]

@LarsSchaaf there is a chance your results are also due to a bug. Let's not despair yet. EDIT: probably not after all. The many-body ACE basis and evaluator do look ok.

[cortner]

@WillBaldwin0 Would you be willing to try again?

[cortner]

analytic transform usage:

trans = ACE.Transforms.AnalyticTransform(
            "r -> exp(- 1.234 * r^2)", "x -> sqrt(- log(x) / 1.234)")

[WillBaldwin0]

I've ran the same code and it seems to work fine, thanks!

[cortner]

To everyone in this discussion - I’ll f you have sensible unit tests to contribute please post them here

[cortner]

is this ready for merging - i.e. does this now have the functionality you need? (we can also add bug-fixes and/or tweaks later of course...)

[LarsSchaaf]

I think it is. Thanks for the analytic transform. I think it's a great addition to the simplicity of testing different transforms going forth.

The pair potential also works for me now and I'm starting to see some advantages using the species dependant cutoffs. I just ran a quick test, where I only changed the inner cutoffs based on the RDF. Specifically, I set the inner cutoffs, such that there are 4 data points below the cutoff for each species pair (to help the pair potential remain repulsive). Looking forward to seeing what we can get out of this added functionality.

[cortner]

Wow, those don't seem to be small advantages. Not bad. I'd be interested to have my intuition confirmed or rejected that adjusting the inner cutoff will give you more than the outer cutoff. (in the sense that you can always just increase the outer cutoff without much penalty)

[LarsSchaaf]

I assume this will heavily depend on the functional form of the transform. If we continue using transforms, where the gradient approaches zero near the outer cutoff, then I would guess the penalty of increasing the outer cutoff will be minimal as you say. However, I feel like we should use a transform whose gradient, by design can never be less than a certain value, with respect to r_out - r_in (I can post more about this later, once I've run some tests). If the gradient near the cutoff is not negligible, I would then expect varying the outer cutoff to have a bigger effect on how well we resolve around our focus distance r0.

[cortner]

That’s interesting - I haven’t really thought about this. It’s another thing against my intuition but I can imagine it may really help in some contexts. I’ll be curious to learn more

[cortner]

this is now released as v0.8.9, please put this version bound into your Project.toml.